Tensor function projections strange behaviour

13 months ago by
Hi everyone, I have been doing some 2D elasticity and need to obtain sigma. For it, lets just deal with getting grad(u) numerically. Apparently projecting on different meshes raises the error. Here is a MWE:

N = 10
mesh = UnitSquareMesh(N,N)
V = VectorFunctionSpace(mesh, 'CG', 1)
u = Function(V)
VT = TensorFunctionSpace(mesh, 'CG', 1)
gradu = project(grad(u),VT)

# This one raises the error
VT = TensorFunctionSpace(UnitSquareMesh(N,N), 'CG', 1)
gradu = project(grad(u), VT)


I found a possible workaround... not too elegant nor fast, but functional:
N = 10
mesh = UnitSquareMesh(N,N)
V = FunctionSpace(mesh, 'CG', 1)
DV = FunctionSpace(mesh, 'DG', 0)
ux = project(u[0], V)
uy = project(u[1], V)
uxdx = project(ux.dx(0), DV)
uxdy = project(ux.dx(1), DV)
uydx = project(uy.dx(0), DV)
uydy = project(uy.dx(1), DV)
gradu = as_matrix([[uxdx, uxdy], [uydx, uydy]])
VT = TensorFunctionSpace(UnitSquareMesh(N,N), 'DG', 0)
gradu = project(gradu, VT)​
 In fact, directly projecting grad(u) also fails, which forced me to project partial derivatives too. 

Best regards.

Community: FEniCS Project
Could you clarify 'if instead of creating u for this I create it to hold the solution values'? It would be helpful to post the precise code that leads to the error so it can be tested unambiguously. 
written 13 months ago by Garth Wells  
Ok, this helped me find a simpler way of exposing the issue. I am updating it, should be up in some minutes. Thanks!
written 13 months ago by Nicolás Barnafi  
Garth, I have had a lot of trouble with tensor functions lately (which have been solved so far so no big deal) and I believe that the core of the issue is the difference between tensor elements and row wise vector elements. I understand how this is directly influenced by RT and BDM being constructed row wise, but maybe using tensor elements as also row wise but projected into the row wise canonical basis can solve the issue (i.e \( \vec{\mathbb P_k} = \vec e_0 \mathbb P_k + \vec e_1 \mathbb P_k \)) . It would be great to hear some comments on this. Best regards!
written 13 months ago by Nicolás Barnafi  
Hi, the MWE seems fine, no errors, solutions match to 1e-14?
written 12 months ago by David Nolte  
Hi David, what version are you using? Just rechecked (2017.1.0) and I still get the error.
written 12 months ago by Nicolás Barnafi  
Oh, I tried it with 2016.2.0.
written 12 months ago by David Nolte  
Ok, confirmed the error with 2017.1, JIT fails. The error.log in the jitfailure-ffc_form_* folder which is created says that some const double J_c0 are being redeclared, and indeed the cpp file has duplicated code... This also happens for scalar functions and their gradients, not only for TensorFunctionSpaces. Don't know what the intended behavior is. Maybe file a bug report?

On the other hand, if you do it in 2 steps, it works:
- project the gradient onto a space on the same mesh
- interpolate the result onto a different mesh
written 12 months ago by David Nolte  
Thanks for rechecking, maybe I had a compilation bug going on super hidden. I thought about interpolating but I need the projection for a numerical convergence study. Thanks!
written 12 months ago by Nicolás Barnafi  

1 Answer

12 months ago by

I hope it could help you.

Kind regards,


# -*- coding: utf-8 -*-
Created on Fri Jul 28 06:59:32 2017

@author: leocosta
from dolfin import*
N = 10
mesh = UnitSquareMesh(N,N)
boundaries = FacetFunctionSizet(mesh, 1)

V = VectorFunctionSpace(mesh, 'CG', 1)

u= TrialFunction(V)
v = TestFunction(V)
E = 1.0
nu = 0.3

lmbda, mu = Constant(E*nu/((1.0 + nu )*(1.0-2.0*nu))) , Constant(E/(2*(1+nu)))

def epsilon(v):
    return 0.5*(grad(v) + grad(v).T)
def sigma(u):
    return lmbda*tr(epsilon(u))*Identity(2) + 2.0*mu*epsilon(u)

class bottom(SubDomain):
    def inside(self,x,on_boundary):
        tol = 1E-14
        return abs(x[1]) < tol and on_boundary

class top(SubDomain):
    def inside(self,x,on_boundary):
        tol = 1E-14
        return abs(x[1]-1.) < tol and on_boundary

Topo = top()
Topo.mark(boundaries, 1)

Bottom.mark(boundaries, 2)

bc = DirichletBC(V, (0.,0.), boundaries, 2)
b = Constant((0., 0.))
s= Constant((0.0 ,1.E-5))

ds = Measure("ds", domain=mesh, subdomain_data=boundaries)

# Variational form
a = inner(sigma(u), grad(v))*dx
L = + dot(b,v)*dx \
    - dot(s,v)*ds(1)

u = Function(V)

solve(a == L, u, bcs= bc,
      solver_parameters={"linear_solver": "lu"},
      form_compiler_parameters={"optimize": True})

d = u.geometric_dimension() 
# Plot stress
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)  # deviatoric stress
von_Mises = sqrt(3./2*inner(s, s))
V2 = FunctionSpace(mesh, 'P', 1)
von_Mises = project(von_Mises, V2)
plot(von_Mises, title='Stress intensity')

Thanks, but the answer is not related to the question, it regards an issue with pointers to different meshes and projections. This solutions has already been presented at the docs.
written 12 months ago by Nicolás Barnafi  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »