]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add and clarify some comments. Add a second set of transform_{co,contra}variant funct...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 3 Sep 2002 15:30:31 +0000 (15:30 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 3 Sep 2002 15:30:31 +0000 (15:30 +0000)
git-svn-id: https://svn.dealii.org/trunk@6349 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/mapping.h
deal.II/deal.II/include/fe/mapping_cartesian.h
deal.II/deal.II/include/fe/mapping_q.h
deal.II/deal.II/include/fe/mapping_q1.h
deal.II/deal.II/include/fe/mapping_q1_eulerian.h
deal.II/deal.II/source/fe/mapping_cartesian.cc
deal.II/deal.II/source/fe/mapping_q.cc
deal.II/deal.II/source/fe/mapping_q1.cc

index e0213bc5f3567a0151e83563a0e8f58a94ce232e..52618daa7fc88b0abf2e8f0c4393ab4364493554 100644 (file)
@@ -18,7 +18,6 @@
 #include <cmath>
 #include <base/point.h>
 #include <base/subscriptor.h>
-#include <base/vector2d.h>
 #include <grid/tria.h>
 #include <dofs/dof_handler.h>
 #include <fe/fe_update_flags.h>
@@ -195,36 +194,97 @@ class Mapping : public Subscriptor
     };
     
                                     /**
-                                     * Tranform a field of covariant
-                                     * vectors. The vector spanned by
-                                     * @p{begin} and @p{end} must
+                                     * Transform a field of covariant
+                                     * vectors. The covariant
+                                     * transformation multiplies a
+                                     * vector from the right with the
+                                     * inverse of the Jacobian matrix
+                                     * of the transformation from
+                                     * unit to real space
+                                     * cell. Alternatively, this is
+                                     * equivalent to a multiplication
+                                     * from the left with the
+                                     * transpose if the inverse
+                                     * matrix.
+                                     *
+                                     * The range of vectors spanned
+                                     * by @p{begin} and @p{end} must
                                      * have as many elements as there
                                      * are quadrature points (not
                                      * tested inside the function).
+                                     * Also note the different
+                                     * convention for parameters
+                                     * compared to the standard C++
+                                     * @p{transform} function: here,
+                                     * first destination, then source
+                                     * are specified, and the number
+                                     * of elements is determined by a
+                                     * range of destination
+                                     * vectors. This convention is
+                                     * due to the way the function is
+                                     * usually called.
                                      *
                                      * The vector @p{src} must
                                      * contain at least as many
-                                     * elements.
+                                     * elements as there are
+                                     * quadrature points.
                                      */
     virtual
     void
     transform_covariant (Tensor<1,dim>          *begin,
                         Tensor<1,dim>          *end,
                         const Tensor<1,dim>    *src,
+                        const InternalDataBase &internal) const = 0;
+
+                                     /**
+                                      * Transform a set of matrices
+                                      * covariantly, i.e. multiply
+                                      * each function in the input
+                                      * range by the Jacobian matrices
+                                      * at the different quadrature
+                                      * points from the left.
+                                      *
+                                      * The same applies as to the
+                                      * function above regarding input
+                                      * and output ranges.
+                                      */
+    virtual
+    void
+    transform_covariant (Tensor<2,dim>          *begin,
+                        Tensor<2,dim>          *end,
+                        const Tensor<2,dim>    *src,
                         const InternalDataBase &internal) const = 0;
     
                                     /**
-                                     * Tranform a field of
+                                     * Transform a field of
                                      * contravariant vectors. The
-                                     * vector spanned by @p{begin}
-                                     * and @p{end} must have as many
-                                     * elements as there are
-                                     * quadrature points (not tested
-                                     * inside the function).
+                                     * contravariant transformation
+                                     * multiplies a vector from the
+                                     * left with the Jacobian matrix
+                                     * of the transformation from
+                                     * unit to real space cell.
+                                     * 
+                                     * The range of vectors spanned
+                                     * by @p{begin} and @p{end} must
+                                     * have as many elements as there
+                                     * are quadrature points (not
+                                     * tested inside the function).
+                                     * Also note the different
+                                     * convention for parameters
+                                     * compared to the standard C++
+                                     * @p{transform} function: here,
+                                     * first destination, then source
+                                     * are specified, and the number
+                                     * of elements is determined by a
+                                     * range of destination
+                                     * vectors. This convention is
+                                     * due to the way the function is
+                                     * usually called.
                                      *
                                      * The vector @p{src} must
                                      * contain at least as many
-                                     * elements.
+                                     * elements as there are
+                                     * quadrature points.
                                      */
     virtual
     void
