A C# code for solving 3D topology optimization problems using SAP2000

SAP2000 is well-known commercial software for analysis and design of structural systems that is equipped with an open application programming interface (OAPI). In this work, a code written in C# able to solve three-dimensional topology optimization problems is presented, where a topology optimization framework was integrated into SAP2000 taking advantage of its OAPI feature. The code is partially based on the 99 and 88 line codes written by Sigmund (Struct Multidiscip Optim 21(2):120–127, 2001) and Andreassen et al. (Struct Multidiscip Optim 43(1):1–16, 2011). The code solves the problem of minimum compliance while through OAPI it takes advantage of all modelling capabilities of SAP2000. The paper covers the theoretical aspects of topology optimization incorporated in the code and provides detailed description of their numerical implementations. Special effort was made to the latter one, describing in detail all numerical aspects of the code, in order to facilitate the reader to understand the code, and therefore being able to further enhance its capabilities. The complete code can be downloaded from GitHub (https://github.com/nikoslagaros/TOCP).


Introduction
Since 1970 structural optimization has been the subject of intensive research and several methodologies for achieving optimized structural designs have been advocated (Gallagher and Zienkiewicz 1973;Haug and Arora 1974;Moses 1974;Sheu and Prager 1968;Spunt 1971).Topology optimization represents a material distribution numerical procedure for synthesizing structural layouts without any preconceived shape.In this study we present the C# code that was developed in the framework of the research project ''OPTARCH: Optimization Driven Architectural Design of Structures'' (H2020-MSCA-RISE-2015) aiming to facilitate the needs of the project.In the rest of the paper this code is abbreviated as TOCP standing for topology optimization computing platform.One of the objectives of the project is to integrate architectural design criteria into shape and topology (S&T) optimization procedures.In this direction, a first goal of the OPTARCH project is to exploit the use of S&T optimization techniques in computer aided architectural design.This scope, among others, will be fulfilled by formulating S&T problems as well as developing procedures for dealing with.Studies where criteria imposed by architects are integrated into topology optimization are rather few (Beghini et al. 2014).According to OPTARCH project, it is anticipated that specially tailored topology optimization design tools will equip the engineer with the capacity to directly define various design alternatives.As a result, it is of importance to develop a platform that will provide a common ground with computer aided design (CAD) tools both for defining the architectural design constraints and for depicting the optimized results.In this direction, TOCP code was developed to become the basis and to serve OPTARCH objectives.Moreover, exchange of ideas is considered as an issue of major importance and therefore the gold model for open access was considered as an essential pillar for better communication and dissemination of OPTARCH's results.Consequently, it was decided that students and newcomers to the field of topology optimization to be able to download TOCP code from GitHub (https://github.com/nikoslagaros/TOCP).During the four years of the project, the code will further be enhanced and new capabilities incorporated will also be provided for free through the above mentioned repository.Among different S&T optimization techniques that exist in the literature, the solid isotropic material with penalization (SIMP) (Bendsøe and Sigmund 1999;Bendsøe and Sigmund 2003) and level-set (Allaire 2002;Allaire et al. 2004) methods are the most popular ones.Although, in the framework of OPTARCH project both methods will be examined with respect to their suitability for the needs of conceptual design for architectural criteria, in the current version of TOCP only SIMP was implemented.
In literature, various free source codes are provided, among others by the pioneering researchers in the field of S&T optimization, like the Matlab codes by Sigmund (2001), Andreassen et al. (2011) and the Scilab ones by Allaire (2017).The applicability of most of these codes is limited to 2D elastic structural systems.Indicatively, the two Matlab codes written by Sigmund (2001) and Andreassen et al. (2011) both for solving 2D topology optimization problems based on SIMP.Allaire and his co-workers released a set of Scilab routines (Allaire 2017) that perform S&T optimization for 2D elastic structures based on the level set method.Talischi et al. (2012) presented a 2D Matlab code for structural topology optimization that includes a general finite element routine based on isoparametric polygonal elements, which can be considered as an extension of linear triangles and bilinear quads.As far as 3D topology optimization problems only limited number of codes can be found, like the 169-line Matlab code by Liu and Tovar (2014) where SIMP is used and the 100-line Python code by Zuo and Xie (2015) that is based on the bidirectional evolutionary structural optimization method.All these codes are limited to implementations of S&T optimization using interpret programming languages (like Matlab, Python and Scilab) while the solution of the system of equations resulted from the finite element (FE) modelling is also part of these codes.More exceptional geometries and structural systems are not easy to deal with these codes, while only linear analysis is support by most of these codes.There are several commercial software applications of topology optimization; however all of them are integrated into mechanical and aerospace engineering mostly oriented analysis and design software, like the OpiStruct by Altair (OptiStruct 2017), the Tosca optimization suite by Dassault Systems (2017) that works with standard FEA software (like ABAQUS and MSC Nastran), Ansys topology optimization solution (2017) and recently by SOLIDWORKS (2017).In this work a C# source code that is a general-purpose, object-oriented programming language, is presented.C# language was developed by Microsoft within its.NET initiative and standardized by ECMA (ECMA-334) and ISO (ISO/IEC 23270:2006), however it is not limited to just Microsoft platforms; C# compilers exist for FreeBSD, Linux and Macintosh computing platforms.C# language was chosen because it is a multi-paradigm programming language encompassing strong typing, imperative, declarative, functional, generic, object-oriented (class-based), and component-oriented programming disciplines.The contribution of the C# code presented herein and the major differences over those presented so far, originates from the additional capabilities and implementation features that are required to be developed in order to make feasible the integration of a topology optimization framework into a civil engineering analysis and design software like SAP2000.In particular, the calculation of the stiffness matrix derivatives are bypassed, generation of the FE mesh, etc.
The basic advantage of the C# code presented here in, is that it provides a fully functional integration of the S&T optimization concept into a commercial structural and analysis software, specially tailored for structural engineering purposes.Therefore, the users of TOCP code will be able to benefit from all modelling and analysis-design capabilities that commercial software provide.SAP2000 (Wilson and Habibullah 2017), that was introduced over 30 years ago, is equipped with an advanced and rather simple to use open application programming interface (OAPI) making ease the communication of external developers with the capabilities of the software.In addition, a second reason for choosing SAP2000 (Wilson and Habibullah 2017) is because it is available for academic purposes almost free of charge.
The rest of the paper is organized as follows: In order to make the manuscript readable by unfamiliar to S&T optimization researchers, a rather extensive theoretical background regarding the optimization problem formulation and the solution methods implemented in TOCP code, are presented in Sect. 2. Section 3 discusses the general framework of SAP2000 OAPI implementation in order to facilitate future developers.Section 4 is devoted to the description of TOCP code, component by component, representing the numerical implementation of the theoretical parts of Sect. 2. Finally, some simple optimization test examples are presented in order to assist the understanding both of the code and also its integration into SAP2000, as well as how to take advantage of a commercial software capabilities and use those in research practice.

