From 603edc13ef3899caafc36b7106ab9af8ca6115f9 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Wed, 27 Sep 2017 16:25:27 +0200 Subject: [PATCH] extend MappingQEulerian to work with multi-grid solution vectors --- include/deal.II/fe/mapping_q_eulerian.h | 30 +- source/fe/mapping_q_eulerian.cc | 33 +- tests/mappings/mapping_q_eulerian_08.cc | 320 ++++++++++++++++++ ...ulerian_08.with_p4est=true.mpirun=1.output | 4 + ...ulerian_08.with_p4est=true.mpirun=3.output | 14 + 5 files changed, 380 insertions(+), 21 deletions(-) create mode 100644 tests/mappings/mapping_q_eulerian_08.cc create mode 100644 tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=1.output create mode 100644 tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=3.output diff --git a/include/deal.II/fe/mapping_q_eulerian.h b/include/deal.II/fe/mapping_q_eulerian.h index 5c1a7e130a..1b3117f153 100644 --- a/include/deal.II/fe/mapping_q_eulerian.h +++ b/include/deal.II/fe/mapping_q_eulerian.h @@ -43,11 +43,11 @@ template class Vector; *

Usage

* * The constructor of this class takes three arguments: the polynomial degree - * of the desire Qp mapping, a reference to the vector that defines the + * of the desired Qp mapping, a reference to the vector that defines the * mapping from the initial configuration to the current configuration, and a * reference to the DoFHandler. The most common case is to use the solution * vector for the problem under consideration as the shift vector. The key - * requirement is that the number of components of the given vector field be + * requirement is that the number of components of the given vector field must be * equal to (or possibly greater than) the number of space dimensions. If * there are more components than space dimensions (for example, if one is * working with a coupled problem where there are additional solution @@ -57,7 +57,7 @@ template class Vector; * the triangulation and associate the desired shift vector to it. * * Typically, the DoFHandler operates on a finite element that is constructed - * as a system element (FESystem) from continuous FE_Q() objects. An example + * as a system element (FESystem) from continuous FE_Q objects. An example * is shown below: * @code * FESystem fe(FE_Q(2), dim, FE_Q(1), 1); @@ -80,7 +80,7 @@ template class Vector; * that whenever you use this object, the given objects still represent valid * data. * - * To enable the use of the MappingQ1Eulerian class also in the context of + * To enable the use of the MappingQEulerian class also in the context of * parallel codes using the PETSc or Trilinos wrapper classes, the type * of the vector can be specified as template parameter VectorType. * @@ -102,10 +102,14 @@ public: * the second argument. The first dim components of this function will be * interpreted as the displacement we use in defining the mapping, relative * to the location of cells of the underlying triangulation. + * @param[in] level. Is the multi-grid level at which the mapping will + * be used. It is mainly used to check if the size of the @p euler_vector + * is consistent with the @p euler_dof_handler . */ MappingQEulerian (const unsigned int degree, const DoFHandler &euler_dof_handler, - const VectorType &euler_vector); + const VectorType &euler_vector, + const unsigned int level = numbers::invalid_unsigned_int); /** * Return the mapped vertices of the cell. For the current class, this @@ -125,15 +129,16 @@ public: Mapping *clone () const; /** - * Always returns @p false because MappingQ1Eulerian does not in general + * Always returns @p false because MappingQEulerian does not in general * preserve vertex locations (unless the translation vector happens to - * provide for zero displacements at vertex locations). + * provide zero displacements at vertex locations). */ virtual bool preserves_vertex_locations () const; /** - * Exception + * Exception which is thrown when the mapping is being evaluated at + * non-active cell. */ DeclException0 (ExcInactiveCell); @@ -167,6 +172,11 @@ protected: private: + /** + * Multigrid level at which the mapping is to be used. + */ + const unsigned int level; + /** * A class derived from MappingQGeneric that provides the generic mapping * with support points on boundary objects so that the corresponding Q3 @@ -202,7 +212,7 @@ private: compute_mapping_support_points(const typename Triangulation::cell_iterator &cell) const; /** - * Always returns @p false because MappingQ1Eulerian does not in general + * Always returns @p false because MappingQEulerianGeneric does not in general * preserve vertex locations (unless the translation vector happens to * provide for zero displacements at vertex locations). */ @@ -220,14 +230,12 @@ private: * Special quadrature rule used to define the support points in the * reference configuration. */ - class SupportQuadrature : public Quadrature { public: /** * Constructor, with an argument defining the desired polynomial degree. */ - SupportQuadrature (const unsigned int map_degree); }; diff --git a/source/fe/mapping_q_eulerian.cc b/source/fe/mapping_q_eulerian.cc index b5a3296689..d06ea777af 100644 --- a/source/fe/mapping_q_eulerian.cc +++ b/source/fe/mapping_q_eulerian.cc @@ -58,11 +58,13 @@ template MappingQEulerian:: MappingQEulerian (const unsigned int degree, const DoFHandler &euler_dof_handler, - const VectorType &euler_vector) + const VectorType &euler_vector, + const unsigned int level) : MappingQ(degree, true), euler_vector(&euler_vector), - euler_dof_handler(&euler_dof_handler) + euler_dof_handler(&euler_dof_handler), + level(level) { // reset the q1 mapping we use for interior cells (and previously // set by the MappingQ constructor) to a MappingQ1Eulerian with the @@ -168,10 +170,13 @@ std::vector > MappingQEulerian::MappingQEulerianGeneric:: compute_mapping_support_points (const typename Triangulation::cell_iterator &cell) const { - // first, basic assertion with respect to vector size, + const bool mg_vector = mapping_q_eulerian.level != numbers::invalid_unsigned_int; - const types::global_dof_index n_dofs = mapping_q_eulerian.euler_dof_handler->n_dofs(); + const types::global_dof_index n_dofs = mg_vector ? + mapping_q_eulerian.euler_dof_handler->n_dofs(mapping_q_eulerian.level) : + mapping_q_eulerian.euler_dof_handler->n_dofs(); const types::global_dof_index vector_size = mapping_q_eulerian.euler_vector->size(); + (void)n_dofs; (void)vector_size; @@ -182,20 +187,19 @@ compute_mapping_support_points (const typename Triangulation::cell typename DoFHandler::cell_iterator dof_cell(*cell, mapping_q_eulerian.euler_dof_handler); - Assert (dof_cell->active() == true, ExcInactiveCell()); + Assert (mg_vector || dof_cell->active() == true, ExcInactiveCell()); // our quadrature rule is chosen so that each quadrature point corresponds - // to a support point in the undeformed configuration. we can then query + // to a support point in the undeformed configuration. We can then query // the given displacement field at these points to determine the shift // vector that maps the support points to the deformed configuration. - // we assume that the given field contains dim displacement components, but + // We assume that the given field contains dim displacement components, but // that there may be other solution components as well (e.g. pressures). // this class therefore assumes that the first dim components represent the // actual shift vector we need, and simply ignores any components after - // that. this implies that the user should order components appropriately, + // that. This implies that the user should order components appropriately, // or create a separate dof handler for the displacements. - const unsigned int n_support_pts = support_quadrature.size(); const unsigned int n_components = mapping_q_eulerian.euler_dof_handler->get_fe(0).n_components(); @@ -205,12 +209,21 @@ compute_mapping_support_points (const typename Triangulation::cell shift_vector(n_support_pts, Vector(n_components)); + std::vector dof_indices(mapping_q_eulerian.euler_dof_handler->get_fe(0).dofs_per_cell); // fill shift vector for each support point using an fe_values object. make // sure that the fe_values variable isn't used simultaneously from different // threads Threads::Mutex::ScopedLock lock(fe_values_mutex); fe_values.reinit(dof_cell); - fe_values.get_function_values(*mapping_q_eulerian.euler_vector, shift_vector); + if (mg_vector) + { + dof_cell->get_mg_dof_indices(dof_indices); + fe_values.get_function_values(*mapping_q_eulerian.euler_vector, + dof_indices, + shift_vector); + } + else + fe_values.get_function_values(*mapping_q_eulerian.euler_vector, shift_vector); // and finally compute the positions of the support points in the deformed // configuration. diff --git a/tests/mappings/mapping_q_eulerian_08.cc b/tests/mappings/mapping_q_eulerian_08.cc new file mode 100644 index 0000000000..2f1bc4fa32 --- /dev/null +++ b/tests/mappings/mapping_q_eulerian_08.cc @@ -0,0 +1,320 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test that MappingQEulerian works in parallel with geometric multigrids. +// We apply a simple linear deformation which can be represented exactly +// at the coarse level. +// +// inspired by _07 + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + + + + +using namespace dealii; + + +template +class Displacement : public Function +{ +public: + Displacement() : + Function(dim) + {} + + double value (const Point &p, + const unsigned int component) const + { + return (0.1 + 2.*component) + (0.2 + 3.*component) * p[component]; + } + + template + Tensor<1,dim, VectorizedArray > + shift_value(const Point> &p_vec) const + { + Tensor<1,dim, VectorizedArray> shift_vec; + Point p; + for (unsigned int v = 0; v < VectorizedArray::n_array_elements; ++v) + { + for (unsigned int d = 0; d < dim; ++d) + p[d] = p_vec[d][v]; + + for (unsigned int d = 0; d < dim; ++d) + shift_vec[d][v] = this->value(p,d); + } + + return shift_vec; + + } +}; + + +template +void test (const unsigned int n_ref = 0) +{ + Displacement displacement_function; + const unsigned int euler_fe_degree=2; + + deallog << "dim=" << dim << std::endl; + MPI_Comm mpi_communicator(MPI_COMM_WORLD); + unsigned int myid = Utilities::MPI::this_mpi_process (mpi_communicator); + unsigned int numproc = Utilities::MPI::n_mpi_processes (mpi_communicator); + + deallog << "numproc=" << numproc << std::endl; + + parallel::distributed::Triangulation + triangulation (mpi_communicator, + Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation::construct_multigrid_hierarchy); + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global(3); + + // do some adaptive refinement + for (unsigned int ref=0; ref::active_cell_iterator cell=triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (cell->is_locally_owned() && + (cell->center().norm() < 0.5 && (cell->level() < 5 || + cell->center().norm() > 0.45) + || + (dim == 2 && cell->center().norm() > 1.2))) + cell->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + } + + FE_Q fe (fe_degree); + DoFHandler dof_handler (triangulation); + dof_handler.distribute_dofs(fe); + dof_handler.distribute_mg_dofs(); + + // quadrature for MatrixFree, not related to the degree in the + // FE displacement field. + QGauss<1> quadrature_formula(n_q_points); + + + FESystem fe_euler(FE_Q(euler_fe_degree), dim); + DoFHandler dof_handler_euler(triangulation); + dof_handler_euler.distribute_dofs(fe_euler); + dof_handler_euler.distribute_mg_dofs (); + + const IndexSet &locally_owned_dofs_euler = dof_handler_euler.locally_owned_dofs (); + IndexSet locally_relevant_dofs_euler; + DoFTools::extract_locally_relevant_dofs (dof_handler_euler,locally_relevant_dofs_euler); + + const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs (); + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs (dof_handler,locally_relevant_dofs); + + // constraints: + ConstraintMatrix constraints_euler; + constraints_euler.reinit (locally_relevant_dofs_euler); + DoFTools::make_hanging_node_constraints (dof_handler_euler, constraints_euler); + constraints_euler.close (); + + ConstraintMatrix constraints; + constraints.reinit (locally_relevant_dofs); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints_euler.close (); + + // MG constraints: + MGConstrainedDoFs mg_constrained_dofs_euler; + mg_constrained_dofs_euler.initialize(dof_handler_euler); + + // Displacement vector + LinearAlgebra::distributed::Vector displacement; + displacement.reinit(locally_owned_dofs_euler, + locally_relevant_dofs_euler, + mpi_communicator); + + VectorTools::project,dim> + (dof_handler_euler, + constraints_euler, + QGauss(n_q_points), + displacement_function, + displacement); + displacement.update_ghost_values(); + + MGTransferMatrixFree mg_transfer_euler(mg_constrained_dofs_euler); + mg_transfer_euler.build(dof_handler); + + // now the core of the test: + const unsigned int max_level = dof_handler.get_triangulation().n_global_levels()-1; + const unsigned int min_level = 0; + MGLevelObject> displacement_level(min_level, max_level); + mg_transfer_euler.interpolate_to_mg(dof_handler_euler,displacement_level, displacement); + + // First, check fine-level only: + { + MappingQEulerian> + euler_fine(euler_fe_degree, dof_handler_euler, displacement); + + + + MatrixFree matrix_free_euler; + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::partition_color; + data.mapping_update_flags = update_values | update_gradients | update_JxW_values | update_quadrature_points; + matrix_free_euler.reinit (euler_fine, dof_handler, constraints, quadrature_formula, data); + + MatrixFree matrix_free; + matrix_free.reinit (dof_handler, constraints, quadrature_formula, data); + + + // test fine-level mapping: + { + FEEvaluation fe_eval_euler(matrix_free_euler); + FEEvaluation fe_eval(matrix_free); + const unsigned int n_cells = matrix_free_euler.n_macro_cells(); + Assert (matrix_free_euler.n_macro_cells() == matrix_free.n_macro_cells(), + ExcInternalError()); + const unsigned int nqp = fe_eval.n_q_points; + for (unsigned int cell=0; cell dist = v1.distance(v2); + for (unsigned int v = 0; v < VectorizedArray::n_array_elements; ++v) + AssertThrow (dist[v] < 1e-8, + ExcMessage("distance: " + std::to_string(dist[v]))); + } + } + } + + } + + // now go through all GMG levels: + std::set dirichlet_boundary_ids; + dirichlet_boundary_ids.insert(0); + MGConstrainedDoFs mg_constrained_dofs; + mg_constrained_dofs.initialize(dof_handler); + mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,dirichlet_boundary_ids); + + for (unsigned int level = min_level; level <= max_level; ++level) + { + typename MatrixFree::AdditionalData mg_additional_data; + mg_additional_data.tasks_parallel_scheme = MatrixFree::AdditionalData::partition_color; + mg_additional_data.level_mg_handler = level; + mg_additional_data.mapping_update_flags = update_values | update_gradients | update_JxW_values | update_quadrature_points; + + ConstraintMatrix level_constraints; + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler, level, + relevant_dofs); + level_constraints.reinit(relevant_dofs); + level_constraints.add_lines(mg_constrained_dofs.get_boundary_indices(level)); + level_constraints.close(); + + MatrixFree mg_level_euler; + + MappingQEulerian> + euler_level(euler_fe_degree, dof_handler_euler, displacement_level[level], level); + + mg_level_euler.reinit (euler_level, dof_handler, level_constraints, quadrature_formula, + mg_additional_data); + + MatrixFree mg_level; + mg_level.reinit (dof_handler, level_constraints, quadrature_formula, + mg_additional_data); + + // go through all cells and quadrature points: + { + FEEvaluation fe_eval_euler(mg_level_euler); + FEEvaluation fe_eval(mg_level); + const unsigned int n_cells = mg_level_euler.n_macro_cells(); + Assert (mg_level_euler.n_macro_cells() == mg_level.n_macro_cells(), + ExcInternalError()); + const unsigned int nqp = fe_eval.n_q_points; + for (unsigned int cell=0; cell dist = v1.distance(v2); + for (unsigned int v = 0; v < VectorizedArray::n_array_elements; ++v) + AssertThrow (dist[v] < 1e-8, + ExcMessage("Level " + std::to_string(level) + + " distance: " + std::to_string(dist[v]))); + } + } + } + } + + deallog << "Ok" << std::endl; +} + + + +int main (int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + MPILogInitAll log; + + test<2> (); +} + + diff --git a/tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=1.output b/tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..047c058b1b --- /dev/null +++ b/tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=1.output @@ -0,0 +1,4 @@ + +DEAL:0::dim=2 +DEAL:0::numproc=1 +DEAL:0::Ok diff --git a/tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=3.output b/tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..d171703b0b --- /dev/null +++ b/tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=3.output @@ -0,0 +1,14 @@ + +DEAL:0::dim=2 +DEAL:0::numproc=3 +DEAL:0::Ok + +DEAL:1::dim=2 +DEAL:1::numproc=3 +DEAL:1::Ok + + +DEAL:2::dim=2 +DEAL:2::numproc=3 +DEAL:2::Ok + -- 2.39.5