Error: "Can't add expressions with different shapes."?

1

Hi!

In the code below u is a scalar function of x and y and U = (0, 0, u(x,y)). The code

dot(grad(U) + grad(U).T, n) gives the error

ufl.log.UFLException: Can't add expressions with different shapes.

I am not quite able to figure out why. Any suggestions would be helpful. Here’s the MWE:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt


from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp



# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


import numpy
import meshio
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *


class boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary


   class disk(SubDomain):
    def inside(self, x, on_boundary):
        return True


mesh = generate_mesh(Circle(Point(0, 0), 3), 50)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
boundary().mark(boundary_markers, 2)
disk().mark(surface_markers, 1)


ds = Measure('ds', subdomain_data=boundary_markers)
dx = Measure('dx', subdomain_data=surface_markers)
n = FacetNormal(mesh)

W1 = FunctionSpace(mesh, "Lagrange", 1)
n = FacetNormal(mesh)
G , mu = 1, 0.1
u_D=Constant(0.0)

bc = DirichletBC(W1, u_D, boundary_markers, 2)

# Define variational problem
u = TrialFunction(W1)
v = TestFunction(W1)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx
L = -f*v*dx

# Compute solution
u = Function(W1)
solve(a == L, u, bc)

# Calculating shear stress S
Ux=0.0
Uy=0.0
Uz=u
U = as_vector((Ux, Uy, u))
Sn = mu/2 * dot(grad(U) + grad(U).T, n)
2

Given \vec{U} \in \mathbb{R}^m and \vec{x} \in \mathbb{R}^d then \nabla \vec{U} \in \mathbb{R}^{m \times d}.

Your mesh looks 2D, so d=2. And you’ve defined \vec{U} such that m=3, so \nabla \vec{U} \in \mathbb{R}^{3 \times 2}. Which is verified:

print(grad(U).ufl_shape)
>>> (3, 2)

So grad(U) + grad(U).T is not a valid operation.

Given your system setup, you should carefully redesign grad or consider a different definition of U to match your problem.

2 Likes
3

Hi Nate, thanks for your reply. I completely agree with your explanation. The code that I posted was copied from another script where the mesh was 2D and was created in GMSH. In that case, while transporting the mesh, I added the following commands to get rid of this exact error:

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    #if prune_z:  # if you remove the "#", the "invalid operation" message crops up
        #out_mesh.prune_z_0()
    return out_mesh

I don’t know how to implement a similar set of command when creating a 2D mesh in the script itself.

4

prune_z_0 in this context will remove all of the z coordinate entries if they’re the same value (essentially you have a 2D mesh).

It’s up to you whether you want to solve a truly 2D problem. Or a 2D manifold in a 3D space problem.

5

My problem is on a 2D manifold in a 3D space. I am a little confused as to how to relay that information to the script itself while creating a in-built 2D geometry. Are there specific commands to be plugged in so that FEniCS recognizes that it’s dealing with a 2D structure embedded in a 3D space?

I guess that whenever I am creating a 2D geometry (in the script itself), FEniCS by default assumes that it is embedded in a 2D space.

I apologize for so many questions. Since I usually import structures from gmsh, I am not quite adept yet at creating in-built structures and defining boundaries in the script itself. But to create MWEs, I need to master these techniques.

6

You can check the geometrical dimension of your mesh by:
print(mesh.geometry().dim()). If this dimension is equal to 2, dolfin thinks of the problem as a pure 2D problem, if it is 3, it assumes that it is embedded in a 3D space (manifold). This is determined by the dimension of the points sent into the mesh constructor. If you are using mshr to generate the mesh (the disk above), it is going to set the geometrical dimension to 2.

1 Like
7

Hi dokken! Thanks for the explanation. Here is what I observed with a couple of geometries created through the script and gmsh.

In gmsh, all 2D geometries are by default embedded in a 3D space.
If the geometry is created through the script, it depends on the dimension of the points as you said.