From a160dcc79962b1a930edf8204035d7930e07e5bf Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 23 Sep 2009 17:48:36 +0000 Subject: [PATCH] Patch by Andrew McBride: implement extractors for symmetric rank 2 tensors. git-svn-id: https://svn.dealii.org/trunk@19506 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_values.h | 671 ++++++++++++++----- deal.II/deal.II/source/fe/fe_values.cc | 681 ++++++++++++++------ deal.II/deal.II/source/fe/fe_values.inst.in | 22 +- 3 files changed, 1001 insertions(+), 373 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 27fbd4d625..172e698a69 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -55,8 +55,10 @@ template class FEValuesBase; * 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 dim elements. + * element. There are extractors for single scalar components, + * vector components consisting of dim elements, + * and second order symmetric tensors consisting + * of (dim*dim + dim)/2 components * * See the description of the @ref vector_valued module for examples how to * use the features of this namespace. @@ -90,7 +92,7 @@ namespace FEValuesExtractors Scalar (const unsigned int component); }; - + /** * Extractor for a vector of * dim components of a @@ -123,10 +125,53 @@ namespace FEValuesExtractors */ 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 + * (dim*dim + dim)/2 + * components of a vector-valued + * element. The value of + * dim 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 + 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 @@ -134,8 +179,10 @@ namespace FEValuesExtractors * 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 dim elements. + * There are classes that present views for single scalar components, + * vector components consisting of dim elements, and + * symmetric second order tensor components + * consisting of (dim*dim + dim)/2 elements * * See the description of the @ref vector_valued module for examples how to * use the features of this namespace. @@ -152,7 +199,7 @@ namespace FEValuesViews * * @ingroup feaccess vector_valued */ - template + template class Scalar { public: @@ -190,7 +237,7 @@ namespace FEValuesViews * invalid object. */ Scalar (); - + /** * Constructor for an object that * represents a single scalar component @@ -208,7 +255,7 @@ namespace FEValuesViews * this function is called. */ Scalar & operator= (const Scalar &); - + /** * Return the value of the vector * component selected by this view, for @@ -330,12 +377,12 @@ namespace FEValuesViews template void get_function_laplacians (const InputVector& fe_function, std::vector& laplacians) const; - - + + private: /** * A reference to the FEValuesBase object - * we operator on. + * we operate on. */ const FEValuesBase &fe_values; @@ -394,7 +441,7 @@ namespace FEValuesViews */ std::vector shape_function_data; }; - + /** @@ -406,7 +453,7 @@ namespace FEValuesViews * * @ingroup feaccess vector_valued */ - template + template class Vector { public: @@ -438,7 +485,7 @@ namespace FEValuesViews * symmetrized gradient is a * SymmetricTensor@<2,spacedim@>. */ - typedef SymmetricTensor<2,spacedim> symmetric_gradient_type; + typedef dealii::SymmetricTensor<2,spacedim> symmetric_gradient_type; /** * A typedef for the type of the @@ -449,7 +496,7 @@ namespace FEValuesViews * course is a scalar. */ typedef double divergence_type; - + /** * A typedef for the type of second * derivatives of the view this class @@ -625,7 +672,7 @@ namespace FEValuesViews template void get_function_symmetric_gradients (const InputVector& fe_function, std::vector& symmetric_gradients) const; - + /** * Return the divergence of the selected * vector components of the finite @@ -648,7 +695,7 @@ namespace FEValuesViews template void get_function_divergences (const InputVector& fe_function, std::vector& divergences) const; - + /** * Return the Hessians of the selected * vector components of the finite @@ -695,11 +742,11 @@ namespace FEValuesViews template void get_function_laplacians (const InputVector& fe_function, std::vector& laplacians) const; - + private: /** * A reference to the FEValuesBase object - * we operator on. + * we operate on. */ const FEValuesBase &fe_values; @@ -717,7 +764,7 @@ namespace FEValuesViews * much cheaper. */ struct ShapeFunctionData - { + { /** * For each pair (shape * function,component within @@ -780,6 +827,251 @@ namespace FEValuesViews */ std::vector shape_function_data; }; + + + template + class SymmetricTensor; + + /** + * A class representing a view to a set of + * (dim*dim + dim)/2 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 + 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 (dim*dim + dim)/2 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 (dim*dim + dim)/2 unique components + * of the finite element representing a symmetric second-order + * tensor, the divergence of + * course is a * Tensor@<1,dim@>. + */ + typedef Tensor<1, spacedim> divergence_type; + + /** + * Default constructor. Creates an + * invalid object. + */ + SymmetricTensor(); + + /** + * Constructor for an object that + * represents (dim*dim + dim)/2 + * 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 &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 + * (dim*dim + dim)/2 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 + * fe_function at the + * quadrature points of the cell, face + * or subface selected the last time + * the reinit 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 + void get_function_values(const InputVector& fe_function, + std::vector& values) const; + + /** + * Return the divergence of the selected + * vector components of the finite + * element function characterized by + * fe_function at the + * quadrature points of the cell, face + * or subface selected the last time + * the reinit 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 + void get_function_divergences(const InputVector& fe_function, + std::vector& divergences) const; + + private: + /** + * A reference to the FEValuesBase object + * we operate on. + */ + const FEValuesBase &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 > 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 shape_function_data; + }; } @@ -804,10 +1096,13 @@ namespace internal { /** * Caches for scalar and - * vector-valued views. + * vector, and symmetric second-order + * tensor valued views. */ std::vector > scalars; std::vector > vectors; + std::vector > + symmetric_second_order_tensors; /** * Constructor. @@ -819,7 +1114,7 @@ namespace internal - + /*!@addtogroup feaccess */ /*@{*/ @@ -839,7 +1134,7 @@ namespace internal * * The purpose of this class is discussed * on the page on @ref UpdateFlagsEssay. - * + * * @author Guido Kanschat, 2000 */ template @@ -979,7 +1274,7 @@ class FEValuesData std::vector > jacobian_grads; /** - * Array of the inverse Jacobian matrices + * Array of the inverse Jacobian matrices * at the quadrature points. */ std::vector > inverse_jacobians; @@ -1001,7 +1296,7 @@ class FEValuesData std::vector > 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). @@ -1027,14 +1322,14 @@ class FEValuesData * 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 @@ -1065,7 +1360,7 @@ class FEValuesData * non-zero components. */ std::vector shape_function_to_row_table; - + /** * Original update flags handed * to the constructor of @@ -1139,7 +1434,7 @@ class FEValuesData * * The functions of this class fall into different cathegories: *
    - *
  • shape_value(), shape_grad(), etc: return one of the values + *
  • 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 @@ -1158,13 +1453,13 @@ class FEValuesData * 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. - * + * *
  • get_function_values(), get_function_gradients(), etc.: Compute a * finite element function or its derivative in quadrature points. * *
  • 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. *