Topology optimization
In order to present the components of TOCP code and make the paper self-contained several theoretical parts of S&T optimization are provided in this section.Structural topology optimization can be considered as a procedure for optimizing the topological arrangement of material into the design domain, eliminating the material volume that is not needed.It can also be seen as the problem of finding the structural layout that best transfers specific loading conditions to supports.It can be employed in order to generate an acceptable initial layout of the structural system, which is then refined by means of shape optimization procedure.Therefore, it can be used to assist designer to define the structural system that satisfies best the operating conditions.Topology optimization procedure operates with gradual ''removal'' of low stressed material.This is treated as typical structural reanalysis problem where small variations are encountered on the stiffness matrix between subsequent optimization steps.

Problem formulation
In topology optimization problem formulations the quantities provided are the domain X, where the optimized layout will be created, the required volume fraction of the optimized layout, the boundary C and loading conditions F. As shown in Fig. 1, the boundary conditions C consist of C o , C s , C t and C u , parts where The surface tractions t are applied at region C t ; regions C s denote the non-optimizable areas; C u represent the support conditions and C o are the geometric boundaries of X.
The problem can be solved using material distribution methods (Bendsøe and Sigmund 1999) for finding the optimum layout of a structural system composed by linearly elastic isotropic material.Therefore, the question under investigation is how to distribute material volume into domain X in order to minimize a specific criterion; compliance C is a commonly used criterion.The distribution of the material volume in domain X is controlled by the density values x distributed over the domain.More specifically, it is controlled by design parameters that are represented by the densities x e assigned to the FE discretization of domain X.The densities (x or x e ) take values in the range [0, 1], where zero denotes no material in the specific element.The mathematical formulation of the topology optimization problem can be expressed as: min x CðxÞ ¼ F T uðxÞ s:t: where C(x) is the compliance of the structure, F is a force vector and uðxÞ is the corresponding global displacement vector.The second expression of Eq. (1) refers to a volume constraint where f is the volume fractions of domain X that the optimized layout should occupy.The third equality of Eq. (1) corresponds to the equilibrium equation where K(x) is the global stiffness matrix of the structural system.The inequality of Eq. (1) denotes the definition set of the density values.

Solid isotropic material with penalization
SIMP method was proposed by Bendsøe and Sigmund (1999) aiming to deal with the problem of Eq. ( 1) and is implemented in TOCP code presented herein.According to SIMP method the finite elements' density values are correlated with the corresponding Young modulus value E, through the following expression: where the parameter p is usually taken equal to 3.This power law correlation is implemented in SIMP to achieve density values closer to the lower and upper bounds of the design variables (i.e.0 and 1).The calculation formula of compliance can be written as follows: Fig. 1 The generalized shape design problem of finding the optimal material distribution in 2D domain CðxÞ ¼ F T uðxÞ ¼ uðxÞ T KðxÞ uðxÞ ð 3Þ Therefore, based on Eq. ( 2), compliance calculation formula can be expressed as: where N is the number of finite elements and u e is the displacement vector in the elements' local coordinate system.Thus, using Eqs.( 2) and (4) the formulation of the topology optimization problem of Eq. ( 1) can now be expressed as follows: min x CðxÞ ¼ uðxÞ T KðxÞ uðxÞ ¼ P N e¼1 ðx e Þ p u T e k 0 e u e s:t: Among others, the optimization problem formulated according to SIMP in Eq. ( 5) can be solved using either the optimality criteria (OC) method or the method of moving asymptotes (MMA), both described in the following sections.

