From 8e6a4484d662ef11edfd9f2335d11c46a00bedf9 Mon Sep 17 00:00:00 2001 From: buerg Date: Wed, 25 Jul 2012 18:42:23 +0000 Subject: [PATCH] Fix introduced bug. git-svn-id: https://svn.dealii.org/trunk@25732 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/source/fe/fe_nedelec.cc | 4797 ++++++++++++++++++++----------- 1 file changed, 3093 insertions(+), 1704 deletions(-) diff --git a/deal.II/source/fe/fe_nedelec.cc b/deal.II/source/fe/fe_nedelec.cc index a2791856c2..e1c05bb37b 100644 --- a/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/source/fe/fe_nedelec.cc @@ -161,6 +161,7 @@ FE_PolyTensor, dim> default: Assert (false, ExcNotImplemented ()); } + } @@ -496,8 +497,6 @@ FE_Nedelec<1>::initialize_restriction () - - // Restriction operator template void @@ -516,9 +515,6 @@ FE_Nedelec::initialize_restriction () n_edge_quadrature_points = edge_quadrature.size (); const unsigned int index = RefinementCase::isotropic_refinement - 1; - const unsigned int deg = this->degree-1; - const std::vector >& - legendre_polynomials = Polynomials::Legendre::generate_complete_basis (deg); switch (dim) { @@ -526,83 +522,232 @@ FE_Nedelec::initialize_restriction () { // First interpolate the shape // functions of the child cells - // to the face shape functions - // of the parent cell. + // to the lowest order shape + // functions of the parent cell. for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; ++q_point) { + const double weight = 2.0 * edge_quadrature.weight (q_point); + if (edge_quadrature_points[q_point] (0) < 0.5) - for (unsigned int i = 0; i <= deg; ++i) - { - const double weight - = 2.0 * edge_quadrature.weight (q_point) - * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0)); - Point quadrature_point (0.0, - 2.0 * edge_quadrature_points[q_point] (0)); - - this->restriction[index][0] (i, dof) += weight - * this->shape_value_component - (dof, - quadrature_point, - 1); - quadrature_point (0) = 1.0; - this->restriction[index][1] (i + this->degree, dof) - += weight * this->shape_value_component (dof, - quadrature_point, - 1); - quadrature_point (0) = quadrature_point (1); - quadrature_point (1) = 0.0; - this->restriction[index][0] (i + 2 * this->degree, dof) - += weight * this->shape_value_component (dof, - quadrature_point, - 0); - quadrature_point (1) = 1.0; - this->restriction[index][2] (i + 3 * this->degree, dof) - += weight * this->shape_value_component (dof, - quadrature_point, - 0); - } + { + Point quadrature_point (0.0, + 2.0 * edge_quadrature_points[q_point] (0)); + + this->restriction[index][0] (0, dof) += weight + * this->shape_value_component + (dof, + quadrature_point, + 1); + quadrature_point (0) = 1.0; + this->restriction[index][1] (this->degree, dof) + += weight * this->shape_value_component (dof, + quadrature_point, + 1); + quadrature_point (0) = quadrature_point (1); + quadrature_point (1) = 0.0; + this->restriction[index][0] (2 * this->degree, dof) + += weight * this->shape_value_component (dof, + quadrature_point, + 0); + quadrature_point (1) = 1.0; + this->restriction[index][2] (3 * this->degree, dof) + += weight * this->shape_value_component (dof, + quadrature_point, + 0); + } else - for (unsigned int i = 0; i <= deg; ++i) - { - const double weight - = 2.0 * edge_quadrature.weight (q_point) - * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0)); - Point quadrature_point (0.0, - 2.0 * edge_quadrature_points[q_point] (0) - - 1.0); - - this->restriction[index][2] (i, dof) += weight - * this->shape_value_component - (dof, - quadrature_point, - 1); - quadrature_point (0) = 1.0; - this->restriction[index][3] (i + this->degree, dof) - += weight * this->shape_value_component (dof, - quadrature_point, - 1); - quadrature_point (0) = quadrature_point (1); - quadrature_point (1) = 0.0; - this->restriction[index][1] (i + 2 * this->degree, dof) - += weight * this->shape_value_component (dof, - quadrature_point, - 0); - quadrature_point (1) = 1.0; - this->restriction[index][3] (i + 3 * this->degree, dof) - += weight * this->shape_value_component (dof, - quadrature_point, - 0); - } + { + Point quadrature_point (0.0, + 2.0 * edge_quadrature_points[q_point] (0) + - 1.0); + + this->restriction[index][2] (0, dof) += weight + * this->shape_value_component + (dof, + quadrature_point, + 1); + quadrature_point (0) = 1.0; + this->restriction[index][3] (this->degree, dof) + += weight * this->shape_value_component (dof, + quadrature_point, + 1); + quadrature_point (0) = quadrature_point (1); + quadrature_point (1) = 0.0; + this->restriction[index][1] (2 * this->degree, dof) + += weight * this->shape_value_component (dof, + quadrature_point, + 0); + quadrature_point (1) = 1.0; + this->restriction[index][3] (3 * this->degree, dof) + += weight * this->shape_value_component (dof, + quadrature_point, + 0); + } } - // Then interpolate the shape functions - // of the child cells to the interior - // shape functions of the parent cell. + // Then project the shape functions + // of the child cells to the higher + // order shape functions of the + // parent cell. if (this->degree > 1) { + const unsigned int deg = this->degree-1; + const std::vector >& + legendre_polynomials + = Polynomials::Legendre::generate_complete_basis (deg); + FullMatrix system_matrix_inv (deg, deg); + + { + FullMatrix assembling_matrix (deg, + n_edge_quadrature_points); + + for (unsigned int q_point = 0; + q_point < n_edge_quadrature_points; ++q_point) + { + const double weight + = std::sqrt (edge_quadrature.weight (q_point)); + + for (unsigned int i = 0; i < deg; ++i) + assembling_matrix (i, q_point) = weight + * legendre_polynomials[i + 1].value + (edge_quadrature_points[q_point] (0)); + } + + FullMatrix system_matrix (deg, deg); + + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_matrix_inv.invert (system_matrix); + } + + FullMatrix solution (this->degree-1, 4); + FullMatrix system_rhs (this->degree-1, 4); + Vector tmp (4); + + for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) + for (unsigned int i = 0; i < 2; ++i) + { + system_rhs = 0.0; + + for (unsigned int q_point = 0; + q_point < n_edge_quadrature_points; ++q_point) + { + const double weight + = edge_quadrature.weight (q_point); + const Point quadrature_point_0 (i, + edge_quadrature_points[q_point] (0)); + const Point quadrature_point_1 + (edge_quadrature_points[q_point] (0), + i); + + if (edge_quadrature_points[q_point] (0) < 0.5) + { + Point quadrature_point_2 (i, + 2.0 * edge_quadrature_points[q_point] (0)); + + tmp (0) = weight + * (2.0 * this->shape_value_component + (dof, quadrature_point_2, 1) + - this->restriction[index][i] + (i * this->degree, dof) + * this->shape_value_component + (i * this->degree, + quadrature_point_0, 1)); + tmp (1) = -1.0 * weight + * this->restriction[index][i + 2] + (i * this->degree, dof) + * this->shape_value_component + (i * this->degree, + quadrature_point_0, 1); + quadrature_point_2 + = Point (2.0 * edge_quadrature_points[q_point] (0), + i); + tmp (2) = weight + * (2.0 * this->shape_value_component + (dof, quadrature_point_2, 0) + - this->restriction[index][2 * i] + ((i + 2) * this->degree, dof) + * this->shape_value_component + ((i + 2) * this->degree, + quadrature_point_1, 0)); + tmp (3) = -1.0 * weight + * this->restriction[index][2 * i + 1] + ((i + 2) * this->degree, dof) + * this->shape_value_component + ((i + 2) * this->degree, + quadrature_point_1, 0); + } + + else + { + tmp (0) = -1.0 * weight + * this->restriction[index][i] + (i * this->degree, dof) + * this->shape_value_component + (i * this->degree, + quadrature_point_0, 1); + + Point quadrature_point_2 (i, + 2.0 * edge_quadrature_points[q_point] (0) + - 1.0); + + tmp (1) = weight + * (2.0 * this->shape_value_component + (dof, quadrature_point_2, 1) + - this->restriction[index][i + 2] + (i * this->degree, dof) + * this->shape_value_component + (i * this->degree, + quadrature_point_0, 1)); + tmp (2) = -1.0 * weight + * this->restriction[index][2 * i] + ((i + 2) * this->degree, dof) + * this->shape_value_component + ((i + 2) * this->degree, + quadrature_point_1, 0); + quadrature_point_2 + = Point (2.0 * edge_quadrature_points[q_point] (0) + - 1.0, i); + tmp (3) = weight + * (2.0 * this->shape_value_component + (dof, quadrature_point_2, 0) + - this->restriction[index][2 * i + 1] + ((i + 2) * this->degree, dof) + * this->shape_value_component + ((i + 2) * this->degree, + quadrature_point_1, 0)); + } + + for (unsigned int j = 0; j < this->degree-1; ++j) + { + const double L_j + = legendre_polynomials[j + 1].value + (edge_quadrature_points[q_point] (0)); + + for (unsigned int k = 0; k < tmp.size (); ++k) + system_rhs (j, k) += tmp (k) * L_j; + } + } + + system_matrix_inv.mmult (solution, system_rhs); + + for (unsigned int j = 0; j < this->degree-1; ++j) + for (unsigned int k = 0; k < 2; ++k) + { + if (std::abs (solution (j, k)) > 1e-14) + this->restriction[index][i + 2 * k] + (i * this->degree + j + 1, dof) + = solution (j, k); + + if (std::abs (solution (j, k + 2)) > 1e-14) + this->restriction[index][2 * i + k] + ((i + 2) * this->degree + j + 1, dof) + = solution (j, k + 2); + } + } + const QGauss quadrature (2 * this->degree); const std::vector >& quadrature_points = quadrature.get_points (); @@ -612,71 +757,185 @@ FE_Nedelec::initialize_restriction () (this->degree); const unsigned int n_boundary_dofs = GeometryInfo::faces_per_cell * this->degree; + const unsigned int& n_quadrature_points = quadrature.size (); - for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) - for (unsigned int q_point = 0; q_point < quadrature.size (); + { + FullMatrix assembling_matrix ((this->degree-1) * this->degree, + n_quadrature_points); + + for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { - const double weight = 2.0 * quadrature.weight (q_point); - Point quadrature_point (2.0 * quadrature_points[q_point] (0), - 2.0 * quadrature_points[q_point] (1)); - unsigned int k; + const double weight + = std::sqrt (quadrature.weight (q_point)); - if (quadrature_points[q_point] (0) < 0.5) + for (unsigned int i = 0; i < this->degree; ++i) { + const double L_i = weight + * legendre_polynomials[i].value + (quadrature_points[q_point] (0)); + + for (unsigned int j = 0; j < this->degree-1; ++j) + assembling_matrix (i * (this->degree-1) + j, q_point) + = L_i * lobatto_polynomials[j + 2].value + (quadrature_points[q_point] (1)); + } + } + + FullMatrix system_matrix (assembling_matrix.m (), + assembling_matrix.m ()); + + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_matrix_inv.reinit (system_matrix.m (), system_matrix.m ()); + system_matrix_inv.invert (system_matrix); + } + + solution.reinit (system_matrix_inv.m (), 8); + system_rhs.reinit (system_matrix_inv.m (), 8); + tmp.reinit (8); + + for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) + { + system_rhs = 0.0; + + for (unsigned int q_point = 0; + q_point < n_quadrature_points; ++q_point) + { + tmp = 0.0; + + if (quadrature_points[q_point] (0) < 0.5) + { + if (quadrature_points[q_point] (1) < 0.5) + { + const Point quadrature_point + (2.0 * quadrature_points[q_point] (0), + 2.0 * quadrature_points[q_point] (1)); + + tmp (0) += 2.0 * this->shape_value_component + (dof, quadrature_point, 0); + tmp (1) += 2.0 * this->shape_value_component + (dof, quadrature_point, 1); + } + + else + { + const Point quadrature_point + (2.0 * quadrature_points[q_point] (0), + 2.0 * quadrature_points[q_point] (1) + - 1.0); + + tmp (4) += 2.0 * this->shape_value_component + (dof, quadrature_point, 0); + tmp (5) += 2.0 * this->shape_value_component + (dof, quadrature_point, 1); + } + } + + else if (quadrature_points[q_point] (1) < 0.5) - k = 0; - + { + const Point quadrature_point + (2.0 * quadrature_points[q_point] (0) + - 1.0, + 2.0 * quadrature_points[q_point] (1)); + + tmp (2) += 2.0 * this->shape_value_component + (dof, quadrature_point, 0); + tmp (3) += 2.0 * this->shape_value_component + (dof, quadrature_point, 1); + } + else { - quadrature_point (1) -= 1.0; - k = 1; + const Point quadrature_point + (2.0 * quadrature_points[q_point] (0) + - 1.0, + 2.0 * quadrature_points[q_point] (1) + - 1.0); + + tmp (6) += 2.0 * this->shape_value_component + (dof, quadrature_point, 0); + tmp (7) += 2.0 * this->shape_value_component + (dof, quadrature_point, 1); } - } - - else - if (quadrature_points[q_point] (1) < 0.5) + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < this->degree; ++j) + { + tmp (2 * i) -= this->restriction[index][i] + (j + 2 * this->degree, dof) + * this->shape_value_component + (j + 2 * this->degree, + quadrature_points[q_point], 0); + tmp (2 * i + 1) -= this->restriction[index][i] + (i * this->degree + j, dof) + * this->shape_value_component + (i * this->degree + j, + quadrature_points[q_point], 1); + tmp (2 * (i + 2)) -= this->restriction[index][i + 2] + (j + 3 * this->degree, dof) + * this->shape_value_component + (j + 3 * this->degree, + quadrature_points[q_point], + 0); + tmp (2 * i + 5) -= this->restriction[index][i + 2] + (i * this->degree + j, dof) + * this->shape_value_component + (i * this->degree + j, + quadrature_points[q_point], 1); + } + + tmp *= quadrature.weight (q_point); + + for (unsigned int i = 0; i < this->degree; ++i) { - quadrature_point (0) -= 1.0; - k = 2; + const double L_i_0 + = legendre_polynomials[i].value + (quadrature_points[q_point] (0)); + const double L_i_1 + = legendre_polynomials[i].value + (quadrature_points[q_point] (1)); + + for (unsigned int j = 0; j < this->degree-1; ++j) + { + const double l_j_0 + = L_i_0 * lobatto_polynomials[j + 2].value + (quadrature_points[q_point] (1)); + const double l_j_1 + = L_i_1 * lobatto_polynomials[j + 2].value + (quadrature_points[q_point] (0)); + + for (unsigned int k = 0; k < 4; ++k) + { + system_rhs (i * (this->degree-1) + j, 2 * k) + += tmp (2 * k) * l_j_0; + system_rhs (i * (this->degree-1) + j, 2 * k + 1) + += tmp (2 * k + 1) * l_j_1; + } + } } - - else + } + + system_matrix_inv.mmult (solution, system_rhs); + + for (unsigned int i = 0; i < this->degree; ++i) + for (unsigned int j = 0; j < this->degree-1; ++j) + for (unsigned int k = 0; k < 4; ++k) { - quadrature_point (0) -= 1.0; - quadrature_point (1) -= 1.0; - k = 3; + if (std::abs (solution (i * (this->degree-1) + j, 2 * k)) + > 1e-14) + this->restriction[index][k] + (i * (this->degree-1) + j + n_boundary_dofs, dof) + = solution (i * (this->degree-1) + j, 2 * k); + + if (std::abs (solution (i * (this->degree-1) + j, 2 * k + 1)) + > 1e-14) + this->restriction[index][k] + (i + (this->degree-1 + j) * this->degree + n_boundary_dofs, + dof) + = solution (i * (this->degree-1) + j, 2 * k + 1); } - - for (unsigned int i = 0; i < this->degree; ++i) - { - const double L_i_0 = weight - * legendre_polynomials[i].value - (quadrature_points[q_point] (0)); - const double L_i_1 = weight - * legendre_polynomials[i].value - (quadrature_points[q_point] (1)); - - for (unsigned int j = 0; j < deg; ++j) - { - this->restriction[index][k] (i * deg + j + n_boundary_dofs, dof) - += L_i_0 * (this->shape_grad_component (dof, quadrature_point, 0)[1] - * legendre_polynomials[j + 1].value - (quadrature_points[q_point] (1)) - + this->shape_value_component (dof, quadrature_point, 0) - * lobatto_polynomials[j + 2].value - (quadrature_points[q_point] (1))); - this->restriction[index][k] (i + (deg + j) * this->degree - + n_boundary_dofs, dof) - += L_i_1 * (this->shape_grad_component (dof, quadrature_point, 1)[0] - * legendre_polynomials[j + 1].value - (quadrature_points[q_point] (0)) - + this->shape_value_component (dof, quadrature_point, 1) - * lobatto_polynomials[j + 2].value - (quadrature_points[q_point] (0))); - } - } - } + } } break; @@ -686,365 +945,1017 @@ FE_Nedelec::initialize_restriction () { // First interpolate the shape // functions of the child cells - // to the edge shape functions - // of the parent cell. + // to the lowest order shape + // functions of the parent cell. for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; ++q_point) { const double weight = 2.0 * edge_quadrature.weight (q_point); - - if (edge_quadrature_points[q_point] (0) < 0.5) - for (unsigned int i = 0; i < this->degree; ++i) - { - const double L_i - = weight * legendre_polynomials[i].value - (edge_quadrature_points[q_point] (0)); - for (unsigned int j = 0; j < 2; ++j) - for (unsigned int k = 0; k < 2; ++k) - { - Point quadrature_point (j, - 2.0 * edge_quadrature_points[q_point] (0), - k); - - this->restriction[index][j + 4 * k] - (i + (j + 4 * k) * this->degree, dof) - += L_i * this->shape_value_component (dof, - quadrature_point, - 1); - quadrature_point - = Point (2.0 * edge_quadrature_points[q_point] (0), - j, k); - this->restriction[index][2 * (j + 2 * k)] - (i + (j + 4 * k + 2) * this->degree, dof) - += L_i * this->shape_value_component (dof, - quadrature_point, - 0); - quadrature_point = Point (j, k, - 2.0 * edge_quadrature_points[q_point] (0)); - this->restriction[index][j + 2 * k] - (i + (j + 2 * (k + 4)) * this->degree, dof) - += L_i * this->shape_value_component (dof, - quadrature_point, - 2); - } - } - - else - for (unsigned int i = 0; i < this->degree; ++i) - { - const double L_i - = weight * legendre_polynomials[i].value - (edge_quadrature_points[q_point] (0)); + if (edge_quadrature_points[q_point] (0) < 0.5) + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < 2; ++j) + { + Point quadrature_point (i, + 2.0 * edge_quadrature_points[q_point] (0), + j); + + this->restriction[index][i + 4 * j] + ((i + 4 * j) * this->degree, dof) + += weight * this->shape_value_component (dof, + quadrature_point, + 1); + quadrature_point + = Point (2.0 * edge_quadrature_points[q_point] (0), + i, j); + this->restriction[index][2 * (i + 2 * j)] + ((i + 4 * j + 2) * this->degree, dof) + += weight * this->shape_value_component (dof, + quadrature_point, + 0); + quadrature_point = Point (i, j, + 2.0 * edge_quadrature_points[q_point] (0)); + this->restriction[index][i + 2 * j] + ((i + 2 * (j + 4)) * this->degree, dof) + += weight * this->shape_value_component (dof, + quadrature_point, + 2); + } - for (unsigned int j = 0; j < 2; ++j) - for (unsigned int k = 0; k < 2; ++k) + else + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < 2; ++j) { - Point quadrature_point (j, + Point quadrature_point (i, 2.0 * edge_quadrature_points[q_point] (0) - - 1.0, k); + - 1.0, j); - this->restriction[index][j + 4 * k + 2] - (i + (j + 4 * k) * this->degree, dof) - += L_i * this->shape_value_component (dof, - quadrature_point, - 1); + this->restriction[index][i + 4 * j + 2] + ((i + 4 * j) * this->degree, dof) + += weight * this->shape_value_component (dof, + quadrature_point, + 1); quadrature_point - = Point (2.0 * edge_quadrature_points[q_point] (0) - 1.0, - j, k); - this->restriction[index][2 * (j + 2 * k) + 1] - (i + (j + 4 * k + 2) * this->degree, dof) - += L_i * this->shape_value_component (dof, - quadrature_point, - 0); - quadrature_point = Point (j, k, + = Point (2.0 * edge_quadrature_points[q_point] (0) + - 1.0, i, j); + this->restriction[index][2 * (i + 2 * j) + 1] + ((i + 4 * j + 2) * this->degree, dof) + += weight * this->shape_value_component (dof, + quadrature_point, + 0); + quadrature_point = Point (i, j, 2.0 * edge_quadrature_points[q_point] (0) - 1.0); - this->restriction[index][j + 2 * (k + 2)] - (i + (j + 2 * (k + 4)) * this->degree, dof) - += L_i * this->shape_value_component (dof, - quadrature_point, - 2); + this->restriction[index][i + 2 * (j + 2)] + ((i + 2 * (j + 4)) * this->degree, dof) + += weight * this->shape_value_component (dof, + quadrature_point, + 2); } - } } - // Then interpolate the shape functions - // of the child cells to the face - // and interior shape functions of - // the parent cell. + // Then project the shape functions + // of the child cells to the higher + // order shape functions of the + // parent cell. if (this->degree > 1) { - const QGauss<2> face_quadrature (2 * this->degree); - const std::vector >& face_quadrature_points - = face_quadrature.get_points (); + const unsigned int deg = this->degree-1; const std::vector >& - lobatto_polynomials - = Polynomials::Lobatto::generate_complete_basis - (this->degree); - const QGauss quadrature (2 * this->degree); - const std::vector >& - quadrature_points = quadrature.get_points (); - const unsigned int n_edge_dofs - = GeometryInfo::lines_per_cell * this->degree; - const unsigned int n_boundary_dofs - = 2 * GeometryInfo::faces_per_cell * deg * this->degree - + n_edge_dofs; - const unsigned int& n_face_quadrature_points - = face_quadrature.size (); - const unsigned int& n_quadrature_points = quadrature.size (); - - // First, the interpolation to - // the face shape functions. - for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) - { - for (unsigned int q_point = 0; q_point < n_face_quadrature_points; - ++q_point) + legendre_polynomials + = Polynomials::Legendre::generate_complete_basis (deg); + FullMatrix system_matrix_inv (deg, deg); + + { + FullMatrix assembling_matrix (deg, + n_edge_quadrature_points); + + for (unsigned int q_point = 0; + q_point < n_edge_quadrature_points; ++q_point) + { + const double weight = std::sqrt (edge_quadrature.weight + (q_point)); + + for (unsigned int i = 0; i < deg; ++i) + assembling_matrix (i, q_point) = weight + * legendre_polynomials[i + 1].value + (edge_quadrature_points[q_point] (0)); + } + + FullMatrix system_matrix (deg, deg); + + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_matrix_inv.invert (system_matrix); + } + + FullMatrix solution (deg, 6); + FullMatrix system_rhs (deg, 6); + Vector tmp (6); + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < 2; ++j) + for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) { - const double weight = 2.0 * face_quadrature.weight (q_point); - - for (unsigned int i = 0; i < 2; ++i) + system_rhs = 0.0; + + for (unsigned int q_point = 0; + q_point < n_edge_quadrature_points; ++q_point) { - Point quadrature_point (i, - 2.0 * face_quadrature_points[q_point] (0), - 2.0 * face_quadrature_points[q_point] (1)); - unsigned int l; - unsigned int m; - - if (face_quadrature_points[q_point] (0) < 0.5) + const double weight = edge_quadrature.weight + (q_point); + const Point quadrature_point_0 (i, + edge_quadrature_points[q_point] (0), + j); + const Point + quadrature_point_1 + (edge_quadrature_points[q_point] (0), i, j); + const Point quadrature_point_2 (i, j, + edge_quadrature_points[q_point] (0)); + + if (edge_quadrature_points[q_point] (0) < 0.5) { - m = 0; - - if (face_quadrature_points[q_point] (1) < 0.5) - l = 0; - - else - { - quadrature_point (2) -= 1.0; - l = 1; - } + Point quadrature_point_3 (i, + 2.0 * edge_quadrature_points[q_point] (0), + j); + + tmp (0) = weight + * (2.0 * this->shape_value_component + (dof, quadrature_point_3, 1) + - this->restriction[index][i + 4 * j] + ((i + 4 * j) * this->degree, + dof) + * this->shape_value_component + ((i + 4 * j) * this->degree, + quadrature_point_0, 1)); + tmp (1) = -1.0 * weight + * this->restriction[index][i + 4 * j + 2] + ((i + 4 * j) * this->degree, + dof) + * this->shape_value_component + ((i + 4 * j) * this->degree, + quadrature_point_0, 1); + quadrature_point_3 + = Point (2.0 * edge_quadrature_points[q_point] (0), + i, j); + tmp (2) = weight + * (2.0 * this->shape_value_component + (dof, quadrature_point_3, 0) + - this->restriction[index][2 * (i + 2 * j)] + ((i + 4 * j + 2) * this->degree, + dof) + * this->shape_value_component + ((i + 4 * j + 2) * this->degree, + quadrature_point_1, 0)); + tmp (3) = -1.0 * weight + * this->restriction[index][2 * (i + 2 * j) + 1] + ((i + 4 * j + 2) * this->degree, + dof) + * this->shape_value_component + ((i + 4 * j + 2) * this->degree, + quadrature_point_1, 0); + quadrature_point_3 = Point (i, j, + 2.0 * edge_quadrature_points[q_point] (0)); + tmp (4) = weight + * (2.0 * this->shape_value_component + (dof, quadrature_point_3, 2) + - this->restriction[index][i + 2 * j] + ((i + 2 * (j + 4)) * this->degree, + dof) + * this->shape_value_component + ((i + 2 * (j + 4)) * this->degree, + quadrature_point_2, 2)); + tmp (5) = -1.0 * weight + * this->restriction[index][i + 2 * (j + 2)] + ((i + 2 * (j + 4)) * this->degree, + dof) + * this->shape_value_component + ((i + 2 * (j + 4)) * this->degree, + quadrature_point_2, 2); } - + else { - quadrature_point (1) -= 1.0; - m = 1; - - if (face_quadrature_points[q_point] (1) < 0.5) - l = 0; - - else - { - quadrature_point (2) -= 1.0; - l = 1; - } + tmp (0) = -1.0 * weight + * this->restriction[index][i + 4 * j] + ((i + 4 * j) * this->degree, + dof) + * this->shape_value_component + ((i + 4 * j) * this->degree, + quadrature_point_0, 1); + + Point quadrature_point_3 (i, + 2.0 * edge_quadrature_points[q_point] (0) + - 1.0, j); + + tmp (1) = weight + * (2.0 * this->shape_value_component + (dof, quadrature_point_3, 1) + - this->restriction[index][i + 4 * j + 2] + ((i + 4 * j) * this->degree, + dof) + * this->shape_value_component + ((i + 4 * j) * this->degree, + quadrature_point_0, 1)); + tmp (2) = -1.0 * weight + * this->restriction[index][2 * (i + 2 * j)] + ((i + 4 * j + 2) * this->degree, + dof) + * this->shape_value_component + ((i + 4 * j + 2) * this->degree, + quadrature_point_1, 0); + quadrature_point_3 + = Point (2.0 * edge_quadrature_points[q_point] (0) + - 1.0, i, j); + tmp (3) = weight + * (2.0 * this->shape_value_component + (dof, quadrature_point_3, 0) + - this->restriction[index][2 * (i + 2 * j) + 1] + ((i + 4 * j + 2) * this->degree, + dof) + * this->shape_value_component + ((i + 4 * j + 2) * this->degree, + quadrature_point_1, 0)); + tmp (4) = -1.0 * weight + * this->restriction[index][i + 2 * j] + ((i + 2 * (j + 4)) * this->degree, + dof) + * this->shape_value_component + ((i + 2 * (j + 4)) * this->degree, + quadrature_point_2, 2); + quadrature_point_3 = Point (i, j, + 2.0 * edge_quadrature_points[q_point] (0) + - 1.0); + tmp (5) = weight + * (2.0 * this->shape_value_component + (dof, quadrature_point_3, 2) + - this->restriction[index][i + 2 * (j + 2)] + ((i + 2 * (j + 4)) * this->degree, + dof) + * this->shape_value_component + ((i + 2 * (j + 4)) * this->degree, + quadrature_point_2, 2)); } - - for (unsigned int j = 0; j < this->degree; ++j) + + for (unsigned int k = 0; k < deg; ++k) { - const double L_j_0 - = weight * legendre_polynomials[j].value - (face_quadrature_points[q_point] (0)); - const double L_j_1 - = weight * legendre_polynomials[j].value - (face_quadrature_points[q_point] (1)); - - for (unsigned int k = 0; k < deg; ++k) - { - const double Le_k_0 - = legendre_polynomials[k + 1].value - (face_quadrature_points[q_point] (0)); - const double Le_k_1 - = legendre_polynomials[k + 1].value - (face_quadrature_points[q_point] (1)); - const double lo_k_0 - = lobatto_polynomials[k + 2].value - (face_quadrature_points[q_point] (0)); - const double lo_k_1 - = lobatto_polynomials[k + 2].value - (face_quadrature_points[q_point] (1)); - - this->restriction[index][i + 2 * (2 * l + m)] - ((2 * i * this->degree + j) * deg + k + n_edge_dofs, dof) - += L_j_0 * (this->shape_grad_component (dof, quadrature_point, 1)[2] - * Le_k_1 - + this->shape_value_component (dof, quadrature_point, 1) - * lo_k_1); - this->restriction[index][i + 2 * (2 * l + m)] - (((2 * i + 1) * deg + k) * this->degree + j + n_edge_dofs, dof) - += L_j_1 * (this->shape_grad_component (dof, quadrature_point, 2)[1] - * Le_k_0 - + this->shape_value_component (dof, quadrature_point, 2) - * lo_k_0); - this->restriction[index][2 * (i + 2 * l) + m] - ((2 * (i + 2) * this->degree + j) * deg + k + n_edge_dofs, dof) - += L_j_0 * (this->shape_grad_component (dof, quadrature_point, 2)[0] - * Le_k_1 - + this->shape_value_component (dof, quadrature_point, 2) - * lo_k_1); - this->restriction[index][2 * (i + 2 * l) + m] - (((2 * i + 5) * deg + k) * this->degree + j + n_edge_dofs, dof) - += L_j_1 * (this->shape_grad_component (dof, quadrature_point, 0)[2] - * Le_k_0 - + this->shape_value_component (dof, quadrature_point, 0) - * lo_k_0); - this->restriction[index][2 * (2 * i + l) + m] - ((2 * (i + 4) * this->degree + j) * deg + k + n_edge_dofs, dof) - += L_j_0 * (this->shape_grad_component (dof, quadrature_point, 0)[1] - * Le_k_1 - + this->shape_value_component (dof, quadrature_point, 0) - * lo_k_1); - this->restriction[index][2 * (2 * i + l) + m] - (((2 * i + 9) * deg + k) * this->degree + j + n_edge_dofs, dof) - += L_j_1 * (this->shape_grad_component (dof, quadrature_point, 1)[0] - * Le_k_0 - + this->shape_value_component (dof, quadrature_point, 1) - * lo_k_0); - } + const double L_k + = legendre_polynomials[k + 1].value + (edge_quadrature_points[q_point] (0)); + + for (unsigned int l = 0; l < tmp.size (); ++l) + system_rhs (k, l) += tmp (l) * L_k; } } + + system_matrix_inv.mmult (solution, system_rhs); + + for (unsigned int k = 0; k < 2; ++k) + for (unsigned int l = 0; l < deg; ++l) + { + if (std::abs (solution (l, k)) > 1e-14) + this->restriction[index][i + 2 * (2 * j + k)] + ((i + 4 * j) * this->degree + l + 1, dof) + = solution (l, k); + + if (std::abs (solution (l, k + 2)) > 1e-14) + this->restriction[index][2 * (i + 2 * j) + k] + ((i + 4 * j + 2) * this->degree + l + 1, dof) + = solution (l, k + 2); + + if (std::abs (solution (l, k + 4)) > 1e-14) + this->restriction[index][i + 2 * (j + 2 * k)] + ((i + 2 * (j + 4)) * this->degree + l + 1, + dof) + = solution (l, k + 4); + } } - - // Then, the interpolation to - // the interior shape functions. - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - { - Point quadrature_point = 2.0 * quadrature_points[q_point]; - unsigned int l; - unsigned int m; - unsigned int n; - - if (quadrature_points[q_point] (0) < 0.5) - { - n = 0; - - if (quadrature_points[q_point] (1) < 0.5) - { - m = 0; - - if (quadrature_points[q_point] (2) < 0.5) - l = 0; - - else - { - quadrature_point (2) -= 1.0; - l = 1; - } - } - - else - { - quadrature_point (1) -= 1.0; - m = 1; - - if (quadrature_points[q_point] (2) < 0.5) - l = 0; - - else - { - quadrature_point (2) -= 1.0; - l = 1; - } - } - } - - else - { - quadrature_point (0) -= 1.0; - n = 1; - - if (quadrature_points[q_point] (1) < 0.5) + + const QGauss<2> face_quadrature (2 * this->degree); + const std::vector >& face_quadrature_points + = face_quadrature.get_points (); + const std::vector >& + lobatto_polynomials + = Polynomials::Lobatto::generate_complete_basis + (this->degree); + const unsigned int n_edge_dofs + = GeometryInfo::lines_per_cell * this->degree; + const unsigned int& n_face_quadrature_points + = face_quadrature.size (); + + { + FullMatrix assembling_matrix + (deg * this->degree, + n_face_quadrature_points); + + for (unsigned int q_point = 0; + q_point < n_face_quadrature_points; ++q_point) + { + const double weight + = std::sqrt (face_quadrature.weight (q_point)); + + for (unsigned int i = 0; i <= deg; ++i) + { + const double L_i = weight + * legendre_polynomials[i].value + (face_quadrature_points[q_point] (0)); + + for (unsigned int j = 0; j < deg; ++j) + assembling_matrix (i * deg + j, q_point) + = L_i * lobatto_polynomials[j + 2].value + (face_quadrature_points[q_point] (1)); + } + } + + FullMatrix system_matrix (assembling_matrix.m (), + assembling_matrix.m ()); + + assembling_matrix.mTmult (system_matrix, + assembling_matrix); + system_matrix_inv.reinit (system_matrix.m (), + system_matrix.m ()); + system_matrix_inv.invert (system_matrix); + } + + solution.reinit (system_matrix_inv.m (), 24); + system_rhs.reinit (system_matrix_inv.m (), 24); + tmp.reinit (24); + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) + { + system_rhs = 0.0; + + for (unsigned int q_point = 0; + q_point < n_face_quadrature_points; ++q_point) + { + tmp = 0.0; + + if (face_quadrature_points[q_point] (0) < 0.5) + { + if (face_quadrature_points[q_point] (1) < 0.5) + { + Point quadrature_point_0 (i, + 2.0 * face_quadrature_points[q_point] (0), + 2.0 * face_quadrature_points[q_point] (1)); + + tmp (0) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 1); + tmp (1) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 2); + quadrature_point_0 + = Point (2.0 * face_quadrature_points[q_point] (0), + i, + 2.0 * face_quadrature_points[q_point] (1)); + tmp (8) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 2); + tmp (9) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 0); + quadrature_point_0 + = Point (2.0 * face_quadrature_points[q_point] (0), + 2.0 * face_quadrature_points[q_point] (1), + i); + tmp (16) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 0); + tmp (17) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 1); + } + + else + { + Point quadrature_point_0 (i, + 2.0 * face_quadrature_points[q_point] (0), + 2.0 * face_quadrature_points[q_point] (1) + - 1.0); + + tmp (2) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 1); + tmp (3) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 2); + quadrature_point_0 + = Point (2.0 * face_quadrature_points[q_point] (0), + i, + 2.0 * face_quadrature_points[q_point] (1) + - 1.0); + tmp (10) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 2); + tmp (11) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 0); + quadrature_point_0 + = Point (2.0 * face_quadrature_points[q_point] (0), + 2.0 * face_quadrature_points[q_point] (1) + - 1.0, i); + tmp (18) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 0); + tmp (19) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 1); + } + } + + else + if (face_quadrature_points[q_point] (1) < 0.5) { - m = 0; - - if (quadrature_points[q_point] (2) < 0.5) - l = 0; - - else - { - quadrature_point (2) -= 1.0; - l = 1; - } + Point quadrature_point_0 (i, + 2.0 * face_quadrature_points[q_point] (0) + - 1.0, + 2.0 * face_quadrature_points[q_point] (1)); + + tmp (4) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 1); + tmp (5) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 2); + quadrature_point_0 + = Point (2.0 * face_quadrature_points[q_point] (0) + - 1.0, i, + 2.0 * face_quadrature_points[q_point] (1)); + tmp (12) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 2); + tmp (13) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 0); + quadrature_point_0 + = Point (2.0 * face_quadrature_points[q_point] (0) + - 1.0, + 2.0 * face_quadrature_points[q_point] (1), + i); + tmp (20) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 0); + tmp (21) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 1); } - + else { - quadrature_point (1) -= 1.0; - m = 1; - + Point quadrature_point_0 (i, + 2.0 * face_quadrature_points[q_point] (0) + - 1.0, + 2.0 * face_quadrature_points[q_point] (1) + - 1.0); + + tmp (6) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 1); + tmp (7) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 2); + quadrature_point_0 + = Point (2.0 * face_quadrature_points[q_point] (0) + - 1.0, i, + 2.0 * face_quadrature_points[q_point] (1) + - 1.0); + tmp (14) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 2); + tmp (15) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 0); + quadrature_point_0 + = Point (2.0 * face_quadrature_points[q_point] (0) + - 1.0, + 2.0 * face_quadrature_points[q_point] (1) + - 1.0, i); + tmp (22) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 0); + tmp (23) += 2.0 * this->shape_value_component + (dof, quadrature_point_0, 1); + } + + const Point quadrature_point_0 (i, + face_quadrature_points[q_point] (0), + face_quadrature_points[q_point] (1)); + const Point quadrature_point_1 + (face_quadrature_points[q_point] (0), + i, + face_quadrature_points[q_point] (1)); + const Point quadrature_point_2 + (face_quadrature_points[q_point] (0), + face_quadrature_points[q_point] (1), + i); + + for (unsigned int j = 0; j < 2; ++j) + for (unsigned int k = 0; k < 2; ++k) + for (unsigned int l = 0; l <= deg; ++l) + { + tmp (2 * (j + 2 * k)) + -= this->restriction[index][i + 2 * (2 * j + k)] + ((i + 4 * j) * this->degree + l, dof) + * this->shape_value_component + ((i + 4 * j) * this->degree + l, + quadrature_point_0, 1); + tmp (2 * (j + 2 * k) + 1) + -= this->restriction[index][i + 2 * (2 * j + k)] + ((i + 2 * (k + 4)) * this->degree + l, + dof) + * this->shape_value_component + ((i + 2 * (k + 4)) * this->degree + l, + quadrature_point_0, 2); + tmp (2 * (j + 2 * (k + 2))) + -= this->restriction[index][2 * (i + 2 * j) + k] + ((2 * (i + 4) + k) * this->degree + l, + dof) + * this->shape_value_component + ((2 * (i + 4) + k) * this->degree + l, + quadrature_point_1, 2); + tmp (2 * (j + 2 * k) + 9) + -= this->restriction[index][2 * (i + 2 * j) + k] + ((i + 4 * j + 2) * this->degree + l, + dof) + * this->shape_value_component + ((i + 4 * j + 2) * this->degree + l, + quadrature_point_1, 0); + tmp (2 * (j + 2 * (k + 4))) + -= this->restriction[index][2 * (2 * i + j) + k] + ((4 * i + j + 2) * this->degree + l, + dof) + * this->shape_value_component + ((4 * i + j + 2) * this->degree + l, + quadrature_point_2, 0); + tmp (2 * (j + 2 * k) + 17) + -= this->restriction[index][2 * (2 * i + j) + k] + ((4 * i + k) * this->degree + l, dof) + * this->shape_value_component + ((4 * i + k) * this->degree + l, + quadrature_point_2, 1); + } + + tmp *= face_quadrature.weight (q_point); + + for (unsigned int j = 0; j <= deg; ++j) + { + const double L_j_0 + = legendre_polynomials[j].value + (face_quadrature_points[q_point] (0)); + const double L_j_1 + = legendre_polynomials[j].value + (face_quadrature_points[q_point] (1)); + + for (unsigned int k = 0; k < deg; ++k) + { + const double l_k_0 + = L_j_0 * lobatto_polynomials[k + 2].value + (face_quadrature_points[q_point] (1)); + const double l_k_1 + = L_j_1 * lobatto_polynomials[k + 2].value + (face_quadrature_points[q_point] (0)); + + for (unsigned int l = 0; l < 4; ++l) + { + system_rhs (j * deg + k, 2 * l) + += tmp (2 * l) * l_k_0; + system_rhs (j * deg + k, 2 * l + 1) + += tmp (2 * l + 1) * l_k_1; + system_rhs (j * deg + k, 2 * (l + 4)) + += tmp (2 * (l + 4)) * l_k_1; + system_rhs (j * deg + k, 2 * l + 9) + += tmp (2 * l + 9) * l_k_0; + system_rhs (j * deg + k, 2 * (l + 8)) + += tmp (2 * (l + 8)) * l_k_0; + system_rhs (j * deg + k, 2 * l + 17) + += tmp (2 * l + 17) * l_k_1; + } + } + } + } + + system_matrix_inv.mmult (solution, system_rhs); + + for (unsigned int j = 0; j < 2; ++j) + for (unsigned int k = 0; k < 2; ++k) + for (unsigned int l = 0; l <= deg; ++l) + for (unsigned int m = 0; m < deg; ++m) + { + if (std::abs (solution (l * deg + m, + 2 * (j + 2 * k))) + > 1e-14) + this->restriction[index][i + 2 * (2 * j + k)] + ((2 * i * this->degree + l) * deg + m + + n_edge_dofs, + dof) = solution (l * deg + m, + 2 * (j + 2 * k)); + + if (std::abs (solution (l * deg + m, + 2 * (j + 2 * k) + 1)) + > 1e-14) + this->restriction[index][i + 2 * (2 * j + k)] + (((2 * i + 1) * deg + m) * this->degree + l + + n_edge_dofs, dof) + = solution (l * deg + m, + 2 * (j + 2 * k) + 1); + + if (std::abs (solution (l * deg + m, + 2 * (j + 2 * (k + 2)))) + > 1e-14) + this->restriction[index][2 * (i + 2 * j) + k] + ((2 * (i + 2) * this->degree + l) * deg + m + + n_edge_dofs, + dof) = solution (l * deg + m, + 2 * (j + 2 * (k + 2))); + + if (std::abs (solution (l * deg + m, + 2 * (j + 2 * k) + 9)) + > 1e-14) + this->restriction[index][2 * (i + 2 * j) + k] + (((2 * i + 5) * deg + m) * this->degree + l + + n_edge_dofs, dof) + = solution (l * deg + m, + 2 * (j + 2 * k) + 9); + + if (std::abs (solution (l * deg + m, + 2 * (j + 2 * (k + 4)))) + > 1e-14) + this->restriction[index][2 * (2 * i + j) + k] + ((2 * (i + 4) * this->degree + l) * deg + m + + n_edge_dofs, + dof) = solution (l * deg + m, + 2 * (j + 2 * (k + 4))); + + if (std::abs (solution (l * deg + m, + 2 * (j + 2 * k) + 17)) + > 1e-14) + this->restriction[index][2 * (2 * i + j) + k] + (((2 * i + 9) * deg + m) * this->degree + l + + n_edge_dofs, dof) + = solution (l * deg + m, + 2 * (j + 2 * k) + 17); + } + } + + const QGauss quadrature (2 * this->degree); + const std::vector >& + quadrature_points = quadrature.get_points (); + const unsigned int n_boundary_dofs + = 2 * GeometryInfo::faces_per_cell * deg * this->degree + + n_edge_dofs; + const unsigned int& n_quadrature_points = quadrature.size (); + + { + FullMatrix + assembling_matrix (deg * deg * this->degree, + n_quadrature_points); + + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + { + const double weight = std::sqrt (quadrature.weight + (q_point)); + + for (unsigned int i = 0; i <= deg; ++i) + { + const double L_i = weight + * legendre_polynomials[i].value + (quadrature_points[q_point] (0)); + + for (unsigned int j = 0; j < deg; ++j) + { + const double l_j + = L_i * lobatto_polynomials[j + 2].value + (quadrature_points[q_point] (1)); + + for (unsigned int k = 0; k < deg; ++k) + assembling_matrix ((i * deg + j) * deg + k, + q_point) + = l_j * lobatto_polynomials[k + 2].value + (quadrature_points[q_point] (2)); + } + } + } + + FullMatrix system_matrix (assembling_matrix.m (), + assembling_matrix.m ()); + + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_matrix_inv.reinit (system_matrix.m (), + system_matrix.m ()); + system_matrix_inv.invert (system_matrix); + } + + solution.reinit (system_matrix_inv.m (), 24); + system_rhs.reinit (system_matrix_inv.m (), 24); + tmp.reinit (24); + + for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) + { + system_rhs = 0.0; + + for (unsigned int q_point = 0; + q_point < n_quadrature_points; ++q_point) + { + tmp = 0.0; + + if (quadrature_points[q_point] (0) < 0.5) + { + if (quadrature_points[q_point] (1) < 0.5) + { if (quadrature_points[q_point] (2) < 0.5) - l = 0; - + { + const Point quadrature_point + (2.0 * quadrature_points[q_point] (0), + 2.0 * quadrature_points[q_point] (1), + 2.0 * quadrature_points[q_point] (2)); + + tmp (0) += 2.0 * this->shape_value_component + (dof, quadrature_point, 0); + tmp (1) += 2.0 * this->shape_value_component + (dof, quadrature_point, 1); + tmp (2) += 2.0 * this->shape_value_component + (dof, quadrature_point, 2); + } + else { - quadrature_point (2) -= 1.0; - l = 1; + const Point quadrature_point + (2.0 * quadrature_points[q_point] (0), + 2.0 * quadrature_points[q_point] (1), + 2.0 * quadrature_points[q_point] (2) + - 1.0); + + tmp (3) += 2.0 * this->shape_value_component + (dof, quadrature_point, 0); + tmp (4) += 2.0 * this->shape_value_component + (dof, quadrature_point, 1); + tmp (5) += 2.0 * this->shape_value_component + (dof, quadrature_point, 2); } } + + else + if (quadrature_points[q_point] (2) < 0.5) + { + const Point quadrature_point + (2.0 * quadrature_points[q_point] (0), + 2.0 * quadrature_points[q_point] (1) + - 1.0, + 2.0 * quadrature_points[q_point] (2)); + + tmp (6) += 2.0 * this->shape_value_component + (dof, quadrature_point, 0); + tmp (7) += 2.0 * this->shape_value_component + (dof, quadrature_point, 1); + tmp (8) += 2.0 * this->shape_value_component + (dof, quadrature_point, 2); + } + + else + { + const Point quadrature_point + (2.0 * quadrature_points[q_point] (0), + 2.0 * quadrature_points[q_point] (1) + - 1.0, + 2.0 * quadrature_points[q_point] (2) + - 1.0); + + tmp (9) += 2.0 * this->shape_value_component + (dof, quadrature_point, 0); + tmp (10) += 2.0 * this->shape_value_component + (dof, quadrature_point, 1); + tmp (11) += 2.0 * this->shape_value_component + (dof, quadrature_point, 2); + } } - - const double weight = 2.0 * quadrature.weight (q_point); - - for (unsigned int i = 0; i < this->degree; ++i) + + else + if (quadrature_points[q_point] (1) < 0.5) + { + if (quadrature_points[q_point] (2) < 0.5) + { + const Point quadrature_point + (2.0 * quadrature_points[q_point] (0) + - 1.0, + 2.0 * quadrature_points[q_point] (1), + 2.0 * quadrature_points[q_point] (2)); + + tmp (12) += 2.0 * this->shape_value_component + (dof, quadrature_point, 0); + tmp (13) += 2.0 * this->shape_value_component + (dof, quadrature_point, 1); + tmp (14) += 2.0 * this->shape_value_component + (dof, quadrature_point, 2); + } + + else + { + const Point quadrature_point + (2.0 * quadrature_points[q_point] (0) + - 1.0, + 2.0 * quadrature_points[q_point] (1), + 2.0 * quadrature_points[q_point] (2) + - 1.0); + + tmp (15) += 2.0 * this->shape_value_component + (dof, quadrature_point, 0); + tmp (16) += 2.0 * this->shape_value_component + (dof, quadrature_point, 1); + tmp (17) += 2.0 * this->shape_value_component + (dof, quadrature_point, 2); + } + } + + else + if (quadrature_points[q_point] (2) < 0.5) + { + const Point quadrature_point + (2.0 * quadrature_points[q_point] (0) + - 1.0, + 2.0 * quadrature_points[q_point] (1) + - 1.0, + 2.0 * quadrature_points[q_point] (2)); + + tmp (18) += 2.0 * this->shape_value_component + (dof, quadrature_point, 0); + tmp (19) += 2.0 * this->shape_value_component + (dof, quadrature_point, 1); + tmp (20) += 2.0 * this->shape_value_component + (dof, quadrature_point, 2); + } + + else + { + const Point quadrature_point + (2.0 * quadrature_points[q_point] (0) + - 1.0, + 2.0 * quadrature_points[q_point] (1) + - 1.0, + 2.0 * quadrature_points[q_point] (2) + - 1.0); + + tmp (21) += 2.0 * this->shape_value_component + (dof, quadrature_point, 0); + tmp (22) += 2.0 * this->shape_value_component + (dof, quadrature_point, 1); + tmp (23) += 2.0 * this->shape_value_component + (dof, quadrature_point, 2); + } + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < 2; ++j) + for (unsigned int k = 0; k < 2; ++k) + for (unsigned int l = 0; l <= deg; ++l) + { + tmp (3 * (i + 2 * (j + 2 * k))) + -= this->restriction[index][2 * (2 * i + j) + k] + ((4 * i + j + 2) * this->degree + l, dof) + * this->shape_value_component + ((4 * i + j + 2) * this->degree + l, + quadrature_points[q_point], 0); + tmp (3 * (i + 2 * (j + 2 * k)) + 1) + -= this->restriction[index][2 * (2 * i + j) + k] + ((4 * i + k) * this->degree + l, dof) + * this->shape_value_component + ((4 * i + k) * this->degree + l, + quadrature_points[q_point], 1); + tmp (3 * (i + 2 * (j + 2 * k)) + 2) + -= this->restriction[index][2 * (2 * i + j) + k] + ((2 * (j + 4) + k) * this->degree + l, + dof) + * this->shape_value_component + ((2 * (j + 4) + k) * this->degree + l, + quadrature_points[q_point], 2); + + for (unsigned int m = 0; m < deg; ++m) + { + tmp (3 * (i + 2 * (j + 2 * k))) + -= this->restriction[index][2 * (2 * i + j) + k] + (((2 * j + 5) * deg + m) + * this->degree + l + n_edge_dofs, + dof) + * this->shape_value_component + (((2 * j + 5) * deg + m) + * this->degree + l + n_edge_dofs, + quadrature_points[q_point], 0); + tmp (3 * (i + 2 * (j + 2 * k))) + -= this->restriction[index][2 * (2 * i + j) + k] + ((2 * (i + 4) * this->degree + l) + * deg + m + n_edge_dofs, dof) + * this->shape_value_component + ((2 * (i + 4) * this->degree + l) + * deg + m + n_edge_dofs, + quadrature_points[q_point], 0); + tmp (3 * (i + 2 * (j + 2 * k)) + 1) + -= this->restriction[index][2 * (2 * i + j) + k] + ((2 * k * this->degree + l) * deg + m + + n_edge_dofs, + dof) + * this->shape_value_component + ((2 * k * this->degree + l) * deg + m + + n_edge_dofs, + quadrature_points[q_point], 1); + tmp (3 * (i + 2 * (j + 2 * k)) + 1) + -= this->restriction[index][2 * (2 * i + j) + k] + (((2 * i + 9) * deg + m) + * this->degree + l + n_edge_dofs, + dof) + * this->shape_value_component + (((2 * i + 9) * deg + m) + * this->degree + l + n_edge_dofs, + quadrature_points[q_point], 1); + tmp (3 * (i + 2 * (j + 2 * k)) + 2) + -= this->restriction[index][2 * (2 * i + j) + k] + (((2 * k + 1) * deg + m) + * this->degree + l + n_edge_dofs, + dof) + * this->shape_value_component + (((2 * k + 1) * deg + m) + * this->degree + l + n_edge_dofs, + quadrature_points[q_point], 2); + tmp (3 * (i + 2 * (j + 2 * k)) + 2) + -= this->restriction[index][2 * (2 * i + j) + k] + ((2 * (j + 2) * this->degree + l) + * deg + m + n_edge_dofs, dof) + * this->shape_value_component + ((2 * (j + 2) * this->degree + l) + * deg + m + n_edge_dofs, + quadrature_points[q_point], 2); + } + } + + tmp *= quadrature.weight (q_point); + + for (unsigned int i = 0; i <= deg; ++i) { const double L_i_0 - = weight * legendre_polynomials[i].value (quadrature_points[q_point] (0)); + = legendre_polynomials[i].value + (quadrature_points[q_point] (0)); const double L_i_1 - = weight * legendre_polynomials[i].value (quadrature_points[q_point] (1)); + = legendre_polynomials[i].value + (quadrature_points[q_point] (1)); const double L_i_2 - = weight * legendre_polynomials[i].value (quadrature_points[q_point] (2)); - + = legendre_polynomials[i].value + (quadrature_points[q_point] (2)); + for (unsigned int j = 0; j < deg; ++j) { - const double Le_j_0 - = legendre_polynomials[j + 1].value (quadrature_points[q_point] (0)); - const double Le_j_1 - = legendre_polynomials[j + 1].value (quadrature_points[q_point] (1)); - const double lo_j_0 - = lobatto_polynomials[j + 2].value (quadrature_points[q_point] (0)); - const double lo_j_1 - = lobatto_polynomials[j + 2].value (quadrature_points[q_point] (1)); - + const double l_j_0 + = L_i_0 * lobatto_polynomials[j + 2].value + (quadrature_points[q_point] (1)); + const double l_j_1 + = L_i_1 * lobatto_polynomials[j + 2].value + (quadrature_points[q_point] (0)); + const double l_j_2 + = L_i_2 * lobatto_polynomials[j + 2].value + (quadrature_points[q_point] (0)); + for (unsigned int k = 0; k < deg; ++k) { - const double L_k - = legendre_polynomials[k + 1].value (quadrature_points[q_point] (2)); + const double l_k_0 + = l_j_0 * lobatto_polynomials[k + 2].value + (quadrature_points[q_point] (2)); const double l_k_1 - = lobatto_polynomials[k + 2].value (quadrature_points[q_point] (1)); + = l_j_1 * lobatto_polynomials[k + 2].value + (quadrature_points[q_point] (2)); const double l_k_2 - = lobatto_polynomials[k + 2].value (quadrature_points[q_point] (2)); - - this->restriction[index][2 * (2 * l + m) + n] - ((i * deg + j) * deg + k + n_boundary_dofs, dof) - += L_i_0 * (this->shape_grad_component (dof, quadrature_point, 0)[2] - * lo_j_1 * L_k - + this->shape_grad_component (dof, quadrature_point, 0)[1] - * Le_j_1 * l_k_2 - + this->shape_value_component (dof, quadrature_point, 0) - * lo_j_1 * l_k_2); - this->restriction[index][2 * (2 * l + m) + n] - ((i + (j + deg) * this->degree) * deg + k + n_boundary_dofs, dof) - += L_i_1 * (this->shape_grad_component (dof, quadrature_point, 1)[2] - * lo_j_0 * L_k - + this->shape_grad_component (dof, quadrature_point, 1)[0] - * Le_j_0 * l_k_2 - + this->shape_value_component (dof, quadrature_point, 1) - * lo_j_0 * l_k_2); - this->restriction[index][2 * (2 * l + m) + n] - (i + (k + (j + 2 * deg) * deg) * this->degree + n_boundary_dofs, dof) - += L_i_2 * (this->shape_grad_component (dof, quadrature_point, 2)[1] - * lo_j_0 - * legendre_polynomials[k + 1].value (quadrature_points[q_point] (1)) - + this->shape_grad_component (dof, quadrature_point, 2)[0] - * Le_j_0 * l_k_1 - + this->shape_value_component (dof, quadrature_point, 2) - * lo_j_0 * l_k_1); + = l_j_2 * lobatto_polynomials[k + 2].value + (quadrature_points[q_point] (1)); + + for (unsigned int l = 0; l < 8; ++l) + { + system_rhs ((i * deg + j) * deg + k, + 3 * l) + += tmp (3 * l) * l_k_0; + system_rhs ((i * deg + j) * deg + k, + 3 * l + 1) + += tmp (3 * l + 1) * l_k_1; + system_rhs ((i * deg + j) * deg + k, + 3 * l + 2) + += tmp (3 * l + 2) * l_k_2; + } } } } } + + system_matrix_inv.mmult (solution, system_rhs); + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < 2; ++j) + for (unsigned int k = 0; k < 2; ++k) + for (unsigned int l = 0; l <= deg; ++l) + for (unsigned int m = 0; m < deg; ++m) + for (unsigned int n = 0; n < deg; ++n) + { + if (std::abs (solution + ((l * deg + m) * deg + n, + 3 * (i + 2 * (j + 2 * k)))) + > 1e-14) + this->restriction[index][2 * (2 * i + j) + k] + ((l * deg + m) * deg + n + n_boundary_dofs, + dof) = solution ((l * deg + m) * deg + n, + 3 * (i + 2 * (j + 2 * k))); + + if (std::abs (solution + ((l * deg + m) * deg + n, + 3 * (i + 2 * (j + 2 * k)) + 1)) + > 1e-14) + this->restriction[index][2 * (2 * i + j) + k] + ((l + (m + deg) * this->degree) * deg + n + + n_boundary_dofs, + dof) = solution ((l * deg + m) * deg + n, + 3 * (i + 2 * (j + 2 * k)) + 1); + + if (std::abs (solution + ((l * deg + m) * deg + n, + 3 * (i + 2 * (j + 2 * k)) + 2)) + > 1e-14) + this->restriction[index][2 * (2 * i + j) + k] + (l + ((m + 2 * deg) * deg + n) * this->degree + + n_boundary_dofs, dof) + = solution ((l * deg + m) * deg + n, + 3 * (i + 2 * (j + 2 * k)) + 2); + } } } @@ -1057,6 +1968,7 @@ FE_Nedelec::initialize_restriction () } + template std::vector FE_Nedelec::get_dpo_vector (const unsigned int degree, bool dg) @@ -1080,7 +1992,6 @@ FE_Nedelec::get_dpo_vector (const unsigned int degree, bool dg) return dpo; } - //--------------------------------------------------------------------------- // Data field initialization //--------------------------------------------------------------------------- @@ -1596,9 +2507,8 @@ const + GeometryInfo::lines_per_face) = 1.0; interpolation_matrix - (i + (j + GeometryInfo::lines_per_face + source_fe.degree - 1) - * source_fe.degree, - i + (j + GeometryInfo::lines_per_face + this->degree-1) * this->degree) + (i + (j + GeometryInfo::lines_per_face) * source_fe.degree, + i + (j + GeometryInfo::lines_per_face) * this->degree) = 1.0; } } @@ -1620,6 +2530,20 @@ FE_Nedelec<1>::get_subface_interpolation_matrix( // In this function we compute the // subface interpolation matrix. + // This is done by a projection- + // based interpolation. Therefore + // we first interpolate the + // shape functions of the higher + // order element on the lowest + // order edge shape functions. + // Then the remaining part of + // the interpolated shape + // functions is projected on the + // higher order edge shape + // functions, the face shape + // functions and the interior + // shape functions (if they all + // exist). template void FE_Nedelec::get_subface_interpolation_matrix( @@ -1627,9 +2551,9 @@ FE_Nedelec::get_subface_interpolation_matrix( const unsigned int subface, FullMatrix& interpolation_matrix) const { - // this is only implemented, if the - // source FE is also a - // Nedelec element + // this is only implemented, if the + // source FE is also a + // Nedelec element typedef FE_Nedelec FEN; typedef FiniteElement FEL; @@ -1643,12 +2567,12 @@ FE_Nedelec::get_subface_interpolation_matrix( ExcDimensionMismatch (interpolation_matrix.n (), this->dofs_per_face)); - // ok, source is a Nedelec element, so - // we will be able to do the work + // ok, source is a Nedelec element, so + // we will be able to do the work const FE_Nedelec &source_fe = dynamic_cast&> (source); - // Make sure, that the element, + // Make sure, that the element, // for which the DoFs should be // constrained is the one with // the higher polynomial degree. @@ -1656,41 +2580,101 @@ FE_Nedelec::get_subface_interpolation_matrix( // also if this assertion is not // satisfied. But the matrices // produced in that case might - // lead to problems in the + // lead to problems in the // hp procedures, which use this - // method. + // method. Assert (this->dofs_per_face <= source_fe.dofs_per_face, typename FEL::ExcInterpolationNotImplemented ()); interpolation_matrix = 0.0; - // Set the degrees of freedom - // by interpolation as usual. + // Perform projection-based interpolation + // as usual. const QGauss<1> edge_quadrature (source_fe.degree); const std::vector >& edge_quadrature_points = edge_quadrature.get_points (); - const std::vector >& - legendre_polynomials - = Polynomials::Legendre::generate_complete_basis (source_fe.degree - 1); const unsigned int& n_edge_quadrature_points = edge_quadrature.size (); switch (dim) { case 2: { - for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; - ++q_point) + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) + for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; + ++q_point) + { + const Point quadrature_point (0.0, + 0.5 * (edge_quadrature_points[q_point] (0) + + subface)); + + interpolation_matrix (0, dof) += 0.5 + * edge_quadrature.weight (q_point) + * this->shape_value_component + (dof, quadrature_point, 1); + } + + if (source_fe.degree > 1) { - const Point quadrature_point (0.0, - 0.5 * (edge_quadrature_points[q_point] (0) - + subface)); - const double weight = 0.5 * edge_quadrature.weight (q_point); - + const std::vector >& + legendre_polynomials + = Polynomials::Legendre::generate_complete_basis (source_fe.degree - 1); + FullMatrix system_matrix_inv (source_fe.degree - 1, + source_fe.degree - 1); + + { + FullMatrix assembling_matrix (source_fe.degree - 1, + n_edge_quadrature_points); + + for (unsigned int q_point = 0; + q_point < n_edge_quadrature_points; ++q_point) + { + const double weight + = std::sqrt (edge_quadrature.weight (q_point)); + + for (unsigned int i = 0; i < source_fe.degree - 1; ++i) + assembling_matrix (i, q_point) = weight + * legendre_polynomials[i + 1].value + (edge_quadrature_points[q_point] (0)); + } + + FullMatrix system_matrix (source_fe.degree - 1, source_fe.degree - 1); + + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_matrix_inv.invert (system_matrix); + } + + Vector solution (source_fe.degree - 1); + Vector system_rhs (source_fe.degree - 1); + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - const double shape_value = weight * this->shape_value_component (dof, quadrature_point, 1); - - for (unsigned int i = 0; i < source_fe.dofs_per_face; ++i) - interpolation_matrix (i, dof) += shape_value - * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0)); + system_rhs = 0.0; + + for (unsigned int q_point = 0; + q_point < n_edge_quadrature_points; ++q_point) + { + const Point quadrature_point_0 (0.0, + 0.5 * (edge_quadrature_points[q_point] (0) + + subface)); + const Point quadrature_point_1 (0.0, + edge_quadrature_points[q_point] (0)); + const double tmp = edge_quadrature.weight (q_point) + * (0.5 * this->shape_value_component + (dof, quadrature_point_0, 1) + - interpolation_matrix (0, + dof) + * source_fe.shape_value_component + (0, quadrature_point_1, 1)); + + for (unsigned int i = 0; i < source_fe.degree - 1; ++i) + system_rhs (i) += tmp + * legendre_polynomials[i + 1].value + (edge_quadrature_points[q_point] (0)); + } + + system_matrix_inv.vmult (solution, system_rhs); + + for (unsigned int i = 0; i < source_fe.degree - 1; ++i) + if (std::abs (solution (i)) > 1e-14) + interpolation_matrix (i + 1, dof) = solution (i); } } @@ -1702,81 +2686,169 @@ FE_Nedelec::get_subface_interpolation_matrix( const double shifts[4][2] = { { 0.0, 0.0 }, { 1.0, 0.0 }, { 0.0, 1.0 }, { 1.0, 1.0 } }; - for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; - ++q_point) - { - const double weight = 0.5 * edge_quadrature.weight (q_point); + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) + for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; + ++q_point) + { + const double weight = 0.5 * edge_quadrature.weight (q_point); - for (unsigned int i = 0; i < 2; ++i) - { - Point - quadrature_point (0.5 * (i + shifts[subface][0]), - 0.5 * (edge_quadrature_points[q_point] (0) - + shifts[subface][1]), - 0.0); - - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) - { - const double shape_value - = weight * this->shape_value_component (this->face_to_cell_index (dof, 4), - quadrature_point, 1); - - for (unsigned int j = 0; j < source_fe.degree; ++j) - interpolation_matrix (i * source_fe.degree + j, dof) - += shape_value - * legendre_polynomials[j].value (edge_quadrature_points[q_point] (0)); - } - - quadrature_point - = Point (0.5 * (edge_quadrature_points[q_point] (0) - + shifts[subface][0]), - 0.5 * (i + shifts[subface][1]), 0.0); - - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) - { - const double shape_value - = weight * this->shape_value_component (this->face_to_cell_index (dof, 4), - quadrature_point, 0); - - for (unsigned int j = 0; j < source_fe.degree; ++j) - interpolation_matrix ((i + 2) * source_fe.degree + j, dof) - += shape_value - * legendre_polynomials[j].value (edge_quadrature_points[q_point] (0)); - } - } - } + for (unsigned int i = 0; i < 2; ++i) + { + Point + quadrature_point (0.5 * (i + shifts[subface][0]), + 0.5 * (edge_quadrature_points[q_point] (0) + + shifts[subface][1]), + 0.0); + + interpolation_matrix (i * source_fe.degree, dof) += weight + * this->shape_value_component + (this->face_to_cell_index (dof, 4), + quadrature_point, + 1); + quadrature_point + = Point (0.5 * (edge_quadrature_points[q_point] (0) + + shifts[subface][0]), + 0.5 * (i + shifts[subface][1]), 0.0); + interpolation_matrix ((i + 2) * source_fe.degree, dof) + += weight * this->shape_value_component + (this->face_to_cell_index (dof, 4), + quadrature_point, 0); + } + } if (source_fe.degree > 1) { - const QGauss<2> quadrature (source_fe.degree); - const std::vector >& - quadrature_points = quadrature.get_points (); const std::vector >& - lobatto_polynomials - = Polynomials::Lobatto::generate_complete_basis - (source_fe.degree); - const unsigned int n_boundary_dofs - = GeometryInfo::lines_per_face * source_fe.degree; - const unsigned int& n_quadrature_points = quadrature.size (); - FullMatrix - system_matrix_inv (source_fe.degree * (source_fe.degree - 1), - source_fe.degree * (source_fe.degree - 1)); + legendre_polynomials + = Polynomials::Legendre::generate_complete_basis (source_fe.degree - 1); + FullMatrix system_matrix_inv (source_fe.degree - 1, + source_fe.degree - 1); { - FullMatrix - assembling_matrix (source_fe.degree * (source_fe.degree - 1), - n_quadrature_points); + FullMatrix assembling_matrix (source_fe.degree - 1, + n_edge_quadrature_points); - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) + for (unsigned int q_point = 0; + q_point < n_edge_quadrature_points; ++q_point) { - const double weight = std::sqrt (quadrature.weight (q_point)); + const double weight + = std::sqrt (edge_quadrature.weight (q_point)); - for (unsigned int i = 0; i < source_fe.degree; ++i) - { - const double L_i = weight - * legendre_polynomials[i].value - (quadrature_points[q_point] (0)); + for (unsigned int i = 0; i < source_fe.degree - 1; ++i) + assembling_matrix (i, q_point) = weight + * legendre_polynomials[i + 1].value + (edge_quadrature_points[q_point] (0)); + } + + FullMatrix system_matrix (source_fe.degree - 1, source_fe.degree - 1); + + assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_matrix_inv.invert (system_matrix); + } + + FullMatrix solution (source_fe.degree - 1, + GeometryInfo::lines_per_face); + FullMatrix system_rhs (source_fe.degree - 1, + GeometryInfo::lines_per_face); + Vector tmp (GeometryInfo::lines_per_face); + + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) + { + system_rhs = 0.0; + + for (unsigned int q_point = 0; + q_point < n_edge_quadrature_points; ++q_point) + { + const double weight = edge_quadrature.weight (q_point); + + for (unsigned int i = 0; i < 2; ++i) + { + Point + quadrature_point_0 + (0.5 * (i + shifts[subface][0]), + 0.5 * (edge_quadrature_points[q_point] (0) + + shifts[subface][1]), 0.0); + Point quadrature_point_1 (i, + edge_quadrature_points[q_point] (0), + 0.0); + + tmp (i) = weight + * (0.5 * this->shape_value_component + (this->face_to_cell_index (dof, 4), + quadrature_point_0, 1) + - interpolation_matrix + (i * source_fe.degree, dof) + * source_fe.shape_value_component + (i * source_fe.degree, + quadrature_point_1, 1)); + quadrature_point_0 + = Point (0.5 * (edge_quadrature_points[q_point] (0) + + shifts[subface][0]), + 0.5 * (i + shifts[subface][1]), + 0.0); + quadrature_point_1 + = Point (edge_quadrature_points[q_point] (0), + i, 0.0); + tmp (i + 2) = weight + * (0.5 * this->shape_value_component + (this->face_to_cell_index (dof, 4), + quadrature_point_0, 0) + - interpolation_matrix + ((i + 2) * source_fe.degree, + dof) + * source_fe.shape_value_component + ((i + 2) * source_fe.degree, + quadrature_point_1, 0)); + } + + for (unsigned int i = 0; i < source_fe.degree - 1; ++i) + { + const double L_i + = legendre_polynomials[i + 1].value + (edge_quadrature_points[q_point] (0)); + + for (unsigned int j = 0; + j < GeometryInfo::lines_per_face; ++j) + system_rhs (i, j) += tmp (j) * L_i; + } + } + + system_matrix_inv.mmult (solution, system_rhs); + + for (unsigned int i = 0; + i < GeometryInfo::lines_per_face; ++i) + for (unsigned int j = 0; j < source_fe.degree - 1; ++j) + if (std::abs (solution (j, i)) > 1e-14) + interpolation_matrix (i * source_fe.degree + j + 1, + dof) = solution (j, i); + } + + const QGauss<2> quadrature (source_fe.degree); + const std::vector >& + quadrature_points = quadrature.get_points (); + const std::vector >& + lobatto_polynomials + = Polynomials::Lobatto::generate_complete_basis + (source_fe.degree); + const unsigned int n_boundary_dofs + = GeometryInfo::lines_per_face * source_fe.degree; + const unsigned int& n_quadrature_points = quadrature.size (); + + { + FullMatrix + assembling_matrix (source_fe.degree * (source_fe.degree - 1), + n_quadrature_points); + + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + { + const double weight = std::sqrt (quadrature.weight (q_point)); + + for (unsigned int i = 0; i < source_fe.degree; ++i) + { + const double L_i = weight + * legendre_polynomials[i].value + (quadrature_points[q_point] (0)); for (unsigned int j = 0; j < source_fe.degree - 1; ++j) assembling_matrix (i * (source_fe.degree - 1) + j, @@ -1790,12 +2862,14 @@ FE_Nedelec::get_subface_interpolation_matrix( assembling_matrix.m ()); assembling_matrix.mTmult (system_matrix, assembling_matrix); + system_matrix_inv.reinit (system_matrix.m (), + system_matrix.m ()); system_matrix_inv.invert (system_matrix); } - FullMatrix solution (system_matrix_inv.m (), 2); - FullMatrix system_rhs (system_matrix_inv.m (), 2); - Vector tmp (2); + solution.reinit (system_matrix_inv.m (), 2); + system_rhs.reinit (system_matrix_inv.m (), 2); + tmp.reinit (2); for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { @@ -1901,21 +2975,16 @@ void FE_Nedelec::interpolate (std::vector&, const std::vector void FE_Nedelec::interpolate (std::vector& local_dofs, const std::vector >& values, unsigned int offset) const { + const unsigned int deg = this->degree-1; + Assert (values.size () == this->generalized_support_points.size (), ExcDimensionMismatch (values.size (), this->generalized_support_points.size ())); @@ -1927,872 +2996,938 @@ FE_Nedelec::interpolate (std::vector& local_dofs, std::fill (local_dofs.begin (), local_dofs.end (), 0.); if (offset < dim) - { - const unsigned int deg = this->degree-1; - const std::vector >& - legendre_polynomials - = Polynomials::Legendre::generate_complete_basis (deg); - const QGauss<1> edge_quadrature (this->degree); - const std::vector >& - edge_quadrature_points = edge_quadrature.get_points (); - const unsigned int& n_edge_points = edge_quadrature.size (); - - switch (dim) - { - case 2: - { + switch (dim) + { + case 2: + { + const QGauss<1> reference_edge_quadrature (this->degree); + const unsigned int& n_edge_points + = reference_edge_quadrature.size (); + // Let us begin with the - // edge degrees of freedom. - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) + // interpolation part. + for (unsigned int i = 0; i < 2; ++i) + { + for (unsigned int q_point = 0; q_point < n_edge_points; + ++q_point) + local_dofs[i * this->degree] + += reference_edge_quadrature.weight (q_point) + * values[q_point + i * n_edge_points] (1); + + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + if (std::abs (local_dofs[i * this->degree]) < 1e-14) + local_dofs[i * this->degree] = 0.0; + } + + if (offset == 0) + for (unsigned int i = 0; i < 2; ++i) { - const double weight = edge_quadrature.weight (q_point); - - for (unsigned int i = 0; i < this->degree; ++i) + for (unsigned int q_point = 0; q_point < n_edge_points; + ++q_point) + local_dofs[(i + 2) * this->degree] + += reference_edge_quadrature.weight (q_point) + * values[q_point + (i + 2) * n_edge_points] (0); + + if (std::abs (local_dofs[(i + 2) * this->degree]) < 1e-14) + local_dofs[(i + 2) * this->degree] = 0.0; + } + + // If the degree is greater + // than 0, then we have still + // some higher order edge + // shape functions to + // consider. + // Here the projection part + // starts. The dof values + // are obtained by solving + // a linear system of + // equations. + if (this->degree > 1) + { + // We start with projection + // on the higher order edge + // shape function. + const std::vector >& + lobatto_polynomials + = Polynomials::Lobatto::generate_complete_basis + (this->degree); + const unsigned int + line_coordinate[GeometryInfo<2>::lines_per_cell] + = {1, 1, 0, 0}; + std::vector > + lobatto_polynomials_grad (this->degree); + + for (unsigned int i = 0; i < lobatto_polynomials_grad.size (); + ++i) + lobatto_polynomials_grad[i] + = lobatto_polynomials[i + 1].derivative (); + + // Set up the system matrix. + // This can be used for all + // edges. + FullMatrix system_matrix (this->degree-1, this->degree-1); + + for (unsigned int i = 0; i < system_matrix.m (); ++i) + for (unsigned int j = 0; j < system_matrix.n (); ++j) + for (unsigned int q_point = 0; q_point < n_edge_points; + ++q_point) + system_matrix (i, j) + += boundary_weights (q_point, j) + * lobatto_polynomials_grad[i + 1].value + (this->generalized_face_support_points[q_point] + (1)); + + FullMatrix system_matrix_inv (this->degree-1, this->degree-1); + + system_matrix_inv.invert (system_matrix); + + Vector system_rhs (system_matrix.m ()); + Vector solution (system_rhs.size ()); + + for (unsigned int line = 0; + line < GeometryInfo::lines_per_cell; ++line) + if ((line < 2) || (offset == 0)) { - const double L_i - = weight * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0)); - - for (unsigned int j = 0; j < 2; ++j) + // Set up the right hand side. + system_rhs = 0; + + for (unsigned int q_point = 0; q_point < n_edge_points; + ++q_point) { - local_dofs[i + j * this->degree] - += L_i * values[q_point + j * n_edge_points] (1); - - if (offset == 0) - local_dofs[i + (j + 2) * this->degree] - += L_i * values[q_point + (j + 2) * n_edge_points] (0); + const double tmp + = values[line * n_edge_points + q_point] + (line_coordinate[line]) + - local_dofs[line * this->degree] + * this->shape_value_component + (line * this->degree, + this->generalized_support_points[line + * n_edge_points + + q_point], + line_coordinate[line]); + + for (unsigned int i = 0; i < system_rhs.size (); + ++i) + system_rhs (i) += boundary_weights (q_point, i) + * tmp; } - } - } - + + system_matrix_inv.vmult (solution, system_rhs); + // Add the computed values // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i < GeometryInfo::lines_per_cell * this->degree; ++i) - if (std::abs (local_dofs[i]) < 1e-14) - local_dofs[i] = 0.0; + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + local_dofs[line * this->degree + i + 1] + = solution (i); + } - // If the degree is greater - // than 0, then we have still - // some interpolations onto - // the interior shape - // functions left. - if (this->degree > 1) - { - const std::vector >& - lobatto_polynomials - = Polynomials::Lobatto::generate_complete_basis - (this->degree); - std::vector > - lobatto_polynomials_grad (this->degree); - - for (unsigned int i = 0; i < lobatto_polynomials_grad.size (); - ++i) - lobatto_polynomials_grad[i] - = lobatto_polynomials[i + 1].derivative (); - - // We set up the system - // matrix and use it for - // both, the horizontal - // and the vertical + // Then we go on to the // interior shape + // functions. Again we + // set up the system + // matrix and use it + // for both, the + // horizontal and the + // vertical, interior + // shape functions. + const QGauss reference_quadrature (this->degree); + const std::vector >& + legendre_polynomials + = Polynomials::Legendre::generate_complete_basis (this->degree-1); + const unsigned int& n_interior_points + = reference_quadrature.size (); + + system_matrix.reinit ((this->degree-1) * this->degree, + (this->degree-1) * this->degree); + system_matrix = 0; + + for (unsigned int i = 0; i < this->degree; ++i) + for (unsigned int j = 0; j < this->degree-1; ++j) + for (unsigned int k = 0; k < this->degree; ++k) + for (unsigned int l = 0; l < this->degree-1; ++l) + for (unsigned int q_point = 0; + q_point < n_interior_points; ++q_point) + system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l) + += reference_quadrature.weight (q_point) + * legendre_polynomials[i].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points] + (0)) + * lobatto_polynomials[j + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points] + (1)) + * lobatto_polynomials_grad[k].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points] + (0)) + * lobatto_polynomials[l + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points] + (1)); + + system_matrix_inv.reinit (system_matrix.m (), + system_matrix.m ()); + system_matrix_inv.invert (system_matrix); + solution.reinit (system_matrix_inv.m ()); + system_rhs.reinit (system_matrix.m ()); + + if (offset == 0) + { + // Set up the right hand side + // for the horizontal shape // functions. - const QGauss reference_quadrature (this->degree); - const unsigned int& n_interior_points - = reference_quadrature.size (); + system_rhs = 0; - FullMatrix system_matrix (deg * this->degree, - deg * this->degree); + for (unsigned int q_point = 0; + q_point < n_interior_points; ++q_point) + { + double tmp + = values[q_point + GeometryInfo::lines_per_cell + * n_edge_points] (0); - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - for (unsigned int k = 0; k < this->degree; ++k) - for (unsigned int l = 0; l < deg; ++l) - for (unsigned int q_point = 0; - q_point < n_interior_points; ++q_point) - system_matrix (i * deg + j, k * deg + l) - += reference_quadrature.weight (q_point) - * legendre_polynomials[i].value + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < this->degree; ++j) + tmp -= local_dofs[(i + 2) * this->degree + j] + * this->shape_value_component + ((i + 2) * this->degree + j, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points], + 0); + + for (unsigned int i = 0; i < this->degree; ++i) + for (unsigned int j = 0; j < this->degree-1; ++j) + system_rhs (i * (this->degree-1) + j) + += reference_quadrature.weight (q_point) * tmp + * lobatto_polynomials_grad[i].value (this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points] (0)) * lobatto_polynomials[j + 2].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points] - (1)) - * lobatto_polynomials_grad[k].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points] - (0)) - * lobatto_polynomials[l + 2].value (this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points] (1)); + } - FullMatrix system_matrix_inv (system_matrix.m (), - system_matrix.m ()); - - system_matrix_inv.invert (system_matrix); - - Vector solution (system_matrix_inv.m ()); - Vector system_rhs (system_matrix.m ()); - - if (offset == 0) - { - // Set up the right hand side - // for the horizontal shape - // functions. - for (unsigned int q_point = 0; - q_point < n_interior_points; ++q_point) - { - double tmp - = values[q_point + GeometryInfo::lines_per_cell - * n_edge_points] (0); - - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp -= local_dofs[(i + 2) * this->degree + j] - * this->shape_value_component - ((i + 2) * this->degree + j, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points], - 0); - - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < this->degree-1; ++j) - system_rhs (i * (this->degree-1) + j) - += reference_quadrature.weight (q_point) * tmp - * lobatto_polynomials_grad[i].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points] - (0)) - * lobatto_polynomials[j + 2].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points] - (1)); - } - - system_matrix_inv.vmult (solution, system_rhs); + system_matrix_inv.vmult (solution, system_rhs); // Add the computed values // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - if (std::abs (solution (i * deg + j)) > 1e-14) - local_dofs[(i + GeometryInfo::lines_per_cell) - * deg + j - + GeometryInfo::lines_per_cell] - = solution (i * deg + j); - } + for (unsigned int i = 0; i < this->degree; ++i) + for (unsigned int j = 0; j < this->degree-1; ++j) + if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14) + local_dofs[(i + GeometryInfo::lines_per_cell) + * (this->degree-1) + j + + GeometryInfo::lines_per_cell] + = solution (i * (this->degree-1) + j); + } // Set up the right hand side // for the vertical shape // functions. - system_rhs = 0; - - for (unsigned int q_point = 0; q_point < n_interior_points; - ++q_point) - { - double tmp - = values[q_point + GeometryInfo::lines_per_cell - * n_edge_points] (1); + system_rhs = 0; - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp -= local_dofs[i * this->degree + j] - * this->shape_value_component - (i * this->degree + j, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points], - 1); + for (unsigned int q_point = 0; q_point < n_interior_points; + ++q_point) + { + double tmp + = values[q_point + GeometryInfo::lines_per_cell + * n_edge_points] (1); - for (unsigned i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - system_rhs (i * deg + j) - += reference_quadrature.weight (q_point) * tmp - * lobatto_polynomials_grad[i].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points] - (1)) - * lobatto_polynomials[j + 2].value - (this->generalized_support_points[q_point + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < this->degree; ++j) + tmp -= local_dofs[i * this->degree + j] + * this->shape_value_component + (i * this->degree + j, + this->generalized_support_points[q_point + GeometryInfo::lines_per_cell - * n_edge_points] - (0)); - } + * n_edge_points], + 1); - system_matrix_inv.vmult (solution, system_rhs); + for (unsigned i = 0; i < this->degree; ++i) + for (unsigned int j = 0; j < this->degree-1; ++j) + system_rhs (i * (this->degree-1) + j) + += reference_quadrature.weight (q_point) * tmp + * lobatto_polynomials_grad[i].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points] + (1)) + * lobatto_polynomials[j + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points] + (0)); + } + + system_matrix_inv.vmult (solution, system_rhs); // Add the computed values // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - if (std::abs (solution (i * deg + j)) > 1e-14) - local_dofs[i + (j + GeometryInfo::lines_per_cell - + deg) * this->degree] - = solution (i * deg + j); - } + for (unsigned int i = 0; i < this->degree; ++i) + for (unsigned int j = 0; j < this->degree-1; ++j) + if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14) + local_dofs[i + (j + GeometryInfo::lines_per_cell + + this->degree-1) * this->degree] + = solution (i * (this->degree-1) + j); + } - break; - } + break; + } + + case 3: + { + const QGauss<1> + reference_edge_quadrature (this->degree); + const unsigned int& + n_edge_points = reference_edge_quadrature.size (); - case 3: - { // Let us begin with the - // edge degrees of freedom. - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) - { - const double weight = edge_quadrature.weight (q_point); - - for (unsigned int i = 0; i < this->degree; ++i) + // interpolation part. + for (unsigned int i = 0; i < 4; ++i) + { + for (unsigned int q_point = 0; q_point < n_edge_points; + ++q_point) + local_dofs[(i + 8) * this->degree] + += reference_edge_quadrature.weight (q_point) + * values[q_point + (i + 8) * n_edge_points] (2); + + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + if (std::abs (local_dofs[(i + 8) * this->degree]) < 1e-14) + local_dofs[(i + 8) * this->degree] = 0.0; + } + + if (offset + 1 < dim) + { + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < 2; ++j) { - const double L_i - = weight * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0)); - - for (unsigned int j = 0; j < 4; ++j) - local_dofs[i + (j + 8) * this->degree] - += L_i * values[q_point + (j + 8) * n_edge_points] (2); - - if (offset + 1 < dim) - { - for (unsigned int j = 0; j < 2; ++j) - for (unsigned int k = 0; k < 2; ++k) - local_dofs[i + (j + 4 * k) * this->degree] - += L_i * values[q_point + (j + 4 * k) * n_edge_points] (1); - - if (offset == 0) - for (unsigned int j = 0; j < 2; ++j) - for (unsigned int k = 0; k < 2; ++k) - local_dofs[i + (j + 4 * k + 2) * this->degree] - += L_i * values[q_point + (j + 4 * k + 2) * n_edge_points] (0); - } - } - } - + for (unsigned int q_point = 0; q_point < n_edge_points; + ++q_point) + local_dofs[(i + 4 * j) * this->degree] + += reference_edge_quadrature.weight (q_point) + * values[q_point + (i + 4 * j) * n_edge_points] + (1); + // Add the computed values // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i < GeometryInfo::lines_per_cell * this->degree; ++i) - if (std::abs (local_dofs[i]) < 1e-14) - local_dofs[i] = 0.0; + if (std::abs (local_dofs[(i + 4 * j) * this->degree]) + < 1e-14) + local_dofs[(i + 4 * j) * this->degree] = 0.0; + } - // If the degree is greater - // than 0, then we have still - // some interpolation to the - // face and interior shape - // functions left. - if (this->degree > 1) - { - const std::vector >& - lobatto_polynomials - = Polynomials::Lobatto::generate_complete_basis - (this->degree); - FullMatrix system_matrix (deg * this->degree, - deg * this->degree); - std::vector > - lobatto_polynomials_grad (this->degree); - - for (unsigned int i = 0; i < lobatto_polynomials_grad.size (); - ++i) - lobatto_polynomials_grad[i] - = lobatto_polynomials[i + 1].derivative (); - - // We set up the system - // matrix and use it for - // both, the horizontal - // and the vertical, shape - // functions. - const unsigned int - n_face_points = n_edge_points * n_edge_points; + if (offset == 0) + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < 2; ++j) + { + for (unsigned int q_point = 0; + q_point < n_edge_points; ++q_point) + local_dofs[(i + 4 * j + 2) * this->degree] + += reference_edge_quadrature.weight (q_point) + * values[q_point + (i + 4 * j + 2) + * n_edge_points] (0); - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - for (unsigned int k = 0; k < this->degree; ++k) - for (unsigned int l = 0; l < deg; ++l) - for (unsigned int q_point = 0; q_point < n_face_points; - ++q_point) - system_matrix (i * deg + j, k * deg + l) - += boundary_weights (q_point + n_edge_points, - 2 * (k * deg + l)) - * legendre_polynomials[i].value - (this->generalized_face_support_points[q_point - + 4 - * n_edge_points] - (0)) - * lobatto_polynomials[j + 2].value - (this->generalized_face_support_points[q_point - + 4 - * n_edge_points] - (1)); + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + if (std::abs (local_dofs[(i + 4 * j + 2) + * this->degree]) < 1e-14) + local_dofs[(i + 4 * j + 2) * this->degree] = 0.0; + } + } - FullMatrix system_matrix_inv (system_matrix.m (), - system_matrix.m ()); - - system_matrix_inv.invert (system_matrix); - - Vector solution (system_matrix.m ()); - Vector system_rhs (system_matrix.m ()); + // If the degree is greater + // than 0, then we have still + // some higher order shape + // functions to consider. + // Here the projection part + // starts. The dof values + // are obtained by solving + // a linear system of + // equations. + if (this->degree > 1) + { + // We start with projection + // on the higher order edge + // shape function. + const std::vector >& + lobatto_polynomials + = Polynomials::Lobatto::generate_complete_basis + (this->degree); + const unsigned int + line_coordinate[GeometryInfo<3>::lines_per_cell] + = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2}; + FullMatrix system_matrix (this->degree-1, this->degree-1); + FullMatrix system_matrix_inv (this->degree-1, this->degree-1); + std::vector > + lobatto_polynomials_grad (this->degree); + + for (unsigned int i = 0; i < lobatto_polynomials_grad.size (); + ++i) + lobatto_polynomials_grad[i] + = lobatto_polynomials[i + 1].derivative (); + + Vector system_rhs (system_matrix.m ()); + Vector solution (system_rhs.size ()); + + // Set up the system matrix. + // This can be used for all + // edges. + for (unsigned int i = 0; i < system_matrix.m (); ++i) + for (unsigned int j = 0; j < system_matrix.n (); ++j) + for (unsigned int q_point = 0; q_point < n_edge_points; + ++q_point) + system_matrix (i, j) + += boundary_weights (q_point, j) + * lobatto_polynomials_grad[i + 1].value + (this->generalized_face_support_points[q_point] + (1)); - for (unsigned int face = 0; - face < GeometryInfo::faces_per_cell; ++face) - { - switch (face) - { - case 0: - { - if (offset + 1 < dim) - { - // Set up the right hand side - // for the horizontal shape - // functions. - system_rhs = 0; + system_matrix_inv.invert (system_matrix); - for (unsigned int q_point = 0; - q_point < n_face_points; ++q_point) - { - double tmp - = values[q_point - + GeometryInfo::lines_per_cell - * n_edge_points] (1); + for (unsigned int line = 0; + line < GeometryInfo::lines_per_cell; ++line) + { + // Set up the right hand side. + system_rhs = 0; - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp - -= local_dofs[4 * i * this->degree - + j] - * this->shape_value_component - (4 * i * this->degree + j, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points], - 1); + if ((((line == 0) || (line == 1) || (line == 4) || + (line == 5)) && (offset + 1 < dim)) || + (((line == 2) || (line == 3) || (line == 6) || + (line == 7)) && (offset == 0)) || (line > 7)) + { + for (unsigned int q_point = 0; q_point < n_edge_points; + ++q_point) + { + double tmp + = values[line * n_edge_points + q_point] + (line_coordinate[line]) + - local_dofs[line * this->degree] + * this->shape_value_component + (line * this->degree, + this->generalized_support_points[line + * this->degree + + q_point], + line_coordinate[line]); + + for (unsigned int i = 0; i < system_rhs.size (); + ++i) + system_rhs (i) + += boundary_weights (q_point, i) * tmp; + } - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - system_rhs (i * deg + j) - += boundary_weights - (q_point + n_edge_points, - 2 * (i * deg + j)) * tmp; - } + system_matrix_inv.vmult (solution, system_rhs); - system_matrix_inv.vmult (solution, system_rhs); + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + local_dofs[line * this->degree + i + 1] + = solution (i); + } + } - // Add the computed values - // to the resulting vector - // only, if they are not - // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - if (std::abs (solution (i * deg + j)) - > 1e-14) - local_dofs[(i - + GeometryInfo::lines_per_cell) - * deg + j - + GeometryInfo::lines_per_cell] - = solution (i * deg + j); - } + // Then we go on to the + // face shape functions. + // Again we set up the + // system matrix and + // use it for both, the + // horizontal and the + // vertical, shape + // functions. + const std::vector >& + legendre_polynomials + = Polynomials::Legendre::generate_complete_basis (this->degree-1); + const unsigned int + n_face_points = n_edge_points * n_edge_points; + + system_matrix.reinit ((this->degree-1) * this->degree, + (this->degree-1) * this->degree); + system_matrix = 0; + + for (unsigned int i = 0; i < this->degree; ++i) + for (unsigned int j = 0; j < this->degree-1; ++j) + for (unsigned int k = 0; k < this->degree; ++k) + for (unsigned int l = 0; l < this->degree-1; ++l) + for (unsigned int q_point = 0; q_point < n_face_points; + ++q_point) + system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l) + += boundary_weights (q_point + n_edge_points, + 2 * (k * (this->degree-1) + l)) + * legendre_polynomials[i].value + (this->generalized_face_support_points[q_point + + 4 + * n_edge_points] + (0)) + * lobatto_polynomials[j + 2].value + (this->generalized_face_support_points[q_point + + 4 + * n_edge_points] + (1)); + + system_matrix_inv.reinit (system_matrix.m (), + system_matrix.n ()); + system_matrix_inv.invert (system_matrix); + solution.reinit (system_matrix.m ()); + system_rhs.reinit (system_matrix.m ()); + for (unsigned int face = 0; + face < GeometryInfo::faces_per_cell; ++face) + { + switch (face) + { + case 0: + { + if (offset + 1 < dim) + { // Set up the right hand side - // for the vertical shape + // for the horizontal shape // functions. - system_rhs = 0; + system_rhs = 0; - for (unsigned int q_point = 0; - q_point < n_face_points; ++q_point) - { - double tmp - = values[q_point - + GeometryInfo::lines_per_cell - * n_edge_points] (2); + for (unsigned int q_point = 0; + q_point < n_face_points; ++q_point) + { + double tmp + = values[q_point + + GeometryInfo::lines_per_cell + * n_edge_points] (1); - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp -= local_dofs[2 * (i + 4) - * this->degree + j] + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < this->degree; ++j) + tmp + -= local_dofs[4 * i * this->degree + + j] * this->shape_value_component - (2 * (i + 4) * this->degree + j, + (4 * i * this->degree + j, this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points], - 2); + 1); - for (unsigned i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - system_rhs (i * deg + j) - += boundary_weights - (q_point + n_edge_points, - 2 * (i * deg + j) + 1) - * tmp; - } + for (unsigned int i = 0; i < this->degree; ++i) + for (unsigned int j = 0; j < this->degree-1; ++j) + system_rhs (i * (this->degree-1) + j) + += boundary_weights + (q_point + n_edge_points, + 2 * (i * (this->degree-1) + j)) * tmp; + } - system_matrix_inv.vmult (solution, system_rhs); + system_matrix_inv.vmult (solution, system_rhs); // Add the computed values // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - if (std::abs (solution (i * deg + j)) > 1e-14) - local_dofs[i + (j + GeometryInfo::lines_per_cell - + deg) - * this->degree] - = solution (i * deg + j); - - break; - } - - case 1: - { - if (offset + 1 < dim) - { - // Set up the right hand side - // for the horizontal shape - // functions. - system_rhs = 0; - - for (unsigned int q_point = 0; - q_point < n_face_points; ++q_point) - { - double tmp - = values[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + n_face_points] (1); + for (unsigned int i = 0; i < this->degree; ++i) + for (unsigned int j = 0; j < this->degree-1; ++j) + if (std::abs (solution (i * (this->degree-1) + j)) + > 1e-14) + local_dofs[(i + + GeometryInfo::lines_per_cell) + * (this->degree-1) + j + + GeometryInfo::lines_per_cell] + = solution (i * (this->degree-1) + j); + } - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp -= local_dofs[(4 * i + 1) - * this->degree + j] - * this->shape_value_component - ((4 * i + 1) * this->degree - + j, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + n_face_points], - 1); + // Set up the right hand side + // for the vertical shape + // functions. + system_rhs = 0; - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - system_rhs (i * deg + j) - += boundary_weights - (q_point + n_edge_points, - 2 * (i * deg + j)) * tmp; - } + for (unsigned int q_point = 0; + q_point < n_face_points; ++q_point) + { + double tmp + = values[q_point + + GeometryInfo::lines_per_cell + * n_edge_points] (2); + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < this->degree; ++j) + tmp -= local_dofs[2 * (i + 4) + * this->degree + j] + * this->shape_value_component + (2 * (i + 4) * this->degree + j, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points], + 2); + + for (unsigned i = 0; i < this->degree; ++i) + for (unsigned int j = 0; j < this->degree-1; ++j) + system_rhs (i * (this->degree-1) + j) + += boundary_weights + (q_point + n_edge_points, + 2 * (i * (this->degree-1) + j) + 1) + * tmp; + } - system_matrix_inv.vmult (solution, system_rhs); + system_matrix_inv.vmult (solution, system_rhs); - // Add the computed values - // to the resulting vector - // only, if they are not - // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - if (std::abs (solution (i * deg + j)) - > 1e-14) - local_dofs[(i + GeometryInfo::lines_per_cell - + 2 * this->degree) * deg + j - + GeometryInfo::lines_per_cell] - = solution (i * deg + j); - } + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + for (unsigned int i = 0; i < this->degree; ++i) + for (unsigned int j = 0; j < this->degree-1; ++j) + if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14) + local_dofs[i + (j + GeometryInfo::lines_per_cell + + this->degree-1) + * this->degree] + = solution (i * (this->degree-1) + j); + + break; + } + case 1: + { + if (offset + 1 < dim) + { // Set up the right hand side - // for the vertical shape + // for the horizontal shape // functions. - system_rhs = 0; + system_rhs = 0; - for (unsigned int q_point = 0; - q_point < n_face_points; ++q_point) - { - double tmp - = values[q_point - + GeometryInfo::lines_per_cell - * n_edge_points + n_face_points] - (2); + for (unsigned int q_point = 0; + q_point < n_face_points; ++q_point) + { + double tmp + = values[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + n_face_points] (1); - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp -= local_dofs[(2 * (i + 4) + 1) - * this->degree + j] - * this->shape_value_component - ((2 * (i + 4) + 1) * this->degree - + j, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + n_face_points], - 2); + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= deg; ++j) + tmp -= local_dofs[(4 * i + 1) + * this->degree + j] + * this->shape_value_component + ((4 * i + 1) * this->degree + + j, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + n_face_points], + 1); - for (unsigned i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - system_rhs (i * deg + j) - += boundary_weights - (q_point + n_edge_points, - 2 * (i * deg + j) + 1) * tmp; - } + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + system_rhs (i * deg + j) + += boundary_weights + (q_point + n_edge_points, + 2 * (i * deg + j)) * tmp; + } - system_matrix_inv.vmult (solution, system_rhs); + system_matrix_inv.vmult (solution, system_rhs); // Add the computed values // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - if (std::abs (solution (i * deg + j)) > 1e-14) - local_dofs[i + (j + GeometryInfo::lines_per_cell - + 3 * deg) - * this->degree] - = solution (i * deg + j); - - break; - } + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + if (std::abs (solution (i * deg + j)) + > 1e-14) + local_dofs[(i + GeometryInfo::lines_per_cell + + 2 * this->degree) * deg + j + + GeometryInfo::lines_per_cell] + = solution (i * deg + j); + } - case 2: - { - if (offset == 0) - { - // Set up the right hand side - // for the horizontal shape - // functions. - system_rhs = 0; + // Set up the right hand side + // for the vertical shape + // functions. + system_rhs = 0; - for (unsigned int q_point = 0; - q_point < n_face_points; ++q_point) - { - double tmp - = values[q_point - + GeometryInfo::lines_per_cell - * n_edge_points + 2 * n_face_points] - (2); + for (unsigned int q_point = 0; + q_point < n_face_points; ++q_point) + { + double tmp + = values[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + n_face_points] + (2); + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= deg; ++j) + tmp -= local_dofs[(2 * (i + 4) + 1) + * this->degree + j] + * this->shape_value_component + ((2 * (i + 4) + 1) * this->degree + + j, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + n_face_points], + 2); - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp -= local_dofs[(i + 8) * this->degree - + j] - * this->shape_value_component - ((i + 8) * this->degree + j, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + 2 - * n_face_points], - 2); - - for (unsigned i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - system_rhs (i * deg + j) - += boundary_weights - (q_point + n_edge_points, - 2 * (i * deg + j)) * tmp; - } + for (unsigned i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + system_rhs (i * deg + j) + += boundary_weights + (q_point + n_edge_points, + 2 * (i * deg + j) + 1) * tmp; + } - system_matrix_inv.vmult (solution, system_rhs); + system_matrix_inv.vmult (solution, system_rhs); - // Add the computed values - // to the resulting vector - // only, if they are not - // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - if (std::abs (solution (i * deg + j)) - > 1e-14) - local_dofs[(i + GeometryInfo::lines_per_cell - + 4 * this->degree) * deg - + j - + GeometryInfo::lines_per_cell] - = solution (i * deg + j); - } + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + if (std::abs (solution (i * deg + j)) > 1e-14) + local_dofs[i + (j + GeometryInfo::lines_per_cell + + 3 * deg) + * this->degree] + = solution (i * deg + j); + + break; + } + case 2: + { + if (offset == 0) + { // Set up the right hand side - // for the vertical shape + // for the horizontal shape // functions. - system_rhs = 0; + system_rhs = 0; - for (unsigned int q_point = 0; - q_point < n_face_points; ++q_point) - { - double tmp - = values[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + 2 * n_face_points] (0); + for (unsigned int q_point = 0; + q_point < n_face_points; ++q_point) + { + double tmp + = values[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + 2 * n_face_points] + (2); - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp -= local_dofs[(4 * i + 2) - * this->degree + j] - * this->shape_value_component - ((4 * i + 2) * this->degree - + j, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + 2 - * n_face_points], - 0); - - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - system_rhs (i * deg + j) - += boundary_weights - (q_point + n_edge_points, - 2 * (i * deg + j) + 1) * tmp; - } + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= deg; ++j) + tmp -= local_dofs[(i + 8) * this->degree + + j] + * this->shape_value_component + ((i + 8) * this->degree + j, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + 2 + * n_face_points], + 2); + + for (unsigned i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + system_rhs (i * deg + j) + += boundary_weights + (q_point + n_edge_points, + 2 * (i * deg + j)) * tmp; + } - system_matrix_inv.vmult (solution, system_rhs); + system_matrix_inv.vmult (solution, system_rhs); // Add the computed values // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - if (std::abs (solution (i * deg + j)) > 1e-14) - local_dofs[i + (j + GeometryInfo::lines_per_cell - + 5 * deg) * this->degree] - = solution (i * deg + j); - - break; - } - - case 3: - { - if (offset == 0) - { - // Set up the right hand side - // for the horizontal shape - // functions. - system_rhs = 0; + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + if (std::abs (solution (i * deg + j)) + > 1e-14) + local_dofs[(i + GeometryInfo::lines_per_cell + + 4 * this->degree) * deg + + j + + GeometryInfo::lines_per_cell] + = solution (i * deg + j); + } - for (unsigned int q_point = 0; - q_point < n_face_points; ++q_point) - { - double tmp - = values[q_point - + GeometryInfo::lines_per_cell - * n_edge_points + 3 * n_face_points] - (2); + // Set up the right hand side + // for the vertical shape + // functions. + system_rhs = 0; - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp -= local_dofs[(i + 10) * this->degree - + j] - * this->shape_value_component - ((i + 10) * this->degree + j, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + 3 - * n_face_points], - 2); - - for (unsigned i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - system_rhs (i * deg + j) - += boundary_weights - (q_point + n_edge_points, - 2 * (i * deg + j)) * tmp; - } + for (unsigned int q_point = 0; + q_point < n_face_points; ++q_point) + { + double tmp + = values[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + 2 * n_face_points] (0); + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= deg; ++j) + tmp -= local_dofs[(4 * i + 2) + * this->degree + j] + * this->shape_value_component + ((4 * i + 2) * this->degree + + j, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + 2 + * n_face_points], + 0); + + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + system_rhs (i * deg + j) + += boundary_weights + (q_point + n_edge_points, + 2 * (i * deg + j) + 1) * tmp; + } - system_matrix_inv.vmult (solution, system_rhs); + system_matrix_inv.vmult (solution, system_rhs); - // Add the computed values - // to the resulting vector - // only, if they are not - // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - if (std::abs (solution (i * deg + j)) - > 1e-14) - local_dofs[(i + GeometryInfo::lines_per_cell - + 6 * this->degree) * deg + j - + GeometryInfo::lines_per_cell] - = solution (i * deg + j); - } + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + if (std::abs (solution (i * deg + j)) > 1e-14) + local_dofs[i + (j + GeometryInfo::lines_per_cell + + 5 * deg) * this->degree] + = solution (i * deg + j); + + break; + } + case 3: + { + if (offset == 0) + { // Set up the right hand side - // for the vertical shape + // for the horizontal shape // functions. - system_rhs = 0; + system_rhs = 0; - for (unsigned int q_point = 0; - q_point < n_face_points; ++q_point) - { - double tmp - = values[q_point - + GeometryInfo::lines_per_cell - * n_edge_points + 3 - * n_face_points] (0); + for (unsigned int q_point = 0; + q_point < n_face_points; ++q_point) + { + double tmp + = values[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + 3 * n_face_points] + (2); - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp -= local_dofs[(4 * i + 3) - * this->degree + j] - * this->shape_value_component - ((4 * i + 3) * this->degree - + j, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + 3 - * n_face_points], - 0); - - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - system_rhs (i * deg + j) - += boundary_weights - (q_point + n_edge_points, - 2 * (i * deg + j) + 1) * tmp; - } + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= deg; ++j) + tmp -= local_dofs[(i + 10) * this->degree + + j] + * this->shape_value_component + ((i + 10) * this->degree + j, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + 3 + * n_face_points], + 2); + + for (unsigned i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + system_rhs (i * deg + j) + += boundary_weights + (q_point + n_edge_points, + 2 * (i * deg + j)) * tmp; + } - system_matrix_inv.vmult (solution, system_rhs); + system_matrix_inv.vmult (solution, system_rhs); // Add the computed values // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - if (std::abs (solution (i * deg + j)) > 1e-14) - local_dofs[i + (j + GeometryInfo::lines_per_cell - + 7 * deg) * this->degree] - = solution (i * deg + j); - - break; - } - - case 4: - { - if (offset + 1 < dim) - { - // Set up the right hand side - // for the horizontal shape - // functions. - if (offset == 0) - { - system_rhs = 0; - - for (unsigned int q_point = 0; - q_point < n_face_points; ++q_point) - { - double tmp - = values[q_point - + GeometryInfo::lines_per_cell - * n_edge_points + 4 - * n_face_points] (0); - - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp -= local_dofs[(i + 2) - * this->degree - + j] - * this->shape_value_component - ((i + 2) * this->degree - + j, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + 4 - * n_face_points], - 0); - - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - system_rhs (i * deg + j) - += boundary_weights - (q_point + n_edge_points, - 2 * (i * deg + j)) * tmp; - } - - system_matrix_inv.vmult - (solution, system_rhs); - - // Add the computed values - // to the resulting vector - // only, if they are not - // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - if (std::abs (solution (i * deg + j)) - > 1e-14) - local_dofs[(i + GeometryInfo::lines_per_cell - + 8 * this->degree) * deg - + j - + GeometryInfo::lines_per_cell] - = solution (i * deg + j); - } - - // Set up the right hand side - // for the vertical shape - // functions. - system_rhs = 0; - - for (unsigned int q_point = 0; - q_point < n_face_points; ++q_point) - { - double tmp - = values[q_point - + GeometryInfo::lines_per_cell - * n_edge_points + 4 - * n_face_points] (1); - - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp -= local_dofs[i * this->degree + j] - * this->shape_value_component - (i * this->degree + j, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + 4 - * n_face_points], - 1); + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + if (std::abs (solution (i * deg + j)) + > 1e-14) + local_dofs[(i + GeometryInfo::lines_per_cell + + 6 * this->degree) * deg + j + + GeometryInfo::lines_per_cell] + = solution (i * deg + j); + } - for (unsigned i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - system_rhs (i * deg + j) - += boundary_weights - (q_point + n_edge_points, - 2 * (i * deg + j) + 1) * tmp; - } + // Set up the right hand side + // for the vertical shape + // functions. + system_rhs = 0; - system_matrix_inv.vmult (solution, system_rhs); + for (unsigned int q_point = 0; + q_point < n_face_points; ++q_point) + { + double tmp + = values[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + 3 + * n_face_points] (0); + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= deg; ++j) + tmp -= local_dofs[(4 * i + 3) + * this->degree + j] + * this->shape_value_component + ((4 * i + 3) * this->degree + + j, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + 3 + * n_face_points], + 0); + + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + system_rhs (i * deg + j) + += boundary_weights + (q_point + n_edge_points, + 2 * (i * deg + j) + 1) * tmp; + } - // Add the computed values - // to the resulting vector - // only, if they are not - // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - if (std::abs (solution (i * deg + j)) - > 1e-14) - local_dofs[i + (j + GeometryInfo::lines_per_cell - + 9 * deg) - * this->degree] - = solution (i * deg + j); - } + system_matrix_inv.vmult (solution, system_rhs); - break; - } + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + if (std::abs (solution (i * deg + j)) > 1e-14) + local_dofs[i + (j + GeometryInfo::lines_per_cell + + 7 * deg) * this->degree] + = solution (i * deg + j); + + break; + } - default: + case 4: + { if (offset + 1 < dim) { // Set up the right hand side @@ -2808,23 +3943,25 @@ FE_Nedelec::interpolate (std::vector& local_dofs, double tmp = values[q_point + GeometryInfo::lines_per_cell - * n_edge_points - + 5 * n_face_points] (0); + * n_edge_points + 4 + * n_face_points] (0); for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp -= local_dofs[(i + 6) - * this->degree + j] + for (unsigned int j = 0; j <= deg; ++j) + tmp -= local_dofs[(i + 2) + * this->degree + + j] * this->shape_value_component - ((i + 6) * this->degree + j, + ((i + 2) * this->degree + + j, this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points - + 5 + + 4 * n_face_points], 0); - for (unsigned int i = 0; i < this->degree; ++i) + for (unsigned int i = 0; i <= deg; ++i) for (unsigned int j = 0; j < deg; ++j) system_rhs (i * deg + j) += boundary_weights @@ -2835,18 +3972,18 @@ FE_Nedelec::interpolate (std::vector& local_dofs, system_matrix_inv.vmult (solution, system_rhs); - // Add the computed values - // to the resulting vector - // only, if they are not - // too small. - for (unsigned int i = 0; i < this->degree; ++i) + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + for (unsigned int i = 0; i <= deg; ++i) for (unsigned int j = 0; j < deg; ++j) if (std::abs (solution (i * deg + j)) > 1e-14) local_dofs[(i + GeometryInfo::lines_per_cell - + 10 * this->degree) - * deg + j - + GeometryInfo::lines_per_cell] + + 8 * this->degree) * deg + + j + + GeometryInfo::lines_per_cell] = solution (i * deg + j); } @@ -2861,23 +3998,22 @@ FE_Nedelec::interpolate (std::vector& local_dofs, double tmp = values[q_point + GeometryInfo::lines_per_cell - * n_edge_points + 5 + * n_edge_points + 4 * n_face_points] (1); for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < this->degree; ++j) - tmp -= local_dofs[(i + 4) - * this->degree + j] + for (unsigned int j = 0; j <= deg; ++j) + tmp -= local_dofs[i * this->degree + j] * this->shape_value_component - ((i + 4) * this->degree + j, + (i * this->degree + j, this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points - + 5 + + 4 * n_face_points], 1); - for (unsigned i = 0; i < this->degree; ++i) + for (unsigned i = 0; i <= deg; ++i) for (unsigned int j = 0; j < deg; ++j) system_rhs (i * deg + j) += boundary_weights @@ -2891,397 +4027,508 @@ FE_Nedelec::interpolate (std::vector& local_dofs, // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i < this->degree; ++i) + for (unsigned int i = 0; i <= deg; ++i) for (unsigned int j = 0; j < deg; ++j) if (std::abs (solution (i * deg + j)) > 1e-14) local_dofs[i + (j + GeometryInfo::lines_per_cell - + 11 * deg) * this->degree] + + 9 * deg) + * this->degree] = solution (i * deg + j); } - } - } - // Finally we project - // the remaining parts - // of the function on - // the interior shape + break; + } + + default: + if (offset + 1 < dim) + { + // Set up the right hand side + // for the horizontal shape // functions. - const QGauss reference_quadrature (this->degree); - const unsigned int& - n_interior_points = reference_quadrature.size (); + if (offset == 0) + { + system_rhs = 0; - // We create the - // system matrix. - system_matrix.reinit (this->degree * deg * deg, - this->degree * deg * deg); - system_matrix = 0; + for (unsigned int q_point = 0; + q_point < n_face_points; ++q_point) + { + double tmp + = values[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + 5 * n_face_points] (0); - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - for (unsigned int k = 0; k < deg; ++k) - for (unsigned int l = 0; l < this->degree; ++l) - for (unsigned int m = 0; m < deg; ++m) - for (unsigned int n = 0; n < deg; ++n) - for (unsigned int q_point = 0; - q_point < n_interior_points; ++q_point) - system_matrix ((i * deg + j) * deg + k, - (l * deg + m) * deg + n) - += reference_quadrature.weight (q_point) - * legendre_polynomials[i].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] - (0)) * lobatto_polynomials[j + 2].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] - (1)) - * lobatto_polynomials[k + 2].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] - (2)) - * lobatto_polynomials_grad[l].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] - (0)) - * lobatto_polynomials[m + 2].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] - (1)) - * lobatto_polynomials[n + 2].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] - (2)); + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= deg; ++j) + tmp -= local_dofs[(i + 6) + * this->degree + j] + * this->shape_value_component + ((i + 6) * this->degree + j, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + 5 + * n_face_points], + 0); + + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + system_rhs (i * deg + j) + += boundary_weights + (q_point + n_edge_points, + 2 * (i * deg + j)) * tmp; + } - system_matrix_inv.reinit (system_matrix.m (), - system_matrix.m ()); - system_matrix_inv.invert (system_matrix); - system_rhs.reinit (system_matrix_inv.m ()); - solution.reinit (system_matrix.m ()); + system_matrix_inv.vmult + (solution, system_rhs); - if (offset + 1 < dim) - { - if (offset == 0) - { - // Set up the right hand side. - system_rhs = 0; + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + if (std::abs (solution (i * deg + j)) + > 1e-14) + local_dofs[(i + GeometryInfo::lines_per_cell + + 10 * this->degree) + * deg + j + + GeometryInfo::lines_per_cell] + = solution (i * deg + j); + } - for (unsigned int q_point = 0; - q_point < n_interior_points; ++q_point) - { - double tmp - = values[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] (0); - - for (unsigned int i = 0; i < this->degree; ++i) + // Set up the right hand side + // for the vertical shape + // functions. + system_rhs = 0; + + for (unsigned int q_point = 0; + q_point < n_face_points; ++q_point) { - for (unsigned int j = 0; j < 2; ++j) - for (unsigned int k = 0; k < 2; ++k) - tmp -= local_dofs[i + (j + 4 * k + 2) - * this->degree] - * this->shape_value_component - (i + (j + 4 * k + 2) - * this->degree, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points], - 0); + double tmp + = values[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + 5 + * n_face_points] (1); - for (unsigned int j = 0; j < deg; ++j) - for (unsigned int k = 0; k < 4; ++k) - tmp -= local_dofs[(i + 2 * (k + 2) - * this->degree - + GeometryInfo::lines_per_cell) - * deg + j - + GeometryInfo::lines_per_cell] + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= deg; ++j) + tmp -= local_dofs[(i + 4) + * this->degree + j] * this->shape_value_component - ((i + 2 * (k + 2) * this->degree - + GeometryInfo::lines_per_cell) - * deg + j - + GeometryInfo::lines_per_cell, + ((i + 4) * this->degree + j, this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points - + GeometryInfo::faces_per_cell + + 5 * n_face_points], - 0); + 1); + + for (unsigned i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + system_rhs (i * deg + j) + += boundary_weights + (q_point + n_edge_points, + 2 * (i * deg + j) + 1) * tmp; } - for (unsigned int i = 0; i < this->degree; ++i) + system_matrix_inv.vmult (solution, system_rhs); + + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + for (unsigned int i = 0; i <= deg; ++i) for (unsigned int j = 0; j < deg; ++j) - for (unsigned int k = 0; k < deg; ++k) - system_rhs ((i * deg + j) * deg + k) - += reference_quadrature.weight (q_point) - * tmp - * lobatto_polynomials_grad[i].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] - (0)) - * lobatto_polynomials[j + 2].value + if (std::abs (solution (i * deg + j)) + > 1e-14) + local_dofs[i + (j + GeometryInfo::lines_per_cell + + 11 * deg) * this->degree] + = solution (i * deg + j); + } + } + } + + // Finally we project + // the remaining parts + // of the function on + // the interior shape + // functions. + const QGauss reference_quadrature (this->degree); + const unsigned int& + n_interior_points = reference_quadrature.size (); + + // We create the + // system matrix. + system_matrix.reinit (this->degree * deg * deg, + this->degree * deg * deg); + system_matrix = 0; + + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + for (unsigned int k = 0; k < deg; ++k) + for (unsigned int l = 0; l <= deg; ++l) + for (unsigned int m = 0; m < deg; ++m) + for (unsigned int n = 0; n < deg; ++n) + for (unsigned int q_point = 0; + q_point < n_interior_points; ++q_point) + system_matrix ((i * deg + j) * deg + k, + (l * deg + m) * deg + n) + += reference_quadrature.weight (q_point) + * legendre_polynomials[i].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (0)) * lobatto_polynomials[j + 2].value (this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * n_face_points] - (1)) + (1)) * lobatto_polynomials[k + 2].value (this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * n_face_points] - (2)); - } - - system_matrix_inv.vmult (solution, system_rhs); + (2)) + * lobatto_polynomials_grad[l].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (0)) + * lobatto_polynomials[m + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (1)) + * lobatto_polynomials[n + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (2)); - // Add the computed values - // to the resulting vector - // only, if they are not - // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - for (unsigned int k = 0; k < deg; ++k) - if (std::abs (solution ((i * deg + j) * deg + k)) - > 1e-14) - local_dofs[((i + 2 - * GeometryInfo::faces_per_cell) - * deg + j - + GeometryInfo::lines_per_cell - + 2 - * GeometryInfo::faces_per_cell) - * deg + k - + GeometryInfo::lines_per_cell] - = solution ((i * deg + j) * deg + k); - } + system_matrix_inv.reinit (system_matrix.m (), + system_matrix.m ()); + system_matrix_inv.invert (system_matrix); + system_rhs.reinit (system_matrix_inv.m ()); + solution.reinit (system_matrix.m ()); + if (offset + 1 < dim) + { + if (offset == 0) + { // Set up the right hand side. - system_rhs = 0; - - for (unsigned int q_point = 0; q_point < n_interior_points; - ++q_point) - { - double tmp - = values[q_point + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] (1); + system_rhs = 0; - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < 2; ++j) + for (unsigned int q_point = 0; + q_point < n_interior_points; ++q_point) + { + double tmp + = values[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] (0); + + for (unsigned int i = 0; i <= deg; ++i) { - for (unsigned int k = 0; k < 2; ++k) - tmp -= local_dofs[i + (4 * j + k) + for (unsigned int j = 0; j < 2; ++j) + for (unsigned int k = 0; k < 2; ++k) + tmp -= local_dofs[i + (j + 4 * k + 2) * this->degree] - * this->shape_value_component - (i + (4 * j + k) * this->degree, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points], - 1); + * this->shape_value_component + (i + (j + 4 * k + 2) + * this->degree, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points], + 0); - for (unsigned int k = 0; k < deg; ++k) - tmp -= local_dofs[(i + 2 * j * this->degree + for (unsigned int j = 0; j < deg; ++j) + for (unsigned int k = 0; k < 4; ++k) + tmp -= local_dofs[(i + 2 * (k + 2) + * this->degree + GeometryInfo::lines_per_cell) - * deg + k - + GeometryInfo::lines_per_cell] - * this->shape_value_component - ((i + 2 * j * this->degree + * deg + j + + GeometryInfo::lines_per_cell] + * this->shape_value_component + ((i + 2 * (k + 2) * this->degree + GeometryInfo::lines_per_cell) - * deg + k - + GeometryInfo::lines_per_cell, - this->generalized_support_points[q_point + * deg + j + + GeometryInfo::lines_per_cell, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points], + 0); + } + + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + for (unsigned int k = 0; k < deg; ++k) + system_rhs ((i * deg + j) * deg + k) + += reference_quadrature.weight (q_point) + * tmp + * lobatto_polynomials_grad[i].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (0)) + * lobatto_polynomials[j + 2].value + (this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell - * n_face_points], - 1) - + local_dofs[i + ((2 * j + 9) * deg + k - + GeometryInfo::lines_per_cell) - * this->degree] - * this->shape_value_component - (i + ((2 * j + 9) * deg + k - + GeometryInfo::lines_per_cell) - * this->degree, - this->generalized_support_points[q_point + * n_face_points] + (1)) + * lobatto_polynomials[k + 2].value + (this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell - * n_face_points], - 1); - } - - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - for (unsigned int k = 0; k < deg; ++k) - system_rhs ((i * deg + j) * deg + k) - += reference_quadrature.weight (q_point) * tmp - * lobatto_polynomials_grad[i].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] - (1)) - * lobatto_polynomials[j + 2].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] - (0)) - * lobatto_polynomials[k + 2].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] - (2)); - } + * n_face_points] + (2)); + } - system_matrix_inv.vmult (solution, system_rhs); + system_matrix_inv.vmult (solution, system_rhs); // Add the computed values // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - for (unsigned int k = 0; k < deg; ++k) - if (std::abs (solution ((i * deg + j) * deg + k)) - > 1e-14) - local_dofs[((i + this->degree + 2 + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + for (unsigned int k = 0; k < deg; ++k) + if (std::abs (solution ((i * deg + j) * deg + k)) + > 1e-14) + local_dofs[((i + 2 * GeometryInfo::faces_per_cell) - * deg + j - + GeometryInfo::lines_per_cell + 2 - * GeometryInfo::faces_per_cell) - * deg + k - + GeometryInfo::lines_per_cell] + * deg + j + + GeometryInfo::lines_per_cell + + 2 + * GeometryInfo::faces_per_cell) + * deg + k + + GeometryInfo::lines_per_cell] = solution ((i * deg + j) * deg + k); - } + } // Set up the right hand side. - system_rhs = 0; + system_rhs = 0; - for (unsigned int q_point = 0; q_point < n_interior_points; - ++q_point) - { - double tmp - = values[q_point + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] (2); + for (unsigned int q_point = 0; q_point < n_interior_points; + ++q_point) + { + double tmp + = values[q_point + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] (1); - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < 4; ++j) - { - tmp -= local_dofs[i + (j + 8) * this->degree] - * this->shape_value_component - (i + (j + 8) * this->degree, - this->generalized_support_points[q_point + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < 2; ++j) + { + for (unsigned int k = 0; k < 2; ++k) + tmp -= local_dofs[i + (4 * j + k) + * this->degree] + * this->shape_value_component + (i + (4 * j + k) * this->degree, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points], + 1); + + for (unsigned int k = 0; k < deg; ++k) + tmp -= local_dofs[(i + 2 * j * this->degree + + GeometryInfo::lines_per_cell) + * deg + k + + GeometryInfo::lines_per_cell] + * this->shape_value_component + ((i + 2 * j * this->degree + + GeometryInfo::lines_per_cell) + * deg + k + + GeometryInfo::lines_per_cell, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points], + 1) + + local_dofs[i + ((2 * j + 9) * deg + k + + GeometryInfo::lines_per_cell) + * this->degree] + * this->shape_value_component + (i + ((2 * j + 9) * deg + k + + GeometryInfo::lines_per_cell) + * this->degree, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points], + 1); + } + + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + for (unsigned int k = 0; k < deg; ++k) + system_rhs ((i * deg + j) * deg + k) + += reference_quadrature.weight (q_point) * tmp + * lobatto_polynomials_grad[i].value + (this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell - * n_face_points], - 2); + * n_face_points] + (1)) + * lobatto_polynomials[j + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (0)) + * lobatto_polynomials[k + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (2)); + } - for (unsigned int k = 0; k < deg; ++k) - tmp -= local_dofs[i + ((2 * j + 1) * deg + k - + GeometryInfo::lines_per_cell) - * this->degree] - * this->shape_value_component - (i + ((2 * j + 1) * deg + k - + GeometryInfo::lines_per_cell) - * this->degree, - this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points], - 2); - } + system_matrix_inv.vmult (solution, system_rhs); - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - for (unsigned int k = 0; k < deg; ++k) - system_rhs ((i * deg + j) * deg + k) - += reference_quadrature.weight (q_point) * tmp - * lobatto_polynomials_grad[i].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] - (2)) - * lobatto_polynomials[j + 2].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points - + GeometryInfo::faces_per_cell - * n_face_points] - (0)) - * lobatto_polynomials[k + 2].value - (this->generalized_support_points[q_point + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + for (unsigned int k = 0; k < deg; ++k) + if (std::abs (solution ((i * deg + j) * deg + k)) + > 1e-14) + local_dofs[((i + this->degree + 2 + * GeometryInfo::faces_per_cell) + * deg + j + + GeometryInfo::lines_per_cell + 2 + * GeometryInfo::faces_per_cell) + * deg + k + + GeometryInfo::lines_per_cell] + = solution ((i * deg + j) * deg + k); + } + + // Set up the right hand side. + system_rhs = 0; + + for (unsigned int q_point = 0; q_point < n_interior_points; + ++q_point) + { + double tmp + = values[q_point + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] (2); + + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < 4; ++j) + { + tmp -= local_dofs[i + (j + 8) * this->degree] + * this->shape_value_component + (i + (j + 8) * this->degree, + this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell - * n_face_points] - (1)); - } + * n_face_points], + 2); - system_matrix_inv.vmult (solution, system_rhs); + for (unsigned int k = 0; k < deg; ++k) + tmp -= local_dofs[i + ((2 * j + 1) * deg + k + + GeometryInfo::lines_per_cell) + * this->degree] + * this->shape_value_component + (i + ((2 * j + 1) * deg + k + + GeometryInfo::lines_per_cell) + * this->degree, + this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points], + 2); + } + + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + for (unsigned int k = 0; k < deg; ++k) + system_rhs ((i * deg + j) * deg + k) + += reference_quadrature.weight (q_point) * tmp + * lobatto_polynomials_grad[i].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (2)) + * lobatto_polynomials[j + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (0)) + * lobatto_polynomials[k + 2].value + (this->generalized_support_points[q_point + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_points] + (1)); + } + + system_matrix_inv.vmult (solution, system_rhs); // Add the computed values // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) - for (unsigned int k = 0; k < deg; ++k) - if (std::abs (solution ((i * deg + j) * deg + k)) - > 1e-14) - local_dofs[i + ((j + 2 - * (deg + GeometryInfo::faces_per_cell)) - * deg + k - + GeometryInfo::lines_per_cell) - * this->degree] - = solution ((i * deg + j) * deg + k); - } + for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int j = 0; j < deg; ++j) + for (unsigned int k = 0; k < deg; ++k) + if (std::abs (solution ((i * deg + j) * deg + k)) + > 1e-14) + local_dofs[i + ((j + 2 + * (deg + GeometryInfo::faces_per_cell)) + * deg + k + + GeometryInfo::lines_per_cell) + * this->degree] + = solution ((i * deg + j) * deg + k); + } - break; - } + break; + } - default: - Assert (false, ExcNotImplemented ()); - } - } + default: + Assert (false, ExcNotImplemented ()); + } } @@ -3297,6 +4544,7 @@ FE_Nedelec::interpolate (std::vector& local_dofs, const VectorSlice > >& values) const { + const unsigned int deg = this->degree-1; Assert (values.size () == this->n_components (), ExcDimensionMismatch (values.size (), this->n_components ())); Assert (values[0].size () == this->generalized_support_points.size (), @@ -3306,56 +4554,53 @@ const ExcDimensionMismatch (local_dofs.size (), this->dofs_per_cell)); std::fill (local_dofs.begin (), local_dofs.end (), 0.0); - const unsigned int deg = this->degree-1; - const std::vector >& - legendre_polynomials - = Polynomials::Legendre::generate_complete_basis (deg); - const QGauss<1> edge_quadrature (this->degree); - const std::vector >& - edge_quadrature_points = edge_quadrature.get_points (); - const unsigned int& n_edge_points = edge_quadrature.size (); - switch (dim) { case 2: { // Let us begin with the - // edge degrees of freedom. - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) - { - const double weight = edge_quadrature.weight (q_point); - - for (unsigned int i = 0; i < this->degree; ++i) - { - const double L_i - = weight * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0)); - - for (unsigned int j = 0; j < 2; ++j) - for (unsigned int k = 0; k < 2; ++k) - local_dofs[i + (j + 2 * k) * this->degree] - += L_i * values[1 - k][q_point + (j + 2 * k) * n_edge_points]; - } - } - + // interpolation part. + const QGauss reference_edge_quadrature (this->degree); + const unsigned int& + n_edge_points = reference_edge_quadrature.size (); + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < 2; ++j) + { + for (unsigned int q_point = 0; q_point < n_edge_points; + ++q_point) + local_dofs[(i + 2 * j) * this->degree] + += reference_edge_quadrature.weight (q_point) + * values[1 - j][q_point + (i + 2 * j) * n_edge_points]; + // Add the computed values // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i < GeometryInfo::lines_per_cell * this->degree; ++i) - if (std::abs (local_dofs[i]) < 1e-14) - local_dofs[i] = 0.0; + if (std::abs (local_dofs[(i + 2 * j) * this->degree]) < 1e-14) + local_dofs[(i + 2 * j) * this->degree] = 0.0; + } // If the degree is greater // than 0, then we have still - // some interpolations onto - // the interior shape - // functions left. - if (this->degree > 1) + // some higher order edge + // shape functions to + // consider. + // Here the projection part + // starts. The dof values + // are obtained by solving + // a linear system of + // equations. + if (this->degree-1 > 1) { + // We start with projection + // on the higher order edge + // shape function. const std::vector >& lobatto_polynomials = Polynomials::Lobatto::generate_complete_basis (this->degree); + FullMatrix system_matrix (this->degree-1, this->degree-1); std::vector > lobatto_polynomials_grad (this->degree); @@ -3364,26 +4609,91 @@ const lobatto_polynomials_grad[i] = lobatto_polynomials[i + 1].derivative (); - // We set up the system - // matrix and use it for - // both, the horizontal - // and the vertical + // Set up the system matrix. + // This can be used for all + // edges. + for (unsigned int i = 0; i < system_matrix.m (); ++i) + for (unsigned int j = 0; j < system_matrix.n (); ++j) + for (unsigned int q_point = 0; q_point < n_edge_points; + ++q_point) + system_matrix (i, j) + += boundary_weights (q_point, j) + * lobatto_polynomials_grad[i + 1].value + (this->generalized_face_support_points[q_point] + (1)); + + FullMatrix system_matrix_inv (this->degree-1, this->degree-1); + + system_matrix_inv.invert (system_matrix); + + const unsigned int + line_coordinate[GeometryInfo<2>::lines_per_cell] + = {1, 1, 0, 0}; + Vector system_rhs (system_matrix.m ()); + Vector solution (system_rhs.size ()); + + for (unsigned int line = 0; + line < GeometryInfo::lines_per_cell; ++line) + { + // Set up the right hand side. + system_rhs = 0; + + for (unsigned int q_point = 0; q_point < n_edge_points; + ++q_point) + { + const double tmp + = values[line_coordinate[line]][line * n_edge_points + + q_point] + - local_dofs[line * this->degree] + * this->shape_value_component + (line * this->degree, + this->generalized_support_points[line + * n_edge_points + + q_point], + line_coordinate[line]); + + for (unsigned int i = 0; i < system_rhs.size (); ++i) + system_rhs (i) += boundary_weights (q_point, i) * tmp; + } + + system_matrix_inv.vmult (solution, system_rhs); + + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + local_dofs[line * this->degree + i + 1] = solution (i); + } + + // Then we go on to the // interior shape - // functions. + // functions. Again we + // set up the system + // matrix and use it + // for both, the + // horizontal and the + // vertical, interior + // shape functions. const QGauss reference_quadrature (this->degree); - const unsigned int& n_interior_points - = reference_quadrature.size (); + const unsigned int& + n_interior_points = reference_quadrature.size (); + const std::vector >& + legendre_polynomials + = Polynomials::Legendre::generate_complete_basis (this->degree-1); - FullMatrix system_matrix (deg * this->degree, - deg * this->degree); + system_matrix.reinit ((this->degree-1) * this->degree, + (this->degree-1) * this->degree); + system_matrix = 0; for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) + for (unsigned int j = 0; j < this->degree-1; ++j) for (unsigned int k = 0; k < this->degree; ++k) - for (unsigned int l = 0; l < deg; ++l) + for (unsigned int l = 0; l < this->degree-1; ++l) for (unsigned int q_point = 0; q_point < n_interior_points; ++q_point) - system_matrix (i * deg + j, k * deg + l) + system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l) += reference_quadrature.weight (q_point) * legendre_polynomials[i].value (this->generalized_support_points[q_point @@ -3406,13 +4716,14 @@ const * n_edge_points] (1)); - FullMatrix system_matrix_inv (system_matrix.m (), - system_matrix.m ()); - + system_matrix_inv.reinit (system_matrix.m (), + system_matrix.m ()); system_matrix_inv.invert (system_matrix); - - Vector solution (system_matrix_inv.m ()); - Vector system_rhs (system_matrix.m ()); + // Set up the right hand side + // for the horizontal shape + // functions. + system_rhs.reinit (system_matrix_inv.m ()); + system_rhs = 0; for (unsigned int q_point = 0; q_point < n_interior_points; ++q_point) @@ -3519,49 +4830,61 @@ const case 3: { // Let us begin with the - // edge degrees of freedom. + // interpolation part. + const QGauss<1> reference_edge_quadrature (this->degree); + const unsigned int& + n_edge_points = reference_edge_quadrature.size (); + for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - const double weight = edge_quadrature.weight (q_point); - - for (unsigned int i = 0; i < this->degree; ++i) - { - const double L_i - = weight * legendre_polynomials[i].value (edge_quadrature_points[q_point] (0)); - - for (unsigned int j = 0; j < 4; ++j) - local_dofs[i + (j + 8) * this->degree] - += L_i * values[2][q_point + (j + 8) * n_edge_points]; - - for (unsigned int j = 0; j < 2; ++j) - for (unsigned int k = 0; k < 2; ++k) - for (unsigned int l = 0; l < 2; ++l) - local_dofs[i + (j + 4 * k + 2 * l) * this->degree] - += L_i * values[1 - l][q_point + (j + 4 * k + 2 * l) * n_edge_points]; - } + for (unsigned int i = 0; i < 4; ++i) + local_dofs[(i + 8) * this->degree] + += reference_edge_quadrature.weight (q_point) + * values[2][q_point + (i + 8) * n_edge_points]; + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < 2; ++j) + for (unsigned int k = 0; k < 2; ++k) + local_dofs[(i + 2 * (2 * j + k)) * this->degree] + += reference_edge_quadrature.weight (q_point) + * values[1 - k][q_point + (i + 2 * (2 * j + k)) + * n_edge_points]; } - + // Add the computed values // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i < GeometryInfo::lines_per_cell * this->degree; ++i) - if (std::abs (local_dofs[i]) < 1e-14) - local_dofs[i] = 0.0; + for (unsigned int i = 0; i < 4; ++i) + if (std::abs (local_dofs[(i + 8) * this->degree]) < 1e-14) + local_dofs[(i + 8) * this->degree] = 0.0; + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < 2; ++j) + for (unsigned int k = 0; k < 2; ++k) + if (std::abs (local_dofs[(i + 2 * (2 * j + k)) * this->degree]) + < 1e-14) + local_dofs[(i + 2 * (2 * j + k)) * this->degree] = 0.0; // If the degree is greater // than 0, then we have still - // some interpolation to the - // face and interior shape - // functions left. + // some higher order shape + // functions to consider. + // Here the projection part + // starts. The dof values + // are obtained by solving + // a linear system of + // equations. if (this->degree > 1) { + // We start with projection + // on the higher order edge + // shape function. const std::vector >& lobatto_polynomials = Polynomials::Lobatto::generate_complete_basis (this->degree); - FullMatrix system_matrix (deg * this->degree, - deg * this->degree); + FullMatrix system_matrix (this->degree-1, this->degree-1); std::vector > lobatto_polynomials_grad (this->degree); @@ -3570,23 +4893,90 @@ const lobatto_polynomials_grad[i] = lobatto_polynomials[i + 1].derivative (); - // We set up the system - // matrix and use it for - // both, the horizontal - // and the vertical, shape - // functions. + // Set up the system matrix. + // This can be used for all + // edges. + for (unsigned int i = 0; i < system_matrix.m (); ++i) + for (unsigned int j = 0; j < system_matrix.n (); ++j) + for (unsigned int q_point = 0; q_point < n_edge_points; + ++q_point) + system_matrix (i, j) + += boundary_weights (q_point, j) + * lobatto_polynomials_grad[i + 1].value + (this->generalized_face_support_points[q_point] + (1)); + + FullMatrix system_matrix_inv (this->degree-1, this->degree-1); + + system_matrix_inv.invert (system_matrix); + const unsigned int - n_face_points = n_edge_points * n_edge_points; + line_coordinate[GeometryInfo<3>::lines_per_cell] + = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2}; + Vector system_rhs (system_matrix.m ()); + Vector solution (system_rhs.size ()); + + for (unsigned int line = 0; + line < GeometryInfo::lines_per_cell; ++line) + { + // Set up the right hand side. + system_rhs = 0; + + for (unsigned int q_point = 0; q_point < this->degree; ++q_point) + { + const double tmp + = values[line_coordinate[line]][line * this->degree + + q_point] + - local_dofs[line * this->degree] + * this->shape_value_component + (line * this->degree, + this->generalized_support_points[line + * this->degree + + q_point], + line_coordinate[line]); + + for (unsigned int i = 0; i < system_rhs.size (); ++i) + system_rhs (i) += boundary_weights (q_point, i) + * tmp; + } + + system_matrix_inv.vmult (solution, system_rhs); + + // Add the computed values + // to the resulting vector + // only, if they are not + // too small. + for (unsigned int i = 0; i < solution.size (); ++i) + if (std::abs (solution (i)) > 1e-14) + local_dofs[line * this->degree + i + 1] = solution (i); + } + + // Then we go on to the + // face shape functions. + // Again we set up the + // system matrix and + // use it for both, the + // horizontal and the + // vertical, shape + // functions. + const std::vector >& + legendre_polynomials + = Polynomials::Legendre::generate_complete_basis (this->degree-1); + const unsigned int n_face_points = n_edge_points * n_edge_points; + + system_matrix.reinit ((this->degree-1) * this->degree, + (this->degree-1) * this->degree); + system_matrix = 0; for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < deg; ++j) + for (unsigned int j = 0; j < this->degree-1; ++j) for (unsigned int k = 0; k < this->degree; ++k) - for (unsigned int l = 0; l < deg; ++l) + for (unsigned int l = 0; l < this->degree-1; ++l) for (unsigned int q_point = 0; q_point < n_face_points; ++q_point) - system_matrix (i * deg + j, k * deg + l) + system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l) += boundary_weights (q_point + n_edge_points, - 2 * (k * deg + l)) + 2 * (k * (this->degree-1) + l)) * legendre_polynomials[i].value (this->generalized_face_support_points[q_point + 4 @@ -3598,11 +4988,12 @@ const * n_edge_points] (1)); - FullMatrix system_matrix_inv (system_matrix.m (), - system_matrix.m ()); - + system_matrix_inv.reinit (system_matrix.m (), + system_matrix.m ()); system_matrix_inv.invert (system_matrix); - + solution.reinit (system_matrix.m ()); + system_rhs.reinit (system_matrix.m ()); + const unsigned int face_coordinates[GeometryInfo<3>::faces_per_cell][2] = {{1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}}; @@ -3610,8 +5001,6 @@ const edge_indices[GeometryInfo<3>::faces_per_cell][GeometryInfo<3>::lines_per_face] = {{0, 4, 8, 10}, {1, 5, 9, 11}, {8, 9, 2, 6}, {10, 11, 3, 7}, {2, 3, 0, 1}, {6, 7, 4, 5}}; - Vector solution (system_matrix.m ()); - Vector system_rhs (system_matrix.m ()); for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) -- 2.39.5