From 964fa6545418a73bd97e59c6cce08853e04729b4 Mon Sep 17 00:00:00 2001 From: Markus Buerg Date: Thu, 20 Sep 2012 09:20:12 +0000 Subject: [PATCH] Make project_boundary_values_curl_conforming work for arbitrary FESystem. git-svn-id: https://svn.dealii.org/trunk@26554 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/numerics/vector_tools.templates.h | 530 +++++++++--------- 1 file changed, 275 insertions(+), 255 deletions(-) diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h index ee7bb166ed..2b7c2ad0cc 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@ -2769,7 +2769,8 @@ namespace VectorTools hp::FEValues<3>& hp_fe_values, const Function<3>& boundary_function, const unsigned int first_vector_component, - std::vector& dof_values) + std::vector& dof_values, + std::vector& dofs_processed) { const double tol = 0.5 * cell->get_fe ().degree * 1e-13 / cell->face (face)->line (line)->diameter (); const unsigned int dim = 3; @@ -2784,6 +2785,7 @@ namespace VectorTools // objects. const FEValues& fe_values = hp_fe_values.get_present_fe_values (); + const FiniteElement& fe = cell->get_fe (); const std::vector< DerivativeForm<1,dim,spacedim> > & jacobians = fe_values.get_jacobians (); const std::vector >& @@ -2791,7 +2793,7 @@ namespace VectorTools std::vector > tangentials (fe_values.n_quadrature_points); std::vector > values (fe_values.n_quadrature_points, - Vector (cell->get_fe ().n_components ())); + Vector (fe.n_components ())); // Get boundary function values // at quadrature points. @@ -2799,8 +2801,27 @@ namespace VectorTools const std::vector >& reference_quadrature_points = fe_values.get_quadrature ().get_points (); - const unsigned int superdegree = cell->get_fe ().degree; - const unsigned int degree = superdegree - 1; + const unsigned int degree = fe.degree - 1; + std::pair base_indices (0, 0); + + if (dynamic_cast*> (&cell->get_fe ()) != 0) + { + unsigned int fe_index = 0; + unsigned int fe_index_old = 0; + unsigned int i = 0; + + for (; i < fe.n_base_elements (); ++i) + { + fe_index_old = fe_index; + fe_index += fe.element_multiplicity (i) * fe.base_element (i).n_components (); + + if (fe_index >= first_vector_component) + break; + } + + base_indices.first = i; + base_indices.second = (first_vector_component - fe_index_old) / fe.base_element (i).n_components (); + } // coordinate directions of // the edges of the face. @@ -2848,22 +2869,33 @@ namespace VectorTools // Compute the degrees of // freedom. - for (unsigned int i = 0; i <= degree; ++i) - dof_values[i + line * superdegree] - += (fe_values.JxW (q_point) - * (values[q_point] (first_vector_component) * tangentials[q_point] (0) - + values[q_point] (first_vector_component + 1) * tangentials[q_point] (1) - + values[q_point] (first_vector_component + 2) * tangentials[q_point] (2)) - * (fe_values[vec].value (cell->get_fe ().face_to_cell_index (i + line * superdegree, - face), - q_point) - * tangentials[q_point]) - / std::sqrt (jacobians[q_point][0][edge_coordinate_direction[face][line]] - * jacobians[q_point][0][edge_coordinate_direction[face][line]] - + jacobians[q_point][1][edge_coordinate_direction[face][line]] - * jacobians[q_point][1][edge_coordinate_direction[face][line]] - + jacobians[q_point][2][edge_coordinate_direction[face][line]] - * jacobians[q_point][2][edge_coordinate_direction[face][line]])); + for (unsigned int i = 0; i < fe.dofs_per_face; ++i) + if (((dynamic_cast*> (&fe) != 0) + && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices) + && (fe.base_element (base_indices.first).face_to_cell_index (line * fe.degree, face) + <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second) + && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second + <= fe.base_element (base_indices.first).face_to_cell_index + ((line + 1) * fe.degree - 1, face))) + || ((dynamic_cast*> (&fe) != 0) && (line * fe.degree <= i) + && (i < (line + 1) * fe.degree))) + { + dof_values[i] + += (fe_values.JxW (q_point) + * (values[q_point] (first_vector_component) * tangentials[q_point] (0) + + values[q_point] (first_vector_component + 1) * tangentials[q_point] (1) + + values[q_point] (first_vector_component + 2) * tangentials[q_point] (2)) + * (fe_values[vec].value (fe.face_to_cell_index (i, face), q_point) * tangentials[q_point]) + / std::sqrt (jacobians[q_point][0][edge_coordinate_direction[face][line]] + * jacobians[q_point][0][edge_coordinate_direction[face][line]] + + jacobians[q_point][1][edge_coordinate_direction[face][line]] + * jacobians[q_point][1][edge_coordinate_direction[face][line]] + + jacobians[q_point][2][edge_coordinate_direction[face][line]] + * jacobians[q_point][2][edge_coordinate_direction[face][line]])); + + if (q_point == 0) + dofs_processed[i] = true; + } } } @@ -2878,7 +2910,8 @@ namespace VectorTools hp::FEValues&, const Function&, const unsigned int, - std::vector&) + std::vector&, + std::vector&) { Assert (false, ExcInternalError ()); } @@ -2894,7 +2927,8 @@ namespace VectorTools hp::FEValues& hp_fe_values, const Function& boundary_function, const unsigned int first_vector_component, - std::vector& dof_values) + std::vector& dof_values, + std::vector& dofs_processed) { const unsigned int spacedim = dim; hp_fe_values.reinit (cell, cell->active_fe_index () @@ -2903,29 +2937,50 @@ namespace VectorTools // objects. const FEValues& fe_values = hp_fe_values.get_present_fe_values (); + const FiniteElement& fe = cell->get_fe (); const std::vector< DerivativeForm<1,dim,spacedim> > & jacobians = fe_values.get_jacobians (); + const std::vector >& + quadrature_points = fe_values.get_quadrature_points (); + const unsigned int degree = fe.degree - 1; + std::pair base_indices (0, 0); + + if (dynamic_cast*> (&cell->get_fe ()) != 0) + { + unsigned int fe_index = 0; + unsigned int fe_index_old = 0; + unsigned int i = 0; + + for (; i < fe.n_base_elements (); ++i) + { + fe_index_old = fe_index; + fe_index += fe.element_multiplicity (i) * fe.base_element (i).n_components (); + + if (fe_index >= first_vector_component) + break; + } + + base_indices.first = i; + base_indices.second = (first_vector_component - fe_index_old) / fe.base_element (i).n_components (); + } std::vector > - values (fe_values.n_quadrature_points, Vector (cell->get_fe ().n_components ())); + values (fe_values.n_quadrature_points, Vector (fe.n_components ())); + + // Get boundary function + // values at quadrature + // points. + boundary_function.vector_value_list (quadrature_points, values); switch (dim) { case 2: { - const std::vector >& - quadrature_points = fe_values.get_quadrature_points (); std::vector > tangentials (fe_values.n_quadrature_points); - // Get boundary function - // values at quadrature - // points. - boundary_function.vector_value_list (quadrature_points, values); - const std::vector >& reference_quadrature_points = fe_values.get_quadrature ().get_points (); - const unsigned int degree = cell->get_fe ().degree - 1; // coordinate directions // of the face. @@ -2970,19 +3025,27 @@ namespace VectorTools /= std::sqrt (tangentials[q_point].square ()); // Compute the degrees // of freedom. - for (unsigned int i = 0; i <= degree; ++i) - dof_values[i] - += fe_values.JxW (q_point) - * (values[q_point] (first_vector_component) - * tangentials[q_point] (0) - + values[q_point] (first_vector_component + 1) - * tangentials[q_point] (1)) - * (fe_values[vec].value (cell->get_fe ().face_to_cell_index (i, face), - q_point) * tangentials[q_point]) - / std::sqrt (jacobians[q_point][0][face_coordinate_direction[face]] - * jacobians[q_point][0][face_coordinate_direction[face]] - + jacobians[q_point][1][face_coordinate_direction[face]] - * jacobians[q_point][1][face_coordinate_direction[face]]); + for (unsigned int i = 0; i < fe.dofs_per_face; ++i) + if (((dynamic_cast*> (&fe) != 0) + && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)) + || (dynamic_cast*> (&fe) != 0)) + { + dof_values[i] + += fe_values.JxW (q_point) + * (values[q_point] (first_vector_component) + * tangentials[q_point] (0) + + values[q_point] (first_vector_component + 1) + * tangentials[q_point] (1)) + * (fe_values[vec].value (fe.face_to_cell_index (i, face), q_point) + * tangentials[q_point]) + / std::sqrt (jacobians[q_point][0][face_coordinate_direction[face]] + * jacobians[q_point][0][face_coordinate_direction[face]] + + jacobians[q_point][1][face_coordinate_direction[face]] + * jacobians[q_point][1][face_coordinate_direction[face]]); + + if (q_point == 0) + dofs_processed[i] = true; + } } break; @@ -2990,19 +3053,9 @@ namespace VectorTools case 3: { - const std::vector >& - quadrature_points = fe_values.get_quadrature_points (); - - // Get boundary function - // values at quadrature - // points. - boundary_function.vector_value_list (quadrature_points, values); - const FEValuesExtractors::Vector vec (first_vector_component); - const unsigned int superdegree = cell->get_fe ().degree; - const unsigned int degree = superdegree - 1; FullMatrix - assembling_matrix (degree * superdegree, + assembling_matrix (degree * fe.degree, dim * fe_values.n_quadrature_points); Vector assembling_vector (assembling_matrix.n ()); Vector cell_rhs (assembling_matrix.m ()); @@ -3058,13 +3111,17 @@ namespace VectorTools for (unsigned int d = 0; d < dim; ++d) tmp[d] = values[q_point] (first_vector_component + d); - - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j <= degree; ++j) - tmp -= dof_values[(i + 2) * superdegree + j] - * fe_values[vec].value (cell->get_fe ().face_to_cell_index - ((i + 2) * superdegree + j, - face), q_point); + + for (unsigned int i = 0; i < fe.dofs_per_face; ++i) + if (((dynamic_cast*> (&fe) != 0) + && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices) + && (fe.base_element (base_indices.first).face_to_cell_index (2 * fe.degree, face) + <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second) + && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second + < fe.base_element (base_indices.first).face_to_cell_index (4 * fe.degree, face))) + || ((dynamic_cast*> (&fe) != 0) && (2 * fe.degree <= i) + && (i < 4 * fe.degree))) + tmp -= dof_values[i] * fe_values[vec].value (fe.face_to_cell_index (i, face), q_point); const double JxW = std::sqrt (fe_values.JxW (q_point) @@ -3097,25 +3154,30 @@ namespace VectorTools // the face. for (unsigned int d = 0; d < dim; ++d) assembling_vector (dim * q_point + d) = JxW * tmp[d]; - - for (unsigned int i = 0; i <= degree; ++i) - for (unsigned int j = 0; j < degree; ++j) + + unsigned int index = 0; + + for (unsigned int i = 0; i < fe.dofs_per_face; ++i) + if (((dynamic_cast*> (&fe) != 0) + && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices) + && (fe.base_element (base_indices.first).face_to_cell_index + (GeometryInfo::lines_per_face * fe.degree, face) + <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second) + && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second + < fe.base_element (base_indices.first).face_to_cell_index + ((degree + GeometryInfo::lines_per_face) * fe.degree, face))) + || ((dynamic_cast*> (&fe) != 0) + && (GeometryInfo::lines_per_face * fe.degree <= i) + && (i < (degree + GeometryInfo::lines_per_face) * fe.degree))) { const Tensor<1, dim> shape_value - = (JxW - * fe_values[vec].value (cell->get_fe () - .face_to_cell_index - ((i + GeometryInfo::lines_per_face) - * degree - + j - + GeometryInfo::lines_per_face, - face), - q_point)); + = (JxW * fe_values[vec].value (fe.face_to_cell_index (i, face), + q_point)); for (unsigned int d = 0; d < dim; ++d) - assembling_matrix (i * degree + j, - dim * q_point + d) - = shape_value[d]; + assembling_matrix (index, dim * q_point + d) = shape_value[d]; + + ++index; } } @@ -3139,11 +3201,27 @@ namespace VectorTools // Store the computed // values. - for (unsigned int i = 0; i <= degree; ++i) - for (unsigned int j = 0; j < degree; ++j) - dof_values[(i + GeometryInfo::lines_per_face) - * degree + j + GeometryInfo::lines_per_face] - = solution (i * degree + j); + { + unsigned int index = 0; + + for (unsigned int i = 0; i < fe.dofs_per_face; ++i) + if (((dynamic_cast*> (&fe) != 0) + && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices) + && (fe.base_element (base_indices.first).face_to_cell_index + (GeometryInfo::lines_per_face * fe.degree, face) + <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second) + && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second + < fe.base_element (base_indices.first).face_to_cell_index + ((degree + GeometryInfo::lines_per_face) * fe.degree, face))) + || ((dynamic_cast*> (&fe) != 0) + && (GeometryInfo::lines_per_face * fe.degree <= i) + && (i < (degree + GeometryInfo::lines_per_face) * fe.degree))) + { + dof_values[i] = solution (index); + dofs_processed[i] = true; + ++index; + } + } // Now we do the // same as above @@ -3159,12 +3237,13 @@ namespace VectorTools for (unsigned int d = 0; d < dim; ++d) tmp[d] = values[q_point] (first_vector_component + d); - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j <= degree; ++j) - tmp - -= dof_values[i * superdegree + j] - * fe_values[vec].value (cell->get_fe ().face_to_cell_index - (i * superdegree + j, face), q_point); + for (unsigned int i = 0; i < fe.dofs_per_face; ++i) + if (((dynamic_cast*> (&fe) != 0) + && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices) + && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second + < fe.base_element (base_indices.first).face_to_cell_index (2 * fe.degree, face))) + || ((dynamic_cast*> (&fe) != 0) && (i < 2 * fe.degree))) + tmp -= dof_values[i] * fe_values[vec].value (fe.face_to_cell_index (i, face), q_point); const double JxW = std::sqrt (fe_values.JxW (q_point) @@ -3189,18 +3268,25 @@ namespace VectorTools for (unsigned int d = 0; d < dim; ++d) assembling_vector (dim * q_point + d) = JxW * tmp[d]; - for (unsigned int i = 0; i < degree; ++i) - for (unsigned int j = 0; j <= degree; ++j) + unsigned int index = 0; + + for (unsigned int i = 0; i < fe.dofs_per_face; ++i) + if (((dynamic_cast*> (&fe) != 0) + && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices) + && (fe.base_element (base_indices.first).face_to_cell_index + ((degree + GeometryInfo::lines_per_face) * fe.degree, face) + <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second)) + || ((dynamic_cast*> (&fe) != 0) + && ((degree + GeometryInfo::lines_per_face) * fe.degree <= i))) { const Tensor<1, dim> shape_value - = (JxW - * fe_values[vec].value (cell->get_fe ().face_to_cell_index - ((i + degree + GeometryInfo::lines_per_face) - * superdegree + j, face), q_point)); + = JxW * fe_values[vec].value (fe.face_to_cell_index (i, face), + q_point); for (unsigned int d = 0; d < dim; ++d) - assembling_matrix (i * superdegree + j, dim * q_point + d) - = shape_value[d]; + assembling_matrix (index, dim * q_point + d) = shape_value[d]; + + ++index; } } @@ -3209,10 +3295,21 @@ namespace VectorTools assembling_matrix.vmult (cell_rhs, assembling_vector); cell_matrix_inv.vmult (solution, cell_rhs); - for (unsigned int i = 0; i < degree; ++i) - for (unsigned int j = 0; j <= degree; ++j) - dof_values[(i + degree + GeometryInfo::lines_per_face) * superdegree + j] - = solution (i * superdegree + j); + unsigned int index = 0; + + for (unsigned int i = 0; i < fe.dofs_per_face; ++i) + if (((dynamic_cast*> (&fe) != 0) + && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices) + && (fe.base_element (base_indices.first).face_to_cell_index + ((degree + GeometryInfo::lines_per_face) * fe.degree, face) + <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second)) + || ((dynamic_cast*> (&fe) != 0) + && ((degree + GeometryInfo::lines_per_face) * fe.degree <= i))) + { + dof_values[i] = solution (index); + dofs_processed[i] = true; + ++index; + } break; } @@ -3280,6 +3377,7 @@ namespace VectorTools update_quadrature_points | update_values); + std::vector dofs_processed (dofs_per_face); std::vector dof_values (dofs_per_face); std::vector face_dof_indices (dofs_per_face); typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); @@ -3306,7 +3404,8 @@ namespace VectorTools // element. If the FE // is a FESystem, we // cannot check this. - if (dynamic_cast*> (&cell->get_fe ()) == 0) { + if (dynamic_cast*> (&cell->get_fe ()) == 0) + { typedef FiniteElement FEL; AssertThrow (dynamic_cast*> (&cell->get_fe ()) != 0, @@ -3314,7 +3413,10 @@ namespace VectorTools } for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + { dof_values[dof] = 0.0; + dofs_processed[dof] = false; + } // Compute the // projection of the @@ -3324,7 +3426,7 @@ namespace VectorTools ::compute_face_projection_curl_conforming (cell, face, fe_face_values, boundary_function, first_vector_component, - dof_values); + dof_values, dofs_processed); cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); @@ -3334,11 +3436,11 @@ namespace VectorTools // if the degree of // freedom is not // already constrained. - const double tol = 1e-13; + const double tol = 0.5 * cell->get_fe ().degree * 1e-13 / cell->face (face)->diameter (); for (unsigned int dof = 0; dof < dofs_per_face; ++dof) - if (constraints.can_store_line (face_dof_indices[dof]) - && !(constraints.is_constrained (face_dof_indices[dof]))) + if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof]) + && !(constraints.is_constrained (face_dof_indices[dof]))) { constraints.add_line (face_dof_indices[dof]); @@ -3354,7 +3456,6 @@ namespace VectorTools { const QGauss reference_edge_quadrature (2 * superdegree); const unsigned int degree = superdegree - 1; - const unsigned int n_dofs = dof_handler.n_dofs (); hp::QCollection edge_quadrature_collection; for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) @@ -3370,13 +3471,6 @@ namespace VectorTools update_JxW_values | update_quadrature_points | update_values); - std::vector computed_constraints (n_dofs); - std::vector projected_dofs (n_dofs); - - for (unsigned int dof = 0; dof < n_dofs; ++dof) { - computed_constraints[dof] = 0.0; - projected_dofs[dof] = -1; - } for (; cell != dof_handler.end (); ++cell) if (cell->at_boundary () && cell->is_locally_owned ()) @@ -3405,62 +3499,21 @@ namespace VectorTools } for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + { dof_values[dof] = 0.0; - - cell->face (face)->get_dof_indices (face_dof_indices, - cell->active_fe_index ()); + dofs_processed[dof] = false; + } // First we compute the // projection on the // edges. for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) - { - // If we have reached - // this edge through - // another cell - // before, we do not - // do here anything - // unless we have a - // good reason, i.e. - // a higher - // polynomial degree. - if (projected_dofs[face_dof_indices[line * superdegree]] - < - (int) degree) - { - // Compute the - // projection of - // the boundary - // function on the - // edge. - internals - ::compute_edge_projection (cell, face, line, - fe_edge_values, - boundary_function, - first_vector_component, - dof_values); - // Mark the - // projected - // degrees of - // freedom. - for (unsigned int dof = line * superdegree; - dof < (line + 1) * superdegree; ++dof) - projected_dofs[face_dof_indices[dof]] = degree; - } - - // If we have - // computed the - // values in a - // previous step of - // the loop, we just - // copy the values in - // the local vector. - else - for (unsigned int dof = line * superdegree; - dof < (line + 1) * superdegree; ++dof) - dof_values[dof] = computed_constraints[face_dof_indices[dof]]; - } + internals + ::compute_edge_projection (cell, face, line, fe_edge_values, + boundary_function, + first_vector_component, + dof_values, dofs_processed); // If there are higher // order shape @@ -3468,47 +3521,32 @@ namespace VectorTools // still some work // left. if (degree > 0) - { - // Compute the - // projection of the - // boundary function - // on the interior of - // the face. - internals - ::compute_face_projection_curl_conforming (cell, face, fe_face_values, - boundary_function, - first_vector_component, - dof_values); - - // Mark the projected - // degrees of - // freedom. - for (unsigned int dof = GeometryInfo::lines_per_face * superdegree; - dof < dofs_per_face; ++dof) - projected_dofs[face_dof_indices[dof]] = degree; - } + internals + ::compute_face_projection_curl_conforming (cell, face, fe_face_values, + boundary_function, + first_vector_component, + dof_values, + dofs_processed); // Store the computed // values in the global // vector. + + cell->face (face)->get_dof_indices (face_dof_indices, + cell->active_fe_index ()); + const double tol = 0.5 * superdegree * 1e-13 / cell->face (face)->diameter (); for (unsigned int dof = 0; dof < dofs_per_face; ++dof) - if (std::abs (computed_constraints[face_dof_indices[dof]] - dof_values[dof]) > tol) - computed_constraints[face_dof_indices[dof]] = dof_values[dof]; - } + if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof]) + && !(constraints.is_constrained (face_dof_indices[dof]))) + { + constraints.add_line (face_dof_indices[dof]); - // Add the computed constraints - // to the constraint matrix, if - // the degree of freedom is not - // already constrained. - for (unsigned int dof = 0; dof < n_dofs; ++dof) - if ((projected_dofs[dof] != -1) && constraints.can_store_line (dof) - && !(constraints.is_constrained (dof))) - { - constraints.add_line (dof); - constraints.set_inhomogeneity (dof, computed_constraints[dof]); - } + if (std::abs (dof_values[dof]) > tol) + constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]); + } + } break; } @@ -3549,6 +3587,7 @@ namespace VectorTools update_JxW_values | update_quadrature_points | update_values); + std::vector dofs_processed; std::vector dof_values; std::vector face_dof_indices; typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); @@ -3584,17 +3623,21 @@ namespace VectorTools } const unsigned int dofs_per_face = cell->get_fe ().dofs_per_face; - + + dofs_processed.resize (dofs_per_face); dof_values.resize (dofs_per_face); for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + { dof_values[dof] = 0.0; + dofs_processed[dof] = false; + } internals ::compute_face_projection_curl_conforming (cell, face, fe_face_values, boundary_function, first_vector_component, - dof_values); + dof_values, dofs_processed); face_dof_indices.resize (dofs_per_face); cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); @@ -3602,8 +3645,8 @@ namespace VectorTools const double tol = 0.5 * cell->get_fe ().degree * 1e-13 / cell->face (face)->diameter (); for (unsigned int dof = 0; dof < dofs_per_face; ++dof) - if (constraints.can_store_line (face_dof_indices[dof]) - && !(constraints.is_constrained (face_dof_indices[dof]))) + if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof]) + && !(constraints.is_constrained (face_dof_indices[dof]))) { constraints.add_line (face_dof_indices[dof]); @@ -3617,7 +3660,6 @@ namespace VectorTools case 3: { - const unsigned int n_dofs = dof_handler.n_dofs (); hp::QCollection edge_quadrature_collection; for (unsigned int i = 0; i < fe_collection.size (); ++i) @@ -3639,13 +3681,6 @@ namespace VectorTools update_JxW_values | update_quadrature_points | update_values); - std::vector computed_constraints (n_dofs); - std::vector projected_dofs (n_dofs); - - for (unsigned int dof = 0; dof < n_dofs; ++dof) { - computed_constraints[dof] = 0.0; - projected_dofs[dof] = -1; - } for (; cell != dof_handler.end (); ++cell) if (cell->at_boundary () && cell->is_locally_owned ()) @@ -3676,68 +3711,53 @@ namespace VectorTools const unsigned int superdegree = cell->get_fe ().degree; const unsigned int degree = superdegree - 1; const unsigned int dofs_per_face = cell->get_fe ().dofs_per_face; - + + dofs_processed.resize (dofs_per_face); dof_values.resize (dofs_per_face); for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + { dof_values[dof] = 0.0; - - face_dof_indices.resize (dofs_per_face); - cell->face (face)->get_dof_indices (face_dof_indices, - cell->active_fe_index ()); + dofs_processed[dof] = false; + } for (unsigned int line = 0; line < GeometryInfo::lines_per_face; ++line) - { - if (projected_dofs[face_dof_indices[line * superdegree]] - < - (int) degree) - { - internals - ::compute_edge_projection (cell, face, line, - fe_edge_values, - boundary_function, - first_vector_component, - dof_values); - - for (unsigned int dof = line * superdegree; - dof < (line + 1) * superdegree; ++dof) - projected_dofs[face_dof_indices[dof]] = degree; - } - - else - for (unsigned int dof = line * superdegree; - dof < (line + 1) * superdegree; ++dof) - dof_values[dof] = computed_constraints[face_dof_indices[dof]]; - } + internals + ::compute_edge_projection (cell, face, line, fe_edge_values, + boundary_function, + first_vector_component, + dof_values, dofs_processed); + // If there are higher + // order shape + // functions, there is + // still some work + // left. if (degree > 0) - { - internals - ::compute_face_projection_curl_conforming (cell, face, fe_face_values, - boundary_function, - first_vector_component, - dof_values); - - for (unsigned int dof = GeometryInfo::lines_per_face * superdegree; - dof < dofs_per_face; ++dof) - projected_dofs[face_dof_indices[dof]] = degree; - } + internals + ::compute_face_projection_curl_conforming (cell, face, fe_face_values, + boundary_function, + first_vector_component, + dof_values, dofs_processed); + + face_dof_indices.resize (dofs_per_face); + cell->face (face)->get_dof_indices (face_dof_indices, + cell->active_fe_index ()); + const double tol = 0.5 * superdegree * 1e-13 / cell->face (face)->diameter (); for (unsigned int dof = 0; dof < dofs_per_face; ++dof) - if (std::abs (computed_constraints[face_dof_indices[dof]] - dof_values[dof]) > tol) - computed_constraints[face_dof_indices[dof]] = dof_values[dof]; - } + if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof]) + && !(constraints.is_constrained (face_dof_indices[dof]))) + { + constraints.add_line (face_dof_indices[dof]); - for (unsigned int dof = 0; dof < n_dofs; ++dof) - if ((projected_dofs[dof] != -1) && constraints.can_store_line (dof) - && !(constraints.is_constrained (dof))) - { - constraints.add_line (dof); - constraints.set_inhomogeneity (dof, computed_constraints[dof]); - } + if (std::abs (dof_values[dof]) > tol) + constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]); + } + } break; } -- 2.39.5