P1isoP2 elements

mesh
1

Is it possible to implement “P1isoP2” elements in fenicsx, where one has a mixed element on a grid and a refined version of that same grid?

I’m interested in trying this in the context of solving the Navier-Stokes equations, with linear elements for the pressure on a regular triangular grid, and linear elements for velocity on a refined triangulation. My naïve approach using two meshes didn’t get very far:

mesh = create_unit_square(MPI.COMM_WORLD, h//2, h//2, CellType.triangle)
mesh_refined = create_unit_square(MPI.COMM_WORLD, h, h, CellType.triangle)

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

Which leads to odd but not surprising errors when attempting to construct the forms:

TypeError: '<' not supported between instances of 'Mesh' and 'Mesh'
2

Hello,

I don’t think this is possible yet, but it may be possible in the (quite distant) future once we’ve looked at supporting macro elements (eg Guzmán-Neilan elements).

3

Hello,
Are there any advances on the topic ?

4

P1 iso P2 elements are now supported, you can create them with:

import basix
import basix.ufl

element = basix.ufl.element(basix.ElementFamily.iso, basix.CellType.triangle, 1)

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10, CellType.triangle)
W = fem.FunctionSpace(mesh, element)
2 Likes
5

Hi,
Thank you for your answer, I didn’t know this was available. However in addition to the velocity I would need to define tensors on the same refined grid which would be of order 0. I tried to pass 0 as order in your code but I got an error “Cannot create a continuous order 0 Lagrange basis function”. It seems that it forces the elements to be continuous. Is there any workaround ?

6

This should work:

element = basix.ufl.element(basix.ElementFamily.P, basix.CellType.triangle, 0, discontinuous=true)
2 Likes
7

Did you mean basix.ElementFamily.iso ?

Also this does not provide a tensor but only a scalar, I wanted to plug this into ufl.TensorElement but the family name ‘iso’ does not exist (that’s actually why I thought it was not implemented). I thought using “ufl.finiteelement.elementlist.register_element” would solve the problem but depending on the order of the element, the sobolev space varies from H1 to L2, but the basix family name remains ‘iso’, so I can’t register both order 1 and order 0 in the same script. Did I miss anything ?

8

ufl.TensorElement and ufl.finiteelement.elementlist.register_element are deprecated and will be removed in the next UFL release, so we recommend moving away from using these.

The functions in basix.ufl should be able to do everything that you used to do with these. You can create a tensor element by passing a shape into the element function, for example:

element = basix.ufl.element(
    basix.ElementFamily.P, basix.CellType.triangle, 0,
    discontinuous=true, shape=(2, 2))

There’s no need to use an iso element here as it would be the same as a DP0 element

1 Like
9

Okay thank you ! In that case I guess that ufl.VectorElement should also be replaced by something like

vector = basix.ufl.element(
basix.ElementFamily.iso,basic.CellType.triangle,1,shape=(2,))

?
Concerning the use of basix.ElementFamily.iso, I don’t understand why this is the same as basix.ElementFamily.P when the degree is 0. For me the latter designates the use of Lagrange elements, so if we consider the order 0 this makes only one degree of freedom per triangle. What I want is Lagrange elements of order 0 but on the refined grid, so this makes 3 degrees of freedom per triangle on the coarse mesh. Is there any way to get that ?

10

Yes, that’s exactly how to create a vector element.

In Basix, the basis functions of an iso element are always continuous across the boundaries between edges on the fine mesh (that are not edges in the coarse mesh), so there’s currently no way to make this element. I’ve opened an issue for this: Add discontinuous iso element · Issue #794 · FEniCS/basix · GitHub

11

Yes.
I’ll let @mscroggs comment on the other part of the question.

1 Like
12

Okay I see, thank you

13

Hi, @mscroggs, I am also trying to use P1-iso-P2 element. But I have difficulties plotting results like P2 element. I get the following error

     97 if V.ufl_element().family_name not in [
     98     "Discontinuous Lagrange",
     99     "Lagrange",
   (...)
    103     "P",
    104 ]:
--> 105     raise RuntimeError(
    106         "Can only create meshes from continuous or discontinuous Lagrange spaces"
    107     )
    109 degree = V.ufl_element().degree
    110 if degree == 0:

RuntimeError: Can only create meshes from continuous or discontinuous Lagrange spaces"

The plotting code is modified from a Taylor-Hood element example.

P2 = element(basix.ElementFamily.iso, msh.topology.cell_name(), 1, shape=(2,))
P1 = element("Lagrange", msh.topology.cell_name(), 1)
mixed = mixed_element([P2, P1])

