I’m trying to create a piecewise vector function in FEniCS; I’m using the legacy version. I was able to run the piecewise scalar variant of it. The tutorial is that of a plate (L x H); this is taken from here.

We want to solve div \sigma=f. The code is

```
from dolfin import *
import fenics as fe
import matplotlib.pyplot as plt
import numpy as np
L = 2.
H = 2.
Nx = 250
Ny = 10
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, "crossed")
x = SpatialCoordinate(mesh)
def eps(v):
return sym(grad(v))
E = Constant(1e5)
nu = Constant(0.3)
model = "plane_stress"
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
lmbda = 2*mu*lmbda/(lmbda+2*mu)
def sigma(v):
return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
#Defining the function space and the functions.
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2, dim=2)
du = TrialFunction(V)
u_ = TestFunction(V)
rho_g1 = Constant(1e-3) #constant to be used later
rho_g2 = Constant(1). #constant to be used later
#Testingpiecewise fn defn
p1 = Expression(("0","c1"), c1=rho_g1, degree=2, domain=mesh)
p2 = Expression(("c2","0"), c2=rho_g2, degree=2, domain=mesh)
f = Expression('x[0] <= 0.2 + DOLFIN_EPS ? p1 : p2', p1=p1, p2=p2, degree=2)
##Testing vector function
#f = Expression(("c2","0"), c2=rho_g2, degree=2, domain=mesh)
print("shape of f=",f.ufl_shape)
print("shape of u_=",u_.ufl_shape)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx
def left(x, on_boundary):
return near(x[0], 0.)
bc = DirichletBC(V, Constant((0.,0.)), left)
u = Function(V, name="Displacement")
solve(a == l, u, bc)
```

There are two sections to the above code; both of them relate to defining the loading f. In the code above, we have commented “#Testing vector function” and left “#Testingpiecewise fn defn” uncommented. The above tests a piecewise vector function as input for the force. The code “print(“shape of f=”,f.ufl_shape)” returns null. If you change the comment from #Testing vector function" to the other, then the code runs perfectly.