SIMP element for topology optimization

[PhiHo-eng]

Hi everyone,

I am busy with the reactivation of the Topology Optimization Application, but I am having trouble implementing the SIMP element based on the small-dispalcement solid element (small_displacement.cpp) of the Structural Mechanics Application. In the legacy TopOpt code, the SIMP element was based on the Solid Mechanics Application and it seems that the implementation "philosophy" of these elements is different. I am specifically having a problem with the Calculate function.

I am developing this in the following repository: MyRepository

The test example that I am using: exmaples/01_Small_Cantilever_Hexahedra/ (01_Small_Cantilever_Hexahedra/)

I am not able to get to the root of this problem, does anyone have en idea and can help me with this?

Thanks in advance!

[loumalouomega]

Taking a fast at your code I have many comments, but the important thing, what is the problem?

(and the approaches are no so different)

[PhiHo-eng]

Thanks for that super fast answer!

This is the error that I get:

Linear-Solver-Factory: Constructing a regular (non-complex) linear-solver 
terminate called after throwing an instance of 'Kratos::Exception'
  what():  Error: attempting to do a GetValue for: CONSTITUTIVE_LAW variable #503075072CONSTITUTIVE_LAW variable #503075072 #503075072 unfortunately the variable is not in the database and the operations is not threadsafe (this function is being called from within a parallel region)

in kratos/containers/data_value_container.h:199:TDataType& DataValueContainer::GetValue(const Variable<TDataType>&) [with TDataType = shared_ptr<ConstitutiveLaw>]
   /home/philipp/opt/kratosDev/applications/StructuralMechanicsApplication/custom_elements/base_solid_element.cpp:236:virtual void BaseSolidElement::InitializeMaterial()
   /home/philipp/opt/kratosDev/applications/StructuralMechanicsApplication/custom_elements/base_solid_element.cpp:67:virtual void BaseSolidElement::Initialize(const ProcessInfo&)

terminate called recursively
Abgebrochen (Speicherabzug geschrieben)

My assumtion is that this error, due to the slightly different functions, is caused in the calclaute() funtion. Would it be possibel for you to make a quick video chat via MSTeams or Zoom in order to give me all the comments you have for me?!

Thanks in advance!

[philbucher]

considering where the error comes from I would say you didn't specify a constitutive law

[PhiHo-eng]

Thank you for checking out the error.

I am not quite sure, what you mean. Do I have to make this specification in the SIMP_element or do I have to specify this in the runTopOpt file?

[philbucher]

usually the constitutive law is set in the solver If I understand correctly you are using the basic structuralmechanics workflow and slightly modify it to make the topopt work no? Does your example run in the regular workflow, aka without your modifications?

[PhiHo-eng]

Thank you again for helping me! I was able to modify the code, so that the example launched.

But I still have a problem, that the objective function is always zero. So it does not change at all. It seems that Calculate() in structure_response_function_utilities.h is either called in the wrong way or I made a mistake by defining it wrong in the small_dispalcement_eimp_element.cpp?

Do you have any idea on this or can give me some tips how I might be able to solve this?

Thanks in advance!

[loumalouomega]

Can you post the code, or link it, please?

[PhiHo-eng]

Thanks for looking into it!

This is the link to the code in my repository: structure_response_function_utilities.h

[loumalouomega]

Comment, I would clean that includes, but nothing related with the topic of this issue

[loumalouomega]

Two notes, use CalculateOnIntegrationPoints, Calculate methods are not implemented (@philbucher we should clean that, and deprecate them)

The method must be implemented in the solid element, because otherwise will return an empty vector

[PhiHo-eng]

I still get a value of zero for the objective function, with the CalculateOnIntergrationPoints. I think this function is not even called.

Do you mean to ues BaseSolidElement as BaseType for the element rather than SmallDisplacement?

[loumalouomega]

You should implement CalculateOnIntergrationPoints for LOCAL_STRAIN_ENERGY. If the implementation is general you can do it in BaseSolidElement, otherwise SmallDisplacement will be OK

[PhiHo-eng]

Still learning Kratos here and sorry that I don't understand. Where should implement CalculateOnIntegrationPoints?

I already tried it, but I am not sure if it is right?!

structure_respone_function_utilities.h

[loumalouomega]

In the element you are using (or in upper classes if it may be used there too). For example:

Take the old implementation: https://github.com/KratosMultiphysics/KratosLegacyApplications/blob/05ea65215990ba8e3999627341bd20c019b84775/TopologyOptimizationApplication/custom_elements/small_displacement_simp_element.cpp#L140

