From: Luca Heltai Date: Sat, 11 Apr 2015 14:48:07 +0000 (+0200) Subject: Added get_position_vector to VectorTools X-Git-Tag: v8.3.0-rc1~267^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=47f5efdcccbf6977ac4c43a593f0836e5bc125f9;p=dealii.git Added get_position_vector to VectorTools Removed references to update_euler_vector from tets. --- diff --git a/doc/news/changes.h b/doc/news/changes.h index 2c3d4c119b..d3ef58962d 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -459,6 +459,20 @@ inconvenience this causes. (Luca Heltai, 2015/04/12) +
  • New: A new VectorTools::get_position_vector() function has been + added to the library that allows one to interpolate the Geometry of + a (possibly curved) triangulation to vector finite element fields + of at least spacedim components. +
    + (Luca Heltai, 2015/04/11) +
  • + +
  • New: A new MappingFEField() class has been added to the library + that generalizes MappingQEulerian to allow arbitrary FiniteElements. +
    + (Marco Tezzele, Luca Heltai, 2015/04/06) +
  • +
  • Changed: The cells of coarse meshes in parallel::distributed::Triangulation used to be ordered in a Cuthill-McKee numbering, which yields very high surface-to-volume ratios of the individual @@ -511,12 +525,6 @@ inconvenience this causes. (Wolfgang Bangerth, 2015/04/08)
    • -
  • New: A new MappingFEField() class has been added to the library - that generalizes MappingQEulerian to allow arbitrary FiniteElements. -
    - (Marco Tezzele, Luca Heltai, 2015/04/06) -
  • -
  • Changed: In the spirit of the changes made to the distinction between Point and Tensor objects discussed above, the first argument to GridTools::shift() has been changed from a Point to a Tensor@<1,dim@>. diff --git a/include/deal.II/fe/mapping_fe_field.h b/include/deal.II/fe/mapping_fe_field.h index 7e065af27c..6ecb446bb5 100644 --- a/include/deal.II/fe/mapping_fe_field.h +++ b/include/deal.II/fe/mapping_fe_field.h @@ -51,7 +51,7 @@ DEAL_II_NAMESPACE_OPEN * components of the FiniteElement to use for the mapping. * * Typically, the DoFHandler operates on a finite element that is - * constructed as a system element (FESystem) from continuous FE_Q() + * constructed as a system element (FESystem()) from continuous FE_Q() * (for iso-parametric discretizations) or FE_Bernstein() (for * iso-geometric discretizations) objects. An example is shown below: * @@ -61,9 +61,10 @@ DEAL_II_NAMESPACE_OPEN * DoFHandler dhq(triangulation); * dhq.distribute_dofs(fesystem); * Vector eulerq(dhq.n_dofs()); + * // Fills the euler vector with information from the Triangulation + * VectorTools::get_position_vector(dhq, eulerq); * const ComponentMask mask(spacedim, true); * MappingFEField map(eulerq, dhq, mask); - * map.update_euler_vector_using_triangulation(eulerq); * @endcode * * @author Luca Heltai, Marco Tezzele 2013, 2015 @@ -89,16 +90,10 @@ public: * euler_vector actually represents a valid geometry (i.e., one with * no inverted cells, or with no zero-volume cells). * - * In order to facilitate the user, we provide an - * update_euler_vector_using_triangulation() method to initialize a - * euler_vector in a way similar to what happens in - * MappingQEulerian, by creating a euler_vector which interpolates - * the underlying Triangulation. - * * If the underlying FiniteElement is a system of FE_Q(), and * euler_vector is initialized using - * update_euler_vector_using_triangulation(), then this class is in - * all respects identical to MappingQ(). + * VectorTools::get_position_vector(), then this class is in all + * respects identical to MappingQ(). * * The optional ComponentMask argument can be used to specify what * components of the FiniteElement to use for the geometrical @@ -126,24 +121,6 @@ public: */ MappingFEField (const MappingFEField &mapping); - /** - * Helper function to fill the given euler vector with information - * coming from the triangulation, by interpolating the geometry of - * the Triangulation associated with the DoFHandler used at - * construction time. - * - * The resulting map is guaranteed to be interpolatory at the - * vertices of the Triangulation. Notice that this may or may not be - * meaningful, depending on the FiniteElement you have used to - * construct this MappingFEField. - * - * If the underlying FiniteElement is a system of FE_Q(), and - * euler_vector is initialized using this function, then this class - * is in all respects identical to MappingQ(). - */ - void update_euler_vector_using_triangulation(VECTOR &vector) const; - - /** * Transforms the point @p p on the unit cell to the point @p p_real on the diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index b4c5c2b041..ee6c85d6f9 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -2188,7 +2188,44 @@ namespace VectorTools const InVector &v, const unsigned int component); //@} + /** + * Geometrical interpolation + */ + //@{ + /** + * Given a DoFHandler containing at least a spacedim vector field, + * this function interpolates the Triangulation at the support + * points of a FE_Q() finite element of the same degree as + * dh->get_fe().degree. + * + * Curved manifold are respected, and the resulting VECTOR will be + * geometrically consistent. + * + * The resulting map is guaranteed to be interpolatory at the + * support points of a FE_Q() finite element of the same degree as + * dh->get_fe().degree. Notice that this may or may not be + * meaningful, depending on the FiniteElement you have distribed in + * dh. + * + * If the underlying finite element is an FE_Q(1)^spacedim, then the + * resulting VECTOR is a finite element field representation of the + * vertices of the Triangulation. + * + * The optional ComponentMask argument can be used to specify what + * components of the FiniteElement to use to describe the + * geometry. If no mask is specified at construction time, then a + * default one is used, i.e., the first spacedim components of the + * FiniteElement are assumed to represent the geometry of the + * problem. + * + * @author Luca Heltai, 2015 + */ + template + void get_position_vector(const DH &dh, + VECTOR &vector, + const ComponentMask mask=ComponentMask()); + //@} /** * Exception diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 0a41091f16..d156bb0ac3 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -6880,6 +6880,98 @@ namespace VectorTools { return compute_mean_value(StaticMappingQ1::mapping, dof, quadrature, v, component); } + + + template + void get_position_vector(const DH &dh, + VECTOR &vector, + const ComponentMask mask) + { + AssertDimension(vector.size(), dh.n_dofs()); + + const unsigned int dim=DH::dimension; + const unsigned int spacedim=DH::space_dimension; + const FiniteElement &fe = dh.get_fe(); + + + // Construct default fe_mask; + const ComponentMask fe_mask(mask.size() ? mask : + ComponentMask(fe.get_nonzero_components(0).size(), true)); + + AssertDimension(fe_mask.size(), fe.get_nonzero_components(0).size()); + + std::vector fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int); + unsigned int size = 0; + for (unsigned int i=0; i quad(fe.get_unit_support_points()); + + MappingQ map_q(fe.degree); + FEValues fe_v(map_q, fe, quad, update_quadrature_points); + std::vector dofs(fe.dofs_per_cell); + + AssertDimension(fe.dofs_per_cell, fe.get_unit_support_points().size()); + Assert(fe.is_primitive(), ExcMessage("FE is not Primitive! This won't work.")); + + for (cell = dh.begin_active(); cell != dh.end(); ++cell) + { + fe_v.reinit(cell); + cell->get_dof_indices(dofs); + const std::vector > &points = fe_v.get_quadrature_points(); + for (unsigned int q = 0; q < points.size(); ++q) + { + unsigned int comp = fe.system_to_component_index(q).first; + if (fe_mask[comp]) + vector(dofs[q]) = points[q][fe_to_real[comp]]; + } + } + } + else + { + // Construct a FiniteElement with FE_Q^spacedim, and call this + // function again. + // + // Once we have this, interpolate with the given finite element + // to get a Mapping which is interpolatory at the support points + // of FE_Q(fe.degree()) + FESystem feq(FE_Q(fe.degree), spacedim); + DH dhq(dh.get_tria()); + dhq.distribute_dofs(feq); + Vector eulerq(dhq.n_dofs()); + const ComponentMask maskq(spacedim, true); + get_position_vector(dhq, eulerq); + + FullMatrix transfer(fe.dofs_per_cell, feq.dofs_per_cell); + const std::vector > &points = feq.get_unit_support_points(); + + // Here construct the interpolation matrix from FE_Q^spacedim to + // the FiniteElement used by euler_dof_handler. + // + // The interpolation matrix is then passed to the + // VectorTools::interpolate() function to generate + for (unsigned int i=0; i::update_internal_dofs ( local_dofs[i] = (*euler_vector)(dof_indices[i]); } - -template -void MappingFEField::update_euler_vector_using_triangulation -(VECTOR &vector) const -{ - AssertDimension(vector.size(), euler_dof_handler->n_dofs()); - if ( fe->has_support_points() ) - { - typename DH::active_cell_iterator cell; - const Quadrature quad(fe->get_unit_support_points()); - - MappingQ map_q(fe->degree); - FEValues fe_v(map_q, *fe, quad, update_quadrature_points); - std::vector dofs(fe->dofs_per_cell); - - AssertDimension(fe->dofs_per_cell, fe->get_unit_support_points().size()); - Assert(fe->is_primitive(), ExcMessage("FE is not Primitive! This won't work.")); - - for (cell = euler_dof_handler->begin_active(); cell != euler_dof_handler->end(); ++cell) - { - fe_v.reinit(cell); - cell->get_dof_indices(dofs); - const std::vector > &points = fe_v.get_quadrature_points(); - for (unsigned int q = 0; q < points.size(); ++q) - { - unsigned int comp = fe->system_to_component_index(q).first; - if (fe_mask[comp]) - vector(dofs[q]) = points[q][fe_to_real[comp]]; - } - } - } - else - { - // Construct a MappingFEField with an FEQ, and construct a - // standard iso-parametric interpolation on the Triangulation. - // - // Once we have this, interpolate with the given finite element - // to get a Mapping which is interpolatory at the support points - // of FE_Q(fe->degree()) - FESystem feq(FE_Q(fe->degree), spacedim); - DH dhq(euler_dof_handler->get_tria()); - dhq.distribute_dofs(feq); - VECTOR eulerq(dhq.n_dofs()); - const ComponentMask maskq(spacedim, true); - MappingFEField newfe(eulerq, dhq, maskq); - - newfe.update_euler_vector_using_triangulation(eulerq); - - FullMatrix transfer(fe->dofs_per_cell, feq.dofs_per_cell); - const std::vector > &points = feq.get_unit_support_points(); - - // Here construct the interpolation matrix from FE_Q^spacedim to - // the FiniteElement used by euler_dof_handler. - // - // The interpolation matrix is then passed to the - // VectorTools::interpolate() function to generate - for (unsigned int i=0; idofs_per_cell; ++i) - { - unsigned int comp_i = fe->system_to_component_index(i).first; - if (fe_mask[comp_i]) - for (unsigned int j=0; jshape_value(i, points[j]); - } - } - VectorTools::interpolate(dhq, *euler_dof_handler, transfer, eulerq, vector); - } -} - - - - // explicit instantiations #include "mapping_fe_field.inst" diff --git a/source/numerics/vector_tools_interpolate.inst.in b/source/numerics/vector_tools_interpolate.inst.in index 5430ee6abc..3c3bad4c85 100644 --- a/source/numerics/vector_tools_interpolate.inst.in +++ b/source/numerics/vector_tools_interpolate.inst.in @@ -40,6 +40,12 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const VEC&, VEC&); + template + void get_position_vector + (const DoFHandler&, + VEC&, + const ComponentMask); + \} #endif } diff --git a/tests/fe/mapping_fe_field_real_to_unit_b1.cc b/tests/fe/mapping_fe_field_real_to_unit_b1.cc index 0463f9b182..6ac4bd9a47 100644 --- a/tests/fe/mapping_fe_field_real_to_unit_b1.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_b1.cc @@ -31,6 +31,7 @@ #include #include #include +#include using namespace dealii; @@ -81,9 +82,9 @@ void test_real_to_unit_cell() Vector eulerq(dhb.n_dofs()); const ComponentMask mask(spacedim, true); + VectorTools::get_position_vector(dhb, eulerq, mask); MappingFEField map(eulerq, dhb, mask); - map.update_euler_vector_using_triangulation(eulerq); typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); diff --git a/tests/fe/mapping_fe_field_real_to_unit_b2_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_b2_curved.cc index 7efb2357d9..a11b2d7989 100644 --- a/tests/fe/mapping_fe_field_real_to_unit_b2_curved.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_b2_curved.cc @@ -37,6 +37,7 @@ #include #include #include +#include using namespace dealii; @@ -99,9 +100,8 @@ void test_real_to_unit_cell() Vector eulerq(dhb.n_dofs()); const ComponentMask mask(spacedim, true); + VectorTools::get_position_vector(dhb, eulerq, mask); MappingFEField map(eulerq, dhb, mask); - - map.update_euler_vector_using_triangulation(eulerq); typename Triangulation::active_cell_iterator diff --git a/tests/fe/mapping_fe_field_real_to_unit_b2_mask.cc b/tests/fe/mapping_fe_field_real_to_unit_b2_mask.cc index f8e4be9d7c..e8edd1d89d 100644 --- a/tests/fe/mapping_fe_field_real_to_unit_b2_mask.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_b2_mask.cc @@ -33,6 +33,7 @@ #include #include #include +#include using namespace dealii; @@ -98,11 +99,9 @@ void test_real_to_unit_cell() ComponentMask mask(spacedim+1, true); mask.set(0, false); + VectorTools::get_position_vector(dhb, eulerq, mask); MappingFEField map(eulerq, dhb, mask); - map.update_euler_vector_using_triangulation(eulerq); - - typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); diff --git a/tests/fe/mapping_fe_field_real_to_unit_b3_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_b3_curved.cc index 1ed55433fe..0df01fc4a3 100644 --- a/tests/fe/mapping_fe_field_real_to_unit_b3_curved.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_b3_curved.cc @@ -37,6 +37,7 @@ #include #include #include +#include using namespace dealii; @@ -99,9 +100,8 @@ void test_real_to_unit_cell() Vector eulerq(dhb.n_dofs()); const ComponentMask mask(spacedim, true); + VectorTools::get_position_vector(dhb, eulerq, mask); MappingFEField map(eulerq, dhb, mask); - - map.update_euler_vector_using_triangulation(eulerq); typename Triangulation::active_cell_iterator diff --git a/tests/fe/mapping_fe_field_real_to_unit_b4_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_b4_curved.cc index 2cdf1cebae..990ea4d28e 100644 --- a/tests/fe/mapping_fe_field_real_to_unit_b4_curved.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_b4_curved.cc @@ -37,6 +37,7 @@ #include #include #include +#include using namespace dealii; @@ -99,11 +100,9 @@ void test_real_to_unit_cell() Vector eulerq(dhb.n_dofs()); const ComponentMask mask(spacedim, true); + VectorTools::get_position_vector(dhb, eulerq, mask); MappingFEField map(eulerq, dhb, mask); - map.update_euler_vector_using_triangulation(eulerq); - - typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); diff --git a/tests/fe/mapping_fe_field_real_to_unit_b5_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_b5_curved.cc index 7886f74919..e6843b985f 100644 --- a/tests/fe/mapping_fe_field_real_to_unit_b5_curved.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_b5_curved.cc @@ -37,6 +37,7 @@ #include #include #include +#include using namespace dealii; @@ -99,11 +100,9 @@ void test_real_to_unit_cell() Vector eulerq(dhb.n_dofs()); const ComponentMask mask(spacedim, true); + VectorTools::get_position_vector(dhb, eulerq, mask); MappingFEField map(eulerq, dhb, mask); - map.update_euler_vector_using_triangulation(eulerq); - - typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); diff --git a/tests/fe/mapping_fe_field_real_to_unit_q1.cc b/tests/fe/mapping_fe_field_real_to_unit_q1.cc index 9ac09034dc..debf15ae94 100644 --- a/tests/fe/mapping_fe_field_real_to_unit_q1.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_q1.cc @@ -30,6 +30,7 @@ #include #include #include +#include using namespace dealii; @@ -80,10 +81,9 @@ void test_real_to_unit_cell() Vector eulerq(dhq.n_dofs()); const ComponentMask mask(spacedim, true); + VectorTools::get_position_vector(dhq, eulerq, mask); MappingFEField map(eulerq, dhq, mask); - map.update_euler_vector_using_triangulation(eulerq); - typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); diff --git a/tests/fe/mapping_fe_field_real_to_unit_q2_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_q2_curved.cc index 0cce6f055b..adfe6f3ed1 100644 --- a/tests/fe/mapping_fe_field_real_to_unit_q2_curved.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_q2_curved.cc @@ -36,6 +36,7 @@ #include #include #include +#include #include #include @@ -101,11 +102,9 @@ void test_real_to_unit_cell() Vector eulerq(dhq.n_dofs()); const ComponentMask mask(spacedim, true); + VectorTools::get_position_vector(dhq, eulerq, mask); MappingFEField map(eulerq, dhq, mask); - map.update_euler_vector_using_triangulation(eulerq); - - typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); diff --git a/tests/fe/mapping_fe_field_real_to_unit_q2_mask.cc b/tests/fe/mapping_fe_field_real_to_unit_q2_mask.cc index bd8ef97191..074c7aa957 100644 --- a/tests/fe/mapping_fe_field_real_to_unit_q2_mask.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_q2_mask.cc @@ -30,6 +30,7 @@ #include #include #include +#include using namespace dealii; @@ -83,10 +84,9 @@ void test_real_to_unit_cell() ComponentMask mask(spacedim+1, true); mask.set(0, false); + VectorTools::get_position_vector(dhq, eulerq, mask); MappingFEField map(eulerq, dhq, mask); - map.update_euler_vector_using_triangulation(eulerq); - typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); diff --git a/tests/fe/mapping_fe_field_real_to_unit_q3_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_q3_curved.cc index e70517993f..54c5afb21f 100644 --- a/tests/fe/mapping_fe_field_real_to_unit_q3_curved.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_q3_curved.cc @@ -36,6 +36,7 @@ #include #include #include +#include #include #include @@ -101,11 +102,9 @@ void test_real_to_unit_cell() Vector eulerq(dhq.n_dofs()); const ComponentMask mask(spacedim, true); + VectorTools::get_position_vector(dhq, eulerq, mask); MappingFEField map(eulerq, dhq, mask); - map.update_euler_vector_using_triangulation(eulerq); - - typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); diff --git a/tests/fe/mapping_fe_field_real_to_unit_q4_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_q4_curved.cc index 3ff565dfe3..2690693320 100644 --- a/tests/fe/mapping_fe_field_real_to_unit_q4_curved.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_q4_curved.cc @@ -36,6 +36,7 @@ #include #include #include +#include using namespace dealii; @@ -98,11 +99,9 @@ void test_real_to_unit_cell() Vector eulerq(dhq.n_dofs()); const ComponentMask mask(spacedim, true); + VectorTools::get_position_vector(dhq, eulerq, mask); MappingFEField map(eulerq, dhq, mask); - map.update_euler_vector_using_triangulation(eulerq); - - typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); diff --git a/tests/fe/mapping_fe_field_real_to_unit_q5_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_q5_curved.cc index 1e04f9527b..669ea368a0 100644 --- a/tests/fe/mapping_fe_field_real_to_unit_q5_curved.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_q5_curved.cc @@ -36,6 +36,7 @@ #include #include #include +#include using namespace dealii; @@ -98,11 +99,9 @@ void test_real_to_unit_cell() Vector eulerq(dhq.n_dofs()); const ComponentMask mask(spacedim, true); + VectorTools::get_position_vector(dhq, eulerq, mask); MappingFEField map(eulerq, dhq, mask); - map.update_euler_vector_using_triangulation(eulerq); - - typename Triangulation::active_cell_iterator cell = triangulation.begin_active();