### How to extract SubMesh and corresponding MeshFunction from a mesh imported from Gmsh

108
views
0
10 weeks ago by
Dear All,

I am solving a 3D coupled fluid and structure problem. Now I have a full mesh and I want divide it into two submesh:
filer.read(*mesh,std::string("mesh"),false);
filer.read(*cell_marker, std::string("subdomains_mark")); //MeshFunction for subdomains_mark with label=1 (fluid) and 2 (solid).
filer.read(*bd_marker, std::string("facet_mark")); // MeshFunction for facet_mark with some labels representing different boundaries.​

then I can get the SubMesh and MeshFuntion by:

fluid_mesh.reset( new SubMesh(*mesh, *cell_mark, 1)); //Here label=1 means taking fluid mesh.
fluid_cell_marker.reset( new MeshFunction<std::size_t>(fluid_mesh, fluid_mesh->topology().dim()));
fluid_bd_marker.reset( new MeshFunction<std::size_t>(fluid_mesh, fluid_mesh->topology().dim()-1));​

But, fluid_cell_marker and fluid_bd_marker are not consistent with cell_marker and bd_marker, i.e., the same boundary has different boundary labels in mesh and submesh. How to solve this problem?

I found a similar post https://www.allanswered.com/post/nqzzq/#lozzv. In this post, the author was trying to impose a dirichlet boundary condition on a boundary. I think the difficulty is that this boundary is labeled as 3 in the full mesh, however we do not know what the label is in submesh.

Thanks!
Leon

Community: FEniCS Project
Can you post your .geo file? IT seems like you either have forgotten to add physical groups or you are using a label twice in GMSH. (For starter try to use the command 'newreg' in GMSH, instead of manually labeling them. This will prevent repetitive registers.) Again, I could be more helpful if I saw your .geo file.
written 10 weeks ago by Nima

Dear Nima,

Here is the script of .geo. If you need all the files, I will upload later today.

Merge "inlet.stl";
Merge "outlet.stl";
Merge "innerwall.stl";
Merge "outerwall.stl";
Merge "inlet_wall.stl";
Merge "outlet_wall.stl";

CreateTopology;
//Surface Loop(5) = {1};
//Surface Loop(6) = {2};
//Surface Loop(7) = {3};
//Surface Loop(8) = {4};
//Volume(9) = {5, 6, 7, 8};

Surface Loop(6) = {3, 1, 2};
Volume(7) = {6};
Surface Loop(8) = {3, 4, 5,6};
Volume(9) = {8};

//+
Physical Surface("innerwall",1) = {3};
//+
Physical Surface("outerwall",2) = {4};
//+
Physical Surface("inlet",3) = {1};
//+
Physical Surface("inlet_wall",4) = {5};
//+
Physical Surface("outlet",5) = {2};
//+
Physical Surface("outlet_wall",6) = {6};

//+
Physical Volume("wall", 2) = {9};
//+
Physical Volume("blood", 1) = {7};


Thank you very much!
Leon

written 10 weeks ago by Leon
One more comment:

When I check the file fluid_cell_marker, I found the labels of cells are assigned randomly.

I am trying to use fluid_cell_marker.data().array('parent_vertex_indices', 0)  to find a map from the full mesh and submesh, then I can mark them manually.
written 10 weeks ago by Leon
I am looking at this on my phone and your .geo script looks ok. (I would choose 7 and 8 for "wall" and "blood", to not have the same markers as innerwall and outerwall, just in case.)
Also I have a question. So when you mesh this region using GMSH it will save a .msh file with the same name as your .geo file. Have you tried converting it using 'dolfin-convert' command?
if you are using linux the command is:

dolfin-convert yourfile.msh yourfile.xml​
this will create 3 files for you.
1- yourfile.xml
2- yourfile_physical_region.xml
3- yourfile_facet_region.xml

then you can read your markers directly from these files by double clicking on them and checking the long list. Also reading from these files is less difficult.
written 10 weeks ago by Nima
Yes, of course. I also use those three .xml files to build a .h5 file and then use HDF5File filer to read the mesh data. For the full mesh, it works well.
My question is about the Submesh, when I tried to extract a submesh from the full mesh by using
fluid_mesh.reset( new SubMesh(*mesh, *cell_mark, 8));

I do not know how to get the corresponding CellFunction and FacetFunction. Although I can use

fluid_cell_marker.reset( new MeshFunction<std::size_t>(fluid_mesh, fluid_mesh->topology().dim()));
fluid_bd_marker.reset( new MeshFunction<std::size_t>(fluid_mesh, fluid_mesh->topology().dim()-1));​

to create such MeshFunctions, the labels in these MeshFunctions are random, not consistent with the full mesh.

written 10 weeks ago by Leon