From: Wolfgang Bangerth Date: Tue, 2 Jul 2013 21:08:25 +0000 (+0000) Subject: Patch by Denis Davydov: Add Tensor extractors for FEValues objects. X-Git-Tag: v8.0.0~197 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d7b66956cd4073f79d2daf9b44ae553a79d46eea;p=dealii.git Patch by Denis Davydov: Add Tensor extractors for FEValues objects. git-svn-id: https://svn.dealii.org/trunk@29921 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 1049ce5462..e2b2a59a36 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -147,6 +147,14 @@ this function.
    +
  1. New: In addition to the FEValuesExtractors::Scalar, +FEValuesExtractors::Vector, and FEValuesExtractors::SymmetricTensor classes, +there are now also fully featured FEValuesExtractors::Tensor extractors +for non-symmetric tensors of rank 2. +
    +(Denis Davydov, 2013/07/02) +
    2. +
  2. New: There are now functions Tensor::component_to_unrolled_index() and Tensor::unrolled_to_component_indices() in the same way as they already exist for the SymmetricTensor class. diff --git a/deal.II/include/deal.II/fe/fe_values.h b/deal.II/include/deal.II/fe/fe_values.h index 42498edbec..4da0529a08 100644 --- a/deal.II/include/deal.II/fe/fe_values.h +++ b/deal.II/include/deal.II/fe/fe_values.h @@ -166,7 +166,7 @@ namespace FEValuesViews * gradient is a * Tensor@<1,dim@>. */ - typedef Tensor<1,spacedim> gradient_type; + typedef dealii::Tensor<1,spacedim> gradient_type; /** * A typedef for the type of second @@ -176,7 +176,7 @@ namespace FEValuesViews * Hessian is a * Tensor@<2,dim@>. */ - typedef Tensor<2,spacedim> hessian_type; + typedef dealii::Tensor<2,spacedim> hessian_type; /** * A structure where for each shape @@ -453,7 +453,7 @@ namespace FEValuesViews * of dim components, the * value type is a Tensor<1,spacedim>. */ - typedef Tensor<1,spacedim> value_type; + typedef dealii::Tensor<1,spacedim> value_type; /** * A typedef for the type of gradients @@ -463,7 +463,7 @@ namespace FEValuesViews * finite element, the gradient is a * Tensor@<2,spacedim@>. */ - typedef Tensor<2,spacedim> gradient_type; + typedef dealii::Tensor<2,spacedim> gradient_type; /** * A typedef for the type of @@ -512,7 +512,7 @@ namespace FEValuesViews * finite element, the Hessian is a * Tensor@<3,dim@>. */ - typedef Tensor<3,spacedim> hessian_type; + typedef dealii::Tensor<3,spacedim> hessian_type; /** * A structure where for each shape @@ -985,7 +985,7 @@ namespace FEValuesViews * definition of the * divergence. */ - typedef Tensor<1, spacedim> divergence_type; + typedef dealii::Tensor<1, spacedim> divergence_type; /** * A structure where for each shape @@ -1198,6 +1198,270 @@ namespace FEValuesViews */ std::vector shape_function_data; }; + + + template + class Tensor; + + /** + * A class representing a view to a set of + * dim*dim components forming a + * second-order tensor from a + * vector-valued finite + * element. Views are discussed in the + * @ref vector_valued module. + * + * This class allows to query the + * value and divergence of + * (components of) shape functions + * and solutions representing + * tensors. The + * divergence of a tensor + * $T_{ij}, 0\le i,j<\text{dim}$ is + * defined as + * $d_i = \sum_j \frac{\partial T_{ji}}{\partial x_j}, + * 0\le i<\text{dim}$. + * + * You get an object of this type if you + * apply a + * FEValuesExtractors::Tensor to + * an FEValues, FEFaceValues or + * FESubfaceValues object. + * + * @ingroup feaccess vector_valued + * + * @author Denis Davydov, 2013 + */ + template + class Tensor<2,dim,spacedim> + { + public: + + /** + * Data type for what you get when you apply an extractor + * of this kind to a vector-valued finite element. + */ + typedef dealii::Tensor<2, spacedim> value_type; + + /** + * Data type for taking the divergence of a tensor: a vector. + */ + typedef dealii::Tensor<1, spacedim> divergence_type; + + /** + * 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[value_type::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[value_type::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) 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; + }; + + /** + * Default constructor. Creates an + * invalid object. + */ + Tensor(); + + + /** + * Constructor for an object that + * represents (dim*dim) + * components of a + * FEValuesBase object (or of one of + * the classes derived from + * FEValuesBase), representing the unique + * components comprising a second-order + * tensor valued variable. + * + * The second argument denotes the + * index of the first component of the + * selected symmetric second order tensor. + */ + Tensor(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. + */ + Tensor &operator=(const Tensor<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) components + * (the unique components of a second-order tensor), + * the return type is a tensor of rank 2. + * + * @param shape_function Number + * of the shape function to be + * evaluated. Note that this + * number runs from zero to + * dofs_per_cell, even in the + * case of an FEFaceValues or + * FESubfaceValues object. + * + * @param q_point Number of + * the quadrature point at which + * function is to be evaluated + */ + 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. + * + * See the general discussion + * of this class for a + * definition of the + * divergence. + * + * @note The meaning of the arguments + * is as documented for the value() + * function. + */ + 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. + * + * 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. + * + * There is no equivalent function such + * as + * FEValuesBase::get_function_divergences + * in the FEValues classes but the + * information can be obtained from + * FEValuesBase::get_function_gradients, + * of course. + * + * See the general discussion + * of this class for a + * definition of the + * divergence. + */ + 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; + + /** + * Store the data about shape + * functions. + */ + std::vector shape_function_data; + }; + } @@ -1229,6 +1493,8 @@ namespace internal std::vector > vectors; std::vector > symmetric_second_order_tensors; + std::vector > + second_order_tensors; /** * Constructor. @@ -2832,6 +3098,21 @@ public: const FEValuesViews::SymmetricTensor<2,dim,spacedim> & operator[] (const FEValuesExtractors::SymmetricTensor<2> &tensor) const; + + /** + * Create a view of the current FEValues + * object that represents a set of + * (dim*dim) scalar components + * (i.e. a 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::Tensor<2,dim,spacedim> & + operator[] (const FEValuesExtractors::Tensor<2> &tensor) const; + //@} /// @name Access to the raw data @@ -3134,6 +3415,7 @@ private: template friend class FEValuesViews::Scalar; template friend class FEValuesViews::Vector; template friend class FEValuesViews::SymmetricTensor; + template friend class FEValuesViews::Tensor; }; @@ -4098,7 +4380,7 @@ namespace FEValuesViews inline dealii::SymmetricTensor<2,1> symmetrize_single_row (const unsigned int n, - const Tensor<1,1> &t) + const dealii::Tensor<1,1> &t) { Assert (n < 1, ExcIndexRange (n, 0, 1)); (void)n; // removes -Wunused-parameter warning in optimized mode @@ -4111,7 +4393,7 @@ namespace FEValuesViews inline dealii::SymmetricTensor<2,2> symmetrize_single_row (const unsigned int n, - const Tensor<1,2> &t) + const dealii::Tensor<1,2> &t) { switch (n) { @@ -4137,7 +4419,7 @@ namespace FEValuesViews inline dealii::SymmetricTensor<2,3> symmetrize_single_row (const unsigned int n, - const Tensor<1,3> &t) + const dealii::Tensor<1,3> &t) { switch (n) { @@ -4318,7 +4600,7 @@ namespace FEValuesViews // b_jj := \dfrac{\partial phi_{ii,jj}}{\partial x_jj}. // again, all other entries of 'b' are // zero - const Tensor<1, spacedim> phi_grad = fe_values.shape_gradients[snc][q_point]; + const dealii::Tensor<1, spacedim> phi_grad = fe_values.shape_gradients[snc][q_point]; divergence_type return_value; return_value[ii] = phi_grad[jj]; @@ -4388,6 +4670,21 @@ operator[] (const FEValuesExtractors::SymmetricTensor<2> &tensor) const return fe_values_views_cache.symmetric_second_order_tensors[tensor.first_tensor_component]; } +template +inline +const FEValuesViews::Tensor<2,dim,spacedim> & +FEValuesBase:: +operator[] (const FEValuesExtractors::Tensor<2> &tensor) const +{ + Assert (tensor.first_tensor_component < + fe_values_views_cache.second_order_tensors.size(), + ExcIndexRange (tensor.first_tensor_component, + 0, fe_values_views_cache.second_order_tensors.size())); + + return fe_values_views_cache.second_order_tensors[tensor.first_tensor_component]; +} + + template diff --git a/deal.II/include/deal.II/fe/fe_values_extractors.h b/deal.II/include/deal.II/fe/fe_values_extractors.h index 5d85836b08..6aed40676c 100644 --- a/deal.II/include/deal.II/fe/fe_values_extractors.h +++ b/deal.II/include/deal.II/fe/fe_values_extractors.h @@ -200,6 +200,61 @@ namespace FEValuesExtractors */ SymmetricTensor (const unsigned int first_tensor_component); }; + + + /** + * Extractor for a (possible non-)symmetric tensor of a + * rank specified by the template + * argument. For a second order + * tensor, this represents a collection of + * (dim*dim) + * components of a vector-valued + * element. The value of dim + * is defined by the FEValues object the + * extractor is applied to. The result of + * applying an object of this type to an + * FEValues, FEFaceValues or + * FESubfaceValues object is of type + * FEValuesViews::Tensor. + * + * 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 Denis Davydov, 2013 + */ + template + struct Tensor + { + /** + * The first component of the tensor + * view. + */ + unsigned int first_tensor_component; + + /** + * Default constructor. Initialize the + * object with an invalid component. This leads to + * an object that can not be used, but it allows + * objects of this kind to be put into arrays that + * require a default constructor upon resizing the + * array, and then later assigning a suitable + * object to each element of the array. + */ + Tensor (); + + /** + * Constructor. Take the first + * component of the selected tensor + * inside the FEValues object as + * argument. + */ + Tensor (const unsigned int first_tensor_component); + }; } @@ -251,6 +306,22 @@ namespace FEValuesExtractors : first_tensor_component (first_tensor_component) {} + + + template + inline + Tensor::Tensor () + : + first_tensor_component(numbers::invalid_unsigned_int) + {} + + + template + inline + Tensor::Tensor (const unsigned int first_tensor_component) + : + first_tensor_component (first_tensor_component) + {} } diff --git a/deal.II/source/fe/fe_values.cc b/deal.II/source/fe/fe_values.cc index b482eba843..7957a09afc 100644 --- a/deal.II/source/fe/fe_values.cc +++ b/deal.II/source/fe/fe_values.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998-2012 by the deal.II authors +// Copyright (C) 1998-2013 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -342,6 +342,104 @@ namespace FEValuesViews } + template + Tensor<2, dim, spacedim>:: + Tensor(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 + dim*dim - 1 + < + fe_values.fe->n_components(), + ExcIndexRange(first_tensor_component + + dim*dim - 1, + 0, + fe_values.fe->n_components())); +//TODO: we'd like to use the fields with the same name as these +// variables from FEValuesData, but they aren't initialized yet +// at the time we get here, so re-create it all + const std::vector shape_function_to_row_table + = make_shape_function_to_row_table (*fe_values.fe); + + for (unsigned int d = 0; d < dim*dim; ++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) + shape_function_data[i].row_index[d] + = shape_function_to_row_table[i*fe_values.fe->n_components()+component]; + 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 < dim*dim; ++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 < dim*dim; ++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; + } + } + } + } + + + + template + Tensor<2, dim, spacedim>::Tensor() + : + fe_values(*static_cast*> (0)), + first_tensor_component(numbers::invalid_unsigned_int) + {} + + + + template + Tensor<2, dim, spacedim> & + Tensor<2, dim, spacedim>::operator=(const Tensor<2, dim, spacedim> &) + { + // we shouldn't be copying these objects + Assert(false, ExcInternalError()); + return *this; + } + namespace internal { @@ -391,9 +489,9 @@ namespace FEValuesViews template void do_function_derivatives (const ::dealii::Vector &dof_values, - const std::vector > > &shape_derivatives, + const std::vector > > &shape_derivatives, const std::vector::ShapeFunctionData> &shape_function_data, - std::vector > &derivatives) + std::vector > &derivatives) { const unsigned int dofs_per_cell = dof_values.size(); const unsigned int n_quadrature_points = dofs_per_cell > 0 ? @@ -401,7 +499,7 @@ namespace FEValuesViews AssertDimension (derivatives.size(), n_quadrature_points); std::fill (derivatives.begin(), derivatives.end(), - Tensor()); + dealii::Tensor()); for (unsigned int shape_function=0; shape_function *shape_derivative_ptr = + const dealii::Tensor *shape_derivative_ptr = &shape_derivatives[shape_function_data[shape_function].row_index][0]; for (unsigned int q_point=0; q_point void do_function_laplacians (const ::dealii::Vector &dof_values, - const std::vector > > &shape_hessians, + const std::vector > > &shape_hessians, const std::vector::ShapeFunctionData> &shape_function_data, std::vector &laplacians) { @@ -442,7 +540,7 @@ namespace FEValuesViews if (value == 0.) continue; - const Tensor<2,spacedim> *shape_hessian_ptr = + const dealii::Tensor<2,spacedim> *shape_hessian_ptr = &shape_hessians[shape_function_data[shape_function].row_index][0]; for (unsigned int q_point=0; q_point &dof_values, const Table<2,double> &shape_values, const std::vector::ShapeFunctionData> &shape_function_data, - std::vector > &values) + std::vector > &values) { const unsigned int dofs_per_cell = dof_values.size(); const unsigned int n_quadrature_points = dofs_per_cell > 0 ? shape_values.n_cols() : values.size(); AssertDimension (values.size(), n_quadrature_points); - std::fill (values.begin(), values.end(), Tensor<1,spacedim>()); + std::fill (values.begin(), values.end(), dealii::Tensor<1,spacedim>()); for (unsigned int shape_function=0; shape_function void do_function_derivatives (const ::dealii::Vector &dof_values, - const std::vector > > &shape_derivatives, + const std::vector > > &shape_derivatives, const std::vector::ShapeFunctionData> &shape_function_data, - std::vector > &derivatives) + std::vector > &derivatives) { const unsigned int dofs_per_cell = dof_values.size(); const unsigned int n_quadrature_points = dofs_per_cell > 0 ? @@ -514,7 +612,7 @@ namespace FEValuesViews AssertDimension (derivatives.size(), n_quadrature_points); std::fill (derivatives.begin(), derivatives.end(), - Tensor()); + dealii::Tensor()); for (unsigned int shape_function=0; shape_function *shape_derivative_ptr = + const dealii::Tensor *shape_derivative_ptr = &shape_derivatives[snc][0]; for (unsigned int q_point=0; q_point *shape_derivative_ptr = + const dealii::Tensor *shape_derivative_ptr = &shape_derivatives[shape_function_data[shape_function]. row_index[d]][0]; for (unsigned int q_point=0; q_point void do_function_symmetric_gradients (const ::dealii::Vector &dof_values, - const std::vector > > &shape_gradients, + const std::vector > > &shape_gradients, const std::vector::ShapeFunctionData> &shape_function_data, std::vector > &symmetric_gradients) { @@ -585,7 +683,7 @@ namespace FEValuesViews { const unsigned int comp = shape_function_data[shape_function].single_nonzero_component_index; - const Tensor<1,spacedim> *shape_gradient_ptr = + const dealii::Tensor<1,spacedim> *shape_gradient_ptr = &shape_gradients[snc][0]; for (unsigned int q_point=0; q_point grad; + dealii::Tensor<2,spacedim> grad; for (unsigned int d=0; d void do_function_divergences (const ::dealii::Vector &dof_values, - const std::vector > > &shape_gradients, + const std::vector > > &shape_gradients, const std::vector::ShapeFunctionData> &shape_function_data, std::vector &divergences) { @@ -637,7 +735,7 @@ namespace FEValuesViews { const unsigned int comp = shape_function_data[shape_function].single_nonzero_component_index; - const Tensor<1,spacedim> *shape_gradient_ptr = &shape_gradients[snc][0]; + const dealii::Tensor<1,spacedim> *shape_gradient_ptr = &shape_gradients[snc][0]; for (unsigned int q_point=0; q_point *shape_gradient_ptr = + const dealii::Tensor<1,spacedim> *shape_gradient_ptr = &shape_gradients[shape_function_data[shape_function]. row_index[d]][0]; for (unsigned int q_point=0; q_point void do_function_curls (const ::dealii::Vector &dof_values, - const std::vector > > &shape_gradients, + const std::vector > > &shape_gradients, const std::vector::ShapeFunctionData> &shape_function_data, std::vector::type> &curls) { @@ -696,7 +794,7 @@ namespace FEValuesViews if (snc != -1) { - const Tensor<1, spacedim> *shape_gradient_ptr = + const dealii::Tensor<1, spacedim> *shape_gradient_ptr = &shape_gradients[snc][0]; switch (shape_function_data[shape_function].single_nonzero_component_index) @@ -720,7 +818,7 @@ namespace FEValuesViews { if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) { - const Tensor<1,spacedim> *shape_gradient_ptr = + const dealii::Tensor<1,spacedim> *shape_gradient_ptr = &shape_gradients[shape_function_data[shape_function].row_index[0]][0]; for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) @@ -729,7 +827,7 @@ namespace FEValuesViews if (shape_function_data[shape_function].is_nonzero_shape_function_component[1]) { - const Tensor<1,spacedim> *shape_gradient_ptr = + const dealii::Tensor<1,spacedim> *shape_gradient_ptr = &shape_gradients[shape_function_data[shape_function].row_index[1]][0]; for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) @@ -758,7 +856,7 @@ namespace FEValuesViews if (snc != -1) { - const Tensor<1, spacedim> *shape_gradient_ptr = &shape_gradients[snc][0]; + const dealii::Tensor<1, spacedim> *shape_gradient_ptr = &shape_gradients[snc][0]; switch (shape_function_data[shape_function].single_nonzero_component_index) { @@ -800,7 +898,7 @@ namespace FEValuesViews { if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) { - const Tensor<1,spacedim> *shape_gradient_ptr = + const dealii::Tensor<1,spacedim> *shape_gradient_ptr = &shape_gradients[shape_function_data[shape_function].row_index[0]][0]; for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) @@ -812,7 +910,7 @@ namespace FEValuesViews if (shape_function_data[shape_function].is_nonzero_shape_function_component[1]) { - const Tensor<1,spacedim> *shape_gradient_ptr = + const dealii::Tensor<1,spacedim> *shape_gradient_ptr = &shape_gradients[shape_function_data[shape_function].row_index[1]][0]; for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) @@ -824,7 +922,7 @@ namespace FEValuesViews if (shape_function_data[shape_function].is_nonzero_shape_function_component[2]) { - const Tensor<1,spacedim> *shape_gradient_ptr = + const dealii::Tensor<1,spacedim> *shape_gradient_ptr = &shape_gradients[shape_function_data[shape_function].row_index[2]][0]; for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) @@ -844,16 +942,17 @@ namespace FEValuesViews template void do_function_laplacians (const ::dealii::Vector &dof_values, - const std::vector > > &shape_hessians, + const std::vector > > &shape_hessians, const std::vector::ShapeFunctionData> &shape_function_data, - std::vector > &laplacians) + std::vector > &laplacians) { const unsigned int dofs_per_cell = dof_values.size(); const unsigned int n_quadrature_points = dofs_per_cell > 0 ? shape_hessians[0].size() : laplacians.size(); AssertDimension (laplacians.size(), n_quadrature_points); - std::fill (laplacians.begin(), laplacians.end(), Tensor<1,spacedim>()); + std::fill (laplacians.begin(), laplacians.end(), + dealii::Tensor<1,spacedim>()); for (unsigned int shape_function=0; shape_function *shape_hessian_ptr = + const dealii::Tensor<2,spacedim> *shape_hessian_ptr = &shape_hessians[snc][0]; for (unsigned int q_point=0; q_point *shape_hessian_ptr = + const dealii::Tensor<2,spacedim> *shape_hessian_ptr = &shape_hessians[shape_function_data[shape_function]. row_index[d]][0]; for (unsigned int q_point=0; q_point void do_function_divergences (const ::dealii::Vector &dof_values, - const std::vector > > &shape_gradients, + const std::vector > > &shape_gradients, const std::vector::ShapeFunctionData> &shape_function_data, - std::vector > &divergences) + std::vector > &divergences) { const unsigned int dofs_per_cell = dof_values.size(); const unsigned int n_quadrature_points = dofs_per_cell > 0 ? shape_gradients[0].size() : divergences.size(); AssertDimension (divergences.size(), n_quadrature_points); - std::fill (divergences.begin(), divergences.end(), Tensor<1,spacedim>()); + std::fill (divergences.begin(), divergences.end(), + dealii::Tensor<1,spacedim>()); for (unsigned int shape_function=0; shape_function *shape_gradient_ptr = + const dealii::Tensor < 1, spacedim> *shape_gradient_ptr = &shape_gradients[snc][0]; const unsigned int ii = dealii::SymmetricTensor<2,spacedim>:: @@ -1017,7 +1117,7 @@ namespace FEValuesViews const unsigned int comp = shape_function_data[shape_function].single_nonzero_component_index; - const Tensor < 1, spacedim> *shape_gradient_ptr = + const dealii::Tensor < 1, spacedim> *shape_gradient_ptr = &shape_gradients[shape_function_data[shape_function]. row_index[d]][0]; for (unsigned int q_point = 0; q_point < n_quadrature_points; @@ -1034,6 +1134,129 @@ namespace FEValuesViews } } + // ---------------------- non-symmetric tensor part ------------------------ + + template + void + do_function_values (const ::dealii::Vector &dof_values, + const Table<2,double> &shape_values, + const std::vector::ShapeFunctionData> &shape_function_data, + std::vector > &values) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = dofs_per_cell > 0 ? + shape_values.n_cols() : values.size(); + AssertDimension (values.size(), n_quadrature_points); + + std::fill (values.begin(), values.end(), + dealii::Tensor<2,spacedim>()); + + for (unsigned int shape_function=0; + shape_function indices = dealii::Tensor<2,spacedim>::unrolled_to_component_indices(comp); + const unsigned int i = indices[0], + j = indices[1]; + + const double *shape_value_ptr = &shape_values(snc,0); + for (unsigned int q_point=0; q_point indices = dealii::Tensor<2,spacedim>::unrolled_to_component_indices(d); + const unsigned int i = indices[0], + j = indices[1]; + + const double *shape_value_ptr = + &shape_values(shape_function_data[shape_function].row_index[d],0); + for (unsigned int q_point=0; q_point + void + do_function_divergences (const ::dealii::Vector &dof_values, + const std::vector > > &shape_gradients, + const std::vector::ShapeFunctionData> &shape_function_data, + std::vector > &divergences) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = dofs_per_cell > 0 ? + shape_gradients[0].size() : divergences.size(); + AssertDimension (divergences.size(), n_quadrature_points); + + std::fill (divergences.begin(), divergences.end(), + dealii::Tensor<1,spacedim>()); + + for (unsigned int shape_function=0; + shape_function *shape_gradient_ptr = + &shape_gradients[snc][0]; + + const TableIndices<2> indices = dealii::Tensor<2,spacedim>::unrolled_to_component_indices(comp); + const unsigned int ii = indices[0]; + const unsigned int jj = indices[1]; + + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point, ++shape_gradient_ptr) + { + divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj]; + + if (ii != jj) + divergences[q_point][jj] += value * (*shape_gradient_ptr)[ii]; + } + } + else + { + for (unsigned int d = 0; + d < dim*dim; ++d) + if (shape_function_data[shape_function].is_nonzero_shape_function_component[d]) + { + Assert (false, ExcNotImplemented()); + } + } + } + } + } // end of namespace internal @@ -1352,6 +1575,53 @@ namespace FEValuesViews internal::do_function_divergences (dof_values, fe_values.shape_gradients, shape_function_data, divergences); } + + template + template + void + Tensor<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(fe_values.present_cell.get() != 0, + ExcMessage("FEValues object is not reinit'ed to any cell")); + AssertDimension(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); + internal::do_function_values + (dof_values, fe_values.shape_values, shape_function_data, values); + } + + + + template + template + void + Tensor<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(fe_values.present_cell.get() != 0, + ExcMessage("FEValues object is not reinit'ed to any cell")); + AssertDimension(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); + internal::do_function_divergences + (dof_values, fe_values.shape_gradients, shape_function_data, divergences); + } } @@ -1431,6 +1701,28 @@ namespace internal dealii::FEValuesViews::SymmetricTensor<2, dim, spacedim > (fe_values, component); } + + + // compute number of symmetric + // tensors in the same way as above + const unsigned int n_second_order_tensors + = (fe.n_components() >= dim*dim ? + fe.n_components() - dim*dim + 1 : + 0); + second_order_tensors.resize(n_second_order_tensors); + for (unsigned int component = 0; component < n_second_order_tensors; ++component) + { +#ifndef DEAL_II_EXPLICIT_DESTRUCTOR_BUG + second_order_tensors[component] + .dealii::FEValuesViews::Tensor<2, dim, spacedim>::~Tensor(); +#else + typedef dealii::FEValuesViews::Tensor<2, dim, spacedim> TensorView; + second_order_tensors[component].TensorView::~TensorView(); +#endif + new (&second_order_tensors[component]) + dealii::FEValuesViews::Tensor<2, dim, spacedim > (fe_values, + component); + } } } } diff --git a/deal.II/source/fe/fe_values.inst.in b/deal.II/source/fe/fe_values.inst.in index fb01f675dc..8298609ef1 100644 --- a/deal.II/source/fe/fe_values.inst.in +++ b/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, 2009, 2010, 2012 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2012, 2013 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -36,6 +36,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS template class Scalar; template class Vector; template class SymmetricTensor<2, deal_II_dimension, deal_II_space_dimension>; + template class Tensor<2, deal_II_dimension, deal_II_space_dimension>; \} #endif #endif @@ -70,11 +71,11 @@ for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimensi template void FEValuesViews::Scalar ::get_function_gradients - (const dealii::VEC&, std::vector >&) const; + (const dealii::VEC&, std::vector >&) const; template void FEValuesViews::Scalar ::get_function_hessians - (const dealii::VEC&, std::vector >&) const; + (const dealii::VEC&, std::vector >&) const; template void FEValuesViews::Scalar ::get_function_laplacians @@ -83,11 +84,11 @@ for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimensi template void FEValuesViews::Vector ::get_function_values - (const dealii::VEC&, std::vector >&) const; + (const dealii::VEC&, std::vector >&) const; template void FEValuesViews::Vector ::get_function_gradients - (const dealii::VEC&, std::vector >&) const; + (const dealii::VEC&, std::vector >&) const; template void FEValuesViews::Vector ::get_function_symmetric_gradients @@ -103,11 +104,11 @@ for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimensi template void FEValuesViews::Vector ::get_function_hessians - (const dealii::VEC&, std::vector >&) const; + (const dealii::VEC&, std::vector >&) const; template void FEValuesViews::Vector ::get_function_laplacians - (const dealii::VEC&, std::vector >&) const; + (const dealii::VEC&, std::vector >&) const; template void FEValuesViews::SymmetricTensor<2, deal_II_dimension, deal_II_space_dimension> @@ -116,7 +117,16 @@ for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimensi template void FEValuesViews::SymmetricTensor<2, deal_II_dimension, deal_II_space_dimension> ::get_function_divergences - (const dealii::VEC&, std::vector >&) const; + (const dealii::VEC&, std::vector >&) const; + + template + void FEValuesViews::Tensor<2, deal_II_dimension, deal_II_space_dimension> + ::get_function_values + (const dealii::VEC&, std::vector >&) const; + template + void FEValuesViews::Tensor<2, deal_II_dimension, deal_II_space_dimension> + ::get_function_divergences + (const dealii::VEC&, std::vector >&) const; #endif #endif } @@ -168,35 +178,35 @@ for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimensi template void FEValuesBase::get_function_gradients - (const VEC&, std::vector > &) const; + (const VEC&, std::vector > &) const; template void FEValuesBase::get_function_gradients (const VEC&, const VectorSlice >&, - std::vector > &) const; + std::vector > &) const; template void FEValuesBase::get_function_gradients - (const VEC&, std::vector > > &) const; + (const VEC&, std::vector > > &) const; template void FEValuesBase::get_function_gradients (const VEC&, const VectorSlice >&, - VectorSlice > > >, bool) const; + VectorSlice > > >, bool) const; template void FEValuesBase::get_function_hessians - (const VEC&, std::vector > &) const; + (const VEC&, std::vector > &) const; template void FEValuesBase::get_function_hessians (const VEC&, const VectorSlice >&, - std::vector > &) const; + std::vector > &) const; template void FEValuesBase::get_function_hessians - (const VEC&, std::vector > > &, bool) const; + (const VEC&, std::vector > > &, bool) const; template void FEValuesBase::get_function_hessians (const VEC&, const VectorSlice >&, - VectorSlice > > >, bool) const; + VectorSlice > > >, bool) const; template void FEValuesBase::get_function_laplacians @@ -251,20 +261,20 @@ for (deal_II_dimension : DIMENSIONS) (const dealii::IndexSet&, std::vector&) const; template void FEValuesViews::Scalar::get_function_gradients - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; template void FEValuesViews::Scalar::get_function_hessians - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; template void FEValuesViews::Scalar::get_function_laplacians (const dealii::IndexSet&, std::vector&) const; template void FEValuesViews::Vector::get_function_values - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; template void FEValuesViews::Vector::get_function_gradients - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; template void FEValuesViews::Vector::get_function_symmetric_gradients (const dealii::IndexSet&, std::vector >&) const; @@ -276,17 +286,24 @@ for (deal_II_dimension : DIMENSIONS) (const dealii::IndexSet&, std::vector&) const; template void FEValuesViews::Vector::get_function_hessians - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; template void FEValuesViews::Vector::get_function_laplacians - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; template void FEValuesViews::SymmetricTensor<2,deal_II_dimension,deal_II_dimension>::get_function_values (const dealii::IndexSet&, std::vector >&) const; template void FEValuesViews::SymmetricTensor<2,deal_II_dimension,deal_II_dimension>::get_function_divergences - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; + + template + void FEValuesViews::Tensor<2,deal_II_dimension,deal_II_dimension>::get_function_values + (const dealii::IndexSet&, std::vector >&) const; + template + void FEValuesViews::Tensor<2,deal_II_dimension,deal_II_dimension>::get_function_divergences + (const dealii::IndexSet&, std::vector >&) const; #if deal_II_dimension != 3 @@ -297,11 +314,11 @@ for (deal_II_dimension : DIMENSIONS) template void FEValuesViews::Scalar ::get_function_gradients - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; template void FEValuesViews::Scalar ::get_function_hessians - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; template void FEValuesViews::Scalar ::get_function_laplacians @@ -310,11 +327,11 @@ for (deal_II_dimension : DIMENSIONS) template void FEValuesViews::Vector ::get_function_values - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; template void FEValuesViews::Vector ::get_function_gradients - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; template void FEValuesViews::Vector ::get_function_symmetric_gradients @@ -326,11 +343,11 @@ for (deal_II_dimension : DIMENSIONS) template void FEValuesViews::Vector ::get_function_hessians - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; template void FEValuesViews::Vector ::get_function_laplacians - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; template void FEValuesViews::SymmetricTensor<2, deal_II_dimension, deal_II_dimension+1> @@ -339,7 +356,16 @@ for (deal_II_dimension : DIMENSIONS) template void FEValuesViews::SymmetricTensor<2, deal_II_dimension, deal_II_dimension+1> ::get_function_divergences - (const dealii::IndexSet&, std::vector >&) const; + (const dealii::IndexSet&, std::vector >&) const; + + template + void FEValuesViews::Tensor<2, deal_II_dimension, deal_II_dimension+1> + ::get_function_values + (const dealii::IndexSet&, std::vector >&) const; + template + void FEValuesViews::Tensor<2, deal_II_dimension, deal_II_dimension+1> + ::get_function_divergences + (const dealii::IndexSet&, std::vector >&) const; #endif #endif @@ -392,35 +418,35 @@ for (deal_II_dimension : DIMENSIONS) template void FEValuesBase::get_function_gradients - (const IndexSet&, std::vector > &) const; + (const IndexSet&, std::vector > &) const; template void FEValuesBase::get_function_gradients (const IndexSet&, const VectorSlice >&, - std::vector > &) const; + std::vector > &) const; template void FEValuesBase::get_function_gradients - (const IndexSet&, std::vector > > &) const; + (const IndexSet&, std::vector > > &) const; template void FEValuesBase::get_function_gradients (const IndexSet&, const VectorSlice >&, - VectorSlice > > >, bool) const; + VectorSlice > > >, bool) const; template void FEValuesBase::get_function_hessians - (const IndexSet&, std::vector > &) const; + (const IndexSet&, std::vector > &) const; template void FEValuesBase::get_function_hessians (const IndexSet&, const VectorSlice >&, - std::vector > &) const; + std::vector > &) const; template void FEValuesBase::get_function_hessians - (const IndexSet&, std::vector > > &, bool) const; + (const IndexSet&, std::vector > > &, bool) const; template void FEValuesBase::get_function_hessians (const IndexSet&, const VectorSlice >&, - VectorSlice > > >, bool) const; + VectorSlice > > >, bool) const; template void FEValuesBase::get_function_laplacians @@ -503,35 +529,35 @@ for (deal_II_dimension : DIMENSIONS) template void FEValuesBase::get_function_gradients - (const IndexSet&, std::vector > &) const; + (const IndexSet&, std::vector > &) const; template void FEValuesBase::get_function_gradients (const IndexSet&, const VectorSlice >&, - std::vector > &) const; + std::vector > &) const; template void FEValuesBase::get_function_gradients - (const IndexSet&, std::vector > > &) const; + (const IndexSet&, std::vector > > &) const; template void FEValuesBase::get_function_gradients (const IndexSet&, const VectorSlice >&, - VectorSlice > > >, bool) const; + VectorSlice > > >, bool) const; template void FEValuesBase::get_function_hessians - (const IndexSet&, std::vector > &) const; + (const IndexSet&, std::vector > &) const; template void FEValuesBase::get_function_hessians (const IndexSet&, const VectorSlice >&, - std::vector > &) const; + std::vector > &) const; template void FEValuesBase::get_function_hessians - (const IndexSet&, std::vector > > &, bool) const; + (const IndexSet&, std::vector > > &, bool) const; template void FEValuesBase::get_function_hessians (const IndexSet&, const VectorSlice >&, - VectorSlice > > >, bool) const; + VectorSlice > > >, bool) const; template void FEValuesBase::get_function_laplacians