]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
More work with regard to the transition of gradients from data type Point<dim> to...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 7 Nov 1998 02:15:57 +0000 (02:15 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 7 Nov 1998 02:15:57 +0000 (02:15 +0000)
git-svn-id: https://svn.dealii.org/trunk@655 0785d39b-7218-0410-832d-ea1e28bc413d

13 files changed:
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_values.h
deal.II/deal.II/include/grid/tria.h
deal.II/deal.II/source/dofs/dof_handler.cc
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.linear.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/error_estimator.cc
deal.II/deal.II/source/numerics/vectors.cc

index 89f92709280b753c2396c0f973a974250f8956e1..69a3bc511d5f3945ad93c18c35f5aa6c015715ea 100644 (file)
@@ -578,8 +578,8 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      * #i#th shape function for the mapping
                                      * from unit to real cell.
                                      */
-    virtual Point<dim> shape_grad_transform (const unsigned int i,
-                                            const Point<dim> &p) const = 0;    
+    virtual Tensor<1,dim> shape_grad_transform (const unsigned int i,
+                                               const Point<dim> &p) const = 0;    
     
                                     /**
                                      * Compute the Jacobian matrix and the
@@ -646,7 +646,7 @@ class FiniteElement : public FiniteElementBase<dim> {
                                 vector<Point<dim> > &q_points,
                                 const bool           compute_q_points,
                                 const dFMatrix      &shape_values_transform,
-                                const vector<vector<Point<dim> > > &shape_grads_transform,
+                                const vector<vector<Tensor<1,dim> > > &shape_grads_transform,
                                 const Boundary<dim> &boundary) const;
 
                                     /**
@@ -768,7 +768,7 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      vector<Point<dim> > &normal_vectors,
                                      const bool           compute_normal_vectors,
                                      const dFMatrix      &shape_values_transform,
-                                     const vector<vector<Point<dim> > > &shape_grads_transform,
+                                     const vector<vector<Tensor<1,dim> > > &shape_grads_transform,
                                      const Boundary<dim> &boundary) const;
 
                                     /**
@@ -813,7 +813,7 @@ class FiniteElement : public FiniteElementBase<dim> {
                                         vector<Point<dim> > &normal_vectors,
                                         const bool           compute_normal_vectors,
                                         const dFMatrix      &shape_values_transform,
-                                        const vector<vector<Point<dim> > > &shape_grads_transform,
+                                        const vector<vector<Tensor<1,dim> > > &shape_grads_transform,
                                         const Boundary<dim> &boundary) const;
 
                                     /**
index d27931ed9625c56f166eb9cc7635c5687ba6fa68..68d41b8ea3a164dcdfba415d2d11cf8334cb6808 100644 (file)
@@ -210,16 +210,16 @@ class FECrissCross : public FiniteElement<dim> {
                                      * information on this function.
                                      */
     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;
 
                                     /**
                                      * Refer to the base class for detailed
                                      * information on this function.
                                      */
     virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
-                                        const Boundary<dim> &boundary,
-                                        vector<Point<dim> > &support_points) const;
+                                         const Boundary<dim> &boundary,
+                                         vector<Point<dim> > &support_points) const;
 
                                     /**
                                      * Refer to the base class for detailed
@@ -248,8 +248,8 @@ class FECrissCross : public FiniteElement<dim> {
                                      * computation of the transformation from
                                      * unit cell to real space cell.
                                      */
-    virtual Point<dim> shape_grad_transform (const unsigned int i,
-                                            const Point<dim> &p) const;
+    virtual Tensor<1,dim> shape_grad_transform (const unsigned int i,
+                                               const Point<dim> &p) const;
 
                                     /**
                                      * Refer to the base class for detailed
@@ -339,7 +339,7 @@ class FECrissCross : public FiniteElement<dim> {
                                 vector<Point<dim> > &q_points,
                                 const bool           compute_q_points,
                                 const dFMatrix      &shape_values_transform,
-                                const vector<vector<Point<dim> > > &shape_grad_transform,
+                                const vector<vector<Tensor<1,dim> > > &shape_grad_transform,
                                 const Boundary<dim> &boundary) const;
 
     DeclException0 (ExcNotUseful);
index 735023da01c908c2d46cdd69f019826aef68e68d..9be39b199fb1775397a083cd1ed134acd2dc9c9b 100644 (file)
@@ -51,8 +51,8 @@ class FELinearMapping : public FiniteElement<dim> {
                                      * computation of the transformation from
                                      * unit cell to real space cell.
                                      */