And put in that element in the CalculateOnIntegrationPoints

[PhiHo-eng]

I implemented CalcualteOnIntegrationPoints as you suggestet:

small_displacement_simp_element.cpp

But that gives me an error while compiling, because CaluclateOnIntegrationPoints is already used in the elment:

[ 90%] Building CXX object applications/TopologyOptimizationApplication/CMakeFiles/KratosTopologyOptimizationApplication.dir/custom_elements/small_displacement_simp_element.cpp.o
In file included from /home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/custom_elements/small_displacement_simp_element.cpp:15:
/home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/custom_elements/small_displacement_simp_element.hpp:274:18: warning: ‘virtual void Kratos::SmallDisplacementSIMPElement::save(Kratos::Serializer&) const’ can be marked override [-Wsuggest-override]
  274 |     virtual void save(Serializer& rSerializer) const;
      |                  ^~~~
/home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/custom_elements/small_displacement_simp_element.hpp:275:18: warning: ‘virtual void Kratos::SmallDisplacementSIMPElement::load(Kratos::Serializer&)’ can be marked override [-Wsuggest-override]
  275 |     virtual void load(Serializer& rSerializer);
      |                  ^~~~
In file included from /home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/custom_elements/small_displacement_simp_element.cpp:17:
/home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/topology_optimization_application.h:125:23: warning: ‘virtual std::string Kratos::KratosTopologyOptimizationApplication::Info() const’ can be marked override [-Wsuggest-override]
  125 |   virtual std::string Info() const
      |                       ^~~~
/home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/topology_optimization_application.h:131:16: warning: ‘virtual void Kratos::KratosTopologyOptimizationApplication::PrintInfo(std::ostream&) const’ can be marked override [-Wsuggest-override]
  131 |   virtual void PrintInfo(std::ostream& rOStream) const
      |                ^~~~~~~~~
/home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/topology_optimization_application.h:138:20: warning: ‘virtual void Kratos::KratosTopologyOptimizationApplication::PrintData(std::ostream&) const’ can be marked override [-Wsuggest-override]
  138 |       virtual void PrintData(std::ostream& rOStream) const
      |                    ^~~~~~~~~
