From c8929dc724c1601d4dfa28766420974c7bd2356a Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 5 May 1998 08:54:09 +0000 Subject: [PATCH] Add some new functions for the restriction of finite elements to subfaces. Fix a small problem. git-svn-id: https://svn.dealii.org/trunk@248 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/fe.cc | 61 +++++++- deal.II/deal.II/source/fe/fe_lib.linear.cc | 173 ++++++++++++++++++++- deal.II/deal.II/source/fe/fe_values.cc | 42 ++--- 3 files changed, 249 insertions(+), 27 deletions(-) diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 04cb987add..602e91e364 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -200,7 +200,7 @@ void FiniteElement<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator & const bool, const Boundary<1> &) const { Assert (false, ExcNotImplemented()); -} +}; @@ -250,6 +250,65 @@ void FiniteElement::fill_fe_face_values (const DoFHandler::cell_iterat + +void FiniteElement<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterator &, + const unsigned int , + const unsigned int , + const vector > &, + const vector > &, + vector &, + const bool , + vector > &, + const bool , + vector &, + const bool , + vector > &, + const bool, + const Boundary<1> &) const { + Assert (false, ExcNotImplemented()); +}; + + + +template +void FiniteElement::fill_fe_subface_values (const DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const vector > &unit_points, + const vector > &global_unit_points, + vector &jacobians, + const bool compute_jacobians, + vector > &q_points, + const bool compute_q_points, + vector &face_jacobi_determinants, + const bool compute_face_jacobians, + vector > &normal_vectors, + const bool compute_normal_vectors, + const Boundary &boundary) const { + Assert (jacobians.size() == unit_points.size(), + ExcWrongFieldDimension(jacobians.size(), unit_points.size())); + Assert (q_points.size() == unit_points.size(), + ExcWrongFieldDimension(q_points.size(), unit_points.size())); + Assert (global_unit_points.size() == unit_points.size(), + ExcWrongFieldDimension(global_unit_points.size(), unit_points.size())); + + static vector > dummy(total_dofs); + fill_fe_values (cell, global_unit_points, + jacobians, compute_jacobians, + dummy, false, + q_points, compute_q_points, + boundary); + + if (compute_face_jacobians) + get_subface_jacobians (cell->face(face_no), subface_no, + unit_points, face_jacobi_determinants); + + if (compute_normal_vectors) + get_normal_vectors (cell, face_no, subface_no, + unit_points, normal_vectors); +}; + + /*------------------------------- Explicit Instantiations -------------*/ template class FiniteElementData<1>; diff --git a/deal.II/deal.II/source/fe/fe_lib.linear.cc b/deal.II/deal.II/source/fe/fe_lib.linear.cc index 739df3b6ad..e242af0205 100644 --- a/deal.II/deal.II/source/fe/fe_lib.linear.cc +++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc @@ -106,6 +106,15 @@ void FELinear<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, +void FELinear<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, + const unsigned int , + const vector > &, + vector &) const { + Assert (false, ExcInternalError()); +}; + + + void FELinear<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, const unsigned int, const Boundary<1> &, @@ -116,6 +125,16 @@ void FELinear<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, +void FELinear<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const unsigned int, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + + FELinear<2>::FELinear () : @@ -369,6 +388,27 @@ void FELinear<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face, +void FELinear<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face, + const unsigned int, + const vector > &unit_points, + vector &face_jacobians) const { + Assert (unit_points.size() == face_jacobians.size(), + ExcWrongFieldDimension (unit_points.size(), face_jacobians.size())); + Assert (face->at_boundary() == false, + ExcBoundaryFaceUsed ()); + + // a linear mapping for a single line + // produces particularly simple + // expressions for the jacobi + // determinant :-) + const double h = sqrt((face->vertex(1) - face->vertex(0)).square()); + fill_n (face_jacobians.begin(), + unit_points.size(), + h/2); +}; + + + void FELinear<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell, const unsigned int face_no, const Boundary<2> &, @@ -398,6 +438,40 @@ void FELinear<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell, +void FELinear<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell, + const unsigned int face_no, + const unsigned int, + const vector > &unit_points, + vector > &normal_vectors) const { + // note, that in 2D the normal vectors to the + // subface have the same direction as that + // for the face + Assert (unit_points.size() == normal_vectors.size(), + ExcWrongFieldDimension (unit_points.size(), normal_vectors.size())); + Assert (cell->face(face_no)->at_boundary() == false, + ExcBoundaryFaceUsed ()); + + const DoFHandler<2>::face_iterator face = cell->face(face_no); + // compute direction of line + const Point<2> line_direction = (face->vertex(1) - face->vertex(0)); + // rotate to the right by 90 degrees + const Point<2> normal_direction(line_direction(1), + -line_direction(0)); + + if (face_no <= 1) + // for sides 0 and 1: return the correctly + // scaled vector + fill (normal_vectors.begin(), normal_vectors.end(), + normal_direction / sqrt(normal_direction.square())); + else + // for sides 2 and 3: scale and invert + // vector + fill (normal_vectors.begin(), normal_vectors.end(), + normal_direction / (-sqrt(normal_direction.square()))); +}; + + + @@ -441,6 +515,15 @@ void FEQuadratic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, +void FEQuadratic<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, + const unsigned int , + const vector > &, + vector &) const { + Assert (false, ExcInternalError()); +}; + + + void FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, const unsigned int, const Boundary<1> &, @@ -451,6 +534,15 @@ void FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, +void FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const unsigned int, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + FEQuadratic<2>::FEQuadratic () : @@ -528,9 +620,22 @@ void FEQuadratic::get_face_ansatz_points (const typename DoFHandler::f template void FEQuadratic::get_face_jacobians (const DoFHandler::face_iterator &, - const Boundary &, - const vector > &, - vector &) const { + const Boundary &, + const vector > &, + vector &) const { + Assert (false, ExcNotImplemented()); +}; + + + +template +void FEQuadratic::get_subface_jacobians (const DoFHandler::face_iterator &face, + const unsigned int , + const vector > &, + vector &) const { + Assert (face->at_boundary() == false, + ExcBoundaryFaceUsed ()); + Assert (false, ExcNotImplemented()); }; @@ -542,7 +647,20 @@ void FEQuadratic::get_normal_vectors (const DoFHandler::cell_iterator const Boundary &, const vector > &, vector > &) const { - Assert (false, ExcInternalError()); + Assert (false, ExcNotImplemented()); +}; + + + +template +void FEQuadratic::get_normal_vectors (const DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int, + const vector > &, + vector > &) const { + Assert (cell->face(face_no)->at_boundary() == false, + ExcBoundaryFaceUsed ()); + Assert (false, ExcNotImplemented()); }; @@ -590,6 +708,15 @@ void FECubic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, +void FECubic<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, + const unsigned int , + const vector > &, + vector &) const { + Assert (false, ExcInternalError()); +}; + + + void FECubic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, const unsigned int, const Boundary<1> &, @@ -600,6 +727,15 @@ void FECubic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, +void FECubic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const unsigned int, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + FECubic<2>::FECubic () : FiniteElement<2> (1, 2, 4) {}; @@ -672,13 +808,40 @@ void FECubic::get_face_jacobians (const DoFHandler::face_iterator &, +template +void FECubic::get_subface_jacobians (const DoFHandler::face_iterator &face, + const unsigned int , + const vector > &, + vector &) const { + Assert (face->at_boundary() == false, + ExcBoundaryFaceUsed ()); + + Assert (false, ExcNotImplemented()); +}; + + + template void FECubic::get_normal_vectors (const DoFHandler::cell_iterator &, const unsigned int, const Boundary &, const vector > &, vector > &) const { - Assert (false, ExcInternalError()); + Assert (false, ExcNotImplemented()); +}; + + + +template +void FECubic::get_normal_vectors (const DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int , + const vector > &, + vector > &) const { + Assert (cell->face(face_no)->at_boundary() == false, + ExcBoundaryFaceUsed ()); + + Assert (false, ExcNotImplemented()); }; diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 462690cf5c..eca2ef337d 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -160,7 +160,7 @@ FEValues::FEValues (const FiniteElement &fe, vector >(quadrature.n_quadrature_points)), unit_quadrature_points(quadrature.get_quad_points()) { - Assert ((update_flags | update_normal_vectors) == false, + Assert ((update_flags & update_normal_vectors) == false, ExcInvalidUpdateFlag()); for (unsigned int i=0; i::FEFaceValuesBase (const unsigned int n_q_points, n_dofs, n_faces_or_subfaces, update_flags), + unit_shape_gradients (n_faces_or_subfaces, + vector > >(n_dofs, + vector >(n_q_points))), unit_face_quadrature_points (n_q_points, Point()), unit_quadrature_points (n_faces_or_subfaces, vector >(n_q_points, Point())), face_jacobi_determinants (n_q_points, 0), normal_vectors (n_q_points) -{ - for (unsigned int i=0; i >(n_q_points)); - unit_quadrature_points[i].resize (n_q_points, - Point()); - }; -}; +{}; @@ -358,10 +353,11 @@ void FEFaceValues::reinit (const typename DoFHandler::cell_iterator &c Assert (update_flags & update_jacobians, ExcCannotInitializeField()); for (unsigned int i=0; i(); - + { + fill_n (shape_gradients[i].begin(), + n_quadrature_points, + Point()); + for (unsigned int j=0; j::reinit (const typename DoFHandler::cell_iterator &c shape_gradients[i][j](s) += (unit_shape_gradients[face_no][i][j](b) * jacobi_matrices[j](b,s)); - }; + }; }; @@ -443,6 +439,9 @@ void FESubfaceValues::reinit (const typename DoFHandler::cell_iterator const unsigned int subface_no, const FiniteElement &fe, const Boundary &boundary) { + Assert (cell->face(face_no)->at_boundary() == false, + ExcReinitCalledWithBoundaryFace()); + present_cell = cell; selected_dataset = face_no*(1<<(dim-1)) + subface_no; // fill jacobi matrices and real @@ -471,11 +470,12 @@ void FESubfaceValues::reinit (const typename DoFHandler::cell_iterator { Assert (update_flags & update_jacobians, ExcCannotInitializeField()); - for (unsigned int i=0; i(); - + for (unsigned int i=0; i()); + for (unsigned int j=0; j::reinit (const typename DoFHandler::cell_iterator shape_gradients[i][j](s) += (unit_shape_gradients[selected_dataset][i][j](b) * jacobi_matrices[j](b,s)); - }; + }; }; -- 2.39.5