It doesn’t work after all, but I don’t understand what the error is. Here are the two pictures. Does anyone understand what is wrong?
How do I find out which points collides with the bb_tree?
So these points from the example:
tol = 0.001 # Avoid hitting the outside of the domain
y = np.linspace(-1 + tol, 1 - tol, 101)
points = np.zeros((3, 101))
points[1] = y
MWE:
#Mesh Erzeugung:(mit GMSH)
from mpi4py import MPI
from dolfinx import mesh
import warnings
warnings.filterwarnings("ignore")
import gmsh
gmsh.initialize()
from dolfinx import fem
from dolfinx import nls
from dolfinx import log
from dolfinx import io
from dolfinx import geometry
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx.fem import Constant, Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc
from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector, create_vector, set_bc
from dolfinx.graph import create_adjacencylist
from dolfinx.geometry import BoundingBoxTree, compute_collisions, compute_colliding_cells
from dolfinx.io import (XDMFFile, cell_perm_gmsh, distribute_entity_data, extract_gmsh_geometry,
extract_gmsh_topology_and_markers, ufl_mesh_from_gmsh)
from dolfinx.mesh import create_mesh, meshtags_from_entities
import ufl
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
import numpy as np
import matplotlib.pyplot as plt
import tqdm.notebook
from mpi4py import MPI
from petsc4py import PETSc
import scipy.constants as c
#Domain-Definieren:
L = 1.
H = 1.
c_x = 0.50
c_y = 0.50
r = 0.01
gdim = 2
rank = MPI.COMM_WORLD.rank
if rank == 0:
rectangle = gmsh.model.occ.addRectangle(0,0,0, L, H, tag=1)
obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)
if rank == 0:
Mesh_domain = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
gmsh.model.occ.synchronize()
#Verbinden aller Mesh-Punkte
Mesh_domain_marker = 1
if rank == 0:
volumes = gmsh.model.getEntities(dim=gdim)
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], Mesh_domain_marker)
gmsh.model.setPhysicalName(volumes[0][0], Mesh_domain_marker, "Mesh_domain")
inlet_marker, outlet_marker, wall_marker, obstacle_marker, wall_top = 2, 3, 4, 5, 6
inflow, outflow, walls, obstacle, top = [], [], [], [], []
if rank == 0: #Walls sind Ground, Obstical ist Leiter
boundaries = gmsh.model.getBoundary(volumes, oriented=False)
for boundary in boundaries:
center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
if np.allclose(center_of_mass, [0, H/2, 0]):
inflow.append(boundary[1])
elif np.allclose(center_of_mass, [L, H/2, 0]):
outflow.append(boundary[1])
elif np.allclose(center_of_mass, [L/2, 0, 0]):
walls.append(boundary[1])
elif np.allclose(center_of_mass, [L/2, H, 0]):
top.append(boundary[1])
else:
obstacle.append(boundary[1])
gmsh.model.addPhysicalGroup(1, walls, wall_marker)
gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")
gmsh.model.addPhysicalGroup(1, top, wall_top)
gmsh.model.setPhysicalName(1, wall_top, "top")
if rank == 0:
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 8)
gmsh.option.setNumber("Mesh.RecombineAll", 2)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.optimize("Netgen")
if MPI.COMM_WORLD.rank == 0:
# Mesh Geometrie
x = extract_gmsh_geometry(gmsh.model)
# Topologie fĂĽr jedes Element vom Mesh
topologies = extract_gmsh_topology_and_markers(gmsh.model)
num_cell_types = len(topologies.keys())
cell_information = {}
cell_dimensions = np.zeros(num_cell_types, dtype=np.int32)
for i, element in enumerate(topologies.keys()):
properties = gmsh.model.mesh.getElementProperties(element)
name, dim, order, num_nodes, local_coords, _ = properties
cell_information[i] = {"id": element, "dim": dim, "num_nodes": num_nodes}
cell_dimensions[i] = dim
# Sort elements by ascending dimension
perm_sort = np.argsort(cell_dimensions)
# Broadcast cell type data and geometric dimension
cell_id = cell_information[perm_sort[-1]]["id"]
tdim = cell_information[perm_sort[-1]]["dim"]
num_nodes = cell_information[perm_sort[-1]]["num_nodes"]
cell_id, num_nodes = MPI.COMM_WORLD.bcast([cell_id, num_nodes], root=0)
if tdim - 1 in cell_dimensions:
num_facet_nodes = MPI.COMM_WORLD.bcast( cell_information[perm_sort[-2]]["num_nodes"], root=0)
gmsh_facet_id = cell_information[perm_sort[-2]]["id"]
marked_facets = np.asarray(topologies[gmsh_facet_id]["topology"], dtype=np.int64)
facet_values = np.asarray(topologies[gmsh_facet_id]["cell_data"], dtype=np.int32)
cells = np.asarray(topologies[cell_id]["topology"], dtype=np.int64)
cell_values = np.asarray(topologies[cell_id]["cell_data"], dtype=np.int32)
else:
cell_id, num_nodes = MPI.COMM_WORLD.bcast([None, None], root=0)
cells, x = np.empty([0, num_nodes], np.int64), np.empty([0, gdim])
cell_values = np.empty((0,), dtype=np.int32)
num_facet_nodes = MPI.COMM_WORLD.bcast(None, root=0)
marked_facets = np.empty((0, num_facet_nodes), dtype=np.int64)
facet_values = np.empty((0,), dtype=np.int32)
# Create distributed mesh
ufl_domain = ufl_mesh_from_gmsh(cell_id, gdim)
gmsh_cell_perm = cell_perm_gmsh(to_type(str(ufl_domain.ufl_cell())), num_nodes)
cells = cells[:, gmsh_cell_perm]
mesh = create_mesh(MPI.COMM_WORLD, cells, x[:, :gdim], ufl_domain)
tdim = mesh.topology.dim
fdim = tdim - 1
# Permute facets from MSH to DOLFINx ordering
facet_type = cell_entity_type(to_type(str(ufl_domain.ufl_cell())), fdim, 0)
gmsh_facet_perm = cell_perm_gmsh(facet_type, num_facet_nodes)
marked_facets = np.asarray(marked_facets[:, gmsh_facet_perm], dtype=np.int64)
local_entities, local_values = distribute_entity_data(mesh, fdim, marked_facets, facet_values)
mesh.topology.create_connectivity(fdim, tdim)
adj = create_adjacencylist(local_entities)
bb_tree = geometry.BoundingBoxTree(mesh, mesh.topology.dim)
# Create DOLFINx MeshTags
ft = meshtags_from_entities(mesh, fdim, adj, np.int32(local_values))
ft.name = "Facet tags"
#Funktionspace und Vektorfunktionspace defenieren
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, s_cg1)
V_E= fem.VectorFunctionSpace(mesh, ("DG", 0))
fdim = mesh.topology.dim - 1
#Konstanten, Funktionen und Randbedingungen definieren:
phi_0 = 80000 #V
h = c_y
u1_0 = Constant(mesh, PETSc.ScalarType(phi_0))
#Boundary conditions definieren
# Walls
u1_nonslip = np.array((0,) *mesh.geometry.dim, dtype=PETSc.ScalarType)
u11_nonslip = np.array((0,) *mesh.geometry.dim, dtype=PETSc.ScalarType)
wall_facets = ft.indices[ft.values == wall_marker]
bcu1_walls = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(V, fdim, wall_facets), V)
# Obstacle
obstacle_facets = ft.indices[ft.values == obstacle_marker]
bcu1_obstacle = dirichletbc(PETSc.ScalarType(8e4), locate_dofs_topological(V, fdim, obstacle_facets), V)
#Boundary conditions for equations:
bcu1 = [bcu1_obstacle, bcu1_walls] # Spannung am Leiter und Ränder
#Laplace
u1 = ufl.TrialFunction(V)
v1 = ufl.TestFunction(V)
k = Constant(mesh, PETSc.ScalarType(0))
f = k/c.epsilon_0
epsilon = c.epsilon_0
a = ufl.inner(ufl.grad(u1),ufl.grad(v1)) * ufl.dx
L = ufl.inner(f , v1) * ufl.dx
u1 = fem.Function(V)
problem1 = fem.petsc.LinearProblem(a, L, bcu1, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
phi1 = problem1.solve()
with io.VTKFile(mesh.comm, "Laplace.pvd", "w") as vtk:
vtk.write([phi1._cpp_object])
with io.XDMFFile(mesh.comm, "Laplace.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(phi1)
###https://jorgensd.github.io/dolfinx-tutorial/chapter1/membrane_code.html#making-curve-plots-throughout-the-domain
tol = 0.001
y = np.linspace(-1 + tol, 1 - tol, 101)
points = np.zeros((3, 101))
points[1] = y
u_values = []
p_values = []
cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions(bb_tree, points.T)
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(mesh, cell_candidates, points.T)
for i, point in enumerate(points.T):
if len(colliding_cells.links(i))>0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
u_values = phi1.eval(points_on_proc, cells)
import matplotlib.pyplot as plt
fig = plt.figure()
plt.plot(points_on_proc[:,1], u_values, "k", linewidth=2, label="Phi")
plt.grid(True)
plt.xlabel("x")
plt.legend()
# If run in parallel as a python file, we save a plot per processor
plt.savefig(f"membrane_rank{MPI.COMM_WORLD.rank:d}.png")
print("finish")
It works in version 0.4.1