V, Q = functionspace(msh, P2), functionspace(msh, P1)

pyvista.start_xvfb()
topology, cell_types, geometry = plot.vtk_mesh(V)
values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
values[:, :len(outcome)] = outcome.x.array.real[0:geometry.shape[0]*2].reshape((geometry.shape[0], len(outcome)))

# Create a point cloud of glyphs
function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
function_grid["u"] = values
glyphs = function_grid.glyph(orient="u", factor=0.5)

# Create a pyvista-grid for the mesh
msh.topology.create_connectivity(msh.topology.dim, msh.topology.dim)

# Create plotter
plotter = pyvista.Plotter()
plotter.add_mesh(function_grid, style="wireframe", color="k")
plotter.add_mesh(glyphs)
plotter.view_xy()
plotter.show("glyphs.png")
14

You would have to interpolate the P1-iso-P2 into an appropriate continuous or discontinuous Lagrange space to visualize it with Pyvista.

15

Thanks @dokken, after interpolate we can view the output now. But I find the results from p1-iso-p2/p1 is not matching the results as p2/p1.

The first problem I had is that create_mesh doesn’t seem to work with iso element?

domain = basix.ufl.element(basix.ElementFamily.iso, "triangle", 1, shape=(2,))

x = [[0., 0.], [1., 0.], [0., 1.], [1., 1.]]
cells = cells = [[0, 3, 1], [0, 2, 3]]
msh = create_mesh(MPI.COMM_WORLD, cells, x, domain)

This will lead to error

munmap_chunk(): invalid pointer
[ca1ee4a86145:30195] *** Process received signal ***
[ca1ee4a86145:30195] Signal: Aborted (6)
[ca1ee4a86145:30195] Signal code:  (-6)
[ca1ee4a86145:30195] [ 0] /usr/lib/x86_64-linux-gnu/libc.so.6(+0x42520)[0x7b457e13a520]
[ca1ee4a86145:30195] [ 1] /usr/lib/x86_64-linux-gnu/libc.so.6(pthread_kill+0x12c)[0x7b457e18e9fc]
[ca1ee4a86145:30195] [ 2] /usr/lib/x86_64-linux-gnu/libc.so.6(raise+0x16)[0x7b457e13a476]
[ca1ee4a86145:30195] [ 3] /usr/lib/x86_64-linux-gnu/libc.so.6(abort+0xd3)[0x7b457e1207f3]
[ca1ee4a86145:30195] [ 4] /usr/lib/x86_64-linux-gnu/libc.so.6(+0x89676)[0x7b457e181676]
[ca1ee4a86145:30195] [ 5] /usr/lib/x86_64-linux-gnu/libc.so.6(+0xa0cfc)[0x7b457e198cfc]
[ca1ee4a86145:30195] [ 6] /usr/lib/x86_64-linux-gnu/libc.so.6(+0xa0fdc)[0x7b457e198fdc]
[ca1ee4a86145:30195] [ 7] /usr/lib/x86_64-linux-gnu/libc.so.6(free+0xba)[0x7b457e19d49a]
[ca1ee4a86145:30195] [ 8] /usr/lib/x86_64-linux-gnu/libmpi.so.40(ompi_comm_free+0x273)[0x7b4570058ea3]
[ca1ee4a86145:30195] [ 9] /usr/lib/x86_64-linux-gnu/libmpi.so.40(PMPI_Comm_free+0x76)[0x7b4570097246]
[ca1ee4a86145:30195] [10] /usr/lib/x86_64-linux-gnu/libdolfinx_real.so.0.9(_ZN7dolfinx3MPI4CommD2Ev+0x19)[0x7b455d71fca9]
[ca1ee4a86145:30195] [11] /usr/lib/petsc/lib/python3/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x215a47)[0x7b455da67a47]
[ca1ee4a86145:30195] [12] /usr/lib/petsc/lib/python3/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x25e4b2)[0x7b455dab04b2]
[ca1ee4a86145:30195] [13] /usr/lib/petsc/lib/python3/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x29f8f9)[0x7b455daf18f9]
[ca1ee4a86145:30195] [14] python3(+0x169b87)[0x5d02e7575b87]
[ca1ee4a86145:30195] [15] python3(+0x1a2751)[0x5d02e75ae751]
[ca1ee4a86145:30195] [16] python3(+0x169bc7)[0x5d02e7575bc7]
[ca1ee4a86145:30195] [17] python3(+0x1a2751)[0x5d02e75ae751]
[ca1ee4a86145:30195] [18] python3(+0x15b797)[0x5d02e7567797]
[ca1ee4a86145:30195] [19] python3(_PyEval_EvalFrameDefault+0x1be5)[0x5d02e7582eb5]
[ca1ee4a86145:30195] [20] python3(+0x25ae56)[0x5d02e7666e56]
[ca1ee4a86145:30195] [21] python3(PyEval_EvalCode+0x86)[0x5d02e7666d26]
[ca1ee4a86145:30195] [22] python3(+0x281ae8)[0x5d02e768dae8]
[ca1ee4a86145:30195] [23] python3(+0x27c2ef)[0x5d02e76882ef]
[ca1ee4a86145:30195] [24] python3(+0x281885)[0x5d02e768d885]
[ca1ee4a86145:30195] [25] python3(_PyRun_SimpleFileObject+0x1a8)[0x5d02e768ce68]
[ca1ee4a86145:30195] [26] python3(_PyRun_AnyFileObject+0x47)[0x5d02e768cb47]
[ca1ee4a86145:30195] [27] python3(Py_RunMain+0x2be)[0x5d02e768102e]
[ca1ee4a86145:30195] [28] python3(Py_BytesMain+0x2d)[0x5d02e765ad6d]
[ca1ee4a86145:30195] [29] /usr/lib/x86_64-linux-gnu/libc.so.6(+0x29d90)[0x7b457e121d90]
[ca1ee4a86145:30195] *** End of error message ***
Aborted (core dumped)

