From: Markus Buerg Date: Tue, 22 Feb 2011 13:48:15 +0000 (+0000) Subject: Function to interpolate Dirichlet boundary values for Raviart-Thomas elements. X-Git-Tag: v8.0.0~4288 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e43172e2030d87024743c3964743069ea4f39a4f;p=dealii.git Function to interpolate Dirichlet boundary values for Raviart-Thomas elements. git-svn-id: https://svn.dealii.org/trunk@23414 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/numerics/vectors.h b/deal.II/include/deal.II/numerics/vectors.h index 3287568555..a148cf6f7d 100644 --- a/deal.II/include/deal.II/numerics/vectors.h +++ b/deal.II/include/deal.II/numerics/vectors.h @@ -1090,6 +1090,92 @@ class VectorTools ConstraintMatrix& constraints, const hp::MappingCollection& mapping_collection = hp::StaticMappingQ1::mapping_collection); + + /** + * Compute constraints that correspond to + * boundary conditions of the form + * $\vec{n}^T\vec{u}=\vec{n}^T\vec{f}$, + * i.e. the normal components of $u$ + * and $f$ shall coincide. + * + * If the ConstraintMatrix @p constraints + * contained values or other + * constraints before, the new ones are + * added or the old ones overwritten, + * if a node of the boundary part to be + * used was already in the list of + * constraints. This is handled by + * using inhomogeneous constraints. Please + * note that when combining adaptive meshes + * and this kind of constraints, the + * Dirichlet conditions should be set + * first, and then completed by hanging + * node constraints, in order to make sure + * that the discretization remains + * consistent. See the discussion on + * conflicting constraints in the + * module on @ref constraints . + * + * This function is explecitly written to + * use with the FE_RaviartThomas elements. + * Thus it throws an exception, if it is + * called with other finite elements. + * + * The second argument of this function + * denotes the first vector component in + * the finite element that corresponds to + * the vector function that you want to + * constrain. Vectors are implicitly + * assumed to have exactly + * dim components that are + * ordered in the same way as we + * usually order the coordinate directions, + * i.e. $x$-, $y$-, and finally + * $z$-component. + * + * The parameter @p boundary_component + * corresponds to the number + * @p boundary_indicator of the face. 255 + * is an illegal value, since it is + * reserved for interior faces. + * + * The last argument is denoted to compute + * the normal vector $\vec{n}$ at the + * boundary points. + * + *

Computing constraints