Optimality criteria method
Calculating the derivative of C(x) represents an important part of both OC and MMA algorithms.In order to avoid calculating the derivative of the displacement vector uðxÞ with respect to design vector (x e ) a zero part is subtracted from C(x): Therefore, the partial derivatives of C(x) become: OC is an iterative search algorithm where the solution vector is updated in very iteration until convergence.First a linear approximation of C(x) is defined close to the design variable vector x k : Then C(x) becomes: where since the derivatives of the compliance (C(x)) can take negative values only: Therefore, in order to maximize the subtracting part of the objective function C(x) only the positive part need to be minimized: The following subproblem can now be formulated as follows: min In order to solve this problem the Lagrangian Duality method is applied where the augmented Lagrangian function is expressed as follows: The minimum value of x e resulted from the solution of the subproblem of Eq. ( 15) is obtained minimizing L(x e , k) with respect to x e and maximizing L(x e , k) with respect to k.In order to calculate x e , the derivatives of L(x e , k) with respect to x e are defined first: and therefore the values of x e are obtained: Since x e takes values in the range [0, 1] and large changes should be avoided, x e is updated according to the following rules: m is the maximum change allowed for x e .Similar to x e , in order to calculate k, the derivatives of L(x, k) in respect to k are defined: The calculation of k is achieved by iteratively choosing values of k for each x e until satisfaction of the following equality:

Method of moving asymptotes
Similar to OC, the method of moving asymptotes is an iterative procedure proposed for solving nonlinear optimization problems, in particular of the following form: min x2R N f 0 ðxÞ s:t: f i ðxÞ f ; for i ¼ 1; 2; . ..; n x j x j x j ; for j ¼ 1; 2; . ..; For every iteration k, the values of functions together with their partial derivatives need to be calculated.MMA rely on first order approximations, where the required estimates for design x k are defined using the values of previous steps.Afterwards, Lagrangian Duality method is implemented in order to derive design x k?1 using these approximations.Specifically functions are approximated through the following expression: where and the second order partial derivatives of f i with respect to x j are calculated through the following expressions: where, as it can be observed for Eq. ( 25) they are positive definite, and their values increase when U j and L j , called ''moving asymptotes'', tend to x j .The implementation of MMA into the topology optimization problem of Eq. ( 5) requires some modifications, i.e.only compliance C(x) needs to be approximated, while since its partial derivatives are always negative then p k = 0. Thus, the topology optimization problem of Eq. ( 15) becomes: x e À L ðkÞ e " # s:t: where Lagrangian Duality method is applied also for solving this problem and the augmented Lagrangian function is expressed as follows: The values of x e and k are obtained by minimizing L(x e , k) with respect to x e and maximizing L(x e , k) with respect to k.In order to calculate x e , the derivative of L(x e , k) with respect to x e is defined first: thus the values of x e are obtained as follows: implementing a maximum movement factor m and a, b as move limits, x e is updated according to the following rules: To obtain k, L(x e , k) is maximized with respect to k: This is achieved using the bisection method.The ''moving limits'' are used to avoid division by zero and can be calculated as follows: are used.According to Liu and Tovar (2014) the moving asymptotes are calculated using an update scheme that is based on three successive iterations.In particular, for k = 1 and 2:

SAP2000 OAPI implementation
In computer programming, an application programming interface (API) refers to a set of subroutine definitions, protocols, and tools for building application software.
In general terms, it's a set of clearly defined methods of communication between various software components.A good API makes it easier to develop a computer program by providing all the building blocks, which are then put together by the programming developer.The SAP2000 open application programming interface allows developers to automate many of the processes required to build, analyse, and design models and to obtain customized analysis and design results.It also allows developers to link SAP2000 with third-party tools, providing a path for two-way exchange of model information with other applications.Based on these capabilities provided by SAP2000 OAPI, in this section details on the TOCP C# implementation of minimum compliance topology optimization framework, interfacing with SAP2000, is presented.Although, the code described herein is based on SAP2000 version 15 OAPI documentation, the corresponding code for the latest version 19 is also provided and can be downloaded from the above referenced repository, users are free to use these codes as long as they acknowledge their source.Before describing how parts of topology optimization described previously are modified in order to take advantage of SAP2000 OAPI, the basic issues related to the use of OAPI are provided in this section along with justification for choosing SAP2000 over other also well-known software for analysis and design of structural systems.

