Fiber direction import from VTU files for anisotropic strain energy models in GetFEM

General
1

Dear all,

I am currently transitioning from FEniCS to GetFEM and planning to implement an anisotropic strain energy model inspired by the works of Holzapfel and Ogden. These types of models typically involve preferred directions associated with fiber and sheet structures in biological tissues.

A generic example of the strain energy function I aim to use can be written as:
\Psi = \Psi_{\text{iso}}(C) + \sum_{i=1}^N \Psi_{\text{aniso}}(C, \mathbf{a}_i)

where \Psi_{\text{iso}}(C) is the isotropic part depending on the right Cauchy-Green tensor C = F^T F, and \Psi_{\text{aniso}} includes the fiber contributions along preferred directions \mathbf{a}_i. These contributions often take the form:
\Psi_{\text{aniso}}(C, \mathbf{a}_i) = k_i \left( \mathbf{a}_i \cdot C \mathbf{a}_i - 1 \right)^2,

and involve terms like \mathbf{a}_i that need to be precomputed and/or interpolated on the mesh.

In my previous FEniCS workflow, I used the Laplace-Dirichlet Rule-Based (LDRB) method to compute fiber and sheet directions (see images). These directions were then saved as vector fields in VTU format (one vector per node or cell).

To avoid rewriting all the code I already developed in FEniCS, I would like to know:
Is there a way in GetFEM to read solutions stored in a VTU file and load them as vector variables or simply variables, so that they can be used later in the computation of quantities such as
\mathbf{a}_i \otimes \mathbf{a}_i or for futher calulus in GWFL?

This would greatly simplify the transition and allow me to reuse the preprocessing steps already implemented in my current pipeline.

Any guidance or example on how to read VTU files (or any alternative approach to inject external directional data into a GetFEM simulation) would be highly appreciated.

Best regards,

radial
radial517Ă—399 173 KB

circ
circ513Ă—375 134 KB

2

Normally you need to have the unit direction vectors at all integration points that you are going to use.

In a 2D model where we import experimentally measured fiber orientations we use something like this:

mimd1 = gf.MeshImData(mim, -1)

# Model
md = gf.Model("real")

...
md.add_im_data("s0", mimd1)         # sin(phi0): sine of the fiber angle
md.add_im_data("c0", mimd1)         # cos(phi0): cosine of the fiber angle

x = np.arange(-L/2, L/2+1e-6, voxel_size)
y = np.arange(-H/2, H/2+1e-6, voxel_size)
fiber_angle_distribution_interp = scipy.interpolate.RegularGridInterpolator((y, x), fiber_angle_distribution)

