Hello,
I have a question about Neumann boundary conditions.
I have a 2D plot in a room where the bottom is the ground with Dirichlet BC = 0 and 2 electrical conductors with a positive voltage and a negative voltage. The surface of the conductors should contain the Neumann boundary condition. I can’t get both conditions to be executed.
Here is the code. Unfortunately it is a bit big.
But only the bottom part is important.
#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.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 import geometry
from dolfinx.mesh import create_mesh, meshtags_from_entities
import ufl
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad
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 =24.
H =18.
c_x1 = 8.
c_y1 = 13.
r1 = 0.01
c_x2 = 16.
c_y2 = 13.
r2 = 0.001
gdim = 2
rank = MPI.COMM_WORLD.rank
if rank == 0:
rectangle = gmsh.model.occ.addRectangle(0,0,0, L, H, tag=1)
obstacle1 = gmsh.model.occ.addDisk(c_x1, c_y1, 0, r1, r1)
obstacle2 = gmsh.model.occ.addDisk(c_x2, c_y2, 0, r2, r2)
if rank == 0:
Mesh_domain = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle1), (gdim, obstacle2)])
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, obstacle1_marker,obstacle2_marker, wall_top = 2, 3, 4, 5, 6, 7
inflow, outflow, walls, obstacle1, obstacle2 , top = [], [], [], [], [], []
if rank == 0:
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])
elif np.allclose(center_of_mass, [c_x1, c_y1, 0]):
obstacle1.append(boundary[1])
elif np.allclose(center_of_mass, [c_x2, c_y2, 0]):
obstacle2.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, obstacle1, obstacle1_marker)
gmsh.model.setPhysicalName(1, obstacle1_marker, "Obstacle1")
gmsh.model.addPhysicalGroup(1, obstacle2, obstacle2_marker)
gmsh.model.setPhysicalName(1, obstacle2_marker, "Obstacle2")
gmsh.model.addPhysicalGroup(1, top, wall_top)
gmsh.model.setPhysicalName(1, wall_top, "top")
if rank == 0:
gmsh.option.setNumber('Mesh.MeshSizeMin', 2e-1)
gmsh.option.setNumber('Mesh.MeshSizeMax', 2e-1)
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
E_cp = -16.5e5
E_cn = 16.5e5
#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]
bc_walls = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(V, fdim, wall_facets), V)
#Boundary conditions for equations:
bcwall = [bc_walls]
#Laplace
u1 = ufl.TrialFunction(V)
v1 = ufl.TestFunction(V)
dObt = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=obstacle1_marker)
dObs = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=obstacle2_marker)
a = ufl.inner(ufl.grad(u1), ufl.grad(v1)) * ufl.dx
x = SpatialCoordinate(mesh)
L = ufl.inner(E_cp, v1) * dObs + ufl.inner(E_cn,v1) * dObt
problem1 = fem.petsc.LinearProblem(a, L, bcwall, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
phi = problem1.solve()
print("phi_max: ", np.max(phi.x.array), " / phi_min: ", np.min(phi.x.array))
print(phi.x.array)
with io.VTKFile(mesh.comm, "Laplace.pvd", "w") as vtk:
vtk.write([phi._cpp_object])
with io.XDMFFile(mesh.comm, "Laplace.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(phi)