From e8c2a621898432a517181353d7148199d0645189 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Wed, 7 Feb 2024 17:27:04 +0000 Subject: [PATCH] Mesh Hierarchy support when working with Netgen mesh 2D (#3314) * Mesh Hierarchy support when working with Netgen mesh --- firedrake/mg/mesh.py | 57 +++-- tests/conftest.py | 16 ++ tests/multigrid/test_netgen_gmg.py | 101 ++++++++ tests/regression/test_netgen.py | 383 ++++++++++++++++++----------- 4 files changed, 394 insertions(+), 163 deletions(-) create mode 100644 tests/multigrid/test_netgen_gmg.py diff --git a/firedrake/mg/mesh.py b/firedrake/mg/mesh.py index e4eb1f67b8..a40fad62f7 100644 --- a/firedrake/mg/mesh.py +++ b/firedrake/mg/mesh.py @@ -9,7 +9,6 @@ from firedrake.cython import mgimpl as impl from .utils import set_level - __all__ = ("HierarchyBase", "MeshHierarchy", "ExtrudedMeshHierarchy", "NonNestedHierarchy", "SemiCoarsenedExtrudedHierarchy") @@ -80,28 +79,58 @@ def __getitem__(self, idx): def MeshHierarchy(mesh, refinement_levels, refinements_per_level=1, + netgen_flags=False, reorder=None, distribution_parameters=None, callbacks=None, mesh_builder=firedrake.Mesh): """Build a hierarchy of meshes by uniformly refining a coarse mesh. - :arg mesh: the coarse :func:`~.Mesh` to refine - :arg refinement_levels: the number of levels of refinement - :arg refinements_per_level: the number of refinements for each - level in the hierarchy. - :arg distribution_parameters: options controlling mesh - distribution, see :func:`~.Mesh` for details. If ``None``, - use the same distribution parameters as were used to - distribute the coarse mesh, otherwise, these options override - the default. - :arg reorder: optional flag indicating whether to reorder the - refined meshes. - :arg callbacks: A 2-tuple of callbacks to call before and + Parameters + ---------- + mesh : MeshGeometry + the coarse mesh to refine + refinement_levels : int + the number of levels of refinement + refinements_per_level : int + the number of refinements for each level in the hierarchy. + netgen_flags : bool, dict + either a bool or a dictionary containing options for Netgen. + If not False the hierachy is constructed using ngsPETSc, if + None hierarchy constructed in a standard manner. + distribution_parameters : dict + options controlling mesh distribution, see :py:func:`.Mesh` + for details. If ``None``, use the same distribution + parameters as were used to distribute the coarse mesh, + otherwise, these options override the default. + reorder : bool + optional flag indicating whether to reorder the + refined meshes. + callbacks : tuple + A 2-tuple of callbacks to call before and after refinement of the DM. The before callback receives the DM to be refined (and the current level), the after callback receives the refined DM (and the current level). - :arg mesh_builder: Function to turn a DM into a ``Mesh``. Used by pyadjoint. + mesh_builder + Function to turn a DM into a ``Mesh``. Used by pyadjoint. + Returns + ------- + A :py:class:`HierarchyBase` object representing the + mesh hierarchy. """ + + if (isinstance(netgen_flags, bool) and netgen_flags) or isinstance(netgen_flags, dict): + try: + from ngsPETSc import NetgenHierarchy + except ImportError: + raise ImportError("Unable to import netgen and ngsPETSc. Please ensure that netgen and ngsPETSc\ + are installed and available to Firedrake. You can do this via \ + firedrake-update --netgen.") + if hasattr(mesh, "netgen_mesh"): + return NetgenHierarchy(mesh, refinement_levels, flags=netgen_flags) + else: + raise RuntimeError("Cannot create a NetgenHierarchy from a mesh that has not been generated by\ + Netgen.") + cdm = mesh.topology_dm cdm.setRefinementUniform(True) dms = [] diff --git a/tests/conftest.py b/tests/conftest.py index 5e0418bd31..57866d20c4 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -20,6 +20,9 @@ def pytest_configure(config): config.addinivalue_line( "markers", "skipplot: mark as skipped if matplotlib is not installed") + config.addinivalue_line( + "markers", + "skipnetgen: mark as skipped if netgen and ngsPETSc is not installed") def pytest_collection_modifyitems(session, config, items): @@ -39,6 +42,15 @@ def pytest_collection_modifyitems(session, config, items): except ImportError: ml_backend = False + try: + import netgen + del netgen + import ngsPETSc + del ngsPETSc + netgen_installed = True + except ImportError: + netgen_installed = False + for item in items: if complex_mode: if item.get_closest_marker("skipcomplex") is not None: @@ -57,6 +69,10 @@ def pytest_collection_modifyitems(session, config, items): if item.get_closest_marker("skipplot") is not None: item.add_marker(pytest.mark.skip(reason="Test cannot be run unless Matplotlib is installed")) + if not netgen_installed: + if item.get_closest_marker("skipnetgen") is not None: + item.add_marker(pytest.mark.skip(reason="Test cannot be run unless Netgen and ngsPETSc are installed")) + @pytest.fixture(scope="module", autouse=True) def check_empty_tape(request): diff --git a/tests/multigrid/test_netgen_gmg.py b/tests/multigrid/test_netgen_gmg.py new file mode 100644 index 0000000000..375d28f202 --- /dev/null +++ b/tests/multigrid/test_netgen_gmg.py @@ -0,0 +1,101 @@ +from firedrake import * +import pytest + + +def create_netgen_mesh_circle(): + from netgen.geom2d import Circle, CSG2d + geo = CSG2d() + + circle = Circle(center=(0, 0), radius=1.0, mat="mat1", bc="circle") + geo.Add(circle) + + ngmesh = geo.GenerateMesh(maxh=0.75) + return ngmesh + + +@pytest.mark.skipcomplex +@pytest.mark.skipnetgen +def test_netgen_mg_circle(): + ngmesh = create_netgen_mesh_circle() + mesh = Mesh(ngmesh) + nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": 3}) + mesh = nh[-1] + + V = FunctionSpace(mesh, "CG", 3) + + u = TrialFunction(V) + v = TestFunction(V) + + a = inner(grad(u), grad(v))*dx + labels = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(codim=1)) if name in ["circle"]] + bcs = DirichletBC(V, zero(), labels) + x, y = SpatialCoordinate(mesh) + + f = 4+0*x + L = f*v*dx + exact = (1-x**2-y**2) + + u = Function(V) + solve(a == L, u, bcs=bcs, solver_parameters={"ksp_type": "cg", + "pc_type": "mg"}) + expect = Function(V).interpolate(exact) + assert (norm(assemble(u - expect)) <= 1e-6) + + +@pytest.mark.skipcomplex +@pytest.mark.skipnetgeun +def test_netgen_mg_circle_non_uniform_degree(): + ngmesh = create_netgen_mesh_circle() + mesh = Mesh(ngmesh) + nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": [1, 2, 3]}) + mesh = nh[-1] + + V = FunctionSpace(mesh, "CG", 3) + + u = TrialFunction(V) + v = TestFunction(V) + + a = inner(grad(u), grad(v))*dx + labels = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(codim=1)) if name in ["circle"]] + bcs = DirichletBC(V, zero(), labels) + x, y = SpatialCoordinate(mesh) + + f = 4+0*x + L = f*v*dx + exact = (1-x**2-y**2) + + u = Function(V) + solve(a == L, u, bcs=bcs, solver_parameters={"ksp_type": "cg", + "pc_type": "mg"}) + expect = Function(V).interpolate(exact) + assert (norm(assemble(u - expect)) <= 1e-6) + + +@pytest.mark.skipcomplex +@pytest.mark.parallel +@pytest.mark.skipnetgen +def test_netgen_mg_circle_parallel(): + ngmesh = create_netgen_mesh_circle() + mesh = Mesh(ngmesh) + nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": 3}) + mesh = nh[-1] + + V = FunctionSpace(mesh, "CG", 3) + + u = TrialFunction(V) + v = TestFunction(V) + + a = inner(grad(u), grad(v))*dx + labels = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(codim=1)) if name in ["circle"]] + bcs = DirichletBC(V, zero(), labels) + x, y = SpatialCoordinate(mesh) + + f = 4+0*x + L = f*v*dx + exact = (1-x**2-y**2) + + u = Function(V) + solve(a == L, u, bcs=bcs, solver_parameters={"ksp_type": "cg", + "pc_type": "mg"}) + expect = Function(V).interpolate(exact) + assert norm(assemble(u - expect)) <= 1e-6 diff --git a/tests/regression/test_netgen.py b/tests/regression/test_netgen.py index 7f17ffc2a8..7560a77db6 100644 --- a/tests/regression/test_netgen.py +++ b/tests/regression/test_netgen.py @@ -1,25 +1,26 @@ from firedrake import * +from firedrake.__future__ import interpolate import numpy as np -import gc from petsc4py import PETSc import pytest printf = lambda msg: PETSc.Sys.Print(msg) -def poisson(h, degree=2): - try: - from netgen.geom2d import SplineGeometry - import netgen - except ImportError: - pytest.skip(reason="Netgen unavailable, skipping Netgen test.") +def square_geometry(h): + from netgen.geom2d import SplineGeometry + geo = SplineGeometry() + geo.AddRectangle((0, 0), (np.pi, np.pi), bc="rect") + ngmesh = geo.GenerateMesh(maxh=h) + return ngmesh + +def poisson(h, degree=2): + import netgen comm = COMM_WORLD # Setting up Netgen geometry and mesh if comm.Get_rank() == 0: - geo = SplineGeometry() - geo.AddRectangle((0, 0), (np.pi, np.pi), bc="rect") - ngmesh = geo.GenerateMesh(maxh=h) + ngmesh = square_geometry(h) labels = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(codim=1)) if name == "rect"] else: ngmesh = netgen.libngpy._meshing.Mesh(2) @@ -31,9 +32,8 @@ def poisson(h, degree=2): V = FunctionSpace(msh, "CG", degree) u = TrialFunction(V) v = TestFunction(V) - f = Function(V) x, y = SpatialCoordinate(msh) - f.interpolate(2*sin(x)*sin(y)) + f = assemble(interpolate(2*sin(x)*sin(y), V)) a = inner(grad(u), grad(v))*dx l = inner(f, v) * dx u = Function(V) @@ -53,11 +53,8 @@ def poisson(h, degree=2): def poisson3D(h, degree=2): - try: - from netgen.csg import CSGeometry, OrthoBrick, Pnt - import netgen - except ImportError: - pytest.skip(reason="Netgen unavailable, skipping Netgen test.") + from netgen.csg import CSGeometry, OrthoBrick, Pnt + import netgen comm = COMM_WORLD # Setting up Netgen geometry and mesh @@ -79,9 +76,8 @@ def poisson3D(h, degree=2): V = FunctionSpace(msh, "CG", degree) u = TrialFunction(V) v = TestFunction(V) - f = Function(V) x, y, z = SpatialCoordinate(msh) - f.interpolate(3*sin(x)*sin(y)*sin(z)) + f = assemble(interpolate(3*sin(x)*sin(y)*sin(z), V)) a = inner(grad(u), grad(v))*dx l = inner(f, v) * dx u = Function(V) @@ -101,6 +97,7 @@ def poisson3D(h, degree=2): return S +@pytest.mark.skipnetgen def test_firedrake_Poisson_netgen(): diff = np.array([poisson(h)[0] for h in [1/2, 1/4, 1/8]]) print("l2 error norms:", diff) @@ -109,6 +106,17 @@ def test_firedrake_Poisson_netgen(): assert (np.array(conv) > 2.8).all() +@pytest.mark.skipnetgen +@pytest.mark.parallel +def test_firedrake_Poisson_netgen_parallel(): + diff = np.array([poisson(h)[0] for h in [1/2, 1/4, 1/8]]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > 2.8).all() + + +@pytest.mark.skipnetgen def test_firedrake_Poisson3D_netgen(): diff = np.array([poisson3D(h) for h in [1, 1/2, 1/4]]) print("l2 error norms:", diff) @@ -117,12 +125,10 @@ def test_firedrake_Poisson3D_netgen(): assert (np.array(conv) > 2.8).all() +@pytest.mark.skipnetgen def test_firedrake_integral_2D_netgen(): - try: - from netgen.geom2d import SplineGeometry - import netgen - except ImportError: - pytest.skip(reason="Netgen unavailable, skipping Netgen test.") + from netgen.geom2d import SplineGeometry + import netgen comm = COMM_WORLD if comm.Get_rank() == 0: @@ -137,18 +143,15 @@ def test_firedrake_integral_2D_netgen(): msh = Mesh(ngmesh) V = FunctionSpace(msh, "CG", 3) x, y = SpatialCoordinate(msh) - f = Function(V).interpolate(x*x+y*y*y+x*y) + f = assemble(interpolate(x*x+y*y*y+x*y, V)) assert abs(assemble(f * dx) - (5/6)) < 1.e-10 +@pytest.mark.skipnetgen def test_firedrake_integral_3D_netgen(): - try: - from netgen.csg import CSGeometry, OrthoBrick, Pnt - import netgen - except ImportError: - pytest.skip(reason="Netgen unavailable, skipping Netgen test.") + from netgen.csg import CSGeometry, OrthoBrick, Pnt + import netgen - # Setting up Netgen geometry and mes comm = COMM_WORLD if comm.Get_rank() == 0: box = OrthoBrick(Pnt(0, 0, 0), Pnt(1, 1, 1)) @@ -165,20 +168,17 @@ def test_firedrake_integral_3D_netgen(): msh = Mesh(ngmesh) V = FunctionSpace(msh, "CG", 3) x, y, z = SpatialCoordinate(msh) - f = Function(V).interpolate(2 * x + 3 * y * y + 4 * z * z * z) + f = assemble(interpolate(2 * x + 3 * y * y + 4 * z * z * z, V)) assert abs(assemble(f * ds) - (2 + 4 + 2 + 5 + 2 + 6)) < 1.e-10 +@pytest.mark.skipnetgen def test_firedrake_integral_ball_netgen(): - try: - from netgen.csg import CSGeometry, Pnt, Sphere - from netgen.meshing import MeshingParameters - from netgen.meshing import MeshingStep - import netgen - except ImportError: - pytest.skip(reason="Netgen unavailable, skipping Netgen test.") - - # Setting up Netgen geometry and mes + from netgen.csg import CSGeometry, Pnt, Sphere + from netgen.meshing import MeshingParameters + from netgen.meshing import MeshingStep + import netgen + comm = COMM_WORLD if comm.Get_rank() == 0: geo = CSGeometry() @@ -191,18 +191,15 @@ def test_firedrake_integral_ball_netgen(): msh = Mesh(ngmesh) V = FunctionSpace(msh, "CG", 3) x, y, z = SpatialCoordinate(msh) - f = Function(V).interpolate(1+0*x) + f = assemble(interpolate(1+0*x, V)) assert abs(assemble(f * dx) - 4*np.pi) < 1.e-2 +@pytest.mark.skipnetgen def test_firedrake_integral_sphere_high_order_netgen(): - try: - from netgen.csg import CSGeometry, Pnt, Sphere - import netgen - except ImportError: - pytest.skip(reason="Netgen unavailable, skipping Netgen test.") + from netgen.csg import CSGeometry, Pnt, Sphere + import netgen - # Setting up Netgen geometry and mes comm = COMM_WORLD if comm.Get_rank() == 0: geo = CSGeometry() @@ -215,113 +212,201 @@ def test_firedrake_integral_sphere_high_order_netgen(): homsh = Mesh(msh.curve_field(4)) V = FunctionSpace(homsh, "CG", 4) x, y, z = SpatialCoordinate(homsh) - f = Function(V).interpolate(1+0*x) + f = assemble(interpolate(1+0*x, V)) assert abs(assemble(f * dx) - (4/3)*np.pi) < 1.e-4 +@pytest.mark.skipnetgen +@pytest.mark.parallel +def test_firedrake_integral_sphere_high_order_netgen_parallel(): + from netgen.csg import CSGeometry, Pnt, Sphere + import netgen + + comm = COMM_WORLD + if comm.Get_rank() == 0: + geo = CSGeometry() + geo.Add(Sphere(Pnt(0, 0, 0), 1).bc("sphere")) + ngmesh = geo.GenerateMesh(maxh=0.7) + else: + ngmesh = netgen.libngpy._meshing.Mesh(3) + + msh = Mesh(ngmesh) + homsh = Mesh(msh.curve_field(2)) + V = FunctionSpace(homsh, "CG", 2) + x, y, z = SpatialCoordinate(homsh) + f = assemble(interpolate(1+0*x, V)) + assert abs(assemble(f * dx) - (4/3)*np.pi) < 1.e-2 + + +@pytest.mark.skipnetgen @pytest.mark.skipcomplex def test_firedrake_Adaptivity_netgen(): - try: - from netgen.geom2d import SplineGeometry - import netgen - except ImportError: - pytest.skip(reason="Netgen unavailable, skipping Netgen test.") - - try: - from slepc4py import SLEPc - except ImportError: - pytest.skip(reason="SLEPc unavailable, skipping adaptive test refinement.") - - gc.collect() - comm = COMM_WORLD + from netgen.occ import WorkPlane, OCCGeometry, Axes + from netgen.occ import X, Z - def Solve(msh, labels): - V = FunctionSpace(msh, "CG", 2) - u = TrialFunction(V) + def solve_poisson(mesh): + V = FunctionSpace(mesh, "CG", 1) + uh = Function(V, name="Solution") v = TestFunction(V) - a = inner(grad(u), grad(v))*dx - m = (u*v)*dx - uh = Function(V) - bc = DirichletBC(V, 0, labels) - A = assemble(a, bcs=bc) - M = assemble(m) - Asc, Msc = A.M.handle, M.M.handle - E = SLEPc.EPS().create() - E.setType(SLEPc.EPS.Type.ARNOLDI) - E.setProblemType(SLEPc.EPS.ProblemType.GHEP) - E.setDimensions(1, SLEPc.DECIDE) - E.setOperators(Asc, Msc) - ST = E.getST() - ST.setType(SLEPc.ST.Type.SINVERT) - PC = ST.getKSP().getPC() - PC.setType("lu") - PC.setFactorSolverType("mumps") - E.setST(ST) - E.solve() - vr, vi = Asc.getVecs() - with uh.dat.vec_wo as vr: - lam = E.getEigenpair(0, vr, vi) - return (lam, uh, V) - - def Mark(msh, uh, lam): - W = FunctionSpace(msh, "DG", 0) + bc = DirichletBC(V, 0, "on_boundary") + f = Constant(1) + F = inner(grad(uh), grad(v))*dx - inner(f, v)*dx + solve(F == 0, uh, bc) + return uh + + def estimate_error(mesh, uh): + W = FunctionSpace(mesh, "DG", 0) + eta_sq = Function(W) w = TestFunction(W) - R_T = lam.real*uh + div(grad(uh)) - n = FacetNormal(V.mesh()) - h = CellDiameter(msh) - R_dT = dot(grad(uh), n) - eta = assemble(h**2*R_T**2*w*dx + (h("+")+h("-"))*(R_dT("+")-R_dT("-"))**2*(w("+")+w("-"))*dS) - frac = .95 - delfrac = .05 - part = .2 - mark = Function(W) - with mark.dat.vec as markedVec: - with eta.dat.vec as etaVec: - sum_eta = etaVec.sum() - if sum_eta < tolerance: - return markedVec - eta_max = etaVec.max()[1] - sct, etaVec0 = PETSc.Scatter.toZero(etaVec) - markedVec0 = etaVec0.duplicate() - sct(etaVec, etaVec0) - if etaVec.getComm().getRank() == 0: - eta = etaVec0.getArray() - marked = np.zeros(eta.size, dtype='bool') - sum_marked_eta = 0. - while sum_marked_eta < part*sum_eta: - new_marked = (~marked) & (eta > frac*eta_max) - sum_marked_eta += sum(eta[new_marked]) - marked += new_marked - frac -= delfrac - markedVec0.getArray()[:] = 1.0*marked[:] - sct(markedVec0, markedVec, mode=PETSc.Scatter.Mode.REVERSE) - return mark - tolerance = 1e-16 - max_iterations = 15 - exact = 3.375610652693620492628**2 - geo = SplineGeometry() - pnts = [(0, 0), (1, 0), (1, 1), - (0, 1), (-1, 1), (-1, 0), - (-1, -1), (0, -1)] - p1, p2, p3, p4, p5, p6, p7, p8 = [geo.AppendPoint(*pnt) for pnt in pnts] - curves = [[["line", p1, p2], "line"], - [["spline3", p2, p3, p4], "curve"], - [["spline3", p4, p5, p6], "curve"], - [["spline3", p6, p7, p8], "curve"], - [["line", p8, p1], "line"]] - [geo.Append(c, bc=bc) for c, bc in curves] - if comm.Get_rank() == 0: - ngmsh = geo.GenerateMesh(maxh=0.2) - labels = [i+1 for i, name in enumerate(ngmsh.GetRegionNames(codim=1)) if name == "line" or name == "curve"] - else: - ngmsh = netgen.libngpy._meshing.Mesh(2) - labels = None - labels = comm.bcast(labels, root=0) - msh = Mesh(ngmsh) + f = Constant(1) + h = CellDiameter(mesh) + n = FacetNormal(mesh) + v = CellVolume(mesh) + + # Compute error indicator cellwise + G = inner(eta_sq / v, w)*dx + G = G - inner(h**2 * (f + div(grad(uh)))**2, w) * dx + G = G - inner(h('+')/2 * jump(grad(uh), n)**2, w('+')) * dS + + # Each cell is an independent 1x1 solve, so Jacobi is exact + sp = {"mat_type": "matfree", + "ksp_type": "richardson", + "pc_type": "jacobi"} + solve(G == 0, eta_sq, solver_parameters=sp) + eta = Function(W) + eta.interpolate(sqrt(eta_sq)) # the above computed eta^2 + + with eta.dat.vec_ro as eta_: + error_est = sqrt(eta_.dot(eta_)) + return (eta, error_est) + + def adapt(mesh, eta): + W = FunctionSpace(mesh, "DG", 0) + markers = Function(W) + with eta.dat.vec_ro as eta_: + eta_max = eta_.max()[1] + + theta = 0.5 + should_refine = conditional(gt(eta, theta*eta_max), 1, 0) + markers.interpolate(should_refine) + + refined_mesh = mesh.refine_marked_elements(markers) + return refined_mesh + + rect1 = WorkPlane(Axes((0, 0, 0), n=Z, h=X)).Rectangle(1, 2).Face() + rect2 = WorkPlane(Axes((0, 1, 0), n=Z, h=X)).Rectangle(2, 1).Face() + L = rect1 + rect2 + + geo = OCCGeometry(L, dim=2) + ngmsh = geo.GenerateMesh(maxh=0.1) + mesh = Mesh(ngmsh) + + max_iterations = 10 + error_estimators = [] + dofs = [] for i in range(max_iterations): - printf("level {}".format(i)) - lam, uh, V = Solve(msh, labels) - mark = Mark(msh, uh, lam) - msh = msh.refine_marked_elements(mark) - File("Sol.pvd").write(uh) - assert (abs(lam-exact) < 1e-2) + uh = solve_poisson(mesh) + (eta, error_est) = estimate_error(mesh, uh) + error_estimators.append(error_est) + dofs.append(uh.function_space().dim()) + mesh = adapt(mesh, eta) + assert error_estimators[-1] < 0.05 + + +@pytest.mark.skipnetgen +@pytest.mark.skipcomplex +@pytest.mark.parallel +def test_firedrake_Adaptivity_netgen_parallel(): + from netgen.occ import WorkPlane, OCCGeometry, Axes + from netgen.occ import X, Z + + def solve_poisson(mesh): + V = FunctionSpace(mesh, "CG", 1) + uh = Function(V, name="Solution") + v = TestFunction(V) + bc = DirichletBC(V, 0, "on_boundary") + f = Constant(1) + F = inner(grad(uh), grad(v))*dx - inner(f, v)*dx + solve(F == 0, uh, bc) + return uh + + def estimate_error(mesh, uh): + W = FunctionSpace(mesh, "DG", 0) + eta_sq = Function(W) + w = TestFunction(W) + f = Constant(1) + h = CellDiameter(mesh) + n = FacetNormal(mesh) + v = CellVolume(mesh) + + # Compute error indicator cellwise + G = inner(eta_sq / v, w)*dx + G = G - inner(h**2 * (f + div(grad(uh)))**2, w) * dx + G = G - inner(h('+')/2 * jump(grad(uh), n)**2, w('+')) * dS + + # Each cell is an independent 1x1 solve, so Jacobi is exact + sp = {"mat_type": "matfree", + "ksp_type": "richardson", + "pc_type": "jacobi"} + solve(G == 0, eta_sq, solver_parameters=sp) + eta = Function(W) + eta.interpolate(sqrt(eta_sq)) # the above computed eta^2 + + with eta.dat.vec_ro as eta_: + error_est = sqrt(eta_.dot(eta_)) + return (eta, error_est) + + def adapt(mesh, eta): + W = FunctionSpace(mesh, "DG", 0) + markers = Function(W) + with eta.dat.vec_ro as eta_: + eta_max = eta_.max()[1] + + theta = 0.5 + should_refine = conditional(gt(eta, theta*eta_max), 1, 0) + markers.interpolate(should_refine) + + refined_mesh = mesh.refine_marked_elements(markers) + return refined_mesh + + rect1 = WorkPlane(Axes((0, 0, 0), n=Z, h=X)).Rectangle(1, 2).Face() + rect2 = WorkPlane(Axes((0, 1, 0), n=Z, h=X)).Rectangle(2, 1).Face() + L = rect1 + rect2 + + geo = OCCGeometry(L, dim=2) + ngmsh = geo.GenerateMesh(maxh=0.1) + mesh = Mesh(ngmsh) + + max_iterations = 10 + error_estimators = [] + dofs = [] + for i in range(max_iterations): + uh = solve_poisson(mesh) + (eta, error_est) = estimate_error(mesh, uh) + error_estimators.append(error_est) + dofs.append(uh.function_space().dim()) + mesh = adapt(mesh, eta) + assert error_estimators[-1] < 0.05 + + +@pytest.mark.skipnetgen +@pytest.mark.parallel +@pytest.mark.skipcomplex +def test_alfeld_stokes_netgen(): + ngmesh = square_geometry(0.75) + labels = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(codim=1)) if name == "rect"] + msh = Mesh(ngmesh, netgen_flags={"split": "Alfeld"}) + V = VectorFunctionSpace(msh, "CG", 2) + Q = FunctionSpace(msh, "DG", 1) + W = V*Q + u, p = TrialFunctions(W) + v, q = TestFunctions(W) + x, y = SpatialCoordinate(msh) + f = assemble(interpolate(as_vector([sin(x)*sin(y), sin(x)*sin(y)]), V)) + a = (inner(grad(u), grad(v)) - div(u)*q - div(v)*p)*dx + l = inner(f, v) * dx + w = Function(W) + bc = DirichletBC(W.sub(0), as_vector([0, 0]), labels) + solve(a == l, w, bcs=bc) + u = w.split()[0] + assert assemble(div(u)*div(u)*dx) < 1e-16