]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move the generation of intergrid constraints to the library.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 14 Apr 2000 17:09:32 +0000 (17:09 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 14 Apr 2000 17:09:32 +0000 (17:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@2728 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_tools.h
deal.II/deal.II/source/dofs/dof_tools.cc

index 5bcb95a0a46bd1c6611c1734047b6c91712db638..985cf0b36d2413b9b485993996cb08f0e99c52de 100644 (file)
@@ -248,8 +248,9 @@ class DoFTools
                                      * considered.
                                      */
     template<int dim>
-    static void make_flux_sparsity_pattern (const DoFHandler<dim> &dof_handler,
-                                           SparsityPattern       &sparsity_pattern);
+    static void
+    make_flux_sparsity_pattern (const DoFHandler<dim> &dof_handler,
+                               SparsityPattern       &sparsity_pattern);
     
                                     /**
                                      * Make up the constraints which
@@ -353,22 +354,35 @@ class DoFTools
                                      * the number of components in the
                                      * finite element used by #dof#.
                                      */
-    template<int dim>
-    static void extract_dofs(const DoFHandler<dim> &dof_handler,
-                            const vector<bool>    &select,
-                            vector<bool>          &selected_dofs);
+    template <int dim>
+    static void
+    extract_dofs(const DoFHandler<dim> &dof_handler,
+                const vector<bool>    &select,
+                vector<bool>          &selected_dofs);
 
                                     /**
                                      * Do the same thing as #extract_dofs#
                                      * for one level of a multi-grid DoF
                                      * numbering.
                                      */
-    template<int dim>
-    static void extract_level_dofs(const unsigned int       level,
-                                  const MGDoFHandler<dim> &dof,
-                                  const vector<bool>      &select,
-                                  vector<bool>            &selected_dofs);
+    template <int dim>
+    static void
+    extract_level_dofs(const unsigned int       level,
+                      const MGDoFHandler<dim> &dof,
+                      const vector<bool>      &select,
+                      vector<bool>            &selected_dofs);
 
+
+    template <int dim>
+    static void
+    compute_intergrid_constraints (const DoFHandler<dim>              &coarse_grid,
+                                  const unsigned int                  coarse_component,
+                                  const DoFHandler<dim>              &fine_grid,
+                                  const unsigned int                  fine_component,
+                                  const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
+                                  ConstraintMatrix                   &constraints);
+
+                                  
                                     /**
                                      * Exception
                                      */
@@ -385,6 +399,10 @@ class DoFTools
                    << "is invalid with respect to the number "
                    << "of components in the finite element "
                    << "(" << arg2 << ")");
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcFiniteElementsDontMatch);
 };
 
 
index e7a841a06d3493a8b53e697112e75a22af3170e5..1a7cceb263bcdc0b251e8f5efddc250f798c780c 100644 (file)
@@ -14,6 +14,7 @@
 
 #include <grid/tria.h>
 #include <grid/tria_iterator.h>
+#include <grid/intergrid_map.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_constraints.h>
@@ -709,6 +710,433 @@ DoFTools::extract_level_dofs(const unsigned int       level,
 }
 
 
