]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Work on the implementation of second derivatives.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 8 Nov 1998 20:36:10 +0000 (20:36 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 8 Nov 1998 20:36:10 +0000 (20:36 +0000)
git-svn-id: https://svn.dealii.org/trunk@657 0785d39b-7218-0410-832d-ea1e28bc413d

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_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_linear_mapping.cc
deal.II/deal.II/source/fe/fe_values.cc
deal.II/deal.II/source/fe/scripts/2d/postprocess

index 69a3bc511d5f3945ad93c18c35f5aa6c015715ea..836f4c892a436c2ead7dd8f285f0e23510d496ad 100644 (file)
@@ -589,12 +589,15 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      * and the given quadrature points on the
                                      * unit cell. The Jacobian matrix is to
                                      * be computed at every quadrature point.
+                                     * The derivative of the jacobian matrix
+                                     * is the derivative with respect to the
+                                     * unit cell coordinates.
                                      * This function has to be in the finite
                                      * element class, since different finite
                                      * elements need different transformations
                                      * of the unit cell to a real cell.
                                      *
-                                     * The computation of the three fields may
+                                     * The computation of these fields may
                                      * share some common code, which is why we
                                      * put it in one function. However, it may
                                      * not always be necessary to really
@@ -607,25 +610,33 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      * of the Jacobi matrix and of the various
                                      * structures to be filled.
                                      *
-                                     * It is provided for the finite element
-                                     * class in one space dimension, but for
-                                     * higher dimensions, it depends on the
-                                     * present fe and needs reimplementation
-                                     * by the user. This is due to the fact
+                                     * This function is provided for
+                                     * the finite element class in
+                                     * one space dimension, but for
+                                     * higher dimensions, it depends
+                                     * on the present fe and needs
+                                     * reimplementation by the
+                                     * user. This is due to the fact
                                      * that the user may want to use
-                                     * iso- or subparametric mappings of the
-                                     * unit cell to the real cell, which
-                                     * makes things much more complicated.
+                                     * iso- or subparametric mappings
+                                     * of the unit cell to the real
+                                     * cell, which makes things much
+                                     * more complicated.
                                      *
-                                     * The #shape_values/grads_transform#
-                                     * arrays store the values and gradients
-                                     * of the transformation basis functions.
-                                     * While this information is not necessary
-                                     * for the computation of the other fields,
-                                     * it allows for significant speedups, since
-                                     * the values and gradients of the transform
-                                     * functions at the quadrature points
-                                     * need not be recomputed each time this
+                                     * The
+                                     * #shape_values/grads_transform#
+                                     * arrays store the values and
+                                     * gradients of the
+                                     * transformation basis
+                                     * functions.  While this
+                                     * information is not necessary
+                                     * for the computation of the
+                                     * other fields, it allows for
+                                     * significant speedups, since
+                                     * the values and gradients of
+                                     * the transform functions at the
+                                     * quadrature points need not be
+                                     * recomputed each time this
                                      * function is called.
                                      *
                                      * The function assumes that the fields
@@ -635,17 +646,18 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      * This function is more or less an
                                      * interface to the #FEValues# class and
                                      * should not be used by users unless
-                                     * absolutely needed.
-                                     */
+                                     * absolutely needed.  */
     virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
-                                const vector<Point<dim> >            &unit_points,
-                                vector<Tensor<2,dim> >               &jacobians,
-                                const bool           compute_jacobians,
-                                vector<Point<dim> > &support_points,
-                                const bool           compute_support_points,
-                                vector<Point<dim> > &q_points,
-                                const bool           compute_q_points,
-                                const dFMatrix      &shape_values_transform,
+                                const vector<Point<dim> > &unit_points,
+                                vector<Tensor<2,dim> >    &jacobians,
+                                const bool                 compute_jacobians,
+                                vector<Tensor<3,dim> >    &jacobians_grad,
+                                const bool                 compute_jacobians_grad,
+                                vector<Point<dim> >       &support_points,
+                                const bool                 compute_support_points,
+                                vector<Point<dim> >       &q_points,
+                                const bool                 compute_q_points,
+                                const dFMatrix            &shape_values_transform,
                                 const vector<vector<Tensor<1,dim> > > &shape_grads_transform,
                                 const Boundary<dim> &boundary) const;
 
