* A namespace in which we declare "extractors", i.e. classes that when used
* as subscripts in operator[] expressions on FEValues, FEFaceValues, and
* FESubfaceValues objects extract certain components of a vector-valued
- * element. There are extractors for single scalar components as well as
- * vector components consisting of <code>dim</code> elements.
+ * element. There are extractors for single scalar components,
+ * vector components consisting of <code>dim</code> elements,
+ * and second order symmetric tensors consisting
+ * of <code>(dim*dim + dim)/2<code> components
*
* See the description of the @ref vector_valued module for examples how to
* use the features of this namespace.
Scalar (const unsigned int component);
};
-
+
/**
* Extractor for a vector of
* <code>dim</code> components of a
*/
Vector (const unsigned int first_vector_component);
};
+
+
+ /**
+ * Extractor for a symmetric tensor
+ * of a rank specified by the
+ * template argument. For a second
+ * order symmetric tensor, this
+ * represents a collection of
+ * <code>(dim*dim + dim)/2</code>
+ * components of a vector-valued
+ * element. The value of
+ * <code>dim</code> is defined by
+ * the FEValues object the
+ * extractor is applied to.
+ *
+ * The concept of
+ * extractors is defined in the
+ * documentation of the namespace
+ * FEValuesExtractors and in the @ref
+ * vector_valued module.
+ *
+ * @ingroup feaccess vector_valued
+ *
+ * @author Andrew McBride, 2009
+ */
+ template <int rank>
+ struct SymmetricTensor
+ {
+ /**
+ * The first component of the tensor
+ * view.
+ */
+ unsigned int first_tensor_component;
+
+ /**
+ * Constructor. Take the first
+ * component of the selected tensor
+ * inside the FEValues object as
+ * argument.
+ */
+ SymmetricTensor (const unsigned int first_tensor_component);
+ };
}
+
/**
* A namespace for "views" on a FEValues, FEFaceValues, or FESubfaceValues
* object. A view represents only a certain part of the whole: whereas the
* derivatives of all components of a vector-valued element, views restrict
* the attention to only a single component or a subset of components.
*
- * There are classes that present views for single scalar components as well
- * as vector components consisting of <code>dim</code> elements.
+ * There are classes that present views for single scalar components,
+ * vector components consisting of <code>dim</code> elements, and
+ * symmetric second order tensor components
+ * consisting of <code>(dim*dim + dim)/2</code> elements
*
* See the description of the @ref vector_valued module for examples how to
* use the features of this namespace.
*
* @ingroup feaccess vector_valued
*/
- template <int dim, int spacedim=dim>
+ template <int dim, int spacedim=dim>
class Scalar
{
public:
* invalid object.
*/
Scalar ();
-
+
/**
* Constructor for an object that
* represents a single scalar component
* this function is called.
*/
Scalar & operator= (const Scalar<dim,spacedim> &);
-
+
/**
* Return the value of the vector
* component selected by this view, for
template <class InputVector>
void get_function_laplacians (const InputVector& fe_function,
std::vector<value_type>& laplacians) const;
-
-
+
+
private:
/**
* A reference to the FEValuesBase object
- * we operator on.
+ * we operate on.
*/
const FEValuesBase<dim,spacedim> &fe_values;
*/
std::vector<ShapeFunctionData> shape_function_data;
};
-
+
/**
*
* @ingroup feaccess vector_valued
*/
- template <int dim, int spacedim=dim>
+ template <int dim, int spacedim=dim>
class Vector
{
public:
* symmetrized gradient is a
* <code>SymmetricTensor@<2,spacedim@></code>.
*/
- typedef SymmetricTensor<2,spacedim> symmetric_gradient_type;
+ typedef dealii::SymmetricTensor<2,spacedim> symmetric_gradient_type;
/**
* A typedef for the type of the
* course is a scalar.
*/
typedef double divergence_type;
-
+
/**
* A typedef for the type of second
* derivatives of the view this class
template <class InputVector>
void get_function_symmetric_gradients (const InputVector& fe_function,
std::vector<symmetric_gradient_type>& symmetric_gradients) const;
-
+
/**
* Return the divergence of the selected
* vector components of the finite
template <class InputVector>
void get_function_divergences (const InputVector& fe_function,
std::vector<divergence_type>& divergences) const;
-
+
/**
* Return the Hessians of the selected
* vector components of the finite
template <class InputVector>
void get_function_laplacians (const InputVector& fe_function,
std::vector<value_type>& laplacians) const;
-
+
private:
/**
* A reference to the FEValuesBase object
- * we operator on.
+ * we operate on.
*/
const FEValuesBase<dim,spacedim> &fe_values;
* much cheaper.
*/
struct ShapeFunctionData
- {
+ {
/**
* For each pair (shape
* function,component within
*/
std::vector<ShapeFunctionData> shape_function_data;
};
+
+
+ template <int rank, int dim, int spacedim = dim>
+ class SymmetricTensor;
+
+ /**
+ * A class representing a view to a set of
+ * <code>(dim*dim + dim)/2</code> components forming a
+ * symmetric second-order tensor from a
+ * vector-valued finite
+ * element. Views are discussed in the
+ * @ref vector_valued module.
+ *
+ * @ingroup feaccess vector_valued
+ *
+ * @author Andrew McBride, 2009
+ */
+ template <int dim, int spacedim>
+ class SymmetricTensor<2,dim,spacedim>
+ {
+ public:
+ /**
+ * A typedef for the data type of
+ * values of the view this class
+ * represents. Since we deal with a set
+ * of <code>(dim*dim + dim)/2</code> components
+ * (i.e. the unique components of a symmetric second-order
+ * tensor), the value type is a SymmetricTensor<2,spacedim>.
+ */
+ typedef dealii::SymmetricTensor<2, spacedim> value_type;
+
+ /**
+ * A typedef for the type of the
+ * divergence of the view this class
+ * represents. Here, for a set of
+ * of <code>(dim*dim + dim)/2</code> unique components
+ * of the finite element representing a symmetric second-order
+ * tensor, the divergence of
+ * course is a * <code>Tensor@<1,dim@></code>.
+ */
+ typedef Tensor<1, spacedim> divergence_type;
+
+ /**
+ * Default constructor. Creates an
+ * invalid object.
+ */
+ SymmetricTensor();
+
+ /**
+ * Constructor for an object that
+ * represents <code>(dim*dim + dim)/2</code>
+ * components of a
+ * FEValuesBase object (or of one of
+ * the classes derived from
+ * FEValuesBase), representing the unique
+ * components comprising a symmetric second-
+ * order tensor valued variable.
+ *
+ * The second argument denotes the
+ * index of the first component of the
+ * selected symmetric second order tensor.
+ */
+ SymmetricTensor(const FEValuesBase<dim, spacedim> &fe_values_base,
+ const unsigned int first_tensor_component);
+
+ /**
+ * Copy operator. This is not a
+ * lightweight object so we don't allow
+ * copying and generate an exception if
+ * this function is called.
+ */
+ SymmetricTensor & operator=(const SymmetricTensor<2, dim, spacedim> &);
+
+ /**
+ * Return the value of the vector
+ * components selected by this view,
+ * for the shape function and
+ * quadrature point selected by the
+ * arguments. Here, since the view
+ * represents a vector-valued part of
+ * the FEValues object with
+ * <code>(dim*dim + dim)/2</code> components
+ * (the unique components of a symmetric second-order tensor),
+ * the return type is a symmetric tensor of rank 2.
+ */
+ value_type
+ value(const unsigned int shape_function,
+ const unsigned int q_point) const;
+
+
+ /**
+ * Return the vector divergence of
+ * the vector components selected by
+ * this view, for the shape function
+ * and quadrature point selected by the
+ * arguments.
+ */
+ divergence_type
+ divergence(const unsigned int shape_function,
+ const unsigned int q_point) const;
+
+ /**
+ * Return the values 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.
+ *
+ * This function is the equivalent of
+ * the
+ * FEValuesBase::get_function_values
+ * function but it only works on the
+ * selected vector components.
+ */
+ template <class InputVector>
+ void get_function_values(const InputVector& fe_function,
+ std::vector<value_type>& values) const;
+
+ /**
+ * Return the divergence 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_gradients
+ * in the FEValues classes but the
+ * information can be obtained from
+ * FEValuesBase::get_function_gradients,
+ * of course.
+ */
+ template <class InputVector>
+ void get_function_divergences(const InputVector& fe_function,
+ std::vector<divergence_type>& divergences) const;
+
+ private:
+ /**
+ * A reference to the FEValuesBase object
+ * we operate on.
+ */
+ const FEValuesBase<dim, spacedim> &fe_values;
+
+ /**
+ * The first component of the vector
+ * this view represents of the
+ * FEValuesBase object.
+ */
+ const unsigned int first_tensor_component;
+
+ /**
+ * The number of unique
+ * components of a symmetric
+ * second-order tensor
+ */
+ static const unsigned int n_independent_components = (dim*dim + dim)/2;
+
+ /**
+ * The index of the first-order
+ * tensor representation of a
+ * symmetric second-order
+ * tensor, stored as the
+ * components of a tensor
+ */
+ std::vector < std::vector<unsigned int> > vector_to_symmetric_tensor_data;
+
+ /**
+ * A structure where for each shape
+ * function we pre-compute a bunch of
+ * data that will make later accesses
+ * much cheaper.
+ */
+ struct ShapeFunctionData
+ {
+ /**
+ * For each pair (shape
+ * function,component within
+ * vector), store whether the
+ * selected vector component may be
+ * nonzero. For primitive shape
+ * functions we know for sure
+ * whether a certain scalar
+ * component of a given shape
+ * function is nonzero, whereas for
+ * non-primitive shape functions
+ * this may not be entirely clear
+ * (e.g. for RT elements it depends
+ * on the shape of a cell).
+ */
+ bool is_nonzero_shape_function_component[n_independent_components];
+
+ /**
+ * For each pair (shape function,
+ * component within vector), store
+ * the row index within the
+ * shape_values, shape_gradients,
+ * and shape_hessians tables (the
+ * column index is the quadrature
+ * point index). If the shape
+ * function is primitive, then we
+ * can get this information from
+ * the shape_function_to_row_table
+ * of the FEValues object;
+ * otherwise, we have to work a bit
+ * harder to compute this
+ * information.
+ */
+ unsigned int row_index[n_independent_components];
+
+ /**
+ * For each shape function say the
+ * following: if only a single
+ * entry in
+ * is_nonzero_shape_function_component
+ * for this shape function is
+ * nonzero, then store the
+ * corresponding value of row_index
+ * and
+ * single_nonzero_component_index
+ * represents the index between 0
+ * and (dim^2 + dim)/2 for which it is
+ * attained. If multiple components
+ * are nonzero, then store -1. If
+ * no components are nonzero then
+ * store -2.
+ */
+ int single_nonzero_component;
+ unsigned int single_nonzero_component_index;
+ };
+
+ /**
+ * Store the data about shape
+ * functions.
+ */
+ std::vector<ShapeFunctionData> shape_function_data;
+ };
}
{
/**
* Caches for scalar and
- * vector-valued views.
+ * vector, and symmetric second-order
+ * tensor valued views.
*/
std::vector<dealii::FEValuesViews::Scalar<dim,spacedim> > scalars;
std::vector<dealii::FEValuesViews::Vector<dim,spacedim> > vectors;
+ std::vector<dealii::FEValuesViews::SymmetricTensor<2,dim,spacedim> >
+ symmetric_second_order_tensors;
/**
* Constructor.
-
+
/*!@addtogroup feaccess */
/*@{*/
*
* The purpose of this class is discussed
* on the page on @ref UpdateFlagsEssay.
- *
+ *
* @author Guido Kanschat, 2000
*/
template <int dim, int spacedim=dim>
std::vector<Tensor<3,spacedim> > jacobian_grads;
/**
- * Array of the inverse Jacobian matrices
+ * Array of the inverse Jacobian matrices
* at the quadrature points.
*/
std::vector<Tensor<2,spacedim> > inverse_jacobians;
std::vector<Point<spacedim> > normal_vectors;
/**
- * List of outward vectors normal to the cell
+ * List of outward vectors normal to the cell
* surface (line) at the quadrature points
* for the codimension 1 case,
* when spacedim=3 (=2).
* filled by Mapping.
*/
HessianVector support_jacobians;
-
+
/**
* Array of the inverse Jacobian of the
* mapping in the support points,
* filled by Mapping.
*/
HessianVector support_inverse_jacobians;
-
+
/**
* Indicate the first row which a
* given shape function occupies
* non-zero components.
*/
std::vector<unsigned int> shape_function_to_row_table;
-
+
/**
* Original update flags handed
* to the constructor of
*
* The functions of this class fall into different cathegories:
* <ul>
- * <li> shape_value(), shape_grad(), etc: return one of the values
+ * <li> shape_value(), shape_grad(), etc: return one of the values
* of this object at a time. These functions are inlined, so this
* is the suggested access to all finite element values. There
* should be no loss in performance with an optimizing compiler. If
* non-zero component, since then the above functions cannot be used,
* and you have to walk over all (or only the non-zero) components of
* the shape function using this set of functions.
- *
+ *
* <li> get_function_values(), get_function_gradients(), etc.: Compute a
* finite element function or its derivative in quadrature points.
*
* <li> reinit: initialize the FEValues object for a certain cell.
* This function is not in the present class but only in the derived
- * classes and has a variable call syntax.
+ * classes and has a variable call syntax.
* See the docs for the derived classes for more information.
* </ul>
*
* give <b>all</b> needed update flags to FEValues.
*
* The mechanisms by which this class works is also discussed
- * on the page on @ref UpdateFlagsEssay.
+ * on the page on @ref UpdateFlagsEssay.
*
* @author Wolfgang Bangerth, 1998, 2003, Guido Kanschat, 2001
*/
*/
const unsigned int dofs_per_cell;
-
+
/**
* Constructor. Set up the array
* sizes with <tt>n_q_points</tt>
const FEValuesViews::Vector<dim,spacedim> &
operator[] (const FEValuesExtractors::Vector &vector) const;
+ /**
+ * Create a view of the current FEValues
+ * object that represents a set of
+ * <code>(dim*dim + dim)/2</code> scalar components
+ * (i.e. a symmetric 2nd order tensor)
+ * of the vector-valued
+ * finite element. The concept of views
+ * is explained in the documentation of
+ * the namespace FEValuesViews and in particular
+ * in the @ref vector_valued module.
+ */
+ const FEValuesViews::SymmetricTensor<2,dim,spacedim> &
+ operator[] (const FEValuesExtractors::SymmetricTensor<2> &tensor) const;
+
//@}
-
+
/// @name ShapeAccess Access to shape function values
//@{
-
+
/**
* Value of a shape function at a
* quadrature point on the cell,
shape_hessian_component (const unsigned int function_no,
const unsigned int point_no,
const unsigned int component) const;
-
+
/**
* @deprecated Wrapper for shape_hessian_component()
*/
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
* it is a vector-valued one, then use
* the other get_function_values()
* function.
- *
+ *
* The function assumes that the
* <tt>values</tt> object already has the
- * correct size.
+ * correct size.
*
* The actual data type of the
* input vector may be either a
/**
* Generate function values from
* an arbitrary vector.
- *
+ *
* This function offers the
* possibility to extract
* function values in quadrature
* it is a vector-valued one, then use
* the other get_function_gradients()
* function.
- *
+ *
* The function assumes that the
* @p gradients object already has the
* right size.
void get_function_gradients (const InputVector &fe_function,
std::vector<std::vector<Tensor<1,spacedim> > > &gradients) const;
-
+
/**
* @deprecated Use
* get_function_gradients() instead.
* the other
* get_function_hessians()
* function.
- *
+ *
* The actual data type of the
* input vector may be either a
* Vector<T>,
get_function_hessians (const InputVector& fe_function,
std::vector<Tensor<2,spacedim> >& hessians) const;
-
+
/**
* Compute the tensor of second
* derivatives of the finite
get_function_2nd_derivatives (const InputVector&,
std::vector<std::vector<Tensor<2,spacedim> > >&,
bool = false) const;
-
+
/**
* Compute the (scalar) Laplacian
* the other
* get_function_laplacians()
* function.
- *
+ *
* The actual data type of the
* input vector may be either a
* Vector<T>,
get_function_laplacians (const InputVector& fe_function,
std::vector<number>& laplacians) const;
-
+
/**
* Compute the (scalar) Laplacian
* of the finite element function
//@}
-
+
/**
* Position of the <tt>i</tt>th
* quadrature point in real space.
* the values returned by JxW().
*/
const std::vector<double> & get_JxW_values () const;
-
+
/**
* Return the Jacobian of the
* transformation at the specified
/**
* Pointer to the array holding
- * the values returned by
+ * the values returned by
* inverse_jacobian().
*/
const std::vector<Tensor<2,spacedim> > & get_inverse_jacobians () const;
* for this object.
*/
UpdateFlags get_update_flags () const;
-
+
/**
* Return a triangulation
* iterator to the current cell.
* @ingroup Exceptions
*/
DeclException0 (ExcFENotPrimitive);
-
+
protected:
/**
* Objects of the FEValues
*
* @author Wolfgang Bangerth, 2003
*/
- class CellIteratorBase
+ class CellIteratorBase
{
public:
/**
* class.
*/
virtual ~CellIteratorBase ();
-
+
/**
* Conversion operator to an
* iterator for
*/
template <typename CI> class CellIterator;
class TriaCellIterator;
-
+
/**
* Store the cell selected last time the
* reinit() function was called. This is
* the extractor classes.
*/
std::auto_ptr<const CellIteratorBase> present_cell;
-
+
/**
* Storage for the mapping object.
*/
const SmartPointer<const Mapping<dim,spacedim> > mapping;
-
+
/**
* Store the finite element for later use.
*/
const SmartPointer<const FiniteElement<dim,spacedim> > fe;
-
+
/**
* Internal data of mapping.
*/
* and also do not implement it.
*/
FEValuesBase (const FEValuesBase &);
-
+
/**
* Copy operator. Since
* objects of this class are not
* objects.
*/
internal::FEValuesViews::Cache<dim,spacedim> fe_values_views_cache;
-
+
/**
* Make the view classes friends of this
* class, since they access internal
*/
template <int, int> friend class FEValuesViews::Scalar;
template <int, int> friend class FEValuesViews::Vector;
+ template <int, int, int> friend class FEValuesViews::SymmetricTensor;
};
* This function implements the initialization routines for
* FEValuesBase, if values in quadrature points of a cell are
* needed. For further documentation see this class.
- *
+ *
* @author Wolfgang Bangerth, 1998, Guido Kanschat, 2001
*/
template <int dim, int spacedim=dim>
* to <code>dim</code>.
*/
static const unsigned int integral_dimension = dim;
-
+
/**
* Constructor. Gets cell
* independent data from mapping
FEValues (const FiniteElement<dim,spacedim> &fe,
const Quadrature<dim> &quadrature,
const UpdateFlags update_flags);
-
+
/**
* Reinitialize the gradients,
* Jacobi determinants, etc for
* stored by this object.
*/
const Quadrature<dim> & get_quadrature () const;
-
+
/**
* Determine an estimate for the
* memory consumption (in bytes)
* FEValues/hpFEValues).
*/
const FEValues<dim,spacedim> & get_present_fe_values () const;
-
+
private:
/**
* Store a copy of the quadrature
*
* @author Wolfgang Bangerth, 1998, Guido Kanschat, 2000, 2001
*/
-template <int dim, int spacedim=dim>
+template <int dim, int spacedim=dim>
class FEFaceValuesBase : public FEValuesBase<dim,spacedim>
{
public:
* is normalized to one.
*/
const Point<dim> & normal_vector (const unsigned int i) const;
-
+
/**
* Boundary form of the
* transformation of the cell at
* <tt>n.ds</tt>.
*/
const Tensor<1,spacedim> & boundary_form (const unsigned int i) const;
-
+
/**
* Return the list of outward normal
* vectors to the cell at the
* stored by this object.
*/
const Quadrature<dim-1> & get_quadrature () const;
-
+
/**
* Determine an estimate for the
* memory consumption (in bytes)
* of this object.
*/
unsigned int memory_consumption () const;
-
+
protected:
/**
*
* @author Wolfgang Bangerth, 1998, Guido Kanschat, 2000, 2001
*/
-template <int dim, int spacedim=dim>
+template <int dim, int spacedim=dim>
class FEFaceValues : public FEFaceValuesBase<dim,spacedim>
{
public:
*/
void reinit (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no);
-
+
/**
* Return a reference to this
* very object.
* Do work common to the two
* constructors.
*/
- void initialize (const UpdateFlags update_flags);
+ void initialize (const UpdateFlags update_flags);
/**
* The reinit() functions do
*
* @author Wolfgang Bangerth, 1998, Guido Kanschat, 2000, 2001
*/
-template <int dim, int spacedim=dim>
+template <int dim, int spacedim=dim>
class FESubfaceValues : public FEFaceValuesBase<dim,spacedim>
{
public:
void reinit (const typename hp::DoFHandler<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int subface_no);
-
+
/**
* Reinitialize the gradients,
* Jacobi determinants, etc for
void reinit (const typename MGDoFHandler<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int subface_no);
-
+
/**
* Reinitialize the gradients,
* Jacobi determinants, etc for
:
first_vector_component (first_vector_component)
{}
+
+ template <int rank>
+ inline
+ SymmetricTensor<rank>::SymmetricTensor (const unsigned int first_tensor_component)
+ :
+ first_tensor_component (first_tensor_component)
+ {}
}
namespace FEValuesViews
{
- template <int dim, int spacedim>
+ template <int dim, int spacedim>
inline
typename Scalar<dim,spacedim>::value_type
Scalar<dim,spacedim>::value (const unsigned int shape_function,
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_values,
- typename FVB::ExcAccessToUninitializedField());
+ typename FVB::ExcAccessToUninitializedField());
// an adaptation of the
// FEValuesBase::shape_value_component
else
return 0;
}
-
- template <int dim, int spacedim>
+
+ template <int dim, int spacedim>
inline
typename Scalar<dim,spacedim>::gradient_type
Scalar<dim,spacedim>::gradient (const unsigned int shape_function,
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());
+ typename FVB::ExcAccessToUninitializedField());
// an adaptation of the
// FEValuesBase::shape_grad_component
- template <int dim, int spacedim>
+ template <int dim, int spacedim>
inline
typename Scalar<dim,spacedim>::hessian_type
Scalar<dim,spacedim>::hessian (const unsigned int shape_function,
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_hessians,
- typename FVB::ExcAccessToUninitializedField());
+ typename FVB::ExcAccessToUninitializedField());
// an adaptation of the
// FEValuesBase::shape_grad_component
- template <int dim, int spacedim>
+ template <int dim, int spacedim>
inline
typename Vector<dim,spacedim>::value_type
Vector<dim,spacedim>::value (const unsigned int shape_function,
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_values,
- typename FVB::ExcAccessToUninitializedField());
+ typename FVB::ExcAccessToUninitializedField());
// same as for the scalar case except
// that we have one more index
else
{
value_type return_value;
- for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int d=0; d<dim; ++d)
if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
return_value[d]
= fe_values.shape_values(shape_function_data[shape_function].row_index[d],q_point);
-
+
return return_value;
}
}
-
- template <int dim, int spacedim>
+
+ template <int dim, int spacedim>
inline
typename Vector<dim,spacedim>::gradient_type
Vector<dim,spacedim>::gradient (const unsigned int shape_function,
const unsigned int q_point) const
{
- typedef FEValuesBase<dim,spacedim> FVB;
+ 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());
+ typename FVB::ExcAccessToUninitializedField());
// same as for the scalar case except
// that we have one more index
else
{
gradient_type return_value;
- for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int d=0; d<dim; ++d)
if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
return_value[d]
= fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point];
- template <int dim, int spacedim>
+ template <int dim, int spacedim>
inline
typename Vector<dim,spacedim>::divergence_type
Vector<dim,spacedim>::divergence (const unsigned int shape_function,
{
// this function works like in
// the case above
- typedef FEValuesBase<dim,spacedim> FVB;
+ 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());
+ typename FVB::ExcAccessToUninitializedField());
// same as for the scalar case except
// that we have one more index
else
{
divergence_type return_value = 0;
- for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int d=0; d<dim; ++d)
if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
return_value
+= fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point][d];
return return_value;
}
}
-
- template <int dim, int spacedim>
+
+ 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
{
// this function works like in
- // the case above
- typedef FEValuesBase<dim,spacedim> FVB;
+ // 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_hessians,
- typename FVB::ExcAccessToUninitializedField());
+ typename FVB::ExcAccessToUninitializedField());
// same as for the scalar case except
// that we have one more index
else
{
hessian_type return_value;
- for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int d=0; d<dim; ++d)
if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
return_value[d]
= fe_values.shape_hessians[shape_function_data[shape_function].row_index[d]][q_point];
* equal to zero.
*/
inline
- SymmetricTensor<2,1>
+ dealii::SymmetricTensor<2,1>
symmetrize_single_row (const unsigned int n,
const Tensor<1,1> &t)
{
const unsigned int dim = 1;
Assert (n < dim, ExcIndexRange (n, 0, dim));
-
+
const double array[1] = { t[0] };
- return SymmetricTensor<2,1>(array);
+ return dealii::SymmetricTensor<2,1>(array);
}
inline
- SymmetricTensor<2,2>
+ dealii::SymmetricTensor<2,2>
symmetrize_single_row (const unsigned int n,
const Tensor<1,2> &t)
{
case 0:
{
const double array[3] = { t[0], 0, t[1]/2 };
- return SymmetricTensor<2,2>(array);
+ return dealii::SymmetricTensor<2,2>(array);
}
case 1:
{
const double array[3] = { 0, t[1], t[0]/2 };
- return SymmetricTensor<2,2>(array);
+ return dealii::SymmetricTensor<2,2>(array);
}
default:
{
Assert (false, ExcIndexRange (n, 0, 2));
- return SymmetricTensor<2,2>();
+ return dealii::SymmetricTensor<2,2>();
}
}
}
inline
- SymmetricTensor<2,3>
+ dealii::SymmetricTensor<2,3>
symmetrize_single_row (const unsigned int n,
const Tensor<1,3> &t)
{
case 0:
{
const double array[6] = { t[0], 0, 0, t[1]/2, t[2]/2, 0 };
- return SymmetricTensor<2,3>(array);
+ return dealii::SymmetricTensor<2,3>(array);
}
case 1:
{
const double array[6] = { 0, t[1], 0, t[0]/2, 0, t[2]/2 };
- return SymmetricTensor<2,3>(array);
+ return dealii::SymmetricTensor<2,3>(array);
}
case 2:
{
const double array[6] = { 0, 0, t[2], 0, t[0]/2, t[1]/2 };
- return SymmetricTensor<2,3>(array);
+ return dealii::SymmetricTensor<2,3>(array);
}
default:
{
Assert (false, ExcIndexRange (n, 0, 3));
- return SymmetricTensor<2,3>();
+ return dealii::SymmetricTensor<2,3>();
}
}
}
}
-
+
template <int dim, int spacedim>
inline
Vector<dim,spacedim>::symmetric_gradient (const unsigned int shape_function,
const unsigned int q_point) const
{
- typedef FEValuesBase<dim,spacedim> FVB;
+ 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());
+ typename FVB::ExcAccessToUninitializedField());
// same as for the scalar case except
// that we have one more index
else
{
gradient_type return_value;
- for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int d=0; d<dim; ++d)
if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
return_value[d]
= fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point];
Assert (scalar.component < fe_values_views_cache.scalars.size(),
ExcIndexRange (scalar.component,
0, fe_values_views_cache.scalars.size()));
-
+
return fe_values_views_cache.scalars[scalar.component];
}
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const FEValuesViews::Vector<dim,spacedim> &
FEValuesBase<dim,spacedim>::
fe_values_views_cache.vectors.size(),
ExcIndexRange (vector.first_vector_component,
0, fe_values_views_cache.vectors.size()));
-
+
return fe_values_views_cache.vectors[vector.first_vector_component];
}
+template <int dim, int spacedim>
+inline
+const FEValuesViews::SymmetricTensor<2,dim,spacedim> &
+FEValuesBase<dim,spacedim>::
+operator[] (const FEValuesExtractors::SymmetricTensor<2> &tensor) const
+{
+ Assert (tensor.first_tensor_component <
+ fe_values_views_cache.symmetric_second_order_tensors.size(),
+ ExcIndexRange (tensor.first_tensor_component,
+ 0, fe_values_views_cache.symmetric_second_order_tensors.size()));
+ return fe_values_views_cache.symmetric_second_order_tensors[tensor.first_tensor_component];
+}
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const double &
FEValuesBase<dim,spacedim>::shape_value (const unsigned int i,
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
double
FEValuesBase<dim,spacedim>::shape_value_component (const unsigned int i,
ExcAccessToUninitializedField());
Assert (component < fe->n_components(),
ExcIndexRange(component, 0, fe->n_components()));
-
+
// if this particular shape
// function is primitive, then we
// can take a short-cut by checking
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const Tensor<1,spacedim> &
FEValuesBase<dim,spacedim>::shape_grad (const unsigned int i,
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
Tensor<1,spacedim>
FEValuesBase<dim,spacedim>::shape_grad_component (const unsigned int i,
ExcAccessToUninitializedField());
Assert (component < fe->n_components(),
ExcIndexRange(component, 0, fe->n_components()));
-
+
// if this particular shape
// function is primitive, then we
// can take a short-cut by checking
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const Tensor<2,spacedim> &
FEValuesBase<dim,spacedim>::shape_hessian (const unsigned int i,
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const Tensor<2,spacedim> &
FEValuesBase<dim,spacedim>::shape_2nd_derivative (const unsigned int i,
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
Tensor<2,spacedim>
FEValuesBase<dim,spacedim>::shape_hessian_component (const unsigned int i,
ExcAccessToUninitializedField());
Assert (component < fe->n_components(),
ExcIndexRange(component, 0, fe->n_components()));
-
+
// if this particular shape
// function is primitive, then we
// can take a short-cut by checking
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
Tensor<2,spacedim>
FEValuesBase<dim,spacedim>::shape_2nd_derivative_component (const unsigned int i,
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
-const FiniteElement<dim,spacedim> &
+const FiniteElement<dim,spacedim> &
FEValuesBase<dim,spacedim>::get_fe () const
{
return *fe;
}
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
-const Mapping<dim,spacedim> &
+const Mapping<dim,spacedim> &
FEValuesBase<dim,spacedim>::get_mapping () const
{
return *mapping;
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
UpdateFlags
FEValuesBase<dim,spacedim>::get_update_flags () const
-template <int dim, int spacedim>
+template <int dim, int spacedim>
const typename Triangulation<dim,spacedim>::cell_iterator
FEValuesBase<dim,spacedim>::get_cell () const
{
}
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const std::vector<Point<spacedim> > &
FEValuesBase<dim,spacedim>::get_quadrature_points () const
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const std::vector<double> &
FEValuesBase<dim,spacedim>::get_JxW_values () const
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const std::vector<Tensor<2,spacedim> >&
FEValuesBase<dim,spacedim>::get_jacobians () const
{
Assert (this->update_flags & update_quadrature_points, ExcAccessToUninitializedField());
Assert (i<this->quadrature_points.size(), ExcIndexRange(i, 0, this->quadrature_points.size()));
-
+
return this->quadrature_points[i];
}
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
double
FEValuesBase<dim,spacedim>::JxW (const unsigned int i) const
{
Assert (this->update_flags & update_JxW_values, ExcAccessToUninitializedField());
Assert (i<this->JxW_values.size(), ExcIndexRange(i, 0, this->JxW_values.size()));
-
+
return this->JxW_values[i];
}
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const Tensor<2,spacedim> &
FEValuesBase<dim,spacedim>::jacobian (const unsigned int i) const
{
Assert (this->update_flags & update_jacobians, ExcAccessToUninitializedField());
Assert (i<this->jacobians.size(), ExcIndexRange(i, 0, this->jacobians.size()));
-
+
return this->jacobians[i];
}
{
Assert (this->update_flags & update_jacobian_grads, ExcAccessToUninitializedField());
Assert (i<this->jacobian_grads.size(), ExcIndexRange(i, 0, this->jacobian_grads.size()));
-
+
return this->jacobian_grads[i];
}
{
Assert (this->update_flags & update_inverse_jacobians, ExcAccessToUninitializedField());
Assert (i<this->inverse_jacobians.size(), ExcIndexRange(i, 0, this->inverse_jacobians.size()));
-
+
return this->inverse_jacobians[i];
}
-template <int dim, int spacedim>
+template <int dim, int spacedim>
template <class InputVector>
inline
void
-template <int dim, int spacedim>
+template <int dim, int spacedim>
template <class InputVector>
inline
void
-template <int dim, int spacedim>
+template <int dim, int spacedim>
template <class InputVector>
inline
void
-template <int dim, int spacedim>
+template <int dim, int spacedim>
template <class InputVector>
inline
void
-template <int dim, int spacedim>
+template <int dim, int spacedim>
template <class InputVector>
inline
void
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const Point<spacedim> &
FEValuesBase<dim,spacedim>::cell_normal_vector (const unsigned int i) const
ExcIndexRange(i, 0, this->cell_normal_vectors.size()));
Assert (this->update_flags & update_cell_normal_vectors,
typename FVB::ExcAccessToUninitializedField());
-
+
return this->cell_normal_vectors[i];
}
/*------------------------ Inline functions: FEValues ----------------------------*/
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const Quadrature<dim> &
-FEValues<dim,spacedim>::get_quadrature () const
+FEValues<dim,spacedim>::get_quadrature () const
{
return quadrature;
}
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const FEValues<dim,spacedim> &
-FEValues<dim,spacedim>::get_present_fe_values () const
+FEValues<dim,spacedim>::get_present_fe_values () const
{
return *this;
}
/*------------------------ Inline functions: FEFaceValuesBase --------------------*/
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const Point<dim> &
FEFaceValuesBase<dim,spacedim>::normal_vector (const unsigned int i) const
ExcIndexRange(i, 0, this->normal_vectors.size()));
Assert (this->update_flags & update_normal_vectors,
typename FVB::ExcAccessToUninitializedField());
-
+
return this->normal_vectors[i];
}
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
unsigned int
FEFaceValuesBase<dim,spacedim>::get_face_index () const
/*------------------------ Inline functions: FE*FaceValues --------------------*/
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const Quadrature<dim-1> &
-FEFaceValuesBase<dim,spacedim>::get_quadrature () const
+FEFaceValuesBase<dim,spacedim>::get_quadrature () const
{
return quadrature;
}
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const FEFaceValues<dim,spacedim> &
-FEFaceValues<dim,spacedim>::get_present_fe_values () const
+FEFaceValues<dim,spacedim>::get_present_fe_values () const
{
return *this;
}
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const FESubfaceValues<dim,spacedim> &
-FESubfaceValues<dim,spacedim>::get_present_fe_values () const
+FESubfaceValues<dim,spacedim>::get_present_fe_values () const
{
return *this;
}
-template <int dim, int spacedim>
+template <int dim, int spacedim>
inline
const Tensor<1,spacedim> &
FEFaceValuesBase<dim,spacedim>::boundary_form (const unsigned int i) const
ExcIndexRange(i, 0, this->boundary_forms.size()));
Assert (this->update_flags & update_boundary_forms,
typename FVB::ExcAccessToUninitializedField());
-
+
return this->boundary_forms[i];
}
namespace FEValuesViews
{
- template <int dim, int spacedim>
+ template <int dim, int spacedim>
Scalar<dim,spacedim>::Scalar (const FEValuesBase<dim,spacedim> &fe_values,
const unsigned int component)
:
const std::vector<unsigned int> shape_function_to_row_table
= make_shape_function_to_row_table (*fe_values.fe);
-
+
for (unsigned int i=0; i<fe_values.fe->dofs_per_cell; ++i)
{
const bool is_primitive = (fe_values.fe->is_primitive() ||
Assert (false, ExcInternalError());
return *this;
}
-
- template <int dim, int spacedim>
+
+ template <int dim, int spacedim>
Vector<dim,spacedim>::Vector (const FEValuesBase<dim,spacedim> &fe_values,
const unsigned int first_vector_component)
:
const std::vector<unsigned int> shape_function_to_row_table
= make_shape_function_to_row_table (*fe_values.fe);
-
+
for (unsigned int d=0; d<dim; ++d)
{
const unsigned int component = first_vector_component + d;
-
+
for (unsigned int i=0; i<fe_values.fe->dofs_per_cell; ++i)
{
const bool is_primitive = (fe_values.fe->is_primitive() ||
}
}
}
-
+
template <int dim, int spacedim>
Vector<dim,spacedim>::Vector ()
}
+ template <int dim, int spacedim>
+ SymmetricTensor<2, dim, spacedim>::
+ SymmetricTensor(const FEValuesBase<dim, spacedim> &fe_values,
+ const unsigned int first_tensor_component)
+ :
+ fe_values(fe_values),
+ first_tensor_component(first_tensor_component),
+ shape_function_data(fe_values.fe->dofs_per_cell)
+ {
+
+ Assert(first_tensor_component + n_independent_components - 1 < fe_values.fe->n_components(),
+ ExcIndexRange(first_tensor_component + n_independent_components - 1, 0,
+ fe_values.fe->n_components()));
+
+ const std::vector<unsigned int> shape_function_to_row_table
+ = make_shape_function_to_row_table(*fe_values.fe);
+
+ for (unsigned int d = 0; d < n_independent_components; ++d) {
+ const unsigned int component = first_tensor_component + d;
+
+ for (unsigned int i = 0; i < fe_values.fe->dofs_per_cell; ++i) {
+ const bool is_primitive = (fe_values.fe->is_primitive() ||
+ fe_values.fe->is_primitive(i));
+
+ if (is_primitive == true)
+ shape_function_data[i].is_nonzero_shape_function_component[d]
+ = (component ==
+ fe_values.fe->system_to_component_index(i).first);
+ else
+ shape_function_data[i].is_nonzero_shape_function_component[d]
+ = (fe_values.fe->get_nonzero_components(i)[component]
+ == true);
+
+ if (shape_function_data[i].is_nonzero_shape_function_component[d]
+ == true) {
+ if (is_primitive == true)
+ shape_function_data[i].row_index[d]
+ = shape_function_to_row_table[i];
+ else
+ shape_function_data[i].row_index[d]
+ = (shape_function_to_row_table[i]
+ +
+ std::count(fe_values.fe->get_nonzero_components(i).begin(),
+ fe_values.fe->get_nonzero_components(i).begin() +
+ component,
+ true));
+ } else
+ shape_function_data[i].row_index[d]
+ = numbers::invalid_unsigned_int;
+ }
+ }
+
+ for (unsigned int i = 0; i < fe_values.fe->dofs_per_cell; ++i) {
+ unsigned int n_nonzero_components = 0;
+ for (unsigned int d = 0; d < n_independent_components; ++d)
+ if (shape_function_data[i].is_nonzero_shape_function_component[d]
+ == true)
+ ++n_nonzero_components;
+
+ if (n_nonzero_components == 0)
+ shape_function_data[i].single_nonzero_component = -2;
+ else if (n_nonzero_components > 1)
+ shape_function_data[i].single_nonzero_component = -1;
+ else {
+ for (unsigned int d = 0; d < n_independent_components; ++d)
+ if (shape_function_data[i].is_nonzero_shape_function_component[d]
+ == true) {
+ shape_function_data[i].single_nonzero_component
+ = shape_function_data[i].row_index[d];
+ shape_function_data[i].single_nonzero_component_index
+ = d;
+ break;
+ }
+ }
+ }
+
+ switch(dim) {
+ case(1):
+ vector_to_symmetric_tensor_data[0][0] = 0;
+ break;
+ case(2):
+ vector_to_symmetric_tensor_data[0][0] = 0;
+ vector_to_symmetric_tensor_data[1][1] = 1;
+ vector_to_symmetric_tensor_data[1][2] = 2;
+ vector_to_symmetric_tensor_data[2][1] = 2;
+ break;
+ case(3):
+ vector_to_symmetric_tensor_data[0][0] = 0;
+ vector_to_symmetric_tensor_data[1][1] = 1;
+ vector_to_symmetric_tensor_data[2][2] = 2;
+ vector_to_symmetric_tensor_data[0][1] = 3;
+ vector_to_symmetric_tensor_data[1][0] = 3;
+ vector_to_symmetric_tensor_data[0][2] = 4;
+ vector_to_symmetric_tensor_data[2][0] = 4;
+ vector_to_symmetric_tensor_data[1][2] = 5;
+ vector_to_symmetric_tensor_data[2][1] = 5;
+ break;
+ default:
+ ;
+ }
+
+
+
+ }
+
+ template <int dim, int spacedim>
+ SymmetricTensor<2, dim, spacedim>::SymmetricTensor()
+ :
+ fe_values(*static_cast<dealii::FEValuesBase<dim, spacedim>*> (0)),
+ first_tensor_component(numbers::invalid_unsigned_int)
+ {}
+
+
+ template <int dim, int spacedim>
+ SymmetricTensor<2, dim, spacedim> &
+ SymmetricTensor<2, dim, spacedim>::operator=(const SymmetricTensor<2, dim, spacedim> &)
+ {
+ // we shouldn't be copying these objects
+ Assert(false, ExcInternalError());
+ return *this;
+ }
+
+
template <int dim, int spacedim>
template <class InputVector>
void
Scalar<dim,spacedim>::
- get_function_values (const InputVector &fe_function,
+ 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());
+ 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,
if (value == 0.)
continue;
- const double * shape_value_ptr =
+ 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 <int dim, int spacedim>
template <class InputVector>
void
Scalar<dim,spacedim>::
- get_function_gradients (const InputVector &fe_function,
+ 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());
+ 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,
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)
if (value == 0.)
continue;
- const Tensor<1,spacedim> * shape_gradient_ptr =
+ 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)
- template <int dim, int spacedim>
+ template <int dim, int spacedim>
template <class InputVector>
void
Scalar<dim,spacedim>::
- get_function_hessians (const InputVector &fe_function,
+ get_function_hessians (const InputVector &fe_function,
std::vector<hessian_type> &hessians) const
{
typedef FEValuesBase<dim,spacedim> FVB;
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);
if (value == 0.)
continue;
- const Tensor<2,spacedim> * shape_hessian_ptr =
+ 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)
template <class InputVector>
void
Scalar<dim,spacedim>::
- get_function_laplacians (const InputVector &fe_function,
+ 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());
+ 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,
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] +=
+ laplacians[q_point] +=
value * trace(fe_values.shape_hessians[row_index][q_point]);
}
}
template <class InputVector>
void
Vector<dim,spacedim>::
- get_function_values (const InputVector &fe_function,
+ 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());
+ 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,
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
}
}
}
-
- template <int dim, int spacedim>
+
+ template <int dim, int spacedim>
template <class InputVector>
void
Vector<dim,spacedim>::
- get_function_gradients (const InputVector &fe_function,
+ 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());
+ 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,
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
if (snc != -1)
{
- const unsigned int comp =
+ const unsigned int comp =
shape_function_data[shape_function].single_nonzero_component_index;
- const Tensor<1,spacedim> * shape_gradient_ptr =
+ 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)
+ 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 =
+ 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)
- template <int dim, int spacedim>
+ template <int dim, int spacedim>
template <class InputVector>
void
Vector<dim,spacedim>::
- get_function_symmetric_gradients (const InputVector &fe_function,
+ 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());
+ 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,
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(),
+ 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
if (snc != -1)
{
- const unsigned int comp =
+ const unsigned int comp =
shape_function_data[shape_function].single_nonzero_component_index;
- const Tensor<1,spacedim> * shape_gradient_ptr =
+ 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]
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)
+ 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];
- template <int dim, int spacedim>
+ template <int dim, int spacedim>
template <class InputVector>
void
Vector<dim,spacedim>::
- get_function_divergences (const InputVector &fe_function,
+ 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());
+ 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,
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
if (snc != -1)
{
- const unsigned int comp =
+ 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];
divergences[q_point] += value * (*shape_gradient_ptr++)[comp];
}
else
- for (unsigned int d=0; d<dim; ++d)
+ 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 =
}
}
}
-
- template <int dim, int spacedim>
+
+ template <int dim, int spacedim>
template <class InputVector>
void
Vector<dim,spacedim>::
- get_function_hessians (const InputVector &fe_function,
+ get_function_hessians (const InputVector &fe_function,
std::vector<hessian_type> &hessians) const
{
typedef FEValuesBase<dim,spacedim> FVB;
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);
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
if (snc != -1)
{
- const unsigned int comp =
+ const unsigned int comp =
shape_function_data[shape_function].single_nonzero_component_index;
- const Tensor<2,spacedim> * shape_hessian_ptr =
+ 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++;
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 =
+ 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)
template <class InputVector>
void
Vector<dim,spacedim>::
- get_function_laplacians (const InputVector &fe_function,
+ 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());
+ 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,
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
{
const unsigned int comp =
shape_function_data[shape_function].single_nonzero_component_index;
- const Tensor<2,spacedim> * shape_hessian_ptr =
+ 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++);
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 =
+ 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)
}
}
}
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ SymmetricTensor<2, 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());
+
+ // the unique components of the
+ // second order tensor stored as
+ // a vector (i.e. a first-order
+ // tensor)
+ typedef Tensor<1, n_independent_components> base_tensor_type;
+
+ std::vector< base_tensor_type > values_in_vector_form(values.size(), base_tensor_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 (value == 0.)
+ continue;
+
+ 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_in_vector_form[q_point][comp] += value * *shape_value_ptr++;
+ } else
+ for (unsigned int d = 0; d < n_independent_components; ++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_in_vector_form[q_point][d] += value * *shape_value_ptr++;
+ }
+ }
+ // copy entries in std::vector to an array as there is no constructor
+ // for a second order tensor that take a std::vector
+ double values_array[n_independent_components];
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+ for (unsigned int d = 0; d < n_independent_components; d++)
+ values_array[d] = values_in_vector_form[q_point][d];
+
+ values[q_point] = dealii::SymmetricTensor<2, dim>(values_array);
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ SymmetricTensor<2, 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 (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];
+ for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
+ ++q_point, ++shape_gradient_ptr) {
+ for (unsigned int j = 0; j < dim; ++j) {
+ const unsigned int vector_component = vector_to_symmetric_tensor_data[comp][j];
+ divergences[q_point][vector_component] += value * (*shape_gradient_ptr)[j];
+ }
+ }
+ } else
+ for (unsigned int d = 0; d < n_independent_components; ++d)
+ if (shape_function_data[shape_function].is_nonzero_shape_function_component[d]) {
+ const unsigned int comp =
+ shape_function_data[shape_function].single_nonzero_component_index;
+
+ 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, ++shape_gradient_ptr) {
+ for (unsigned int j = 0; j < dim; ++j) {
+ const unsigned int vector_component = vector_to_symmetric_tensor_data[comp][j];
+ divergences[q_point][vector_component] += value * (*shape_gradient_ptr++)[j];
+ }
+ }
+ }
+ }
+ }
}
dealii::FEValuesViews::Scalar<dim,spacedim>(fe_values,
component);
}
-
+
const unsigned int n_vectors = (fe.n_components() >= dim ?
fe.n_components()-dim+1 :
0);
new (&vectors[component])
dealii::FEValuesViews::Vector<dim,spacedim>(fe_values,
component);
- }
+ }
+
+ const unsigned int n_symmetric_second_order_tensors = (fe.n_components() >= 0.5*(dim*dim + dim) ?
+ fe.n_components() - 0.5*(dim*dim + dim) + 1 :
+ 0);
+ symmetric_second_order_tensors.resize(n_symmetric_second_order_tensors);
+ for (unsigned int component = 0; component < n_symmetric_second_order_tensors; ++component)
+ {
+ typedef dealii::FEValuesViews::SymmetricTensor<2, dim, spacedim> SymmetricTensorView;
+ symmetric_second_order_tensors[component].SymmetricTensorView::~SymmetricTensorView();
+ new (&symmetric_second_order_tensors[component])
+ dealii::FEValuesViews::SymmetricTensor<2, dim, spacedim > (fe_values,
+ component);
+ }
}
}
}
* interface. See there for a
* description of the use of
* these classes.
- *
+ *
* @author Wolfgang Bangerth, 2003
*/
template <int dim, int spacedim>
* this class.
*/
CellIterator (const CI &cell);
-
+
/**
* Conversion operator to an
* iterator for
*/
virtual
operator const typename Triangulation<dim,spacedim>::cell_iterator () const;
-
+
/**
* Return the number of
* degrees of freedom the DoF
*/
const CI cell;
};
-
+
/**
* Implementation of a derived
* iterator that allows to
* extract degree of freedom
* information.
- *
+ *
* @author Wolfgang Bangerth, 2003
*/
template <int dim, int spacedim>
* this class.
*/
TriaCellIterator (const typename Triangulation<dim,spacedim>::cell_iterator &cell);
-
+
/**
* Conversion operator to an
* iterator for
*/
static const char * const message_string;
};
-
+
// functions
unsigned int n_nonzero_shape_components = 0;
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- n_nonzero_shape_components += fe.n_nonzero_components (i);
+ n_nonzero_shape_components += fe.n_nonzero_components (i);
Assert (n_nonzero_shape_components >= fe.dofs_per_cell,
ExcInternalError());
if (flags & update_hessians)
this->shape_hessians.resize (n_nonzero_shape_components,
std::vector<Tensor<2,spacedim> > (n_quadrature_points));
-
+
if (flags & update_quadrature_points)
this->quadrature_points.resize(n_quadrature_points);
// each quadrature point
Assert (values.size() == n_quadrature_points,
ExcDimensionMismatch(values.size(), n_quadrature_points));
-
+
// initialize with zero
std::fill_n (values.begin(), n_quadrature_points, 0);
-
+
// add up contributions of trial
// functions. note that here we
// deal with scalar finite
const InputVector& fe_function,
std::vector<Vector<number> >& values) const
{
-//TODO: Find out how to do this assertion.
+//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(),
for (unsigned i=0;i<values.size();++i)
Assert (values[i].size() == n_components,
ExcDimensionMismatch(values[i].size(), n_components));
-
+
Assert (this->update_flags & update_values, 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<values.size();++i)
std::fill_n (values[i].begin(), values[i].size(), 0);
// function can't be used because
// it doesn't return us a pointer
// to the data).
- const unsigned int
- row = fe->is_primitive() ?
+ const unsigned int
+ row = fe->is_primitive() ?
shape_func : this->shape_function_to_row_table[shape_func];
const double *shape_value_ptr = &this->shape_values(row, 0);
// 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
// multiple of dofs_per_cell such
// that an integer number of
// number of components of the
// finite element
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));
// 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);
if (fe->is_primitive(shape_func))
{
- const unsigned int
- row = fe->is_primitive() ?
+ const unsigned int
+ row = fe->is_primitive() ?
shape_func : this->shape_function_to_row_table[shape_func];
const double *shape_value_ptr = &this->shape_values(row, 0);
- const unsigned int comp = fe->system_to_component_index(shape_func).first
+ 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)
values[point](comp) += value * *shape_value_ptr++;
const double *shape_value_ptr = &this->shape_values(row, 0);
const unsigned int comp = c + mc * n_components;
-
+
for (unsigned int point=0; point<n_quadrature_points; ++point)
values[point](comp) += value * *shape_value_ptr++;
}
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
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);
if (fe->is_primitive(shape_func))
{
- const unsigned int
- row = fe->is_primitive() ?
+ const unsigned int
+ row = fe->is_primitive() ?
shape_func : this->shape_function_to_row_table[shape_func];
const double *shape_value_ptr = &this->shape_values(row, 0);
- const unsigned int comp = fe->system_to_component_index(shape_func).first
+ const unsigned int comp = fe->system_to_component_index(shape_func).first
+ mc * n_components;
if (quadrature_points_fastest)
const double *shape_value_ptr = &this->shape_values(row, 0);
const unsigned int comp = c + mc * n_components;
-
+
if (quadrature_points_fastest)
for (unsigned int point=0; point<n_quadrature_points; ++point)
values[comp][point] += value * *shape_value_ptr++;
if (value == 0.)
continue;
- const Tensor<1,spacedim> *shape_gradient_ptr
+ 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++;
// each quadrature point
Assert (gradients.size() == n_quadrature_points,
ExcDimensionMismatch(gradients.size(), n_quadrature_points));
-
+
// initialize with zero
std::fill_n (gradients.begin(), n_quadrature_points, Tensor<1,spacedim>());
-
+
// add up contributions of trial
// functions. note that here we
// deal with scalar finite
if (value == 0.)
continue;
- const Tensor<1,spacedim> *shape_gradient_ptr
+ 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++;
if (fe->is_primitive(shape_func))
{
- const unsigned int
- row = fe->is_primitive() ?
+ const unsigned int
+ row = fe->is_primitive() ?
shape_func : this->shape_function_to_row_table[shape_func];
- const Tensor<1,spacedim> *shape_gradient_ptr
+ const Tensor<1,spacedim> *shape_gradient_ptr
= &this->shape_gradients[row][0];
const unsigned int comp = fe->system_to_component_index(shape_func).first;
for (unsigned int point=0; point<n_quadrature_points; ++point)
fe->get_nonzero_components(shape_func).begin()+c,
true));
- const Tensor<1,spacedim> *shape_gradient_ptr
+ const Tensor<1,spacedim> *shape_gradient_ptr
= &this->shape_gradients[row][0];
for (unsigned int point=0; point<n_quadrature_points; ++point)
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
// 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)
// 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<gradients.size();++i)
std::fill_n (gradients[i].begin(), gradients[i].size(), Tensor<1,spacedim>());
if (fe->is_primitive(shape_func))
{
- const unsigned int
- row = fe->is_primitive() ?
+ const unsigned int
+ row = fe->is_primitive() ?
shape_func : this->shape_function_to_row_table[shape_func];
- const Tensor<1,spacedim> *shape_gradient_ptr
+ const Tensor<1,spacedim> *shape_gradient_ptr
= &this->shape_gradients[row][0];
const unsigned int comp = fe->system_to_component_index(shape_func).first
+ mc * n_components;
fe->get_nonzero_components(shape_func).begin()+c,
true));
- const Tensor<1,spacedim> *shape_gradient_ptr
+ const Tensor<1,spacedim> *shape_gradient_ptr
= &this->shape_gradients[row][0];
const unsigned int comp = c + mc * n_components;
if (value == 0.)
continue;
- const Tensor<2,spacedim> *shape_hessians_ptr
+ 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++;
// each quadrature point
Assert (hessians.size() == n_quadrature_points,
ExcDimensionMismatch(hessians.size(), n_quadrature_points));
-
+
// initialize with zero
std::fill_n (hessians.begin(), n_quadrature_points, Tensor<2,spacedim>());
-
+
// add up contributions of trial
// functions. note that here we
// deal with scalar finite
if (value == 0.)
continue;
- const Tensor<2,spacedim> *shape_hessians_ptr
+ 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++;
if (fe->is_primitive(shape_func))
{
- const unsigned int
- row = fe->is_primitive() ?
+ const unsigned int
+ row = fe->is_primitive() ?
shape_func : this->shape_function_to_row_table[shape_func];
- const Tensor<2,spacedim> *shape_hessian_ptr
+ const Tensor<2,spacedim> *shape_hessian_ptr
= &this->shape_hessians[row][0];
const unsigned int comp = fe->system_to_component_index(shape_func).first;
fe->get_nonzero_components(shape_func).begin()+c,
true));
- const Tensor<2,spacedim> *shape_hessian_ptr
+ const Tensor<2,spacedim> *shape_hessian_ptr
= &this->shape_hessians[row][0];
if (quadrature_points_fastest)
bool quadrature_points_fastest) const
{
Assert (this->update_flags & update_second_derivatives, ExcAccessToUninitializedField());
-
+
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
// 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)
// element, we need this number for
// loops filling all components
const unsigned int component_multiple = result_components / n_components;
-
+
// initialize with zero
for (unsigned i=0;i<hessians.size();++i)
std::fill_n (hessians[i].begin(), hessians[i].size(), Tensor<2,spacedim>());
if (fe->is_primitive(shape_func))
{
- const unsigned int
- row = fe->is_primitive() ?
+ const unsigned int
+ row = fe->is_primitive() ?
shape_func : this->shape_function_to_row_table[shape_func];
- const Tensor<2,spacedim> *shape_hessian_ptr
+ 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;
fe->get_nonzero_components(shape_func).begin()+c,
true));
- const Tensor<2,spacedim> *shape_hessian_ptr
+ const Tensor<2,spacedim> *shape_hessian_ptr
= &this->shape_hessians[row][0];
const unsigned int comp = c + mc * n_components;
if (value == 0.)
continue;
- const Tensor<2,spacedim> *shape_hessian_ptr
+ 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++);
// 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
if (value == 0.)
continue;
- const Tensor<2,spacedim> *shape_hessian_ptr
+ 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++);
const InputVector& fe_function,
std::vector<Vector<number> >& laplacians) const
{
-//TODO: Find out how to do this assertion.
+//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(),
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);
if (fe->is_primitive(shape_func))
{
- const unsigned int
- row = fe->is_primitive() ?
+ const unsigned int
+ row = fe->is_primitive() ?
shape_func : this->shape_function_to_row_table[shape_func];
- const Tensor<2,spacedim> *shape_hessian_ptr
+ 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)
fe->get_nonzero_components(shape_func).begin()+c,
true));
- const Tensor<2,spacedim> *shape_hessian_ptr
+ const Tensor<2,spacedim> *shape_hessian_ptr
= &this->shape_hessians[row][0];
for (unsigned int point=0; point<n_quadrature_points; ++point)
// 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
// 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));
// 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);
if (fe->is_primitive(shape_func))
{
- const unsigned int
- row = fe->is_primitive() ?
+ 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
+ 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++);
fe->get_nonzero_components(shape_func).begin()+c,
true));
- const Tensor<2,spacedim> *shape_hessian_ptr
+ 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++);
}
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
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);
if (fe->is_primitive(shape_func))
{
- const unsigned int
- row = fe->is_primitive() ?
+ const unsigned int
+ row = fe->is_primitive() ?
shape_func : this->shape_function_to_row_table[shape_func];
- const Tensor<2,spacedim> *shape_hessian_ptr
+ const Tensor<2,spacedim> *shape_hessian_ptr
= &this->shape_hessians[row][0];
- const unsigned int comp = fe->system_to_component_index(shape_func).first
+ 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)
fe->get_nonzero_components(shape_func).begin()+c,
true));
- const Tensor<2,spacedim> *shape_hessian_ptr
+ 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++);
UpdateFlags
FEValuesBase<dim,spacedim>::compute_update_flags (const UpdateFlags update_flags) const
{
-
+
// first find out which objects
// need to be recomputed on each
// cell we visit. this we have to
template <int dim, int spacedim>
inline
void
-FEValuesBase<dim,spacedim>::check_cell_similarity
+FEValuesBase<dim,spacedim>::check_cell_similarity
(const typename Triangulation<dim,spacedim>::cell_iterator &cell)
{
// case that there has not been any cell
else
cell_similarity = (cell->is_translation_of
(static_cast<const typename Triangulation<dim,spacedim>::cell_iterator>(*this->present_cell))
- ?
- CellSimilarity::translation
- :
+ ?
+ CellSimilarity::translation
+ :
CellSimilarity::none);
-
+
// TODO: here, one could implement other
// checks for similarity, e.g. for
// children of a parallelogram.
template <int dim, int spacedim>
-CellSimilarity::Similarity
+CellSimilarity::Similarity
FEValuesBase<dim,spacedim>::get_cell_similarity () const
{
return cell_similarity;
const UpdateFlags flags = this->compute_update_flags (update_flags);
-
+
// then get objects into which the
// FE and the Mapping can store
// intermediate data used across
// destroyed and also that this
// object gets destroyed in the
// destruction of this class
- this->present_cell.reset
+ this->present_cell.reset
(new typename FEValuesBase<dim,spacedim>::template
CellIterator<typename DoFHandler<dim,spacedim>::cell_iterator> (cell));
// destroyed and also that this
// object gets destroyed in the
// destruction of this class
- this->present_cell.reset
+ this->present_cell.reset
(new typename FEValuesBase<dim,spacedim>::template
CellIterator<typename hp::DoFHandler<dim,spacedim>::cell_iterator> (cell));
// this cell, are the same
//TODO: This was documented out ith the repository. Why?
-
+
// Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
// static_cast<const FiniteElementData<dim>&>(cell->get_fe()),
// typename FEValuesBase<dim,spacedim>::ExcFEDontMatch());
// destroyed and also that this
// object gets destroyed in the
// destruction of this class
- this->present_cell.reset
+ this->present_cell.reset
(new typename FEValuesBase<dim,spacedim>::template
CellIterator<typename MGDoFHandler<dim,spacedim>::cell_iterator> (cell));
// destroyed and also that this
// object gets destroyed in the
// destruction of this class
- this->present_cell.reset
+ this->present_cell.reset
(new typename FEValuesBase<dim,spacedim>::TriaCellIterator (cell));
// this was the part of the work
// that is dependent on the actual
FEFaceValuesBase<dim,spacedim>::FEFaceValuesBase (const unsigned int n_q_points,
const unsigned int dofs_per_cell,
const UpdateFlags,
- const Mapping<dim,spacedim> &mapping,
+ const Mapping<dim,spacedim> &mapping,
const FiniteElement<dim,spacedim> &fe,
const Quadrature<dim-1>& quadrature)
:
FEFaceValues<dim,spacedim>::initialize (const UpdateFlags update_flags)
{
const UpdateFlags flags = this->compute_update_flags (update_flags);
-
+
// then get objects into which the
// FE and the Mapping can store
// intermediate data used across
Assert (face_no < GeometryInfo<dim>::faces_per_cell,
ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
-
+
// set new cell. auto_ptr will take
// care that old object gets
// destroyed and also that this
// object gets destroyed in the
// destruction of this class
- this->present_cell.reset
+ this->present_cell.reset
(new typename FEValuesBase<dim,spacedim>::template
CellIterator<typename DoFHandler<dim,spacedim>::cell_iterator> (cell));
-
+
// this was the part of the work
// that is dependent on the actual
// data type of the iterator. now
Assert (face_no < GeometryInfo<dim>::faces_per_cell,
ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
-
+
// set new cell. auto_ptr will take
// care that old object gets
// destroyed and also that this
// object gets destroyed in the
// destruction of this class
- this->present_cell.reset
+ this->present_cell.reset
(new typename FEValuesBase<dim,spacedim>::template
CellIterator<typename hp::DoFHandler<dim,spacedim>::cell_iterator> (cell));
Assert (face_no < GeometryInfo<dim>::faces_per_cell,
ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
-
+
// set new cell. auto_ptr will take
// care that old object gets
// destroyed and also that this
// object gets destroyed in the
// destruction of this class
- this->present_cell.reset
+ this->present_cell.reset
(new typename FEValuesBase<dim,spacedim>::template
CellIterator<typename MGDoFHandler<dim,spacedim>::cell_iterator> (cell));
// destroyed and also that this
// object gets destroyed in the
// destruction of this class
- this->present_cell.reset
+ this->present_cell.reset
(new typename FEValuesBase<dim,spacedim>::TriaCellIterator (cell));
// this was the part of the work
// available)
const typename Triangulation<dim,spacedim>::cell_iterator cell=*this->present_cell;
this->present_face_index=cell->face_index(face_no);
-
+
this->get_mapping().fill_fe_face_values(*this->present_cell, face_no,
this->quadrature,
*this->mapping_data,
FESubfaceValues<dim,spacedim>::initialize (const UpdateFlags update_flags)
{
const UpdateFlags flags = this->compute_update_flags (update_flags);
-
+
// then get objects into which the
// FE and the Mapping can store
// intermediate data used across
// destroyed and also that this
// object gets destroyed in the
// destruction of this class
- this->present_cell.reset
+ this->present_cell.reset
(new typename FEValuesBase<dim,spacedim>::template
CellIterator<typename DoFHandler<dim,spacedim>::cell_iterator> (cell));
do_reinit (face_no, subface_no);
}
-
+
template <int dim, int spacedim>
void FESubfaceValues<dim,spacedim>::reinit (const typename MGDoFHandler<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
// destroyed and also that this
// object gets destroyed in the
// destruction of this class
- this->present_cell.reset
+ this->present_cell.reset
(new typename FEValuesBase<dim,spacedim>::template
CellIterator<typename MGDoFHandler<dim,spacedim>::cell_iterator> (cell));
ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
Assert (subface_no < cell->face(face_no)->n_children(),
ExcIndexRange (subface_no, 0, cell->face(face_no)->n_children()));
-
+
// set new cell. auto_ptr will take
// care that old object gets
// destroyed and also that this
// object gets destroyed in the
// destruction of this class
- this->present_cell.reset
+ this->present_cell.reset
(new typename FEValuesBase<dim,spacedim>::TriaCellIterator (cell));
// this was the part of the work
ExcInternalError());
this->present_face_index=subface_index;
}
-
+
// now ask the mapping and the finite element
// to do the actual work
this->get_mapping().fill_fe_subface_values(*this->present_cell,
this->JxW_values,
this->boundary_forms,
this->normal_vectors);
-
+
this->get_fe().fill_fe_subface_values(this->get_mapping(),
*this->present_cell,
face_no, subface_no,
template class FEValuesData<deal_II_dimension,deal_II_dimension+1>;
template class FEValuesBase<deal_II_dimension,deal_II_dimension+1>;
template class FEValues<deal_II_dimension,deal_II_dimension+1>;
-template class FEValuesBase<deal_II_dimension,deal_II_dimension+1>::
-CellIterator<DoFHandler<deal_II_dimension,deal_II_dimension+1>::cell_iterator>;
+template class FEValuesBase<deal_II_dimension,deal_II_dimension+1>::
+CellIterator<DoFHandler<deal_II_dimension,deal_II_dimension+1>::cell_iterator>;
//template class FEValuesBase<deal_II_dimension,deal_II_dimension+1>::
// CellIterator<MGDoFHandler<deal_II_dimension,deal_II_dimension+1>::cell_iterator>;