I have just read the gmsh manual and verified that the physical tags to
surfaces are indeed getting assigned properly. The relevant lines are 25-31
in the msh file provided above.
I have arranged these lines in a table and verified that the physical tags
are getting assigned exactly as I prescribed in the geo file. The picture
shows the first 9 space-separated entries of lines 25-31 of cube.msh. I
have inserted an additional column commenting about the plane's location.
[image: cube_surface_tags.png]
According to the msh file, both x normal surfaces have tag 2, both y normal
surfaces have tag 3 and z normal surfaces have tags 4 and 5 -- exactly as
in the geo file.
So, I think the issue is with the procedure I am employing to read the
mesh. I am attaching the main.cc (provided in the MWE above) here again,
inline. I am unable to understand why the tags for y normal and z normal
surfaces get swapped when I read and print them using this code. I used a
similar code for 2D meshes and it worked fine. Are there any tricky things
to consider additionally in 3D?
Vachan
----------------------------------------------------------------------------
#include <fstream>
#include <string>
#include <deal.II/base/mpi.h>
#include <deal.II/base/utilities.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/grid/grid_in.h>
using namespace dealii;
int main(int argc, char **argv)
{
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
constexpr int dim = 3;
MPI_Comm mpi_comm(MPI_COMM_WORLD);
parallel::distributed::Triangulation<dim> triang(mpi_comm);
GridIn<dim> grid_in;
grid_in.attach_triangulation(triang);
std::ifstream file("cube.msh");
grid_in.read_msh(file);
const int periodic_dir = std::stoi(argv[1]); // first command line
argument
for(auto cell: triang.active_cell_iterators()){
const auto face_1 = cell->face(2*periodic_dir);
const auto face_2 = cell->face(2*periodic_dir+1);
if(face_1->at_boundary()){
std::cout << "Found face1 of local id " << 2*periodic_dir << "
at boundary with bid "
<< face_1->boundary_id() << "\n";
}
if(face_2->at_boundary()){
std::cout << "Found face2 of local id " << 2*periodic_dir+1 <<
" at boundary with bid "
<< face_2->boundary_id() << "\n";
}
}
return 0;
}
On Wednesday, May 5, 2021 at 11:45:21 AM UTC+5:30 [email protected]
wrote:
> Daniel,
>
> Thank you for responding.
>
> I am not familiar with the msh file format. But I see the 'Entities'
> section. Definitely, some info is being assigned, but I am not sure if it
> is correct.
>
> I am attaching a minimal working example containing 4 files: main.cc,
> cube.geo, cube.msh and CMakeLists.txt. The executable (named a.out in cmake
> file) takes one command line argument for the direction. I am using gmsh
> 4.3.0.
>
> The code reads 'cube.msh' and loops over all faces with local ids
> 2*direction and 2*direction+1. If any of them is on the boundary, it prints
> the boundary id. The geo file is supposed to assign boundary id 2 to x
> direction faces, 3 to y direction faces and 4,5 to z- (normal in negative z
> direction) and z+ faces.
>
> The output for x faces is as expected, but the boundary ids assigned to y
> normal and z normal faces are somehow swapped.
>
> Please let me know if any additional information is required.
>
> On Tuesday, May 4, 2021 at 8:24:18 PM UTC+5:30 [email protected] wrote:
>
>> Vachan,
>>
>> It is difficult to say what goes wrong here with the information you
>> provided. How does the msh file created by gmsh look like. Is the boundary
>> information correct there?
>> Which gmsh version are you using? Can you post a complete minimal example
>> of your code?
>>
>> Best,
>> Daniel
>>
>> Am Di., 4. Mai 2021 um 05:56 Uhr schrieb [email protected] <
>> [email protected]>:
>>
>>> Dear all,
>>>
>>> I have encountered a strange problem while dealing with physical
>>> boundary ids. I have a simple cube mesh (annexure 1) which I am exporting
>>> to msh4 format for reading through GridIn. The mesh is constructed by
>>> extruding a square 2D mesh in x-y plane to 3D.
>>>
>>> I was facing problems with using collect_periodic_faces()
>>> <https://www.dealii.org/8.5.0/doxygen/deal.II/namespaceGridTools.html#add6c6052ecc070b8f6ee1489b55897a7>
>>> on
>>> this mesh. So I mimicked the function by writing my own test function
>>> (annexure 2). It loops over faces with ids 2*direction and 2*direction+1
>>> (with direction=2) which lie on a boundary and prints the boundary id.
>>>
>>> According to the construction of the mesh (see geo file) face1 must have
>>> boundary id 4 and face 2 must have boundary id 5 in the test function
>>> (coincidentally, matching with the cell-local id). However, when I execute
>>> the test, the output shows both faces have a boundary id 0.
>>>
>>> More interestingly, if I uncomment the lines in geo file and set
>>> boundary ids for all boundaries (in the trailing lines), then the code
>>> execution shows boundary ids of y-normal and z-normal faces swapped; i.e.;
>>> face1 and face2 show boundary id 3 which was assigned to y-normal faces
>>> (see the output below).
>>>
>>> Found face1 of local id 4 at boundary with bid 3
>>> Found face2 of local id 5 at boundary with bid 3
>>>
>>> If I instead define face1 and face2 with direction=1 using
>>>
>>> const auto face_1 = cell->face(2);
>>> const auto face_2 = cell->face(3);
>>>
>>> in the test function and modify the local id in print statements, then
>>> the boundary ids printed are 4 and 5 which were assigned to z-normal faces.
>>>
>>> Found face1 of local id 2 at boundary with bid 4
>>> Found face2 of local id 3 at boundary with bid 5
>>>
>>> The x normal faces lying on boundary have the correct boundary id. So it
>>> looks like y and z axes are getting swapped?
>>>
>>> Am I doing something wrong? Has anyone experienced this before? I am
>>> unable to find an explanation for this behaviour.
>>>
>>> Thanking in anticipation
>>> Vachan
>>>
>>> *Annexure 1*
>>> ------------------------------------------------------------------------
>>> // cube.geo file
>>> SetFactory("OpenCASCADE");
>>>
>>> // Simple cube geometry
>>> lc=0.1;
>>> l=1;
>>>
>>> // 2D geometry
>>> p1=newp; Point(p1) = {0,0,0,lc};
>>> p2=newp; Point(p2) = {l,0,0,lc};
>>> p3=newp; Point(p3) = {l,l,0,lc};
>>> p4=newp; Point(p4) = {0,l,0,lc};
>>>
>>> l1=newl; Line(l1) = {p1,p2};
>>> l2=newl; Line(l2) = {p2,p3};
>>> l3=newl; Line(l3) = {p3,p4};
>>> l4=newl; Line(l4) = {p4,p1};
>>>
>>> Transfinite Curve{l1,l2,l3,l4} = 10;
>>>
>>> ll1=newll; Curve Loop(ll1) = {l1,l2,l3,l4};
>>> s1=news; Plane Surface(s1) = {ll1};
>>> Transfinite Surface{s1} = {p1,p2,p3,p4};
>>> Recombine Surface{s1};
>>>
>>> // 3D extrusion
>>> out[] = Extrude {0,0,l} {
>>> Surface{s1}; Layers{10}; Recombine;
>>> };
>>>
>>> // Physical entities
>>> Physical Volume(1) = {out[1]};
>>> // Physical Surface(2) = {out[3],out[5]}; // x dir faces
>>> // Physical Surface(3) = {out[2],out[4]}; // y dir faces
>>> Physical Surface(4) = {s1}; // z- dir face
>>> Physical Surface(5) = {out[0]}; // z+ dir face
>>>
>>> *Annexure 2*
>>> --------------------------------------------------------------
>>> // snippet from test function mimicking collect_periodic_faces()
>>> for(auto cell: problem.triang.active_cell_iterators()){
>>> const auto face_1 = cell->face(4);
>>> const auto face_2 = cell->face(5);
>>> if(face_1->at_boundary()){
>>> std::cout << "Found face1 of local id 4 at boundary with bid "
>>> << face_1->boundary_id()
>>> << "\n";
>>> }
>>> if(face_2->at_boundary()){
>>> std::cout << "Found face2 of local id 5 at boundary with bid "
>>> << face_2->boundary_id()
>>> << "\n";
>>> }
>>> }
>>>
>>> --
>>> The deal.II project is located at http://www.dealii.org/
>>> For mailing list/forum options, see
>>> https://groups.google.com/d/forum/dealii?hl=en
>>> ---
>>> You received this message because you are subscribed to the Google
>>> Groups "deal.II User Group" group.
>>> To unsubscribe from this group and stop receiving emails from it, send
>>> an email to [email protected].
>>> To view this discussion on the web visit
>>> https://groups.google.com/d/msgid/dealii/97ddd5ab-18be-4ce2-8ecf-b0d77cd70f9dn%40googlegroups.com
>>>
>>> <https://groups.google.com/d/msgid/dealii/97ddd5ab-18be-4ce2-8ecf-b0d77cd70f9dn%40googlegroups.com?utm_medium=email&utm_source=footer>
>>> .
>>>
>>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/9197bb5d-7484-4316-9949-e7a6cd70b38fn%40googlegroups.com.
