From: bangerth Date: Mon, 19 Aug 2013 03:15:49 +0000 (+0000) Subject: Make it run, finally.. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5d907a85c0e5f9d5dc02f697e0ab93fc3f6817d6;p=dealii-svn.git Make it run, finally.. git-svn-id: https://svn.dealii.org/trunk@30342 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/source/fe/fe_q_base.cc b/deal.II/source/fe/fe_q_base.cc index 25f5a9f134..0fa2931eb9 100644 --- a/deal.II/source/fe/fe_q_base.cc +++ b/deal.II/source/fe/fe_q_base.cc @@ -459,18 +459,17 @@ FE_Q_Base:: get_interpolation_matrix (const FiniteElement &x_source_fe, FullMatrix &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 *source_fe = dynamic_cast*>(&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(this->degree+1); @@ -530,6 +529,16 @@ get_interpolation_matrix (const FiniteElement &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, diff --git a/deal.II/source/numerics/solution_transfer.cc b/deal.II/source/numerics/solution_transfer.cc index fef10bdc17..9e65552c14 100644 --- a/deal.II/source/numerics/solution_transfer.cc +++ b/deal.II/source/numerics/solution_transfer.cc @@ -571,24 +571,33 @@ interpolate (const std::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 &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; diff --git a/tests/hp/solution_transfer_03.cc b/tests/hp/solution_transfer_03.cc index 82a08810ab..c8d9075566 100644 --- a/tests/hp/solution_transfer_03.cc +++ b/tests/hp/solution_transfer_03.cc @@ -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, 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, 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 new_solution(dof_handler.n_dofs()); solultion_trans.interpolate(solution, new_solution);