This is my mesh:
In fact using this mesh the facet_markers only include the global boundary:
msh, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)
But still I can access the internal subdomain boundaries by
tdim = msh.topology.dim
fdim = tdim - 1
boundary_facets = mesh.exterior_facet_indices(msh.topology)
This is an example of what I am doing:
import sys
import pyvista
import math
from dolfinx.io import XDMFFile, gmshio
from dolfinx.io.gmshio import model_to_mesh
from dolfinx.plot import vtk_mesh
sys.path.insert(0, 'meshes')
#import gmsh
import matplotlib.pyplot as plt
import numpy as np
import numpy.linalg as npla
from scipy.sparse.linalg import spsolve
import scipy.sparse as sps
import ufl
from dolfinx import fem
from dolfinx import default_real_type, mesh
from ufl import FacetNormal, SpatialCoordinate, ds, dx, grad, CellDiameter, inner, dot, as_ufl #Unified form language pieces
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement, TensorElement,
as_vector, div, dot, dS,ds, dx, inner, outer, lhs, grad, nabla_grad, rhs, sym, jump, avg)
from dolfinx.fem import (Expression, Constant, Function, FunctionSpace,
assemble_scalar, dirichletbc, form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
import numpy
import gmsh
from mpi4py import MPI
from petsc4py import PETSc
#import functions1
import dolfinx
from dolfinx.io import XDMFFile
#from dolfinx.plot import create_vtk_mesh
import dolfinx.io
from dolfinx import io
currentmesh = "mesh10" # Change N and the filename for different refinements
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gmsh.open(currentmesh + ".msh")
model = gmsh.model
model.setCurrent(currentmesh)
print('Model ' + model.getCurrent() + ' (' + str(model.getDimension()) + 'D)')
# Define mesh
msh, cell_markers, facet_markers = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)
n = FacetNormal(msh)
degree=1
V = FunctionSpace(msh, ("CG", degree))
tdim = msh.topology.dim
fdim = tdim - 1
u, v = TrialFunction(V), TestFunction(V)
msh.topology.create_connectivity(fdim, tdim)
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
(2, lambda x: np.isclose(x[0], 1)),
(3, lambda x: np.isclose(x[1], 0)),
(4, lambda x: np.isclose(x[1], 1))]
facet_indices1, facet_markers1 = [], []
fdim = msh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(msh, fdim, locator)
facet_indices1.append(facets)
facet_markers1.append(np.full_like(facets, marker))
facet_indices1 = np.hstack(facet_indices1).astype(np.int32)
facet_markers1 = np.hstack(facet_markers1).astype(np.int32)
sorted_facets = np.argsort(facet_indices1)
facet_tag = meshtags(msh, fdim, facet_indices1[sorted_facets], facet_markers1[sorted_facets])
msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
ds = Measure("ds", domain=msh, subdomain_data=facet_tag)
class BoundaryCondition():
def __init__(self, type, marker, values):
self._type = type
if type == "Dirichlet":
u_D = Function(V)
u_D.interpolate(values)
facets = facet_tag.find(marker)
dofs = locate_dofs_topological(V, fdim, facets)
self._bc = dirichletbc(u_D, dofs)
elif type == "Neumann":
self._bc = inner(values, v) * ds(marker)
elif type == "Robin":
self._bc = values[0] * inner(u-values[1], v)* ds(marker)
else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc
@property
def type(self):
return self._type
u_1 = Function(V)
u_1.interpolate(lambda x: 0*x[1]**0)
boundary_conditions = [BoundaryCondition("Dirichlet", 1, u_1),
BoundaryCondition("Dirichlet", 2, u_1),
BoundaryCondition("Dirichlet", 3, u_1),
BoundaryCondition("Dirichlet", 4, u_1),
]
#####################################################################################################################
jump_grad_u = dot(grad(u)('+'),n('+'))-dot(grad(u)('-'),n('-'))
jump_grad_v = dot(grad(v)('+'),n('+'))-dot(grad(v)('-'),n('-'))
jump_u = u('+')*n('+')+u('-')*n('-')
jump_v = v('+')*n('+')+v('-')*n('-')
epsilon=10.0**(-2)
alpha = dolfinx.fem.Constant(msh, 1.0)
h = CellDiameter(msh)
dS= Measure("dS", domain=msh, subdomain_data=facet_markers)
###################################################################################
a = epsilon*dot(grad(u), grad(v))*dx+ u * v * dx- epsilon*inner(avg(grad(v)), jump_u) * dS + epsilon*inner(jump_v, avg(grad(u))) * dS+epsilon*(alpha / avg(h))* epsilon* inner(jump_u, jump_v) * dS#+epsilon*(alpha / avg(h))* epsilon* inner(jump_u, jump_v) * dS
LL= dolfinx.fem.Constant(msh, 1.0) * v * dx
F=a-LL
bcs = []
for condition in boundary_conditions:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
else:
F += condition.bc
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
I want only integrate over internal subdomain boundaries.