X = md.interpolation("X(1)", mimd1, -1) # x-coordinate of all integration points in mim
Y = md.interpolation("X(2)", mimd1, -1) # y-coordinate of all integration points in mim
PHI0 = fiber_angle_distribution_interp(np.stack((Y,X), axis=1)
md.set_variable("s0", np.sin(PHI0))
md.set_variable("c0", np.cos(PHI0))

Alternatively, you can also export the vectors X and Y (and Z in 3D), and ask FEniCS to interpolate the quantity at these points.

But in general if you have some code that calculates the fiber directions in FEniCS it should be rather trivial to port to GetFEM.

3

First of all, thank you for your previous explanation, I would try to implement it. Indeed, these are unit vectors at each Gauss point. you can see the getfem version results (not getting the same result like in Fenics still have to fixe it)

I would like to go a bit deeper with my current question. In some physical problems, especially in time-dependent simulations, the initial time steps are not necessarily significant for the final results. These early steps are mainly used to reach an intrinsic steady or equilibrium state before coupling the system with a more complex physical problem.

In this context, I would like to know whether it is possible, within GetFEM, to stop the simulation, save the current solution state (check-pointing), and then use this state as an initial or boundary condition for another problem so reused vtu data.

Additionally, I would highly appreciate if you could share any documentation, notes, or resources related to anisotropic hyperelastic material models in GetFEM.

Thank you very much for your support and time.

Best regards,

radial
radial499Ă—437 192 KB

cir
cir450Ă—396 150 KB

4

this is an excerpt from one of my codes that uses check-pointing:

    if rank == 0:
      output = (mfu, md.variable("u"), "Displacements",
                mfx, md.variable("x"), "x",
                mfx, md.variable("x")-md.variable("x0"), "dx",
                ...)
      mfout.export_to_vtu(f"{resultspath}/optimization_step_{i}.vtu", *output)
      np.save(f"{resultspath}/x_{i}.npy", md.variable("x"))
      np.save(f"{resultspath}/state_{i}.npy", md.from_variables())

and this is the corresponding snippet for resuming a run:

if initstate > 0:
  md.to_variables(np.load(resultspath+f"/state_{initstate}.npy"))
  ...
else:
  ... do some pre-simulations

... actual simulation

The vtu format is a format that is meant for post-processing, it is not meant for importing data to a model, and such a use case is not supported in GetFEM.

5

Regarding anisotropic hyperelasticity, there is nothing specific to GetFEM that you really need to learn. You can either type in your formula for the elastic strain energy density, as a function of \nabla u (0-th order expression), using GWFL syntax, or the expression for the virtual work (1st order expression). Normally I choose the second, and I type in the virtual work expression with GWFL syntax in the format P(u):Grad(Test_u).

Just learn how to write the P(u) closed form expression in a form as simple as possible, to avoid unnecessary computations. In the literature, you will see many intermediate quantities computed to go from \nabla u to P (1st PK), but they are not really necessary for the implementation, you only need the deformation gradient F and the direction unit vector T.

In general, I would recommend you to try to learn the formalism I use in my hyperelasticity notes. They include also a case for orthotropic hyperelasticity.

Otherwise, we have a series of publications like this one, which involve anisotropy, plasticity, and higher order strains, and they are all implemented in GetFEM based on GWFL expressions only.

6

Thank you very much for the checkpoint control code and the valuable advice you provided. I’ve implemented everything using the GWFL formulation, and I personally find it more straightforward to work with the first Piola-Kirchhoff tensor through the principle of virtual work.

I’ve written a code for the anisotropic case and I’m currently in the process of testing it. If it performs well, I’ll be happy to share it.

Best regards.

7

Dear Konstantinos,

First, thank you for the suggestions you provided earlier they worked great!
I was able to implement a staggered scheme coupling three different models.

Now I have a new question regarding visualization in a deformed configuration.

I have a scalar field ul defined on a scalar MeshFem:

mf1 = gf.MeshFem(m, 1)
mf1.set_fem(gf.Fem("FEM_PK(3,2)"))

and I also have a displacement field u defined on a vectorial MeshFem:

mf2 = gf.MeshFem(m, m.dim())
mf2.set_fem(gf.Fem("FEM_PK(3,2)"))

I would like to visualize the scalar field ul in the deformed configuration that is, as if it were “attached” to the mesh deformed by u.

Is there a way to:

  • either project ul onto mf2, or
  • associate the scalar field ul with the displacement u during export,
    so that both the displacement and the scalar field appear together in the same VTU file (or any other visualization format), with the mesh deformed by u?

Thanks in advance for your support!

8

just export them in the same file and use the WarpByVector filter in Paraview.

That’s exactly what is done in he code snippet I gave you before:

      output = (mfu, md.variable("u"), "Displacements",
                mfx, md.variable("x"), "x",
                mfx, md.variable("x")-md.variable("x0"), "dx",
                ...)
      mfout.export_to_vtu(f"{resultspath}/optimization_step_{i}.vtu", *output)

this is equivalent to

      mfout.export_to_vtu(f"{resultspath}/optimization_step_{i}.vtu",
                          mfu, md.variable("u"), "Displacements",
                          mfx, md.variable("x"), "x",
                          mfx, md.variable("x")-md.variable("x0"), "dx",
                          ...)

mfx is a scalar mesh_fem and mfu is a vector mesh_fem, both defined on the same mesh. Scalar variable “x”, displacement variable “u”.

1 Like
9

Dear @Konstantinos.Poulios
Thanks very much you helps, I was able to do the check point controls see my strategy and it works.
check-pointing

 state = {
    "variables": md.from_variables(),
    "Previous_ui": md.variable("Previous_ui").copy(),
    "Previous_vi": md.variable("Previous_vi").copy(),
    "Previous_ul": md.variable("Previous_ul").copy(),
    "Previous_vl": md.variable("Previous_vl").copy(),
    }
np.save(f"ELEC/full_state.npy", state)
np.save(f"ELEC/time.npy", t)
np.save(f"ELEC/n.npy", n)

resuming

# Load check-point data
state = np.load(f"ELEC/full_state.npy", allow_pickle=True).item()
md2.to_variables(state["variables"])
md2.set_variable("Previous_ui", state["Previous_ui"])
md2.set_variable("Previous_vi", state["Previous_vi"])
md2.set_variable("Previous_ul", state["Previous_ul"])
md2.set_variable("Previous_vl", state["Previous_vl"])
t = float(np.load("ELEC/time.npy"))
n = int(np.load("ELEC/n.npy"))

please I have another question for the contact. can you please tell me how I can fix the following warning?

Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-7.97854
Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-7.78622
Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-5.08899
Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-8.39849
Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-8.13806
Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-8.19618
Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-6.72738
Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-6.95948
Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-8.36482
Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-7.37768
Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-8.49248
Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-5.68726
Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-7.71223
Level 1 Warning in getfem_contact_and_friction_common.cc, line 47: Inverted element !-5.89841

Thanks in advance,
Best.

10

These warnings are basically an indication that your Newton solution is diverging.

To avoid this situation consider:

  • Performing smaller load steps
  • Trying different augmentation parameter if your contact model is based on an augmented Lagrangian
  • Using line search or changing the line search settings
  • Changing the side of master/slave
  • Changing the order of the Lagrange multiplier field
  • …
1 Like