I thought that this would overcomplicate the code. Also, I can’t provide the exact mesh as it involves sensitive data.
import basix
import pyvista as pv
import numpy as np
from dolfinx.mesh import create_mesh, exterior_facet_indices, meshtags
from dolfinx.fem import *
from dolfinx import fem
from mpi4py import MPI
from scipy.spatial import cKDTree
import ufl
def get_original_index_mapping(mesh, original_points):
# get new coordinates of nodes from the DOLFINx mesh
new_coordinates = mesh.geometry.x
# build a KDTree for the new coordinates
tree = cKDTree(new_coordinates)
# query the KDTree to find the nearest points in new coordinates for each original point
_, nearest_indices = tree.query(original_points)
# create a mapping dictionary
original_to_new_mapping = [nearest_indices[i].astype(int) for i in range(len(original_points))]
return np.array(original_to_new_mapping)
def find_closest_node(points, target_node):
# convert the target_node to a numpy array (if it's a list or tuple)
target_node = np.array(target_node)
# calculate Euclidean distances between the target node and all other points
distances = np.linalg.norm(points - target_node, axis=1)
# find the index of the minimum distance
closest_node_index = np.argmin(distances)
return closest_node_index
# read mesh with pyvista and get nodes and cells
mesh_init = pv.read('mesh.vtk')
nodes = np.asarray(mesh_init.points)
cells = np.array([mesh_init.get_cell(i).point_ids for i in range(mesh_init.n_cells)])
set_conductivities = np.array([0.0, 0.3, 1/10, 0.7, 1/50, 1/40])
conductivities = set_conductivities[mesh_init.cell_data['material']]
# read electrode coordinates and get corresponding nodes
electrodes = np.loadtxt('electrodes/electrodes.txt')
# last coordinates belong to reference electrode
ref_electrode = electrodes[-1]
electrodes = electrodes[:-1]
e_idx = []
for e in electrodes:
idx = np.where(np.all(nodes==e, axis=1))[0]
if idx.size==1:
e_idx.append(int(idx))
else:
idx = find_closest_node(nodes, e)
e_idx.append(idx)
electrode_nodes = np.array(e_idx).reshape(4,16)
ref_electrode_node = np.where(np.all(nodes==ref_electrode,axis=1))[0]
if len(ref_electrode_node)==0:
for i, node in enumerate(nodes):
if np.allclose(node, ref_electrode):
ref_electrode_node = [i]
break
# set source, sink, curr
ref_node = ref_electrode_node[0]
curr = 1.
sink = electrode_nodes[0,0]
source = electrode_nodes[0,1]
source_coord = nodes[source].reshape(1,3)
sink_coord = nodes[sink].reshape(1,3)
ground_coord = nodes[ref_node].reshape(1,3)
# create mesh and compute all connections
c_el = ufl.Mesh(basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(nodes.shape[1],)))
mesh = create_mesh(MPI.COMM_WORLD, cells, nodes, c_el)
mesh.topology.create_connectivity(mesh.topology.dim, 2)
mesh.topology.create_connectivity(mesh.topology.dim-1, 1)
mesh.topology.create_connectivity(mesh.topology.dim-2, 0)
mesh.topology.create_connectivity(0, mesh.topology.dim-1)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
node_map = get_original_index_mapping(mesh, nodes) # from dolfinx mesh to orignal; return idx of node in dolfinx mesh
# define function space
V = fem.functionspace(mesh, ("CG", 1))
# define DG0 function space for conductivity
DG0 = fem.functionspace(mesh, ("DG", 0))
sigma = fem.Function(DG0)
sigma.x.array[:] = conductivities[mesh.topology.original_cell_index]
# find corresponding dof (source, sink, ref_nod and electrodes are just ids)
source_dof = node_map[source]
sink_dof = node_map[sink]
ground_dof = node_map[ref_node]
electrodes_dof = []
for e in electrode_nodes:
electrodes_dof.append(node_map[e])
source_facet_indices = mesh.topology.connectivity(mesh.topology.dim - 1, 0).links(source_dof)
sink_facet_indices = mesh.topology.connectivity(mesh.topology.dim - 1, 0).links(sink_dof)
source_facet_indices = [i for i in source_facet_indices if i in boundary_facets]
sink_facet_indices = [i for i in sink_facet_indices if i in boundary_facets]
# define boundary conditions for potential at source, sink, and ground
bc_source = fem.dirichletbc(np.array(curr, dtype=np.float64), np.asarray(source_dof.astype(np.int32)).reshape(1), V)
bc_sink = fem.dirichletbc(np.array(-curr, dtype=np.float64), np.asarray(sink_dof.astype(np.int32)).reshape(1), V)
bc_ground = fem.dirichletbc(np.array(0., dtype=np.float64), np.asarray(ground_dof.astype(np.int32)).reshape(1), V)
# define the trial and test functions for the weak form
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# define Poisson's equation weak form
a = ufl.dot(sigma *ufl.grad(u), ufl.grad(v)) * ufl.dx
f = fem.Constant(mesh, 0.)
L = f * v * ufl.dx
# use provided facet IDs for Neumann boundary condition
facet_ids = np.concatenate([source_facet_indices, sink_facet_indices])
facet_tags = np.concatenate([np.full(len(source_facet_indices), 1, dtype=np.int32), np.full(len(sink_facet_indices), 2, dtype=np.int32)])
facet_tag = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, facet_ids, facet_tags)
# create a measure for the boundary
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag)
# add Neumann boundary condition to the weak form
L += curr * v * ds(1)
L += -curr * v * ds(2)
# solve the linear system
problem = LinearProblem(a, L, bcs=[bc_sink, bc_source],
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
uh = problem.solve()
I don’t get any error message, the values don’t match the excepted ones.
