When using CG degree 1, the result corresponds to degree 2

1

Hello everyone, I’m trying to solve a very simple problem from linear elasticity, where a square mesh of 2x2 elements has a different shear modulus at each element (the shear modulus is piecewise constant, DG degree 0). Then apply a shear to the mesh, and compute the resulting strain.

My problem is that, even though I set the elements to CG with degree 1, I obtain what seems to be the result of degree 2, i.e., instead of linear it looks quadratic. You can appreciate in the attached image the strain plotted (in Paraview) along the diagonal of the mesh. In brown is the linear result I expected, and in orange, it is shown the result from Fenics (quadratic-like).

I would like to know how can I use linear interpolation, and why Fenics is still using quadratic. You can see my code bellow.

Thanks in advance.

Screenshot from 2020-03-30 02-47-40
Screenshot from 2020-03-30 02-47-40646×658 7.67 KB

The code to reproduce this is the following:

from __future__ import print_function
from dolfin import *
import numpy as np
import matplotlib.tri as tri
from mshr import *
import matplotlib.pyplot as plt

mesh = RectangleMesh.create( [Point(0., 0.),Point(2, 2)], [2,2], CellType.Type.quadrilateral) 

nu = Constant(0.3)    

# piecewise constant shear modulus
V = FunctionSpace(mesh, "DG", 0)
G = Function(V)  
G.vector().vec().setValueLocal(0, 1)
G.vector().vec().setValueLocal(1, 6)
G.vector().vec().setValueLocal(2, 0.3)
G.vector().vec().setValueLocal(3, 4)

# isotropic stifness tensor
a = 2*G/(1-nu)    
C = as_matrix([[a,a*nu,0],[a*nu,a,0],[0,0,G]])  
   
def eps(v):
return sym(grad(v))

def strain2voigt(e):
"""e is a 2nd-order tensor, returns its Voigt vectorial representation"""
return as_vector([e[0,0],e[1,1],2*e[0,1]])

def voigt2stress(s):
"""s is a stress-like vector (no 2 factor on last component)returns its tensorial representation"""
return as_tensor([[s[0],s[2]],[s[2],s[1]]])

def sigma(v):
return voigt2stress(dot(C,strain2voigt(eps(v))))   

# use linear shape functions
fem_order = 1

# define problem
f = Constant((0,0))
V = VectorFunctionSpace(mesh, 'CG', fem_order)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx

# apply a shear deformation of eps_xy=0.25 to the whole boundadry
eps_xy = 0.25
def all_boundary(x, on_boundary): return on_boundary
inflow_profile = ('{}*x[1]'.format(eps_xy), '{}*x[0]'.format(eps_xy))                             
dirichelt_shear_BCs  = DirichletBC(V, Expression(inflow_profile, degree=2), all_boundary)    

# solve problem
u = Function(V, name="Displacement")
solve(a == l, u, dirichelt_shear_BCs, form_compiler_parameters={"optimize": True})