@@ -758,7 +770,9 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      const vector<Point<dim-1> > &unit_points,
                                      const vector<Point<dim> >   &global_unit_points,
                                      vector<Tensor<2,dim> >      &jacobians,
-                                     const bool           compute_jacobians,
+                                     const bool                   compute_jacobians,
+                                     vector<Tensor<3,dim> >      &jacobians_grad,
+                                     const bool                   compute_jacobians_grad,
                                      vector<Point<dim> > &support_points,
                                      const bool           compute_support_points,
                                      vector<Point<dim> > &q_points,
@@ -805,7 +819,9 @@ class FiniteElement : public FiniteElementBase<dim> {
                                         const vector<Point<dim-1> > &unit_points,
                                         const vector<Point<dim> >   &global_unit_points,
                                         vector<Tensor<2,dim> >      &jacobians,
-                                        const bool           compute_jacobians,
+                                        const bool                   compute_jacobians,
+                                        vector<Tensor<3,dim> >      &jacobians_grad,
+                                        const bool           compute_jacobians_grad,
                                         vector<Point<dim> > &q_points,
                                         const bool           compute_q_points,
                                         vector<double>      &face_jacobi_determinants,
index 68d41b8ea3a164dcdfba415d2d11cf8334cb6808..6e7d2c493195be41918fed3ef8d20d1c8baf8666 100644 (file)
@@ -333,12 +333,14 @@ class FECrissCross : public FiniteElement<dim> {
     virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
                                 const vector<Point<dim> >            &unit_points,
                                 vector<Tensor<2,dim> >               &jacobians,
-                                const bool           compute_jacobians,
-                                vector<Point<dim> > &support_points,
-                                const bool           compute_support_points,
-                                vector<Point<dim> > &q_points,
-                                const bool           compute_q_points,
-                                const dFMatrix      &shape_values_transform,
+                                const bool              compute_jacobians,
+                                vector<Tensor<3,dim> > &jacobians_grad,
+                                const bool              compute_jacobians_grad,
+                                vector<Point<dim> >    &support_points,
+                                const bool              compute_support_points,
+                                vector<Point<dim> >    &q_points,
+                                const bool              compute_q_points,
+                                const dFMatrix         &shape_values_transform,
                                 const vector<vector<Tensor<1,dim> > > &shape_grad_transform,
                                 const Boundary<dim> &boundary) const;
 
index 9be39b199fb1775397a083cd1ed134acd2dc9c9b..0b0af49d049eceeb2b66e317af20c18e1214ed98 100644 (file)
@@ -136,12 +136,14 @@ class FELinearMapping : public FiniteElement<dim> {
     virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
                                 const vector<Point<dim> >            &unit_points,
                                 vector<Tensor<2,dim> >               &jacobians,
-                                const bool           compute_jacobians,
-                                vector<Point<dim> > &support_points,
-                                const bool           compute_support_points,
-                                vector<Point<dim> > &q_points,
-                                const bool           compute_q_points,
-                                const dFMatrix      &shape_values_transform,
+                                const bool              compute_jacobians,
+                                vector<Tensor<3,dim> > &jacobians_grad,
+                                const bool              compute_jacobians_grad,
+                                vector<Point<dim> >    &support_points,
+                                const bool              compute_support_points,
+                                vector<Point<dim> >    &q_points,
+                                const bool              compute_q_points,
+                                const dFMatrix         &shape_values_transform,
                                 const vector<vector<Tensor<1,dim> > > &shape_grad_transform,
                                 const Boundary<dim> &boundary) const;
 };
index 9e61b96271680773ab7af6e8129a485af8df10ff..793211f6992445f60b81300a2ed8d5132a3716ca 100644 (file)
@@ -71,7 +71,7 @@ enum UpdateFlags {
                                        * Update the second derivatives of the
                                        * shape functions on the real cell.
                                        */
-      update_second_derivatives
+      update_second_derivatives = 64
 };
 
 
index 303c126a246261c32f9d842263ffbda63e2d66a9..35c15b2f76a3a1a234d41aa5f74da4428516006b 100644 (file)
@@ -148,9 +148,9 @@ template <int dim> class Quadrature;
  *    a reference to a whole field. Usually these fields contain
  *    the values of all trial functions at all quadrature points.
  *
