]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Make it run, finally..
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 19 Aug 2013 03:15:49 +0000 (03:15 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 19 Aug 2013 03:15:49 +0000 (03:15 +0000)
git-svn-id: https://svn.dealii.org/trunk@30342 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/source/fe/fe_q_base.cc
deal.II/source/numerics/solution_transfer.cc
tests/hp/solution_transfer_03.cc

index 25f5a9f13493b9e851f301984aefea5a9ccceeba..0fa2931eb9a520f671e3c847af8a7e88f9130a09 100644 (file)
@@ -459,18 +459,17 @@ FE_Q_Base<POLY,dim,spacedim>::
 get_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
                           FullMatrix<double>       &interpolation_matrix) const
 {
-  Assert (interpolation_matrix.m() == this->dofs_per_cell,
-          ExcDimensionMismatch (interpolation_matrix.m(),
-                                this->dofs_per_cell));
-  Assert (interpolation_matrix.n() == x_source_fe.dofs_per_cell,
-          ExcDimensionMismatch (interpolation_matrix.m(),
-                                x_source_fe.dofs_per_cell));
-
   // go through the list of elements we can interpolate from
   if (const FE_Q_Base<POLY,dim,spacedim> *source_fe
       = dynamic_cast<const FE_Q_Base<POLY,dim,spacedim>*>(&x_source_fe))
     {
       // ok, source is a Q element, so we will be able to do the work
+      Assert (interpolation_matrix.m() == this->dofs_per_cell,
+              ExcDimensionMismatch (interpolation_matrix.m(),
+                                    this->dofs_per_cell));
+      Assert (interpolation_matrix.n() == x_source_fe.dofs_per_cell,
+              ExcDimensionMismatch (interpolation_matrix.m(),
+                                    x_source_fe.dofs_per_cell));
 
       // only evaluate Q dofs
       const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(this->degree+1);
@@ -530,6 +529,16 @@ get_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
       // element represents a function that is constant zero and has no
       // degrees of freedom, so the interpolation is simply a multiplication
       // with a n_dofs x 0 matrix. there is nothing to do here
+
+      // we would like to verify that the number of rows and columns of
+      // the matrix equals this->dofs_per_cell and zero. unfortunately,
+      // whenever we do FullMatrix::reinit(m,0), it sets both rows and
+      // columns to zero, instead of m and zero. thus, only test the
+      // number of columns
+      Assert (interpolation_matrix.n() == x_source_fe.dofs_per_cell,
+              ExcDimensionMismatch (interpolation_matrix.m(),
+                                    x_source_fe.dofs_per_cell));
+
     }
   else
     AssertThrow (false,
index fef10bdc175905103389a3faec2b8ce20b2f5864..9e65552c147d1a9f4bedcb3f4c5d6af871398c85 100644 (file)
@@ -571,24 +571,33 @@ interpolate (const std::vector<VECTOR> &all_in,
                          // back to size zero. so if we get here and that is
                          // the wrong size, then this may be because the elements
                          // haven't implemented the correct function yet
+                         //
+                         // there is one wrinkle. we would like to only error out if
+                         // the size of the matrix is 0 times 0 but at least one
+                         // of the elements has more than one dof per cell. the
+                         // problem is that if you reinit a matrix to 4x0, it automatically
+                         // sets its size to 0x0. so we can only execute the following
+                         // test if *both* elements have dofs_per_cell>0, not if *at
+                         // least one* have.
                          Assert (! ((interpolation_hp(new_fe_index,old_index).m() == 0)
                                     &&
                                     (interpolation_hp(new_fe_index,old_index).n() == 0)
                                     &&
                                     ((dof_handler->get_fe()[new_fe_index].dofs_per_cell > 0)
-                                     ||
+                                     &&
                                      (dof_handler->get_fe()[old_index].dofs_per_cell > 0))),
                                  ExcMessage ("The interpolation between two different "
                                              "elements you are trying to use here has "
                                              "not been implemented for this pair of "
-                                             "elememts!"));
+                                             "elements!"));
 
-                          AssertDimension (local_values.size(),
-                                           interpolation_hp(new_fe_index,old_index).m());
                           // simple case where all children have the
                           // same FE index: just interpolate to their FE
                           // first and then use the standard routines
-                          interpolation_hp(new_fe_index,old_index).vmult (local_values, tmp);
+                         if (tmp.size() > 0)
+                           interpolation_hp(new_fe_index,old_index).vmult (local_values, tmp);
+                         else
+                           local_values = 0;
                         }
 
                       if (cell->has_children() == false)
@@ -612,16 +621,17 @@ interpolate (const std::vector<VECTOR> &all_in,
                                               (interpolation_hp(c_index,old_index).n() == 0)
                                               &&
                                               ((dof_handler->get_fe()[c_index].dofs_per_cell > 0)
-                                               ||
+                                               &&
                                                (dof_handler->get_fe()[old_index].dofs_per_cell > 0))),
                                            ExcMessage ("The interpolation between two different "
                                                        "elements you are trying to use here has "
                                                        "not been implemented for this pair of "
-                                                       "elememts!"));
-                                   
-                                    AssertDimension (local_values.size(),
-                                                     interpolation_hp(c_index,old_index).m());
-                                    interpolation_hp(c_index,old_index).vmult (local_values, tmp);
+                                                       "elements!"));
+
+                                   if (tmp.size() > 0)
+                                     interpolation_hp(c_index,old_index).vmult (local_values, tmp);
+                                    else
+                                      local_values = 0;
                                   }
                                 else
                                   local_values = tmp;
index 82a08810ab1a53a2328d875615bc532f8aae5eca..c8d90755666e17ccb606b9340f7aff0989a166ad 100644 (file)
@@ -57,13 +57,11 @@ int main()
   triangulation.refine_global (1);
 
 
-  hp::DoFHandler<2> dof_handler(triangulation);
   hp::FECollection<2>     fe_collection;
-
-
   fe_collection.push_back(FE_Q<2>(1));
   fe_collection.push_back(FE_Nothing<2>());
 
+  hp::DoFHandler<2> dof_handler(triangulation);
 
   // Assign FE to cells
   hp::DoFHandler<2>::active_cell_iterator cell;
@@ -95,18 +93,18 @@ int main()
   data_out.write_vtu (deallog.get_file_stream());
 
 
+  // Interpoalte solution
+  SolutionTransfer<2, Vector<double>, hp::DoFHandler<2> > solultion_trans(dof_handler);
+  solultion_trans.prepare_for_coarsening_and_refinement(solution);
+
+  triangulation.execute_coarsening_and_refinement ();
   // Assign FE_Q_ to all cells
   cell = dof_handler.begin_active();
   for( ; cell != endc; ++cell){
     cell->set_active_fe_index(0);          
   }
-
-
-  // Interpoalte solution
-  SolutionTransfer<2, Vector<double>, hp::DoFHandler<2> > solultion_trans(dof_handler);
-  solultion_trans.prepare_for_coarsening_and_refinement(solution);
-  triangulation.execute_coarsening_and_refinement ();
   dof_handler.distribute_dofs (fe_collection);
+  
   Vector<double> new_solution(dof_handler.n_dofs());
   solultion_trans.interpolate(solution, new_solution);
 

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.