From e18cab567976676e1c91bdb9eead3d43d57180b7 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 22 Jul 1998 08:26:16 +0000 Subject: [PATCH] Separate the linear mapping from the rest of the finite element stuff. This is a Good Thing (TM) since most otehr parts of the library don't need specifics about a finite element, but only the general FiniteElement class. git-svn-id: https://svn.dealii.org/trunk@452 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe.h | 137 ------- deal.II/deal.II/include/fe/fe_lib.lagrange.h | 2 +- .../deal.II/include/fe/fe_linear_mapping.h | 154 ++++++++ deal.II/deal.II/source/fe/fe.cc | 359 +---------------- .../deal.II/source/fe/fe_linear_mapping.cc | 374 ++++++++++++++++++ 5 files changed, 530 insertions(+), 496 deletions(-) create mode 100644 deal.II/deal.II/include/fe/fe_linear_mapping.h create mode 100644 deal.II/deal.II/source/fe/fe_linear_mapping.cc diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 0dcc0001f4..fbe50a070f 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -1041,143 +1041,6 @@ 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, - GeometryInfo::vertices_per_cell) {}; - - /** - * 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. - */ - virtual double shape_value_transform (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. - */ - virtual Point shape_grad_transform (const unsigned int i, - const Point &p) 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. - * - * 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 dFMatrix &shape_values_transform, - const vector > > &shape_grad_transform, - const Boundary &boundary) 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 608dd3345f..22da61caf8 100644 --- a/deal.II/deal.II/include/fe/fe_lib.lagrange.h +++ b/deal.II/deal.II/include/fe/fe_lib.lagrange.h @@ -5,7 +5,7 @@ /*---------------------------- fe_lib.h ---------------------------*/ -#include +#include diff --git a/deal.II/deal.II/include/fe/fe_linear_mapping.h b/deal.II/deal.II/include/fe/fe_linear_mapping.h new file mode 100644 index 0000000000..86f2935a98 --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_linear_mapping.h @@ -0,0 +1,154 @@ +/*---------------------------- fe_linear_mapping.h ---------------------------*/ +/* $Id$ */ +#ifndef __fe_linear_mapping_H +#define __fe_linear_mapping_H +/*---------------------------- fe_linear_mapping.h ---------------------------*/ + + +#include + + + + +/** + * 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, + GeometryInfo::vertices_per_cell) {}; + + /** + * 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. + */ + virtual double shape_value_transform (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. + */ + virtual Point shape_grad_transform (const unsigned int i, + const Point &p) 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. + * + * 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 dFMatrix &shape_values_transform, + const vector > > &shape_grad_transform, + const Boundary &boundary) const; +}; + + + + +/*---------------------------- fe_linear_mapping.h ---------------------------*/ +/* end of #ifndef __fe_linear_mapping_H */ +#endif +/*---------------------------- fe_linear_mapping.h ---------------------------*/ diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index ece34091c5..449de51e7c 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -434,368 +434,11 @@ void FiniteElement::get_ansatz_points (const DoFHandler::cell_iterator -/*---------------------------- FELinearMapping ----------------------------------*/ - -#if deal_II_dimension == 1 - -template <> -inline -double -FELinearMapping<1>::shape_value_transform (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>::shape_grad_transform(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 FELinearMapping<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, - const Boundary<1> &, - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FELinearMapping<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, - const unsigned int , - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FELinearMapping<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const Boundary<1> &, - const vector > &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FELinearMapping<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 dFMatrix &shape_values_transform, - const vector > > &shape_gradients_transform, - 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, - shape_values_transform, shape_gradients_transform, - boundary); -}; - - - -#endif - - - -#if deal_II_dimension == 2 - -template <> -inline -double -FELinearMapping<2>::shape_value_transform (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>::shape_grad_transform (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 dFMatrix &shape_values_transform, - const vector > > &shape_grad_transform, - 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 = shape_grad_transform[s][l]; - for (unsigned int i=0; i; template class FiniteElement; -template class FELinearMapping; + diff --git a/deal.II/deal.II/source/fe/fe_linear_mapping.cc b/deal.II/deal.II/source/fe/fe_linear_mapping.cc new file mode 100644 index 0000000000..4f86319dc6 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_linear_mapping.cc @@ -0,0 +1,374 @@ +/* $Id$ */ + +#include +#include +#include +#include +#include + + + + + +/*---------------------------- FELinearMapping ----------------------------------*/ + +#if deal_II_dimension == 1 + +template <> +inline +double +FELinearMapping<1>::shape_value_transform (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>::shape_grad_transform(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 FELinearMapping<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, + const Boundary<1> &, + const vector > &, + vector &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FELinearMapping<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, + const unsigned int , + const vector > &, + vector &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FELinearMapping<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const Boundary<1> &, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FELinearMapping<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 dFMatrix &shape_values_transform, + const vector > > &shape_gradients_transform, + 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, + shape_values_transform, shape_gradients_transform, + boundary); +}; + + + +#endif + + + +#if deal_II_dimension == 2 + +template <> +inline +double +FELinearMapping<2>::shape_value_transform (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>::shape_grad_transform (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 dFMatrix &shape_values_transform, + const vector > > &shape_grad_transform, + 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 = shape_grad_transform[s][l]; + for (unsigned int i=0; i; -- 2.39.5