@@ -233,6 +293,32 @@ class Mapping : public Subscriptor
                             const Tensor<1,dim>    *src,
                             const InternalDataBase &internal) const = 0;
 
+                                     /**
+                                      * Transform a set of matrices
+                                      * contravariantly, i.e. multiply
+                                      * each function in the input
+                                      * range by the inverse Jacobian
+                                      * matrices at the different
+                                      * quadrature points from the
+                                      * right. Note that here it is no
+                                      * more equivalent to
+                                      * multiplication with the
+                                      * transpose of the inverse
+                                      * matrices from the left, unless
+                                      * the matrices to be multiplied
+                                      * with are symmetric.
+                                      *
+                                      * The same applies as to the
+                                      * function above regarding input
+                                      * and output ranges.
+                                      */
+    virtual
+    void
+    transform_contravariant (Tensor<2,dim>          *begin,
+                             Tensor<2,dim>          *end,
+                             const Tensor<2,dim>    *src,
+                             const InternalDataBase &internal) const = 0;
+    
                                     /**
                                      * Indicate fields to be updated
                                      * in the constructor of
@@ -352,7 +438,7 @@ class Mapping : public Subscriptor
     fill_fe_values (const typename DoFHandler<dim>::cell_iterator &cell,
                    const Quadrature<dim>                         &quadrature,
                    InternalDataBase                              &internal,
-                   std::vector<Point<dim> >             &quadrature_points,
+                   std::vector<Point<dim> >                      &quadrature_points,
                    std::vector<double>                           &JxW_values) const = 0;
 
                                     /**
index 0027696eb5939fe4cc6558ac0b139e2193d081da..134943693bcecd00d70ab1f11e81aa2cb804ac50 100644 (file)
@@ -108,6 +108,16 @@ class MappingCartesian : public Mapping<dim>
     transform_covariant (Tensor<1,dim>          *begin,
                         Tensor<1,dim>          *end,
                         const Tensor<1,dim>    *src,
+                        const typename Mapping<dim>::InternalDataBase &internal) const;
+
+                                    /**
+                                     * Implementation of the interface in
+                                     * @ref{Mapping}.
+                                     */
+    virtual void
+    transform_covariant (Tensor<2,dim>          *begin,
+                        Tensor<2,dim>          *end,
+                        const Tensor<2,dim>    *src,
                         const typename Mapping<dim>::InternalDataBase &internal) const;
     
                                     /**
@@ -118,6 +128,16 @@ class MappingCartesian : public Mapping<dim>
     transform_contravariant (Tensor<1,dim>          *begin,
                             Tensor<1,dim>          *end,
                             const Tensor<1,dim>    *src,
+                            const typename Mapping<dim>::InternalDataBase &internal) const;
+    
+                                    /**
+                                     * Implementation of the interface in
+                                     * @ref{Mapping}.
+                                     */
+    virtual void
+    transform_contravariant (Tensor<2,dim>          *begin,
+                            Tensor<2,dim>          *end,
+                            const Tensor<2,dim>    *src,
                             const typename Mapping<dim>::InternalDataBase &internal) const;
 
                                     /**
@@ -205,12 +225,12 @@ class MappingCartesian : public Mapping<dim>
                                          *
                                          * Filled once.
                                          */
-        vector2d<Tensor<1,dim> > unit_tangentials;
+        Table<2,Tensor<1,dim> > unit_tangentials;
        
                                         /**
                                          * Auxiliary vectors for internal use.
                                          */
