Serendipity Elements

1

Is it possible to use serendipity-type elements in fenicsx on a quad mesh? I have a code working with standard polynomial elements (“P”), but changing the element family to “S” results in a PETSc error:

[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run 
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(50176059) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 50176059) - process 0
2

Please supply a minimal code example where this is happening.

I’ve tested using S, with the Poisson demo, and it works as expected.

3

My original code was fairly complex, so distilling down a minimal example was non-trivial, but here it is:

import numpy as np
import ufl

from dolfinx import log, plot, fem
from dolfinx.mesh import CellType, create_unit_square, compute_boundary_facets

from mpi4py import MPI
from petsc4py import PETSc

try:
    import pyvista as pv
    import pyvistaqt as pvqt

    have_pyvista = True
    if pv.OFF_SCREEN:
        pv.start_xvfb(wait=0.5)
except ModuleNotFoundError:
    print("pyvista is required to visualise the solution")
    have_pyvista = False
    
have_pyvista = False

# Save all logging to file
log.set_output_file("log.txt")

#
# Next, various model parameters are defined::
def solve_dh(h, dt=1.0e-2, T=8, decoupled=False, k=2):
    domain = create_unit_square(MPI.COMM_WORLD, h, h, CellType.quadrilateral)

    P2 = ufl.FiniteElement("S", domain.ufl_cell(), k)
    H = fem.FunctionSpace(domain, P2)

    x = ufl.SpatialCoordinate(domain)
    t = fem.Constant(domain, 0.0)

    def phi_e(x):
        return np.sin(np.pi * x[0]) * np.cos(np.pi * x[1]) ** 2

    def f_3():
        return 3 * ufl.pi ** 2 * ufl.cos(ufl.pi * x[1]) ** 2 * ufl.sin(ufl.pi * x[0]) \
               - 2 * ufl.pi ** 2 * ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]) ** 2

    # geometric boundary
    def boundary_D(x):
        return np.logical_or(np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1)),
                             np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1)))

    dofs_D = fem.locate_dofs_geometrical(H, boundary_D)

    tdim = domain.topology.dim
    fdim = tdim - 1
    domain.topology.create_connectivity(fdim, tdim)

    phi_e_bc = fem.Function(H)
    phi_e_bc.interpolate(phi_e)

    # Weak statement of the equations
    # terms for Eq. 13
    phi_e = ufl.TrialFunction(H)
    psi = ufl.TestFunction(H)
    # Eq. 13
    F3 = ufl.inner(ufl.grad(phi_e), ufl.grad(psi)) * ufl.dx \
         - f_3() * psi * ufl.dx

    a3 = fem.form(ufl.lhs(F3))
    L3 = fem.form(ufl.rhs(F3))
    A3 = fem.assemble_matrix(a3, bcs=[fem.dirichletbc(phi_e_bc, dofs_D)])
    A3.assemble()
    b3 = fem.assemble_vector(L3)

    phi_e_m1, phi_e_m2 = fem.Function(H), fem.Function(H)
    phi_e = fem.Function(H)

    # numbered by step, not equation
    solver1 = PETSc.KSP().create(domain.comm)
    solver1.setOperators(A3)
    solver1.setType(PETSc.KSP.Type.PREONLY)
    pc1 = solver1.getPC()
    pc1.setType(PETSc.PC.Type.LU)

    # interpolate exact solution at initial time
    t.value = 0.0
    phi_e_m2.interpolate(phi_e)

    t.value += dt
    phi_e_m1.interpolate(phi_e)

    # Prepare viewer for plotting solution during the computation
    if have_pyvista:
        topology, cell_types = plot.create_vtk_topology(domain, domain.topology.dim)
        grid = pv.UnstructuredGrid(topology, cell_types, domain.geometry.x)
        # glyph_grid = pv.UnstructuredGrid(topology, cell_types, mesh.geometry.x)

        grid.point_data["phi_e"] = phi_e.compute_point_values().real
        p = pvqt.BackgroundPlotter(title="Plot", auto_update=True, shape=(2, 2))

        p.subplot(1, 1)
        phi_e_grid = grid.copy()
        phi_e_grid.set_active_scalars("phi_e")
        p.add_mesh(phi_e_grid, clim=[-1, 1])
        #glyph_actor = p.add_mesh(glyphs)
        p.view_xy(True)

        p.subplot(0, 0)
        p.add_mesh(grid, color="k", style="wireframe")
        #glyph_actor = p.add_mesh(glyphs)
        p.view_xy(True)

        p.add_text(f"time: {t}", font_size=10, name="timelabel")

    # time stepping
    print(f"starting timestepping...")
    while t.value < T:
        t.value += dt

        phi_e_bc.interpolate(phi_e)

        A3 = fem.assemble_matrix(a3, bcs=[fem.dirichletbc(phi_e_bc, dofs_D)])
        A3.assemble()

        solver1.setOperators(A3)

        with b3.localForm() as loc_b:
            loc_b.set(0)

        fem.assemble_vector(b3, L3)

        # Apply Dirichlet boundary condition to the vector
        fem.apply_lifting(b3, [a3], [[fem.dirichletbc(phi_e_bc, dofs_D)]])
        b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        fem.set_bc(b3, [fem.dirichletbc(phi_e_bc, dofs_D)])

        r = solver1.solve(b3, phi_e.vector)
        phi_e.x.scatter_forward()

        # Compute L2 error and error at nodes
        error_e_L2 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form((phi_e - phi_e_bc) ** 2 * ufl.dx)), op=MPI.SUM))

        # Update the plot window
        if have_pyvista:
            p.subplot(1, 1)
            p.add_text(f"phi_e L2: {error_e_L2:.3e}", font_size=8, name="phi_e_L2", position=(10, 10))

            p.add_text(f"time: {t.value:.2e}", font_size=10, name="timelabel")

            phi_e_grid.point_data["phi_e"] = phi_e.compute_point_values().real

            p.app.processEvents()

    print(f"completed solve with h = {h}")
    return None


