Creating a piecewise vector function

Legacy DOLFIN
1

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.

2

Looks like this was just a syntax issue. When I defined the functions component-wise, the code ran smoothly. The part that was changed was “#Testingpiecewise fn defn”. It was changed to:

#piecewise fn defn, component wise.
p11 = Expression('0', degree=2, domain=mesh)
p12 = Expression('c1', c1=rho_g1, degree=2, domain=mesh)
p21 = Expression('c2', c2=rho_g2, degree=2, domain=mesh)
p22 = Expression('0', degree=2, domain=mesh)
f = Expression(('x[0] <= 0.2 + DOLFIN_EPS ? p11 : p12', 'x[0]<=0.2+DOLFIN_EPS ? p21 : p22'), p11=p11, p12=p12, p21=p21, p22=p22, degree=2, domain=mesh)
1 Like