From db9151cfd967c224095ff145a2d5dfae98f34aa7 Mon Sep 17 00:00:00 2001 From: Markus Buerg Date: Fri, 2 Jul 2010 13:44:24 +0000 Subject: [PATCH] Documentation and style changes. git-svn-id: https://svn.dealii.org/trunk@21443 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/fe/fe_q_hierarchical.h | 30 + .../deal.II/source/fe/fe_q_hierarchical.cc | 934 +++++++++++++++++- 2 files changed, 963 insertions(+), 1 deletion(-) diff --git a/deal.II/deal.II/include/fe/fe_q_hierarchical.h b/deal.II/deal.II/include/fe/fe_q_hierarchical.h index 7dbfa15a50..70b547eb15 100644 --- a/deal.II/deal.II/include/fe/fe_q_hierarchical.h +++ b/deal.II/deal.II/include/fe/fe_q_hierarchical.h @@ -322,6 +322,36 @@ class FE_Q_Hierarchical : public FE_Poly,dim> virtual std::vector > hp_vertex_dof_identities (const FiniteElement &fe_other) const; + + /** + * Return the matrix interpolating from a face of one + * element to the face of the neighboring element. The + * size of the matrix is then source.dofs_per_face + * times this->dofs_per_face. + * + * Derived elements will have to implement this function. + * They may only provide interpolation matrices for certain + * source finite elements, for example those from the same + * family. If they don't implement interpolation from a given + * element, then they must throw an exception of type + * FiniteElement::ExcInterpolationNotImplemented. + */ + virtual void get_face_interpolation_matrix (const FiniteElement& source, FullMatrix& matrix) const; + + /** + * Return the matrix interpolating from a face of one element + * to the subface of the neighboring element. The size of + * the matrix is then source.dofs_per_face times + * this->dofs_per_face. + * + * Derived elements will have to implement this function. + * They may only provide interpolation matrices for certain + * source finite elements, for example those from the same + * family. If they don't implement interpolation from a given + * element, then they must throw an exception of type + * ExcInterpolationNotImplemented. + */ + virtual void get_subface_interpolation_matrix (const FiniteElement& source, const unsigned int subface, FullMatrix& matrix) const; /** * Determine an estimate for the diff --git a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc index 03a75d22c2..bbc5a613ef 100644 --- a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc +++ b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc @@ -11,8 +11,13 @@ // //--------------------------------------------------------------------------- +#include +#include #include #include +#include +#include +#include #include #include @@ -126,7 +131,7 @@ template bool FE_Q_Hierarchical::hp_constraints_are_implemented () const { - return false; + return true; } @@ -552,6 +557,933 @@ void FE_Q_Hierarchical<1>::initialize_unit_face_support_points () // no faces in 1d, so nothing to do } + +template <> +void FE_Q_Hierarchical<1>:: +get_face_interpolation_matrix (const FiniteElement<1> &/*x_source_fe*/, + FullMatrix &/*interpolation_matrix*/) const +{ + Assert (false, + FiniteElement<1>:: + ExcInterpolationNotImplemented ()); +} + + +template <> +void +FE_Q_Hierarchical<1>:: +get_subface_interpolation_matrix (const FiniteElement<1> &/*x_source_fe*/, + const unsigned int /*subface*/, + FullMatrix &/*interpolation_matrix*/) const +{ + Assert (false, + FiniteElement<1>:: + ExcInterpolationNotImplemented ()); +} + +#endif + + +#if deal_II_dimension > 1 + +template +void +FE_Q_Hierarchical:: +get_face_interpolation_matrix (const FiniteElement &x_source_fe, + FullMatrix &interpolation_matrix) const +{ + // this is only implemented, if the + // source FE is also a + // Q_Hierarchical element + typedef FE_Q_Hierarchical FEQHierarchical; + typedef FiniteElement FEL; + AssertThrow ((x_source_fe.get_name().find ("FE_Q_Hierarchical<") == 0) + || + (dynamic_cast(&x_source_fe) != 0), + typename FEL:: + ExcInterpolationNotImplemented()); + + Assert (interpolation_matrix.n() == this->dofs_per_face, + ExcDimensionMismatch (interpolation_matrix.n(), + this->dofs_per_face)); + Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face, + ExcDimensionMismatch (interpolation_matrix.m(), + x_source_fe.dofs_per_face)); + + // ok, source is a Q_Hierarchical element, so + // we will be able to do the work + const FE_Q_Hierarchical &source_fe + = dynamic_cast&>(x_source_fe); + + // Make sure, that the element, + // for which the DoFs should be + // constrained is the one with + // the higher polynomial degree. + // Actually the procedure will work + // also if this assertion is not + // satisfied. But the matrices + // produced in that case might + // lead to problems in the + // hp procedures, which use this + // method. + Assert (this->dofs_per_face <= source_fe.dofs_per_face, + typename FEL:: + ExcInterpolationNotImplemented ()); + interpolation_matrix = 0; + + switch (dim) { + case 2: { + for (unsigned int i = 0; i < this->dofs_per_face; ++i) + interpolation_matrix (i, i) = 1; + + break; + } + + case 3: { + for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_face; ++i) + interpolation_matrix (i, i) = 1; + + for (unsigned int i = 0; i < this->degree - 1; ++i) { + for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face; ++j) + interpolation_matrix ( + i + j * (x_source_fe.degree - 1) + GeometryInfo<3>::vertices_per_face, + i + j * (this->degree - 1) + GeometryInfo<3>::vertices_per_face) = 1; + + for (unsigned int j = 0; j < this->degree - 1; ++j) + interpolation_matrix ( + (i + GeometryInfo<3>::lines_per_face) * (x_source_fe.degree - 1) + j + + GeometryInfo<3>::vertices_per_face, + (i + GeometryInfo<3>::lines_per_face) * (this->degree - 1) + j + + GeometryInfo<3>::vertices_per_face) = 1; + } + } + } +} + + + +template +void +FE_Q_Hierarchical:: +get_subface_interpolation_matrix (const FiniteElement &x_source_fe, + const unsigned int subface, + FullMatrix &interpolation_matrix) const +{ + // this is only implemented, if the + // source FE is also a + // Q_Hierarchical element + typedef FE_Q_Hierarchical FEQHierarchical; + typedef FiniteElement FEL; + AssertThrow ((x_source_fe.get_name().find ("FE_Q_Hierarchical<") == 0) + || + (dynamic_cast(&x_source_fe) != 0), + typename FEL:: + ExcInterpolationNotImplemented()); + + Assert (interpolation_matrix.n() == this->dofs_per_face, + ExcDimensionMismatch (interpolation_matrix.n(), + this->dofs_per_face)); + Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face, + ExcDimensionMismatch (interpolation_matrix.m(), + x_source_fe.dofs_per_face)); + + // ok, source is a Q_Hierarchical element, so + // we will be able to do the work + const FE_Q_Hierarchical &source_fe + = dynamic_cast&>(x_source_fe); + + // Make sure, that the element, + // for which the DoFs should be + // constrained is the one with + // the higher polynomial degree. + // Actually the procedure will work + // also if this assertion is not + // satisfied. But the matrices + // produced in that case might + // lead to problems in the + // hp procedures, which use this + // method. + Assert (this->dofs_per_face <= source_fe.dofs_per_face, + typename FEL:: + ExcInterpolationNotImplemented ()); + + switch (dim) { + case 2: { + switch (subface) { + case 0: { + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 0), Point ())) > 1e-14) + interpolation_matrix (0, dof) = this->shape_value (this->face_to_cell_index (dof, 0), Point ()); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 0), Point (0.0, 0.5))) > 1e-14) + interpolation_matrix (1, dof) = this->shape_value (this->face_to_cell_index (dof, 0), Point (0.0, 0.5)); + } + + if (this->degree > 1) { + double weight; + QGauss reference_edge_quadrature (this->degree); + Quadrature edge_quadrature = QProjector::project_to_child (reference_edge_quadrature, subface); + const unsigned int n_edge_points = edge_quadrature.size (); + FullMatrix assembling_matrix (this->degree - 1, dim * n_edge_points); + FullMatrix system_matrix (this->degree - 1, this->degree - 1); + FullMatrix system_matrix_inv (system_matrix.m (), system_matrix.m ()); + Point point; + std::vector > quadrature_points = edge_quadrature.get_points (); + Tensor<1, dim> grad; + Vector assembling_vector (dim * n_edge_points); + Vector solution (system_matrix.m ()); + Vector system_rhs (system_matrix.m ()); + + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (0.0, 2.0 * quadrature_points[q_point] (0)); + weight = std::sqrt (edge_quadrature.weight (q_point)); + + for (unsigned int i = 0; i < system_matrix.m (); ++i) { + grad = weight * this->shape_grad (i + 4, point); + + for (unsigned int d = 0; d < dim; ++d) + assembling_matrix (i, dim * q_point + d) = grad[d]; + } + } + + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_matrix_inv.invert (system_matrix); + + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (0.0, 2.0 * quadrature_points[q_point] (0)); + grad = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_grad (this->face_to_cell_index (dof, 0), Point (0.0, quadrature_points[q_point] (0))) - interpolation_matrix (dof, 0) * source_fe.shape_grad (0, point) - interpolation_matrix (dof, 1) * source_fe.shape_grad (1, point)); + + for (unsigned int d = 0; d < dim; ++d) + assembling_vector (dim * q_point + d) = grad[d]; + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 2, dof) = solution (i); + } + } + + break; + } + + case 1: { + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 0), Point (0.0, 0.5))) > 1e-14) + interpolation_matrix (0, dof) = this->shape_value (this->face_to_cell_index (dof, 0), Point (0.0, 0.5)); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 0), Point (0.0, 1.0))) > 1e-14) + interpolation_matrix (1, dof) = this->shape_value (this->face_to_cell_index (dof, 0), Point (0.0, 1.0)); + } + + if (this->degree > 1) { + double weight; + QGauss reference_edge_quadrature (this->degree); + Quadrature edge_quadrature = QProjector::project_to_child (reference_edge_quadrature, subface); + const unsigned int n_edge_points = edge_quadrature.size (); + FullMatrix assembling_matrix (this->degree - 1, dim * n_edge_points); + FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); + FullMatrix system_matrix_inv (system_matrix.m (), system_matrix.m ()); + Point point; + std::vector > quadrature_points = edge_quadrature.get_points (); + Tensor<1, dim> grad; + Vector assembling_vector (dim * n_edge_points); + Vector solution (system_matrix.m ()); + Vector system_rhs (system_matrix.m ()); + + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (0.0, 2.0 * quadrature_points[q_point] (0) - 1.0); + weight = std::sqrt (edge_quadrature.weight (q_point)); + + for (unsigned int i = 0; i < system_matrix.m (); ++i) { + grad = weight * this->shape_grad (i + 4, point); + + for (unsigned int d = 0; d < dim; ++d) + assembling_matrix (i, dim * q_point + d) = grad[d]; + } + } + + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_matrix_inv.invert (system_matrix); + + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (0.0, 2.0 * quadrature_points[q_point] (0) - 1.0); + grad = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_grad (this->face_to_cell_index (dof, 0), Point (0.0, quadrature_points[q_point] (0))) - interpolation_matrix (0, dof) * source_fe.shape_grad (0, point) - interpolation_matrix (1, dof) * source_fe.shape_grad (1, point)); + + for (unsigned int d = 0; d < dim; ++d) + assembling_vector (dim * q_point + d) = grad[d]; + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 2, dof) = solution (i); + } + } + } + } + + break; + } + + case 3: { + switch (subface) { + case 0: { + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point ())) > 1e-14) + interpolation_matrix (0, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point ()); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.0, 0.0))) > 1e-14) + interpolation_matrix (1, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.0, 0.0)); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 0.5, 0.0))) > 1e-14) + interpolation_matrix (2, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 0.5, 0.0)); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0))) > 1e-14) + interpolation_matrix (3, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0)); + } + + if (this->degree > 1) { + double weight; + QGauss reference_edge_quadrature (this->degree + 1); + Quadrature edge_quadrature = QProjector::project_to_child (reference_edge_quadrature, 0); + const unsigned int n_edge_points = reference_edge_quadrature.size (); + QGauss reference_quadrature (this->degree); + Quadrature quadrature = QProjector::project_to_child (reference_quadrature, subface); + const unsigned int n_quadrature_points = quadrature.size (); + FullMatrix assembling_matrix (this->degree - 1, n_edge_points); + FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); + FullMatrix system_matrix_inv (assembling_matrix.m (), assembling_matrix.m ()); + Point point; + PreconditionJacobi > precondition; + std::vector > edge_quadrature_points = edge_quadrature.get_points (); + std::vector > quadrature_points = quadrature.get_points (); + Tensor<1, dim> grad; + Vector assembling_vector (n_edge_points); + Vector solution (assembling_matrix.m ()); + Vector system_rhs (assembling_matrix.m ()); + + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (0.0, 2.0 * edge_quadrature_points[q_point] (0), 0.0); + weight = std::sqrt (edge_quadrature.weight (q_point)); + + for (unsigned int i = 0; i < system_matrix.m (); ++i) + assembling_matrix (i, q_point) = weight * this->shape_value (i + 8, point); + } + + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_matrix_inv.invert (system_matrix); + + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) + switch (line) { + case 0: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (0.0, 2.0 * edge_quadrature_points[q_point] (0), 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, edge_quadrature_points[q_point] (0), 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 4, dof) = solution (i); + + break; + } + + case 1: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (1.0, 2.0 * edge_quadrature_points[q_point] (0), 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, edge_quadrature_points[q_point] (0), 0.0)) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + source_fe.degree + 3, dof) = solution (i); + + break; + } + + case 2: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (2.0 * edge_quadrature_points[q_point] (0), 0.0, 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_points[q_point] (0), 0.0, 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 2 * source_fe.degree + 2, dof) = solution (i); + + break; + } + + case 3: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (2.0 * edge_quadrature_points[q_point] (0), 1.0, 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_points[q_point] (0), 0.5, 0.0)) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 3 * source_fe.degree + 1, dof) = solution (i); + } + } + + assembling_matrix.reinit ((this->degree - 1) * (this->degree - 1), dim * n_quadrature_points); + assembling_vector.reinit (assembling_matrix.n ()); + + for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { + point = Point (2.0 * quadrature_points[q_point] (0), 2.0 * quadrature_points[q_point] (1), 0.0); + grad = this->shape_grad (this->face_to_cell_index (dof, 4), Point (quadrature_points[q_point] (0), quadrature_points[q_point] (1), 0.0)); + + for (unsigned int vertex = 0; vertex < GeometryInfo::vertices_per_face; ++vertex) + grad -= interpolation_matrix (vertex, dof) * source_fe.shape_grad (vertex, point); + + for (unsigned int line = 0; line < GeometryInfo::lines_per_face; ++line) + for (unsigned int i = 0; i < this->degree - 1; ++i) + grad -= interpolation_matrix (i + line * (source_fe.degree - 1) + 4, dof) * source_fe.shape_grad (i + line * (source_fe.degree - 1) + 8, point); + + weight = std::sqrt (quadrature.weight (q_point)); + grad *= weight; + + for (unsigned int d = 0; d < dim; ++d) + assembling_vector (dim * q_point + d) = grad[d]; + + for (unsigned int i = 0; i < assembling_matrix.m (); ++i) { + grad = weight * this->shape_grad (i + 4 * ((this->degree + 2) * (this->degree - 1) + 2), point); + + for (unsigned int d = 0; d < dim; ++d) + assembling_matrix (i, dim * q_point + d) = grad[d]; + } + } + + system_matrix.reinit (assembling_matrix.m (), assembling_matrix.n ()); + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_rhs.reinit (system_matrix.m ()); + assembling_matrix.vmult (system_rhs, assembling_vector); + solution.reinit (system_matrix.m ()); + + SolverControl solver_control (system_matrix.m (), 1e-13, false, false); + SolverCG<> cg (solver_control); + + precondition.initialize (system_matrix); + cg.solve (system_matrix, solution, system_rhs, precondition); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 4 * source_fe.degree, dof) = solution (i); + } + } + + break; + } + + case 1: { + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point ())) > 1e-14) + interpolation_matrix (0, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.0, 0.0)); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.0, 0.0))) > 1e-14) + interpolation_matrix (1, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (1.0, 0.0, 0.0)); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 0.5, 0.0))) > 1e-14) + interpolation_matrix (2, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0)); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0))) > 1e-14) + interpolation_matrix (3, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (1.0, 0.5, 0.0)); + } + + if (this->degree > 1) { + double weight; + QGauss reference_edge_quadrature (this->degree + 1); + Quadrature edge_quadrature_x = QProjector::project_to_child (reference_edge_quadrature, 1); + Quadrature edge_quadrature_y = QProjector::project_to_child (reference_edge_quadrature, 0); + const unsigned int n_edge_points = reference_edge_quadrature.size (); + QGauss reference_quadrature (this->degree); + Quadrature quadrature = QProjector::project_to_child (reference_quadrature, subface); + const unsigned int n_quadrature_points = quadrature.size (); + FullMatrix assembling_matrix (this->degree - 1, n_edge_points); + FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); + FullMatrix system_matrix_inv (assembling_matrix.m (), assembling_matrix.m ()); + Point point; + PreconditionJacobi > precondition; + std::vector > edge_quadrature_x_points = edge_quadrature_x.get_points (); + std::vector > edge_quadrature_y_points = edge_quadrature_y.get_points (); + std::vector > quadrature_points = quadrature.get_points (); + Tensor<1, dim> grad; + Vector assembling_vector (n_edge_points); + Vector solution (assembling_matrix.m ()); + Vector system_rhs (assembling_matrix.m ()); + + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (0.0, 2.0 * edge_quadrature_y_points[q_point] (0), 0.0); + weight = std::sqrt (edge_quadrature_y.weight (q_point)); + + for (unsigned int i = 0; i < system_matrix.m (); ++i) + assembling_matrix (i, q_point) = weight * this->shape_value (i + 8, point); + } + + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_matrix_inv.invert (system_matrix); + + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) + switch (line) { + case 0: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (0.0, 2.0 * edge_quadrature_y_points[q_point] (0), 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature_y.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, edge_quadrature_y_points[q_point] (0), 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 4, dof) = solution (i); + + break; + } + + case 1: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (1.0, 2.0 * edge_quadrature_y_points[q_point] (0), 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature_y.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (1.0, edge_quadrature_y_points[q_point] (0), 0.0)) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + source_fe.degree + 3, dof) = solution (i); + + break; + } + + case 2: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (2.0 * edge_quadrature_x_points[q_point] (0) - 1.0, 0.0, 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature_x.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_x_points[q_point] (0), 0.0, 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 2 * source_fe.degree + 2, dof) = solution (i); + + break; + } + + case 3: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (2.0 * edge_quadrature_x_points[q_point] (0) - 1.0, 1.0, 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature_x.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_x_points[q_point] (0), 0.5, 0.0)) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 3 * source_fe.degree + 1, dof) = solution (i); + } + } + + assembling_matrix.reinit ((this->degree - 1) * (this->degree - 1), dim * n_quadrature_points); + assembling_vector.reinit (assembling_matrix.n ()); + + for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { + point = Point (2.0 * quadrature_points[q_point] (0) - 1.0, 2.0 * quadrature_points[q_point] (1), 0.0); + grad = this->shape_grad (this->face_to_cell_index (dof, 4), Point (quadrature_points[q_point] (0), quadrature_points[q_point] (1), 0.0)); + + for (unsigned int vertex = 0; vertex < GeometryInfo::vertices_per_face; ++vertex) + grad -= interpolation_matrix (vertex, dof) * source_fe.shape_grad (vertex, point); + + for (unsigned int line = 0; line < GeometryInfo::lines_per_face; ++line) + for (unsigned int i = 0; i < this->degree - 1; ++i) + grad -= interpolation_matrix (i + line * (source_fe.degree - 1) + 4, dof) * source_fe.shape_grad (i + line * (source_fe.degree - 1) + 8, point); + + weight = std::sqrt (quadrature.weight (q_point)); + grad *= weight; + + for (unsigned int d = 0; d < dim; ++d) + assembling_vector (dim * q_point + d) = grad[d]; + + for (unsigned int i = 0; i < assembling_matrix.m (); ++i) { + grad = weight * this->shape_grad (i + 4 * ((this->degree + 2) * (this->degree - 1) + 2), point); + + for (unsigned int d = 0; d < dim; ++d) + assembling_matrix (i, dim * q_point + d) = grad[d]; + } + } + + system_matrix.reinit (assembling_matrix.m (), assembling_matrix.n ()); + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_rhs.reinit (system_matrix.m ()); + assembling_matrix.vmult (system_rhs, assembling_vector); + solution.reinit (system_matrix.m ()); + + SolverControl solver_control (system_matrix.m (), 1e-13, false, false); + SolverCG<> cg (solver_control); + + precondition.initialize (system_matrix); + cg.solve (system_matrix, solution, system_rhs, precondition); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 4 * source_fe.degree, dof) = solution (i); + } + } + + break; + } + + case 2: { + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point ())) > 1e-14) + interpolation_matrix (0, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 0.5, 0.0)); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.0, 0.0))) > 1e-14) + interpolation_matrix (1, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0)); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 0.5, 0.0))) > 1e-14) + interpolation_matrix (2, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 1.0, 0.0)); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0))) > 1e-14) + interpolation_matrix (3, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 1.0, 0.0)); + } + + if (this->degree > 1) { + double weight; + QGauss reference_edge_quadrature (this->degree + 1); + Quadrature edge_quadrature_x = QProjector::project_to_child (reference_edge_quadrature, 0); + Quadrature edge_quadrature_y = QProjector::project_to_child (reference_edge_quadrature, 1); + const unsigned int n_edge_points = reference_edge_quadrature.size (); + QGauss reference_quadrature (this->degree); + Quadrature quadrature = QProjector::project_to_child (reference_quadrature, subface); + const unsigned int n_quadrature_points = quadrature.size (); + FullMatrix assembling_matrix (this->degree - 1, n_edge_points); + FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); + FullMatrix system_matrix_inv (assembling_matrix.m (), assembling_matrix.m ()); + Point point; + PreconditionJacobi > precondition; + std::vector > edge_quadrature_x_points = edge_quadrature_x.get_points (); + std::vector > edge_quadrature_y_points = edge_quadrature_y.get_points (); + std::vector > quadrature_points = quadrature.get_points (); + Tensor<1, dim> grad; + Vector assembling_vector (n_edge_points); + Vector solution (assembling_matrix.m ()); + Vector system_rhs (assembling_matrix.m ()); + + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (0.0, 2.0 * edge_quadrature_y_points[q_point] (0) - 1.0, 0.0); + weight = std::sqrt (edge_quadrature_y.weight (q_point)); + + for (unsigned int i = 0; i < system_matrix.m (); ++i) + assembling_matrix (i, q_point) = weight * this->shape_value (i + 8, point); + } + + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_matrix_inv.invert (system_matrix); + + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) + switch (line) { + case 0: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (0.0, 2.0 * edge_quadrature_y_points[q_point] (0) - 1.0, 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature_y.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, edge_quadrature_y_points[q_point] (0), 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 4, dof) = solution (i); + + break; + } + + case 1: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (1.0, 2.0 * edge_quadrature_y_points[q_point] (0) - 1.0, 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature_y.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, edge_quadrature_y_points[q_point] (0), 0.0)) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + source_fe.degree + 3, dof) = solution (i); + + break; + } + + case 2: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (2.0 * edge_quadrature_x_points[q_point] (0), 0.0, 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature_x.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_x_points[q_point] (0), 0.5, 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 2 * source_fe.degree + 2, dof) = solution (i); + + break; + } + + case 3: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (2.0 * edge_quadrature_x_points[q_point] (0), 1.0, 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature_x.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_x_points[q_point] (0), 1.0, 0.0)) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 3 * source_fe.degree + 1, dof) = solution (i); + } + } + + assembling_matrix.reinit ((this->degree - 1) * (this->degree - 1), dim * n_quadrature_points); + assembling_vector.reinit (assembling_matrix.n ()); + + for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { + point = Point (2.0 * quadrature_points[q_point] (0), 2.0 * quadrature_points[q_point] (1) - 1.0, 0.0); + grad = this->shape_grad (this->face_to_cell_index (dof, 4), Point (quadrature_points[q_point] (0), quadrature_points[q_point] (1), 0.0)); + + for (unsigned int vertex = 0; vertex < GeometryInfo::vertices_per_face; ++vertex) + grad -= interpolation_matrix (vertex, dof) * source_fe.shape_grad (vertex, point); + + for (unsigned int line = 0; line < GeometryInfo::lines_per_face; ++line) + for (unsigned int i = 0; i < this->degree - 1; ++i) + grad -= interpolation_matrix (i + line * (source_fe.degree - 1) + 4, dof) * source_fe.shape_grad (i + line * (source_fe.degree - 1) + 8, point); + + weight = std::sqrt (quadrature.weight (q_point)); + grad *= weight; + + for (unsigned int d = 0; d < dim; ++d) + assembling_vector (dim * q_point + d) = grad[d]; + + for (unsigned int i = 0; i < assembling_matrix.m (); ++i) { + grad = weight * this->shape_grad (i + 4 * ((this->degree + 2) * (this->degree - 1) + 2), point); + + for (unsigned int d = 0; d < dim; ++d) + assembling_matrix (i, dim * q_point + d) = grad[d]; + } + } + + system_matrix.reinit (assembling_matrix.m (), assembling_matrix.n ()); + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_rhs.reinit (system_matrix.m ()); + assembling_matrix.vmult (system_rhs, assembling_vector); + solution.reinit (system_matrix.m ()); + + SolverControl solver_control (system_matrix.m (), 1e-13, false, false); + SolverCG<> cg (solver_control); + + precondition.initialize (system_matrix); + cg.solve (system_matrix, solution, system_rhs, precondition); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 4 * source_fe.degree, dof) = solution (i); + } + } + + break; + } + + case 3: { + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point ())) > 1e-14) + interpolation_matrix (0, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0)); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.0, 0.0))) > 1e-14) + interpolation_matrix (1, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (1.0, 0.5, 0.0)); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 0.5, 0.0))) > 1e-14) + interpolation_matrix (2, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 1.0, 0.0)); + + if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0))) > 1e-14) + interpolation_matrix (3, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (1.0, 1.0, 0.0)); + } + + if (this->degree > 1) { + double weight; + QGauss reference_edge_quadrature (this->degree + 1); + Quadrature edge_quadrature = QProjector::project_to_child (reference_edge_quadrature, 1); + const unsigned int n_edge_points = reference_edge_quadrature.size (); + QGauss reference_quadrature (this->degree); + Quadrature quadrature = QProjector::project_to_child (reference_quadrature, subface); + const unsigned int n_quadrature_points = quadrature.size (); + FullMatrix assembling_matrix (this->degree - 1, n_edge_points); + FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); + FullMatrix system_matrix_inv (assembling_matrix.m (), assembling_matrix.m ()); + Point point; + PreconditionJacobi > precondition; + std::vector > edge_quadrature_points = edge_quadrature.get_points (); + std::vector > quadrature_points = quadrature.get_points (); + Tensor<1, dim> grad; + Vector assembling_vector (n_edge_points); + Vector solution (assembling_matrix.m ()); + Vector system_rhs (assembling_matrix.m ()); + + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (0.0, 2.0 * edge_quadrature_points[q_point] (0) - 1.0, 0.0); + weight = std::sqrt (edge_quadrature.weight (q_point)); + + for (unsigned int i = 0; i < system_matrix.m (); ++i) + assembling_matrix (i, q_point) = weight * this->shape_value (i + 8, point); + } + + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_matrix_inv.invert (system_matrix); + + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) + switch (line) { + case 0: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (0.0, 2.0 * edge_quadrature_points[q_point] (0) - 1.0, 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, edge_quadrature_points[q_point] (0), 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 4, dof) = solution (i); + + break; + } + + case 1: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (1.0, 2.0 * edge_quadrature_points[q_point] (0) - 1.0, 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (1.0, edge_quadrature_points[q_point] (0), 0.0)) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + source_fe.degree + 3, dof) = solution (i); + + break; + } + + case 2: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (2.0 * edge_quadrature_points[q_point] (0) - 1.0, 0.0, 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_points[q_point] (0), 0.5, 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 2 * source_fe.degree + 2, dof) = solution (i); + + break; + } + + case 3: { + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { + point = Point (2.0 * edge_quadrature_points[q_point] (0) - 1.0, 1.0, 0.0); + assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_points[q_point] (0), 1.0, 0.0)) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); + } + + assembling_matrix.vmult (system_rhs, assembling_vector); + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 3 * source_fe.degree + 1, dof) = solution (i); + } + } + + assembling_matrix.reinit ((this->degree - 1) * (this->degree - 1), dim * n_quadrature_points); + assembling_vector.reinit (assembling_matrix.n ()); + + for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { + point = Point (2.0 * quadrature_points[q_point] (0) - 1.0, 2.0 * quadrature_points[q_point] (1) - 1.0, 0.0); + grad = this->shape_grad (this->face_to_cell_index (dof, 4), Point (quadrature_points[q_point] (0), quadrature_points[q_point] (1), 0.0)); + + for (unsigned int vertex = 0; vertex < GeometryInfo::vertices_per_face; ++vertex) + grad -= interpolation_matrix (vertex, dof) * source_fe.shape_grad (vertex, point); + + for (unsigned int line = 0; line < GeometryInfo::lines_per_face; ++line) + for (unsigned int i = 0; i < this->degree - 1; ++i) + grad -= interpolation_matrix (i + line * (source_fe.degree - 1) + 4, dof) * source_fe.shape_grad (i + line * (source_fe.degree - 1) + 8, point); + + weight = std::sqrt (quadrature.weight (q_point)); + grad *= weight; + + for (unsigned int d = 0; d < dim; ++d) + assembling_vector (dim * q_point + d) = grad[d]; + + for (unsigned int i = 0; i < assembling_matrix.m (); ++i) { + grad = weight * this->shape_grad (i + 4 * ((this->degree + 2) * (this->degree - 1) + 2), point); + + for (unsigned int d = 0; d < dim; ++d) + assembling_matrix (i, dim * q_point + d) = grad[d]; + } + } + + system_matrix.reinit (assembling_matrix.m (), assembling_matrix.n ()); + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_rhs.reinit (system_matrix.m ()); + assembling_matrix.vmult (system_rhs, assembling_vector); + solution.reinit (system_matrix.m ()); + + SolverControl solver_control (system_matrix.m (), 1e-13, false, false); + SolverCG<> cg (solver_control); + + precondition.initialize (system_matrix); + cg.solve (system_matrix, solution, system_rhs, precondition); + + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 4 * source_fe.degree, dof) = solution (i); + } + } + } + } + } + } +} + #endif -- 2.39.5