Set integration rule to closed/vertex

General
1

I want to set the integration rule to closed for all terms of my problem and not only the mass matrix. I have a MWE on the unit interval with P1 elements and a term f = x

from dolfinx import mesh, fem
import ufl
from mpi4py import MPI

nx = 2
domain = mesh.create_unit_interval(MPI.COMM_WORLD, nx)
V = fem.FunctionSpace(domain, ("CG", 1))

x = ufl.SpatialCoordinate( domain )
function = x[0]

v = ufl.TestFunction(V)
b = function*v*ufl.dx

linear_form   = fem.form(b)
b = fem.petsc.create_vector(linear_form)
fem.petsc.assemble_vector(b, linear_form)

print( "b = \n", b.getValues(range(nx+1)) )

Exact integration gives b = [0.04166667, 0.25, 0.20833333] which is what this code outputs.
Closed integration gives b = [0.0, 0.25, 0.25]

I tried setting the differential like it is mentioned in Mass lumping in time dependent problem :

dxL = dx(scheme='vertex', degree=1,
         metadata={'representation': 'quadrature',
                   'degree': 1})

but this outputs an error on my machine:

AttributeError: 'numpy.ndarray' object has no attribute 'integral_type'

I am using dolfinx 0.4.1 downloaded through apt

2

Thanks for reporting this.
I’ve made a fix for this at: Fix interval integrals with vertex rules by jorgensd · Pull Request #562 · FEniCS/ffcx · GitHub
In general, I would use the following syntax:

from dolfinx import mesh, fem
import ufl
from mpi4py import MPI

nx = 2
domain = mesh.create_unit_interval(MPI.COMM_WORLD, nx)
V = fem.FunctionSpace(domain, ("CG", 1))

x = ufl.SpatialCoordinate(domain)
function = x[0]

v = ufl.TestFunction(V)
dx = ufl.Measure(
    "dx", metadata={"quadrature_rule": "vertex"})
b = function*v*dx

linear_form = fem.form(b)
b = fem.petsc.create_vector(linear_form)
fem.petsc.assemble_vector(b, linear_form)

print("b = \n", b.getValues(range(nx+1)))
1 Like
3

Thanks for fixing this @dokken. I tried it out using the docker nightly build and it works for line and triangle elements. It does not work for quadrilateral elements, which is what I ultimately need. Here is a MWE with f = xy

from dolfinx import mesh, fem
import ufl
from mpi4py import MPI

N = 1
domain = mesh.create_unit_square(MPI.COMM_WORLD, N, N,
                                 mesh.CellType.quadrilateral
                                 )
V = fem.FunctionSpace(domain, ("CG", 1))

x = ufl.SpatialCoordinate(domain)

v = ufl.TestFunction(V)
dx = ufl.Measure("dx",
                 #metadata={"quadrature_rule": "vertex"}
                 )
b = x[0]*x[1]*v*dx

linear_form   = fem.form( b )

b = fem.petsc.create_vector(linear_form)
fem.petsc.assemble_vector(b, linear_form)

print( "b = \n", b.getValues( range(b.size) ) )

Exact integration gives b = [0.02777778, 0.05555556, 0.05555556, 0.11111111] which is what this code prints. Closed integration gives b = [0.0, 0.0, 0.0, 0.25].

When I uncomment the line #metadata={"quadrature_rule": "vertex"} , I get the following error:

UnboundLocalError: local variable 'points' referenced before assignment
4

Another PR has been opened: Add quads and hexes for vertex scheme by jorgensd · Pull Request #564 · FEniCS/ffcx · GitHub

1 Like
5

I encountered the same problem for boundary integrals. Here is a MWE:

from dolfinx import mesh, fem
import ufl
from mpi4py import MPI
import numpy as np

N = 1
domain = mesh.create_unit_square(MPI.COMM_WORLD, N, N)
V = fem.FunctionSpace(domain, ("CG", 1))
v = ufl.TestFunction(V)

x = ufl.SpatialCoordinate(domain)

ds = ufl.Measure(
        "ds",
        metadata={
            #"quadrature_rule":"vertex",
            "quadrature_rule":"custom",
            "quadrature_points":np.array([[0.0], [1.0]]),
            "quadrature_weights":np.array([1.0 / 2.0, 1.0 / 2.0]),
            }
        )

b = x[0]*v*ds

linear_form   = fem.form( b )

b = fem.petsc.create_vector(linear_form)
fem.petsc.assemble_vector(b, linear_form)

print( "b = \n", b.getValues( range(b.size) ) )

Uncommenting the vertex quadrature rule line and commenting the custom quadrature rules lines gives an AssertionError. I am working with the nightly version.

6

Fix proposed at: Support facet integrals with vertex scheme by jorgensd · Pull Request #566 · FEniCS/ffcx · GitHub

1 Like