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.

1 Like
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.

1 Like
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
1 Like
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

1 Like
5

I obtained stress without problem. Thanks again @dokken

6

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

gradu = nabla_grad(u)
#u1xp = project(gradu[0,1], W)
#uxp = project(gradu[0], WW)
grad_up = project(gradu, WWW)

I getting the follows error:

Traceback (most recent call last):
File “/home/ex1.py”, line 134, in
grad_up = project(gradu, WWW)
^^^^^^^^^^^^^^^^^^^
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

*** https://fenicsproject.discourse.group/

*** 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
author: Saulo Rodrigo Medrado
'''
# <----- 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
	
		gradu = nabla_grad(u)		
		#u1xp 		= project(gradu[0,1], W)
		#uxp 		= project(gradu[0], WW)
		#grad_up = project(gradu, WWW)

# <----- export graphics ----->		
#u_file    = File("Data_Paraview_2D/approx_u.pvd")    << u
#gradu_file    = File("Data_Paraview_2D/approx_grad_u.pvd")    << grad_u
7

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

1 Like