From 66d61dd68f3d5e4da84686c0ead37b1a4fee50bf Mon Sep 17 00:00:00 2001 From: buerg Date: Fri, 3 Jun 2011 13:59:58 +0000 Subject: [PATCH] Extended VectorTools::compute_no_normal_flux_constraints to hp::DoFHandler. git-svn-id: https://svn.dealii.org/trunk@23787 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/numerics/vectors.templates.h | 88 +++++++++++++++---- deal.II/source/numerics/vectors.inst.in | 7 ++ 2 files changed, 76 insertions(+), 19 deletions(-) diff --git a/deal.II/include/deal.II/numerics/vectors.templates.h b/deal.II/include/deal.II/numerics/vectors.templates.h index 7f3fe5fc5b..c37ddb6b90 100644 --- a/deal.II/include/deal.II/numerics/vectors.templates.h +++ b/deal.II/include/deal.II/numerics/vectors.templates.h @@ -41,7 +41,7 @@ #include #include #include -#include +#include #include #include #include @@ -2447,6 +2447,7 @@ namespace internals { const unsigned int first_vector_component, std::vector& dof_values) { + const double tol = 0.5 * cell->get_fe ().degree * 1e-13 / cell->face (face)->line (line)->diameter (); const unsigned int dim = 3; hp_fe_values.reinit @@ -2504,17 +2505,18 @@ namespace internals { Point shifted_reference_point_1 = reference_quadrature_points[q_point]; Point shifted_reference_point_2 = reference_quadrature_points[q_point]; - shifted_reference_point_1 (edge_coordinate_direction[face][line]) += 1e-13; - shifted_reference_point_2 (edge_coordinate_direction[face][line]) -= 1e-13; + shifted_reference_point_1 (edge_coordinate_direction[face][line]) += tol; + shifted_reference_point_2 (edge_coordinate_direction[face][line]) -= tol; tangentials[q_point] - = (2e13 * + = (0.5 * (fe_values.get_mapping () .transform_unit_to_real_cell (cell, shifted_reference_point_1) - fe_values.get_mapping () .transform_unit_to_real_cell (cell, - shifted_reference_point_2))); + shifted_reference_point_2)) + / tol); tangentials[q_point] /= std::sqrt (tangentials[q_point].square ()); @@ -2675,6 +2677,7 @@ namespace internals { { case 2: { + const double tol = 0.5 * cell->get_fe ().degree * 1e-13 / cell->face (face)->diameter (); const std::vector >& quadrature_points = fe_values.get_quadrature_points (); std::vector > @@ -2715,18 +2718,19 @@ namespace internals { = reference_quadrature_points[q_point]; shifted_reference_point_1 (face_coordinate_direction[face]) - += 1e-13; + += tol; shifted_reference_point_2 (face_coordinate_direction[face]) - -= 1e-13; + -= tol; tangentials[q_point] - = 2e13 + = 0.5 * (fe_values.get_mapping () .transform_unit_to_real_cell (cell, shifted_reference_point_1) - fe_values.get_mapping () .transform_unit_to_real_cell (cell, - shifted_reference_point_2)); + shifted_reference_point_2) + / tol); tangentials[q_point] /= std::sqrt (tangentials[q_point].square ()); // Compute the mean @@ -3169,6 +3173,13 @@ project_boundary_values_curl_conforming (const DoFHandler& dof_handler, for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) if (cell->face (face)->boundary_indicator () == boundary_component) { + // if the FE is a + // FE_Nothing object + // there is no work to + // do + if (dynamic_cast*> (&cell->get_fe ()) != 0) + return; + // this is only // implemented, if the // FE is a Nedelec @@ -3198,12 +3209,14 @@ project_boundary_values_curl_conforming (const DoFHandler& dof_handler, // if the degree of // freedom is not // already constrained. + const double tol = 0.5 * superdegree * 1e-13 / cell->face (face)->diameter (); + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) if (!(constraints.is_constrained (face_dof_indices[dof]))) { constraints.add_line (face_dof_indices[dof]); - if (std::abs (dof_values[dof]) > 1e-14) + if (std::abs (dof_values[dof]) > tol) constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]); } } @@ -3234,14 +3247,23 @@ project_boundary_values_curl_conforming (const DoFHandler& dof_handler, std::vector computed_constraints (n_dofs); std::vector projected_dofs (n_dofs); - for (unsigned int dof = 0; dof < n_dofs; ++dof) + 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 ()) for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) if (cell->face (face)->boundary_indicator () == boundary_component) { + // if the FE is a + // FE_Nothing object + // there is no work to + // do + if (dynamic_cast*> (&cell->get_fe ()) != 0) + return; + // this is only // implemented, if the // FE is a Nedelec @@ -3337,8 +3359,10 @@ project_boundary_values_curl_conforming (const DoFHandler& dof_handler, // Store the computed // values in the global // vector. + 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 (dof_values[dof]) > 1e-14) + if (std::abs (computed_constraints[face_dof_indices[dof]] - dof_values[dof]) > tol) computed_constraints[face_dof_indices[dof]] = dof_values[dof]; } @@ -3405,6 +3429,13 @@ project_boundary_values_curl_conforming (const hp::DoFHandler& dof_handler, for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) if (cell->face (face)->boundary_indicator () == boundary_component) { + // if the FE is a + // FE_Nothing object + // there is no work to + // do + if (dynamic_cast*> (&cell->get_fe ()) != 0) + return; + // This is only // implemented, if the // FE is a Nedelec @@ -3435,12 +3466,14 @@ project_boundary_values_curl_conforming (const hp::DoFHandler& dof_handler, cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); + 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.is_constrained (face_dof_indices[dof]))) { constraints.add_line (face_dof_indices[dof]); - if (std::abs (dof_values[dof]) > 1e-14) + if (std::abs (dof_values[dof]) > tol) constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]); } } @@ -3475,14 +3508,23 @@ project_boundary_values_curl_conforming (const hp::DoFHandler& dof_handler, std::vector computed_constraints (n_dofs); std::vector projected_dofs (n_dofs); - for (unsigned int dof = 0; dof < n_dofs; ++dof) + 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 ()) for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) if (cell->face (face)->boundary_indicator () == boundary_component) { + // if the FE is a + // FE_Nothing object + // there is no work to + // do + if (dynamic_cast*> (&cell->get_fe ()) != 0) + return; + // This is only // implemented, if the // FE is a Nedelec @@ -3548,8 +3590,10 @@ project_boundary_values_curl_conforming (const hp::DoFHandler& dof_handler, projected_dofs[face_dof_indices[dof]] = degree; } + 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 (dof_values[dof]) > 1e-14) + if (std::abs (computed_constraints[face_dof_indices[dof]] - dof_values[dof]) > tol) computed_constraints[face_dof_indices[dof]] = dof_values[dof]; } @@ -3811,6 +3855,13 @@ VectorTools::project_boundary_values_div_conforming (const DoFHandler& dof_ for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) if (cell->face (face)->boundary_indicator () == boundary_component) { + // if the FE is a + // FE_Nothing object + // there is no work to + // do + if (dynamic_cast*> (&cell->get_fe ()) != 0) + return; + // This is only // implemented, if the // FE is a Raviart-Thomas @@ -4072,10 +4123,7 @@ compute_no_normal_flux_constraints (const DH &dof_handler, "to imposing Dirichlet values on the vector-valued " "quantity.")); - const FiniteElement &fe = dof_handler.get_fe(); - - std::vector face_dofs (fe.dofs_per_face); - std::vector > dof_locations (fe.dofs_per_face); + std::vector face_dofs; // have a map that stores normal // vectors for each vector-dof @@ -4120,10 +4168,12 @@ compute_no_normal_flux_constraints (const DH &dof_handler, if (boundary_ids.find(cell->face(face_no)->boundary_indicator()) != boundary_ids.end()) { + const FiniteElement& fe = cell->get_fe (); typename DH::face_iterator face = cell->face(face_no); // get the indices of the // dofs on this cell... + face_dofs.resize (fe.dofs_per_face); face->get_dof_indices (face_dofs, cell->active_fe_index()); // ...and the normal diff --git a/deal.II/source/numerics/vectors.inst.in b/deal.II/source/numerics/vectors.inst.in index 2d3b5e8b95..472f8c1ee2 100644 --- a/deal.II/source/numerics/vectors.inst.in +++ b/deal.II/source/numerics/vectors.inst.in @@ -690,6 +690,13 @@ VectorTools::compute_no_normal_flux_constraints (const DoFHandler &boundary_ids, ConstraintMatrix &constraints, const Mapping &mapping); +template +void +VectorTools::compute_no_normal_flux_constraints (const hp::DoFHandler &dof_handler, + const unsigned int first_vector_component, + const std::set &boundary_ids, + ConstraintMatrix &constraints, + const Mapping &mapping); #endif -- 2.39.5