]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fixing Jacobians.
authorLuca Heltai <luca.heltai@sissa.it>
Thu, 3 Mar 2016 18:52:56 +0000 (19:52 +0100)
committerLuca Heltai <luca.heltai@sissa.it>
Wed, 6 Apr 2016 11:00:14 +0000 (13:00 +0200)
include/deal.II/fe/mapping_manifold.h
source/fe/mapping_manifold.cc

index 52371c38e0bed706d181f01070393462e32a24d9..4a16e5866f5beb0fea696f09a4de75ed4c0afeb8 100644 (file)
@@ -205,17 +205,30 @@ public:
     void
     compute_manifold_quadrature_weights(const Quadrature<dim> &quadrature);
 
+    /**
+     * Store vertices internally.
+     */
+    void
+    store_vertices(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
+
     /**
      * Return an estimate (in bytes) or the memory consumption of this object.
      */
     virtual std::size_t memory_consumption () const;
 
     /**
-     * Store the current cell.
+     * The current cell vertices.
+     *
+     * Computed each.
+     */
+    mutable std::vector<Point<spacedim> > vertices;
+
+    /**
+     * The current cell.
      *
      * Computed each.
      */
-    typename Triangulation<dim,spacedim>::cell_iterator current_cell;
+    mutable typename Triangulation<dim,spacedim>::cell_iterator cell;
 
     /**
      * The actual quadrature on the reference cell.
@@ -331,11 +344,16 @@ public:
     //  */
     // mutable typename Triangulation<dim,spacedim>::cell_iterator cell_of_current_support_points;
 
-    // /**
-    //  * The determinant of the Jacobian in each quadrature point. Filled if
-    //  * #update_volume_elements.
-    //  */
-    // mutable std::vector<double> volume_elements;
+    /**
+     * The determinant of the Jacobian in each quadrature point. Filled if
+     * #update_volume_elements.
+     */
+    mutable std::vector<double> volume_elements;
+
+    /**
+     * A Q1 Finite element, to compute weights.
+     */
+    const FE_Q<dim> fe_q;
   };
 
 
@@ -399,63 +417,6 @@ protected:
    */
   const FE_Q<dim,spacedim> fe_q;
 
-  /**
-   * A table of weights by which we multiply the locations of the support
-   * points on the perimeter of a quad to get the location of interior support
-   * points.
-   *
-   * Sizes: support_point_weights_on_quad.size()= number of inner
-   * unit_support_points support_point_weights_on_quad[i].size()= number of
-   * outer unit_support_points, i.e.  unit_support_points on the boundary of
-   * the quad
-   *
-   * For the definition of this vector see equation (8) of the `mapping'
-   * report.
-   */
-  Table<2,double> support_point_weights_on_quad;
-
-  /**
-   * A table of weights by which we multiply the locations of the support
-   * points on the perimeter of a hex to get the location of interior support
-   * points.
-   *
-   * For the definition of this vector see equation (8) of the `mapping'
-   * report.
-   */
-  Table<2,double> support_point_weights_on_hex;
-
-  /**
-   * Return the locations of support points for the mapping. For example, for
-   * $Q_1$ mappings these are the vertices, and for higher order polynomial
-   * mappings they are the vertices plus interior points on edges, faces, and
-   * the cell interior that are placed in consultation with the Manifold
-   * description of the domain and its boundary. However, other classes may
-   * override this function differently. In particular, the MappingQ1Eulerian
-   * class does exactly this by not computing the support points from the
-   * geometry of the current cell but instead evaluating an externally given
-   * displacement field in addition to the geometry of the cell.
-   *
-   * The default implementation of this function is appropriate for most
-   * cases. It takes the locations of support points on the boundary of the
-   * cell from the underlying manifold. Interior support points (ie. support
-   * points in quads for 2d, in hexes for 3d) are then computed using the
-   * solution of a Laplace equation with the position of the outer support
-   * points as boundary values, in order to make the transformation as smooth
-   * as possible.
-   *
-   * The function works its way from the vertices (which it takes from the
-   * given cell) via the support points on the line (for which it calls the
-   * add_line_support_points() function) and the support points on the quad
-   * faces (in 3d, for which it calls the add_quad_support_points() function).
-   * It then adds interior support points that are either computed by
-   * interpolation from the surrounding points using weights computed by
-   * solving a Laplace equation, or if dim<spacedim, it asks the underlying
-   * manifold for the locations of interior points.
-   */
-  // virtual
-  // std::vector<Point<spacedim> >
-  // compute_mapping_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
-
   /**
    * Transforms the point @p p on the real cell to the corresponding point on
    * the unit cell @p cell by a Newton iteration.
@@ -464,13 +425,6 @@ protected:
   // transform_real_to_unit_cell_internal (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   //                                       const Point<spacedim> &p,
   //                                       const Point<dim> &initial_p_unit) const;
-
-  /**
-   * Make MappingQ a friend since it needs to call the fill_fe_values()
-   * functions on its MappingManifold(1) sub-object.
-   */
-  template <int, int> friend class MappingQ;
-
 };
 
 
