--- /dev/null
+Extend VectorTools::project() function with MatrixFree quadrature data to
+optionally take the finite element component index.
+<br>
+(Denis Davydov, 2018/02/07)
* @endcode
* where <code>qp_data</code> is a an object of type Table<2, VectorizedArray<double> >,
* which stores quadrature point data.
+ *
+ * @p fe_component allow to additionally specify which component of @p data
+ * to use in case it was constructed with an <code>std::vector<const DoFHandler<dim>*></code>
+ * in will be internally used in constructor of FEEvaluation object.
*/
template <int dim, typename VectorType>
void project (std::shared_ptr<const MatrixFree<dim,typename VectorType::value_type> > data,
const ConstraintMatrix &constraints,
const unsigned int n_q_points_1d,
const std::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> &func,
- VectorType &vec_result);
+ VectorType &vec_result,
+ const unsigned int fe_component = 0);
/**
* Same as above but for <code>n_q_points_1d = matrix_free.get_dof_handler().get_fe().degree+1</code>.
void project (std::shared_ptr<const MatrixFree<dim,typename VectorType::value_type> > data,
const ConstraintMatrix &constraints,
const std::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> &func,
- VectorType &vec_result);
+ VectorType &vec_result,
+ const unsigned int fe_component = 0);
/**
* Compute Dirichlet boundary conditions. This function makes up a map of
void project_parallel (std::shared_ptr<const MatrixFree<dim,typename VectorType::value_type> > matrix_free,
const ConstraintMatrix &constraints,
const std::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> &func,
- VectorType &vec_result)
+ VectorType &vec_result,
+ const unsigned int fe_component)
{
- const DoFHandler<dim,spacedim> &dof = matrix_free->get_dof_handler();
+ const DoFHandler<dim,spacedim> &dof = matrix_free->get_dof_handler(fe_component);
typedef typename VectorType::value_type Number;
Assert (dof.get_fe(0).n_components() == 1,
typedef MatrixFreeOperators::MassOperator<dim, fe_degree, n_q_points_1d, 1, LinearAlgebra::distributed::Vector<Number> > MatrixType;
MatrixType mass_matrix;
- mass_matrix.initialize(matrix_free);
+ mass_matrix.initialize(matrix_free, {{fe_component}});
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);
+ matrix_free->initialize_dof_vector(vec,fe_component);
+ matrix_free->initialize_dof_vector(rhs,fe_component);
+ matrix_free->initialize_dof_vector(inhomogeneities,fe_component);
constraints.distribute(inhomogeneities);
inhomogeneities*=-1.;
// assemble right hand side:
{
- FEEvaluation<dim,fe_degree,n_q_points_1d,1,Number> fe_eval(*matrix_free);
+ FEEvaluation<dim,fe_degree,n_q_points_1d,1,Number> fe_eval(*matrix_free, fe_component);
const unsigned int n_cells = matrix_free->n_macro_cells();
const unsigned int n_q_points = fe_eval.n_q_points;
const ConstraintMatrix &constraints,
const unsigned int n_q_points_1d,
const std::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> &func,
- VectorType &vec_result)
+ VectorType &vec_result,
+ const unsigned int fe_component)
{
- (void) n_q_points_1d;
- const unsigned int fe_degree = matrix_free->get_dof_handler().get_fe().degree;
-
- Assert (fe_degree+1 == n_q_points_1d,
- ExcNotImplemented());
+ const unsigned int fe_degree = matrix_free->get_dof_handler(fe_component).get_fe().degree;
- switch (fe_degree)
- {
- case 1:
- project_parallel<dim,VectorType,dim,1,2> (matrix_free,constraints,func,vec_result);
- break;
- case 2:
- project_parallel<dim,VectorType,dim,2,3> (matrix_free,constraints,func,vec_result);
- break;
- case 3:
- project_parallel<dim,VectorType,dim,3,4> (matrix_free,constraints,func,vec_result);
- break;
- default:
- project_parallel<dim,VectorType,dim,-1,0> (matrix_free,constraints,func,vec_result);
- }
+ if (fe_degree+1 == n_q_points_1d)
+ switch (fe_degree)
+ {
+ case 1:
+ project_parallel<dim,VectorType,dim,1,2> (matrix_free,constraints,func,vec_result,fe_component);
+ break;
+ case 2:
+ project_parallel<dim,VectorType,dim,2,3> (matrix_free,constraints,func,vec_result,fe_component);
+ break;
+ case 3:
+ project_parallel<dim,VectorType,dim,3,4> (matrix_free,constraints,func,vec_result,fe_component);
+ break;
+ default:
+ project_parallel<dim,VectorType,dim,-1,0> (matrix_free,constraints,func,vec_result,fe_component);
+ }
+ else
+ project_parallel<dim,VectorType,dim,-1,0> (matrix_free,constraints,func,vec_result,fe_component);
}
void project (std::shared_ptr<const MatrixFree<dim,typename VectorType::value_type> > matrix_free,
const ConstraintMatrix &constraints,
const std::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> &func,
- VectorType &vec_result)
+ VectorType &vec_result,
+ const unsigned int fe_component)
{
project (matrix_free,
constraints,
- matrix_free->get_dof_handler().get_fe().degree+1,
+ matrix_free->get_dof_handler(fe_component).get_fe().degree+1,
func,
- vec_result);
+ vec_result,
+ fe_component);
}
std::shared_ptr<const MatrixFree<deal_II_dimension,VEC::value_type> > matrix_free,
const ConstraintMatrix &constraints,
const std::function< VectorizedArray<VEC::value_type> (const unsigned int, const unsigned int)> &,
- VEC &);
+ VEC &,
+ const unsigned int);
template
void project<deal_II_dimension, VEC>(
const ConstraintMatrix &constraints,
const unsigned int,
const std::function< VectorizedArray<VEC::value_type> (const unsigned int, const unsigned int)> &,
- VEC &);
+ VEC &,
+ const unsigned int);
\}
#include <fstream>
#include <vector>
+/**
+ * A function of polynomial degree @p q with n_components.
+ */
template <int dim>
class F : public Function<dim>
{
};
+/**
+ * A general function to perform projection of polynomial function of degrees [0,p]
+ * onto the FE space given by @p fe and @p triangulation .
+ */
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)
+ const std::vector<const FiniteElement<dim>*> &fes,
+ const unsigned int p,
+ const unsigned int fe_index = 0)
{
- AssertThrow(fe.n_components() ==1,
+ // we only project scalar value function
+ AssertThrow(fes[fe_index]->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);
+ std::vector<std::shared_ptr<DoFHandler<dim>>> dof_handlers(fes.size());
+ std::vector<IndexSet> locally_relevant_dofs(fes.size());
+ std::vector<ConstraintMatrix> constraints(fes.size());
- ConstraintMatrix constraints;
- constraints.reinit (locally_relevant_dofs);
- DoFTools::make_hanging_node_constraints (dof_handler, constraints);
- constraints.close ();
+ std::vector<const DoFHandler<dim>*> dof_handlers_mf(fes.size());
+ std::vector<const ConstraintMatrix *> constraints_mf(fes.size());
+ for (unsigned int i = 0; i < fes.size(); ++i)
+ {
+ dof_handlers[i] = std::make_shared<DoFHandler<dim>>(triangulation);
+ dof_handlers[i]->distribute_dofs (*fes[i]);
+ deallog << "n_dofs=" << dof_handlers[i]->n_dofs() << std::endl;
+
+ DoFTools::extract_locally_relevant_dofs (*dof_handlers[i],
+ locally_relevant_dofs[i]);
+
+ constraints[i].reinit (locally_relevant_dofs[i]);
+ DoFTools::make_hanging_node_constraints (*dof_handlers[i], constraints[i]);
+ constraints[i].close ();
+
+ dof_handlers_mf[i] = dof_handlers[i].get();
+ constraints_mf[i] = &constraints[i];
+ }
QGauss<1> quadrature_formula_1d(n_q_points_1d);
additional_data.tasks_parallel_scheme = MatrixFree<dim,double>::AdditionalData::partition_color;
additional_data.mapping_update_flags = update_values | update_JxW_values | update_quadrature_points;
std::shared_ptr<MatrixFree<dim,double> > data(new MatrixFree<dim,double> ());
- data->reinit (dof_handler, constraints, quadrature_formula_1d, additional_data);
+ data->reinit (dof_handlers_mf, constraints_mf, quadrature_formula_1d, additional_data);
for (unsigned int q=0; q<=p; ++q)
{
// setup quadrature data:
- F<dim> function(q, fe.n_components());
+ F<dim> function(q, fes[fe_index]->n_components());
// initialize a quadrature data
{
- FEEvaluation<dim,fe_degree,n_q_points_1d,1,double> fe_eval(*data);
+ FEEvaluation<dim,fe_degree,n_q_points_1d,1,double> fe_eval(*data,fe_index);
const unsigned int n_cells = data->n_macro_cells();
const unsigned int n_q_points = fe_eval.n_q_points;
}
LinearAlgebra::distributed::Vector<double> field;
- data->initialize_dof_vector(field);
+ data->initialize_dof_vector(field,fe_index);
VectorTools::project<dim,LinearAlgebra::distributed::Vector<double> >
(data,
- constraints,
+ constraints[fe_index],
n_q_points_1d,
[=] (const unsigned int cell, const unsigned int q) -> VectorizedArray<double> { return qp_data(cell,q); },
- field);
+ field,
+ fe_index);
field.update_ghost_values();
double L2_norm = 0.;
{
QGauss<dim> quadrature_formula_error(std::max(p,q)+1);
- FEValues<dim> fe_values (fe, quadrature_formula_error,
+ FEValues<dim> fe_values (*fes[fe_index], quadrature_formula_error,
update_values | update_quadrature_points | update_JxW_values);
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int dofs_per_cell = fes[fe_index]->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();
+ cell = dof_handlers[fe_index]->begin_active(),
+ endc = dof_handlers[fe_index]->end();
for (; cell!=endc; ++cell)
if (cell->is_locally_owned())
{
L2_norm = Utilities::MPI::sum(L2_norm, MPI_COMM_WORLD);
L2_norm = std::sqrt(L2_norm);
- deallog << fe.get_name() << ", P_" << q
+ deallog << fes[fe_index]->get_name() << ", P_" << q
<< ", rel. error=" << L2_norm / field_l2_norm
<< std::endl;
GridGenerator::hyper_cube (triangulation);
triangulation.refine_global (3);
- do_project<fe_degree, n_q_points_1d, dim> (triangulation, fe, p);
+ do_project<fe_degree, n_q_points_1d, dim> (triangulation, {{&fe}}, p);
}
triangulation.execute_coarsening_and_refinement ();
triangulation.refine_global (1);
- do_project<fe_degree, n_q_points_1d, dim> (triangulation, fe, p);
+ do_project<fe_degree, n_q_points_1d, dim> (triangulation, {{&fe}}, p);
}
+
+// same as above but for multiple fes
+template <int fe_degree, int n_q_points_1d, int dim>
+void test_with_hanging_nodes (const std::vector<const FiniteElement<dim>*> &fes,
+ const unsigned int p,
+ const unsigned int fe_index)
+{
+ 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, fes, p, fe_index);
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 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 and specified fe_component
+
+#include "../tests.h"
+#include "project_parallel_qpmf_common.h"
+
+template <int dim>
+void test ()
+{
+ FE_Q<dim> fe1(1);
+ FE_Q<dim> fe2(2);
+ FE_Q<dim> fe4(4);
+ test_with_hanging_nodes<1, 2,dim> ({{&fe1,&fe2}}, 1, 0); // take first
+ test_with_hanging_nodes<2, 3,dim> ({{&fe1,&fe2}}, 2, 1); // take second
+ test_with_hanging_nodes<2, 6,dim> ({{&fe2,&fe4}}, 2, 0); // take first with non-standard n_q_points_1d
+}
+
+int main(int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv,
+ numbers::invalid_unsigned_int);
+
+ initlog();
+
+ test<2>();
+ test<3>();
+}
--- /dev/null
+
+DEAL::n_cells=76
+DEAL::n_dofs=97
+DEAL::n_dofs=349
+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=97
+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=349
+DEAL::n_dofs=1309
+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=568
+DEAL::n_dofs=827
+DEAL::n_dofs=5559
+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=827
+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=5559
+DEAL::n_dofs=40319
+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
projects / dealii.git / commitdiff