*/
template <class InputVector>
void get_function_gradients (const InputVector& fe_function,
- std::vector<gradient_type>& values) const;
+ std::vector<gradient_type>& gradients) const;
/**
* Return the Hessians of the selected
*/
template <class InputVector>
void get_function_hessians (const InputVector& fe_function,
- std::vector<hessian_type>& values) const;
+ std::vector<hessian_type>& hessians) const;
+
+ /**
+ * Return the Laplacians of the
+ * selected scalar component of
+ * the finite element function
+ * characterized by
+ * <tt>fe_function</tt> at the
+ * quadrature points of the
+ * cell, face or subface
+ * selected the last time the
+ * <tt>reinit</tt> function of
+ * the FEValues object was
+ * called, at the quadrature
+ * points. The Laplacians are
+ * the trace of the Hessians.
+ *
+ * This function is the
+ * equivalent of the
+ * FEValuesBase::get_function_laplacians
+ * function but it only works
+ * on the selected scalar
+ * component.
+ */
+ template <class InputVector>
+ void get_function_laplacians (const InputVector& fe_function,
+ std::vector<value_type>& laplacians) const;
private:
* finite element, the Hessian is a
* <code>Tensor@<3,dim@></code>.
*/
- typedef Tensor<3,dim> hessian_type;
+ typedef Tensor<3,spacedim> hessian_type;
/**
* Default constructor. Creates an
divergence_type
divergence (const unsigned int shape_function,
const unsigned int q_point) const;
-
+
/**
* Return the Hessian (the tensor of
* rank 2 of all second derivatives) of
*/
template <class InputVector>
void get_function_gradients (const InputVector& fe_function,
- std::vector<gradient_type>& values) const;
+ std::vector<gradient_type>& gradients) const;
/**
* Return the symmetrized gradients of
*/
template <class InputVector>
void get_function_symmetric_gradients (const InputVector& fe_function,
- std::vector<symmetric_gradient_type>& values) const;
+ std::vector<symmetric_gradient_type>& symmetric_gradients) const;
/**
* Return the divergence of the selected
*/
template <class InputVector>
void get_function_divergences (const InputVector& fe_function,
- std::vector<divergence_type>& values) const;
+ std::vector<divergence_type>& divergences) const;
/**
* Return the Hessians of the selected
*/
template <class InputVector>
void get_function_hessians (const InputVector& fe_function,
- std::vector<hessian_type>& values) const;
+ std::vector<hessian_type>& hessians) const;
+
+ /**
+ * Return the Laplacians of the
+ * selected vector components
+ * of the finite element
+ * function characterized by
+ * <tt>fe_function</tt> at the
+ * quadrature points of the
+ * cell, face or subface
+ * selected the last time the
+ * <tt>reinit</tt> function of
+ * the FEValues object was
+ * called, at the quadrature
+ * points. The Laplacians are
+ * the trace of the Hessians.
+ *
+ * This function is the
+ * equivalent of the
+ * FEValuesBase::get_function_laplacians
+ * function but it only works
+ * on the selected vector
+ * components.
+ */
+ template <class InputVector>
+ void get_function_laplacians (const InputVector& fe_function,
+ std::vector<value_type>& laplacians) const;
private:
/**
*/
std::vector<Point<spacedim> > normal_vectors;
- /**
+ /**
* List of outward vectors normal to the cell
* surface (line) at the quadrature points
* for the codimension 1 case,
* BlockVector<T>, or one
* of the sequential PETSc or
* Trilinos vector wrapper
- * classes..It represents a
+ * classes. It represents a
* global vector of DoF values
* associated with the DofHandler
* object with which this
* FEValues object was last
* initialized.
*
- * The output are the second derivatives
- * of the function represented by
- * these DoF values, as computed
- * in real space (as opposed to
- * on the unit cell).
+ * The output are the second
+ * derivatives of the function
+ * represented by these DoF
+ * values, as computed in real
+ * space (as opposed to on the
+ * unit cell).
*/
template <class InputVector>
void
void get_function_hessians (
const InputVector& fe_function,
const VectorSlice<const std::vector<unsigned int> >& indices,
- std::vector<Tensor<2,spacedim> >& result) const;
+ std::vector<Tensor<2,spacedim> >& hessians) const;
/**
* Access to the second
void get_function_hessians (
const InputVector& fe_function,
const VectorSlice<const std::vector<unsigned int> >& indices,
- std::vector<std::vector<Tensor<2,spacedim> > >& result,
+ std::vector<std::vector<Tensor<2,spacedim> > >& hessians,
bool quadrature_points_fastest = false) const;
/**
bool = false) const;
+ /**
+ * Compute the (scalar) Laplacian
+ * of the finite element function
+ * characterized by @p
+ * fe_function restricted to @p
+ * cell at the quadrature
+ * points. The Laplacian output
+ * vector is equivalent to
+ * getting
+ * <tt>trace(hessians)</tt>,
+ * where <tt>hessian</tt> would
+ * be the output of the
+ * get_function_hessians()
+ * function.
+ *
+ * The function assumes that the
+ * @p laplacians object
+ * already has the correct size.
+ *
+ * This function may only be used if the
+ * finite element in use is a scalar one,
+ * i.e. has only one vector component. If
+ * it is a vector-valued one, then use
+ * the other
+ * get_function_laplacians()
+ * function.
+ *
+ * The actual data type of the
+ * input vector may be either a
+ * Vector<T>,
+ * BlockVector<T>, or one
+ * of the sequential PETSc or
+ * Trilinos vector wrapper
+ * classes. It represents a
+ * global vector of DoF values
+ * associated with the DofHandler
+ * object with which this
+ * FEValues object was last
+ * initialized.
+ *
+ * The output are the traces of
+ * the second derivatives
+ * (i.e. Laplacians) of the
+ * function represented by these
+ * DoF values, as computed in
+ * real space (as opposed to on
+ * the unit cell).
+ */
+ template <class InputVector, typename number>
+ void
+ get_function_laplacians (const InputVector& fe_function,
+ std::vector<number>& laplacians) const;
+
+
+ /**
+ * Compute the (scalar) Laplacian
+ * of the finite element function
+ * characterized by @p
+ * fe_function restricted to @p
+ * cell at the quadrature
+ * points. The Laplacian output
+ * vector is equivalent to
+ * getting
+ * <tt>trace(hessians)</tt>, with
+ * <tt>hessian</tt> corresponding
+ * to the output of the
+ * get_function_hessians()
+ * function.
+ *
+ * The function assumes that the
+ * @p laplacians object
+ * already has the correct size.
+ *
+ * This function does the same as
+ * the other one with the same
+ * name, but applies to
+ * vector-valued finite elements.
+ *
+ * The actual data type of the
+ * input vector may be either a
+ * Vector<T>,
+ * BlockVector<T>, or one
+ * of the sequential PETSc or
+ * Trilinos vector wrapper
+ * classes. It represents a
+ * global vector of DoF values
+ * associated with the DofHandler
+ * object with which this
+ * FEValues object was last
+ * initialized.
+ *
+ * The output are the traces of
+ * the second derivatives (i.e.
+ * Laplacians) of the function
+ * represented by these DoF
+ * values, as computed in real
+ * space (as opposed to on the
+ * unit cell).
+ */
+ template <class InputVector, typename number>
+ void
+ get_function_laplacians (const InputVector &fe_function,
+ std::vector<Vector<number> > &laplacians) const;
+
+ /**
+ * Access to the second
+ * derivatives of a function with
+ * more flexibility. see
+ * get_function_values() with
+ * corresponding arguments.
+ */
+ template <class InputVector, typename number>
+ void get_function_laplacians (
+ const InputVector& fe_function,
+ const VectorSlice<const std::vector<unsigned int> >& indices,
+ std::vector<number>& laplacians) const;
+
+ /**
+ * Access to the second
+ * derivatives of a function with
+ * more flexibility. see
+ * get_function_values() with
+ * corresponding arguments.
+ */
+ template <class InputVector, typename number>
+ void get_function_laplacians (
+ const InputVector& fe_function,
+ const VectorSlice<const std::vector<unsigned int> >& indices,
+ std::vector<Vector<number> >& laplacians) const;
+
+ /**
+ * Access to the second
+ * derivatives of a function with
+ * more flexibility. see
+ * get_function_values() with
+ * corresponding arguments.
+ */
+ template <class InputVector, typename number>
+ void get_function_laplacians (
+ const InputVector& fe_function,
+ const VectorSlice<const std::vector<unsigned int> >& indices,
+ std::vector<std::vector<number> >& laplacians,
+ bool quadrature_points_fastest = false) const;
+
+
//@}
/**
* about degrees of
* freedom. These functions are,
* above all, the
- * <tt>get_function_value/gradients/hessians</tt>
+ * <tt>get_function_value/gradients/hessians/laplacians</tt>
* functions. If you want to call
* these functions, you have to
* call the @p reinit variants
* function was called.
*/
unsigned int present_face_index;
-
+
/**
* Store a copy of the quadrature
* formula here.
- template <int dim, int spacedim>
- template <class InputVector>
- inline
- void
- Scalar<dim,spacedim>::
- get_function_values (const InputVector &fe_function,
- std::vector<value_type> &values) const
- {
- typedef FEValuesBase<dim,spacedim> FVB;
- Assert (fe_values.update_flags & update_values,
- typename FVB::ExcAccessToUninitializedField());
- Assert (values.size() == fe_values.n_quadrature_points,
- ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points));
- Assert (fe_values.present_cell.get() != 0,
- ExcMessage ("FEValues object is not reinit'ed to any cell"));
- Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
- ExcDimensionMismatch(fe_function.size(),
- fe_values.present_cell->n_dofs_for_dof_handler()));
-
- // get function values of dofs
- // on this cell
- dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
- fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
-
- std::fill (values.begin(), values.end(), value_type());
-
- for (unsigned int shape_function=0;
- shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
- if (shape_function_data[shape_function].is_nonzero_shape_function_component)
- for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
- values[q_point] +=
- dof_values(shape_function) *
- fe_values.shape_values(shape_function_data[shape_function]
- .row_index,
- q_point);
- }
-
-
-
-
- template <int dim, int spacedim>
- template <class InputVector>
- inline
- void
- Scalar<dim,spacedim>::
- get_function_gradients (const InputVector &fe_function,
- std::vector<gradient_type> &values) const
- {
- typedef FEValuesBase<dim,spacedim> FVB;
- Assert (fe_values.update_flags & update_gradients,
- typename FVB::ExcAccessToUninitializedField());
- Assert (values.size() == fe_values.n_quadrature_points,
- ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points));
- Assert (fe_values.present_cell.get() != 0,
- ExcMessage ("FEValues object is not reinit'ed to any cell"));
- Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
- ExcDimensionMismatch(fe_function.size(),
- fe_values.present_cell->n_dofs_for_dof_handler()));
-
- // get function values of dofs
- // on this cell
- dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
- fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
-
- std::fill (values.begin(), values.end(), gradient_type());
-
- for (unsigned int shape_function=0;
- shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
- if (shape_function_data[shape_function].is_nonzero_shape_function_component)
- for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
- values[q_point] +=
- dof_values(shape_function) *
- fe_values.shape_gradients[shape_function_data[shape_function]
- .row_index][q_point];
- }
-
-
-
- template <int dim, int spacedim>
- template <class InputVector>
- inline
- void
- Scalar<dim,spacedim>::
- get_function_hessians (const InputVector &fe_function,
- std::vector<hessian_type> &values) const
- {
- typedef FEValuesBase<dim,spacedim> FVB;
- Assert (fe_values.update_flags & update_hessians,
- typename FVB::ExcAccessToUninitializedField());
- Assert (values.size() == fe_values.n_quadrature_points,
- ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points));
- Assert (fe_values.present_cell.get() != 0,
- ExcMessage ("FEValues object is not reinit'ed to any cell"));
- Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
- ExcDimensionMismatch(fe_function.size(),
- fe_values.present_cell->n_dofs_for_dof_handler()));
-
- // get function values of dofs
- // on this cell
- dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
- fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
-
- std::fill (values.begin(), values.end(), hessian_type());
-
- for (unsigned int shape_function=0;
- shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
- for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
- if (shape_function_data[shape_function].is_nonzero_shape_function_component)
- values[q_point] +=
- dof_values(shape_function) *
- fe_values.shape_hessians[shape_function_data[shape_function]
- .row_index][q_point];
- }
-
-
-
template <int dim, int spacedim>
inline
typename Vector<dim,spacedim>::value_type
return symmetrize(return_value);
}
}
-
-
-
- template <int dim, int spacedim>
- template <class InputVector>
- inline
- void
- Vector<dim,spacedim>::
- get_function_values (const InputVector &fe_function,
- std::vector<value_type> &values) const
- {
- typedef FEValuesBase<dim,spacedim> FVB;
- Assert (fe_values.update_flags & update_values,
- typename FVB::ExcAccessToUninitializedField());
- Assert (values.size() == fe_values.n_quadrature_points,
- ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points));
- Assert (fe_values.present_cell.get() != 0,
- ExcMessage ("FEValues object is not reinit'ed to any cell"));
- Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
- ExcDimensionMismatch(fe_function.size(),
- fe_values.present_cell->n_dofs_for_dof_handler()));
-
- // get function values of dofs
- // on this cell
- dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
- fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
-
- std::fill (values.begin(), values.end(), value_type());
-
- for (unsigned int shape_function=0;
- shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
- {
- const int snc = shape_function_data[shape_function].single_nonzero_component;
-
- if (snc == -2)
- // shape function is zero for the
- // selected components
- continue;
- else if (snc != -1)
- for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
- values[q_point][shape_function_data[shape_function].single_nonzero_component_index]
- += dof_values(shape_function) * fe_values.shape_values(snc,q_point);
- else
- for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
- for (unsigned int d=0; d<dim; ++d)
- if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
- values[q_point][d]
- += (dof_values(shape_function) *
- fe_values.shape_values(shape_function_data[shape_function].row_index[d],q_point));
- }
- }
-
-
-
-
- template <int dim, int spacedim>
- template <class InputVector>
- inline
- void
- Vector<dim,spacedim>::
- get_function_gradients (const InputVector &fe_function,
- std::vector<gradient_type> &values) const
- {
- typedef FEValuesBase<dim,spacedim> FVB;
- Assert (fe_values.update_flags & update_gradients,
- typename FVB::ExcAccessToUninitializedField());
- Assert (values.size() == fe_values.n_quadrature_points,
- ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points));
- Assert (fe_values.present_cell.get() != 0,
- ExcMessage ("FEValues object is not reinit'ed to any cell"));
- Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
- ExcDimensionMismatch(fe_function.size(),
- fe_values.present_cell->n_dofs_for_dof_handler()));
-
- // get function values of dofs
- // on this cell
- dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
- fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
-
- std::fill (values.begin(), values.end(), gradient_type());
-
- for (unsigned int shape_function=0;
- shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
- {
- const int snc = shape_function_data[shape_function].single_nonzero_component;
-
- if (snc == -2)
- // shape function is zero for the
- // selected components
- continue;
- else if (snc != -1)
- for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
- values[q_point][shape_function_data[shape_function].single_nonzero_component_index]
- += dof_values(shape_function) * fe_values.shape_gradients[snc][q_point];
- else
- for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
- for (unsigned int d=0; d<dim; ++d)
- if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
- values[q_point][d]
- += (dof_values(shape_function) *
- fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point]);
- }
- }
-
-
-
- template <int dim, int spacedim>
- template <class InputVector>
- inline
- void
- Vector<dim,spacedim>::
- get_function_symmetric_gradients (const InputVector &fe_function,
- std::vector<symmetric_gradient_type> &values) const
- {
- typedef FEValuesBase<dim,spacedim> FVB;
- Assert (fe_values.update_flags & update_gradients,
- typename FVB::ExcAccessToUninitializedField());
- Assert (values.size() == fe_values.n_quadrature_points,
- ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points));
- Assert (fe_values.present_cell.get() != 0,
- ExcMessage ("FEValues object is not reinit'ed to any cell"));
- Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
- ExcDimensionMismatch(fe_function.size(),
- fe_values.present_cell->n_dofs_for_dof_handler()));
-
- // get function values of dofs
- // on this cell
- dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
- fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
-
- std::fill (values.begin(), values.end(), symmetric_gradient_type());
-
- for (unsigned int shape_function=0;
- shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
- {
- const int snc = shape_function_data[shape_function].single_nonzero_component;
-
- if (snc == -2)
- // shape function is zero for the
- // selected components
- continue;
- else if (snc != -1)
- for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
- values[q_point]
- += dof_values(shape_function) *
- symmetrize_single_row (shape_function_data[shape_function].single_nonzero_component_index,
- fe_values.shape_gradients[snc][q_point]);
- else
- for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
- {
- gradient_type grad;
- for (unsigned int d=0; d<dim; ++d)
- if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
- grad[d] = (dof_values(shape_function) *
- fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point]);
- values[q_point] += symmetrize(grad);
- }
- }
- }
-
-
-
- template <int dim, int spacedim>
- template <class InputVector>
- inline
- void
- Vector<dim,spacedim>::
- get_function_divergences (const InputVector &fe_function,
- std::vector<divergence_type> &values) const
- {
- typedef FEValuesBase<dim,spacedim> FVB;
- Assert (fe_values.update_flags & update_gradients,
- typename FVB::ExcAccessToUninitializedField());
- Assert (values.size() == fe_values.n_quadrature_points,
- ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points));
- Assert (fe_values.present_cell.get() != 0,
- ExcMessage ("FEValues object is not reinit'ed to any cell"));
- Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
- ExcDimensionMismatch(fe_function.size(),
- fe_values.present_cell->n_dofs_for_dof_handler()));
-
- // get function values of dofs
- // on this cell
- dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
- fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
-
- std::fill (values.begin(), values.end(), divergence_type());
-
- for (unsigned int shape_function=0;
- shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
- {
- const int snc = shape_function_data[shape_function].single_nonzero_component;
-
- if (snc == -2)
- // shape function is zero for the
- // selected components
- continue;
- else if (snc != -1)
- for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
- values[q_point]
- += dof_values(shape_function) *
- fe_values.shape_gradients[snc][q_point][shape_function_data[shape_function].single_nonzero_component_index];
- else
- for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
- {
- for (unsigned int d=0; d<dim; ++d)
- if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
- values[q_point] +=
- dof_values(shape_function) *
- fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point][d];
- }
- }
- }
-
-
-
- template <int dim, int spacedim>
- template <class InputVector>
- inline
- void
- Vector<dim,spacedim>::
- get_function_hessians (const InputVector &fe_function,
- std::vector<hessian_type> &values) const
- {
- typedef FEValuesBase<dim,spacedim> FVB;
- Assert (fe_values.update_flags & update_hessians,
- typename FVB::ExcAccessToUninitializedField());
- Assert (values.size() == fe_values.n_quadrature_points,
- ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points));
- Assert (fe_values.present_cell.get() != 0,
- ExcMessage ("FEValues object is not reinit'ed to any cell"));
- Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
- ExcDimensionMismatch(fe_function.size(),
- fe_values.present_cell->n_dofs_for_dof_handler()));
-
- // get function values of dofs
- // on this cell
- dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
- fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
-
- std::fill (values.begin(), values.end(), hessian_type());
-
- for (unsigned int shape_function=0;
- shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
- {
- const int snc = shape_function_data[shape_function].single_nonzero_component;
-
- if (snc == -2)
- // shape function is zero for the
- // selected components
- continue;
- else if (snc != -1)
- for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
- values[q_point][shape_function_data[shape_function].single_nonzero_component_index]
- += dof_values(shape_function) * fe_values.shape_hessians[snc][q_point];
- else
- for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
- for (unsigned int d=0; d<dim; ++d)
- if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
- values[q_point][d]
- += (dof_values(shape_function) *
- fe_values.shape_hessians[shape_function_data[shape_function].row_index[d]][q_point]);
- }
- }
}
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
Tensor<1,spacedim>
FEValuesBase<dim,spacedim>::shape_grad_component (const unsigned int i,
{
template <int dim, int spacedim>
Scalar<dim,spacedim>::Scalar (const FEValuesBase<dim,spacedim> &fe_values,
- const unsigned int component)
+ const unsigned int component)
:
fe_values (fe_values),
component (component),
Assert (false, ExcInternalError());
return *this;
}
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ Scalar<dim,spacedim>::
+ get_function_values (const InputVector &fe_function,
+ std::vector<value_type> &values) const
+ {
+ typedef FEValuesBase<dim,spacedim> FVB;
+ Assert (fe_values.update_flags & update_values,
+ typename FVB::ExcAccessToUninitializedField());
+ Assert (values.size() == fe_values.n_quadrature_points,
+ ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points));
+ Assert (fe_values.present_cell.get() != 0,
+ ExcMessage ("FEValues object is not reinit'ed to any cell"));
+ Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
+ ExcDimensionMismatch(fe_function.size(),
+ fe_values.present_cell->n_dofs_for_dof_handler()));
+
+ // get function values of dofs
+ // on this cell
+ dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
+ fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+ std::fill (values.begin(), values.end(), value_type());
+
+ for (unsigned int shape_function=0;
+ shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component)
+ {
+ const double value = dof_values(shape_function);
+ const double * shape_value_ptr =
+ &fe_values.shape_values(shape_function_data[shape_function].row_index, 0);
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ values[q_point] += value * *shape_value_ptr++;
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ Scalar<dim,spacedim>::
+ get_function_gradients (const InputVector &fe_function,
+ std::vector<gradient_type> &gradients) const
+ {
+ typedef FEValuesBase<dim,spacedim> FVB;
+ Assert (fe_values.update_flags & update_gradients,
+ typename FVB::ExcAccessToUninitializedField());
+ Assert (gradients.size() == fe_values.n_quadrature_points,
+ ExcDimensionMismatch(gradients.size(), fe_values.n_quadrature_points));
+ Assert (fe_values.present_cell.get() != 0,
+ ExcMessage ("FEValues object is not reinit'ed to any cell"));
+ Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
+ ExcDimensionMismatch(fe_function.size(),
+ fe_values.present_cell->n_dofs_for_dof_handler()));
+
+ // get function values of dofs
+ // on this cell
+ dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
+ fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+ std::fill (gradients.begin(), gradients.end(), gradient_type());
+
+ for (unsigned int shape_function=0;
+ shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component)
+ {
+ const double value = dof_values(shape_function);
+ const Tensor<1,spacedim> * shape_gradient_ptr =
+ &fe_values.shape_gradients[shape_function_data[shape_function].
+ row_index][0];
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ gradients[q_point] += value * *shape_gradient_ptr++;
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ Scalar<dim,spacedim>::
+ get_function_hessians (const InputVector &fe_function,
+ std::vector<hessian_type> &hessians) const
+ {
+ typedef FEValuesBase<dim,spacedim> FVB;
+ Assert (fe_values.update_flags & update_hessians,
+ typename FVB::ExcAccessToUninitializedField());
+ Assert (hessians.size() == fe_values.n_quadrature_points,
+ ExcDimensionMismatch(hessians.size(), fe_values.n_quadrature_points));
+ Assert (fe_values.present_cell.get() != 0,
+ ExcMessage ("FEValues object is not reinit'ed to any cell"));
+ Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
+ ExcDimensionMismatch(fe_function.size(),
+ fe_values.present_cell->n_dofs_for_dof_handler()));
+
+ // get function values of dofs
+ // on this cell
+ dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
+ fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+ std::fill (hessians.begin(), hessians.end(), hessian_type());
+
+ for (unsigned int shape_function=0;
+ shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component)
+ {
+ const double value = dof_values(shape_function);
+ const Tensor<2,spacedim> * shape_hessian_ptr =
+ &fe_values.shape_hessians[shape_function_data[shape_function].
+ row_index][0];
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ hessians[q_point] += value * *shape_hessian_ptr++;
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ Scalar<dim,spacedim>::
+ get_function_laplacians (const InputVector &fe_function,
+ std::vector<value_type> &laplacians) const
+ {
+ typedef FEValuesBase<dim,spacedim> FVB;
+ Assert (fe_values.update_flags & update_hessians,
+ typename FVB::ExcAccessToUninitializedField());
+ Assert (laplacians.size() == fe_values.n_quadrature_points,
+ ExcDimensionMismatch(laplacians.size(), fe_values.n_quadrature_points));
+ Assert (fe_values.present_cell.get() != 0,
+ ExcMessage ("FEValues object is not reinit'ed to any cell"));
+ Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
+ ExcDimensionMismatch(fe_function.size(),
+ fe_values.present_cell->n_dofs_for_dof_handler()));
+
+ // get function values of dofs
+ // on this cell
+ dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
+ fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+ std::fill (laplacians.begin(), laplacians.end(), value_type());
+
+ for (unsigned int shape_function=0;
+ shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component)
+ {
+ const double value = dof_values(shape_function);
+ const unsigned int row_index = shape_function_data[shape_function].row_index;
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ laplacians[q_point] +=
+ value * trace(fe_values.shape_hessians[row_index][q_point]);
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ Vector<dim,spacedim>::
+ get_function_values (const InputVector &fe_function,
+ std::vector<value_type> &values) const
+ {
+ typedef FEValuesBase<dim,spacedim> FVB;
+ Assert (fe_values.update_flags & update_values,
+ typename FVB::ExcAccessToUninitializedField());
+ Assert (values.size() == fe_values.n_quadrature_points,
+ ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points));
+ Assert (fe_values.present_cell.get() != 0,
+ ExcMessage ("FEValues object is not reinit'ed to any cell"));
+ Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
+ ExcDimensionMismatch(fe_function.size(),
+ fe_values.present_cell->n_dofs_for_dof_handler()));
+
+ // get function values of dofs
+ // on this cell
+ dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
+ fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+ std::fill (values.begin(), values.end(), value_type());
+
+ for (unsigned int shape_function=0;
+ shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
+ {
+ const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+ if (snc == -2)
+ // shape function is zero for the
+ // selected components
+ continue;
+
+ const double value = dof_values(shape_function);
+ if (snc != -1)
+ {
+ const unsigned int comp =
+ shape_function_data[shape_function].single_nonzero_component_index;
+ const double * shape_value_ptr = &fe_values.shape_values(snc,0);
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ values[q_point][comp] += value * *shape_value_ptr++;
+ }
+ else
+ for (unsigned int d=0; d<dim; ++d)
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
+ {
+ const double * shape_value_ptr =
+ &fe_values.shape_values(shape_function_data[shape_function].row_index[d],0);
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ values[q_point][d] += value * *shape_value_ptr++;
+ }
+ }
+ }
+
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ Vector<dim,spacedim>::
+ get_function_gradients (const InputVector &fe_function,
+ std::vector<gradient_type> &gradients) const
+ {
+ typedef FEValuesBase<dim,spacedim> FVB;
+ Assert (fe_values.update_flags & update_gradients,
+ typename FVB::ExcAccessToUninitializedField());
+ Assert (gradients.size() == fe_values.n_quadrature_points,
+ ExcDimensionMismatch(gradients.size(), fe_values.n_quadrature_points));
+ Assert (fe_values.present_cell.get() != 0,
+ ExcMessage ("FEValues object is not reinit'ed to any cell"));
+ Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
+ ExcDimensionMismatch(fe_function.size(),
+ fe_values.present_cell->n_dofs_for_dof_handler()));
+
+ // get function values of dofs
+ // on this cell
+ dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
+ fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+ std::fill (gradients.begin(), gradients.end(), gradient_type());
+
+ for (unsigned int shape_function=0;
+ shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
+ {
+ const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+ if (snc == -2)
+ // shape function is zero for the
+ // selected components
+ continue;
+
+ const double value = dof_values(shape_function);
+ if (snc != -1)
+ {
+ const unsigned int comp =
+ shape_function_data[shape_function].single_nonzero_component_index;
+ const Tensor<1,spacedim> * shape_gradient_ptr =
+ &fe_values.shape_gradients[snc][0];
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ gradients[q_point][comp] += value * *shape_gradient_ptr++;
+ }
+ else
+ for (unsigned int d=0; d<dim; ++d)
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
+ {
+ const Tensor<1,spacedim> * shape_gradient_ptr =
+ &fe_values.shape_gradients[shape_function_data[shape_function].
+ row_index[d]][0];
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ gradients[q_point][d] += value * *shape_gradient_ptr++;
+ }
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ Vector<dim,spacedim>::
+ get_function_symmetric_gradients (const InputVector &fe_function,
+ std::vector<symmetric_gradient_type> &symmetric_gradients) const
+ {
+ typedef FEValuesBase<dim,spacedim> FVB;
+ Assert (fe_values.update_flags & update_gradients,
+ typename FVB::ExcAccessToUninitializedField());
+ Assert (symmetric_gradients.size() == fe_values.n_quadrature_points,
+ ExcDimensionMismatch(symmetric_gradients.size(), fe_values.n_quadrature_points));
+ Assert (fe_values.present_cell.get() != 0,
+ ExcMessage ("FEValues object is not reinit'ed to any cell"));
+ Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
+ ExcDimensionMismatch(fe_function.size(),
+ fe_values.present_cell->n_dofs_for_dof_handler()));
+
+ // get function values of dofs
+ // on this cell
+ dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
+ fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+ std::fill (symmetric_gradients.begin(), symmetric_gradients.end(),
+ symmetric_gradient_type());
+
+ for (unsigned int shape_function=0;
+ shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
+ {
+ const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+ if (snc == -2)
+ // shape function is zero for the
+ // selected components
+ continue;
+
+ const double value = dof_values(shape_function);
+ if (snc != -1)
+ {
+ const unsigned int comp =
+ shape_function_data[shape_function].single_nonzero_component_index;
+ const Tensor<1,spacedim> * shape_gradient_ptr =
+ &fe_values.shape_gradients[snc][0];
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ symmetric_gradients[q_point]
+ += value * symmetrize_single_row (comp,*shape_gradient_ptr++);
+ }
+ else
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ {
+ gradient_type grad;
+ for (unsigned int d=0; d<dim; ++d)
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
+ grad[d] = value *
+ fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point];
+ symmetric_gradients[q_point] += symmetrize(grad);
+ }
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ Vector<dim,spacedim>::
+ get_function_divergences (const InputVector &fe_function,
+ std::vector<divergence_type> &divergences) const
+ {
+ typedef FEValuesBase<dim,spacedim> FVB;
+ Assert (fe_values.update_flags & update_gradients,
+ typename FVB::ExcAccessToUninitializedField());
+ Assert (divergences.size() == fe_values.n_quadrature_points,
+ ExcDimensionMismatch(divergences.size(), fe_values.n_quadrature_points));
+ Assert (fe_values.present_cell.get() != 0,
+ ExcMessage ("FEValues object is not reinit'ed to any cell"));
+ Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
+ ExcDimensionMismatch(fe_function.size(),
+ fe_values.present_cell->n_dofs_for_dof_handler()));
+
+ // get function values of dofs
+ // on this cell
+ dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
+ fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+ std::fill (divergences.begin(), divergences.end(), divergence_type());
+
+ for (unsigned int shape_function=0;
+ shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
+ {
+ const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+ if (snc == -2)
+ // shape function is zero for the
+ // selected components
+ continue;
+
+ const double value = dof_values(shape_function);
+ if (snc != -1)
+ {
+ const unsigned int comp =
+ shape_function_data[shape_function].single_nonzero_component_index;
+ const Tensor<1,spacedim> * shape_gradient_ptr =
+ &fe_values.shape_gradients[snc][0];
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ divergences[q_point] += value * (*shape_gradient_ptr++)[comp];
+ }
+ else
+ for (unsigned int d=0; d<dim; ++d)
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
+ {
+ const Tensor<1,spacedim> * shape_gradient_ptr =
+ &fe_values.shape_gradients[shape_function_data[shape_function].
+ row_index[d]][0];
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ divergences[q_point] += value * (*shape_gradient_ptr++)[d];
+ }
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ Vector<dim,spacedim>::
+ get_function_hessians (const InputVector &fe_function,
+ std::vector<hessian_type> &hessians) const
+ {
+ typedef FEValuesBase<dim,spacedim> FVB;
+ Assert (fe_values.update_flags & update_hessians,
+ typename FVB::ExcAccessToUninitializedField());
+ Assert (hessians.size() == fe_values.n_quadrature_points,
+ ExcDimensionMismatch(hessians.size(), fe_values.n_quadrature_points));
+ Assert (fe_values.present_cell.get() != 0,
+ ExcMessage ("FEValues object is not reinit'ed to any cell"));
+ Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
+ ExcDimensionMismatch(fe_function.size(),
+ fe_values.present_cell->n_dofs_for_dof_handler()));
+
+ // get function values of dofs
+ // on this cell
+ dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
+ fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+ std::fill (hessians.begin(), hessians.end(), hessian_type());
+
+ for (unsigned int shape_function=0;
+ shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
+ {
+ const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+ if (snc == -2)
+ // shape function is zero for the
+ // selected components
+ continue;
+
+ const double value = dof_values(shape_function);
+ if (snc != -1)
+ {
+ const unsigned int comp =
+ shape_function_data[shape_function].single_nonzero_component_index;
+ const Tensor<2,spacedim> * shape_hessian_ptr =
+ &fe_values.shape_hessians[snc][0];
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ hessians[q_point][comp] += value * *shape_hessian_ptr++;
+ }
+ else
+ for (unsigned int d=0; d<dim; ++d)
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
+ {
+ const Tensor<2,spacedim> * shape_hessian_ptr =
+ &fe_values.shape_hessians[shape_function_data[shape_function].
+ row_index[d]][0];
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ hessians[q_point][d] += value * *shape_hessian_ptr++;
+ }
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ Vector<dim,spacedim>::
+ get_function_laplacians (const InputVector &fe_function,
+ std::vector<value_type> &laplacians) const
+ {
+ typedef FEValuesBase<dim,spacedim> FVB;
+ Assert (fe_values.update_flags & update_hessians,
+ typename FVB::ExcAccessToUninitializedField());
+ Assert (laplacians.size() == fe_values.n_quadrature_points,
+ ExcDimensionMismatch(laplacians.size(), fe_values.n_quadrature_points));
+ Assert (fe_values.present_cell.get() != 0,
+ ExcMessage ("FEValues object is not reinit'ed to any cell"));
+ Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
+ ExcDimensionMismatch(fe_function.size(),
+ fe_values.present_cell->n_dofs_for_dof_handler()));
+
+ // get function values of dofs
+ // on this cell
+ dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
+ fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+ std::fill (laplacians.begin(), laplacians.end(), value_type());
+
+ for (unsigned int shape_function=0;
+ shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
+ {
+ const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+ if (snc == -2)
+ // shape function is zero for the
+ // selected components
+ continue;
+
+ const double value = dof_values(shape_function);
+ if (snc != -1)
+ {
+ const unsigned int comp =
+ shape_function_data[shape_function].single_nonzero_component_index;
+ const Tensor<2,spacedim> * shape_hessian_ptr =
+ &fe_values.shape_hessians[snc][0];
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ laplacians[q_point][comp] += value * trace(*shape_hessian_ptr++);
+ }
+ else
+ for (unsigned int d=0; d<dim; ++d)
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
+ {
+ const Tensor<2,spacedim> * shape_hessian_ptr =
+ &fe_values.shape_hessians[shape_function_data[shape_function].
+ row_index[d]][0];
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ laplacians[q_point][d] += value * trace(*shape_hessian_ptr++);
+ }
+ }
+ }
}
// deal with scalar finite
// elements, so no need to check
// for non-primitivity of shape
- // functions
+ // functions. in order to increase
+ // the speed of this function, we
+ // directly access the data in the
+ // shape_values array, and
+ // increment pointers for accessing
+ // the data. this saves some lookup
+ // time and indexing. moreover, the
+ // order of the loops is such that
+ // we can access the shape_values
+ // data stored contiguously (which
+ // is also advantageous because
+ // access to dof_values is
+ // generally more expensive than
+ // access to the std::vector values
+ // - so we do the cheaper operation
+ // in the innermost loop)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
{
const double value = dof_values(shape_func);
// deal with scalar finite
// elements, so no need to check
// for non-primitivity of shape
- // functions
+ // functions. in order to increase
+ // the speed of this function, we
+ // directly access the data in the
+ // shape_values array, and
+ // increment pointers for accessing
+ // the data. this saves some lookup
+ // time and indexing. moreover, the
+ // order of the loops is such that
+ // we can access the shape_values
+ // data stored contiguously (which
+ // is also advantageous because
+ // access to the global vector
+ // fe_function is more expensive
+ // than access to the small
+ // std::vector values - so we do
+ // the cheaper operation in the
+ // innermost loop)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
{
const double value = fe_function(indices[shape_func]);
// shape function is primitive or
// not. if it is, then set its only
// non-zero component, otherwise
- // loop over components
+ // loop over components. in order
+ // to increase the speed of this
+ // function, we directly access the
+ // data in the shape_values array,
+ // and increment pointers for
+ // accessing the data. this saves
+ // some lookup time and
+ // indexing. moreover, in order of
+ // the loops is such that we can
+ // access the shape_values data
+ // stored contiguously (which is
+ // also advantageous because access
+ // to the global vector fe_function
+ // is more expensive than access to
+ // the small std::vector values -
+ // so we do the cheaper operation
+ // in the innermost loop)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
{
const double value = dof_values(shape_func);
if (fe->is_primitive(shape_func))
{
+ // compared to the scalar
+ // functions, finding the correct
+ // index in the shape_value table
+ // is more involved, since we have
+ // to find the row in shape_values
+ // that corresponds to the present
+ // shape_func. this is done
+ // manually in the same way as in
+ // shape_value_component() (that
+ // function can't be used because
+ // it doesn't return us a pointer
+ // to the data).
const unsigned int
row = fe->is_primitive() ?
shape_func : this->shape_function_to_row_table[shape_func];
for (unsigned int point=0; point<n_quadrature_points; ++point)
values[point](comp) += value * *shape_value_ptr++;
}
+ // non-primitive case (vector-valued
+ // element)
else
for (unsigned int c=0; c<n_components; ++c)
{
const double *shape_value_ptr = &this->shape_values(row, 0);
const unsigned int comp = fe->system_to_component_index(shape_func).first
+ mc * n_components;
+
if (quadrature_points_fastest)
for (unsigned int point=0; point<n_quadrature_points; ++point)
values[comp][point] += value * *shape_value_ptr++;
// deal with scalar finite
// elements, so no need to check
// for non-primitivity of shape
- // functions
+ // functions. in order to increase
+ // the speed of this function, we
+ // directly access the data in the
+ // shape_gradients array, and
+ // increment pointers for accessing
+ // the data. this saves some lookup
+ // time and indexing. moreover, the
+ // order of the loops is such that
+ // we can access the
+ // shape_gradients data stored
+ // contiguously (which is also
+ // advantageous because access to
+ // the vector dof_values is
+ // gerenally more expensive than
+ // access to the std::vector
+ // gradients - so we do the cheaper
+ // operation in the innermost loop)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
{
const double value = dof_values(shape_func);
- const Tensor<1,spacedim> *shape_gradient_ptr = &this->shape_gradients[shape_func][0];
+ const Tensor<1,spacedim> *shape_gradient_ptr
+ = &this->shape_gradients[shape_func][0];
+
for (unsigned int point=0; point<n_quadrature_points; ++point)
gradients[point] += value * *shape_gradient_ptr++;
}
// deal with scalar finite
// elements, so no need to check
// for non-primitivity of shape
- // functions
+ // functions. in order to increase
+ // the speed of this function, we
+ // directly access the data in the
+ // shape_gradients array, and
+ // increment pointers for accessing
+ // the data. this saves some lookup
+ // time and indexing. moreover, the
+ // order of the loops is such that
+ // we can access the
+ // shape_gradients data stored
+ // contiguously (which is also
+ // advantageous because access to
+ // the global vector fe_function is
+ // more expensive than access to
+ // the small std::vector gradients
+ // - so we do the cheaper operation
+ // in the innermost loop)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
{
const double value = fe_function(indices[shape_func]);
- const Tensor<1,spacedim> *shape_gradient_ptr = &this->shape_gradients[shape_func][0];
+ const Tensor<1,spacedim> *shape_gradient_ptr
+ = &this->shape_gradients[shape_func][0];
+
for (unsigned int point=0; point<n_quadrature_points; ++point)
gradients[point] += value * *shape_gradient_ptr++;
}
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
{
const double value = dof_values(shape_func);
- const Tensor<2,spacedim> *shape_hessians_ptr = &this->shape_hessians[shape_func][0];
+ const Tensor<2,spacedim> *shape_hessians_ptr
+ = &this->shape_hessians[shape_func][0];
+
for (unsigned int point=0; point<n_quadrature_points; ++point)
hessians[point] += value * *shape_hessians_ptr++;
}
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
{
const double value = fe_function(indices[shape_func]);
- const Tensor<2,spacedim> *shape_hessians_ptr = &this->shape_hessians[shape_func][0];
+ const Tensor<2,spacedim> *shape_hessians_ptr
+ = &this->shape_hessians[shape_func][0];
+
for (unsigned int point=0; point<n_quadrature_points; ++point)
hessians[point] += value * *shape_hessians_ptr++;
}
+template <int dim, int spacedim>
+template <class InputVector, typename number>
+void FEValuesBase<dim,spacedim>::get_function_laplacians (
+ const InputVector &fe_function,
+ std::vector<number> &laplacians) const
+{
+ Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField());
+ Assert (fe->n_components() == 1,
+ ExcDimensionMismatch(fe->n_components(), 1));
+ Assert (laplacians.size() == n_quadrature_points,
+ ExcDimensionMismatch(laplacians.size(), n_quadrature_points));
+ Assert (present_cell.get() != 0,
+ ExcMessage ("FEValues object is not reinit'ed to any cell"));
+ Assert (fe_function.size() == present_cell->n_dofs_for_dof_handler(),
+ ExcDimensionMismatch(fe_function.size(),
+ present_cell->n_dofs_for_dof_handler()));
+
+ // get function values of dofs
+ // on this cell
+ Vector<typename InputVector::value_type> dof_values (dofs_per_cell);
+ present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+ // initialize with zero
+ std::fill_n (laplacians.begin(), n_quadrature_points, 0);
+
+ // 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. note that the
+ // laplacian is the trace of the
+ // hessian, so we use a pointer to
+ // the hessians and get their trace
+ for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+ {
+ const double value = dof_values(shape_func);
+ const Tensor<2,spacedim> *shape_hessian_ptr
+ = &this->shape_hessians[shape_func][0];
+
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ laplacians[point] += value * trace(*shape_hessian_ptr++);
+ }
+}
+
+
+
+template <int dim, int spacedim>
+template <class InputVector, typename number>
+void FEValuesBase<dim,spacedim>::get_function_laplacians (
+ const InputVector& fe_function,
+ const VectorSlice<const std::vector<unsigned int> >& indices,
+ std::vector<number> &laplacians) const
+{
+ Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField());
+ // This function fills a single
+ // component only
+ Assert (fe->n_components() == 1,
+ ExcDimensionMismatch(fe->n_components(), 1));
+ // 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 (laplacians.size() == n_quadrature_points,
+ ExcDimensionMismatch(laplacians.size(), n_quadrature_points));
+
+ // initialize with zero
+ std::fill_n (laplacians.begin(), n_quadrature_points, 0);
+
+ // 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. note that the
+ // laplacian is the trace of the
+ // hessian, so we use a pointer to
+ // the hessians and get their trace
+ for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+ {
+ const double value = fe_function(indices[shape_func]);
+ const Tensor<2,spacedim> *shape_hessian_ptr
+ = &this->shape_hessians[shape_func][0];
+
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ laplacians[point] += value * trace(*shape_hessian_ptr++);
+ }
+}
+
+
+
+template <int dim, int spacedim>
+template <class InputVector, typename number>
+void FEValuesBase<dim,spacedim>::get_function_laplacians (
+ const InputVector& fe_function,
+ std::vector<Vector<number> >& laplacians) const
+{
+//TODO: Find out how to do this assertion.
+ // This vector must correspond to a
+ // complete discretization
+// Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(),
+// ExcDimensionMismatch(fe_function.size(),
+// present_cell->get_dof_handler().n_dofs()));
+ // One entry per quadrature point
+ Assert (present_cell.get() != 0,
+ ExcMessage ("FEValues object is not reinit'ed to any cell"));
+ Assert (laplacians.size() == n_quadrature_points,
+ ExcDimensionMismatch(laplacians.size(), n_quadrature_points));
+
+ const unsigned int n_components = fe->n_components();
+ // Assert that we can write all
+ // components into the result
+ // vectors
+ for (unsigned i=0;i<laplacians.size();++i)
+ Assert (laplacians[i].size() == n_components,
+ ExcDimensionMismatch(laplacians[i].size(), n_components));
+
+ Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField());
+ Assert (fe_function.size() == present_cell->n_dofs_for_dof_handler(),
+ ExcDimensionMismatch(fe_function.size(), present_cell->n_dofs_for_dof_handler()));
+
+ // get function values of dofs
+ // on this cell
+ Vector<typename InputVector::value_type> dof_values (dofs_per_cell);
+ present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+ // initialize with zero
+ for (unsigned i=0;i<laplacians.size();++i)
+ std::fill_n (laplacians[i].begin(), laplacians[i].size(), 0);
+
+ // 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 shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+ {
+ const double value = dof_values(shape_func);
+
+ if (fe->is_primitive(shape_func))
+ {
+ const unsigned int
+ row = fe->is_primitive() ?
+ shape_func : this->shape_function_to_row_table[shape_func];
+
+ const Tensor<2,spacedim> *shape_hessian_ptr
+ = &this->shape_hessians[row][0];
+ const unsigned int comp = fe->system_to_component_index(shape_func).first;
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ laplacians[point](comp) += value * trace(*shape_hessian_ptr++);
+ }
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ if (fe->get_nonzero_components(shape_func)[c] == false)
+ continue;
+
+ const unsigned int
+ row = (this->shape_function_to_row_table[shape_func]
+ +
+ std::count (fe->get_nonzero_components(shape_func).begin(),
+ fe->get_nonzero_components(shape_func).begin()+c,
+ true));
+
+ const Tensor<2,spacedim> *shape_hessian_ptr
+ = &this->shape_hessians[row][0];
+
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ laplacians[point](c) += value * trace(*shape_hessian_ptr++);
+ }
+ }
+}
+
+
+
+template <int dim, int spacedim>
+template <class InputVector, typename number>
+void FEValuesBase<dim,spacedim>::get_function_laplacians (
+ const InputVector& fe_function,
+ const VectorSlice<const std::vector<unsigned int> >& indices,
+ std::vector<Vector<number> >& laplacians) const
+{
+ // One value per quadrature point
+ Assert (n_quadrature_points == laplacians.size(),
+ ExcDimensionMismatch(laplacians.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() * n_components / dofs_per_cell;
+
+ for (unsigned i=0;i<laplacians.size();++i)
+ Assert (laplacians[i].size() == result_components,
+ ExcDimensionMismatch(laplacians[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_hessians, ExcAccessToUninitializedField());
+
+ // initialize with zero
+ for (unsigned i=0;i<laplacians.size();++i)
+ std::fill_n (laplacians[i].begin(), laplacians[i].size(), 0);
+
+ // 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 shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+ {
+ const double value = fe_function(indices[shape_func+mc*dofs_per_cell]);
+
+ if (fe->is_primitive(shape_func))
+ {
+ const unsigned int
+ row = fe->is_primitive() ?
+ shape_func : this->shape_function_to_row_table[shape_func];
+
+ const Tensor<2,spacedim> *shape_hessian_ptr
+ = &this->shape_hessians[row][0];
+ const unsigned int comp = fe->system_to_component_index(shape_func).first
+ + mc * n_components;
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ laplacians[point](comp) += value * trace(*shape_hessian_ptr++);
+ }
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ if (fe->get_nonzero_components(shape_func)[c] == false)
+ continue;
+
+ const unsigned int
+ row = (this->shape_function_to_row_table[shape_func]
+ +
+ std::count (fe->get_nonzero_components(shape_func).begin(),
+ fe->get_nonzero_components(shape_func).begin()+c,
+ true));
+
+ const Tensor<2,spacedim> *shape_hessian_ptr
+ = &this->shape_hessians[row][0];
+ const unsigned int comp = c + mc * n_components;
+
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ laplacians[point](comp) += value * trace(*shape_hessian_ptr++);
+ }
+ }
+}
+
+
+
+template <int dim, int spacedim>
+template <class InputVector, typename number>
+void FEValuesBase<dim,spacedim>::get_function_laplacians (
+ const InputVector& fe_function,
+ const VectorSlice<const std::vector<unsigned int> >& indices,
+ std::vector<std::vector<number> >& laplacians,
+ bool quadrature_points_fastest) const
+{
+ 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() * n_components / dofs_per_cell;
+
+ // Check if the value argument is
+ // initialized to the correct sizes
+ if (quadrature_points_fastest)
+ {
+ Assert (laplacians.size() == result_components,
+ ExcDimensionMismatch(laplacians.size(), result_components));
+ for (unsigned i=0;i<laplacians.size();++i)
+ Assert (laplacians[i].size() == n_quadrature_points,
+ ExcDimensionMismatch(laplacians[i].size(), n_quadrature_points));
+ }
+ else
+ {
+ Assert(laplacians.size() == n_quadrature_points,
+ ExcDimensionMismatch(laplacians.size(), n_quadrature_points));
+ for (unsigned i=0;i<laplacians.size();++i)
+ Assert (laplacians[i].size() == result_components,
+ ExcDimensionMismatch(laplacians[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_hessians, ExcAccessToUninitializedField());
+
+ // initialize with zero
+ for (unsigned i=0;i<laplacians.size();++i)
+ std::fill_n (laplacians[i].begin(), laplacians[i].size(), 0);
+
+ // 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 shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+ {
+ const double value = fe_function(indices[shape_func+mc*dofs_per_cell]);
+
+ if (fe->is_primitive(shape_func))
+ {
+ const unsigned int
+ row = fe->is_primitive() ?
+ shape_func : this->shape_function_to_row_table[shape_func];
+
+ const Tensor<2,spacedim> *shape_hessian_ptr
+ = &this->shape_hessians[row][0];
+ const unsigned int comp = fe->system_to_component_index(shape_func).first
+ + mc * n_components;
+ if (quadrature_points_fastest)
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ laplacians[comp][point] += value * trace(*shape_hessian_ptr++);
+ else
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ laplacians[point][comp] += value * trace(*shape_hessian_ptr++);
+ }
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ if (fe->get_nonzero_components(shape_func)[c] == false)
+ continue;
+
+ const unsigned int
+ row = (this->shape_function_to_row_table[shape_func]
+ +
+ std::count (fe->get_nonzero_components(shape_func).begin(),
+ fe->get_nonzero_components(shape_func).begin()+c,
+ true));
+
+ const Tensor<2,spacedim> *shape_hessian_ptr
+ = &this->shape_hessians[row][0];
+ const unsigned int comp = c + mc * n_components;
+
+ if (quadrature_points_fastest)
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ laplacians[comp][point] += value * trace(*shape_hessian_ptr++);
+ else
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ laplacians[point][comp] += value * trace(*shape_hessian_ptr++);
+ }
+ }
+}
+
+
+
template <int dim, int spacedim>
const std::vector<Point<spacedim> > &
FEValuesBase<dim,spacedim>::get_cell_normal_vectors () const