ns = [2, 4, 8, 16, 32, 64]
errors = [solve_dh(h=N, T=2, dt=1.0/(10*N), decoupled=True, k=2) for N in ns]

This is just a poisson problem as well, but with time-stepping. In this example the boundary conditions do not vary in time, but in general they may. My original problem involves a coupled system of which this is the simplest equation.

I have observed that the failure mode is not always the same as what I stated originally. Occasionally, the apparent heap corruption will be caught by the OS (macOS) rather than by PETSc. When it occurs in the computation also varies substantially.

4

It seems like you are using a slightly out of date dolfinx.
I’ve made some changes to your code, including using locate_dofs_topological for the boundary conditions, as you cannot use locate_dofs_geometrical on higher order Serendipity spaces, as the dofs are defined through intergral moments.
Here is my altered version, that runs for me:

import numpy as np
import ufl

from dolfinx import log, fem
from dolfinx import __version__ as version
from dolfinx.common import git_commit_hash
from dolfinx.mesh import CellType, create_unit_square, locate_entities_boundary

from mpi4py import MPI
from petsc4py import PETSc


# Save all logging to file
log.set_output_file("log.txt")

#
# Next, various model parameters are defined::

print(f"Running DOLFINx version {version}, git commit hash {git_commit_hash}")


def solve_dh(h, dt=1.0e-2, T=8, decoupled=False, k=2):
    domain = create_unit_square(MPI.COMM_WORLD, h, h, CellType.quadrilateral)

    P2 = ufl.FiniteElement("S", domain.ufl_cell(), k)
    H = fem.FunctionSpace(domain, P2)

    x = ufl.SpatialCoordinate(domain)
    t = fem.Constant(domain, 0.0)

    def phi_e(x):
        return np.sin(np.pi * x[0]) * np.cos(np.pi * x[1]) ** 2

    def f_3():
        return 3 * ufl.pi ** 2 * ufl.cos(ufl.pi * x[1]) ** 2 * ufl.sin(ufl.pi * x[0]) \
            - 2 * ufl.pi ** 2 * \
            ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]) ** 2

    # geometric boundary
    def boundary_D(x):
        return np.logical_or(np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1)),
                             np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1)))

    facets_D = locate_entities_boundary(
        domain, domain.topology.dim-1, boundary_D)
    dofs_D = fem.locate_dofs_topological(H, domain.topology.dim-1, facets_D)

    tdim = domain.topology.dim
    fdim = tdim - 1
    domain.topology.create_connectivity(fdim, tdim)

    phi_e_bc = fem.Function(H)
    phi_e_bc.interpolate(phi_e)

    # Weak statement of the equations
    # terms for Eq. 13
    phi_e = ufl.TrialFunction(H)
    psi = ufl.TestFunction(H)
    # Eq. 13
    F3 = ufl.inner(ufl.grad(phi_e), ufl.grad(psi)) * ufl.dx \
        - f_3() * psi * ufl.dx

    a3 = fem.form(ufl.lhs(F3))
    L3 = fem.form(ufl.rhs(F3))
    A3 = fem.petsc.assemble_matrix(a3, bcs=[fem.dirichletbc(phi_e_bc, dofs_D)])
    A3.assemble()
    b3 = fem.petsc.assemble_vector(L3)

    phi_e_m1, phi_e_m2 = fem.Function(H), fem.Function(H)
    phi_e = fem.Function(H)

    # numbered by step, not equation
    solver1 = PETSc.KSP().create(domain.comm)
    solver1.setOperators(A3)
    solver1.setType(PETSc.KSP.Type.PREONLY)
    pc1 = solver1.getPC()
    pc1.setType(PETSc.PC.Type.LU)

    # interpolate exact solution at initial time
    t.value = 0.0
    phi_e_m2.interpolate(phi_e)

    t.value += dt
    phi_e_m1.interpolate(phi_e)

    # time stepping
    print(f"starting timestepping...")
    while t.value < T:
        t.value += dt

        phi_e_bc.interpolate(phi_e)
        A3.zeroEntries()
        fem.petsc.assemble_matrix(
            A3, a3, bcs=[fem.dirichletbc(phi_e_bc, dofs_D)])
        A3.assemble()

        solver1.setOperators(A3)

        with b3.localForm() as loc_b:
            loc_b.set(0)

        fem.petsc.assemble_vector(b3, L3)

        # Apply Dirichlet boundary condition to the vector
        fem.apply_lifting(b3, [a3], [[fem.dirichletbc(phi_e_bc, dofs_D)]])
        b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,
                       mode=PETSc.ScatterMode.REVERSE)
        fem.set_bc(b3, [fem.dirichletbc(phi_e_bc, dofs_D)])

        r = solver1.solve(b3, phi_e.vector)
        phi_e.x.scatter_forward()

        # Compute L2 error and error at nodes
        error_e_L2 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(
            fem.form((phi_e - phi_e_bc) ** 2 * ufl.dx)), op=MPI.SUM))

    print(f"completed solve with h = {h}")
    return None


