Sorry, I thought I worked around the problem but actually I did not. The
mapping fails again when I work with a refined mesh, for both the cases of
a mesh read in from gmsh and a mesh generated with Deal.
I try to explain better what I mean:
1) I generate a grid, using:
*GridGenerator::subdivided_hyper_rectangle(triangulation,repetitions,p1,p2)*
;
2) I refine one time the grid with: *triangulation.refine_global(1)* ;
3) I run over the cells and I try to get the mapped position of a given
point, which I set to be in the centre of the considered cell:
*mapping.transform_real_to_unit_cell*.
Since the point I am testing is exactly in the centre of the cell I expect
the mapped position, in the unit cell, to be x = 0.5; y = 0.5.
However, for some "random" cell, again, the mapping fails.
Here I am not dealing anymore with the digits of accuracy.
The numbers I am dealing with are the coordinates of the vertexes provided
by the GridGenerator and the coordinates of a testing point, which have
just a couple of digits after comma. Moreover, the spacing between vertexes
for the cells I am considering is of the order of 0.1.
I attach a very simple example, where I consider a uniform grid of 5 x 5
elements (repetitions[0]=repetitions[1]=5), refined one time, so 100
elements in total.
For cell with index 96 of my example, which vertexes are:
vertex(0): 0.8 0.8
vertex(1): 0.9 0.8
vertex(2): 0.8 0.9
vertex(3): 0.9 0.9
I test the point (0.85, 0.85) and the mapped position in the unit cell
results to be: (0.5, 0.46875) instead of (0.5, 0.5).
You may notice that, by running the example with refineglobal(false), the
mapping on an equivalent grid made of 10 x 10 elements, with no refinements
(repetitions[0]=repetitions[1]=10), works just fine.
I don't see any difference between the two cases.
What is missing here? Can you help me please?
Thank you very much,
Best,
Giovanni
On Saturday, July 22, 2017 at 8:57:57 PM UTC+2, Wolfgang Bangerth wrote:
>
> On 07/22/2017 04:42 AM, Giovanni Di Ilio wrote:
> >
> > Yes, I agree, the issue is definitely related to round-off error. But I
> still
> > don't understand why the mapping fails, if the vertex coordinates read
> in from
> > my .msh file are apparently stored as doubles...
> > Anyway, yes, the problem is easily overcome by "manipulating" the values
> > before mapping.
>
> The numbers you consider differ in the last two or three bits when stored
> as
> double precision. Any computations you do with them, for example in the
> mapping, will then necessarily be completely off.
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> <javascript:>
> 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 RefineTest
{
public:
static constexpr int dim = 2 ;
public:
RefineTest(bool _refine)
:
triangulation(),
dof_handler(triangulation),
refine(_refine)
{
make_grid();
test_mapping();
}
private:
void make_grid()
{
if(refine)
{
std::vector< unsigned int > repetitions(dim);
repetitions[0] = repetitions[1] = 5;
std::array<double,dim> L = {1.0, 1.0};
GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, Point<dim>(0.0,0.0), Point<dim>(L[0],L[1]), true);
triangulation.refine_global(1);
}
else
{
std::vector< unsigned int > repetitions(dim);
repetitions[0] = repetitions[1] = 10;
std::array<double,dim> L = {1.0, 1.0};
GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, Point<dim>(0.0,0.0), Point<dim>(L[0],L[1]), true);
triangulation.refine_global(0);
}
}
void test_mapping()
{
auto cell = dof_handler.begin_active();
auto endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
auto cellN = cell->index();
auto P1 = cell->vertex(0);
auto P2 = cell->vertex(1);
auto P3 = cell->vertex(2);
Point<dim> test;
test(0) = (P1[0]+P2[0])/2.0;
test(1) = (P2[1]+P3[1])/2.0;
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);
if (abs(test[0]-dp_real[0]) > 1e-8 || abs(test[1]-dp_real[1]) > 1e-8)
{
std::cout << " " << std::endl;
std::cout << "ERROR" << std::endl;
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 refine;
};
int main(int argc, char *argv[])
{
RefineTest refineglobal(true);
return 0;
}