Then I also can’t get good results with p1-iso-p2 element, while change the element to P2 will work out just fine. I attach my complete code below for your review.

import time
import pyvista
import numpy as np

import ufl
import basix
import dolfinx
import dolfinx.fem.petsc
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block

from dolfinx import cpp as _cpp
from dolfinx import fem
from dolfinx import plot
from dolfinx.fem import (Constant, Function, functionspace, dirichletbc,
                         extract_function_spaces, form, Expression,
                         locate_dofs_geometrical, locate_dofs_topological,
                         create_interpolation_data, assemble_scalar)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_rectangle, locate_entities,
                          locate_entities_boundary, meshtags, DiagonalType, refine,
                          create_mesh)
from dolfinx import geometry

from ufl import div, dx, grad, inner, Measure
from basix.ufl import element, mixed_element

from mpi4py import MPI
from petsc4py import PETSc

domain = basix.ufl.element("Lagrange", "triangle", 2, shape=(2,))
# domain = basix.ufl.element(basix.ElementFamily.iso, "triangle", 1, shape=(2,))

# x = [[0., 0.], [1., 0.], [0., 1.], [1., 1.]]
# cells = cells = [[0, 3, 1], [0, 2, 3]]
# msh = create_mesh(MPI.COMM_WORLD, cells, x, domain)
# for i in range(2):
#     msh.topology.create_entities(1)
#     msh, _, _ = refine(msh)

msh = create_rectangle(MPI.COMM_WORLD,
                        [np.array([0, 0]), np.array([1., 1.])],
                        [8, 8],
                        CellType.triangle, ghost_mode=GhostMode.none,
                        diagonal=DiagonalType.left_right)

# P2 = element(basix.ElementFamily.iso, msh.topology.cell_name(), 1, shape=(2,))
P2 = element("Lagrange", msh.topology.cell_name(), 2, shape=(2,))
P1 = element("Lagrange", msh.topology.cell_name(), 1)
V, Q = functionspace(msh, P2), functionspace(msh, P1)
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)

f = Constant(msh, (PETSc.ScalarType(0.0), PETSc.ScalarType(0.0)))

a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx], [inner(div(u), q) * dx, Constant(msh, PETSc.ScalarType(0))*p*q*dx]])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx])  # type: ignore

#-----------------SET THE BOUNDARY CONDITIONS FOR THE PROBLEM----------------------------------------------
def inlet_profile_expression(x):
    return np.stack((x[0]*0, x[0]**2-x[0]))

# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0)

# Function to mark the lid (y = 1)
def inlet(x):
    return np.isclose(x[1], 1.0)
    
# +
# No-slip condition on boundaries where x = 0, x = 1, and y = 0
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType)  # type: ignore
facets = locate_entities_boundary(msh, 1, noslip_boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, 1, facets), V)

# Driving (lid) velocity condition on top boundary (y = 1)
inlet_profile = Function(V)
inlet_profile.interpolate(inlet_profile_expression)
facets = locate_entities_boundary(msh, 1, inlet)
bc1 = dirichletbc(inlet_profile, locate_dofs_topological(V, 1, facets))

