]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Convert part of the library from Point and dFMatrix to Tensor where necessary.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 6 Nov 1998 19:05:39 +0000 (19:05 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 6 Nov 1998 19:05:39 +0000 (19:05 +0000)
git-svn-id: https://svn.dealii.org/trunk@652 0785d39b-7218-0410-832d-ea1e28bc413d

18 files changed:
deal.II/deal.II/Attic/examples/poisson/equation.cc
deal.II/deal.II/Todo
deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_lib.criss_cross.h
deal.II/deal.II/include/fe/fe_lib.lagrange.h
deal.II/deal.II/include/fe/fe_linear_mapping.h
deal.II/deal.II/include/fe/fe_update_flags.h
deal.II/deal.II/include/fe/fe_values.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_lib.criss_cross.cc
deal.II/deal.II/source/fe/fe_lib.cubic.cc
deal.II/deal.II/source/fe/fe_lib.linear.cc
deal.II/deal.II/source/fe/fe_lib.quadratic.cc
deal.II/deal.II/source/fe/fe_lib.quartic.cc
deal.II/deal.II/source/fe/fe_linear_mapping.cc
deal.II/deal.II/source/fe/fe_values.cc
deal.II/deal.II/source/numerics/matrices.cc
tests/big-tests/poisson/equation.cc

index 7844de59d5ab6951eb18bdf71b790c3ca91c2b9d..e30bfefa7677925b7c708998665549615a707194 100644 (file)
@@ -36,7 +36,7 @@ void PoissonEquation<2>::assemble (dFMatrix            &cell_matrix,
                                   dVector             &rhs,
                                   const FEValues<2>   &fe_values,
                                   const Triangulation<2>::cell_iterator &) const {
-  const vector<vector<Point<2> > >&gradients = fe_values.get_shape_grads ();
+  const vector<vector<Tensor<1,2> > >&gradients = fe_values.get_shape_grads ();
   const dFMatrix       &values    = fe_values.get_shape_values ();
   vector<double>        rhs_values (fe_values.n_quadrature_points);
   const vector<double> &weights   = fe_values.get_JxW_values ();
index 53360beabe8be456749c4efd5ff8696e57ba3929..5ceb7ec3972ecd8d2d4b0b951d847092716decdf 100644 (file)
@@ -175,3 +175,11 @@ Review DoFCellAccessor<dim>::get_interpolated_dof_values: it may be
 Move the exceptions in the Tensor package (tensor_base.h and tensor.h)
   back into the Tensor class, if the compiler allows to do so. Also
   rename them back (i.e. drop the initial Tensor* from the name).
+
+
+
+Think about the determinant function in the tensor package. Is it
+  useful, can it be generalized?
+
+
+Implement FELinearMapping::fill_fe_values->jacobians
\ No newline at end of file
index b48770180b20c107d63acfc0437b738abd6abcb9..89f92709280b753c2396c0f973a974250f8956e1 100644 (file)
@@ -550,8 +550,8 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      * function at the point #p#.
                                      * #p# is a point on the reference element,
                                      */
-    virtual Point<dim> shape_grad (const unsigned int  i,
-                                  const Point<dim>   &p) const = 0;
+    virtual Tensor<1,dim> shape_grad (const unsigned int  i,
+                                     const Point<dim>   &p) const = 0;
 
                                     /**
                                      * Return the tensor of second derivatives
@@ -639,7 +639,7 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      */
     virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
                                 const vector<Point<dim> >            &unit_points,
-                                vector<dFMatrix>    &jacobians,
+                                vector<Tensor<2,dim> >               &jacobians,
                                 const bool           compute_jacobians,
                                 vector<Point<dim> > &support_points,
                                 const bool           compute_support_points,
@@ -757,7 +757,7 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      const unsigned int           face_no,
                                      const vector<Point<dim-1> > &unit_points,
                                      const vector<Point<dim> >   &global_unit_points,
-                                     vector<dFMatrix>    &jacobians,
+                                     vector<Tensor<2,dim> >      &jacobians,
                                      const bool           compute_jacobians,
                                      vector<Point<dim> > &support_points,
                                      const bool           compute_support_points,
@@ -804,7 +804,7 @@ class FiniteElement : public FiniteElementBase<dim> {
                                         const unsigned int           subface_no,
                                         const vector<Point<dim-1> > &unit_points,
                                         const vector<Point<dim> >   &global_unit_points,
-                                        vector<dFMatrix>    &jacobians,
+                                        vector<Tensor<2,dim> >      &jacobians,
                                         const bool           compute_jacobians,
                                         vector<Point<dim> > &q_points,
                                         const bool           compute_q_points,
@@ -863,8 +863,8 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      * overwritten function has to be provided.
                                      */
     virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
-                                   const Boundary<dim> &boundary,
-                                   vector<Point<dim> > &support_points) const;
+                                    const Boundary<dim> &boundary,
+                                    vector<Point<dim> > &support_points) const;
     
                                     /**
                                      * Compute the off-points of the finite
@@ -907,8 +907,8 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      * when called with #dim==1#.
                                      */
     virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
-                                        const Boundary<dim> &boundary,
-                                        vector<Point<dim> > &support_points) const =0;
+                                         const Boundary<dim> &boundary,
+                                         vector<Point<dim> > &support_points) const =0;
 
                                     /**
                                      * This is the second separated function
index 51da24e2c06fe56ce03583f617cf417d994226ff..d27931ed9625c56f166eb9cc7635c5687ba6fa68 100644 (file)
@@ -181,8 +181,8 @@ class FECrissCross : public FiniteElement<dim> {
                                      * Return the gradient of the #i#th shape
                                      * function at point #p# on the unit cell.
                                      */
-    virtual Point<dim> shape_grad(const unsigned int i,
-                                 const Point<dim>& p) const;
+    virtual Tensor<1,dim> shape_grad(const unsigned int i,
+                                    const Point<dim>& p) const;
 
                                     /**
                                      * Return the tensor of second derivatives
@@ -197,7 +197,7 @@ class FECrissCross : public FiniteElement<dim> {
                                      * need to do some evaluation by hand.
                                      */
     virtual Tensor<2,dim> shape_grad_grad (const unsigned int  i,
-                                          const Point<dim>   &p) const = 0;
+                                          const Point<dim>   &p) const;
 
                                     /**
                                      * Refer to the base class for detailed
@@ -332,7 +332,7 @@ class FECrissCross : public FiniteElement<dim> {
                                      */
     virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
                                 const vector<Point<dim> >            &unit_points,
