]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
New interpolate, back_interpolate and interpolation_difference functions with Constra...
authorhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 13 Jun 2002 12:49:10 +0000 (12:49 +0000)
committerhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 13 Jun 2002 12:49:10 +0000 (12:49 +0000)
git-svn-id: https://svn.dealii.org/trunk@6092 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_tools.h
deal.II/deal.II/source/fe/fe_tools.cc
deal.II/doc/news/2002/c-3-4.html

index f01d0bb95efc0027176674cdded96af7a3642e16..94a6a8c14e1a686ce38001a4551a4df47b770758 100644 (file)
@@ -20,6 +20,8 @@ template <int dim> class FiniteElement;
 template <int dim> class DoFHandler;
 template <typename number> class Vector;
 template <int dim> class FE_Q;
+class ConstraintMatrix;
+
 
 #include <base/config.h>
 #include <base/exceptions.h>
@@ -127,7 +129,12 @@ class FETools
                                      * at the DoF values on the
                                      * discontinuities.
                                      *
-                                     * Note that the resulting output
+                                     * Note that for continuous
+                                     * elements on grids with hanging
+                                     * nodes (i.e. locally refined
+                                     * grids) this function does not
+                                     * give the expected output.
+                                     * Indeed, the resulting output
                                      * vector does not necessarily
                                      * respect continuity
                                      * requirements at hanging nodes:
@@ -144,12 +151,17 @@ class FETools
                                      * triangulation, although it is
                                      * of course Q1 on each cell.
                                      *
-                                     * To make the field conforming
-                                     * again, you have to have a
-                                     * hanging node constraints
-                                     * object, and call its
-                                     * @p{distribute} function on the
-                                     * result of this function.
+                                     * For this case (continuous
+                                     * elements on grids with hanging
+                                     * nodes), please use the
+                                     * @p{interpolate} function with
+                                     * an additional
+                                     * @p{ConstraintMatrix} argument,
+                                     * see below, or make the field
+                                     * conforming yourself by calling
+                                     * the @p{distribute} function of
+                                     * your hanging node constraints
+                                     * object.
                                      */
     template <int dim, typename number>
     static void interpolate (const DoFHandler<dim> &dof1,
@@ -157,6 +169,47 @@ class FETools
                             const DoFHandler<dim> &dof2,
                             Vector<number>        &u2);
     
