Polynomial level of function space changes the grid?

Hey,
I am trying to solve the Laplace equation outside a sphere, where the boundary condition is Neumann over all of the sphere except for small patch where the potential is constant.
not sure if the problem is in the visualization, but when I use the V = fem.functionspace(msh, ('CG', 1)) the output looks correct
image

but when I use V = fem.functionspace(msh, ('CG', 2)) something weird is happened
image

I tried to use the second-order because the flux calculation of with the first-order didn’t matched the theory.

the outline of the code before the solver:

  1. creating a big sphere and a small one so I can try to simulate the decay of the function like 1/r at long distance.
  2. refine the mesh near the point of the patch.
  3. find the dof which will be the “patch” with the Dirichlet bc (I assumed the problem is here but I tried to take only the dofs on the inner boundery)
import dolfinx.io.gmshio
import numpy as np
from mpi4py import MPI
import gmsh
import ufl
from dolfinx import fem, io, plot, mesh
from dolfinx.fem import dirichletbc, locate_dofs_geometrical, FunctionSpace, Constant
from dolfinx.fem.petsc import LinearProblem
from ufl import TrialFunction, TestFunction, grad, inner, ds, dx
import pyvista
from dolfinx.plot import vtk_mesh as create_vtk_mesh
from utils import *

shell_volume_tag = 1
inner_boundary_tag = 2
outer_boundary_tag = 3
REFINEMENT_RADIUS = 0.1

def generate_spherical_shell_mesh(a, R, inner_mesh_size, outer_mesh_size):
    gmsh.initialize()
    gmsh.model.add("spherical_shell")
    # Create the inner sphere (volume)
    inner_sphere = gmsh.model.occ.addSphere(0, 0, 0, a)
    gmsh.model.occ.synchronize()
    # Create the outer sphere (volume)
    outer_sphere = gmsh.model.occ.addSphere(0, 0, 0, R)
    gmsh.model.occ.synchronize()
    # Create the spherical shell by cutting the inner sphere from the outer sphere
    shell = gmsh.model.occ.cut([(3, outer_sphere)], [(3, inner_sphere)])
    gmsh.model.occ.synchronize()
    # Get the volume tag of the spherical shell
    shell_volume = shell[0][0][1]
    # Get the surfaces of the spherical shell
    shell_surfaces = gmsh.model.occ.getEntities(dim=2)
    gmsh.model.occ.synchronize()
    # Add physical groups
    gmsh.model.addPhysicalGroup(3, [shell_volume], tag=shell_volume_tag)  # Spherical shell volume
    gmsh.model.setPhysicalName(3, 1, "ShellVolume")
    gmsh.model.addPhysicalGroup(2, [inner_sphere], tag=inner_boundary_tag)  # Inner boundary
    gmsh.model.setPhysicalName(2, 1, "InnerBoundary")
    gmsh.model.addPhysicalGroup(2, [outer_sphere], tag=outer_boundary_tag)  # Outer boundary
    gmsh.model.setPhysicalName(2, 2, "OuterBoundary")

    gmsh.model.occ.synchronize()
    # Set mesh sizes on the surfaces
    gmsh.model.mesh.setSize(
        gmsh.model.getBoundary([(2, inner_sphere)], recursive=True),
        inner_mesh_size
    )
    gmsh.model.mesh.setSize(
        gmsh.model.getBoundary([(2, outer_sphere)], recursive=True),
        outer_mesh_size
    )

    # Mesh generation
    gmsh.model.mesh.generate(3)
    # Convert to DOLFINx mesh
    msh, cell_tags, facet_tags = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)
    gmsh.finalize()

    return msh, cell_tags, facet_tags


def refine_mesh(msh, facet_tags, points, number_of_iterations):
    print('Start refinment')
    for receptor_center in points:
        dim = msh.topology.dim
        # ------------------------------
        # Refinement in DOLFINx
        # ------------------------------
        # Create connectivity for edges (dim-2 to dim-1)
        msh.topology.create_connectivity(dim - 1, dim)  # Create facet-cell connectivity
        msh.topology.create_connectivity(dim - 2, dim)  # Create edge-facet connectivity (for refinement)

        refinement_radius = REFINEMENT_RADIUS

        # Define the receptor area function
        def receptor_point(x):
            dist_squared = (x[0] - receptor_center[0]) ** 2 + (x[1] - receptor_center[1]) ** 2 + (
                    x[2] - receptor_center[2]) ** 2
            return dist_squared < refinement_radius ** 2

        for n in range(number_of_iterations):  # How many times to refine.
            refinement_radius = 0.8 * refinement_radius
            # Locate edges near the receptor, dim=1 is the edge dim.
            edges = mesh.locate_entities(msh, 1, receptor_point)

            # Ensure mesh has connectivity between facets and cells
            msh.topology.create_entities(dim - 1)
            f_to_c = msh.topology.connectivity(dim - 1, dim)

            facet_indices = []

            # Find facets on the boundary (those connected to only one cell)
            for f in range(msh.topology.index_map(dim - 1).size_local):
                if len(f_to_c.links(f)) == 1:
                    facet_indices.append(f)

            # Create meshtag for refinement
            meshtag = mesh.meshtags(
                msh,
                dim - 1,
                np.array(facet_indices, dtype=np.int32),
                facet_tags.values,
            )

            # Refine the mesh using mesh.refine_plaza
            msh, parent_cell, parent_facet = mesh.refine_plaza(
                msh, edges, redistribute=False, option=mesh.RefinementOption.parent_cell_and_facet
            )
            msh.topology.create_entities(dim - 1)

            # Transfer meshtags to the refined mesh
            facet_tags = mesh.transfer_meshtag(meshtag, msh, parent_cell, parent_facet)

            msh.topology.create_connectivity(msh.topology.dim - 1, msh.topology.dim)
    print('Finished refinment')
    return msh, facet_tags


