Gradient projection - Shapes do not match

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)

ux = uxp.compute_vertex_values(mesh)
uy = uyp.compute_vertex_values(mesh)
uxy = uyp.compute_vertex_values(mesh)
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>
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>.


1 Like

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.

1 Like

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)
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 = assemble(a)
b = assemble(L)

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

V = TensorFunctionSpace(mesh, 'DG', 1)
ux = uxp.compute_vertex_values(mesh) #vector composed by numbers
uy = uyp.compute_vertex_values(mesh) #vector composed by numbers
uxy = uyp.compute_vertex_values(mesh) #vector composed by numbers
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
1 Like

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)


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

1 Like

I obtained stress without problem. Thanks again @dokken

Hello everyone. I have a similar problem. I can project some terms of the solution, but the grad(u) projections don’t work.

W = FunctionSpace(mesh, ‘DG’, order)
WW = VectorFunctionSpace(mesh, ‘DG’, order )
WWW = TensorFunctionSpace(mesh, ‘DG’, order )

u1p = project(u[0], W)
sigh12p = project(sig[0,1],W) #to plot in paraview

I getting the follows error:

Traceback (most recent call last):
File “/home/ex1.py”, line 134, in
^^^^^^^^^^^^^^^^^^^
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/projection.py”, line 132, in project
A, b = assemble_system(a, L, bcs=bcs,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py”, line 436, in assemble_system
assembler = cpp.fem.SystemAssembler(A_dolfin_form, b_dolfin_form, bcs)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

*** -------------------------------------------------------------------------
*** Error: Unable to assemble system.
*** Reason: expected a linear form for L.
*** Where: This error was encountered inside SystemAssembler.cpp.
*** Process: 0

*** DOLFIN version: 2019.2.0.13.dev0

Can someone help me please? My mininall code is:

'''
This code solves a mixed FEM for the Navier-Stokes equations
in symmetric-pseudostress-velocity-vorticity formulation
'''
# <----- preamble ----->
from fenics import *
import numpy as np
from numpy import linalg as LA
import sympy as sp # required by symbolic expressions
from mshr import *

# create classes for defining boundary
def gamma(x,on_boundary):
return on_boundary

# <----- model parameters and auxiliary symbolic expressions ----->
Id   = Identity(2)
x, y = sp.symbols('x[0] x[1]')

# <--------------------------------------------------------------------->
# <----- Data parameters ----->
nu =   1.0

# <----- boundary condition  ----->
ue_c1 = pow(x,2)*pow((x-1.),2)*sp.sin(y)
ue_c2 = 2.*x*(x-1.)*(2*x-1.)*sp.cos(y)

ue        = Expression((sp.printing.ccode(ue_c1), sp.printing.ccode(ue_c2) ), degree = 5)
# <----- RHS  ----->
ff1 = 4*pow(x, 4)*pow(x - 1.0, 3)*pow(sp.sin(y), 2) + 4*pow(x, 3)*pow(x - 1.0, 4)*pow(sp.sin(y), 2) \
- 2.0*pow(x, 3)*pow(x - 1.0, 3)*(2*x - 1.0)*pow(sp.sin(y), 2)	\
+ 2.0*pow(x, 3)*pow(x - 1.0, 3)*(2*x - 1.0)*pow(sp.cos(y), 2) 	\
+ 1.0*pow(x, 2)*pow(x - 1.0, 2)*sp.sin(y) - 2.0*pow(x, 2)*sp.sin(y) 	\
- 4.0*x*(2*x - 2.0)*sp.sin(y) - 2.0*pow(x - 1.0, 2)*sp.sin(y) \
- 3.1415926535897931*sp.exp(3.1415926535897931*y)*sp.sin(3.1415926535897931*x)

ff2 =  4.0*pow(x, 3)*pow(x - 1.0, 3)*sp.sin(y)*sp.cos(y) 	\
+ 6.0*pow(x, 3)*pow(x - 1.0, 2)*(2*x - 1.0)*sp.sin(y)*sp.cos(y) \
+ 6.0*pow(x, 2)*pow(x - 1.0, 3)*(2*x - 1.0)*sp.sin(y)*sp.cos(y) 	\
- 32.0*pow(x, 2)*pow(x - 1.0, 2)*pow(x - 0.5, 2)*sp.sin(y)*sp.cos(y) 	\
+ 2.0*x*(x - 1.0)*(2*x - 1.0)*sp.cos(y) - 8.0*x*sp.cos(y) - 8.0*(x - 1.0)*sp.cos(y) 	\
- 4.0*(2*x - 1.0)*sp.cos(y) + 3.1415926535897931*sp.exp(3.1415926535897931*y)*sp.cos(3.1415926535897931*x)

ff        = Expression(( sp.printing.ccode(ff1), sp.printing.ccode(ff2) ), degree = 5)

# <--------------------------------------------------------------------->
# <----- solver  ----->
# <------ lists to hold errors, meshsize, etc. ------->
err_rel0    = 10.0
tol         = 1e-6
itermax     = 20

NN = [4, 8, 16, 32]
for i in NN:
domain_vertices = [Point( 0.0, 0.0),
Point( 1.0, 0.0),
Point( 1.0, 1.0),
Point( 0.0, 1.0)]
domain  = Polygon(domain_vertices)
mesh		= generate_mesh(domain, i)

Om			= assemble( Constant(1.0)*dx(mesh) )
n				= FacetNormal(mesh)
s_e			= as_vector((-n[1],n[0]))
h     	= mesh.hmax()
nelem 	= mesh.num_cells()

# <---------Polynomial order degree ----->
order = 0

# <-------- define function space -------->
Hsig  = VectorElement("BDM",mesh.ufl_cell(), order+1)  # tensor BDM_k+1
Hu    = VectorElement("DG" ,mesh.ufl_cell(), order)     # vector P_k discontinuous
L     = FiniteElement("R"  ,mesh.ufl_cell(), 0)          # Lagrange multiplier space
Xh    = FunctionSpace(mesh,MixedElement([Hsig,Hu,L]))
dim   = Xh.dim()

# <----- Newton's method ----->
iter    = 0.
err_rel = err_rel0
# <----- Initialization ----->
u0 = Constant((0,0))
# <----- Kernel ----->
while (err_rel > tol) and (iter< itermax) :
# <-----  1. solve step ----->
# <----- define trial and test functions ----->
(sig, u, l1) = TrialFunctions(Xh)
(tau, v, t1) = TestFunctions(Xh)

# <---- involved forms ----->
A   = (1/nu)*inner(dev(sig),dev(tau))*dx
An  = (1/nu)*inner(outer(u0,u) + outer(u,u0),dev(tau))*dx
B   = dot(u,div(tau))*dx
Bt  = dot(v,div(sig))*dx
M1  = l1*tr(tau)*dx
M2  = t1*tr(sig)*dx
AA  = An + A + B + Bt + M1 + M2

F1  = 1/nu*inner(outer(u0,u0),dev(tau))*dx
F2  =  dot(tau*n,ue)*ds
Ff2 =  - dot(ff,v)*dx
FF  = F1 + F2 + Ff2

# <---- function to solve for ---->
sol = Function(Xh)
solve(AA == FF, sol, solver_parameters = {'linear_solver':'mumps'} ) #solver_parameters = {'linear_solver':'mumps,umfpack''}
( sig, u, l1) = sol.split()
solaux = np.reshape(sol.vector()[:], (Xh.dim(),1))
coeffm = solaux
if iter==0:
err_rel = err_rel0
else:
err_rel = LA.norm(coeffm - coeff_old)/LA.norm(coeffm)
# <----- 2. update data ----->
coeff_old = coeffm
u0        = u
iter = iter + 1.0

print('it.Newton',iter)
print('number of elements: ', nelem, 'DOFs:',dim,  'err_rel', err_rel )

# <------ post-processing ------->
W  = FunctionSpace(mesh, 'DG', order)
WW = VectorFunctionSpace(mesh, 'DG', order )
WWW = TensorFunctionSpace(mesh, 'DG', order )

u1p 		= project(u[0], W)
sigh12p = project(sig[0,1],W) #to plot in paraview


Your u is in a Vector DG-0 space, hence grad_u will always be zero.