-        vector2d<Tensor<1,dim> > aux;
+        Table<2,Tensor<1,dim> > aux;
     };
     
                                     /**
index a935afb8060d16814b002aca2105d32d8fcfb568..dbc26e4a1eca0050b786b67b3e10975d57b3f5c4 100644 (file)
@@ -15,7 +15,7 @@
 
 
 #include <base/config.h>
-#include <base/vector2d.h>
+#include <base/table.h>
 #include <fe/mapping_q1.h>
 
 template <int dim> class TensorProductPolynomials;
@@ -79,8 +79,8 @@ class MappingQ : public MappingQ1<dim>
       const Point<dim> &p) const;
 
                                     /**
-                                     * Implementation of the interface in
-                                     * @ref{Mapping}.
+                                     * Implementation of the
+                                     * interface in @ref{Mapping}.
                                      */
     virtual void
     transform_covariant (Tensor<1,dim>          *begin,
@@ -89,8 +89,18 @@ class MappingQ : public MappingQ1<dim>
                         const typename Mapping<dim>::InternalDataBase &internal) const;
     
                                     /**
-                                     * Implementation of the interface in
-                                     * @ref{Mapping}.
+                                     * Implementation of the
+                                     * interface in @ref{Mapping}.
+                                     */
+    virtual void
+    transform_covariant (Tensor<2,dim>          *begin,
+                        Tensor<2,dim>          *end,
+                        const Tensor<2,dim>    *src,
+                        const typename Mapping<dim>::InternalDataBase &internal) const;
+    
+                                    /**
+                                     * Implementation of the
+                                     * interface in @ref{Mapping}.
                                      */
     virtual void
     transform_contravariant (Tensor<1,dim>          *begin,
@@ -98,6 +108,16 @@ class MappingQ : public MappingQ1<dim>
                             const Tensor<1,dim>    *src,
                             const typename Mapping<dim>::InternalDataBase &internal) const;    
 
+                                    /**
+                                     * Implementation of the
+                                     * interface in @ref{Mapping}.
+                                     */
+    virtual void
+    transform_contravariant (Tensor<2,dim>          *begin,
+                            Tensor<2,dim>          *end,
+                            const Tensor<2,dim>    *src,
+                            const typename Mapping<dim>::InternalDataBase &internal) const;    
+    
                                     /**
                                      * Return the degree of the
                                      * mapping, i.e. the value which
@@ -304,7 +324,7 @@ class MappingQ : public MappingQ1<dim>
                                      * this vector is computed.
                                      */
     void
-    set_laplace_on_quad_vector(vector2d<double> &loqvs) const;
+    set_laplace_on_quad_vector(Table<2,double> &loqvs) const;
     
                                     /**
                                      * This function is needed by the
@@ -317,7 +337,7 @@ class MappingQ : public MappingQ1<dim>
                                      * @p{degree>2} this vector is
                                      * computed.
                                      */
-    void set_laplace_on_hex_vector(vector2d<double> &lohvs) const;
+    void set_laplace_on_hex_vector(Table<2,double> &lohvs) const;
     
                                     /**
                                      * Computes the @p{laplace_on_quad(hex)_vector}.
@@ -327,7 +347,7 @@ class MappingQ : public MappingQ1<dim>
                                      * functions if the data is not
                                      * yet hardcoded.
                                      */
-    void compute_laplace_vector(vector2d<double> &lvs) const;
+    void compute_laplace_vector(Table<2,double> &lvs) const;
 
                                     /**
                                      * Takes a
@@ -344,7 +364,7 @@ class MappingQ : public MappingQ1<dim>
                                      * @p{n_inner} computed inner
                                      * points are appended.
                                      */
