]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add the functionality inverse_jacobian to the FEValues class, which gives the inverse...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 13 Aug 2008 22:05:24 +0000 (22:05 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 13 Aug 2008 22:05:24 +0000 (22:05 +0000)
git-svn-id: https://svn.dealii.org/trunk@16539 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_update_flags.h
deal.II/deal.II/include/fe/fe_values.h
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/source/fe/fe_values.cc
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 5d6711e033f661f9cfb193db150324b399d74721..027ee86bd1d3f42b532b3459360c7229610251ca 100644 (file)
@@ -165,6 +165,15 @@ enum UpdateFlags
                                        * transformation.
                                        */
       update_jacobian_grads               = 0x0100,
+                                      //! Volume element
+                                      /**
+                                       * Compute the inverse 
+                                       * Jacobian of the
+                                       * transformation from the
+                                       * reference cell to the real
+                                       * cell.
+                                       */
+      update_inverse_jacobians            = 0x0200,
                                       //! Covariant transformation
                                       /**
                                        * Compute all values the
@@ -174,9 +183,9 @@ enum UpdateFlags
                                        * mappings like
                                        * MappingCartesian this may be
                                        * simpler than
-                                       * #update_jacobians.
+                                       * #update_inverse_jacobians.
                                        */
-      update_covariant_transformation     = 0x0200,
+      update_covariant_transformation     = 0x0400,
                                       //! Contravariant transformation
                                       /**
                                        * Compute all values the
@@ -188,14 +197,14 @@ enum UpdateFlags
                                        * simpler than
                                        * #update_jacobians.
                                        */
-      update_contravariant_transformation = 0x0400,
+      update_contravariant_transformation = 0x0800,
                                       //! Shape function values of transformation
                                       /**
                                        * Compute the shape function
                                        * values of the transformation
                                        * defined by the Mapping.
                                        */
-      update_transformation_values        = 0x0800,
+      update_transformation_values        = 0x1000,
                                       //! Shape function gradients of transformation
                                       /**
                                        * Compute the shape function
@@ -203,7 +212,7 @@ enum UpdateFlags
                                        * transformation defined by
                                        * the Mapping.
                                        */
-      update_transformation_gradients     = 0x1000,
+      update_transformation_gradients     = 0x2000,
 
                                       /**
                                        * Compute the JxW values
@@ -214,25 +223,25 @@ enum UpdateFlags
                                        * is used in conjunction with
                                        * H_div subspaces like RT and ABF.
                                        */
-      update_cell_JxW_values              = 0x2000,
+      update_cell_JxW_values              = 0x4000,
                                       /**
                                        * Update the location of the
                                        * mapped generalized support
                                        * points of the element.
                                        */
-      update_support_points               = 0x4000,
+      update_support_points               = 0x8000,
                                       /**
                                        * Update the Jacobian of the
                                        * mapping in generalized
                                        * support points.
                                        */
-      update_support_jacobians            = 0x8000,
+      update_support_jacobians            = 0x10000,
                                       /**
                                        * Update the inverse Jacobian
                                        * of the mapping in
                                        * generalized support points.
                                        */
-      update_support_inverse_jacobians    = 0x10000,
+      update_support_inverse_jacobians    = 0x20000,
                                       /**
                                        * @deprecated Update
                                        * quadrature points
index 54edbc94e98c798b9330da057559a95d016c5a5b..cc843b40ca1096d69a83253b823f22f95b350e1a 100644 (file)
@@ -562,13 +562,19 @@ class FEValuesData
                                      * quadrature points.
                                      */
     std::vector<Tensor<2,dim> > jacobians;
-    
+
                                     /**
                                      * Array of the derivatives of the Jacobian
                                      * matrices at the quadrature points.
                                      */
     std::vector<Tensor<3,dim> > jacobian_grads;
-    
+
+                                    /**
+                                     * Array of the inverse Jacobian matrices 
+                                     * at the quadrature points.
+                                     */
+    std::vector<Tensor<2,dim> > inverse_jacobians;
+
                                     /**
                                      * Store an array of weights
                                      * times the Jacobi determinant
@@ -1633,7 +1639,22 @@ class FEValuesBase : protected FEValuesData<dim>,
                                      * jacobian_grads().
                                      */
     const std::vector<Tensor<3,dim> > & get_jacobian_grads () const;
-    
+
+                                    /**
+                                     * Return the inverse Jacobian of the
+                                     * transformation at the specified
+                                     * quadrature point, i.e.
+                                     * $J_{ij}=d\hat x_i/dx_j$
+                                     */
+    const Tensor<2,dim> & inverse_jacobian (const unsigned int quadrature_point) const;
+
+                                    /**
+                                     * Pointer to the array holding
+                                     * the values returned by 
+                                     * inverse_jacobian().
+                                     */
+    const std::vector<Tensor<2,dim> > & get_inverse_jacobians () const;
+
                                     /**
                                      * Constant reference to the
                                      * selected mapping object.
@@ -4038,6 +4059,17 @@ FEValuesBase<dim>::get_jacobian_grads () const
 
 
 
+template <int dim>
+inline
+const std::vector<Tensor<2,dim> >&
+FEValuesBase<dim>::get_inverse_jacobians () const
+{
+  Assert (this->update_flags & update_inverse_jacobians, ExcAccessToUninitializedField());
+  return this->inverse_jacobians;
+}
+
+
+
 template <int dim>
 inline
 const Point<dim> &
@@ -4091,6 +4123,19 @@ FEValuesBase<dim>::jacobian_grad (const unsigned int i) const
 
 
 
+template <int dim>
+inline
+const Tensor<2,dim> &
+FEValuesBase<dim>::inverse_jacobian (const unsigned int i) const
+{
+  Assert (this->update_flags & update_inverse_jacobians, ExcAccessToUninitializedField());
+  Assert (i<this->inverse_jacobians.size(), ExcIndexRange(i, 0, this->inverse_jacobians.size()));
+  
+  return this->inverse_jacobians[i];
+}
+
+
+
 template <int dim>
 template <class InputVector>
 inline
index 88c3992a00c6b69cbdb400e7a059e0a14f54365c..efc90904bfcab0fe8dc5d882296cb6c331cc3de9 100644 (file)
@@ -484,9 +484,11 @@ class Mapping : public Subscriptor
                                      * matrices needed to transform
                                      * vector-valued functions,
                                      * namely
-                                     * @p jacobians
-                                     * and the derivatives
-                                     * @p jacobian_grads.
+                                     * @p jacobians,
+                                     * the derivatives
+                                     * @p jacobian_grads,
+                                     * and the inverse operations in
+                                     * @p inverse_jacobians.
                                      */
     virtual void
     fill_fe_values (const typename Triangulation<dim>::cell_iterator &cell,
@@ -495,7 +497,8 @@ class Mapping : public Subscriptor
                    std::vector<Point<dim> >                      &quadrature_points,
                    std::vector<double>                           &JxW_values,
                    std::vector<Tensor<2,dim> >                   &jacobians,
-                   std::vector<Tensor<3,dim> >                   &jacobian_grads) const = 0;
+                   std::vector<Tensor<3,dim> >                   &jacobian_grads,
+                   std::vector<Tensor<2,dim> >                   &inverse_jacobians) const = 0;
 
                                     /**
                                      * Performs the same as @p fill_fe_values
index d9bf66f544ed614244cd2ca4951eda2bf8be45c7..2ef9a7d79ff68437644d57760637d8589988e4f8 100644 (file)
@@ -61,7 +61,8 @@ class MappingCartesian : public Mapping<dim>
                    std::vector<Point<dim> >        &quadrature_points,
                    std::vector<double>             &JxW_values,
                    std::vector<Tensor<2,dim> >     &jacobians,
-                   std::vector<Tensor<3,dim> >     &jacobian_grads) const ;
+                   std::vector<Tensor<3,dim> >     &jacobian_grads,
+                   std::vector<Tensor<2,dim> >     &inverse_jacobians) const ;
 
     virtual void
     fill_fe_face_values (const typename Triangulation<dim>::cell_iterator &cell,
index 2f903a3e4f34b98e931e476beda9f9105f442d96..53b3290ba643c00a1e210a92496d1a3420b5ac5b 100644 (file)
@@ -233,7 +233,8 @@ class MappingQ : public MappingQ1<dim>
                    typename std::vector<Point<dim> >             &quadrature_points,
                    std::vector<double>                  &JxW_values,
                    std::vector<Tensor<2,dim> >          &jacobians,
-                   std::vector<Tensor<3,dim> >          &jacobian_grads) const ;
+                   std::vector<Tensor<3,dim> >          &jacobian_grads,
+                   std::vector<Tensor<2,dim> >          &inverse_jacobians) const ;
 
                                     /**
                                      * Implementation of the interface in
index c600375b5eb937cddbdddf0ef43aa544ef401cf5..293609b83c354d178bba9f15ef285ff4a7e2b921 100644 (file)
@@ -309,7 +309,8 @@ class MappingQ1 : public Mapping<dim>
                    typename std::vector<Point<dim> >       &quadrature_points,
                    std::vector<double>             &JxW_values,
                    std::vector<Tensor<2,dim> >     &jacobians,
-                   std::vector<Tensor<3,dim> >     &jacobian_grads) const ;
+                   std::vector<Tensor<3,dim> >     &jacobian_grads,
+                   std::vector<Tensor<2,dim> >     &inverse_jacobians) const ;
 
                                     /**
                                      * Implementation of the interface in
index 991cd04fd8c6178aca4722920cc4ef184aee7253..8c619ec4c8762b7aa11b87e134e06f42267ee824 100644 (file)
@@ -314,6 +314,9 @@ FEValuesData<dim>::initialize (const unsigned int        n_quadrature_points,
   if (flags & update_jacobian_grads)
     this->jacobian_grads.resize(n_quadrature_points);
 
+  if (flags & update_inverse_jacobians)
+    this->inverse_jacobians.resize(n_quadrature_points);
+
   if (flags & update_boundary_forms)
     this->boundary_forms.resize(n_quadrature_points);
 
@@ -1386,7 +1389,8 @@ void FEValues<dim>::do_reinit ()
                                     this->quadrature_points,
                                     this->JxW_values,
                                     this->jacobians,
-                                    this->jacobian_grads);
+                                    this->jacobian_grads,
+                                    this->inverse_jacobians);
   
   this->get_fe().fill_fe_values(this->get_mapping(),
                                *this->present_cell,
index 323b5afc24ed23c5482a3bcfe0e183e0097abf62..b36b23d0047bd37160a3d54e2173aefdce4c7158 100644 (file)
@@ -303,7 +303,8 @@ fill_fe_values (const typename Triangulation<dim>::cell_iterator& cell,
                 std::vector<Point<dim> >& quadrature_points,
                 std::vector<double>& JxW_values,
                std::vector<Tensor<2,dim> >& jacobians,
-               std::vector<Tensor<3,dim> >& jacobian_grads) const
+               std::vector<Tensor<3,dim> >& jacobian_grads,
+               std::vector<Tensor<2,dim> >& inverse_jacobians) const
 {
                                   // convert data object to internal
                                   // data for this class. fails with
@@ -346,6 +347,16 @@ fill_fe_values (const typename Triangulation<dim>::cell_iterator& cell,
   if (data.current_update_flags() & update_jacobian_grads)
     for (unsigned int i=0; i<jacobian_grads.size();++i)
       jacobian_grads[i]=Tensor<3,dim>();
+                                  // "compute" inverse Jacobian at the 
+                                  // quadrature points, which are all 
+                                  // the same
+  if (data.current_update_flags() & update_inverse_jacobians)
+    for (unsigned int i=0; i<inverse_jacobians.size();++i)
+      {
+       inverse_jacobians[i]=Tensor<2,dim>();
+       for (unsigned int j=0; j<dim; ++j)
+         inverse_jacobians[j][j]=1./data.length[j];
+      }
 }
 
 
index 66e4899d8793fbbdcd27b2917826f878e9dbacbd..c4053122f652db4d25d65b800124a61c1e8fd1c0 100644 (file)
@@ -307,7 +307,8 @@ MappingQ<dim>::fill_fe_values (const typename Triangulation<dim>::cell_iterator
                               std::vector<Point<dim> >             &quadrature_points,
                               std::vector<double>                  &JxW_values,
                               std::vector<Tensor<2,dim> >          &jacobians,
-                              std::vector<Tensor<3,dim> >          &jacobian_grads) const
+                              std::vector<Tensor<3,dim> >          &jacobian_grads,
+                              std::vector<Tensor<2,dim> >          &inverse_jacobians) const
 {
                                   // convert data object to internal
                                   // data for this class. fails with
@@ -334,7 +335,7 @@ MappingQ<dim>::fill_fe_values (const typename Triangulation<dim>::cell_iterator
   
   MappingQ1<dim>::fill_fe_values(cell, q, *p_data,
                                 quadrature_points, JxW_values,
-                                jacobians, jacobian_grads);
+                                jacobians, jacobian_grads, inverse_jacobians);
 }
 
 
index 877748f94713bfa19980b48d2119146de65faf50..fec6d69eb458df854d6faa82979c4ecc5259892b 100644 (file)
@@ -343,7 +343,8 @@ MappingQ1<dim>::update_once (const UpdateFlags in) const
            | update_boundary_forms
            | update_normal_vectors
            | update_jacobians
-           | update_jacobian_grads))
+           | update_jacobian_grads
+           | update_inverse_jacobians))
     out |= update_transformation_gradients;
 
   return out;
@@ -365,7 +366,8 @@ MappingQ1<dim>::update_each (const UpdateFlags in) const
                                      | update_boundary_forms
                                      | update_normal_vectors
                                      | update_jacobians
-                                     | update_jacobian_grads));
+                                     | update_jacobian_grads
+                                     | update_inverse_jacobians));
 
                                   // add a few flags. note that some
                                   // flags appear in both conditions
@@ -373,10 +375,10 @@ MappingQ1<dim>::update_each (const UpdateFlags in) const
                                   // operations. this leads to some
                                   // circular logic. the only way to
                                   // treat this is to iterate. since
-                                  // there are 3 if-clauses in the
+                                  // there are 4 if-clauses in the
                                   // loop, it will take at most 3
                                   // iterations to converge. do them:
-  for (unsigned int i=0; i<3; ++i)
+  for (unsigned int i=0; i<4; ++i)
     {
                                       // The following is a little incorrect:
                                       // If not applied on a face,
@@ -398,6 +400,9 @@ MappingQ1<dim>::update_each (const UpdateFlags in) const
                 | update_boundary_forms
                 | update_normal_vectors))
        out |= update_contravariant_transformation;
+      
+      if (out & (update_inverse_jacobians))
+       out |= update_covariant_transformation;
 
                                       // The contravariant transformation
                                       // is a Piola transformation, which
@@ -672,12 +677,13 @@ MappingQ1<dim>::compute_mapping_support_points(
 template <int dim>
 void
 MappingQ1<dim>::fill_fe_values (const typename Triangulation<dim>::cell_iterator &cell,
-                               const Quadrature<dim>                &q,
-                               typename Mapping<dim>::InternalDataBase      &mapping_data,
+                               const Quadrature<dim>                     &q,
+                               typename Mapping<dim>::InternalDataBase   &mapping_data,
                                std::vector<Point<dim> >                  &quadrature_points,
                                std::vector<double>                       &JxW_values,
                                std::vector<Tensor<2,dim> >               &jacobians,
-                               std::vector<Tensor<3,dim> >               &jacobian_grads) const
+                               std::vector<Tensor<3,dim> >               &jacobian_grads,
+                               std::vector<Tensor<2,dim> >               &inverse_jacobians) const
 {
   InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
   Assert(data_ptr!=0, ExcInternalError());
@@ -733,6 +739,16 @@ MappingQ1<dim>::fill_fe_values (const typename Triangulation<dim>::cell_iterator
                  += (data.second_derivative(point+DataSetDescriptor::cell (), k)[j][l]
                      *
                      data.mapping_support_points[k][i]);
+    }
+                                  // copy values from InternalData to vector
+                                  // given by reference
+  if (update_flags & update_inverse_jacobians)
+    {      
+      Assert (inverse_jacobians.size() == n_q_points,
+             ExcDimensionMismatch(inverse_jacobians.size(), n_q_points));
+      for (unsigned int point=0; point<n_q_points; ++point)
+       inverse_jacobians[point]
+         = transpose(data.covariant[point]);
     }
 }
 

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.