Why choosing SAP2000
Over the years, SAP2000 has proven to be one of the most integrated, generalpurpose structural software on engineering practice today.It is considered as one of the easiest, most productive software for structural analysis and design needs (from a simple small 2D static frame analysis to a large complex 3D nonlinear dynamic analysis).The interface of SAP2000 allows creating structural models rapidly and intuitively.Complex models can be generated and meshed using built in templates.This is the reason of less importance we have decided to use the specific software to integrate with the S&T optimization concept.Worth mentioning also that in case of academic purposes needs, multiple licenses for unlimited usage are provided almost for free.
However, the main reason for choosing SAP2000 is due to the fact that through OAPI it's modelling, analysis-design capabilities are fully accessible.By means of its OAPI, integrated design code features can be used to automatically generate wind, wave, bridge, and seismic loads with automatic steel and concrete design code checks per European, US, Canadian, Chinese and other international design standards.In addition, various analytical techniques that allow for step-by-step large deformation analysis, Eigen and Ritz analyses based on stiffness of nonlinear cases, catenary cable analysis, material nonlinear analysis with fibre hinges, buckling analysis, progressive collapse analysis, or plasticity and nonlinear segmental construction analysis along with specialised modelling capabilities using multilayered nonlinear shell element, velocity-dependent dampers, or base isolators, are also accessible through OAPI.

SAP2000 OAPI basic features
In order to access SAP2000 from an external application using its OAPI the following steps need to be followed: The first step is to reference SAP2000 from the application.Specifically, as soon as a new visual studio (VS) C# console application project was created a reference (COM) is added in the project to SAP2000.EXE.Next an instance of the SAP2000 object needs to be created (also known as instantiating the object) within the developed application.This is accomplished as follows: where the command of first line #L1 creates the object variable and the second line #L2 creates an instance of the SAP2000 object.Now that an instance of the SAP2000 object has been created in the application, SAP2000 can be started using the following command: L3 Object.ApplicationStart(SAP2000v15.eUnits.kN_m_C,true,""); At this point an existing model can be opened, or a new one to be created and perform any action required.In general, OAPI commands are accessed through Object.SapModel.It might be helpful to define a SapModel object so that the OAPI commands are accessed through Model instead of Object.SapModel.For example, a new SAP2000 is initiated with the following steps, a Model object is created first: and then the model is initialized as follows: L6 int ret = Model.InitializeNewModel(Sap2000.eUnits.kip_in_F); When finished with a model, it might be needed to close SAP2000 application.This can be accomplished using the following command: L7 Object.ApplicationExit(false); As a last step, Model and Object objects should be set to null.This is accomplished as follows: Setting the objects to null is a very important step.It breaks the connection between developer's application and SAP2000 and frees up system resources.If the objects are not set to null, the SAP2000 application will not completely close (it will still be seen running in Windows task manager).Note that the line numbering of this section (lines are labelled as L1, L2 etc.) is indicative and not related to the numbering of the methods that will be used in the following parts of the study.123

