Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
Ziemnono opened this issue Jan 28, 2022 · 2 comments
Closed

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

Ziemnono opened this issue Jan 28, 2022 · 2 comments
Assignees
Labels
help wanted Extra attention is needed question Further information is requested

Comments

@Ziemnono
Copy link
Collaborator

Ziemnono commented Jan 28, 2022

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.

  • What did I do wrong?
  • If you can think about a better solution, do not hesitate !
  • If the issue is unclear tell me where or we can also have a meeting to discuss it.

Thank you and have a nice day.

@Ziemnono Ziemnono added help wanted Extra attention is needed question Further information is requested labels Jan 28, 2022
@jnbrunet
Copy link
Member

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
Copy link
Collaborator Author

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!

@Sidaty1 Sidaty1 closed this as completed Mar 19, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
help wanted Extra attention is needed question Further information is requested
Projects
None yet
Development

No branches or pull requests

3 participants