* @@ -1183,7 +1478,7 @@ class FEValuesData * give all 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 */ @@ -1208,7 +1503,7 @@ class FEValuesBase : protected FEValuesData, */ const unsigned int dofs_per_cell; - + /** * Constructor. Set up the array * sizes with n_q_points @@ -1264,12 +1559,26 @@ class FEValuesBase : protected FEValuesData, const FEValuesViews::Vector & operator[] (const FEValuesExtractors::Vector &vector) const; + /** + * Create a view of the current FEValues + * object that represents a set of + * (dim*dim + dim)/2 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, @@ -1470,7 +1779,7 @@ class FEValuesBase : protected FEValuesData, shape_hessian_component (const unsigned int function_no, const unsigned int point_no, const unsigned int component) const; - + /** * @deprecated Wrapper for shape_hessian_component() */ @@ -1478,12 +1787,12 @@ class FEValuesBase : protected FEValuesData, 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 @@ -1512,10 +1821,10 @@ class FEValuesBase : protected FEValuesData, * it is a vector-valued one, then use * the other get_function_values() * function. - * + * * The function assumes that the * values object already has the - * correct size. + * correct size. * * The actual data type of the * input vector may be either a @@ -1564,7 +1873,7 @@ class FEValuesBase : protected FEValuesData, /** * Generate function values from * an arbitrary vector. - * + * * This function offers the * possibility to extract * function values in quadrature @@ -1721,7 +2030,7 @@ class FEValuesBase : protected FEValuesData, * 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. @@ -1799,7 +2108,7 @@ class FEValuesBase : protected FEValuesData, void get_function_gradients (const InputVector &fe_function, std::vector > > &gradients) const; - + /** * @deprecated Use * get_function_gradients() instead. @@ -1869,7 +2178,7 @@ class FEValuesBase : protected FEValuesData, * the other * get_function_hessians() * function. - * + * * The actual data type of the * input vector may be either a * Vector<T>, @@ -1895,7 +2204,7 @@ class FEValuesBase : protected FEValuesData, get_function_hessians (const InputVector& fe_function, std::vector >& hessians) const; - + /** * Compute the tensor of second * derivatives of the finite @@ -1980,7 +2289,7 @@ class FEValuesBase : protected FEValuesData, get_function_2nd_derivatives (const InputVector&, std::vector > >&, bool = false) const; - + /** * Compute the (scalar) Laplacian @@ -2008,7 +2317,7 @@ class FEValuesBase : protected FEValuesData, * the other * get_function_laplacians() * function. - * + * * The actual data type of the * input vector may be either a * Vector<T>, @@ -2035,7 +2344,7 @@ class FEValuesBase : protected FEValuesData, get_function_laplacians (const InputVector& fe_function, std::vector& laplacians) const; - + /** * Compute the (scalar) Laplacian * of the finite element function @@ -2128,7 +2437,7 @@ class FEValuesBase : protected FEValuesData, //@} - + /** * Position of the ith * quadrature point in real space. @@ -2171,7 +2480,7 @@ class FEValuesBase : protected FEValuesData, * the values returned by JxW(). */ const std::vector & get_JxW_values () const; - + /** * Return the Jacobian of the * transformation at the specified @@ -2212,7 +2521,7 @@ class FEValuesBase : protected FEValuesData, /** * Pointer to the array holding - * the values returned by + * the values returned by * inverse_jacobian(). */ const std::vector > & get_inverse_jacobians () const; @@ -2235,7 +2544,7 @@ class FEValuesBase : protected FEValuesData, * for this object. */ UpdateFlags get_update_flags () const; - + /** * Return a triangulation * iterator to the current cell. @@ -2324,7 +2633,7 @@ class FEValuesBase : protected FEValuesData, * @ingroup Exceptions */ DeclException0 (ExcFENotPrimitive); - + protected: /** * Objects of the FEValues @@ -2385,7 +2694,7 @@ class FEValuesBase : protected FEValuesData, * * @author Wolfgang Bangerth, 2003 */ - class CellIteratorBase + class CellIteratorBase { public: /** @@ -2395,7 +2704,7 @@ class FEValuesBase : protected FEValuesData, * class. */ virtual ~CellIteratorBase (); - + /** * Conversion operator to an * iterator for @@ -2568,7 +2877,7 @@ class FEValuesBase : protected FEValuesData, */ template class CellIterator; class TriaCellIterator; - + /** * Store the cell selected last time the * reinit() function was called. This is @@ -2578,18 +2887,18 @@ class FEValuesBase : protected FEValuesData, * the extractor classes. */ std::auto_ptr present_cell; - + /** * Storage for the mapping object. */ const SmartPointer > mapping; - + /** * Store the finite element for later use. */ const SmartPointer > fe; - + /** * Internal data of mapping. */ @@ -2649,7 +2958,7 @@ class FEValuesBase : protected FEValuesData, * and also do not implement it. */ FEValuesBase (const FEValuesBase &); - + /** * Copy operator. Since * objects of this class are not @@ -2663,7 +2972,7 @@ class FEValuesBase : protected FEValuesData, * objects. */ internal::FEValuesViews::Cache fe_values_views_cache; - + /** * Make the view classes friends of this * class, since they access internal @@ -2671,6 +2980,7 @@ class FEValuesBase : protected FEValuesData, */ template friend class FEValuesViews::Scalar; template friend class FEValuesViews::Vector; + template friend class FEValuesViews::SymmetricTensor; }; @@ -2681,7 +2991,7 @@ class FEValuesBase : protected FEValuesData, * 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 @@ -2707,7 +3017,7 @@ class FEValues : public FEValuesBase * to dim. */ static const unsigned int integral_dimension = dim; - + /** * Constructor. Gets cell * independent data from mapping @@ -2727,7 +3037,7 @@ class FEValues : public FEValuesBase FEValues (const FiniteElement &fe, const Quadrature &quadrature, const UpdateFlags update_flags); - + /** * Reinitialize the gradients, * Jacobi determinants, etc for @@ -2807,7 +3117,7 @@ class FEValues : public FEValuesBase * stored by this object. */ const Quadrature & get_quadrature () const; - + /** * Determine an estimate for the * memory consumption (in bytes) @@ -2839,7 +3149,7 @@ class FEValues : public FEValuesBase * FEValues/hpFEValues). */ const FEValues & get_present_fe_values () const; - + private: /** * Store a copy of the quadrature @@ -2889,7 +3199,7 @@ class FEValues : public FEValuesBase * * @author Wolfgang Bangerth, 1998, Guido Kanschat, 2000, 2001 */ -template +template class FEFaceValuesBase : public FEValuesBase { public: @@ -2933,7 +3243,7 @@ class FEFaceValuesBase : public FEValuesBase * is normalized to one. */ const Point & normal_vector (const unsigned int i) const; - + /** * Boundary form of the * transformation of the cell at @@ -2952,7 +3262,7 @@ class FEFaceValuesBase : public FEValuesBase * n.ds. */ const Tensor<1,spacedim> & boundary_form (const unsigned int i) const; - + /** * Return the list of outward normal * vectors to the cell at the @@ -2981,14 +3291,14 @@ class FEFaceValuesBase : public FEValuesBase * stored by this object. */ const Quadrature & get_quadrature () const; - + /** * Determine an estimate for the * memory consumption (in bytes) * of this object. */ unsigned int memory_consumption () const; - + protected: /** @@ -3020,7 +3330,7 @@ class FEFaceValuesBase : public FEValuesBase * * @author Wolfgang Bangerth, 1998, Guido Kanschat, 2000, 2001 */ -template +template class FEFaceValues : public FEFaceValuesBase { public: @@ -3130,7 +3440,7 @@ class FEFaceValues : public FEFaceValuesBase */ void reinit (const typename Triangulation::cell_iterator &cell, const unsigned int face_no); - + /** * Return a reference to this * very object. @@ -3161,7 +3471,7 @@ class FEFaceValues : public FEFaceValuesBase * Do work common to the two * constructors. */ - void initialize (const UpdateFlags update_flags); + void initialize (const UpdateFlags update_flags); /** * The reinit() functions do @@ -3196,7 +3506,7 @@ class FEFaceValues : public FEFaceValuesBase * * @author Wolfgang Bangerth, 1998, Guido Kanschat, 2000, 2001 */ -template +template class FESubfaceValues : public FEFaceValuesBase { public: @@ -3273,7 +3583,7 @@ class FESubfaceValues : public FEFaceValuesBase void reinit (const typename hp::DoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no); - + /** * Reinitialize the gradients, * Jacobi determinants, etc for @@ -3290,7 +3600,7 @@ class FESubfaceValues : public FEFaceValuesBase void reinit (const typename MGDoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no); - + /** * Reinitialize the gradients, * Jacobi determinants, etc for @@ -3406,6 +3716,13 @@ namespace FEValuesExtractors : first_vector_component (first_vector_component) {} + + template + inline + SymmetricTensor::SymmetricTensor (const unsigned int first_tensor_component) + : + first_tensor_component (first_tensor_component) + {} } @@ -3413,7 +3730,7 @@ namespace FEValuesExtractors namespace FEValuesViews { - template + template inline typename Scalar::value_type Scalar::value (const unsigned int shape_function, @@ -3423,7 +3740,7 @@ namespace FEValuesViews 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 @@ -3438,11 +3755,11 @@ namespace FEValuesViews else return 0; } - - template + + template inline typename Scalar::gradient_type Scalar::gradient (const unsigned int shape_function, @@ -3452,7 +3769,7 @@ namespace FEValuesViews 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 @@ -3469,7 +3786,7 @@ namespace FEValuesViews - template + template inline typename Scalar::hessian_type Scalar::hessian (const unsigned int shape_function, @@ -3479,7 +3796,7 @@ namespace FEValuesViews 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 @@ -3495,7 +3812,7 @@ namespace FEValuesViews - template + template inline typename Vector::value_type Vector::value (const unsigned int shape_function, @@ -3505,7 +3822,7 @@ namespace FEValuesViews 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 @@ -3522,28 +3839,28 @@ namespace FEValuesViews else { value_type return_value; - for (unsigned int d=0; d + + template inline typename Vector::gradient_type Vector::gradient (const unsigned int shape_function, const unsigned int q_point) const { - typedef FEValuesBase FVB; + typedef FEValuesBase 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 @@ -3560,7 +3877,7 @@ namespace FEValuesViews else { gradient_type return_value; - for (unsigned int d=0; d + template inline typename Vector::divergence_type Vector::divergence (const unsigned int shape_function, @@ -3579,11 +3896,11 @@ namespace FEValuesViews { // this function works like in // the case above - typedef FEValuesBase FVB; + typedef FEValuesBase 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 @@ -3596,7 +3913,7 @@ namespace FEValuesViews else { divergence_type return_value = 0; - for (unsigned int d=0; d + + template inline typename Vector::hessian_type Vector::hessian (const unsigned int shape_function, const unsigned int q_point) const { // this function works like in - // the case above - typedef FEValuesBase FVB; + // the case above + typedef FEValuesBase 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 @@ -3636,7 +3953,7 @@ namespace FEValuesViews else { hessian_type return_value; - for (unsigned int d=0; d + 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) { @@ -3677,24 +3994,24 @@ namespace FEValuesViews 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) { @@ -3703,27 +4020,27 @@ namespace FEValuesViews 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 inline @@ -3731,11 +4048,11 @@ namespace FEValuesViews Vector::symmetric_gradient (const unsigned int shape_function, const unsigned int q_point) const { - typedef FEValuesBase FVB; + typedef FEValuesBase 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 @@ -3748,7 +4065,7 @@ namespace FEValuesViews else { gradient_type return_value; - for (unsigned int d=0; d +template inline const FEValuesViews::Vector & FEValuesBase:: @@ -3789,13 +4106,25 @@ operator[] (const FEValuesExtractors::Vector &vector) const 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 +inline +const FEValuesViews::SymmetricTensor<2,dim,spacedim> & +FEValuesBase:: +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 +template inline const double & FEValuesBase::shape_value (const unsigned int i, @@ -3827,7 +4156,7 @@ FEValuesBase::shape_value (const unsigned int i, -template +template inline double FEValuesBase::shape_value_component (const unsigned int i, @@ -3840,7 +4169,7 @@ FEValuesBase::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 @@ -3891,7 +4220,7 @@ FEValuesBase::shape_value_component (const unsigned int i, -template +template inline const Tensor<1,spacedim> & FEValuesBase::shape_grad (const unsigned int i, @@ -3927,7 +4256,7 @@ FEValuesBase::shape_grad (const unsigned int i, -template +template inline Tensor<1,spacedim> FEValuesBase::shape_grad_component (const unsigned int i, @@ -3940,7 +4269,7 @@ FEValuesBase::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 @@ -3991,7 +4320,7 @@ FEValuesBase::shape_grad_component (const unsigned int i, -template +template inline const Tensor<2,spacedim> & FEValuesBase::shape_hessian (const unsigned int i, @@ -4027,7 +4356,7 @@ FEValuesBase::shape_hessian (const unsigned int i, -template +template inline const Tensor<2,spacedim> & FEValuesBase::shape_2nd_derivative (const unsigned int i, @@ -4038,7 +4367,7 @@ FEValuesBase::shape_2nd_derivative (const unsigned int i, -template +template inline Tensor<2,spacedim> FEValuesBase::shape_hessian_component (const unsigned int i, @@ -4051,7 +4380,7 @@ FEValuesBase::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 @@ -4102,7 +4431,7 @@ FEValuesBase::shape_hessian_component (const unsigned int i, -template +template inline Tensor<2,spacedim> FEValuesBase::shape_2nd_derivative_component (const unsigned int i, @@ -4114,18 +4443,18 @@ FEValuesBase::shape_2nd_derivative_component (const unsigned int i -template +template inline -const FiniteElement & +const FiniteElement & FEValuesBase::get_fe () const { return *fe; } -template +template inline -const Mapping & +const Mapping & FEValuesBase::get_mapping () const { return *mapping; @@ -4133,7 +4462,7 @@ FEValuesBase::get_mapping () const -template +template inline UpdateFlags FEValuesBase::get_update_flags () const @@ -4143,7 +4472,7 @@ FEValuesBase::get_update_flags () const -template +template const typename Triangulation::cell_iterator FEValuesBase::get_cell () const { @@ -4151,7 +4480,7 @@ FEValuesBase::get_cell () const } -template +template inline const std::vector > & FEValuesBase::get_quadrature_points () const @@ -4162,7 +4491,7 @@ FEValuesBase::get_quadrature_points () const -template +template inline const std::vector & FEValuesBase::get_JxW_values () const @@ -4173,7 +4502,7 @@ FEValuesBase::get_JxW_values () const -template +template inline const std::vector >& FEValuesBase::get_jacobians () const @@ -4212,34 +4541,34 @@ FEValuesBase::quadrature_point (const unsigned int i) const { Assert (this->update_flags & update_quadrature_points, ExcAccessToUninitializedField()); Assert (iquadrature_points.size(), ExcIndexRange(i, 0, this->quadrature_points.size())); - + return this->quadrature_points[i]; } -template +template inline double FEValuesBase::JxW (const unsigned int i) const { Assert (this->update_flags & update_JxW_values, ExcAccessToUninitializedField()); Assert (iJxW_values.size(), ExcIndexRange(i, 0, this->JxW_values.size())); - + return this->JxW_values[i]; } -template +template inline const Tensor<2,spacedim> & FEValuesBase::jacobian (const unsigned int i) const { Assert (this->update_flags & update_jacobians, ExcAccessToUninitializedField()); Assert (ijacobians.size(), ExcIndexRange(i, 0, this->jacobians.size())); - + return this->jacobians[i]; } @@ -4252,7 +4581,7 @@ FEValuesBase::jacobian_grad (const unsigned int i) const { Assert (this->update_flags & update_jacobian_grads, ExcAccessToUninitializedField()); Assert (ijacobian_grads.size(), ExcIndexRange(i, 0, this->jacobian_grads.size())); - + return this->jacobian_grads[i]; } @@ -4265,7 +4594,7 @@ FEValuesBase::inverse_jacobian (const unsigned int i) const { Assert (this->update_flags & update_inverse_jacobians, ExcAccessToUninitializedField()); Assert (iinverse_jacobians.size(), ExcIndexRange(i, 0, this->inverse_jacobians.size())); - + return this->inverse_jacobians[i]; } @@ -4283,7 +4612,7 @@ FEValuesBase::get_function_grads (const InputVector &fe_ -template +template template inline void @@ -4297,7 +4626,7 @@ FEValuesBase::get_function_grads ( -template +template template inline void @@ -4310,7 +4639,7 @@ get_function_grads (const InputVector &fe_function, -template +template template inline void @@ -4325,7 +4654,7 @@ FEValuesBase::get_function_grads ( -template +template template inline void @@ -4338,7 +4667,7 @@ get_function_2nd_derivatives (const InputVector &fe_function, -template +template template inline void @@ -4352,7 +4681,7 @@ get_function_2nd_derivatives (const InputVector &fe_func -template +template inline const Point & FEValuesBase::cell_normal_vector (const unsigned int i) const @@ -4362,7 +4691,7 @@ FEValuesBase::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]; } @@ -4372,20 +4701,20 @@ FEValuesBase::cell_normal_vector (const unsigned int i) const /*------------------------ Inline functions: FEValues ----------------------------*/ -template +template inline const Quadrature & -FEValues::get_quadrature () const +FEValues::get_quadrature () const { return quadrature; } -template +template inline const FEValues & -FEValues::get_present_fe_values () const +FEValues::get_present_fe_values () const { return *this; } @@ -4394,7 +4723,7 @@ FEValues::get_present_fe_values () const /*------------------------ Inline functions: FEFaceValuesBase --------------------*/ -template +template inline const Point & FEFaceValuesBase::normal_vector (const unsigned int i) const @@ -4404,12 +4733,12 @@ FEFaceValuesBase::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 +template inline unsigned int FEFaceValuesBase::get_face_index () const @@ -4420,37 +4749,37 @@ FEFaceValuesBase::get_face_index () const /*------------------------ Inline functions: FE*FaceValues --------------------*/ -template +template inline const Quadrature & -FEFaceValuesBase::get_quadrature () const +FEFaceValuesBase::get_quadrature () const { return quadrature; } -template +template inline const FEFaceValues & -FEFaceValues::get_present_fe_values () const +FEFaceValues::get_present_fe_values () const { return *this; } -template +template inline const FESubfaceValues & -FESubfaceValues::get_present_fe_values () const +FESubfaceValues::get_present_fe_values () const { return *this; } -template +template inline const Tensor<1,spacedim> & FEFaceValuesBase::boundary_form (const unsigned int i) const @@ -4460,7 +4789,7 @@ FEFaceValuesBase::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]; } diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 333157f2a3..8a3a581771 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -55,7 +55,7 @@ namespace namespace FEValuesViews { - template + template Scalar::Scalar (const FEValuesBase &fe_values, const unsigned int component) : @@ -68,7 +68,7 @@ namespace FEValuesViews const std::vector shape_function_to_row_table = make_shape_function_to_row_table (*fe_values.fe); - + for (unsigned int i=0; idofs_per_cell; ++i) { const bool is_primitive = (fe_values.fe->is_primitive() || @@ -119,10 +119,10 @@ namespace FEValuesViews Assert (false, ExcInternalError()); return *this; } - - template + + template Vector::Vector (const FEValuesBase &fe_values, const unsigned int first_vector_component) : @@ -136,11 +136,11 @@ namespace FEValuesViews const std::vector shape_function_to_row_table = make_shape_function_to_row_table (*fe_values.fe); - + for (unsigned int d=0; ddofs_per_cell; ++i) { const bool is_primitive = (fe_values.fe->is_primitive() || @@ -203,7 +203,7 @@ namespace FEValuesViews } } } - + template Vector::Vector () @@ -224,17 +224,140 @@ namespace FEValuesViews } + template + SymmetricTensor<2, dim, spacedim>:: + SymmetricTensor(const FEValuesBase &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 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 + SymmetricTensor<2, dim, spacedim>::SymmetricTensor() + : + fe_values(*static_cast*> (0)), + first_tensor_component(numbers::invalid_unsigned_int) + {} + + + template + 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 template void Scalar:: - get_function_values (const InputVector &fe_function, + get_function_values (const InputVector &fe_function, std::vector &values) const { typedef FEValuesBase 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, @@ -258,25 +381,25 @@ namespace FEValuesViews 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 + + template template void Scalar:: - get_function_gradients (const InputVector &fe_function, + get_function_gradients (const InputVector &fe_function, std::vector &gradients) const { typedef FEValuesBase 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, @@ -291,7 +414,7 @@ namespace FEValuesViews 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_functiondofs_per_cell; ++shape_function) if (shape_function_data[shape_function].is_nonzero_shape_function_component) @@ -300,7 +423,7 @@ namespace FEValuesViews 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 + template template void Scalar:: - get_function_hessians (const InputVector &fe_function, + get_function_hessians (const InputVector &fe_function, std::vector &hessians) const { typedef FEValuesBase FVB; @@ -327,7 +450,7 @@ namespace FEValuesViews 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 dof_values (fe_values.dofs_per_cell); @@ -343,7 +466,7 @@ namespace FEValuesViews 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 void Scalar:: - get_function_laplacians (const InputVector &fe_function, + get_function_laplacians (const InputVector &fe_function, std::vector &laplacians) const { typedef FEValuesBase 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, @@ -388,7 +511,7 @@ namespace FEValuesViews const unsigned int row_index = shape_function_data[shape_function].row_index; for (unsigned int q_point=0; q_point void Vector:: - get_function_values (const InputVector &fe_function, + get_function_values (const InputVector &fe_function, std::vector &values) const { typedef FEValuesBase 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, @@ -424,7 +547,7 @@ namespace FEValuesViews shape_functiondofs_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 @@ -453,20 +576,20 @@ namespace FEValuesViews } } } - - template + + template template void Vector:: - get_function_gradients (const InputVector &fe_function, + get_function_gradients (const InputVector &fe_function, std::vector &gradients) const { typedef FEValuesBase 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, @@ -481,12 +604,12 @@ namespace FEValuesViews 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_functiondofs_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 @@ -498,18 +621,18 @@ namespace FEValuesViews 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 * 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 + template template void Vector:: - get_function_symmetric_gradients (const InputVector &fe_function, + get_function_symmetric_gradients (const InputVector &fe_function, std::vector &symmetric_gradients) const { typedef FEValuesBase 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, @@ -543,14 +666,14 @@ namespace FEValuesViews dealii::Vector 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_functiondofs_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 @@ -562,9 +685,9 @@ namespace FEValuesViews 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 + template template void Vector:: - get_function_divergences (const InputVector &fe_function, + get_function_divergences (const InputVector &fe_function, std::vector &divergences) const { typedef FEValuesBase 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, @@ -609,12 +732,12 @@ namespace FEValuesViews 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_functiondofs_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 @@ -626,7 +749,7 @@ namespace FEValuesViews 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]; @@ -634,7 +757,7 @@ namespace FEValuesViews divergences[q_point] += value * (*shape_gradient_ptr++)[comp]; } else - for (unsigned int d=0; d * shape_gradient_ptr = @@ -645,14 +768,14 @@ namespace FEValuesViews } } } - - template + + template template void Vector:: - get_function_hessians (const InputVector &fe_function, + get_function_hessians (const InputVector &fe_function, std::vector &hessians) const { typedef FEValuesBase FVB; @@ -665,7 +788,7 @@ namespace FEValuesViews 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 dof_values (fe_values.dofs_per_cell); @@ -677,7 +800,7 @@ namespace FEValuesViews shape_functiondofs_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 @@ -689,9 +812,9 @@ namespace FEValuesViews 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 * 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 void Vector:: - get_function_laplacians (const InputVector &fe_function, + get_function_laplacians (const InputVector &fe_function, std::vector &laplacians) const { typedef FEValuesBase 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, @@ -740,7 +863,7 @@ namespace FEValuesViews shape_functiondofs_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 @@ -754,7 +877,7 @@ namespace FEValuesViews { 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 * 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 + template + void + SymmetricTensor<2, dim, spacedim>:: + get_function_values(const InputVector &fe_function, + std::vector &values) const + { + typedef FEValuesBase 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 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 + template + void + SymmetricTensor<2, dim, spacedim>:: + get_function_divergences(const InputVector &fe_function, + std::vector &divergences) const + { + typedef FEValuesBase 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 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]; + } + } + } + } + } } @@ -801,7 +1071,7 @@ namespace internal dealii::FEValuesViews::Scalar(fe_values, component); } - + const unsigned int n_vectors = (fe.n_components() >= dim ? fe.n_components()-dim+1 : 0); @@ -813,7 +1083,20 @@ namespace internal new (&vectors[component]) dealii::FEValuesViews::Vector(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); + } } } } @@ -835,7 +1118,7 @@ FEValuesBase::CellIteratorBase::~CellIteratorBase () * interface. See there for a * description of the use of * these classes. - * + * * @author Wolfgang Bangerth, 2003 */ template @@ -849,7 +1132,7 @@ class FEValuesBase::CellIterator : public FEValuesBase::CellIterator : public FEValuesBase::cell_iterator () const; - + /** * Return the number of * degrees of freedom the DoF @@ -1020,7 +1303,7 @@ class FEValuesBase::CellIterator : public FEValuesBase::CellIterator : public FEValuesBase @@ -1075,7 +1358,7 @@ class FEValuesBase::TriaCellIterator : public FEValuesBase::cell_iterator &cell); - + /** * Conversion operator to an * iterator for @@ -1263,7 +1546,7 @@ class FEValuesBase::TriaCellIterator : public FEValuesBase::initialize (const unsigned int n_quadrature_p // functions unsigned int n_nonzero_shape_components = 0; for (unsigned int i=0; i= fe.dofs_per_cell, ExcInternalError()); @@ -1663,7 +1946,7 @@ FEValuesData::initialize (const unsigned int n_quadrature_p if (flags & update_hessians) this->shape_hessians.resize (n_nonzero_shape_components, std::vector > (n_quadrature_points)); - + if (flags & update_quadrature_points) this->quadrature_points.resize(n_quadrature_points); @@ -1818,10 +2101,10 @@ void FEValuesBase::get_function_values ( // 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 @@ -1864,7 +2147,7 @@ void FEValuesBase::get_function_values ( const InputVector& fe_function, std::vector >& 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(), @@ -1883,16 +2166,16 @@ void FEValuesBase::get_function_values ( for (unsigned i=0;iupdate_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 dof_values (dofs_per_cell); present_cell->get_interpolated_dof_values(fe_function, dof_values); - + // initialize with zero for (unsigned i=0;i::get_function_values ( // 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); @@ -1983,9 +2266,9 @@ void FEValuesBase::get_function_values ( // 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 @@ -1999,7 +2282,7 @@ void FEValuesBase::get_function_values ( // number of components of the // finite element const unsigned int result_components = indices.size() * n_components / dofs_per_cell; - + for (unsigned i=0;i::get_function_values ( // 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::get_function_values ( 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::get_function_values ( const double *shape_value_ptr = &this->shape_values(row, 0); const unsigned int comp = c + mc * n_components; - + for (unsigned int point=0; point::get_function_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 @@ -2107,15 +2390,15 @@ void FEValuesBase::get_function_values ( 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::get_function_values ( 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) @@ -2165,7 +2448,7 @@ void FEValuesBase::get_function_values ( 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::get_function_gradients ( 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::get_function_gradients ( // 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 @@ -2293,7 +2576,7 @@ void FEValuesBase::get_function_gradients ( 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::get_function_gradients ( 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::get_function_gradients ( 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::get_function_gradients ( 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 @@ -2404,7 +2687,7 @@ void FEValuesBase::get_function_gradients ( // 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) @@ -2429,9 +2712,9 @@ void FEValuesBase::get_function_gradients ( // 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()); @@ -2451,10 +2734,10 @@ void FEValuesBase::get_function_gradients ( 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; @@ -2479,7 +2762,7 @@ void FEValuesBase::get_function_gradients ( 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; @@ -2533,7 +2816,7 @@ get_function_hessians (const InputVector &fe_function, 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::get_function_hessians ( // 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 @@ -2577,7 +2860,7 @@ void FEValuesBase::get_function_hessians ( 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; pointis_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; @@ -2657,7 +2940,7 @@ get_function_hessians (const InputVector &fe_function, 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) @@ -2681,9 +2964,9 @@ void FEValuesBase::get_function_hessians ( 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 @@ -2697,7 +2980,7 @@ void FEValuesBase::get_function_hessians ( // 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) @@ -2722,7 +3005,7 @@ void FEValuesBase::get_function_hessians ( // 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()); @@ -2742,11 +3025,11 @@ void FEValuesBase::get_function_hessians ( 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; @@ -2771,7 +3054,7 @@ void FEValuesBase::get_function_hessians ( 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; @@ -2827,7 +3110,7 @@ void FEValuesBase::get_function_laplacians ( 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::get_function_laplacians ( // 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 @@ -2874,7 +3157,7 @@ void FEValuesBase::get_function_laplacians ( 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::get_function_laplacians ( const InputVector& fe_function, std::vector >& 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(), @@ -2908,16 +3191,16 @@ void FEValuesBase::get_function_laplacians ( for (unsigned i=0;iupdate_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 dof_values (dofs_per_cell); present_cell->get_interpolated_dof_values(fe_function, dof_values); - + // initialize with zero for (unsigned i=0;i::get_function_laplacians ( 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::get_function_laplacians ( 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::get_function_laplacians ( // 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 @@ -2996,7 +3279,7 @@ void FEValuesBase::get_function_laplacians ( // number of components of the // finite element const unsigned int result_components = indices.size() * n_components / dofs_per_cell; - + for (unsigned i=0;i::get_function_laplacians ( // 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::get_function_laplacians ( 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::get_function_laplacians ( 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::get_function_laplacians ( bool quadrature_points_fastest) const { const unsigned int n_components = fe->n_components(); - + // Size of indices must be a // multiple of dofs_per_cell such // that an integer number of @@ -3106,15 +3389,15 @@ void FEValuesBase::get_function_laplacians ( 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::get_function_laplacians ( 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::get_function_laplacians ( 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 UpdateFlags FEValuesBase::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 @@ -3243,7 +3526,7 @@ FEValuesBase::compute_update_flags (const UpdateFlags update_flags template inline void -FEValuesBase::check_cell_similarity +FEValuesBase::check_cell_similarity (const typename Triangulation::cell_iterator &cell) { // case that there has not been any cell @@ -3260,11 +3543,11 @@ FEValuesBase::check_cell_similarity else cell_similarity = (cell->is_translation_of (static_cast::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. @@ -3273,7 +3556,7 @@ FEValuesBase::check_cell_similarity template -CellSimilarity::Similarity +CellSimilarity::Similarity FEValuesBase::get_cell_similarity () const { return cell_similarity; @@ -3346,7 +3629,7 @@ FEValues::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 @@ -3381,7 +3664,7 @@ FEValues::reinit (const typename DoFHandler::cell_it // 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::template CellIterator::cell_iterator> (cell)); @@ -3415,7 +3698,7 @@ FEValues::reinit (const typename hp::DoFHandler::cel // 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::template CellIterator::cell_iterator> (cell)); @@ -3460,7 +3743,7 @@ FEValues::reinit (const typename MGDoFHandler::cell_ // this cell, are the same //TODO: This was documented out ith the repository. Why? - + // Assert (static_cast&>(*this->fe) == // static_cast&>(cell->get_fe()), // typename FEValuesBase::ExcFEDontMatch()); @@ -3472,7 +3755,7 @@ FEValues::reinit (const typename MGDoFHandler::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::template CellIterator::cell_iterator> (cell)); @@ -3498,7 +3781,7 @@ void FEValues::reinit (const typename Triangulation: // 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::TriaCellIterator (cell)); // this was the part of the work // that is dependent on the actual @@ -3554,7 +3837,7 @@ template FEFaceValuesBase::FEFaceValuesBase (const unsigned int n_q_points, const unsigned int dofs_per_cell, const UpdateFlags, - const Mapping &mapping, + const Mapping &mapping, const FiniteElement &fe, const Quadrature& quadrature) : @@ -3649,7 +3932,7 @@ void FEFaceValues::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 @@ -3677,16 +3960,16 @@ void FEFaceValues::reinit (const typename DoFHandler Assert (face_no < GeometryInfo::faces_per_cell, ExcIndexRange (face_no, 0, GeometryInfo::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::template CellIterator::cell_iterator> (cell)); - + // this was the part of the work // that is dependent on the actual // data type of the iterator. now @@ -3713,13 +3996,13 @@ void FEFaceValues::reinit (const typename hp::DoFHandler::faces_per_cell, ExcIndexRange (face_no, 0, GeometryInfo::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::template CellIterator::cell_iterator> (cell)); @@ -3769,13 +4052,13 @@ void FEFaceValues::reinit (const typename MGDoFHandler::faces_per_cell, ExcIndexRange (face_no, 0, GeometryInfo::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::template CellIterator::cell_iterator> (cell)); @@ -3800,7 +4083,7 @@ void FEFaceValues::reinit (const typename Triangulationpresent_cell.reset + this->present_cell.reset (new typename FEValuesBase::TriaCellIterator (cell)); // this was the part of the work @@ -3821,7 +4104,7 @@ void FEFaceValues::do_reinit (const unsigned int face_no) // available) const typename Triangulation::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, @@ -3892,7 +4175,7 @@ void FESubfaceValues::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 @@ -3947,7 +4230,7 @@ void FESubfaceValues::reinit (const typename DoFHandlerpresent_cell.reset + this->present_cell.reset (new typename FEValuesBase::template CellIterator::cell_iterator> (cell)); @@ -4001,7 +4284,7 @@ void FESubfaceValues::reinit (const typename hp::DoFHandler void FESubfaceValues::reinit (const typename MGDoFHandler::cell_iterator &cell, const unsigned int face_no, @@ -4025,7 +4308,7 @@ void FESubfaceValues::reinit (const typename MGDoFHandlerpresent_cell.reset + this->present_cell.reset (new typename FEValuesBase::template CellIterator::cell_iterator> (cell)); @@ -4048,13 +4331,13 @@ void FESubfaceValues::reinit (const typename Triangulation::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::TriaCellIterator (cell)); // this was the part of the work @@ -4139,7 +4422,7 @@ void FESubfaceValues::do_reinit (const unsigned int face_no, 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, @@ -4150,7 +4433,7 @@ void FESubfaceValues::do_reinit (const unsigned int face_no, 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, @@ -4184,8 +4467,8 @@ template class FESubfaceValues; template class FEValuesData; template class FEValuesBase; template class FEValues; -template class FEValuesBase:: -CellIterator::cell_iterator>; +template class FEValuesBase:: +CellIterator::cell_iterator>; //template class FEValuesBase:: // CellIterator::cell_iterator>; diff --git a/deal.II/deal.II/source/fe/fe_values.inst.in b/deal.II/deal.II/source/fe/fe_values.inst.in index 5009baf9a5..de06ad83ff 100644 --- a/deal.II/deal.II/source/fe/fe_values.inst.in +++ b/deal.II/deal.II/source/fe/fe_values.inst.in @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -37,7 +37,7 @@ for (VEC : SERIAL_VECTORS) (const dealii::VEC&, std::vector >&) const; template void FEValuesViews::Vector::get_function_symmetric_gradients - (const dealii::VEC&, std::vector >&) const; + (const dealii::VEC&, std::vector >&) const; template void FEValuesViews::Vector::get_function_divergences (const dealii::VEC&, std::vector&) const; @@ -48,6 +48,13 @@ for (VEC : SERIAL_VECTORS) void FEValuesViews::Vector::get_function_laplacians (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::SymmetricTensor<2,deal_II_dimension>::get_function_values + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::SymmetricTensor<2,deal_II_dimension>::get_function_divergences + (const dealii::VEC&, std::vector >&) const; + #if deal_II_dimension != 3 template @@ -78,7 +85,7 @@ for (VEC : SERIAL_VECTORS) template void FEValuesViews::Vector ::get_function_symmetric_gradients - (const dealii::VEC&, std::vector >&) const; + (const dealii::VEC&, std::vector >&) const; template void FEValuesViews::Vector ::get_function_divergences @@ -92,6 +99,15 @@ for (VEC : SERIAL_VECTORS) ::get_function_laplacians (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::SymmetricTensor<2,deal_II_dimension+1> + ::get_function_values + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::SymmetricTensor<2,deal_II_dimension+1> + ::get_function_divergences + (const dealii::VEC&, std::vector >&) const; + #endif } -- 2.39.5