From: bangerth Date: Thu, 10 Aug 2006 20:41:57 +0000 (+0000) Subject: Implement VectorTools::interpolate for hp DoFHandlers X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2ef1e36a32d488cc68b90a0a51244c508bc86b6f;p=dealii-svn.git Implement VectorTools::interpolate for hp DoFHandlers git-svn-id: https://svn.dealii.org/trunk@13643 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index ca7561820a..b2876ed662 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -411,13 +411,17 @@ class VectorTools * make the result continuous * again. * + * The template argument DH + * may either be of type DoFHandler or + * hp::DoFHandler. + * * See the general documentation * of this class for further * information. */ - template + template static void interpolate (const Mapping &mapping, - const DoFHandler &dof, + const DH &dof, const Function &function, VECTOR &vec); @@ -426,11 +430,11 @@ class VectorTools * function above with * mapping=MappingQ1@@(). */ - template - static void interpolate (const DoFHandler &dof, + template + static void interpolate (const DH &dof, const Function &function, VECTOR &vec); - + /** * Interpolate different finite * element spaces. The diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index efb99a2b99..4d6334d0c5 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -44,37 +44,45 @@ //TODO[RH]: Use StaticQ1Mapping where appropriate -template +template void VectorTools::interpolate (const Mapping &mapping, - const DoFHandler &dof, + const DH &dof, const Function &function, VECTOR &vec) { Assert (dof.get_fe().n_components() == function.n_components, ExcComponentMismatch()); - const FiniteElement &fe = dof.get_fe(); - const unsigned int n_components = fe.n_components(); - const bool fe_is_system = (n_components != 1); + const hp::FECollection fe (dof.get_fe()); + const unsigned int n_components = fe.n_components(); + const bool fe_is_system = (n_components != 1); - typename DoFHandler::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); + typename DH::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); // For FESystems many of the - // unit_support_points will - // appear multiply, as a point - // may be unit_support_point - // for several of the components - // of the system. - // The following is rather - // complicated as it is - // avoided to evaluate - // the vectorfunction multiply at - // the same point on a cell. - const std::vector > & - unit_support_points = fe.get_unit_support_points(); - Assert (unit_support_points.size() != 0, - ExcNonInterpolatingFE()); + // unit_support_points will appear + // multiply, as a point may be + // unit_support_point for several of the + // components of the system. The following + // is rather complicated, but at least + // attempts to avoid evaluating the + // vectorfunction multiple times at the + // same point on a cell. + // + // note that we have to set up all of the + // following arrays for each of the + // elements in the FECollection (which + // means only once if this is for a regular + // DoFHandler) + std::vector > > unit_support_points (fe.size()); + for (unsigned int fe_index=0; fe_index &mapping, // the following vector collects all rep dofs. // the position of a rep dof within this vector // is called rep index. - std::vector dofs_of_rep_points; + std::vector > dofs_of_rep_points(fe.size()); // the following table converts a dof i // to the rep index. - std::vector dof_to_rep_index_table; - unsigned int n_rep_points=0; - for (unsigned int i=0; i > dof_to_rep_index_table(fe.size()); + + std::vector n_rep_points (fe.size(), 0); + + for (unsigned int fe_index=0; fe_index0; --j) - if (unit_support_points[i] - == unit_support_points[dofs_of_rep_points[j-1]]) - { - dof_to_rep_index_table.push_back(j-1); - representative=false; - break; - } - - if (representative) + for (unsigned int i=0; i0; --j) + if (unit_support_points[fe_index][i] + == unit_support_points[fe_index][dofs_of_rep_points[fe_index][j-1]]) + { + dof_to_rep_index_table[fe_index].push_back(j-1); + representative=false; + break; + } + + if (representative) + { + // rep_index=dofs_of_rep_points.size() + dof_to_rep_index_table[fe_index].push_back(dofs_of_rep_points[fe_index].size()); + // dofs_of_rep_points[rep_index]=i + dofs_of_rep_points[fe_index].push_back(i); + ++n_rep_points[fe_index]; + } } - } - Assert(dofs_of_rep_points.size()==n_rep_points, ExcInternalError()); - Assert(dof_to_rep_index_table.size()==fe.dofs_per_cell, ExcInternalError()); - - std::vector dofs_on_cell (fe.dofs_per_cell); - std::vector > rep_points (n_rep_points); + Assert(dofs_of_rep_points[fe_index].size()==n_rep_points[fe_index], + ExcInternalError()); + Assert(dof_to_rep_index_table[fe_index].size()==fe[fe_index].dofs_per_cell, + ExcInternalError()); + } + + const unsigned int max_rep_points = *std::max_element (n_rep_points.begin(), + n_rep_points.end()); + std::vector dofs_on_cell (fe.max_dofs_per_cell()); + std::vector > rep_points (max_rep_points); + // get space for the values of the // function at the rep support points. // // have two versions, one for system fe // and one for scalar ones, to take the // more efficient one respectively - std::vector function_values_scalar (n_rep_points); - std::vector > function_values_system (n_rep_points, - Vector(fe.n_components())); + std::vector > function_values_scalar(fe.size()); + std::vector > > function_values_system(fe.size()); // Make a quadrature rule from support points // to feed it into FEValues - Quadrature support_quadrature(unit_support_points); + hp::QCollection support_quadrature; + for (unsigned int fe_index=0; fe_index(unit_support_points[fe_index])); // Transformed support points are computed by // FEValues - FEValues fe_values (mapping, fe, support_quadrature, update_q_points); + hp::MappingCollection mapping_collection (mapping); + + hp::FEValues fe_values (mapping_collection, + fe, support_quadrature, update_q_points); for (; cell!=endc; ++cell) { + const unsigned int fe_index = cell->active_fe_index(); + // for each cell: // get location of finite element // support_points fe_values.reinit(cell); const std::vector >& support_points = - fe_values.get_quadrature_points(); + fe_values.get_present_fe_values().get_quadrature_points(); // pick out the representative // support points - for (unsigned int j=0; jget_dof_indices (dofs_on_cell); @@ -174,39 +200,45 @@ void VectorTools::interpolate (const Mapping &mapping, // get function values at // these points. Here: get // all components - function.vector_value_list (rep_points, function_values_system); + function_values_system[fe_index] + .resize (n_rep_points[fe_index], + Vector (fe[fe_index].n_components())); + function.vector_value_list (rep_points, + function_values_system[fe_index]); // distribute the function // values to the global // vector - for (unsigned int i=0; i -void VectorTools::interpolate (const DoFHandler &dof, +template +void VectorTools::interpolate (const DH &dof, const Function &function, VECTOR &vec) { diff --git a/deal.II/deal.II/source/numerics/vectors.instance.h b/deal.II/deal.II/source/numerics/vectors.instance.h index e83be170f4..538dd5f66c 100644 --- a/deal.II/deal.II/source/numerics/vectors.instance.h +++ b/deal.II/deal.II/source/numerics/vectors.instance.h @@ -27,6 +27,18 @@ void VectorTools::interpolate const Function&, VEC&); +template +void VectorTools::interpolate +(const Mapping&, + const hp::DoFHandler&, + const Function&, + VEC&); +template +void VectorTools::interpolate +(const hp::DoFHandler&, + const Function&, + VEC&); + // Should these be instantiated for every combination of two types? template void VectorTools::interpolate