Mechanical equilibrium of a single finite element

General
1

Dear community. I would like to know if there is way to compute the nodal forces on the nodes of one finite element that mechanically balance the mechanical deformation of such an element, i.e., such that the virtual power of such nodal forces equals the strain energy of the element. In other words, I wonder how I could implement the following integral,

\mathbf{f}^{n}=\int_{\mathcal{B}e}\sigma^t \cdot\nabla N^n dV

where \nabla N^n dV Is the gradient of the shape function of the node n; and \mathcal{B}e the volume domain of the single element.

To this end, I guess that I need a way to locate an individual finite element within a larger mesh, to obtain the derivative of the nodal shape functions, and to compute the integral on that element. Then, do the same for the rest of nodes of the element.

Many thanks.

2

Assuming uh is the deformation vector

def sigma(u):
    return lmbda * ufl.nabla_div(u) * Identity(len(u)) + 2*mu*epsilon(u)

v = TestFunction(V)
f = assemble(inner(sigma(uh), grad(v))*dx)
3

Thank you, dokken.
How can I get the displacement vector uh?
I guess it must contain only the solved nodal degrees of freedom. I have tried u.vector().array() using FEniCSx, but I get the error “‘petsc4py.PETSc.Vec’ object is not callable”.
Also, which Function from dolfinx.fem.assemble should I use?

4

I would assume that you are solving a linear elasticity problem to get your deformation field uh, otherwise how would you define \sigma?

Im not sure what you mean by this.

You did not state what version of DOLFIN you were running, so I was assuming the legacy version. The equivalent for DOLFINx would be

def sigma(u):
    return lmbda * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*epsilon(u)

v = ufl.TestFunction(V)
f = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(ufl.inner(sigma(uh), ufl.grad(v))*ufl.dx))
print(f.get_local())
5

No, I work with finite strains and the deformation gradient (\mathbf{F}). More precisely, it is not the Cauchy stress tensor the one I want to use, but the Eshelby one, but this does not affect the implementation. I define u as a Function from a Function Space, and then, \mathbf{F} = \mathbf{I} + \text{ufl.grad}(\mathbf{u}). With this implementation, I aim to compute the material (configurational) nodal forces.
Hence, is this u solution function what I should use in the assemble you proposed?

Since you wrote uh, I am not sure if that is simply my Function \mathbf{u} already solved, or just the nodal d.o.f.s ordered in a different manner.

That for the assemble works, thanks. However, with print(f.get_local()) I get the following error: AttributeError: 'petsc4py.PETSc.Vec' object has no attribute 'get_local'

6

Sorry, mixing legacy and dolfinx, replace

With f.array()

7

As you didn’t post a u in your initial code, I used uh to illustrate the numerical displacement function (stemming from one process or another). uh should therefore be replaced with your u

8

When I do this:

i.e.,
f = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(ufl.inner(eshelby_material_stress(F), ufl.grad(v_u))*dx))
print(f.array())

I still get this error: TypeError: 'numpy.ndarray' object is not callable

9

print(f.array) then.

This error states that f.array is a numpy array, and does not have a call function, i.e. f.array().

10

This seems to work. Thank you for the explanations.

Then, the only thing I still need to wrap up this implementation would be to get a vector with the coordinates of the nodes for which the nodal forces in the vector f.array correspond to. I would appreciate your help again here.

11

V.tabulate_dof_coordinates(), see

12

Thanks. This works.

And for the coordinates of the d.o.f.s in the deformed configuration?

13

You can just sum the coordinate of the unperturbed state with the deformation (as you can get the displacement at every node in the exact same way).

14

Now everything seems to be working. Many thanks.

15

One additional concern related to the nodal forces:
How could I write the forces on the nodes, i.e., the vector \mathbf{f}, to my .xdmf file to visualize them on Paraview?

16

Rewrite this to:

def sigma(u):
    return lmbda * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*epsilon(u)

v = ufl.TestFunction(V)
f = dolfinx.fem.Function(V)
dolfinx.fem.petsc.assemble_vector(f.vector, dolfinx.fem.form(ufl.inner(sigma(uh), ufl.grad(v))*ufl.dx))

and then save f to file as you would usually do with DOLFINx