Gmsh, boundaries, Output,

[Noya]

Hello everybody,

I have a problem. I defined my geometry and mesh in gmsh.

SetFactory("OpenCASCADE");
L = 2;
R = 0.2;
r = 0.1;
Box(1) = {-L/2,-L/2,-L/2,L,L,L};  //tag=1
//+
Torus(2) = {0, 0, 0, R, r, Pi};  //tag=2
//+
Torus(3) = {0, 0, 0, R, r, Pi};  //tag=3
//+
Rotate {{1, 0, 0}, {0, 0, 0}, Pi/2}{
 Volume{2}; 
}
//+
Rotate {{1, 0, 0}, {0, 0, 0}, -Pi/2}{
 Volume{3}; 
}
Physical Volume(1) = {2, 3};	//Torus
//+
Physical Volume(2) = {1};	//surrounding
//+
Physical Surface(3) = {13, 16};    //internal boundary
//+
Physical Surface(4) = {1, 2, 3, 4, 5, 6};    //external boundary
Mesh 3;

But when I look at the integrals, the integral over the internal boundarie is 0.
What is the reason for this?

from fenics import*
import meshio
# Konvertiere Mesh zu XDMF für die Benutzung in FEniCS
msh = meshio.read("torus.msh")
tet_data = msh.cell_data_dict["gmsh:physical"]["tetra"]
tri_data = msh.cell_data_dict["gmsh:physical"]["triangle"]

meshio.write("mesht.xdmf",
    meshio.Mesh(points=msh.points,
        cells={"tetra": msh.cells_dict["tetra"]}
    )
)
meshio.write("randt.xdmf",
    meshio.Mesh(points=msh.points,
        cells={"triangle": msh.cells_dict["triangle"]},
        cell_data={"bnd_marker": [msh.cell_data_dict["gmsh:physical"]["triangle"]]}
    )
)
meshio.write("gebiett.xdmf",
     meshio.Mesh(
        points=msh.points,
        cells={"tetra": msh.cells_dict["tetra"]},
        cell_data={"dom_marker": [msh.get_cell_data("gmsh:physical", "tetra")]}
    )
)


mesh = Mesh()
with XDMFFile("mesht.xdmf") as mshfile:
    mshfile.read(mesh)



mft = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
with XDMFFile("gebiett.xdmf") as domfile:
    domfile.read(mft)
File("mft.pvd") << mft
dx = Measure('dx', domain=mesh, subdomain_data=mft)

submesh_torus = SubMesh(mesh, mft, 1)
submesh_OT = SubMesh(mesh, mft, 2)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("randt.xdmf") as bndfile:
    bndfile.read(mvc, "bnd_marker")
mff = MeshFunction('size_t', mesh, mvc)
for i in range(mff.size()):
    if mff[i] not in (3, 4):
        mff[i] = 0
File("mff.pvd") << mff
dS = Measure('dS', domain=mesh, subdomain_data=mff) # internal facets
ds = Measure("ds", domain=mesh, subdomain_data=mff) #external facets

print(assemble(Constant(1)*dx))
print(assemble(Constant(1)*dx(1)))
print(assemble(Constant(1)*dx(2)))
print(assemble(Constant(1)*ds(4)))
print(assemble(Constant(1)*dS(3)))

8.031561967647798
0.03156196764777432
8.000000000000007
23.999999999999968
0.0

I have oriented myself on this Mesh, Submesh , Torus, Boundary - #3 by conpierce8 but the outputs of “mesht”, “randt” and “gebiett” were different from the xdmf files created there if you ignore the “refinement” of the mesh.

left my mesht right mesh

left my gebiett right domains

left my randt right boundaries

But why? Can someone please have a look and give me a tip to fix it?

Thank you and best regards,
noya

[dokken]

This is because you need to use Boolean Fragments in your geo file. Currently, the box and the torus is meshed has individual meshes.
In the reply

Mesh, Submesh , Torus, Boundary