Specific features of SAP2000 OAPI used in topology optimization implementation
Most of the functions of SAP2000 OAPI that are used in this study in order to integrate the topology optimization implementation into SAP2000 are described in this section.GetNameList (ret = Model.SolidObj.GetNameList();), this function retrieves the names of all defined solid objects.SetLocalAxes (ret = Model.SolidObj.SetLocalAxes();), this function sets the local axes angles for solid objects.AddMaterial (ret = Model.PropMaterial.AddMate- rial();) that is a function that adds a new material property to the model based on the pre-defined material properties.GetMPIsotropic (ret = Model.PropMate- rial.GetMPIsotropic();) that is a function that retrieves the mechanical properties for a material with an isotropic directional symmetry type.SetMPIsotropic (ret = Model.PropMaterial.SetMPIsotropic();) that is a function that sets the material directional symmetry type to isotropic, and assigns the isotropic mechanical properties.SetProp (ret = Model.PropSolid.Set- Prop();), this function assigns a solid property to solid objects.ClearSelection (ret = Model.SelectObj.ClearSelection();), this function deselects all objects in the model.CoordinateRange (ret = Model.SelectObj.Coordi- nateRange();), this function selects or deselects objects inside the box defined by the X Min , X Max , Y Min , Y Max , Z Min and Z Max coordinates.GetSelected (ret = Model.SelectObj.GetSelected();), this function retrieves a list of selected objects.SetModelIsLocked (ret = Model.SetModelIsLocked();) this function locks or unlocks the model.SolidJointForce (ret = Model.Re- sults.SolidJointForce();), this function reports the joint forces for the point elements at every corner of the specified solid elements.JointDispl (ret = Model.Results.JointDispl();), this function reports the joint displacements for the specified point elements.The displacements reported by this function are relative displacements.GetProperty (ret = Model.AreaObj.GetProperty();), this function retrieves the area property assigned to an area object.GetType (ret = Model.PropArea.Get- Type();), this function retrieves the property type for the specified area property.Delete (ret = Model.SolidObj.Delete();), this function deletes solid objects.SetRunCaseFlag (ret = Model.Analyze.SetRunCaseFlag();), this function sets the run flag for load cases.RunAnalysis (ret = Model.Analyze.RunAnal- ysis();), this function runs the analysis.The analysis model is automatically created as part of this function.DeselectAllCasesAndCombosForOutput (ret = Model.Results.Setup.DeselectAllCasesAnCombosForOut- put();), this function deselects all load cases and response combinations for output.SetCaseSelectedForOutput (ret = Model.Results.Setup.SetCaseS- electedForOutput();), this function sets a load case selected for output flag.
The current study is accompanied by four visual studio C# projects, all four VS projects correspond to different integrations of the S&T optimization concept into SAP2000 software using its OAPI.The first two projects labelled as TOCP1 and TOCP2 refer to the implementations of OC and MMA algorithms, respectively, and the mesh of the optimizable domain is provided by the user; while the latter two projects labelled as TOCP3 and TOCP4 denote also the implementations of the two algorithms however the mesh of the optimizable domain is generated automatically and the user needs to provide also the representative element's discretization nlx 9 nly 9 nlz.Loading and boundary conditions are denoted by the TOCP user through SAP2000 interface when describing the domain to be optimized.
In this section, the C# code corresponding to the implementation of OC is outlined while description of the numerical implementation of MMA is also provided.Once this code is understood the rest ones can easily be comprehended.The term method that is used in the current and next sections refers to a group of statements that together perform a task.In this section ten methods that are used in OC and MMA visual studio projects are described, together with the methods written for implementing OC and MMA algorithms.All these methods are described in detail below together with a flowchart at the end to better understand their use, while in order to facilitate the understanding of the TOCP code structure the flowchart of Fig. 2 is provided.

ReadModel method
ReadModel refers to the method that is used first by the four projects.Taking advantage of the modelling capabilities of SAP2000, a structural system can be modelled using 1D, 2D and/or 3D finite elements.Using an instance of this method several arrays are returned containing the nodes, the finite elements (solid, frame and area ones) of the model and the load cases considered.Part of the ReadModel method is provided in the code below, where indicatively in line #30 the 3D finite elements (solid ones) used in the model are retrieved using function GetNameList and in line #34 each solid element's local axes angles are set the same as the global ones.Similar procedure is implemented both for 1D and 2D finite elements and their corresponding nodes.

CreateMatSolid method
In order to define 3D finite elements in SAP2000, solid sections need to be defined first.Each solid section property, containing among others the material characteristics, is assigned to groups of solid finite elements.According to SIMP, the elements' modulus of elasticity are penalized through density values x e using the formula of Eq. ( 2).This is carried out by the CreateMatSolid method, where a number of material classes are generated first (for example 101 different material classes as described in lines #72 to #77 of the following code); initially, having identical properties and they are stored in array MatName.This is described in the following part of CreateMatSolid method: In subsequent part of CreateMatSolid method (i.e.lines #79 to #87), the modulus of elasticity is modified according to the following formula (see Line #84): In order to avoid section properties having modulus of elasticity values equal to zero, that might generate numerical instability issues during the analysis, the first material property corresponding to the lowest modulus of elasticity value becomes equal to E 0 ¼ 1:0E -07 Â E 0 (where E 0 corresponds to the first material class while E 0 is the initial value of the modulus of elasticity).This is described in lines #79 and #80 of the code provided below along with the generation of the rest of the material properties (i.e.lines #81 to #87): In the last part of the code of CreateMatSolid method (i.e.lines #89 to #96) each material property defined according to Eq. ( 34) is assigned to a solid section property.Each solid section is labelled as OptiSolidProp#, where # stands for an id value that denotes the modulus of elasticity value of the specific section (i.e.OptiSolidProp0.01stands for the solid section property whose modulus of elasticity is equal to E i ¼ 0:01 Â E 0 ).This is described in the code below: 89 90 91 92 93 94 95 96 metr = 0.00; for (i = 0; i <= 100; i++) { temp_name = "OptiSolidProp" + string.Format("{0:N2}", metr); ret = Model.PropSolid.SetProp(temp_name,MatName[i],0.0,0.0,0.0,true,-1,"","");metr = Math.Round(metr, 2) + 0.01; metr = Math.Round(metr, 2); } Using an instance of CreateMatSolid method, modulus of elasticity values that will be used in the following steps of the topology optimization procedure are assigned to solid section properties.