ns = [2, 4, 8, 16, 32]
errors = [solve_dh(h=N, T=2, dt=1.0/(10*N), decoupled=True, k=2) for N in ns]

and returns

Running DOLFINx version 0.3.1.0, git commit hash 2f63364f479899b1500490d9f27869a83d6a5d8e
starting timestepping...
completed solve with h = 2
starting timestepping...
completed solve with h = 4
starting timestepping...
completed solve with h = 8
starting timestepping...
completed solve with h = 16
starting timestepping...
completed solve with h = 32
starting timestepping...
5

Thank you, indeed using locate_dofs_topological appears to avoid any issues with out of bounds memory accesses, but appears to locate the wrong boundary nodes when using "S" elements. ("P" gave me the same results as before).

My dolfinx was out of date by a couple months; I thought I had updated all the fenicsx packages recently, but must have missed dolfinx. Now that I’m up-to-date I find that some functions I had been flying on have moved or changed; in particular, compute_point_values() on Functions, which I used to visualize the solutions in PyVista, appears to have moved. Do you know if there is any replacement?

6

See: Using the method compute_point_values() · Issue #2136 · FEniCS/dolfinx · GitHub

If you can supply a minimal example which shows that it finds the wrong boundary nodes I can have a look.

7

With the updated dolfinx the topological boundaries appear to be identified correctly.

