From 672743006c9b0c6846ab280efcef2d026fe78c81 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 22 Jul 2019 18:30:47 +0200 Subject: [PATCH] Enable constructor of MappingFEField on MG levels --- include/deal.II/fe/mapping_fe_field.h | 36 ++++++- source/fe/mapping_fe_field.cc | 136 ++++++++++++++++++++++---- 2 files changed, 153 insertions(+), 19 deletions(-) diff --git a/include/deal.II/fe/mapping_fe_field.h b/include/deal.II/fe/mapping_fe_field.h index ef7cd3c6bc..ef93982531 100644 --- a/include/deal.II/fe/mapping_fe_field.h +++ b/include/deal.II/fe/mapping_fe_field.h @@ -19,6 +19,7 @@ #include +#include #include #include @@ -118,6 +119,31 @@ public: const VectorType & euler_vector, const ComponentMask & mask = ComponentMask()); + /** + * Constructor taking vectors on the multigrid levels rather than the active + * cells only. The vector of vectors is expected to have as many entries as + * there are levels in the triangulation and provide valid data on each + * level, i.e., be of compatible length DoFHandler::n_dofs(level). Apart + * from the source of the vectors, the same arguments as in the other + * constructor need to be provided. + */ + MappingFEField(const DoFHandlerType & euler_dof_handler, + const std::vector &euler_vector, + const ComponentMask & mask = ComponentMask()); + + /** + * Constructor taking vectors on the multigrid levels rather than the active + * cells only, variant with MGLevelObject instead of std::vector. The vector + * of vectors is expected to have as many entries as there are levels in the + * triangulation and provide valid data on each level, i.e., be of + * compatible length DoFHandler::n_dofs(level). Apart from the source of the + * vectors, the same arguments as in the other constructor need to be + * provided. + */ + MappingFEField(const DoFHandlerType & euler_dof_handler, + const MGLevelObject &euler_vector, + const ComponentMask & mask = ComponentMask()); + /** * Copy constructor. */ @@ -504,12 +530,18 @@ private: * @} */ + /** + * Specifies whether we access unknowns on the active dofs (with a single + * Euler vector) or on the level dofs (via a vector of Euler vectors). + */ + const bool uses_level_dofs; /** * Reference to the vector of shifts. */ - SmartPointer> + std::vector< + SmartPointer>> euler_vector; /** diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc index fa22390109..5e0420a0c8 100644 --- a/source/fe/mapping_fe_field.cc +++ b/source/fe/mapping_fe_field.cc @@ -240,7 +240,8 @@ MappingFEField::MappingFEField( const DoFHandlerType &euler_dof_handler, const VectorType & euler_vector, const ComponentMask & mask) - : euler_vector(&euler_vector) + : uses_level_dofs(false) + , euler_vector({&euler_vector}) , euler_dof_handler(&euler_dof_handler) , fe_mask(mask.size() ? mask : @@ -262,10 +263,92 @@ MappingFEField::MappingFEField( } + +template +MappingFEField::MappingFEField( + const DoFHandlerType & euler_dof_handler, + const std::vector &euler_vector, + const ComponentMask & mask) + : uses_level_dofs(true) + , euler_dof_handler(&euler_dof_handler) + , fe_mask(mask.size() ? + mask : + ComponentMask( + euler_dof_handler.get_fe().get_nonzero_components(0).size(), + true)) + , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int) + , fe_values(this->euler_dof_handler->get_fe(), + get_vertex_quadrature(), + update_values) +{ + unsigned int size = 0; + for (unsigned int i = 0; i < fe_mask.size(); ++i) + { + if (fe_mask[i]) + fe_to_real[i] = size++; + } + AssertDimension(size, spacedim); + + Assert(euler_dof_handler.has_level_dofs(), + ExcMessage("The underlying did not call distribute_mg_dofs(). " + "In this case, the construction via level vectors does not " + "make sense.")); + AssertDimension(euler_vector.size(), + euler_dof_handler.get_triangulation().n_global_levels()); + this->euler_vector.clear(); + this->euler_vector.resize(euler_vector.size()); + for (unsigned int i = 0; i < euler_vector.size(); ++i) + this->euler_vector[i] = &euler_vector[i]; +} + + + +template +MappingFEField::MappingFEField( + const DoFHandlerType & euler_dof_handler, + const MGLevelObject &euler_vector, + const ComponentMask & mask) + : uses_level_dofs(true) + , euler_dof_handler(&euler_dof_handler) + , fe_mask(mask.size() ? + mask : + ComponentMask( + euler_dof_handler.get_fe().get_nonzero_components(0).size(), + true)) + , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int) + , fe_values(this->euler_dof_handler->get_fe(), + get_vertex_quadrature(), + update_values) +{ + unsigned int size = 0; + for (unsigned int i = 0; i < fe_mask.size(); ++i) + { + if (fe_mask[i]) + fe_to_real[i] = size++; + } + AssertDimension(size, spacedim); + + Assert(euler_dof_handler.has_level_dofs(), + ExcMessage("The underlying did not call distribute_mg_dofs(). " + "In this case, the construction via level vectors does not " + "make sense.")); + AssertDimension(euler_vector.max_level() + 1 - euler_vector.min_level(), + euler_dof_handler.get_triangulation().n_global_levels()); + this->euler_vector.clear(); + this->euler_vector.resize( + euler_dof_handler.get_triangulation().n_global_levels()); + for (unsigned int i = euler_vector.min_level(); i <= euler_vector.max_level(); + ++i) + this->euler_vector[i - euler_vector.min_level()] = &euler_vector[i]; +} + + + template MappingFEField::MappingFEField( const MappingFEField &mapping) - : euler_vector(mapping.euler_vector) + : uses_level_dofs(mapping.uses_level_dofs) + , euler_vector(mapping.euler_vector) , euler_dof_handler(mapping.euler_dof_handler) , fe_mask(mapping.fe_mask) , fe_to_real(mapping.fe_to_real) @@ -311,28 +394,42 @@ MappingFEField::get_vertices( const typename DoFHandler::cell_iterator dof_cell( *cell, euler_dof_handler); - Assert(dof_cell->active() == true, ExcInactiveCell()); + Assert(uses_level_dofs || dof_cell->active() == true, ExcInactiveCell()); AssertDimension(GeometryInfo::vertices_per_cell, fe_values.n_quadrature_points); AssertDimension(fe_to_real.size(), euler_dof_handler->get_fe().n_components()); - std::vector> values( - fe_values.n_quadrature_points, - Vector( - euler_dof_handler->get_fe().n_components())); - { std::lock_guard lock(fe_values_mutex); fe_values.reinit(dof_cell); - fe_values.get_function_values(*euler_vector, values); } + const unsigned int dofs_per_cell = euler_dof_handler->get_fe().dofs_per_cell; + std::vector dof_indices(dofs_per_cell); + if (uses_level_dofs) + dof_cell->get_mg_dof_indices(dof_indices); + else + dof_cell->get_dof_indices(dof_indices); + + const VectorType &vector = + uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0]; std::array, GeometryInfo::vertices_per_cell> vertices; - - for (unsigned int i = 0; i < GeometryInfo::vertices_per_cell; ++i) - for (unsigned int j = 0; j < fe_to_real.size(); ++j) - if (fe_to_real[j] != numbers::invalid_unsigned_int) - vertices[i][fe_to_real[j]] = values[i][j]; + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const unsigned int comp = fe_to_real + [euler_dof_handler->get_fe().system_to_component_index(i).first]; + if (comp != numbers::invalid_unsigned_int) + { + typename VectorType::value_type value = + internal::ElementAccess::get(vector, dof_indices[i]); + if (euler_dof_handler->get_fe().is_primitive(i)) + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; + ++v) + vertices[v][comp] += fe_values.shape_value(i, v) * value; + else + Assert(false, ExcNotImplemented()); + } + } return vertices; } @@ -2349,13 +2446,18 @@ MappingFEField::update_internal_dofs( ExcMessage("euler_dof_handler is empty")); typename DoFHandlerType::cell_iterator dof_cell(*cell, euler_dof_handler); - Assert(dof_cell->active() == true, ExcInactiveCell()); + Assert(uses_level_dofs || dof_cell->active() == true, ExcInactiveCell()); - dof_cell->get_dof_indices(data.local_dof_indices); + if (uses_level_dofs) + dof_cell->get_mg_dof_indices(data.local_dof_indices); + else + dof_cell->get_dof_indices(data.local_dof_indices); + const VectorType &vector = + uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0]; for (unsigned int i = 0; i < data.local_dof_values.size(); ++i) data.local_dof_values[i] = - internal::ElementAccess::get(*euler_vector, + internal::ElementAccess::get(vector, data.local_dof_indices[i]); } -- 2.39.5