* and the given quadrature points on the
* unit cell. The Jacobian matrix is to
* be computed at every quadrature point.
+ * The derivative of the jacobian matrix
+ * is the derivative with respect to the
+ * unit cell coordinates.
* This function has to be in the finite
* element class, since different finite
* elements need different transformations
* of the unit cell to a real cell.
*
- * The computation of the three fields may
+ * The computation of these fields may
* share some common code, which is why we
* put it in one function. However, it may
* not always be necessary to really
* of the Jacobi matrix and of the various
* structures to be filled.
*
- * It is provided for the finite element
- * class in one space dimension, but for
- * higher dimensions, it depends on the
- * present fe and needs reimplementation
- * by the user. This is due to the fact
+ * This function is provided for
+ * the finite element class in
+ * one space dimension, but for
+ * higher dimensions, it depends
+ * on the present fe and needs
+ * reimplementation by the
+ * user. This is due to the fact
* that the user may want to use
- * iso- or subparametric mappings of the
- * unit cell to the real cell, which
- * makes things much more complicated.
+ * iso- or subparametric mappings
+ * of the 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
+ * 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
* This function is more or less an
* interface to the #FEValues# class and
* should not be used by users unless
- * absolutely needed.
- */
+ * absolutely needed. */
virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
- const vector<Point<dim> > &unit_points,
- vector<Tensor<2,dim> > &jacobians,
- const bool compute_jacobians,
- vector<Point<dim> > &support_points,
- const bool compute_support_points,
- vector<Point<dim> > &q_points,
- const bool compute_q_points,
- const dFMatrix &shape_values_transform,
+ const vector<Point<dim> > &unit_points,
+ vector<Tensor<2,dim> > &jacobians,
+ const bool compute_jacobians,
+ vector<Tensor<3,dim> > &jacobians_grad,
+ const bool compute_jacobians_grad,
+ vector<Point<dim> > &support_points,
+ const bool compute_support_points,
+ vector<Point<dim> > &q_points,
+ const bool compute_q_points,
+ const dFMatrix &shape_values_transform,
const vector<vector<Tensor<1,dim> > > &shape_grads_transform,
const Boundary<dim> &boundary) const;
const vector<Point<dim-1> > &unit_points,
const vector<Point<dim> > &global_unit_points,
vector<Tensor<2,dim> > &jacobians,
- const bool compute_jacobians,
+ const bool compute_jacobians,
+ vector<Tensor<3,dim> > &jacobians_grad,
+ const bool compute_jacobians_grad,
vector<Point<dim> > &support_points,
const bool compute_support_points,
vector<Point<dim> > &q_points,
const vector<Point<dim-1> > &unit_points,
const vector<Point<dim> > &global_unit_points,
vector<Tensor<2,dim> > &jacobians,
- const bool compute_jacobians,
+ const bool compute_jacobians,
+ vector<Tensor<3,dim> > &jacobians_grad,
+ const bool compute_jacobians_grad,
vector<Point<dim> > &q_points,
const bool compute_q_points,
vector<double> &face_jacobi_determinants,
virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
vector<Tensor<2,dim> > &jacobians,
- const bool compute_jacobians,
- vector<Point<dim> > &support_points,
- const bool compute_support_points,
- vector<Point<dim> > &q_points,
- const bool compute_q_points,
- const dFMatrix &shape_values_transform,
+ const bool compute_jacobians,
+ vector<Tensor<3,dim> > &jacobians_grad,
+ const bool compute_jacobians_grad,
+ vector<Point<dim> > &support_points,
+ const bool compute_support_points,
+ vector<Point<dim> > &q_points,
+ const bool compute_q_points,
+ const dFMatrix &shape_values_transform,
const vector<vector<Tensor<1,dim> > > &shape_grad_transform,
const Boundary<dim> &boundary) const;
virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
vector<Tensor<2,dim> > &jacobians,
- const bool compute_jacobians,
- vector<Point<dim> > &support_points,
- const bool compute_support_points,
- vector<Point<dim> > &q_points,
- const bool compute_q_points,
- const dFMatrix &shape_values_transform,
+ const bool compute_jacobians,
+ vector<Tensor<3,dim> > &jacobians_grad,
+ const bool compute_jacobians_grad,
+ vector<Point<dim> > &support_points,
+ const bool compute_support_points,
+ vector<Point<dim> > &q_points,
+ const bool compute_q_points,
+ const dFMatrix &shape_values_transform,
const vector<vector<Tensor<1,dim> > > &shape_grad_transform,
const Boundary<dim> &boundary) const;
};
* Update the second derivatives of the
* shape functions on the real cell.
*/
- update_second_derivatives
+ update_second_derivatives = 64
};
* a reference to a whole field. Usually these fields contain
* the values of all trial functions at all quadrature points.
*
- * \item #get_function_values#, #get_function_gradients#: these
- * two functions offer a simple way to avoid the detour of the
- * trial functions, if you have a finite solution (resp. the
+ * \item #get_function_values#, #get_function_grads#, #...#: these
+ * functions offer a simple way to avoid the detour of the
+ * trial functions, if you have a finite element solution (resp. the
* vector of values associated with the different trial functions.)
* Then you may want to get information from the restriction of
* the finite element function to a certain cell, e.g. the values
* you pass it a vector holding the finite element solution and the
* functions return the values or gradients of the finite element
* function restricted to the cell which was given last time the
- * #reinit# function was given.
+ * #reinit# function was given. The same applies for the functions
+ * returning higher derivatives of the solution.
*
* Though possible in principle, these functions do not call the
* #reinit# function, you have to do so yourself beforehand. On the
* documentation for the matrix itself.
*/
const vector<vector<Tensor<1,dim> > > & get_shape_grads () const;
-
+
/**
* Return the gradients of the finite
* element function characterized
void get_function_grads (const dVector &fe_function,
vector<Tensor<1,dim> > &gradients) const;
+ /**
+ * Return the 2nd derivatives of
+ * the #i#th shape function at
+ * the #j# quadrature point. If
+ * you want to get the derivatives
+ * in one of the coordinate
+ * directions, use the
+ * appropriate function of the
+ * #Tensor# class to extract one
+ * component. Since only a
+ * reference to the derivatives'
+ * values is returned, there
+ * should be no major performance
+ * drawback. The function
+ * returns the derivatives on the
+ * real element, not the
+ * reference element.
+ */
+ const Tensor<2,dim> & shape_2nd_derivative (const unsigned int function,
+ const unsigned int quadrature_point) const;
+
+ /**
+ * Return a pointer to the
+ * matrix holding all 2nd
+ * derivatives of shape functions
+ * at all integration points, on
+ * the present cell. For the
+ * format of this matrix, see the
+ * documentation for the matrix
+ * itself.
+ */
+ const vector<vector<Tensor<2,dim> > > & get_shape_2nd_derivatives () const;
+
+ /**
+ * Return the tensor of second
+ * derivatives of the finite
+ * element function characterized
+ * by #fe_function# restricted to
+ * #cell# at the quadrature points.
+ *
+ * The function assumes that the
+ * #second_derivatives# object already has
+ * the right size.
+ */
+ void get_function_2nd_derivatives (const dVector &fe_function,
+ vector<Tensor<2,dim> > &second_derivatives) const;
+
/**
* Return the position of the #i#th
* quadrature point in real space.
+template <int dim>
+inline
+const vector<vector<Tensor<2,dim> > > &
+FEValuesBase<dim>::get_shape_2nd_derivatives () const {
+ Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField());
+ return shape_2nd_derivatives;
+};
+
+
+
template <int dim>
inline
const vector<Point<dim> > &
const vector<Point<dim> > &,
vector<Tensor<2,dim> > &,
const bool,
+ vector<Tensor<3,dim> > &,
+ const bool,
vector<Point<dim> > &,
const bool,
vector<Point<dim> > &,
const vector<Point<dim> > &global_unit_points,
vector<Tensor<2,dim> > &jacobians,
const bool compute_jacobians,
+ vector<Tensor<3,dim> > &jacobians_grad,
+ const bool compute_jacobians_grad,
vector<Point<dim> > &support_points,
const bool compute_support_points,
vector<Point<dim> > &q_points,
vector<Point<dim> > dummy(total_dofs);
fill_fe_values (cell, global_unit_points,
jacobians, compute_jacobians,
+ jacobians_grad, compute_jacobians_grad,
dummy, false,
q_points, compute_q_points,
shape_values_transform, shape_gradients_transform,
const vector<Point<dim> > &global_unit_points,
vector<Tensor<2,dim> > &jacobians,
const bool compute_jacobians,
+ vector<Tensor<3,dim> > &jacobians_grad,
+ const bool compute_jacobians_grad,
vector<Point<dim> > &q_points,
const bool compute_q_points,
vector<double> &face_jacobi_determinants,
vector<Point<dim> > dummy(total_dofs);
fill_fe_values (cell, global_unit_points,
jacobians, compute_jacobians,
+ jacobians_grad, compute_jacobians_grad,
dummy, false,
q_points, compute_q_points,
shape_values_transform, shape_gradients_transform,
void FECrissCross<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
vector<Tensor<2,dim> > &jacobians,
- const bool compute_jacobians,
- vector<Point<dim> > &support_points,
+ const bool compute_jacobians,
+ vector<Tensor<3,dim> > &jacobians_grad,
+ const bool compute_jacobians_grad,
+ vector<Point<dim> > &support_points,
const bool,
- vector<Point<dim> > &q_points,
- const bool compute_q_points,
- const dFMatrix &shape_values_transform,
+ vector<Point<dim> > &q_points,
+ const bool compute_q_points,
+ const dFMatrix &shape_values_transform,
const vector<vector<Tensor<1,dim> > > &/*shape_grad_transform*/,
const Boundary<dim> &boundary) const {
Assert (jacobians.size() == unit_points.size(),
// not implemented at present
Assert (false, ExcNotImplemented());
};
+
+
+ if (compute_jacobians_grad)
+ switch (dim)
+ {
+ case 1:
+ {
+ // derivative of the
+ // jacobian is always zero
+ // for a linear mapping in 1d
+ for (unsigned int point=0; point<n_points; ++point)
+ jacobians_grad[point][0][0][0] = 0;
+ break;
+ };
+
+ case 2:
+ {
+ for (unsigned int point=0; point<n_points; ++point)
+ {
+ const double xi = unit_points[point](0);
+ const double eta= unit_points[point](1);
+
+ const double t2 = vertices[1](0)*eta;
+ const double t4 = vertices[3](0)*vertices[2](1);
+ const double t6 = vertices[0](0)*vertices[2](1);
+ const double t8 = vertices[0](0)*vertices[3](1);
+ const double t10 = vertices[3](0)*xi;
+ const double t13 = vertices[2](0)*xi;
+ const double t16 = vertices[3](0)*vertices[1](1);
+ const double t18 = vertices[0](0)*vertices[1](1);
+ const double t19 = vertices[3](0)*vertices[0](1);
+ const double t20 = -t2*vertices[3](1)-t4*eta-t6*xi+t8*xi-
+ t10*vertices[0](1)+t10*vertices[1](1)+
+ t13*vertices[0](1)-t13*vertices[1](1)+t16
+ *eta+t18+t19;
+ const double t23 = vertices[1](0)*vertices[3](1);
+ const double t26 = vertices[2](0)*eta;
+ const double t29 = vertices[1](0)*vertices[0](1);
+ const double t30 = vertices[1](0)*vertices[2](1);
+ const double t32 = -t16-t18*eta+t6*eta-t23*xi+t2*vertices[0](1)-
+ t26*vertices[0](1)+t26*vertices[3](1)-
+ t8-t29+t23+t30
+ *xi;
+ const double t33 = t20+t32;
+ const double t34 = 1/t33;
+ const double t35 = (vertices[0](1)-vertices[1](1)+
+ vertices[2](1)-vertices[3](1))*t34;
+ const double t41 = t33*t33;
+ const double t42 = 1/t41;
+ const double t43 = (-vertices[0](1)+vertices[0](1)*xi-
+ vertices[1](1)*xi+vertices[2](1)*xi+
+ vertices[3](1)-vertices[3](1)*xi)*t42;
+ const double t44 = vertices[2](0)*vertices[0](1);
+ const double t46 = -t6+t8-t19+t16+t44-
+ vertices[2](0)*vertices[1](1)-
+ t23+t30;
+ const double t50 = (vertices[0](0)-vertices[1](0)+
+ vertices[2](0)-vertices[3](0))*t34;
+ const double t54 = (-vertices[0](0)+vertices[0](0)*xi-
+ vertices[1](0)*xi+t13+
+ vertices[3](0)-t10)*t42;
+ const double t62 = (-vertices[0](1)+vertices[0](1)*eta+
+ vertices[1](1)-vertices[1](1)*eta+
+ vertices[2](1)*eta-
+ vertices[3](1)*eta)*t42;
+ const double t67 = (-vertices[0](0)+vertices[0](0)*eta+
+ vertices[1](0)-t2+t26-vertices[3](0)*eta)*t42;
+ const double t70 = -t23-t4+t16-t18+t6+t29-t44+
+ vertices[2](0)*vertices[3](1);
+ jacobians_grad[point][0][0][0] = t35-t43*t46;
+ jacobians_grad[point][0][0][1] = -t50+t54*t46;
+ jacobians_grad[point][0][1][0] = t62*t46;
+ jacobians_grad[point][0][1][1] = -t67*t46;
+ jacobians_grad[point][1][0][0] = -t43*t70;
+ jacobians_grad[point][1][0][1] = t54*t70;
+ jacobians_grad[point][1][1][0] = -t35+t62*t70;
+ jacobians_grad[point][1][1][1] = t50-t67*t70;
+ };
+ break;
+
+ };
+
+ default:
+ // not implemented at present
+ Assert (false, ExcNotImplemented());
+ };
};
#endif
void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
vector<Tensor<2,dim> > &jacobians,
- const bool compute_jacobians,
+ const bool compute_jacobians,
+ vector<Tensor<3,dim> > &jacobians_grad,
+ const bool compute_jacobians_grad,
vector<Point<dim> > &support_points,
const bool compute_support_points,
vector<Point<dim> > &q_points,
is multiplied to the unit cell gradient *from the right*! be very careful
with these things.
- The following little program tests the correct behaviour:
+ The following little program tests the correct behaviour. The program can
+ also be found in the <tests> directory.
-------------------------------------------
#include <grid/tria.h>
// not implemented at present
Assert (false, ExcNotImplemented());
};
+
+ if (compute_jacobians_grad)
+ switch (dim)
+ {
+ case 1:
+ {
+ // derivative of the
+ // jacobian is always zero
+ // for a linear mapping in 1d
+ for (unsigned int point=0; point<n_points; ++point)
+ jacobians_grad[point][0][0][0] = 0;
+ break;
+ };
+
+ case 2:
+ {
+ for (unsigned int point=0; point<n_points; ++point)
+ {
+ const double xi = unit_points[point](0);
+ const double eta= unit_points[point](1);
+
+ const double t2 = vertices[1](0)*eta;
+ const double t4 = vertices[3](0)*vertices[2](1);
+ const double t6 = vertices[0](0)*vertices[2](1);
+ const double t8 = vertices[0](0)*vertices[3](1);
+ const double t10 = vertices[3](0)*xi;
+ const double t13 = vertices[2](0)*xi;
+ const double t16 = vertices[3](0)*vertices[1](1);
+ const double t18 = vertices[0](0)*vertices[1](1);
+ const double t19 = vertices[3](0)*vertices[0](1);
+ const double t20 = -t2*vertices[3](1)-t4*eta-t6*xi+t8*xi-
+ t10*vertices[0](1)+t10*vertices[1](1)+
+ t13*vertices[0](1)-t13*vertices[1](1)+t16
+ *eta+t18+t19;
+ const double t23 = vertices[1](0)*vertices[3](1);
+ const double t26 = vertices[2](0)*eta;
+ const double t29 = vertices[1](0)*vertices[0](1);
+ const double t30 = vertices[1](0)*vertices[2](1);
+ const double t32 = -t16-t18*eta+t6*eta-t23*xi+t2*vertices[0](1)-
+ t26*vertices[0](1)+t26*vertices[3](1)-
+ t8-t29+t23+t30
+ *xi;
+ const double t33 = t20+t32;
+ const double t34 = 1/t33;
+ const double t35 = (vertices[0](1)-vertices[1](1)+
+ vertices[2](1)-vertices[3](1))*t34;
+ const double t41 = t33*t33;
+ const double t42 = 1/t41;
+ const double t43 = (-vertices[0](1)+vertices[0](1)*xi-
+ vertices[1](1)*xi+vertices[2](1)*xi+
+ vertices[3](1)-vertices[3](1)*xi)*t42;
+ const double t44 = vertices[2](0)*vertices[0](1);
+ const double t46 = -t6+t8-t19+t16+t44-
+ vertices[2](0)*vertices[1](1)-
+ t23+t30;
+ const double t50 = (vertices[0](0)-vertices[1](0)+
+ vertices[2](0)-vertices[3](0))*t34;
+ const double t54 = (-vertices[0](0)+vertices[0](0)*xi-
+ vertices[1](0)*xi+t13+
+ vertices[3](0)-t10)*t42;
+ const double t62 = (-vertices[0](1)+vertices[0](1)*eta+
+ vertices[1](1)-vertices[1](1)*eta+
+ vertices[2](1)*eta-
+ vertices[3](1)*eta)*t42;
+ const double t67 = (-vertices[0](0)+vertices[0](0)*eta+
+ vertices[1](0)-t2+t26-vertices[3](0)*eta)*t42;
+ const double t70 = -t23-t4+t16-t18+t6+t29-t44+
+ vertices[2](0)*vertices[3](1);
+ jacobians_grad[point][0][0][0] = t35-t43*t46;
+ jacobians_grad[point][0][0][1] = -t50+t54*t46;
+ jacobians_grad[point][0][1][0] = t62*t46;
+ jacobians_grad[point][0][1][1] = -t67*t46;
+ jacobians_grad[point][1][0][0] = -t43*t70;
+ jacobians_grad[point][1][0][1] = t54*t70;
+ jacobians_grad[point][1][1][0] = -t35+t62*t70;
+ jacobians_grad[point][1][1][1] = t50-t67*t70;
+ };
+ break;
+
+ };
+
+ default:
+ // not implemented at present
+ Assert (false, ExcNotImplemented());
+ };
+
+
if (compute_support_points)
n_transform_functions (n_transform_functions),
shape_values (n_values_arrays, dFMatrix(n_dofs, n_q_points)),
shape_gradients (n_dofs, vector<Tensor<1,dim> >(n_q_points)),
+ shape_2nd_derivatives (n_dofs, vector<Tensor<2,dim> >(n_q_points)),
weights (n_q_points, 0),
JxW_values (n_q_points, 0),
quadrature_points (n_q_points, Point<dim>()),
+template <int dim>
+const Tensor<2,dim> &
+FEValuesBase<dim>::shape_2nd_derivative (const unsigned int i,
+ const unsigned int j) const {
+ Assert (i<shape_2nd_derivatives.size(),
+ ExcInvalidIndex (i, shape_2nd_derivatives.size()));
+ Assert (j<shape_2nd_derivatives[i].size(),
+ ExcInvalidIndex (j, shape_2nd_derivatives[i].size()));
+ Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField());
+
+ return shape_2nd_derivatives[i][j];
+};
+
+
+
+template <int dim>
+void FEValuesBase<dim>::get_function_2nd_derivatives (const dVector &fe_function,
+ vector<Tensor<2,dim> > &second_derivatives) const {
+ Assert (second_derivatives.size() == n_quadrature_points,
+ ExcWrongVectorSize(second_derivatives.size(), n_quadrature_points));
+
+ // get function values of dofs
+ // on this cell
+ dVector dof_values (total_dofs);
+ present_cell->get_dof_values (fe_function, dof_values);
+
+ // initialize with zero
+ fill_n (second_derivatives.begin(), n_quadrature_points, Tensor<2,dim>());
+
+ // add up contributions of trial
+ // functions
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
+ {
+ Tensor<2,dim> tmp(shape_2nd_derivatives[shape_func][point]);
+ tmp *= dof_values(shape_func);
+ second_derivatives[point] += tmp;
+ };
+};
+
+
+
template <int dim>
const Point<dim> &
FEValuesBase<dim>::quadrature_point (const unsigned int i) const {
present_cell = cell;
// fill jacobi matrices and real
// quadrature points
- if ((update_flags & update_jacobians) ||
- (update_flags & update_JxW_values)||
- (update_flags & update_q_points) ||
- (update_flags & update_gradients) ||
+ if ((update_flags & update_jacobians) ||
+ (update_flags & update_JxW_values) ||
+ (update_flags & update_q_points) ||
+ (update_flags & update_gradients) ||
+ (update_flags & update_second_derivatives) ||
(update_flags & update_support_points))
fe->fill_fe_values (cell,
unit_quadrature_points,
update_JxW_values |
update_gradients |
update_second_derivatives),
-// jacobi_matrices_grad,
-// update_flags & update_second_derivatives,
+ jacobi_matrices_grad,
+ update_flags & update_second_derivatives,
support_points,
update_flags & update_support_points,
quadrature_points,
selected_dataset = face_no;
// fill jacobi matrices and real
// quadrature points
- if ((update_flags & update_jacobians) ||
- (update_flags & update_JxW_values)||
- (update_flags & update_q_points) ||
- (update_flags & update_gradients) ||
- (update_flags & update_support_points) ||
+ if ((update_flags & update_jacobians) ||
+ (update_flags & update_JxW_values) ||
+ (update_flags & update_q_points) ||
+ (update_flags & update_gradients) ||
+ (update_flags & update_second_derivatives) ||
+ (update_flags & update_support_points) ||
(update_flags & update_JxW_values))
fe->fill_fe_face_values (cell,
- face_no,
- unit_face_quadrature_points,
- unit_quadrature_points[face_no],
- jacobi_matrices,
- update_flags & (update_jacobians |
- update_gradients |
- update_JxW_values),
- support_points,
- update_flags & update_support_points,
- quadrature_points,
- update_flags & update_q_points,
- face_jacobi_determinants,
- update_flags & update_JxW_values,
- normal_vectors,
- update_flags & update_normal_vectors,
- shape_values_transform[face_no],
- unit_shape_gradients_transform[face_no],
- boundary);
+ face_no,
+ unit_face_quadrature_points,
+ unit_quadrature_points[face_no],
+ jacobi_matrices,
+ update_flags & (update_jacobians |
+ update_gradients |
+ update_JxW_values |
+ update_second_derivatives),
+ jacobi_matrices_grad,
+ update_flags & update_second_derivatives,
+ support_points,
+ update_flags & update_support_points,
+ quadrature_points,
+ update_flags & update_q_points,
+ face_jacobi_determinants,
+ 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
// requested
selected_dataset = face_no*(1<<(dim-1)) + subface_no;
// fill jacobi matrices and real
// quadrature points
- if ((update_flags & update_jacobians) ||
- (update_flags & update_JxW_values)||
- (update_flags & update_q_points) ||
- (update_flags & update_gradients) ||
+ if ((update_flags & update_jacobians) ||
+ (update_flags & update_JxW_values) ||
+ (update_flags & update_q_points) ||
+ (update_flags & update_gradients) ||
+ (update_flags & update_second_derivatives) ||
(update_flags & update_JxW_values))
fe->fill_fe_subface_values (cell,
- face_no,
- subface_no,
- unit_face_quadrature_points,
- unit_quadrature_points[selected_dataset],
- jacobi_matrices,
- update_flags & (update_jacobians |
- update_gradients |
- update_JxW_values),
- quadrature_points,
- update_flags & update_q_points,
- face_jacobi_determinants,
- update_flags & update_JxW_values,
- normal_vectors,
- update_flags & update_normal_vectors,
- shape_values_transform[selected_dataset],
- unit_shape_gradients_transform[selected_dataset],
- boundary);
+ face_no,
+ subface_no,
+ unit_face_quadrature_points,
+ unit_quadrature_points[selected_dataset],
+ jacobi_matrices,
+ update_flags & (update_jacobians |
+ update_gradients |
+ update_JxW_values|
+ update_second_derivatives),
+ jacobi_matrices_grad,
+ update_flags & update_second_derivatives,
+ quadrature_points,
+ update_flags & update_q_points,
+ face_jacobi_determinants,
+ 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
// requested
perl -pi -e 's/^\s*t/const double t/g;' *2d_inverse_jacobian_grad
perl -pi -e 's/x\[(\d)\]/vertices[$1](0)/g;' *2d_inverse_jacobian_grad
perl -pi -e 's/y\[(\d)\]/vertices[$1](1)/g;' *2d_inverse_jacobian_grad
+perl -pi -e 's/inverseJacobian/jacobians_grad[point]/g;' *2d_inverse_jacobian_grad
perl -pi -e 's/([^;])\n/$1/g;' *2d_shape_grad
perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[0\] = (.*);/case $1: return Point<2>($2,/g;' *2d_shape_grad