]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move the remaining pieces of transform_r_to_u() from MappingQ to MappingQGeneric.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 1 Oct 2015 15:16:51 +0000 (10:16 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 1 Oct 2015 20:16:34 +0000 (15:16 -0500)
include/deal.II/fe/mapping.h
include/deal.II/fe/mapping_q1_eulerian.h
include/deal.II/fe/mapping_q_eulerian.h
source/fe/mapping_q.cc
source/fe/mapping_q_eulerian.cc
source/fe/mapping_q_generic.cc

index be9671888c1b73f9bf8b6e0dae436bc182b7882d..3a99fbc5422339417bbf935aecb170f8fb16c67b 100644 (file)
@@ -369,7 +369,7 @@ public:
    * equals the given point @p p.  If this is the case then this function
    * throws an exception of type Mapping::ExcTransformationFailed . Whether
    * the given point @p p lies outside the cell can therefore be determined by
-   * checking whether the return reference coordinates lie inside or outside
+   * checking whether the returned reference coordinates lie inside or outside
    * the reference cell (e.g., using GeometryInfo::is_inside_unit_cell()) or
    * whether the exception mentioned above has been thrown.
    *
index bd2e3078a5f1fbf8ce2fd2ab5850b827719437ba..fc17a1734a8d4825789ccc8cb6a4d177c6a1cbe1 100644 (file)
@@ -108,8 +108,7 @@ public:
    */
   virtual
   std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
-  get_vertices
-  (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
+  get_vertices (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
 
   /**
    * Return a pointer to a copy of the present object. The caller of this copy
index c81e62a68821d368d3e929f026eb5e0421aeb297..f4bd91316b074866b2bf18e9401c573a00deb4b0 100644 (file)
@@ -116,6 +116,16 @@ public:
                     const VECTOR                   &euler_vector,
                     const DoFHandler<dim,spacedim> &euler_dof_handler) DEAL_II_DEPRECATED;
 
+  /**
+   * Return the mapped vertices of the cell. For the current class, this function does
+   * not use the support points from the geometry of the current cell but
+   * instead evaluates an externally given displacement field in addition to
+   * the geometry of the cell.
+   */
+  virtual
+  std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+  get_vertices (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
+
   /**
    * Return a pointer to a copy of the present object. The caller of this copy
    * then assumes ownership of it.
index 400d2998d6282f04765b01376776cccf8b64fcaf..e2162e64f6eec691e63a0537619a3bcc8c462e59 100644 (file)
@@ -445,60 +445,14 @@ 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 to get a good starting
-  // point, the idea being that this is cheaper than trying to start with the
-  // real mapping and likely also more robust.
-  //
-  // that said, this doesn't always work: there are cases where the point is
-  // outside the cell and the inverse mapping doesn't converge. in that case,
-  // use the center point of the cell as a starting point if we are to go on
-  // using the full mapping, or just propagate up the exception if we had no
-  // intention of continuing with the full mapping
-  Point<dim> initial_p_unit;
-  try
-    {
-      initial_p_unit = q1_mapping->transform_real_to_unit_cell(cell, p);
-    }
-  catch (const typename Mapping<dim,spacedim>::ExcTransformationFailed &)
-    {
-      // mirror the conditions of the code below to determine if we need to
-      // use an arbitrary starting point or if we just need to rethrow the
-      // exception
-      if (cell->has_boundary_lines()
-          ||
-          use_mapping_q_on_all_cells
-          ||
-          (dim!=spacedim) )
-        {
-          for (unsigned int d=0; d<dim; ++d)
-            initial_p_unit[d] = 0.5;
-        }
-      else
-        throw;
-    }
-
-  // then a Newton iteration based on the full MappingQ if we need this. note
-  // that for interior cells with dim==spacedim, the mapping used is in fact a
-  // Q1 mapping, so there is nothing we need to do unless the iteration above
-  // failed
   if (cell->has_boundary_lines()
       ||
       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
-      initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
-
-      return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
-    }
+    return MappingQGeneric<dim,spacedim>::transform_real_to_unit_cell(cell, p);
   else
-    return initial_p_unit;
+    return q1_mapping->transform_real_to_unit_cell(cell, p);
 }
 
 
index c7a91efb09562e9a11e0b05ce1a317f68eaaf187..51b86544983a447bb303f6a337d4719460b3721f 100644 (file)
@@ -126,13 +126,32 @@ SupportQuadrature (const unsigned int map_degree)
 
 // .... COMPUTE MAPPING SUPPORT POINTS
 
+template <int dim, class EulerVectorType, int spacedim>
+std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+MappingQEulerian<dim, EulerVectorType, spacedim>::
+get_vertices
+(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
+{
+  // get the vertices as the first 2^dim mapping support points
+  std::vector<Point<spacedim> > a;
+  compute_mapping_support_points(cell, a);
+
+  std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertex_locations;
+  std::copy (a.begin(),
+             a.begin()+GeometryInfo<dim>::vertices_per_cell,
+             vertex_locations.begin());
+
+  return vertex_locations;
+}
+
+
+
 template <int dim, class EulerVectorType, int spacedim>
 void
 MappingQEulerian<dim, EulerVectorType, spacedim>::
 compute_mapping_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                                 std::vector<Point<spacedim> > &a) const
 {
-
   // first, basic assertion with respect to vector size,
 
   const types::global_dof_index n_dofs  = euler_dof_handler->n_dofs();
index 034f2b4edb7741126d88bcc12dada25308e4b5aa..946a655788b4de7c72dc98e6b1b70bc5b660bc4b 100644 (file)
@@ -31,6 +31,7 @@
 #include <deal.II/fe/fe.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping_q_generic.h>
+#include <deal.II/fe/mapping_q1.h>
 
 #include <cmath>
 #include <algorithm>
@@ -1859,15 +1860,49 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
     }
 
 
-  // Find the initial value for the Newton iteration by a normal
-  // projection to the least square plane determined by the vertices
-  // of the cell
-  std::vector<Point<spacedim> > a;
-  compute_mapping_support_points (cell,a);
-  Assert(a.size() == GeometryInfo<dim>::vertices_per_cell,
-         ExcInternalError());
-  const Point<dim> initial_p_unit =
-    internal::MappingQ1::transform_real_to_unit_cell_initial_guess<dim,spacedim>(a,p);
+  Point<dim> initial_p_unit;
+  if (polynomial_degree == 1)
+    {
+      // Find the initial value for the Newton iteration by a normal
+      // projection to the least square plane determined by the vertices
+      // of the cell
+      std::vector<Point<spacedim> > a;
+      compute_mapping_support_points (cell,a);
+      Assert(a.size() == GeometryInfo<dim>::vertices_per_cell,
+             ExcInternalError());
+      initial_p_unit = internal::MappingQ1::transform_real_to_unit_cell_initial_guess<dim,spacedim>(a,p);
+    }
+  else
+    {
+      try
+        {
+          // Find the initial value for the Newton iteration by a normal
+          // projection to the least square plane determined by the vertices
+          // of the cell
+          //
+          // we do this by first getting all support points, then
+          // throwing away all but the vertices, and finally calling
+          // the same function as above
+          std::vector<Point<spacedim> > a;
+          compute_mapping_support_points (cell,a);
+          a.resize(GeometryInfo<dim>::vertices_per_cell);
+          initial_p_unit = internal::MappingQ1::transform_real_to_unit_cell_initial_guess<dim,spacedim>(a,p);
+        }
+      catch (const typename Mapping<dim,spacedim>::ExcTransformationFailed &)
+        {
+          for (unsigned int d=0; d<dim; ++d)
+            initial_p_unit[d] = 0.5;
+        }
+
+      // 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
+      initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
+    }
 
   // perform the Newton iteration and return the result. note that
   // this statement may throw an exception, which we simply pass up to

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.