-    virtual Point<dim> shape_grad_transform (const unsigned int i,
-                                            const Point<dim> &p) const;
+    virtual Tensor<1,dim> shape_grad_transform (const unsigned int i,
+                                               const Point<dim> &p) const;
 
                                     /**
                                      * Refer to the base class for detailed
@@ -142,7 +142,7 @@ class FELinearMapping : public FiniteElement<dim> {
                                 vector<Point<dim> > &q_points,
                                 const bool           compute_q_points,
                                 const dFMatrix      &shape_values_transform,
-                                const vector<vector<Point<dim> > > &shape_grad_transform,
+                                const vector<vector<Tensor<1,dim> > > &shape_grad_transform,
                                 const Boundary<dim> &boundary) const;
 };
 
index 12fa99453f11eedb7a673eb044c1112e80f79c57..303c126a246261c32f9d842263ffbda63e2d66a9 100644 (file)
@@ -297,7 +297,7 @@ class FEValuesBase {
                                      * function at the #j# quadrature point.
                                      * If you want to get the derivative in
                                      * one of the coordinate directions, use
-                                     * the appropriate function of the #Point#
+                                     * the appropriate function of the #Tensor#
                                      * class to extract one component. Since
                                      * only a reference to the gradient's value
                                      * is returned, there should be no major
@@ -327,8 +327,8 @@ class FEValuesBase {
                                      * #gradients# object already has the
                                      * right size.
                                      */
-    void get_function_grads (const dVector       &fe_function,
-                            vector<Point<dim> > &gradients) const;
+    void get_function_grads (const dVector          &fe_function,
+                            vector<Tensor<1,dim> > &gradients) const;
 
                                     /**
                                      * Return the position of the #i#th
@@ -704,7 +704,7 @@ class FEValues : public FEValuesBase<dim> {
                                      * Analogous to the #shape_values_transform#
                                      * array of the base class.
                                      */
-    vector<vector<Point<dim> > >   unit_shape_gradients_transform;
+    vector<vector<Tensor<1,dim> > > unit_shape_gradients_transform;
     
                                     /**
                                      * Array of quadrature points in the unit
@@ -861,7 +861,7 @@ class FEFaceValuesBase : public FEValuesBase<dim> {
                                      * Analogous to the #shape_values_transform#
                                      * array of the base class.
                                      */
-    vector<vector<vector<Point<dim> > > > unit_shape_gradients_transform;
+    vector<vector<vector<Tensor<1,dim> > > > unit_shape_gradients_transform;
 
                                     /**
                                      * Array of quadrature points on the
index f67067033334ac13f3f8b149ff28a51d1ea474b3..1cb011c71fb23da072e9584e653ef9b869c2acae 100644 (file)
@@ -8,7 +8,7 @@
 #include <vector>
 #include <base/point.h>
 #include <grid/geometry_info.h>
-
+#include <base/subscriptor.h>
 
 
 //forward declaration needed
@@ -1144,7 +1144,7 @@ class TriaDimensionInfo<2> {
  *   @author Wolfgang Bangerth, 1998
  */
 template <int dim>