CreateOptiSolidName method
The use of an instance of the CreateOptiSolidName method creates the array OptiSolidName containing the ids of the model's solid elements whose density values will be the design variables of the topology optimization problem.The optimizable solid finite elements are generated into a rectangular domain (optimizable domain) that is defined by the user providing the bounds of its three dimensions (i.e.xmin, xmax, ymin, ymax, zmin and zmax, see Fig. 3).In addition, the dimension of the quadratic finite elements needs also to be provided by the user.
CreateOptiSolidName method based on the nodal coordinates of the solid finite elements, in line #119 selects which elements will be optimized using OAPI function Model.SelectObj.CoordinateRange.Next, in line #124 the method obtains the id (i.e.ObjNameSel) of the selected solid element and stores it in the proper location of array OptiSolidName (see line #134).The arrangement of the solid elements' ids in the array OptiSolidName is the same with that used in (Liu and Tovar 2014).An example of the numbering implementation for a mesh 5 9 395 of optimizable solid finite elements (along x, y and z axes, respectively) can be seen in Fig. 4

CreateMapDist method
In order to avoid the checkerboard problem filtering needs to be applied both to new density values and derivatives.In every iteration the filter is applied to each solid finite element based on its neighbouring elements.To facilitate this process an instance of the CreateMapDist method needs to be used first before initiating the iterations.This method is used for creating the array Map that contains the ids of each solid finite element's neighbouring elements.The first three iterative loops (that begin in lines #175, #177 and #179) are implemented in order to identify in line #181 each solid finite element (e_int) while the next three iterative loops (that begin in lines #189, #191 and #193) are employed in order to find in line #195 all its neighbouring located into a sphere having radius equal to rmin.In addition, CreateMapDist method calculates the distances H e;i ¼ max 0; rmin À Dðe; iÞ ð Þ for each solid finite element (see line #196), where D(e,i) is the distance between the solid finite element e and every of its neighbouring solid finite elements i, and stores them in array Dist (see line #198).
The procedure of finding neighbouring elements that is implemented with the formulas of lines #181 and #195 is based on the elements' names that Sigmund (2001) proposed and Liu and Tovar (2014) also applied in their 3D code and is currently adopted in CreateOptiSolidName method.

UpdateModel method
In this part of the VS project, the solid section properties with the corresponding modulus of elasticity values are assigned to the optimizable solid finite elements.This is performed with the use of an instance of UpdateModel method as described in the lines of the code below: The 3rd power penalization of the density value according to SIMP is calculated in line #217, over the total number of solid finite elements, whose density values represent variables of the problem.In line #218, this density value is converted into a string variable and is stored in the temporary variable temp_string composed by the prefix ''OptiSolidProp'' and followed by the value (i.e.OptiSolidProp0.01).
In line #219, a solid section property is assigned to every optimizable solid finite element based on the name that is stored in temporal variable temp_string.

CalcCompliance method
In this part of the code, the compliance of the structure is calculated.In line #341 OAPI function SolidJointForce is used over the solid elements in order to obtain from SAP2000 joint forces F1, F2, F3, M1, M2 and M3 for every node and they are stored in arrays.In the next for loop, function JointDispl is used in order to obtain the displacements U1, U2, U3, R1, R2 and R3 also for every node.The cross product between the arrays of displacements and forces over the nodes of each element results into the total compliance of the structural system (see line #347).Worth mentioning that the calculation of compliance is not a required part for the implementation of the optimization proccess, however, it is provided as a measure for the progress of the search procedure.

CalcDerivatives method
The calculation of the partial derivatives with respect to x e is the most important part of the iterative steps of the optimization procedure.This is performed with an instance of CalcDerivatives method, whose input argument is the array xkfil.According to the theoretical part provided previously the partial derivatives are calculated using Eq. ( 8).However, the finite elements' stiffness matrices k e ðxÞ are not available through the OAPI of SAP2000.In order to overcome this difficulty a modified expression is proposed in this study.According to SIMP changing the modulus of elasticity of the optimizable solid finite elements, results into modified stiffness matrices: and Thus, Eq. ( 8) can be rewritten as follows: oCðxÞ ox e ¼ Àu e ðxÞ T ok e ðx e Þ ox e u e ðxÞ ¼ À u e ðxÞ T 3 k e ðx e Þ x e u e ðxÞ CalcDerivatives method is used to calculate the work produced by the optimizable solid finite elements based on the procedure described in the previous method (CalcCompliance) and then to calculate the partial derivatives according to Eq. ( 37).The method returns the array der containing the derivatives arranged in a similar fashion with the rest ones that have been created already.

FilterDer method
Given that the derivatives have been calculated, their values need to be filtered in order to avoid the checkerboard problem.
Finally, FilterDer method returns the array derfil[i] containing the filtered derivatives.For each value of k defined as the average value of the corresponding interval, a new density value (variable xnew[i]) is calculated according to Eq. ( 19), using the unfiltered density value of the previous iteration (variable xk[i]).
This is described in the code that follows: In the code of lines #596 to #608 the constraint described in Eq. ( 21) is verified.If the sum (variable metrx) of filtered density values (x e ) over all optimizable solid finite elements is greater than the required volume fraction value (V), then density values need to be reduced, and therefore k (variable lamda) that lies in the interval 0:5ða þ bÞ; b ð Þneeds to be increased.Accordingly, if the sum of density values (x e ) is less than the required volume fraction value, density threshold value 0.6 then the procedure described in lines #782 to #784 is performed where its solid section property is updated based on its density value, otherwise the element is deleted (see line #788).When this threshold checked is performed for all optimizable solid finite elements, the optimized layout can be seen through SAP2000 environment.

