### Tensor function projections strange behaviour

372

views

-1

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:

I found a possible workaround... not too elegant nor fast, but functional:

```
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

### 1 Answer

0

I hope it could help you.

Kind regards,

Leo

```
# -*- 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
boundaries.set_all(0)
Topo = top()
Topo.mark(boundaries, 1)
Bottom=bottom()
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})
plot(u)
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')
interactive()
```

written
12 months ago by
Nicolás Barnafi

Please login to add an answer/comment or follow this question.

`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

?