- *  \item #get_function_values#, #get_function_gradients#: these
- *    two functions offer a simple way to avoid the detour of the
- *    trial functions, if you have a finite solution (resp. the
+ *  \item #get_function_values#, #get_function_grads#, #...#: these
+ *    functions offer a simple way to avoid the detour of the
+ *    trial functions, if you have a finite element solution (resp. the
  *    vector of values associated with the different trial functions.)
  *    Then you may want to get information from the restriction of
  *    the finite element function to a certain cell, e.g. the values
@@ -159,7 +159,8 @@ template <int dim> class Quadrature;
  *    you pass it a vector holding the finite element solution and the
  *    functions return the values or gradients of the finite element
  *    function restricted to the cell which was given last time the
- *    #reinit# function was given.
+ *    #reinit# function was given. The same applies for the functions
+ *    returning higher derivatives of the solution.
  *   
  *    Though possible in principle, these functions do not call the
  *    #reinit# function, you have to do so yourself beforehand. On the
@@ -316,7 +317,7 @@ class FEValuesBase {
                                      * documentation for the matrix itself.
                                      */
     const vector<vector<Tensor<1,dim> > > & get_shape_grads () const;
-    
+
                                     /**
                                      * Return the gradients of the finite
                                      * element function characterized
@@ -330,6 +331,53 @@ class FEValuesBase {
     void get_function_grads (const dVector          &fe_function,
                             vector<Tensor<1,dim> > &gradients) const;
 
+                                    /**
+                                     * Return the 2nd derivatives of
+                                     * the #i#th shape function at
+                                     * the #j# quadrature point. If
+                                     * you want to get the derivatives
+                                     * in one of the coordinate
+                                     * directions, use the
+                                     * appropriate function of the
+                                     * #Tensor# class to extract one
+                                     * component. Since only a
+                                     * reference to the derivatives'
+                                     * values is returned, there
+                                     * should be no major performance
+                                     * drawback.  The function
+                                     * returns the derivatives on the
+                                     * real element, not the
+                                     * reference element.
+                                     */
+    const Tensor<2,dim> & shape_2nd_derivative (const unsigned int function,
+                                               const unsigned int quadrature_point) const;
+
+                                    /**
+                                     * Return a pointer to the
+                                     * matrix holding all 2nd
+                                     * derivatives of shape functions
+                                     * at all integration points, on
+                                     * the present cell.  For the
+                                     * format of this matrix, see the
+                                     * documentation for the matrix
+                                     * itself.
+                                     */
+    const vector<vector<Tensor<2,dim> > > & get_shape_2nd_derivatives () const;
+    
+                                    /**
+                                     * Return the tensor of second
+                                     * derivatives of the finite
+                                     * element function characterized
+                                     * by #fe_function# restricted to
+                                     * #cell# at the quadrature points.
+                                     *
+                                     * The function assumes that the
+                                     * #second_derivatives# object already has
+                                     * the right size.
+                                     */
+    void get_function_2nd_derivatives (const dVector          &fe_function,
+                                      vector<Tensor<2,dim> > &second_derivatives) const;
+
                                     /**
                                      * Return the position of the #i#th
                                      * quadrature point in real space.
@@ -1155,6 +1203,16 @@ FEValuesBase<dim>::get_shape_grads () const {
 
 
 
+template <int dim>
+inline
+const vector<vector<Tensor<2,dim> > > &
+FEValuesBase<dim>::get_shape_2nd_derivatives () const {
+  Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField());
+  return shape_2nd_derivatives;
+};
+
+
+
 template <int dim>
 inline
 const vector<Point<dim> > &
index 57151da0e95cdac968f6ceeb19ae4ce5eaf8644a..7fe883f092edd593e8f4736bb554ffdd5e417b91 100644 (file)
@@ -318,6 +318,8 @@ void FiniteElement<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &,
                                         const vector<Point<dim> > &,
                                         vector<Tensor<2,dim> > &,
                                         const bool,
+                                        vector<Tensor<3,dim> > &,
+                                        const bool,
                                         vector<Point<dim> > &,
                                         const bool,
                                         vector<Point<dim> > &,
@@ -337,6 +339,8 @@ void FiniteElement<dim>::fill_fe_face_values (const DoFHandler<dim>::cell_iterat
                                              const vector<Point<dim> > &global_unit_points,
                                              vector<Tensor<2,dim> >    &jacobians,
                                              const bool           compute_jacobians,
+                                             vector<Tensor<3,dim> >    &jacobians_grad,
+                                             const bool           compute_jacobians_grad,
                                              vector<Point<dim> > &support_points,
                                              const bool           compute_support_points,
                                              vector<Point<dim> > &q_points,
@@ -360,6 +364,7 @@ void FiniteElement<dim>::fill_fe_face_values (const DoFHandler<dim>::cell_iterat
   vector<Point<dim> > dummy(total_dofs);
   fill_fe_values (cell, global_unit_points,
                  jacobians, compute_jacobians,
+                 jacobians_grad, compute_jacobians_grad,
                  dummy, false,
                  q_points, compute_q_points,
                  shape_values_transform, shape_gradients_transform,
@@ -388,6 +393,8 @@ void FiniteElement<dim>::fill_fe_subface_values (const DoFHandler<dim>::cell_ite
                                                 const vector<Point<dim> > &global_unit_points,
                                                 vector<Tensor<2,dim> >    &jacobians,
                                                 const bool           compute_jacobians,
+                                                vector<Tensor<3,dim> >    &jacobians_grad,
+                                                const bool           compute_jacobians_grad,
                                                 vector<Point<dim> > &q_points,
                                                 const bool           compute_q_points,
                                                 vector<double>      &face_jacobi_determinants,
@@ -407,6 +414,7 @@ void FiniteElement<dim>::fill_fe_subface_values (const DoFHandler<dim>::cell_ite
   vector<Point<dim> > dummy(total_dofs);
   fill_fe_values (cell, global_unit_points,
                  jacobians, compute_jacobians,
+                 jacobians_grad, compute_jacobians_grad,
                  dummy, false,
                  q_points, compute_q_points,
                  shape_values_transform, shape_gradients_transform,
index 0af836516f4eae63c514508e96ff48ef52eda336..de29c860edc3dd913d9cf9b7213a07f38aba5d13 100644 (file)
@@ -928,12 +928,14 @@ template <int dim>
 void FECrissCross<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
                                        const vector<Point<dim> >            &unit_points,
                                        vector<Tensor<2,dim> >               &jacobians,
-                                       const bool           compute_jacobians,
-                                       vector<Point<dim> > &support_points,
+                                       const bool              compute_jacobians,
+                                       vector<Tensor<3,dim> > &jacobians_grad,
+                                       const bool              compute_jacobians_grad,
+                                       vector<Point<dim> >    &support_points,
                                        const bool,
-                                       vector<Point<dim> > &q_points,
-                                       const bool           compute_q_points,
-                                       const dFMatrix      &shape_values_transform,
+                                       vector<Point<dim> >    &q_points,
+                                       const bool              compute_q_points,
+                                       const dFMatrix         &shape_values_transform,
                                        const vector<vector<Tensor<1,dim> > > &/*shape_grad_transform*/,
                                        const Boundary<dim> &boundary) const {
   Assert (jacobians.size() == unit_points.size(),
@@ -1051,6 +1053,92 @@ void FECrissCross<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &ce
                                               // not implemented at present
              Assert (false, ExcNotImplemented());
       };
+
+
+  if (compute_jacobians_grad)
+    switch (dim) 
+      {
+       case 1:
+       {
+                                          // derivative of the
+                                          // jacobian is always zero
+                                          // for a linear mapping in 1d
+         for (unsigned int point=0; point<n_points; ++point)
+           jacobians_grad[point][0][0][0] = 0;
+         break;
+       };
+
+       case 2:
+       {
+         for (unsigned int point=0; point<n_points; ++point)
+           {
+             const double xi = unit_points[point](0);
+             const double eta= unit_points[point](1);
+       
+             const double t2 = vertices[1](0)*eta;
+             const double t4 = vertices[3](0)*vertices[2](1);
+             const double t6 = vertices[0](0)*vertices[2](1);
+             const double t8 = vertices[0](0)*vertices[3](1);
+             const double t10 = vertices[3](0)*xi;
+             const double t13 = vertices[2](0)*xi;
+             const double t16 = vertices[3](0)*vertices[1](1);
+             const double t18 = vertices[0](0)*vertices[1](1);
+             const double t19 = vertices[3](0)*vertices[0](1);
+             const double t20 = -t2*vertices[3](1)-t4*eta-t6*xi+t8*xi-
+                                t10*vertices[0](1)+t10*vertices[1](1)+
+                                t13*vertices[0](1)-t13*vertices[1](1)+t16
+                                *eta+t18+t19;
+             const double t23 = vertices[1](0)*vertices[3](1);
+             const double t26 = vertices[2](0)*eta;
+             const double t29 = vertices[1](0)*vertices[0](1);
+             const double t30 = vertices[1](0)*vertices[2](1);
+             const double t32 = -t16-t18*eta+t6*eta-t23*xi+t2*vertices[0](1)-
+                                t26*vertices[0](1)+t26*vertices[3](1)-
+                                t8-t29+t23+t30
+                                *xi;
+             const double t33 = t20+t32;
+             const double t34 = 1/t33;
+             const double t35 = (vertices[0](1)-vertices[1](1)+
+                                 vertices[2](1)-vertices[3](1))*t34;
+             const double t41 = t33*t33;
+             const double t42 = 1/t41;
+             const double t43 = (-vertices[0](1)+vertices[0](1)*xi-
+                                 vertices[1](1)*xi+vertices[2](1)*xi+
+                                 vertices[3](1)-vertices[3](1)*xi)*t42;
+             const double t44 = vertices[2](0)*vertices[0](1);
+             const double t46 = -t6+t8-t19+t16+t44-
+                                vertices[2](0)*vertices[1](1)-
+                                t23+t30;
+             const double t50 = (vertices[0](0)-vertices[1](0)+
+                                 vertices[2](0)-vertices[3](0))*t34;
+             const double t54 = (-vertices[0](0)+vertices[0](0)*xi-
+                                 vertices[1](0)*xi+t13+
+                                 vertices[3](0)-t10)*t42;
+             const double t62 = (-vertices[0](1)+vertices[0](1)*eta+
+                                 vertices[1](1)-vertices[1](1)*eta+
+                                 vertices[2](1)*eta-
+                                 vertices[3](1)*eta)*t42;
+             const double t67 = (-vertices[0](0)+vertices[0](0)*eta+
+                                 vertices[1](0)-t2+t26-vertices[3](0)*eta)*t42;
+             const double t70 = -t23-t4+t16-t18+t6+t29-t44+
+                                vertices[2](0)*vertices[3](1);
+             jacobians_grad[point][0][0][0] = t35-t43*t46;
+             jacobians_grad[point][0][0][1] = -t50+t54*t46;
+             jacobians_grad[point][0][1][0] = t62*t46;
+             jacobians_grad[point][0][1][1] = -t67*t46;
+             jacobians_grad[point][1][0][0] = -t43*t70;
+             jacobians_grad[point][1][0][1] = t54*t70;
+             jacobians_grad[point][1][1][0] = -t35+t62*t70;
+             jacobians_grad[point][1][1][1] = t50-t67*t70;
+           };
+         break;
+         
+       };
+       
+       default:
+                                              // not implemented at present
+             Assert (false, ExcNotImplemented());
+      };             
 };
 
 #endif
index f9694d653f297f39f12a2b78483a0f0ba56ee0aa..eb0dd4c6d5a6ebe07d53d22ac77d3e618e5dbb78 100644 (file)
@@ -280,7 +280,9 @@ template <int dim>
 void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
                                           const vector<Point<dim> >            &unit_points,
                                           vector<Tensor<2,dim> >               &jacobians,
-                                          const bool           compute_jacobians,
+                                          const bool              compute_jacobians,
+                                          vector<Tensor<3,dim> > &jacobians_grad,
+                                          const bool              compute_jacobians_grad,
                                           vector<Point<dim> > &support_points,
                                           const bool           compute_support_points,
                                           vector<Point<dim> > &q_points,
@@ -379,7 +381,8 @@ void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator
   is multiplied to the unit cell gradient *from the right*! be very careful
   with these things.
 
-  The following little program tests the correct behaviour:
+  The following little program tests the correct behaviour. The program can
+  also be found in the <tests> directory.
 
   -------------------------------------------
   #include <grid/tria.h>
@@ -488,7 +491,94 @@ void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator
                                               // not implemented at present
              Assert (false, ExcNotImplemented());
       };