+
+
+template <int dim>
+void
+DoFTools::compute_intergrid_constraints (const DoFHandler<dim>              &coarse_grid,
+                                        const unsigned int                  coarse_component,
+                                        const DoFHandler<dim>              &fine_grid,
+                                        const unsigned int                  fine_component,
+                                        const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
+                                        ConstraintMatrix                   &constraints)
+{
+                                  // aliases to the finite elements
+                                  // used by the dof handlers:
+  const FiniteElement<dim> &coarse_fe = coarse_grid.get_fe(),
+                          &fine_fe   = fine_grid.get_fe();
+
+                                  // global numbers of dofs
+  const unsigned int n_coarse_dofs = coarse_grid.n_dofs(),
+                    n_fine_dofs   = fine_grid.n_dofs();
+
+                                  // local numbers of dofs
+  const unsigned int coarse_dofs_per_cell = coarse_fe.dofs_per_cell,
+                    fine_dofs_per_cell   = fine_fe.dofs_per_cell;
+
+                                  // alias the number of dofs per
+                                  // cell belonging to the
+                                  // coarse_component which is to be
+                                  // the restriction of the fine
+                                  // grid:
+  const unsigned int coarse_dofs_per_cell_component
+    = coarse_fe.base_element(coarse_fe.component_to_base(coarse_component)).dofs_per_cell;
+  
+  
+                                  // check whether component numbers
+                                  // are valid
+  Assert (coarse_component < coarse_fe.n_components(),
+         ExcInvalidComponent (coarse_component, coarse_fe.n_components()));
+  Assert (fine_component < fine_fe.n_components(),
+         ExcInvalidComponent (fine_component, fine_fe.n_components()));
+                                  // check whether respective finite
+                                  // elements are equal
+  Assert (coarse_fe.base_element (coarse_fe.component_to_base(coarse_component))
+         ==
+         fine_fe.base_element (fine_fe.component_to_base(fine_component)),
+         ExcFiniteElementsDontMatch());
+  
+
+/*
+ * From here on: the term `parameter' refers to the selected component
+ * on the coarse grid and its analogon on the fine grid. The naming of
+ * variables containing this term is due to the fact that
+ * `selected_component' is longer, but also due to the fact that the
+ * code of this function was initially written for a program where the
+ * component which we wanted to match between grids was actually the
+ * `parameter' variable.
+ *
+ * Likewise, the terms `parameter grid' and `state grid' refer to the
+ * coarse and fine grids, respectively.
+ *
+ * Changing the names of variables would in principle be a good idea,
+ * but would not make things simpler and would be another source of
+ * errors. If anyone feels like doing so: patches would be welcome!
+ */
+
+
+  
+                                  // set up vectors of cell-local
+                                  // data; each vector represents one
+                                  // degree of freedom of the
+                                  // coarse-grid variable in the
+                                  // fine-grid element
+  vector<Vector<double> > parameter_dofs (coarse_dofs_per_cell_component,
+                                         Vector<double>(fine_dofs_per_cell));
+                                  // for each coarse dof: find its
+                                  // position within the fine element
+                                  // and set this value to one in the
+                                  // respective vector (all other values
+                                  // are zero by construction)
+  for (unsigned int local_coarse_dof=0;
+       local_coarse_dof<coarse_dofs_per_cell_component;
+       ++local_coarse_dof)
+    parameter_dofs[local_coarse_dof]
+      (fine_fe.component_to_system_index (fine_component,local_coarse_dof)) = 1.;
+
+
+                                  // find out how many DoFs there are
+                                  // on the grids belonging to the
+                                  // components we want to match
+  unsigned int n_parameters_on_fine_grid=0;
+  if (true)
+    {
+                                      // have a flag for each dof on
+                                      // the fine grid and set it
+                                      // to true if this is an
+                                      // interesting dof. finally count
+                                      // how many true's there
+      vector<bool>          dof_is_interesting (fine_grid.n_dofs(), false);
+      vector<unsigned int>  local_dof_indices (fine_fe.dofs_per_cell);
+      
+      for (DoFHandler<dim>::active_cell_iterator cell=fine_grid.begin_active();
+          cell!=fine_grid.end(); ++cell)
+       {
+         cell->get_dof_indices (local_dof_indices);
+         for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
+           if (fine_fe.system_to_component_index(i).first == fine_component)
+             dof_is_interesting[local_dof_indices[i]] = true;
+       };
+
+      n_parameters_on_fine_grid = count (dof_is_interesting.begin(),
+                                           dof_is_interesting.end(),
+                                           true);
+    };
+
+                                  // vector to hold the representation of
+                                  // a single degree of freedom on the
+                                  // coarse grid (for the selected fe)
+                                  // on the fine grid
+  Vector<double> global_parameter_representation (n_fine_dofs);
+
+                                  // store the weights with which a dof
+                                  // on the parameter grid contributes
+                                  // to a dof on the fine grid. see the
+                                  // long doc below for more info
+                                  //
+                                  // allocate as many rows as there are
+                                  // parameter dofs on the coarse grid
+                                  // and as many columns as there are
+                                  // parameter dofs on the fine grid.
+                                  //
+                                  // weight_mapping is used to map the
+                                  // global (fine grid) parameter dof
+                                  // indices to the columns
+  FullMatrix<double> weights (n_coarse_dofs,
+                             n_parameters_on_fine_grid);
+                                  // this is this mapping. there is one
+                                  // entry for each dof on the fine grid;
+                                  // if it is a parameter dof, then its
+                                  // value is the column in weights for
+                                  // that parameter dof, if it is any
+                                  // other dof, then its value is -1,
+                                  // indicating an error
+  vector<int> weight_mapping (n_fine_dofs, -1);
+
+                                  // set up this mapping
+  if (true)
+    {
+                                      // list of local dof numbers
+                                      // belonging to the parameter
+                                      // part of the FE on the fine grid
+      vector<unsigned int> parameter_dofs;
+      for (unsigned int i=0; i<coarse_dofs_per_cell; ++i)
+       parameter_dofs.push_back(fine_fe.component_to_system_index(coarse_component,i));
+
+      vector<unsigned int> local_dof_indices(fine_fe.dofs_per_cell);
+      unsigned int next_free_index=0;
+      for (typename DoFHandler<dim>::active_cell_iterator cell=fine_grid.begin_active();
+          cell != fine_grid.end(); ++cell)
+       {
+         cell->get_dof_indices (local_dof_indices);
+         for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
+                                            // if this DoF is a
+                                            // parameter dof and has
+                                            // not yet been numbered,
+                                            // then do so
+           if ((fine_fe.system_to_component_index(i).first == fine_component) &&
+               (weight_mapping[local_dof_indices[i]] == -1))
+             {
+               weight_mapping[local_dof_indices[i]] = next_free_index;
+               ++next_free_index;
+             };
+       };
+
+      Assert (next_free_index == n_parameters_on_fine_grid,
+             ExcInternalError());
+    };
+  
+  
+                                  // for each cell on the parameter grid:
+                                  // find out which degrees of freedom on the
+                                  // fine grid correspond in which way to
+                                  // the degrees of freedom on the parameter
+                                  // grid
+                                  //
+                                  // since for continuous FEs some
+                                  // dofs exist on more than one
+                                  // cell, we have to track which
+                                  // ones were already visited. the
+                                  // problem is that if we visit a
+                                  // dof first on one cell and
+                                  // compute its weight with respect
+                                  // to some global dofs to be
+                                  // non-zero, and later visit the
+                                  // dof again on another cell and
+                                  // (since we are on another cell)
+                                  // recompute the weights with
+                                  // respect to the same dofs as
+                                  // above to be zero now, we have to
+                                  // preserve them. we therefore
+                                  // overwrite all weights if they
+                                  // are nonzero and do not enforce
+                                  // zero weights since that might be
+                                  // only due to the fact that we are
+                                  // on another cell.
+                                  //
+                                  // example:
+                                  // coarse grid
+                                  //  |     |     |
+                                  //  *-----*-----*
+                                  //  | cell|cell |
+                                  //  |  1  |  2  |
+                                  //  |     |     |
+                                  //  0-----1-----*
+                                  //
+                                  // fine grid
+                                  //  |  |  |  |  |
+                                  //  *--*--*--*--*
+                                  //  |  |  |  |  |
+                                  //  *--*--*--*--*
+                                  //  |  |  |  |  |
+                                  //  *--x--y--*--*
+                                  //
+                                  // when on cell 1, we compute the
+                                  // weights of dof 'x' to be 1/2
+                                  // from parameter dofs 0 and 1,
+                                  // respectively. however, when
+                                  // later we are on cell 2, we again
+                                  // compute the prolongation of
+                                  // shape function 1 restricted to
+                                  // cell 2 to the globla grid and
+                                  // find that the weight of global
+                                  // dof 'x' now is zero. however, we
+                                  // should not overwrite the old
+                                  // value.
+                                  //
+                                  // we therefore always only set
+                                  // nonzero values. why adding up is
+                                  // not useful: dof 'y' would get
+                                  // weight 1 from parameter dof 1 on
+                                  // both cells 1 and 2, but the
+                                  // correct weight is nevertheless
+                                  // only 1.
+  typename DoFHandler<dim>::active_cell_iterator cell;
+  vector<unsigned int> parameter_dof_indices (coarse_fe.dofs_per_cell);
+  
+  for (cell=coarse_grid.begin_active(); cell!=coarse_grid.end(); ++cell)
+    {
+                                      // get the global indices of the
+                                      // parameter dofs on this parameter
+                                      // grid cell
+      cell->get_dof_indices (parameter_dof_indices);
+
+                                      // loop over all dofs on this
+                                      // cell and check whether they
+                                      // are interesting for us
+      for (unsigned int local_parameter_dof=0;
+          local_parameter_dof<coarse_fe.dofs_per_cell;
+          ++local_parameter_dof)
+       if (coarse_fe.system_to_component_index(local_parameter_dof).first
+           ==
+           coarse_component)
+         {
+           global_parameter_representation.clear ();
+           
+                                            // distribute the representation of
+                                            // #local_parameter_dof# on the
+                                            // parameter grid cell #cell# to
+                                            // the global data space
+           coarse_to_fine_grid_map[cell]->
+             set_dof_values_by_interpolation (parameter_dofs[local_parameter_dof],
+                                              global_parameter_representation);
+                                            // now that we've got the global
+                                            // representation of each parameter
+                                            // dof, we've only got to clobber the
+                                            // non-zero entries in that vector and
+                                            // store the result
+                                            //
+                                            // what we have learned: if entry #i#
+                                            // of the global vector holds the value
+                                            // #v[i]#, then this is the weight with
+                                            // which the present dof contributes
+                                            // to #i#. there may be several such
+                                            // #i#s and their weights' sum should
+                                            // be one. Then, #v[i]# should be
+                                            // equal to #\sum_j w_{ij} p[j]# with
+                                            // #p[j]# be the values of the degrees
+                                            // of freedom on the coarse grid. we
+                                            // can thus compute constraints which
+                                            // link the degrees of freedom #v[i]#
+                                            // on the fine grid to those on the
+                                            // coarse grid, #p[j]#. Now to use
+                                            // these as real constraints, rather
+                                            // than as additional equations, we
+                                            // have to identify representants
+                                            // among the #i# for each #j#. this will
+                                            // be done by simply taking the first
+                                            // #i# for which #w_{ij}==1#.
+           for (unsigned int i=0; i<global_parameter_representation.size(); ++i)
+                                              // set this weight if it belongs
+                                              // to a parameter dof.
+             if (weight_mapping[i] != -1)
+               {
+                                                  // only overwrite old
+                                                  // value if not by
+                                                  // zero
+                 if (global_parameter_representation(i) != 0)
+                   {
+                     const unsigned int wi = parameter_dof_indices[local_parameter_dof],
+                                        wj = weight_mapping[i];
+                     weights(wi,wj) = global_parameter_representation(i);
+                   };
+               }
+             else
+               Assert (global_parameter_representation(i) == 0,
+                       ExcInternalError());
+         };
+    };
+  
+                                  // ok, now we have all weights for each
+                                  // dof on the fine grid. if in debug
+                                  // mode lets see if everything went smooth,
+                                  // i.e. each dof has sum of weights one
+                                  //
+                                  // in other words this means that
+                                  // if the sum of all shape
+                                  // functions on the parameter grid
+                                  // is one (which is always the
+                                  // case), then the representation
+                                  // on the state grid should be as
+                                  // well (division of unity)
+                                  //
+                                  // if the parameter grid has more
+                                  // than one component, then the
+                                  // respective dofs of the other
+                                  // components have sum of weights
+                                  // zero, of course. we do not
+                                  // explicitely ask which component
+                                  // a dof belongs to, but this at
+                                  // least tests some errors
+#ifdef DEBUG
+  for (unsigned int col=0; col<weights.n(); ++col)
+    {
+      double sum=0;
+      for (unsigned int row=0; row<weights.m(); ++row)
+       sum += weights(row,col);
+      Assert ((sum==1) ||
+             ((coarse_fe.n_components()>1) && (sum==0)), ExcInternalError());
+    };
+#endif
+
+
+                                  // now we know that the weights in each
+                                  // row constitute a constraint. enter this
+                                  // into the constraints object
+                                  //
+                                  // first task: for each parameter dof on
+                                  // the parameter grid, find a representant
+                                  // on the fine, global grid. this is
+                                  // possible since we use conforming
+                                  // finite element. we take this representant
+                                  // to be the first element in this row with
+                                  // weight identical to one
+  vector<int> representants(weights.m(), -1);
+  for (unsigned int parameter_dof=0; parameter_dof<weights.m(); ++parameter_dof)
+    {
+      unsigned int column=0;
+      for (; column<weights.n(); ++column)
+       if (weights(parameter_dof,column) == 1)
+         break;
+      Assert (column < weights.n(), ExcInternalError());
+
+                                      // now we know in which column of
+                                      // weights the representant is, but
+                                      // we don't know its global index. get
+                                      // it using the inverse operation of
+                                      // the weight_mapping
+      unsigned int global_dof=0;
+      for (; global_dof<weight_mapping.size(); ++global_dof)
+       if (weight_mapping[global_dof] == static_cast<int>(column))
+         break;
+      Assert (global_dof < weight_mapping.size(), ExcInternalError());
+
+                                      // now enter the representants global
+                                      // index into our list
+      representants[parameter_dof] = global_dof;
+    };
+  
+  for (unsigned int global_dof=0; global_dof<n_fine_dofs; ++global_dof)
+    if (weight_mapping[global_dof] != -1)
+                                      // this global dof is a parameter
+                                      // dof, so it may carry a constraint
+                                      // note that for each global dof,
+                                      // the sum of weights shall be one,
+                                      // so we can find out whether this
+                                      // dof is constrained in the following
+                                      // way: if the only weight in this row
+                                      // is a one, and the representant for
+                                      // the parameter dof of the line in
+                                      // which this one is is the present
+                                      // dof, then we consider this dof
+                                      // to be unconstrained. otherwise,
+                                      // all other dofs are constrained
+      {
+       unsigned int first_used_row=0;
+       for (; first_used_row<weights.m(); ++first_used_row)
+         if (weights(first_used_row,weight_mapping[global_dof]) != 0)
+           break;
+
+       if ((weights(first_used_row,weight_mapping[global_dof]) == 1) &&
+           (representants[first_used_row] == static_cast<int>(global_dof)))
+                                          // dof unconstrained
+         continue;
+
+                                        // otherwise enter all constraints
+       constraints.add_line (global_dof);
+       for (unsigned int row=first_used_row; row<weights.m(); ++row)
+         if (weights(row,weight_mapping[global_dof]) != 0)
+           constraints.add_entry (global_dof,
+                                  representants[row],
+                                  weights(row,weight_mapping[global_dof]));
+      };
+};
+
+
+
+
+
+
 // explicit instantiations
 #if deal_II_dimension > 1
 template void
@@ -755,3 +1183,12 @@ template void DoFTools::extract_level_dofs(unsigned int level,
                                           const vector<bool>& local_select,
                                           vector<bool>& selected_dofs);
 
+
+template
+void
+DoFTools::compute_intergrid_constraints (const DoFHandler<deal_II_dimension> &,
+                                        const unsigned int                   ,
+                                        const DoFHandler<deal_II_dimension> &,
+                                        const unsigned int                   ,
+                                        const InterGridMap<DoFHandler,deal_II_dimension> &,
+                                        ConstraintMatrix                    &);

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.