Dear teacher:

- I am a learner in FEniCS and Gmsh. Now I want to integral my function in Lines. So I give two domains, this is my Gmsh code:

//Gmsh Code:

Geometry.OCCTargetUnit = “M”;

mesh1 = 0.04;

mesh2 = 0.1;

L = 1;

rate = 1/4;

//Subdomain1

p0 = newp; Point(p0) = {0.1,0,0,mesh1};

p1 = newp; Point(p1) = {0.45,rate*L/2/2,0,mesh1};
p2 = newp; Point(p2) = {0.1,rate*L/2/2,0,mesh1};

p3 = newp; Point(p3) = {0.1,-rate

*L/2/2,0,mesh1};*

p4 = newp; Point(p4) = {0.45,-rateL/2/2,0,mesh1};

p4 = newp; Point(p4) = {0.45,-rate

l1 = newl; Line(l1) = {p1,p2};

l2 = newl; Line(l2) = {p2,p3};

l3 = newl; Line(l3) = {p3,p4};

l4 = newl; Line(l4) = {p4,p1};

ll1 = newll; Line loop(ll1) = {l1,l2,l3,l4};

s1 = news; Surface(s1) = {ll1};

//Subdomain2

p10 = newp; Point(p10) = {0,0,0,mesh2};

p11 = newp; Point(p11) = {L/2,rate*L/2,0,mesh2};
p12 = newp; Point(p12) = {-L/2,rate*L/2,0,mesh2};

p13 = newp; Point(p13) = {-L/2,-rate

*L/2,0,mesh2};*

p14 = newp; Point(p14) = {L/2,-rateL/2,0,mesh2};

p14 = newp; Point(p14) = {L/2,-rate

l11 = newl; Line(l11) = {p11,p12};

l12 = newl; Line(l12) = {p12,p13};

l13 = newl; Line(l13) = {p13,p14};

l14 = newl; Line(l14) = {p14,p11};

ll11 = newll; Line loop(ll11) = {l11,l12,l13,l14};

s11 = news; Surface(s11) = {ll11, ll1};

Physical Surface(1) = {s11,s1};

Physical Line(1) = {l1};

Physical Line(2) = {l11};

- Now I want to integral in Lines, So I first put “mesh.geo” into FEniCS,the fenics code is:

#fenics Code:

import meshio

from dolfin import *

import matplotlib.pyplot as plt

import numpy

#Mesh

mesh_from_file = meshio.read(“111.msh”)

def create_mesh(mesh, cell_type, prune_z=False):

cells = numpy.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])

cell_data = numpy.hstack([mesh.cell_data_dict[“gmsh:physical”][key]

for key in mesh.cell_data_dict[“gmsh:physical”].keys() if key==cell_type])

# Remove z-coordinates from mesh if we have a 2D cell and all points have the same third coordinate

points= mesh.points

if prune_z:

points = points[:,:2]

mesh_new = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={“name_to_read”:[cell_data]})

return mesh_new

line_mesh = create_mesh(mesh_from_file, “line”, prune_z=True)

meshio.write(“facet_mesh.xdmf”, line_mesh)

triangle_mesh = create_mesh(mesh_from_file, “triangle”, prune_z=True)

meshio.write(“mesh.xdmf”, triangle_mesh)

mesh = Mesh()

with XDMFFile(“mesh.xdmf”) as infile:

infile.read(mesh)

mvc = MeshValueCollection(“size_t”, mesh, 1)

with XDMFFile(“facet_mesh.xdmf”) as infile:

infile.read(mvc, “name_to_read”)

mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

ds1 = Measure(“ds”, domain=mesh, subdomain_data=mf)#边界的识别

plot(mesh,title=‘Mesh’)

#Integral

print(‘ds1(1)=’,assemble(1*ds1(1)))
print(‘ds1(2)=’,assemble(1*ds1(2)))

- My Question is:When I want to integral in Inner Line(ds1(1)), the answer is 0;

When I want to integral in Outer Line(ds1(2)), the answer is 1;It seems that I failed to identify the inner line 1 in fenics. So if I want to integral in ds(1), what can I do?

Thanks for your answer!