-                                vector<dFMatrix>    &jacobians,
+                                vector<Tensor<2,dim> >               &jacobians,
                                 const bool           compute_jacobians,
                                 vector<Point<dim> > &support_points,
                                 const bool           compute_support_points,
index cc1bf1ed8c026bfb1a3ecd8c70431d3e09a8849b..e443652c2a5196d3ce0b6029bb2a82418ca1f42c 100644 (file)
@@ -44,8 +44,8 @@ class FELinear : public FELinearMapping<dim> {
                                      * Return the gradient of the #i#th shape
                                      * function at point #p# on the unit cell.
                                      */
-    virtual Point<dim> shape_grad(const unsigned int i,
-                                 const Point<dim>& p) const;
+    virtual Tensor<1,dim> shape_grad(const unsigned int i,
+                                    const Point<dim>& p) const;
 
                                     /**
                                      * Return the tensor of second derivatives
@@ -56,7 +56,7 @@ class FELinear : public FELinearMapping<dim> {
                                      * derivatives on the unit cell are zero.
                                      */
     virtual Tensor<2,dim> shape_grad_grad (const unsigned int  i,
-                                          const Point<dim>   &p) const = 0;
+                                          const Point<dim>   &p) const;
 
                                     /**
                                      * Refer to the base class for detailed
@@ -118,8 +118,8 @@ class FEQuadraticSub : public FELinearMapping<dim> {
                                      * Return the gradient of the #i#th shape
                                      * function at point #p# on the unit cell.
                                      */
-    virtual Point<dim> shape_grad(const unsigned int i,
-                                 const Point<dim>& p) const;
+    virtual Tensor<1,dim> shape_grad(const unsigned int i,
+                                    const Point<dim>& p) const;
 
                                     /**
                                      * Return the tensor of second derivatives
@@ -127,7 +127,7 @@ class FEQuadraticSub : public FELinearMapping<dim> {
                                      * point #p# on the unit cell.
                                      */
     virtual Tensor<2,dim> shape_grad_grad (const unsigned int  i,
-                                          const Point<dim>   &p) const = 0;
+                                          const Point<dim>   &p) const;
 
                                     /**
                                      * Refer to the base class for detailed
@@ -208,8 +208,8 @@ class FECubicSub : public FELinearMapping<dim> {
                                      * Return the gradient of the #i#th shape
                                      * function at point #p# on the unit cell.
                                      */
-    virtual Point<dim> shape_grad(const unsigned int i,
-                                 const Point<dim>& p) const;
+    virtual Tensor<1,dim> shape_grad(const unsigned int i,
+                                    const Point<dim>& p) const;
 
                                     /**
                                      * Return the tensor of second derivatives
@@ -217,7 +217,7 @@ class FECubicSub : public FELinearMapping<dim> {
                                      * point #p# on the unit cell.
                                      */
     virtual Tensor<2,dim> shape_grad_grad (const unsigned int  i,
-                                          const Point<dim>   &p) const = 0;
+                                          const Point<dim>   &p) const;
 
                                     /**
                                      * Refer to the base class for detailed
@@ -299,8 +299,8 @@ class FEQuarticSub : public FELinearMapping<dim> {
                                      * Return the gradient of the #i#th shape
                                      * function at point #p# on the unit cell.
                                      */
-    virtual Point<dim> shape_grad(const unsigned int i,
-                                 const Point<dim>& p) const;
+    virtual Tensor<1,dim> shape_grad(const unsigned int i,
+                                    const Point<dim>& p) const;
 
                                     /**
                                      * Return the tensor of second derivatives
@@ -308,7 +308,7 @@ class FEQuarticSub : public FELinearMapping<dim> {
                                      * point #p# on the unit cell.
                                      */
     virtual Tensor<2,dim> shape_grad_grad (const unsigned int  i,
-                                          const Point<dim>   &p) const = 0;
+                                          const Point<dim>   &p) const;
 
                                     /**
                                      * Refer to the base class for detailed
index 09476f8169653772f55f32358a8ba035c128c974..735023da01c908c2d46cdd69f019826aef68e68d 100644 (file)
@@ -135,7 +135,7 @@ class FELinearMapping : public FiniteElement<dim> {
                                      */
     virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
                                 const vector<Point<dim> >            &unit_points,
-                                vector<dFMatrix>    &jacobians,
+                                vector<Tensor<2,dim> >               &jacobians,
                                 const bool           compute_jacobians,
                                 vector<Point<dim> > &support_points,
                                 const bool           compute_support_points,
index 46d729bb595c7c797b99fa0592f97dccc5ccad4f..9e61b96271680773ab7af6e8129a485af8df10ff 100644 (file)
@@ -65,7 +65,13 @@ enum UpdateFlags {
                                        * an error, since normal vectors are
                                        * not useful in that case.
                                        */
-      update_normal_vectors = 32
+      update_normal_vectors = 32,
+
+                                      /**
+                                       * Update the second derivatives of the
+                                       * shape functions on the real cell.
+                                       */
+      update_second_derivatives
 };
 
 
index bf0768e08ec70cb2361658162851f469a1078e2e..12fa99453f11eedb7a673eb044c1112e80f79c57 100644 (file)
@@ -118,6 +118,22 @@ template <int dim> class Quadrature;
  *  #FELinearMapping<dim>::fill_fe_values# function, where also a small test
  *  program is presented.
  *
+ *  The derivatives of the Jacobi matrices at the quadrature points with respect
+ *  to unit cell coordinates is stored in the field names
+ *  #jacobi_matrices_grad#. Since the gradient of a shape function is given by
+ *  $\partial_i \phi = \sum_k  \hat\partial_k \hat\phi  J_{ki}$, where $\hat$
+ *  denotes differentiation on the unit cell, the second derivative of a
+ *  function is given by
+ *  $\partial_j \partial i \phi
+ *   =
+ *   \hat\partial_l [ (\hat \partial_k \hat\phi) J_{ki} ] J_{lj}
+ *   =
+ *   (\hat\partial_k \hat\partial_l \hat\phi) J_{ki} J_{lj}
+ *   +
+ *   (\hat \partial_l \hat\phi) (\hat\partial_l J_{ki}) J_{lj}$.
+ *  While we already have access to the Jacobian matrix, the derivatives are
+ *  stored in the named field.
+ *
  *  
  *  \subsection{Member functions}
  *
