@DEAL_II_EXPAND_PETSC_MPI_BLOCKVECTOR_REAL@;
}
+REAL_NONBLOCK_VECTORS := { Vector<double>;
+ Vector<float> ;
+
+ LinearAlgebra::Vector<double>;
+ LinearAlgebra::Vector<float> ;
+
+ LinearAlgebra::distributed::Vector<double>;
+ LinearAlgebra::distributed::Vector<float> ;
+
+ @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@;
<h3>Specific improvements</h3>
<ol>
+ <li> New: Add VectorTools::project() to do L2 projection
+ of scalar-valued quadrature point data in parallel.
+ <br>
+ (Denis Davydov, 2016/10/28)
+ <li>
<li> Fixed: Increased precision of timesteps in DataOutInterface::write_pvd_record().
<br>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/point.h>
+#include <deal.II/base/std_cxx11/function.h>
#include <deal.II/dofs/function_map.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/hp/dof_handler.h>
template <int dim, typename Number> class Function;
template <int dim> class Quadrature;
template <int dim> class QGauss;
+template <int dim, typename number> class MatrixFree;
template <typename number> class Vector;
template <typename number> class FullMatrix;
Quadrature<dim-1>(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<dim>::active_cell_iterator & cell, const unsigned int q) -> double
+ * { return qp_data.get_data(cell)[q]->density; },
+ * field);
+ * @endcode
+ * where <code>qp_data</code> is a CellDataStorage object, which stores
+ * quadrature point data.
+ */
+ template <int dim, typename VectorType, int spacedim>
+ void project (const Mapping<dim, spacedim> &mapping,
+ const DoFHandler<dim,spacedim> &dof,
+ const ConstraintMatrix &constraints,
+ const Quadrature<dim> &quadrature,
+ const std_cxx11::function< typename VectorType::value_type (const typename DoFHandler<dim, spacedim>::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<double>
+ * { return qp_data(cell,q); },
+ * field);
+ * @endcode
+ * where <code>qp_data</code> is a an object of type Table<2, VectorizedArray<double> >,
+ * which stores quadrature point data.
+ */
+ template <int dim, typename VectorType>
+ void project (const MatrixFree<dim,typename VectorType::value_type> &data,
+ const ConstraintMatrix &constraints,
+ const unsigned int n_q_points_1d,
+ const std_cxx11::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> func,
+ VectorType &vec_result);
+
+ /**
+ * Same as above but for <code>n_q_points_1d = matrix_free.get_dof_handler().get_fe().degree+1</code>.
+ */
+ template <int dim, typename VectorType>
+ void project (const MatrixFree<dim,typename VectorType::value_type> &data,
+ const ConstraintMatrix &constraints,
+ const std_cxx11::function< VectorizedArray<typename VectorType::value_type> (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
#include <deal.II/base/function.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/qprojector.h>
+#include <deal.II/distributed/tria_base.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/la_vector.h>
#include <deal.II/lac/vector_memory.h>
#include <deal.II/lac/filtered_matrix.h>
#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/operators.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_boundary.h>
for (unsigned int i=0; i<vec.size(); ++i)
vec_result(i) = vec(i);
}
+
+
+
+ template <int dim, typename VectorType, int spacedim, int fe_degree>
+ void project_parallel (const Mapping<dim, spacedim> &mapping,
+ const DoFHandler<dim,spacedim> &dof,
+ const ConstraintMatrix &constraints,
+ const Quadrature<dim> &quadrature,
+ const std_cxx11::function< typename VectorType::value_type (const typename DoFHandler<dim, spacedim>::active_cell_iterator &, const unsigned int)> func,
+ VectorType &vec_result)
+ {
+ const parallel::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::Triangulation<dim,spacedim>*>(&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<dim,Number>::AdditionalData additional_data;
+ additional_data.tasks_parallel_scheme =
+ MatrixFree<dim,Number>::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<dim, Number> matrix_free;
+ matrix_free.reinit (mapping, dof, constraints,
+ QGauss<1>(fe_degree+1), additional_data);
+ typedef MatrixFreeOperators::MassOperator<dim, fe_degree, 1, Number> MatrixType;
+ MatrixType mass_matrix;
+ mass_matrix.initialize(matrix_free);
+ mass_matrix.compute_diagonal();
+
+ typedef LinearAlgebra::distributed::Vector<Number> 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<dim> 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<Number> cell_rhs (dofs_per_cell);
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+ typename DoFHandler<dim, spacedim>::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_point<n_q_points; ++q_point)
+ {
+ const double val_q = func(cell, q_point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+ val_q *
+ fe_values.JxW(q_point));
+ }
+
+ cell->get_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<LinearAlgebra::distributed::Vector<Number> > cg(control);
+ typename PreconditionChebyshev<MatrixType, LocalVectorType>::AdditionalData data;
+ data.matrix_diagonal_inverse = mass_matrix.get_matrix_diagonal_inverse();
+ PreconditionChebyshev<MatrixType, LocalVectorType> 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 <int dim, typename VectorType, int fe_degree, int n_q_points_1d>
+ void project_parallel (const MatrixFree<dim,typename VectorType::value_type> &matrix_free,
+ const ConstraintMatrix &constraints,
+ const std_cxx11::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> func,
+ VectorType &vec_result)
+ {
+ const DoFHandler<dim> &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<dim, fe_degree, 1, Number> MatrixType;
+ MatrixType mass_matrix;
+ mass_matrix.initialize(matrix_free);
+ mass_matrix.compute_diagonal();
+
+ typedef LinearAlgebra::distributed::Vector<Number> 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<dim,fe_degree,n_q_points_1d,1,Number> 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<n_cells; ++cell)
+ {
+ fe_eval.reinit(cell);
+ for (unsigned int q= 0; q < n_q_points; ++q)
+ fe_eval.submit_value(func(cell,q), q);
+
+ fe_eval.integrate (true,false);
+ fe_eval.distribute_local_to_global (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<LinearAlgebra::distributed::Vector<Number> > cg(control);
+ typename PreconditionChebyshev<MatrixType, LocalVectorType>::AdditionalData data;
+ data.matrix_diagonal_inverse = mass_matrix.get_matrix_diagonal_inverse();
+ PreconditionChebyshev<MatrixType, LocalVectorType> 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 <int dim, typename VectorType, int spacedim>
+ void project (const Mapping<dim, spacedim> &mapping,
+ const DoFHandler<dim,spacedim> &dof,
+ const ConstraintMatrix &constraints,
+ const Quadrature<dim> &quadrature,
+ const std_cxx11::function< typename VectorType::value_type (const typename DoFHandler<dim, spacedim>::active_cell_iterator &, const unsigned int)> func,
+ VectorType &vec_result)
+ {
+ switch (dof.get_fe().degree)
+ {
+ case 1:
+ project_parallel<dim, VectorType,spacedim,1> (mapping,dof,constraints,quadrature,func,vec_result);
+ break;
+ case 2:
+ project_parallel<dim, VectorType,spacedim,2> (mapping,dof,constraints,quadrature,func,vec_result);
+ break;
+ case 3:
+ project_parallel<dim, VectorType,spacedim,3> (mapping,dof,constraints,quadrature,func,vec_result);
+ break;
+ case 4:
+ project_parallel<dim, VectorType,spacedim,4> (mapping,dof,constraints,quadrature,func,vec_result);
+ break;
+ case 5:
+ project_parallel<dim, VectorType,spacedim,5> (mapping,dof,constraints,quadrature,func,vec_result);
+ break;
+ case 6:
+ project_parallel<dim, VectorType,spacedim,6> (mapping,dof,constraints,quadrature,func,vec_result);
+ break;
+ case 7:
+ project_parallel<dim, VectorType,spacedim,7> (mapping,dof,constraints,quadrature,func,vec_result);
+ break;
+ case 8:
+ project_parallel<dim, VectorType,spacedim,8> (mapping,dof,constraints,quadrature,func,vec_result);
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
}
+
+ template <int dim, typename VectorType>
+ void project (const MatrixFree<dim,typename VectorType::value_type> &matrix_free,
+ const ConstraintMatrix &constraints,
+ const unsigned int n_q_points_1d,
+ const std_cxx11::function< VectorizedArray<typename VectorType::value_type> (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<dim, VectorType,1,2> (matrix_free,constraints,func,vec_result);
+ break;
+ case 2:
+ project_parallel<dim, VectorType,2,3> (matrix_free,constraints,func,vec_result);
+ break;
+ case 3:
+ project_parallel<dim, VectorType,3,4> (matrix_free,constraints,func,vec_result);
+ break;
+ case 4:
+ project_parallel<dim, VectorType,4,5> (matrix_free,constraints,func,vec_result);
+ break;
+ case 5:
+ project_parallel<dim, VectorType,5,6> (matrix_free,constraints,func,vec_result);
+ break;
+ case 6:
+ project_parallel<dim, VectorType,6,7> (matrix_free,constraints,func,vec_result);
+ break;
+ case 7:
+ project_parallel<dim, VectorType,7,8> (matrix_free,constraints,func,vec_result);
+ break;
+ case 8:
+ project_parallel<dim, VectorType,8,9> (matrix_free,constraints,func,vec_result);
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ }
+
+
+
+ template <int dim, typename VectorType>
+ void project (const MatrixFree<dim,typename VectorType::value_type> &matrix_free,
+ const ConstraintMatrix &constraints,
+ const std_cxx11::function< VectorizedArray<typename VectorType::value_type> (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 <int dim, typename VectorType, int spacedim>
void project (const Mapping<dim, spacedim> &mapping,
const DoFHandler<dim,spacedim> &dof,
\}
#endif
}
+
+for (VEC: REAL_NONBLOCK_VECTORS; deal_II_dimension : DIMENSIONS)
+{
+ namespace VectorTools \{
+
+ template
+ void project<deal_II_dimension,VEC,deal_II_dimension>(
+ const Mapping<deal_II_dimension, deal_II_dimension> &,
+ const DoFHandler<deal_II_dimension,deal_II_dimension> &,
+ const ConstraintMatrix &,
+ const Quadrature<deal_II_dimension> &,
+ const std_cxx11::function<VEC::value_type (const DoFHandler<deal_II_dimension, deal_II_dimension>::active_cell_iterator &, const unsigned int)>,
+ VEC &);
+
+ template
+ void project<deal_II_dimension, VEC>(
+ const MatrixFree<deal_II_dimension,VEC::value_type> &matrix_free,
+ const ConstraintMatrix &constraints,
+ const std_cxx11::function< VectorizedArray<VEC::value_type> (const unsigned int, const unsigned int)>,
+ VEC &);
+
+ template
+ void project<deal_II_dimension, VEC>(
+ const MatrixFree<deal_II_dimension,VEC::value_type> &matrix_free,
+ const ConstraintMatrix &constraints,
+ const unsigned int,
+ const std_cxx11::function< VectorizedArray<VEC::value_type> (const unsigned int, const unsigned int)>,
+ VEC &);
+
+ \}
+
+}
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.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_tools.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/grid_refinement.h>
+
+#include <deal.II/lac/generic_linear_algebra.h>
+#include <deal.II/base/quadrature_point_data.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+
+#include <fstream>
+#include <vector>
+
+template <int dim>
+class F : public Function<dim>
+{
+public:
+ F (const unsigned int q,
+ const unsigned int n_components)
+ :
+ Function<dim>(n_components),
+ q(q)
+ {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const
+ {
+ Assert ((component == 0) && (this->n_components == 1),
+ ExcInternalError());
+ double val = 0;
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int i=0; i<=q; ++i)
+ val += (d+1)*(i+1)*std::pow (p[d], 1.*i);
+ return val;
+ }
+
+ VectorizedArray<double> value(const Point<dim, VectorizedArray<double>> &p_vec) const
+ {
+ VectorizedArray<double> res = make_vectorized_array (0.);
+ Point<dim> p;
+ for (unsigned int v = 0; v < VectorizedArray<double>::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 <typename VectorType, int dim>
+void do_project (const parallel::distributed::Triangulation<dim> &triangulation,
+ const FiniteElement<dim> &fe,
+ const unsigned int p)
+{
+ AssertThrow(fe.n_components() ==1,
+ ExcNotImplemented());
+ DoFHandler<dim> 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<dim> function(q, fe.n_components());
+ VectorType field_relevant;
+ double field_l2_norm = 0.;
+ {
+ QGauss<dim> quadrature_formula(p+2);
+ CellDataStorage<typename parallel::distributed::Triangulation<dim,dim>::cell_iterator,QData> qp_manager;
+ {
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_quadrature_points);
+
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ typename DoFHandler<dim>::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<std::shared_ptr<QData> > qpd = qp_manager.get_data(cell);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ qpd[q]->density = function.value(fe_values.quadrature_point(q));
+ }
+ }
+
+ VectorType field(dof_handler.locally_owned_dofs(),
+ MPI_COMM_WORLD);
+ VectorTools::project
+ (MappingQGeneric<dim>(1),
+ dof_handler,
+ constraints,
+ quadrature_formula,
+ [=] (const typename DoFHandler<dim>::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<dim> quadrature_formula_error(std::max(p,q)+1);
+ FEValues<dim> 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<double> values (n_q_points);
+
+ typename DoFHandler<dim>::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<n_q_points; ++q_point)
+ L2_norm += Utilities::fixed_power<2>(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 <typename VectorType, int dim>
+void test_no_hanging_nodes (const FiniteElement<dim> &fe,
+ const unsigned int p)
+{
+ parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+ GridGenerator::hyper_cube (triangulation);
+ triangulation.refine_global (3);
+
+ do_project<VectorType> (triangulation, fe, p);
+}
+
+
+
+// same test as above, but this time with a mesh that has hanging nodes
+template <typename VectorType, int dim>
+void test_with_hanging_nodes (const FiniteElement<dim> &fe,
+ const unsigned int p)
+{
+ parallel::distributed::Triangulation<dim> 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<VectorType> (triangulation, fe, p);
+}
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <int dim>
+void test ()
+{
+ for (unsigned int p=1; p<7-dim; ++p)
+ test_no_hanging_nodes<LA::MPI::Vector> (FE_Q<dim>(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>();
+}
+
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <int dim>
+void test ()
+{
+ for (unsigned int p=1; p<7-dim; ++p)
+ test_with_hanging_nodes<LA::MPI::Vector> (FE_Q<dim>(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>();
+}
+
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.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_tools.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/grid_refinement.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+
+#include <fstream>
+#include <vector>
+
+template <int dim>
+class F : public Function<dim>
+{
+public:
+ F (const unsigned int q,
+ const unsigned int n_components)
+ :
+ Function<dim>(n_components),
+ q(q)
+ {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const
+ {
+ Assert ((component == 0) && (this->n_components == 1),
+ ExcInternalError());
+ double val = 0;
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int i=0; i<=q; ++i)
+ val += (d+1)*(i+1)*std::pow (p[d], 1.*i);
+ return val;
+ }
+
+ VectorizedArray<double> value(const Point<dim, VectorizedArray<double>> &p_vec) const
+ {
+ VectorizedArray<double> res = make_vectorized_array (0.);
+ Point<dim> p;
+ for (unsigned int v = 0; v < VectorizedArray<double>::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 <int fe_degree, int n_q_points_1d, int dim>
+void do_project (const parallel::distributed::Triangulation<dim> &triangulation,
+ const FiniteElement<dim> &fe,
+ const unsigned int p)
+{
+ AssertThrow(fe.n_components() ==1,
+ ExcNotImplemented());
+ DoFHandler<dim> 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<double> > qp_data;
+
+ typename MatrixFree<dim,double>::AdditionalData additional_data;
+ additional_data.tasks_parallel_scheme = MatrixFree<dim,double>::AdditionalData::partition_color;
+ additional_data.mpi_communicator = MPI_COMM_WORLD;
+ additional_data.mapping_update_flags = update_values | update_JxW_values | update_quadrature_points;
+ MatrixFree<dim,double> data;
+ data.reinit (dof_handler, constraints, quadrature_formula_1d, additional_data);
+
+ for (unsigned int q=0; q<=p; ++q)
+ {
+ // setup quadrature data:
+ F<dim> function(q, fe.n_components());
+
+ // initialize a quadrature data
+ {
+ FEEvaluation<dim,fe_degree,n_q_points_1d,1,double> 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<n_cells; ++cell)
+ {
+ fe_eval.reinit(cell);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ qp_data(cell,q) = function.value(fe_eval.quadrature_point(q));
+ }
+ }
+
+ LinearAlgebra::distributed::Vector<double> field;
+ data.initialize_dof_vector(field);
+ VectorTools::project
+ (data,
+ constraints,
+ n_q_points_1d,
+ [=] (const unsigned int cell, const unsigned int q) -> VectorizedArray<double> { 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<dim> quadrature_formula_error(std::max(p,q)+1);
+ FEValues<dim> 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<double> values (n_q_points);
+
+ typename DoFHandler<dim>::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<n_q_points; ++q_point)
+ L2_norm += Utilities::fixed_power<2>(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 <int fe_degree, int n_q_points_1d, int dim>
+void test_no_hanging_nodes (const FiniteElement<dim> &fe,
+ const unsigned int p)
+{
+ parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+ GridGenerator::hyper_cube (triangulation);
+ triangulation.refine_global (3);
+
+ do_project<fe_degree, n_q_points_1d, dim> (triangulation, fe, p);
+}
+
+
+
+// same test as above, but this time with a mesh that has hanging nodes
+template <int fe_degree, int n_q_points_1d, int dim>
+void test_with_hanging_nodes (const FiniteElement<dim> &fe,
+ const unsigned int p)
+{
+ parallel::distributed::Triangulation<dim> 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<fe_degree, n_q_points_1d, dim> (triangulation, fe, p);
+}
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <int dim>
+void test ()
+{
+ test_no_hanging_nodes<1, 2,dim> (FE_Q<dim>(1), 1);
+ test_no_hanging_nodes<2, 3,dim> (FE_Q<dim>(2), 2);
+ test_no_hanging_nodes<3, 4,dim> (FE_Q<dim>(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>();
+}
+
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <int dim>
+void test ()
+{
+ test_with_hanging_nodes<1, 2,dim> (FE_Q<dim>(1), 1);
+ test_with_hanging_nodes<2, 3,dim> (FE_Q<dim>(2), 2);
+ test_with_hanging_nodes<3, 4,dim> (FE_Q<dim>(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>();
+}
+
--- /dev/null
+
+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
projects / dealii.git / commitdiff