Loss of Information During Projection

3 months ago by
Hello all,

Normally I use projection function very frequently to visualize and check the code in each iteration. However, I noticed that projection causes much of loss information in the functionspace. For example I have (3,3) tensonfunctionspace and I am calculating  its normal tensor by   $n=\frac{A}{sqrt\left(A:A\right)}$n=Asqrt(A:A)   . Normally its norm should equal to 1. Then  I am checking it by N= project(n,V) where V is a first order Lagrange TensorFunctionSpace. When I check specific points to see what does it calculates such as N(0,0,0), it gives a tensor which's norm is much higher than 1.

So, is it only visual thing in projection or does it really calculates it higher than 1 also for UFL function ? Does the projection function calculates the function values on nodes or quadrature points ? What is the exact meaning if I got unexpected results when I check values by projection ?

Thanks a lot.

Community: FEniCS Project

Here is the sample code on elasticity short code. I calculated normal tensor from strain and getting its norm in a specific point at (10,1,1). The norm should be equal to 1.0 normally, but it calculates as 1.19. There is %19 error. The norm also is not equal to 1.0 in any other different location.
from __future__ import print_function
from fenics import *
%matplotlib inline

# Scaled variables
L = 10; W = 5
mu = 1
rho = 0.01
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 5, 5)
V = VectorFunctionSpace(mesh, 'Lagrange', 2)
Vepspn = TensorFunctionSpace(mesh, "Lagrange",1, shape=(3,3))
# Define boundary condition
tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

# Define strain and stress

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
    #return sym(nabla_grad(u))

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = Function(V)
d = 3  # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
T = Constant((0, 0, 0))


Pi=Psi*dx-dot(f, u)*dx - dot(T, u)*ds

F=derivative(Pi,u,v)  #To get linear form ?? fenics-ufl.pdf p:26


eps = project(epsilon(u),Vepspn)

normal = project(eps / sqrt(inner(eps,eps)),Vepspn)
norm_sqr = 0.0
sample = normal(10,1,1)
for i in range(len(sample)):
    norm_sqr = norm_sqr + (sample[i])**2
norm = sqrt(norm_sqr)


written 3 months ago by Christian  

2 Answers

3 months ago by
K D  
The set of tensors that you have defined is not a tensor space (roughly they are a manifold).  They cannot be represented by standard shape functions that are designed for spaces, and the error is to be expected.

To see the issue, consider the simpler case of unit vectors with linear interpolation in 1D (the vectors can be in 2D or 3D, the dependence is only on one spatial variable).  If you have unit vectors at two adjacent nodes, the interpolation between them will not satisfy that the vector field is unit anywhere (except in the trivial constant case).  People in liquid crystals and micromagnetics have to spend a lot of effort to get around this issue when they apply FEM with a unit vector constraint.

Hello Dayal,

Thanks for your response.

Actually I could not get the meaning of the first paragraph. Why they are not tensor space ? But I do projection to map them to tensor space. Is not it meaningful ? What have should I done actually ?

About the second paragraph of your comment, the same situation is valid also for stress calculation. Stresses cannot be same for adjacent nodes. What does FEniCs do then ? Taking the avarage of two nodes or does it override one to other ?


written 3 months ago by Christian  
An important requirement of a vector space is that if I take 2 vectors and add them, the result should also be in the space.  If I add 2 unit vectors, I don't get another unit vector.  Hence unit vectors are not a vector space.  The analogous thing is true for your normalized tensors, i.e. adding 2 tensors normalized as you have does not give another normalized tensor.

It's simply not possible to accurately represent such objects using standard FEM function spaces.  How to deal the issue often varies with the physical problem.

About the stress question: stress tensors can form a space, because it is meaningful to have any stress tensor.  They may need to be symmetric, but symmetric tensors fortunately are a tensor space.  So projecting them is pretty easy.  As the other response stated, L2 projections are used.
written 3 months ago by K D  
Another way how to see it is that, assuming your $A$A is polynomial, resulting  $\frac{A}{\left|A\right|}$A|A|   is not a polynomial. Hence it can't be exactly represented by a polynomial FE field.
written 3 months ago by Jan Blechta  
That's a nice and simple explanation.
written 3 months ago by K D  
3 months ago by
project is  $L^2$L2 -orthogonal projection. (Without a minimal working example it is hard to say more about your specific problem.)
Please login to add an answer/comment or follow this question.

Similar posts:
Search »