Permalink
Browse files

Merge pull request #3773 from kronbichler/mg_transfer_mf_legendre

Allow for FE_DGQLegendre in MGTransferMatrixFree
  • Loading branch information...
2 parents 2fc43b3 + dc37703 commit 7ed5c238d8d12f24a7c23c1caa3ccab00e45ec5d @kronbichler kronbichler committed on GitHub Jan 12, 2017
@@ -20,7 +20,6 @@
#include <deal.II/base/mg_level_object.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/lac/la_parallel_vector.h>
-#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
DEAL_II_NAMESPACE_OPEN
@@ -88,10 +87,17 @@ namespace internal
unsigned int n_child_cell_dofs;
/**
+ * Holds the numbering between the numbering of degrees of freedom in
+ * the finite element and the lexicographic numbering needed for the
+ * tensor product application.
+ */
+ std::vector<unsigned int> lexicographic_numbering;
+
+ /**
* Holds the one-dimensional embedding (prolongation) matrix from mother
- * element to the children.
+ * element to all the children.
*/
- internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info;
+ AlignedVector<VectorizedArray<Number> > prolongation_matrix_1d;
};
@@ -23,7 +23,6 @@
#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/base/mg_level_object.h>
#include <deal.II/multigrid/mg_transfer.h>
-#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/dofs/dof_handler.h>
@@ -188,9 +187,9 @@ class MGTransferMatrixFree : public MGLevelGlobalTransfer<LinearAlgebra::distrib
/**
* Holds the one-dimensional embedding (prolongation) matrix from mother
- * element to the children.
+ * element to all the children.
*/
- internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info;
+ AlignedVector<VectorizedArray<Number> > prolongation_matrix_1d;
/**
* Holds the temporary values for the tensor evaluation
@@ -17,6 +17,7 @@
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_tools.h>
+#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/multigrid/mg_transfer_internal.h>
DEAL_II_NAMESPACE_OPEN
@@ -414,35 +415,30 @@ namespace internal
renumbering[fe.dofs_per_cell-fe.dofs_per_vertex] = fe.dofs_per_vertex;
}
- // step 1.3: create a 1D quadrature formula from the finite element that
- // collects the support points of the basis functions on the two children.
+ // step 1.3: create a dummy 1D quadrature formula to extract the
+ // lexicographic numbering for the elements
std::vector<Point<1> > basic_support_points = fe.get_unit_support_points();
Assert(fe.dofs_per_vertex == 0 || fe.dofs_per_vertex == 1,
ExcNotImplemented());
- std::vector<Point<1> > points_refined(fe.dofs_per_vertex > 0 ?
+ const unsigned int shift = fe.dofs_per_cell - fe.dofs_per_vertex;
+ const unsigned int n_child_dofs_1d = (fe.dofs_per_vertex > 0 ?
(2 * fe.dofs_per_cell - 1) :
(2 * fe.dofs_per_cell));
- const unsigned int shift = fe.dofs_per_cell - fe.dofs_per_vertex;
- for (unsigned int c=0; c<GeometryInfo<1>::max_children_per_cell; ++c)
- for (unsigned int j=0; j<basic_support_points.size(); ++j)
- points_refined[shift*c+j][0] =
- c*0.5 + 0.5 * basic_support_points[renumbering[j]][0];
-
- const unsigned int n_child_dofs_1d = points_refined.size();
elem_info.n_child_cell_dofs = elem_info.n_components*Utilities::fixed_power<dim>(n_child_dofs_1d);
+ const Quadrature<1> dummy_quadrature(std::vector<Point<1> >(1, Point<1>()));
+ internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info;
+ shape_info.reinit(dummy_quadrature, mg_dof.get_fe(), 0);
+ elem_info.lexicographic_numbering = shape_info.lexicographic_numbering;
- // step 1.4: evaluate the polynomials and store the data in ShapeInfo
- const Quadrature<1> quadrature(points_refined);
- elem_info.shape_info.reinit(quadrature, mg_dof.get_fe(), 0);
+ // step 1.4: get the 1d prolongation matrix and combine from both children
+ elem_info.prolongation_matrix_1d.resize(fe.dofs_per_cell*n_child_dofs_1d);
for (unsigned int c=0; c<GeometryInfo<1>::max_children_per_cell; ++c)
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
- Assert(std::abs(elem_info.shape_info.shape_values[i*n_child_dofs_1d+j+c*shift][0] -
- fe.get_prolongation_matrix(c)(renumbering[j],renumbering[i]))
- < std::max(2.*(double)std::numeric_limits<Number>::epsilon(),1e-12),
- ExcInternalError());
+ elem_info.prolongation_matrix_1d[i*n_child_dofs_1d+j+c*shift] =
+ fe.get_prolongation_matrix(c)(renumbering[j],renumbering[i]);
}
@@ -567,7 +563,7 @@ namespace internal
ghosted_level_dofs.push_back(local_dof_indices[i]);
add_child_indices<dim>(c, fe.dofs_per_cell - fe.dofs_per_vertex,
- fe.degree, elem_info.shape_info.lexicographic_numbering,
+ fe.degree, elem_info.lexicographic_numbering,
local_dof_indices,
&next_indices[start_index]);
@@ -598,7 +594,7 @@ namespace internal
if (mg_constrained_dofs != 0)
for (unsigned int i=0; i<mg_dof.get_fe().dofs_per_cell; ++i)
if (mg_constrained_dofs->is_boundary_index(level,
- local_dof_indices[elem_info.shape_info.lexicographic_numbering[i]]))
+ local_dof_indices[elem_info.lexicographic_numbering[i]]))
dirichlet_indices[level][child_index].push_back(i);
}
}
@@ -626,14 +622,14 @@ namespace internal
global_level_dof_indices_l0.resize(start_index+elem_info.n_child_cell_dofs,
numbers::invalid_dof_index);
add_child_indices<dim>(0, fe.dofs_per_cell - fe.dofs_per_vertex,
- fe.degree, elem_info.shape_info.lexicographic_numbering,
+ fe.degree, elem_info.lexicographic_numbering,
local_dof_indices,
&global_level_dof_indices_l0[start_index]);
dirichlet_indices[0].push_back(std::vector<unsigned short>());
if (mg_constrained_dofs != 0)
for (unsigned int i=0; i<mg_dof.get_fe().dofs_per_cell; ++i)
- if (mg_constrained_dofs->is_boundary_index(0, local_dof_indices[elem_info.shape_info.lexicographic_numbering[i]]))
+ if (mg_constrained_dofs->is_boundary_index(0, local_dof_indices[elem_info.lexicographic_numbering[i]]))
dirichlet_indices[0].back().push_back(i);
}
}
@@ -27,7 +27,6 @@
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/multigrid/mg_transfer_internal.h>
-#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <algorithm>
@@ -85,7 +84,7 @@ void MGTransferMatrixFree<dim,Number>::clear ()
level_dof_indices.clear();
parent_child_connect.clear();
n_owned_level_cells.clear();
- shape_info = internal::MatrixFreeFunctions::ShapeInfo<Number>();
+ prolongation_matrix_1d.clear();
evaluation_data.clear();
weights_on_refined.clear();
}
@@ -116,7 +115,7 @@ void MGTransferMatrixFree<dim,Number>::build
element_is_continuous = elem_info.element_is_continuous;
n_components = elem_info.n_components;
n_child_cell_dofs = elem_info.n_child_cell_dofs;
- shape_info = elem_info.shape_info;
+ prolongation_matrix_1d = elem_info.prolongation_matrix_1d;
// reshuffle into aligned vector of vectorized arrays
@@ -387,16 +386,14 @@ ::do_prolongate_add (const unsigned int to_level,
}
// perform tensorized operation
- Assert(shape_info.element_type ==
- internal::MatrixFreeFunctions::tensor_symmetric, ExcNotImplemented());
if (element_is_continuous)
{
- AssertDimension(shape_info.shape_val_evenodd.size(),
- (degree+1)*(degree+1));
- typedef internal::EvaluatorTensorProduct<internal::evaluate_evenodd,dim,degree,2*degree+1,VectorizedArray<Number> > Evaluator;
- Evaluator evaluator(shape_info.shape_val_evenodd,
- shape_info.shape_val_evenodd,
- shape_info.shape_val_evenodd);
+ AssertDimension(prolongation_matrix_1d.size(),
+ (2*degree+1)*(degree+1));
+ typedef internal::EvaluatorTensorProduct<internal::evaluate_general,dim,degree,2*degree+1,VectorizedArray<Number> > Evaluator;
+ Evaluator evaluator(prolongation_matrix_1d,
+ prolongation_matrix_1d,
+ prolongation_matrix_1d);
perform_tensorized_op<dim,Evaluator,Number,true>(evaluator,
n_child_cell_dofs,
n_components,
@@ -407,12 +404,12 @@ ::do_prolongate_add (const unsigned int to_level,
}
else
{
- AssertDimension(shape_info.shape_val_evenodd.size(),
- (degree+1)*(degree+1));
- typedef internal::EvaluatorTensorProduct<internal::evaluate_evenodd,dim,degree,2*degree+2,VectorizedArray<Number> > Evaluator;
- Evaluator evaluator(shape_info.shape_val_evenodd,
- shape_info.shape_val_evenodd,
- shape_info.shape_val_evenodd);
+ AssertDimension(prolongation_matrix_1d.size(),
+ 2*(degree+1)*(degree+1));
+ typedef internal::EvaluatorTensorProduct<internal::evaluate_general,dim,degree,2*(degree+1),VectorizedArray<Number> > Evaluator;
+ Evaluator evaluator(prolongation_matrix_1d,
+ prolongation_matrix_1d,
+ prolongation_matrix_1d);
perform_tensorized_op<dim,Evaluator,Number,true>(evaluator,
n_child_cell_dofs,
n_components,
@@ -464,16 +461,14 @@ ::do_restrict_add (const unsigned int from_level,
}
// perform tensorized operation
- Assert(shape_info.element_type ==
- internal::MatrixFreeFunctions::tensor_symmetric, ExcNotImplemented());
if (element_is_continuous)
{
- AssertDimension(shape_info.shape_val_evenodd.size(),
- (degree+1)*(degree+1));
- typedef internal::EvaluatorTensorProduct<internal::evaluate_evenodd,dim,degree,2*degree+1,VectorizedArray<Number> > Evaluator;
- Evaluator evaluator(shape_info.shape_val_evenodd,
- shape_info.shape_val_evenodd,
- shape_info.shape_val_evenodd);
+ AssertDimension(prolongation_matrix_1d.size(),
+ (2*degree+1)*(degree+1));
+ typedef internal::EvaluatorTensorProduct<internal::evaluate_general,dim,degree,2*degree+1,VectorizedArray<Number> > Evaluator;
+ Evaluator evaluator(prolongation_matrix_1d,
+ prolongation_matrix_1d,
+ prolongation_matrix_1d);
weight_dofs_on_child<dim,degree,Number>(&weights_on_refined[from_level-1][(cell/vec_size)*three_to_dim],
n_components,
&evaluation_data[0]);
@@ -484,12 +479,12 @@ ::do_restrict_add (const unsigned int from_level,
}
else
{
- AssertDimension(shape_info.shape_val_evenodd.size(),
- (degree+1)*(degree+1));
- typedef internal::EvaluatorTensorProduct<internal::evaluate_evenodd,dim,degree,2*degree+2,VectorizedArray<Number> > Evaluator;
- Evaluator evaluator(shape_info.shape_val_evenodd,
- shape_info.shape_val_evenodd,
- shape_info.shape_val_evenodd);
+ AssertDimension(prolongation_matrix_1d.size(),
+ 2*(degree+1)*(degree+1));
+ typedef internal::EvaluatorTensorProduct<internal::evaluate_general,dim,degree,2*(degree+1),VectorizedArray<Number> > Evaluator;
+ Evaluator evaluator(prolongation_matrix_1d,
+ prolongation_matrix_1d,
+ prolongation_matrix_1d);
perform_tensorized_op<dim,Evaluator,Number,false>(evaluator,
n_child_cell_dofs,
n_components,
@@ -534,7 +529,7 @@ MGTransferMatrixFree<dim,Number>::memory_consumption() const
memory += MemoryConsumption::memory_consumption(level_dof_indices);
memory += MemoryConsumption::memory_consumption(parent_child_connect);
memory += MemoryConsumption::memory_consumption(n_owned_level_cells);
- memory += shape_info.memory_consumption();
+ memory += MemoryConsumption::memory_consumption(prolongation_matrix_1d);
memory += MemoryConsumption::memory_consumption(evaluation_data);
memory += MemoryConsumption::memory_consumption(weights_on_refined);
memory += MemoryConsumption::memory_consumption(dirichlet_indices);
Oops, something went wrong.

0 comments on commit 7ed5c23

Please sign in to comment.