+
   
+  if (compute_jacobians_grad)
+    switch (dim) 
+      {
+       case 1:
+       {
+                                          // derivative of the
+                                          // jacobian is always zero
+                                          // for a linear mapping in 1d
+         for (unsigned int point=0; point<n_points; ++point)
+           jacobians_grad[point][0][0][0] = 0;
+         break;
+       };
+
+       case 2:
+       {
+         for (unsigned int point=0; point<n_points; ++point)
+           {
+             const double xi = unit_points[point](0);
+             const double eta= unit_points[point](1);
+       
+             const double t2 = vertices[1](0)*eta;
+             const double t4 = vertices[3](0)*vertices[2](1);
+             const double t6 = vertices[0](0)*vertices[2](1);
+             const double t8 = vertices[0](0)*vertices[3](1);
+             const double t10 = vertices[3](0)*xi;
+             const double t13 = vertices[2](0)*xi;
+             const double t16 = vertices[3](0)*vertices[1](1);
+             const double t18 = vertices[0](0)*vertices[1](1);
+             const double t19 = vertices[3](0)*vertices[0](1);
+             const double t20 = -t2*vertices[3](1)-t4*eta-t6*xi+t8*xi-
+                                t10*vertices[0](1)+t10*vertices[1](1)+
+                                t13*vertices[0](1)-t13*vertices[1](1)+t16
+                                *eta+t18+t19;
+             const double t23 = vertices[1](0)*vertices[3](1);
+             const double t26 = vertices[2](0)*eta;
+             const double t29 = vertices[1](0)*vertices[0](1);
+             const double t30 = vertices[1](0)*vertices[2](1);
+             const double t32 = -t16-t18*eta+t6*eta-t23*xi+t2*vertices[0](1)-
+                                t26*vertices[0](1)+t26*vertices[3](1)-
+                                t8-t29+t23+t30
+                                *xi;
+             const double t33 = t20+t32;
+             const double t34 = 1/t33;
+             const double t35 = (vertices[0](1)-vertices[1](1)+
+                                 vertices[2](1)-vertices[3](1))*t34;
+             const double t41 = t33*t33;
+             const double t42 = 1/t41;
+             const double t43 = (-vertices[0](1)+vertices[0](1)*xi-
+                                 vertices[1](1)*xi+vertices[2](1)*xi+
+                                 vertices[3](1)-vertices[3](1)*xi)*t42;
+             const double t44 = vertices[2](0)*vertices[0](1);
+             const double t46 = -t6+t8-t19+t16+t44-
+                                vertices[2](0)*vertices[1](1)-
+                                t23+t30;
+             const double t50 = (vertices[0](0)-vertices[1](0)+
+                                 vertices[2](0)-vertices[3](0))*t34;
+             const double t54 = (-vertices[0](0)+vertices[0](0)*xi-
+                                 vertices[1](0)*xi+t13+
+                                 vertices[3](0)-t10)*t42;
+             const double t62 = (-vertices[0](1)+vertices[0](1)*eta+
+                                 vertices[1](1)-vertices[1](1)*eta+
+                                 vertices[2](1)*eta-
+                                 vertices[3](1)*eta)*t42;
+             const double t67 = (-vertices[0](0)+vertices[0](0)*eta+
+                                 vertices[1](0)-t2+t26-vertices[3](0)*eta)*t42;
+             const double t70 = -t23-t4+t16-t18+t6+t29-t44+
+                                vertices[2](0)*vertices[3](1);
+             jacobians_grad[point][0][0][0] = t35-t43*t46;
+             jacobians_grad[point][0][0][1] = -t50+t54*t46;
+             jacobians_grad[point][0][1][0] = t62*t46;
+             jacobians_grad[point][0][1][1] = -t67*t46;
+             jacobians_grad[point][1][0][0] = -t43*t70;
+             jacobians_grad[point][1][0][1] = t54*t70;
+             jacobians_grad[point][1][1][0] = -t35+t62*t70;
+             jacobians_grad[point][1][1][1] = t50-t67*t70;
+           };
+         break;
+         
+       };
+       
+       default:
+                                              // not implemented at present
+             Assert (false, ExcNotImplemented());
+      };
+       
+             
   
     
   if (compute_support_points)
