I am currently trying to solve the poisson equation for neumann boundary conditions and I seem to get seemingly incorrect answers. I am inexperienced with differential equation and FEM and I would hope that there could be some light shed on this topic.
I am trying to solve a simple zero neumann boundary condition poisson equation with a dirac delta source in a circle.
In order for a neumann boundary condition PDE to be well defined, there are two conditions that need to be met.
- The integral of the source must be equal to the integral of the boundary function g ( referring to the neumann poisson). This is due to the physics of the problem.
- Either the integral of the solution u must be 0 or some point must be known. This is due to fact that the solution is up to a constant and therefore some point must be fixed to a value.
According to the neumann poisson post, if you enforce the solution of u to be zero you sort of meet both conditions described. Unfortunately following the recipe on thsi post is not possible on fenicsx so instead I attempt to remove the nullspace manually from the linear system.
my code looks like this
import dolfinx
import dolfinx.fem as fem
import numpy as np
import ufl
from mpi4py import MPI
from petsc4py.PETSc import ScalarType as default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from dolfinx.fem import petsc
# Define group names mapping
GROUP_NAMES = {"Cancellous Bone": 1, "Cortical Bone": 2, "Muscle": 3, "Fat": 4, "Skin": 5, "Boundary": 6}
# Define conductivity constants
ANISOTROPY_RATIO = 5
CONDUCTIVITY = {
"Cancellous Bone": 0.075,
"Cortical Bone": 0.02,
"Fat": 0.0379,
"Skin": 4.55E-4,
"Muscle": np.diag(
[
0.2455,
ANISOTROPY_RATIO * 0.2455
]),
}
from mpi4py import MPI
from dolfinx import mesh,plot
from dolfinx.fem import FunctionSpace,Function,form
from dolfinx.fem.petsc import create_vector,assemble_matrix,assemble_vector
from ufl import inner,TrialFunction,TestFunction,ds,dx,grad,exp,sin,SpatialCoordinate
from petsc4py import PETSc
# To solve pure neumann conditions PDE, we constraint the solutions
# to sum to zero as in Neurodec
# Alternative options include using a lagrange multiplier in legacy dolfin
class ConstrainedLinearProblem():
def __init__(self, a, L, V):
self.a = a
self.L = L
self.V = V
def solve(self):
A = assemble_matrix(form(self.a))
A.assemble()
b=create_vector(form(self.L))
with b.localForm() as b_loc:
b_loc.set(0)
assemble_vector(b,form(self.L))
# Solution Function
uh = Function(self.V)
# Create Krylov solver
solver = PETSc.KSP().create(A.getComm())
solver.setOperators(A)
# Create vector that spans the null space
nullspace = PETSc.NullSpace().create(constant=True,comm=MPI.COMM_WORLD)
A.setNullSpace(nullspace)
# orthogonalize b with respect to the nullspace ensures that
# b does not contain any component in the nullspace
nullspace.remove(b)
# Finally we are able to solve our linear system ::
solver.solve(b,uh.vector)
return uh
class FixedPointProblem():
def __init__(self, a, L, V, bcs):
self.a = a
self.L = L
self.V = V
self.bc = bcs
def solve(self):
A = petsc.assemble_matrix(fem.form(self.a), bcs=[self.bc])
A.assemble()
b = petsc.assemble_vector(fem.form(self.L))
fem.apply_lifting(b, [fem.form(self.a)], [[self.bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.set_bc(b, [self.bc])
solver = PETSc.KSP().create(MPI.COMM_WORLD)
solver.setOperators(A)
# solver.setType(PETSc.KSP.Type.PREONLY)
# solver.getPC().setType(PETSc.PC.Type.LU)
solver.setFromOptions()
uh = fem.Function(self.V)
solver.solve(b, uh.vector)
uh.x.scatter_forward()
return uh
class DiskFEMModel:
def __init__(self, msh_file_path: str):
self.mesh, self.cell_markers, self.facet_markers = dolfinx.io.gmshio.read_from_msh(msh_file_path, MPI.COMM_WORLD, gdim=2)
self.mesh_geometry = self.mesh.geometry
self.mesh_topology = self.mesh.topology
self.build_model()
def build_model(self):
# Build out index references
self.tree = dolfinx.geometry.bb_tree(self.mesh, self.mesh.topology.dim)
self.midpoints = dolfinx.geometry.create_midpoint_tree(self.mesh, self.mesh.topology.dim, np.arange(self.cell_markers.values.size, dtype=np.int32))
# Link the 1D and 3D topology
self.mesh_topology.create_connectivity(self.mesh.topology.dim, 0)
self.cell_to_vertex = self.mesh_topology.connectivity(self.mesh.topology.dim, 0)
# Define function spaces
element = ufl.TensorElement("DG", self.mesh.ufl_cell(), degree=0, shape=(2, 2))
self.V_tensor = fem.functionspace(self.mesh, element)
#self.V_tensor = fem.TensorFunctionSpace(self.mesh, ("DG", 1), shape=(3, 3))
self.V_scalar = fem.FunctionSpace(self.mesh, ("CG", 1))
self.V_source = fem.FunctionSpace(self.mesh, ("CG", 1))
# Assign conductivity
self.build_conductivity_map()
def build_conductivity_map(self):
self.sigma_anisotropic = fem.Function(self.V_tensor)
num_cells = len(self.mesh.topology.connectivity(2, 0))
with self.sigma_anisotropic.vector.localForm() as lf:
for i in range(num_cells):
lf.setValuesBlocked([i], CONDUCTIVITY['Muscle'].flatten())
def extend_array(self, x):
zeros_column = np.zeros((x.shape[0], 1))
extended_array = np.hstack((x, zeros_column))
return extended_array
# Maybe create a more generic way of doing this instead of repeating code.
def evaluate_solution_at_points(self, points: np.array):
points = self.extend_array(points)
cell_ids = dolfinx.geometry.compute_closest_entity(self.tree, self.midpoints, self.mesh, points).squeeze()
return self.uh.eval(points, cell_ids)
def evaluate_f_at_points(self, points: np.array):
points = self.extend_array(points)
cell_ids = dolfinx.geometry.compute_closest_entity(self.tree, self.midpoints, self.mesh, points).squeeze()
return self.source_function.eval(points, cell_ids)
def point_to_cell_id(self, point: np.array):
# Get the closest cell ID to the point
cell_id = dolfinx.geometry.compute_closest_entity(self.tree, self.midpoints, self.mesh, point).squeeze()
# Check if point is in the cell
# vertex_indices = self.cell_to_vertex.links(cell_id)
# vertices = self.mesh.geometry.x[vertex_indices]
# min_bounds = vertices.min(axis=0)
# max_bounds = vertices.max(axis=0)
# assert np.all(point >= min_bounds) and np.all(point <= max_bounds), "Point co-ordinate's are not inside the volume."
return cell_id
def point_to_vertex_locations(self, point: np.array):
cell_id = self.point_to_cell_id(point)
vertex_indices = self.cell_to_vertex.links(cell_id)
return np.array([self.mesh_geometry.x[vertex_index] for vertex_index in vertex_indices])
def point_to_cell_material(self, point: np.array):
cell_id = self.point_to_cell_id(point)
material_marker = self.cell_markers.values[cell_id]
material_name = {v: k for k, v in GROUP_NAMES.items()}[material_marker]
return material_name
def assign_source_to_point(self, point: np.array, source_value: float = 1.0):
# Initialize the source function
self.source_function = fem.Function(self.V_source)
with self.source_function.vector.localForm() as local:
local.set(0.0)
def gaussian_nd(x, x_mu, sigma):
"""
Calculate the n-dimensional Gaussian function for a given mean vector and scalar variance.
:param x: ndarray, where each row represents a coordinate in n-dimensional space.
:param x_mu: ndarray, the mean vector in n-dimensional space.
:param sigma: float, the standard deviation (sqrt of variance).
:return: ndarray, values of the Gaussian function at each point in x.
"""
# Ensure x and x_mu are numpy arrays for broadcasting
# x = np.array(x)
# x_mu = np.array(x_mu)
x = x.T
# Calculate the squared Euclidean distance from each point in x to the mean
squared_dist = np.sum((x - x_mu)**2, axis=-1)
# Compute the Gaussian function
return 1 / ((2 * np.pi * sigma**2) ** (x.shape[-1] / 2)) * np.exp(-squared_dist / (2 * sigma**2))
self.source_function.interpolate(lambda x: gaussian_nd(x, point, 0.1))
integral = fem.assemble_scalar(fem.form(self.source_function * dx))
print(integral)
#self.source_function.vector.axpy(-integral, self.source_function.vector)
u_constant = fem.Constant(self.mesh, default_scalar_type(1.))
volume = fem.assemble_scalar(fem.form(u_constant * dx))
#self.source_function.interpolate(lambda x: gaussian_nd(x, point, 0.1)-integral/volume)
self.source_integral = fem.assemble_scalar(fem.form(self.source_function * dx))
def fix_cell_value(self, point, value = 1.0):
cell_id = self.point_to_cell_id(point)
vertex_indices = self.cell_to_vertex.links(cell_id)
u_D = fem.Function(self.V_scalar)
with u_D.vector.localForm() as loc:
loc.set(0)
u_D.vector[vertex_indices[0]] = value
u_D.vector[vertex_indices[1]] = value
u_D.vector[vertex_indices[2]] = value
boundary_dofs = vertex_indices
bc = fem.dirichletbc(u_D, boundary_dofs)
return bc
def solve_for_point(self, point: np.array, source_value: float = 1.0):
# Create the point source
self.assign_source_to_point(point, source_value)
# Define the test and trial functions of weak form
u = ufl.TrialFunction(self.V_scalar)
v = ufl.TestFunction(self.V_scalar)
g = fem.Constant(self.mesh, default_scalar_type(self.source_integral)) # Neumann boundary condition
# Define the variational problem
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
#a = ufl.dot(ufl.dot(self.sigma_anisotropic, ufl.grad(u)), ufl.grad(v)) * ufl.dx
L = self.source_function * v * ufl.dx - g * v * ufl.ds
# Apply Dirichlet boundary condition
def on_boundary(x):
return np.isclose(np.sqrt(x[0]**2 + x[1]**2), 1)
#boundary_facets = facet_indices[np.where(facet_values == GROUP_NAMES["Boundary"])]
boundary_dofs = fem.locate_dofs_geometrical(self.V_scalar, on_boundary)
u_bc = fem.Function(self.V_scalar)
with u_bc.vector.localForm() as loc:
loc.set(5.0)
boundary_condition = fem.dirichletbc(u_bc, boundary_dofs)
# Pass to dolfinx and solve
print(point)
# bcs = self.fix_cell_value(point[0])
problem = ConstrainedLinearProblem(a, L, self.V_scalar)
#problem = LinearProblem(a, L, petsc_options={"ksp_type": "preonly", "pc_type": "lu"} )
#problem = FixedPointProblem(a, L, self.V_scalar, bcs=bcs)
self.uh = problem.solve()
return np.array(self.uh.vector)
# if __name__ == "__main__":
# model = DiskFEMModel("cylinders.msh")
# point = np.expand_dims(np.array([[20, 10, 10]], dtype=np.float64), 0)
# cell_material = model.point_to_cell_material(point[0])
# print(cell_material, "\n", model.point_to_conductivity(point[0]))
# uh = model.solve_for_point(point)
# print(uh)
# print(np.nonzero(uh))
Pay attention at the ConstrainedProblem definition that essentially adds the removal of the nullspace from a regular linear problem.
The main topic of my question is the following. If I use this solution to solve the PDE my results look strange. Here follows a source function and the solution just by ensuring the integral of u is zero using the nullspace trick.
solution: image0sol.png - Google Drive
Now these are the results when I enforce the source function to have integral 0 as per the 1st condition for the pure neumann condition poisson.
solution: image1sol.png - Google Drive
I am not sure if both results are correct. I can understand just enforcing the 2nd condition can be enough to have a well-posed problem, but what does this mean physically and why do the results look so different.
Another question is how does this solution differ from the solution using fenics and a lagrangian constraint.
Thank you