How to apply quadratic element type

1

Hello,

I recently found a code for the cantilever beam with uniform distributed load problem. It works well with CG element type. I want it run with quadratic element type. So I tried to edit some codes inside.
The code is here:

from future import print_function
from dolfin import *
from fenics import *
import matplotlib.pyplot as plt
import math

L = 25.
H = 1.
Nx = 250
Ny = 10
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, “crossed”)

def eps(v):
return sym(grad(v))

E = Constant(1e5)
nu = Constant(0.3)

lmbda = Enu/(1+nu)/(1-2nu)
mu = E/2/(1+nu)
lmbda = 2mulmbda/(lmbda+2*mu)

def sigma(v):
return lmbdatr(eps(v))Identity(2) + 2.0mueps(v)

rho_g = 1e-3
f = Constant((0, -rho_g))

fe = VectorElement(“Quadrature”, cell=mesh.ufl_cell(), degree = 1, quad_scheme = “default”)
V = VectorFunctionSpace(mesh, fe, degree = 1)
du = TrialFunction(V)
u_ = TestFunction(V)
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)

plot(1e3*u, mode=“displacement”)
plt.title(“2D elastic bending model”)
plt.xlabel(“length”)
plt.ylabel(“bending”)
plt.show()

print(“Maximal deflection:”, -u(L,H/2.)[1])
print(“Beam theory deflection:”, float(3rho_gL4/2/E/H3))

Vsig = TensorFunctionSpace(mesh, “DG”, degree=0)
sig = Function(Vsig, name=“Stress”)
sig.assign(project(sigma(u), Vsig))
print(“Stress at (0,H):”, sig(0, H))

error = abs(-u(L,H/2.)[1] - float(3rho_gL4/2/E/H3))
print(error)

But the an err occurred:

runcell(0, ‘/Users/pigdtt12345/Desktop/Bill/FEniCS/Training/2D_elasticity_simple.py’)
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Symmetric part of tensor with rank != 2 is undefined.
Traceback (most recent call last):

File “/Users/pigdtt12345/Desktop/Bill/FEniCS/Training/2D_elasticity_simple.py”, line 34, in
a = inner(sigma(du), eps(u_))*dx

File “/Users/pigdtt12345/Desktop/Bill/FEniCS/Training/2D_elasticity_simple.py”, line 25, in sigma
return lmbdatr(eps(v))Identity(2) + 2.0mueps(v)

File “/Users/pigdtt12345/Desktop/Bill/FEniCS/Training/2D_elasticity_simple.py”, line 14, in eps
return sym(grad(v))

File “/Users/pigdtt12345/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/operators.py”, line 318, in sym
return Sym(A)

File “/Users/pigdtt12345/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/tensoralgebra.py”, line 457, in new
error(“Symmetric part of tensor with rank != 2 is undefined.”)

File “/Users/pigdtt12345/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/log.py”, line 172, in error
raise self._exception_type(self._format_raw(*message))

UFLException: Symmetric part of tensor with rank != 2 is undefined.

Could some please help me with it? I just want it work with “Quadrature” element type. I am sorry if I make some format error here. I am a new member with this qa forum.
Thanks a lot in advance.

2

If what you want is to use quadratic shape functions, I think you just need to set degree=2.

4

Hi,
Thanks for your help. I have tried it before.
I directly tried

V = VectorFunctionSpace(mesh, “Quadrature”, degree = 2)

but the error comes out too. The error said: Missing quad_scheme in quadrature element. So I checked online to solve how to fix it, and I found this page that has same issue.

https://fenicsproject.org/qa/9635/quadrature-element/

But when I try it, it has the exception as the initial problem as “Symmetric part of tensor with rank != 2 is undefined.”

So here I am totally lost what I can do next to fix this bug. Sorry that I am beginner and this is a bit out of my knowledge.
Thanks in advance for your help.

6

Additionally for using quadrature elements, you might have to change the representation for the form compiler. The default is uflacs. This can be done using something like:

parameters["form_compiler"]["representation"] = 'quadrature'
# And to suppress `depecration` warnings 
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)

see here

7

Hello,
you are mixing “quadratic” interpolation (implictly Lagrange P2 elements) with “quadrature” elements which correspond to fields evaluated at quadrature points only and from which one cannot compute gradients and so on.
For quadratic elements, just use V = VectorFunctionSpace(mesh, "CG", 2)

3 Likes
8

Ah! I see! Thank you so much!

9

Hello,

I put it at the top and the error occurs again. It said, “Symmetric part of tensor with rank != 2 is undefined.”. But I think I mixed the quadratic interpolation with quadrature elements. But thank you for your help very much too! Sorry for taking up your time.

11

Yeah. I just figured that you wanted to use quadratic and not quadrature elements. Both

deg = 2
V = VectorFunctionSpace(mesh, "CG", deg)

and

W = VectorElement("CG", mesh.ufl_cell(), deg)
V = FunctionSpace(mesh, W)

are equivalent!!

2 Likes
12

Yeah, because actually I have just got familiar with FEniCS and FEM. I am sorry I did not make my question clear. Thank you for your added up information. This also helps me a lot.