Re: [deal.II] Gmsh unexpected behaviour - axes directions swapped

[vachanpo...@gmail.com]

Daniel,

As you rightly said, the documentation very clearly states how the 
direction parameter would be used. But the presence of a rotation in the 
mapping physical and reference cell was unexpected for me considering that 
the mesh was Cartesian. So the tricky part here is how the mesh is being 
read.. May be one small comment on this too would be helpful.

Thanks for the clarification and help!


On Friday, May 7, 2021 at 6:40:48 PM UTC+5:30 d.arnd...@gmail.com wrote:

> Vachan,
>
> The documentation for the second version says
>
> "Instead of defining a 'first' and 'second' boundary with the help of two 
> boundary_ids this function defines a 'left' boundary as all faces with 
> local face index 2*direction and boundary indicator b_id and, similarly, 
> a 'right' boundary consisting of all face with local face index 
> 2*direction+1 and boundary indicator b_id. Faces with coordinates only 
> differing in the direction component are identified."
>
> So there is an assumption about the local face index and the direction 
> parameter that is also used for identifying faces. This implies that a case 
> like yours where the cell is rotated is not supported by this overload. 
> Just as you said, it's not enough for the cell to be in standard 
> orientation.
>
> I'm sorry that this caused so much trouble but I'm glad that the first 
> version (that I would probably always recommend using anyway) works for you.
>
> Best,
> Daniel
>
>
>
> Am Fr., 7. Mai 2021 um 02:58 Uhr schrieb vachanpo...@gmail.com <
> vachanpo...@gmail.com>:
>
>> I think I have found an explanation.
>>
>> Since the mesh generated is a cartesian mesh, I was assuming that the 
>> faces of physical and reference cells have normals in the same direction 
>> and the mapping between then would be a simple stretching. But when I 
>> printed out vectors connecting face centers (0 to 1, 2 to 3 and 4 to 5), I 
>> found that this was not the case. The physical and reference cells have a 
>> rotation too included in the mapping.
>>
>> That was the reason why version 1 
>> 
>>  of 
>> collect_periodic_faces() was working and version 2 
>> 
>>  was 
>> not: since the direction parameter is also used to access faces in the 
>> latter version. This would produce an exception because there is a rotation 
>> of axes between physical and reference space and the boundary id wouldn't 
>> match.
>>
>> This has provoked another question. Version 2 of collect_periodic_faces() 
>> says 'This compatibility version of collect_periodic_faces() only works on 
>> grids with cells in standard orientation'. But is that sufficient? In my 
>> case, the mesh, after reading through GridIn, was in standard orientation. 
>> I had verified this by printing face orientation, flip and rotation. But 
>> since the boundary ids were set according to physical axes directions and 
>> these are different from reference axes directions, the output was 
>> unexpected to me.
>>
>> Any comments on whatever I concluded so far would be greatly helpful to 
>> me!
>>
>> On Thursday, May 6, 2021 at 10:41:21 AM UTC+5:30 vachanpo...@gmail.com 
>> wrote:
>>
>>> 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 
>>> #include 
>>>
>>> #include 
>>> #include 
>>> #include 
>>> #include 
>>> #include 
>>>
>>> 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 trian

Re: [deal.II] Gmsh unexpected behaviour - axes directions swapped

[Daniel Arndt]

Vachan,

The documentation for the second version says

"Instead of defining a 'first' and 'second' boundary with the help of two
boundary_ids this function defines a 'left' boundary as all faces with
local face index 2*direction and boundary indicator b_id and, similarly, a
'right' boundary consisting of all face with local face index 2*direction+1
and boundary indicator b_id. Faces with coordinates only differing in the
direction component are identified."

So there is an assumption about the local face index and the direction
parameter that is also used for identifying faces. This implies that a case
like yours where the cell is rotated is not supported by this overload.
Just as you said, it's not enough for the cell to be in standard
orientation.

I'm sorry that this caused so much trouble but I'm glad that the first
version (that I would probably always recommend using anyway) works for you.

Best,
Daniel



Am Fr., 7. Mai 2021 um 02:58 Uhr schrieb vachanpo...@gmail.com <
vachanpotluri1...@gmail.com>:

> I think I have found an explanation.
>
> Since the mesh generated is a cartesian mesh, I was assuming that the
> faces of physical and reference cells have normals in the same direction
> and the mapping between then would be a simple stretching. But when I
> printed out vectors connecting face centers (0 to 1, 2 to 3 and 4 to 5), I
> found that this was not the case. The physical and reference cells have a
> rotation too included in the mapping.
>
> That was the reason why version 1
> 
>  of
> collect_periodic_faces() was working and version 2
> 
>  was
> not: since the direction parameter is also used to access faces in the
> latter version. This would produce an exception because there is a rotation
> of axes between physical and reference space and the boundary id wouldn't
> match.
>
> This has provoked another question. Version 2 of collect_periodic_faces()
> says 'This compatibility version of collect_periodic_faces() only works on
> grids with cells in standard orientation'. But is that sufficient? In my
> case, the mesh, after reading through GridIn, was in standard orientation.
> I had verified this by printing face orientation, flip and rotation. But
> since the boundary ids were set according to physical axes directions and
> these are different from reference axes directions, the output was
> unexpected to me.
>
> Any comments on whatever I concluded so far would be greatly helpful to me!
>
> On Thursday, May 6, 2021 at 10:41:21 AM UTC+5:30 vachanpo...@gmail.com
> wrote:
>
>> 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 
>> #include 
>>
>> #include 
>> #include 
>> #include 
>> #include 
>> #include 
>>
>> 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 triang(mpi_comm);
>>
>> GridIn 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-