* derivatives.
*/
typedef std::vector<std::vector<Tensor<2,dim> > > GradGradVector;
-
+
/**
* Store the values of the shape
* functions at the quadrature
* Destructor.
*/
~FEValuesBase ();
+ /// @name ShapeAccess Access to shape function values
+ //@{
+
/**
* Value of a shape function at a
* quadrature point on the cell,
const unsigned int point_no,
const unsigned int component) const;
+ /**
+ * Compute the gradient of the
+ * @p ith shape function at the
+ * @p j quadrature point with
+ * respect to real cell
+ * coordinates. If you want to
+ * get the derivative in one of
+ * the coordinate directions, use
+ * the appropriate function of
+ * the Tensor class to
+ * extract one component. Since
+ * only a reference to the
+ * gradient's value is returned,
+ * there should be no major
+ * performance drawback.
+ *
+ * If the shape function is
+ * vector-valued, then this
+ * returns the only non-zero
+ * component. If the shape
+ * function has more than one
+ * non-zero component (i.e. it is
+ * not primitive), then throw an
+ * exception of type
+ * ExcShapeFunctionNotPrimitive. In
+ * that case, use the
+ * shape_grad_component()
+ * function.
+ */
+ const Tensor<1,dim> &
+ shape_grad (const unsigned int function,
+ const unsigned int quadrature_point) const;
+
+ /**
+ * Return one vector component of
+ * the gradient of a shape function
+ * at a quadrature point. If the
+ * finite element is scalar, then
+ * only component zero is allowed
+ * and the return value equals
+ * that of the shape_grad()
+ * function. If the finite
+ * element is vector valued but
+ * all shape functions are
+ * primitive (i.e. they are
+ * non-zero in only one
+ * component), then the value
+ * returned by shape_grad()
+ * equals that of this function
+ * for exactly one
+ * component. This function is
+ * therefore only of greater
+ * interest if the shape function
+ * is not primitive, but then it
+ * is necessary since the other
+ * function cannot be used.
+ */
+ Tensor<1,dim>
+ shape_grad_component (const unsigned int function_no,
+ const unsigned int point_no,
+ const unsigned int component) const;
+
+ /**
+ * Second derivatives of
+ * the @p function_noth shape function at
+ * the @p point_noth quadrature point
+ * with respect to real cell
+ * coordinates. If you want to
+ * get the derivatives in one of
+ * the coordinate directions, use
+ * the appropriate function of
+ * the @p Tensor class to
+ * extract one component. Since
+ * only a reference to the
+ * derivative values is returned,
+ * there should be no major
+ * performance drawback.
+ *
+ * If the shape function is
+ * vector-valued, then this
+ * returns the only non-zero
+ * component. If the shape
+ * function has more than one
+ * non-zero component (i.e. it is
+ * not primitive), then throw an
+ * exception of type
+ * @p ExcShapeFunctionNotPrimitive. In
+ * that case, use the
+ * shape_grad_grad_component()
+ * function.
+ */
+ const Tensor<2,dim> &
+ shape_2nd_derivative (const unsigned int function_no,
+ const unsigned int point_no) const;
+
+
+ /**
+ * Return one vector component of
+ * the gradient of a shape
+ * function at a quadrature
+ * point. If the finite element
+ * is scalar, then only component
+ * zero is allowed and the return
+ * value equals that of the
+ * @p shape_2nd_derivative
+ * function. If the finite
+ * element is vector valued but
+ * all shape functions are
+ * primitive (i.e. they are
+ * non-zero in only one
+ * component), then the value
+ * returned by
+ * @p shape_2nd_derivative
+ * equals that of this function
+ * for exactly one
+ * component. This function is
+ * therefore only of greater
+ * interest if the shape function
+ * is not primitive, but then it
+ * is necessary since the other
+ * function cannot be used.
+ */
+ Tensor<2,dim>
+ shape_2nd_derivative_component (const unsigned int function_no,
+ const unsigned int point_no,
+ const unsigned int component) const;
+
+
+ //@}
+ /// @name FunctionAccess Access to values of global finite element functions
+ //@{
+
/**
* Returns the values of the
* finite element function
const std::vector<unsigned int>& indices,
std::vector<Vector<number> >& values) const;
- /**
- * Compute the gradient of the
- * @p ith shape function at the
- * @p j quadrature point with
- * respect to real cell
- * coordinates. If you want to
- * get the derivative in one of
- * the coordinate directions, use
- * the appropriate function of
- * the Tensor class to
- * extract one component. Since
- * only a reference to the
- * gradient's value is returned,
- * there should be no major
- * performance drawback.
- *
- * If the shape function is
- * vector-valued, then this
- * returns the only non-zero
- * component. If the shape
- * function has more than one
- * non-zero component (i.e. it is
- * not primitive), then throw an
- * exception of type
- * ExcShapeFunctionNotPrimitive. In
- * that case, use the
- * shape_grad_component()
- * function.
- */
- const Tensor<1,dim> &
- shape_grad (const unsigned int function,
- const unsigned int quadrature_point) const;
-
- /**
- * Return one vector component of
- * the gradient of a shape function
- * at a quadrature point. If the
- * finite element is scalar, then
- * only component zero is allowed
- * and the return value equals
- * that of the shape_grad()
- * function. If the finite
- * element is vector valued but
- * all shape functions are
- * primitive (i.e. they are
- * non-zero in only one
- * component), then the value
- * returned by shape_grad()
- * equals that of this function
- * for exactly one
- * component. This function is
- * therefore only of greater
- * interest if the shape function
- * is not primitive, but then it
- * is necessary since the other
- * function cannot be used.
- */
- Tensor<1,dim>
- shape_grad_component (const unsigned int function_no,
- const unsigned int point_no,
- const unsigned int component) const;
-
/**
* Compute the gradients of the finite
* element function characterized
void get_function_grads (const InputVector &fe_function,
std::vector<std::vector<Tensor<1,dim> > > &gradients) const;
- /**
- * Second derivatives of
- * the @p function_noth shape function at
- * the @p point_noth quadrature point
- * with respect to real cell
- * coordinates. If you want to
- * get the derivatives in one of
- * the coordinate directions, use
- * the appropriate function of
- * the @p Tensor class to
- * extract one component. Since
- * only a reference to the
- * derivative values is returned,
- * there should be no major
- * performance drawback.
- *
- * If the shape function is
- * vector-valued, then this
- * returns the only non-zero
- * component. If the shape
- * function has more than one
- * non-zero component (i.e. it is
- * not primitive), then throw an
- * exception of type
- * @p ExcShapeFunctionNotPrimitive. In
- * that case, use the
- * shape_grad_grad_component()
- * function.
+ /**
+ * Function gradient access with
+ * more flexibility. see
+ * get_function_values() with
+ * corresponding arguments.
*/
- const Tensor<2,dim> &
- shape_2nd_derivative (const unsigned int function_no,
- const unsigned int point_no) const;
-
+ template <class InputVector>
+ void get_function_grads (const InputVector& fe_function,
+ const std::vector<unsigned int>& indices,
+ std::vector<Tensor<1,dim> >& gradients) const;
/**
- * Return one vector component of
- * the gradient of a shape
- * function at a quadrature
- * point. If the finite element
- * is scalar, then only component
- * zero is allowed and the return
- * value equals that of the
- * @p shape_2nd_derivative
- * function. If the finite
- * element is vector valued but
- * all shape functions are
- * primitive (i.e. they are
- * non-zero in only one
- * component), then the value
- * returned by
- * @p shape_2nd_derivative
- * equals that of this function
- * for exactly one
- * component. This function is
- * therefore only of greater
- * interest if the shape function
- * is not primitive, but then it
- * is necessary since the other
- * function cannot be used.
+ * Function gradient access with
+ * more flexibility. see
+ * get_function_values() with
+ * corresponding arguments.
*/
- Tensor<2,dim>
- shape_2nd_derivative_component (const unsigned int function_no,
- const unsigned int point_no,
- const unsigned int component) const;
-
+ template <class InputVector>
+ void get_function_grads (const InputVector& fe_function,
+ const std::vector<unsigned int>& indices,
+ std::vector<std::vector<Tensor<1,dim> > >& gradients) const;
+
/**
* Compute the tensor of second
* derivatives of the finite
void
get_function_2nd_derivatives (const InputVector &fe_function,
std::vector<std::vector<Tensor<2,dim> > > &second_derivatives) const;
+ //@}
/**
* Position of the @p ith
}
+template <int dim>
+template <class InputVector>
+void FEValuesBase<dim>::get_function_grads (
+ const InputVector& fe_function,
+ const std::vector<unsigned int>& indices,
+ std::vector<Tensor<1,dim> > &values) const
+{
+ Assert (this->update_flags & update_gradients, ExcAccessToUninitializedField());
+ // This function fills a single
+ // component only
+ Assert (fe->n_components() == 1,
+ ExcWrongNoOfComponents());
+ // One index for each dof
+ Assert (indices.size() == dofs_per_cell,
+ ExcDimensionMismatch(indices.size(), dofs_per_cell));
+ // This vector has one entry for
+ // each quadrature point
+ Assert (values.size() == n_quadrature_points,
+ ExcWrongVectorSize(values.size(), n_quadrature_points));
+
+ // initialize with zero
+ std::fill_n (values.begin(), n_quadrature_points, Tensor<1,dim>());
+
+ // add up contributions of trial
+ // functions. note that here we
+ // deal with scalar finite
+ // elements, so no need to check
+ // for non-primitivity of shape
+ // functions
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+ values[point] += (fe_function(indices[shape_func]) *
+ this->shape_grad(shape_func, point));
+}
+
+
+
template <int dim>
template <class InputVector>
+template <int dim>
+template <class InputVector>
+void FEValuesBase<dim>::get_function_grads (
+ const InputVector& fe_function,
+ const std::vector<unsigned int>& indices,
+ std::vector<std::vector<Tensor<1,dim> > >& values) const
+{
+ // One value per quadrature point
+ Assert (n_quadrature_points == values.size(),
+ ExcWrongVectorSize(values.size(), n_quadrature_points));
+
+ const unsigned int n_components = fe->n_components();
+
+ // Size of indices must be a
+ // multiple of dofs_per_cell such
+ // that an integer number of
+ // function values is generated in
+ // each point.
+ Assert (indices.size() % dofs_per_cell == 0,
+ ExcNotMultiple(indices.size(), dofs_per_cell));
+
+ // The number of components of the
+ // result may be a multiple of the
+ // number of components of the
+ // finite element
+ const unsigned int result_components = indices.size() / dofs_per_cell;
+
+ for (unsigned i=0;i<values.size();++i)
+ Assert (values[i].size() == result_components,
+ ExcDimensionMismatch(values[i].size(), result_components));
+
+ // If the result has more
+ // components than the finite
+ // element, we need this number for
+ // loops filling all components
+ const unsigned int component_multiple = result_components / n_components;
+
+ Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
+
+ // initialize with zero
+ for (unsigned i=0;i<values.size();++i)
+ std::fill_n (values[i].begin(), values[i].size(), Tensor<1,dim>());
+
+ // add up contributions of trial
+ // functions. now check whether the
+ // shape function is primitive or
+ // not. if it is, then set its only
+ // non-zero component, otherwise
+ // loop over components
+ for (unsigned int mc = 0; mc < component_multiple; ++mc)
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+ if (fe->is_primitive(shape_func))
+ values[point][fe->system_to_component_index(shape_func).first
+ +mc * n_components]
+ += fe_function(indices[shape_func+mc*dofs_per_cell])
+ * shape_grad(shape_func, point);
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ values[point][c] += (fe_function(indices[shape_func]) *
+ shape_grad_component(shape_func, point, c));
+}
+
+
+
template <int dim>
template <class InputVector>
void