From dbe54da52a8372cbfac56f035f004c93c098a956 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 27 Apr 1998 18:03:45 +0000 Subject: [PATCH] Still changes for the support of integration along faces. git-svn-id: https://svn.dealii.org/trunk@205 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe.h | 34 ++++++-- deal.II/deal.II/include/fe/fe_lib.lagrange.h | 37 +++++++++ deal.II/deal.II/include/fe/fe_update_flags.h | 10 ++- deal.II/deal.II/include/fe/fe_values.h | 44 +++++++++- deal.II/deal.II/source/fe/fe.cc | 16 +++- deal.II/deal.II/source/fe/fe_lib.linear.cc | 84 ++++++++++++++++++++ deal.II/deal.II/source/fe/fe_values.cc | 13 +++ 7 files changed, 224 insertions(+), 14 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 188a130e5c..142b1b7cd4 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -546,15 +546,16 @@ class FiniteElement : public FiniteElementBase { * performed by the specific finite element * class. * - * Two fields remain to be finite element + * Three fields remain to be finite element * specific in this standard implementation: * The jacobi determinants of the * transformation from the unit face to the - * real face and the ansatz points. For - * these two fields, there exist pure - * virtual functions, #get_face_jacobians# - * and the #get_face_ansatz_points# - * function. + * real face, the ansatz points + * and the outward normal vectors. For + * these fields, there exist pure + * virtual functions, #get_face_jacobians#, + * #get_face_ansatz_points# and + * #get_normal_vectors#. * * Though there is a standard * implementation, there @@ -587,6 +588,8 @@ class FiniteElement : public FiniteElementBase { 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; /** @@ -644,6 +647,25 @@ class FiniteElement : public FiniteElementBase { const Boundary &boundary, const vector > &unit_points, vector &face_jacobi_determinants) const =0; + + /** + * Compute the normal vectors to the cell + * at the quadrature points. See the + * documentation for the #fill_fe_face_values# + * function for more details. The function + * must guarantee that the length of the + * vectors be one. + * + * Since any implementation for one + * dimension would be senseless, all + * derived classes should throw an error + * when called with #dim==1#. + */ + virtual void get_normal_vectors (const DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const Boundary &boundary, + const vector > &unit_points, + vector > &normal_vectors) const =0; /** * Exception diff --git a/deal.II/deal.II/include/fe/fe_lib.lagrange.h b/deal.II/deal.II/include/fe/fe_lib.lagrange.h index ab435d964e..b340857559 100644 --- a/deal.II/deal.II/include/fe/fe_lib.lagrange.h +++ b/deal.II/deal.II/include/fe/fe_lib.lagrange.h @@ -83,6 +83,23 @@ class FELinear : public FiniteElement { const Boundary &boundary, const vector > &unit_points, vector &face_jacobi_determinants) const; + + /** + * For linear finite elements, this function + * is particularly simple since all normal + * vectors are equal and can easiliy be + * computed from the direction of the face + * without using the transformation (Jacobi) + * matrix, at least for two dimensions. + * + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_normal_vectors (const DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const Boundary &boundary, + const vector > &unit_points, + vector > &normal_vectors) const; }; @@ -149,6 +166,16 @@ class FEQuadratic : public FiniteElement { const Boundary &boundary, const vector > &unit_points, vector &face_jacobi_determinants) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_normal_vectors (const DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const Boundary &boundary, + const vector > &unit_points, + vector > &normal_vectors) const; }; @@ -215,6 +242,16 @@ class FECubic : public FiniteElement { const Boundary &boundary, const vector > &unit_points, vector &face_jacobi_determinants) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_normal_vectors (const DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const Boundary &boundary, + const vector > &unit_points, + vector > &normal_vectors) const; }; diff --git a/deal.II/deal.II/include/fe/fe_update_flags.h b/deal.II/deal.II/include/fe/fe_update_flags.h index c6562c858e..1a464f2878 100644 --- a/deal.II/deal.II/include/fe/fe_update_flags.h +++ b/deal.II/deal.II/include/fe/fe_update_flags.h @@ -47,12 +47,18 @@ enum UpdateFlags { * on which the ansatz functions are * located. */ - update_ansatz_points = 16 + update_ansatz_points = 16, + /** + * Update the outward normal vectors + * to the face relative to this cell. + * This flag is only evaluated by + * the #FEFaceValues# class. + */ + update_normal_vectors = 32 }; - /*---------------------------- fe_update_flags.h ---------------------------*/ /* end of #ifndef __fe_update_flags_H */ #endif diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index e4946a8299..0715c95ad8 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -206,7 +206,7 @@ class FEValues { * segments, but higher order elements * may use other ways.) */ - void reinit (const DoFHandler::cell_iterator &, + void reinit (const typename DoFHandler::cell_iterator &, const FiniteElement &, const Boundary &); @@ -386,6 +386,14 @@ class FEValues { for the transformation of the unit face to the real face (needed to compute the weight factors for integration along faces). These two concepts have to be carefully separated. + + Finally, we will often need the outward normal to a cell at the quadrature + points. While this could in principle be easily done using the Jacobi + matrices at the quadrature points and the normal vectors to the unit cell + (also easily derived, since they have an appealingly easy form for the unit + cell ;-), it is more efficiently done by the finite element class itself. + For example for (bi-, tri-)linear mappings the normal vector is readily + available without compicated matrix-vector-multiplications. */ template class FEFaceValues { @@ -533,6 +541,21 @@ class FEFaceValues { * alike cells. */ const vector & get_JxW_values () const; + + /** + * Return the outward normal vector to + * the cell at the #i#th quadrature + * point. The length of the vector + * is normalized to one. + */ + const Point & normal_vector (const unsigned int i) const; + + /** + * Return the list of outward normal + * vectors to the cell at the + * quadrature points. + */ + const vector > & get_normal_vectors () const; /** * Reinitialize the gradients, Jacobi @@ -553,7 +576,7 @@ class FEFaceValues { * segments, but higher order elements * may use other ways.) */ - void reinit (const DoFHandler::cell_iterator &cell, + void reinit (const typename DoFHandler::cell_iterator &cell, const unsigned int face_no, const FiniteElement &fe, const Boundary &boundary); @@ -707,6 +730,13 @@ class FEFaceValues { * compute the JxW values. */ vector face_jacobi_determinants; + + /** + * List of outward normal vectors at the + * quadrature points. This field is filled + * in by the finite element class. + */ + vector > normal_vectors; /** * Store which fields are to be updated by @@ -833,6 +863,16 @@ FEFaceValues::get_JxW_values () const { +template +inline +const vector > & +FEFaceValues::get_normal_vectors () const { + Assert (update_flags & update_normal_vectors, ExcAccessToUninitializedField()); + return normal_vectors; +}; + + + /*---------------------------- fe_values.h ---------------------------*/ /* end of #ifndef __fe_values_H */ diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index e08fda1ac9..fea46fe848 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -196,6 +196,8 @@ void FiniteElement<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator & const bool , vector &, const bool , + vector > &, + const bool, const Boundary<1> &) const { Assert (false, ExcNotImplemented()); } @@ -215,6 +217,8 @@ void FiniteElement::fill_fe_face_values (const DoFHandler::cell_iterat 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())); @@ -232,17 +236,21 @@ void FiniteElement::fill_fe_face_values (const DoFHandler::cell_iterat q_points, compute_q_points, boundary); - cout << "Global unit points:\n"; - for (unsigned int p=0; pface(face_no), boundary, ansatz_points); if (compute_face_jacobians) get_face_jacobians (cell->face(face_no), boundary, unit_points, face_jacobi_determinants); + + if (compute_normal_vectors) + get_normal_vectors (cell, face_no, boundary, + unit_points, normal_vectors); + cout << "Global unit points:\n"; + for (unsigned int p=0; p::get_face_jacobians (const DoFHandler<1>::face_iterator &, +void FELinear<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const Boundary<1> &, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + + + + FELinear<2>::FELinear () : FiniteElement<2> (1, 0, 0) { @@ -357,6 +369,35 @@ void FELinear<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face, +void FELinear<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell, + const unsigned int face_no, + const Boundary<2> &, + const vector > &unit_points, + vector > &normal_vectors) const { + Assert (unit_points.size() == normal_vectors.size(), + ExcWrongFieldDimension (unit_points.size(), normal_vectors.size())); + + 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()))); +}; + + + @@ -400,6 +441,18 @@ void FEQuadratic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, +void FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const Boundary<1> &, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + + + + FEQuadratic<2>::FEQuadratic () : FiniteElement<2> (1, 1, 1) { @@ -483,6 +536,17 @@ void FEQuadratic::get_face_jacobians (const DoFHandler::face_iterator +template +void FEQuadratic::get_normal_vectors (const DoFHandler::cell_iterator &, + const unsigned int, + const Boundary &, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + + @@ -526,6 +590,16 @@ void FECubic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, +void FECubic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const Boundary<1> &, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + + FECubic<2>::FECubic () : FiniteElement<2> (1, 2, 4) {}; @@ -598,6 +672,16 @@ void FECubic::get_face_jacobians (const DoFHandler::face_iterator &, +template +void FECubic::get_normal_vectors (const DoFHandler::cell_iterator &, + const unsigned int, + const Boundary &, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + // explicit instantiations diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index a332c178a4..3fd49eb326 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -180,6 +180,7 @@ FEFaceValues::FEFaceValues (const FiniteElement &fe, ansatz_points (fe.total_dofs, Point()), jacobi_matrices (quadrature.n_quadrature_points,dFMatrix(dim,dim)), face_jacobi_determinants (quadrature.n_quadrature_points,0), + normal_vectors (quadrature.n_quadrature_points,Point()), update_flags (update_flags), selected_face(0) { @@ -298,6 +299,16 @@ const Point & FEFaceValues::ansatz_point (const unsigned int i) const +template +const Point & FEFaceValues::normal_vector (const unsigned int i) const { + Assert (i double FEFaceValues::JxW (const unsigned int i) const { Assert (i::reinit (const typename DoFHandler::cell_iterator &c update_flags & update_q_points, face_jacobi_determinants, update_flags & update_JxW_values, + normal_vectors, + update_flags & update_normal_vectors, boundary); // compute gradients on real element if -- 2.39.5