Adaptive mesh on curved domain with refine command

General
,
1

Dear community.

I’m trying to replicate the adaptive mesh adaptive refinement from Adaptive mesh refinement with NetGen and DOLFINx but using the command refine instead of netgen. The code is below.

I generate the mesh with

SetFactory("OpenCASCADE");
//+
N=2; // MESH RESOLUTION

h=1/N; // MESH SIZE BASED ON N
Mesh.CharacteristicLengthMin=h;
Mesh.CharacteristicLengthMax=h;//+//+
//+
Circle(1) = {0, 0, 0, 1, 0, 3*Pi/2.};
//+
Point(3) = {0, 0, 0, 1.0};
//+
Line(2) = {1, 3};
//+
Line(3) = {3, 2};
//+
Curve Loop(1) = {2, 3, -1};
//+
Plane Surface(1) = {1};
//+
Physical Curve("1", 1) = {2, 1, 3};
//+
Physical Surface("1", 1) = {1};

Mesh 2;
Save "circle_mesh.msh"

then run

import dolfinx
import dolfinx.fem.petsc
import numpy as np
import pyvista
import ufl
from dolfinx.cpp.mesh import create_cell_partitioner
from dolfinx.mesh import (
    RefinementOption,
    refine,
    GhostMode,
)
from mpi4py import MPI
from packaging.version import Version
from petsc4py import PETSc
from slepc4py import SLEPc


def write_frame(plotter: pyvista.Plotter, uh_r: dolfinx.fem.Function):
    # Scale uh_r to be consistent between refinement steps, as it can be multiplied by -1
    uh_r_min = msh.comm.allreduce(uh_r.x.array.min(), op=MPI.MIN)
    uh_r_max = msh.comm.allreduce(uh_r.x.array.max(), op=MPI.MAX)
    uh_sign = np.sign(uh_r_min)
    if np.isclose(uh_sign, 0):
        uh_sign = np.sign(uh_r_max)
    assert not np.isclose(uh_sign, 0), "uh_r has zero values, cannot determine sign."
    uh_r.x.array[:] *= uh_sign

    # Update plot with refined mesh
    grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(msh))
    curved_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(uh_r.function_space))
    curved_grid.point_data["u"] = uh_r.x.array
    curved_grid = curved_grid.tessellate()
    curved_actor = plotter.add_mesh(
        curved_grid,
        show_edges=False,
    )

    actor = plotter.add_mesh(grid, style="wireframe", color="black")
    plotter.view_xy()
    plotter.write_frame()
    plotter.remove_actor(actor)
    plotter.remove_actor(curved_actor)

def solve(
    mesh: dolfinx.mesh.Mesh,
    facet_tags: dolfinx.mesh.MeshTags,
) -> tuple[float, dolfinx.fem.Function, dolfinx.fem.Function]:
    # We define the lhs and rhs bilinear forms
    V = dolfinx.fem.functionspace(mesh, ("Lagrange", 3))
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
    m = ufl.inner(u, v) * ufl.dx

    tdim = msh.topology.dim
    fdim = tdim - 1

    msh.topology.create_connectivity(fdim, tdim)
    boundary_facets =dolfinx.mesh.exterior_facet_indices(msh.topology)

    mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
    boundary_dofs = dolfinx.fem.locate_dofs_topological(
        V, mesh.topology.dim - 1, boundary_facets
    )

    # We create a zero boundary condition for these dofs to be in the suitable space, and
    # set up the discrete matrices `A` and `M`
    bc = dolfinx.fem.dirichletbc(0.0, boundary_dofs, V)
    A = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a), bcs=[bc])
    A.assemble()
    if Version(dolfinx.__version__) < Version("0.10.0"):
        diag_kwargs = {"diagonal": 0.0}
    else:
        diag_kwargs = {"diag": 0.0}

    M = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(m), bcs=[bc], **diag_kwargs)
    M.assemble()

    # Next, we define the SLEPc Eigenvalue Problem Solver (EPS), and set up to use a shift
    # and invert (SINVERT) spectral transformation where the preconditioner factorisation
    # is computed using [MUMPS](https://mumps-solver.org/index.php).

    E = SLEPc.EPS().create(mesh.comm)
    E.setType(SLEPc.EPS.Type.ARNOLDI)
    E.setProblemType(SLEPc.EPS.ProblemType.GHEP)
    E.setDimensions(1, SLEPc.DECIDE)
    E.setOperators(A, M)
    ST = E.getST()
    ST.setType(SLEPc.ST.Type.SINVERT)
    PC = ST.getKSP().getPC()
    PC.setType("lu")
    PC.setFactorSolverType("mumps")
    E.setST(ST)
    E.solve()
    assert E.getConvergedReason() >= 0, "Eigenvalue solver did not converge"

    # We get the real and imaginary parts of the first eigenvector along with the eigenvalue.
    uh_r = dolfinx.fem.Function(V)
    uh_i = dolfinx.fem.Function(V)
    lam = E.getEigenpair(0, uh_r.x.petsc_vec, uh_i.x.petsc_vec)
    E.destroy()
    uh_r.x.scatter_forward()
    uh_i.x.scatter_forward()
    return (lam, uh_r, uh_i)


