mimesis-inria / caribou

Multi-physics computation library
GNU Lesser General Public License v3.0
29 stars 17 forks source link

Did not succeed in templating with Elements a material. #100

Closed Ziemnono closed 2 years ago

Ziemnono commented 2 years ago

Hi @Sidaty1 and @jnbrunet,

For this project, I created a FEniCS-features branch where you should be able to checkout and compile (hopefully). All the concerning files are in the directory /src/SofaCaribou/FEniCS. As you see, in the file /src/SofaCaribou/FEniCS/Forcefield/HyperelasticForcefield_FEniCS.inl, we use the function

integral->tabulate_tensor_float64(nodal_forces.data(), coefficients.data(), constants, current_nodes_position.data(), nullptr, nullptr);

in addForce and assemble_stiffness to generate the residual vector and the stiffness matrix.

The function comes from const ufcx_integral *integral = form_SaintVenantKirchhoff_Hexa_Order2_J->integrals(ufcx_integral_type::cell)[0]; that should be changed manually depending on the chosen element (Hexa of order 2 in this case).

This is quite uncomfortable for the user (myself included). To solve this problem, I thought about doing something similar to what you were doing by creating a material called /src/SofaCaribou/FEniCS/Material/SaintVenantKirchhoffMaterial_FEniCS and calling a function called calculate_nodal_forces instead of the PK2_stress for example.

This calculate_nodal_forces function needs to know the size of the R vector (NumberOfNodePerElement, Dimension) that depends on the element type and the interpolation degree.

I tried to create a HyperelasticMaterial_FEniCS class in /src/SofaCaribou/FEniCS/Material/ to be able to template the Element type and be able to use the NumberOfNodePerElement variable for example. Once done, I am adding it to the factory in /src/SofaCaribou/Material/HyperelasticMaterial.cpp

Once I am running the example ./Validation/fenics/cantilever_beam_SOFA_FENiCS_Hexa_quadratic.py. I obtain the following message: [ERROR] [ObjectFactory] Component already registered: SaintVenantKirchhoffMaterial_FEniCS<Vec3d> This means that the SaintVenantKirchhoffMaterial_FEniCS is only templated using DataTypes while I wanted to template it with Element, DataTypes.

Thank you and have a nice day.

jnbrunet commented 2 years ago

Hey @Ziemnono !

You even created a github project, this is great ! Let's hope it takes more momentum than the last time I tried it :-)

I am super close to have a Windows workflow which will allow us to produce nightly binaries for Windows users. I just want to finish this once and for all, and then I'm 100% with this new project.

For the templated-based material we will have to manage the instantiation of the good template from the scene definition. For example, the following:

meca.addObject('CaribouTopology', template='Hexahedron')
meca.addObject('CaribouMass', density=mass_density, template='Hexahedron')
meca.addObject('SaintVenantKirchhoffMaterial', young_modulus=youngModulus, poisson_ratio=poissonRatio)
meca.addObject('HyperelasticForcefield', template="Hexahedron")

would become

meca.addObject('CaribouTopology', template='Hexahedron')
meca.addObject('CaribouMass', density=mass_density, template='Hexahedron')
meca.addObject('SaintVenantKirchhoffMaterial', template="Hexahedron", young_modulus=youngModulus, poisson_ratio=poissonRatio)
meca.addObject('HyperelasticForcefield', template="Hexahedron")

I will try to put some time into it next week. Should I push directly in your branch? I can also use pull requests which can be done on any target branch.

/cc @hugtalbot and @fredroy

Ziemnono commented 2 years ago

Hi @jnbrunet , Thank you for responding so fast. Yes that was the idea that I have in mind, we could even change to:

meca.addObject('CaribouTopology', template='Hexahedron')
meca.addObject('CaribouMass', density=mass_density, template='Hexahedron')
meca.addObject('SaintVenantKirchhoffMaterial', template="Hexahedron", young_modulus=youngModulus, poisson_ratio=poissonRatio)
meca.addObject('HyperelasticForcefield')

Where the HyperelasticForcefield could deduce its template from the material. If you can make a PR it would be perfect, I will be working in the company for 2 weeks so I will not have too much time to dig into those problems ! Thank you!