Problem with orthonormal vectors and outer product

variational formulation
1

Hello everyone,

I’m having an issue when I try to use an orthonormal basis of vectors. I will be using these vectors as fiber directions in a nonlinear solid mechanics problem. For example, let’s consider a 2D domain, where I have defined two orthonormal directions \vec f and \vec s. Then I can build a tensor using these vectors as:

G = \gamma_f \vec f \otimes \vec f + \gamma_s \vec s \otimes \vec s

If \gamma_f = \gamma_s = 1, G should be equal to the identity. The problem is that I’m not getting that with Fenics, which is causing me trouble because, for example, if I don’t put any loads to a problem I still get a non-zero residual. I track the error to this. Here is a MEW:

import dolfin as d

mesh = d.RectangleMesh(d.Point(-5, -5), d.Point(5, 5), 5, 5)

V = d.VectorFunctionSpace(mesh, 'DG', 2)

N1 = d.Expression(('cos(x[0])/sqrt(pow(cos(x[0]), 2) + pow(sin(x[1]), 2))', 'sin(x[1])/sqrt(pow(cos(x[0]), 2) + pow(sin(x[1]), 2))'), element = V.ufl_element())
N2 = d.Expression(('-sin(x[1])/sqrt(pow(cos(x[0]), 2) + pow(sin(x[1]), 2))', 'cos(x[0])/sqrt(pow(cos(x[0]), 2) + pow(sin(x[1]), 2))'), element = V.ufl_element())

N1 = d.project(N1, V)
N2 = d.project(N2, V)
d.File('testN1.pvd') << N1
d.File('testN2.pvd') << N2
   
T = d.TensorFunctionSpace(mesh, 'DG', 2)    
I_df = d.outer(N1, N1) + d.outer(N2, N2)
F = d.dot(I_df, I_df)

I_df_ = d.project(I_df, T)
F = d.project(F, T)

print(I_df_.vector().get_local().reshape([-1,2,2]))


import numpy as np
for i in range(N1.vector().get_local().reshape([-1,2]).shape[0]):
    n1 = N1.vector().get_local().reshape([-1,2])[i]
    n2 = N2.vector().get_local().reshape([-1,2])[i]
    if np.abs(np.dot(n1,n2)) > 1e-12:   # Check that the vectors are perpendicular
        print(i)
    if np.abs(np.dot(n1, n1) - 1) > 1e-12:
        print(i)
    if np.abs(np.dot(n2, n2) - 1) > 1e-12:
        print(i)
    I_np = np.outer(n1,n1) + np.outer(n2,n2)

The output is

[[[ 1.07626026e+00 -5.71851283e-16]
  [-5.71851283e-16  1.07626026e+00]]

 [[ 7.36428348e-01 -1.67397705e-17]
  [-1.67397705e-17  7.36428348e-01]]

 [[ 7.55745073e-01 -1.67397705e-17]
  [-1.67397705e-17  7.55745073e-01]]

 ...

 [[ 1.07723900e+00 -1.52486776e-16]
  [-1.52486776e-16  1.07723900e+00]]

 [[ 1.04033445e+00  1.93728740e-17]
  [ 1.93728740e-17  1.04033445e+00]]

 [[ 1.07945641e+00  2.66692601e-17]
  [ 2.66692601e-17  1.07945641e+00]]]

which is not the identity anywhere. I checked that the vectors at the nodes were orthonormal in the nodes and they are. I’m using a DG element of second degree because that’s what I need to use in further applications. If I refine the mesh it does get better, however, I won’t be able to do that in a bigger problem. I’m thinking that when fenics projects is evaluating the vector functions in the quadrature points, where, because of the interpolation, they might not be orthonormal. Is there a way to enforce orthonormality everywhere in the functions? I’m using docker with the development version of fenics.

Thank you in advance!

2

I would consider instead defining \mathbf{N}_1 and \mathbf{N}_2 using a UFL SpatialCoordinate. See the following example:

import dolfin as d

mesh = d.RectangleMesh(d.Point(-5, -5), d.Point(5, 5), 5, 5)

#V = d.VectorFunctionSpace(mesh, 'DG', 2)
#N1 = d.Expression(('cos(x[0])/sqrt(pow(cos(x[0]), 2) + pow(sin(x[1]), 2))', 'sin(x[1])/sqrt(pow(cos(x[0]), 2) + pow(sin(x[1]), 2))'), element = V.ufl_element())
#N2 = d.Expression(('-sin(x[1])/sqrt(pow(cos(x[0]), 2) + pow(sin(x[1]), 2))', 'cos(x[0])/sqrt(pow(cos(x[0]), 2) + pow(sin(x[1]), 2))'), element = V.ufl_element())
x = d.SpatialCoordinate(mesh)
from ufl import sin, cos, sqrt
N1 = d.as_vector((cos(x[0])/sqrt(pow(cos(x[0]), 2) + pow(sin(x[1]), 2)),
                  sin(x[1])/sqrt(pow(cos(x[0]), 2) + pow(sin(x[1]), 2))))
N2 = d.as_vector((-sin(x[1])/sqrt(pow(cos(x[0]), 2) + pow(sin(x[1]), 2)),
                  cos(x[0])/sqrt(pow(cos(x[0]), 2) + pow(sin(x[1]), 2))))

#N1 = d.project(N1, V)
#N2 = d.project(N2, V)
#d.File('testN1.pvd') << N1
#d.File('testN2.pvd') << N2
   
#T = d.TensorFunctionSpace(mesh, 'DG', 2)    
I_df = d.outer(N1, N1) + d.outer(N2, N2)
F = d.dot(I_df, I_df)

#I_df_ = d.project(I_df, T)
#F = d.project(F, T)

# Verify that $L^2$ norm of difference between F and identity is machine zero;
# because equality holds pointwise for UFL definition of N1, N2, we don't need
# high-order quadrature, even though the integrand involves trig functions:
d.dx = d.dx(metadata={"quadrature_degree":1})
print(sqrt(d.assemble(((F - d.Identity(2))**2)*d.dx)))

The problem with what you’re doing now is that there are approximation errors associated with interpolating and projecting things into finite element spaces.

4

Thanks for the answer, it indeed fixes the issue for the MWE. However, in my real problem, I have as input a fiber field that can’t be computed analytically. I understand that my problem came from interpolating and projecting things into the finite element spaces. So, is there a way to tell Fenics to compute I_df at the nodes (where I know the vectors are orthonormal) and then interpolate this tensor instead of interpolating the vectors and computing the outer product at the quadrature point?

5

You could create a tensor-valued UserExpression subclass, then interpolate an instance of that.