From: Wolfgang Bangerth Date: Thu, 16 Jul 1998 16:08:17 +0000 (+0000) Subject: Implement a faster scheme to compute the FEValues by storing the values and gradients... X-Git-Tag: v8.0.0~22807 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=436954d00e87184d0a24cb3aa5590311038fecb5;p=dealii.git Implement a faster scheme to compute the FEValues by storing the values and gradients of the transformation basis fnctions at teh quadrature points only once at the initialization rather than on each cell. git-svn-id: https://svn.dealii.org/trunk@447 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 cfe3502f39..87c4b073e3 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -53,28 +53,34 @@ struct FiniteElementData<1> { */ const unsigned int total_dofs; + /** + * Number of basis functions used for the + * transformation from unit cell to real + * cell. For a linear mapping, this number + * equals the number of vertices. + */ + const unsigned int n_transform_functions; + /** * Default constructor. Constructs an element * which is not so useful. Checking * #total_dofs# is therefore a good way to * check if something went wrong. */ - FiniteElementData () : - dofs_per_vertex(0), - dofs_per_line(0), - dofs_per_face(0), - total_dofs(0) {}; + FiniteElementData (); /** * A more useful version to construct * an object of this type. */ FiniteElementData (const unsigned int dofs_per_vertex, - const unsigned int dofs_per_line) : + const unsigned int dofs_per_line, + const unsigned int n_transform_functions) : dofs_per_vertex(dofs_per_vertex), dofs_per_line(dofs_per_line), dofs_per_face(dofs_per_vertex), - total_dofs (2*dofs_per_vertex+dofs_per_line) {}; + total_dofs (2*dofs_per_vertex+dofs_per_line), + n_transform_functions(n_transform_functions) {}; /** * Comparison operator. It is not clear to @@ -82,6 +88,11 @@ struct FiniteElementData<1> { * this one explicitely. */ bool operator == (const FiniteElementData<1> &) const; + + /** + * Exception + */ + DeclException0 (ExcInternalError); }; @@ -128,18 +139,21 @@ struct FiniteElementData<2> { */ const unsigned int total_dofs; + /** + * Number of basis functions used for the + * transformation from unit cell to real + * cell. For a linear mapping, this number + * equals the number of vertices. + */ + const unsigned int n_transform_functions; + /** * Default constructor. Constructs an element * which is not so useful. Checking * #total_dofs# is therefore a good way to * check if something went wrong. */ - FiniteElementData () : - dofs_per_vertex(0), - dofs_per_line(0), - dofs_per_quad(0), - dofs_per_face(0), - total_dofs(0) {}; + FiniteElementData (); /** * A more useful version to construct @@ -147,7 +161,8 @@ struct FiniteElementData<2> { */ FiniteElementData (const unsigned int dofs_per_vertex, const unsigned int dofs_per_line, - const unsigned int dofs_per_quad) : + const unsigned int dofs_per_quad, + const unsigned int n_transform_functions) : dofs_per_vertex(dofs_per_vertex), dofs_per_line(dofs_per_line), dofs_per_quad(dofs_per_quad), @@ -155,7 +170,8 @@ struct FiniteElementData<2> { dofs_per_line), total_dofs (4*dofs_per_vertex+ 4*dofs_per_line+ - dofs_per_quad) {}; + dofs_per_quad), + n_transform_functions (n_transform_functions) {}; /** * Comparison operator. It is not clear to @@ -163,6 +179,11 @@ struct FiniteElementData<2> { * this one explicitely. */ bool operator == (const FiniteElementData<2> &) const; + + /** + * Exception + */ + DeclException0 (ExcInternalError); }; @@ -202,7 +223,8 @@ struct FiniteElementBase : public FiniteElementData { */ FiniteElementBase (const unsigned int dofs_per_vertex, const unsigned int dofs_per_line, - const unsigned int dofs_per_quad=0); + const unsigned int dofs_per_quad, + const unsigned int n_transform_functions); /** * Return a readonly reference to the @@ -266,10 +288,6 @@ struct FiniteElementBase : public FiniteElementData { << "The interface matrix has a size of " << arg1 << "x" << arg2 << ", which is not reasonable in the present dimension."); - /** - * Exception - */ - DeclException0 (ExcInternalError); protected: /** @@ -474,10 +492,12 @@ class FiniteElement : public FiniteElementBase { */ FiniteElement (const unsigned int dofs_per_vertex, const unsigned int dofs_per_line, - const unsigned int dofs_per_quad=0) : + const unsigned int dofs_per_quad, + const unsigned int n_transform_functions) : FiniteElementBase (dofs_per_vertex, dofs_per_line, - dofs_per_quad) {}; + dofs_per_quad, + n_transform_functions) {}; /** * Destructor. Only declared to have a @@ -492,7 +512,7 @@ class FiniteElement : public FiniteElementBase { * #p# is a point on the reference element. */ virtual double shape_value (const unsigned int i, - const Point& p) const = 0; + const Point &p) const = 0; /** * Return the gradient of the #i#th shape @@ -500,8 +520,28 @@ class FiniteElement : public FiniteElementBase { * #p# is a point on the reference element, */ virtual Point shape_grad (const unsigned int i, - const Point& p) const = 0; + const Point &p) const = 0; + /** + * Return the value of the #i#th shape + * function of the transformation mapping + * from unit cell to real cell. For + * isoparametric elements, this function + * is the same as the ansatz functions, + * but for sublinear or other mappings, + * they differ. + */ + virtual double shape_value_transform (const unsigned int i, + const Point &p) const = 0; + + /** + * Same as above: return gradient of the + * #i#th shape function for the mapping + * from unit to real cell. + */ + virtual Point shape_grad_transform (const unsigned int i, + const Point &p) const = 0; + /** * Compute the Jacobian matrix and the * quadrature points as well as the ansatz @@ -538,6 +578,17 @@ class FiniteElement : public FiniteElementBase { * unit cell to the real cell, which * makes things much more complicated. * + * The #shape_values/grads_transform# + * arrays store the values and gradients + * of the transformation basis functions. + * While this information is not necessary + * for the computation of the other fields, + * it allows for significant speedups, since + * the values and gradients of the transform + * functions at the quadrature points + * need not be recomputed each time this + * function is called. + * * The function assumes that the fields * already have the right number of * elements. @@ -555,6 +606,8 @@ class FiniteElement : public FiniteElementBase { const bool compute_ansatz_points, vector > &q_points, const bool compute_q_points, + const dFMatrix &shape_values_transform, + const vector > > &shape_grads_transform, const Boundary &boundary) const; /** @@ -675,6 +728,8 @@ class FiniteElement : public FiniteElementBase { const bool compute_face_jacobians, vector > &normal_vectors, const bool compute_normal_vectors, + const dFMatrix &shape_values_transform, + const vector > > &shape_grads_transform, const Boundary &boundary) const; /** @@ -718,6 +773,8 @@ class FiniteElement : public FiniteElementBase { const bool compute_face_jacobians, vector > &normal_vectors, const bool compute_normal_vectors, + const dFMatrix &shape_values_transform, + const vector > > &shape_grads_transform, const Boundary &boundary) const; /** @@ -984,7 +1041,30 @@ class FELinearMapping : public FiniteElement { const unsigned int dofs_per_quad=0) : FiniteElement (dofs_per_vertex, dofs_per_line, - dofs_per_quad) {}; + 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 @@ -1073,30 +1153,9 @@ class FELinearMapping : public FiniteElement { 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; - - 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; }; diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 7b153329a3..fc156f0f1f 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -57,6 +57,32 @@ template class Quadrature; * the derived classes and the #get_values# function for the exact usage of * this variable. * + * For many of the actual computations done by the #fill_fe_*# functions of + * the #FiniteElement# class and its decendants, the values and gradients of + * the transformation functions are needed. For example, for the computation + * of the real location of a quadrature point from the location on the unit + * cell, the values are needed, while for the computation of the Jacobian + * matrix the gradient is needed. While for linear elements the transformation + * functions coincide with the ansatz functions, this does not hold for higher + * order elements with subparametric mappings and for other types of elements + * such as non-conforming ones, etc, such that the precomputed values and + * gradients of the ansatz functions (#unit_shape_values# and + * #unit_shape_grads#) cannot be used for the present purpose. + * In principle, these values could be computed each time the #fill_fe_*# + * function is called; however, this computation is highly redundant, since + * only the values on the unit cell and only at the quadrature points are + * needed, i.e. they are the same for each cell that #fill_fe_*# is called. + * Therefore, two additional arrays, #unit_shape_values_transform# and + * #unit_shape_grads_transform# are provided, which are filled upon construction + * of an object of this type, which the actual finite element may or may not + * use. Later on, the #fill_fe_*# functions are passed pointers to these + * arrays, which they may use to extract the values and gradients of the + * transform functions. If a concrete finite element choses not to use this + * field, it shall set its field #n_transform_functions# to zero. + * + * The #unit_shape_grads_transform# array is provided by the derived classes + * to allow for the inclusion of multiple faces, etc. + * * * \subsection{Definitions} * @@ -141,7 +167,6 @@ template class Quadrature; * not assume that someone who wants to get the #JxW_values# must set the * #update_jacobians# flag besides the #update_JxW_values# flag. * - * * It is also forbidden that the constructor of a class set the * #update_jacobians# flag if the user specifies #update_JxW_values#. This is * since derived classes may be able to compute the #JxW_values# field without @@ -169,6 +194,14 @@ class FEValuesBase { */ const unsigned int total_dofs; + /** + * Number of basis functions for the + * transformation from the unit cell + * to the real cell. See the docs for + * more information on this field. + */ + const unsigned int n_transform_functions; + /** * Constructor. Set up the array sizes * with #n_q_points# quadrature points, @@ -186,6 +219,7 @@ class FEValuesBase { FEValuesBase (const unsigned int n_q_points, const unsigned int n_ansatz_points, const unsigned int n_dofs, + const unsigned int n_transform_functions, const unsigned int n_values_array, const UpdateFlags update_flags); @@ -470,6 +504,27 @@ class FEValuesBase { */ vector jacobi_matrices; + /** + * Store the values of the basis functions + * of the transformation from unit cell + * to real cell at the quadrature points. + * + * This field stores some data which is not + * really needed for the assemblage of + * matrices and vectors but makes that + * operation much faster. Each time the + * #FEValues::reinit# function calls + * the #FiniteElemenet::fill_fe_values# + * function, this and the next array are + * passed. The #fill_fe_values# function + * may or may not use this field. + * + * The element #(i,j)# denotes the value + * of the #i#th transfer basis function + * at the #j#th quadrature point. + */ + vector shape_values_transform; + /** * Store which of the data sets in the * #shape_values# array is presently @@ -484,7 +539,7 @@ class FEValuesBase { * Store which fields are to be updated by * the reinit function. */ - UpdateFlags update_flags; + UpdateFlags update_flags; /** * Store the cell selected last time @@ -575,7 +630,14 @@ class FEValues : public FEValuesBase { */ vector > > unit_shape_gradients; - + /** + * Gradients of the basis + * functions of the transformation. + * Analogous to the #shape_values_transform# + * array of the base class. + */ + vector > > unit_shape_gradients_transform; + /** * Array of quadrature points in the unit * cell. This array is set up upon @@ -680,6 +742,7 @@ class FEFaceValuesBase : public FEValuesBase { FEFaceValuesBase (const unsigned int n_q_points, const unsigned int n_ansatz_points, const unsigned int n_dofs, + const unsigned int n_transform_functions, const unsigned int n_faces_or_subfaces, const UpdateFlags update_flags); @@ -712,6 +775,14 @@ class FEFaceValuesBase : public FEValuesBase { * #unit_shape_gradients[face][dof][q_point]# */ vector > > > unit_shape_gradients; + + /** + * Gradients of the basis + * functions of the transformation. + * Analogous to the #shape_values_transform# + * array of the base class. + */ + vector > > > unit_shape_gradients_transform; /** * Array of quadrature points on the diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index fefe9e5745..cc8e4a672e 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -14,6 +14,17 @@ #if deal_II_dimension == 1 +FiniteElementData<1>::FiniteElementData () : + dofs_per_vertex(0), + dofs_per_line(0), + dofs_per_face(0), + total_dofs(0), + n_transform_functions(0) { + Assert (false, ExcInternalError()); +}; + + + bool FiniteElementData<1>::operator== (const FiniteElementData<1> &f) const { return ((dofs_per_vertex == f.dofs_per_vertex) && (dofs_per_line == f.dofs_per_line) && @@ -26,6 +37,19 @@ bool FiniteElementData<1>::operator== (const FiniteElementData<1> &f) const { #if deal_II_dimension == 2 + +FiniteElementData<2>::FiniteElementData () : + dofs_per_vertex(0), + dofs_per_line(0), + dofs_per_quad(0), + dofs_per_face(0), + total_dofs(0), + n_transform_functions(0) { + Assert (false, ExcInternalError()); +}; + + + bool FiniteElementData<2>::operator== (const FiniteElementData<2> &f) const { return ((dofs_per_vertex == f.dofs_per_vertex) && (dofs_per_line == f.dofs_per_line) && @@ -45,9 +69,11 @@ bool FiniteElementData<2>::operator== (const FiniteElementData<2> &f) const { template <> FiniteElementBase<1>::FiniteElementBase (const unsigned int dofs_per_vertex, const unsigned int dofs_per_line, - const unsigned int dofs_per_quad) : + const unsigned int dofs_per_quad, + const unsigned int n_transform_funcs) : FiniteElementData<1> (dofs_per_vertex, - dofs_per_line) + dofs_per_line, + n_transform_functs) { Assert (dofs_per_quad==0, ExcInternalError()); @@ -69,10 +95,12 @@ FiniteElementBase<1>::FiniteElementBase (const unsigned int dofs_per_vertex, template <> FiniteElementBase<2>::FiniteElementBase (const unsigned int dofs_per_vertex, const unsigned int dofs_per_line, - const unsigned int dofs_per_quad) : + const unsigned int dofs_per_quad, + const unsigned int n_transform_funcs) : FiniteElementData<2> (dofs_per_vertex, dofs_per_line, - dofs_per_quad) + dofs_per_quad, + n_transform_funcs) { const unsigned int dim=2; for (unsigned int i=0; i::children_per_cell; ++i) @@ -153,6 +181,8 @@ void FiniteElement<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, 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<1> &boundary) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); @@ -198,6 +228,8 @@ void FiniteElement<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator & const bool , vector > &, const bool, + const dFMatrix &, + const vector > > &, const Boundary<1> &) const { Assert (false, ExcNotImplemented()); }; @@ -217,6 +249,8 @@ void FiniteElement<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterato const bool , vector > &, const bool, + const dFMatrix &, + const vector > > &, const Boundary<1> &) const { Assert (false, ExcNotImplemented()); }; @@ -260,6 +294,8 @@ void FiniteElement::fill_fe_values (const DoFHandler::cell_iterator &, const bool, vector > &, const bool, + const dFMatrix &, + const vector > > &, const Boundary &) const { Assert (false, ExcPureFunctionCalled()); }; @@ -281,6 +317,8 @@ void FiniteElement::fill_fe_face_values (const DoFHandler::cell_iterat const bool compute_face_jacobians, vector > &normal_vectors, const bool compute_normal_vectors, + const dFMatrix &shape_values_transform, + const vector > > &shape_gradients_transform, const Boundary &boundary) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); @@ -296,6 +334,7 @@ void FiniteElement::fill_fe_face_values (const DoFHandler::cell_iterat jacobians, compute_jacobians, dummy, false, q_points, compute_q_points, + shape_values_transform, shape_gradients_transform, boundary); if (compute_ansatz_points) @@ -327,6 +366,8 @@ void FiniteElement::fill_fe_subface_values (const DoFHandler::cell_ite const bool compute_face_jacobians, vector > &normal_vectors, const bool compute_normal_vectors, + const dFMatrix &shape_values_transform, + const vector > > &shape_gradients_transform, const Boundary &boundary) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); @@ -340,6 +381,7 @@ void FiniteElement::fill_fe_subface_values (const DoFHandler::cell_ite jacobians, compute_jacobians, dummy, false, q_points, compute_q_points, + shape_values_transform, shape_gradients_transform, boundary); if (compute_face_jacobians) @@ -362,13 +404,16 @@ void FiniteElement::get_ansatz_points (const DoFHandler::cell_iterator + +/*---------------------------- FELinearMapping ----------------------------------*/ + #if deal_II_dimension == 1 template <> inline double -FELinearMapping<1>::linear_shape_value(const unsigned int i, - const Point<1> &p) const +FELinearMapping<1>::shape_value_transform (const unsigned int i, + const Point<1> &p) const { Assert((i<2), ExcInvalidIndex(i)); const double xi = p(0); @@ -385,8 +430,8 @@ FELinearMapping<1>::linear_shape_value(const unsigned int i, template <> inline Point<1> -FELinearMapping<1>::linear_shape_grad(const unsigned int i, - const Point<1>&) const +FELinearMapping<1>::shape_grad_transform(const unsigned int i, + const Point<1>&) const { Assert((i<2), ExcInvalidIndex(i)); switch (i) @@ -400,7 +445,7 @@ FELinearMapping<1>::linear_shape_grad(const unsigned int i, template <> -void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, +void FELinearMapping<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, const Boundary<1> &, const vector > &, vector &) const { @@ -410,7 +455,7 @@ void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator & template <> -void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, +void FELinearMapping<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, const unsigned int , const vector > &, vector &) const { @@ -420,7 +465,7 @@ void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterato template <> -void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, +void FELinearMapping<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, const unsigned int, const Boundary<1> &, const vector > &, @@ -431,7 +476,7 @@ void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator & template <> -void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, +void FELinearMapping<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, const unsigned int, const unsigned int, const vector > &, @@ -449,12 +494,16 @@ void FELinearMapping<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cel 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, boundary); + q_points, compute_q_points, + shape_values_transform, shape_gradients_transform, + boundary); }; @@ -468,8 +517,8 @@ void FELinearMapping<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cel template <> inline double -FELinearMapping<2>::linear_shape_value (const unsigned int i, - const Point<2>& p) const +FELinearMapping<2>::shape_value_transform (const unsigned int i, + const Point<2>& p) const { Assert((i<4), ExcInvalidIndex(i)); switch (i) @@ -487,8 +536,8 @@ FELinearMapping<2>::linear_shape_value (const unsigned int i, template <> inline Point<2> -FELinearMapping<2>::linear_shape_grad (const unsigned int i, - const Point<2>& p) const +FELinearMapping<2>::shape_grad_transform (const unsigned int i, + const Point<2>& p) const { Assert((i<4), ExcInvalidIndex(i)); switch (i) @@ -629,6 +678,8 @@ void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator 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())); @@ -665,7 +716,7 @@ void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator // use a subparametric mapping for (unsigned int j=0; j::vertices_per_cell; ++j) for (unsigned int l=0; l::fill_fe_values (const DoFHandler::cell_iterator { // we want the linear transform, // so use that function - const Point gradient = linear_shape_grad (s, unit_points[l]); + const Point gradient = shape_grad_transform[s][l]; for (unsigned int i=0; i::fill_fe_values (const DoFHandler::cell_iterator /*------------------------------- Explicit Instantiations -------------*/ -template class FiniteElementData; template class FiniteElementBase; template class FiniteElement; template class FELinearMapping; diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index afa2c70791..0a8c85e5de 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -17,10 +17,12 @@ template FEValuesBase::FEValuesBase (const unsigned int n_q_points, const unsigned int n_ansatz_points, const unsigned int n_dofs, + const unsigned int n_transform_functions, const unsigned int n_values_arrays, const UpdateFlags update_flags) : n_quadrature_points (n_q_points), total_dofs (n_dofs), + n_transform_functions (n_transform_functions), shape_values (n_values_arrays, dFMatrix(n_dofs, n_q_points)), shape_gradients (n_dofs, vector >(n_q_points)), weights (n_q_points, 0), @@ -28,6 +30,9 @@ FEValuesBase::FEValuesBase (const unsigned int n_q_points, quadrature_points (n_q_points, Point()), ansatz_points (n_ansatz_points, Point()), jacobi_matrices (n_q_points, dFMatrix(dim,dim)), + shape_values_transform (n_values_arrays, + dFMatrix(n_transform_functions, + n_quadrature_points)), selected_dataset (0), update_flags (update_flags) {}; @@ -154,10 +159,13 @@ FEValues::FEValues (const FiniteElement &fe, FEValuesBase (quadrature.n_quadrature_points, fe.total_dofs, fe.total_dofs, + fe.n_transform_functions, 1, update_flags), unit_shape_gradients(fe.total_dofs, vector >(quadrature.n_quadrature_points)), + unit_shape_gradients_transform(fe.n_transform_functions, + vector >(quadrature.n_quadrature_points)), unit_quadrature_points(quadrature.get_quad_points()) { Assert ((update_flags & update_normal_vectors) == false, @@ -171,6 +179,15 @@ FEValues::FEValues (const FiniteElement &fe, = fe.shape_grad(i, unit_quadrature_points[j]); }; + for (unsigned int i=0; i::reinit (const typename DoFHandler::cell_iterator &cell, update_flags & update_ansatz_points, quadrature_points, update_flags & update_q_points, + shape_values_transform[0], unit_shape_gradients_transform, boundary); // compute gradients on real element if @@ -244,16 +262,21 @@ template FEFaceValuesBase::FEFaceValuesBase (const unsigned int n_q_points, const unsigned int n_ansatz_points, const unsigned int n_dofs, + const unsigned int n_transform_functions, const unsigned int n_faces_or_subfaces, const UpdateFlags update_flags) : FEValuesBase (n_q_points, n_ansatz_points, n_dofs, + n_transform_functions, n_faces_or_subfaces, update_flags), unit_shape_gradients (n_faces_or_subfaces, vector > >(n_dofs, vector >(n_q_points))), + unit_shape_gradients_transform (n_faces_or_subfaces, + vector > >(n_transform_functions, + vector >(n_q_points))), unit_face_quadrature_points (n_q_points, Point()), unit_quadrature_points (n_faces_or_subfaces, vector >(n_q_points, Point())), @@ -287,7 +310,8 @@ FEFaceValues::FEFaceValues (const FiniteElement &fe, FEFaceValuesBase (quadrature.n_quadrature_points, fe.dofs_per_face, fe.total_dofs, - 2*dim, + fe.n_transform_functions, + GeometryInfo::faces_per_cell, update_flags) { unit_face_quadrature_points = quadrature.get_quad_points(); @@ -298,10 +322,10 @@ FEFaceValues::FEFaceValues (const FiniteElement &fe, // of the space with #dim# dimensions. // the points are still on the unit // cell, not on the real cell. - for (unsigned int face=0; face<2*dim; ++face) + for (unsigned int face=0; face::faces_per_cell; ++face) QProjector::project_to_face (quadrature, face, unit_quadrature_points[face]); - for (unsigned int face=0; face<2*dim; ++face) + for (unsigned int face=0; face::faces_per_cell; ++face) for (unsigned int i=0; i::FEFaceValues (const FiniteElement &fe, unit_shape_gradients[face][i][j] = fe.shape_grad(i, unit_quadrature_points[face][j]); }; + + for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int i=0; i::reinit (const typename DoFHandler::cell_iterator &c update_flags & update_JxW_values, normal_vectors, update_flags & update_normal_vectors, + shape_values_transform[face_no], + unit_shape_gradients_transform[face_no], boundary); // compute gradients on real element if @@ -392,7 +428,8 @@ FESubfaceValues::FESubfaceValues (const FiniteElement &fe, FEFaceValuesBase (quadrature.n_quadrature_points, 0, fe.total_dofs, - 2*dim*(1<<(dim-1)), + fe.n_transform_functions, + GeometryInfo::faces_per_cell * GeometryInfo::subfaces_per_face, update_flags) { Assert ((update_flags & update_ansatz_points) == false, @@ -406,21 +443,31 @@ FESubfaceValues::FESubfaceValues (const FiniteElement &fe, // of the space with #dim# dimensions. // the points are still on the unit // cell, not on the real cell. - for (unsigned int face=0; face<2*dim; ++face) - for (unsigned int subface=0; subface<(1<<(dim-1)); ++subface) + for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int subface=0; subface::subfaces_per_face; ++subface) QProjector::project_to_subface (quadrature, face, subface, unit_quadrature_points[face*(1<<(dim-1))+subface]); - for (unsigned int face=0; face<2*dim; ++face) - for (unsigned int subface=0; subface<(1<<(dim-1)); ++subface) + for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int subface=0; subface::subfaces_per_face; ++subface) for (unsigned int i=0; i::subfaces_per_face+subface](i,j) + = fe.shape_value(i, unit_quadrature_points[face*GeometryInfo::subfaces_per_face+subface][j]); + unit_shape_gradients[face*GeometryInfo::subfaces_per_face+subface][i][j] + = fe.shape_grad(i, unit_quadrature_points[face*GeometryInfo::subfaces_per_face+subface][j]); + }; + for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int subface=0; subface::subfaces_per_face; ++subface) + for (unsigned int i=0; i::subfaces_per_face+subface] (i,j) + = fe.shape_value_transform (i, unit_quadrature_points[face*GeometryInfo::subfaces_per_face+subface][j]); + unit_shape_gradients_transform[face*GeometryInfo::subfaces_per_face+subface][i][j] + = fe.shape_grad_transform(i, unit_quadrature_points[face*GeometryInfo::subfaces_per_face+subface][j]); }; }; @@ -459,6 +506,8 @@ void FESubfaceValues::reinit (const typename DoFHandler::cell_iterator update_flags & update_JxW_values, normal_vectors, update_flags & update_normal_vectors, + shape_values_transform[selected_dataset], + unit_shape_gradients_transform[selected_dataset], boundary); // compute gradients on real element if