file_results = XDMFFile("./output.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)

# compute strain
T = TensorFunctionSpace(mesh, 'CG', fem_order)
strain = project(eps(u), T)
[e11, e12, e21, e22] = strain.split(True)
S = project(e12, FunctionSpace(mesh, 'CG', fem_order))
file_results.write(S, 0.)
2

To visualize other elements than CG 1 in paraview, you Need to use XDMFFile.write_checkpoint, see: Loading xdmf data back in

3

But what I’m trying to visualize is precisely CG 1 so then I shouldn’t need that right?

4

The gradient of a CG 1 Function is a DG 0 function.

5

But why I see quadratic interpolation in Paraview?

6

Which line are you plotting over? I cant reproduce your plots (also you have changed the labels in the code compared to the figure you are showing, so this is clearly not the code to reproduce the behavior).

7

The provided code creates the same quadratic interpolation in the output file, but with the label f_90.

The brown curve is added to paraview manually by me, and corresponds to the solution from another FEM solver when using linear shape functions.

8

To me what one obtains is a 2D function (which is natural as you use a 2D mesh), so you are clearly plotting your function over a line. Which line are you plotting the data over?

Additionally, as I said previously, you should not save the strain as a CG 1 function, but a DG 0 function (using write_checkpoint).

9

I’m plotting along the diagonal (from bottom left to top right).

I should use XDMFFile.write_checkpoint to visualize other elements than linear, but actually I want to visualize linear. Yet, it looks not linear. Either way, I have used write_checkpoint and there’s no difference in the visualization in Paraview, in this case.

I try the following now: I assemble the system with quadratic shape functions and compute the strain with linear shape functions, which is something more consistent (see code below, I set fem_order=2, and later to 1 for computing strains). But I still see that the strain is not linear in this Paraview plot (see below). So it leads me to the following questions, about why the result is not linear: am I using Fenics wrongly? is Fenics doing the right thing but writing the output file wrongly? Or is Paraview doing its own interpolation?

Screenshot from 2020-03-30 17-26-49
Screenshot from 2020-03-30 17-26-49571×640 5.73 KB

The data for this plot is generated with this code (the same as before but changing the degree of the shape functions):

from __future__ import print_function
from dolfin import *
import numpy as np
import matplotlib.tri as tri
from mshr import *
import matplotlib.pyplot as plt

mesh = RectangleMesh.create( [Point(0., 0.),Point(2, 2)], [2,2], CellType.Type.quadrilateral) 

nu = Constant(0.3)    

# piecewise constant shear modulus
V = FunctionSpace(mesh, "DG", 0)
G = Function(V)  
G.vector().vec().setValueLocal(0, 1)
G.vector().vec().setValueLocal(1, 6)
G.vector().vec().setValueLocal(2, 0.3)
G.vector().vec().setValueLocal(3, 4)

# isotropic stifness tensor
a = 2*G/(1-nu)    
C = as_matrix([[a,a*nu,0],[a*nu,a,0],[0,0,G]])  
   
def eps(v):
    return sym(grad(v))

def strain2voigt(e):
    """e is a 2nd-order tensor, returns its Voigt vectorial representation"""
    return as_vector([e[0,0],e[1,1],2*e[0,1]])

def voigt2stress(s):
    """s is a stress-like vector (no 2 factor on last component)returns its tensorial representation"""
    return as_tensor([[s[0],s[2]],[s[2],s[1]]])

def sigma(v):
    return voigt2stress(dot(C,strain2voigt(eps(v))))   

# use linear shape functions
fem_order = 2

# define problem
f = Constant((0,0))
V = VectorFunctionSpace(mesh, 'CG', fem_order)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx

# apply a shear deformation of eps_xy=0.25 to the whole boundadry
eps_xy = 0.25
def all_boundary(x, on_boundary): return on_boundary
inflow_profile = ('{}*x[1]'.format(eps_xy), '{}*x[0]'.format(eps_xy))                             
dirichelt_shear_BCs  = DirichletBC(V, Expression(inflow_profile, degree=2), all_boundary)    

# solve problem
u = Function(V, name="Displacement")
solve(a == l, u, dirichelt_shear_BCs, form_compiler_parameters={"optimize": True})

file_results = XDMFFile("./output.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)

# compute strain
T = TensorFunctionSpace(mesh, 'CG', fem_order-1)
strain = project(eps(u), T)
[e11, e12, e21, e22] = strain.split(True)
S = project(e12, FunctionSpace(mesh, 'CG', fem_order-1))
file_results.write(S, 0.)
10

There are two things that you are ignoring. The first one is that even if you would like to have a linear stress, you do not want a continuous lagrange representation.
Therefore you need to use write_checkpoint. Using a function space that is linear discontinuous you obtain the attached image.

T = TensorFunctionSpace(mesh, 'DQ', fem_order-1)
strain = project(eps(u), T)
[e11, e12, e21, e22] = strain.split(True)
file_results = XDMFFile("./output2.xdmf")
file_results.write_checkpoint(e12, "e12",0.)

DQ
DQ3775×2096 419 KB

If you plot this over the line from 0,0.2 to 2,0.2, you obtain the linear stress you are looking for.

11

Thank you, this solves my problem well enough. But I cannot understand why Fenics keeps using a non-linear interpolation.

12

There is not any nonlinear interpolation. What write_checkpoint does that it writes the degrees of freedom to file, as a finite element function, see xdmf documentation. As you clearly see from plotting over the line I suggested, the solution is linear in directions aligned with the coordinate system.
The basis functions for a linear CG/DQ quadrilateral can be found on slide two in this lecture note.

13

I see your point now clearly, thank you.