Hi @Noya, You may want to try creating your mesh using an external meshing tool, e.g. Gmsh. The advantage of this approach is that you can mark the subdomains and boundaries using Gmsh, save the markers to files, and load them into FEniCS MeshFunctions, which can be used to define integration measures and boundary conditions. Using your current code (i.e. with BoxMesh), you obtain a mesh that doesn’t respect the internal boundary of the torus, i.e. you will have elements that pass through the t…

This is done with the command:

@conpierce8

gmsh.model.occ.fragment([(3,3)], [(3,1),(3,2)])

I’ve also shown how to use it with the geo files in: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #62 by dokken

You can realize this mistake by looking at your output. You get the volume 8 by integrating dx(2) which is the whole box.

[Noya]

Ah okay.
originally I had BooleanFragments{ Volume{1}, Volume{2}, Volume{3}; Delete; }{} but that didn’t change anything, so I deleted it again.
This means for me to add between the geometry and the physical volumes and surfaces the following

BooleanFragments{ Volume{1}; Delete;}{ Volume{2,3}; Delete; };

?
But then I’ve got the error:

line 44, in
bndfile.read(mvc, “bnd_marker”)
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** fenics-support@googlegroups.com
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to find entity in map.
*** Reason: Error reading MeshValueCollection.
*** Where: This error was encountered inside HDF5File.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

This error is the same as in your link.
In your solution you used

v() = BooleanFragments{ Volume{1}; Delete;}{ Volume{2,3,4}; Delete; };
Physical Volume("air", 33) ={#v()};
Physical Volume("hBN", 32) ={#v()-1};

What means or what do {#v()}; ,{#v()-1}; ? I couldn’t find it in the gmsh documentation. Can you explain me that?

Edit : I found it in tutorial/t16.geo · gmsh_4_8_4 · gmsh / gmsh · GitLab but I don’t know if I understand it correctly.

The tag of the volume has changed by “BooleanFragments”, so we need to access it programmatically.

with {#v()} I access the whole cube? Is this correct?
But how I access on my torus?

best regards,
noya

[Noya]

hello,

I think I understand the problem above.
There are more meshes and they are overlapping each other!
So I edit my geo-file and have

SetFactory("OpenCASCADE");

L = 2;
R = 0.2;
r = 0.1;
hmin = r/2;
hmax = L/10;

Box(1) = {-L/2,-L/2,-L/2,L,L,L};  //tag=1
//+
Torus(2) = {0, 0, 0, R, r, Pi};  //tag=2
//+
Torus(3) = {0, 0, 0, R, r, Pi};  //tag=3
//+
Rotate {{1, 0, 0}, {0, 0, 0}, Pi/2}{
  Volume{2}; 
}
//+
Rotate {{1, 0, 0}, {0, 0, 0}, -Pi/2}{
  Volume{3}; 
}
BooleanFragments{ Volume{1,2,3}; Delete;};
Physical Volume(1) = {2,3};	//Torus
//+
Physical Volume(2) = {1};	//surrounding
//+
Physical Surface(3) = {13,16};    //boundary torus
//+
Physical Surface(4) = {1, 2, 3, 4, 5, 6};    //boundary cube
Mesh 3;

I have no error but the integral over the cube without the torus print(assemble(Constant(1)*dx(2))) is still \approx 8

What do I wrong?

Best regards,
noya

[dokken]

Please read through the gmsh documentation.
Here is a working gmsh code making the boolean fragments you would like:

SetFactory("OpenCASCADE");

L = 2;
R = 0.3;
r = 0.2;

Box(1) = {-L/2,-L/2,-L/2,L,L,L};  //tag=1
Torus(2) = {0, 0, 0, R, r, 2*Pi};  //tag=2


v() = BooleanFragments{ Volume{1, 2}; Delete; }{};
Physical Volume(22) = {v[0]};
Physical Volume(23) = {v[1]};

BooleanUnion{ Volume{2}; Delete; }{ Volume{3}; Delete; }

[Noya]

Thank you very much.
I think I got the geo-file right.