]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Improve an algorithm in MappingQ::transform_real_to_unit_cell.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 30 May 2012 21:11:51 +0000 (21:11 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 30 May 2012 21:11:51 +0000 (21:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@25578 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/source/fe/mapping_q.cc
deal.II/source/fe/mapping_q1.cc
tests/fe/mapping_real_to_unit_q4_curved.cc
tests/fe/mapping_real_to_unit_q4_curved_codim.cc [new file with mode: 0644]
tests/fe/mapping_real_to_unit_q4_curved_codim/cmp/generic [new file with mode: 0644]

index be820c9c3b7b1bd38e6761d53b6f5b98bee2d6bb..8a4f303a513d7edce221df3e0fbd55a77bdcc2e2 100644 (file)
@@ -201,13 +201,18 @@ enabled due to a missing include file in file
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Fixed: A problem in MappingQ::transform_real_to_unit_cell
+that sometimes led the algorithm in this function to abort.
+<br>
+(Wolfgang Bangerth 2012/05/30)
+
 <li> New: The function DataOutInterface::write_pvd_record can be used
 to provide Paraview with metadata that describes which time in a
 simulation a particular output file corresponds to.
 <br>
 (Marco Engelhard 2012/05/30)
 
-<li> Fixed: Bug in 3d with hanging nodes in GridTools::find_cells_adjacent_to_vertex()
+<li> Fixed: A bug in 3d with hanging nodes in GridTools::find_cells_adjacent_to_vertex()
 that caused find_active_cell_around_point() to fail in those cases.
 <br>
 (Timo Heister, Wolfgang Bangerth 2012/05/30)
index a6e9222bf2ef5749f38635968d72d644a87c5e29..0beb8330d2042c4148d2466245e9f2329bf00345 100644 (file)
@@ -1445,9 +1445,10 @@ MappingQ<dim,spacedim>::
 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
@@ -1459,6 +1460,19 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
       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);
 
index b3d683b0b54441afa6732db83b4b6b1b22dd2263..87485d368c518939075e8d1e445910dee0890327 100644 (file)
@@ -1933,10 +1933,8 @@ transform_real_to_unit_cell_internal_codim1
                                    // eps. Only if the first condition
                                    // failed, loop will have been
                                    // increased and tested, and thus
-                                   // havereached the limit.
-  AssertThrow(loop<loop_limit, ExcTransformationFailed());
-
-
+                                   // have reached the limit.
+  AssertThrow (loop<loop_limit, ExcTransformationFailed());
 }
 
 
index 3401d6c108638ea389e41c8362ddd21d0664e2d7..9c54286dd2016655bd8665ecd4dc0f3e9519d3d1 100644 (file)
@@ -120,7 +120,6 @@ main()
 
 // the following is not currently implemented
 //  test_real_to_unit_cell<1,3>();
-//  test_real_to_unit_cell<2,3>();
 
   return 0;
 }
diff --git a/tests/fe/mapping_real_to_unit_q4_curved_codim.cc b/tests/fe/mapping_real_to_unit_q4_curved_codim.cc
new file mode 100644 (file)
index 0000000..2349ba7
--- /dev/null
@@ -0,0 +1,121 @@
+//-----------------------------------------------------------------------------
+//    $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;
+}
+
+
+
diff --git a/tests/fe/mapping_real_to_unit_q4_curved_codim/cmp/generic b/tests/fe/mapping_real_to_unit_q4_curved_codim/cmp/generic
new file mode 100644 (file)
index 0000000..4e66a74
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::dim=2, spacedim=3
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.