]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add another test for MappingQEulerian by comparing to manually deformed tria 5252/head
authorDenis Davydov <davydden@gmail.com>
Mon, 16 Oct 2017 11:06:00 +0000 (13:06 +0200)
committerDenis Davydov <davydden@gmail.com>
Mon, 16 Oct 2017 13:13:47 +0000 (15:13 +0200)
tests/mappings/mapping_q_eulerian_09.cc [new file with mode: 0644]
tests/mappings/mapping_q_eulerian_09.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/mappings/mapping_q_eulerian_09.with_p4est=true.mpirun=3.output [new file with mode: 0644]

diff --git a/tests/mappings/mapping_q_eulerian_09.cc b/tests/mappings/mapping_q_eulerian_09.cc
new file mode 100644 (file)
index 0000000..3082feb
--- /dev/null
@@ -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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/thread_management.h>
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/fe/mapping_q_eulerian.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/operators.h>
+#include <deal.II/numerics/vector_tools.h>
+
+// Displacement to be applied
+template <int dim>
+class Displacement : public Function<dim>
+{
+public:
+  Displacement() :
+    Function<dim>(dim)
+  {}
+
+  double value (const Point<dim> &p,
+                const unsigned int component) const
+  {
+    return p[component]*std::exp(-2.*std::abs(p[component]));
+  }
+};
+
+
+template <int dim, int fe_degree=2, int n_q_points=fe_degree+1, typename NumberType=double>
+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<dim>
+  triangulation (mpi_communicator,
+                 Triangulation<dim>::limit_level_difference_at_vertices,
+                 parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+
+  GridGenerator::subdivided_hyper_cube (triangulation, 8, -5, 5);
+
+  FE_Q<dim> fe (fe_degree);
+  DoFHandler<dim> dof_handler (triangulation);
+  dof_handler.distribute_dofs(fe);
+  dof_handler.distribute_mg_dofs();
+
+  const unsigned int euler_fe_degree=1;
+  FESystem<dim> fe_euler(FE_Q<dim>(euler_fe_degree), dim);
+  DoFHandler<dim> 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<NumberType> displacement;
+  displacement.reinit(locally_owned_dofs_euler,
+                      locally_relevant_dofs_euler,
+                      mpi_communicator);
+  displacement = 0.;
+
+  Displacement<dim> 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<dim,LinearAlgebra::distributed::Vector<NumberType>>
+      euler_mapping(euler_fe_degree, dof_handler_euler, displacement);
+
+  std::shared_ptr<MatrixFree<dim,NumberType>> matrix_free(new MatrixFree<dim,NumberType> ());
+  typename MatrixFree<dim,NumberType>::AdditionalData data;
+  data.tasks_parallel_scheme = MatrixFree<dim,NumberType>::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<NumberType> 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<dim,fe_degree, n_q_points, 1, LinearAlgebra::distributed::Vector<NumberType>> mf_mass;
+  MatrixFreeOperators::MassOperator<dim,fe_degree, n_q_points, 1, LinearAlgebra::distributed::Vector<NumberType>> 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<dim>::cell_iterator cell=dof_handler.begin_active(),
+                                            endc=dof_handler.end();
+    std::vector<bool> vertex_touched (triangulation.n_vertices(), false);
+    for (cell=dof_handler.begin_active(); cell!=endc; ++cell)
+      for (unsigned int vertex_no=0;
+           vertex_no<GeometryInfo<dim>::vertices_per_cell; ++vertex_no)
+        if (vertex_touched[cell->vertex_index(vertex_no)] == false)
+          {
+            Point<dim> &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 (file)
index 0000000..a0dcd70
--- /dev/null
@@ -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 (file)
index 0000000..a0dcd70
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.