bcs = [bc0, bc1]

# Assemble LHS matrix and RHS vector
a,L = form(a),form(L)

A = assemble_matrix_block(a, bcs=bcs)
A.assemble()
b = assemble_vector_block(L, a, bcs=bcs)

# Create a solver
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("preonly")

# Set the solver type to MUMPS (LU solver) and configure MUMPS to
# handle pressure nullspace
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")
pc.setFactorSetUpSolverType()
pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)  # For pressure nullspace
pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)  # For pressure nullspace

# Create a block vector (x) to store the full solution, and solve
x = A.createVecLeft()
ksp.solve(b, x)

# Create Functions and scatter x solution
u, p = Function(V), Function(Q)
offset = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
u.x.array[:offset] = x.array_r[:offset]
p.x.array[: (len(x.array_r) - offset)] = x.array_r[offset:]

# P2iso = element(basix.ElementFamily.iso, msh.topology.cell_name(), 1, shape=(2,))
# P2lag = element("Lagrange", "triangle", 1, shape=(2,))
# Vlag = functionspace(msh, P2lag)
# Viso = functionspace(msh, P2iso)
# cells_lag_map = Vlag.mesh.topology.index_map(Vlag.mesh.topology.dim)
# num_cells_lag = cells_lag_map.size_local + cells_lag_map.num_ghosts
# cells_lag = np.arange(num_cells_lag, dtype=np.int32)

# ulag = Function(Vlag)
# interpolation_data = create_interpolation_data(Vlag, Viso, cells_lag)
# ulag.interpolate_nonmatching(u, cells_lag, interpolation_data=interpolation_data)
# ulag.x.scatter_forward()

# pyvista.start_xvfb()
# topology, cell_types, geometry = plot.vtk_mesh(Vlag)
# values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
# values[:, :len(u)] = u.x.array.real[0:geometry.shape[0]*2].reshape((geometry.shape[0], len(u)))

pyvista.start_xvfb()
topology, cell_types, geometry = plot.vtk_mesh(V)
values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
values[:, :len(u)] = u.x.array.real[0:geometry.shape[0]*2].reshape((geometry.shape[0], len(u)))

# Create a point cloud of glyphs
function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
function_grid["u"] = values
glyphs = function_grid.glyph(orient="u", factor=0.5)

# Create a pyvista-grid for the mesh
msh.topology.create_connectivity(msh.topology.dim, msh.topology.dim)

# Create plotter
plotter = pyvista.Plotter(off_screen=True)
plotter.add_mesh(function_grid, style="wireframe", color="k")
plotter.add_mesh(glyphs)
plotter.view_xy()
fig_as_array = plotter.screenshot("glyphs.png")

Results from iso

glyphs_iso
glyphs_iso1024×768 20.7 KB

Results from P2

glyphs
glyphs1024×768 128 KB

16

The current design of DOLFINx doesn’t allow for other elements than Lagrange as the coordinate element.

I guess to visualize a P1-iso-P2 element, one should really refine the grid and interpolate onto an appropriate space, for instance (DG, 2):

from mpi4py import MPI
import ufl
import basix.ufl
import dolfinx
import numpy as np

msh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 1, 1)



iso = basix.ufl.element(basix.ElementFamily.iso, msh.topology.cell_name(), 1)
V = dolfinx.fem.functionspace(msh, iso)
u = dolfinx.fem.Function(V)


visualizing_mesh, _,_ = dolfinx.mesh.refine(msh)
Q = dolfinx.fem.functionspace(visualizing_mesh, ("DG", 2))
q = dolfinx.fem.Function(Q)
num_cells_out = visualizing_mesh.topology.index_map(2).size_local
ipd = dolfinx.fem.create_interpolation_data(Q, V, np.arange(num_cells_out, dtype=np.int32), padding=1e-10)


with dolfinx.io.XDMFFile(msh.comm,"output.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh)
with dolfinx.io.VTXWriter(visualizing_mesh.comm, "output.bp", [q]) as bp:
    for i in range(len(u.x.array)):
        u.x.array[:] = 0.0
        u.x.array[i] = 1.0

        q.x.array[:] = 0.0
        q.interpolate_nonmatching(u, np.arange(num_cells_out, dtype=np.int32), interpolation_data=ipd)
        q.x.scatter_forward()


        bp.write(i)

Maybe @mscroggs has a clever idea for visualizing them.

17

You can interpolation in to (DG, 1) on the refined mesh as P1-iso-P2 is linear on each sub-triangle. Other than that, your suggestion is the only way I can think of of visualising this