From: hartmann Date: Thu, 13 Jun 2002 12:49:10 +0000 (+0000) Subject: New interpolate, back_interpolate and interpolation_difference functions with Constra... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=38e6977d6f6cd694bb9942085fc7cf18c90a1d60;p=dealii-svn.git New interpolate, back_interpolate and interpolation_difference functions with ConstraintMatrix arguments. git-svn-id: https://svn.dealii.org/trunk@6092 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_tools.h b/deal.II/deal.II/include/fe/fe_tools.h index f01d0bb95e..94a6a8c14e 100644 --- a/deal.II/deal.II/include/fe/fe_tools.h +++ b/deal.II/deal.II/include/fe/fe_tools.h @@ -20,6 +20,8 @@ template class FiniteElement; template class DoFHandler; template class Vector; template class FE_Q; +class ConstraintMatrix; + #include #include @@ -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 static void interpolate (const DoFHandler &dof1, @@ -157,6 +169,47 @@ class FETools const DoFHandler &dof2, Vector &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 + static void interpolate (const DoFHandler &dof1, + const Vector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints, + Vector &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 static void back_interpolate (const DoFHandler &dof1, const Vector &u1, const FiniteElement &fe2, + Vector &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 + static void back_interpolate (const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const Vector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, Vector &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 static void interpolation_difference(const DoFHandler &dof1, @@ -201,6 +299,31 @@ class FETools const FiniteElement &fe2, Vector &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 + static void interpolation_difference(const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const Vector &z1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + Vector &z1_difference); + /** * Gives the patchwise * extrapolation of a @p{dof1} @@ -306,6 +429,11 @@ class FETools * Exception */ DeclException0 (ExcTriangulationMismatch); + + /** + * Exception + */ + DeclException0 (ExcHangingNodesNotAllowed); /** * Exception diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index e02a3db509..28d615f943 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -22,6 +22,7 @@ #include #include #include +#include #include #include @@ -122,11 +123,21 @@ void FETools::get_interpolation_difference_matrix(const FiniteElement &fe1, } +template +void FETools::interpolate(const DoFHandler &dof1, + const Vector &u1, + const DoFHandler &dof2, + Vector &u2) +{ + interpolate(dof1, u1, dof2, ConstraintMatrix(), u2); +} + template void FETools::interpolate(const DoFHandler &dof1, const Vector &u1, const DoFHandler &dof2, + const ConstraintMatrix &constraints, Vector &u2) { Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch()); @@ -135,6 +146,16 @@ void FETools::interpolate(const DoFHandler &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 &dof1, std::vector index_multiplicity(dof2.n_dofs(),0); std::vector dofs (dofs_per_cell2); u2.clear (); - + for (; cell1!=endc1, cell2!=endc2; ++cell1, ++cell2) { + if (hanging_nodes_not_allowed) + for (unsigned int face=0; face::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 &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 @@ -185,30 +218,83 @@ void FETools::back_interpolate(const DoFHandler &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 u1_local(dofs_per_cell1); Vector u1_int_local(dofs_per_cell1); + typename DoFHandler::active_cell_iterator cell = dof1.begin_active(), + endc = dof1.end(); + FullMatrix interpolation_matrix(dofs_per_cell1, dofs_per_cell1); FETools::get_back_interpolation_matrix(dof1.get_fe(), fe2, interpolation_matrix); - - typename DoFHandler::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::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 +void FETools::back_interpolate(const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const Vector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + Vector &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 u2(dof2.n_dofs()); + interpolate(dof1, u1, dof2, constraints2, u2); + interpolate(dof2, u2, dof1, constraints1, u1_interpolated); + } +} + + template void FETools::interpolation_difference(const DoFHandler &dof1, const Vector &u1, @@ -219,30 +305,69 @@ void FETools::interpolation_difference(const DoFHandler &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 u1_local(dofs_per_cell1); - Vector 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 u1_local(dofs_per_cell); + Vector u1_diff_local(dofs_per_cell); - FullMatrix difference_matrix(dofs_per_cell1, dofs_per_cell1); + FullMatrix difference_matrix(dofs_per_cell, dofs_per_cell); FETools::get_interpolation_difference_matrix(dof1.get_fe(), fe2, difference_matrix); - typename DoFHandler::active_cell_iterator cell1 = dof1.begin_active(), - endc1 = dof1.end(); + typename DoFHandler::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::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 +void FETools::interpolation_difference(const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const Vector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + Vector &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 void FETools::extrapolate(const DoFHandler &dof1, const Vector &u1, @@ -732,18 +857,41 @@ void FETools::interpolate const DoFHandler &, Vector &); template +void FETools::interpolate +(const DoFHandler &, + const Vector &, + const DoFHandler &, + const ConstraintMatrix &, + Vector &); +template void FETools::back_interpolate (const DoFHandler &, const Vector &, const FiniteElement &, Vector &); template +void FETools::back_interpolate +(const DoFHandler &, + const ConstraintMatrix &, + const Vector &, + const DoFHandler &, + const ConstraintMatrix &, + Vector &); +template void FETools::interpolation_difference (const DoFHandler &, const Vector &, const FiniteElement &, Vector &); template +void FETools::interpolation_difference +(const DoFHandler &, + const ConstraintMatrix &, + const Vector &, + const DoFHandler &, + const ConstraintMatrix &, + Vector &); +template void FETools::extrapolate (const DoFHandler &, const Vector &, diff --git a/deal.II/doc/news/2002/c-3-4.html b/deal.II/doc/news/2002/c-3-4.html index b14d5b9d51..50c6aa2f5a 100644 --- a/deal.II/doc/news/2002/c-3-4.html +++ b/deal.II/doc/news/2002/c-3-4.html @@ -105,6 +105,22 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

deal.II

    +
  1. New: The existing FETools::interpolate, FETools::back_interpolate and FETools::interpolate_difference 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 ConstraintMatrix 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. +
    + (RH 2002/06/13) +

    +
  2. Change: The constructor of DoFHandler now takes a reference to a const Triangulation.