-    void apply_laplace_vector(const vector2d<double>   &lvs,
+    void apply_laplace_vector(const Table<2,double>   &lvs,
                              std::vector<Point<dim> > &a) const;
     
                                     /**
@@ -426,7 +446,7 @@ class MappingQ : public MappingQ1<dim>
                                      *   unit_support_points on the
                                      *   boundary of the quad
                                      */
-    vector2d<double> laplace_on_quad_vector;
+    Table<2,double> laplace_on_quad_vector;
     
                                     /**
                                      * Needed by the
@@ -434,7 +454,7 @@ class MappingQ : public MappingQ1<dim>
                                      * (for @p{dim==3}). Filled by the
                                      * constructor.
                                      */
-    vector2d<double> laplace_on_hex_vector;
+    Table<2,double> laplace_on_hex_vector;
 
                                     /**
                                      * Exception.
@@ -514,11 +534,11 @@ template<> void MappingQ<1>::compute_shapes_virtual (
   const std::vector<Point<1> > &unit_points,
   MappingQ1<1>::InternalData   &data) const;
 template <> void MappingQ<1>::set_laplace_on_quad_vector(
-  vector2d<double> &) const;
+  Table<2,double> &) const;
 template <> void MappingQ<3>::set_laplace_on_hex_vector(
-  vector2d<double> &lohvs) const;
+  Table<2,double> &lohvs) const;
 template <> void MappingQ<1>::compute_laplace_vector(
-  vector2d<double> &) const;
+  Table<2,double> &) const;
 template <> void MappingQ<1>::add_line_support_points (
   const Triangulation<1>::cell_iterator &,
   std::vector<Point<1> > &) const;
index 24237bb452cc591612233045ad1e5393b0eb1fb4..6aa96558432ac2bc8216def46516d2d98c7f2859 100644 (file)
@@ -15,6 +15,7 @@
 
 
 #include <base/config.h>
+#include <base/table.h>
 #include <cmath>
 #include <fe/mapping.h>
 
@@ -82,9 +83,29 @@ class MappingQ1 : public Mapping<dim>
                                      * @ref{Mapping}.
                                      */
     virtual void
+    transform_covariant (Tensor<2,dim>          *begin,
+                        Tensor<2,dim>          *end,
+                        const Tensor<2,dim>    *src,
+                        const typename Mapping<dim>::InternalDataBase &internal) const;
+    
+                                    /**
+                                     * Implementation of the interface in
+                                     * @ref{Mapping}.
+                                     */
+    virtual void
     transform_contravariant (Tensor<1,dim>          *begin,
                             Tensor<1,dim>          *end,
                             const Tensor<1,dim>    *src,
+                            const typename Mapping<dim>::InternalDataBase &internal) const;
+    
+                                    /**
+                                     * Implementation of the interface in
+                                     * @ref{Mapping}.
+                                     */
+    virtual void
+    transform_contravariant (Tensor<2,dim>          *begin,
+                            Tensor<2,dim>          *end,
+                            const Tensor<2,dim>    *src,
                             const typename Mapping<dim>::InternalDataBase &internal) const;
     
                                     /**
@@ -237,7 +258,14 @@ class MappingQ1 : public Mapping<dim>
        
                                         /**
                                          * Tensors of covariant
-                                         * transformation.
+                                         * transformation at each of
+                                         * the quadrature points. The
+                                         * matrix stored is the
+                                         * inverse of the Jacobian
+                                         * matrix, which itself is
+                                         * stored in the
+                                         * @p{contravariant} field of
+                                         * this structure.
                                          *
                                          * Computed on each cell.
                                          */
@@ -245,7 +273,12 @@ class MappingQ1 : public Mapping<dim>
        
                                         /**
                                          * Tensors of covariant
-                                         * transformation.
+                                         * transformation at each of
+                                         * the quadrature points. The
+                                         * contravariant matrix is
+                                         * the Jacobian of the
+                                         * transformation,
+                                         * i.e. $J_ij=dx_i/d\hat x_j$.
                                          *
                                          * Computed on each cell.
                                          */
@@ -259,12 +292,12 @@ class MappingQ1 : public Mapping<dim>
                                          *
                                          * Filled once.
                                          */
-        vector2d<Tensor<1,dim> > unit_tangentials;
+        Table<2,Tensor<1,dim> > unit_tangentials;
        
                                         /**
                                          * Auxuliary vectors for internal use.
                                          */
-        vector2d<Tensor<1,dim> > aux;
+        Table<2,Tensor<1,dim> > aux;
 
                                         /**
                                          * Stores the support points of
index 003cc106154458ebb4d308a016b562b03b1e51f2..c2d0260095dd63e7174b3e3c15c8dc1eba20935d 100644 (file)
@@ -95,8 +95,8 @@ class MappingQ1Eulerian : public MappingQ1<dim>
                                      * can be initialized by
                                      * @p{DoFObjectAccessor::set_dof_values()}.
                                      */
-    MappingQ1Eulerian ( const Vector<double>  &euler_transform_vectors,
-                       const DoFHandler<dim> &shiftmap_dof_handler);
+    MappingQ1Eulerian (const Vector<double>  &euler_transform_vectors,
+                       const DoFHandler<dim> &shiftmap_dof_handler);
 
 
                                     /**
index cc90555189a0cf1022456a8551f02e81edd238f5..7fa2d74579879c65f1109899cc7ef486f5f2367d 100644 (file)
@@ -460,6 +460,33 @@ MappingCartesian<dim>::transform_covariant (Tensor<1,dim>       *begin,
 }
 
 
+
+template <int dim>
+void
+MappingCartesian<dim>::transform_covariant (Tensor<2,dim>       *begin,
+                                           Tensor<2,dim>       *end,
+                                           const Tensor<2,dim> *src,
+                                           const typename Mapping<dim>::InternalDataBase &mapping_data) const
+{
+  const InternalData &data = dynamic_cast<const InternalData&> (mapping_data);
+
+  Assert (data.update_flags & update_covariant_transformation,
+         typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+  
+                                  // simply scale by inverse Jacobian
+                                  // (which is diagonal here)
+  while (begin!=end)
+    {
+      for (unsigned int d=0; d<dim; ++d)
+        for (unsigned int p=0; p<dim; ++p)
+       (*begin)[d][p] = (*src)[d][p] / data.length[p];
+      begin++;
+      src++;
+    }
+}
+
+
+
 template <int dim>
 void
 MappingCartesian<dim>::transform_contravariant (Tensor<1,dim>       *begin,
@@ -488,6 +515,37 @@ MappingCartesian<dim>::transform_contravariant (Tensor<1,dim>       *begin,
 }
 
 
+
+template <int dim>
+void
+MappingCartesian<dim>::transform_contravariant (Tensor<2,dim>       *begin,
+                                               Tensor<2,dim>       *end,
+                                               const Tensor<2,dim> *src,
+  const typename Mapping<dim>::InternalDataBase &mapping_data) const
+{
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
+  const InternalData &data = dynamic_cast<const InternalData&> (mapping_data);
+
+  Assert (data.update_flags & update_contravariant_transformation,
+         typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+
+                                  // simply scale by Jacobian
+                                  // (which is diagonal here)
+  while (begin!=end)
+    {
+      for (unsigned int d=0; d<dim; ++d)
+        for (unsigned int p=0; p<dim; ++p)
+       (*begin)[d][p] = data.length[d] * (*src)[d][p];
+      begin++;
+      src++;
+    }
+}
+
+
+
 template <int dim>
 Point<dim> MappingCartesian<dim>::transform_unit_to_real_cell (
   const typename Triangulation<dim>::cell_iterator cell,
index d55bee4885d4980f2950c06cdf8a0c73ce3442f0..b0d5ac2ab7a799662cc506802fb8655c80df89bd 100644 (file)
@@ -61,8 +61,6 @@ MappingQ<dim>::InternalData::memory_consumption () const
 // cells are scaled linearly
 template<>
 MappingQ<1>::MappingQ (const unsigned int):
-               laplace_on_quad_vector(0),
-               laplace_on_hex_vector(0),
                degree(1),
                n_inner(0),
                n_outer(0),
@@ -100,9 +98,8 @@ static number power(const number x, const unsigned int y)
 
 
 template<int dim>
-MappingQ<dim>::MappingQ (const unsigned int p):
-               laplace_on_quad_vector(0),
-               laplace_on_hex_vector(0),
+MappingQ<dim>::MappingQ (const unsigned int p)
+                :
                degree(p),
                n_inner(power(degree-1, dim)),
                n_outer((dim==2) ? 4+4*(degree-1)
@@ -409,7 +406,7 @@ MappingQ<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_iter
 
 template <>
 void
-MappingQ<1>::set_laplace_on_quad_vector(vector2d<double> &) const
+MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const
 {
   Assert(false, ExcInternalError());
 }
@@ -418,7 +415,7 @@ MappingQ<1>::set_laplace_on_quad_vector(vector2d<double> &) const
 
 template <int dim>
 void
-MappingQ<dim>::set_laplace_on_quad_vector(vector2d<double> &loqvs) const
+MappingQ<dim>::set_laplace_on_quad_vector(Table<2,double> &loqvs) const
 {
   Assert(degree>1, ExcInternalError());
   const unsigned int n_inner_2d=(degree-1)*(degree-1);
@@ -490,7 +487,7 @@ MappingQ<dim>::set_laplace_on_quad_vector(vector2d<double> &loqvs) const
 
 template <>
 void
-MappingQ<3>::set_laplace_on_hex_vector(vector2d<double> &lohvs) const
+MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const
 {
   Assert(degree>1, ExcInternalError());
 
@@ -536,7 +533,7 @@ MappingQ<3>::set_laplace_on_hex_vector(vector2d<double> &lohvs) const
 
 template <int dim>
 void
-MappingQ<dim>::set_laplace_on_hex_vector(vector2d<double> &) const
+MappingQ<dim>::set_laplace_on_hex_vector(Table<2,double> &) const
 {
   Assert(false, ExcInternalError());
 }
@@ -548,7 +545,7 @@ MappingQ<dim>::set_laplace_on_hex_vector(vector2d<double> &) const
 
 template <>
 void
-MappingQ<1>::compute_laplace_vector(vector2d<double> &) const
+MappingQ<1>::compute_laplace_vector(Table<2,double> &) const
 {
   Assert(false, ExcInternalError());
 }
@@ -558,7 +555,7 @@ MappingQ<1>::compute_laplace_vector(vector2d<double> &) const
 
 template <int dim>
 void
-MappingQ<dim>::compute_laplace_vector(vector2d<double> &lvs) const
+MappingQ<dim>::compute_laplace_vector(Table<2,double> &lvs) const
 {
   Assert(lvs.n_rows()==0, ExcInternalError());
   Assert(dim==2 || dim==3, ExcNotImplemented());
@@ -617,7 +614,7 @@ MappingQ<dim>::compute_laplace_vector(vector2d<double> &lvs) const
 
 template <int dim>
 void
-MappingQ<dim>::apply_laplace_vector(const vector2d<double> &lvs,
+MappingQ<dim>::apply_laplace_vector(const Table<2,double> &lvs,
                                    std::vector<Point<dim> > &a) const
 {
   Assert(lvs.n_rows()!=0, ExcLaplaceVectorNotSet(degree));
@@ -1086,6 +1083,40 @@ MappingQ<dim>::transform_covariant (Tensor<1,dim>       *begin,
 }
 
 
+
+template <int dim>
+void
+MappingQ<dim>::transform_covariant (Tensor<2,dim>       *begin,
+                                   Tensor<2,dim>       *end,
+                                   const Tensor<2,dim> *src,
+                                   const typename Mapping<dim>::InternalDataBase &mapping_data) const
+{
+  const typename MappingQ1<dim>::InternalData *q1_data =
+    dynamic_cast<const typename MappingQ1<dim>::InternalData *> (&mapping_data);
+  Assert(q1_data!=0, ExcInternalError());
+  
+  typename std::vector<Tensor<2,dim> >::const_iterator tensor;
+
+  if (q1_data->is_mapping_q1_data)
+    tensor = q1_data->covariant.begin();
+  else
+    {
+      const InternalData *data = dynamic_cast<const InternalData *> (q1_data);
+      Assert(data!=0, ExcInternalError());
+
+      if (data->use_mapping_q1_on_current_cell)
+       tensor = data->mapping_q1_data.covariant.begin();
+      else
+       tensor = data->covariant.begin();
+    }
+  while (begin!=end)
+    {
+      contract (*(begin++), *(src++), *(tensor++));
+    }
+}
+
+
+
 template <int dim>
 void
 MappingQ<dim>::transform_contravariant (Tensor<1,dim>       *begin,
@@ -1118,6 +1149,40 @@ MappingQ<dim>::transform_contravariant (Tensor<1,dim>       *begin,
 }
 
 
+
+template <int dim>
+void
+MappingQ<dim>::transform_contravariant (Tensor<2,dim>       *begin,
+                                       Tensor<2,dim>       *end,
+                                       const Tensor<2,dim> *src,
+                                       const typename Mapping<dim>::InternalDataBase &mapping_data) const
+{
+  const typename MappingQ1<dim>::InternalData *q1_data =
+    dynamic_cast<const typename MappingQ1<dim>::InternalData *> (&mapping_data);
+  Assert(q1_data!=0, ExcInternalError());
+  
+  typename std::vector<Tensor<2,dim> >::const_iterator tensor;
+
+  if (q1_data->is_mapping_q1_data)
+    tensor = q1_data->contravariant.begin();
+  else
+    {
+      const InternalData *data = dynamic_cast<const InternalData *> (q1_data);
+      Assert(data!=0, ExcInternalError());
+
+      if (data->use_mapping_q1_on_current_cell)
+       tensor = data->mapping_q1_data.contravariant.begin();
+      else
+       tensor = data->contravariant.begin();    
+    }
+  while (begin!=end)
+    {
+      contract (*(begin++), *(tensor++), *(src++));
+    }
+}
+
+
+
 template <int dim>
 Point<dim> MappingQ<dim>::transform_unit_to_real_cell (
   const typename Triangulation<dim>::cell_iterator cell,
index bc31b613961eae82e9b689e03d4e3be2c3372fd5..8a31a3129b1dde3c172ac93754de64adcdf5e640 100644 (file)
@@ -467,47 +467,45 @@ MappingQ1<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterator &cel
 //           ExcDimensionMismatch(covariant_grads.size(), npts));
     }
   
-  std::vector<Point<dim> > &a=data.mapping_support_points;
-  
-                                  // store all Lagrangian
-                                  // support points in a
-  if (a.size()==0
-      || (&cell->get_triangulation() !=
-         &data.cell_of_current_support_points->get_triangulation())
-      || (cell!=data.cell_of_current_support_points))
+                                  // if necessary, recompute the
+                                  // support points of the
+                                  // transformation of this cell
+  if ((data.mapping_support_points.size() == 0)
+      ||
+      (&cell->get_triangulation() !=
+       &data.cell_of_current_support_points->get_triangulation())
+      ||
+      (cell != data.cell_of_current_support_points))
     {
-      compute_mapping_support_points(cell, a);
-      data.cell_of_current_support_points=cell;
+      compute_mapping_support_points(cell, data.mapping_support_points);
+      data.cell_of_current_support_points = cell;
     }
+
+                                   // first compute quadrature points
+  if (update_flags & update_q_points)
+    for (unsigned int point=0; point<npts; ++point)
+      for (unsigned int k=0; k<data.n_shape_functions; ++k)
+        quadrature_points[point]
+          += data.shape(point+offset,k) * data.mapping_support_points[k];
   
-  for (unsigned int point=0; point<npts; ++point)
-    {
-                                      // First, compute function
-                                      // values and derivatives
+                                   // then Jacobians
+  if (update_flags & update_contravariant_transformation)
+    for (unsigned int point=0; point<npts; ++point)
       for (unsigned int k=0; k<data.n_shape_functions; ++k)
-       {
-         if (update_flags & update_q_points)
-           {
-             quadrature_points[point]
-               += data.shape(point+offset,k) * a[k];
-           }
-         if (update_flags & update_contravariant_transformation)
-           {
-             for (unsigned int i=0; i<dim; ++i)
-               for (unsigned int j=0; j<dim; ++j)
-                 data.contravariant[point][i][j]
-                   += data.derivative(point+offset, k)[j]
-                   * a[k][i];
-           }
-       }
-
-                                      // Invert contravariant for
+        for (unsigned int i=0; i<dim; ++i)
+          for (unsigned int j=0; j<dim; ++j)
+            data.contravariant[point][i][j]
+              += (data.derivative(point+offset, k)[j]
+                  *
+                  data.mapping_support_points[k][i]);
+
+                                      // invert contravariant for
                                       // covariant transformation
                                       // matrices
-      if (update_flags & update_covariant_transformation)
-       data.covariant[point]
-         = invert(data.contravariant[point]);
-    }
+  if (update_flags & update_covariant_transformation)
+    for (unsigned int point=0; point<npts; ++point)
+      data.covariant[point] = invert(data.contravariant[point]);
+
   data.first_cell = false;
 }
 
@@ -798,7 +796,37 @@ MappingQ1<dim>::transform_covariant (Tensor<1,dim>       *begin,
 //  Assert (dst.size() == data.contravariant.size(),
 //       ExcDimensionMismatch(dst.size() + src_offset, data.contravariant.size()));
 
-  typename std::vector<Tensor<2,dim> >::const_iterator tensor = data.covariant.begin();
+  typename std::vector<Tensor<2,dim> >::const_iterator
+    tensor = data.covariant.begin();
+  
+  while (begin!=end)
+    {
+      contract (*(begin++), *(src++), *(tensor++));
+    }
+}
+
+
+
+template <int dim>
+void
+MappingQ1<dim>::transform_covariant (Tensor<2,dim>       *begin,
+                                    Tensor<2,dim>       *end,
+                                    const Tensor<2,dim> *src,
+                                    const typename Mapping<dim>::InternalDataBase &mapping_data) const
+{
+  const InternalData *data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
+  Assert(data_ptr!=0, ExcInternalError());
+  const InternalData &data=*data_ptr;
+
+  Assert (data.update_flags & update_covariant_transformation,
+         typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+
+//TODO: [GK] Can we do a similar assertion?  
+//  Assert (dst.size() == data.contravariant.size(),
+//       ExcDimensionMismatch(dst.size() + src_offset, data.contravariant.size()));
+
+  typename std::vector<Tensor<2,dim> >::const_iterator
+    tensor = data.covariant.begin();
   
   while (begin!=end)
     {
@@ -807,6 +835,7 @@ MappingQ1<dim>::transform_covariant (Tensor<1,dim>       *begin,
 }
 
 
+
 template <int dim>
 void
 MappingQ1<dim>::transform_contravariant (Tensor<1,dim>       *begin,
@@ -821,7 +850,8 @@ MappingQ1<dim>::transform_contravariant (Tensor<1,dim>       *begin,
   Assert (data.update_flags & update_contravariant_transformation,
          typename FEValuesBase<dim>::ExcAccessToUninitializedField());
 
-  typename std::vector<Tensor<2,dim> >::const_iterator tensor = data.contravariant.begin();
+  typename std::vector<Tensor<2,dim> >::const_iterator
+    tensor = data.contravariant.begin();
   
   while (begin!=end)
     {
@@ -829,6 +859,32 @@ MappingQ1<dim>::transform_contravariant (Tensor<1,dim>       *begin,
     }
 }
 
+
+template <int dim>
+void
+MappingQ1<dim>::transform_contravariant (Tensor<2,dim>       *begin,
+                                        Tensor<2,dim>       *end,
+                                        const Tensor<2,dim> *src,
+                                        const typename Mapping<dim>::InternalDataBase &mapping_data) const
+{
+  const InternalData* data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
+  Assert(data_ptr!=0, ExcInternalError());
+  const InternalData &data=*data_ptr;
+
+  Assert (data.update_flags & update_contravariant_transformation,
+         typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+
+  typename std::vector<Tensor<2,dim> >::const_iterator
+    tensor = data.contravariant.begin();
+  
+  while (begin!=end)
+    {
+      contract (*(begin++), *(tensor++), *(src++));
+    }
+}
+
+
+
 template <int dim>
 Point<dim> MappingQ1<dim>::transform_unit_to_real_cell (
   const typename Triangulation<dim>::cell_iterator cell,

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.