index b194e4b83ed398593f1812f8bb8e5b7544f2eed5..69a8ea3578ae4a1ad58ee627925a3aaa26e2202d 100644 (file)
@@ -28,6 +28,7 @@ FEValuesBase<dim>::FEValuesBase (const unsigned int n_q_points,
                n_transform_functions (n_transform_functions),
                shape_values (n_values_arrays, dFMatrix(n_dofs, n_q_points)),
                shape_gradients (n_dofs, vector<Tensor<1,dim> >(n_q_points)),
+               shape_2nd_derivatives (n_dofs, vector<Tensor<2,dim> >(n_q_points)),
                weights (n_q_points, 0),
                JxW_values (n_q_points, 0),
                quadrature_points (n_q_points, Point<dim>()),
@@ -127,6 +128,48 @@ void FEValuesBase<dim>::get_function_grads (const dVector          &fe_function,
 
 
 
+template <int dim>
+const Tensor<2,dim> &
+FEValuesBase<dim>::shape_2nd_derivative (const unsigned int i,
+                                        const unsigned int j) const {
+  Assert (i<shape_2nd_derivatives.size(),
+         ExcInvalidIndex (i, shape_2nd_derivatives.size()));
+  Assert (j<shape_2nd_derivatives[i].size(),
+         ExcInvalidIndex (j, shape_2nd_derivatives[i].size()));
+  Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField());
+
+  return shape_2nd_derivatives[i][j];
+};
+
+
+
+template <int dim>
+void FEValuesBase<dim>::get_function_2nd_derivatives (const dVector &fe_function,
+                                                     vector<Tensor<2,dim> > &second_derivatives) const {
+  Assert (second_derivatives.size() == n_quadrature_points,
+         ExcWrongVectorSize(second_derivatives.size(), n_quadrature_points));
+
+                                  // get function values of dofs
+                                  // on this cell
+  dVector dof_values (total_dofs);
+  present_cell->get_dof_values (fe_function, dof_values);
+
+                                  // initialize with zero
+  fill_n (second_derivatives.begin(), n_quadrature_points, Tensor<2,dim>());
+
+                                  // add up contributions of trial
+                                  // functions
+  for (unsigned int point=0; point<n_quadrature_points; ++point)
+    for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
+      {
+       Tensor<2,dim> tmp(shape_2nd_derivatives[shape_func][point]);
+       tmp *= dof_values(shape_func);
+       second_derivatives[point] += tmp;
+      };
+};
+
+
+
 template <int dim>
 const Point<dim> &
 FEValuesBase<dim>::quadrature_point (const unsigned int i) const {
@@ -218,10 +261,11 @@ void FEValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell,
   present_cell = cell;
                                   // fill jacobi matrices and real
                                   // quadrature points
-  if ((update_flags & update_jacobians) ||
-      (update_flags & update_JxW_values)||
-      (update_flags & update_q_points)  ||
-      (update_flags & update_gradients) ||
+  if ((update_flags & update_jacobians)          ||
+      (update_flags & update_JxW_values)         ||
+      (update_flags & update_q_points)           ||
+      (update_flags & update_gradients)          ||
+      (update_flags & update_second_derivatives) ||
       (update_flags & update_support_points))
     fe->fill_fe_values (cell,
                        unit_quadrature_points,
@@ -230,8 +274,8 @@ void FEValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell,
                                        update_JxW_values |
                                        update_gradients  |
                                        update_second_derivatives),
-//                     jacobi_matrices_grad,
-//                     update_flags & update_second_derivatives,
+                       jacobi_matrices_grad,
+                       update_flags & update_second_derivatives,
                        support_points,
                        update_flags & update_support_points,
                        quadrature_points,
@@ -407,31 +451,35 @@ void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &c
   selected_dataset = face_no;
                                   // fill jacobi matrices and real
                                   // quadrature points
-  if ((update_flags & update_jacobians) ||
-      (update_flags & update_JxW_values)||
-      (update_flags & update_q_points)  ||
-      (update_flags & update_gradients) ||
-      (update_flags & update_support_points) ||
+  if ((update_flags & update_jacobians)          ||
+      (update_flags & update_JxW_values)         ||
+      (update_flags & update_q_points)           ||
+      (update_flags & update_gradients)          ||
+      (update_flags & update_second_derivatives) ||
+      (update_flags & update_support_points)     ||
       (update_flags & update_JxW_values))
     fe->fill_fe_face_values (cell,
-                           face_no,
-                           unit_face_quadrature_points,
-                           unit_quadrature_points[face_no],
-                           jacobi_matrices,
-                           update_flags & (update_jacobians |
-                                           update_gradients |
-                                           update_JxW_values),
-                           support_points,
-                           update_flags & update_support_points,
-                           quadrature_points,
-                           update_flags & update_q_points,
-                           face_jacobi_determinants,
-                           update_flags & update_JxW_values,
-                           normal_vectors,
-                           update_flags & update_normal_vectors,
-                           shape_values_transform[face_no],
-                           unit_shape_gradients_transform[face_no],
-                           boundary);
+                            face_no,
+                            unit_face_quadrature_points,
+                            unit_quadrature_points[face_no],
+                            jacobi_matrices,
+                            update_flags & (update_jacobians |
+                                            update_gradients |
+                                            update_JxW_values |
+                                            update_second_derivatives),
+                            jacobi_matrices_grad,
+                            update_flags & update_second_derivatives,
+                            support_points,
+                            update_flags & update_support_points,
+                            quadrature_points,
+                            update_flags & update_q_points,
+                            face_jacobi_determinants,
+                            update_flags & update_JxW_values,
+                            normal_vectors,
+                            update_flags & update_normal_vectors,
+                            shape_values_transform[face_no],
+                            unit_shape_gradients_transform[face_no],
+                            boundary);
 
                                   // compute gradients on real element if
                                   // requested
@@ -579,29 +627,33 @@ void FESubfaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator
   selected_dataset = face_no*(1<<(dim-1)) + subface_no;
                                   // fill jacobi matrices and real
                                   // quadrature points
-  if ((update_flags & update_jacobians) ||
-      (update_flags & update_JxW_values)||
-      (update_flags & update_q_points)  ||
-      (update_flags & update_gradients) ||
+  if ((update_flags & update_jacobians)          ||
+      (update_flags & update_JxW_values)         ||
+      (update_flags & update_q_points)           ||
+      (update_flags & update_gradients)          ||
+      (update_flags & update_second_derivatives) ||
       (update_flags & update_JxW_values))
     fe->fill_fe_subface_values (cell,
-                              face_no,
-                              subface_no,
-                              unit_face_quadrature_points,
-                              unit_quadrature_points[selected_dataset],
-                              jacobi_matrices,
-                              update_flags & (update_jacobians |
-                                              update_gradients |
-                                              update_JxW_values),
-                              quadrature_points,
-                              update_flags & update_q_points,
-                              face_jacobi_determinants,
-                              update_flags & update_JxW_values,
-                              normal_vectors,
-                              update_flags & update_normal_vectors,
-                              shape_values_transform[selected_dataset],
-                              unit_shape_gradients_transform[selected_dataset],
-                              boundary);
+                               face_no,
+                               subface_no,
+                               unit_face_quadrature_points,
+                               unit_quadrature_points[selected_dataset],
+                               jacobi_matrices,
+                               update_flags & (update_jacobians |
+                                               update_gradients |
+                                               update_JxW_values|
+                                               update_second_derivatives),
+                               jacobi_matrices_grad,
+                               update_flags & update_second_derivatives,
+                               quadrature_points,
+                               update_flags & update_q_points,
+                               face_jacobi_determinants,
+                               update_flags & update_JxW_values,
+                               normal_vectors,
+                               update_flags & update_normal_vectors,
+                               shape_values_transform[selected_dataset],
+                               unit_shape_gradients_transform[selected_dataset],
+                               boundary);
 
                                   // compute gradients on real element if
                                   // requested
index 57425e87fc311fa01cc2bbba78c11205afb70f20..e7034ba43079e9aa8f45c77db1c82d6418be85e7 100644 (file)
@@ -22,6 +22,7 @@ perl -pi -e 's/\[(\d)\]\[(\d)\] =/($1,$2) =/g;' *2d_inverse_jacobian
 perl -pi -e 's/^\s*t/const double t/g;' *2d_inverse_jacobian_grad
 perl -pi -e 's/x\[(\d)\]/vertices[$1](0)/g;' *2d_inverse_jacobian_grad
 perl -pi -e 's/y\[(\d)\]/vertices[$1](1)/g;' *2d_inverse_jacobian_grad
+perl -pi -e 's/inverseJacobian/jacobians_grad[point]/g;' *2d_inverse_jacobian_grad
 
 perl -pi -e 's/([^;])\n/$1/g;' *2d_shape_grad
 perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[0\] = (.*);/case $1: return Point<2>($2,/g;' *2d_shape_grad

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.