From 6b3aa8205ee56466c9ca276a362fc385c54e7436 Mon Sep 17 00:00:00 2001 From: buerg Date: Sun, 22 Jul 2012 21:07:53 +0000 Subject: [PATCH] Restore comment and fix projection on faces. git-svn-id: https://svn.dealii.org/trunk@25718 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/source/fe/fe_nedelec.cc | 183 ++++++++++++++++++++++---------- 1 file changed, 125 insertions(+), 58 deletions(-) diff --git a/deal.II/source/fe/fe_nedelec.cc b/deal.II/source/fe/fe_nedelec.cc index b1da9ccd6e..a2791856c2 100644 --- a/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/source/fe/fe_nedelec.cc @@ -16,6 +16,10 @@ #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 @@ -1053,7 +1057,6 @@ FE_Nedelec::initialize_restriction () } - template std::vector FE_Nedelec::get_dpo_vector (const unsigned int degree, bool dg) @@ -1752,59 +1755,127 @@ FE_Nedelec::get_subface_interpolation_matrix( lobatto_polynomials = Polynomials::Lobatto::generate_complete_basis (source_fe.degree); - const unsigned int n_edge_dofs + const unsigned int n_boundary_dofs = GeometryInfo::lines_per_face * source_fe.degree; const unsigned int& n_quadrature_points = quadrature.size (); + FullMatrix + system_matrix_inv (source_fe.degree * (source_fe.degree - 1), + source_fe.degree * (source_fe.degree - 1)); - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) + { + 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.invert (system_matrix); + } + + FullMatrix solution (system_matrix_inv.m (), 2); + FullMatrix system_rhs (system_matrix_inv.m (), 2); + Vector tmp (2); + + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { - 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) + system_rhs = 0.0; + + for (unsigned int q_point = 0; + q_point < n_quadrature_points; ++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); - + 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); + for (unsigned int i = 0; i < source_fe.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)); - + 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 < source_fe.degree - 1; ++j) { - 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_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)); } } } + + 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); + } } } @@ -1897,10 +1968,9 @@ FE_Nedelec::interpolate (std::vector& local_dofs, // 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) - if (std::abs (local_dofs[i + e * this->degree]) < 1e-14) - local_dofs[i + e * this->degree] = 0.0; + for (unsigned int i = 0; i < GeometryInfo::lines_per_cell * this->degree; ++i) + if (std::abs (local_dofs[i]) < 1e-14) + local_dofs[i] = 0.0; // If the degree is greater // than 0, then we have still @@ -2116,10 +2186,9 @@ FE_Nedelec::interpolate (std::vector& local_dofs, // 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) - if (std::abs (local_dofs[i + e * this->degree]) < 1e-14) - local_dofs[i + e * this->degree] = 0.0; + for (unsigned int i = 0; i < GeometryInfo::lines_per_cell * this->degree; ++i) + if (std::abs (local_dofs[i]) < 1e-14) + local_dofs[i] = 0.0; // If the degree is greater // than 0, then we have still @@ -3272,10 +3341,9 @@ const // 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) - if (std::abs (local_dofs[i + e * this->degree]) < 1e-14) - local_dofs[i + e * this->degree] = 0.0; + for (unsigned int i = 0; i < GeometryInfo::lines_per_cell * this->degree; ++i) + if (std::abs (local_dofs[i]) < 1e-14) + local_dofs[i] = 0.0; // If the degree is greater // than 0, then we have still @@ -3477,10 +3545,9 @@ const // 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) - if (std::abs (local_dofs[i + e * this->degree]) < 1e-14) - local_dofs[i + e * this->degree] = 0.0; + for (unsigned int i = 0; i < GeometryInfo::lines_per_cell * this->degree; ++i) + if (std::abs (local_dofs[i]) < 1e-14) + local_dofs[i] = 0.0; // If the degree is greater // than 0, then we have still -- 2.39.5