From: Denis Davydov Date: Mon, 16 Oct 2017 11:06:00 +0000 (+0200) Subject: add another test for MappingQEulerian by comparing to manually deformed tria X-Git-Tag: v9.0.0-rc1~945^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F5252%2Fhead;p=dealii.git add another test for MappingQEulerian by comparing to manually deformed tria --- diff --git a/tests/mappings/mapping_q_eulerian_09.cc b/tests/mappings/mapping_q_eulerian_09.cc new file mode 100644 index 0000000000..3082feb03a --- /dev/null +++ b/tests/mappings/mapping_q_eulerian_09.cc @@ -0,0 +1,208 @@ +// --------------------------------------------------------------------- +// +// 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 matrix-free Laplace and Mass operators with MappingQEulerian +// by comparing to the results obtained from a deformed mesh. +// As a displacement function use: x exp (-2|x|) + + +#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 + +// Displacement to be applied +template +class Displacement : public Function +{ +public: + Displacement() : + Function(dim) + {} + + double value (const Point &p, + const unsigned int component) const + { + return p[component]*std::exp(-2.*std::abs(p[component])); + } +}; + + +template +void test () +{ + 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); + + parallel::distributed::Triangulation + triangulation (mpi_communicator, + Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation::construct_multigrid_hierarchy); + + GridGenerator::subdivided_hyper_cube (triangulation, 8, -5, 5); + + FE_Q fe (fe_degree); + DoFHandler dof_handler (triangulation); + dof_handler.distribute_dofs(fe); + dof_handler.distribute_mg_dofs(); + + const unsigned int euler_fe_degree=1; + 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 (); + + // IndexSets and constraints + 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 (); + + // Displacement vector + LinearAlgebra::distributed::Vector displacement; + displacement.reinit(locally_owned_dofs_euler, + locally_relevant_dofs_euler, + mpi_communicator); + displacement = 0.; + + Displacement displacement_function; + + // first, move via mapping: + VectorTools::interpolate(dof_handler_euler, + displacement_function, + displacement); + displacement.compress(VectorOperation::insert); + displacement.update_ghost_values(); + + // The core : compute matrix-vector products with mass and laplace operators + QGauss<1> quadrature_formula(n_q_points); + + MappingQEulerian> + euler_mapping(euler_fe_degree, dof_handler_euler, displacement); + + std::shared_ptr> matrix_free(new MatrixFree ()); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::partition_color; + data.mapping_update_flags = update_values | update_gradients | update_JxW_values; + matrix_free->reinit (euler_mapping, dof_handler, constraints, quadrature_formula, data); + + LinearAlgebra::distributed::Vector src, dst1, dst2, dst3, dst4; + matrix_free->initialize_dof_vector(src); + matrix_free->initialize_dof_vector(dst1); + matrix_free->initialize_dof_vector(dst2); + matrix_free->initialize_dof_vector(dst3); + matrix_free->initialize_dof_vector(dst4); + + for (unsigned int i = 0; i < src.local_size(); ++i) + src.local_element(i) = (NumberType)Testing::rand()/RAND_MAX; + + MatrixFreeOperators::MassOperator> mf_mass; + MatrixFreeOperators::MassOperator> mf_laplace; + mf_mass.initialize(matrix_free); + mf_laplace.initialize(matrix_free); + + mf_mass.vmult(dst1,src); + mf_laplace.vmult(dst3,src); + + // now move manually + displacement = 0.; + displacement.update_ghost_values(); + { + typename DoFHandler::cell_iterator cell=dof_handler.begin_active(), + endc=dof_handler.end(); + std::vector vertex_touched (triangulation.n_vertices(), false); + for (cell=dof_handler.begin_active(); cell!=endc; ++cell) + for (unsigned int vertex_no=0; + vertex_no::vertices_per_cell; ++vertex_no) + if (vertex_touched[cell->vertex_index(vertex_no)] == false) + { + Point &v = cell->vertex(vertex_no); + Tensor<1,dim> d; + for (unsigned int i = 0; i < dim; ++i) + d[i] = displacement_function.value(v, i); + + v+= d; + vertex_touched[cell->vertex_index(vertex_no)] = true; + } + } + // minimize the data that is re-computed + data.initialize_indices = false; + matrix_free->reinit (euler_mapping, dof_handler, constraints, quadrature_formula, data); + + mf_mass.vmult(dst2,src); + mf_laplace.vmult(dst4,src); + + deallog << "Mass operator: " << std::endl + << "l2_norm: " << std::abs(dst1.l2_norm() - dst2.l2_norm()) << std::endl + << "l1_norm: " << std::abs(dst1.l1_norm() - dst2.l1_norm()) << std::endl + << "linfty_norm: " << std::abs(dst1.linfty_norm() - dst2.linfty_norm()) << std::endl + << "Laplace operator: " << std::endl + << "l2_norm: " << std::abs(dst3.l2_norm() - dst4.l2_norm()) << std::endl + << "l1_norm: " << std::abs(dst3.l1_norm() - dst4.l1_norm()) << std::endl + << "linfty_norm: " << std::abs(dst3.linfty_norm() - dst4.linfty_norm()) << std::endl; + +} + +int main (int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads()); + + mpi_initlog(); + + deallog.push("2d"); + test<2,2>(); + deallog.pop(); + + deallog.push("3d"); + test<3,2>(); + deallog.pop(); +} diff --git a/tests/mappings/mapping_q_eulerian_09.with_p4est=true.mpirun=1.output b/tests/mappings/mapping_q_eulerian_09.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..a0dcd70b75 --- /dev/null +++ b/tests/mappings/mapping_q_eulerian_09.with_p4est=true.mpirun=1.output @@ -0,0 +1,17 @@ + +DEAL:2d::Mass operator: +DEAL:2d::l2_norm: 0.00000 +DEAL:2d::l1_norm: 0.00000 +DEAL:2d::linfty_norm: 0.00000 +DEAL:2d::Laplace operator: +DEAL:2d::l2_norm: 0.00000 +DEAL:2d::l1_norm: 0.00000 +DEAL:2d::linfty_norm: 0.00000 +DEAL:3d::Mass operator: +DEAL:3d::l2_norm: 0.00000 +DEAL:3d::l1_norm: 0.00000 +DEAL:3d::linfty_norm: 0.00000 +DEAL:3d::Laplace operator: +DEAL:3d::l2_norm: 0.00000 +DEAL:3d::l1_norm: 0.00000 +DEAL:3d::linfty_norm: 0.00000 diff --git a/tests/mappings/mapping_q_eulerian_09.with_p4est=true.mpirun=3.output b/tests/mappings/mapping_q_eulerian_09.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..a0dcd70b75 --- /dev/null +++ b/tests/mappings/mapping_q_eulerian_09.with_p4est=true.mpirun=3.output @@ -0,0 +1,17 @@ + +DEAL:2d::Mass operator: +DEAL:2d::l2_norm: 0.00000 +DEAL:2d::l1_norm: 0.00000 +DEAL:2d::linfty_norm: 0.00000 +DEAL:2d::Laplace operator: +DEAL:2d::l2_norm: 0.00000 +DEAL:2d::l1_norm: 0.00000 +DEAL:2d::linfty_norm: 0.00000 +DEAL:3d::Mass operator: +DEAL:3d::l2_norm: 0.00000 +DEAL:3d::l1_norm: 0.00000 +DEAL:3d::linfty_norm: 0.00000 +DEAL:3d::Laplace operator: +DEAL:3d::l2_norm: 0.00000 +DEAL:3d::l1_norm: 0.00000 +DEAL:3d::linfty_norm: 0.00000