From 073a3aea2f6c33d18a63253fbd2da2506c44d969 Mon Sep 17 00:00:00 2001 From: bangerth Date: Sat, 8 Feb 2014 13:56:26 +0000 Subject: [PATCH] Make DerivativeApproximation a namespace, rather than a class with all static members. git-svn-id: https://svn.dealii.org/trunk@32438 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 7 +- .../numerics/derivative_approximation.h | 318 +-- .../numerics/derivative_approximation.cc | 1709 ++++++++++------- .../numerics/derivative_approximation.inst.in | 36 +- 4 files changed, 1005 insertions(+), 1065 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 0126d3df68..0b5acf04ed 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -124,6 +124,11 @@ inconvenience this causes.

Specific improvements

    +
  1. Changed: DerivativeApproximation used to be a class that only had + static members. It is now a namespace. +
    + (Wolfgang Bangerth, 2014/02/08) +
  2. New: ThreadLocalStorage::clear() clears out all objects allocated on the current and all other threads.
    @@ -172,7 +177,7 @@ inconvenience this causes. (Matthias Maier, 2014/02/01)
  3. New/fixed: The ParameterHandler::print_parameters_section - method not worked for XML output. There is now a flag + method did not work for XML output. There is now a flag include_top_level_elements which prints all higher subsection elements, default is false. For XML output setting this flag to true is required diff --git a/deal.II/include/deal.II/numerics/derivative_approximation.h b/deal.II/include/deal.II/numerics/derivative_approximation.h index cc43eeb479..660365f7b1 100644 --- a/deal.II/include/deal.II/numerics/derivative_approximation.h +++ b/deal.II/include/deal.II/numerics/derivative_approximation.h @@ -39,7 +39,7 @@ namespace hp /** - * This class computes a cell-wise approximation of the norm of a + * This namespace provides functions that compute a cell-wise approximation of the norm of a * derivative of a finite element field by taking difference quotients * between neighboring cells. This is a rather simple but efficient * form to get an error indicator, since it can be computed with @@ -85,7 +85,7 @@ namespace hp * @code * -------------------------------------------------------- * An error occurred in line <749> of file in function - * static void DerivativeApproximation::approximate(const Mapping&, const DH&, const InputVector&, unsigned int, const + * void DerivativeApproximation::approximate(const Mapping&, const DH&, const InputVector&, unsigned int, const * std::pair&, Vector&) [with DerivativeDescription = DerivativeApproximation::Gradient<3>, int * dim = 3, DH = DoFHandler, InputVector = Vector] * The violated condition was: @@ -158,7 +158,7 @@ namespace hp * The formulae for the computation of approximations to the gradient * and to the tensor of second derivatives shown above are very much * alike. The basic difference is that in one case the finite - * difference quotiont is a scalar, while in the other case it is a + * difference quotient is a scalar, while in the other case it is a * vector. For higher derivatives, this would be a tensor of even * higher rank. We then have to form the outer product of this * difference quotient with the distance vector $y_{KK'}$, symmetrize @@ -185,9 +185,8 @@ namespace hp * @ingroup numerics * @author Wolfgang Bangerth, 2000 */ -class DerivativeApproximation +namespace DerivativeApproximation { -public: /** * This function is used to obtain an approximation of the gradient. Pass it * the DoF handler object that describes the finite element field, a nodal @@ -204,7 +203,7 @@ public: * locally relevant unknowns. */ template class DH, class InputVector, int spacedim> - static void + void approximate_gradient (const Mapping &mapping, const DH &dof, const InputVector &solution, @@ -216,7 +215,7 @@ public: * mapping=MappingQ1@(). */ template class DH, class InputVector, int spacedim> - static void + void approximate_gradient (const DH &dof, const InputVector &solution, Vector &derivative_norm, @@ -240,7 +239,7 @@ public: * locally relevant unknowns. */ template class DH, class InputVector, int spacedim> - static void + void approximate_second_derivative (const Mapping &mapping, const DH &dof, const InputVector &solution, @@ -252,7 +251,7 @@ public: * mapping=MappingQ1@(). */ template class DH, class InputVector, int spacedim> - static void + void approximate_second_derivative (const DH &dof, const InputVector &solution, Vector &derivative_norm, @@ -272,7 +271,7 @@ public: * locally relevant unknowns. */ template class DH, class InputVector, int order, int spacedim> - static void + void approximate_derivative_tensor (const Mapping &mapping, const DH &dof, const InputVector &solution, @@ -284,7 +283,7 @@ public: * Same as above, with mapping=MappingQ1@(). */ template class DH, class InputVector, int order, int spacedim> - static void + void approximate_derivative_tensor (const DH &dof, const InputVector &solution, const typename DH::active_cell_iterator &cell, @@ -295,8 +294,8 @@ public: * Return the norm of the derivative. */ template - static double - derivative_norm(const Tensor &derivative); + double + derivative_norm (const Tensor &derivative); /** * Exception @@ -309,299 +308,8 @@ public: * Exception */ DeclException0 (ExcInsufficientDirections); +} -private: - - /** - * The following class is used to describe the data needed to compute the - * finite difference approximation to the gradient on a cell. See the - * general documentation of this class for more information on - * implementational details. - * - * @author Wolfgang Bangerth, 2000 - */ - template - class Gradient - { - 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<1,dim> Derivative; - - /** - * Likewise declare the data type that holds the derivative projected to a - * certain directions. - */ - typedef double 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 gradient, we - * choose the Euclidian norm of the gradient vector. - */ - 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 first derivatives, no such thing is necessary, so this function - * is a no-op. - */ - static void symmetrize (Derivative &derivative_tensor); - }; - - - - /** - * The following class is used to describe the data needed to compute the - * finite difference approximation to the second derivatives on a cell. See - * the general documentation of this class for more information on - * implementational details. - * - * @author Wolfgang Bangerth, 2000 - */ - template - class SecondDerivative - { - 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<2,dim> Derivative; - - /** - * Likewise declare the data type that holds the derivative projected to a - * certain directions. - */ - typedef Tensor<1,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 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 the range of indices on which a certain - * thread shall operate. - */ - typedef std::pair IndexInterval; - - /** - * Kind of the main function of this class. It is called by the public entry - * points to this class with the correct template first argument and then - * simply calls the @p approximate function, after setting up several - * threads and doing some administration that is independent of the actual - * derivative to be computed. - * - * The @p component argument denotes which component of the solution vector - * we are to work on. - */ - template class DH, class InputVector, int spacedim> - static void - approximate_derivative (const Mapping &mapping, - const DH &dof, - const InputVector &solution, - const unsigned int component, - Vector &derivative_norm); - - /** - * Compute the derivative approximation on a given cell. Fill the @p - * derivative_norm vector with the norm of the computed derivative tensors - * on the cell. - */ - template class DH, class InputVector, int spacedim> - static void - approximate (SynchronousIterators::active_cell_iterator, - Vector::iterator> > const &cell, - const Mapping &mapping, - const DH &dof, - const InputVector &solution, - const unsigned int component); - - /** - * Compute the derivative approximation on one cell. This computes the full - * derivative tensor. - */ - template class DH, class InputVector, int spacedim> - 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); -}; - - -/* -------------- declaration of explicit specializations ------------- */ - -template <> -double -DerivativeApproximation::SecondDerivative<1>::derivative_norm (const Derivative &d); - -template <> -double -DerivativeApproximation::SecondDerivative<2>::derivative_norm (const Derivative &d); - -template <> -double -DerivativeApproximation::SecondDerivative<3>::derivative_norm (const Derivative &d); DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/numerics/derivative_approximation.cc b/deal.II/source/numerics/derivative_approximation.cc index 76ffcdbac4..75a4578557 100644 --- a/deal.II/source/numerics/derivative_approximation.cc +++ b/deal.II/source/numerics/derivative_approximation.cc @@ -45,832 +45,1083 @@ DEAL_II_NAMESPACE_OPEN -template -static inline T sqr (const T t) +namespace { - return t*t; + template + inline T sqr (const T t) + { + return t*t; + } } +// --------------- First the classes and functions that describe individual +// --------------- derivatives +namespace DerivativeApproximation +{ + namespace internal + { + /** + * The following class is used to describe the data needed to compute the + * finite difference approximation to the gradient on a cell. See the + * general documentation of this class for more information on + * implementational details. + * + * @author Wolfgang Bangerth, 2000 + */ + template + class Gradient + { + 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<1,dim> Derivative; + + /** + * Likewise declare the data type that holds the derivative projected to a + * certain directions. + */ + typedef double 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 gradient, we + * choose the Euclidian norm of the gradient vector. + */ + 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 first derivatives, no such thing is necessary, so this function + * is a no-op. + */ + static void symmetrize (Derivative &derivative_tensor); + }; -// static variables -template -const UpdateFlags DerivativeApproximation::Gradient::update_flags = update_values; - -template -const UpdateFlags DerivativeApproximation::SecondDerivative::update_flags = update_gradients; + // static variables + template + const UpdateFlags Gradient::update_flags = update_values; -template -const UpdateFlags DerivativeApproximation::ThirdDerivative::update_flags = update_hessians; + template + template + inline + typename Gradient::ProjectedDerivative + Gradient:: + 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_values (solution, values); + return values[0]; + } + else + { + std::vector > values + (1, Vector(fe_values.get_fe().n_components())); + fe_values.get_function_values (solution, values); + return values[0](component); + } + } -// Dummy structures and dummy function used for WorkStream -namespace internal -{ - namespace Assembler - { - struct Scratch - { - Scratch() {} - }; - struct CopyData + template + inline + double + Gradient::derivative_norm (const Derivative &d) { - CopyData() {} - }; - } -} + double s = 0; + for (unsigned int i=0; i -template -inline -typename DerivativeApproximation::Gradient::ProjectedDerivative -DerivativeApproximation::Gradient:: -get_projected_derivative (const FEValues &fe_values, - const InputVector &solution, - const unsigned int component) -{ - if (fe_values.get_fe().n_components() == 1) + template + inline + void + Gradient::symmetrize (Derivative &) { - std::vector values (1); - fe_values.get_function_values (solution, values); - return values[0]; + // nothing to do here } - else + + + + /** + * The following class is used to describe the data needed to compute the + * finite difference approximation to the second derivatives on a cell. See + * the general documentation of this class for more information on + * implementational details. + * + * @author Wolfgang Bangerth, 2000 + */ + template + class SecondDerivative { - std::vector > values - (1, Vector(fe_values.get_fe().n_components())); - fe_values.get_function_values (solution, values); - return values[0](component); - } -} + 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<2,dim> Derivative; + + /** + * Likewise declare the data type that holds the derivative projected to a + * certain directions. + */ + typedef Tensor<1,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 + const UpdateFlags SecondDerivative::update_flags = update_gradients; -template -inline -double -DerivativeApproximation::Gradient::derivative_norm (const Derivative &d) -{ - double s = 0; - for (unsigned int i=0; i + template + inline + typename SecondDerivative::ProjectedDerivative + SecondDerivative:: + 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_gradients (solution, values); + return values[0]; + } + else + { + std::vector > values + (1, std::vector(fe_values.get_fe().n_components())); + fe_values.get_function_gradients (solution, values); + return values[0][component]; + }; + } -template -inline -void -DerivativeApproximation::Gradient::symmetrize (Derivative &) -{ - // nothing to do here -} + template <> + inline + double + SecondDerivative<1>:: + derivative_norm (const Derivative &d) + { + return std::fabs (d[0][0]); + } -template -template -inline -typename DerivativeApproximation::SecondDerivative::ProjectedDerivative -DerivativeApproximation::SecondDerivative:: -get_projected_derivative (const FEValues &fe_values, - const InputVector &solution, - const unsigned int component) -{ - if (fe_values.get_fe().n_components() == 1) + template <> + inline + double + SecondDerivative<2>:: + derivative_norm (const Derivative &d) { - std::vector values (1); - fe_values.get_function_gradients (solution, values); - return values[0]; + // note that d should be a + // symmetric 2x2 tensor, so the + // eigenvalues are: + // + // 1/2(a+b\pm\sqrt((a-b)^2+4c^2)) + // + // if the d_11=a, d_22=b, + // d_12=d_21=c + const double radicand = dealii::sqr(d[0][0] - d[1][1]) + + 4*dealii::sqr(d[0][1]); + const double eigenvalues[2] + = { 0.5*(d[0][0] + d[1][1] + std::sqrt(radicand)), + 0.5*(d[0][0] + d[1][1] - std::sqrt(radicand)) + }; + + return std::max (std::fabs (eigenvalues[0]), + std::fabs (eigenvalues[1])); } - else + + + + template <> + inline + double + SecondDerivative<3>:: + derivative_norm (const Derivative &d) { - std::vector > values - (1, std::vector(fe_values.get_fe().n_components())); - fe_values.get_function_gradients (solution, values); - return values[0][component]; - }; -} + /* + compute the three eigenvalues of the tensor @p{d} and take the + largest. one could use the following maple script to generate C + code: + + with(linalg); + readlib(C); + A:=matrix(3,3,[[a00,a01,a02],[a01,a11,a12],[a02,a12,a22]]); + E:=eigenvals(A); + EE:=vector(3,[E[1],E[2],E[3]]); + C(EE); + + Unfortunately, with both optimized and non-optimized output, at some + places the code `sqrt(-1.0)' is emitted, and I don't know what + Maple intends to do with it. This happens both with Maple4 and + Maple5. + + Fortunately, Roger Young provided the following Fortran code, which + is transcribed below to C. The code uses an algorithm that uses the + invariants of a symmetric matrix. (The translated algorithm is + augmented by a test for R>0, since R==0 indicates that all three + eigenvalues are equal.) + + + PROGRAM MAIN + + C FIND EIGENVALUES OF REAL SYMMETRIC MATRIX + C (ROGER YOUNG, 2001) + + IMPLICIT NONE + + REAL*8 A11, A12, A13, A22, A23, A33 + REAL*8 I1, J2, J3, AM + REAL*8 S11, S12, S13, S22, S23, S33 + REAL*8 SS12, SS23, SS13 + REAL*8 R,R3, XX,YY, THETA + REAL*8 A1,A2,A3 + REAL*8 PI + PARAMETER (PI=3.141592653587932384D0) + REAL*8 A,B,C, TOL + PARAMETER (TOL=1.D-14) + + C DEFINE A TEST MATRIX + + A11 = -1.D0 + A12 = 5.D0 + A13 = 3.D0 + A22 = -2.D0 + A23 = 0.5D0 + A33 = 4.D0 + + + I1 = A11 + A22 + A33 + AM = I1/3.D0 + + S11 = A11 - AM + S22 = A22 - AM + S33 = A33 - AM + S12 = A12 + S13 = A13 + S23 = A23 + + SS12 = S12*S12 + SS23 = S23*S23 + SS13 = S13*S13 + + J2 = S11*S11 + S22*S22 + S33*S33 + J2 = J2 + 2.D0*(SS12 + SS23 + SS13) + J2 = J2/2.D0 + + J3 = S11**3 + S22**3 + S33**3 + J3 = J3 + 3.D0*S11*(SS12 + SS13) + J3 = J3 + 3.D0*S22*(SS12 + SS23) + J3 = J3 + 3.D0*S33*(SS13 + SS23) + J3 = J3 + 6.D0*S12*S23*S13 + J3 = J3/3.D0 + + R = SQRT(4.D0*J2/3.D0) + R3 = R*R*R + XX = 4.D0*J3/R3 + + YY = 1.D0 - DABS(XX) + IF(YY.LE.0.D0)THEN + IF(YY.GT.(-TOL))THEN + WRITE(6,*)'Equal roots: XX= ',XX + A = -(XX/DABS(XX))*SQRT(J2/3.D0) + B = AM + A + C = AM - 2.D0*A + WRITE(6,*)B,' (twice) ',C + STOP + ELSE + WRITE(6,*)'Error: XX= ',XX + STOP + ENDIF + ENDIF + + THETA = (ACOS(XX))/3.D0 + + A1 = AM + R*COS(THETA) + A2 = AM + R*COS(THETA + 2.D0*PI/3.D0) + A3 = AM + R*COS(THETA + 4.D0*PI/3.D0) + + WRITE(6,*)A1,A2,A3 + + STOP + END + + */ + + const double am = trace(d) / 3.; + + // s := d - trace(d) I + Tensor<2,3> s = d; + for (unsigned int i=0; i<3; ++i) + s[i][i] -= am; + + const double ss01 = s[0][1] * s[0][1], + ss12 = s[1][2] * s[1][2], + ss02 = s[0][2] * s[0][2]; + + const double J2 = (s[0][0]*s[0][0] + s[1][1]*s[1][1] + s[2][2]*s[2][2] + + 2 * (ss01 + ss02 + ss12)) / 2.; + const double J3 = (std::pow(s[0][0],3) + std::pow(s[1][1],3) + std::pow(s[2][2],3) + + 3. * s[0][0] * (ss01 + ss02) + + 3. * s[1][1] * (ss01 + ss12) + + 3. * s[2][2] * (ss02 + ss12) + + 6. * s[0][1] * s[0][2] * s[1][2]) / 3.; + + const double R = std::sqrt (4. * J2 / 3.); + + double EE[3] = { 0, 0, 0 }; + // the eigenvalues are away from + // @p{am} in the order of R. thus, + // if R< -1e-14, ExcInternalError()); + + if (YY < 0) + { + // two roots are equal + const double a = (XX>0 ? -1. : 1.) * R / 2; + EE[0] = EE[1] = am + a; + EE[2] = am - 2.*a; + } + else + { + const double theta = std::acos(XX) / 3.; + EE[0] = am + R*std::cos(theta); + EE[1] = am + R*std::cos(theta + 2./3.*numbers::PI); + EE[2] = am + R*std::cos(theta + 4./3.*numbers::PI); + }; + }; + return std::max (std::fabs (EE[0]), + std::max (std::fabs (EE[1]), + std::fabs (EE[2]))); + } -template <> -inline -double -DerivativeApproximation::SecondDerivative<1>:: -derivative_norm (const Derivative &d) -{ - return std::fabs (d[0][0]); -} + template + inline + double + SecondDerivative:: + derivative_norm (const Derivative &) + { + // computing the spectral norm is + // not so simple in general. it is + // feasible for dim==3 as shown + // above, since then there are + // still closed form expressions of + // the roots of the characteristic + // polynomial, and they can easily + // be computed using + // maple. however, for higher + // dimensions, some other method + // needs to be employed. maybe some + // steps of the power method would + // suffice? + Assert (false, ExcNotImplemented()); + return 0; + } -template <> -inline -double -DerivativeApproximation::SecondDerivative<2>:: -derivative_norm (const Derivative &d) -{ - // note that d should be a - // symmetric 2x2 tensor, so the - // eigenvalues are: - // - // 1/2(a+b\pm\sqrt((a-b)^2+4c^2)) - // - // if the d_11=a, d_22=b, - // d_12=d_21=c - const double radicand = dealii::sqr(d[0][0] - d[1][1]) + - 4*dealii::sqr(d[0][1]); - const double eigenvalues[2] - = { 0.5*(d[0][0] + d[1][1] + std::sqrt(radicand)), - 0.5*(d[0][0] + d[1][1] - std::sqrt(radicand)) - }; - return std::max (std::fabs (eigenvalues[0]), - std::fabs (eigenvalues[1])); -} + template + inline + void + SecondDerivative::symmetrize (Derivative &d) + { + // symmetrize non-diagonal entries + for (unsigned int i=0; i -inline -double -DerivativeApproximation::SecondDerivative<3>:: -derivative_norm (const Derivative &d) -{ - /* - compute the three eigenvalues of the tensor @p{d} and take the - largest. one could use the following maple script to generate C - code: - - with(linalg); - readlib(C); - A:=matrix(3,3,[[a00,a01,a02],[a01,a11,a12],[a02,a12,a22]]); - E:=eigenvals(A); - EE:=vector(3,[E[1],E[2],E[3]]); - C(EE); - - Unfortunately, with both optimized and non-optimized output, at some - places the code `sqrt(-1.0)' is emitted, and I don't know what - Maple intends to do with it. This happens both with Maple4 and - Maple5. - - Fortunately, Roger Young provided the following Fortran code, which - is transcribed below to C. The code uses an algorithm that uses the - invariants of a symmetric matrix. (The translated algorithm is - augmented by a test for R>0, since R==0 indicates that all three - eigenvalues are equal.) - - - PROGRAM MAIN - - C FIND EIGENVALUES OF REAL SYMMETRIC MATRIX - C (ROGER YOUNG, 2001) - - IMPLICIT NONE - - REAL*8 A11, A12, A13, A22, A23, A33 - REAL*8 I1, J2, J3, AM - REAL*8 S11, S12, S13, S22, S23, S33 - REAL*8 SS12, SS23, SS13 - REAL*8 R,R3, XX,YY, THETA - REAL*8 A1,A2,A3 - REAL*8 PI - PARAMETER (PI=3.141592653587932384D0) - REAL*8 A,B,C, TOL - PARAMETER (TOL=1.D-14) - - C DEFINE A TEST MATRIX - - A11 = -1.D0 - A12 = 5.D0 - A13 = 3.D0 - A22 = -2.D0 - A23 = 0.5D0 - A33 = 4.D0 - - - I1 = A11 + A22 + A33 - AM = I1/3.D0 - - S11 = A11 - AM - S22 = A22 - AM - S33 = A33 - AM - S12 = A12 - S13 = A13 - S23 = A23 - - SS12 = S12*S12 - SS23 = S23*S23 - SS13 = S13*S13 - - J2 = S11*S11 + S22*S22 + S33*S33 - J2 = J2 + 2.D0*(SS12 + SS23 + SS13) - J2 = J2/2.D0 - - J3 = S11**3 + S22**3 + S33**3 - J3 = J3 + 3.D0*S11*(SS12 + SS13) - J3 = J3 + 3.D0*S22*(SS12 + SS23) - J3 = J3 + 3.D0*S33*(SS13 + SS23) - J3 = J3 + 6.D0*S12*S23*S13 - J3 = J3/3.D0 - - R = SQRT(4.D0*J2/3.D0) - R3 = R*R*R - XX = 4.D0*J3/R3 - - YY = 1.D0 - DABS(XX) - IF(YY.LE.0.D0)THEN - IF(YY.GT.(-TOL))THEN - WRITE(6,*)'Equal roots: XX= ',XX - A = -(XX/DABS(XX))*SQRT(J2/3.D0) - B = AM + A - C = AM - 2.D0*A - WRITE(6,*)B,' (twice) ',C - STOP - ELSE - WRITE(6,*)'Error: XX= ',XX - STOP - ENDIF - ENDIF - - THETA = (ACOS(XX))/3.D0 - - A1 = AM + R*COS(THETA) - A2 = AM + R*COS(THETA + 2.D0*PI/3.D0) - A3 = AM + R*COS(THETA + 4.D0*PI/3.D0) - - WRITE(6,*)A1,A2,A3 - - STOP - END - - */ - - const double am = trace(d) / 3.; - - // s := d - trace(d) I - Tensor<2,3> s = d; - for (unsigned int i=0; i<3; ++i) - s[i][i] -= am; - - const double ss01 = s[0][1] * s[0][1], - ss12 = s[1][2] * s[1][2], - ss02 = s[0][2] * s[0][2]; - - const double J2 = (s[0][0]*s[0][0] + s[1][1]*s[1][1] + s[2][2]*s[2][2] - + 2 * (ss01 + ss02 + ss12)) / 2.; - const double J3 = (std::pow(s[0][0],3) + std::pow(s[1][1],3) + std::pow(s[2][2],3) - + 3. * s[0][0] * (ss01 + ss02) - + 3. * s[1][1] * (ss01 + ss12) - + 3. * s[2][2] * (ss02 + ss12) - + 6. * s[0][1] * s[0][2] * s[1][2]) / 3.; - - const double R = std::sqrt (4. * J2 / 3.); - - double EE[3] = { 0, 0, 0 }; - // the eigenvalues are away from - // @p{am} in the order of R. thus, - // if R< + class ThirdDerivative { - // at least two eigenvalues are - // distinct - const double R3 = R*R*R; - const double XX = 4. * J3 / R3; - const double YY = 1. - std::fabs(XX); + 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); + }; - Assert (YY > -1e-14, ExcInternalError()); + template + const UpdateFlags ThirdDerivative::update_flags = update_hessians; - if (YY < 0) + + template + template + inline + typename ThirdDerivative::ProjectedDerivative + ThirdDerivative:: + get_projected_derivative (const FEValues &fe_values, + const InputVector &solution, + const unsigned int component) { - // two roots are equal - const double a = (XX>0 ? -1. : 1.) * R / 2; - EE[0] = EE[1] = am + a; - EE[2] = am - 2.*a; + if (fe_values.get_fe().n_components() == 1) + { + std::vector values (1); + fe_values.get_function_hessians (solution, values); + return values[0]; + } + else + { + std::vector > values + (1, std::vector(fe_values.get_fe().n_components())); + fe_values.get_function_hessians (solution, values); + return values[0][component]; + }; } - else - { - const double theta = std::acos(XX) / 3.; - EE[0] = am + R*std::cos(theta); - EE[1] = am + R*std::cos(theta + 2./3.*numbers::PI); - EE[2] = am + R*std::cos(theta + 4./3.*numbers::PI); - }; - }; - - return std::max (std::fabs (EE[0]), - std::max (std::fabs (EE[1]), - std::fabs (EE[2]))); -} -template -inline -double -DerivativeApproximation::SecondDerivative:: -derivative_norm (const Derivative &) -{ - // computing the spectral norm is - // not so simple in general. it is - // feasible for dim==3 as shown - // above, since then there are - // still closed form expressions of - // the roots of the characteristic - // polynomial, and they can easily - // be computed using - // maple. however, for higher - // dimensions, some other method - // needs to be employed. maybe some - // steps of the power method would - // suffice? - Assert (false, ExcNotImplemented()); - return 0; -} + template <> + inline + double + ThirdDerivative<1>:: + derivative_norm (const Derivative &d) + { + return std::fabs (d[0][0][0]); + } -template -inline -void -DerivativeApproximation::SecondDerivative::symmetrize (Derivative &d) -{ - // symmetrize non-diagonal entries - for (unsigned int i=0; i + inline + double + ThirdDerivative:: + derivative_norm (const Derivative &d) + { + // return the Frobenius-norm. this is a + // member function of Tensor + return d.norm(); + } -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) + template + inline + void + ThirdDerivative::symmetrize (Derivative &d) { - std::vector values (1); - fe_values.get_function_hessians (solution, values); - return values[0]; + // 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 DerivativeSelector { - std::vector > values - (1, std::vector(fe_values.get_fe().n_components())); - fe_values.get_function_hessians (solution, values); - return values[0][component]; + 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 <> -inline -double -DerivativeApproximation::ThirdDerivative<1>:: -derivative_norm (const Derivative &d) -{ - return std::fabs (d[0][0][0]); -} + template + class DerivativeSelector<2,dim> + { + public: + typedef SecondDerivative DerivDescr; + }; + template + class DerivativeSelector<3,dim> + { + public: -template -inline -double -DerivativeApproximation::ThirdDerivative:: -derivative_norm (const Derivative &d) -{ - // return the Frobenius-norm. this is a - // member function of Tensor - return d.norm(); + typedef ThirdDerivative DerivDescr; + }; + } } - -template -inline -void -DerivativeApproximation::ThirdDerivative::symmetrize (Derivative &d) +// Dummy structures and dummy function used for WorkStream +namespace DerivativeApproximation { - // 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, int spacedim> -void -DerivativeApproximation:: -approximate_gradient (const Mapping &mapping, - const DH &dof_handler, - const InputVector &solution, - Vector &derivative_norm, - const unsigned int component) +namespace DerivativeApproximation { - approximate_derivative,dim> (mapping, - dof_handler, - solution, - component, - derivative_norm); -} + namespace internal + { + /** + * Compute the derivative approximation on one cell. This computes the full + * derivative tensor. + */ + template class DH, class InputVector, int spacedim> + void + 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 + // single quadratures, mappings, + // and finite elements. if we have + // an hp DoFHandler, + // dof_handler.get_fe() returns a + // collection of which we do a + // shallow copy instead + const hp::QCollection q_collection (midpoint_rule); + const hp::FECollection fe_collection(dof_handler.get_fe()); + const hp::MappingCollection mapping_collection (mapping); + + hp::FEValues x_fe_midpoint_value (mapping_collection, fe_collection, + q_collection, + DerivativeDescription::update_flags | + update_quadrature_points); + + // matrix Y=sum_i y_i y_i^T + Tensor<2,dim> Y; + + + // vector to hold iterators to all + // active neighbors of a cell + // reserve the maximal number of + // active neighbors + std::vector::active_cell_iterator> active_neighbors; + active_neighbors.reserve (GeometryInfo::faces_per_cell * + GeometryInfo::max_children_per_face); + + // vector + // g=sum_i y_i (f(x+y_i)-f(x))/|y_i| + // or related type for higher + // derivatives + typename DerivativeDescription::Derivative projected_derivative; + + // reinit fe values object... + x_fe_midpoint_value.reinit (cell); + const FEValues &fe_midpoint_value + = x_fe_midpoint_value.get_present_fe_values(); + + // ...and get the value of the + // projected derivative... + const typename DerivativeDescription::ProjectedDerivative + this_midpoint_value + = DerivativeDescription::get_projected_derivative (fe_midpoint_value, + solution, + component); + // ...and the place where it lives + const Point this_center = fe_midpoint_value.quadrature_point(0); + + // loop over all neighbors and + // accumulate the difference + // quotients from them. note + // that things get a bit more + // complicated if the neighbor + // is more refined than the + // present one + // + // to make processing simpler, + // first collect all neighbor + // cells in a vector, and then + // collect the data from them + GridTools::get_active_neighbors >(cell, active_neighbors); + + // now loop over all active + // neighbors and collect the + // data we need + typename std::vector::active_cell_iterator>::const_iterator + neighbor_ptr = active_neighbors.begin(); + for (; neighbor_ptr!=active_neighbors.end(); ++neighbor_ptr) + { + const typename DH::active_cell_iterator + neighbor = *neighbor_ptr; + + // reinit fe values object... + x_fe_midpoint_value.reinit (neighbor); + const FEValues &fe_midpoint_value + = x_fe_midpoint_value.get_present_fe_values(); + + // ...and get the value of the + // solution... + const typename DerivativeDescription::ProjectedDerivative + neighbor_midpoint_value + = DerivativeDescription::get_projected_derivative (fe_midpoint_value, + solution, component); + + // ...and the place where it lives + const Point + neighbor_center = fe_midpoint_value.quadrature_point(0); + + + // vector for the + // normalized + // direction between + // the centers of two + // cells + Point y = neighbor_center - this_center; + const double distance = std::sqrt(y.square()); + // normalize y + y /= distance; + // *** note that unlike in + // the docs, y denotes the + // normalized vector + // connecting the centers + // of the two cells, rather + // than the normal + // difference! *** + + // add up the + // contribution of + // this cell to Y + for (unsigned int i=0; i Y_inverse = invert(Y); + + contract (derivative, Y_inverse, projected_derivative); + + // finally symmetrize the derivative + DerivativeDescription::symmetrize (derivative); + } -template class DH, class InputVector, int spacedim> -void -DerivativeApproximation:: -approximate_gradient (const DH &dof_handler, - const InputVector &solution, - Vector &derivative_norm, - const unsigned int component) -{ - approximate_derivative,dim> (StaticMappingQ1::mapping, - dof_handler, - solution, - component, - derivative_norm); -} + /** + * Compute the derivative approximation on a given cell. Fill the @p + * derivative_norm vector with the norm of the computed derivative tensors + * on the cell. + */ + template class DH, class InputVector, int spacedim> + void + approximate (SynchronousIterators::active_cell_iterator,Vector::iterator> > const &cell, + const Mapping &mapping, + const DH &dof_handler, + const InputVector &solution, + const unsigned int component) + { + // if the cell is not locally owned, then there is nothing to do + if (std_cxx1x::get<0>(cell.iterators)->is_locally_owned() == false) + *std_cxx1x::get<1>(cell.iterators) = 0; + else + { + typename DerivativeDescription::Derivative derivative; + // call the function doing the actual + // work on this cell + approximate_cell + (mapping,dof_handler,solution,component,std_cxx1x::get<0>(cell.iterators),derivative); + // evaluate the norm and fill the vector + //*derivative_norm_on_this_cell + *std_cxx1x::get<1>(cell.iterators) = DerivativeDescription::derivative_norm (derivative); + } + } -template class DH, class InputVector, int spacedim> -void -DerivativeApproximation:: -approximate_second_derivative (const Mapping &mapping, - const DH &dof_handler, - const InputVector &solution, - Vector &derivative_norm, - const unsigned int component) -{ - approximate_derivative,dim> (mapping, - dof_handler, - solution, - component, - derivative_norm); -} + /** + * Kind of the main function of this class. It is called by the public entry + * points to this class with the correct template first argument and then + * simply calls the @p approximate function, after setting up several + * threads and doing some administration that is independent of the actual + * derivative to be computed. + * + * The @p component argument denotes which component of the solution vector + * we are to work on. + */ + template class DH, class InputVector, int spacedim> + void + approximate_derivative (const Mapping &mapping, + const DH &dof_handler, + const InputVector &solution, + const unsigned int component, + Vector &derivative_norm) + { + Assert (derivative_norm.size() == dof_handler.get_tria().n_active_cells(), + ExcInvalidVectorLength (derivative_norm.size(), + dof_handler.get_tria().n_active_cells())); + Assert (component < dof_handler.get_fe().n_components(), + ExcIndexRange (component, 0, dof_handler.get_fe().n_components())); + + typedef std_cxx1x::tuple::active_cell_iterator,Vector::iterator> + Iterators; + SynchronousIterators begin(Iterators(dof_handler.begin_active(), + derivative_norm.begin())), + end(Iterators(dof_handler.end(), + derivative_norm.end())); + + // There is no need for a copier because there is no conflict between threads + // to write in derivative_norm. Scratch and CopyData are also useless. + WorkStream::run(begin, + end, + static_cast const &, + Assembler::Scratch const &, Assembler::CopyData &)> > + (std_cxx1x::bind(&approximate, + std_cxx1x::_1, + std_cxx1x::cref(mapping), + std_cxx1x::cref(dof_handler), + std_cxx1x::cref(solution),component)), + std_cxx1x::function (), + internal::Assembler::Scratch (),internal::Assembler::CopyData ()); + } -template class DH, class InputVector, int spacedim> -void -DerivativeApproximation:: -approximate_second_derivative (const DH &dof_handler, - const InputVector &solution, - Vector &derivative_norm, - const unsigned int component) -{ - approximate_derivative,dim> (StaticMappingQ1::mapping, - dof_handler, - solution, - component, - derivative_norm); -} + } // namespace internal +} // namespace DerivativeApproximation -template class DH, class InputVector, int order, int spacedim> -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) + +// ------------------------ finally for the public interface of this namespace + +namespace DerivativeApproximation { - approximate_cell::DerivDescr,dim,DH,InputVector> - (mapping, - dof, - solution, - component, - cell, - derivative); -} + template class DH, class InputVector, int spacedim> + void + approximate_gradient (const Mapping &mapping, + const DH &dof_handler, + const InputVector &solution, + Vector &derivative_norm, + const unsigned int component) + { + internal::approximate_derivative,dim> (mapping, + dof_handler, + solution, + component, + derivative_norm); + } + template class DH, class InputVector, int spacedim> + void + approximate_gradient (const DH &dof_handler, + const InputVector &solution, + Vector &derivative_norm, + const unsigned int component) + { + internal::approximate_derivative,dim> (StaticMappingQ1::mapping, + dof_handler, + solution, + component, + derivative_norm); + } -template class DH, class InputVector, int order, int spacedim> -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 - approximate_derivative_tensor - (StaticMappingQ1::mapping, - dof, - solution, - cell, - derivative, - component); -} + template class DH, class InputVector, int spacedim> + void + approximate_second_derivative (const Mapping &mapping, + const DH &dof_handler, + const InputVector &solution, + Vector &derivative_norm, + const unsigned int component) + { + internal::approximate_derivative,dim> (mapping, + dof_handler, + solution, + component, + derivative_norm); + } -template class DH, class InputVector, int spacedim> -void -DerivativeApproximation:: -approximate_derivative (const Mapping &mapping, - const DH &dof_handler, - const InputVector &solution, - const unsigned int component, - Vector &derivative_norm) -{ - Assert (derivative_norm.size() == dof_handler.get_tria().n_active_cells(), - ExcInvalidVectorLength (derivative_norm.size(), - dof_handler.get_tria().n_active_cells())); - Assert (component < dof_handler.get_fe().n_components(), - ExcIndexRange (component, 0, dof_handler.get_fe().n_components())); - - typedef std_cxx1x::tuple::active_cell_iterator,Vector::iterator> - Iterators; - SynchronousIterators begin(Iterators(dof_handler.begin_active(), - derivative_norm.begin())), - end(Iterators(dof_handler.end(), - derivative_norm.end())); - - // There is no need for a copier because there is no conflict between threads - // to write in derivative_norm. Scratch and CopyData are also useless. - WorkStream::run(begin, - end, - static_cast const &, - internal::Assembler::Scratch const &,internal::Assembler::CopyData &)> > - (std_cxx1x::bind(&DerivativeApproximation::template approximate, - std_cxx1x::_1, - std_cxx1x::cref(mapping), - std_cxx1x::cref(dof_handler), - std_cxx1x::cref(solution),component)), - std_cxx1x::function (), - internal::Assembler::Scratch (),internal::Assembler::CopyData ()); -} + template class DH, class InputVector, int spacedim> + void + approximate_second_derivative (const DH &dof_handler, + const InputVector &solution, + Vector &derivative_norm, + const unsigned int component) + { + internal::approximate_derivative,dim> (StaticMappingQ1::mapping, + dof_handler, + solution, + component, + derivative_norm); + } + template class DH, class InputVector, int order, int spacedim> + 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) + { + internal::approximate_cell::DerivDescr,dim,DH,InputVector> + (mapping, + dof, + solution, + component, + cell, + derivative); + } -template class DH, class InputVector, int spacedim> -void -DerivativeApproximation::approximate (SynchronousIterators::active_cell_iterator,Vector::iterator> > const &cell, - const Mapping &mapping, - const DH &dof_handler, - const InputVector &solution, - const unsigned int component) -{ - // if the cell is not locally owned, then there is nothing to do - if (std_cxx1x::get<0>(cell.iterators)->is_locally_owned() == false) - *std_cxx1x::get<1>(cell.iterators) = 0; - else - { - typename DerivativeDescription::Derivative derivative; - // call the function doing the actual - // work on this cell - DerivativeApproximation::template approximate_cell - (mapping,dof_handler,solution,component,std_cxx1x::get<0>(cell.iterators),derivative); - // evaluate the norm and fill the vector - //*derivative_norm_on_this_cell - *std_cxx1x::get<1>(cell.iterators) = DerivativeDescription::derivative_norm (derivative); - } -} + template class DH, class InputVector, int order, int spacedim> + void + 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 + approximate_derivative_tensor + (StaticMappingQ1::mapping, + dof, + solution, + cell, + derivative, + component); + } -template class DH, class InputVector, int spacedim> -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 - // single quadratures, mappings, - // and finite elements. if we have - // an hp DoFHandler, - // dof_handler.get_fe() returns a - // collection of which we do a - // shallow copy instead - const hp::QCollection q_collection (midpoint_rule); - const hp::FECollection fe_collection(dof_handler.get_fe()); - const hp::MappingCollection mapping_collection (mapping); - - hp::FEValues x_fe_midpoint_value (mapping_collection, fe_collection, - q_collection, - DerivativeDescription::update_flags | - update_quadrature_points); - - // matrix Y=sum_i y_i y_i^T - Tensor<2,dim> Y; - - - // vector to hold iterators to all - // active neighbors of a cell - // reserve the maximal number of - // active neighbors - std::vector::active_cell_iterator> active_neighbors; - active_neighbors.reserve (GeometryInfo::faces_per_cell * - GeometryInfo::max_children_per_face); - - // vector - // g=sum_i y_i (f(x+y_i)-f(x))/|y_i| - // or related type for higher - // derivatives - typename DerivativeDescription::Derivative projected_derivative; - - // reinit fe values object... - x_fe_midpoint_value.reinit (cell); - const FEValues &fe_midpoint_value - = x_fe_midpoint_value.get_present_fe_values(); - - // ...and get the value of the - // projected derivative... - const typename DerivativeDescription::ProjectedDerivative - this_midpoint_value - = DerivativeDescription::get_projected_derivative (fe_midpoint_value, - solution, - component); - // ...and the place where it lives - const Point this_center = fe_midpoint_value.quadrature_point(0); - - // loop over all neighbors and - // accumulate the difference - // quotients from them. note - // that things get a bit more - // complicated if the neighbor - // is more refined than the - // present one - // - // to make processing simpler, - // first collect all neighbor - // cells in a vector, and then - // collect the data from them - GridTools::get_active_neighbors >(cell, active_neighbors); - - // now loop over all active - // neighbors and collect the - // data we need - typename std::vector::active_cell_iterator>::const_iterator - neighbor_ptr = active_neighbors.begin(); - for (; neighbor_ptr!=active_neighbors.end(); ++neighbor_ptr) - { - const typename DH::active_cell_iterator - neighbor = *neighbor_ptr; - - // reinit fe values object... - x_fe_midpoint_value.reinit (neighbor); - const FEValues &fe_midpoint_value - = x_fe_midpoint_value.get_present_fe_values(); - - // ...and get the value of the - // solution... - const typename DerivativeDescription::ProjectedDerivative - neighbor_midpoint_value - = DerivativeDescription::get_projected_derivative (fe_midpoint_value, - solution, component); - - // ...and the place where it lives - const Point - neighbor_center = fe_midpoint_value.quadrature_point(0); - - - // vector for the - // normalized - // direction between - // the centers of two - // cells - Point y = neighbor_center - this_center; - const double distance = std::sqrt(y.square()); - // normalize y - y /= distance; - // *** note that unlike in - // the docs, y denotes the - // normalized vector - // connecting the centers - // of the two cells, rather - // than the normal - // difference! *** - - // add up the - // contribution of - // this cell to Y - for (unsigned int i=0; i Y_inverse = invert(Y); - - contract (derivative, Y_inverse, projected_derivative); - - // finally symmetrize the derivative - DerivativeDescription::symmetrize (derivative); -} -template -double -DerivativeApproximation:: -derivative_norm(const Tensor &derivative) -{ - return DerivativeSelector::DerivDescr::derivative_norm(derivative); -} + template + double + derivative_norm (const Tensor &derivative) + { + return internal::DerivativeSelector::DerivDescr::derivative_norm(derivative); + } + +} // --------------------------- explicit instantiations --------------------- diff --git a/deal.II/source/numerics/derivative_approximation.inst.in b/deal.II/source/numerics/derivative_approximation.inst.in index 6e7ed876ab..430c8a095f 100644 --- a/deal.II/source/numerics/derivative_approximation.inst.in +++ b/deal.II/source/numerics/derivative_approximation.inst.in @@ -17,9 +17,10 @@ for (deal_II_dimension : DIMENSIONS ; VEC : SERIAL_VECTORS ; DH : DOFHANDLER_TEMPLATES) { + namespace DerivativeApproximation + \{ template void -DerivativeApproximation:: approximate_gradient (const Mapping &mapping, const DH &dof_handler, @@ -29,7 +30,6 @@ approximate_gradient template void -DerivativeApproximation:: approximate_gradient (const DH &dof_handler, const VEC &solution, @@ -38,7 +38,6 @@ approximate_gradient template void -DerivativeApproximation:: approximate_second_derivative (const Mapping &mapping, const DH &dof_handler, @@ -48,7 +47,6 @@ approximate_second_derivative template void -DerivativeApproximation:: approximate_second_derivative (const DH &dof_handler, const VEC &solution, @@ -57,7 +55,6 @@ approximate_second_derivative template void -DerivativeApproximation:: approximate_derivative_tensor (const Mapping & mapping, const DH &dof_handler, @@ -68,7 +65,6 @@ approximate_derivative_tensor template void -DerivativeApproximation:: approximate_derivative_tensor (const Mapping & mapping, const DH &dof_handler, @@ -79,7 +75,6 @@ approximate_derivative_tensor template void -DerivativeApproximation:: approximate_derivative_tensor (const Mapping & mapping, const DH &dof_handler, @@ -90,7 +85,6 @@ approximate_derivative_tensor template void -DerivativeApproximation:: approximate_derivative_tensor (const DH &dof_handler, const VEC &solution, @@ -100,7 +94,6 @@ approximate_derivative_tensor template void -DerivativeApproximation:: approximate_derivative_tensor (const DH &dof_handler, const VEC &solution, @@ -110,7 +103,6 @@ approximate_derivative_tensor template void -DerivativeApproximation:: approximate_derivative_tensor (const DH &dof_handler, const VEC &solution, @@ -118,40 +110,24 @@ approximate_derivative_tensor Tensor<3,deal_II_dimension> &derivative, const unsigned int component); + \} } for (deal_II_dimension : DIMENSIONS) { + namespace DerivativeApproximation + \{ 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 -// -// on AIX, the linker is unhappy about some missing symbols. they -// should really be there, but explicitly instantiating them will also -// not hurt -template -const UpdateFlags -DerivativeApproximation::Gradient::update_flags; - -template -const UpdateFlags -DerivativeApproximation::SecondDerivative::update_flags; -template -const UpdateFlags -DerivativeApproximation::ThirdDerivative::update_flags; + \} } -- 2.39.5