-class Triangulation : public TriaDimensionInfo<dim> {
+class Triangulation : public TriaDimensionInfo<dim>, public Subscriptor {
   public:
                                     // insert these definitions for gcc2.8,
                                     // since it can't inherit typedefs (I
index fb090012be9d7329ddbb323db8d010c06c8ac1ba..5a79b3586887caf8b4198bac1110235fe3568a1a 100644 (file)
 template <int dim>
 DoFHandler<dim>::DoFHandler (Triangulation<dim> *tria) :
                tria(tria),
-               used_dofs (0) {};
+               used_dofs (0)
+{
+  tria->subscribe ();
+};
 
 
 template <int dim>
-DoFHandler<dim>::~DoFHandler ()
-{}
+DoFHandler<dim>::~DoFHandler () {
+  tria->unsubscribe ();
+}
 
 
 
index 0833482a4a9b4ca49952e1529b939dd810bdd6ff..57151da0e95cdac968f6ceeb19ae4ce5eaf8644a 100644 (file)
@@ -185,7 +185,7 @@ void FiniteElement<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
                                       vector<Point<1> > &q_points,
                                       const bool         compute_q_points,
                                       const dFMatrix      &,
-                                      const vector<vector<Point<1> > > &,
+                                      const vector<vector<Tensor<1,1> > > &,
                                       const Boundary<1> &boundary) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
@@ -235,7 +235,7 @@ void FiniteElement<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &
                                            vector<Point<1> >       &,
                                            const bool,
                                            const dFMatrix          &,
-                                           const vector<vector<Point<1> > > &,
+                                           const vector<vector<Tensor<1,1> > > &,
                                            const Boundary<1>       &) const {
   Assert (false, ExcNotImplemented());
 };
@@ -256,7 +256,7 @@ void FiniteElement<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterato
                                               vector<Point<1> >       &,
                                               const bool,
                                               const dFMatrix          &,
-                                              const vector<vector<Point<1> > > &,
+                                              const vector<vector<Tensor<1,1> > > &,
                                               const Boundary<1>       &) const {
   Assert (false, ExcNotImplemented());
 };
@@ -323,7 +323,7 @@ void FiniteElement<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &,
                                         vector<Point<dim> > &,
                                         const bool,
                                         const dFMatrix      &,
-                                        const vector<vector<Point<dim> > > &,
+                                        const vector<vector<Tensor<1,dim> > > &,
                                         const Boundary<dim> &) const {
   Assert (false, ExcPureFunctionCalled());
 };
@@ -346,7 +346,7 @@ void FiniteElement<dim>::fill_fe_face_values (const DoFHandler<dim>::cell_iterat
                                              vector<Point<dim> > &normal_vectors,
                                              const bool           compute_normal_vectors,
                                              const dFMatrix      &shape_values_transform,
-                                             const vector<vector<Point<dim> > > &shape_gradients_transform,
+                                             const vector<vector<Tensor<1,dim> > > &shape_gradients_transform,
                                              const Boundary<dim> &boundary) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
@@ -395,7 +395,7 @@ void FiniteElement<dim>::fill_fe_subface_values (const DoFHandler<dim>::cell_ite
                                                 vector<Point<dim> > &normal_vectors,
                                                 const bool           compute_normal_vectors,
                                                 const dFMatrix      &shape_values_transform,
-                                                const vector<vector<Point<dim> > > &shape_gradients_transform,
+                                                const vector<vector<Tensor<1,dim> > > &shape_gradients_transform,
                                                 const Boundary<dim> &boundary) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
