*/
typedef double divergence_type;
+ /**
+ * A typedef for the type of the curl
+ * of the view this class
+ * represents. Here, for a set of
+ * <code>dim</code> components of the
+ * finite element, the curl is of type
+ * <code>Tensor@<1, 1@></code> (i.e. a
+ * scalar) in 2d, and of type
+ * <code>Tensor@<1, 3@></code> (i.e. a
+ * vector) in 3d.
+ */
+ typedef Tensor<1, (spacedim == 3)? 3 : 1> curl_type;
+
/**
* A typedef for the type of second
* derivatives of the view this class
divergence (const unsigned int shape_function,
const unsigned int q_point) const;
+ /**
+ * Return the vector curl of
+ * the vector components selected by
+ * this view, for the shape function
+ * and quadrature point selected by the
+ * arguments.
+ */
+ curl_type
+ curl (const unsigned int shape_function,
+ const unsigned int q_point) const;
+
/**
* Return the Hessian (the tensor of
* rank 2 of all second derivatives) of
*
* There is no equivalent function such
* as
- * FEValuesBase::get_function_gradients
+ * FEValuesBase::get_function_divergences
* in the FEValues classes but the
* information can be obtained from
* FEValuesBase::get_function_gradients,
void get_function_divergences (const InputVector& fe_function,
std::vector<divergence_type>& divergences) const;
+ /**
+ * Return the curl 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.
+ *
+ * There is no equivalent function such
+ * as
+ * FEValuesBase::get_function_curls
+ * in the FEValues classes but the
+ * information can be obtained from
+ * FEValuesBase::get_function_gradients,
+ * of course.
+ */
+ template <class InputVector>
+ void get_function_curls (const InputVector& fe_function,
+ std::vector<curl_type>& curls) const;
+
/**
* Return the Hessians of the selected
* vector components of the finite
template <int dim, int spacedim>
inline
+ typename Vector<dim,spacedim>::curl_type
+ Vector<dim,spacedim>::curl (const unsigned int shape_function,
+ const unsigned int q_point) const
+ {
+ // this function works like in the case
+ // above
+ typedef FEValuesBase<dim,spacedim> FVB;
+
+ Assert (shape_function < fe_values.fe->dofs_per_cell,
+ ExcIndexRange (shape_function, 0, fe_values.fe->dofs_per_cell));
+ Assert (fe_values.update_flags & update_gradients,
+ typename FVB::ExcAccessToUninitializedField());
+ // same as for the scalar case except
+ // that we have one more index
+ const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+ if (snc == -2)
+ return curl_type ();
+
+ else
+ switch (dim) {
+ case 1:
+ return curl_type ();
+
+ case 2: {
+ if (snc != -1) {
+ curl_type return_value;
+
+ switch (shape_function_data[shape_function].single_nonzero_component_index) {
+ case 0: {
+ return_value[0] = -1.0 * fe_values.shape_gradients[snc][q_point][1];
+ return return_value;
+ }
+
+ default: {
+ return_value[0] = fe_values.shape_gradients[snc][q_point][2];
+ return return_value;
+ }
+ }
+ }
+
+ else {
+ curl_type return_value;
+
+ return_value[0] = 0.0;
+
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[0])
+ return_value[0] -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][1];
+
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[1])
+ return_value[0] += fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][0];
+
+ return return_value;
+ }
+ }
+
+ case 3: {
+ if (snc != -1) {
+ curl_type return_value;
+
+ switch (shape_function_data[shape_function].single_nonzero_component_index) {
+ case 0: {
+ return_value[0] = 0;
+ return_value[1] = fe_values.shape_gradients[snc][q_point][2];
+ return_value[2] = -1.0 * fe_values.shape_gradients[snc][q_point][1];
+ return return_value;
+ }
+
+ case 1: {
+ return_value[0] = -1.0 * fe_values.shape_gradients[snc][q_point][2];
+ return_value[1] = 0;
+ return_value[2] = fe_values.shape_gradients[snc][q_point][0];
+ return return_value;
+ }
+
+ default: {
+ return_value[0] = fe_values.shape_gradients[snc][q_point][1];
+ return_value[1] = -1.0 * fe_values.shape_gradients[snc][q_point][0];
+ return_value[2] = 0;
+ return return_value;
+ }
+ }
+ }
+
+ else {
+ curl_type return_value;
+
+ for (unsigned int i = 0; i < dim; ++i)
+ return_value[i] = 0.0;
+
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) {
+ return_value[1] += fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][2];
+ return_value[2] -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][1];
+ }
+
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[1]) {
+ return_value[0] -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][2];
+ return_value[2] += fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][0];
+ }
+
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[2]) {
+ return_value[0] += fe_values.shape_gradients[shape_function_data[shape_function].row_index[2]][q_point][1];
+ return_value[1] -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[2]][q_point][0];
+ }
+
+ return return_value;
+ }
+ }
+ }
+ }
+
+
+ template <int dim, int spacedim>
+ inline
typename Vector<dim,spacedim>::hessian_type
Vector<dim,spacedim>::hessian (const unsigned int shape_function,
const unsigned int q_point) const
}
}
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ Vector<dim,spacedim>::get_function_curls (const InputVector &fe_function,
+ std::vector<curl_type> &curls) const {
+ typedef FEValuesBase<dim,spacedim> FVB;
+
+ Assert (fe_values.update_flags & update_gradients,
+ typename FVB::ExcAccessToUninitializedField());
+ Assert (curls.size() == fe_values.n_quadrature_points,
+ ExcDimensionMismatch (curls.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 (curls.begin (), curls.end (), curl_type ());
+
+ switch (dim) {
+ case 2: {
+ 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 (value == 0.)
+ continue;
+
+ 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];
+
+ switch (shape_function_data[shape_function].single_nonzero_component_index) {
+ case 0: {
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+ curls[q_point][0] = -1.0 * value * (*shape_gradient_ptr++)[1];
+
+ break;
+ }
+
+ default:
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+ curls[q_point][0] = value * (*shape_gradient_ptr)[0];
+ }
+ }
+
+ else {
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+ curls[q_point][0] = 0;
+
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) {
+ const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][0];
+
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+ curls[q_point][0] -= value * (*shape_gradient_ptr++)[1];
+ }
+
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[1]) {
+ const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][0];
+
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+ curls[q_point][0] += value * (*shape_gradient_ptr++)[0];
+ }
+ }
+ }
+ break;
+ }
+
+ case 3: {
+ 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 (value == 0.)
+ continue;
+
+ 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];
+
+ switch (shape_function_data[shape_function].single_nonzero_component_index) {
+ case 0: {
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+ curls[q_point][0] = 0;
+ curls[q_point][1] = value * (*shape_gradient_ptr)[2];
+ curls[q_point][2] = -1.0 * value * (*shape_gradient_ptr++)[1];
+ }
+
+ break;
+ }
+
+ case 1: {
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+ curls[q_point][0] = -1.0 * value * (*shape_gradient_ptr)[2];
+ curls[q_point][1] = 0;
+ curls[q_point][2] = value * (*shape_gradient_ptr++)[0];
+ }
+
+ break;
+ }
+
+ default:
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+ curls[q_point][0] = value * (*shape_gradient_ptr)[1];
+ curls[q_point][1] = -1.0 * value * (*shape_gradient_ptr++)[0];
+ curls[q_point][2] = 0;
+ }
+ }
+ }
+
+ else {
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+ for (unsigned int d = 0; d < dim; ++d)
+ curls[q_point][d] = 0;
+
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) {
+ const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][0];
+
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+ curls[q_point][1] += value * (*shape_gradient_ptr++)[2];
+ curls[q_point][2] -= value * (*shape_gradient_ptr++)[1];
+ }
+ }
+
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[1]) {
+ const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][0];
+
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+ curls[q_point][0] -= value * (*shape_gradient_ptr++)[2];
+ curls[q_point][2] += value * (*shape_gradient_ptr++)[0];
+ }
+ }
+
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[2]) {
+ const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[2]][0];
+
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+ curls[q_point][0] += value * (*shape_gradient_ptr++)[1];
+ curls[q_point][1] -= value * (*shape_gradient_ptr++)[0];
+ }
+ }
+ }
+ }
+ }
+ }
+ }
- template <int dim, int spacedim>
+ template <int dim, int spacedim>
template <class InputVector>
void
Vector<dim,spacedim>::