From 57479aea58f4c636e8a134fe1192388bc61d8feb Mon Sep 17 00:00:00 2001 From: buerg Date: Sun, 22 Jul 2012 15:21:35 +0000 Subject: [PATCH] Replaced projection-based interpolation by explicit interpolation formulas wherever possible. git-svn-id: https://svn.dealii.org/trunk@25717 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/source/fe/fe_nedelec.cc | 4892 +++++++++++-------------------- 1 file changed, 1718 insertions(+), 3174 deletions(-) diff --git a/deal.II/source/fe/fe_nedelec.cc b/deal.II/source/fe/fe_nedelec.cc index e1c05bb37b..b1da9ccd6e 100644 --- a/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/source/fe/fe_nedelec.cc @@ -16,10 +16,6 @@ #include #include -//TODO: implement the adjust_quad_dof_index_for_face_orientation_table and -//adjust_line_dof_index_for_line_orientation_table fields, and write tests -//similar to bits/face_orientation_and_fe_q_* - DEAL_II_NAMESPACE_OPEN @@ -161,7 +157,6 @@ FE_PolyTensor, dim> default: Assert (false, ExcNotImplemented ()); } - } @@ -497,6 +492,8 @@ FE_Nedelec<1>::initialize_restriction () + + // Restriction operator template void @@ -515,6 +512,9 @@ 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) { @@ -522,232 +522,83 @@ FE_Nedelec::initialize_restriction () { // First interpolate the shape // functions of the child cells - // to the lowest order shape - // functions of the parent cell. + // to the face 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) - { - 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); - } + 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); + } else - { - 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); - } + 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); + } } - // Then project the shape functions - // of the child cells to the higher - // order shape functions of the - // parent cell. + // Then interpolate the shape functions + // of the child cells to the interior + // 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 (); @@ -757,185 +608,71 @@ 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 (); - { - FullMatrix assembling_matrix ((this->degree-1) * this->degree, - n_quadrature_points); - - for (unsigned int q_point = 0; q_point < n_quadrature_points; + for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) + for (unsigned int q_point = 0; q_point < quadrature.size (); ++q_point) { - const double weight - = std::sqrt (quadrature.weight (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; - for (unsigned int i = 0; i < this->degree; ++i) + if (quadrature_points[q_point] (0) < 0.5) { - 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) - { - 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); - } - + k = 0; + else { - 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); - } - - 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); + quadrature_point (1) -= 1.0; + k = 1; } - - tmp *= quadrature.weight (q_point); - - for (unsigned int i = 0; i < this->degree; ++i) + } + + else + if (quadrature_points[q_point] (1) < 0.5) { - 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; - } - } + quadrature_point (0) -= 1.0; + k = 2; } - } - - 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) + + else { - 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); + quadrature_point (0) -= 1.0; + quadrature_point (1) -= 1.0; + k = 3; } - } + + 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; @@ -945,292 +682,96 @@ FE_Nedelec::initialize_restriction () { // First interpolate the shape // functions of the child cells - // to the lowest order shape - // functions of the parent cell. + // to the edge 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 < 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 i = 0; i < this->degree; ++i) + { + const double L_i + = weight * legendre_polynomials[i].value + (edge_quadrature_points[q_point] (0)); - else - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < 2; ++j) + 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)); + + for (unsigned int j = 0; j < 2; ++j) + for (unsigned int k = 0; k < 2; ++k) { - Point quadrature_point (i, + Point quadrature_point (j, 2.0 * edge_quadrature_points[q_point] (0) - - 1.0, j); + - 1.0, k); - this->restriction[index][i + 4 * j + 2] - ((i + 4 * j) * this->degree, dof) - += weight * this->shape_value_component (dof, - quadrature_point, - 1); + this->restriction[index][j + 4 * k + 2] + (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) - - 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, + = 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, 2.0 * edge_quadrature_points[q_point] (0) - 1.0); - this->restriction[index][i + 2 * (j + 2)] - ((i + 2 * (j + 4)) * this->degree, dof) - += weight * this->shape_value_component (dof, - quadrature_point, - 2); + 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); } + } } - // Then project the shape functions - // of the child cells to the higher - // order shape functions of the - // parent cell. + // Then interpolate the shape functions + // of the child cells to the face + // and interior 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 (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) - { - 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), - 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) - { - 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 - { - 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 k = 0; k < deg; ++k) - { - 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); - } - } - const QGauss<2> face_quadrature (2 * this->degree); const std::vector >& face_quadrature_points = face_quadrature.get_points (); @@ -1238,724 +779,268 @@ FE_Nedelec::initialize_restriction () 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 (); - - { - 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) + 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) + { + const double weight = 2.0 * face_quadrature.weight (q_point); + + for (unsigned int i = 0; i < 2; ++i) + { + 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) { - 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); + m = 0; + + if (face_quadrature_points[q_point] (1) < 0.5) + l = 0; + + else + { + quadrature_point (2) -= 1.0; + l = 1; + } } - + else { - 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); + 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; + } } - - 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) + + for (unsigned int j = 0; j < this->degree; ++j) { - 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) + 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); + } + } + } + } + + // Then, the interpolation to + // the interior shape functions. + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) { - tmp = 0.0; - + 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 { - 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); + 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 { - 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); + quadrature_point (2) -= 1.0; + l = 1; } } - - 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); - } } - + 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) + { + quadrature_point (0) -= 1.0; + n = 1; + + if (quadrature_points[q_point] (1) < 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); + m = 0; + + if (quadrature_points[q_point] (2) < 0.5) + l = 0; + + else + { + quadrature_point (2) -= 1.0; + l = 1; + } } - + 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); + 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; + } } - - 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 weight = 2.0 * quadrature.weight (q_point); + + for (unsigned int i = 0; i < this->degree; ++i) { const double L_i_0 - = legendre_polynomials[i].value - (quadrature_points[q_point] (0)); + = weight * legendre_polynomials[i].value (quadrature_points[q_point] (0)); const double L_i_1 - = legendre_polynomials[i].value - (quadrature_points[q_point] (1)); + = weight * legendre_polynomials[i].value (quadrature_points[q_point] (1)); const double L_i_2 - = legendre_polynomials[i].value - (quadrature_points[q_point] (2)); - + = weight * legendre_polynomials[i].value (quadrature_points[q_point] (2)); + for (unsigned int j = 0; j < deg; ++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)); - const double l_j_2 - = L_i_2 * lobatto_polynomials[j + 2].value - (quadrature_points[q_point] (0)); - + 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)); + for (unsigned int k = 0; k < deg; ++k) { - const double l_k_0 - = l_j_0 * lobatto_polynomials[k + 2].value - (quadrature_points[q_point] (2)); + const double L_k + = legendre_polynomials[k + 1].value (quadrature_points[q_point] (2)); const double l_k_1 - = l_j_1 * lobatto_polynomials[k + 2].value - (quadrature_points[q_point] (2)); + = lobatto_polynomials[k + 2].value (quadrature_points[q_point] (1)); const double l_k_2 - = 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; - } + = 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); } } } } - - 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); - } } } @@ -1992,6 +1077,7 @@ FE_Nedelec::get_dpo_vector (const unsigned int degree, bool dg) return dpo; } + //--------------------------------------------------------------------------- // Data field initialization //--------------------------------------------------------------------------- @@ -2507,8 +1593,9 @@ const + GeometryInfo::lines_per_face) = 1.0; interpolation_matrix - (i + (j + GeometryInfo::lines_per_face) * source_fe.degree, - i + (j + GeometryInfo::lines_per_face) * this->degree) + (i + (j + GeometryInfo::lines_per_face + source_fe.degree - 1) + * source_fe.degree, + i + (j + GeometryInfo::lines_per_face + this->degree-1) * this->degree) = 1.0; } } @@ -2530,20 +1617,6 @@ 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( @@ -2551,9 +1624,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; @@ -2567,12 +1640,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. @@ -2580,101 +1653,41 @@ 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; - // Perform projection-based interpolation - // as usual. + // Set the degrees of freedom + // by 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 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) + for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; + ++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); - + 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); + 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 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); + 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)); } } @@ -2686,143 +1699,52 @@ 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 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); - - 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) + for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; + ++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); + const double weight = 0.5 * edge_quadrature.weight (q_point); - 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) + 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 (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 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; - } + 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)); } - - 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); } + } + if (source_fe.degree > 1) + { const QGauss<2> quadrature (source_fe.degree); const std::vector >& quadrature_points = quadrature.get_points (); @@ -2830,126 +1752,59 @@ FE_Nedelec::get_subface_interpolation_matrix( lobatto_polynomials = Polynomials::Lobatto::generate_complete_basis (source_fe.degree); - const unsigned int n_boundary_dofs + const unsigned int n_edge_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, - 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 (), 2); - system_rhs.reinit (system_matrix_inv.m (), 2); - tmp.reinit (2); - - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) { - system_rhs = 0.0; - - for (unsigned int q_point = 0; - q_point < n_quadrature_points; ++q_point) + const double weight = 0.5 * quadrature.weight (q_point); + const Point + quadrature_point + (0.5 * (quadrature_points[q_point] (0) + + shifts[subface][0]), + 0.5 * (quadrature_points[q_point] (1) + + shifts[subface][1]), 0.0); + + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - Point - quadrature_point - (0.5 * (quadrature_points[q_point] (0) - + shifts[subface][0]), - 0.5 * (quadrature_points[q_point] (1) - + shifts[subface][1]), 0.0); - tmp (0) = 0.5 * this->shape_value_component - (this->face_to_cell_index (dof, 4), - quadrature_point, 0); - tmp (1) = 0.5 * this->shape_value_component - (this->face_to_cell_index (dof, 4), - quadrature_point, 1); - quadrature_point - = Point (quadrature_points[q_point] (0), - quadrature_points[q_point] (1), 0.0); - - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < source_fe.degree; ++j) - { - tmp (0) -= interpolation_matrix - ((i + 2) * source_fe.degree + j, dof) - * source_fe.shape_value_component - ((i + 2) * source_fe.degree + j, - quadrature_point, 0); - tmp (1) -= interpolation_matrix - (i * source_fe.degree + j, dof) - * source_fe.shape_value_component - (i * source_fe.degree + j, - quadrature_point, 1); - } - - tmp *= quadrature.weight (q_point); - + const double shape_grad_0 + = this->shape_grad_component (this->face_to_cell_index (dof, 4), + quadrature_point, 0)[1]; + const double shape_grad_1 + = this->shape_grad_component (this->face_to_cell_index (dof, 4), + quadrature_point, 1)[0]; + const double shape_value_0 + = this->shape_value_component (this->face_to_cell_index (dof, 4), + quadrature_point, 0); + const double shape_value_1 + = this->shape_value_component (this->face_to_cell_index (dof, 4), + quadrature_point, 1); + for (unsigned int i = 0; i < source_fe.degree; ++i) { - 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)); - + 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 < source_fe.degree - 1; ++j) { - system_rhs (i * (source_fe.degree - 1) + j, 0) - += tmp (0) * L_i_0 - * lobatto_polynomials[j + 2].value - (quadrature_points[q_point] (1)); - system_rhs (i * (source_fe.degree - 1) + j, 1) - += tmp (1) * L_i_1 - * lobatto_polynomials[j + 2].value - (quadrature_points[q_point] (0)); + interpolation_matrix (i * (source_fe.degree - 1) + j + n_edge_dofs, dof) + += L_i_0 * (shape_grad_0 + * legendre_polynomials[j + 1].value (quadrature_points[q_point] (1)) + + shape_value_0 + * lobatto_polynomials[j + 2].value (quadrature_points[q_point] (1))); + interpolation_matrix (i + (j + source_fe.degree - 1) * source_fe.degree + + n_edge_dofs, dof) + += L_i_1 * (shape_grad_1 + * legendre_polynomials[j + 1].value (quadrature_points[q_point] (0)) + + shape_value_1 + * lobatto_polynomials[j + 2].value (quadrature_points[q_point] (0))); } } } - - system_matrix_inv.mmult (solution, system_rhs); - - for (unsigned int i = 0; i < source_fe.degree; ++i) - for (unsigned int j = 0; j < source_fe.degree - 1; ++j) - { - if (std::abs (solution (i * (source_fe.degree - 1) + j, 0)) - > 1e-14) - interpolation_matrix (i * (source_fe.degree - 1) - + j + n_boundary_dofs, dof) - = solution (i * (source_fe.degree - 1) + j, 0); - - if (std::abs (solution (i * (source_fe.degree - 1) + j, 1)) - > 1e-14) - interpolation_matrix (i + (j + source_fe.degree - 1) - * source_fe.degree - + n_boundary_dofs, dof) - = solution (i * (source_fe.degree - 1) + j, 1); - } } } @@ -2975,16 +1830,21 @@ 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 ())); @@ -2996,122 +1856,157 @@ FE_Nedelec::interpolate (std::vector& local_dofs, std::fill (local_dofs.begin (), local_dofs.end (), 0.); if (offset < dim) - switch (dim) - { - case 2: - { - const QGauss<1> reference_edge_quadrature (this->degree); - const unsigned int& n_edge_points - = reference_edge_quadrature.size (); - + { + 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 - // 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); - + // 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) + { + 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); + } + } + } + // 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; - } + for (unsigned int e = 0; e < GeometryInfo::lines_per_cell; ++e) + for (unsigned int i = 0; i < this->degree; ++i) + if (std::abs (local_dofs[i + e * this->degree]) < 1e-14) + local_dofs[i + e * this->degree] = 0.0; - if (offset == 0) - for (unsigned int i = 0; i < 2; ++i) + // If the degree is greater + // than 0, then we have still + // some interpolations onto + // the interior shape + // functions left. + if (this->degree > 1) { - 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); + 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 + // interior shape + // functions. + const QGauss reference_quadrature (this->degree); + const unsigned int& n_interior_points + = reference_quadrature.size (); - if (std::abs (local_dofs[(i + 2) * this->degree]) < 1e-14) - local_dofs[(i + 2) * this->degree] = 0.0; - } + FullMatrix system_matrix (deg * this->degree, + deg * this->degree); - // 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); + 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 + (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.invert (system_matrix); + FullMatrix system_matrix_inv (system_matrix.m (), + system_matrix.m ()); + + system_matrix_inv.invert (system_matrix); - Vector system_rhs (system_matrix.m ()); - Vector solution (system_rhs.size ()); + Vector solution (system_matrix_inv.m ()); + Vector system_rhs (system_matrix.m ()); - for (unsigned int line = 0; - line < GeometryInfo::lines_per_cell; ++line) - if ((line < 2) || (offset == 0)) + if (offset == 0) { - // Set up the right hand side. - system_rhs = 0; - - for (unsigned int q_point = 0; q_point < n_edge_points; - ++q_point) + // Set up the right hand side + // for the horizontal shape + // functions. + for (unsigned int q_point = 0; + q_point < n_interior_points; ++q_point) { - 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; + 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); @@ -3120,814 +2015,715 @@ FE_Nedelec::interpolate (std::vector& local_dofs, // 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); + 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 - // 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; + // 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 < 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 + 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 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 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_grad[k].value - (this->generalized_support_points[q_point - + GeometryInfo::lines_per_cell - * n_edge_points] - (0)) - * lobatto_polynomials[l + 2].value + * lobatto_polynomials[j + 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. - 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] (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)); - } + (0)); + } - 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 < 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); - - 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 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)); - } + 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); + } - system_matrix_inv.vmult (solution, system_rhs); + break; + } + 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) + { + 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); + } + } + } + // Add the computed values // to the resulting vector // only, if they are not // too small. + for (unsigned int e = 0; e < GeometryInfo::lines_per_cell; ++e) 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; - } + if (std::abs (local_dofs[i + e * this->degree]) < 1e-14) + local_dofs[i + e * this->degree] = 0.0; - case 3: - { - const QGauss<1> - reference_edge_quadrature (this->degree); - const unsigned int& - n_edge_points = reference_edge_quadrature.size (); + // 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; - // Let us begin with the - // 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); + 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 + 8) * this->degree]) < 1e-14) - local_dofs[(i + 8) * 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 (offset + 1 < dim) - { - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j < 2; ++j) + for (unsigned int face = 0; + face < GeometryInfo::faces_per_cell; ++face) { - 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); + switch (face) + { + case 0: + { + if (offset + 1 < dim) + { + // Set up the right hand side + // for the horizontal shape + // functions. + system_rhs = 0; - // Add the computed values - // to the resulting vector - // only, if they are not - // too small. - if (std::abs (local_dofs[(i + 4 * j) * this->degree]) - < 1e-14) - local_dofs[(i + 4 * j) * this->degree] = 0.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] (1); - 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 < 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); - // 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; - } - } + 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; + } - // 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)); + system_matrix_inv.vmult (solution, system_rhs); - system_matrix_inv.invert (system_matrix); + // 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 line = 0; - line < GeometryInfo::lines_per_cell; ++line) - { - // Set up the right hand side. - system_rhs = 0; + // Set up the right hand side + // for the vertical shape + // functions. + system_rhs = 0; - 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 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 < 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 < 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 < 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 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] (1); - - 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); + 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) - system_rhs (i * (this->degree-1) + j) - += boundary_weights - (q_point + n_edge_points, - 2 * (i * (this->degree-1) + j)) * tmp; - } + break; + } - system_matrix_inv.vmult (solution, system_rhs); + case 1: + { + if (offset + 1 < dim) + { + // Set up the right hand side + // for the horizontal shape + // functions. + system_rhs = 0; - // 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 - + GeometryInfo::lines_per_cell) - * (this->degree-1) + j - + GeometryInfo::lines_per_cell] - = solution (i * (this->degree-1) + 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 + + n_face_points] (1); - // 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[(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 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; - } + 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 < 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; - } + // 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); + } - case 1: - { - if (offset + 1 < dim) - { // Set up the right hand side - // for the horizontal shape + // for the vertical 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] (1); + 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[(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 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 <= 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; - } + 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; + } - 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 <= 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); - } - - // 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 + 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 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); + 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); - // 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; - } + break; + } - case 2: - { - if (offset == 0) - { - // Set up the right hand side - // for the horizontal shape - // functions. - system_rhs = 0; + case 2: + { + if (offset == 0) + { + // 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 + 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 + 2 * n_face_points] + (2); - 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; - } + 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; + } - 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 <= 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); - } + // 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); + } // Set up the right hand side // for the vertical 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 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; - } + 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 < 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; + } - 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 <= 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; - } + 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); - case 3: - { - if (offset == 0) - { - // Set up the right hand side - // for the horizontal shape - // functions. - system_rhs = 0; + break; + } - 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); + 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 < 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; - } + 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); - system_matrix_inv.vmult (solution, system_rhs); + 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; + } - // 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 - + 6 * this->degree) * deg + j - + GeometryInfo::lines_per_cell] - = solution (i * deg + j); - } + 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); + } // Set up the right hand side // for the vertical 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 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; - } + 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 < 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; + } - 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 <= 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; - } + 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); - case 4: - { + 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 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; + } + + 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 + + 9 * deg) + * this->degree] + = solution (i * deg + j); + } + + break; + } + + default: if (offset + 1 < dim) { // Set up the right hand side @@ -3943,25 +2739,23 @@ FE_Nedelec::interpolate (std::vector& local_dofs, double tmp = values[q_point + GeometryInfo::lines_per_cell - * n_edge_points + 4 - * n_face_points] (0); + * n_edge_points + + 5 * n_face_points] (0); for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j <= deg; ++j) - tmp -= local_dofs[(i + 2) - * this->degree - + j] + for (unsigned int j = 0; j < this->degree; ++j) + tmp -= local_dofs[(i + 6) + * this->degree + j] * this->shape_value_component - ((i + 2) * this->degree - + j, + ((i + 6) * this->degree + j, this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points - + 4 + + 5 * n_face_points], 0); - for (unsigned int i = 0; i <= deg; ++i) + for (unsigned int i = 0; i < this->degree; ++i) for (unsigned int j = 0; j < deg; ++j) system_rhs (i * deg + j) += boundary_weights @@ -3972,18 +2766,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 <= deg; ++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 - + 8 * this->degree) * deg - + j - + GeometryInfo::lines_per_cell] + + 10 * this->degree) + * deg + j + + GeometryInfo::lines_per_cell] = solution (i * deg + j); } @@ -3998,22 +2792,23 @@ FE_Nedelec::interpolate (std::vector& local_dofs, double tmp = values[q_point + GeometryInfo::lines_per_cell - * n_edge_points + 4 + * n_edge_points + 5 * n_face_points] (1); for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j <= deg; ++j) - tmp -= local_dofs[i * this->degree + j] + for (unsigned int j = 0; j < this->degree; ++j) + tmp -= local_dofs[(i + 4) + * this->degree + j] * this->shape_value_component - (i * this->degree + j, + ((i + 4) * this->degree + j, this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points - + 4 + + 5 * n_face_points], 1); - for (unsigned i = 0; i <= deg; ++i) + for (unsigned i = 0; i < this->degree; ++i) for (unsigned int j = 0; j < deg; ++j) system_rhs (i * deg + j) += boundary_weights @@ -4027,508 +2822,397 @@ FE_Nedelec::interpolate (std::vector& local_dofs, // to the resulting vector // only, if they are not // too small. - for (unsigned int i = 0; i <= deg; ++i) + 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] + + 11 * deg) * this->degree] = solution (i * deg + j); } + } + } - break; - } - - default: - if (offset + 1 < dim) - { - // Set up the right hand side - // for the horizontal shape + // Finally we project + // the remaining parts + // of the function on + // the interior 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 - + 5 * n_face_points] (0); - - 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); + const QGauss reference_quadrature (this->degree); + const unsigned int& + n_interior_points = reference_quadrature.size (); - 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; - } + // We create the + // system matrix. + system_matrix.reinit (this->degree * deg * deg, + this->degree * deg * deg); + system_matrix = 0; - 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) + 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)); - // 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); - } + 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 ()); - // Set up the right hand side - // for the vertical shape - // functions. - system_rhs = 0; + 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_face_points; ++q_point) + 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) { - double tmp - = values[q_point - + GeometryInfo::lines_per_cell - * n_edge_points + 5 - * n_face_points] (1); - - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j <= deg; ++j) - tmp -= local_dofs[(i + 4) - * this->degree + j] + 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) * this->degree + j, + (i + (j + 4 * k + 2) + * this->degree, this->generalized_support_points[q_point + GeometryInfo::lines_per_cell * n_edge_points - + 5 + + GeometryInfo::faces_per_cell * n_face_points], - 1); + 0); - 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 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] + * this->shape_value_component + ((i + 2 * (k + 2) * this->degree + + GeometryInfo::lines_per_cell) + * 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); } - 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 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 - + 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 + 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)) + (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)); + (2)); + } - 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); + + // 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); + } - if (offset + 1 < dim) - { - if (offset == 0) - { // 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] (0); - - for (unsigned int i = 0; i <= deg; ++i) + 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 < 2; ++j) { - for (unsigned int j = 0; j < 2; ++j) - for (unsigned int k = 0; k < 2; ++k) - tmp -= local_dofs[i + (j + 4 * k + 2) + for (unsigned int k = 0; k < 2; ++k) + tmp -= local_dofs[i + (4 * j + k) * 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); - - 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] - * this->shape_value_component - ((i + 2 * (k + 2) * this->degree - + GeometryInfo::lines_per_cell) - * 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); - } + * 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 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 + 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)) - * lobatto_polynomials[k + 2].value - (this->generalized_support_points[q_point + * 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] - (2)); - } - - 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) - 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); - } - - // 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); - - 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); + * 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) - 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] - (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)); - } + 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)); + } - 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 <= 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); - } + 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 + * 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); + system_rhs = 0; - 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], - 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] (2); - for (unsigned int k = 0; k < deg; ++k) - tmp -= local_dofs[i + ((2 * j + 1) * deg + k - + GeometryInfo::lines_per_cell) - * this->degree] + 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 + ((2 * j + 1) * deg + k - + GeometryInfo::lines_per_cell) - * this->degree, + (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], 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)); - } + 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 < 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 + + GeometryInfo::lines_per_cell + * n_edge_points + + GeometryInfo::faces_per_cell + * n_face_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 <= 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); - } + 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); + } - break; - } + break; + } - default: - Assert (false, ExcNotImplemented ()); - } + default: + Assert (false, ExcNotImplemented ()); + } + } } @@ -4544,7 +3228,6 @@ 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 (), @@ -4554,53 +3237,57 @@ 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 - // 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]; - + // 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]; + } + } + // Add the computed values // to the resulting vector // only, if they are not // too small. - if (std::abs (local_dofs[(i + 2 * j) * this->degree]) < 1e-14) - local_dofs[(i + 2 * j) * this->degree] = 0.0; - } + for (unsigned int e = 0; e < GeometryInfo::lines_per_cell; ++e) + for (unsigned int i = 0; i < this->degree; ++i) + if (std::abs (local_dofs[i + e * this->degree]) < 1e-14) + local_dofs[i + e * 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 > 1) + // some interpolations onto + // the interior shape + // functions left. + 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 (this->degree-1, this->degree-1); std::vector > lobatto_polynomials_grad (this->degree); @@ -4609,91 +3296,26 @@ const lobatto_polynomials_grad[i] = lobatto_polynomials[i + 1].derivative (); - // 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 + // We set up the system + // matrix and use it for + // both, the horizontal + // and the vertical // interior shape - // functions. Again we - // set up the system - // matrix and use it - // for both, the - // horizontal and the - // vertical, interior - // shape functions. + // functions. const QGauss reference_quadrature (this->degree); - const unsigned int& - n_interior_points = reference_quadrature.size (); - 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; + FullMatrix system_matrix (deg * this->degree, + deg * this->degree); for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < this->degree-1; ++j) + for (unsigned int j = 0; j < deg; ++j) for (unsigned int k = 0; k < this->degree; ++k) - for (unsigned int l = 0; l < this->degree-1; ++l) + for (unsigned int l = 0; l < deg; ++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) + system_matrix (i * deg + j, k * deg + l) += reference_quadrature.weight (q_point) * legendre_polynomials[i].value (this->generalized_support_points[q_point @@ -4716,14 +3338,13 @@ const * n_edge_points] (1)); - system_matrix_inv.reinit (system_matrix.m (), - system_matrix.m ()); + FullMatrix system_matrix_inv (system_matrix.m (), + system_matrix.m ()); + system_matrix_inv.invert (system_matrix); - // Set up the right hand side - // for the horizontal shape - // functions. - system_rhs.reinit (system_matrix_inv.m ()); - system_rhs = 0; + + Vector solution (system_matrix_inv.m ()); + Vector system_rhs (system_matrix.m ()); for (unsigned int q_point = 0; q_point < n_interior_points; ++q_point) @@ -4830,61 +3451,50 @@ const case 3: { // Let us begin with the - // interpolation part. - const QGauss<1> reference_edge_quadrature (this->degree); - const unsigned int& - n_edge_points = reference_edge_quadrature.size (); - + // edge degrees of freedom. for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - 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]; + 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]; + } } - + // Add the computed values // to the resulting vector // only, if they are not // too small. - 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; + for (unsigned int e = 0; e < GeometryInfo::lines_per_cell; ++e) + for (unsigned int i = 0; i < this->degree; ++i) + if (std::abs (local_dofs[i + e * this->degree]) < 1e-14) + local_dofs[i + e * this->degree] = 0.0; // 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. + // some interpolation to the + // face and interior shape + // functions left. 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 (this->degree-1, this->degree-1); + FullMatrix system_matrix (deg * this->degree, + deg * this->degree); std::vector > lobatto_polynomials_grad (this->degree); @@ -4893,90 +3503,23 @@ const lobatto_polynomials_grad[i] = lobatto_polynomials[i + 1].derivative (); - // 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<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 + // 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; + const unsigned int + n_face_points = n_edge_points * n_edge_points; for (unsigned int i = 0; i < this->degree; ++i) - for (unsigned int j = 0; j < this->degree-1; ++j) + for (unsigned int j = 0; j < deg; ++j) for (unsigned int k = 0; k < this->degree; ++k) - for (unsigned int l = 0; l < this->degree-1; ++l) + for (unsigned int l = 0; l < deg; ++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) + system_matrix (i * deg + j, k * deg + l) += boundary_weights (q_point + n_edge_points, - 2 * (k * (this->degree-1) + l)) + 2 * (k * deg + l)) * legendre_polynomials[i].value (this->generalized_face_support_points[q_point + 4 @@ -4988,12 +3531,11 @@ const * n_edge_points] (1)); - system_matrix_inv.reinit (system_matrix.m (), - system_matrix.m ()); + FullMatrix system_matrix_inv (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}}; @@ -5001,6 +3543,8 @@ 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