@@ -289,8 +305,8 @@ class FEValuesBase {
                                      * The function returns the gradient on the
                                      * real element, not the reference element.
                                      */
-    const Point<dim> & shape_grad (const unsigned int function,
-                                  const unsigned int quadrature_point) const;
+    const Tensor<1,dim> & shape_grad (const unsigned int function,
+                                     const unsigned int quadrature_point) const;
 
                                     /** 
                                      * Return a pointer to the matrix holding
@@ -299,7 +315,7 @@ class FEValuesBase {
                                      * For the format of this matrix, see the
                                      * documentation for the matrix itself.
                                      */
-    const vector<vector<Point<dim> > > & get_shape_grads () const;
+    const vector<vector<Tensor<1,dim> > > & get_shape_grads () const;
     
                                     /**
                                      * Return the gradients of the finite
@@ -458,7 +474,18 @@ class FEValuesBase {
                                      * gradients on the real element, rather
                                      * than on the reference element.
                                      */
-    vector<vector<Point<dim> > >  shape_gradients;
+    vector<vector<Tensor<1,dim> > >  shape_gradients;
+
+                                    /**
+                                     * Store the 2nd derivatives of the shape
+                                     * functions at the quadrature points.
+                                     *
+                                     * This field is reset each time
+                                     * #reinit# is called and contains the
+                                     * gradients on the real element, rather
+                                     * than on the reference element.
+                                     */
+    vector<vector<Tensor<2,dim> > >  shape_2nd_derivatives;
 
                                     /**
                                      * Store an array of the weights of the
@@ -516,8 +543,21 @@ class FEValuesBase {
                                      * transformation of the gradients to the
                                      * real cell.
                                      */
-    vector<dFMatrix>     jacobi_matrices;
+    vector<Tensor<2,dim> > jacobi_matrices;
 
+                                    /**
+                                     * Store the derivatives of the jacobi
+                                     * matrices. If #J[j][k]# is the jacobi
+                                     * matrix, then the index #[i][j][k]#
+                                     * of this field denotes the derivation
+                                     * of #J[j][k]# with respect to the
+                                     * #i#th variable.
+                                     *
+                                     * The same general remarks apply as for
+                                     * #jacobi_matrices#.
+                                     */
+    vector<Tensor<3,dim> > jacobi_matrices_grad;
+    
                                     /**
                                      * Store the values of the basis functions
                                      * of the transformation from unit cell
@@ -646,7 +686,17 @@ class FEValues : public FEValuesBase<dim> {
                                      * of the object and contains the gradients
                                      * on the reference element.
                                      */
-    vector<vector<Point<dim> > >   unit_shape_gradients;
+    vector<vector<Tensor<1,dim> > > unit_shape_gradients;
+
+                                    /**
+                                     * Store the 2nd derivatives of the shape
+                                     * functions at the quadrature points on
+                                     * the unit cell.
+                                     * This field is set up upon construction
+                                     * of the object and contains the
+                                     * derivatives on the reference element.
+                                     */
+    vector<vector<Tensor<2,dim> > > unit_shape_2nd_derivatives;
     
                                     /**
                                      * Gradients of the basis
@@ -793,8 +843,18 @@ class FEFaceValuesBase : public FEValuesBase<dim> {
                                      * subface, with indices like that:
                                      * #unit_shape_gradients[face][dof][q_point]#
                                      */
-    vector<vector<vector<Point<dim> > > > unit_shape_gradients;
+    vector<vector<vector<Tensor<1,dim> > > > unit_shape_gradients;
     
+                                    /**
+                                     * Store the 2nd derivatives of the shape
+                                     * functions at the quadrature points on
+                                     * the unit cell for each face.
+                                     * This field is set up upon construction
+                                     * of the object and contains the
+                                     * derivatives on the reference element.
+                                     */
+    vector<vector<vector<Tensor<2,dim> > > > unit_shape_2nd_derivatives;
+
                                     /**
                                      * Gradients of the basis
                                      * functions of the transformation.
@@ -1087,7 +1147,7 @@ const dFMatrix & FEValuesBase<dim>::get_shape_values () const {
 
 template <int dim>
 inline
-const vector<vector<Point<dim> > > &
+const vector<vector<Tensor<1,dim> > > &
 FEValuesBase<dim>::get_shape_grads () const {
   Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
   return shape_gradients;
index f3379ebea65f98bde50e9f1fab48c487d1b80c92..0833482a4a9b4ca49952e1529b939dd810bdd6ff 100644 (file)
@@ -178,7 +178,7 @@ bool FiniteElementBase<dim>::operator == (const FiniteElementBase<dim> &f) const
 template <>
 void FiniteElement<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
                                       const vector<Point<1> > &unit_points,
-                                      vector<dFMatrix>  &jacobians,
+                                      vector<Tensor<2,1> >    &jacobians,
                                       const bool         compute_jacobians,
                                       vector<Point<1> > &support_points,
                                       const bool         compute_support_points,
@@ -224,7 +224,7 @@ void FiniteElement<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &
                                            const unsigned int       ,
                                            const vector<Point<0> > &,
                                            const vector<Point<1> > &,
-                                           vector<dFMatrix>        &,
+                                           vector<Tensor<2,1> >    &,
                                            const bool               ,
                                            vector<Point<1> >       &,
                                            const bool               ,
@@ -247,7 +247,7 @@ void FiniteElement<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterato
                                               const unsigned int       ,
                                               const vector<Point<0> > &,
                                               const vector<Point<1> > &,
-                                              vector<dFMatrix>        &,
+                                              vector<Tensor<2,1> >    &,
                                               const bool               ,
                                               vector<Point<1> >       &,
                                               const bool               ,
@@ -316,7 +316,7 @@ void FiniteElement<1>::get_support_points (const DoFHandler<1>::cell_iterator &c
 template <int dim>
 void FiniteElement<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &,
                                         const vector<Point<dim> > &,
-                                        vector<dFMatrix> &,
+                                        vector<Tensor<2,dim> > &,
                                         const bool,
                                         vector<Point<dim> > &,
                                         const bool,
@@ -335,7 +335,7 @@ void FiniteElement<dim>::fill_fe_face_values (const DoFHandler<dim>::cell_iterat
                                              const unsigned int           face_no,
                                              const vector<Point<dim-1> > &unit_points,
                                              const vector<Point<dim> > &global_unit_points,
-                                             vector<dFMatrix>    &jacobians,
+                                             vector<Tensor<2,dim> >    &jacobians,
                                              const bool           compute_jacobians,
                                              vector<Point<dim> > &support_points,
                                              const bool           compute_support_points,
@@ -386,7 +386,7 @@ void FiniteElement<dim>::fill_fe_subface_values (const DoFHandler<dim>::cell_ite
                                                 const unsigned int           subface_no,
                                                 const vector<Point<dim-1> > &unit_points,
                                                 const vector<Point<dim> > &global_unit_points,
-                                                vector<dFMatrix>    &jacobians,
+                                                vector<Tensor<2,dim> >    &jacobians,
                                                 const bool           compute_jacobians,
                                                 vector<Point<dim> > &q_points,
                                                 const bool           compute_q_points,
index cb46024d834421f023ceebf7bcc3b3101a155c7e..5cace7e1d7a9cd3b06f1552d7d3b1050b4671baf 100644 (file)
@@ -323,7 +323,7 @@ double FECrissCross<1>::shape_value (const unsigned int, const Point<1> &) const
 
 
 template <>
-Point<1> FECrissCross<1>::shape_grad (const unsigned int, const Point<1> &) const {
+Tensor<1,1> FECrissCross<1>::shape_grad (const unsigned int, const Point<1> &) const {
   Assert (false, ExcNotUseful());
   return Point<1>();
 };
@@ -538,7 +538,7 @@ double FECrissCross<2>::shape_value (const unsigned int i,
 
 template <>
 inline
-Point<2> FECrissCross<2>::shape_grad (const unsigned int i, const Point<2> &p) const {
+Tensor<1,2> FECrissCross<2>::shape_grad (const unsigned int i, const Point<2> &p) const {
   Assert((i<total_dofs), ExcInvalidIndex(i));
 
   const double x = p(0),
@@ -927,7 +927,7 @@ void FECrissCross<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &ce
 template <int dim>
 void FECrissCross<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
                                        const vector<Point<dim> >            &unit_points,
-                                       vector<dFMatrix>    &jacobians,
+                                       vector<Tensor<2,dim> >               &jacobians,
                                        const bool           compute_jacobians,
                                        vector<Point<dim> > &support_points,
                                        const bool,
@@ -1023,16 +1023,16 @@ void FECrissCross<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &ce
                           vertices[0](1)+t28-t22*xi;
        const double t38 = 1/(t24+t36);
 
-       jacobians[point](0,0) = (-vertices[0](1)+vertices[0](1)*xi-
-                                vertices[1](1)*xi+vertices[2](1)*xi+
-                                vertices[3](1)-vertices[3](1)*xi)*t38;
-       jacobians[point](0,1) = -(-vertices[0](0)+vertices[0](0)*xi-
-                                 vertices[1](0)*xi+t8+vertices[3](0)-t16)*t38;
-       jacobians[point](1,0) = -(-vertices[0](1)+vertices[0](1)*eta+
-                                 vertices[1](1)-vertices[1](1)*eta+
-                                 vertices[2](1)*eta-vertices[3](1)*eta)*t38;
-       jacobians[point](1,1) = (-vertices[0](0)+vertices[0](0)*eta+
-                                vertices[1](0)-t10+t31-vertices[3](0)*eta)*t38;
+       jacobians[point][0][0] = (-vertices[0](1)+vertices[0](1)*xi-
+                                 vertices[1](1)*xi+vertices[2](1)*xi+
+                                 vertices[3](1)-vertices[3](1)*xi)*t38;
+       jacobians[point][0][1] = -(-vertices[0](0)+vertices[0](0)*xi-
+                                  vertices[1](0)*xi+t8+vertices[3](0)-t16)*t38;
+       jacobians[point][1][0] = -(-vertices[0](1)+vertices[0](1)*eta+
+                                  vertices[1](1)-vertices[1](1)*eta+
+                                  vertices[2](1)*eta-vertices[3](1)*eta)*t38;
+       jacobians[point][1][1] = (-vertices[0](0)+vertices[0](0)*eta+
+                                 vertices[1](0)-t10+t31-vertices[3](0)*eta)*t38;
       };
 };
 
index 1216637d6343a7c6f9c0ad07e1c86c54b4a868d8..9c62a710d803f6cc329ecf31eba2a8a8dd7758d8 100644 (file)
@@ -107,12 +107,17 @@ FECubicSub<1>::shape_value (const unsigned int i,
 
 
 template <>
-Point<1>
+Tensor<1,1>
 FECubicSub<1>::shape_grad (const unsigned int i,
                           const Point<1>    &p) const
 {
   Assert((i<total_dofs), ExcInvalidIndex(i));
   const double xi = p(0);
+                                  // originally, the return type of the
+                                  // function was Point<dim>, so we
+                                  // still construct it as that. it should
+                                  // make no difference in practice,
+                                  // however
   switch (i)
     {
       case 0: return Point<1>(-27.0/2.0*xi*xi+18.0*xi-11.0/2.0);
@@ -721,7 +726,7 @@ xi*xi)*eta*eta*eta;
 
 
 template <>
-Point<2>
+Tensor<1,2>
 FECubicSub<2>::shape_grad (const unsigned int i,
                               const Point<2>    &p) const
 {
@@ -729,6 +734,11 @@ FECubicSub<2>::shape_grad (const unsigned int i,
 
   const double xi = p(0),
               eta= p(1);
+                                  // originally, the return type of the
+                                  // function was Point<dim>, so we
+                                  // still construct it as that. it should
+                                  // make no difference in practice,
+                                  // however
   switch (i)
     {
       case 0: return Point<2>(-11.0/2.0+18.0*xi-27.0/2.0*xi*xi+(121.0/4.0-99.0*xi+297.0/4.0*xi*xi)*eta+(-99.0/2.0+162.0*xi-243.0/2.0*xi*xi)*eta*eta+(99.0/4.0-81.0*xi+243.0/4.0*xi*xi)*eta*eta*eta,
index 33babbfa5c0c9872120da14736df6007f03a4719..08d1580d4da71648461a23fbbef3b671ef98c7f5 100644 (file)
@@ -62,11 +62,16 @@ FELinear<1>::shape_value(const unsigned int i,
 
 template <>
 inline
-Point<1>
+Tensor<1,1>
 FELinear<1>::shape_grad(const unsigned int i,
                        const Point<1>&) const
 {
   Assert((i<total_dofs), ExcInvalidIndex(i));
+                                  // originally, the return type of the
+                                  // function was Point<dim>, so we
+                                  // still construct it as that. it should
+                                  // make no difference in practice,
+                                  // however
   switch (i)
     {
     case 0: return Point<1>(-1.);
@@ -215,11 +220,16 @@ FELinear<2>::shape_value (const unsigned int i,
 
 template <>
 inline
-Point<2>
+Tensor<1,2>
 FELinear<2>::shape_grad (const unsigned int i,
                         const Point<2>& p) const
 {
   Assert((i<total_dofs), ExcInvalidIndex(i));
+                                  // originally, the return type of the
+                                  // function was Point<dim>, so we
+                                  // still construct it as that. it should
+                                  // make no difference in practice,
+                                  // however
   switch (i)
     {
       case 0: return Point<2> (p(1)-1., p(0)-1.);
index 49236de6d6c4e71e0aa197055ae2ffb3edbc20ab..33c3a4c7e317927dc2c7a98e7043c231c3d6dc4f 100644 (file)
@@ -118,12 +118,17 @@ FEQuadraticSub<1>::shape_value(const unsigned int i,
 
 
 template <>
-Point<1>
+Tensor<1,1>
 FEQuadraticSub<1>::shape_grad(const unsigned int i,
                              const Point<1>    &p) const
 {
   Assert((i<total_dofs), ExcInvalidIndex(i));
   const double xi = p(0);
+                                  // originally, the return type of the
+                                  // function was Point<dim>, so we
+                                  // still construct it as that. it should
+                                  // make no difference in practice,
+                                  // however
   switch (i)
     {
       case 0: return Point<1>(-3+4*xi);
@@ -632,7 +637,7 @@ FEQuadraticSub<2>::shape_value (const unsigned int i,
 
 
 template <>
-Point<2>
+Tensor<1,2>
 FEQuadraticSub<2>::shape_grad (const unsigned int i,
                               const Point<2>    &p) const
 {
@@ -640,6 +645,11 @@ FEQuadraticSub<2>::shape_grad (const unsigned int i,
 
   const double xi = p(0),
               eta= p(1);
+                                  // originally, the return type of the
+                                  // function was Point<dim>, so we
+                                  // still construct it as that. it should
+                                  // make no difference in practice,
+                                  // however
   switch (i)
     {
       case 0: return Point<2>(-(2*xi-1)*(1-eta)*(2*eta-1)+2*(1-xi)*(1-eta)*(2*eta-1),
index 7cd303f6e1bd613dd241fb8ea48239b80f963162..a6cf7ffd3560cf5d9fffa9998ffc40e0df9dd671 100644 (file)
@@ -75,12 +75,17 @@ FEQuarticSub<1>::shape_value(const unsigned int i,
 
 
 template <>
-Point<1>
+Tensor<1,1>
 FEQuarticSub<1>::shape_grad(const unsigned int i,
                            const Point<1>    &p) const
 {
   Assert((i<total_dofs), ExcInvalidIndex(i));
   const double xi = p(0);
+                                  // originally, the return type of the
+                                  // function was Point<dim>, so we
+                                  // still construct it as that. it should
+                                  // make no difference in practice,
+                                  // however
   switch (i)
     {
       case 0: return Point<1>(128.0/3.0*xi*xi*xi-80.0*xi*xi+140.0/3.0*xi-25.0/3.0);
@@ -1089,7 +1094,7 @@ FEQuarticSub<2>::shape_value (const unsigned int i,
 
 
 template <>
-Point<2>
+Tensor<1,2>
 FEQuarticSub<2>::shape_grad (const unsigned int i,
                             const Point<2>    &p) const
 {
@@ -1097,6 +1102,11 @@ FEQuarticSub<2>::shape_grad (const unsigned int i,
 
   const double xi = p(0),
               eta= p(1);
+                                  // originally, the return type of the
+                                  // function was Point<dim>, so we
+                                  // still construct it as that. it should
+                                  // make no difference in practice,
+                                  // however
   switch (i)
     {
       case 0: return Point<2>(-25.0/3.0+140.0/3.0*xi-80.0*xi*xi+128.0/3.0*xi*xi*xi+(625.0/9.0-3500.0/9.0*xi+2000.0/3.0*xi*xi-3200.0/9.0*xi*xi*xi)*eta+(-1750.0/9.0+9800.0/9.0*xi-5600.0/3.0*xi*xi+8960.0/9.0*xi*xi*xi)*eta*eta+(2000.0/9.0-11200.0/9.0*xi+6400.0/3.0*xi*xi-10240.0/9.0*xi*xi*xi)*eta*eta*eta+(-800.0/9.0+4480.0/9.0*xi-2560.0/3.0*xi*xi+4096.0/9.0*xi*xi*xi)*eta*eta*eta*eta,
index a6eada6302a429b889002f58b93767065570e609..e2728ad7635cfc73b8d452a7fffffb19a15ca7c9 100644 (file)
@@ -279,7 +279,7 @@ void FELinearMapping<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator
 template <int dim>
 void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
                                           const vector<Point<dim> >            &unit_points,
-                                          vector<dFMatrix>    &jacobians,
+                                          vector<Tensor<2,dim> >               &jacobians,
                                           const bool           compute_jacobians,
                                           vector<Point<dim> > &support_points,
                                           const bool           compute_support_points,
@@ -429,40 +429,53 @@ void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator
 
   if (compute_jacobians)
     for (unsigned int point=0; point<n_points; ++point)
-      {
-       const double xi = unit_points[point](0);
-       const double eta= unit_points[point](1);
+      switch (dim)
+       {
+         case 2:
+         {
+           
+           const double xi = unit_points[point](0);
+           const double eta= unit_points[point](1);
        
-       const double t6 = vertices[0](0)*vertices[3](1);
-       const double t8 = vertices[2](0)*xi;
-       const double t10 = vertices[1](0)*eta;
-       const double t12 = vertices[3](0)*vertices[1](1);
-       const double t16 = vertices[3](0)*xi;
-       const double t20 = vertices[0](0)*vertices[1](1);
-       const double t22 = vertices[0](0)*vertices[2](1);
-       const double t24 = t6*xi-t8*vertices[1](1)-t10*vertices[3](1)+
-                          t12*eta-vertices[3](0)*vertices[2](1)*eta-
-                          t16*vertices[0](1)+t16*vertices[1](1)-t12+
-                          vertices[3](0)*vertices[0](1)-t20*eta+t22*eta;
-       const double t28 = vertices[1](0)*vertices[3](1);
-       const double t31 = vertices[2](0)*eta;
-       const double t36 = t8*vertices[0](1)+vertices[1](0)*vertices[2](1)*xi-
-                          t28*xi+t10*vertices[0](1)-t31*vertices[0](1)+
-                          t31*vertices[3](1)+t20-t6-vertices[1](0)*
-                          vertices[0](1)+t28-t22*xi;
-       const double t38 = 1/(t24+t36);
-
-       jacobians[point](0,0) = (-vertices[0](1)+vertices[0](1)*xi-
-                                vertices[1](1)*xi+vertices[2](1)*xi+
-                                vertices[3](1)-vertices[3](1)*xi)*t38;
-       jacobians[point](0,1) = -(-vertices[0](0)+vertices[0](0)*xi-
-                                 vertices[1](0)*xi+t8+vertices[3](0)-t16)*t38;
-       jacobians[point](1,0) = -(-vertices[0](1)+vertices[0](1)*eta+
-                                 vertices[1](1)-vertices[1](1)*eta+
-                                 vertices[2](1)*eta-vertices[3](1)*eta)*t38;
-       jacobians[point](1,1) = (-vertices[0](0)+vertices[0](0)*eta+
-                                vertices[1](0)-t10+t31-vertices[3](0)*eta)*t38;
-      };
+           const double t6 = vertices[0](0)*vertices[3](1);
+           const double t8 = vertices[2](0)*xi;
+           const double t10 = vertices[1](0)*eta;
+           const double t12 = vertices[3](0)*vertices[1](1);
+           const double t16 = vertices[3](0)*xi;
+           const double t20 = vertices[0](0)*vertices[1](1);
+           const double t22 = vertices[0](0)*vertices[2](1);
+           const double t24 = t6*xi-t8*vertices[1](1)-t10*vertices[3](1)+
+                              t12*eta-vertices[3](0)*vertices[2](1)*eta-
+                              t16*vertices[0](1)+t16*vertices[1](1)-t12+
+                              vertices[3](0)*vertices[0](1)-t20*eta+t22*eta;
+           const double t28 = vertices[1](0)*vertices[3](1);
+           const double t31 = vertices[2](0)*eta;
+           const double t36 = t8*vertices[0](1)+vertices[1](0)*vertices[2](1)*xi-
+                              t28*xi+t10*vertices[0](1)-t31*vertices[0](1)+
+                              t31*vertices[3](1)+t20-t6-vertices[1](0)*
+                              vertices[0](1)+t28-t22*xi;
+           const double t38 = 1/(t24+t36);
+
+           jacobians[point][0][0] = (-vertices[0](1)+vertices[0](1)*xi-
+                                     vertices[1](1)*xi+vertices[2](1)*xi+
+                                     vertices[3](1)-vertices[3](1)*xi)*t38;
+           jacobians[point][0][1] = -(-vertices[0](0)+vertices[0](0)*xi-
+                                      vertices[1](0)*xi+t8+vertices[3](0)-t16)*t38;
+           jacobians[point][1][0] = -(-vertices[0](1)+vertices[0](1)*eta+
+                                      vertices[1](1)-vertices[1](1)*eta+
+                                      vertices[2](1)*eta-vertices[3](1)*eta)*t38;
+           jacobians[point][1][1] = (-vertices[0](0)+vertices[0](0)*eta+
+                                     vertices[1](0)-t10+t31-vertices[3](0)*eta)*t38;
+           
+           break;
+         };
+
+         default:
+                                                // not implemented at present,
+                                                // because of the changes above
+               Assert (false, ExcNotImplemented());
+       };
+  
   
     
   if (compute_support_points)
index cec0e6e9ad9c678951184cc049a89f2d1bbc14d9..b2ec9d3e8ce4fdb8fc00aba2095be1a8104f74d7 100644 (file)
@@ -27,12 +27,13 @@ FEValuesBase<dim>::FEValuesBase (const unsigned int n_q_points,
                total_dofs (n_dofs),
                n_transform_functions (n_transform_functions),
                shape_values (n_values_arrays, dFMatrix(n_dofs, n_q_points)),
-               shape_gradients (n_dofs, vector<Point<dim> >(n_q_points)),
+               shape_gradients (n_dofs, vector<Tensor<1,dim> >(n_q_points)),
                weights (n_q_points, 0),
                JxW_values (n_q_points, 0),
                quadrature_points (n_q_points, Point<dim>()),
                support_points (n_support_points, Point<dim>()),
-               jacobi_matrices (n_q_points, dFMatrix(dim,dim)),
+               jacobi_matrices (n_q_points, Tensor<2,dim>()),
+               jacobi_matrices_grad (n_q_points, Tensor<3,dim>()),
                shape_values_transform (n_values_arrays,
                                        dFMatrix(n_transform_functions,
                                                 n_quadrature_points)),
@@ -85,7 +86,7 @@ void FEValuesBase<dim>::get_function_values (const dVector  &fe_function,
 
 
 template <int dim>
-const Point<dim> &
+const Tensor<1,dim> &
 FEValuesBase<dim>::shape_grad (const unsigned int i,
                               const unsigned int j) const {
   Assert (i<shape_gradients.size(),
@@ -117,8 +118,11 @@ void FEValuesBase<dim>::get_function_grads (const dVector       &fe_function,
                                   // functions
   for (unsigned int point=0; point<n_quadrature_points; ++point)
     for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
-      gradients[point] += (dof_values(shape_func) *
-                          shape_gradients[shape_func][point]);
+      {
+       Tensor<1,dim> tmp(shape_gradients[shape_func][point]);
+       tmp *= dof_values(shape_func);
+       gradients[point] += tmp;
+      };
 };
 
 
@@ -168,7 +172,9 @@ FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
                                   update_flags,
                                   fe),
                unit_shape_gradients(fe.total_dofs,
-                                    vector<Point<dim> >(quadrature.n_quadrature_points)),
+                                    vector<Tensor<1,dim> >(quadrature.n_quadrature_points)),
+               unit_shape_2nd_derivatives(fe.total_dofs,
+                                          vector<Tensor<2,dim> >(quadrature.n_quadrature_points)),
                unit_shape_gradients_transform(fe.n_transform_functions,
                                               vector<Point<dim> >(quadrature.n_quadrature_points)),
                unit_quadrature_points(quadrature.get_quad_points())
@@ -182,6 +188,8 @@ FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
        shape_values[0](i,j) = fe.shape_value(i, unit_quadrature_points[j]);
        unit_shape_gradients[i][j]
          = fe.shape_grad(i, unit_quadrature_points[j]);
+       unit_shape_2nd_derivatives[i][j]
+         = fe.shape_grad_grad(i, unit_quadrature_points[j]);
       };
 
   for (unsigned int i=0; i<n_transform_functions; ++i)
@@ -218,7 +226,10 @@ void FEValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell,
                        jacobi_matrices,
                        update_flags & (update_jacobians  |
                                        update_JxW_values |
-                                       update_gradients),
+                                       update_gradients  |
+                                       update_second_derivatives),
+//                     jacobi_matrices_grad,
+//                     update_flags & update_second_derivatives,
                        support_points,
                        update_flags & update_support_points,
                        quadrature_points,
@@ -233,16 +244,41 @@ void FEValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell,
       for (unsigned int j=0; j<n_quadrature_points; ++j)
        for (unsigned int s=0; s<dim; ++s)
          {
-           shape_gradients[i][j](s) = 0;
+           shape_gradients[i][j][s] = 0;
            
                                             // (grad psi)_s =
                                             // (grad_{\xi\eta})_b J_{bs}
                                             // with J_{bs}=(d\xi_b)/(dx_s)
            for (unsigned int b=0; b<dim; ++b)
-             shape_gradients[i][j](s)
+             shape_gradients[i][j][s]
                +=
-               unit_shape_gradients[i][j](b) * jacobi_matrices[j](b,s);
+               unit_shape_gradients[i][j][b] * jacobi_matrices[j][b][s];
          };
+
+  
+  if (update_flags & update_second_derivatives)
+    for (unsigned int i=0; i<fe->total_dofs; ++i)
+      for (unsigned int j=0; j<n_quadrature_points; ++j)
+       {
+         static Tensor<2,dim> tmp1, tmp2;
+
+                                          // tmp1 := (d_k d_l phi) J_lj
+         contract (tmp1, unit_shape_2nd_derivatives[i][j], jacobi_matrices[j]);
+                                          // tmp2 := tmp1_kj J_ki
+         contract (tmp2, tmp1, 1, jacobi_matrices[j], 1);
+
+
+                                          // second part:
+                                          // tmp1 := (d_k J_lj) (d_l phi)
+         contract (tmp1, jacobi_matrices_grad[j], 2, unit_shape_gradients[i][j]);
+                                          // tmp1_kj J_ki
+         contract (shape_2nd_derivatives[i][j],
+                   jacobi_matrices[j], 1,
+                   tmp1, 1);
+
+                                          // add up first contribution
+         shape_2nd_derivatives[i][j] += tmp2;
+       };
   
   
                                   // compute Jacobi determinants in
@@ -252,7 +288,7 @@ void FEValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell,
                                   // determinant
   if (update_flags & update_JxW_values) 
     for (unsigned int i=0; i<n_quadrature_points; ++i)
-      JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
+      JxW_values[i] = weights[i] / determinant(jacobi_matrices[i]);
 };
 
 
@@ -278,8 +314,11 @@ FEFaceValuesBase<dim>::FEFaceValuesBase (const unsigned int n_q_points,
                                   update_flags,
                                   fe),
                unit_shape_gradients (n_faces_or_subfaces,
-                                     vector<vector<Point<dim> > >(n_dofs,
-                                                                  vector<Point<dim> >(n_q_points))),
+                                     vector<vector<Tensor<1,dim> > >(n_dofs,
+                                                                  vector<Tensor<1,dim> >(n_q_points))),
+               unit_shape_2nd_derivatives(n_dofs,
+                                          vector<vector<Tensor<2,dim> > >(n_dofs,
+                                                                          vector<Tensor<2,dim> >(n_q_points))),
                unit_shape_gradients_transform (n_faces_or_subfaces,
                                                vector<vector<Point<dim> > >(n_transform_functions,
                                                                             vector<Point<dim> >(n_q_points))),
@@ -340,6 +379,8 @@ FEFaceValues<dim>::FEFaceValues (const FiniteElement<dim> &fe,
            = fe.shape_value(i, unit_quadrature_points[face][j]);
          unit_shape_gradients[face][i][j]
            = fe.shape_grad(i, unit_quadrature_points[face][j]);
+         unit_shape_2nd_derivatives[face][i][j]
+           = fe.shape_grad_grad(i, unit_quadrature_points[face][j]);
        };
 
   for (unsigned int i=0; i<n_transform_functions; ++i)
@@ -403,10 +444,37 @@ void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &c
                                             // (grad_{\xi\eta})_b J_{bs}
                                             // with J_{bs}=(d\xi_b)/(dx_s)
            for (unsigned int b=0; b<dim; ++b)
-             shape_gradients[i][j](s)
-               += (unit_shape_gradients[face_no][i][j](b) *
-                   jacobi_matrices[j](b,s));
+             shape_gradients[i][j][s]
+               += (unit_shape_gradients[face_no][i][j][b] *
+                   jacobi_matrices[j][b][s]);
       };
+
+
+  if (update_flags & update_second_derivatives)
+    for (unsigned int i=0; i<fe->total_dofs; ++i)
+      for (unsigned int j=0; j<n_quadrature_points; ++j)
+       {
+         static Tensor<2,dim> tmp1, tmp2;
+
+                                          // tmp1 := (d_k d_l phi) J_lj
+         contract (tmp1, unit_shape_2nd_derivatives[face_no][i][j], jacobi_matrices[j]);
+                                          // tmp2 := tmp1_kj J_ki
+         contract (tmp2, tmp1, 1, jacobi_matrices[j], 1);
+
+
+                                          // second part:
+                                          // tmp1 := (d_k J_lj) (d_l phi)
+         contract (tmp1,
+                   jacobi_matrices_grad[j], 2,
+                   unit_shape_gradients[face_no][i][j]);
+                                          // tmp1_kj J_ki
+         contract (shape_2nd_derivatives[i][j],
+                   jacobi_matrices[j], 1,
+                   tmp1, 1);
+
+                                          // add up first contribution
+         shape_2nd_derivatives[i][j] += tmp2;
+       };
   
   
                                   // compute Jacobi determinants in
@@ -424,7 +492,7 @@ void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &c
 
 
 
-/*------------------------------- FEFaceValues -------------------------------*/
+/*------------------------------- FESubFaceValues -------------------------------*/
 
 
 template <int dim>
@@ -469,6 +537,10 @@ FESubfaceValues<dim>::FESubfaceValues (const FiniteElement<dim> &fe,
              = fe.shape_grad(i, unit_quadrature_points[face *
                                                       GeometryInfo<dim>::
                                                       subfaces_per_face+subface][j]);
+           unit_shape_2nd_derivatives[face*GeometryInfo<dim>::subfaces_per_face+subface][i][j]
+             = fe.shape_grad_grad(i, unit_quadrature_points[face *
+                                                           GeometryInfo<dim>::
+                                                           subfaces_per_face+subface][j]);
          };
   for (unsigned int i=0; i<n_transform_functions; ++i)
     for (unsigned int j=0; j<n_quadrature_points; ++j)
@@ -542,10 +614,38 @@ void FESubfaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator
                                             // (grad_{\xi\eta})_b J_{bs}
                                             // with J_{bs}=(d\xi_b)/(dx_s)
            for (unsigned int b=0; b<dim; ++b)
-             shape_gradients[i][j](s)
-               += (unit_shape_gradients[selected_dataset][i][j](b) *
-                   jacobi_matrices[j](b,s));
+             shape_gradients[i][j][s]
+               += (unit_shape_gradients[selected_dataset][i][j][b] *
+                   jacobi_matrices[j][b][s]);
       };
