From d8942f8304870adb3ca6172df16fc66c3861374a Mon Sep 17 00:00:00 2001 From: leicht Date: Fri, 4 Aug 2006 05:02:02 +0000 Subject: [PATCH] extend DerivativeApproximation to third order derivatives, enable returning the full derivative tensor git-svn-id: https://svn.dealii.org/trunk@13591 0785d39b-7218-0410-832d-ea1e28bc413d --- .../numerics/derivative_approximation.h | 190 ++++++++++ .../numerics/derivative_approximation.cc | 353 +++++++++++++++--- 2 files changed, 485 insertions(+), 58 deletions(-) diff --git a/deal.II/deal.II/include/numerics/derivative_approximation.h b/deal.II/deal.II/include/numerics/derivative_approximation.h index 291db5f75f..62098de3d6 100644 --- a/deal.II/deal.II/include/numerics/derivative_approximation.h +++ b/deal.II/deal.II/include/numerics/derivative_approximation.h @@ -1,3 +1,4 @@ + //--------------------------------------------------------------------------- // $Id$ // Version: $Name$ @@ -265,6 +266,53 @@ class DerivativeApproximation Vector &derivative_norm, const unsigned int component = 0); + /** + * This function calculates the + * order-th order approximate + * derivative and returns the full tensor + * for a single cell. + * + * The last parameter denotes the + * solution component, for which + * the gradient is to be + * computed. It defaults to the + * first component. For + * scalar elements, this is the only + * valid choice; for vector-valued ones, + * any component between zero and the + * number of vector components can be + * given here. + */ + + template class DH, class InputVector, int order> + static void + approximate_derivative_tensor (const Mapping &mapping, + const DH &dof, + const InputVector &solution, + const typename DH::active_cell_iterator &cell, + Tensor &derivative, + const unsigned int component = 0); + + /** + * Same as above, with + * mapping=MappingQ1@(). + */ + + template class DH, class InputVector, int order> + static void + approximate_derivative_tensor (const DH &dof, + const InputVector &solution, + const typename DH::active_cell_iterator &cell, + Tensor &derivative, + const unsigned int component = 0); + + /** + * Return the norm of the derivative. + */ + template + static double + derivative_norm(const Tensor &derivative); + /** * Exception */ @@ -449,6 +497,130 @@ class DerivativeApproximation */ static void symmetrize (Derivative &derivative_tensor); }; + + template + class ThirdDerivative + { + public: + /** + * Declare which data fields have + * to be updated for the function + * @p get_projected_derivative + * to work. + */ + static const UpdateFlags update_flags; + + /** + * Declare the data type which + * holds the derivative described + * by this class. + */ + typedef Tensor<3,dim> Derivative; + + /** + * Likewise declare the data type + * that holds the derivative + * projected to a certain + * directions. + */ + typedef Tensor<2,dim> ProjectedDerivative; + + /** + * Given an FEValues object + * initialized to a cell, and a + * solution vector, extract the + * desired derivative at the + * first quadrature point (which + * is the only one, as we only + * evaluate the finite element + * field at the center of each + * cell). + */ + template + static ProjectedDerivative + get_projected_derivative (const FEValues &fe_values, + const InputVector &solution, + const unsigned int component); + + /** + * Return the norm of the + * derivative object. Here, for + * the (symmetric) tensor of + * second derivatives, we choose + * the absolute value of the + * largest eigenvalue, which is + * the matrix norm associated to + * the $l_2$ norm of vectors. It + * is also the largest value of + * the curvature of the solution. + */ + static double derivative_norm (const Derivative &d); + + /** + * If for the present derivative + * order, symmetrization of the + * derivative tensor is + * necessary, then do so on the + * argument. + * + * For the second derivatives, + * each entry of the tensor is + * set to the mean of its value + * and the value of the transpose + * element. + * + * Note that this function + * actually modifies its + * argument. + */ + static void symmetrize (Derivative &derivative_tensor); + }; + + template + class DerivativeSelector + { + public: + /** + * typedef to select the + * DerivativeDescription corresponding + * to the orderth + * derivative. In this general template + * we set an unvalid typedef to void, + * the real typedefs have to be + * specialized. + */ + typedef void DerivDescr; + + }; + + template + class DerivativeSelector<1,dim> + { + public: + + typedef Gradient DerivDescr; + }; + + template + class DerivativeSelector<2,dim> + { + public: + + typedef SecondDerivative DerivDescr; + }; + + template + class DerivativeSelector<3,dim> + { + public: + + typedef ThirdDerivative DerivDescr; + }; + + + + + private: /** * Convenience typedef denoting @@ -491,6 +663,9 @@ class DerivativeApproximation * approximation on the cells in * the range given by the third * parameter. + * Fill the @p derivative_norm vector with + * the norm of the computed derivative + * tensors on each cell. */ template class DH, class InputVector> @@ -501,6 +676,21 @@ class DerivativeApproximation const unsigned int component, const IndexInterval &index_interval, Vector &derivative_norm); + + /** + * Compute the derivative approximation on + * one cell. This computes the full + * derivative tensor. + */ + template class DH, class InputVector> + static void + approximate_cell (const Mapping &mapping, + const DH &dof, + const InputVector &solution, + const unsigned int component, + const typename DH::active_cell_iterator &cell, + typename DerivativeDescription::Derivative &derivative); }; diff --git a/deal.II/deal.II/source/numerics/derivative_approximation.cc b/deal.II/deal.II/source/numerics/derivative_approximation.cc index f33de9df38..8bffc3207f 100644 --- a/deal.II/deal.II/source/numerics/derivative_approximation.cc +++ b/deal.II/deal.II/source/numerics/derivative_approximation.cc @@ -50,6 +50,9 @@ const UpdateFlags DerivativeApproximation::Gradient::update_flags = update_ template const UpdateFlags DerivativeApproximation::SecondDerivative::update_flags = update_gradients; +template +const UpdateFlags DerivativeApproximation::ThirdDerivative::update_flags = update_second_derivatives; + template @@ -388,6 +391,113 @@ DerivativeApproximation::SecondDerivative::symmetrize (Derivative &d) } +template +template +inline +typename DerivativeApproximation::ThirdDerivative::ProjectedDerivative +DerivativeApproximation::ThirdDerivative:: +get_projected_derivative (const FEValues &fe_values, + const InputVector &solution, + const unsigned int component) +{ + if (fe_values.get_fe().n_components() == 1) + { + std::vector values (1); + fe_values.get_function_2nd_derivatives (solution, values); + return values[0]; + } + else + { + std::vector > values + (1, std::vector(fe_values.get_fe().n_components())); + fe_values.get_function_2nd_derivatives (solution, values); + return values[0][component]; + }; +} + + +#if deal_II_dimension == 1 + +template <> +inline +double +DerivativeApproximation::ThirdDerivative<1>:: +derivative_norm (const Derivative &d) +{ + return std::fabs (d[0][0][0]); +} + +#endif + + +template +inline +double +DerivativeApproximation::ThirdDerivative:: +derivative_norm (const Derivative &d) +{ + // return the Frobenius-norm. this is a + // member function of Tensor + return d.norm(); +} + + +template +inline +void +DerivativeApproximation::ThirdDerivative::symmetrize (Derivative &d) +{ + // symmetrize non-diagonal entries + + // first do it in the case, that i,j,k are + // pairwise different (which can onlky happen + // in dim >= 3) + for (unsigned int i=0; i class DH, class InputVector> @@ -458,6 +568,49 @@ approximate_second_derivative (const DH &dof_handler, } +template class DH, class InputVector, int order> +void +DerivativeApproximation:: +approximate_derivative_tensor (const Mapping &mapping, + const DH &dof, + const InputVector &solution, + const typename DH::active_cell_iterator &cell, + Tensor &derivative, + const unsigned int component) +{ + approximate_cell::DerivDescr,dim,DH,InputVector> + (mapping, + dof, + solution, + component, + cell, + derivative); +} + + + +template class DH, class InputVector, int order> +void +DerivativeApproximation:: +approximate_derivative_tensor (const DH &dof, + const InputVector &solution, + const typename DH::active_cell_iterator &cell, + Tensor &derivative, + const unsigned int component) +{ + // just call the respective function with Q1 + // mapping + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + approximate_derivative_tensor (StaticMappingQ1::mapping, + dof, + solution, + cell, + derivative, + component); +} + + + namespace WorkAround { // gcc 2.95 is not happy if we take the address of a template function as in @@ -538,6 +691,51 @@ DerivativeApproximation::approximate (const Mapping &mapping, const IndexInterval &index_interval, Vector &derivative_norm) { + // iterators over all cells and the + // respective entries in the output + // vector: + Vector::iterator + derivative_norm_on_this_cell + = derivative_norm.begin() + index_interval.first; + + typename DH::active_cell_iterator cell, endc; + cell = endc = dof_handler.begin_active(); + // (static_cast to avoid warnings + // about unsigned always >=0) + std::advance (cell, static_cast(index_interval.first)); + std::advance (endc, static_cast(index_interval.second)); + + for (; cell!=endc; ++cell, ++derivative_norm_on_this_cell) + { + typename DerivativeDescription::Derivative derivative; + // call the function doing the actual + // work on this cell + DerivativeApproximation:: + template approximate_cell + (mapping, + dof_handler, + solution, + component, + cell, + derivative); + // evaluate the norm and fill the vector + *derivative_norm_on_this_cell + = DerivativeDescription::derivative_norm (derivative); + } +} + + +template class DH, class InputVector> +void +DerivativeApproximation:: +approximate_cell (const Mapping &mapping, + const DH &dof_handler, + const InputVector &solution, + const unsigned int component, + const typename DH::active_cell_iterator &cell, + typename DerivativeDescription::Derivative &derivative) +{ QMidpoint midpoint_rule; // create collection objects from @@ -559,19 +757,6 @@ DerivativeApproximation::approximate (const Mapping &mapping, // matrix Y=sum_i y_i y_i^T Tensor<2,dim> Y; - // iterators over all cells and the - // respective entries in the output - // vector: - Vector::iterator - derivative_norm_on_this_cell - = derivative_norm.begin() + index_interval.first; - - typename DH::active_cell_iterator cell, endc; - cell = endc = dof_handler.begin_active(); - // (static_cast to avoid warnings - // about unsigned always >=0) - std::advance (cell, static_cast(index_interval.first)); - std::advance (endc, static_cast(index_interval.second)); // vector to hold iterators to all // active neighbors of a cell @@ -581,9 +766,6 @@ DerivativeApproximation::approximate (const Mapping &mapping, active_neighbors.reserve (GeometryInfo::faces_per_cell * GeometryInfo::subfaces_per_face); - for (; cell!=endc; ++cell, ++derivative_norm_on_this_cell) - { - Y.clear (); // vector // g=sum_i y_i (f(x+y_i)-f(x))/|y_i| // or related type for higher @@ -665,7 +847,7 @@ DerivativeApproximation::approximate (const Mapping &mapping, for (unsigned int c=0; cn_children(); ++c) for (unsigned int f=0; f::faces_per_cell; ++f) if (neighbor->child(c)->neighbor(f) == cell) - active_neighbors.push_back (neighbor->child(c)); + active_neighbors.push_back (neighbor->child(c)); }; }; @@ -751,60 +933,95 @@ DerivativeApproximation::approximate (const Mapping &mapping, // compute Y^-1 g const Tensor<2,dim> Y_inverse = invert(Y); - typename DerivativeDescription::Derivative derivative; contract (derivative, Y_inverse, projected_derivative); // finally symmetrize the derivative DerivativeDescription::symmetrize (derivative); +} - *derivative_norm_on_this_cell - = DerivativeDescription::derivative_norm (derivative); - } + +template +double +DerivativeApproximation:: +derivative_norm(const Tensor &derivative) +{ + return DerivativeSelector::DerivDescr::derivative_norm(derivative); } // -------------------------------------------------------------------- // explicit instantiations -#define INSTANTIATE(InputVector,DH) \ -template \ -void \ -DerivativeApproximation:: \ -approximate_gradient \ -(const Mapping &mapping, \ - const DH &dof_handler, \ - const InputVector &solution, \ - Vector &derivative_norm, \ - const unsigned int component); \ - \ -template \ -void \ -DerivativeApproximation:: \ -approximate_gradient \ -(const DH &dof_handler, \ - const InputVector &solution, \ - Vector &derivative_norm, \ - const unsigned int component); \ - \ -template \ -void \ -DerivativeApproximation:: \ -approximate_second_derivative \ -(const Mapping &mapping, \ - const DH &dof_handler, \ - const InputVector &solution, \ - Vector &derivative_norm, \ - const unsigned int component); \ - \ -template \ -void \ -DerivativeApproximation:: \ -approximate_second_derivative \ -(const DH &dof_handler, \ - const InputVector &solution, \ - Vector &derivative_norm, \ +#define INSTANTIATE(InputVector,DH) \ +template \ +void \ +DerivativeApproximation:: \ +approximate_gradient \ +(const Mapping &mapping, \ + const DH &dof_handler, \ + const InputVector &solution, \ + Vector &derivative_norm, \ + const unsigned int component); \ + \ +template \ +void \ +DerivativeApproximation:: \ +approximate_gradient \ +(const DH &dof_handler, \ + const InputVector &solution, \ + Vector &derivative_norm, \ + const unsigned int component); \ + \ +template \ +void \ +DerivativeApproximation:: \ +approximate_second_derivative \ +(const Mapping &mapping, \ + const DH &dof_handler, \ + const InputVector &solution, \ + Vector &derivative_norm, \ + const unsigned int component); \ + \ +template \ +void \ +DerivativeApproximation:: \ +approximate_second_derivative \ +(const DH &dof_handler, \ + const InputVector &solution, \ + Vector &derivative_norm, \ + const unsigned int component); \ + \ +template \ +void \ +DerivativeApproximation:: \ +approximate_derivative_tensor \ +(const DH &dof_handler, \ + const InputVector &solution, \ + const DH::active_cell_iterator &cell,\ + Tensor<1,deal_II_dimension> &derivative, \ + const unsigned int component); \ + \ +template \ +void \ +DerivativeApproximation:: \ +approximate_derivative_tensor \ +(const DH &dof_handler, \ + const InputVector &solution, \ + const DH::active_cell_iterator &cell,\ + Tensor<2,deal_II_dimension> &derivative, \ + const unsigned int component); \ + \ +template \ +void \ +DerivativeApproximation:: \ +approximate_derivative_tensor \ +(const DH &dof_handler, \ + const InputVector &solution, \ + const DH::active_cell_iterator &cell,\ + Tensor<3,deal_II_dimension> &derivative, \ const unsigned int component) + INSTANTIATE(Vector, DoFHandler); INSTANTIATE(Vector, DoFHandler); INSTANTIATE(BlockVector, DoFHandler); @@ -823,6 +1040,22 @@ INSTANTIATE(PETScWrappers::Vector, hp::DoFHandler); INSTANTIATE(PETScWrappers::BlockVector, hp::DoFHandler); #endif +template +double +DerivativeApproximation:: +derivative_norm(const Tensor<1,deal_II_dimension> &derivative); + +template +double +DerivativeApproximation:: +derivative_norm(const Tensor<2,deal_II_dimension> &derivative); + +template +double +DerivativeApproximation:: +derivative_norm(const Tensor<3,deal_II_dimension> &derivative); + + // static variables // @@ -836,3 +1069,7 @@ DerivativeApproximation::Gradient::update_flags; template const UpdateFlags DerivativeApproximation::SecondDerivative::update_flags; + +template +const UpdateFlags +DerivativeApproximation::ThirdDerivative::update_flags; -- 2.39.5