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.