+                                    /**
+                                     * Gives the interpolation of a
+                                     * the @p{dof1}-function @p{u1}
+                                     * to a @p{dof2}-function
+                                     * @p{u2}. @p{dof1} and @p{dof2}
+                                     * need to be @ref{DoFHandler}s
+                                     * based on the same
+                                     * triangulation.
+                                     * @p{constraints} is a hanging
+                                     * node constraints object
+                                     * corresponding to
+                                     * @p{dof2}. This object is
+                                     * particular important when
+                                     * interpolating onto continuous
+                                     * elements on grids with hanging
+                                     * nodes (locally refined grids).
+                                     *
+                                     * If the elements @p{fe1} and @p{fe2}
+                                     * are either both continuous or
+                                     * both discontinuous then this
+                                     * interpolation is the usual point
+                                     * interpolation. The same is true
+                                     * if @p{fe1} is a continuous and
+                                     * @p{fe2} is a discontinuous finite
+                                     * element. For the case that @p{fe1}
+                                     * is a discontinuous and @p{fe2} is
+                                     * a continuous finite element
+                                     * there is no point interpolation
+                                     * defined at the discontinuities.
+                                     * Therefore the meanvalue is taken
+                                     * at the DoF values on the
+                                     * discontinuities.
+                                     */
+    template <int dim, typename number>
+    static void interpolate (const DoFHandler<dim>  &dof1,
+                            const Vector<number>   &u1,
+                            const DoFHandler<dim>  &dof2,
+                            const ConstraintMatrix &constraints,
+                            Vector<number>         &u2);
+    
+
                                     /**
                                      * Gives the interpolation of the
                                      * @p{fe1}-function @p{u1} to a
@@ -165,23 +218,63 @@ class FETools
                                      * @p{fe1}-function named
                                      * @p{u1_interpolated}.
                                      *
-                                     * Note, that this function only
-                                     * makes sense if the finite
-                                     * element space due to @p{fe1}
-                                     * is not a subset of the finite
-                                     * element space due to @p{fe2},
-                                     * as if it were a subset then
-                                     * @p{u1_interpolated} would be
-                                     * equal to @p{u1}.
-                                     *
                                      * Note, that this function does
                                      * not work on continuous
-                                     * elements at hanging nodes.
+                                     * elements at hanging nodes. For
+                                     * that case use the
+                                     * @p{back_interpolate} function,
+                                     * below, that takes an
+                                     * additional
+                                     * @p{ConstraintMatrix} object.
+                                     *
+                                     * Furthermore note, that for the
+                                     * specific case when the finite
+                                     * element space corresponding to
+                                     * @p{fe1} is a subset of the
+                                     * finite element space
+                                     * corresponding to @p{fe2}, this
+                                     * function is simply an identity
+                                     * mapping.
                                      */
     template <int dim, typename number>
     static void back_interpolate (const DoFHandler<dim>    &dof1,
                                  const Vector<number>     &u1,
                                  const FiniteElement<dim> &fe2,
+                                 Vector<number>           &u1_interpolated);
+
+                                    /**
+                                     * Gives the interpolation of the
+                                     * @p{dof1}-function @p{u1} to a
+                                     * @p{dof2}-function, and
+                                     * interpolates this to a second
+                                     * @p{dof1}-function named
+                                     * @p{u1_interpolated}.
+                                     * @p{constraints1} and
+                                     * @p{constraints2} are the
+                                     * hanging node constraints
+                                     * corresponding to @p{dof1} and
+                                     * @p{dof2}, respectively. These
+                                     * objects are particular
+                                     * important when continuous
+                                     * elements on grids with hanging
+                                     * nodes (locally refined grids)
+                                     * are involved.
+                                     *
+                                     * Furthermore note, that for the
+                                     * specific case when the finite
+                                     * element space corresponding to
+                                     * @p{dof1} is a subset of the
+                                     * finite element space
+                                     * corresponding to @p{dof2}, this
+                                     * function is simply an identity
+                                     * mapping.
+                                     */
+    template <int dim, typename number>
+    static void back_interpolate (const DoFHandler<dim>    &dof1,
+                                 const ConstraintMatrix   &constraints1,
+                                 const Vector<number>     &u1,
+                                 const DoFHandler<dim>    &dof2,
+                                 const ConstraintMatrix   &constraints2,
                                  Vector<number>           &u1_interpolated);
 
                                     /**
@@ -193,7 +286,12 @@ class FETools
                                      *
                                      * Note, that this function does
                                      * not work on continuous
-                                     * elements at hanging nodes.
+                                     * elements at hanging nodes. For
+                                     * that case use the
+                                     * @p{interpolation_difference}
+                                     * function, below, that takes an
+                                     * additional
+                                     * @p{ConstraintMatrix} object.
                                      */
     template <int dim, typename number>
     static void interpolation_difference(const DoFHandler<dim> &dof1,
@@ -201,6 +299,31 @@ class FETools
                                         const FiniteElement<dim> &fe2,
                                         Vector<number> &z1_difference);    
     
