When I used stl file to create 3D mesh, which contained 10 boundary layers. When the boundary layers were not saved, it works well. But since the mesh file contained the boundary layers, it seemed not work. My code is below, I don’t know why it doesn’t work:
import gmsh
import sys
import os
import math
gmsh.initialize(sys.argv)
path = os.path.dirname(os.path.abspath(file))
gmsh.merge(os.path.join(path, ‘xuluwen.stl’))
gmsh.model.mesh.classifySurfaces(math.pi, True, True)
gmsh.model.mesh.createGeometry()
extrude = True
if extrude:
gmsh.option.setNumber(‘Geometry.ExtrudeReturnLateralEntities’, 0)
sur=gmsh.model.getEntities(2)
surface = [s[1] for s in sur]
print(surface)
e = gmsh.model.geo.extrudeBoundaryLayer(gmsh.model.getEntities(2), [10], [-1],True)
print(e)
# get “top” surfaces created by extrusion
top_ent = [s for s in e if s[0] == 2] #boundary surfaces(interior)
top_surf = [s[1] for s in top_ent]
print(top_surf)
# get boundary of top surfaces, i.e. boundaries of holes
gmsh.model.geo.synchronize()
bnd_ent = gmsh.model.getBoundary(top_ent) #boundary curves
bnd_curv = [c[1] for c in bnd_ent]
print(bnd_curv)
# create plane surfaces filling the holes
loops = gmsh.model.geo.addCurveLoops(bnd_curv)
print(loops)
for l in loops:
top_surf.append(gmsh.model.geo.addPlaneSurface([l]))
print(top_surf)
sur_boundary=gmsh.model.getEntities(2)
surface_boundary = [s[1] for s in sur_boundary]
print(surface_boundary)
splenic_marker=1
smv_marker=2
portal_marker=3
wall_marker=4
domain_marker=5
gmsh.model.geo.addPhysicalGroup(2,[211,42],splenic_marker)
gmsh.model.setPhysicalName(2,splenic_marker,‘splenic’)
gmsh.model.geo.addPhysicalGroup(2,[212,56,99],smv_marker)
gmsh.model.setPhysicalName(2,smv_marker,‘smv’)
gmsh.model.geo.addPhysicalGroup(2,[213,133,177],portal_marker)
gmsh.model.setPhysicalName(2,portal_marker,‘portal’)
gmsh.model.geo.addPhysicalGroup(2,surface,wall_marker)
gmsh.model.setPhysicalName(2,wall_marker,‘wall’)
gmsh.model.geo.addPhysicalGroup(3,[1,2,3,4,5,6,7,8,9,10,11],domain_marker)
gmsh.model.setPhysicalName(3,domain_marker,‘domain’)
gmsh.model.geo.addVolume([gmsh.model.geo.addSurfaceLoop(top_surf)])
gmsh.model.geo.synchronize()
gmsh.option.setNumber(‘Mesh.Algorithm’, 1)
gmsh.option.setNumber(‘Mesh.MeshSizeFactor’, 0.1)
gmsh.model.mesh.generate(3)
gmsh.write(‘msh_generate.msh’)
gmsh.finalize()
import numpy as np
import matplotlib.pyplot as plt
import tqdm.autonotebook
from mpi4py import MPI
from petsc4py import PETSc
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, create_matrix, set_bc)
from dolfinx.graph import adjacencylist
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)
from dolfinx.mesh import create_mesh, meshtags_from_entities
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement, as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
gdim = 3
mesh_comm = MPI.COMM_WORLD
model_rank = 0
mesh,ct, ft = gmshio.read_from_msh(“msh_generate.msh”, MPI.COMM_WORLD, 0, gdim=3)