def apply_receptor_bc(V, receptor_centers, receptor_radius, phi_receptor, facet_tags):
    """
    Applies Dirichlet boundary conditions at receptor locations on the inner boundary.
    Parameters:
        V (FunctionSpace): Function space.
        receptor_centers (list of tuples): Coordinates of receptor centers.
        receptor_radius (float): Radius of receptors.
        phi_receptor (float): Potential at receptors.
        facet_tags (MeshTags): Facet tags of the mesh.
    Returns:
        list: List of DirichletBC objects.
    """
    bcs = []
    # Locate all facets on the inner boundary
    inner_facets = facet_tags.find(inner_boundary_tag)

    for i, receptor_center in enumerate(receptor_centers):
        print(f"Applying boundary condition for receptor {i + 1} at {receptor_center}")

        def receptor_area(x):
            dist_squared = (x[0] - receptor_center[0])**2 + (x[1] - receptor_center[1])**2 + (x[2] - receptor_center[2])**2
            return dist_squared < receptor_radius**2

        # Locate DoFs on the inner boundary within the receptor area
        receptor_dofs = fem.locate_dofs_geometrical(V, receptor_area)
        # take only the DOF on the inner sphere
        receptor_dofs = np.array(list(set(receptor_dofs) & set(inner_facets)))

        if len(receptor_dofs) > 0:
            print(f"DoFs found for receptor {i + 1}: {len(receptor_dofs)}")
            bc_receptor = dirichletbc(phi_receptor, receptor_dofs, V)
            bcs.append(bc_receptor)
        else:
            print(f"No DoFs found for receptor {i + 1}, check mesh refinement.")

    return bcs


def main():
    # Mesh parameters
    a = 1.0  # Cell radius
    R = 10.0  # Outer boundary radius
    inner_mesh_size = 0.05  # Refined mesh near receptor
    outer_mesh_size = 1.0

    # Receptor parameters
    receptor_radius = 0.01  # Physical radius of the receptor on the sphere
    receptor_centers = [(0.0, 0.0, a)]  # Receptor at the north pole
    phi_receptor = 1.0

    number_of_refinment_near_point = 5

    msh, cell_tags, facet_tags = generate_spherical_shell_mesh(a, R, inner_mesh_size, outer_mesh_size)
    msh, facet_tags = refine_mesh(msh, facet_tags, receptor_centers, number_of_refinment_near_point)

    V = fem.functionspace(msh, ('CG', 1))
    # Define measures with facet_tags
    inner_tag = inner_boundary_tag
    outer_tag = outer_boundary_tag
    ds = ufl.Measure("ds", domain=msh, subdomain_data=facet_tags)
    dx = ufl.dx(domain=msh)

    # Define trial and test functions
    u = TrialFunction(V)
    v = TestFunction(V)

    # Define Dirichlet condition at the receptor area
    bcs = apply_receptor_bc(V, receptor_centers, receptor_radius, phi_receptor, facet_tags)

    # Define Neumann condition (zero flux) on the inner boundary
    # Neumann conditions with zero flux should be naturally handled, but we include it explicitly for clarity
    g_N = Constant(msh, 0.0)

    # Define the variational forms
    a_form = inner(grad(u), grad(v)) * dx + (1 / R) * u * v * ds(outer_tag)
    L_form = g_N * v * ds(inner_tag) + Constant(msh, 0.0) * v * dx

    # Set up and solve the linear problem
    problem = LinearProblem(
        a_form,
        L_form,
        bcs=bcs,
        petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
    )
    uh = problem.solve()

    # flux calculation
    n = ufl.FacetNormal(msh)
    flux_form_inner = inner(grad(uh), n) * ds(inner_tag)
    total_flux_inside = fem.assemble_scalar(fem.form(flux_form_inner))
    MPI.COMM_WORLD.barrier()
    # Compute the analytical flux for comparison
    flux_inner_exact = 4 * receptor_radius * phi_receptor
    print(f"Numerical flux over the inner sphere: {total_flux_inside:.6f}")
    print(f"Analytical flux over the inner sphere: {flux_inner_exact:.6f}")
    print(f"Relative error: {abs(total_flux_inside - flux_inner_exact) / abs(flux_inner_exact):.6e}")

    # Visualization using PyVista
    import pyvista
    # Convert mesh and solution to PyVista format
    cells, types, x = plot.vtk_mesh(V)
    grid = pyvista.UnstructuredGrid(cells, types, x)
    grid.point_data["uh"] = uh.x.array.real
    grid.set_active_scalars("uh")
    # Highlight receptor area
    # Points where the potential is equal to phi_receptor are considered receptor points
    receptor_points = grid.points[np.isclose(grid.point_data["uh"], phi_receptor, atol=1e-3)]
    # Plot the mesh and solution
    plotter = pyvista.Plotter()
    plotter.add_mesh(grid, show_edges=True, scalars="uh", cmap="coolwarm")
    plotter.add_points(receptor_points, color="red", point_size=15, label="Receptor")
    plotter.add_scalar_bar(title="Potential")
    plotter.add_legend()
    plotter.show()

if __name__ == "__main__":
    main()