+
+  if (update_flags & update_second_derivatives)
+    for (unsigned int i=0; i<fe->total_dofs; ++i)
+      for (unsigned int j=0; j<n_quadrature_points; ++j)
+       {
+         static Tensor<2,dim> tmp1, tmp2;
+
+                                          // tmp1 := (d_k d_l phi) J_lj
+         contract (tmp1,
+                   unit_shape_2nd_derivatives[selected_dataset][i][j],
+                   jacobi_matrices[j]);
+                                          // tmp2 := tmp1_kj J_ki
+         contract (tmp2, tmp1, 1, jacobi_matrices[j], 1);
+
+
+                                          // second part:
+                                          // tmp1 := (d_k J_lj) (d_l phi)
+         contract (tmp1,
+                   jacobi_matrices_grad[j], 2,
+                   unit_shape_gradients[selected_dataset][i][j]);
+                                          // tmp1_kj J_ki
+         contract (shape_2nd_derivatives[i][j],
+                   jacobi_matrices[j], 1,
+                   tmp1, 1);
+
+                                          // add up first contribution
+         shape_2nd_derivatives[i][j] += tmp2;
+       };
   
   
                                   // compute Jacobi determinants in
index 19e5d6cc7bca9b5c7c64bfde67064123bc52c823..9075fc1e9fe52efc1c0ce63bd5c3daa687645226 100644 (file)
@@ -694,7 +694,7 @@ void LaplaceMatrix<dim>::assemble (dFMatrix            &cell_matrix,
   Assert (cell_matrix.all_zero(), ExcObjectNotEmpty());
   Assert (rhs.all_zero(), ExcObjectNotEmpty());
 
-  const vector<vector<Point<dim> > >&gradients = fe_values.get_shape_grads ();
+  const vector<vector<Tensor<1,dim> > >&gradients = fe_values.get_shape_grads ();
   const dFMatrix       &values    = fe_values.get_shape_values ();
   vector<double>        rhs_values(fe_values.n_quadrature_points);
   const vector<double> &weights   = fe_values.get_JxW_values ();
@@ -748,7 +748,7 @@ void LaplaceMatrix<dim>::assemble (dFMatrix            &cell_matrix,
          ExcWrongSize(cell_matrix.m(), total_dofs));
   Assert (cell_matrix.all_zero(), ExcObjectNotEmpty());
   
-  const vector<vector<Point<dim> > >&gradients = fe_values.get_shape_grads ();
+  const vector<vector<Tensor<1,dim> > >&gradients = fe_values.get_shape_grads ();
   const vector<double> &weights   = fe_values.get_JxW_values ();
    
   if (coefficient != 0)
index 7844de59d5ab6951eb18bdf71b790c3ca91c2b9d..e30bfefa7677925b7c708998665549615a707194 100644 (file)
@@ -36,7 +36,7 @@ void PoissonEquation<2>::assemble (dFMatrix            &cell_matrix,
                                   dVector             &rhs,
                                   const FEValues<2>   &fe_values,
                                   const Triangulation<2>::cell_iterator &) const {
-  const vector<vector<Point<2> > >&gradients = fe_values.get_shape_grads ();
+  const vector<vector<Tensor<1,2> > >&gradients = fe_values.get_shape_grads ();
   const dFMatrix       &values    = fe_values.get_shape_values ();
   vector<double>        rhs_values (fe_values.n_quadrature_points);
   const vector<double> &weights   = fe_values.get_JxW_values ();

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.