Gradient projection - Shapes do not match

General
1

Hey, there.

I am currently working on obtaining stresses in my FEM code. To do this, I need to compute gradient of displacements, u, and obtain certain components of the matrix. u is a 2x1 vector, where u = u(x,y), meaning grad(u) must be a 2x2 matrix.

u is working properly, I can plot my displacements with no problem, but I’m struggling to obtain this components of the gradient. I add a simple small code which, I believe, should be enough to verify my problem. Please feel free to ask for more details about my code if needed.

V = VectorFunctionSpace(mesh, 'CG',1)
u = TrialFunction(V)
solve (a==L,u,bcs)

gradu = nabla_grad(u) #in theory, should be a 2x2 matriz
uxp = project(gradu[0,0],V)
ux = uxp.compute_vertex_values(mesh)
uyp = project(gradu[1,1],V)
uy = uyp.compute_vertex_values(mesh)
uxyp = project(gradu[0,1],V)
uxy = uyp.compute_vertex_values(mesh) 
uyxp = project(gradu[1,0],V)
uyx = uyp.compute_vertex_values(mesh) 

sigmaxx = lmbda * (ux + uy) + 2*mu*ux 
sigmayy = lmbda * (ux + uy) + 2*mu*uy
sigmaxy = mu * (uxy + uyx)

The problem I’m facing is:

Shapes do not match: <Argument id=139735787611328> and <Indexed id=139735787816816>.
Traceback (most recent call last):
  File "Dirichlet_Robin_M2.py", line 145, in <module>
    uxp = project(gradu[0,0],V)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/projection.py", line 129, in project
    L = ufl.inner(w, v) * dx
  File "/usr/lib/python3/dist-packages/ufl/operators.py", line 158, in inner
    return Inner(a, b)
  File "/usr/lib/python3/dist-packages/ufl/tensoralgebra.py", line 147, in __new__
    error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))
  File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shapes do not match: <Argument id=139735787611328> and <Indexed id=139735787816816>.

Thanks in advance.

2

This is not a reproducible example, as

  1. u is a TrialFunction and can thus not be sent into solve
  2. missing import statements and variable definitions.

I would suggest you do something like:

from dolfin import *
mesh = UnitSquareMesh(10,10)
V = VectorFunctionSpace(mesh, 'CG',1)
u= Function(V)

then

which clearly means that you need the space you project gradu into should be a tensor function space, not a vector space.
It should Also be a DG space, as the derivatives are only continuous in the interior of a cell, not cross element.

3

Hello @dokken , thank you for your answer.

  1. I don’t understand why u shouldn’t be sent into solve. In fact, as I mentioned, this gradient implementation is created to obtain stress, but if I do not want stresses, the code works just fine and I can obtain u_x and u_y without any problem. So, this shouldn’t be the problem, I believe. But I would like to learn if there is any issue with this.

  2. I do understand that V needs to be a TensorFunctionSpace if I want to project a gradient. So, I did the change:

gradu = nabla_grad(u) #in theory, should be a 2x2 matriz
V = TensorFunctionSpace(mesh, 'DG', 1)
uxp = project(gradu[0,0],V)
ux = uxp.compute_vertex_values(mesh)

But still I’m facing the same shape error. I did exactly the same implementation with antiplane problem, meaning u had only one component, and working with the gradient was just fine.

Additionally, I’m adding a minimal (but longer) code. There would be no problem in adding full codes, also to generate the mesh. I just feel it’s just too much.

Thanks.


msh2xdmf('malla_modo2.msh',2,directory='.')

mesh, mesh_function, association_table = import_mesh(
    prefix='malla_modo2', 
    dim=2,
    directory='.')

print(mesh)

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

 # Boundary definitions, Im not adding just not to make a giant code 