Test examples
In order to present the capabilities and to provide more information about the use of TOCP source code several test examples are considered in this section.Initially a couple of simple test examples such as the cantilever and simply supported beams are examined.In the second group of tests indicative issues related to conceptual design are integrated into the optimum layout design of moment resisting frames (MRFs) that represent an important representative of the structural members in case of building structures.For all optimization runs that will follow the value of p used in Eq. ( 2) was set equal to 3 and density filtering was applied having filter size r min ¼ ffiffi ffi 3 p .

Simple examples
In the collection of simple examples a cantilever beam and two simply supported ones (beam-A and -B) are studied.For all optimization runs the value of the volume fraction (volfrac) used to obtain the optimized domains was set equal to 30%.
The dimensionless lengths and height, of the cantilever beam of Fig. 5a are 1.2 along the longitudinal axis x, 0.5 along the transverse axis y and 0.8 along the vertical axis z.For the finite element discretization of the cantilever beam shown in Fig. 5a, the following parameters were used: nlx = 12, nly = 5 and nlz = 8 resulting into 480 cubic solid elements where each edge of the basic cubic solid is equal to 0.1 dimensionless length size.All nodes on the right side of the cantilever beam were restrained, while a uniformly distributed load is applied along the upper part of the left side at the vertical direction (i.e.axis z).Both OC and MMA algorithms were used in order to solve the topology optimization problem converging to similar results, both methods converged based on the maximum number of iterations criterion that was set to 150 finite element analyses.The optimized domain obtained by OC is shown in Fig. 5b, a similar one was obtained by MMA as well.
For the simply supported beam-A of Fig. 6a the corresponding dimensionless lengths and height are 3.0 along the longitudinal axis x, 0.8 along the transverse axis y and 1.0 along the vertical axis z.For the finite element discretization of the simply supported beam-A shown in Fig. 6a, the following parameters were used: nlx = 30, nly = 8 and nlz = 10 resulting into 2400 cubic solid elements where each edge of the basic cubic solid is equal to 0.1 dimensionless length size.The nodes of the two lower extreme edges were restrained, while a uniformly distributed load is applied along the centre of the beam at the vertical direction i.e. axis z (as shown in Fig. 6a, corresponding to loading case 1).Both OC and MMA algorithms were used in order to solve the topology optimization problem converging again to similar results; however, OC method converged earlier when the stopping criterion of the maximum difference of the density values between two successive iterations was satisfied at the 90th iterations while MMA converged when the maximum number of iterations criterion was achieved that was set to 150 finite element analyses.The optimized domain obtained by OC is shown in Fig. 6b, a similar one was also obtained by MMA.For the simply supported beam-A a second loading case was also considered, where additional uniformly distributed loads are applied along the upper extreme edges of the beam at the vertical direction (i.e.axis z, see Fig. 7a), the optimized domain obtained by OC for loading case 2 is shown in Fig. 7b, a similar one was obtained also by MMA algorithm.
The last simple example is the second simply supported beam-B, where the dimensionless lengths and height are equal to 1.2 along the longitudinal axis x, 0.5 along the transverse axis y and 0.8 along the vertical axis z.For the finite element discretization of the simply supported beam-B shown in Fig. 8, the following parameters were used: nlx = 24, nly = 10 and nlz = 16 resulting into 3840 cubic solid elements where each edge of the basic cubic solid is equal to 0.05 dimensionless length size.The nodes of the lower extreme edges were restrained, while a uniformly distributed load is applied along the upper left edge at the horizontal direction of axis x.Similar to the previous examples both OC and MMA algorithms were used in order to solve the topology optimization problem converging, however, to rather different results both after satisfying the maximum number of iterations criterion that was set to 150 finite element analyses.In particular, OC resulted into the optimized domain presented in Fig. 9a, while MMA algorithm's result is shown in Fig. 9b.

