From 80331fe3482552b8227efc16d8c80efddd71cf09 Mon Sep 17 00:00:00 2001 From: wolf Date: Fri, 6 Nov 1998 19:05:39 +0000 Subject: [PATCH] Convert part of the library from Point and dFMatrix to Tensor where necessary. git-svn-id: https://svn.dealii.org/trunk@652 0785d39b-7218-0410-832d-ea1e28bc413d --- .../Attic/examples/poisson/equation.cc | 2 +- deal.II/deal.II/Todo | 8 + deal.II/deal.II/include/fe/fe.h | 18 +-- .../deal.II/include/fe/fe_lib.criss_cross.h | 8 +- deal.II/deal.II/include/fe/fe_lib.lagrange.h | 24 +-- .../deal.II/include/fe/fe_linear_mapping.h | 2 +- deal.II/deal.II/include/fe/fe_update_flags.h | 8 +- deal.II/deal.II/include/fe/fe_values.h | 76 +++++++++- deal.II/deal.II/source/fe/fe.cc | 12 +- .../deal.II/source/fe/fe_lib.criss_cross.cc | 26 ++-- deal.II/deal.II/source/fe/fe_lib.cubic.cc | 14 +- deal.II/deal.II/source/fe/fe_lib.linear.cc | 14 +- deal.II/deal.II/source/fe/fe_lib.quadratic.cc | 14 +- deal.II/deal.II/source/fe/fe_lib.quartic.cc | 14 +- .../deal.II/source/fe/fe_linear_mapping.cc | 81 +++++----- deal.II/deal.II/source/fe/fe_values.cc | 140 +++++++++++++++--- deal.II/deal.II/source/numerics/matrices.cc | 4 +- tests/big-tests/poisson/equation.cc | 2 +- 18 files changed, 347 insertions(+), 120 deletions(-) diff --git a/deal.II/deal.II/Attic/examples/poisson/equation.cc b/deal.II/deal.II/Attic/examples/poisson/equation.cc index 7844de59d5..e30bfefa76 100644 --- a/deal.II/deal.II/Attic/examples/poisson/equation.cc +++ b/deal.II/deal.II/Attic/examples/poisson/equation.cc @@ -36,7 +36,7 @@ void PoissonEquation<2>::assemble (dFMatrix &cell_matrix, dVector &rhs, const FEValues<2> &fe_values, const Triangulation<2>::cell_iterator &) const { - const vector > >&gradients = fe_values.get_shape_grads (); + const vector > >&gradients = fe_values.get_shape_grads (); const dFMatrix &values = fe_values.get_shape_values (); vector rhs_values (fe_values.n_quadrature_points); const vector &weights = fe_values.get_JxW_values (); diff --git a/deal.II/deal.II/Todo b/deal.II/deal.II/Todo index 53360beabe..5ceb7ec397 100644 --- a/deal.II/deal.II/Todo +++ b/deal.II/deal.II/Todo @@ -175,3 +175,11 @@ Review DoFCellAccessor::get_interpolated_dof_values: it may be Move the exceptions in the Tensor package (tensor_base.h and tensor.h) back into the Tensor class, if the compiler allows to do so. Also rename them back (i.e. drop the initial Tensor* from the name). + + + +Think about the determinant function in the tensor package. Is it + useful, can it be generalized? + + +Implement FELinearMapping::fill_fe_values->jacobians \ No newline at end of file diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index b48770180b..89f9270928 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -550,8 +550,8 @@ class FiniteElement : public FiniteElementBase { * function at the point #p#. * #p# is a point on the reference element, */ - virtual Point shape_grad (const unsigned int i, - const Point &p) const = 0; + virtual Tensor<1,dim> shape_grad (const unsigned int i, + const Point &p) const = 0; /** * Return the tensor of second derivatives @@ -639,7 +639,7 @@ class FiniteElement : public FiniteElementBase { */ virtual void fill_fe_values (const DoFHandler::cell_iterator &cell, const vector > &unit_points, - vector &jacobians, + vector > &jacobians, const bool compute_jacobians, vector > &support_points, const bool compute_support_points, @@ -757,7 +757,7 @@ class FiniteElement : public FiniteElementBase { const unsigned int face_no, const vector > &unit_points, const vector > &global_unit_points, - vector &jacobians, + vector > &jacobians, const bool compute_jacobians, vector > &support_points, const bool compute_support_points, @@ -804,7 +804,7 @@ class FiniteElement : public FiniteElementBase { const unsigned int subface_no, const vector > &unit_points, const vector > &global_unit_points, - vector &jacobians, + vector > &jacobians, const bool compute_jacobians, vector > &q_points, const bool compute_q_points, @@ -863,8 +863,8 @@ class FiniteElement : public FiniteElementBase { * overwritten function has to be provided. */ virtual void get_support_points (const DoFHandler::cell_iterator &cell, - const Boundary &boundary, - vector > &support_points) const; + const Boundary &boundary, + vector > &support_points) const; /** * Compute the off-points of the finite @@ -907,8 +907,8 @@ class FiniteElement : public FiniteElementBase { * when called with #dim==1#. */ virtual void get_face_support_points (const DoFHandler::face_iterator &face, - const Boundary &boundary, - vector > &support_points) const =0; + const Boundary &boundary, + vector > &support_points) const =0; /** * This is the second separated function diff --git a/deal.II/deal.II/include/fe/fe_lib.criss_cross.h b/deal.II/deal.II/include/fe/fe_lib.criss_cross.h index 51da24e2c0..d27931ed96 100644 --- a/deal.II/deal.II/include/fe/fe_lib.criss_cross.h +++ b/deal.II/deal.II/include/fe/fe_lib.criss_cross.h @@ -181,8 +181,8 @@ class FECrissCross : public FiniteElement { * Return the gradient of the #i#th shape * function at point #p# on the unit cell. */ - virtual Point shape_grad(const unsigned int i, - const Point& p) const; + virtual Tensor<1,dim> shape_grad(const unsigned int i, + const Point& p) const; /** * Return the tensor of second derivatives @@ -197,7 +197,7 @@ class FECrissCross : public FiniteElement { * need to do some evaluation by hand. */ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, - const Point &p) const = 0; + const Point &p) const; /** * Refer to the base class for detailed @@ -332,7 +332,7 @@ class FECrissCross : public FiniteElement { */ virtual void fill_fe_values (const DoFHandler::cell_iterator &cell, const vector > &unit_points, - vector &jacobians, + vector > &jacobians, const bool compute_jacobians, vector > &support_points, const bool compute_support_points, 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 cc1bf1ed8c..e443652c2a 100644 --- a/deal.II/deal.II/include/fe/fe_lib.lagrange.h +++ b/deal.II/deal.II/include/fe/fe_lib.lagrange.h @@ -44,8 +44,8 @@ class FELinear : public FELinearMapping { * Return the gradient of the #i#th shape * function at point #p# on the unit cell. */ - virtual Point shape_grad(const unsigned int i, - const Point& p) const; + virtual Tensor<1,dim> shape_grad(const unsigned int i, + const Point& p) const; /** * Return the tensor of second derivatives @@ -56,7 +56,7 @@ class FELinear : public FELinearMapping { * derivatives on the unit cell are zero. */ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, - const Point &p) const = 0; + const Point &p) const; /** * Refer to the base class for detailed @@ -118,8 +118,8 @@ class FEQuadraticSub : public FELinearMapping { * Return the gradient of the #i#th shape * function at point #p# on the unit cell. */ - virtual Point shape_grad(const unsigned int i, - const Point& p) const; + virtual Tensor<1,dim> shape_grad(const unsigned int i, + const Point& p) const; /** * Return the tensor of second derivatives @@ -127,7 +127,7 @@ class FEQuadraticSub : public FELinearMapping { * point #p# on the unit cell. */ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, - const Point &p) const = 0; + const Point &p) const; /** * Refer to the base class for detailed @@ -208,8 +208,8 @@ class FECubicSub : public FELinearMapping { * Return the gradient of the #i#th shape * function at point #p# on the unit cell. */ - virtual Point shape_grad(const unsigned int i, - const Point& p) const; + virtual Tensor<1,dim> shape_grad(const unsigned int i, + const Point& p) const; /** * Return the tensor of second derivatives @@ -217,7 +217,7 @@ class FECubicSub : public FELinearMapping { * point #p# on the unit cell. */ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, - const Point &p) const = 0; + const Point &p) const; /** * Refer to the base class for detailed @@ -299,8 +299,8 @@ class FEQuarticSub : public FELinearMapping { * Return the gradient of the #i#th shape * function at point #p# on the unit cell. */ - virtual Point shape_grad(const unsigned int i, - const Point& p) const; + virtual Tensor<1,dim> shape_grad(const unsigned int i, + const Point& p) const; /** * Return the tensor of second derivatives @@ -308,7 +308,7 @@ class FEQuarticSub : public FELinearMapping { * point #p# on the unit cell. */ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, - const Point &p) const = 0; + const Point &p) const; /** * Refer to the base class for detailed diff --git a/deal.II/deal.II/include/fe/fe_linear_mapping.h b/deal.II/deal.II/include/fe/fe_linear_mapping.h index 09476f8169..735023da01 100644 --- a/deal.II/deal.II/include/fe/fe_linear_mapping.h +++ b/deal.II/deal.II/include/fe/fe_linear_mapping.h @@ -135,7 +135,7 @@ class FELinearMapping : public FiniteElement { */ virtual void fill_fe_values (const DoFHandler::cell_iterator &cell, const vector > &unit_points, - vector &jacobians, + vector > &jacobians, const bool compute_jacobians, vector > &support_points, const bool compute_support_points, diff --git a/deal.II/deal.II/include/fe/fe_update_flags.h b/deal.II/deal.II/include/fe/fe_update_flags.h index 46d729bb59..9e61b96271 100644 --- a/deal.II/deal.II/include/fe/fe_update_flags.h +++ b/deal.II/deal.II/include/fe/fe_update_flags.h @@ -65,7 +65,13 @@ enum UpdateFlags { * an error, since normal vectors are * not useful in that case. */ - update_normal_vectors = 32 + update_normal_vectors = 32, + + /** + * Update the second derivatives of the + * shape functions on the real cell. + */ + update_second_derivatives }; diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index bf0768e08e..12fa99453f 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -118,6 +118,22 @@ template class Quadrature; * #FELinearMapping::fill_fe_values# function, where also a small test * program is presented. * + * The derivatives of the Jacobi matrices at the quadrature points with respect + * to unit cell coordinates is stored in the field names + * #jacobi_matrices_grad#. Since the gradient of a shape function is given by + * $\partial_i \phi = \sum_k \hat\partial_k \hat\phi J_{ki}$, where $\hat$ + * denotes differentiation on the unit cell, the second derivative of a + * function is given by + * $\partial_j \partial i \phi + * = + * \hat\partial_l [ (\hat \partial_k \hat\phi) J_{ki} ] J_{lj} + * = + * (\hat\partial_k \hat\partial_l \hat\phi) J_{ki} J_{lj} + * + + * (\hat \partial_l \hat\phi) (\hat\partial_l J_{ki}) J_{lj}$. + * While we already have access to the Jacobian matrix, the derivatives are + * stored in the named field. + * * * \subsection{Member functions} * @@ -289,8 +305,8 @@ class FEValuesBase { * The function returns the gradient on the * real element, not the reference element. */ - const Point & shape_grad (const unsigned int function, - const unsigned int quadrature_point) const; + const Tensor<1,dim> & shape_grad (const unsigned int function, + const unsigned int quadrature_point) const; /** * Return a pointer to the matrix holding @@ -299,7 +315,7 @@ class FEValuesBase { * For the format of this matrix, see the * documentation for the matrix itself. */ - const vector > > & get_shape_grads () const; + const vector > > & get_shape_grads () const; /** * Return the gradients of the finite @@ -458,7 +474,18 @@ class FEValuesBase { * gradients on the real element, rather * than on the reference element. */ - vector > > shape_gradients; + vector > > shape_gradients; + + /** + * Store the 2nd derivatives of the shape + * functions at the quadrature points. + * + * This field is reset each time + * #reinit# is called and contains the + * gradients on the real element, rather + * than on the reference element. + */ + vector > > shape_2nd_derivatives; /** * Store an array of the weights of the @@ -516,8 +543,21 @@ class FEValuesBase { * transformation of the gradients to the * real cell. */ - vector jacobi_matrices; + vector > jacobi_matrices; + /** + * Store the derivatives of the jacobi + * matrices. If #J[j][k]# is the jacobi + * matrix, then the index #[i][j][k]# + * of this field denotes the derivation + * of #J[j][k]# with respect to the + * #i#th variable. + * + * The same general remarks apply as for + * #jacobi_matrices#. + */ + vector > jacobi_matrices_grad; + /** * Store the values of the basis functions * of the transformation from unit cell @@ -646,7 +686,17 @@ class FEValues : public FEValuesBase { * of the object and contains the gradients * on the reference element. */ - vector > > unit_shape_gradients; + vector > > unit_shape_gradients; + + /** + * Store the 2nd derivatives of the shape + * functions at the quadrature points on + * the unit cell. + * This field is set up upon construction + * of the object and contains the + * derivatives on the reference element. + */ + vector > > unit_shape_2nd_derivatives; /** * Gradients of the basis @@ -793,8 +843,18 @@ class FEFaceValuesBase : public FEValuesBase { * subface, with indices like that: * #unit_shape_gradients[face][dof][q_point]# */ - vector > > > unit_shape_gradients; + vector > > > unit_shape_gradients; + /** + * Store the 2nd derivatives of the shape + * functions at the quadrature points on + * the unit cell for each face. + * This field is set up upon construction + * of the object and contains the + * derivatives on the reference element. + */ + vector > > > unit_shape_2nd_derivatives; + /** * Gradients of the basis * functions of the transformation. @@ -1087,7 +1147,7 @@ const dFMatrix & FEValuesBase::get_shape_values () const { template inline -const vector > > & +const vector > > & FEValuesBase::get_shape_grads () const { Assert (update_flags & update_gradients, ExcAccessToUninitializedField()); return shape_gradients; diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index f3379ebea6..0833482a4a 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -178,7 +178,7 @@ bool FiniteElementBase::operator == (const FiniteElementBase &f) const template <> void FiniteElement<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, const vector > &unit_points, - vector &jacobians, + vector > &jacobians, const bool compute_jacobians, vector > &support_points, const bool compute_support_points, @@ -224,7 +224,7 @@ void FiniteElement<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator & const unsigned int , const vector > &, const vector > &, - vector &, + vector > &, const bool , vector > &, const bool , @@ -247,7 +247,7 @@ void FiniteElement<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterato const unsigned int , const vector > &, const vector > &, - vector &, + vector > &, const bool , vector > &, const bool , @@ -316,7 +316,7 @@ void FiniteElement<1>::get_support_points (const DoFHandler<1>::cell_iterator &c template void FiniteElement::fill_fe_values (const DoFHandler::cell_iterator &, const vector > &, - vector &, + vector > &, const bool, vector > &, const bool, @@ -335,7 +335,7 @@ void FiniteElement::fill_fe_face_values (const DoFHandler::cell_iterat const unsigned int face_no, const vector > &unit_points, const vector > &global_unit_points, - vector &jacobians, + vector > &jacobians, const bool compute_jacobians, vector > &support_points, const bool compute_support_points, @@ -386,7 +386,7 @@ void FiniteElement::fill_fe_subface_values (const DoFHandler::cell_ite const unsigned int subface_no, const vector > &unit_points, const vector > &global_unit_points, - vector &jacobians, + vector > &jacobians, const bool compute_jacobians, vector > &q_points, const bool compute_q_points, diff --git a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc index cb46024d83..5cace7e1d7 100644 --- a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc +++ b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc @@ -323,7 +323,7 @@ double FECrissCross<1>::shape_value (const unsigned int, const Point<1> &) const template <> -Point<1> FECrissCross<1>::shape_grad (const unsigned int, const Point<1> &) const { +Tensor<1,1> FECrissCross<1>::shape_grad (const unsigned int, const Point<1> &) const { Assert (false, ExcNotUseful()); return Point<1>(); }; @@ -538,7 +538,7 @@ double FECrissCross<2>::shape_value (const unsigned int i, template <> inline -Point<2> FECrissCross<2>::shape_grad (const unsigned int i, const Point<2> &p) const { +Tensor<1,2> FECrissCross<2>::shape_grad (const unsigned int i, const Point<2> &p) const { Assert((i::get_normal_vectors (const DoFHandler<2>::cell_iterator &ce template void FECrissCross::fill_fe_values (const DoFHandler::cell_iterator &cell, const vector > &unit_points, - vector &jacobians, + vector > &jacobians, const bool compute_jacobians, vector > &support_points, const bool, @@ -1023,16 +1023,16 @@ void FECrissCross::fill_fe_values (const DoFHandler::cell_iterator &ce vertices[0](1)+t28-t22*xi; const double t38 = 1/(t24+t36); - jacobians[point](0,0) = (-vertices[0](1)+vertices[0](1)*xi- - vertices[1](1)*xi+vertices[2](1)*xi+ - vertices[3](1)-vertices[3](1)*xi)*t38; - jacobians[point](0,1) = -(-vertices[0](0)+vertices[0](0)*xi- - vertices[1](0)*xi+t8+vertices[3](0)-t16)*t38; - jacobians[point](1,0) = -(-vertices[0](1)+vertices[0](1)*eta+ - vertices[1](1)-vertices[1](1)*eta+ - vertices[2](1)*eta-vertices[3](1)*eta)*t38; - jacobians[point](1,1) = (-vertices[0](0)+vertices[0](0)*eta+ - vertices[1](0)-t10+t31-vertices[3](0)*eta)*t38; + jacobians[point][0][0] = (-vertices[0](1)+vertices[0](1)*xi- + vertices[1](1)*xi+vertices[2](1)*xi+ + vertices[3](1)-vertices[3](1)*xi)*t38; + jacobians[point][0][1] = -(-vertices[0](0)+vertices[0](0)*xi- + vertices[1](0)*xi+t8+vertices[3](0)-t16)*t38; + jacobians[point][1][0] = -(-vertices[0](1)+vertices[0](1)*eta+ + vertices[1](1)-vertices[1](1)*eta+ + vertices[2](1)*eta-vertices[3](1)*eta)*t38; + jacobians[point][1][1] = (-vertices[0](0)+vertices[0](0)*eta+ + vertices[1](0)-t10+t31-vertices[3](0)*eta)*t38; }; }; 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 1216637d63..9c62a710d8 100644 --- a/deal.II/deal.II/source/fe/fe_lib.cubic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.cubic.cc @@ -107,12 +107,17 @@ FECubicSub<1>::shape_value (const unsigned int i, template <> -Point<1> +Tensor<1,1> FECubicSub<1>::shape_grad (const unsigned int i, const Point<1> &p) const { Assert((i, so we + // still construct it as that. it should + // make no difference in practice, + // however switch (i) { case 0: return Point<1>(-27.0/2.0*xi*xi+18.0*xi-11.0/2.0); @@ -721,7 +726,7 @@ xi*xi)*eta*eta*eta; template <> -Point<2> +Tensor<1,2> FECubicSub<2>::shape_grad (const unsigned int i, const Point<2> &p) const { @@ -729,6 +734,11 @@ FECubicSub<2>::shape_grad (const unsigned int i, const double xi = p(0), eta= p(1); + // originally, the return type of the + // function was Point, so we + // still construct it as that. it should + // make no difference in practice, + // however switch (i) { case 0: return Point<2>(-11.0/2.0+18.0*xi-27.0/2.0*xi*xi+(121.0/4.0-99.0*xi+297.0/4.0*xi*xi)*eta+(-99.0/2.0+162.0*xi-243.0/2.0*xi*xi)*eta*eta+(99.0/4.0-81.0*xi+243.0/4.0*xi*xi)*eta*eta*eta, diff --git a/deal.II/deal.II/source/fe/fe_lib.linear.cc b/deal.II/deal.II/source/fe/fe_lib.linear.cc index 33babbfa5c..08d1580d4d 100644 --- a/deal.II/deal.II/source/fe/fe_lib.linear.cc +++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc @@ -62,11 +62,16 @@ FELinear<1>::shape_value(const unsigned int i, template <> inline -Point<1> +Tensor<1,1> FELinear<1>::shape_grad(const unsigned int i, const Point<1>&) const { Assert((i, so we + // still construct it as that. it should + // make no difference in practice, + // however switch (i) { case 0: return Point<1>(-1.); @@ -215,11 +220,16 @@ FELinear<2>::shape_value (const unsigned int i, template <> inline -Point<2> +Tensor<1,2> FELinear<2>::shape_grad (const unsigned int i, const Point<2>& p) const { Assert((i, so we + // still construct it as that. it should + // make no difference in practice, + // however switch (i) { case 0: return Point<2> (p(1)-1., p(0)-1.); 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 49236de6d6..33c3a4c7e3 100644 --- a/deal.II/deal.II/source/fe/fe_lib.quadratic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.quadratic.cc @@ -118,12 +118,17 @@ FEQuadraticSub<1>::shape_value(const unsigned int i, template <> -Point<1> +Tensor<1,1> FEQuadraticSub<1>::shape_grad(const unsigned int i, const Point<1> &p) const { Assert((i, so we + // still construct it as that. it should + // make no difference in practice, + // however switch (i) { case 0: return Point<1>(-3+4*xi); @@ -632,7 +637,7 @@ FEQuadraticSub<2>::shape_value (const unsigned int i, template <> -Point<2> +Tensor<1,2> FEQuadraticSub<2>::shape_grad (const unsigned int i, const Point<2> &p) const { @@ -640,6 +645,11 @@ FEQuadraticSub<2>::shape_grad (const unsigned int i, const double xi = p(0), eta= p(1); + // originally, the return type of the + // function was Point, so we + // still construct it as that. it should + // make no difference in practice, + // however switch (i) { case 0: return Point<2>(-(2*xi-1)*(1-eta)*(2*eta-1)+2*(1-xi)*(1-eta)*(2*eta-1), 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 7cd303f6e1..a6cf7ffd35 100644 --- a/deal.II/deal.II/source/fe/fe_lib.quartic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.quartic.cc @@ -75,12 +75,17 @@ FEQuarticSub<1>::shape_value(const unsigned int i, template <> -Point<1> +Tensor<1,1> FEQuarticSub<1>::shape_grad(const unsigned int i, const Point<1> &p) const { Assert((i, so we + // still construct it as that. it should + // make no difference in practice, + // however switch (i) { case 0: return Point<1>(128.0/3.0*xi*xi*xi-80.0*xi*xi+140.0/3.0*xi-25.0/3.0); @@ -1089,7 +1094,7 @@ FEQuarticSub<2>::shape_value (const unsigned int i, template <> -Point<2> +Tensor<1,2> FEQuarticSub<2>::shape_grad (const unsigned int i, const Point<2> &p) const { @@ -1097,6 +1102,11 @@ FEQuarticSub<2>::shape_grad (const unsigned int i, const double xi = p(0), eta= p(1); + // originally, the return type of the + // function was Point, so we + // still construct it as that. it should + // make no difference in practice, + // however switch (i) { case 0: return Point<2>(-25.0/3.0+140.0/3.0*xi-80.0*xi*xi+128.0/3.0*xi*xi*xi+(625.0/9.0-3500.0/9.0*xi+2000.0/3.0*xi*xi-3200.0/9.0*xi*xi*xi)*eta+(-1750.0/9.0+9800.0/9.0*xi-5600.0/3.0*xi*xi+8960.0/9.0*xi*xi*xi)*eta*eta+(2000.0/9.0-11200.0/9.0*xi+6400.0/3.0*xi*xi-10240.0/9.0*xi*xi*xi)*eta*eta*eta+(-800.0/9.0+4480.0/9.0*xi-2560.0/3.0*xi*xi+4096.0/9.0*xi*xi*xi)*eta*eta*eta*eta, diff --git a/deal.II/deal.II/source/fe/fe_linear_mapping.cc b/deal.II/deal.II/source/fe/fe_linear_mapping.cc index a6eada6302..e2728ad763 100644 --- a/deal.II/deal.II/source/fe/fe_linear_mapping.cc +++ b/deal.II/deal.II/source/fe/fe_linear_mapping.cc @@ -279,7 +279,7 @@ void FELinearMapping<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator template void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator &cell, const vector > &unit_points, - vector &jacobians, + vector > &jacobians, const bool compute_jacobians, vector > &support_points, const bool compute_support_points, @@ -429,40 +429,53 @@ void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator if (compute_jacobians) for (unsigned int point=0; point::FEValuesBase (const unsigned int 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)), + shape_gradients (n_dofs, vector >(n_q_points)), weights (n_q_points, 0), JxW_values (n_q_points, 0), quadrature_points (n_q_points, Point()), support_points (n_support_points, Point()), - jacobi_matrices (n_q_points, dFMatrix(dim,dim)), + jacobi_matrices (n_q_points, Tensor<2,dim>()), + jacobi_matrices_grad (n_q_points, Tensor<3,dim>()), shape_values_transform (n_values_arrays, dFMatrix(n_transform_functions, n_quadrature_points)), @@ -85,7 +86,7 @@ void FEValuesBase::get_function_values (const dVector &fe_function, template -const Point & +const Tensor<1,dim> & FEValuesBase::shape_grad (const unsigned int i, const unsigned int j) const { Assert (i::get_function_grads (const dVector &fe_function, // functions for (unsigned int point=0; point tmp(shape_gradients[shape_func][point]); + tmp *= dof_values(shape_func); + gradients[point] += tmp; + }; }; @@ -168,7 +172,9 @@ FEValues::FEValues (const FiniteElement &fe, update_flags, fe), unit_shape_gradients(fe.total_dofs, - vector >(quadrature.n_quadrature_points)), + vector >(quadrature.n_quadrature_points)), + unit_shape_2nd_derivatives(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()) @@ -182,6 +188,8 @@ FEValues::FEValues (const FiniteElement &fe, shape_values[0](i,j) = fe.shape_value(i, unit_quadrature_points[j]); unit_shape_gradients[i][j] = fe.shape_grad(i, unit_quadrature_points[j]); + unit_shape_2nd_derivatives[i][j] + = fe.shape_grad_grad(i, unit_quadrature_points[j]); }; for (unsigned int i=0; i::reinit (const typename DoFHandler::cell_iterator &cell, jacobi_matrices, update_flags & (update_jacobians | update_JxW_values | - update_gradients), + update_gradients | + update_second_derivatives), +// jacobi_matrices_grad, +// update_flags & update_second_derivatives, support_points, update_flags & update_support_points, quadrature_points, @@ -233,16 +244,41 @@ void FEValues::reinit (const typename DoFHandler::cell_iterator &cell, for (unsigned int j=0; jtotal_dofs; ++i) + for (unsigned int j=0; j tmp1, tmp2; + + // tmp1 := (d_k d_l phi) J_lj + contract (tmp1, unit_shape_2nd_derivatives[i][j], jacobi_matrices[j]); + // tmp2 := tmp1_kj J_ki + contract (tmp2, tmp1, 1, jacobi_matrices[j], 1); + + + // second part: + // tmp1 := (d_k J_lj) (d_l phi) + contract (tmp1, jacobi_matrices_grad[j], 2, unit_shape_gradients[i][j]); + // tmp1_kj J_ki + contract (shape_2nd_derivatives[i][j], + jacobi_matrices[j], 1, + tmp1, 1); + + // add up first contribution + shape_2nd_derivatives[i][j] += tmp2; + }; // compute Jacobi determinants in @@ -252,7 +288,7 @@ void FEValues::reinit (const typename DoFHandler::cell_iterator &cell, // determinant if (update_flags & update_JxW_values) for (unsigned int i=0; i::FEFaceValuesBase (const unsigned int n_q_points, update_flags, fe), unit_shape_gradients (n_faces_or_subfaces, - vector > >(n_dofs, - vector >(n_q_points))), + vector > >(n_dofs, + vector >(n_q_points))), + unit_shape_2nd_derivatives(n_dofs, + vector > >(n_dofs, + vector >(n_q_points))), unit_shape_gradients_transform (n_faces_or_subfaces, vector > >(n_transform_functions, vector >(n_q_points))), @@ -340,6 +379,8 @@ FEFaceValues::FEFaceValues (const FiniteElement &fe, = fe.shape_value(i, unit_quadrature_points[face][j]); unit_shape_gradients[face][i][j] = fe.shape_grad(i, unit_quadrature_points[face][j]); + unit_shape_2nd_derivatives[face][i][j] + = fe.shape_grad_grad(i, unit_quadrature_points[face][j]); }; for (unsigned int i=0; i::reinit (const typename DoFHandler::cell_iterator &c // (grad_{\xi\eta})_b J_{bs} // with J_{bs}=(d\xi_b)/(dx_s) for (unsigned int b=0; btotal_dofs; ++i) + for (unsigned int j=0; j tmp1, tmp2; + + // tmp1 := (d_k d_l phi) J_lj + contract (tmp1, unit_shape_2nd_derivatives[face_no][i][j], jacobi_matrices[j]); + // tmp2 := tmp1_kj J_ki + contract (tmp2, tmp1, 1, jacobi_matrices[j], 1); + + + // second part: + // tmp1 := (d_k J_lj) (d_l phi) + contract (tmp1, + jacobi_matrices_grad[j], 2, + unit_shape_gradients[face_no][i][j]); + // tmp1_kj J_ki + contract (shape_2nd_derivatives[i][j], + jacobi_matrices[j], 1, + tmp1, 1); + + // add up first contribution + shape_2nd_derivatives[i][j] += tmp2; + }; // compute Jacobi determinants in @@ -424,7 +492,7 @@ void FEFaceValues::reinit (const typename DoFHandler::cell_iterator &c -/*------------------------------- FEFaceValues -------------------------------*/ +/*------------------------------- FESubFaceValues -------------------------------*/ template @@ -469,6 +537,10 @@ FESubfaceValues::FESubfaceValues (const FiniteElement &fe, = fe.shape_grad(i, unit_quadrature_points[face * GeometryInfo:: subfaces_per_face+subface][j]); + unit_shape_2nd_derivatives[face*GeometryInfo::subfaces_per_face+subface][i][j] + = fe.shape_grad_grad(i, unit_quadrature_points[face * + GeometryInfo:: + subfaces_per_face+subface][j]); }; for (unsigned int i=0; i::reinit (const typename DoFHandler::cell_iterator // (grad_{\xi\eta})_b J_{bs} // with J_{bs}=(d\xi_b)/(dx_s) for (unsigned int b=0; btotal_dofs; ++i) + for (unsigned int j=0; j tmp1, tmp2; + + // tmp1 := (d_k d_l phi) J_lj + contract (tmp1, + unit_shape_2nd_derivatives[selected_dataset][i][j], + jacobi_matrices[j]); + // tmp2 := tmp1_kj J_ki + contract (tmp2, tmp1, 1, jacobi_matrices[j], 1); + + + // second part: + // tmp1 := (d_k J_lj) (d_l phi) + contract (tmp1, + jacobi_matrices_grad[j], 2, + unit_shape_gradients[selected_dataset][i][j]); + // tmp1_kj J_ki + contract (shape_2nd_derivatives[i][j], + jacobi_matrices[j], 1, + tmp1, 1); + + // add up first contribution + shape_2nd_derivatives[i][j] += tmp2; + }; // compute Jacobi determinants in diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index 19e5d6cc7b..9075fc1e9f 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -694,7 +694,7 @@ void LaplaceMatrix::assemble (dFMatrix &cell_matrix, Assert (cell_matrix.all_zero(), ExcObjectNotEmpty()); Assert (rhs.all_zero(), ExcObjectNotEmpty()); - const vector > >&gradients = fe_values.get_shape_grads (); + const vector > >&gradients = fe_values.get_shape_grads (); const dFMatrix &values = fe_values.get_shape_values (); vector rhs_values(fe_values.n_quadrature_points); const vector &weights = fe_values.get_JxW_values (); @@ -748,7 +748,7 @@ void LaplaceMatrix::assemble (dFMatrix &cell_matrix, ExcWrongSize(cell_matrix.m(), total_dofs)); Assert (cell_matrix.all_zero(), ExcObjectNotEmpty()); - const vector > >&gradients = fe_values.get_shape_grads (); + const vector > >&gradients = fe_values.get_shape_grads (); const vector &weights = fe_values.get_JxW_values (); if (coefficient != 0) diff --git a/tests/big-tests/poisson/equation.cc b/tests/big-tests/poisson/equation.cc index 7844de59d5..e30bfefa76 100644 --- a/tests/big-tests/poisson/equation.cc +++ b/tests/big-tests/poisson/equation.cc @@ -36,7 +36,7 @@ void PoissonEquation<2>::assemble (dFMatrix &cell_matrix, dVector &rhs, const FEValues<2> &fe_values, const Triangulation<2>::cell_iterator &) const { - const vector > >&gradients = fe_values.get_shape_grads (); + const vector > >&gradients = fe_values.get_shape_grads (); const dFMatrix &values = fe_values.get_shape_values (); vector rhs_values (fe_values.n_quadrature_points); const vector &weights = fe_values.get_JxW_values (); -- 2.39.5