A situation regarding the non integrability of quadrilateral mesh

General
1

Dear all!
I encountered a mistake.
When I calculate an integral ：

fv=(div(grad(u_n_1)))**2*dx

When the grid is divided into triangles, code is OK.
But the grid is divided into quadrilateral , will get :

ValueError: ReferenceGrad can only wrap a reference frame type!

This leaves me very puzzled.

code are as follows:

from mpi4py import MPI
import numpy as np
from petsc4py import PETSc
import ufl 
from dolfinx import default_scalar_type, fem, la 
from dolfinx.cpp.fem import compute_integration_domains 
from dolfinx.mesh import (
    GhostMode,
    compute_incident_entities,
    create_box,
    create_rectangle,
    create_submesh,
    create_unit_cube,
    create_unit_interval,
    create_unit_square,
    locate_entities,
    locate_entities_boundary,
    meshtags,
)
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
from dolfinx import default_real_type, fem, io, mesh,default_scalar_type
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from ufl import (CellDiameter, FacetNormal, TestFunction, TrialFunction, avg,
                 conditional, div, dot, dS, ds, dx, grad, gt, inner, outer,
                  TrialFunctions,TestFunctions,curl,cross)
import ufl
from dolfinx.io import gmshio,XDMFFile
from basix.ufl import element,mixed_element 
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, 
                                 create_vector, set_bc,create_matrix,)
import time
import os
import dolfinx 
startT=time.time()

def coeff_expr(x): 
    values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType) 
    values[0]=x[0] 
    values[1]=x[1] 
    return values  
 
def coeff_expr2(x): 
    values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType) 
    values[0]=x[0]**2
    values[1]=x[1]**2
    return values  


comm = MPI.COMM_WORLD 
msh = create_unit_square(comm, 2,  2, ghost_mode=GhostMode.none,cell_type=mesh.CellType.triangle)   
#mesh.CellType.quadrilateral  
#CellType.triangle

U_space=fem.functionspace(msh,("RT",2))

u_n_1=fem.Function(U_space)
u_n_2=fem.Function(U_space)
u_D = fem.Function(U_space)

comm=msh.comm
fv=(div(grad(u_n_1)))**2*dx
fv=fem.form(fv)
FV=np.sqrt(comm.allreduce(fem.assemble_scalar(fv), MPI.SUM))
print('FV is',FV)

Any comments will be appreciated！Thanks!

2

This seems like a bug in UFL, as I can reproduce it with this pure ufl code:

from ufl import (
    FunctionSpace,
    Mesh,
    Coefficient,
    div,
    dot,
    dx,
    grad,
    triangle,
    quadrilateral, classes
)
from ufl.finiteelement import FiniteElement
from ufl.pullback import contravariant_piola, identity_pullback
from ufl.sobolevspace import H1, HDiv

cell = quadrilateral
RT = FiniteElement("Raviart-Thomas", cell, 1, (2,), contravariant_piola, HDiv)
domain = Mesh(FiniteElement("Lagrange", cell, 1, (2,), identity_pullback, H1))
space = FunctionSpace(domain, RT)

u = Coefficient(space)
J = dot(div(grad(u)), div(grad(u))) * dx

from ufl.algorithms import compute_form_data
from ufl.algorithms import compute_form_data
compute_form_data(J,
        do_apply_function_pullbacks=True,
        do_apply_integral_scaling=False,
        do_apply_geometry_lowering=False,
        preserve_geometry_types=(classes.Jacobian,),
        do_apply_restrictions=False,
        do_append_everywhere_integrals=False, )

Seems like there is an issue in apply_function_pullbacks.
I’ve made an issue at:

3

Proposed fix submitted at:

Tested locally with

from mpi4py import MPI
from dolfinx.mesh import create_unit_square, GhostMode, CellType
import numpy as np
from dolfinx import fem
from ufl import div, grad, dx, inner, SpatialCoordinate, as_vector


comm = MPI.COMM_WORLD 

ct = CellType.quadrilateral

msh = create_unit_square(comm, 2,  2, ghost_mode=GhostMode.none,cell_type=ct)   

U_space=fem.functionspace(msh,("RT",2))

u =fem.Function(U_space)

def f(x):
    alpha = 2.0
    beta = -3.0
    return (alpha * x[0]**2, -beta * x[0]*x[1]**2)

u.interpolate(f)

comm=msh.comm
fv=inner(div(grad(u)), div(grad(u)))*dx
fv=fem.form(fv)

x = SpatialCoordinate(msh)
u_ex = div(grad(as_vector(f(x))))
diff = div(grad(u)) - u_ex

FV=np.sqrt(comm.allreduce(fem.assemble_scalar(fv), MPI.SUM))
print('FV is',FV)

FV_ref = np.sqrt(comm.allreduce(fem.assemble_scalar(fem.form(inner(u_ex,u_ex)*dx)), MPI.SUM))
print("FV_ref is",FV_ref)

diff_form = inner(diff, diff) * dx
error = np.sqrt(comm.allreduce(fem.assemble_scalar(fem.form(diff_form)), MPI.SUM))
print("Error is", error)

yielding

FV is 5.291502622129174
FV_ref is 5.29150262212918
Error is 1.4378792021287633e-14
1 Like