def mark_cells(uh_r: dolfinx.fem.Function, lam: float):
    mesh = uh_r.function_space.mesh
    W = dolfinx.fem.functionspace(mesh, ("DG", 0))
    w = ufl.TestFunction(W)
    eta_squared = dolfinx.fem.Function(W)
    f = dolfinx.fem.Constant(mesh, 1.0)
    h = dolfinx.fem.Function(W)
    h.x.array[:] = mesh.h(mesh.topology.dim, np.arange(len(h.x.array), dtype=np.int32))
    n = ufl.FacetNormal(mesh)

    G = (  # compute cellwise error estimator
        ufl.inner(h**2 * (f + ufl.div(ufl.grad(uh_r))) ** 2, w) * ufl.dx
        + ufl.inner(h("+") / 2 * ufl.jump(ufl.grad(uh_r), n) ** 2, w("+")) * ufl.dS
        + ufl.inner(h("-") / 2 * ufl.jump(ufl.grad(uh_r), n) ** 2, w("-")) * ufl.dS
    )
    dolfinx.fem.petsc.assemble_vector(eta_squared.x.petsc_vec, dolfinx.fem.form(G))
    eta = dolfinx.fem.Function(W)
    eta.x.array[:] = np.sqrt(eta_squared.x.array[:])

    eta_max = eta.x.petsc_vec.max()[1]

    theta = 0.5
    should_refine = ufl.conditional(ufl.gt(eta, theta * eta_max), 1, 0)
    markers = dolfinx.fem.Function(W)
    ip = W.element.interpolation_points
    if Version(dolfinx.__version__) < Version("0.10.0"):
        ip = ip()
    markers.interpolate(dolfinx.fem.Expression(should_refine, ip))
    return np.flatnonzero(np.isclose(markers.x.array.astype(np.int32), 1))


max_iterations = 15
exact = 3.375610652693620492628**2
termination_criteria = 1e-5
msh, cell_tags, facet_tags = dolfinx.io.gmshio.read_from_msh(
        "circle_mesh.msh", MPI.COMM_SELF, gdim=2)

plotter = pyvista.Plotter()
plotter.open_gif("amr.gif", fps=1)
for i in range(max_iterations):
    lam, uh_r, _ = solve(msh, facet_tags)

    relative_error = (lam - exact) / abs(exact)
    PETSc.Sys.Print(
        f"Iteration {i + 1}/{max_iterations}, {lam=:.5e}, {exact=:.5e}, {relative_error=:.2e}"
    )

    cells_to_mark = mark_cells(uh_r, lam)
    local_edges = dolfinx.mesh.compute_incident_entities(
        msh.topology, cells_to_mark, 2, 1
    )

    msh, parent_cell, parent_facet = refine(
        msh,
        local_edges,
        create_cell_partitioner(GhostMode.shared_facet),
        RefinementOption.parent_cell_and_facet,
    )

    write_frame(plotter, uh_r)
    if relative_error < termination_criteria:
        PETSc.Sys.Print(f"Converged in {i + 1} iterations.")
        break

plotter.close()

However, the mesh is not refined at the curved boundaries. What I’m missing?. Am I forced to use the Netgen refinement for this type of domain?. Moreover, if I try to use Mesh.ElementOrder = 2;, I get a SEGV after the first adaptive iteration.

Thanks in advance.

amr

2

The built in refinement routine, that uses the plaza algorithm, does not at the moment support higher order geometries.

Additionally, the built in routine does not carry the same information as the netgen meshing tool. The netgen meshing object carries information about the CSG boundary, meaning that when refining this is taken into account.

So for now you are «forced» to use netgen to achieve this behavior.
Long-term, I’ll hopefully have a gmsh interface for this. However it is not the first thing on my list of priorities.

1 Like
3

Thank you @dokken!. I marked the post above as the solution of the principal issue regarding the curved domain. On the other hand, is the SEGV error coming from adapting a degree>1 mesh related to the same problem?

4

Yes, I don’t think we have support for higher order mesh refinement with plaza yet. Should of course rather throw an error than a segfault.

1 Like