Projection on Quadrilateral

1

Hello everyone,

I tried to compared results of 2D triangular and quadrilateral element results in mixed formulation. It seems like projected stress results by triangular elements are smooth but projected stress results by quadrilateral fluctuated. It may be both about the usage of DG0 in the mixed space on quadrilateral mesh or it may be about the projection function which may not be worked well for quadrilateral elements. Does anyone face such a problem ? Is there any suggestion ?

Thanks.

2

Projection is just solving a linear equation inner(u,v)*dx==inner(f,v)*dx, so in itself that should differ much for quadrilaterals and triangles. Could you supply a minimal example?

3

Hello dokken,

Here is the one of the simple example for 2D problem. I extracted a stress component and if you compare the results of T2 and Q1 at the corners (legend setting -50 to 120) there appears weird distribution for quadrilateral element.

From my other observations, it is not about the projection but I could not figure out the reason.

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"

ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

Type ="T2"

inc_load =np.array([ 0. ,0.1,0.2,0.3,0.4, 0.5])

E, nu = 200, 0.49
Element = "CG"

if Type=="T2":
    order_u=2
    mesh = RectangleMesh(Point(0., 0.), Point(1, 1), 64,64, "right/left")

if Type=="Q1":
    order_u=1
    mesh=UnitSquareMesh.create(64, 64, CellType.Type.quadrilateral)

mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

V2 = VectorFunctionSpace(mesh, Element, order_u)
du = TrialFunction(V2)          
v  = TestFunction(V2)            
u  = Function(V2)                

quad_degree=3
metadata = {"quadrature_degree": quad_degree, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)

# Kinematics
d = 3
grad_u = grad(u)
grad_u3 = as_tensor([[grad_u[0, 0], grad_u[0, 1], 0.0 ], [grad_u[1, 0], grad_u[1, 1],0.0], [0.0, 0.0,0.0]])

I = Identity(d)             # Identity tensor
F = I + grad_u3             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor
J=det(F)
C= variable(C)


I1 = tr(C)
I2 = 1/2*((tr(C))**2-tr(C*C))
I3 = det(C)

tol =1e-14

class clamped_boundary_up(SubDomain): #Region definition to be fixed
    def inside(self,x,on_boundary):
        return on_boundary and near(x[1],0,tol)
    
class clamped_boundary_down(SubDomain): #Region definition to be fixed
    def inside(self,x,on_boundary):
        return on_boundary and near(x[1],H,tol)
    

up_bc = clamped_boundary_up()
down_bc = clamped_boundary_down()  
boundary_part = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_part.set_all(0)
up_bc.mark(boundary_part, 2)

Tr = Expression(('pr', '0.0'),degree=1,pr=0)
bc1 = DirichletBC(V2, Tr, up_bc)
bc2 = DirichletBC(V2, Constant((0, 0)), down_bc)
bc = [bc1, bc2]

psi_iso = (lmbda/4)*(I3-1)-((lmbda/2)+mu)*ln((I3)**(0.5))+(mu/2)*(I1-3)

psi = psi_iso


S_ =2*diff(psi_iso,C)
Cauchy = (F*S_*F.T)/J
vert = as_vector((1.0, 0.0, 0.0))
Cauchy_vert = inner(Cauchy,outer(vert,vert))

# Total potential energy
Pi = psi*dxm 

F = derivative(Pi, u, v)

# Compute Jacobian of F
J = derivative(F, u, du)

problem = NonlinearVariationalProblem(F, u, bc, J=J) 
solver = NonlinearVariationalSolver(problem)


file_results = XDMFFile("Result.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

P1_1 = FunctionSpace(mesh, "CG", 1)
sigma = Function (P1_1)

for t in inc_load:
    Tr.pr=t
    solver.solve()
    sigma.assign(project(Cauchy_vert,P1_1,form_compiler_parameters={"quadrature_degree":quad_degree}))
    file_results.write(u, t)
    file_results.write(sigma, t)
4

You are kind of comparing apples and oranges, as you are using a CG-2 space for the triangles, while you are using a CG-1 space for the quadrilaterals. The difference in numbers of degrees of freedom is 33282 to 8450, which is quite different.

As you are using half the number of elements, you also need to slightly increase the quadrature order (to 5 for instance). Using a CG-2 space for both problems gives similar results.