/home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/custom_elements/small_displacement_simp_element.cpp:205:6: error: redefinition of ‘void Kratos::SmallDisplacementSIMPElement::CalculateOnIntegrationPoints(const Kratos::Variable<double>&, std::vector<double>&, const Kratos::ProcessInfo&)’
  205 | void SmallDisplacementSIMPElement::CalculateOnIntegrationPoints(const Variable<double>& rVariable, std::vector<double>& rOutput,
      |      ^~~~~~~~~~~~~~~~~~~~~~~~~~~~
/home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/custom_elements/small_displacement_simp_element.cpp:132:6: note: ‘virtual void Kratos::SmallDisplacementSIMPElement::CalculateOnIntegrationPoints(const Kratos::Variable<double>&, std::vector<double>&, const Kratos::ProcessInfo&)’ previously defined here
  132 | void SmallDisplacementSIMPElement::CalculateOnIntegrationPoints(const Variable<double> &rVariable, std::vector<double>& rOutput, const ProcessInfo &rCurrentProcessInfo)
      |      ^~~~~~~~~~~~~~~~~~~~~~~~~~~~
make[2]: *** [applications/TopologyOptimizationApplication/CMakeFiles/KratosTopologyOptimizationApplication.dir/build.make:102: applications/TopologyOptimizationApplication/CMakeFiles/KratosTopologyOptimizationApplication.dir/custom_elements/small_displacement_simp_element.cpp.o] Fehler 1
make[1]: *** [CMakeFiles/Makefile2:583: applications/TopologyOptimizationApplication/CMakeFiles/KratosTopologyOptimizationApplication.dir/all] Fehler 2
make: *** [Makefile:130: all] Fehler 2

But it got me thinking, that I am not sure if I should calculate on the Integrationpoints, which would give me the strain of the given point (please correct me if I am wrong). But I need the displacement of the element in order to calculate the strain energy (compliance), which for the topology optimization is defined by u^T K u, or at element level ue^T Ke ue.

[loumalouomega]

You should add it as a case in the existing implementation then

[PhiHo-eng]

I do not know what I am doing wrong here and have lost the overview. I need to get the element displacements to calculate the "compliance" for the topopt. This should be all in CalculateOnIntegrationPoints in my SIMPelement?

I tried to do so, but it is still not calculating this (or even going into the function):

small_displacement_simp_element.cpp

[loumalouomega]

GetValueOnIntegrationPoints was removed recently you should move that code to CalculateOnIntegrationPoints

[PhiHo-eng]

Thank you for your help here and your fast responses, it is much appreciated! Unfortunately it is still not working for me and it appears not even to enter the function CalculateOnIntegrationPoints. Can that be? I really do not understand what I am doing wrong here.

[philbucher]

How do you call CalculateOnIntegrationPoints? This is not called in the "standard" workflow => should be called by you when you want to calculate sth

[loumalouomega]

How do you call CalculateOnIntegrationPoints? This is not called in the "standard" workflow => should be called by you when you want to calculate sth

It is inside some utilities

[loumalouomega]

Thank you for your help here and your fast responses, it is much appreciated! Unfortunately it is still not working for me and it appears not even to enter the function CalculateOnIntegrationPoints. Can that be? I really do not understand what I am doing wrong here.

Can you link the current code, so we can take a look

[PhiHo-eng]

CalculationOnIntergrationPoint is called here:

structural_response_function_utilities.h

and this is the element:

small_displacement_simp_element.cpp

And this is the Repository: Repository

[loumalouomega]

Does it enter here? https://github.com/PhiHo-eng/Kratos/blob/ebcb50568eab7cf1a899644e924415d913a74e5e/applications/TopologyOptimizationApplication/custom_elements/small_displacement_simp_element.cpp#L205

[PhiHo-eng]

Right, I think that is my problem. I do not think it is calling this function. Can this be?

[loumalouomega]

The code:

std::vector<double> Out;
element_i->CalculateOnIntegrationPoints(LOCAL_STRAIN_ENERGY, Out, ConstProcessInfo);
Global_Strain_Energy += element_i->GetValue(LOCAL_STRAIN_ENERGY);

Is wrong, the element does not have LOCAL_STRAIN_ENERGY, it is in the Out vector (one value per GP), if this is a GP value.

Maybe I am lost. @philbucher Calculate should set the value in the element, or it is a duplication of the CalculateOnIntegrationPoints?

Sorry, I think i guided you wrongly

[PhiHo-eng]

I am definitely lost. Thanks for the answers nonetheless. It is my understanding that LOCAL_STRAIN_ENERGY is not a GP value. The differences in the old Solid App and the Structural App are just enough to throw me off. Is it correct that the calculate command previously existed? Or is this what I need to implement in my element? I would appreciate any further tips!

[philbucher]

It is my understanding that LOCAL_STRAIN_ENERGY is not a GP value

If it is a GP value I would use aka implement CalculateOnIntegrationPoints. If not use Calculate

Once you implemented it you call it. What I still don't get is if it enters the corresponding function or not. Can you check?

[PhiHo-eng]

I implemented Calculate now:

small_displacement_simp_element.cpp

structure_response_function_utilities.h

But I think it is still not calling the function. The value LOCAL_STRAIN_ENERGY remains zero.

[philbucher]

are you sure this ModelPart has those elements? How do you create them?

[PhiHo-eng]

I think so, shouldn't it take the element-type from the description in the .mdpa file?

The element is created in the simp_element.cpp right (please correct me if I am wrong). I tried to keep the element as close as possible to the form of the old TopOpt Version.

[philbucher]

I think so, shouldn't it take the element-type from the description in the .mdpa file?

The element is created in the simp_element.cpp right (please correct me if I am wrong). I tried to keep the element as close as possible to the form of the old TopOpt Version.

So you have this element in the mdpa file?

I would make sure they are in this ModelPart, seems like sth weird is happening

[PhiHo-eng]

In Python I used the command "model_part.NumberOfElements()" to find that there were 8640 elements, the same number in my model. So it is reading these in from the MDPA file. Is there a way to check which element type they were read in as? Or can I assume that we are good?

[philbucher]

I would assume that it works but apparently it doesn't :)

More serious, can you debug a but your StructureResponseFunctionUtilities? I meant to check the ModelPart that you use in there. E.g. check the number of elements and print the elements (e.g. with cout)

Regarding the mdpa, what name do you use for the elements block?

[PhiHo-eng]

I just went through and it goes into structure_response_function_utilities.h and into the the loop where local strain energy is to be caclualted and summed. It runs trough the loop as many times as there are elements (which is a good sign), althought when when we are in the loop, it does not go into my new element. This is strange because the model is analyzed, i.e. I can open the postprocessing file in GiD and the structural response makes sense. The name for the elements block is SmallDisplacementSIMPElement3D8N.

[loumalouomega]

But after puting std:::cout does it enter in the function?

[PhiHo-eng]

No it seems, that it never even enters the function.

[loumalouomega]

Can you check that you enter the constructor of the element? (to ensure that is even created)

[PhiHo-eng]

Do you mean here?

[loumalouomega]

In that one and the previous one too (all the creates)

[PhiHo-eng]

I made a std::cout inside the creation of the elements, and it never gets printed. So this means the elements don't even get created?

[loumalouomega]

Yep, can you share the input files?

[PhiHo-eng]

The input files are in here.

[loumalouomega]

Strange, should work.

Can you try the following python code?: (to check it enters in the constructor)

NOTE: Applications imports are missing

current_model = KratosMultiphysics.Model()
mp = current_model.CreateModelPart("solid_part")

#create nodes
mp.CreateNewNode(1, 0.00000,  1.00000,  1.00000)
mp.CreateNewNode(2, 0.16500,  0.74500,  0.70200)
mp.CreateNewNode(3, 0.27300,  0.75000,  0.23000)
mp.CreateNewNode(4, 0.78800,  0.69300,  0.64400)
mp.CreateNewNode(5, 0.32000,  0.18600,  0.64300)
mp.CreateNewNode(6, 0.00000,  1.00000,  0.00000)
mp.CreateNewNode(7, 0.00000,  0.00000,  1.00000)
mp.CreateNewNode(8, 1.00000,  1.00000,  1.00000)
mp.CreateNewNode(9, 0.67700,  0.30500,  0.68300)
mp.CreateNewNode(10, 0.24900,  0.34200,  0.19200)
mp.CreateNewNode(11, 0.85000,  0.64900,  0.26300)
mp.CreateNewNode(12, 0.82600,  0.28800,  0.28800)
mp.CreateNewNode(13, 0.00000,  0.00000,  0.00000)
mp.CreateNewNode(14, 1.00000,  1.00000,  0.00000)
mp.CreateNewNode(15, 1.00000,  0.00000,  1.00000)
mp.CreateNewNode(16, 1.00000,  0.00000,  0.00000)

#create Element
mp.CreateNewElement("SmallDisplacementSIMPElement3D8N", 1,[10,5,2,3,13,7,1,6], mp.GetProperties()[1])
mp.CreateNewElement("SmallDisplacementSIMPElement3D8N", 2,[12,9,5,10,16,15,7,13], mp.GetProperties()[1])
mp.CreateNewElement("SmallDisplacementSIMPElement3D8N", 3,[12,11,3,10,9,4,2,5], mp.GetProperties()[1])
mp.CreateNewElement("SmallDisplacementSIMPElement3D8N", 4,[9,4,2,5,15,8,1,7], mp.GetProperties()[1])
mp.CreateNewElement("SmallDisplacementSIMPElement3D8N", 5,[4,11,3,2,8,14,6,1], mp.GetProperties()[1])
mp.CreateNewElement("SmallDisplacementSIMPElement3D8N", 6,[11,4,9,12,14,8,15,16], mp.GetProperties()[1])
mp.CreateNewElement("SmallDisplacementSIMPElement3D8N", 7,[11,12,10,3,14,16,13,6], mp.GetProperties()[1])

[PhiHo-eng]

I tried the code and it seems to work, I don't get an error.

This is the code I used Code.

[philbucher]

can you then also try to call Calculate from python in this small example to see if it enters there too?

[PhiHo-eng]

Also with the "self"-made elements, the function Calculate() is not called! It loops over the elements in the structure_response_function_utlities.h but i wont enter into the Calculate() -function.

[philbucher]

one thing I noticed is that your Calculate method is missing the override hence you might not actually override the fct! => this means that you end up calling a different function

[philbucher]

you should add the override, if it no longer compiles you found the error => then fix the fct signature until it compiles

[PhiHo-eng]

Thank you for the input. I was able to compile it, but it still doesn't work. But I think I might have found the error. It could be in here: topology_optimization_application.h and topology_optimization_application.cpp. I based the elements on small_displacement.h and .cpp, so in the beginning I uesd those elements to get the Application running. Now these elements are still being created in the topology_optimization_application.h and topology_optimization_application.cpp. So (and please correct me if I am wrong) the element will be read as SmallDispacementSimpElement3D8N but it will be saved as a Element of SmallDispacement, and that would be why the function is never being called. I now tried to implement the SmallDispalcementSIMPElements (as it was done in the original Application) as new Elements but I get the following error:

[ 89%] Linking CXX shared library libKratosTopologyOptimizationCore.so
/usr/bin/ld: CMakeFiles/KratosTopologyOptimizationCore.dir/topology_optimization_application.cpp.o: in function `Kratos::KratosTopologyOptimizationApplication::KratosTopologyOptimizationApplication()':
/home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/topology_optimization_application.cpp:68: undefined reference to `Kratos::SmallDisplacementSIMPElement::SmallDisplacementSIMPElement(unsigned long, std::shared_ptr<Kratos::Geometry<Kratos::Node<3ul, Kratos::Dof<double> > > >)'
/usr/bin/ld: /home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/topology_optimization_application.cpp:68: undefined reference to `Kratos::SmallDisplacementSIMPElement::SmallDisplacementSIMPElement(unsigned long, std::shared_ptr<Kratos::Geometry<Kratos::Node<3ul, Kratos::Dof<double> > > >)'
/usr/bin/ld: /home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/topology_optimization_application.cpp:68: undefined reference to `Kratos::SmallDisplacementSIMPElement::SmallDisplacementSIMPElement(unsigned long, std::shared_ptr<Kratos::Geometry<Kratos::Node<3ul, Kratos::Dof<double> > > >)'
/usr/bin/ld: /home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/topology_optimization_application.cpp:68: undefined reference to `Kratos::SmallDisplacementSIMPElement::~SmallDisplacementSIMPElement()'
/usr/bin/ld: /home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/topology_optimization_application.cpp:68: undefined reference to `Kratos::SmallDisplacementSIMPElement::~SmallDisplacementSIMPElement()'
/usr/bin/ld: CMakeFiles/KratosTopologyOptimizationCore.dir/topology_optimization_application.cpp.o: in function `Kratos::SmallDisplacementSIMPElement::SmallDisplacementSIMPElement()':
/home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/custom_elements/small_displacement_simp_element.hpp:224: undefined reference to `vtable for Kratos::SmallDisplacementSIMPElement'
/usr/bin/ld: /home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/custom_elements/small_displacement_simp_element.hpp:224: undefined reference to `vtable for Kratos::SmallDisplacementSIMPElement'
/usr/bin/ld: CMakeFiles/KratosTopologyOptimizationCore.dir/topology_optimization_application.cpp.o: in function `Kratos::KratosTopologyOptimizationApplication::~KratosTopologyOptimizationApplication()':
/home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/topology_optimization_application.h:95: undefined reference to `Kratos::SmallDisplacementSIMPElement::~SmallDisplacementSIMPElement()'
/usr/bin/ld: /home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/topology_optimization_application.h:95: undefined reference to `Kratos::SmallDisplacementSIMPElement::~SmallDisplacementSIMPElement()'
/usr/bin/ld: /home/philipp/opt/kratosDev/applications/TopologyOptimizationApplication/topology_optimization_application.h:95: undefined reference to `Kratos::SmallDisplacementSIMPElement::~SmallDisplacementSIMPElement()'
/usr/bin/ld: CMakeFiles/KratosTopologyOptimizationCore.dir/topology_optimization_application.cpp.o: in function `void Kratos::Serializer::Register<Kratos::SmallDisplacementSIMPElement>(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Kratos::SmallDisplacementSIMPElement const&)':
/home/philipp/opt/kratosDev/kratos/includes/serializer.h:205: undefined reference to `typeinfo for Kratos::SmallDisplacementSIMPElement'
collect2: error: ld returned 1 exit status
make[2]: *** [applications/TopologyOptimizationApplication/CMakeFiles/KratosTopologyOptimizationCore.dir/build.make:89: applications/TopologyOptimizationApplication/libKratosTopologyOptimizationCore.so] Fehler 1
make[1]: *** [CMakeFiles/Makefile2:613: applications/TopologyOptimizationApplication/CMakeFiles/KratosTopologyOptimizationCore.dir/all] Fehler 2
make: *** [Makefile:130: all] Fehler 2
philipp@philipp-IdeaPad-5-15IIL05:~/opt/kratosDev/scripts$ 

Any idea how I might solve this and also if this could even be the source of my problem?