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
but when I use V = fem.functionspace(msh, ('CG', 2)) something weird is happened
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:
- 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.
- refine the mesh near the point of the patch.
- 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()