When I do assemble(1.0*dS(domain=mesh,subdomain_id=2, subdomain_data=boundary_markers)) I also get 0 instead of the area of the surface. When I don’t specify domain=mesh I get an error.
For a MWE, here’s the file.geo that one compiles with gmsh file.geo -3:
Mesh.RandomFactor=1.0e-4;
Mesh.Algorithm = 7;
Mesh.Algorithm3D= 7;
l = 4;
w = 4;
t = 3;
l_diss = 4;
w_diss = 4;
t_diss = 1;
res_xy = 3e-01;
res_z = 1e-01;
Point(1) = {0.0, 0.0, 0.0, res_xy};
Point(2) = {l, 0.0, 0.0, res_xy};
Point(3) = {l, w, 0.0, res_xy};
Point(4) = {0.0, w, 0.0, res_xy};
Point(5) = {0.0, 0.0, t, res_xy};
Point(6) = {l, 0.0, t, res_xy};
Point(7) = {l, w, t, res_xy};
Point(8) = {0.0, w, t, res_xy};
Point(9) = {0.0, 0.0, t + t_diss, res_xy};
Point(10) = {l_diss, 0.0, t + t_diss, res_xy};
Point(11) = {l_diss, w_diss, t+t_diss, res_xy};
Point(12) = {0.0, w_diss, t + t_diss, res_xy};
Line(100) = {1,2};
Line(101) = {2,3};
Line(102) = {3,4};
Line(103) = {4,1};
Line(200) = {5,6};
Line(201) = {6,7};
Line(202) = {7,8};
Line(203) = {8,5};
Line(300) = {9,10};
Line(301) = {10,11};
Line(302) = {11,12};
Line(303) = {12,9};
Line(501) = {1,5};
Line(502) = {2,6};
Line(503) = {3,7};
Line(504) = {4,8};
Line(601) = {5,9};
Line(602) = {6,10};
Line(603) = {7,11};
Line(604) = {8,12};
Line Loop(1000) = {100,101,102,103};
Line Loop(2000) = {200,201,202,203};
Line Loop(3000) = {300,301,302,303};
Line Loop(5000) = {100,502,-200,-501};
Line Loop(5100) = {101,503,-201,-502};
Line Loop(5200) = {102,504,-202,-503};
Line Loop(5300) = {103,501,-203,-504};
Line Loop(6000) = {200,602,-300,-601};
Line Loop(6100) = {201,603,-301,-602};
Line Loop(6200) = {202,604,-302,-603};
Line Loop(6300) = {203,601,-303,-604};
Plane Surface(1001) = {1000};
Plane Surface(2001) = {2000};
Plane Surface(5001) = {5000};
Plane Surface(5101) = {5100};
Plane Surface(5201) = {5200};
Plane Surface(5301) = {5300};
Plane Surface(3001) = {3000};
Plane Surface(6001) = {6000};
Plane Surface(6101) = {6100};
Plane Surface(6201) = {6200};
Plane Surface(6301) = {6300};
Surface Loop(10000) = {1001,2001,5001,5101,5201,5301,3001,6001,6101,6201,6301};
Volume(10001) = {10000};
Field[1] = MathEvalAniso;
Field[1].m11 = "1/(10e-1)^2";
Field[1].m12 = "0.0";
Field[1].m13 = "0.0";
Field[1].m22 = "1/(10e-1)^2";
Field[1].m23 = "0.0";
Field[1].m33 = "1/(5e-1)^2";
Background Field = 1;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthFromCurvature = 0;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
Physical Surface("bottom_0") = {1001};
Physical Surface("top_0") = {2001};
Physical Surface("top_1") = {3001};
Physical Volume("All") = {10001};
Then the Python file:
import meshio
msh = meshio.read("simple_gmsh_mesh.msh")
for key in msh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:physical"][key]
elif key == "tetra":
tetra_data = msh.cell_data_dict["gmsh:physical"][key]
for cell in msh.cells:
if cell.type == "tetra":
tetra_cells = cell.data
elif cell.type == "triangle":
triangle_cells = cell.data
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},
cell_data={"name_to_read":[tetra_data]})
triangle_mesh =meshio.Mesh(points=msh.points,
cells=[("triangle", triangle_cells)],
cell_data={"name_to_read":[triangle_data]})
meshio.write("file2.xdmf", tetra_mesh)
meshio.write("mf2.xdmf", triangle_mesh)
mesh = Mesh()
with XDMFFile("file2.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf2.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc2 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("file2.xdmf") as infile:
infile.read(mvc2, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
dx = Measure("dx", domain=mesh,subdomain_data=cf)
ds = Measure("ds", domain=mesh, subdomain_data=mf)
subdomains = cf
boundary_parts = mf
V = FunctionSpace(mesh, 'CG', 1)
T_bott, T_r, tol, blokus_height, top_topus = 273+60.0, 273+10.0, 1E-14, 3, 1
class Bottom(SubDomain):
def inside(self, x , on_boundary):
return on_boundary and near(x[2], 0.0, tol)
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[2] >= 0.0 + tol and x[2] <= blokus_height - tol
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[2] >= 3.63e-3 + tol and x[2] <= top_topus - tol
class top_blokus(SubDomain):
def inside(self, x , on_boundary):
return on_boundary and near(x[2], blokus_height, tol)
class top_dis(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], blokus_height + top_topus, tol)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(0)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
subdomain_1 = top_blokus()
subdomain_1.mark(boundary_markers, 1)
subdomain_2 = top_dis()
subdomain_2.mark(boundary_markers, 2)
materials = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
subdomain_10 = Omega_0()
subdomain_11 = Omega_1()
subdomain_10.mark(materials, 10)
subdomain_11.mark(materials, 11)
k_0 = 1.0
k_1 = 2.0
t_f = 10
num_steps = 10
dt = t_f / num_steps
class K(UserExpression):
def __init__(self, materials, k_0, k_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.k_0 = k_0
self.k_1 = k_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.k_0
else:
values[0] = self.k_1
def value_shape(self):
return ()
kappa = K(materials, k_0, k_1, degree=1)
u_D0 = Constant(T_bott)
bc0 = DirichletBC(V, u_D0, Bottom())
u_0 = Expression('T_r', T_r = T_r, degree = 1)
u_n = interpolate(u_0, V)
u = TrialFunction(V)
v = TestFunction(V)
F = dt * dot(kappa*grad(u), grad(v)) * dx + v *(u-u_n)*dx
a, L = lhs(F), rhs(F)
u = Function(V)
t = 0
for n in range(num_steps):
u_D0.t = t
t += dt
solve(a == L, u, bcs = bc0)
u_n.assign(u)
print('the temp at the origin is ', u([0,0,0]))
print('the temp at the top is ', u([0,0,blokus_height + top_topus]))
print('debug line: ', assemble(avg(u)*dS(domain=mesh, subdomain_id=2, subdomain_data=boundary_markers)))
This is much longer than 20 lines… I do hope it’s not too much.
