]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new interpolation functions; documentation
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Sat, 24 Sep 2005 11:26:56 +0000 (11:26 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Sat, 24 Sep 2005 11:26:56 +0000 (11:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@11526 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_raviart_thomas.h
deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc

index 9c0115d5eddf6ed734a57677777e06f0d9bf9e87..19040d23ed673daec0a84c7582c9af2bcd1d456e 100644 (file)
@@ -90,6 +90,10 @@ template <int dim> class FECollection;
  *
  * <h3>Notes on the implementation of derived classes</h3>
  *
+ * The following sections list the information to be provided by
+ * derived classes, depending on the space dimension. They are
+ * followed by a list of functions helping to generate these values.
+ *
  * <h4>Finite elements in one dimension</h4>
  *
  * Finite elements in one dimension need only set the #restriction
@@ -213,6 +217,43 @@ template <int dim> class FECollection;
  * constraints that are entered more than once (as is necessary for the case
  * above), it insists that the weights are exactly the same.
  *
+ * <h4>Helper functions</h4>
+ *
+ * Construction of a finite element and computation of the matrices
+ * described above may be a tedious task, in particular if it has to
+ * be performed for several space dimensions. Therefore, some
+ * functions in FETools have been provided to help with these tasks.
+ *
+ * First, aready the basis of the shape function space may be
+ * difficult to implement for arbitrary order and dimension. On the
+ * other hand, if the @ref GlossNodes "node values" are given, then
+ * the duality relation between node functionals and basis functions
+ * defines the basis. As a result, the shape function space may be
+ * defined with arbitrary "raw" basis functions, such that the actual
+ * finite element basis is computed from linear combinations of
+ * them. The coefficients of these combinations are determined by the
+ * duality of node values.
+ *
+ * Using this matrix allows the construction of the basis of shape
+ * functions in two steps.
+ * <ol>
+ *
+ * <li>Define the space of shape functions using an arbitrary basis
+ * <i>w<sub>j</sub></i> and compute the matrix <i>M</i> of node
+ * functionals <i>N<sub>i</sub></i> applied to these basis functions,
+ * such that its entries are <i>m<sub>ij</sub> =
+ * N<sub>i</sub>(w<sub>j</sub>)</i>.
+ *
+ * <li>Compute the basis <i>v<sub>j</sub></i> of the finite element
+ * shape function space by applying <i>M<sup>-1</sup></i> to the basis
+ * <i>w<sub>j</sub></i>.
+ * </ol>
+ *
+ * The function computing the matrix <i>M</i> for you is
+ * FETools::compute_node_matrix(). It relies on the existence of
+ * #generalized_support_points and implementation of interpolate()
+ * with VectorSlice argument.
+ *
  * @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1998, 2000, 2001, 2005
  */
 template <int dim>
index d01d299ac0d24c92722610037d83be84e7698749..bfb805ff94f82c508b6709dc36d353b322590660 100644 (file)
@@ -334,6 +334,15 @@ class FE_RaviartThomas : public FiniteElement<dim>
                            typename Mapping<dim>::InternalDataBase      &fe_internal,
                            FEValuesData<dim>& data) const;
 
+//     virtual void interpolate(std::vector<double>&                local_dofs,
+//                          const std::vector<double>& values) const;
+//     virtual void interpolate(std::vector<double>&                local_dofs,
+//                          const std::vector<Vector<double> >& values,
+//                          unsigned int offset = 0) const;
+    
+//     virtual void interpolate(
+//       std::vector<double>& local_dofs,
+//       const VectorSlice<const std::vector<std::vector<double> > >& values) const;
   private:
                                     /**
                                      * The order of the
@@ -630,6 +639,14 @@ class FE_RaviartThomasNodal
                                      */
     virtual bool has_support_on_face (const unsigned int shape_index,
                                      const unsigned int face_index) const;    
+    virtual void interpolate(std::vector<double>&                local_dofs,
+                            const std::vector<double>& values) const;
+    virtual void interpolate(std::vector<double>&                local_dofs,
+                            const std::vector<Vector<double> >& values,
+                            unsigned int offset = 0) const;    
+    virtual void interpolate(
+      std::vector<double>& local_dofs,
+      const VectorSlice<const std::vector<std::vector<double> > >& values) const;
   private:
                                     /**
                                      * Only for internal use. Its
index e22c31beb40f614f377e2ce7bd259b4e3bbf1410..08fd82fd0d8ea5aae4167f8319086e4e64ff976f 100644 (file)
@@ -52,7 +52,7 @@ FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal (const unsigned int deg)
                                   // basis functions
   initialize_unit_support_points (deg);
   initialize_node_matrix();
-
+  
   for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
     this->prolongation[i].reinit (this->dofs_per_cell,
                                  this->dofs_per_cell);
@@ -96,6 +96,111 @@ FE_RaviartThomasNodal<dim>::has_support_on_face (unsigned int, unsigned int) con
 }
 
 
+template <int dim>
+void
+FE_RaviartThomasNodal<dim>::interpolate(
+  std::vector<double>&,
+  const std::vector<double>&) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+template <int dim>
+void
+FE_RaviartThomasNodal<dim>::interpolate(
+  std::vector<double>&    local_dofs,
+  const std::vector<Vector<double> >& values,
+  unsigned int offset) const
+{
+  Assert (values.size() == generalized_support_points.size(),
+         ExcDimensionMismatch(values.size(), generalized_support_points.size()));
+  Assert (local_dofs.size() == this->dofs_per_cell,
+         ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+  Assert (values[0].size() >= offset+this->n_components(),
+         ExcDimensionMismatch(values[0].size(),offset+this->n_components()));
+
+                                  // First do interpolation on
+                                  // faces. There, the component
+                                  // evaluated depends on the face
+                                  // direction and orientation.
+  unsigned int fbase = 0;
+  unsigned int f=0;
+  for (;f<GeometryInfo<dim>::faces_per_cell;
+       ++f, fbase+=this->dofs_per_face)
+    {
+      for (unsigned int i=0;i<this->dofs_per_face;++i)
+       {
+         local_dofs[fbase+i] = values[fbase+i](offset+GeometryInfo<dim>::unit_normal_direction[f]);
+       }
+    }
+
+                                  // The remaining points form dim
+                                  // chunks, one for each component.
+  const unsigned int istep = (this->dofs_per_cell - fbase) / dim;
+  Assert ((this->dofs_per_cell - fbase) % dim == 0, ExcInternalError());
+
+  f = 0;
+  while (fbase < this->dofs_per_cell)
+    {
+      for (unsigned int i=0;i<istep;++i)
+       {
+         local_dofs[fbase+i] = values[fbase+i](offset+f);
+       }
+      fbase+=istep;
+      ++f;
+    }
+  Assert (fbase == this->dofs_per_cell, ExcInternalError());
+}
+
+
+
+template <int dim>
+void
+FE_RaviartThomasNodal<dim>::interpolate(
+  std::vector<double>& local_dofs,
+  const VectorSlice<const std::vector<std::vector<double> > >& values) const
+{
+  Assert (values.size() == this->n_components(),
+         ExcDimensionMismatch(values.size(), this->n_components()));
+  Assert (values[0].size() == generalized_support_points.size(),
+         ExcDimensionMismatch(values.size(), generalized_support_points.size()));
+  Assert (local_dofs.size() == this->dofs_per_cell,
+         ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+                                  // First do interpolation on
+                                  // faces. There, the component
+                                  // evaluated depends on the face
+                                  // direction and orientation.
+  unsigned int fbase = 0;
+  unsigned int f=0;
+  for (;f<GeometryInfo<dim>::faces_per_cell;
+       ++f, fbase+=this->dofs_per_face)
+    {
+      for (unsigned int i=0;i<this->dofs_per_face;++i)
+       {
+         local_dofs[fbase+i] = values[GeometryInfo<dim>::unit_normal_direction[f]][fbase+i];
+       }
+    }
+                                  // The remaining points form dim
+                                  // chunks, one for each component.
+  const unsigned int istep = (this->dofs_per_cell - fbase) / dim;
+  Assert ((this->dofs_per_cell - fbase) % dim == 0, ExcInternalError());
+
+  f = 0;
+  while (fbase < this->dofs_per_cell)
+    {
+      for (unsigned int i=0;i<istep;++i)
+       {
+         local_dofs[fbase+i] = values[f][fbase+i];
+       }
+      fbase+=istep;
+      ++f;
+    }
+  Assert (fbase == this->dofs_per_cell, ExcInternalError());
+}
+
+
+
 
 template <int dim>
 FiniteElement<dim> *
@@ -199,8 +304,8 @@ template <int dim>
 void
 FE_RaviartThomasNodal<dim>::initialize_unit_support_points (const unsigned int deg)
 {
-  this->unit_support_points.resize (this->dofs_per_cell);
-  this->unit_face_support_points.resize (this->dofs_per_face);
+  this->generalized_support_points.resize (this->dofs_per_cell);
+  this->generalized_face_support_points.resize (this->dofs_per_face);
   
   unsigned int current = 0;
 
@@ -216,11 +321,11 @@ FE_RaviartThomasNodal<dim>::initialize_unit_support_points (const unsigned int d
       Assert (face_points.n_quadrature_points == this->dofs_per_face,
              ExcInternalError());
       for (unsigned int k=0;k<this->dofs_per_face;++k)
-       this->unit_face_support_points[k] = face_points.point(k);
+       this->generalized_face_support_points[k] = face_points.point(k);
       Quadrature<dim> faces = QProjector<dim>::project_to_all_faces(face_points);
       for (unsigned int k=0;k<//faces.n_quadrature_points
                          this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;++k)
-       this->unit_support_points[k] = faces.point(k);
+       this->generalized_support_points[k] = faces.point(k);
 
       current = this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
     }
@@ -242,7 +347,7 @@ FE_RaviartThomasNodal<dim>::initialize_unit_support_points (const unsigned int d
                                                       ((d==1) ? low : high),
                                                       ((d==2) ? low : high));
       for (unsigned int k=0;k<quadrature->n_quadrature_points;++k)
-       this->unit_support_points[current++] = quadrature->point(k);
+       this->generalized_support_points[current++] = quadrature->point(k);
       delete quadrature;
     }
   Assert (current == this->dofs_per_cell, ExcInternalError());
@@ -254,44 +359,15 @@ void
 FE_RaviartThomasNodal<dim>::initialize_node_matrix ()
 {
   const unsigned int n_dofs = this->dofs_per_cell;
-
                                   // We use an auxiliary matrix in
                                   // this function. Therefore,
                                   // inverse_node_matrix is still
                                   // empty and shape_value_component
                                   // returns the 'raw' shape values.
-  FullMatrix<double> N(n_dofs, n_dofs);
-
-                                  // The curent node functional index
-  unsigned int current = 0;
-
-                                  // For each face and all quadrature
-                                  // points on it, the node value is
-                                  // the normal component of the
-                                  // shape function, possibly
-                                  // pointing in negative direction.
-  for (unsigned int face = 0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-    for (unsigned int k=0;k<this->dofs_per_face;++k)
-      {
-       for (unsigned int i=0;i<n_dofs;++i)
-         N(current,i) = this->shape_value_component(
-           i, this->unit_support_points[current],
-           GeometryInfo< dim >::unit_normal_direction[face]);
-       ++current;
-      }
-                                  // Interior degrees of freedom in each direction
-  const unsigned int n_cell = (n_dofs - current) / dim;
-  
-  for (unsigned int d=0;d<dim;++d)
-    for (unsigned int k=0;k<n_cell;++k)
-      {
-       for (unsigned int i=0;i<n_dofs;++i)
-         N(current,i) = this->shape_value_component(i, this->unit_support_points[current], d);
-       ++current;
-      }
-  Assert (current == n_dofs, ExcInternalError());
+  FullMatrix<double> M(n_dofs, n_dofs);
+  FETools::compute_node_matrix(M, *this);
   this->inverse_node_matrix.reinit(n_dofs, n_dofs);
-  this->inverse_node_matrix.invert(N);
+  this->inverse_node_matrix.invert(M);
 }
 
 

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.