I have some questions regarding VectorFunctionSpace(): in the parameters I can define which requirements the base functions have to fulfill, e.g. "CG". When I choose a function to be solved on CG2, I expect the solved function to be continuous after the first differentiation. Is this behavior implemented in FEniCS? When I differentiate the displacement function "u" by building the strain tensor while “u” is defined on CG2, the projection on a new VectorFunctionSpace(mesh, 'DG', 1, dim=3) produces some questions.
As I understand it, the displacement function "u" is CG2 so its differentiation "strain" should still be continuous. Does projecting a continuous function on a discontinuous function space destroy its continuity or simple does not enforce it? When I plot the result of strain along some elements I get a graph which is not continuous at all but I can't draw any conclusions from this because I don't know if this is because CG2 elements don't behave this way or if project() is destroying continuity. Is there any other way to prove or disprove if CG2 elements behave as described above?
some minimal python example (should run without editing on FEniCS 1.3.0):
from dolfin import *
E, nu, nuu = 3.8e6, 0., 0.38
mu, muu, lambada = E/(2.*(1.+nu)), E/(2.*(1.+nuu)), (E*nuu) / ((1.+nuu)*(1.-2.*nuu))
mesh = BoxMesh(.0,.0,.0, .3,.1,.1, 3,1,1)
Space = VectorFunctionSpace(mesh, 'CG', 2, dim=3)
du = TrialFunction(Space)
delu = TestFunction(Space)
u = Function(Space)
plot_Space = VectorFunctionSpace(mesh, 'DG', 1, dim=3)
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0, 10.e-6) and on_boundary
left = Left()
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.3, 10.e-6) and on_boundary
right = Right()
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
ds = Measure('ds')[boundaries]
bcs = DirichletBC(Space, Constant((0.,0.,0.)), boundaries, 1)
traction = Expression(('0.','0.','trac'), trac = -200/(0.1*0.1))
i,j,k,l = indices(4)
delta = Identity(3)
epsilon = as_tensor(1./2.*(du[i].dx(j)+du[j].dx(i)), (i,j))
stiffness = as_tensor(lambada*delta[i,j]*delta[k,l] + muu*(delta[i,k]*delta[j,l]+delta[i,l]*delta[j,k]), (i,j,k,l))
sigma = as_tensor(stiffness[i,j,k,l]*epsilon[k,l], (i,j))
VOL = sigma[i,j] * (1./2.*(delu[i].dx(j)+delu[j].dx(i))) * dx
BDR = traction[j] * delu[j] * ds(2)
solve(VOL == BDR, u, bcs)
#----POST-PROCESSOR----
import numpy as np
vec = as_tensor((1.,0.,0.))
eps = as_tensor(1./2.*(u[i].dx(j)+u[j].dx(i)) * vec[j], (i,))
epsmesh = project(eps, plot_Space)
EPS = np.array([])
X = np.linspace(0., 0.3, (3*300)+1)
for x in X:
EPS = np.append(EPS, epsmesh((x,0.,0.))[0])
#----PLOT----
import matplotlib.pyplot as pylab
pylab.ylabel('strain')
pylab.xlabel('length')
pylab.grid(True)
pylab.plot(X, EPS, marker=',', markersize=9, linestyle='-', linewidth=2)
pylab.show()