Thank you for the references regarding compute_point_values, interpolation into a Lagrange function space is working well instead. One follow-up question: my mesh is a quad mesh, but visualizing it with PyVista via

topology, cell_types, geometry = plot.create_vtk_mesh(V)
grid = pv.UnstructuredGrid(topology, cell_types, geometry)
p.subplot(0, 0)
p.add_mesh(grid, color="k", style="wireframe")

shows a triangulated mesh. Why would that be?

Finally, when I used 2nd-order serendipity elements on a VectorFunctionSpace defined as

P2 = ufl.VectorElement("S", mesh.ufl_cell(), 2)
P1 = ufl.FiniteElement("P", mesh.ufl_cell(), 1)
TH = ufl.MixedElement([P2, P1])
W = fem.FunctionSpace(mesh, TH)

I encountered a new error when interpolating the boundary:

w_bc = fem.Function(W)
w_bc.sub(0).interpolate(exact_solution.u)
w_bc.sub(1).interpolate(exact_solution.p)

  File "/Users/bosmacs/Repos/dolfinx/python/dolfinx/fem/function.py", line 336, in interpolate
    _interpolate(u, cells)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/functools.py", line 888, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/Users/bosmacs/Repos/dolfinx/python/dolfinx/fem/function.py", line 320, in _interpolate
    self._cpp_object.interpolate(u, cells)
RuntimeError: Interpolation data has the wrong shape.

Is this the wrong way to interpolate a subspace of a mixed space? Is works with the "P" elements I had before.

Thanks again for all your help.

8

This is due to how pyvista handles higher order edges, see: VTK_LAGRANGE_HEXAHEDRON not supported with showing edges · Issue #947 · pyvista/pyvista · GitHub

As you haven’t defined what exact_solution.u is in your question, I cannot help you with the error message.

9

@bosmacs I had a look at interpolating into serendipity, and it works if you interpolate into each component of the space:

import dolfinx
from mpi4py import MPI
import ufl
ct = dolfinx.mesh.CellType.quadrilateral

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, cell_type=ct)
el = ufl.VectorElement("S", cell=mesh.ufl_cell(), degree=2)

V = dolfinx.fem.FunctionSpace(mesh, el)
u = dolfinx.fem.Function(V)
u.sub(0).interpolate(lambda x: x[0])
u.sub(1).interpolate(lambda x: x[1])
print(u.x.array)

V_CG = dolfinx.fem.VectorFunctionSpace(mesh, ("CG", 2))
v = dolfinx.fem.Function(V_CG)
v.interpolate(u)

with dolfinx.io.XDMFFile(mesh.comm, "u.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(v)
10

Thanks, interpolating component-wise does work:

w_bc.sub(0).sub(0).interpolate(lambda x: exact_solution.u(x)[0])
w_bc.sub(0).sub(1).interpolate(lambda x: exact_solution.u(x)[1])
w_bc.sub(1).interpolate(exact_solution.p)

where exact_solution is an instance of

@dataclass
class ExactSolution:
    t: ufl.Constant
    x: ufl.Measure

    def u(self, x):
        return np.array(
            (-(1. - np.cos(2. * np.pi * x[0])) * np.sin(2 * np.pi * x[1]) * np.exp(-self.t.value),
            np.sin(2. * np.pi * x[0]) * (1. - np.cos(2 * np.pi * x[1])) * np.exp(-self.t.value))
        )

    def p(self, x):
        return np.sin(2 * np.pi * x[0]) * np.cos(2 * np.pi * x[1])

It seems like a bug that component-wise interpolation is required for this particular choice of element – should I open an issue for it?

11

Ive already made an issue: Interpolation into vector-Serendipty not working for full space · Issue #2155 · FEniCS/dolfinx · GitHub

12

Does locate_dofs_geometrical work for higher-order Lagrange (CG) elements?

13

Yes, as higher order Lagrange DOFs are defined as point evaluations locate_dofs_geometrical is well defined

1 Like