]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Clean up a good amount of code. In particular, do not try to access fe_indices of...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 25 Jan 2014 04:11:46 +0000 (04:11 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 25 Jan 2014 04:11:46 +0000 (04:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@32310 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/source/numerics/solution_transfer.cc

index 54433972bb81d81a8a0230840e73531f597868a2..d15d6768dd7f15efd7f28ec58dd48c80109c22ef 100644 (file)
@@ -234,9 +234,9 @@ void
 SolutionTransfer<dim, VECTOR, DH>::
 prepare_for_coarsening_and_refinement(const std::vector<VECTOR> &all_in)
 {
-  Assert(prepared_for!=pure_refinement, ExcAlreadyPrepForRef());
-  Assert(prepared_for!=coarsening_and_refinement,
-         ExcAlreadyPrepForCoarseAndRef());
+  Assert (prepared_for!=pure_refinement, ExcAlreadyPrepForRef());
+  Assert (prepared_for!=coarsening_and_refinement,
+          ExcAlreadyPrepForCoarseAndRef());
 
   const unsigned int in_size=all_in.size();
   Assert(in_size!=0, ExcNoInVectorsGiven());
@@ -244,7 +244,7 @@ prepare_for_coarsening_and_refinement(const std::vector<VECTOR> &all_in)
   clear();
 
   const unsigned int n_active_cells = dof_handler->get_tria().n_active_cells();
-  n_dofs_old=dof_handler->n_dofs();
+  n_dofs_old = dof_handler->n_dofs();
 
   for (unsigned int i=0; i<in_size; ++i)
     {
@@ -253,14 +253,12 @@ prepare_for_coarsening_and_refinement(const std::vector<VECTOR> &all_in)
     }
 
   // first count the number
-  // of cells that'll be coarsened
+  // of cells that will be coarsened
   // and that'll stay or be refined
   unsigned int n_cells_to_coarsen=0;
   unsigned int n_cells_to_stay_or_refine=0;
-  typename DH::active_cell_iterator
-  act_cell = dof_handler->begin_active(),
-  endc = dof_handler->end();
-  for (; act_cell!=endc; ++act_cell)
+  for (typename DH::active_cell_iterator act_cell = dof_handler->begin_active();
+      act_cell!=dof_handler->end(); ++act_cell)
     {
       if (act_cell->coarsen_flag_set())
         ++n_cells_to_coarsen;
@@ -271,20 +269,16 @@ prepare_for_coarsening_and_refinement(const std::vector<VECTOR> &all_in)
          ExcInternalError());
 
   unsigned int n_coarsen_fathers=0;
-  typename DH::cell_iterator
-  cell=dof_handler->begin();
-  for (; cell!=endc; ++cell)
+  for (typename DH::cell_iterator cell=dof_handler->begin();
+      cell!=dof_handler->end(); ++cell)
     if (!cell->active() && cell->child(0)->coarsen_flag_set())
       ++n_coarsen_fathers;
-
-  if (n_cells_to_coarsen)
-    Assert(n_cells_to_coarsen>=2*n_coarsen_fathers, ExcInternalError());
+  Assert(n_cells_to_coarsen>=2*n_coarsen_fathers, ExcInternalError());
 
   // allocate the needed memory. initialize
   // the following arrays in an efficient
   // way, without copying much
-  std::vector<std::vector<types::global_dof_index> >
-  (n_cells_to_stay_or_refine)
+  std::vector<std::vector<types::global_dof_index> >(n_cells_to_stay_or_refine)
   .swap(indices_on_cell);
 
   std::vector<std::vector<Vector<typename VECTOR::value_type> > >
@@ -293,19 +287,21 @@ prepare_for_coarsening_and_refinement(const std::vector<VECTOR> &all_in)
   .swap(dof_values_on_cell);
 
   typename VECTOR::value_type zero_val = typename VECTOR::value_type();
+
   Table<2,FullMatrix<double> > interpolation_hp;
   std::vector<std::vector<bool> > restriction_is_additive;
+
   internal::extract_interpolation_matrices (*dof_handler, interpolation_hp);
   internal::restriction_additive (dof_handler->get_fe(), restriction_is_additive);
-  Vector<typename VECTOR::value_type> tmp, tmp2;
 
   // we need counters for
   // the 'to_stay_or_refine' cells 'n_sr' and
   // the 'coarsen_fathers' cells 'n_cf',
   unsigned int n_sr=0, n_cf=0;
-  cell = dof_handler->begin();
-  for (; cell!=endc; ++cell)
+  for (typename DH::cell_iterator cell=dof_handler->begin();
+      cell!=dof_handler->end(); ++cell)
     {
+      // CASE 1: active cell that remains as it is
       if (cell->active() && !cell->coarsen_flag_set())
         {
           const unsigned int dofs_per_cell=cell->get_fe().dofs_per_cell;
@@ -316,9 +312,11 @@ prepare_for_coarsening_and_refinement(const std::vector<VECTOR> &all_in)
           // interpolating to the children
           cell->get_dof_indices(indices_on_cell[n_sr]);
           cell_map[std::make_pair(cell->level(), cell->index())]
-            = Pointerstruct(&indices_on_cell[n_sr],cell->active_fe_index());
+            = Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index());
           ++n_sr;
         }
+
+      // CASE 2: cell is inactive but will become active
       else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
         {
           // note that if one child has the
@@ -327,127 +325,51 @@ prepare_for_coarsening_and_refinement(const std::vector<VECTOR> &all_in)
           // worked correctly
           for (unsigned int i=1; i<cell->n_children(); ++i)
             Assert(cell->child(i)->coarsen_flag_set(),
-                   ExcTriaPrepCoarseningNotCalledBefore());
-          const unsigned int dofs_per_cell=cell->get_fe().dofs_per_cell;
-
-          std::vector<Vector<typename VECTOR::value_type> >(in_size,
-                                                            Vector<typename VECTOR::value_type>(dofs_per_cell))
-          .swap(dof_values_on_cell[n_cf]);
+                   ExcMessage("It looks like you didn't call "
+                              "Triangulation::prepare_coarsening_and_refinement before "
+                              "calling the current function. This can't work."));
+
+          // we will need to interpolate from the children of this cell
+          // to the current one. in the hp context, this also means
+          // we need to figure out which finite element space to interpolate
+          // to since that is not implied by the global FE as in the non-hp
+          // case.
+          bool different_fe_on_children = false;
+          for (unsigned int child=1; child<cell->n_children(); ++child)
+            if (cell->child(child)->active_fe_index()
+                != cell->child(0)->active_fe_index())
+              {
+                different_fe_on_children = true;
+                break;
+              }
 
-          unsigned int fe_index = cell->active_fe_index();
+          // take FE index from the child with most
+          // degrees of freedom locally
           unsigned int most_general_child = 0;
-          bool different_elements = false;
-          for (unsigned int child=0; child<cell->n_children(); ++child)
-            {
-              if (cell->child(child)->active_fe_index() != fe_index)
-                different_elements = true;
-              // take FE index from the child with most
-              // degrees of freedom locally
+          if (different_fe_on_children == true)
+            for (unsigned int child=1; child<cell->n_children(); ++child)
               if (cell->child(child)->get_fe().dofs_per_cell >
                   cell->child(most_general_child)->get_fe().dofs_per_cell)
                 most_general_child = child;
-            }
-          if (different_elements == true)
-            fe_index = cell->child(most_general_child)->active_fe_index();
+          const unsigned int target_fe_index = cell->child(most_general_child)->active_fe_index();
+
+          const unsigned int dofs_per_cell=cell->get_dof_handler().get_fe()[target_fe_index].dofs_per_cell;
 
+          std::vector<Vector<typename VECTOR::value_type> >(in_size,
+                                                            Vector<typename VECTOR::value_type>(dofs_per_cell))
+          .swap(dof_values_on_cell[n_cf]);
+
+
+          // store the data of each of the input vectors. get this data
+          // as interpolated onto a finite element space that encompasses
+          // that of all the children. note that cell->get_interpolated_dof_values
+          // already does all of the interpolations between spaces
           for (unsigned int j=0; j<in_size; ++j)
-            {
-              // store the data of each of
-              // the input vectors
-              if (different_elements == false)
-                cell->get_interpolated_dof_values(all_in[j],
-                                                  dof_values_on_cell[n_cf][j]);
-              else
-                {
-                  // if we have different elements, first
-                  // interpolate the children's contribution to
-                  // the most general FE on the child
-                  // level. Then we manually write the
-                  // interpolation operation to the coarser
-                  // level
-                  dof_values_on_cell[n_cf][j].reinit (cell->child(most_general_child)->get_fe().dofs_per_cell);
-                  const unsigned int fe_ind_general =
-                    cell->child(most_general_child)->active_fe_index();
-                  for (unsigned int child=0; child<cell->n_children(); ++child)
-                    {
-                      tmp.reinit (cell->child(child)->get_fe().dofs_per_cell,
-                                  true);
-                      cell->child(child)->get_dof_values (all_in[j],
-                                                          tmp);
-                      const unsigned int child_ind =
-                        cell->child(child)->active_fe_index();
-                      if (child_ind != fe_ind_general)
-                        {
-                          tmp2.reinit (cell->child(most_general_child)->get_fe().dofs_per_cell,
-                                       true);
-
-                          // if the
-                          // matrix
-                          // has size
-                          // 0 and
-                          // the
-                          // corresponding
-                          // elements
-                          // have
-                          // more
-                          // than
-                          // zero
-                          // DoFs,
-                          // then
-                          // this
-                          // means
-                          // that the
-                          // internal::extract_interpolation_matrices
-                          // function
-                          // above
-                          // couldn't
-                          // get the
-                          // interpolation
-                          // matrix
-                          // between
-                          // this
-                          // pair of
-                          // elements. since
-                          // we need
-                          // it here,
-                          // this is
-                          // an error
-                          if ((interpolation_hp(fe_ind_general, child_ind).m() == 0)
-                              &&
-                              (interpolation_hp(fe_ind_general, child_ind).n() == 0)
-                              &&
-                              ((tmp2.size() > 0) || (tmp.size() > 0)))
-                            Assert (false,
-                                    ExcMessage (std::string("No interpolation from element ")
-                                                + cell->child(child)->get_fe().get_name()
-                                                + " to element "
-                                                + cell->child(most_general_child)->get_fe().get_name()
-                                                + " is available, but it is needed here."));
-                          interpolation_hp(fe_ind_general, child_ind).vmult (tmp2, tmp);
-                        }
-                      else
-                        tmp2.swap (tmp);
-
-                      // now do the interpolation
-                      // operation. FullMatrix only wants us
-                      // to call vmult if the matrix size is
-                      // actually non-zero, so check that
-                      // case
-                      const unsigned int dofs_per_cell = tmp2.size();
-                      tmp.reinit (dofs_per_cell, true);
-                      if (dofs_per_cell > 0)
-                        cell->child(most_general_child)->get_fe().
-                        get_restriction_matrix(child, cell->refinement_case()).vmult (tmp, tmp2);
-                      for (unsigned int i=0; i<dofs_per_cell; ++i)
-                        if (restriction_is_additive[fe_ind_general][i])
-                          dof_values_on_cell[n_cf][j](i) += tmp(i);
-                        else if (tmp(i) != zero_val)
-                          dof_values_on_cell[n_cf][j](i) = tmp(i);
-                    }
-                }
-            }
+            cell->get_interpolated_dof_values(all_in[j],
+                                              dof_values_on_cell[n_cf][j],
+                                              target_fe_index);
           cell_map[std::make_pair(cell->level(), cell->index())]
-            = Pointerstruct(&dof_values_on_cell[n_cf], fe_index);
+            = Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index);
           ++n_cf;
         }
     }

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.