index 5cace7e1d7a9cd3b06f1552d7d3b1050b4671baf..0af836516f4eae63c514508e96ff48ef52eda336 100644 (file)
@@ -347,8 +347,8 @@ void FECrissCross<1>::get_unit_support_points (vector<Point<1> >&) const {
 
 template <>
 void FECrissCross<1>::get_support_points (const DoFHandler<1>::cell_iterator &,
-                                        const Boundary<1> &,
-                                        vector<Point<1> > &) const {
+                                         const Boundary<1> &,
+                                         vector<Point<1> > &) const {
   Assert (false, ExcNotUseful());
 };
 
@@ -356,8 +356,8 @@ void FECrissCross<1>::get_support_points (const DoFHandler<1>::cell_iterator &,
 
 template <>
 void FECrissCross<1>::get_face_support_points (const DoFHandler<1>::face_iterator &,
-                                             const Boundary<1> &,
-                                             vector<Point<1> > &) const {
+                                              const Boundary<1> &,
+                                              vector<Point<1> > &) const {
   Assert (false, ExcNotUseful());
 };
 
@@ -382,8 +382,8 @@ double  FECrissCross<1>::shape_value_transform (const unsigned int,
 
 
 template <>
-Point<1> FECrissCross<1>::shape_grad_transform (const unsigned int,
-                                               const Point<1> &) const {
+Tensor<1,1> FECrissCross<1>::shape_grad_transform (const unsigned int,
+                                                  const Point<1> &) const {
   Assert (false, ExcNotUseful());
   return Point<1>();
 };
@@ -442,7 +442,7 @@ void FECrissCross<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &,
                                      vector<Point<1> > &,
                                      const bool,
                                      const dFMatrix      &,
-                                     const vector<vector<Point<1> > > &,
+                                     const vector<vector<Tensor<1,1> > > &,
                                      const Boundary<1> &) const {
   Assert (false, ExcNotUseful());
 };
@@ -594,8 +594,8 @@ void FECrissCross<2>::get_unit_support_points (vector<Point<2> > &unit_points) c
 
 template <>
 void FECrissCross<2>::get_support_points (const DoFHandler<2>::cell_iterator &cell,
-                                        const Boundary<2> &,
-                                        vector<Point<2> > &support_points) const {
+                                         const Boundary<2> &,
+                                         vector<Point<2> > &support_points) const {
   const unsigned int dim = 2;
   
   Assert (support_points.size() == total_dofs,
@@ -632,8 +632,8 @@ void FECrissCross<2>::get_support_points (const DoFHandler<2>::cell_iterator &ce
 
 template <>
 void FECrissCross<2>::get_face_support_points (const DoFHandler<2>::face_iterator &face,
-                                             const Boundary<2> &,
-                                             vector<Point<2> > &support_points) const {
+                                              const Boundary<2> &,
+                                              vector<Point<2> > &support_points) const {
   const unsigned int dim = 2;
   
   Assert ((support_points.size() == dofs_per_face) &&
@@ -793,7 +793,7 @@ void FECrissCross<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator
 template <>
 inline
 double FECrissCross<2>::shape_value_transform (const unsigned int i,
-                                               const Point<2> &p) const {
+                                              const Point<2> &p) const {
                                   // use an isoparametric ansatz
   return shape_value(i,p);
 };
@@ -801,8 +801,8 @@ double FECrissCross<2>::shape_value_transform (const unsigned int i,
 
 
 template <>
-Point<2> FECrissCross<2>::shape_grad_transform (const unsigned int i,
-                                               const Point<2> &p) const {
+Tensor<1,2> FECrissCross<2>::shape_grad_transform (const unsigned int i,
+                                                  const Point<2> &p) const {
                                   // use an isoparametric ansatz
   return shape_grad(i,p);  
 };
@@ -934,7 +934,7 @@ void FECrissCross<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &ce
                                        vector<Point<dim> > &q_points,
                                        const bool           compute_q_points,
                                        const dFMatrix      &shape_values_transform,
-                                       const vector<vector<Point<dim> > > &/*shape_grad_transform*/,
+                                       const vector<vector<Tensor<1,dim> > > &/*shape_grad_transform*/,
                                        const Boundary<dim> &boundary) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
@@ -999,40 +999,57 @@ void FECrissCross<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &ce
     vertices[l] = cell->vertex(l);
 
   if (compute_jacobians)
-    for (unsigned int point=0; point<n_points; ++point)
+    switch (dim)
       {
-       const double xi = unit_points[point](0);
-       const double eta= unit_points[point](1);
+       case 1:
+             for (unsigned int point=0; point<n_points; ++point)
+               jacobians[point][0][0] = 1./(vertices[1](0)-vertices[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 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
+             Assert (false, ExcNotImplemented());
       };
 };
 
index 08d1580d4da71648461a23fbbef3b671ef98c7f5..04ef2edfdf2fad68d74da08608a0ced7b3b74930 100644 (file)
@@ -379,9 +379,10 @@ void FELinear<2>::get_unit_support_points (vector<Point<2> > &unit_points) const
 
 
 template <int dim>
-void FELinear<dim>::get_support_points (const typename DoFHandler<dim>::cell_iterator &cell,
-                                      const Boundary<dim>  &,
-                                      vector<Point<dim> >  &support_points) const {
+void
+FELinear<dim>::get_support_points (const typename DoFHandler<dim>::cell_iterator &cell,
+                                  const Boundary<dim>  &,
+                                  vector<Point<dim> >  &support_points) const {
   Assert (support_points.size() == total_dofs,
          ExcWrongFieldDimension (support_points.size(), total_dofs));
   
@@ -392,9 +393,10 @@ void FELinear<dim>::get_support_points (const typename DoFHandler<dim>::cell_ite
 
 
 template <int dim>
-void FELinear<dim>::get_face_support_points (const typename DoFHandler<dim>::face_iterator &face,
-                                           const Boundary<dim>  &,
-                                           vector<Point<dim> >  &support_points) const {
+void
+FELinear<dim>::get_face_support_points (const typename DoFHandler<dim>::face_iterator &face,
+                                       const Boundary<dim>  &,
+                                       vector<Point<dim> >  &support_points) const {
   Assert ((support_points.size() == dofs_per_face) &&
          (support_points.size() == GeometryInfo<dim>::vertices_per_face),
          ExcWrongFieldDimension (support_points.size(),
index e2728ad7635cfc73b8d452a7fffffb19a15ca7c9..f9694d653f297f39f12a2b78483a0f0ba56ee0aa 100644 (file)
@@ -36,7 +36,7 @@ FELinearMapping<1>::shape_value_transform (const unsigned int i,
 
 template <>
 inline
-Point<1>
+Tensor<1,1>
 FELinearMapping<1>::shape_grad_transform(const unsigned int i,
                                         const Point<1>&) const
 {
@@ -102,7 +102,7 @@ void FELinearMapping<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cel
                                         vector<Point<1> > &q_points,
                                         const bool         compute_q_points,
                                         const dFMatrix      &shape_values_transform,
-                                        const vector<vector<Point<1> > > &shape_gradients_transform,
+                                        const vector<vector<Tensor<1,1> > > &shape_gradients_transform,
                                         const Boundary<1> &boundary) const {
                                   // simply pass down
   FiniteElement<1>::fill_fe_values (cell, unit_points,
@@ -142,7 +142,7 @@ FELinearMapping<2>::shape_value_transform (const unsigned int i,
 
 template <>
 inline
-Point<2>
+Tensor<1,2>
 FELinearMapping<2>::shape_grad_transform (const unsigned int i,
                                          const Point<2>& p) const
 {
@@ -286,7 +286,7 @@ void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator
                                           vector<Point<dim> > &q_points,
                                           const bool           compute_q_points,
                                           const dFMatrix      &shape_values_transform,
-                                          const vector<vector<Point<dim> > > &/*shape_grad_transform*/,
+                                          const vector<vector<Tensor<1,dim> > > &/*shape_grad_transform*/,
                                           const Boundary<dim> &boundary) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
@@ -379,11 +379,19 @@ 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 (you have to find
-  out the right include files, I tested it within a whole project with far
-  more include files than necessary):
+  The following little program tests the correct behaviour:
 
   -------------------------------------------
+  #include <grid/tria.h>
+  #include <grid/tria_boundary.h>
+  #include <grid/dof.h>
+  #include <fe/fe_values.h>
+  #include <fe/fe_lib.lagrange.h>
+  #include <base/quadrature_lib.h>
+  #include <grid/tria_iterator.h>
+  #include <grid/dof_accessor.h>
+  #include <lac/dvector.h>
+
   int main () {
     Triangulation<2> tria;
     tria.create_hypercube (0,1);
@@ -428,53 +436,58 @@ void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator
 */
 
   if (compute_jacobians)
-    for (unsigned int point=0; point<n_points; ++point)
-      switch (dim)
+    switch (dim)
+      {
+       case 1:
+             for (unsigned int point=0; point<n_points; ++point)
+               jacobians[point][0][0] = 1./(vertices[1](0)-vertices[0](0));
+             break;
+               
+       case 2:
        {
-         case 2:
-         {
-           
-           const double xi = unit_points[point](0);
-           const double eta= unit_points[point](1);
+         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 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());
+             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
+             Assert (false, ExcNotImplemented());
+      };
   
   
     
index b2ec9d3e8ce4fdb8fc00aba2095be1a8104f74d7..b194e4b83ed398593f1812f8bb8e5b7544f2eed5 100644 (file)
@@ -101,8 +101,8 @@ FEValuesBase<dim>::shape_grad (const unsigned int i,
 
 
 template <int dim>
-void FEValuesBase<dim>::get_function_grads (const dVector       &fe_function,
-                                           vector<Point<dim> > &gradients) const {
+void FEValuesBase<dim>::get_function_grads (const dVector          &fe_function,
+                                           vector<Tensor<1,dim> > &gradients) const {
   Assert (gradients.size() == n_quadrature_points,
          ExcWrongVectorSize(gradients.size(), n_quadrature_points));
 
@@ -112,7 +112,7 @@ void FEValuesBase<dim>::get_function_grads (const dVector       &fe_function,
   present_cell->get_dof_values (fe_function, dof_values);
 
                                   // initialize with zero
-  fill_n (gradients.begin(), n_quadrature_points, Point<dim>());
+  fill_n (gradients.begin(), n_quadrature_points, Tensor<1,dim>());
 
                                   // add up contributions of trial
                                   // functions
@@ -128,7 +128,8 @@ void FEValuesBase<dim>::get_function_grads (const dVector       &fe_function,
 
 
 template <int dim>
-const Point<dim> & FEValuesBase<dim>::quadrature_point (const unsigned int i) const {
+const Point<dim> &
+FEValuesBase<dim>::quadrature_point (const unsigned int i) const {
   Assert (i<quadrature_points.size(), ExcInvalidIndex(i,quadrature_points.size()));
   Assert (update_flags & update_q_points, ExcAccessToUninitializedField());
   
@@ -138,7 +139,8 @@ const Point<dim> & FEValuesBase<dim>::quadrature_point (const unsigned int i) co
 
 
 template <int dim>
-const Point<dim> & FEValuesBase<dim>::support_point (const unsigned int i) const {
+const Point<dim> &
+FEValuesBase<dim>::support_point (const unsigned int i) const {
   Assert (i<support_points.size(), ExcInvalidIndex(i, support_points.size()));
   Assert (update_flags & update_support_points, ExcAccessToUninitializedField());
   
@@ -176,7 +178,7 @@ FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
                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)),
+                                              vector<Tensor<1,dim> >(quadrature.n_quadrature_points)),
                unit_quadrature_points(quadrature.get_quad_points())
 {
   Assert ((update_flags & update_normal_vectors) == false,
@@ -320,8 +322,8 @@ FEFaceValuesBase<dim>::FEFaceValuesBase (const unsigned int n_q_points,
                                           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))),
+                                               vector<vector<Tensor<1,dim> > >(n_transform_functions,
+                                                                               vector<Tensor<1,dim> >(n_q_points))),
                unit_face_quadrature_points (n_q_points, Point<dim-1>()),
                unit_quadrature_points (n_faces_or_subfaces,
                                        vector<Point<dim> >(n_q_points, Point<dim>())),
@@ -332,7 +334,8 @@ FEFaceValuesBase<dim>::FEFaceValuesBase (const unsigned int n_q_points,
 
 
 template <int dim>
-const Point<dim> & FEFaceValuesBase<dim>::normal_vector (const unsigned int i) const {
+const Point<dim> &
+FEFaceValuesBase<dim>::normal_vector (const unsigned int i) const {
   Assert (i<normal_vectors.size(), ExcInvalidIndex(i, normal_vectors.size()));
   Assert (update_flags & update_normal_vectors,
          ExcAccessToUninitializedField());
@@ -437,7 +440,7 @@ void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &c
       {
        fill_n (shape_gradients[i].begin(),
                n_quadrature_points,
-               Point<dim>());
+               Tensor<1,dim>());
        for (unsigned int j=0; j<n_quadrature_points; ++j) 
          for (unsigned int s=0; s<dim; ++s)
                                             // (grad psi)_s =
@@ -607,7 +610,7 @@ void FESubfaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator
       {
        fill_n (shape_gradients[i].begin(),
                n_quadrature_points,
-               Point<dim>());
+               Tensor<1,dim>());
        for (unsigned int j=0; j<n_quadrature_points; ++j) 
          for (unsigned int s=0; s<dim; ++s)
                                             // (grad psi)_s =
index f0102a9f6026e0b6ab1a2d72e32d8fd6e45e9716..b25f34be0e526f2333c9fad835480c2b434d675f 100644 (file)
@@ -246,7 +246,7 @@ integrate_over_regular_face (const active_cell_iterator &cell,
                                   //
                                   // let psi be a short name for
                                   // [a grad u_h]
-  vector<Point<dim> > psi(n_q_points);
+  vector<Tensor<1,dim> > psi(n_q_points);
   fe_face_values_cell.get_function_grads (solution, psi);
   
   
@@ -280,14 +280,14 @@ integrate_over_regular_face (const active_cell_iterator &cell,
                                       // get a list of the gradients of
                                       // the finite element solution
                                       // restricted to the neighbor cell
-      vector<Point<dim> > neighbor_psi (n_q_points);
+      vector<Tensor<1,dim> > neighbor_psi (n_q_points);
       fe_face_values_neighbor.get_function_grads (solution, neighbor_psi);
       
                                       // compute the jump in the gradients
       transform (psi.begin(), psi.end(),
                 neighbor_psi.begin(),
                 psi.begin(),
-                minus<Point<dim> >());
+                minus<Tensor<1,dim> >());
     };
   
   
@@ -397,7 +397,7 @@ integrate_over_irregular_face (const active_cell_iterator &cell,
                                   //
                                   // let psi be a short name for
                                   // [a grad u_h]
-  vector<Point<dim> > psi(n_q_points);
+  vector<Tensor<1,dim> > psi(n_q_points);
 
                                   // store which number #cell# has in the
                                   // list of neighbors of #neighbor#
@@ -434,7 +434,7 @@ integrate_over_irregular_face (const active_cell_iterator &cell,
                                       // restrict the finite element on the
                                       // neighbor cell to the common #subface#.
                                       // store the gradient in #neighbor_psi#
-      vector<Point<dim> > neighbor_psi (n_q_points);
+      vector<Tensor<1,dim> > neighbor_psi (n_q_points);
       fe_face_values.reinit (neighbor_child, neighbor_neighbor, boundary);
       fe_face_values.get_function_grads (solution, neighbor_psi);
       
@@ -442,7 +442,7 @@ integrate_over_irregular_face (const active_cell_iterator &cell,
       transform (psi.begin(), psi.end(),
                 neighbor_psi.begin(),
                 psi.begin(),
-                minus<Point<dim> >());
+                minus<Tensor<1,dim> >());
 
                                       // next we have to multiply this with
                                       // the normal vector. Since we have
index 9f21f6718dcded7a620da36eb4f83a2791cf84f7..3e7ce4b779186861aaf954d92b6cc9f0992a6f4e 100644 (file)
@@ -32,8 +32,8 @@ inline double sqr (const double x) {
 
 
 template <int dim>
-inline double sqr_point (const Point<dim> &p) {
-  return p.square();
+inline double sqr_point (const Tensor<1,dim> &p) {
+  return p * p;
 };
 
 
@@ -507,7 +507,7 @@ void VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
                                             // same procedure as above, but now
                                             // psi is a vector of gradients
            const unsigned int n_q_points = q.n_quadrature_points;
-           vector<Point<dim> >   psi (n_q_points);
+           vector<Tensor<1,dim> >   psi (n_q_points);
 
                                             // in praxi: first compute
                                             // exact fe_function vector
@@ -518,13 +518,13 @@ void VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
                                             // fe_function
            if (true) 
              {
-               vector<Point<dim> > function_grads (n_q_points, Point<dim>());
+               vector<Tensor<1,dim> > function_grads (n_q_points, Tensor<1,dim>());
                fe_values.get_function_grads (fe_function, function_grads);
 
                transform (psi.begin(), psi.end(),
                           function_grads.begin(),
                           psi.begin(),
-                          minus<Point<dim> >());
+                          minus<Tensor<1,dim> >());
              };
                                             // take square of integrand
            vector<double> psi_square (psi.size(), 0.0);

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.