MRF design test example
Moment-resisting frames are rectilinear assemblages of beams and columns, with the beams rigidly connected to the columns.Resistance to lateral loads is provided primarily by rigid frame action, i.e. by the development of bending moments and shear forces in the frame members and joints.By virtue of the rigid beam-column connections, a moment frame cannot displace laterally without bending the beams or columns depending on the geometry of the connection.The bending rigidity and strength of the frame members is therefore the primary source of lateral stiffness and strength for the entire frame.The 1994 Northridge earthquake revealed a common flaw in the construction, and building design codes were revised to strengthen them.
In the context of OPTARCH research project, among others, topology optimization problems are formulated for deriving multiple alternatives.In this direction, topology optimization is used as a tool to design aesthetically acceptable layouts of MRFs used in the design of high-rise buildings.In order to present the integration of topology optimization formulations in the conceptual  10a is employed; where the domain and finite element mesh discretization is also depicted.The dimensionless lengths and height of the MRF are equal to 2.0 along the longitudinal axis x, 0.3 along the transverse axis y and 9.0 along the vertical axis z.Thus, the first step in the TOCP based design procedure of MRFs is to define the shear wall (optimizable domain), which is the equivalent of MRFs in case of RC structures that through the topology optimization procedure is transformed into an MRF.For the optimization runs of MRF test example various values of the volume fraction were used to obtain variants of optimized domains.For the finite element discretization of the MRF shown in Fig. 10a, the following parameters were used in all cases examined, i.e. nlx = 20, nly = 3 and nlz = 90 resulting into 5400 cubic solid elements.Each edge of the basic cubic solid is equal to 0.1 dimensionless length size and the bottom extreme edges of the domain were fully fixed.For the MRF test example two loading cases have been applied.Among others, these two loading cases included, several issues concerning the conceptual design of According to the first case horizontally loads were applied on the left side of MRF (at the level of the three storeys of the structural system, loading case-1) and in the second one on its both sides (loading case-2).For the first loading case, OC algorithm was adopted and the volume fraction used to obtain the optimized domain was set equal to 40%.The optimized domain obtained for loading case-1 by OC is shown in Fig. 10b.For loading case-2, shown in Fig. 11a, both search algorithms were implemented, and OC resulted in the optimized domain of Fig. 11b where the volume fraction was set equal to 55% whereas MMA algorithm resulted in the optimized domain of Fig. 11c where the volume fraction was set equal to 47%.For all optimization runs performed for this test example the maximum number of iterations criterion was set equal to 400 finite element analyses.In this work we present a code implementing the integration of a minimum compliance shape and topology optimization framework into the SAP2000 structural analysis and design software.The code is written in C# programming language and exploiting the open application programming interface provided by SAP2000, it provides the possibility to take advantage of all its modelling, analysis and design capabilities.In addition, for reasons related to the integration of the Topology Optimization problem into SAP2000, further capabilities and implementation features were developed.The paper covers all theoretical aspects of topology optimization incorporated in the code and provides detailed description of their numerical implementation.For purposes of understanding, several simple test examples are also provided in the current study.TOCP has no restrictions in its application with respect to the problem size; it can be used to solve more realistic engineering problems having hundreds of thousands of elements or even more according to the hardware capabilities.
In particular, four visual studio projects implementing variants of a simple minimum compliance based topology optimization frameworks, are available and can be used for educational purposes.The codes can easily be extended to include additional optimization criteria, like specific displacement or frequency minimization etc., also to handle problems requiring nonlinear and/or dynamic analyses, and to implement continuation on the penalization parameter.Using the four projects in C# are not independent, since the finite element analysis part required by topology optimization is performed by SAP2000; this is considered as the basic advantages of the TOCP codes provided herein.This is because, the analysis and design features that the commercial software provides will enhance significantly the capabilities of future researchers that will extend the TOCP codes' capacities.Furthermore, the user needs not to care for issues like sparsity in the assembly of the global stiffness matrix, acceleration of the finite element equations solution, modelling and analysis issues.Since they will benefit from the acceleration capability and memory management that advanced commercial software have inherently implemented in its computational core.Worth mentioning that in case of academic purposes usage the specific software provides multiple licenses almost for free.
The authors would be happy to receive suggested improvements that can be implemented in the public domain TOCP codes (please address to the first author by his e-mail address nlagaros@central.ntua.gr).

Fig. 5
Fig. 5 Simple test examples: cantilever beam: a domain-mesh discretization and b optimized layout obtained by OC (volume fraction equal to 30%)

Fig. 7
Fig. 7 Simple test examples: simply supported beam-A (loading case 2): a domain-mesh discretization and b optimized layout obtained by OC (volume fraction equal to 30%)

Fig. 9
Fig. 9 Simple test examples: simply supported beam-B, 3D and side views of the optimized domain obtained by (volume fraction equal to 30%): a OC and b MMA

Fig. 10
Fig. 10 MRF test example: a domain-mesh discretization and b optimized layout obtained by OC (volume fraction equal to 40%)

Fig. 11
Fig. 11 MRF test example: a domain-mesh discretization and b optimized layout obtained by OC (volume fraction equal to 55%) and c optimized layout obtained by MMA (volume fraction equal to 47%) .
The input arguments of FilterDer method are the arrays der, xkfil, Map and Dist defined previously by an instance of CreateMapDist method.For each optimizable solid finite element, the ids of its neighbouring elements are retrieved through array Map, the corresponding distances Hei are obtained through array Dist and then the summations of the following expression are calculated: X In order to avoid the checkerboard problem filtering is performed over the optimizable solid elements' density values x e according to the following expression: