Here is the code:
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 (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
L = 24.
H = 18.
c_x1 = 8.
c_y1 = 13.
r1 = 0.001
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()
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])
else:
obstacle2.append(boundary[1])
obstacle1.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.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:
x = extract_gmsh_geometry(gmsh.model)
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
perm_sort = np.argsort(cell_dimensions)
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)
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
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)
ft = meshtags_from_entities(mesh, fdim, adj, np.int32(local_values))
ft.name = "Facet tags"
print(ft.indices[ft.values == obstacle1_marker])
print(ft.indices[ft.values == obstacle2_marker])
I have deleted some of it I can’t get it much smaller