Second elasticity tensor does not show minor symmetries

General
1

Ahoy!

I need the fourth order second elasticity tensor, i.e. the derivative of the second Piola-Kirchhoff stress tensor (S) w.r.t. the right Cauchy-Green strain tensor (C). I’m projecting each component into a function space and saving the result at a point in a txt file. There follows a minimal working example. The point is that the minor symmetries do not show up, which is wrong, because both S and C are symmetric tensors in their own rights. Does anyone know what is the issue?

# Minimal working example for saving elasticity tensor at a point

from dolfin import *
import ufl_legacy as ufl
import numpy as np
import os

L, H, W = 1.0, 0.2, 0.3
mesh = BoxMesh(Point(0,0,0), Point(W,H,L), 10,10,10)

lower_facet = CompiledSubDomain("near(x[1], 0)")
left_facet = CompiledSubDomain("near(x[2], 0)")
right_facet = CompiledSubDomain("near(x[2], L)", L=L)
volume_1 = CompiledSubDomain("x[2]<=L*0.5", L=L)
volume_2 = CompiledSubDomain("x[2]>L*0.5", L=L)

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundary_markers.set_all(0)
lower_facet.mark(boundary_markers, 2)
left_facet.mark(boundary_markers, 3)
right_facet.mark(boundary_markers, 6)

dx = Measure("dx", domain=mesh, metadata={"quadrature_degree": 2})
ds = Measure("ds", domain=mesh, subdomain_data=boundary_markers)

# Constitutive model

class Neo_Hookean:

    def __init__(self, E, v):

        self.E = E
        self.v = v
        self.mu = Constant(self.E/(2*(1+self.v)))
        self.lmbda = Constant(self.v*self.E/((1+self.v)*(1-2*self.v)))

    def strain_energy(self, C):

        I1_C = ufl.tr(C)
        I2_C = ufl.det(C)
        J = ufl.sqrt(I2_C)
        
        psi_1 = (self.mu/2)*(I1_C - 3)

        ln_J = ufl.ln(J)
        psi_2 = -(self.mu*ln_J)+((self.lmbda*0.5)*((ln_J)**2))

        return psi_1+psi_2

    def first_piolaStress(self, u):

        I = Identity(3)
        F = grad(u)+I
        C = (F.T)*F
        C = variable(C)
        
        W = self.strain_energy(C)
        S = 2*diff(W,C)

        return F*S
    
    def second_elasticityTensor(self, u):

        I = Identity(3)
        F = grad(u)+I
        C = (F.T)*F
        C = variable(C)
        
        W = self.strain_energy(C)
        S = 2*diff(W,C)

        return diff(S, C)

constitutive_model = Neo_Hookean(1E6,0.3)
traction_vectors = {2: Constant([0.0, 0.0, 0.0]), 6: Constant([0.0, 0.0, 1E5])}

# Sets the forms

V = VectorFunctionSpace(mesh, "Lagrange", 2)
u_trial = TrialFunction(V)
v = TestFunction(V)
u_solution = Function(V)

bc = [DirichletBC(V, Constant((0.0, 0.0, 0.0)), boundary_markers, 3)]

external_work = 0.0

internal_work = inner(constitutive_model.first_piolaStress(u_solution),grad(v))*dx

for physical_group, traction in traction_vectors.items():
    external_work += dot(traction,v)*ds(physical_group)

# Solver
residual_form = internal_work-external_work
residual_derivative = derivative(residual_form , u_solution, u_trial)

Res = NonlinearVariationalProblem(residual_form, u_solution, bc, J=residual_derivative)
solver = NonlinearVariationalSolver(Res)
solver.solve()

# Gets the fourth order tensor
W_index = FunctionSpace(mesh, "CG", 1)
C_array = np.zeros((3,3,3,3))
for i in range(3):
    for j in range(3):
        for k in range(3):
            for l in range(3):
                C_array[i,j,k,l] = project(constitutive_model.second_elasticityTensor(
                u_solution)[i,j,k,l], W_index, form_compiler_parameters={
                "quadrature_degree": 2})(Point(0.5*W, 0.5*H, 0.5*L))

txt_file = open(os.getcwd()+"//tests//minimal_working_examples_forum//dS_dC_array.txt", "w")
txt_file.write(str(C_array))
txt_file.close()
2

As recently discussed in Issue with dolfinx.fem.Function.interpolate - #2 by dokken, projecting discontinuous data in to continuous spaces can introduce artifacts:

3

Is it possible just to evaluate the elasticity tensor in a Gauss point? Or is there any other work around?

I added the following code to save the second Piola-Kirchhoff stress tensor

S_array = np.zeros((3,3))
for i in range(3):
    for j in range(3):
        S_array[i,j] = project(constitutive_model.second_piolaStress(
        u_solution)[i,j], W_index, form_compiler_parameters={
        "quadrature_degree": 2})(Point(0.5*W, 0.5*H, 0.5*L))

txt_file = open(os.getcwd()+"//tests//minimal_working_examples_forum//S_stress.txt", "w")
txt_file.write(str(S_array))
txt_file.close()

yet this works and the S file is indeed symmetric

4

Use a quadrature space for the tensor.

5

Does the whole solution need to be defined on a quadrature space? If not, do I project the tensor into a quadrature space? If so, wouldn’t it fall onto the same drawback of projection? Or does the quadrature space simply evaluate the tensor at the Gauss points?

6

Projection onto a quadrature space would only evaluate the data at the relevant quadrature points.

See for instance