Dear Jean-Paul, yes, exactly, if the entries were parsed as you mentioned everything would work fine. But unfortunately I guess this is not the case. I post here the simplified problem. Now the mesh described in the .msh file is made of two cells: the one with index 0 which is wrong and one of its neighbours, with index 1, which is ok. The vertexes of cell 0 are the following:
x1 = 5.399999999999241 , y1 = 0.4999999999988218 x2 = 5.499999999999598 , y2 = 0.499999999998822 x4 = 5.399999999999546 , y4 = 0.59999999999853 x5 = 5.499999999999902 , y5 = 0.59999999999853 For this representation the mapping on cell 0 will fail. So, let's just consider the coordinate y1. If we approximate the value to 0.49999999999882*2* the mapping will be ok. However, if we try one of the following, as example: y1 = 0.49999999999882*1* y1 = 0.49999999999882 *21*interesting enough, the mapping will fail again. Another test: I add say two random digits to each vertex coordinate: x1 = 5.399999999999241*21* , y1 = 0.4999999999988218*21* x2 = 5.499999999999598*21* , y2 = 0.499999999998822*21* x4 = 5.399999999999546*21* , y4 = 0.59999999999853*21* x5 = 5.499999999999902*21* , y5 = 0.59999999999853*21* Now it works again...this is mind-blowing. Just something wrong occurs when parsing the values. Hope this is helpful. Please let me know if I can provide some more useful information. Thank you very much again, Best, Giovanni On Thursday, July 20, 2017 at 9:08:34 PM UTC+2, Jean-Paul Pelteret wrote: > > Dear Giovanni, > > Its good to hear that you managed to solve / work around this issue! Would > it be possible for you to post the reduced (problematic) mesh so that we > can see if there's anything that we can learn from it. I'm guessing that > you expected lines with entries such as > 6 0.1999999999991087 0 0 > to be parsed as > 6 0.2 0 0 > ? Did you ultimately have to output a grid that had the coordinates > recorded with a greater precision? > > Thanks, > Jean-Paul > > On Thursday, July 20, 2017 at 7:47:04 PM UTC+2, Giovanni Di Ilio wrote: >> >> Dear Wolfgang, >> >> thanks for your quick reply and for your suggestion. >> I simplified the case to only two cells, the one that was supposed to be >> "wrong" and one of its neighbours, that was supposed to be ok, for >> comparison. >> By watching at the coordinates of the vertices in the file .msh I found >> out that the problem was related to some rounding error. I guess there is a >> bug in the parser. I just converted the coordinate values into double and >> now the mapping works fine everywhere. >> >> Thank you very much again, >> >> Best, >> Giovanni >> >> On Thursday, July 20, 2017 at 4:32:43 PM UTC+2, Wolfgang Bangerth wrote: >>> >>> On 07/20/2017 05:53 AM, Giovanni Di Ilio wrote: >>> > >>> > 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. >>> >>> It's often difficult to debug these cases with such a large mesh. But >>> since >>> you know which cell the problem appears in, take the test.msh file and >>> remove >>> (by hand) all other cells and all of the vertices not needed. This way >>> you >>> should end up with a testcase that has only one cell and that should be >>> much >>> easier to figure out. >>> >>> Best >>> W. >>> >>> -- >>> ------------------------------------------------------------------------ >>> Wolfgang Bangerth email: [email protected] >>> www: >>> http://www.math.colostate.edu/~bangerth/ >>> >>> -- 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]. 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()
{
if(use_hypercube)
{
std::vector< unsigned int > repetitions(dim);
repetitions[0] = 2;
repetitions[1] = 1;
GridGenerator::subdivided_hyper_rectangle(triangulation, repetitions, Point<dim>(5.4,0.5), Point<dim>(5.6,0.6), true);
}
else
{
GridIn<dim> gridin;
gridin.attach_triangulation(triangulation);
std::ifstream f("test.msh");
gridin.read_msh(f);
}
}
void test_mapping()
{
auto cell = dof_handler.begin_active();
auto endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
Point<dim> test;
auto cellN = cell->index();
if (cellN == 0)
{
test(0) = 5.45;
test(1) = 0.55;
}
else if (cellN == 1)
{
test(0) = 5.55;
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 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 gmsh(false);
return 0;
}
test.msh
Description: Mesh model