+ * + * To compute the constraints we use + * interpolation operator proposed + * in Brezzi, Fortin (Mixed and Hybrid + * (Finite Element Methods, Springer, + * 1991) on every face located at the + * boundary. + * + * @ingroup constraints + */ + template + static void project_boundary_values_div_conforming (const DoFHandler& dof_handler, + const unsigned int first_vector_component, + const Function& boundary_function, + const unsigned char boundary_component, + ConstraintMatrix& constraints, + const Mapping& mapping = StaticMappingQ1::mapping); + + /** + * Same as above for the hp-namespace. + * + * @ingroup constraints + */ + template + static void project_boundary_values_div_conforming (const hp::DoFHandler& dof_handler, + const unsigned int first_vector_component, + const Function& boundary_function, + const unsigned char boundary_component, + ConstraintMatrix& constraints, + const hp::MappingCollection& mapping_collection = hp::StaticMappingQ1::mapping_collection); + + /** * Compute the constraints that * correspond to boundary conditions of diff --git a/deal.II/include/deal.II/numerics/vectors.templates.h b/deal.II/include/deal.II/numerics/vectors.templates.h index 6074cc5147..5ba6cd452f 100644 --- a/deal.II/include/deal.II/numerics/vectors.templates.h +++ b/deal.II/include/deal.II/numerics/vectors.templates.h @@ -37,6 +37,7 @@ #include #include #include +#include #include #include #include @@ -2650,12 +2651,12 @@ namespace internals { // faces. template void - compute_face_projection (const cell_iterator& cell, - const unsigned int face, - hp::FEValues& hp_fe_values, - const Function& boundary_function, - const unsigned int first_vector_component, - std::vector& dof_values) + compute_face_projection_curl_conforming (const cell_iterator& cell, + const unsigned int face, + hp::FEValues& hp_fe_values, + const Function& boundary_function, + const unsigned int first_vector_component, + std::vector& dof_values) { hp_fe_values.reinit (cell, cell->active_fe_index () * GeometryInfo::faces_per_cell + face); @@ -3183,9 +3184,10 @@ project_boundary_values_curl_conforming (const DoFHandler& dof_handler, // boundary function on // the edge. internals::VectorTools - ::compute_face_projection (cell, face, fe_face_values, - boundary_function, - first_vector_component, dof_values); + ::compute_face_projection_curl_conforming (cell, face, fe_face_values, + boundary_function, + first_vector_component, + dof_values); cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); @@ -3318,10 +3320,10 @@ project_boundary_values_curl_conforming (const DoFHandler& dof_handler, // on the interior of // the face. internals::VectorTools - ::compute_face_projection (cell, face, fe_face_values, - boundary_function, - first_vector_component, - dof_values); + ::compute_face_projection_curl_conforming (cell, face, fe_face_values, + boundary_function, + first_vector_component, + dof_values); // Mark the projected // degrees of @@ -3349,6 +3351,8 @@ project_boundary_values_curl_conforming (const DoFHandler& dof_handler, constraints.add_line (dof); constraints.set_inhomogeneity (dof, computed_constraints[dof]); } + + break; } default: @@ -3416,10 +3420,10 @@ project_boundary_values_curl_conforming (const hp::DoFHandler& dof_handler, dof_values[dof] = 0.0; internals::VectorTools - ::compute_face_projection (cell, face, fe_face_values, - boundary_function, - first_vector_component, - dof_values); + ::compute_face_projection_curl_conforming (cell, face, fe_face_values, + boundary_function, + first_vector_component, + dof_values); face_dof_indices.resize (dofs_per_face); cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); @@ -3521,10 +3525,10 @@ project_boundary_values_curl_conforming (const hp::DoFHandler& dof_handler, if (degree > 0) { internals::VectorTools - ::compute_face_projection (cell, face, fe_face_values, - boundary_function, - first_vector_component, - dof_values); + ::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) @@ -3552,6 +3556,470 @@ project_boundary_values_curl_conforming (const hp::DoFHandler& dof_handler, } +namespace internals { + namespace VectorTools { + // This function computes the + // projection of the boundary + // function on the boundary + // in 2d. + template + void + compute_face_projection_div_conforming (const cell_iterator& cell, + const unsigned int face, + const FEFaceValues<2>& fe_values, + const unsigned int first_vector_component, + const Function<2>& boundary_function, + const std::vector >& jacobians, + ConstraintMatrix& constraints) + { + // Compute the intergral over + // the product of the normal + // components of the boundary + // function times the normal + // components of the shape + // functions supported on the + // boundary. + const FEValuesExtractors::Vector vec (first_vector_component); + const FiniteElement<2>& fe = cell->get_fe (); + const std::vector >& normals = fe_values.get_normal_vectors (); + const unsigned int + face_coordinate_direction[GeometryInfo<2>::faces_per_cell] = {1, 1, 0, 0}; + std::vector > + values (fe_values.n_quadrature_points, Vector (2)); + Vector dof_values (fe.dofs_per_face); + + // Get the values of the + // boundary function at the + // quadrature points. + { + const std::vector >& + quadrature_points = fe_values.get_quadrature_points (); + + boundary_function.vector_value_list (quadrature_points, values); + } + + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) + { + double tmp = 0.0; + + for (unsigned int d = 0; d < 2; ++d) + tmp += normals[q_point][d] * values[q_point] (d); + + tmp *= fe_values.JxW (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) + dof_values (i) += tmp * (normals[q_point] + * fe_values[vec].value (fe.face_to_cell_index (i, face), q_point)); + } + + std::vector face_dof_indices (fe.dofs_per_face); + + cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); + + // Copy the computed values + // in the ConstraintMatrix only, + // if the degree of freedom is + // not already constrained. + for (unsigned int i = 0; i < fe.dofs_per_face; ++i) + if (!(constraints.is_constrained (face_dof_indices[i]))) + { + constraints.add_line (face_dof_indices[i]); + + if (std::abs (dof_values (i)) > 1e-14) + constraints.set_inhomogeneity (face_dof_indices[i], dof_values (i)); + } + } + + // dummy implementation of above + // function for all other + // dimensions + template + void + compute_face_projection_div_conforming (const cell_iterator&, + const unsigned int, + const FEFaceValues&, + const unsigned int, + const Function&, + const std::vector >&, + ConstraintMatrix&) + { + Assert (false, ExcNotImplemented ()); + } + + // This function computes the + // projection of the boundary + // function on the boundary + // in 3d. + template + void + compute_face_projection_div_conforming (const cell_iterator& cell, + const unsigned int face, + const FEFaceValues<3>& fe_values, + const unsigned int first_vector_component, + const Function<3>& boundary_function, + const std::vector >& jacobians, + std::vector& dof_values, + std::vector& projected_dofs) + { + // Compute the intergral over + // the product of the normal + // components of the boundary + // function times the normal + // components of the shape + // functions supported on the + // boundary. + const FEValuesExtractors::Vector vec (first_vector_component); + const FiniteElement<3>& fe = cell->get_fe (); + const std::vector >& normals = fe_values.get_normal_vectors (); + const unsigned int + face_coordinate_directions[GeometryInfo<3>::faces_per_cell][2] = {{1, 2}, + {1, 2}, + {2, 0}, + {2, 0}, + {0, 1}, + {0, 1}}; + std::vector > + values (fe_values.n_quadrature_points, Vector (3)); + Vector dof_values_local (fe.dofs_per_face); + + { + const std::vector >& + quadrature_points = fe_values.get_quadrature_points (); + + boundary_function.vector_value_list (quadrature_points, values); + } + + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) + { + double tmp = 0.0; + + for (unsigned int d = 0; d < 3; ++d) + tmp += normals[q_point][d] * values[q_point] (d); + + tmp *= fe_values.JxW (q_point) + * std::sqrt ((jacobians[q_point][0][face_coordinate_directions[face][0]] + * jacobians[q_point][0][face_coordinate_directions[face][0]] + + jacobians[q_point][1][face_coordinate_directions[face][0]] + * jacobians[q_point][1][face_coordinate_directions[face][0]] + + jacobians[q_point][2][face_coordinate_directions[face][0]] + * jacobians[q_point][2][face_coordinate_directions[face][0]]) + * (jacobians[q_point][0][face_coordinate_directions[face][1]] + * jacobians[q_point][0][face_coordinate_directions[face][1]] + + jacobians[q_point][1][face_coordinate_directions[face][1]] + * jacobians[q_point][1][face_coordinate_directions[face][1]] + + jacobians[q_point][2][face_coordinate_directions[face][1]] + * jacobians[q_point][2][face_coordinate_directions[face][1]])); + + for (unsigned int i = 0; i < fe.dofs_per_face; ++i) + dof_values_local (i) += tmp * (normals[q_point] + * fe_values[vec].value (fe.face_to_cell_index (i, face), q_point)); + } + + std::vector face_dof_indices (fe.dofs_per_face); + + cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); + + for (unsigned int i = 0; i < fe.dofs_per_face; ++i) + if (projected_dofs[face_dof_indices[i]] < fe.degree) + { + dof_values[face_dof_indices[i]] = dof_values_local (i); + projected_dofs[face_dof_indices[i]] = fe.degree; + } + } + + // dummy implementation of above + // function for all other + // dimensions + template + void + compute_face_projection_div_conforming (const cell_iterator&, + const unsigned int, + const FEFaceValues&, + const unsigned int, + const Function&, + const std::vector >&, + std::vector&, + std::vector&) + { + Assert (false, ExcNotImplemented ()); + } + } +} + + +template +void +VectorTools::project_boundary_values_div_conforming (const DoFHandler& dof_handler, + const unsigned int first_vector_component, + const Function& boundary_function, + const unsigned char boundary_component, + ConstraintMatrix& constraints, + const Mapping& mapping) +{ + // Interpolate the normal components + // of the boundary functions. Since + // the Raviart-Thomas elements are + // constructed from a Lagrangian + // basis, it suffices to compute + // the integral over the product + // of the normal components of the + // boundary function times the + // normal components of the shape + // functions supported on the + // boundary. + const FiniteElement& fe = dof_handler.get_fe (); + QGauss face_quadrature (fe.degree + 1); + FEFaceValues fe_face_values (mapping, fe, face_quadrature, update_JxW_values | + update_normal_vectors | + update_quadrature_points | + update_values); + hp::FECollection fe_collection (fe); + hp::MappingCollection mapping_collection (mapping); + hp::QCollection quadrature_collection; + + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + quadrature_collection.push_back (QProjector::project_to_face (face_quadrature, + face)); + + hp::FEValues fe_values (mapping_collection, fe_collection, quadrature_collection, + update_jacobians); + + switch (dim) + { + case 2: + { + for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); + 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) + { + // this is only + // implemented, if the + // FE is a Raviart-Thomas + // element + typedef FiniteElement FEL; + + AssertThrow (dynamic_cast*> (&cell->get_fe ()) != 0, + typename FEL::ExcInterpolationNotImplemented ()); + fe_values.reinit (cell, face + cell->active_fe_index () + * GeometryInfo::faces_per_cell); + + const std::vector >& + jacobians = fe_values.get_present_fe_values ().get_jacobians (); + + fe_face_values.reinit (cell, face); + internals::VectorTools::compute_face_projection_div_conforming (cell, face, + fe_face_values, + first_vector_component, + boundary_function, + jacobians, + constraints); + } + + break; + } + + case 3: + { + // In three dimensions the + // edges between two faces + // are treated twice. + // Therefore we store the + // computed values in a + // vector and copy them over + // in the ConstraintMatrix + // after all values have been + // computed. + // If we have two values for + // one edge, we choose the one, + // which was computed with the + // higher order element. + // If both elements are of the + // same order, we just keep the + // first value and do not + // compute a second one. + const unsigned int& n_dofs = dof_handler.n_dofs (); + std::vector dof_values (n_dofs); + std::vector projected_dofs (n_dofs); + + for (unsigned int dof = 0; dof < n_dofs; ++dof) + projected_dofs[dof] = 0; + + for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); + 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) + { + // this is only + // implemented, if the + // FE is a Raviart-Thomas + // element + typedef FiniteElement FEL; + + AssertThrow (dynamic_cast*> (&cell->get_fe ()) != 0, + typename FEL::ExcInterpolationNotImplemented ()); + fe_values.reinit (cell, face + cell->active_fe_index () + * GeometryInfo::faces_per_cell); + + const std::vector >& + jacobians = fe_values.get_present_fe_values ().get_jacobians (); + + fe_face_values.reinit (cell, face); + internals::VectorTools::compute_face_projection_div_conforming (cell, face, + fe_face_values, + first_vector_component, + boundary_function, + jacobians, dof_values, + projected_dofs); + } + + for (unsigned int dof = 0; dof < n_dofs; ++dof) + if ((projected_dofs[dof] != 0) && !(constraints.is_constrained (dof))) + { + constraints.add_line (dof); + + if (std::abs (dof_values[dof]) > 1e-14) + constraints.set_inhomogeneity (dof, dof_values[dof]); + } + + break; + } + + default: + Assert (false, ExcNotImplemented ()); + } +} + + +template +void +VectorTools::project_boundary_values_div_conforming (const hp::DoFHandler& dof_handler, + const unsigned int first_vector_component, + const Function& boundary_function, + const unsigned char boundary_component, + ConstraintMatrix& constraints, + const hp::MappingCollection& mapping_collection) +{ + const hp::FECollection& fe_collection = dof_handler.get_fe (); + hp::QCollection face_quadrature_collection; + hp::QCollection quadrature_collection; + + for (unsigned int i = 0; i < fe_collection.size (); ++i) + { + const QGauss quadrature (fe_collection[i].degree + 1); + + face_quadrature_collection.push_back (quadrature); + + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + quadrature_collection.push_back (QProjector::project_to_face (quadrature, + face)); + } + + hp::FEFaceValues fe_face_values (mapping_collection, fe_collection, + face_quadrature_collection, update_JxW_values | + update_normal_vectors | + update_quadrature_points | + update_values); + hp::FEValues fe_values (mapping_collection, fe_collection, quadrature_collection, + update_jacobians); + + switch (dim) + { + case 2: + { + for (typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); + 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) + { + // this is only + // implemented, if the + // FE is a Raviart-Thomas + // element + typedef FiniteElement FEL; + + AssertThrow (dynamic_cast*> (&cell->get_fe ()) != 0, + typename FEL::ExcInterpolationNotImplemented ()); + fe_values.reinit (cell, face + cell->active_fe_index () + * GeometryInfo::faces_per_cell); + + const std::vector >& + jacobians = fe_values.get_present_fe_values ().get_jacobians (); + + fe_face_values.reinit (cell, face); + internals::VectorTools::compute_face_projection_div_conforming (cell, face, + fe_face_values.get_present_fe_values (), + first_vector_component, + boundary_function, + jacobians, + constraints); + } + + break; + } + + case 3: + { + const unsigned int& n_dofs = dof_handler.n_dofs (); + std::vector dof_values (n_dofs); + std::vector projected_dofs (n_dofs); + + for (unsigned int dof = 0; dof < n_dofs; ++dof) + projected_dofs[dof] = 0; + + for (typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); + 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) + { + // this is only + // implemented, if the + // FE is a Raviart-Thomas + // element + typedef FiniteElement FEL; + + AssertThrow (dynamic_cast*> (&cell->get_fe ()) != 0, + typename FEL::ExcInterpolationNotImplemented ()); + fe_values.reinit (cell, face + cell->active_fe_index () + * GeometryInfo::faces_per_cell); + + const std::vector >& + jacobians = fe_values.get_present_fe_values ().get_jacobians (); + + fe_face_values.reinit (cell, face); + internals::VectorTools::compute_face_projection_div_conforming (cell, face, + fe_face_values.get_present_fe_values (), + first_vector_component, + boundary_function, + jacobians, dof_values, + projected_dofs); + } + + for (unsigned int dof = 0; dof < n_dofs; ++dof) + if ((projected_dofs[dof] != 0) && !(constraints.is_constrained (dof))) + { + constraints.add_line (dof); + + if (std::abs (dof_values[dof]) > 1e-14) + constraints.set_inhomogeneity (dof, dof_values[dof]); + } + + break; + } + + default: + Assert (false, ExcNotImplemented ()); + } +} + template class DH, int spacedim> void diff --git a/deal.II/source/numerics/vectors.inst.in b/deal.II/source/numerics/vectors.inst.in index 6bf0f7c76d..5026d7464c 100644 --- a/deal.II/source/numerics/vectors.inst.in +++ b/deal.II/source/numerics/vectors.inst.in @@ -652,6 +652,22 @@ void VectorTools::project_boundary_values_curl_conforming ConstraintMatrix&, const hp::MappingCollection&); template +void VectorTools::project_boundary_values_div_conforming +(const DoFHandler&, + const unsigned int, + const Function&, + const unsigned char, + ConstraintMatrix&, + const Mapping&); +template +void VectorTools::project_boundary_values_div_conforming +(const hp::DoFHandler&, + const unsigned int, + const Function&, + const unsigned char, + ConstraintMatrix&, + const hp::MappingCollection&); +template void VectorTools::compute_no_normal_flux_constraints (const DoFHandler &dof_handler, const unsigned int first_vector_component,