const VectorSlice<const std::vector<unsigned int> >& indices,
std::vector<Vector<number> >& values) const;
+
+ /**
+ * Generate vector function
+ * values from an arbitrary
+ * vector.
+ *
+ * This function offers the
+ * possibility to extract
+ * function values in quadrature
+ * points from vectors not
+ * corresponding to a whole
+ * discretization.
+ *
+ * The length of the vector
+ * <tt>indices</tt> may even be a
+ * multiple of the number of dofs
+ * per cell. Then, the vectors in
+ * <tt>value</tt> should allow
+ * for the same multiple of the
+ * components of the finite
+ * element.
+ *
+ * Depending on the last
+ * argument, the outer vector of
+ * <tt>values</tt> has either the
+ * length of the quadrature rule
+ * (<tt>quadrature_points_fastest
+ * == false</tt>) or the length
+ * of components to be filled
+ * <tt>quadrature_points_fastest
+ * == true</tt>. If <tt>p</tt> is
+ * the xurrent quadrature point
+ * number and <tt>i</tt> is the
+ * vector component of the
+ * solution desired, the access
+ * to <tt>values</tt> is
+ * <tt>values[p][i]</tt> if
+ * <tt>quadrature_points_fastest
+ * == false</tt>, and
+ * <tt>values[i][p]</tt>
+ * otherwise.
+ *
+ * You may want to use this
+ * function, if you want to
+ * access just a single block
+ * from a BlockVector, if you
+ * have a multi-level vector or
+ * if you already have a local
+ * representation of your finite
+ * element data.
+ *
+ * Since this function allows for
+ * fairly general combinations of
+ * argument sizes, be aware that
+ * the checks on the arguments
+ * may not detect errors.
+ */
+ template <class InputVector, typename number>
+ void get_function_values (const InputVector& fe_function,
+ const VectorSlice<const std::vector<unsigned int> >& indices,
+ std::vector<std::vector<number> >& values,
+ bool quadrature_points_fastest) const;
+
/**
* Compute the gradients of the finite
* element function characterized
template <class InputVector>
void get_function_grads (const InputVector& fe_function,
const VectorSlice<const std::vector<unsigned int> >& indices,
- std::vector<std::vector<Tensor<1,dim> > >& gradients) const;
+ std::vector<std::vector<Tensor<1,dim> > >& gradients,
+ bool quadrature_points_fastest = false) const;
/**
* Compute the tensor of second
template <class InputVector>
void
get_function_2nd_derivatives (const InputVector &fe_function,
- std::vector<std::vector<Tensor<2,dim> > > &second_derivatives) const;
+ std::vector<std::vector<Tensor<2,dim> > > &second_derivatives,
+ bool quadrature_points_fastest = false) const;
//@}
/**
// result may be a multiple of the
// number of components of the
// finite element
- const unsigned int result_components = indices.size() / dofs_per_cell;
+ const unsigned int result_components = indices.size() * n_components / dofs_per_cell;
for (unsigned i=0;i<values.size();++i)
Assert (values[i].size() == result_components,
* shape_value(shape_func, point));
else
for (unsigned int c=0; c<n_components; ++c)
- values[point](c) += (fe_function(indices[shape_func]) *
- shape_value_component(shape_func, point, c));
+ values[point](c+mc*n_components)
+ += (fe_function(indices[shape_func])
+ * shape_value_component(shape_func, point, c));
+}
+
+
+
+template <int dim>
+template <class InputVector, typename number>
+void FEValuesBase<dim>::get_function_values (
+ const InputVector& fe_function,
+ const VectorSlice<const std::vector<unsigned int> >& indices,
+ std::vector<std::vector<number> >& values,
+ 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 (values.size() == result_components,
+ ExcDimensionMismatch(values.size(), result_components));
+ for (unsigned i=0;i<values.size();++i)
+ Assert (values[i].size() == n_quadrature_points,
+ ExcDimensionMismatch(values[i].size(), n_quadrature_points));
+ }
+ else
+ {
+ Assert(values.size() == n_quadrature_points,
+ ExcDimensionMismatch(values.size(), n_quadrature_points));
+ 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(), 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
+ if (quadrature_points_fastest)
+ 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[fe->system_to_component_index(shape_func).first
+ +mc * n_components][point]
+ += (fe_function(indices[shape_func+mc*dofs_per_cell])
+ * shape_value(shape_func, point));
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ values[c+mc*n_components][point]
+ += (fe_function(indices[shape_func])
+ * shape_value_component(shape_func, point, c));
+ else
+ 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_value(shape_func, point));
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ values[point][c+mc*n_components]
+ += (fe_function(indices[shape_func])
+ * shape_value_component(shape_func, point, c));
}
{
Tensor<1,dim> tmp = this->shape_grad(shape_func,point);
tmp *= dof_values(shape_func);
- gradients[point][fe->system_to_component_index(shape_func).first]
+ gradients[fe->system_to_component_index(shape_func).first][point]
+= tmp;
}
else
{
Tensor<1,dim> tmp = this->shape_grad_component(shape_func,point,c);
tmp *= dof_values(shape_func);
- gradients[point][c] += tmp;
- };
+ gradients[c][point] += tmp;
+ }
}
void FEValuesBase<dim>::get_function_grads (
const InputVector& fe_function,
const VectorSlice<const std::vector<unsigned int> >& indices,
- std::vector<std::vector<Tensor<1,dim> > >& values) const
+ std::vector<std::vector<Tensor<1,dim> > >& values,
+ bool quadrature_points_fastest) const
{
- // One value per quadrature point
- Assert (n_quadrature_points == values.size(),
- ExcDimensionMismatch(values.size(), n_quadrature_points));
-
const unsigned int n_components = fe->n_components();
// Size of indices must be a
// result may be a multiple of the
// number of components of the
// finite element
- const unsigned int result_components = indices.size() / dofs_per_cell;
+ const unsigned int result_components = indices.size() * n_components / dofs_per_cell;
- for (unsigned i=0;i<values.size();++i)
- Assert (values[i].size() == result_components,
- ExcDimensionMismatch(values[i].size(), result_components));
+ // Check if the value argument is
+ // initialized to the correct sizes
+ if (quadrature_points_fastest)
+ {
+ Assert (values.size() == result_components,
+ ExcDimensionMismatch(values.size(), result_components));
+ for (unsigned i=0;i<values.size();++i)
+ Assert (values[i].size() == n_quadrature_points,
+ ExcDimensionMismatch(values[i].size(), n_quadrature_points));
+ }
+ else
+ {
+ Assert(values.size() == n_quadrature_points,
+ ExcDimensionMismatch(values.size(), n_quadrature_points));
+ 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
// 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));
+ if (quadrature_points_fastest)
+ 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));
+ else
+ 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[fe->system_to_component_index(shape_func).first
+ +mc * n_components][point]
+ += 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[c][point] += (fe_function(indices[shape_func]) *
+ shape_grad_component(shape_func, point, c));
}
void
FEValuesBase<dim>::
get_function_2nd_derivatives (const InputVector &fe_function,
- std::vector<std::vector<Tensor<2,dim> > > &second_derivs) const
+ std::vector<std::vector<Tensor<2,dim> > > &second_derivs,
+ bool quadrature_points_fastest) const
{
Assert (n_quadrature_points == second_derivs.size(),
ExcDimensionMismatch(second_derivs.size(), n_quadrature_points));
// add up contributions of trial
// functions
- 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))
- {
- Tensor<2,dim> tmp(shape_2nd_derivative(shape_func,point));
- tmp *= dof_values(shape_func);
- second_derivs[point][fe->system_to_component_index(shape_func).first]
- += tmp;
- }
- else
- for (unsigned int c=0; c<n_components; ++c)
+ if (quadrature_points_fastest)
+ 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))
+ {
+ Tensor<2,dim> tmp(shape_2nd_derivative(shape_func,point));
+ tmp *= dof_values(shape_func);
+ second_derivs[point][fe->system_to_component_index(shape_func).first]
+ += tmp;
+ }
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ Tensor<2,dim> tmp = this->shape_2nd_derivative_component(shape_func,point,c);
+ tmp *= dof_values(shape_func);
+ second_derivs[point][c] += tmp;
+ }
+ else
+ 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))
{
- Tensor<2,dim> tmp = this->shape_2nd_derivative_component(shape_func,point,c);
+ Tensor<2,dim> tmp(shape_2nd_derivative(shape_func,point));
tmp *= dof_values(shape_func);
- second_derivs[point][c] += tmp;
- };
+ second_derivs[fe->system_to_component_index(shape_func).first][point]
+ += tmp;
+ }
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ Tensor<2,dim> tmp = this->shape_2nd_derivative_component(shape_func,point,c);
+ tmp *= dof_values(shape_func);
+ second_derivs[c][point] += tmp;
+ }
}