Mesh radius for higher order meshes

mesh
1

I want to weakly impose dirichlet boundary condition but it seems that ufl.Circumradius is not supported for higher order meshes. I get the following error (dolfinx 0.9.0, conda build):

ValueError: Circumradius only makes sense for affine simplex cells

Is there any workaround in this situation? The tutorial model is hopefully a good MWE:

from dolfinx import fem, mesh, plot, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
import numpy
from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
from ufl import (Circumradius, FacetNormal, SpatialCoordinate, TrialFunction, TestFunction,
                 div, dx, ds, grad, inner)

N = 8

# domain = mesh.create_unit_square(MPI.COMM_WORLD, N, N)
# create a second order mesh with gmsh
gdim = 2
gmsh.initialize()
gmsh.clear()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
box = gmsh.model.occ.add_rectangle(0, 0, 0, 1, 1)
gmsh.model.occ.synchronize()
    
all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

all_edges = gmsh.model.getEntities(gdim-1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  

gmsh.option.setNumber("Mesh.ElementOrder", 2)
gmsh.model.mesh.generate()

domain, ct, ft = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim)


V = fem.functionspace(domain, ("Lagrange", 1))

uD = fem.Function(V)
x = SpatialCoordinate(domain)
u_ex =  1 + x[0]**2 + 2 * x[1]**2
uD.interpolate(fem.Expression(u_ex, V.element.interpolation_points()))
f = -div(grad(u_ex))

u = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(domain)
h = 2 * Circumradius(domain)
alpha = fem.Constant(domain, default_scalar_type(10))
a = inner(grad(u), grad(v)) * dx - inner(n, grad(u)) * v * ds
a += - inner(n, grad(v)) * u * ds + alpha / h * inner(u, v) * ds
L = inner(f, v) * dx 
L += - inner(n, grad(v)) * uD * ds + alpha / h * inner(uD, v) * ds

problem = LinearProblem(a, L)
uh = problem.solve()

Thanks!

2

You can get area/volume of each cell with a DG-0 function, then compute its circumradius as r = ufl.sqrt(area). See Volume of an element - #4 by dokken

2 Likes