transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
const Point<spacedim> &p) const
{
- // first a Newton iteration based
- // on a Q1 mapping
- Point<dim> p_unit = MappingQ1<dim,spacedim>::transform_real_to_unit_cell(cell, p);
+ // first a Newton iteration based on a Q1
+ // mapping.
+ Point<dim> p_unit =
+ MappingQ1<dim,spacedim>::transform_real_to_unit_cell(cell, p);
// then a Newton iteration based on the
// full MappingQ if we need this. note that
use_mapping_q_on_all_cells ||
(dim!=spacedim) )
{
+ // use the full mapping. in case the
+ // function above should have given us
+ // something back that lies outside the
+ // unit cell (that might happen because
+ // we may have given a point 'p' that
+ // lies inside the cell with the higher
+ // order mapping, but outside the
+ // Q1-mapped reference cell), then
+ // project it back into the reference
+ // cell in hopes that this gives a
+ // better starting point to the
+ // following iteration
+ p_unit = GeometryInfo<dim>::project_to_unit_cell(p_unit);
const Quadrature<dim> point_quadrature(p_unit);
--- /dev/null
+//-----------------------------------------------------------------------------
+// $Id: data_out_base_pvd.cc 25569 2012-05-30 12:53:31Z bangerth $
+// Version: $Name$
+//
+// Copyright (C) 2006, 2007, 2010, 2012 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//-----------------------------------------------------------------------------
+
+
+// on a somewhat deformed cube, verify that if we push forward a bunch
+// of points from the reference to the real cell and then call
+// Mapping::transform_unit_to_real_cell that we get the same point as
+// we had in the beginning.
+//
+// like in the _q4_straight test, we use a Q4 mapping but this time we
+// actually curve one boundary of the cell which ensures that the
+// mapping is really higher order than just Q1
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/fe/mapping_q.h>
+
+
+template<int dim, int spacedim>
+void test_real_to_unit_cell()
+{
+ deallog << "dim=" << dim << ", spacedim=" << spacedim << std::endl;
+
+ // define a boundary that fits the
+ // the vertices of the hyper cube
+ // we're going to create below
+ HyperBallBoundary<dim,spacedim> boundary (Point<spacedim>(),
+ std::sqrt(1.*dim));
+
+ Triangulation<dim, spacedim> triangulation;
+ GridGenerator::hyper_cube (triangulation, -1, 1);
+
+ // set the boundary indicator for
+ // one face of the single cell
+ triangulation.set_boundary (1, boundary);
+ triangulation.begin_active()->face(0)->set_boundary_indicator (1);
+
+ const unsigned int n_points = 5;
+ std::vector< Point<dim> > unit_points(Utilities::fixed_power<dim>(n_points));
+
+ switch (dim)
+ {
+ case 1:
+ for (unsigned int x=0; x<n_points;++x)
+ unit_points[x][0] = double(x)/double(n_points);
+ break;
+
+ case 2:
+ for (unsigned int x=0; x<n_points;++x)
+ for (unsigned int y=0; y<n_points;++y)
+ {
+ unit_points[y * n_points + x][0] = double(x)/double(n_points);
+ unit_points[y * n_points + x][1] = double(y)/double(n_points);
+ }
+ break;
+
+ case 3:
+ for (unsigned int x=0; x<n_points;++x)
+ for (unsigned int y=0; y<n_points;++y)
+ for (unsigned int z=0; z<n_points;++z)
+ {
+ unit_points[z * n_points + y * n_points + x][0] = double(x)/double(n_points);
+ unit_points[z * n_points + y * n_points + x][1] = double(y)/double(n_points);
+ unit_points[z * n_points + y * n_points + x][2] = double(z)/double(n_points);
+ }
+ break;
+ }
+
+
+ MappingQ< dim, spacedim > map(4);
+
+ // work with this cell (unlike the
+ // _q1 test where we move vertices)
+ typename Triangulation<dim, spacedim >::active_cell_iterator
+ cell = triangulation.begin_active();
+ for (unsigned int i=0; i<unit_points.size();++i)
+ {
+ // for each of the points,
+ // verify that if we apply
+ // the forward map and then
+ // pull back that we get
+ // the same point again
+ const Point<spacedim> p = map.transform_unit_to_real_cell(cell,unit_points[i]);
+ const Point<dim> p_unit = map.transform_real_to_unit_cell(cell,p);
+ Assert (unit_points[i].distance(p_unit) < 1e-10,
+ ExcInternalError());
+ }
+ deallog << "OK" << std::endl;
+}
+
+
+int
+main()
+{
+ std::ofstream logfile ("mapping_real_to_unit_q4_curved_codim/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+
+ test_real_to_unit_cell<2,3>();
+
+ return 0;
+}
+
+
+
projects / dealii-svn.git / commitdiff