From: Markus Buerg Date: Mon, 12 Jul 2010 18:56:38 +0000 (+0000) Subject: Documentation and style changes. X-Git-Tag: v8.0.0~5837 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=38c5743b875b6e09e04636a4298e30cd5ed44597;p=dealii.git Documentation and style changes. git-svn-id: https://svn.dealii.org/trunk@21471 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc index bbc5a613ef..3a845c6586 100644 --- a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc +++ b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc @@ -11,13 +11,8 @@ // //--------------------------------------------------------------------------- -#include -#include #include #include -#include -#include -#include #include #include @@ -711,121 +706,67 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, case 2: { switch (subface) { case 0: { - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 0), Point ())) > 1e-14) - interpolation_matrix (0, dof) = this->shape_value (this->face_to_cell_index (dof, 0), Point ()); - - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 0), Point (0.0, 0.5))) > 1e-14) - interpolation_matrix (1, dof) = this->shape_value (this->face_to_cell_index (dof, 0), Point (0.0, 0.5)); + interpolation_matrix (0, 0) = 1.0; + interpolation_matrix (1, 0) = 0.5; + interpolation_matrix (1, 1) = 0.5; + + for (unsigned int dof = 2; dof < this->dofs_per_face;) { + interpolation_matrix (1, dof) = -1.0; + dof = dof + 2; + } + + int faculty_i = 1; + int faculty_ij; + int faculty_j; + + for (int i = 2; i < (int) this->dofs_per_face; ++i) { + interpolation_matrix (i, i) = std::pow (0.5, i); + faculty_i *= i; + faculty_j = faculty_i; + faculty_ij = 1; + + for (int j = i + 1; j < (int) this->dofs_per_face; ++j) { + faculty_ij *= j - i; + faculty_j *= j; + + if ((i + j) & 1) + interpolation_matrix (i, j) = -1.0 * std::pow (0.5, j) * faculty_j / (faculty_i * faculty_ij); + + else + interpolation_matrix (i, j) = std::pow (0.5, j) * faculty_j / (faculty_i * faculty_ij); + } } - - if (this->degree > 1) { - double weight; - QGauss reference_edge_quadrature (this->degree); - Quadrature edge_quadrature = QProjector::project_to_child (reference_edge_quadrature, subface); - const unsigned int n_edge_points = edge_quadrature.size (); - FullMatrix assembling_matrix (this->degree - 1, dim * n_edge_points); - FullMatrix system_matrix (this->degree - 1, this->degree - 1); - FullMatrix system_matrix_inv (system_matrix.m (), system_matrix.m ()); - Point point; - std::vector > quadrature_points = edge_quadrature.get_points (); - Tensor<1, dim> grad; - Vector assembling_vector (dim * n_edge_points); - Vector solution (system_matrix.m ()); - Vector system_rhs (system_matrix.m ()); - - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (0.0, 2.0 * quadrature_points[q_point] (0)); - weight = std::sqrt (edge_quadrature.weight (q_point)); - - for (unsigned int i = 0; i < system_matrix.m (); ++i) { - grad = weight * this->shape_grad (i + 4, point); - - for (unsigned int d = 0; d < dim; ++d) - assembling_matrix (i, dim * q_point + d) = grad[d]; - } - } - - assembling_matrix.mTmult (system_matrix, assembling_matrix); - system_matrix_inv.invert (system_matrix); - - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (0.0, 2.0 * quadrature_points[q_point] (0)); - grad = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_grad (this->face_to_cell_index (dof, 0), Point (0.0, quadrature_points[q_point] (0))) - interpolation_matrix (dof, 0) * source_fe.shape_grad (0, point) - interpolation_matrix (dof, 1) * source_fe.shape_grad (1, point)); - - for (unsigned int d = 0; d < dim; ++d) - assembling_vector (dim * q_point + d) = grad[d]; - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 2, dof) = solution (i); - } - } break; } case 1: { - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 0), Point (0.0, 0.5))) > 1e-14) - interpolation_matrix (0, dof) = this->shape_value (this->face_to_cell_index (dof, 0), Point (0.0, 0.5)); - - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 0), Point (0.0, 1.0))) > 1e-14) - interpolation_matrix (1, dof) = this->shape_value (this->face_to_cell_index (dof, 0), Point (0.0, 1.0)); + interpolation_matrix (0, 0) = 0.5; + interpolation_matrix (0, 1) = 0.5; + + for (unsigned int dof = 2; dof < this->dofs_per_face;) { + interpolation_matrix (0, dof) = -1.0; + dof = dof + 2; + } + + interpolation_matrix (1, 1) = 1.0; + + int faculty_i = 1; + int faculty_ij; + int faculty_j; + + for (int i = 2; i < (int) this->dofs_per_face; ++i) { + interpolation_matrix (i, i) = std::pow (0.5, i); + faculty_i *= i; + faculty_j = faculty_i; + faculty_ij = 1; + + for (int j = i + 1; j < (int) this->dofs_per_face; ++j) { + faculty_ij *= j - i; + faculty_j *= j; + interpolation_matrix (i, j) = std::pow (0.5, j) * faculty_j / (faculty_i * faculty_ij); + } } - - if (this->degree > 1) { - double weight; - QGauss reference_edge_quadrature (this->degree); - Quadrature edge_quadrature = QProjector::project_to_child (reference_edge_quadrature, subface); - const unsigned int n_edge_points = edge_quadrature.size (); - FullMatrix assembling_matrix (this->degree - 1, dim * n_edge_points); - FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); - FullMatrix system_matrix_inv (system_matrix.m (), system_matrix.m ()); - Point point; - std::vector > quadrature_points = edge_quadrature.get_points (); - Tensor<1, dim> grad; - Vector assembling_vector (dim * n_edge_points); - Vector solution (system_matrix.m ()); - Vector system_rhs (system_matrix.m ()); - - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (0.0, 2.0 * quadrature_points[q_point] (0) - 1.0); - weight = std::sqrt (edge_quadrature.weight (q_point)); - - for (unsigned int i = 0; i < system_matrix.m (); ++i) { - grad = weight * this->shape_grad (i + 4, point); - - for (unsigned int d = 0; d < dim; ++d) - assembling_matrix (i, dim * q_point + d) = grad[d]; - } - } - - assembling_matrix.mTmult (system_matrix, assembling_matrix); - system_matrix_inv.invert (system_matrix); - - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (0.0, 2.0 * quadrature_points[q_point] (0) - 1.0); - grad = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_grad (this->face_to_cell_index (dof, 0), Point (0.0, quadrature_points[q_point] (0))) - interpolation_matrix (0, dof) * source_fe.shape_grad (0, point) - interpolation_matrix (1, dof) * source_fe.shape_grad (1, point)); - - for (unsigned int d = 0; d < dim; ++d) - assembling_vector (dim * q_point + d) = grad[d]; - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 2, dof) = solution (i); - } - } } } @@ -835,649 +776,511 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, case 3: { switch (subface) { case 0: { - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point ())) > 1e-14) - interpolation_matrix (0, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point ()); - - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.0, 0.0))) > 1e-14) - interpolation_matrix (1, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.0, 0.0)); + interpolation_matrix (0, 0) = 1.0; + interpolation_matrix (1, 0) = 0.5; + interpolation_matrix (1, 1) = 0.5; + interpolation_matrix (2, 0) = 0.5; + interpolation_matrix (2, 2) = 0.5; + + for (unsigned int i = 0; i < this->degree - 1;) { + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) + interpolation_matrix (3, i + line * (this->degree - 1) + 4) = -0.5; - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 0.5, 0.0))) > 1e-14) - interpolation_matrix (2, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 0.5, 0.0)); + for (unsigned int j = 0; j < this->degree - 1;) { + interpolation_matrix (3, i + (j + 4) * this->degree - j) = 1.0; + j = j + 2; + } - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0))) > 1e-14) - interpolation_matrix (3, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0)); + interpolation_matrix (1, i + 2 * (this->degree + 1)) = -1.0; + interpolation_matrix (2, i + 4) = -1.0; + i = i + 2; } - if (this->degree > 1) { - double weight; - QGauss reference_edge_quadrature (this->degree + 1); - Quadrature edge_quadrature = QProjector::project_to_child (reference_edge_quadrature, 0); - const unsigned int n_edge_points = reference_edge_quadrature.size (); - QGauss reference_quadrature (this->degree); - Quadrature quadrature = QProjector::project_to_child (reference_quadrature, subface); - const unsigned int n_quadrature_points = quadrature.size (); - FullMatrix assembling_matrix (this->degree - 1, n_edge_points); - FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); - FullMatrix system_matrix_inv (assembling_matrix.m (), assembling_matrix.m ()); - Point point; - PreconditionJacobi > precondition; - std::vector > edge_quadrature_points = edge_quadrature.get_points (); - std::vector > quadrature_points = quadrature.get_points (); - Tensor<1, dim> grad; - Vector assembling_vector (n_edge_points); - Vector solution (assembling_matrix.m ()); - Vector system_rhs (assembling_matrix.m ()); + for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_face; ++vertex) + interpolation_matrix (3, vertex) = 0.25; + + double tmp; + int faculty_i = 1; + int faculty_ij; + int faculty_j; + int faculty_k; + int faculty_kl; + int faculty_l; + + for (int i = 2; i <= (int) this->degree; ++i) { + tmp = std::pow (0.5, i); + interpolation_matrix (i + 2, i + 2) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, i + 2 * this->degree) = tmp; + tmp *= 0.5; + interpolation_matrix (i + source_fe.degree + 1, i + 2) = tmp; + interpolation_matrix (i + source_fe.degree + 1, i + this->degree + 1) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, i + 2 * this->degree) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, i + 3 * this->degree - 1) = tmp; + tmp *= -2.0; + + for (unsigned int j = 0; j < this->degree - 1;) { + interpolation_matrix (i + source_fe.degree + 1, (i + 2) * this->degree + j + 2 - i) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, i + (j + 4) * this->degree - j - 2) = tmp; + j = j + 2; + } - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (0.0, 2.0 * edge_quadrature_points[q_point] (0), 0.0); - weight = std::sqrt (edge_quadrature.weight (q_point)); + faculty_k = 1; + + for (int j = 2; j <= (int) this->degree; ++j) { + interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j); + faculty_k *= j; + faculty_kl = 1; + faculty_l = faculty_k; + + for (int k = j + 1; k < (int) this->degree; ++k) { + faculty_kl *= k - j; + faculty_l *= k; + + if ((j + k) & 1) + interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = -1.0 * std::pow (0.5, i + k) * faculty_l / (faculty_k * faculty_kl); - for (unsigned int i = 0; i < system_matrix.m (); ++i) - assembling_matrix (i, q_point) = weight * this->shape_value (i + 8, point); + else + interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = std::pow (0.5, i + k) * faculty_l / (faculty_k * faculty_kl); } + } + + faculty_i *= i; + faculty_j = faculty_i; + faculty_ij = 1; + + for (int j = i + 1; j <= (int) this->degree; ++j) { + faculty_ij *= j - i; + faculty_j *= j; - assembling_matrix.mTmult (system_matrix, assembling_matrix); - system_matrix_inv.invert (system_matrix); - - for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) - switch (line) { - case 0: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (0.0, 2.0 * edge_quadrature_points[q_point] (0), 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, edge_quadrature_points[q_point] (0), 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 4, dof) = solution (i); - - break; - } - - case 1: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (1.0, 2.0 * edge_quadrature_points[q_point] (0), 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, edge_quadrature_points[q_point] (0), 0.0)) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + source_fe.degree + 3, dof) = solution (i); - - break; - } - - case 2: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (2.0 * edge_quadrature_points[q_point] (0), 0.0, 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_points[q_point] (0), 0.0, 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 2 * source_fe.degree + 2, dof) = solution (i); - - break; - } + if ((i + j) & 1) { + tmp = -1.0 * std::pow (0.5, j) * faculty_j / (faculty_i * faculty_ij); + interpolation_matrix (i + 2, j + 2) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp; + faculty_k = 1; + + for (int k = 2; k <= (int) this->degree; ++k) { + interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k); + faculty_k *= k; + faculty_l = faculty_k; + faculty_kl = 1; - case 3: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (2.0 * edge_quadrature_points[q_point] (0), 1.0, 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_points[q_point] (0), 0.5, 0.0)) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 3 * source_fe.degree + 1, dof) = solution (i); + for (int l = k + 1; l <= (int) this->degree; ++l) { + faculty_kl *= l - k; + faculty_l *= l; + + if ((k + l) & 1) + interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = -1.0 * tmp * std::pow (0.5, l) * faculty_l / (faculty_k * faculty_kl); + + else + interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = tmp * std::pow (0.5, l) * faculty_l / (faculty_k * faculty_kl); } } - - assembling_matrix.reinit ((this->degree - 1) * (this->degree - 1), dim * n_quadrature_points); - assembling_vector.reinit (assembling_matrix.n ()); - - for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { - point = Point (2.0 * quadrature_points[q_point] (0), 2.0 * quadrature_points[q_point] (1), 0.0); - grad = this->shape_grad (this->face_to_cell_index (dof, 4), Point (quadrature_points[q_point] (0), quadrature_points[q_point] (1), 0.0)); - for (unsigned int vertex = 0; vertex < GeometryInfo::vertices_per_face; ++vertex) - grad -= interpolation_matrix (vertex, dof) * source_fe.shape_grad (vertex, point); + tmp *= 0.5; + interpolation_matrix (i + source_fe.degree + 1, j + 2) = tmp; + interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, j + 2 * this->degree) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp; + tmp *= -2.0; - for (unsigned int line = 0; line < GeometryInfo::lines_per_face; ++line) - for (unsigned int i = 0; i < this->degree - 1; ++i) - grad -= interpolation_matrix (i + line * (source_fe.degree - 1) + 4, dof) * source_fe.shape_grad (i + line * (source_fe.degree - 1) + 8, point); - - weight = std::sqrt (quadrature.weight (q_point)); - grad *= weight; - - for (unsigned int d = 0; d < dim; ++d) - assembling_vector (dim * q_point + d) = grad[d]; - - for (unsigned int i = 0; i < assembling_matrix.m (); ++i) { - grad = weight * this->shape_grad (i + 4 * ((this->degree + 2) * (this->degree - 1) + 2), point); - - for (unsigned int d = 0; d < dim; ++d) - assembling_matrix (i, dim * q_point + d) = grad[d]; + for (unsigned int k = 0; k < this->degree - 1;) { + interpolation_matrix (i + source_fe.degree + 1, (j + 2) * this->degree + k + 2 - j) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, j + (k + 4) * this->degree - k - 2) = tmp; + k = k + 2; } } - system_matrix.reinit (assembling_matrix.m (), assembling_matrix.n ()); - assembling_matrix.mTmult (system_matrix, assembling_matrix); - system_rhs.reinit (system_matrix.m ()); - assembling_matrix.vmult (system_rhs, assembling_vector); - solution.reinit (system_matrix.m ()); - - SolverControl solver_control (system_matrix.m (), 1e-13, false, false); - SolverCG<> cg (solver_control); - - precondition.initialize (system_matrix); - cg.solve (system_matrix, solution, system_rhs, precondition); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 4 * source_fe.degree, dof) = solution (i); - } - } - - break; - } - - case 1: { - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point ())) > 1e-14) - interpolation_matrix (0, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.0, 0.0)); - - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.0, 0.0))) > 1e-14) - interpolation_matrix (1, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (1.0, 0.0, 0.0)); - - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 0.5, 0.0))) > 1e-14) - interpolation_matrix (2, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0)); - - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0))) > 1e-14) - interpolation_matrix (3, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (1.0, 0.5, 0.0)); - } - - if (this->degree > 1) { - double weight; - QGauss reference_edge_quadrature (this->degree + 1); - Quadrature edge_quadrature_x = QProjector::project_to_child (reference_edge_quadrature, 1); - Quadrature edge_quadrature_y = QProjector::project_to_child (reference_edge_quadrature, 0); - const unsigned int n_edge_points = reference_edge_quadrature.size (); - QGauss reference_quadrature (this->degree); - Quadrature quadrature = QProjector::project_to_child (reference_quadrature, subface); - const unsigned int n_quadrature_points = quadrature.size (); - FullMatrix assembling_matrix (this->degree - 1, n_edge_points); - FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); - FullMatrix system_matrix_inv (assembling_matrix.m (), assembling_matrix.m ()); - Point point; - PreconditionJacobi > precondition; - std::vector > edge_quadrature_x_points = edge_quadrature_x.get_points (); - std::vector > edge_quadrature_y_points = edge_quadrature_y.get_points (); - std::vector > quadrature_points = quadrature.get_points (); - Tensor<1, dim> grad; - Vector assembling_vector (n_edge_points); - Vector solution (assembling_matrix.m ()); - Vector system_rhs (assembling_matrix.m ()); - - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (0.0, 2.0 * edge_quadrature_y_points[q_point] (0), 0.0); - weight = std::sqrt (edge_quadrature_y.weight (q_point)); + else { + tmp = std::pow (0.5, j) * faculty_j / (faculty_i * faculty_ij); + interpolation_matrix (i + 2, j + 2) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp; + faculty_k = 1; - for (unsigned int i = 0; i < system_matrix.m (); ++i) - assembling_matrix (i, q_point) = weight * this->shape_value (i + 8, point); - } - - assembling_matrix.mTmult (system_matrix, assembling_matrix); - system_matrix_inv.invert (system_matrix); - - for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) - switch (line) { - case 0: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (0.0, 2.0 * edge_quadrature_y_points[q_point] (0), 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature_y.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, edge_quadrature_y_points[q_point] (0), 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 4, dof) = solution (i); - - break; - } - - case 1: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (1.0, 2.0 * edge_quadrature_y_points[q_point] (0), 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature_y.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (1.0, edge_quadrature_y_points[q_point] (0), 0.0)) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + source_fe.degree + 3, dof) = solution (i); - - break; - } - - case 2: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (2.0 * edge_quadrature_x_points[q_point] (0) - 1.0, 0.0, 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature_x.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_x_points[q_point] (0), 0.0, 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 2 * source_fe.degree + 2, dof) = solution (i); - - break; - } + for (int k = 2; k <= (int) this->degree; ++k) { + interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k); + faculty_k *= k; + faculty_l = faculty_k; + faculty_kl = 1; - case 3: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (2.0 * edge_quadrature_x_points[q_point] (0) - 1.0, 1.0, 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature_x.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_x_points[q_point] (0), 0.5, 0.0)) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 3 * source_fe.degree + 1, dof) = solution (i); + for (int l = k + 1; l <= (int) this->degree; ++l) { + faculty_kl *= l - k; + faculty_l *= l; + + if ((k + l) & 1) + interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = -1.0 * tmp * std::pow (0.5, l) * faculty_l / (faculty_k * faculty_kl); + + else + interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = tmp * std::pow (0.5, l) * faculty_l / (faculty_k * faculty_kl); } } - - assembling_matrix.reinit ((this->degree - 1) * (this->degree - 1), dim * n_quadrature_points); - assembling_vector.reinit (assembling_matrix.n ()); - - for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { - point = Point (2.0 * quadrature_points[q_point] (0) - 1.0, 2.0 * quadrature_points[q_point] (1), 0.0); - grad = this->shape_grad (this->face_to_cell_index (dof, 4), Point (quadrature_points[q_point] (0), quadrature_points[q_point] (1), 0.0)); - for (unsigned int vertex = 0; vertex < GeometryInfo::vertices_per_face; ++vertex) - grad -= interpolation_matrix (vertex, dof) * source_fe.shape_grad (vertex, point); + tmp *= 0.5; + interpolation_matrix (i + source_fe.degree + 1, j + 2) = tmp; + interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, j + 2 * this->degree) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp; + tmp *= -2.0; - for (unsigned int line = 0; line < GeometryInfo::lines_per_face; ++line) - for (unsigned int i = 0; i < this->degree - 1; ++i) - grad -= interpolation_matrix (i + line * (source_fe.degree - 1) + 4, dof) * source_fe.shape_grad (i + line * (source_fe.degree - 1) + 8, point); - - weight = std::sqrt (quadrature.weight (q_point)); - grad *= weight; - - for (unsigned int d = 0; d < dim; ++d) - assembling_vector (dim * q_point + d) = grad[d]; - - for (unsigned int i = 0; i < assembling_matrix.m (); ++i) { - grad = weight * this->shape_grad (i + 4 * ((this->degree + 2) * (this->degree - 1) + 2), point); - - for (unsigned int d = 0; d < dim; ++d) - assembling_matrix (i, dim * q_point + d) = grad[d]; + for (unsigned int k = 0; k < this->degree - 1;) { + interpolation_matrix (i + source_fe.degree + 1, (j + 2) * this->degree + k + 2 - j) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, j + (k + 4) * this->degree - k - 2) = tmp; + k = k + 2; } } - - system_matrix.reinit (assembling_matrix.m (), assembling_matrix.n ()); - assembling_matrix.mTmult (system_matrix, assembling_matrix); - system_rhs.reinit (system_matrix.m ()); - assembling_matrix.vmult (system_rhs, assembling_vector); - solution.reinit (system_matrix.m ()); - - SolverControl solver_control (system_matrix.m (), 1e-13, false, false); - SolverCG<> cg (solver_control); - - precondition.initialize (system_matrix); - cg.solve (system_matrix, solution, system_rhs, precondition); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 4 * source_fe.degree, dof) = solution (i); } } break; } - case 2: { - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point ())) > 1e-14) - interpolation_matrix (0, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 0.5, 0.0)); - - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.0, 0.0))) > 1e-14) - interpolation_matrix (1, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0)); - - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 0.5, 0.0))) > 1e-14) - interpolation_matrix (2, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 1.0, 0.0)); - - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0))) > 1e-14) - interpolation_matrix (3, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 1.0, 0.0)); + case 1: { + interpolation_matrix (0, 0) = 0.5; + interpolation_matrix (0, 1) = 0.5; + interpolation_matrix (1, 1) = 1.0; + interpolation_matrix (3, 1) = 0.5; + interpolation_matrix (3, 3) = 0.5; + + for (unsigned int i = 0; i < this->degree - 1;) { + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) + interpolation_matrix (2, i + line * (this->degree - 1) + 4) = -0.5; + + for (unsigned int j = 0; j < this->degree - 1;) { + interpolation_matrix (2, i + (j + 4) * this->degree - j) = 1.0; + j = j + 2; + } + + interpolation_matrix (0, i + 2 * (this->degree + 1)) = -1.0; + interpolation_matrix (3, i + this->degree + 3) = -1.0; + i = i + 2; } - if (this->degree > 1) { - double weight; - QGauss reference_edge_quadrature (this->degree + 1); - Quadrature edge_quadrature_x = QProjector::project_to_child (reference_edge_quadrature, 0); - Quadrature edge_quadrature_y = QProjector::project_to_child (reference_edge_quadrature, 1); - const unsigned int n_edge_points = reference_edge_quadrature.size (); - QGauss reference_quadrature (this->degree); - Quadrature quadrature = QProjector::project_to_child (reference_quadrature, subface); - const unsigned int n_quadrature_points = quadrature.size (); - FullMatrix assembling_matrix (this->degree - 1, n_edge_points); - FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); - FullMatrix system_matrix_inv (assembling_matrix.m (), assembling_matrix.m ()); - Point point; - PreconditionJacobi > precondition; - std::vector > edge_quadrature_x_points = edge_quadrature_x.get_points (); - std::vector > edge_quadrature_y_points = edge_quadrature_y.get_points (); - std::vector > quadrature_points = quadrature.get_points (); - Tensor<1, dim> grad; - Vector assembling_vector (n_edge_points); - Vector solution (assembling_matrix.m ()); - Vector system_rhs (assembling_matrix.m ()); + for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_face; ++vertex) + interpolation_matrix (2, vertex) = 0.25; + + double tmp; + int faculty_i = 1; + int faculty_ij; + int faculty_j; + int faculty_k; + int faculty_kl; + int faculty_l; + + for (int i = 2; i <= (int) this->degree; ++i) { + tmp = std::pow (0.5, i + 1); + interpolation_matrix (i + 2, i + 2) = tmp; + interpolation_matrix (i + 2, i + this->degree + 1) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, i + 2 * this->degree) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, i + 3 * this->degree - 1) = tmp; + tmp *= -2.0; - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (0.0, 2.0 * edge_quadrature_y_points[q_point] (0) - 1.0, 0.0); - weight = std::sqrt (edge_quadrature_y.weight (q_point)); - - for (unsigned int i = 0; i < system_matrix.m (); ++i) - assembling_matrix (i, q_point) = weight * this->shape_value (i + 8, point); - } - - assembling_matrix.mTmult (system_matrix, assembling_matrix); - system_matrix_inv.invert (system_matrix); + for (unsigned int j = 0; j < this->degree - 1;) { + interpolation_matrix (i + 2, j + (i + 2) * this->degree + 2 - i) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, i + (j + 4) * this->degree - j - 2) = tmp; + j = j + 2; + } + + tmp *= - 1.0; + interpolation_matrix (i + source_fe.degree + 1, i + this->degree + 1) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, i + 2 * this->degree) = tmp; + faculty_i *= i; + faculty_j = faculty_i; + faculty_ij = 1; + + for (int j = i + 1; j <= (int) this->degree; ++j) { + faculty_ij *= j - i; + faculty_j *= j; + tmp = std::pow (0.5, j) * faculty_j / (faculty_i * faculty_ij); + interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp; + faculty_k = 1; - for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) - switch (line) { - case 0: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (0.0, 2.0 * edge_quadrature_y_points[q_point] (0) - 1.0, 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature_y.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, edge_quadrature_y_points[q_point] (0), 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 4, dof) = solution (i); - - break; - } - - case 1: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (1.0, 2.0 * edge_quadrature_y_points[q_point] (0) - 1.0, 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature_y.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, edge_quadrature_y_points[q_point] (0), 0.0)) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + source_fe.degree + 3, dof) = solution (i); - - break; - } + for (int k = 2; k <= (int) this->degree; ++k) { + interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k); + faculty_k *= k; + faculty_l = faculty_k; + faculty_kl = 1; + + for (int l = k + 1; l <= (int) this->degree; ++l) { + faculty_kl *= l - k; + faculty_l *= l; - case 2: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (2.0 * edge_quadrature_x_points[q_point] (0), 0.0, 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature_x.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_x_points[q_point] (0), 0.5, 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 2 * source_fe.degree + 2, dof) = solution (i); - - break; - } + if ((k + l) & 1) + interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = -1.0 * tmp * std::pow (0.5, l) * faculty_l / (faculty_k * faculty_kl); - case 3: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (2.0 * edge_quadrature_x_points[q_point] (0), 1.0, 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature_x.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_x_points[q_point] (0), 1.0, 0.0)) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 3 * source_fe.degree + 1, dof) = solution (i); - } + else + interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = tmp * std::pow (0.5, l) * faculty_l / (faculty_k * faculty_kl); } + } - assembling_matrix.reinit ((this->degree - 1) * (this->degree - 1), dim * n_quadrature_points); - assembling_vector.reinit (assembling_matrix.n ()); + tmp *= -1.0; - for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { - point = Point (2.0 * quadrature_points[q_point] (0), 2.0 * quadrature_points[q_point] (1) - 1.0, 0.0); - grad = this->shape_grad (this->face_to_cell_index (dof, 4), Point (quadrature_points[q_point] (0), quadrature_points[q_point] (1), 0.0)); - - for (unsigned int vertex = 0; vertex < GeometryInfo::vertices_per_face; ++vertex) - grad -= interpolation_matrix (vertex, dof) * source_fe.shape_grad (vertex, point); - - for (unsigned int line = 0; line < GeometryInfo::lines_per_face; ++line) - for (unsigned int i = 0; i < this->degree - 1; ++i) - grad -= interpolation_matrix (i + line * (source_fe.degree - 1) + 4, dof) * source_fe.shape_grad (i + line * (source_fe.degree - 1) + 8, point); - - weight = std::sqrt (quadrature.weight (q_point)); - grad *= weight; - - for (unsigned int d = 0; d < dim; ++d) - assembling_vector (dim * q_point + d) = grad[d]; - - for (unsigned int i = 0; i < assembling_matrix.m (); ++i) { - grad = weight * this->shape_grad (i + 4 * ((this->degree + 2) * (this->degree - 1) + 2), point); - - for (unsigned int d = 0; d < dim; ++d) - assembling_matrix (i, dim * q_point + d) = grad[d]; - } + for (unsigned int k = 0; k < this->degree - 1;) { + interpolation_matrix (i + 3 * source_fe.degree - 1, j + (k + 4) * this->degree - k - 2) = tmp; + k = k + 2; } - system_matrix.reinit (assembling_matrix.m (), assembling_matrix.n ()); - assembling_matrix.mTmult (system_matrix, assembling_matrix); - system_rhs.reinit (system_matrix.m ()); - assembling_matrix.vmult (system_rhs, assembling_vector); - solution.reinit (system_matrix.m ()); - - SolverControl solver_control (system_matrix.m (), 1e-13, false, false); - SolverCG<> cg (solver_control); + tmp *= -0.5; + interpolation_matrix (i + 3 * source_fe.degree - 1, j + 2 * this->degree) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp; - precondition.initialize (system_matrix); - cg.solve (system_matrix, solution, system_rhs, precondition); + if ((i + j) & 1) + tmp *= -1.0; + + interpolation_matrix (i + 2, j + 2) = tmp; + interpolation_matrix (i + 2, j + this->degree + 1) = tmp; + interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = 2.0 * tmp; + tmp *= -2.0; - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 4 * source_fe.degree, dof) = solution (i); + for (unsigned int k = 0; k < this->degree - 1;) { + interpolation_matrix (i + 2, k + (j + 2) * this->degree + 2 - j) = tmp; + k = k + 2; + } } - } - - break; - } - - case 3: { - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point ())) > 1e-14) - interpolation_matrix (0, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0)); - - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.0, 0.0))) > 1e-14) - interpolation_matrix (1, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (1.0, 0.5, 0.0)); - - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.0, 0.5, 0.0))) > 1e-14) - interpolation_matrix (2, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 1.0, 0.0)); - if (std::abs (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, 0.5, 0.0))) > 1e-14) - interpolation_matrix (3, dof) = this->shape_value (this->face_to_cell_index (dof, 4), Point (1.0, 1.0, 0.0)); - } - - if (this->degree > 1) { - double weight; - QGauss reference_edge_quadrature (this->degree + 1); - Quadrature edge_quadrature = QProjector::project_to_child (reference_edge_quadrature, 1); - const unsigned int n_edge_points = reference_edge_quadrature.size (); - QGauss reference_quadrature (this->degree); - Quadrature quadrature = QProjector::project_to_child (reference_quadrature, subface); - const unsigned int n_quadrature_points = quadrature.size (); - FullMatrix assembling_matrix (this->degree - 1, n_edge_points); - FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); - FullMatrix system_matrix_inv (assembling_matrix.m (), assembling_matrix.m ()); - Point point; - PreconditionJacobi > precondition; - std::vector > edge_quadrature_points = edge_quadrature.get_points (); - std::vector > quadrature_points = quadrature.get_points (); - Tensor<1, dim> grad; - Vector assembling_vector (n_edge_points); - Vector solution (assembling_matrix.m ()); - Vector system_rhs (assembling_matrix.m ()); + faculty_k = 1; - for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (0.0, 2.0 * edge_quadrature_points[q_point] (0) - 1.0, 0.0); - weight = std::sqrt (edge_quadrature.weight (q_point)); - - for (unsigned int i = 0; i < system_matrix.m (); ++i) - assembling_matrix (i, q_point) = weight * this->shape_value (i + 8, point); - } - - assembling_matrix.mTmult (system_matrix, assembling_matrix); - system_matrix_inv.invert (system_matrix); - - for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) - switch (line) { - case 0: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (0.0, 2.0 * edge_quadrature_points[q_point] (0) - 1.0, 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (0.5, edge_quadrature_points[q_point] (0), 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 4, dof) = solution (i); - - break; - } - - case 1: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (1.0, 2.0 * edge_quadrature_points[q_point] (0) - 1.0, 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (1.0, edge_quadrature_points[q_point] (0), 0.0)) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + source_fe.degree + 3, dof) = solution (i); - - break; - } - - case 2: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (2.0 * edge_quadrature_points[q_point] (0) - 1.0, 0.0, 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_points[q_point] (0), 0.5, 0.0)) - interpolation_matrix (0, dof) * source_fe.shape_value (0, point) - interpolation_matrix (1, dof) * source_fe.shape_value (1, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 2 * source_fe.degree + 2, dof) = solution (i); - - break; - } - - case 3: { - for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) { - point = Point (2.0 * edge_quadrature_points[q_point] (0) - 1.0, 1.0, 0.0); - assembling_vector (q_point) = std::sqrt (edge_quadrature.weight (q_point)) * (this->shape_value (this->face_to_cell_index (dof, 4), Point (edge_quadrature_points[q_point] (0), 1.0, 0.0)) - interpolation_matrix (2, dof) * source_fe.shape_value (2, point) - interpolation_matrix (3, dof) * source_fe.shape_value (3, point)); - } - - assembling_matrix.vmult (system_rhs, assembling_vector); - system_matrix_inv.vmult (solution, system_rhs); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 3 * source_fe.degree + 1, dof) = solution (i); - } - } - - assembling_matrix.reinit ((this->degree - 1) * (this->degree - 1), dim * n_quadrature_points); - assembling_vector.reinit (assembling_matrix.n ()); + for (int j = 2; j <= (int) this->degree; ++j) { + interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j); + faculty_k *= j; + faculty_l = faculty_k; + faculty_kl = 1; - for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { - point = Point (2.0 * quadrature_points[q_point] (0) - 1.0, 2.0 * quadrature_points[q_point] (1) - 1.0, 0.0); - grad = this->shape_grad (this->face_to_cell_index (dof, 4), Point (quadrature_points[q_point] (0), quadrature_points[q_point] (1), 0.0)); + for (int k = j + 1; k <= (int) this->degree; ++k) { + faculty_kl *= k - j; + faculty_l *= k; - for (unsigned int vertex = 0; vertex < GeometryInfo::vertices_per_face; ++vertex) - grad -= interpolation_matrix (vertex, dof) * source_fe.shape_grad (vertex, point); + if ((j + k) & 1) + interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = -1.0 * std::pow (0.5, i + k) * faculty_l / (faculty_k * faculty_kl); - for (unsigned int line = 0; line < GeometryInfo::lines_per_face; ++line) - for (unsigned int i = 0; i < this->degree - 1; ++i) - grad -= interpolation_matrix (i + line * (source_fe.degree - 1) + 4, dof) * source_fe.shape_grad (i + line * (source_fe.degree - 1) + 8, point); - - weight = std::sqrt (quadrature.weight (q_point)); - grad *= weight; - - for (unsigned int d = 0; d < dim; ++d) - assembling_vector (dim * q_point + d) = grad[d]; - - for (unsigned int i = 0; i < assembling_matrix.m (); ++i) { - grad = weight * this->shape_grad (i + 4 * ((this->degree + 2) * (this->degree - 1) + 2), point); - - for (unsigned int d = 0; d < dim; ++d) - assembling_matrix (i, dim * q_point + d) = grad[d]; - } + else + interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = std::pow (0.5, i + k) * faculty_l / (faculty_k * faculty_kl); } - - system_matrix.reinit (assembling_matrix.m (), assembling_matrix.n ()); - assembling_matrix.mTmult (system_matrix, assembling_matrix); - system_rhs.reinit (system_matrix.m ()); - assembling_matrix.vmult (system_rhs, assembling_vector); - solution.reinit (system_matrix.m ()); - - SolverControl solver_control (system_matrix.m (), 1e-13, false, false); - SolverCG<> cg (solver_control); - - precondition.initialize (system_matrix); - cg.solve (system_matrix, solution, system_rhs, precondition); - - for (unsigned int i = 0; i < solution.size (); ++i) - if (std::abs (solution (i)) > 1e-14) - interpolation_matrix (i + 4 * source_fe.degree, dof) = solution (i); } } + + break; + } + + case 2: { + interpolation_matrix (0, 0) = 0.5; + interpolation_matrix (0, 2) = 0.5; + interpolation_matrix (2, 2) = 1.0; + interpolation_matrix (3, 2) = 0.5; + interpolation_matrix (3, 3) = 0.5; + + for (unsigned int i = 0; i < this->degree - 1;) { + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) + interpolation_matrix (1, i + line * (this->degree - 1) + 4) = -0.5; + + for (unsigned int j = 0; j < this->degree - 1;) { + interpolation_matrix (1, i + (j + 4) * this->degree - j) = 1.0; + j = j + 2; + } + + interpolation_matrix (0, i + 4) = -1.0; + interpolation_matrix (3, i + 3 * this->degree + 1) = -1.0; + i = i + 2; + } + + for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_face; ++vertex) + interpolation_matrix (1, vertex) = 0.25; + + double tmp; + int faculty_i = 1; + int faculty_ij; + int faculty_j; + int faculty_k; + int faculty_kl; + int faculty_l; + + for (int i = 2; i <= (int) this->degree; ++i) { + tmp = std::pow (0.5, i); + interpolation_matrix (i + 2, i + 2) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, i + 3 * this->degree - 1) = tmp; + tmp *= 0.5; + interpolation_matrix (i + source_fe.degree + 1, i + 2) = tmp; + interpolation_matrix (i + source_fe.degree + 1, i + this->degree + 1) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, i + 2 * this->degree) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, i + 3 * this->degree - 1) = tmp; + tmp *= -2.0; + + for (unsigned int j = 0; j < this->degree - 1;) { + interpolation_matrix (i + source_fe.degree + 1, j + (i + 2) * this->degree + 2 - i) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, i + (j + 4) * this->degree - j - 2) = tmp; + j = j + 2; + } + + faculty_k = 1; + + for (int j = 2; j <= (int) this->degree; ++j) { + interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j); + faculty_k *= j; + faculty_l = faculty_k; + faculty_kl = 1; + + for (int k = j + 1; k <= (int) this->degree; ++k) { + faculty_kl *= k - j; + faculty_l *= k; + interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = std::pow (0.5, i + k) * faculty_l / (faculty_k * faculty_kl); + } + } + + faculty_i *= i; + faculty_j = faculty_i; + faculty_ij = 1; + + for (int j = i + 1; j <= (int) this->degree; ++j) { + faculty_ij *= j - i; + faculty_j *= j; + tmp = std::pow (0.5, j) * faculty_j / (faculty_i * faculty_ij); + interpolation_matrix (i + 2, j + 2) = tmp; + tmp *= -1.0; + + for (unsigned int k = 0; k < this->degree - 1;) { + interpolation_matrix (i + source_fe.degree + 1, k + (j + 2) * this->degree + 2 - j) = tmp; + k = k + 2; + } + + tmp *= -0.5; + interpolation_matrix (i + source_fe.degree + 1, j + 2) = tmp; + interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = tmp; + + if ((i + j) & 1) + tmp *= -1.0; + + interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, j + 3 * this->degree - 1) = tmp; + tmp *= 2.0; + interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp; + faculty_k = 1; + + for (int k = 2; k <= (int) this->degree; ++k) { + interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k); + faculty_k *= k; + faculty_l = faculty_k; + faculty_kl = 1; + + for (int l = k + 1; l <= (int) this->degree; ++l) { + faculty_kl *= l - k; + faculty_l *= l; + interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = tmp * std::pow (0.5, l) * faculty_l / (faculty_k * faculty_kl); + } + } + + tmp *= -1.0; + + for (unsigned int k = 0; k < this->degree - 1;) { + interpolation_matrix (i + 2 * source_fe.degree, j + (k + 4) * this->degree - k - 2) = tmp; + k = k + 2; + } + } + } + + break; + } + + case 3: { + for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_face; ++vertex) + interpolation_matrix (0, vertex) = 0.25; + + for (unsigned int i = 0; i < this->degree - 1;) { + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) + interpolation_matrix (0, i + line * (this->degree - 1) + 4) = -0.5; + + for (unsigned int j = 0; j < this->degree - 1;) { + interpolation_matrix (0, i + (j + 4) * this->degree - j) = 1.0; + j = j + 2; + } + + interpolation_matrix (1, i + 4) = -1.0; + interpolation_matrix (2, i + 3 * this->degree + 1) = -1.0; + i = i + 2; + } + + interpolation_matrix (1, 0) = 0.5; + interpolation_matrix (1, 1) = 0.5; + interpolation_matrix (2, 2) = 0.5; + interpolation_matrix (2, 3) = 0.5; + interpolation_matrix (3, 3) = 1.0; + + double tmp; + int faculty_i = 1; + int faculty_ij; + int faculty_j; + int faculty_k; + int faculty_kl; + int faculty_l; + + for (int i = 2; i <= (int) this->degree; ++i) { + tmp = std::pow (0.5, i + 1); + interpolation_matrix (i + 2, i + 2) = tmp; + interpolation_matrix (i + 2, i + this->degree + 1) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, i + 2 * this->degree) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, i + 3 * this->degree - 1) = tmp; + tmp *= -2.0; + + for (unsigned int j = 0; j < this->degree - 1;) { + interpolation_matrix (i + 2, j + (i + 2) * this->degree + 2 - i) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, i + (j + 4) * this->degree - 2) = tmp; + j = j + 2; + } + + tmp *= -1.0; + interpolation_matrix (i + source_fe.degree + 1, i + this->degree + 1) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, i + 3 * this->degree - 1) = tmp; + faculty_k = 1; + + for (int j = 2; j <= (int) this->degree; ++j) { + interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j); + faculty_k *= j; + faculty_l = faculty_k; + faculty_kl = 1; + + for (int k = j + 1; k <= (int) this->degree; ++k) { + faculty_kl *= k - j; + faculty_l *= k; + interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = std::pow (0.5, i + k) * faculty_l / (faculty_k * faculty_kl); + } + } + + faculty_i *= i; + faculty_j = faculty_i; + faculty_ij = 1; + + for (int j = i + 1; j <= (int) this->degree; ++j) { + faculty_ij *= j - i; + faculty_j *= j; + tmp = std::pow (0.5, j + 1) * faculty_j / (faculty_i * faculty_ij); + interpolation_matrix (i + 2, j + 2) = tmp; + interpolation_matrix (i + 2, j + this->degree + 1) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, j + 3 * this->degree - 1) = tmp; + tmp *= 2.0; + interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = tmp; + interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp; + faculty_k = 1; + + for (int k = 2; k <= (int) this->degree; ++k) { + interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k); + faculty_k *= k; + faculty_l = faculty_k; + faculty_kl = 1; + + for (int l = k + 1; l <= (int) this->degree; ++l) { + faculty_kl *= l - k; + faculty_l *= l; + interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = tmp * std::pow (0.5, l) * faculty_l / (faculty_k * faculty_kl); + } + } + + tmp *= -1.0; + + for (unsigned int k = 0; k < this->degree - 1;) { + interpolation_matrix (i + 2, k + (j + 2) * this->degree + 2 - j) = tmp; + interpolation_matrix (i + 2 * source_fe.degree, j + (k + 4) * this->degree - 2) = tmp; + k = k + 2; + } + } + } } } }