dVector &rhs,
const FEValues<2> &fe_values,
const Triangulation<2>::cell_iterator &) const {
- const vector<vector<Point<2> > >&gradients = fe_values.get_shape_grads ();
+ const vector<vector<Tensor<1,2> > >&gradients = fe_values.get_shape_grads ();
const dFMatrix &values = fe_values.get_shape_values ();
vector<double> rhs_values (fe_values.n_quadrature_points);
const vector<double> &weights = fe_values.get_JxW_values ();
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
* function at the point #p#.
* #p# is a point on the reference element,
*/
- virtual Point<dim> shape_grad (const unsigned int i,
- const Point<dim> &p) const = 0;
+ virtual Tensor<1,dim> shape_grad (const unsigned int i,
+ const Point<dim> &p) const = 0;
/**
* Return the tensor of second derivatives
*/
virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
- vector<dFMatrix> &jacobians,
+ vector<Tensor<2,dim> > &jacobians,
const bool compute_jacobians,
vector<Point<dim> > &support_points,
const bool compute_support_points,
const unsigned int face_no,
const vector<Point<dim-1> > &unit_points,
const vector<Point<dim> > &global_unit_points,
- vector<dFMatrix> &jacobians,
+ vector<Tensor<2,dim> > &jacobians,
const bool compute_jacobians,
vector<Point<dim> > &support_points,
const bool compute_support_points,
const unsigned int subface_no,
const vector<Point<dim-1> > &unit_points,
const vector<Point<dim> > &global_unit_points,
- vector<dFMatrix> &jacobians,
+ vector<Tensor<2,dim> > &jacobians,
const bool compute_jacobians,
vector<Point<dim> > &q_points,
const bool compute_q_points,
* overwritten function has to be provided.
*/
virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
- const Boundary<dim> &boundary,
- vector<Point<dim> > &support_points) const;
+ const Boundary<dim> &boundary,
+ vector<Point<dim> > &support_points) const;
/**
* Compute the off-points of the finite
* when called with #dim==1#.
*/
virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
- const Boundary<dim> &boundary,
- vector<Point<dim> > &support_points) const =0;
+ const Boundary<dim> &boundary,
+ vector<Point<dim> > &support_points) const =0;
/**
* This is the second separated function
* Return the gradient of the #i#th shape
* function at point #p# on the unit cell.
*/
- virtual Point<dim> shape_grad(const unsigned int i,
- const Point<dim>& p) const;
+ virtual Tensor<1,dim> shape_grad(const unsigned int i,
+ const Point<dim>& p) const;
/**
* Return the tensor of second derivatives
* need to do some evaluation by hand.
*/
virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
- const Point<dim> &p) const = 0;
+ const Point<dim> &p) const;
/**
* Refer to the base class for detailed
*/
virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
- vector<dFMatrix> &jacobians,
+ vector<Tensor<2,dim> > &jacobians,
const bool compute_jacobians,
vector<Point<dim> > &support_points,
const bool compute_support_points,
* Return the gradient of the #i#th shape
* function at point #p# on the unit cell.
*/
- virtual Point<dim> shape_grad(const unsigned int i,
- const Point<dim>& p) const;
+ virtual Tensor<1,dim> shape_grad(const unsigned int i,
+ const Point<dim>& p) const;
/**
* Return the tensor of second derivatives
* derivatives on the unit cell are zero.
*/
virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
- const Point<dim> &p) const = 0;
+ const Point<dim> &p) const;
/**
* Refer to the base class for detailed
* Return the gradient of the #i#th shape
* function at point #p# on the unit cell.
*/
- virtual Point<dim> shape_grad(const unsigned int i,
- const Point<dim>& p) const;
+ virtual Tensor<1,dim> shape_grad(const unsigned int i,
+ const Point<dim>& p) const;
/**
* Return the tensor of second derivatives
* point #p# on the unit cell.
*/
virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
- const Point<dim> &p) const = 0;
+ const Point<dim> &p) const;
/**
* Refer to the base class for detailed
* Return the gradient of the #i#th shape
* function at point #p# on the unit cell.
*/
- virtual Point<dim> shape_grad(const unsigned int i,
- const Point<dim>& p) const;
+ virtual Tensor<1,dim> shape_grad(const unsigned int i,
+ const Point<dim>& p) const;
/**
* Return the tensor of second derivatives
* point #p# on the unit cell.
*/
virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
- const Point<dim> &p) const = 0;
+ const Point<dim> &p) const;
/**
* Refer to the base class for detailed
* Return the gradient of the #i#th shape
* function at point #p# on the unit cell.
*/
- virtual Point<dim> shape_grad(const unsigned int i,
- const Point<dim>& p) const;
+ virtual Tensor<1,dim> shape_grad(const unsigned int i,
+ const Point<dim>& p) const;
/**
* Return the tensor of second derivatives
* point #p# on the unit cell.
*/
virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
- const Point<dim> &p) const = 0;
+ const Point<dim> &p) const;
/**
* Refer to the base class for detailed
*/
virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
- vector<dFMatrix> &jacobians,
+ vector<Tensor<2,dim> > &jacobians,
const bool compute_jacobians,
vector<Point<dim> > &support_points,
const bool compute_support_points,
* 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
};
* #FELinearMapping<dim>::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}
*
* The function returns the gradient on the
* real element, not the reference element.
*/
- const Point<dim> & 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
* For the format of this matrix, see the
* documentation for the matrix itself.
*/
- const vector<vector<Point<dim> > > & get_shape_grads () const;
+ const vector<vector<Tensor<1,dim> > > & get_shape_grads () const;
/**
* Return the gradients of the finite
* gradients on the real element, rather
* than on the reference element.
*/
- vector<vector<Point<dim> > > shape_gradients;
+ vector<vector<Tensor<1,dim> > > 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<vector<Tensor<2,dim> > > shape_2nd_derivatives;
/**
* Store an array of the weights of the
* transformation of the gradients to the
* real cell.
*/
- vector<dFMatrix> jacobi_matrices;
+ vector<Tensor<2,dim> > 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<Tensor<3,dim> > jacobi_matrices_grad;
+
/**
* Store the values of the basis functions
* of the transformation from unit cell
* of the object and contains the gradients
* on the reference element.
*/
- vector<vector<Point<dim> > > unit_shape_gradients;
+ vector<vector<Tensor<1,dim> > > 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<vector<Tensor<2,dim> > > unit_shape_2nd_derivatives;
/**
* Gradients of the basis
* subface, with indices like that:
* #unit_shape_gradients[face][dof][q_point]#
*/
- vector<vector<vector<Point<dim> > > > unit_shape_gradients;
+ vector<vector<vector<Tensor<1,dim> > > > 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<vector<vector<Tensor<2,dim> > > > unit_shape_2nd_derivatives;
+
/**
* Gradients of the basis
* functions of the transformation.
template <int dim>
inline
-const vector<vector<Point<dim> > > &
+const vector<vector<Tensor<1,dim> > > &
FEValuesBase<dim>::get_shape_grads () const {
Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
return shape_gradients;
template <>
void FiniteElement<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
const vector<Point<1> > &unit_points,
- vector<dFMatrix> &jacobians,
+ vector<Tensor<2,1> > &jacobians,
const bool compute_jacobians,
vector<Point<1> > &support_points,
const bool compute_support_points,
const unsigned int ,
const vector<Point<0> > &,
const vector<Point<1> > &,
- vector<dFMatrix> &,
+ vector<Tensor<2,1> > &,
const bool ,
vector<Point<1> > &,
const bool ,
const unsigned int ,
const vector<Point<0> > &,
const vector<Point<1> > &,
- vector<dFMatrix> &,
+ vector<Tensor<2,1> > &,
const bool ,
vector<Point<1> > &,
const bool ,
template <int dim>
void FiniteElement<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &,
const vector<Point<dim> > &,
- vector<dFMatrix> &,
+ vector<Tensor<2,dim> > &,
const bool,
vector<Point<dim> > &,
const bool,
const unsigned int face_no,
const vector<Point<dim-1> > &unit_points,
const vector<Point<dim> > &global_unit_points,
- vector<dFMatrix> &jacobians,
+ vector<Tensor<2,dim> > &jacobians,
const bool compute_jacobians,
vector<Point<dim> > &support_points,
const bool compute_support_points,
const unsigned int subface_no,
const vector<Point<dim-1> > &unit_points,
const vector<Point<dim> > &global_unit_points,
- vector<dFMatrix> &jacobians,
+ vector<Tensor<2,dim> > &jacobians,
const bool compute_jacobians,
vector<Point<dim> > &q_points,
const bool compute_q_points,
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>();
};
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<total_dofs), ExcInvalidIndex(i));
const double x = p(0),
template <int dim>
void FECrissCross<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
- vector<dFMatrix> &jacobians,
+ vector<Tensor<2,dim> > &jacobians,
const bool compute_jacobians,
vector<Point<dim> > &support_points,
const bool,
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;
};
};
template <>
-Point<1>
+Tensor<1,1>
FECubicSub<1>::shape_grad (const unsigned int i,
const Point<1> &p) const
{
Assert((i<total_dofs), ExcInvalidIndex(i));
const double xi = p(0);
+ // originally, the return type of the
+ // function was Point<dim>, 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);
template <>
-Point<2>
+Tensor<1,2>
FECubicSub<2>::shape_grad (const unsigned int i,
const Point<2> &p) const
{
const double xi = p(0),
eta= p(1);
+ // originally, the return type of the
+ // function was Point<dim>, 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,
template <>
inline
-Point<1>
+Tensor<1,1>
FELinear<1>::shape_grad(const unsigned int i,
const Point<1>&) const
{
Assert((i<total_dofs), ExcInvalidIndex(i));
+ // originally, the return type of the
+ // function was Point<dim>, so we
+ // still construct it as that. it should
+ // make no difference in practice,
+ // however
switch (i)
{
case 0: return Point<1>(-1.);
template <>
inline
-Point<2>
+Tensor<1,2>
FELinear<2>::shape_grad (const unsigned int i,
const Point<2>& p) const
{
Assert((i<total_dofs), ExcInvalidIndex(i));
+ // originally, the return type of the
+ // function was Point<dim>, 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.);
template <>
-Point<1>
+Tensor<1,1>
FEQuadraticSub<1>::shape_grad(const unsigned int i,
const Point<1> &p) const
{
Assert((i<total_dofs), ExcInvalidIndex(i));
const double xi = p(0);
+ // originally, the return type of the
+ // function was Point<dim>, 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);
template <>
-Point<2>
+Tensor<1,2>
FEQuadraticSub<2>::shape_grad (const unsigned int i,
const Point<2> &p) const
{
const double xi = p(0),
eta= p(1);
+ // originally, the return type of the
+ // function was Point<dim>, 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),
template <>
-Point<1>
+Tensor<1,1>
FEQuarticSub<1>::shape_grad(const unsigned int i,
const Point<1> &p) const
{
Assert((i<total_dofs), ExcInvalidIndex(i));
const double xi = p(0);
+ // originally, the return type of the
+ // function was Point<dim>, 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);
template <>
-Point<2>
+Tensor<1,2>
FEQuarticSub<2>::shape_grad (const unsigned int i,
const Point<2> &p) const
{
const double xi = p(0),
eta= p(1);
+ // originally, the return type of the
+ // function was Point<dim>, 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,
template <int dim>
void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
- vector<dFMatrix> &jacobians,
+ vector<Tensor<2,dim> > &jacobians,
const bool compute_jacobians,
vector<Point<dim> > &support_points,
const bool compute_support_points,
if (compute_jacobians)
for (unsigned int point=0; point<n_points; ++point)
- {
- const double xi = unit_points[point](0);
- const double eta= unit_points[point](1);
+ switch (dim)
+ {
+ case 2:
+ {
+
+ const double xi = unit_points[point](0);
+ const double eta= unit_points[point](1);
- const double t6 = vertices[0](0)*vertices[3](1);
- const double t8 = vertices[2](0)*xi;
- const double t10 = vertices[1](0)*eta;
- const double t12 = vertices[3](0)*vertices[1](1);
- const double t16 = vertices[3](0)*xi;
- const double t20 = vertices[0](0)*vertices[1](1);
- const double t22 = vertices[0](0)*vertices[2](1);
- const double t24 = t6*xi-t8*vertices[1](1)-t10*vertices[3](1)+
- t12*eta-vertices[3](0)*vertices[2](1)*eta-
- t16*vertices[0](1)+t16*vertices[1](1)-t12+
- vertices[3](0)*vertices[0](1)-t20*eta+t22*eta;
- const double t28 = vertices[1](0)*vertices[3](1);
- const double t31 = vertices[2](0)*eta;
- const double t36 = t8*vertices[0](1)+vertices[1](0)*vertices[2](1)*xi-
- t28*xi+t10*vertices[0](1)-t31*vertices[0](1)+
- t31*vertices[3](1)+t20-t6-vertices[1](0)*
- 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;
- };
+ const double t6 = vertices[0](0)*vertices[3](1);
+ const double t8 = vertices[2](0)*xi;
+ const double t10 = vertices[1](0)*eta;
+ const double t12 = vertices[3](0)*vertices[1](1);
+ const double t16 = vertices[3](0)*xi;
+ const double t20 = vertices[0](0)*vertices[1](1);
+ const double t22 = vertices[0](0)*vertices[2](1);
+ const double t24 = t6*xi-t8*vertices[1](1)-t10*vertices[3](1)+
+ t12*eta-vertices[3](0)*vertices[2](1)*eta-
+ t16*vertices[0](1)+t16*vertices[1](1)-t12+
+ vertices[3](0)*vertices[0](1)-t20*eta+t22*eta;
+ const double t28 = vertices[1](0)*vertices[3](1);
+ const double t31 = vertices[2](0)*eta;
+ const double t36 = t8*vertices[0](1)+vertices[1](0)*vertices[2](1)*xi-
+ t28*xi+t10*vertices[0](1)-t31*vertices[0](1)+
+ t31*vertices[3](1)+t20-t6-vertices[1](0)*
+ 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;
+
+ break;
+ };
+
+ default:
+ // not implemented at present,
+ // because of the changes above
+ Assert (false, ExcNotImplemented());
+ };
+
if (compute_support_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<Point<dim> >(n_q_points)),
+ shape_gradients (n_dofs, vector<Tensor<1,dim> >(n_q_points)),
weights (n_q_points, 0),
JxW_values (n_q_points, 0),
quadrature_points (n_q_points, Point<dim>()),
support_points (n_support_points, Point<dim>()),
- 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)),
template <int dim>
-const Point<dim> &
+const Tensor<1,dim> &
FEValuesBase<dim>::shape_grad (const unsigned int i,
const unsigned int j) const {
Assert (i<shape_gradients.size(),
// functions
for (unsigned int point=0; point<n_quadrature_points; ++point)
for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
- gradients[point] += (dof_values(shape_func) *
- shape_gradients[shape_func][point]);
+ {
+ Tensor<1,dim> tmp(shape_gradients[shape_func][point]);
+ tmp *= dof_values(shape_func);
+ gradients[point] += tmp;
+ };
};
update_flags,
fe),
unit_shape_gradients(fe.total_dofs,
- vector<Point<dim> >(quadrature.n_quadrature_points)),
+ vector<Tensor<1,dim> >(quadrature.n_quadrature_points)),
+ unit_shape_2nd_derivatives(fe.total_dofs,
+ vector<Tensor<2,dim> >(quadrature.n_quadrature_points)),
unit_shape_gradients_transform(fe.n_transform_functions,
vector<Point<dim> >(quadrature.n_quadrature_points)),
unit_quadrature_points(quadrature.get_quad_points())
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<n_transform_functions; ++i)
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,
for (unsigned int j=0; j<n_quadrature_points; ++j)
for (unsigned int s=0; s<dim; ++s)
{
- shape_gradients[i][j](s) = 0;
+ shape_gradients[i][j][s] = 0;
// (grad psi)_s =
// (grad_{\xi\eta})_b J_{bs}
// with J_{bs}=(d\xi_b)/(dx_s)
for (unsigned int b=0; b<dim; ++b)
- shape_gradients[i][j](s)
+ shape_gradients[i][j][s]
+=
- unit_shape_gradients[i][j](b) * jacobi_matrices[j](b,s);
+ unit_shape_gradients[i][j][b] * jacobi_matrices[j][b][s];
};
+
+
+ if (update_flags & update_second_derivatives)
+ for (unsigned int i=0; i<fe->total_dofs; ++i)
+ for (unsigned int j=0; j<n_quadrature_points; ++j)
+ {
+ static Tensor<2,dim> 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
// determinant
if (update_flags & update_JxW_values)
for (unsigned int i=0; i<n_quadrature_points; ++i)
- JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
+ JxW_values[i] = weights[i] / determinant(jacobi_matrices[i]);
};
update_flags,
fe),
unit_shape_gradients (n_faces_or_subfaces,
- vector<vector<Point<dim> > >(n_dofs,
- vector<Point<dim> >(n_q_points))),
+ vector<vector<Tensor<1,dim> > >(n_dofs,
+ vector<Tensor<1,dim> >(n_q_points))),
+ unit_shape_2nd_derivatives(n_dofs,
+ vector<vector<Tensor<2,dim> > >(n_dofs,
+ vector<Tensor<2,dim> >(n_q_points))),
unit_shape_gradients_transform (n_faces_or_subfaces,
vector<vector<Point<dim> > >(n_transform_functions,
vector<Point<dim> >(n_q_points))),
= 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<n_transform_functions; ++i)
// (grad_{\xi\eta})_b J_{bs}
// with J_{bs}=(d\xi_b)/(dx_s)
for (unsigned int b=0; b<dim; ++b)
- shape_gradients[i][j](s)
- += (unit_shape_gradients[face_no][i][j](b) *
- jacobi_matrices[j](b,s));
+ shape_gradients[i][j][s]
+ += (unit_shape_gradients[face_no][i][j][b] *
+ jacobi_matrices[j][b][s]);
};
+
+
+ if (update_flags & update_second_derivatives)
+ for (unsigned int i=0; i<fe->total_dofs; ++i)
+ for (unsigned int j=0; j<n_quadrature_points; ++j)
+ {
+ static Tensor<2,dim> 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
-/*------------------------------- FEFaceValues -------------------------------*/
+/*------------------------------- FESubFaceValues -------------------------------*/
template <int dim>
= fe.shape_grad(i, unit_quadrature_points[face *
GeometryInfo<dim>::
subfaces_per_face+subface][j]);
+ unit_shape_2nd_derivatives[face*GeometryInfo<dim>::subfaces_per_face+subface][i][j]
+ = fe.shape_grad_grad(i, unit_quadrature_points[face *
+ GeometryInfo<dim>::
+ subfaces_per_face+subface][j]);
};
for (unsigned int i=0; i<n_transform_functions; ++i)
for (unsigned int j=0; j<n_quadrature_points; ++j)
// (grad_{\xi\eta})_b J_{bs}
// with J_{bs}=(d\xi_b)/(dx_s)
for (unsigned int b=0; b<dim; ++b)
- shape_gradients[i][j](s)
- += (unit_shape_gradients[selected_dataset][i][j](b) *
- jacobi_matrices[j](b,s));
+ shape_gradients[i][j][s]
+ += (unit_shape_gradients[selected_dataset][i][j][b] *
+ jacobi_matrices[j][b][s]);
};
+
+ if (update_flags & update_second_derivatives)
+ for (unsigned int i=0; i<fe->total_dofs; ++i)
+ for (unsigned int j=0; j<n_quadrature_points; ++j)
+ {
+ static Tensor<2,dim> 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
Assert (cell_matrix.all_zero(), ExcObjectNotEmpty());
Assert (rhs.all_zero(), ExcObjectNotEmpty());
- const vector<vector<Point<dim> > >&gradients = fe_values.get_shape_grads ();
+ const vector<vector<Tensor<1,dim> > >&gradients = fe_values.get_shape_grads ();
const dFMatrix &values = fe_values.get_shape_values ();
vector<double> rhs_values(fe_values.n_quadrature_points);
const vector<double> &weights = fe_values.get_JxW_values ();
ExcWrongSize(cell_matrix.m(), total_dofs));
Assert (cell_matrix.all_zero(), ExcObjectNotEmpty());
- const vector<vector<Point<dim> > >&gradients = fe_values.get_shape_grads ();
+ const vector<vector<Tensor<1,dim> > >&gradients = fe_values.get_shape_grads ();
const vector<double> &weights = fe_values.get_JxW_values ();
if (coefficient != 0)
dVector &rhs,
const FEValues<2> &fe_values,
const Triangulation<2>::cell_iterator &) const {
- const vector<vector<Point<2> > >&gradients = fe_values.get_shape_grads ();
+ const vector<vector<Tensor<1,2> > >&gradients = fe_values.get_shape_grads ();
const dFMatrix &values = fe_values.get_shape_values ();
vector<double> rhs_values (fe_values.n_quadrature_points);
const vector<double> &weights = fe_values.get_JxW_values ();