[deal.II] Mapping points of mesh generated by gmsh

[Giovanni Di Ilio]

Dear Deal.II community,

I am facing a problem with the mapping of a point on the real cell to the 
unit cell when the mesh is read in from a .msh file.

So what I have is a .geo file (which contains "physical lines" and 
"physical surfaces", as indicated in Step49). I generate the .msh file and 
I read in using the GridIn class. The result is a uniform structured mesh.

Then, since I want to compute the shape functions, I need to know the 
position of the mapped points on the unit cell. What happens is that for 
some cell (only few cells) the function *transform_real_to_unit_cell* gives 
me a wrong result, that is the coordinates of some mapped point are not 
correct.
If I generate the same identical mesh using *subdivided_hyper_rectangle* 
instead, everything works well.

I checked the .msh file and everything looks fine for me. What concerns me 
the most is that this occurs only for certain "random" cells of the domain.
I don't see any apparent explanation to this. For sure I am missing 
something but I was not able to figure out the solution to this issue.

I provide a small code that I am testing as well as the .msh and .geo files 
I am using.
The code computes the mapped point of a generic real point which is placed 
in the center of a "sick" cell (id=354). I would expect that the 
coordinates of the unit point are x=0.5, y=0.5. However, what I get is 
x=0.5, y=0.0. The same apply if you pick an other real point within this 
cell.

May you please help me? I would really appreciate it.
Thank you very much in advance,

Best,
Giovanni


-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/mapping_q.h>

#include <fstream>
#include <iostream>

using namespace dealii;

class GmeshTest
{

public:
    static constexpr int dim  = 2 ;

public:

    GmeshTest(bool _use_hypercube)
    :
        triangulation(),
        dof_handler(triangulation),
        use_hypercube(_use_hypercube)
    {
        make_grid(); 
        test_mapping();
    }

private:
    void make_grid()
    {

        std::string outmesh_name;
        if(use_hypercube)
        {
            std::vector< unsigned  int > repetitions(dim);
            repetitions[0] = 60;
            repetitions[1] = 20;

            std::array<double,dim> L = {6.0, 2.0};
            GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, Point<dim>(0.0,0.0), Point<dim>(L[0],L[1]), true);
            outmesh_name="mesh_hypercube.vtk";
        }
        else
        {
            GridIn<dim> gridin;
            gridin.attach_triangulation(triangulation);
            std::ifstream f("test.msh"); 
            gridin.read_msh(f);
            outmesh_name="mesh_gmsh.vtk";
        }

        std::ofstream out(outmesh_name);
        GridOut grid_out;
        grid_out.write_vtk(triangulation, out);
    
    }

    void test_mapping()
    {
        auto cell = dof_handler.begin_active();
        auto endc = dof_handler.end();
        for (; cell!=endc; ++cell)
        {
            auto cellN = cell->index();
            if (cellN == 354)
            {
                Point<dim> test;
                test(0) = 5.45;
                test(1) = 0.55;

                MappingQ<dim> mapping(1);

                Point<dim> dp = mapping.transform_real_to_unit_cell(cell,test);
                Point<dim> dp_real = mapping. transform_unit_to_real_cell(cell,dp);

                std::cout << "cell=" << cell->index() << std::endl;
                std::cout << "cell vertex(0): " << cell->vertex(0) << std::endl;
                std::cout << "cell vertex(1): " << cell->vertex(1) << std::endl;
                std::cout << "cell vertex(2): " << cell->vertex(2) << std::endl;
                std::cout << "cell vertex(3): " << cell->vertex(3) << std::endl;
                std::cout << "=================================================" << std::endl;
                std::cout << "test_point=   " << test << "  mapped_point=   " << dp <<" backmapping: "<<dp_real<< std::endl;
                std::cout << "=================================================" << std::endl;
            }
        }
    }

private:
    Triangulation<dim>      triangulation;
    DoFHandler<dim>         dof_handler;
    const bool              use_hypercube;
};

int main(int argc, char *argv[])
{

    GmeshTest hypercube(true);
    GmeshTest gmsh(false);

    return 0;
}

test.geo
Description: Binary data

test.msh
Description: Mesh model