+                                    /**
+                                     * Gives $(Id-I_h)z1$ for a given
+                                     * @p{dof1}-function @p{z1}, where $I_h$
+                                     * is the interpolation from @p{fe1}
+                                     * to @p{fe2}. $(Id-I_h)z1$ is
+                                     * denoted by @p{z1_difference}.
+                                     * @p{constraints1} and
+                                     * @p{constraints2} are the
+                                     * hanging node constraints
+                                     * corresponding to @p{dof1} and
+                                     * @p{dof2}, respectively. These
+                                     * objects are particular
+                                     * important when continuous
+                                     * elements on grids with hanging
+                                     * nodes (locally refined grids)
+                                     * are involved.
+                                     */
+    template <int dim, typename number>
+    static void interpolation_difference(const DoFHandler<dim> &dof1,
+                                        const ConstraintMatrix &constraints1,
+                                        const Vector<number> &z1,
+                                        const DoFHandler<dim> &dof2,
+                                        const ConstraintMatrix &constraints2,
+                                        Vector<number> &z1_difference);    
+
                                     /**
                                      * Gives the patchwise
                                      * extrapolation of a @p{dof1}
@@ -306,6 +429,11 @@ class FETools
                                      * Exception
                                      */
     DeclException0 (ExcTriangulationMismatch);
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcHangingNodesNotAllowed);
   
                                     /**
                                      * Exception
index e02a3db5093e68fa49a5e6f01d05c8ed8fc4fb20..28d615f943a9e495a9692f2d21ed2294fa311447 100644 (file)
@@ -22,6 +22,7 @@
 #include <fe/fe_q.h>
 #include <fe/fe_values.h>
 #include <fe/mapping_q1.h>
+#include <dofs/dof_constraints.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
 
@@ -122,11 +123,21 @@ void FETools::get_interpolation_difference_matrix(const FiniteElement<dim> &fe1,
 }
 
 
+template <int dim, typename number>
+void FETools::interpolate(const DoFHandler<dim> &dof1,
+                         const Vector<number> &u1,
+                         const DoFHandler<dim> &dof2,
+                         Vector<number> &u2)
+{
+  interpolate(dof1, u1, dof2, ConstraintMatrix(), u2);
+}
 
+  
 template <int dim, typename number>
 void FETools::interpolate(const DoFHandler<dim> &dof1,
                          const Vector<number> &u1,
                          const DoFHandler<dim> &dof2,
+                         const ConstraintMatrix &constraints,
                          Vector<number> &u2)
 {
   Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch());
@@ -135,6 +146,16 @@ void FETools::interpolate(const DoFHandler<dim> &dof1,
   Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));  
   Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs()));
 
+                                  // for continuous elements on grids
+                                  // with hanging nodes we need
+                                  // hanging node
+                                  // constraints. Consequentely, if
+                                  // there are no constraints then
+                                  // hanging nodes are not allowed.
+  const bool hanging_nodes_not_allowed=
+    (dof2.get_fe().dofs_per_vertex != 0)
+    && (constraints.n_constraints()==0);
+  
   const unsigned int dofs_per_cell1=dof1.get_fe().dofs_per_cell;
   const unsigned int dofs_per_cell2=dof2.get_fe().dofs_per_cell;
   
@@ -154,9 +175,15 @@ void FETools::interpolate(const DoFHandler<dim> &dof1,
   std::vector<unsigned int> index_multiplicity(dof2.n_dofs(),0);
   std::vector<unsigned int> dofs (dofs_per_cell2);
   u2.clear ();
-  
+
   for (; cell1!=endc1, cell2!=endc2; ++cell1, ++cell2) 
     {
+      if (hanging_nodes_not_allowed)
+       for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+         Assert (cell1->at_boundary(face) ||
+                 cell1->neighbor(face)->level() == cell1->level(),
+                 ExcHangingNodesNotAllowed());
+      
       cell1->get_dof_values(u1, u1_local);
       interpolation_matrix.vmult(u2_local, u1_local);
       cell2->get_dof_indices(dofs);
@@ -166,13 +193,19 @@ void FETools::interpolate(const DoFHandler<dim> &dof1,
          ++index_multiplicity[dofs[i]];
        }
     }
+
+                                  // when a discontinuous element is
+                                  // interpolated to a continuous
+                                  // one, we take the mean values.
   for (unsigned int i=0; i<dof2.n_dofs(); ++i)
     {
       Assert(index_multiplicity[i]!=0, ExcInternalError());
       u2(i) /= index_multiplicity[i];
     }
-}
 
+                                  // Apply hanging node constraints.
+  constraints.distribute(u2);
+}
 
 
 template <int dim, typename number>
@@ -185,30 +218,83 @@ void FETools::back_interpolate(const DoFHandler<dim> &dof1,
         ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components()));
   Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));  
   Assert(u1_interpolated.size()==dof1.n_dofs(),
-        ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs()));  
+        ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs()));
+  
+                                  // For continuous elements on grids
+                                  // with hanging nodes we need
+                                  // hnaging node
+                                  // constraints. Consequently, when
+                                  // the elements are continuous no
+                                  // hanging node constraints are
+                                  // allowed.
+  const bool hanging_nodes_not_allowed=
+    (dof1.get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0);
 
   const unsigned int dofs_per_cell1=dof1.get_fe().dofs_per_cell;
 
   Vector<number> u1_local(dofs_per_cell1);
   Vector<number> u1_int_local(dofs_per_cell1);
   
+  typename DoFHandler<dim>::active_cell_iterator cell = dof1.begin_active(),
+                                                endc = dof1.end();
+
   FullMatrix<number> interpolation_matrix(dofs_per_cell1, dofs_per_cell1);
   FETools::get_back_interpolation_matrix(dof1.get_fe(), fe2,
                                         interpolation_matrix);
-
-  typename DoFHandler<dim>::active_cell_iterator cell1 = dof1.begin_active(),
-                                                endc1 = dof1.end();
-  
-  for (; cell1!=endc1; ++cell1) 
+  for (; cell!=endc; ++cell) 
     {
-      cell1->get_dof_values(u1, u1_local);
+      if (hanging_nodes_not_allowed)
+       for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+         Assert (cell->at_boundary(face) ||
+                 cell->neighbor(face)->level() == cell->level(),
+                 ExcHangingNodesNotAllowed());
+
+      cell->get_dof_values(u1, u1_local);
       interpolation_matrix.vmult(u1_int_local, u1_local);
-      cell1->set_dof_values(u1_int_local, u1_interpolated);
+      cell->set_dof_values(u1_int_local, u1_interpolated);
     }
 }
 
 
+  
+template <int dim, typename number>
+void FETools::back_interpolate(const DoFHandler<dim> &dof1,
+                              const ConstraintMatrix &constraints1,
+                              const Vector<number> &u1,
+                              const DoFHandler<dim> &dof2,
+                              const ConstraintMatrix &constraints2,
+                              Vector<number> &u1_interpolated)
+{
+                                  // For discontinuous elements
+                                  // without constraints take the
+                                  // simpler version of the
+                                  // back_interpolate function.
+  if (dof1.get_fe().dofs_per_vertex==0 && dof2.get_fe().dofs_per_vertex==0
+      && constraints1.n_constraints()==0 && constraints2.n_constraints()==0)
+    back_interpolate(dof1, u1, dof2.get_fe(), u1_interpolated);
+  else
+    {
+      Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(),
+            ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components()));
+      Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));  
+      Assert(u1_interpolated.size()==dof1.n_dofs(),
+            ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs()));
+      
+                                      // For continuous elements
+                                      // first interpolate to dof2,
+                                      // taking into account
+                                      // constraints2, and then
+                                      // interpolate back to dof1
+                                      // taking into account
+                                      // constraints1
+      Vector<double> u2(dof2.n_dofs());
+      interpolate(dof1, u1, dof2, constraints2, u2);
+      interpolate(dof2, u2, dof1, constraints1, u1_interpolated);
+    }
+}
 
+
+  
 template <int dim, typename number>
 void FETools::interpolation_difference(const DoFHandler<dim> &dof1,
                                       const Vector<number> &u1,
@@ -219,30 +305,69 @@ void FETools::interpolation_difference(const DoFHandler<dim> &dof1,
         ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components()));
   Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));  
   Assert(u1_difference.size()==dof1.n_dofs(),
-        ExcDimensionMismatch(u1_difference.size(), dof1.n_dofs()));  
-
-  const unsigned int dofs_per_cell1=dof1.get_fe().dofs_per_cell;
-
-  Vector<number> u1_local(dofs_per_cell1);
-  Vector<number> u1_diff_local(dofs_per_cell1);
+        ExcDimensionMismatch(u1_difference.size(), dof1.n_dofs()));
+  
+                                  // For continuous elements on grids
+                                  // with hanging nodes we need
+                                  // hnaging node
+                                  // constraints. Consequently, when
+                                  // the elements are continuous no
+                                  // hanging node constraints are
+                                  // allowed.
+  const bool hanging_nodes_not_allowed=
+    (dof1.get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0);
+
+  const unsigned int dofs_per_cell=dof1.get_fe().dofs_per_cell;
+
+  Vector<number> u1_local(dofs_per_cell);
+  Vector<number> u1_diff_local(dofs_per_cell);
   
-  FullMatrix<number> difference_matrix(dofs_per_cell1, dofs_per_cell1);
+  FullMatrix<number> difference_matrix(dofs_per_cell, dofs_per_cell);
   FETools::get_interpolation_difference_matrix(dof1.get_fe(), fe2,
                                               difference_matrix);
   
-  typename DoFHandler<dim>::active_cell_iterator cell1 = dof1.begin_active(),
-                                                endc1 = dof1.end();
+  typename DoFHandler<dim>::active_cell_iterator cell = dof1.begin_active(),
+                                                endc = dof1.end();
   
-  for (; cell1!=endc1; ++cell1) 
+  for (; cell!=endc; ++cell)
     {
-      cell1->get_dof_values(u1, u1_local);
+      if (hanging_nodes_not_allowed)
+       for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+         Assert (cell->at_boundary(face) ||
+                 cell->neighbor(face)->level() == cell->level(),
+                 ExcHangingNodesNotAllowed());
+
+      cell->get_dof_values(u1, u1_local);
       difference_matrix.vmult(u1_diff_local, u1_local);
-      cell1->set_dof_values(u1_diff_local, u1_difference);
+      cell->set_dof_values(u1_diff_local, u1_difference);
     }
 }
 
 
 
+template <int dim, typename number>
+void FETools::interpolation_difference(const DoFHandler<dim> &dof1,
+                                      const ConstraintMatrix &constraints1,
+                                      const Vector<number> &u1,
+                                      const DoFHandler<dim> &dof2,
+                                      const ConstraintMatrix &constraints2,
+                                      Vector<number> &u1_difference)
+{
+                                  // For discontinuous elements
+                                  // without constraints take the
+                                  // cheaper version of the
+                                  // interpolation_difference function.
+  if (dof1.get_fe().dofs_per_vertex==0 && dof2.get_fe().dofs_per_vertex==0
+      && constraints1.n_constraints()==0 && constraints2.n_constraints()==0)
+    interpolation_difference(dof1, u1, dof2.get_fe(), u1_difference);
+  else
+    {
+      back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference);
+      u1_difference.sadd(-1, u1);
+    }
+}
+
+  
 template <int dim, typename number>
 void FETools::extrapolate(const DoFHandler<dim> &dof1,
                          const Vector<number> &u1,
@@ -732,18 +857,41 @@ void FETools::interpolate<deal_II_dimension>
  const DoFHandler<deal_II_dimension> &,
  Vector<double> &);
 template
+void FETools::interpolate<deal_II_dimension>
+(const DoFHandler<deal_II_dimension> &,
+ const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &,
+ const ConstraintMatrix &,
+ Vector<double> &);
+template
 void FETools::back_interpolate<deal_II_dimension>
 (const DoFHandler<deal_II_dimension> &,
  const Vector<double> &,
  const FiniteElement<deal_II_dimension> &,
  Vector<double> &);
 template
+void FETools::back_interpolate<deal_II_dimension>
+(const DoFHandler<deal_II_dimension> &,
+ const ConstraintMatrix &, 
+ const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &,
+ const ConstraintMatrix &,
+ Vector<double> &);
+template
 void FETools::interpolation_difference<deal_II_dimension>
 (const DoFHandler<deal_II_dimension> &,
  const Vector<double> &,
  const FiniteElement<deal_II_dimension> &,
  Vector<double> &);
 template
+void FETools::interpolation_difference<deal_II_dimension>
+(const DoFHandler<deal_II_dimension> &,
+ const ConstraintMatrix &,
+ const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &,
+ const ConstraintMatrix &,
+ Vector<double> &);
+template
 void FETools::extrapolate<deal_II_dimension>
 (const DoFHandler<deal_II_dimension> &,
  const Vector<double> &,
index b14d5b9d51602629be0ca8082652921f8299d07e..50c6aa2f5a6e142f6b9cd0a8d1b762710b696697 100644 (file)
@@ -105,6 +105,22 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p> New: The existing <code
+       class="member">FETools::interpolate</code>, <code
+       class="member">FETools::back_interpolate</code> and <code
+       class="member">FETools::interpolate_difference</code> functions
+       do not take into account hanging node constraints. Therefore,
+       they have been working for continuous elements on globally
+       refined grids and on discontinuous elements, only. Now, there
+       are new functions with the same names but additional <code
+       class="class">ConstraintMatrix</code> arguments for the hanging
+       node constraints. These new functions now work for continuous
+       as well as for discontinuous elements. But, the old functions
+       are still supported.
+       <br>
+       (RH 2002/06/13)
+       </p>
+
   <li> <p> Change: The constructor of <code
        class="class">DoFHandler</code> now takes a reference to a
        <tt>const</tt> <code class="class">Triangulation</code>.

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.