disp_x0 = Expression('x[0] * (-0.5 * x[0] / pow(x[0]*x[0]+x[1]*x[1],0.5) - pow(x[0]*x[0]+x[1]*x[1],0.5) * ((2 * std::log(pow(x[0]*x[0]+x[1]*x[1],0.5)) + 1) * cos(2 * atan2(x[1], x[0])) + (-2 * atan2(x[1], x[0]) + 2 * pi) * sin(2 * atan2(x[1], x[0])) + 2 * std::log(pow(x[0]*x[0]+x[1]*x[1],0.5)) - 1) / (4 * pi)) / pow(x[0]*x[0]+x[1]*x[1],0.5) - x[1] * (0.5 * x[1] / pow(x[0]*x[0]+x[1]*x[1],0.5) + pow(x[0]*x[0]+x[1]*x[1],0.5) * ((2 - 2 * cos(2 * atan2(x[1], x[0]))) * (2 * atan2(x[1], x[0]) - 2 * pi) - (2 * std::log(pow(x[0]*x[0]+x[1]*x[1],0.5)) + 1) * sin(2 * atan2(x[1], x[0]))) / (4 * pi)) / pow(x[0]*x[0]+x[1]*x[1],0.5)',degree=3)
disp_x1 = Expression('x[0] * (0.5 * x[1] / pow(x[0]*x[0]+x[1]*x[1],0.5) + pow(x[0]*x[0]+x[1]*x[1],0.5) * ((2 - 2 * cos(2 * atan2(x[1], x[0]))) * (2 * atan2(x[1], x[0]) - 2 * pi) - (2 * std::log(pow(x[0]*x[0]+x[1]*x[1],0.5)) + 1) * sin(2 * atan2(x[1], x[0]))) / (4 * pi)) / pow(x[0]*x[0]+x[1]*x[1],0.5) + x[1] * (-0.5 * x[0] / pow(x[0]*x[0]+x[1]*x[1],0.5) - pow(x[0]*x[0]+x[1]*x[1],0.5) * ((2 * std::log(pow(x[0]*x[0]+x[1]*x[1],0.5)) + 1) * cos(2 * atan2(x[1], x[0])) + (-2 * atan2(x[1], x[0]) + 2 * pi) * sin(2 * atan2(x[1], x[0])) + 2 * std::log(pow(x[0]*x[0]+x[1]*x[1],0.5)) - 1) / (4 * pi)) / pow(x[0]*x[0]+x[1]*x[1],0.5)',degree=3)

r = Constant((1.0, 1.0)) 
g = Constant((0.0, 0.0))
f = Constant((0.0, 0.0))

disp_null = Constant(0)

bcd0 = DirichletBC(V.sub(0),  disp_x0, boundary3, method='topological')
bcd1 = DirichletBC(V.sub(1),  disp_x1, boundary3, method='topological')

bcd2 = DirichletBC(V.sub(0),  disp_null, boundary1, method='topological')
bcd3 = DirichletBC(V.sub(1),  disp_null, boundary1, method='topological')

bcd4 = DirichletBC(V.sub(0),  disp_x0, boundary4, method='topological')
bcd5 = DirichletBC(V.sub(1),  disp_x1, boundary4, method='topological')

bcd6 = DirichletBC(V.sub(0),  disp_x0, boundary5, method='topological')
bcd7 = DirichletBC(V.sub(1),  disp_x1, boundary5, method='topological')

bcs = [bcd1, bcd2, bcd3, bcd4, bcd5, bcd6, bcd7, bcd0]

ds = Measure('ds', domain = mesh, subdomain_data = Boundary_markers)
dx = Measure('dx', domain = mesh, subdomain_data = Boundary_markers)

u = TrialFunction(V)
v = TestFunction(V)

L = inner(f,v)*dx + inner(g,v)*ds
a = inner(grad(u), grad(v)) * dx + (r[0]*v[0]*u[0]+r[1]*v[1]*u[1]) * ds(1)

A = assemble(a)
b = assemble(L)

u = Function(V, name='Displacement')
solve (a==L,u,bcs)

gradu = nabla_grad(u) #in theory, should be a 2x2 matriz
V = TensorFunctionSpace(mesh, 'DG', 1)
uxp = project(gradu[0,0],V)
ux = uxp.compute_vertex_values(mesh) #vector composed by numbers
uyp = project(gradu[1,1],V)
uy = uyp.compute_vertex_values(mesh) #vector composed by numbers
uxyp = project(gradu[0,1],V)
uxy = uyp.compute_vertex_values(mesh) #vector composed by numbers
uyxp = project(gradu[1,0],V)
uyx = uyp.compute_vertex_values(mesh) #vector composed by numbers
x = mesh.coordinates()[:,0]
y = mesh.coordinates()[:,1]
sigmaxx = lmbda * (ux + uy) + 2*mu*ux
4

As you haven’t supplied a complete code, I cannot really guide you any further than saying that a TrialFunction should not be input into the solve(a==L,u), but a function u whose dofs you want to populate.

I misread, i thought you were trying to extract all of gradu.

Since you only want a single component a normal function space is sufficient, ie

gradu = nabla_grad(u) #in theory, should be a 2x2 matriz
V = FunctionSpace(mesh, 'DG', 1)
uxp = project(gradu[0,0],V)

Note that this is not well-defined as the derivative of a continuous function is discontinuous, and you will get different values for each node is shared with multiple cells

5

I obtained stress without problem. Thanks again @dokken