From: Denis Davydov Date: Sat, 29 Oct 2016 14:40:31 +0000 (+0200) Subject: add VectorTools::project() for scalar-valued quadrature data X-Git-Tag: v8.5.0-rc1~519^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9588242b94fc463fb64f8b52a50d381734f82095;p=dealii.git add VectorTools::project() for scalar-valued quadrature data --- diff --git a/cmake/config/template-arguments.in b/cmake/config/template-arguments.in index f2f37a3ecd..7830833e12 100644 --- a/cmake/config/template-arguments.in +++ b/cmake/config/template-arguments.in @@ -69,6 +69,21 @@ REAL_SERIAL_VECTORS := { Vector; @DEAL_II_EXPAND_PETSC_MPI_BLOCKVECTOR_REAL@; } +REAL_NONBLOCK_VECTORS := { Vector; + Vector ; + + LinearAlgebra::Vector; + LinearAlgebra::Vector ; + + LinearAlgebra::distributed::Vector; + LinearAlgebra::distributed::Vector ; + + @DEAL_II_EXPAND_TRILINOS_VECTOR@; + @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@; + @DEAL_II_EXPAND_PETSC_VECTOR_REAL@; + @DEAL_II_EXPAND_PETSC_MPI_VECTOR_REAL@; + } + EXTERNAL_SEQUENTIAL_VECTORS := { @DEAL_II_EXPAND_TRILINOS_VECTOR@; @DEAL_II_EXPAND_TRILINOS_BLOCKVECTOR@; @DEAL_II_EXPAND_PETSC_VECTOR@; diff --git a/doc/news/changes.h b/doc/news/changes.h index 7c697cdc63..bfe83b245a 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -405,6 +405,11 @@ inconvenience this causes.

Specific improvements

    +
  1. New: Add VectorTools::project() to do L2 projection + of scalar-valued quadrature point data in parallel. +
    + (Denis Davydov, 2016/10/28) +
  3. Fixed: Increased precision of timesteps in DataOutInterface::write_pvd_record().
    diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index 4b24acebca..7071137f61 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include @@ -35,6 +36,7 @@ DEAL_II_NAMESPACE_OPEN template class Function; template class Quadrature; template class QGauss; +template class MatrixFree; template class Vector; template class FullMatrix; @@ -856,6 +858,71 @@ namespace VectorTools Quadrature(0)), const bool project_to_boundary_first = false); + /** + * The same as above for projection of scalar-valued quadrature data. + * The user provided function should return a value at the quadrature point + * based on the cell iterator and quadrature number and of course should be + * consistent with the provided @p quadrature object, which will be used + * to assemble the right-hand-side. + * + * This function can be used with lambdas: + * @code + * VectorTools::project + * (mapping, + * dof_handler, + * constraints, + * quadrature_formula, + * [=] (const typename DoFHandler::active_cell_iterator & cell, const unsigned int q) -> double + * { return qp_data.get_data(cell)[q]->density; }, + * field); + * @endcode + * where qp_data is a CellDataStorage object, which stores + * quadrature point data. + */ + template + void project (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const std_cxx11::function< typename VectorType::value_type (const typename DoFHandler::active_cell_iterator &, const unsigned int)> func, + VectorType &vec_result); + + /** + * The same as above for projection of scalar-valued MatrixFree quadrature + * data. + * The user provided function @p func should return a VectorizedArray value + * at the quadrature point based on the cell number and quadrature number and + * should be consistent with the @p n_q_points_1d. + * + * This function can be used with lambdas: + * @code + * VectorTools::project + * (matrix_free_data, + * constraints, + * 3, + * [=] (const unsigned int cell, const unsigned int q) -> VectorizedArray + * { return qp_data(cell,q); }, + * field); + * @endcode + * where qp_data is a an object of type Table<2, VectorizedArray >, + * which stores quadrature point data. + */ + template + void project (const MatrixFree &data, + const ConstraintMatrix &constraints, + const unsigned int n_q_points_1d, + const std_cxx11::function< VectorizedArray (const unsigned int, const unsigned int)> func, + VectorType &vec_result); + + /** + * Same as above but for n_q_points_1d = matrix_free.get_dof_handler().get_fe().degree+1. + */ + template + void project (const MatrixFree &data, + const ConstraintMatrix &constraints, + const std_cxx11::function< VectorizedArray (const unsigned int, const unsigned int)> func, + VectorType &vec_result); + /** * Compute Dirichlet boundary conditions. This function makes up a map of * degrees of freedom subject to Dirichlet boundary conditions and the diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index eae1a44941..0040f2cacd 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include @@ -37,6 +38,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -881,9 +885,282 @@ namespace VectorTools for (unsigned int i=0; i + void project_parallel (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const std_cxx11::function< typename VectorType::value_type (const typename DoFHandler::active_cell_iterator &, const unsigned int)> func, + VectorType &vec_result) + { + const parallel::Triangulation *parallel_tria = + dynamic_cast*>(&dof.get_tria()); + + typedef typename VectorType::value_type Number; + Assert (dof.get_fe().n_components() == 1, + ExcDimensionMismatch(dof.get_fe().n_components(), + 1)); + Assert (vec_result.size() == dof.n_dofs(), + ExcDimensionMismatch (vec_result.size(), dof.n_dofs())); + Assert (dof.get_fe().degree == fe_degree, + ExcDimensionMismatch(fe_degree, dof.get_fe().degree)); + + // set up mass matrix and right hand side + typename MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_color; + if (parallel_tria !=0) + additional_data.mpi_communicator = parallel_tria->get_communicator(); + additional_data.mapping_update_flags = (update_values | update_JxW_values); + MatrixFree matrix_free; + matrix_free.reinit (mapping, dof, constraints, + QGauss<1>(fe_degree+1), additional_data); + typedef MatrixFreeOperators::MassOperator MatrixType; + MatrixType mass_matrix; + mass_matrix.initialize(matrix_free); + mass_matrix.compute_diagonal(); + + typedef LinearAlgebra::distributed::Vector LocalVectorType; + LocalVectorType vec, rhs, inhomogeneities; + matrix_free.initialize_dof_vector(vec); + matrix_free.initialize_dof_vector(rhs); + matrix_free.initialize_dof_vector(inhomogeneities); + constraints.distribute(inhomogeneities); + inhomogeneities*=-1.; + + // assemble right hand side: + { + FEValues fe_values (mapping, dof.get_fe(), quadrature, + update_values | update_JxW_values); + + const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; + const unsigned int n_q_points = quadrature.size(); + Vector cell_rhs (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + typename DoFHandler::active_cell_iterator + cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + cell_rhs = 0; + fe_values.reinit (cell); + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (cell_rhs, + local_dof_indices, + rhs); + } + rhs.compress (VectorOperation::add); + } + + mass_matrix.vmult_add(rhs, inhomogeneities); + + //now invert the matrix + ReductionControl control(rhs.size(), 0., 1e-12, false, false); + SolverCG > cg(control); + typename PreconditionChebyshev::AdditionalData data; + data.matrix_diagonal_inverse = mass_matrix.get_matrix_diagonal_inverse(); + PreconditionChebyshev preconditioner; + preconditioner.initialize(mass_matrix, data); + cg.solve (mass_matrix, vec, rhs, preconditioner); + vec+=inhomogeneities; + + constraints.distribute (vec); + + const IndexSet locally_owned_dofs = dof.locally_owned_dofs(); + IndexSet::ElementIterator it = locally_owned_dofs.begin(); + for (; it!=locally_owned_dofs.end(); ++it) + vec_result(*it) = vec(*it); + vec_result.compress(VectorOperation::insert); + } + + + + template + void project_parallel (const MatrixFree &matrix_free, + const ConstraintMatrix &constraints, + const std_cxx11::function< VectorizedArray (const unsigned int, const unsigned int)> func, + VectorType &vec_result) + { + const DoFHandler &dof = matrix_free.get_dof_handler(); + + typedef typename VectorType::value_type Number; + Assert (dof.get_fe().n_components() == 1, + ExcDimensionMismatch(dof.get_fe().n_components(), + 1)); + Assert (vec_result.size() == dof.n_dofs(), + ExcDimensionMismatch (vec_result.size(), dof.n_dofs())); + Assert (dof.get_fe().degree == fe_degree, + ExcDimensionMismatch(fe_degree, dof.get_fe().degree)); + + typedef MatrixFreeOperators::MassOperator MatrixType; + MatrixType mass_matrix; + mass_matrix.initialize(matrix_free); + mass_matrix.compute_diagonal(); + + typedef LinearAlgebra::distributed::Vector LocalVectorType; + LocalVectorType vec, rhs, inhomogeneities; + matrix_free.initialize_dof_vector(vec); + matrix_free.initialize_dof_vector(rhs); + matrix_free.initialize_dof_vector(inhomogeneities); + constraints.distribute(inhomogeneities); + inhomogeneities*=-1.; + + // assemble right hand side: + { + FEEvaluation fe_eval(matrix_free); + const unsigned int n_cells = matrix_free.n_macro_cells(); + const unsigned int n_q_points = fe_eval.n_q_points; + + for (unsigned int cell=0; cell > cg(control); + typename PreconditionChebyshev::AdditionalData data; + data.matrix_diagonal_inverse = mass_matrix.get_matrix_diagonal_inverse(); + PreconditionChebyshev preconditioner; + preconditioner.initialize(mass_matrix, data); + cg.solve (mass_matrix, vec, rhs, preconditioner); + vec+=inhomogeneities; + + constraints.distribute (vec); + + const IndexSet locally_owned_dofs = dof.locally_owned_dofs(); + IndexSet::ElementIterator it = locally_owned_dofs.begin(); + for (; it!=locally_owned_dofs.end(); ++it) + vec_result(*it) = vec(*it); + vec_result.compress(VectorOperation::insert); + } + } + + + + template + void project (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const std_cxx11::function< typename VectorType::value_type (const typename DoFHandler::active_cell_iterator &, const unsigned int)> func, + VectorType &vec_result) + { + switch (dof.get_fe().degree) + { + case 1: + project_parallel (mapping,dof,constraints,quadrature,func,vec_result); + break; + case 2: + project_parallel (mapping,dof,constraints,quadrature,func,vec_result); + break; + case 3: + project_parallel (mapping,dof,constraints,quadrature,func,vec_result); + break; + case 4: + project_parallel (mapping,dof,constraints,quadrature,func,vec_result); + break; + case 5: + project_parallel (mapping,dof,constraints,quadrature,func,vec_result); + break; + case 6: + project_parallel (mapping,dof,constraints,quadrature,func,vec_result); + break; + case 7: + project_parallel (mapping,dof,constraints,quadrature,func,vec_result); + break; + case 8: + project_parallel (mapping,dof,constraints,quadrature,func,vec_result); + break; + default: + Assert (false, ExcNotImplemented()); + } } + + template + void project (const MatrixFree &matrix_free, + const ConstraintMatrix &constraints, + const unsigned int n_q_points_1d, + const std_cxx11::function< VectorizedArray (const unsigned int, const unsigned int)> func, + VectorType &vec_result) + { + const unsigned int fe_degree = matrix_free.get_dof_handler().get_fe().degree; + + Assert (fe_degree+1 == n_q_points_1d, + ExcNotImplemented()); + + switch (fe_degree) + { + case 1: + project_parallel (matrix_free,constraints,func,vec_result); + break; + case 2: + project_parallel (matrix_free,constraints,func,vec_result); + break; + case 3: + project_parallel (matrix_free,constraints,func,vec_result); + break; + case 4: + project_parallel (matrix_free,constraints,func,vec_result); + break; + case 5: + project_parallel (matrix_free,constraints,func,vec_result); + break; + case 6: + project_parallel (matrix_free,constraints,func,vec_result); + break; + case 7: + project_parallel (matrix_free,constraints,func,vec_result); + break; + case 8: + project_parallel (matrix_free,constraints,func,vec_result); + break; + default: + Assert (false, ExcNotImplemented()); + } + } + + + + template + void project (const MatrixFree &matrix_free, + const ConstraintMatrix &constraints, + const std_cxx11::function< VectorizedArray (const unsigned int, const unsigned int)> func, + VectorType &vec_result) + { + project (matrix_free, + constraints, + matrix_free.get_dof_handler().get_fe().degree+1, + func, + vec_result); + } + + + template void project (const Mapping &mapping, const DoFHandler &dof, diff --git a/source/numerics/vector_tools_project.inst.in b/source/numerics/vector_tools_project.inst.in index 27e1397b54..1c2d2f4141 100644 --- a/source/numerics/vector_tools_project.inst.in +++ b/source/numerics/vector_tools_project.inst.in @@ -67,3 +67,36 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens \} #endif } + +for (VEC: REAL_NONBLOCK_VECTORS; deal_II_dimension : DIMENSIONS) +{ + namespace VectorTools \{ + + template + void project( + const Mapping &, + const DoFHandler &, + const ConstraintMatrix &, + const Quadrature &, + const std_cxx11::function::active_cell_iterator &, const unsigned int)>, + VEC &); + + template + void project( + const MatrixFree &matrix_free, + const ConstraintMatrix &constraints, + const std_cxx11::function< VectorizedArray (const unsigned int, const unsigned int)>, + VEC &); + + template + void project( + const MatrixFree &matrix_free, + const ConstraintMatrix &constraints, + const unsigned int, + const std_cxx11::function< VectorizedArray (const unsigned int, const unsigned int)>, + VEC &); + + \} + +} + diff --git a/tests/numerics/project_parallel_qp_common.h b/tests/numerics/project_parallel_qp_common.h new file mode 100644 index 0000000000..4d9cb56b24 --- /dev/null +++ b/tests/numerics/project_parallel_qp_common.h @@ -0,0 +1,245 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +// common framework to check whether an element of polynomial order p can +// represent functions of order q for projection from quadrature points. + +#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 + +template +class F : public Function +{ +public: + F (const unsigned int q, + const unsigned int n_components) + : + Function(n_components), + q(q) + {} + + virtual double value (const Point &p, + const unsigned int component = 0) const + { + Assert ((component == 0) && (this->n_components == 1), + ExcInternalError()); + double val = 0; + for (unsigned int d=0; d value(const Point> &p_vec) const + { + VectorizedArray res = make_vectorized_array (0.); + 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]; + res[v] = value(p); + } + return res; + } + + +private: + const unsigned int q; +}; + +struct QData +{ + double density; +}; + + +template +void do_project (const parallel::distributed::Triangulation &triangulation, + const FiniteElement &fe, + const unsigned int p) +{ + AssertThrow(fe.n_components() ==1, + ExcNotImplemented()); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs (fe); + + deallog << "n_cells=" << triangulation.n_global_active_cells() << std::endl; + deallog << "n_dofs=" << dof_handler.n_dofs() << std::endl; + + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs (dof_handler, + locally_relevant_dofs); + + ConstraintMatrix constraints; + constraints.reinit (locally_relevant_dofs); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + for (unsigned int q=0; q<=p; ++q) + { + // setup quadrature data: + F function(q, fe.n_components()); + VectorType field_relevant; + double field_l2_norm = 0.; + { + QGauss quadrature_formula(p+2); + CellDataStorage::cell_iterator,QData> qp_manager; + { + FEValues fe_values (fe, quadrature_formula, + update_quadrature_points); + + const unsigned int n_q_points = quadrature_formula.size(); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + fe_values.reinit (cell); + qp_manager.initialize(cell, quadrature_formula.size()); + std::vector > qpd = qp_manager.get_data(cell); + for (unsigned int q=0; qdensity = function.value(fe_values.quadrature_point(q)); + } + } + + VectorType field(dof_handler.locally_owned_dofs(), + MPI_COMM_WORLD); + VectorTools::project + (MappingQGeneric(1), + dof_handler, + constraints, + quadrature_formula, + [=] (const typename DoFHandler::active_cell_iterator & cell, const unsigned int q) -> double { return qp_manager.get_data(cell)[q]->density; }, + field); + + field_relevant.reinit(dof_handler.locally_owned_dofs(), + locally_relevant_dofs, + MPI_COMM_WORLD); + field_relevant = field; + + field_l2_norm = field.l2_norm(); + } + + // L2 norm of the difference between FE field and the function + double L2_norm = 0.; + { + QGauss quadrature_formula_error(std::max(p,q)+1); + FEValues fe_values (fe, quadrature_formula_error, + update_values | update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula_error.size(); + std::vector values (n_q_points); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + fe_values.reinit (cell); + fe_values.get_function_values(field_relevant, + values); + + for (unsigned int q_point=0; q_point(values[q_point] - function.value(fe_values.quadrature_point(q_point))) * + fe_values.JxW(q_point); + + } + } + + L2_norm = Utilities::MPI::sum(L2_norm, MPI_COMM_WORLD); + L2_norm = std::sqrt(L2_norm); + + deallog << fe.get_name() << ", P_" << q + << ", rel. error=" << L2_norm / field_l2_norm + << std::endl; + + if (q<=p) + if (L2_norm > 1e-10*field_l2_norm) + deallog << "Projection failed with relative error " + << L2_norm / field_l2_norm + << std::endl; + } +} + + + +// check the given element of polynomial order p. the last parameter, if +// given, denotes a gap in convergence order; for example, the Nedelec element +// of polynomial degree p has normal components of degree p-1 and therefore +// can only represent polynomials of degree p-1 exactly. the gap is then 1. +template +void test_no_hanging_nodes (const FiniteElement &fe, + const unsigned int p) +{ + parallel::distributed::Triangulation triangulation(MPI_COMM_WORLD); + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (3); + + do_project (triangulation, fe, p); +} + + + +// same test as above, but this time with a mesh that has hanging nodes +template +void test_with_hanging_nodes (const FiniteElement &fe, + const unsigned int p) +{ + parallel::distributed::Triangulation triangulation(MPI_COMM_WORLD); + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (2); + if (triangulation.begin_active()->is_locally_owned()) + triangulation.begin_active()->set_refine_flag(); + triangulation.execute_coarsening_and_refinement (); + triangulation.refine_global (1); + + do_project (triangulation, fe, p); +} + diff --git a/tests/numerics/project_parallel_qp_q_01.cc b/tests/numerics/project_parallel_qp_q_01.cc new file mode 100644 index 0000000000..06027d6824 --- /dev/null +++ b/tests/numerics/project_parallel_qp_q_01.cc @@ -0,0 +1,45 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +// check VectorTools::project_parallel() of quadrature data for +// FE_Q on a mesh without hanging nodes. + +#include "project_parallel_qp_common.h" + +namespace LA +{ + using namespace ::LinearAlgebraTrilinos; +} + +template +void test () +{ + for (unsigned int p=1; p<7-dim; ++p) + test_no_hanging_nodes (FE_Q(p), p); +} + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + test<2>(); + test<3>(); +} + diff --git a/tests/numerics/project_parallel_qp_q_01.mpirun=3.with_cxx11=on.with_trilinos=true.output b/tests/numerics/project_parallel_qp_q_01.mpirun=3.with_cxx11=on.with_trilinos=true.output new file mode 100644 index 0000000000..d338fe1f2b --- /dev/null +++ b/tests/numerics/project_parallel_qp_q_01.mpirun=3.with_cxx11=on.with_trilinos=true.output @@ -0,0 +1,38 @@ + +DEAL::n_cells=64 +DEAL::n_dofs=81 +DEAL::FE_Q<2>(1), P_0, rel. error=0 +DEAL::FE_Q<2>(1), P_1, rel. error=0 +DEAL::n_cells=64 +DEAL::n_dofs=289 +DEAL::FE_Q<2>(2), P_0, rel. error=0 +DEAL::FE_Q<2>(2), P_1, rel. error=0 +DEAL::FE_Q<2>(2), P_2, rel. error=0 +DEAL::n_cells=64 +DEAL::n_dofs=625 +DEAL::FE_Q<2>(3), P_0, rel. error=0 +DEAL::FE_Q<2>(3), P_1, rel. error=0 +DEAL::FE_Q<2>(3), P_2, rel. error=0 +DEAL::FE_Q<2>(3), P_3, rel. error=0 +DEAL::n_cells=64 +DEAL::n_dofs=1089 +DEAL::FE_Q<2>(4), P_0, rel. error=0 +DEAL::FE_Q<2>(4), P_1, rel. error=0 +DEAL::FE_Q<2>(4), P_2, rel. error=0 +DEAL::FE_Q<2>(4), P_3, rel. error=0 +DEAL::FE_Q<2>(4), P_4, rel. error=0 +DEAL::n_cells=512 +DEAL::n_dofs=729 +DEAL::FE_Q<3>(1), P_0, rel. error=0 +DEAL::FE_Q<3>(1), P_1, rel. error=0 +DEAL::n_cells=512 +DEAL::n_dofs=4913 +DEAL::FE_Q<3>(2), P_0, rel. error=0 +DEAL::FE_Q<3>(2), P_1, rel. error=0 +DEAL::FE_Q<3>(2), P_2, rel. error=0 +DEAL::n_cells=512 +DEAL::n_dofs=15625 +DEAL::FE_Q<3>(3), P_0, rel. error=0 +DEAL::FE_Q<3>(3), P_1, rel. error=0 +DEAL::FE_Q<3>(3), P_2, rel. error=0 +DEAL::FE_Q<3>(3), P_3, rel. error=0 diff --git a/tests/numerics/project_parallel_qp_q_02.cc b/tests/numerics/project_parallel_qp_q_02.cc new file mode 100644 index 0000000000..230793a99f --- /dev/null +++ b/tests/numerics/project_parallel_qp_q_02.cc @@ -0,0 +1,45 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +// check VectorTools::project_parallel() of quadrature data for +// FE_Q on a mesh with hanging nodes. + +#include "project_parallel_qp_common.h" + +namespace LA +{ + using namespace ::LinearAlgebraTrilinos; +} + +template +void test () +{ + for (unsigned int p=1; p<7-dim; ++p) + test_with_hanging_nodes (FE_Q(p), p); +} + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + test<2>(); + test<3>(); +} + diff --git a/tests/numerics/project_parallel_qp_q_02.mpirun=3.with_cxx11=on.with_trilinos=true.output b/tests/numerics/project_parallel_qp_q_02.mpirun=3.with_cxx11=on.with_trilinos=true.output new file mode 100644 index 0000000000..5a79d5559c --- /dev/null +++ b/tests/numerics/project_parallel_qp_q_02.mpirun=3.with_cxx11=on.with_trilinos=true.output @@ -0,0 +1,38 @@ + +DEAL::n_cells=76 +DEAL::n_dofs=97 +DEAL::FE_Q<2>(1), P_0, rel. error=0 +DEAL::FE_Q<2>(1), P_1, rel. error=0 +DEAL::n_cells=76 +DEAL::n_dofs=349 +DEAL::FE_Q<2>(2), P_0, rel. error=0 +DEAL::FE_Q<2>(2), P_1, rel. error=0 +DEAL::FE_Q<2>(2), P_2, rel. error=0 +DEAL::n_cells=76 +DEAL::n_dofs=753 +DEAL::FE_Q<2>(3), P_0, rel. error=0 +DEAL::FE_Q<2>(3), P_1, rel. error=0 +DEAL::FE_Q<2>(3), P_2, rel. error=0 +DEAL::FE_Q<2>(3), P_3, rel. error=0 +DEAL::n_cells=76 +DEAL::n_dofs=1309 +DEAL::FE_Q<2>(4), P_0, rel. error=0 +DEAL::FE_Q<2>(4), P_1, rel. error=0 +DEAL::FE_Q<2>(4), P_2, rel. error=0 +DEAL::FE_Q<2>(4), P_3, rel. error=0 +DEAL::FE_Q<2>(4), P_4, rel. error=0 +DEAL::n_cells=568 +DEAL::n_dofs=827 +DEAL::FE_Q<3>(1), P_0, rel. error=0 +DEAL::FE_Q<3>(1), P_1, rel. error=0 +DEAL::n_cells=568 +DEAL::n_dofs=5559 +DEAL::FE_Q<3>(2), P_0, rel. error=0 +DEAL::FE_Q<3>(2), P_1, rel. error=0 +DEAL::FE_Q<3>(2), P_2, rel. error=0 +DEAL::n_cells=568 +DEAL::n_dofs=17587 +DEAL::FE_Q<3>(3), P_0, rel. error=0 +DEAL::FE_Q<3>(3), P_1, rel. error=0 +DEAL::FE_Q<3>(3), P_2, rel. error=0 +DEAL::FE_Q<3>(3), P_3, rel. error=0 diff --git a/tests/numerics/project_parallel_qpmf_common.h b/tests/numerics/project_parallel_qpmf_common.h new file mode 100644 index 0000000000..cd832d1072 --- /dev/null +++ b/tests/numerics/project_parallel_qpmf_common.h @@ -0,0 +1,236 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +// common framework to check whether an element of polynomial order p can +// represent functions of order q for projection from quadrature points of +// matrix-free approach. + +#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 + +template +class F : public Function +{ +public: + F (const unsigned int q, + const unsigned int n_components) + : + Function(n_components), + q(q) + {} + + virtual double value (const Point &p, + const unsigned int component = 0) const + { + Assert ((component == 0) && (this->n_components == 1), + ExcInternalError()); + double val = 0; + for (unsigned int d=0; d value(const Point> &p_vec) const + { + VectorizedArray res = make_vectorized_array (0.); + 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]; + res[v] = value(p); + } + return res; + } + + +private: + const unsigned int q; +}; + + +template +void do_project (const parallel::distributed::Triangulation &triangulation, + const FiniteElement &fe, + const unsigned int p) +{ + AssertThrow(fe.n_components() ==1, + ExcNotImplemented()); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs (fe); + + deallog << "n_cells=" << triangulation.n_global_active_cells() << std::endl; + deallog << "n_dofs=" << dof_handler.n_dofs() << std::endl; + + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs (dof_handler, + locally_relevant_dofs); + + ConstraintMatrix constraints; + constraints.reinit (locally_relevant_dofs); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + QGauss<1> quadrature_formula_1d(n_q_points_1d); + + Table<2, VectorizedArray > qp_data; + + typename MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = MatrixFree::AdditionalData::partition_color; + additional_data.mpi_communicator = MPI_COMM_WORLD; + additional_data.mapping_update_flags = update_values | update_JxW_values | update_quadrature_points; + MatrixFree data; + data.reinit (dof_handler, constraints, quadrature_formula_1d, additional_data); + + for (unsigned int q=0; q<=p; ++q) + { + // setup quadrature data: + F function(q, fe.n_components()); + + // initialize a quadrature data + { + FEEvaluation fe_eval(data); + const unsigned int n_cells = data.n_macro_cells(); + const unsigned int n_q_points = fe_eval.n_q_points; + + qp_data.reinit(n_cells, n_q_points); + for (unsigned int cell=0; cell field; + data.initialize_dof_vector(field); + VectorTools::project + (data, + constraints, + n_q_points_1d, + [=] (const unsigned int cell, const unsigned int q) -> VectorizedArray { return qp_data(cell,q); }, + field); + + field.update_ghost_values(); + + const double field_l2_norm = field.l2_norm(); + + // L2 norm of the difference between FE field and the function + double L2_norm = 0.; + { + QGauss quadrature_formula_error(std::max(p,q)+1); + FEValues fe_values (fe, quadrature_formula_error, + update_values | update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula_error.size(); + std::vector values (n_q_points); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + fe_values.reinit (cell); + fe_values.get_function_values(field, + values); + + for (unsigned int q_point=0; q_point(values[q_point] - function.value(fe_values.quadrature_point(q_point))) * + fe_values.JxW(q_point); + + } + } + + L2_norm = Utilities::MPI::sum(L2_norm, MPI_COMM_WORLD); + L2_norm = std::sqrt(L2_norm); + + deallog << fe.get_name() << ", P_" << q + << ", rel. error=" << L2_norm / field_l2_norm + << std::endl; + + if (q<=p) + if (L2_norm > 1e-10*field_l2_norm) + deallog << "Projection failed with relative error " + << L2_norm / field_l2_norm + << std::endl; + } +} + + + +// check the given element of polynomial order p. the last parameter, if +// given, denotes a gap in convergence order; for example, the Nedelec element +// of polynomial degree p has normal components of degree p-1 and therefore +// can only represent polynomials of degree p-1 exactly. the gap is then 1. +template +void test_no_hanging_nodes (const FiniteElement &fe, + const unsigned int p) +{ + parallel::distributed::Triangulation triangulation(MPI_COMM_WORLD); + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (3); + + do_project (triangulation, fe, p); +} + + + +// same test as above, but this time with a mesh that has hanging nodes +template +void test_with_hanging_nodes (const FiniteElement &fe, + const unsigned int p) +{ + parallel::distributed::Triangulation triangulation(MPI_COMM_WORLD); + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (2); + if (triangulation.begin_active()->is_locally_owned()) + triangulation.begin_active()->set_refine_flag(); + triangulation.execute_coarsening_and_refinement (); + triangulation.refine_global (1); + + do_project (triangulation, fe, p); +} + diff --git a/tests/numerics/project_parallel_qpmf_q_01.cc b/tests/numerics/project_parallel_qpmf_q_01.cc new file mode 100644 index 0000000000..3979347037 --- /dev/null +++ b/tests/numerics/project_parallel_qpmf_q_01.cc @@ -0,0 +1,43 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +// check VectorTools::project_parallel() for matrix-free quadrature data for +// FE_Q on a mesh without hanging nodes. + + +#include "project_parallel_qpmf_common.h" + +template +void test () +{ + test_no_hanging_nodes<1, 2,dim> (FE_Q(1), 1); + test_no_hanging_nodes<2, 3,dim> (FE_Q(2), 2); + test_no_hanging_nodes<3, 4,dim> (FE_Q(3), 3); +} + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, + numbers::invalid_unsigned_int); + + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + test<2>(); + test<3>(); +} + diff --git a/tests/numerics/project_parallel_qpmf_q_01.mpirun=3.with_cxx11=on.output b/tests/numerics/project_parallel_qpmf_q_01.mpirun=3.with_cxx11=on.output new file mode 100644 index 0000000000..ff5f4b5435 --- /dev/null +++ b/tests/numerics/project_parallel_qpmf_q_01.mpirun=3.with_cxx11=on.output @@ -0,0 +1,31 @@ + +DEAL::n_cells=64 +DEAL::n_dofs=81 +DEAL::FE_Q<2>(1), P_0, rel. error=0 +DEAL::FE_Q<2>(1), P_1, rel. error=0 +DEAL::n_cells=64 +DEAL::n_dofs=289 +DEAL::FE_Q<2>(2), P_0, rel. error=0 +DEAL::FE_Q<2>(2), P_1, rel. error=0 +DEAL::FE_Q<2>(2), P_2, rel. error=0 +DEAL::n_cells=64 +DEAL::n_dofs=625 +DEAL::FE_Q<2>(3), P_0, rel. error=0 +DEAL::FE_Q<2>(3), P_1, rel. error=0 +DEAL::FE_Q<2>(3), P_2, rel. error=0 +DEAL::FE_Q<2>(3), P_3, rel. error=0 +DEAL::n_cells=512 +DEAL::n_dofs=729 +DEAL::FE_Q<3>(1), P_0, rel. error=0 +DEAL::FE_Q<3>(1), P_1, rel. error=0 +DEAL::n_cells=512 +DEAL::n_dofs=4913 +DEAL::FE_Q<3>(2), P_0, rel. error=0 +DEAL::FE_Q<3>(2), P_1, rel. error=0 +DEAL::FE_Q<3>(2), P_2, rel. error=0 +DEAL::n_cells=512 +DEAL::n_dofs=15625 +DEAL::FE_Q<3>(3), P_0, rel. error=0 +DEAL::FE_Q<3>(3), P_1, rel. error=0 +DEAL::FE_Q<3>(3), P_2, rel. error=0 +DEAL::FE_Q<3>(3), P_3, rel. error=0 diff --git a/tests/numerics/project_parallel_qpmf_q_02.cc b/tests/numerics/project_parallel_qpmf_q_02.cc new file mode 100644 index 0000000000..577098070c --- /dev/null +++ b/tests/numerics/project_parallel_qpmf_q_02.cc @@ -0,0 +1,42 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +// check VectorTools::project_parallel() for matrix-free quadrature data for +// FE_Q on a mesh with hanging nodes. + +#include "project_parallel_qpmf_common.h" + +template +void test () +{ + test_with_hanging_nodes<1, 2,dim> (FE_Q(1), 1); + test_with_hanging_nodes<2, 3,dim> (FE_Q(2), 2); + test_with_hanging_nodes<3, 4,dim> (FE_Q(3), 3); +} + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, + numbers::invalid_unsigned_int); + + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + test<2>(); + test<3>(); +} + diff --git a/tests/numerics/project_parallel_qpmf_q_02.mpirun=3.with_cxx11=on.output b/tests/numerics/project_parallel_qpmf_q_02.mpirun=3.with_cxx11=on.output new file mode 100644 index 0000000000..78c248a3be --- /dev/null +++ b/tests/numerics/project_parallel_qpmf_q_02.mpirun=3.with_cxx11=on.output @@ -0,0 +1,31 @@ + +DEAL::n_cells=76 +DEAL::n_dofs=97 +DEAL::FE_Q<2>(1), P_0, rel. error=0 +DEAL::FE_Q<2>(1), P_1, rel. error=0 +DEAL::n_cells=76 +DEAL::n_dofs=349 +DEAL::FE_Q<2>(2), P_0, rel. error=0 +DEAL::FE_Q<2>(2), P_1, rel. error=0 +DEAL::FE_Q<2>(2), P_2, rel. error=0 +DEAL::n_cells=76 +DEAL::n_dofs=753 +DEAL::FE_Q<2>(3), P_0, rel. error=0 +DEAL::FE_Q<2>(3), P_1, rel. error=0 +DEAL::FE_Q<2>(3), P_2, rel. error=0 +DEAL::FE_Q<2>(3), P_3, rel. error=0 +DEAL::n_cells=568 +DEAL::n_dofs=827 +DEAL::FE_Q<3>(1), P_0, rel. error=0 +DEAL::FE_Q<3>(1), P_1, rel. error=0 +DEAL::n_cells=568 +DEAL::n_dofs=5559 +DEAL::FE_Q<3>(2), P_0, rel. error=0 +DEAL::FE_Q<3>(2), P_1, rel. error=0 +DEAL::FE_Q<3>(2), P_2, rel. error=0 +DEAL::n_cells=568 +DEAL::n_dofs=17587 +DEAL::FE_Q<3>(3), P_0, rel. error=0 +DEAL::FE_Q<3>(3), P_1, rel. error=0 +DEAL::FE_Q<3>(3), P_2, rel. error=0 +DEAL::FE_Q<3>(3), P_3, rel. error=0