@@ -481,17 +435,16 @@ protected:
 
 #ifndef DOXYGEN
 
-// template<int dim, int spacedim>
-// inline
-// const Tensor<1,dim> &
-// MappingManifold<dim,spacedim>::InternalData::derivative (const unsigned int qpoint,
-//                                                          const unsigned int shape_nr) const
-// {
-//   Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(),
-//          ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
-//                        shape_derivatives.size()));
-//   return shape_derivatives [qpoint*n_shape_functions + shape_nr];
-// }
+template<int dim, int spacedim>
+inline
+void
+MappingManifold<dim,spacedim>::InternalData::store_vertices (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
+{
+  vertices.resize(GeometryInfo<dim>::vertices_per_cell);
+  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+    vertices[i] = cell->vertex(i);
+  this->cell = cell;
+}
 
 
 
index d6216c6779763b95bee61d5a1a3eb397410ea8b5..dfa7449188cd718a8616c90466dad7672f5c5015 100644 (file)
@@ -41,7 +41,8 @@
 DEAL_II_NAMESPACE_OPEN
 
 template<int dim, int spacedim>
-MappingManifold<dim,spacedim>::InternalData::InternalData ()
+MappingManifold<dim,spacedim>::InternalData::InternalData () :
+  fe_q(1)
 {}
 
 
@@ -81,7 +82,7 @@ initialize (const UpdateFlags      update_flags,
 
   // see if we need the (transformation) shape function values
   // and/or gradients and resize the necessary arrays
-  if (this->update_each & update_quadrature_points)
+  if (this->update_each & (update_quadrature_points | update_contravariant_transformation) )
     compute_manifold_quadrature_weights(q);
 
   if (this->update_each & update_covariant_transformation)
@@ -178,8 +179,7 @@ initialize_face (const UpdateFlags      update_flags,
 
 
 template<int dim, int spacedim>
-MappingManifold<dim,spacedim>::MappingManifold ()
-  :
+MappingManifold<dim,spacedim>::MappingManifold () :
   fe_q(1)
   // support_point_weights_on_quad (compute_support_point_weights_on_quad<dim>(this->polynomial_degree)),
   // support_point_weights_on_hex (compute_support_point_weights_on_hex<dim>(this->polynomial_degree)),
@@ -189,8 +189,7 @@ MappingManifold<dim,spacedim>::MappingManifold ()
 
 
 template<int dim, int spacedim>
-MappingManifold<dim,spacedim>::MappingManifold (const MappingManifold<dim,spacedim> &mapping)
-  :
+MappingManifold<dim,spacedim>::MappingManifold (const MappingManifold<dim,spacedim> &mapping) :
   fe_q(1)
 {}
 
@@ -359,28 +358,23 @@ namespace internal
      */
     template <int dim, int spacedim>
     void
-    maybe_compute_q_points (const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell,
-                            const typename QProjector<dim>::DataSetDescriptor                 data_set,
+    maybe_compute_q_points (const typename QProjector<dim>::DataSetDescriptor                 data_set,
                             const typename dealii::MappingManifold<dim,spacedim>::InternalData      &data,
                             std::vector<Point<spacedim> >                                     &quadrature_points)
     {
       const UpdateFlags update_flags = data.update_each;
 
+      AssertDimension(data.vertices.size(), GeometryInfo<dim>::vertices_per_cell);
 
       if (update_flags & update_quadrature_points)
         {
-
           AssertDimension(quadrature_points.size(),
                           data.cell_manifold_quadrature_weights.size());
 
-          std::vector<Point<spacedim> > vertices;
-          for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-            vertices.push_back(cell->vertex(v));
-
           for (unsigned int point=0; point<quadrature_points.size(); ++point)
             {
-              quadrature_points[point] = cell->get_manifold().
-                                         get_new_point(Quadrature<spacedim>(vertices,
+              quadrature_points[point] = data.cell->get_manifold().
+                                         get_new_point(Quadrature<spacedim>(data.vertices,
                                                        data.cell_manifold_quadrature_weights[point]));
             }
         }
@@ -394,8 +388,7 @@ namespace internal
      */
     template <int dim, int spacedim>
     void
-    maybe_update_Jacobians (const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell,
-                            const typename dealii::QProjector<dim>::DataSetDescriptor          data_set,
+    maybe_update_Jacobians (const typename dealii::QProjector<dim>::DataSetDescriptor          data_set,
                             const typename dealii::MappingManifold<dim,spacedim>::InternalData      &data)
     {
       const UpdateFlags update_flags = data.update_each;
@@ -407,13 +400,22 @@ namespace internal
           std::fill(data.contravariant.begin(), data.contravariant.end(),
                     DerivativeForm<1,dim,spacedim>());
 
+          // Cache of weights used to compute points on the reference cell
+          std::vector<double> weights(GeometryInfo<dim>::vertices_per_cell);
 
+          AssertDimension(GeometryInfo<dim>::vertices_per_cell,
+                          data.vertices.size());
           for (unsigned int point=0; point<n_q_points; ++point)
             {
               // Start by figuring out how to compute the direction in
               // the reference space:
               const Point<dim> &p = data.quad.point(point+data_set);
 
+              // And get its image on the manifold:
+              const Point<spacedim> P = data.cell->get_manifold().
+                                        get_new_point(Quadrature<spacedim>(data.vertices,
+                                                      data.cell_manifold_quadrature_weights[point+data_set]));
+
               // Always get the maximum length from the point to the
               // boundary of the reference element, to compute the
               // tangent vectors from the Manifold object
@@ -429,7 +431,17 @@ namespace internal
                   // which is positive, if the coordinate is < .5,
                   double L = ai > .5 ? -ai: 1-ai;
 
-                  data.contravariant[point][i] = cell->get_manifold().get_tangent_vector(p, np)/L;
+                  // Get the weights to compute the np point in real space
+                  for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
+                    weights[j] = data.fe_q.shape_value(j, np);
+
+                  Point<spacedim> NP=data.cell->get_manifold().
+                                     get_new_point(Quadrature<spacedim>(data.vertices, weights));
+
+                  Tensor<1,spacedim> T = data.cell->get_manifold().get_tangent_vector(P, NP);
+
+                  for (unsigned int d=0; d<spacedim; ++d)
+                    data.contravariant[point][i][d] = T[d]/L;
                 }
             }
 
@@ -900,14 +912,14 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 
   const unsigned int n_q_points=quadrature.size();
 
-  internal::maybe_compute_q_points<dim,spacedim> (cell,
-                                                  QProjector<dim>::DataSetDescriptor::cell (),
+  data.store_vertices(cell);
+
+  internal::maybe_compute_q_points<dim,spacedim> (QProjector<dim>::DataSetDescriptor::cell (),
                                                   data,
                                                   output_data.quadrature_points);
 
-//   internal::maybe_update_Jacobians<dim,spacedim> (cell_similarity,
-//                                                   QProjector<dim>::DataSetDescriptor::cell (),
-//                                                   data);
+  internal::maybe_update_Jacobians<dim,spacedim> (QProjector<dim>::DataSetDescriptor::cell (),
+                                                  data);
 
 //   const UpdateFlags update_flags = data.update_each;
 //   const std::vector<double> &weights=quadrature.get_weights();

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.