From: Wolfgang Bangerth Date: Mon, 13 Jul 1998 14:01:58 +0000 (+0000) Subject: Separate the linear mapping from the elements implemented as of yet. X-Git-Tag: v8.0.0~22811 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d2ed99bdfc5dd1f55ecf7d5e15a95cddba46ea02;p=dealii.git Separate the linear mapping from the elements implemented as of yet. git-svn-id: https://svn.dealii.org/trunk@443 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index e04a12f35d..cfe3502f39 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -966,6 +966,141 @@ class FiniteElement : public FiniteElementBase { +/** + * Abstract base class for concrete finite elements which use a + * (bi-,tri-)linear mapping from the unit cell to the real cell. Some + * functions can be singled out from these elements and are collected + * in this one. + */ +template +class FELinearMapping : public FiniteElement { + public: + /** + * Constructor. Simply passes through + * its arguments to the base class. + */ + FELinearMapping (const unsigned int dofs_per_vertex, + const unsigned int dofs_per_line, + const unsigned int dofs_per_quad=0) : + FiniteElement (dofs_per_vertex, + dofs_per_line, + dofs_per_quad) {}; + + /** + * Refer to the base class for detailed + * information on this function. + * + * In two spatial dimensions, this function + * simply returns the length of the face. + */ + virtual void get_face_jacobians (const DoFHandler::face_iterator &face, + const Boundary &boundary, + const vector > &unit_points, + vector &face_jacobi_determinants) const; + + /** + * Refer to the base class for detailed + * information on this function. + * + * In two spatial dimensions, this function + * simply returns half the length of the + * whole face. + */ + virtual void get_subface_jacobians (const DoFHandler::face_iterator &face, + const unsigned int subface_no, + const vector > &unit_points, + vector &face_jacobi_determinants) const; + + /** + * Return the normal vectors to the + * face with number #face_no# of #cell#. + * + * 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; + + /** + * Return the normal vectors to the + * subface with number #subface_no# of + * the face with number #face_no# of #cell#. + * + * 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 unsigned int subface_no, + const vector > &unit_points, + vector > &normal_vectors) const; + + /** + * Refer to the base class for detailed + * information on this function. + * + * For one dimensional elements, this + * function simply passes through to + * the one implemented in the base class. + * For higher dimensional finite elements + * we use linear mappings and therefore + * the boundary object is ignored since + * the boundary is approximated using + * piecewise multilinear boundary segments. + */ + virtual void fill_fe_values (const DoFHandler::cell_iterator &cell, + const vector > &unit_points, + vector &jacobians, + const bool compute_jacobians, + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points, + const Boundary &boundary) const; + + protected: + /** + * Return the value of the #i#th shape + * function at point #p# on the unit cell. + * Here, the (bi-)linear basis functions + * are meant, which are used for the + * computation of the transformation from + * unit cell to real space cell. + */ + double linear_shape_value(const unsigned int i, + const Point& p) const; + + /** + * Return the gradient of the #i#th shape + * function at point #p# on the unit cell. + * Here, the (bi-)linear basis functions + * are meant, which are used for the + * computation of the transformation from + * unit cell to real space cell. + */ + Point linear_shape_grad(const unsigned int i, + const Point& p) const; +}; + + + /*---------------------------- fe.h ---------------------------*/ /* end of #ifndef __fe_H */ 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 d11198b085..45777de5c8 100644 --- a/deal.II/deal.II/include/fe/fe_lib.lagrange.h +++ b/deal.II/deal.II/include/fe/fe_lib.lagrange.h @@ -25,7 +25,7 @@ * @author Wolfgang Bangerth, 1998 */ template -class FELinear : public FiniteElement { +class FELinear : public FELinearMapping { public: /** * Constructor @@ -46,29 +46,6 @@ class FELinear : public FiniteElement { virtual Point shape_grad(const unsigned int i, const Point& p) const; - /** - * Refer to the base class for detailed - * information on this function. - * - * For one dimensional elements, this - * function simply passes through to - * the one implemented in the base class. - * For higher dimensional finite elements - * we use linear mappings and therefore - * the boundary object is ignored since - * the boundary is approximated using - * piecewise straight boundary segments. - */ - virtual void fill_fe_values (const DoFHandler::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - const Boundary &boundary) const; - /** * Refer to the base class for detailed * information on this function. @@ -85,72 +62,6 @@ class FELinear : public FiniteElement { const Boundary &boundary, vector > &ansatz_points) const; - /** - * Refer to the base class for detailed - * information on this function. - * - * In two spatial dimensions, this function - * simply returns the length of the face. - */ - virtual void get_face_jacobians (const DoFHandler::face_iterator &face, - const Boundary &boundary, - const vector > &unit_points, - vector &face_jacobi_determinants) const; - - /** - * Refer to the base class for detailed - * information on this function. - * - * In two spatial dimensions, this function - * simply returns half the length of the - * whole face. - */ - virtual void get_subface_jacobians (const DoFHandler::face_iterator &face, - const unsigned int subface_no, - const vector > &unit_points, - vector &face_jacobi_determinants) const; - - /** - * Return the normal vectors to the - * face with number #face_no# of #cell#. - * - * 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; - - /** - * Return the normal vectors to the - * subface with number #subface_no# of - * the face with number #face_no# of #cell#. - * - * 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 unsigned int subface_no, - const vector > &unit_points, - vector > &normal_vectors) const; - /** * Refer to the base class for detailed * information on this function. @@ -169,7 +80,7 @@ class FELinear : public FiniteElement { * to the real cell is implemented. */ template -class FEQuadraticSub : public FiniteElement { +class FEQuadraticSub : public FELinearMapping { public: /** * Constructor @@ -190,24 +101,6 @@ class FEQuadraticSub : public FiniteElement { virtual Point shape_grad(const unsigned int i, const Point& p) const; - /** - * Refer to the base class for detailed - * information on this function. - * - * For one dimensional elements, this - * function simply passes through to - * the one implemented in the base class. - */ - virtual void fill_fe_values (const DoFHandler::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - const Boundary &boundary) const; - /** * Refer to the base class for detailed * information on this function. @@ -224,44 +117,6 @@ class FEQuadraticSub : public FiniteElement { const Boundary &boundary, vector > &ansatz_points) const; - /** - * Refer to the base class for detailed - * information on this function. - */ - virtual void get_face_jacobians (const DoFHandler::face_iterator &face, - 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_subface_jacobians (const DoFHandler::face_iterator &face, - const unsigned int subface_no, - 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; - - /** - * Refer to the base class for detailed - * information on this function. - */ - virtual void get_normal_vectors (const DoFHandler::cell_iterator &cell, - const unsigned int subface_no, - const unsigned int face_no, - const vector > &unit_points, - vector > &normal_vectors) const; - /** * Refer to the base class for detailed * information on this function. @@ -269,29 +124,6 @@ class FEQuadraticSub : public FiniteElement { virtual void get_local_mass_matrix (const DoFHandler::cell_iterator &cell, const Boundary &boundary, dFMatrix &local_mass_matrix) const; - - private: - /** - * Return the value of the #i#th shape - * function at point #p# on the unit cell. - * Here, the (bi-)linear basis functions - * are meant, which are used for the - * computation of the transformation from - * unit cell to real space cell. - */ - double linear_shape_value(const unsigned int i, - const Point& p) const; - - /** - * Return the gradient of the #i#th shape - * function at point #p# on the unit cell. - * Here, the (bi-)linear basis functions - * are meant, which are used for the - * computation of the transformation from - * unit cell to real space cell. - */ - Point linear_shape_grad(const unsigned int i, - const Point& p) const; }; @@ -322,7 +154,7 @@ class FEQuadraticSub : public FiniteElement { * freedom. */ template -class FECubicSub : public FiniteElement { +class FECubicSub : public FELinearMapping { public: /** * Constructor @@ -343,24 +175,6 @@ class FECubicSub : public FiniteElement { virtual Point shape_grad(const unsigned int i, const Point& p) const; - /** - * Refer to the base class for detailed - * information on this function. - * - * For one dimensional elements, this - * function simply passes through to - * the one implemented in the base class. - */ - virtual void fill_fe_values (const DoFHandler::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - const Boundary &boundary) const; - /** * Refer to the base class for detailed * information on this function. @@ -377,44 +191,6 @@ class FECubicSub : public FiniteElement { const Boundary &boundary, vector > &ansatz_points) const; - /** - * Refer to the base class for detailed - * information on this function. - */ - virtual void get_face_jacobians (const DoFHandler::face_iterator &face, - 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_subface_jacobians (const DoFHandler::face_iterator &face, - const unsigned int subface_no, - 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; - - /** - * Refer to the base class for detailed - * information on this function. - */ - virtual void get_normal_vectors (const DoFHandler::cell_iterator &cell, - const unsigned int subface_no, - const unsigned int face_no, - const vector > &unit_points, - vector > &normal_vectors) const; - /** * Refer to the base class for detailed * information on this function. @@ -422,29 +198,6 @@ class FECubicSub : public FiniteElement { virtual void get_local_mass_matrix (const DoFHandler::cell_iterator &cell, const Boundary &boundary, dFMatrix &local_mass_matrix) const; - - private: - /** - * Return the value of the #i#th shape - * function at point #p# on the unit cell. - * Here, the (bi-)linear basis functions - * are meant, which are used for the - * computation of the transformation from - * unit cell to real space cell. - */ - double linear_shape_value(const unsigned int i, - const Point& p) const; - - /** - * Return the gradient of the #i#th shape - * function at point #p# on the unit cell. - * Here, the (bi-)linear basis functions - * are meant, which are used for the - * computation of the transformation from - * unit cell to real space cell. - */ - Point linear_shape_grad(const unsigned int i, - const Point& p) const; }; @@ -476,7 +229,7 @@ class FECubicSub : public FiniteElement { * freedom. */ template -class FEQuarticSub : public FiniteElement { +class FEQuarticSub : public FELinearMapping { public: /** * Constructor @@ -497,24 +250,6 @@ class FEQuarticSub : public FiniteElement { virtual Point shape_grad(const unsigned int i, const Point& p) const; - /** - * Refer to the base class for detailed - * information on this function. - * - * For one dimensional elements, this - * function simply passes through to - * the one implemented in the base class. - */ - virtual void fill_fe_values (const DoFHandler::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - const Boundary &boundary) const; - /** * Refer to the base class for detailed * information on this function. @@ -531,44 +266,6 @@ class FEQuarticSub : public FiniteElement { const Boundary &boundary, vector > &ansatz_points) const; - /** - * Refer to the base class for detailed - * information on this function. - */ - virtual void get_face_jacobians (const DoFHandler::face_iterator &face, - 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_subface_jacobians (const DoFHandler::face_iterator &face, - const unsigned int subface_no, - 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; - - /** - * Refer to the base class for detailed - * information on this function. - */ - virtual void get_normal_vectors (const DoFHandler::cell_iterator &cell, - const unsigned int subface_no, - const unsigned int face_no, - const vector > &unit_points, - vector > &normal_vectors) const; - /** * Refer to the base class for detailed * information on this function. @@ -576,29 +273,6 @@ class FEQuarticSub : public FiniteElement { virtual void get_local_mass_matrix (const DoFHandler::cell_iterator &cell, const Boundary &boundary, dFMatrix &local_mass_matrix) const; - - private: - /** - * Return the value of the #i#th shape - * function at point #p# on the unit cell. - * Here, the (bi-)linear basis functions - * are meant, which are used for the - * computation of the transformation from - * unit cell to real space cell. - */ - double linear_shape_value(const unsigned int i, - const Point& p) const; - - /** - * Return the gradient of the #i#th shape - * function at point #p# on the unit cell. - * Here, the (bi-)linear basis functions - * are meant, which are used for the - * computation of the transformation from - * unit cell to real space cell. - */ - Point linear_shape_grad(const unsigned int i, - const Point& p) const; }; diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index e64d0b4710..b8e20795d4 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -362,10 +362,363 @@ void FiniteElement::get_ansatz_points (const DoFHandler::cell_iterator +#if deal_II_dimension == 1 + +template <> +inline +double +FELinearMapping<1>::linear_shape_value(const unsigned int i, + const Point<1> &p) const +{ + Assert((i<2), ExcInvalidIndex(i)); + const double xi = p(0); + switch (i) + { + case 0: return 1.-xi; + case 1: return xi; + } + return 0.; +}; + + + +template <> +inline +Point<1> +FELinearMapping<1>::linear_shape_grad(const unsigned int i, + const Point<1>&) const +{ + Assert((i<2), ExcInvalidIndex(i)); + switch (i) + { + case 0: return Point<1>(-1.); + case 1: return Point<1>(1.); + } + return Point<1>(); +}; + + + +template <> +void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, + const Boundary<1> &, + const vector > &, + vector &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, + const unsigned int , + const vector > &, + vector &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const Boundary<1> &, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const unsigned int, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + +template <> +void FELinearMapping<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, + const vector > &unit_points, + vector &jacobians, + const bool compute_jacobians, + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points, + const Boundary<1> &boundary) const { + // simply pass down + FiniteElement<1>::fill_fe_values (cell, unit_points, + jacobians, compute_jacobians, + ansatz_points, compute_ansatz_points, + q_points, compute_q_points, boundary); +}; + + + +#endif + + + +#if deal_II_dimension == 2 + +template <> +inline +double +FELinearMapping<2>::linear_shape_value (const unsigned int i, + const Point<2>& p) const +{ + Assert((i<4), ExcInvalidIndex(i)); + switch (i) + { + case 0: return (1.-p(0)) * (1.-p(1)); + case 1: return p(0) * (1.-p(1)); + case 2: return p(0) * p(1); + case 3: return (1.-p(0)) * p(1); + } + return 0.; +}; + + + +template <> +inline +Point<2> +FELinearMapping<2>::linear_shape_grad (const unsigned int i, + const Point<2>& p) const +{ + Assert((i<4), ExcInvalidIndex(i)); + switch (i) + { + case 0: return Point<2> (p(1)-1., p(0)-1.); + case 1: return Point<2> (1.-p(1), -p(0)); + case 2: return Point<2> (p(1), p(0)); + case 3: return Point<2> (-p(1), 1.-p(0)); + } + return Point<2> (); +}; + + + +template <> +void FELinearMapping<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face, + const Boundary<2> &, + const vector > &unit_points, + vector &face_jacobians) const { + // more or less copied from the linear + // finite element + Assert (unit_points.size() == face_jacobians.size(), + ExcWrongFieldDimension (unit_points.size(), face_jacobians.size())); + + // 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); +}; + + + +template <> +void FELinearMapping<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face, + const unsigned int , + const vector > &unit_points, + vector &face_jacobians) const { + // more or less copied from the linear + // finite element + 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); +}; + + + +template <> +void FELinearMapping<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 { + // more or less copied from the linear + // finite element + 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()))); +}; + + + +template <> +void FELinearMapping<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 { + // more or less copied from the linear + // finite element + // 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()))); +}; + +#endif + + + +template +void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator &cell, + const vector > &unit_points, + vector &jacobians, + const bool compute_jacobians, + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points, + 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 (ansatz_points.size() == total_dofs, + ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); + + + unsigned int n_points=unit_points.size(); + + Point vertices[GeometryInfo::vertices_per_cell]; + for (unsigned int l=0; l::vertices_per_cell; ++l) + vertices[l] = cell->vertex(l); + + + if (compute_q_points) + { + // initialize points to zero + for (unsigned int i=0; i (); + + // note: let x_l be the vector of the + // lth quadrature point in real space and + // xi_l that on the unit cell, let further + // p_j be the vector of the jth vertex + // of the cell in real space and + // N_j(xi_l) be the value of the associated + // basis function at xi_l, then + // x_l(xi_l) = sum_j p_j N_j(xi_l) + // + // Here, N_j is the *linear* basis function, + // not that of the finite element, since we + // use a subparametric mapping + for (unsigned int j=0; j::vertices_per_cell; ++j) + for (unsigned int l=0; l::vertices_per_cell; ++s) + { + // we want the linear transform, + // so use that function + const Point gradient = linear_shape_grad (s, unit_points[l]); + for (unsigned int i=0; i; template class FiniteElementBase; template class FiniteElement; +template class FELinearMapping; diff --git a/deal.II/deal.II/source/fe/fe_lib.cubic.cc b/deal.II/deal.II/source/fe/fe_lib.cubic.cc index 82730253c3..2263f03d14 100644 --- a/deal.II/deal.II/source/fe/fe_lib.cubic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.cubic.cc @@ -303,7 +303,7 @@ template <> FECubicSub<1>::FECubicSub () : - FiniteElement<1> (1, 2) { + FELinearMapping<1> (1, 2) { prolongation[0](0,0) = 1.0; prolongation[0](0,1) = 0.0; prolongation[0](0,2) = 0.0; @@ -340,25 +340,6 @@ FECubicSub<1>::FECubicSub () : -template <> -void FECubicSub<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - const Boundary<1> &boundary) const { - // simply pass down - FiniteElement<1>::fill_fe_values (cell, unit_points, - jacobians, compute_jacobians, - ansatz_points, compute_ansatz_points, - q_points, compute_q_points, boundary); -}; - - - template <> double FECubicSub<1>::shape_value(const unsigned int i, @@ -378,24 +359,6 @@ FECubicSub<1>::shape_value(const unsigned int i, -template <> -inline -double -FECubicSub<1>::linear_shape_value(const unsigned int i, - const Point<1> &p) const -{ - Assert((i<2), ExcInvalidIndex(i)); - const double xi = p(0); - switch (i) - { - case 0: return 1.-xi; - case 1: return xi; - } - return 0.; -}; - - - template <> Point<1> FECubicSub<1>::shape_grad(const unsigned int i, @@ -415,23 +378,6 @@ FECubicSub<1>::shape_grad(const unsigned int i, -template <> -inline -Point<1> -FECubicSub<1>::linear_shape_grad(const unsigned int i, - const Point<1>&) const -{ - Assert((i<2), ExcInvalidIndex(i)); - switch (i) - { - case 0: return Point<1>(-1.); - case 1: return Point<1>(1.); - } - return Point<1>(); -}; - - - template <> void FECubicSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell, const Boundary<1> &boundary, @@ -450,48 +396,6 @@ void FECubicSub<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_i -template <> -void FECubicSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, - const Boundary<1> &, - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FECubicSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, - const unsigned int , - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FECubicSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const Boundary<1> &, - const vector > &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FECubicSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const unsigned int, - const vector > &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - - - template <> void FECubicSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell, const Boundary<1> &, @@ -537,7 +441,7 @@ void FECubicSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &c template <> FECubicSub<2>::FECubicSub () : - FiniteElement<2> (1, 2, 4) + FELinearMapping<2> (1, 2, 4) { interface_constraints(0,0) = -1.0/16.0; interface_constraints(0,1) = -1.0/16.0; @@ -1023,25 +927,6 @@ xi*xi)*eta*eta*eta; -template <> -inline -double -FECubicSub<2>::linear_shape_value (const unsigned int i, - const Point<2>& p) const -{ - Assert((i<4), ExcInvalidIndex(i)); - switch (i) - { - case 0: return (1.-p(0)) * (1.-p(1)); - case 1: return p(0) * (1.-p(1)); - case 2: return p(0) * p(1); - case 3: return (1.-p(0)) * p(1); - } - return 0.; -}; - - - template <> Point<2> FECubicSub<2>::shape_grad (const unsigned int i, @@ -1091,25 +976,6 @@ FECubicSub<2>::shape_grad (const unsigned int i, -template <> -inline -Point<2> -FECubicSub<2>::linear_shape_grad (const unsigned int i, - const Point<2>& p) const -{ - Assert((i<4), ExcInvalidIndex(i)); - switch (i) - { - case 0: return Point<2> (p(1)-1., p(0)-1.); - case 1: return Point<2> (1.-p(1), -p(0)); - case 2: return Point<2> (p(1), p(0)); - case 3: return Point<2> (-p(1), 1.-p(0)); - } - return Point<2> (); -}; - - - template <> void FECubicSub<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cell, const Boundary<2> &, @@ -1763,220 +1629,12 @@ void FECubicSub<2>::get_face_ansatz_points (const typename DoFHandler<2>::face_i -template <> -void FECubicSub<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face, - const Boundary<2> &, - const vector > &unit_points, - vector &face_jacobians) const { - // more or less copied from the linear - // finite element - Assert (unit_points.size() == face_jacobians.size(), - ExcWrongFieldDimension (unit_points.size(), face_jacobians.size())); - - // 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); -}; - - - -template <> -void FECubicSub<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face, - const unsigned int , - const vector > &unit_points, - vector &face_jacobians) const { - // more or less copied from the linear - // finite element - 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); -}; - - - -template <> -void FECubicSub<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 { - // more or less copied from the linear - // finite element - 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()))); -}; - - - -template <> -void FECubicSub<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 { - // more or less copied from the linear - // finite element - // 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()))); -}; - #endif -template -void FECubicSub::fill_fe_values (const DoFHandler::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - 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 (ansatz_points.size() == total_dofs, - ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); - - - unsigned int n_points=unit_points.size(); - - Point vertices[GeometryInfo::vertices_per_cell]; - for (unsigned int l=0; l::vertices_per_cell; ++l) - vertices[l] = cell->vertex(l); - - - if (compute_q_points) - { - // initialize points to zero - for (unsigned int i=0; i (); - - // note: let x_l be the vector of the - // lth quadrature point in real space and - // xi_l that on the unit cell, let further - // p_j be the vector of the jth vertex - // of the cell in real space and - // N_j(xi_l) be the value of the associated - // basis function at xi_l, then - // x_l(xi_l) = sum_j p_j N_j(xi_l) - // - // Here, N_j is the *linear* basis function, - // not that of the finite element, since we - // use a subparametric mapping - for (unsigned int j=0; j::vertices_per_cell; ++j) - for (unsigned int l=0; l::vertices_per_cell; ++s) - { - // we want the linear transform, - // so use that function - const Point gradient = linear_shape_grad (s, unit_points[l]); - for (unsigned int i=0; i FELinear<1>::FELinear () : - FiniteElement<1> (1, 0) + FELinearMapping<1> (1, 0) { // for restriction and prolongation matrices: // note that we do not add up all the @@ -79,25 +79,6 @@ FELinear<1>::shape_grad(const unsigned int i, -template <> -void FELinear<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - const Boundary<1> &boundary) const { - // simply pass down - FiniteElement<1>::fill_fe_values (cell, unit_points, - jacobians, compute_jacobians, - ansatz_points, compute_ansatz_points, - q_points, compute_q_points, boundary); -}; - - - template <> void FELinear<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell, const Boundary<1> &boundary, @@ -116,48 +97,6 @@ void FELinear<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_ite -template <> -void FELinear<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, - const Boundary<1> &, - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FELinear<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, - const unsigned int , - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FELinear<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const Boundary<1> &, - const vector > &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FELinear<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const unsigned int, - const vector > &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - - - template <> void FELinear<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell, const Boundary<1> &, @@ -183,7 +122,7 @@ void FELinear<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cel template <> FELinear<2>::FELinear () : - FiniteElement<2> (1, 0, 0) + FELinearMapping<2> (1, 0, 0) { interface_constraints(0,0) = 1./2.; interface_constraints(0,1) = 1./2.; @@ -413,209 +352,10 @@ void FELinear<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cel -template <> -void FELinear<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face, - const Boundary<2> &, - const vector > &unit_points, - vector &face_jacobians) const { - Assert (unit_points.size() == face_jacobians.size(), - ExcWrongFieldDimension (unit_points.size(), face_jacobians.size())); - - // 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); -}; - - - -template <> -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); -}; - - - -template <> -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()))); -}; - - - -template <> -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()))); -}; - #endif -template -void FELinear::fill_fe_values (const DoFHandler::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - 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 (ansatz_points.size() == total_dofs, - ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); - - unsigned int n_points=unit_points.size(); - - Point vertices[GeometryInfo::vertices_per_cell]; - for (unsigned int l=0; l::vertices_per_cell; ++l) - vertices[l] = cell->vertex(l); - - - if (compute_q_points) - { - // initialize points to zero - for (unsigned int i=0; i (); - - // note: let x_l be the vector of the - // lth quadrature point in real space and - // xi_l that on the unit cell, let further - // p_j be the vector of the jth vertex - // of the cell in real space and - // N_j(xi_l) be the value of the associated - // basis function at xi_l, then - // x_l(xi_l) = sum_j p_j N_j(xi_l) - for (unsigned int j=0; j::vertices_per_cell; ++j) - for (unsigned int l=0; l::shape_value(j, unit_points[l]); - }; - - -/* jacobi matrices: compute d(x)/d(xi) and invert this - Let M(l) be the inverse of J at the quadrature point l, then - M_{ij}(l) = sum_s p_i(s) d(N_s(l))/d(xi_j) - where p_i(s) is the i-th coordinate of the s-th vertex vector, - N_s(l) is the value of the s-th vertex shape function at the - quadrature point l. - - We could therefore write: - l=0..n_points-1 - i=0..dim-1 - j=0..dim-1 - M_{ij}(l) = 0 - s=0..n_vertices - M_{ij}(l) += p_i(s) d(N_s(l))/d(xi_j) - - However, we rewrite the loops to only compute the gradient once for - each integration point and basis function. -*/ - if (compute_jacobians) - { - dFMatrix M(dim,dim); - for (unsigned int l=0; l::vertices_per_cell; ++s) - { - // we want a linear transform and - // if we prepend the class name in - // front of the #shape_grad#, we - // need not use virtual function - // calls. - const Point gradient - = FELinear::shape_grad (s, unit_points[l]); - for (unsigned int i=0; i void FELinear::get_ansatz_points (const typename DoFHandler::cell_iterator &cell, const Boundary &, diff --git a/deal.II/deal.II/source/fe/fe_lib.quadratic.cc b/deal.II/deal.II/source/fe/fe_lib.quadratic.cc index 52884f1df1..2f10e9d1a8 100644 --- a/deal.II/deal.II/source/fe/fe_lib.quadratic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.quadratic.cc @@ -15,7 +15,7 @@ template <> FEQuadraticSub<1>::FEQuadraticSub () : - FiniteElement<1> (1, 1) { + FELinearMapping<1> (1, 1) { /* Get the prolongation matrices by the following little maple script: @@ -62,25 +62,6 @@ FEQuadraticSub<1>::FEQuadraticSub () : -template <> -void FEQuadraticSub<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - const Boundary<1> &boundary) const { - // simply pass down - FiniteElement<1>::fill_fe_values (cell, unit_points, - jacobians, compute_jacobians, - ansatz_points, compute_ansatz_points, - q_points, compute_q_points, boundary); -}; - - - template <> double FEQuadraticSub<1>::shape_value(const unsigned int i, @@ -99,24 +80,6 @@ FEQuadraticSub<1>::shape_value(const unsigned int i, -template <> -inline -double -FEQuadraticSub<1>::linear_shape_value(const unsigned int i, - const Point<1> &p) const -{ - Assert((i<2), ExcInvalidIndex(i)); - const double xi = p(0); - switch (i) - { - case 0: return 1.-xi; - case 1: return xi; - } - return 0.; -}; - - - template <> Point<1> FEQuadraticSub<1>::shape_grad(const unsigned int i, @@ -135,23 +98,6 @@ FEQuadraticSub<1>::shape_grad(const unsigned int i, -template <> -inline -Point<1> -FEQuadraticSub<1>::linear_shape_grad(const unsigned int i, - const Point<1>&) const -{ - Assert((i<2), ExcInvalidIndex(i)); - switch (i) - { - case 0: return Point<1>(-1.); - case 1: return Point<1>(1.); - } - return Point<1>(); -}; - - - template <> void FEQuadraticSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell, const Boundary<1> &boundary, @@ -170,48 +116,6 @@ void FEQuadraticSub<1>::get_face_ansatz_points (const typename DoFHandler<1>::fa -template <> -void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, - const Boundary<1> &, - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, - const unsigned int , - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const Boundary<1> &, - const vector > &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const unsigned int, - const vector > &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - - - template <> void FEQuadraticSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell, const Boundary<1> &, @@ -239,7 +143,7 @@ void FEQuadraticSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterato template <> FEQuadraticSub<2>::FEQuadraticSub () : - FiniteElement<2> (1, 1, 1) + FELinearMapping<2> (1, 1, 1) { interface_constraints(0,2) = 1.0; interface_constraints(1,0) = 3./8.; @@ -641,25 +545,6 @@ FEQuadraticSub<2>::shape_value (const unsigned int i, -template <> -inline -double -FEQuadraticSub<2>::linear_shape_value (const unsigned int i, - const Point<2>& p) const -{ - Assert((i<4), ExcInvalidIndex(i)); - switch (i) - { - case 0: return (1.-p(0)) * (1.-p(1)); - case 1: return p(0) * (1.-p(1)); - case 2: return p(0) * p(1); - case 3: return (1.-p(0)) * p(1); - } - return 0.; -}; - - - template <> Point<2> FEQuadraticSub<2>::shape_grad (const unsigned int i, @@ -695,25 +580,6 @@ FEQuadraticSub<2>::shape_grad (const unsigned int i, -template <> -inline -Point<2> -FEQuadraticSub<2>::linear_shape_grad (const unsigned int i, - const Point<2>& p) const -{ - Assert((i<4), ExcInvalidIndex(i)); - switch (i) - { - case 0: return Point<2> (p(1)-1., p(0)-1.); - case 1: return Point<2> (1.-p(1), -p(0)); - case 2: return Point<2> (p(1), p(0)); - case 3: return Point<2> (-p(1), 1.-p(0)); - } - return Point<2> (); -}; - - - template <> void FEQuadraticSub<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cell, const Boundary<2> &, @@ -1033,220 +899,12 @@ void FEQuadraticSub<2>::get_face_ansatz_points (const typename DoFHandler<2>::fa -template <> -void FEQuadraticSub<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face, - const Boundary<2> &, - const vector > &unit_points, - vector &face_jacobians) const { - // more or less copied from the linear - // finite element - Assert (unit_points.size() == face_jacobians.size(), - ExcWrongFieldDimension (unit_points.size(), face_jacobians.size())); - - // 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); -}; - - - -template <> -void FEQuadraticSub<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face, - const unsigned int , - const vector > &unit_points, - vector &face_jacobians) const { - // more or less copied from the linear - // finite element - 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); -}; - - - -template <> -void FEQuadraticSub<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 { - // more or less copied from the linear - // finite element - 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()))); -}; - - - -template <> -void FEQuadraticSub<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 { - // more or less copied from the linear - // finite element - // 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()))); -}; - #endif -template -void FEQuadraticSub::fill_fe_values (const DoFHandler::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - 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 (ansatz_points.size() == total_dofs, - ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); - - - unsigned int n_points=unit_points.size(); - - Point vertices[GeometryInfo::vertices_per_cell]; - for (unsigned int l=0; l::vertices_per_cell; ++l) - vertices[l] = cell->vertex(l); - - - if (compute_q_points) - { - // initialize points to zero - for (unsigned int i=0; i (); - - // note: let x_l be the vector of the - // lth quadrature point in real space and - // xi_l that on the unit cell, let further - // p_j be the vector of the jth vertex - // of the cell in real space and - // N_j(xi_l) be the value of the associated - // basis function at xi_l, then - // x_l(xi_l) = sum_j p_j N_j(xi_l) - // - // Here, N_j is the *linear* basis function, - // not that of the finite element, since we - // use a subparametric mapping - for (unsigned int j=0; j::vertices_per_cell; ++j) - for (unsigned int l=0; l::vertices_per_cell; ++s) - { - // we want the linear transform, - // so use that function - const Point gradient = linear_shape_grad (s, unit_points[l]); - for (unsigned int i=0; i; diff --git a/deal.II/deal.II/source/fe/fe_lib.quartic.cc b/deal.II/deal.II/source/fe/fe_lib.quartic.cc index fb87596789..71ad977592 100644 --- a/deal.II/deal.II/source/fe/fe_lib.quartic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.quartic.cc @@ -320,7 +320,7 @@ template <> FEQuarticSub<1>::FEQuarticSub () : - FiniteElement<1> (1, 3) { + FELinearMapping<1> (1, 3) { prolongation[0](0,0) = 1.0; prolongation[0](1,3) = 1.0; prolongation[0](2,0) = 35.0/128.0; @@ -351,25 +351,6 @@ FEQuarticSub<1>::FEQuarticSub () : -template <> -void FEQuarticSub<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - const Boundary<1> &boundary) const { - // simply pass down - FiniteElement<1>::fill_fe_values (cell, unit_points, - jacobians, compute_jacobians, - ansatz_points, compute_ansatz_points, - q_points, compute_q_points, boundary); -}; - - - template <> double FEQuarticSub<1>::shape_value(const unsigned int i, @@ -390,24 +371,6 @@ FEQuarticSub<1>::shape_value(const unsigned int i, -template <> -inline -double -FEQuarticSub<1>::linear_shape_value(const unsigned int i, - const Point<1> &p) const -{ - Assert((i<2), ExcInvalidIndex(i)); - const double xi = p(0); - switch (i) - { - case 0: return 1.-xi; - case 1: return xi; - } - return 0.; -}; - - - template <> Point<1> FEQuarticSub<1>::shape_grad(const unsigned int i, @@ -428,23 +391,6 @@ FEQuarticSub<1>::shape_grad(const unsigned int i, -template <> -inline -Point<1> -FEQuarticSub<1>::linear_shape_grad(const unsigned int i, - const Point<1>&) const -{ - Assert((i<2), ExcInvalidIndex(i)); - switch (i) - { - case 0: return Point<1>(-1.); - case 1: return Point<1>(1.); - } - return Point<1>(); -}; - - - template <> void FEQuarticSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell, const Boundary<1> &boundary, @@ -463,48 +409,6 @@ void FEQuarticSub<1>::get_face_ansatz_points (const typename DoFHandler<1>::face -template <> -void FEQuarticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, - const Boundary<1> &, - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FEQuarticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, - const unsigned int , - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FEQuarticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const Boundary<1> &, - const vector > &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FEQuarticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const unsigned int, - const vector > &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - - - template <> void FEQuarticSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell, const Boundary<1> &, @@ -561,7 +465,7 @@ void FEQuarticSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator template <> FEQuarticSub<2>::FEQuarticSub () : - FiniteElement<2> (1, 3, 9) + FELinearMapping<2> (1, 3, 9) { interface_constraints(0,3) = 1.0; interface_constraints(1,0) = 35.0/128.0; @@ -1414,25 +1318,6 @@ FEQuarticSub<2>::shape_value (const unsigned int i, -template <> -inline -double -FEQuarticSub<2>::linear_shape_value (const unsigned int i, - const Point<2>& p) const -{ - Assert((i<4), ExcInvalidIndex(i)); - switch (i) - { - case 0: return (1.-p(0)) * (1.-p(1)); - case 1: return p(0) * (1.-p(1)); - case 2: return p(0) * p(1); - case 3: return (1.-p(0)) * p(1); - } - return 0.; -}; - - - template <> Point<2> FEQuarticSub<2>::shape_grad (const unsigned int i, @@ -1500,25 +1385,6 @@ FEQuarticSub<2>::shape_grad (const unsigned int i, -template <> -inline -Point<2> -FEQuarticSub<2>::linear_shape_grad (const unsigned int i, - const Point<2>& p) const -{ - Assert((i<4), ExcInvalidIndex(i)); - switch (i) - { - case 0: return Point<2> (p(1)-1., p(0)-1.); - case 1: return Point<2> (1.-p(1), -p(0)); - case 2: return Point<2> (p(1), p(0)); - case 3: return Point<2> (-p(1), 1.-p(0)); - } - return Point<2> (); -}; - - - template <> void FEQuarticSub<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cell, const Boundary<2> &, @@ -2912,219 +2778,12 @@ void FEQuarticSub<2>::get_face_ansatz_points (const typename DoFHandler<2>::face -template <> -void FEQuarticSub<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face, - const Boundary<2> &, - const vector > &unit_points, - vector &face_jacobians) const { - // more or less copied from the linear - // finite element - Assert (unit_points.size() == face_jacobians.size(), - ExcWrongFieldDimension (unit_points.size(), face_jacobians.size())); - - // 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); -}; - - - -template <> -void FEQuarticSub<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face, - const unsigned int , - const vector > &unit_points, - vector &face_jacobians) const { - // more or less copied from the linear - // finite element - 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); -}; - - - -template <> -void FEQuarticSub<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 { - // more or less copied from the linear - // finite element - 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()))); -}; - - - -template <> -void FEQuarticSub<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 { - // more or less copied from the linear - // finite element - // 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()))); -}; - #endif -template -void FEQuarticSub::fill_fe_values (const DoFHandler::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - 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 (ansatz_points.size() == total_dofs, - ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); - - - unsigned int n_points=unit_points.size(); - - Point vertices[GeometryInfo::vertices_per_cell]; - for (unsigned int l=0; l::vertices_per_cell; ++l) - vertices[l] = cell->vertex(l); - - - if (compute_q_points) - { - // initialize points to zero - for (unsigned int i=0; i (); - - // note: let x_l be the vector of the - // lth quadrature point in real space and - // xi_l that on the unit cell, let further - // p_j be the vector of the jth vertex - // of the cell in real space and - // N_j(xi_l) be the value of the associated - // basis function at xi_l, then - // x_l(xi_l) = sum_j p_j N_j(xi_l) - // - // Here, N_j is the *linear* basis function, - // not that of the finite element, since we - // use a subparametric mapping - for (unsigned int j=0; j::vertices_per_cell; ++j) - for (unsigned int l=0; l::vertices_per_cell; ++s) - { - // we want the linear transform, - // so use that function - const Point gradient = linear_shape_grad (s, unit_points[l]); - for (unsigned int i=0; i