QuOc

 

anisotropies.h

Go to the documentation of this file.
00001 /* Copyright (c) 2001-2010 AG Rumpf, INS, Universitaet Bonn                      */
00002 /*                                                                               */
00003 /* The contents of this file are subject to the terms of the Common Development  */
00004 /* and Distribution License Version 1.0 (the "License"); you may not use         */
00005 /* this file except in compliance with the License. You may obtain a copy of     */
00006 /* the License at http://www.sun.com/cddl/                                       */
00007 /*                                                                               */
00008 /* Software distributed under the License is distributed on an "AS IS" basis,    */
00009 /* WITHOUT WARRANTY OF ANY KIND, either express or implied.                      */
00010 /*                                                                               */
00011 #ifndef __ANISOTROPIES_H
00012 #define __ANISOTROPIES_H
00013 
00014 #include <aol.h>
00015 #include <polarCoords.h>
00016 #include <scalarArray.h>
00017 
00018 namespace qc {
00019 
00020 //! (Hopefully) efficient graph interface for 3d-anisotropies to be used in the 2d-graph-case.
00021 /*! \author Nemitz */
00022 template <typename RealType, typename Imp>
00023 class Anisotropy3dGraphInterface {
00024 public:
00025   Anisotropy3dGraphInterface(  )  { _graphFlag = false; }
00026 
00027   //! call this method if you want to use the graph-methods.
00028   void setGraphFlag() { _graphFlag = true; }
00029   bool isGraphCase() const { return( _graphFlag ); }
00030 
00031   //! This method is for the graph-case. Here the normal is (-grad u, 1).
00032   //! Notice: To use this method you have to set the graph-flag before!!
00033   RealType gammaNorm ( const aol::Vec2<RealType> &z ) const {
00034     if ( !_graphFlag )
00035       throw aol::Exception ( "\nYou have called the gammaNorm-Method with a 2d-vector.\nSet the graph-flag before!\n", __FILE__, __LINE__ );
00036     aol::Vec3<RealType> z3 ( z[0], z[1], -1 ); 
00037     return this->asImp().gammaNorm ( z3 );
00038   }
00039 
00040   //! This method is for the graph-case. Here the normal is (-grad u, 1).
00041   //! Notice: To use this method you have to set the graph-flag before!!
00042   void gammaFirstDerivative ( const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
00043     if ( !_graphFlag ) {
00044       throw aol::Exception ( "\nYou have called the gammaFirstDerivative-Method with a 2d-vector.\nSet the graph-flag before!\n", __FILE__, __LINE__ );
00045     }
00046     aol::Vec3<RealType> z3 ( z[0], z[1], -1 );
00047     aol::Vec3<RealType> v3;
00048     this->asImp().gammaFirstDerivative ( z3, v3 );
00049     v[0] =  v3[0];     // copy the first two arguments of v3 
00050     v[1] =  v3[1];    
00051   }
00052 
00053   void gammaFirstDerivative ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
00054     gammaFirstDerivative( z, v );
00055   }
00056   void gammaFirstDerivative ( const qc::Element &, const int &, const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
00057     gammaFirstDerivative( z, v );
00058   }
00059   
00060   
00061   void implicitPart ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z,
00062                       aol::Matrix22<RealType> &mat ) const {
00063     implicitPart ( z, mat );
00064   }
00065   void implicitPart ( const qc::Element &, const int &, const aol::Vec2<RealType> &z,
00066                       aol::Matrix22<RealType> &mat ) const {
00067                       implicitPart ( z, mat );
00068   }
00069                       
00070   void explicitPart ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z,
00071                       aol::Matrix22<RealType> &mat ) const {
00072     explicitPart ( z, mat );
00073   }
00074 
00075   //! second derivative, implicit part of the p-regularized willmore flow (graph case)
00076   void implicitPart ( const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
00077     aol::Vec3<RealType> z3 ( z[0], z[1], -1 ); 
00078     aol::Matrix33<RealType> mat33;
00079     this->asImp().implicitPart( z3, mat33 );
00080     mat[0][0] = mat33[0][0];
00081     mat[0][1] = mat33[0][1];
00082     mat[1][0] = mat33[1][0];
00083     mat[1][1] = mat33[1][1];
00084   }
00085   
00086   //! second derivative, explicit part of the p-regularized willmore flow (graph case)
00087   void explicitPart ( const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
00088     aol::Vec3<RealType> z3 ( z[0], z[1], -1 ); 
00089     aol::Matrix33<RealType> mat33;
00090     this->asImp().implicitPart( z3, mat33 );
00091     mat[0][0] = mat33[0][0];
00092     mat[0][1] = mat33[0][1];
00093     mat[1][0] = mat33[1][0];
00094     mat[1][1] = mat33[1][1];
00095   }
00096   
00097   
00098   
00099 protected:
00100   // barton-nackman
00101   inline Imp& asImp() { return static_cast<Imp&> ( *this ); }
00102   inline const Imp& asImp() const { return static_cast<const Imp&> ( *this ); }
00103 
00104   bool  _graphFlag;             //! set this to use the gamma-methods with a 2d-vector-argument (graph-case)
00105 };
00106 
00107 
00108 
00109 //! General class, that applies a linear transformation (2x2-matrix) to a given anisotropy.
00110 template <typename RealType, typename AnisoType>
00111 class LinearTransformedAnisotropy2d {
00112 private:
00113   AnisoType &_anisotropy;
00114   aol::Matrix22<RealType> _mat;              //! the linear transformation matrix
00115   aol::Matrix22<RealType> _matTransposed;
00116 
00117 public:
00118   LinearTransformedAnisotropy2d ( AnisoType &Anisotropy, aol::Matrix22<RealType> Mat ) :
00119       _anisotropy ( Anisotropy ), _mat ( Mat ), _matTransposed ( Mat ) { _matTransposed.transpose(); }
00120 
00121 
00122   //! Evaluating the function itself, i.e. the regularization of
00123   //! \f$ \gamma(z) = \sqrt{z_1^2 + \epsilon^2} + \sqrt{z_1^2 + \epsilon^2} \f$.
00124   RealType gammaNorm ( const aol::Vec2<RealType> &z ) const {
00125     aol::Vec2<RealType> v;
00126     _mat.mult ( z, v );
00127     return ( _anisotropy.gammaNorm ( v ) );
00128   }
00129 
00130   //! evaluating the first derivative: \f$ ... \f$
00131   void gammaFirstDerivative ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
00132     gammaFirstDerivative ( z, v );
00133   }
00134 
00135   //! evaluating the first derivative: \f$ ... \f$
00136   void gammaFirstDerivative ( const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
00137     aol::Vec2<RealType> temp1 ( 0., 0. );
00138     aol::Vec2<RealType> temp2 ( 0., 0. );
00139     _mat.mult ( z, temp1 );
00140     _anisotropy.gammaFirstDerivative ( temp1, temp2 );
00141     _matTransposed.mult ( temp2, v );
00142 //     _anisotropy.gammaFirstDerivative ( temp1, v );      // HACK
00143   }
00144 };
00145 
00146 
00147 
00148 //! Class, that generates a 3d-anisotropy by rotating a 2d-anisotropy
00149 //! around the z-axis. For the derivatives, the 2d-anisotropy has to
00150 //! provide partial derivatives wrt r and phi in the methods
00151 //! partialDerivativesWRTPhiR( const z, pd& )
00152 template <typename RealType, typename AnisoType>
00153 class Rotated3dAnisotropy : public Anisotropy3dGraphInterface<RealType, Rotated3dAnisotropy<RealType, AnisoType> > {
00154 public:  
00155   using Anisotropy3dGraphInterface<RealType, Rotated3dAnisotropy<RealType, AnisoType> >::gammaNorm;
00156   using Anisotropy3dGraphInterface<RealType, Rotated3dAnisotropy<RealType, AnisoType> >::gammaFirstDerivative;
00157 
00158 private:
00159   const AnisoType &_anisotropy;
00160 
00161   const RealType _delta;     //! Regularization-parameter for the radius
00162   RealType _deltaSqr;        //! uh, don't know ...
00163 
00164 public:
00165   Rotated3dAnisotropy ( const AnisoType &Anisotropy, const RealType delta ) :
00166       Anisotropy3dGraphInterface<RealType, Rotated3dAnisotropy<RealType, AnisoType> > ( ),
00167       _anisotropy ( Anisotropy ),
00168       _delta ( delta ) { _deltaSqr = _delta * _delta; }
00169 
00170   //! Evaluating the function itself
00171   RealType gammaNorm ( const aol::Vec3<RealType> &z ) const {
00172     aol::Vec2<RealType> v ( sqrt ( z[0]*z[0] + z[1]*z[1] + _deltaSqr ), z[2] );
00173     return ( _anisotropy.gammaNorm ( v ) );
00174   }
00175 
00176   //! evaluating the first derivative: \f$ ... \f$
00177   void gammaFirstDerivative ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v )
00178   const {    gammaFirstDerivative ( z, v );  }
00179 
00180   //! evaluating the first derivative: \f$ ... \f$
00181   void gammaFirstDerivative ( const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
00182     RealType r   = sqrt ( z[0]*z[0] + z[1]*z[1] + _deltaSqr );
00183 
00184     // compute the partial derivatives of the 2d anisotropy
00185     aol::Vec2<RealType> z2Arg( r, z[2] );
00186     aol::Vec2<RealType> pdGamma;
00187     
00188     _anisotropy.gammaFirstDerivative ( z2Arg, pdGamma );
00189 
00190     RealType drdx = z[0] / r;       //! partial derivatives of r(x,y,z)
00191     RealType drdy = z[1] / r;
00192                 
00193     v[0] = pdGamma[0] * drdx; 
00194     v[1] = pdGamma[0] * drdy; 
00195     v[2] = pdGamma[1];
00196   }             
00197 };
00198 
00199 
00200 //! Class for blending two different anisotropies dependent on the image value, e.g. value>0.5 => use
00201 //! anisotropy1 else use anisotropy2 and use a linear blending in an interval between.
00202 template <typename ConfiguratorType, typename AnisoType1, typename AnisoType2>
00203 class LinearBlendedAnisotropy2dGraph {
00204 
00205 public:
00206   typedef typename ConfiguratorType::RealType RealType;
00207   typedef typename ConfiguratorType::ElementType ElementType;
00208   
00209 private:
00210   AnisoType1 &_anisotropy1;       // for image values < _xMin
00211   AnisoType2 &_anisotropy2;       // for image values > _xMax
00212   const aol::DiscreteFunctionDefault<ConfiguratorType> *_discFunc;
00213   ConfiguratorType _config;
00214   qc::GridDefinition _grid;
00215   RealType _xMin;
00216   RealType _xMax;
00217   RealType _delta;
00218   RealType _deltaSqr;
00219   mutable qc::ScalarArray<RealType, qc::QC_2D> _testImg;
00220 
00221 public:
00222   LinearBlendedAnisotropy2dGraph ( AnisoType1 &Anisotropy1, AnisoType2 &Anisotropy2, 
00223                               const typename ConfiguratorType::InitType &Grid, RealType XMin, RealType XMax, RealType delta ) : 
00224       _anisotropy1 ( Anisotropy1 ), 
00225       _anisotropy2 ( Anisotropy2 ), 
00226       _discFunc( NULL ), 
00227       _config( Grid ), _grid( Grid ),
00228       _xMin( XMin ), _xMax( XMax ),
00229       _delta ( delta ), _deltaSqr ( delta*delta ),
00230       _testImg( Grid ) { _testImg.setAll( -1. ); 
00231   }
00232 
00233   
00234   void setImageReference ( const aol::Vector<RealType> &Image ) {
00235     if ( _discFunc )
00236       delete _discFunc;
00237     _discFunc = new aol::DiscreteFunctionDefault<ConfiguratorType>( this->_config.getInitializer(), Image );
00238   }
00239   
00240   void saveImage() const {
00241     _testImg.save( "BlendingTest.bz2", 9 );
00242   }
00243   
00244   //! this is the blending function.
00245   RealType eta( const RealType x ) const {
00246     if ( x < _xMin ) return 1.;
00247     if ( x > _xMax ) return 0.;
00248     return ( _xMax - x ) / ( _xMax - _xMin );
00249   }
00250   
00251   
00252   //! Evaluating the function itself, i.e. the blending
00253   //! \f$ \gamma(z) = \eta \gamma_1(z) + (1-\eta)\gamma_2(z) \f$.
00254   RealType gammaNorm ( const ElementType &El, const int QuadPoint, const aol::Vec2<RealType> &z ) const {
00255     if (!_discFunc) throw aol::Exception ( "\nSet Image first!\n", __FILE__, __LINE__ );
00256     aol::Vec3<RealType> z3 ( z[0], z[1], -1. );
00257     aol::Vec3<RealType> minusZ3( -z[0], -z[1], 1. );
00258     
00259     // get the weight of the blending function
00260     RealType x = _discFunc->evaluateAtQuadPoint( El, QuadPoint );
00261     RealType etaX = eta(x);
00262     
00263     // here is the linear blending
00264     if ( etaX == 0. ) return _anisotropy2.gammaNorm( z3 );
00265     else if ( etaX == 1. ) return _anisotropy1.gammaNorm( minusZ3 );
00266     else return ( 1. - etaX ) * _anisotropy2.gammaNorm( z3 ) + etaX * _anisotropy1.gammaNorm( minusZ3 );
00267   }
00268 
00269   
00270   //! evaluating the first derivative: \f$ ... \f$
00271   void gammaFirstDerivative ( const ElementType &El, const int QuadPoint,
00272                               const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
00273     aol::Vec3<RealType> z3 ( z[0], z[1], -1. );
00274     aol::Vec3<RealType> minusZ3( -z[0], -z[1], 1. );
00275     
00276     // get the weight of the blending function
00277     RealType x = _discFunc->evaluateAtQuadPoint( El, QuadPoint );
00278     RealType etaX = eta(x);
00279     
00280     aol::Vec3<RealType> resultAniso1, resultAniso2, v3;
00281     if ( etaX == 0. ) {
00282       _anisotropy2.gammaFirstDerivative( z3, resultAniso2 );
00283       resultAniso1.setZero();
00284     } else if ( etaX == 1. ) {
00285       _anisotropy1.gammaFirstDerivative( minusZ3, resultAniso1 );
00286       resultAniso2.setZero();
00287     } else {
00288       _anisotropy1.gammaFirstDerivative( minusZ3, resultAniso1 );
00289       _anisotropy2.gammaFirstDerivative( z3, resultAniso2 );
00290     }
00291     
00292     // here is the linear blending
00293     v3 = ( 1. - eta(x) ) * resultAniso2 - eta(x) * resultAniso1;
00294     
00295     v[0] = v3[0];     // copy the first two arguments of v3 (negative sign from the inner derivative)
00296     v[1] = v3[1];
00297   }
00298 
00299   
00300   //! this method computes \f$ gamma_zz \f$, it's called implicitPart just for technical reasons.
00301   void implicitPart ( const ElementType &El, const int QuadPoint, 
00302                       const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
00303     aol::Vec3<RealType> z3 ( z[0], z[1], -1. );
00304     aol::Vec3<RealType> minusZ3( -z[0], -z[1], 1. );
00305     
00306     // get the weight of the blending function
00307     RealType x = _discFunc->evaluateAtQuadPoint( El, QuadPoint );
00308     RealType etaX = eta(x);
00309     
00310     aol::Matrix33<RealType> resultAniso1, resultAniso2;
00311     if ( etaX == 0. ) {
00312       _anisotropy2.implicitPart( z3, resultAniso2 );
00313       resultAniso1.setZero();
00314     } else if ( etaX == 1. ) {
00315       _anisotropy1.implicitPart( minusZ3, resultAniso1 );
00316       resultAniso2.setZero();
00317     } else {
00318       _anisotropy1.implicitPart( minusZ3, resultAniso1 );
00319       _anisotropy2.implicitPart( z3, resultAniso2 );
00320     }
00321     
00322     resultAniso2 *= ( 1. - eta(x) );                        // weight the two anisotropies
00323     resultAniso1 *= eta(x);
00324     
00325     mat[0][0] = resultAniso1[0][0] + resultAniso2[0][0];    // linear blending
00326     mat[0][1] = resultAniso1[0][1] + resultAniso2[0][1]; 
00327     mat[1][0] = resultAniso1[1][0] + resultAniso2[1][0]; 
00328     mat[1][1] = resultAniso1[1][1] + resultAniso2[1][1]; 
00329   }
00330   
00331 };
00332 
00333 //! Class for blending two different anisotropies dependent on the image value, e.g. value>0.5 => use
00334 //! anisotropy1 else use anisotropy2 and use a linear blending in an interval between.
00335 template <typename ConfiguratorType, typename AnisoType1, typename AnisoType2>
00336 class LinearBlended3dAnisotropy {
00337 
00338 public:
00339   typedef typename ConfiguratorType::RealType RealType;
00340   typedef typename ConfiguratorType::ElementType ElementType;
00341   
00342 private:
00343   AnisoType1 &_anisotropy1;       // for image values < _xMin
00344   AnisoType2 &_anisotropy2;       // for image values > _xMax
00345   const aol::DiscreteFunctionDefault<ConfiguratorType> *_discFunc;    // either use an image for local blending values
00346   RealType _eta;                                                      // or one fixed value
00347   bool _blendingValueIsSet;                                           // flag whether user has defined at least one
00348   ConfiguratorType _config;
00349   qc::GridDefinition _grid;
00350   RealType _xMin;
00351   RealType _xMax;
00352   RealType _delta;
00353   RealType _deltaSqr;
00354   mutable qc::ScalarArray<RealType, qc::QC_3D> _testImg;
00355 
00356 public:
00357   LinearBlended3dAnisotropy ( AnisoType1 &Anisotropy1, AnisoType2 &Anisotropy2, 
00358                               const typename ConfiguratorType::InitType &Grid, RealType XMin, RealType XMax, RealType delta ) : 
00359       _anisotropy1 ( Anisotropy1 ), 
00360       _anisotropy2 ( Anisotropy2 ), 
00361       _discFunc( NULL ), 
00362       _config( Grid ), _grid( Grid ),
00363       _xMin( XMin ), _xMax( XMax ),
00364       _delta ( delta ), _deltaSqr ( delta*delta ),
00365       _testImg( Grid ) { _testImg.setAll( -1. ); 
00366   }
00367 
00368   
00369   void setImageReference ( const aol::Vector<RealType> &Image ) {
00370     if ( _discFunc )
00371       delete _discFunc;
00372     _discFunc = new aol::DiscreteFunctionDefault<ConfiguratorType>( this->_config.getInitializer(), Image );
00373     this->reset();
00374   }
00375     
00376   void setBlendingValue ( const RealType Eta ) {
00377     _eta = Eta;
00378     _blendingValueIsSet = true;
00379   }
00380   
00381   void saveImage() const {
00382     _testImg.save( "BlendingTest.bz2", 9 );
00383   }
00384   
00385   //! this is the blending function.
00386   RealType eta( const RealType x ) const {
00387     if ( x < _xMin ) return 1.;
00388     if ( x > _xMax ) return 0.;
00389     return ( _xMax - x ) / ( _xMax - _xMin );
00390   }
00391   
00392   
00393   //! Evaluating the function itself, i.e. the blending
00394   //! \f$ \gamma(z) = \eta \gamma_1(z) + (1-\eta)\gamma_2(z) \f$.
00395   RealType gammaNorm ( const ElementType &El, const int QuadPoint, const aol::Vec3<RealType> &z3 ) const {
00396     if (!_discFunc) throw aol::Exception ( "\nLinearBlended3dAnisotropy::Set Image first!\n", __FILE__, __LINE__ );
00397     
00398     // get the weight of the blending function
00399     RealType x = _discFunc->evaluateAtQuadPoint( El, QuadPoint );
00400     _eta = eta(x);
00401     
00402     // this flag is here, because otherwise it might happen, that a blending-image is defined, but
00403     // only operators that call the gammaNorm-function without Element and QuadPoint arguments
00404     // are used. Then _eta would be undefined.
00405     _blendingValueIsSet = true;
00406     
00407     // call the linear blending
00408     gammaNorm( z3 );
00409   }
00410 
00411   RealType gammaNorm ( const aol::Vec3<RealType> &z3 ) const {
00412     if ( !_blendingValueIsSet ) throw aol::Exception ( "\nLinearBlended3dAnisotropy: Set Blending-value first!\n", __FILE__, __LINE__ );
00413     aol::Vec3<RealType> minusZ3( -z3[0], -z3[1], -z3[2] );
00414     
00415     // here is the linear blending
00416     if ( _eta == 0. ) return _anisotropy2.gammaNorm( z3 );
00417     else if ( _eta == 1. ) return _anisotropy1.gammaNorm( minusZ3 );
00418     else return ( 1. - _eta ) * _anisotropy2.gammaNorm( z3 ) + _eta * _anisotropy1.gammaNorm( minusZ3 );
00419   }
00420   
00421   
00422   //! evaluating the first derivative: \f$ ... \f$
00423   void gammaFirstDerivative ( const ElementType &El, const int QuadPoint,
00424                               const aol::Vec3<RealType> &z3, aol::Vec3<RealType> &v ) const {
00425     // get the weight of the blending function
00426     RealType x = _discFunc->evaluateAtQuadPoint( El, QuadPoint );
00427     _eta = eta(x);
00428     
00429     _blendingValueIsSet = true;
00430     
00431     gammaFirstDerivative( z3, v );
00432   }
00433   
00434   //! evaluating the first derivative: \f$ ... \f$
00435   void gammaFirstDerivative ( const aol::Vec3<RealType> &z3, aol::Vec3<RealType> &v ) const {
00436     if ( !_blendingValueIsSet ) throw aol::Exception ( "\nLinearBlended3dAnisotropy: Set Blending-value first!\n", __FILE__, __LINE__ );
00437     aol::Vec3<RealType> minusZ3( -z3[0], -z3[1], -z3[2] );
00438     
00439     aol::Vec3<RealType> resultAniso1, resultAniso2;
00440     if ( _eta == 0. ) {
00441       _anisotropy2.gammaFirstDerivative( z3, resultAniso2 );
00442       resultAniso1.setZero();
00443     } else if ( _eta == 1. ) {
00444       _anisotropy1.gammaFirstDerivative( minusZ3, resultAniso1 );
00445       resultAniso2.setZero();
00446     } else {
00447       _anisotropy1.gammaFirstDerivative( minusZ3, resultAniso1 );
00448       _anisotropy2.gammaFirstDerivative( z3, resultAniso2 );
00449     }
00450     
00451     // here is the linear blending
00452     v = ( 1. - _eta ) * resultAniso2 - _eta * resultAniso1;
00453   }
00454    
00455 };
00456 
00457 
00458 //! Class for the simple L1-norm-anisotropy in 2d.
00459 template <typename RealType>
00460 class L1Norm2d {
00461 private:
00462   RealType _epsilon;
00463   RealType _epsSqr;
00464 
00465 public:
00466   L1Norm2d ( RealType epsilon ) : _epsilon ( epsilon ), _epsSqr ( epsilon*epsilon ) { }
00467 
00468 
00469   //! Evaluating the function itself, i.e. the regularization of
00470   //! \f$ \gamma(z) = \sqrt{z_0^2 + \epsilon^2} + \sqrt{z_1^2 + \epsilon^2} \f$.
00471   RealType gammaNorm ( const aol::Vec2<RealType> &z ) const {
00472     return ( sqrt ( z[0]*z[0] + _epsSqr ) + sqrt ( z[1]*z[1] + _epsSqr ) );
00473   }
00474 
00475   //! evaluating the first derivative: \f$ ... \f$
00476   void gammaFirstDerivative ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
00477     gammaFirstDerivative ( z, v );
00478   }
00479 
00480   //! evaluating the first derivative: \f$ ... \f$
00481   void gammaFirstDerivative ( const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
00482     RealType norm = gammaNorm ( z );
00483     v[0] = z[0] / norm;
00484     v[1] = z[1] / norm;
00485   }
00486 };
00487 
00488 
00489 //! Class for the simple parallelogram-anisotropy in 2d.
00490 template <typename RealType>
00491 class Parallelogram2d {
00492 private:
00493   RealType _epsilon;
00494   RealType _epsSqr;
00495   RealType _alpha;        // angle of x-axis with one side
00496   RealType _beta;         // angle of x-axis with the other side
00497   aol::Vec2<RealType> p;  // vector to one corner
00498   aol::Vec2<RealType> q;  // vector to another corner (not -p)
00499 
00500 public:
00501   Parallelogram2d ( RealType Epsilon, RealType Alpha, RealType Beta )
00502     : _epsilon ( Epsilon ),
00503       _epsSqr ( Epsilon*Epsilon ),
00504       _alpha ( Alpha ),
00505       _beta ( Beta )
00506   {
00507     p[0] = cos(_alpha)*cos(_beta) + sin(_alpha);    // the vector p
00508     p[1] = cos(_alpha)*sin(_beta);
00509     q[0] = cos(_alpha) - sin(_alpha)*cos(_beta);    // the vector q
00510     q[1] = - sin(_alpha)*sin(_beta);
00511   }
00512 
00513   //! Evaluating the function itself
00514   RealType gammaNorm ( const aol::Vec2<RealType> &z ) const {
00515      
00516     const RealType temp1 =  aol::Sqr( p*z );    // p*z
00517     const RealType temp2 =  aol::Sqr( q*z );    // q*z
00518 
00519     return sqrt(temp1+_epsSqr) + sqrt(temp2+_epsSqr);
00520   }
00521 
00522 };
00523 
00524 
00525 //! Class for the simple L1-norm-anisotropy in 2d.
00526 template <typename RealType>
00527 class LInfNorm2d {
00528 private:
00529   RealType _epsilon;
00530   RealType _epsSqr;
00531 
00532 public:
00533   LInfNorm2d ( RealType epsilon ) : _epsilon ( epsilon ), _epsSqr ( epsilon*epsilon ) { }
00534 
00535 
00536   //! Evaluating the function itself, i.e. the regularization of
00537   //! \f$ \gamma(z) = \sqrt{z_0^2 + \epsilon^2} + \sqrt{z_1^2 + \epsilon^2} \f$.
00538   RealType gammaNorm ( const aol::Vec2<RealType> &z ) const {
00539     return ( 0.5 * ( sqrt ( aol::Sqr( z[0] - z[1] ) + _epsSqr ) + sqrt ( aol::Sqr( z[0] + z[1] ) + _epsSqr ) ) );
00540   }
00541 
00542   //! evaluating the first derivative: \f$ ... \f$
00543   void gammaFirstDerivative ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
00544     gammaFirstDerivative ( z, v );
00545   }
00546 
00547   //! evaluating the first derivative: \f$ ... \f$
00548   void gammaFirstDerivative ( const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
00549     v[0] = 0.5 * ( aol::signum1at0( z[0] - z[1] ) + aol::signum1at0( z[0] + z[1] ) );
00550     v[1] = 0.5 * ( aol::signum1at0( z[0] + z[1] ) - aol::signum1at0( z[0] - z[1] ) );
00551   }
00552 };
00553 
00554 //! Class for the simple LInf-norm-anisotropy in 3d: 
00555 //! \f$ \gamma(z) = \max \{ |z_1|,|z_2|,|z_3| \} \f$. While the norm itself is simple, 
00556 //! the regularized norm is really awful and is therefore not given in TeX here!
00557 template <typename RealType>
00558 class LInfNorm3d : public Anisotropy3dGraphInterface<RealType, LInfNorm3d<RealType> > {
00559 public:  
00560   using Anisotropy3dGraphInterface<RealType, LInfNorm3d<RealType> >::gammaNorm;
00561   using Anisotropy3dGraphInterface<RealType, LInfNorm3d<RealType> >::gammaFirstDerivative;
00562   using Anisotropy3dGraphInterface<RealType, LInfNorm3d<RealType> >::implicitPart;
00563   
00564 private:
00565   RealType _delta;
00566   RealType _deltaSqr;
00567 
00568 public:
00569   LInfNorm3d ( RealType delta )  :
00570     Anisotropy3dGraphInterface<RealType, LInfNorm3d<RealType> > ( ), _delta ( delta ), _deltaSqr ( delta*delta ) { }
00571 
00572   //! Evaluating the function itself, i.e. the regularization of
00573   //! \f$ \gamma(z) = \max \{ |z_1|,|z_2|,|z_3| \} \f$. Code is generated by Maple and 
00574   //! thus not very readable, but should be correct (otherwise blame Mapele ;-)
00575   RealType gammaNorm ( const aol::Vec3<RealType> &z ) const {
00576       RealType t1 = z[0]*z[0];            RealType t3 = 2.0*z[0]*z[1];            
00577       RealType t4 = z[1]*z[1];            RealType t5 = _deltaSqr;         
00578       RealType t7 = sqrt(t1+t3+t4+t5);    RealType t9 = sqrt(t1-t3+t4+t5);
00579       RealType t10 = 2.0*z[0];            RealType t12 = aol::Sqr( t7+t9+t10 );   
00580       RealType t14 = sqrt(t12+t5);        RealType t17 = aol::Sqr( t7+t9-t10 );
00581       RealType t19 = sqrt(t17+t5);
00582       
00583       return ( 0.25*t14+0.25*t19 );
00584   }
00585 
00586   //! evaluating the first derivative: \f$ ... \f$
00587   void gammaFirstDerivative ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
00588     gammaFirstDerivative ( z, v );
00589   }
00590 
00591   //! evaluating the first derivative: \f$ ... \f$
00592   void gammaFirstDerivative ( const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
00593       RealType t1 = z[0]*z[0];             RealType t3 = 2.0*z[0]*z[1];           
00594       RealType t4 = z[1]*z[1];             RealType t5 = _deltaSqr;             
00595       RealType t7 = sqrt(t1+t3+t4+t5);     RealType t9 = sqrt(t1-t3+t4+t5);
00596       RealType t10 = 2.0*z[2];             RealType t11 = t7+t9+t10;              
00597       RealType t12 = t11*t11;              RealType t14 = sqrt(t12+t5);         
00598       RealType t15 = 1/t14;                RealType t16 = t15*t11;
00599       RealType t19 = 2.0/t7*(z[0]+z[1]);   RealType t20 = 1/t9;                   
00600       RealType t21 = z[0]-z[1];            RealType t23 = t19+2.0*t20*t21;      
00601       RealType t26 = t7+t9-t10;            RealType t27 = t26*t26;
00602       RealType t29 = sqrt(t27+t5);         RealType t30 = 1/t29;                  
00603       RealType t31 = t30*t26;              RealType t36 = t19-2.0*t20*t21;      
00604       RealType t42 = 4.0*t7;               RealType t43 = 4.0*t9;
00605       RealType t44 = 8.0*z[2];             
00606       
00607       v[0] = 0.125*t16*t23+0.125*t31*t23;
00608       v[1] = 0.125*t16*t36+0.125*t31*t36;
00609       v[2] = 0.125*t15*(t42+t43+t44)+0.125*t30*(-t42-t43+t44);
00610   }
00611   
00612     
00613   //! this method computes \f$ gamma_zz \f$, it's called implicitPart just for technical reasons.
00614   void implicitPart ( const aol::Vec3<RealType> &z, aol::Matrix33<RealType> &mat ) const {
00615       RealType t1 = z[0]*z[0];             RealType t3 = 2.0*z[0]*z[1];               
00616       RealType t4 = z[1]*z[1];             RealType t5 = _deltaSqr;              
00617       RealType t6 = t1+t3+t4+t5;           RealType t7 = sqrt(t6);
00618       RealType t8 = t1-t3+t4+t5;           RealType t9 = sqrt(t8);                    
00619       RealType t10 = 2.0*z[2];             RealType t11 = t7+t9+t10;            
00620       RealType t12 = t11*t11;              RealType t13 = t12+t5;
00621       RealType t14 = sqrt(t13);            RealType t16 = 1/t14/t13;                  
00622       RealType t17 = t16*t12;              RealType t18 = 1/t7;                 
00623       RealType t19 = z[0]+z[1];            RealType t20 = 2.0*t18*t19;
00624       RealType t21 = 1/t9;                 RealType t22 = z[0]-z[1];                  
00625       RealType t24 = t20+2.0*t21*t22;      RealType t25 = t24*t24/4.0;          
00626       RealType t28 = 1/t14;                RealType t31 = t28*t11;
00627       RealType t36 = 1/t7/t6*t19*t19;      RealType t38 = 1/t9/t8;                    
00628       RealType t42 = -t36+t18-t38*t22*t22+t21;
00629       RealType t44 = 0.25*t31*t42;         RealType t45 = t7+t9-t10;                  
00630       RealType t46 = t45*t45;              RealType t47 = t46+t5;               
00631       RealType t48 = sqrt(t47);            RealType t50 = 1/t48/t47;
00632       RealType t51 = t50*t46;              RealType t54 = 1/t48;                      
00633       RealType t57 = t54*t45;              RealType t59 = 0.25*t57*t42;         
00634       RealType t62 = t20-2.0*t21*t22;      RealType t63 = t24*t62/4.0;
00635       RealType t66 = t28*t62/2.0;          RealType t72 = -t36+t18+t38*t22*t22-t21;   
00636       RealType t77 = t54*t62/2.0;
00637       RealType t82 = -0.25*t17*t63+0.125*t66*t24+0.25*t31*t72-0.25*t51*t63+0.125*t77*t24+0.25*t57*t72;
00638       RealType t83 = t16*t11;              RealType t84 = 4.0*t7;                 
00639       RealType t85 = 4.0*t9;               RealType t86 = 8.0*z[2];             
00640       RealType t87 = t84+t85+t86;          RealType t90 = 0.625E-1*t83*t24*t87;
00641       RealType t93 = t50*t45;              RealType t94 = -t84-t85+t86;           
00642       RealType t97 = 0.625E-1*t93*t24*t94; RealType t101 = t62*t62/4.0;         
00643       RealType t113 = 0.625E-1*t83*t62*t87;  
00644       RealType t117 = 0.625E-1*t93*t62*t94;
00645       RealType t130 = t87*t87;             RealType t134 = t94*t94;
00646       
00647       mat[0][0] = -0.25*t17*t25+0.25*t28*t25+t44-0.25*t51*t25+0.25*t54*t25+t59;
00648       mat[0][1] = t82;
00649       mat[0][2] = -t90+0.25*t28*t24-t97-0.25*t54*t24;
00650       mat[1][0] = t82;
00651       mat[1][1] = -0.25*t17*t101+0.25*t28*t101+t44-0.25*t51*t101+0.25*t54*t101+t59;
00652       mat[1][2] = -t113+0.5*t66-t117-0.5*t77;
00653       mat[2][0] = -t90+0.25*t28*t24-t97-0.25*t54*t24;
00654       mat[2][1] = -t113+0.25*t28*t62-t117-0.25*t54*t62;
00655       mat[2][2] = -0.625E-1*t16*t130+0.1E1*t28-0.625E-1*t50*t134+0.1E1*t54;
00656     
00657 #ifdef DEBUG    
00658     for (int i=0; i<3; i++)                                          
00659       for (int j=0; j<3; j++)                                        
00660         if ( aol::isNaN( mat[i][j] ) ) cerr << "\n******* InfNorm3D: Found NAN in HessMat!! ******* ";
00661 #endif
00662                                                                          
00663   } 
00664   
00665 };
00666 
00667 
00668 
00669 //! Class for the simple L1-norm-anisotropy in 3d (Wulff shape is a cube).
00670 template <typename RealType>
00671 class L1Norm3d : public Anisotropy3dGraphInterface<RealType, L1Norm3d<RealType> > {
00672 public:  
00673   using Anisotropy3dGraphInterface<RealType, L1Norm3d<RealType> >::gammaNorm;
00674   using Anisotropy3dGraphInterface<RealType, L1Norm3d<RealType> >::gammaFirstDerivative;
00675   using Anisotropy3dGraphInterface<RealType, L1Norm3d<RealType> >::implicitPart;
00676   
00677 private:
00678   RealType _epsilon;
00679   RealType _epsSqr;
00680 
00681 public:
00682   L1Norm3d ( RealType epsilon ) :
00683       Anisotropy3dGraphInterface<RealType, L1Norm3d<RealType> > ( ), _epsilon ( epsilon ), _epsSqr ( epsilon*epsilon ) { }
00684 
00685   //! Evaluating the function itself, i.e. the regularization of
00686   //! \f$ \gamma(z) = \sqrt{z_0^2 + \epsilon^2} + \sqrt{z_1^2 + \epsilon^2} + \sqrt{z_2^2 + \epsilon^2} \f$.
00687   RealType gammaNorm ( const aol::Vec3<RealType> &z ) const {
00688     return ( sqrt ( z[0]*z[0] + _epsSqr ) + sqrt ( z[1]*z[1] + _epsSqr ) + sqrt ( z[2]*z[2] + _epsSqr ) );
00689   }
00690 
00691   //! evaluating the first derivative: \f$ \gamma_z = \frac{z}{|z|_{\varepsilon}} \f$
00692   void gammaFirstDerivative ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
00693     gammaFirstDerivative ( z, v );
00694   }
00695 
00696   //! evaluating the first derivative: \f$ \gamma_z = \frac{z}{|z|_{\varepsilon}} \f$
00697   void gammaFirstDerivative ( const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
00698     v[0] = z[0] / sqrt ( z[0]*z[0] + _epsSqr );
00699     v[1] = z[1] / sqrt ( z[1]*z[1] + _epsSqr );
00700     v[2] = z[2] / sqrt ( z[2]*z[2] + _epsSqr );
00701   }
00702   
00703   //! this method computes \f$ gamma_zz \f$, it's called implicitPart just for technical reasons.
00704   void implicitPart ( const aol::Vec3<RealType> &z, aol::Matrix33<RealType> &mat ) const {
00705     mat[0][0] = _epsSqr / pow( z[0]*z[0] + _epsSqr, 1.5 );
00706     mat[1][1] = _epsSqr / pow( z[1]*z[1] + _epsSqr, 1.5 );
00707     mat[2][2] = _epsSqr / pow( z[2]*z[2] + _epsSqr, 1.5 );
00708     mat[0][1] = mat[0][2] = mat[1][0] = mat[1][2] = mat[2][0] = mat[2][1] = 0.;
00709     
00710 #ifdef DEBUG    
00711     for (int i=0; i<3; i++)                                          
00712       for (int j=0; j<3; j++)                                        
00713         if ( aol::isNaN( mat[i][j] ) ) cerr << "\n******* Found NAN in HessMat!! ******* ";
00714 #endif
00715                                                                          
00716   }           
00717   
00718 };
00719 
00720 //! Class for the double cone anisotropy in 3d (Frank diagram is a cylinder aligned to the z-axis).
00721 template <typename RealType>
00722 class DoubleCone3d : public Anisotropy3dGraphInterface<RealType, DoubleCone3d<RealType> > {
00723 public:  
00724   using Anisotropy3dGraphInterface<RealType, DoubleCone3d<RealType> >::gammaNorm;
00725   using Anisotropy3dGraphInterface<RealType, DoubleCone3d<RealType> >::gammaFirstDerivative;
00726   using Anisotropy3dGraphInterface<RealType, DoubleCone3d<RealType> >::implicitPart;
00727 
00728 private:
00729   RealType _delta;      //! Regularization parameter
00730   RealType _deltaSqr;
00731 
00732 public:
00733   DoubleCone3d ( RealType delta ) :
00734     Anisotropy3dGraphInterface<RealType, DoubleCone3d<RealType> > ( ), _delta ( delta ), _deltaSqr ( delta*delta ) { }
00735 
00736   //! Evaluating the function itself, i.e. the regularization of
00737   //! \f$ \gamma(z) = \sqrt{z_0^2 + \epsilon^2} + \sqrt{z_1^2 + \epsilon^2} + \sqrt{z_2^2 + \epsilon^2} \f$.
00738   RealType gammaNorm ( const aol::Vec3<RealType> &z ) const {
00739     return ( 0.5 * ( sqrt( aol::Sqr( sqrt ( z[0]*z[0] + z[1]*z[1] ) + z[2] ) + _deltaSqr ) 
00740                    + sqrt( aol::Sqr( sqrt ( z[0]*z[0] + z[1]*z[1] ) - z[2] ) + _deltaSqr ) ) );
00741   }
00742 
00743   //! evaluating the first derivative: \f$ \gamma_z = \frac{z}{|z|_{\varepsilon}} \f$
00744   void gammaFirstDerivative ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
00745     gammaFirstDerivative ( z, v );
00746   }
00747 
00748   //! evaluating the first derivative: \f$ \gamma_z = \frac{z}{|z|_{\varepsilon}} \f$
00749   void gammaFirstDerivative ( const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
00750     const RealType t1 = z[0]*z[0];
00751     const RealType t2 = z[1]*z[1];
00752     const RealType t3 = _delta*_delta;
00753     const RealType t5 = sqrt(t1+t2+t3);
00754     const RealType t6 = t5+z[2];
00755     const RealType t7 = t6*t6;
00756     const RealType t9 = sqrt(t7+t3);
00757     const RealType t10 = 1/t9;
00758     const RealType t11 = t10*t6;
00759     const RealType t12 = 1/t5;
00760     const RealType t13 = t12*z[0];
00761     const RealType t16 = t5-z[2];
00762     const RealType t17 = t16*t16;
00763     const RealType t19 = sqrt(t17+t3);
00764     const RealType t20 = 1/t19;
00765     const RealType t21 = t20*t16;
00766     const RealType t25 = t12*z[1];
00767     v[0] = 0.5*t11*t13+0.5*t21*t13;
00768     v[1] = 0.5*t11*t25+0.5*t21*t25;
00769     v[2] = 0.5*t10*t6-0.5*t20*t16;
00770   }
00771   
00772   //! this method computes \f$ gamma_zz \f$, it's called implicitPart just for technical reasons.
00773   void implicitPart ( const aol::Vec3<RealType> &z, aol::Matrix33<RealType> &mat ) const {
00774     // maple part
00775     RealType x = z[0]; RealType y = z[1];
00776     
00777     const RealType t1 = x*x;
00778     const RealType t2 = y*y;
00779     const RealType t3 = _delta*_delta;
00780     const RealType t4 = t1+t2+t3;
00781     const RealType t5 = sqrt(t4);
00782     const RealType t6 = t5+z[2];
00783     const RealType t7 = t6*t6;
00784     const RealType t8 = t7+t3;
00785     const RealType t9 = sqrt(t8);
00786     const RealType t11 = 1/t9/t8;
00787     const RealType t12 = t11*t7;
00788     const RealType t13 = 1/t4;
00789     const RealType t14 = t13*t1;
00790     const RealType t17 = 1/t9;
00791     const RealType t18 = t17*t13;
00792     const RealType t21 = t17*t6;
00793     const RealType t23 = 1/t5/t4;
00794     const RealType t24 = t23*t1;
00795     const RealType t27 = 1/t5;
00796     const RealType t29 = 0.5*t21*t27;
00797     const RealType t30 = t5-z[2];
00798     const RealType t31 = t30*t30;
00799     const RealType t32 = t31+t3;
00800     const RealType t33 = sqrt(t32);
00801     const RealType t35 = 1/t33/t32;
00802     const RealType t36 = t35*t31;
00803     const RealType t39 = 1/t33;
00804     const RealType t40 = t39*t13;
00805     const RealType t43 = t39*t30;
00806     const RealType t47 = 0.5*t43*t27;
00807     const RealType t50 = t13*x*y;
00808     const RealType t53 = x*y;
00809     const RealType t57 = t23*x*y;
00810     const RealType t66 = -0.5*t12*t50+0.5*t18*t53-0.5*t21*t57-0.5*t36*t50+0.5*t40*t53-0.5*t43*t57;
00811     const RealType t68 = t27*x;
00812     const RealType t71 = 0.5*t11*t6*t6*t68;
00813     const RealType t72 = t17*t27;
00814     const RealType t73 = t72*x;
00815     const RealType t78 = -0.5*t35*t30*t30*t68;
00816     const RealType t79 = t39*t27;
00817     const RealType t80 = t79*x;
00818     const RealType t83 = t13*t2;
00819     const RealType t88 = t23*t2;
00820     const RealType t102 = 0.5*t11*t6*t6*t27*y;
00821     const RealType t103 = t72*y;
00822     const RealType t109 = -0.5*t35*t30*t30*t27*y;
00823     const RealType t110 = t79*y;
00824     mat[0][0] = -0.5*t12*t14+0.5*t18*t1-0.5*t21*t24+t29-0.5*t36*t14+0.5*t40*t1-0.5*t43*t24+t47;
00825     mat[0][1] = t66;
00826     mat[0][2] = -t71+0.5*t73-t78-0.5*t80;
00827     mat[1][0] = t66;
00828     mat[1][1] = -0.5*t12*t83+0.5*t18*t2-0.5*t21*t88+t29-0.5*t36*t83+0.5*t40*t2-0.5*t43*t88+t47;
00829     mat[1][2] = -t102+0.5*t103-t109-0.5*t110;
00830     mat[2][0] = -t71+0.5*t73-t78-0.5*t80;
00831     mat[2][1] = -t102+0.5*t103-t109-0.5*t110;
00832     mat[2][2] = -0.5*t11*t6*t6+0.5*t17-0.5*t35*t30*t30+0.5*t39;
00833   }
00834   
00835   
00836 };
00837 
00838 //! Class for a hexagon in 2d.
00839 template <typename RealType>
00840 class Hexagon2d {
00841 
00842 private:
00843   RealType _delta;          //! Regularization parameter
00844   RealType _deltaSqr;
00845   RealType _alpha;          //! the angle of one hexagon-site and the x-axis (or y-axis)
00846   aol::Vec2<RealType> _Vm;  //! V1-V2   
00847   aol::Vec2<RealType> _Vp;  //! V1+V2
00848   aol::Vec2<RealType> _V3;  //! +-V1,+-V2,+-V3 are the corners of the hexagon in the r-z-plane (r=sqrt(x^2+y^2))
00849   
00850 
00851 public:
00852   Hexagon2d ( RealType delta, RealType alpha ) : _delta ( delta ), _deltaSqr ( delta*delta ), _alpha( alpha )
00853   {
00854 //     _Vm[0] = 0.5;                         // regelmaessiges Hexagon, alpha=60 degrees
00855 //     _Vm[1] = -0.8660254037844386467326525391730029923565;
00856 //     _Vp[0] = 1.5;
00857 //     _Vp[1] = 0.8660254037844386467326525391730029923565;
00858 //     _V3[0] = -0.5;
00859 //     _V3[1] = 0.8660254037844386467326525391730029923565;
00860     
00861 //     _Vm[0] = cos( _alpha );        // gleichseitiges Hexagon
00862 //     _Vm[1] = -sin( _alpha );
00863 //     _Vp[0] = 1. + cos( _alpha );
00864 //     _Vp[1] = sin( _alpha );
00865 //     _V3[0] = -0.5;
00866 //     _V3[1] = sin( _alpha );
00867     
00868     _Vm[0] = cos( _alpha ) / sin ( _alpha );  // Die Normalen auf die oberen Seitenflaechen nehmen den gleichen Wert an
00869     _Vm[1] = -1.;
00870     _Vp[0] = ( 2. - cos( _alpha ) ) / sin( _alpha );
00871     _Vp[1] = 1.;
00872     _V3[0] = (cos( _alpha ) - 1.) / sin( _alpha );
00873     _V3[1] = 1.;
00874   }
00875 
00876   //! Evaluating the function itself, i.e. the regularization of
00877   //! \f$ \gamma(z) = \sqrt{z_0^2 + \epsilon^2} + \sqrt{z_1^2 + \epsilon^2} + \sqrt{z_2^2 + \epsilon^2} \f$.
00878   RealType gammaNorm ( const aol::Vec2<RealType> &z ) const {
00879     const RealType t1 = _Vm[0]*_Vm[0];
00880     const RealType t2 = z[0]*z[0];
00881     const RealType t8 = _Vm[1]*_Vm[1];
00882     const RealType t9 = z[1]*z[1];
00883     const RealType t11 = _delta*_delta;
00884     const RealType t13 = sqrt(t1*t2+2.0*_Vm[0]*z[0]*_Vm[1]*z[1]+t8*t9+t11);
00885     const RealType t14 = 0.5*t13;
00886     const RealType t15 = _Vp[0]*_Vp[0];
00887     const RealType t21 = _Vp[1]*_Vp[1];
00888     const RealType t24 = sqrt(t15*t2+2.0*_Vp[0]*z[0]*_Vp[1]*z[1]+t21*t9+t11);
00889     const RealType t25 = 0.5*t24;
00890     const RealType t26 = _V3[0]*z[0];
00891     const RealType t27 = _V3[1]*z[1];
00892     const RealType t29 = pow(t14+t25-t26-t27,2.0);
00893     const RealType t31 = sqrt(t29+t11);
00894     const RealType t34 = pow(t14+t25+t26+t27,2.0);
00895     const RealType t36 = sqrt(t34+t11);
00896      
00897     return( 0.5*t31+0.5*t36 );
00898   }
00899 
00900   //! evaluating the first derivative: \f$ \gamma_z = \frac{z}{|z|_{\varepsilon}} \f$
00901   void gammaFirstDerivative ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
00902     gammaFirstDerivative ( z, v );
00903   }
00904 
00905   //! evaluating the first derivative: \f$ \gamma_z = \frac{z}{|z|_{\varepsilon}} \f$
00906   void gammaFirstDerivative ( const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
00907     const RealType t1 = _Vm[0]*_Vm[0];
00908     const RealType t2 = z[0]*z[0];
00909     const RealType t4 = _Vm[0]*z[0];
00910     const RealType t8 = _Vm[1]*_Vm[1];
00911     const RealType t9 = z[1]*z[1];
00912     const RealType t11 = _delta*_delta;
00913     const RealType t13 = sqrt(t1*t2+2.0*t4*_Vm[1]*z[1]+t8*t9+t11);
00914     const RealType t14 = 0.5*t13;
00915     const RealType t15 = _Vp[0]*_Vp[0];
00916     const RealType t17 = _Vp[0]*z[0];
00917     const RealType t21 = _Vp[1]*_Vp[1];
00918     const RealType t24 = sqrt(t15*t2+2.0*t17*_Vp[1]*z[1]+t21*t9+t11);
00919     const RealType t25 = 0.5*t24;
00920     const RealType t26 = _V3[0]*z[0];
00921     const RealType t27 = _V3[1]*z[1];
00922     const RealType t28 = t14+t25-t26-t27;
00923     const RealType t29 = t28*t28;
00924     const RealType t31 = sqrt(t29+t11);
00925     const RealType t33 = 1/t31*t28;
00926     const RealType t34 = 1/t13;
00927     const RealType t40 = 0.5*t34*(t1*z[0]+_Vm[0]*_Vm[1]*z[1]);
00928     const RealType t41 = 1/t24;
00929     const RealType t47 = 0.5*t41*(t15*z[0]+_Vp[0]*_Vp[1]*z[1]);
00930     const RealType t51 = t14+t25+t26+t27;
00931     const RealType t52 = t51*t51;
00932     const RealType t54 = sqrt(t52+t11);
00933     const RealType t56 = 1/t54*t51;
00934     const RealType t65 = 0.5*t34*(t4*_Vm[1]+t8*z[1]);
00935     const RealType t70 = 0.5*t41*(t17*_Vp[1]+t21*z[1]);
00936     
00937     v[0] = 0.5*t33*(t40+t47-_V3[0])+0.5*t56*(t40+t47+_V3[0]);
00938     v[1] = 0.5*t33*(t65+t70-_V3[1])+0.5*t56*(t65+t70+_V3[1]);
00939   }
00940  
00941 };
00942 
00943 
00944 //! Class for an rotated hexagon in 3d.
00945 template <typename RealType>
00946 class RotatedHexagon3d : public Anisotropy3dGraphInterface<RealType, RotatedHexagon3d<RealType> > {
00947 
00948 public:  
00949   using Anisotropy3dGraphInterface<RealType, RotatedHexagon3d<RealType> >::gammaNorm;
00950   using Anisotropy3dGraphInterface<RealType, RotatedHexagon3d<RealType> >::gammaFirstDerivative;
00951   using Anisotropy3dGraphInterface<RealType, RotatedHexagon3d<RealType> >::implicitPart;
00952   
00953 private:
00954   RealType _delta;          //! Regularization parameter
00955   RealType _deltaSqr;
00956   RealType _alpha;          //! the angle of one hexagon-site and the x-axis (or y-axis)
00957   aol::Vec2<RealType> _Vm;  //! V1-V2   
00958   aol::Vec2<RealType> _Vp;  //! V1+V2
00959   aol::Vec2<RealType> _V3;  //! +-V1,+-V2,+-V3 are the corners of the hexagon in the r-z-plane (r=sqrt(x^2+y^2))
00960   
00961 
00962 public:
00963   RotatedHexagon3d ( RealType delta, RealType alpha ) :
00964     Anisotropy3dGraphInterface<RealType, RotatedHexagon3d<RealType> > ( ), _delta ( delta ), _deltaSqr ( delta*delta ), _alpha( alpha )
00965   {
00966 /*    _Vm[0] = cos( _alpha );        // gleichseitiges Hexagon
00967     _Vm[1] = -sin( _alpha );
00968     _Vp[0] = 1. + cos( _alpha );
00969     _Vp[1] = sin( _alpha );
00970     _V3[0] = -0.5;
00971     _V3[1] = sin( _alpha ); */
00972     
00973     _Vm[0] = cos( _alpha ) / sin ( _alpha );  // Die Normalen auf die oberen Seitenflaechen nehmen den gleichen Wert an
00974     _Vm[1] = -1.;
00975     _Vp[0] = ( 2. - cos( _alpha ) ) / sin( _alpha );
00976     _Vp[1] = 1.;
00977     _V3[0] = (cos( _alpha ) - 1.) / sin( _alpha );
00978     _V3[1] = 1.;
00979    
00980 //     _Vm[0] = 1.;        // zeichnerisch vorgegebene Punkte, so dass die Parallele zur r-Achse erheblich billiger
00981 //     _Vm[1] = -1.;       // ist als die Seitenwaende. Die Seitenwaende haben zur r-Achse einen Winkel von 45 Grad.
00982 //     _Vp[0] = 5.;
00983 //     _Vp[1] = 1.;
00984 //     _V3[0] = -2.;
00985 //     _V3[1] = 1.;
00986     
00987 //     _Vm[0] = 2.5;        // zeichnerisch vorgegebene Punkte, so dass die Parallele zur r-Achse erheblich teurer
00988 //     _Vm[1] = -2.5;       // ist als die Seitenwaende. Die Seitenwaende haben zur r-Achse einen Winkel von 45 Grad.
00989 //     _Vp[0] = 3.5;
00990 //     _Vp[1] = 2.5;
00991 //     _V3[0] = -0.5;
00992 //     _V3[1] = 2.5;
00993     
00994   }
00995 
00996   //! Evaluating the function itself, i.e. the regularization of
00997   //! \f$ \gamma(z) = \sqrt{z_0^2 + \epsilon^2} + \sqrt{z_1^2 + \epsilon^2} + \sqrt{z_2^2 + \epsilon^2} \f$.
00998   RealType gammaNorm ( const aol::Vec3<RealType> &z ) const {
00999     const RealType t1 = z[0]*z[0];   
01000     const RealType t2 = z[1]*z[1];   
01001     const RealType t3 = _delta*_delta;   
01002     const RealType t5 = sqrt(t1+t2+t3);   
01003     const RealType t9 = aol::Sqr(_Vm[0]*t5+_Vm[1]*z[2]);   
01004     const RealType t11 = sqrt(t9+t3);   
01005     const RealType t12 = 0.5*t11;   
01006     const RealType t16 = aol::Sqr(_Vp[0]*t5+_Vp[1]*z[2]);   
01007     const RealType t18 = sqrt(t16+t3);   
01008     const RealType t19 = 0.5*t18;   
01009     const RealType t20 = _V3[0]*t5;   
01010     const RealType t21 = _V3[1]*z[2];   
01011     const RealType t23 = aol::Sqr(t12+t19-t20-t21);   
01012     const RealType t25 = sqrt(t23+t3);   
01013     const RealType t28 = aol::Sqr(t12+t19+t20+t21);   
01014     const RealType t30 = sqrt(t28+t3);   
01015     return( 0.5*t25+0.5*t30 );   
01016   }                                                                  
01017                                                                      
01018   //! evaluating the first derivative: \f$ \gamma_z = \frac{z}{|z|_{\varepsilon}} \f$
01019   void gammaFirstDerivative ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
01020     gammaFirstDerivative ( z, v );                                   
01021   }                                                                  
01022                                                                      
01023   //! evaluating the first derivative: \f$ \gamma_z = \frac{z}{|z|_{\varepsilon}} \f$
01024   void gammaFirstDerivative ( const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
01025     const RealType t1 = z[0]*z[0];          
01026     const RealType t2 = z[1]*z[1];          
01027     const RealType t3 = _delta*_delta;        
01028     const RealType t5 = sqrt(t1+t2+t3);        
01029     const RealType t8 = _Vm[0]*t5+_Vm[1]*z[2];        
01030     const RealType t9 = t8*t8;        
01031     const RealType t11 = sqrt(t9+t3);        
01032     const RealType t12 = 0.5*t11;        
01033     const RealType t15 = _Vp[0]*t5+_Vp[1]*z[2];        
01034     const RealType t16 = t15*t15;        
01035     const RealType t18 = sqrt(t16+t3);        
01036     const RealType t19 = 0.5*t18;        
01037     const RealType t20 = _V3[0]*t5;        
01038     const RealType t21 = _V3[1]*z[2];        
01039     const RealType t22 = t12+t19-t20-t21;        
01040     const RealType t23 = t22*t22;        
01041     const RealType t25 = sqrt(t23+t3);        
01042     const RealType t27 = 1/t25*t22;        
01043     const RealType t29 = 1/t11*t8;        
01044     const RealType t30 = 1/t5;        
01045     const RealType t31 = _Vm[0]*t30;        
01046     const RealType t34 = 0.5*t29*t31*z[0];        
01047     const RealType t36 = 1/t18*t15;        
01048     const RealType t37 = _Vp[0]*t30;        
01049     const RealType t40 = 0.5*t36*t37*z[0];        
01050     const RealType t41 = _V3[0]*t30;        
01051     const RealType t42 = t41*z[0];        
01052     const RealType t46 = t12+t19+t20+t21;        
01053     const RealType t47 = t46*t46;        
01054     const RealType t49 = sqrt(t47+t3);        
01055     const RealType t51 = 1/t49*t46;        
01056     const RealType t58 = 0.5*t29*t31*z[1];        
01057     const RealType t61 = 0.5*t36*t37*z[1];        
01058     const RealType t62 = t41*z[1];        
01059     const RealType t71 = 0.5*t29*_Vm[1];        
01060     const RealType t73 = 0.5*t36*_Vp[1];        
01061     v[0] = 0.5*t27*(t34+t40-t42)+0.5*t51*(t34+t40+t42);        
01062     v[1] = 0.5*t27*(t58+t61-t62)+0.5*t51*(t58+t61+t62);        
01063     v[2] = 0.5*t27*(t71+t73-_V3[1])+0.5*t51*(t71+t73+_V3[1]);      
01064                                                                      
01065 //       for (int j=0; j<3; j++)      
01066 //         cerr << v[j] << ", ";                                  
01067 //         if ( aol::isNaN( v[j] ) ) cerr << "\n******* Found NAN in Grad!! ******* ";
01068   }                                                                  
01069                                                                      
01070    //! this method computes \f$ gamma_zz \f$, it's called implicitPart just for technical reasons.
01071   void implicitPart ( const aol::Vec3<RealType> &z, aol::Matrix33<RealType> &mat ) const {
01072                                                                      
01073     const RealType t1 = z[0]*z[0];              const RealType t2 = z[1]*z[1];                const RealType t3 = _delta*_delta;
01074     const RealType t4 = t1+t2+t3;               const RealType t5 = sqrt(t4);                 const RealType t8 = _Vm[0]*t5+_Vm[1]*z[2];
01075     const RealType t9 = t8*t8;                  const RealType t10 = t9+t3;                   const RealType t11 = sqrt(t10);
01076     const RealType t12 = 0.5*t11;               const RealType t15 = _Vp[0]*t5+_Vp[1]*z[2];   const RealType t16 = t15*t15;
01077     const RealType t17 = t16+t3;                const RealType t18 = sqrt(t17);               const RealType t19 = 0.5*t18;
01078     const RealType t20 = _V3[0]*t5;             const RealType t21 = _V3[1]*z[2];             const RealType t22 = t12+t19-t20-t21;
01079     const RealType t23 = t22*t22;               const RealType t24 = t23+t3;                  const RealType t25 = sqrt(t24);
01080     const RealType t28 = 1/t25/t24*t23;         const RealType t29 = 1/t11;                   const RealType t30 = t29*t8;
01081     const RealType t31 = 1/t5;                  const RealType t32 = _Vm[0]*t31;              const RealType t33 = t32*z[0];
01082     const RealType t35 = 0.5*t30*t33;           const RealType t36 = 1/t18;                   const RealType t37 = t36*t15;
01083     const RealType t38 = _Vp[0]*t31;            const RealType t39 = t38*z[0];                const RealType t41 = 0.5*t37*t39;
01084     const RealType t42 = _V3[0]*t31;            const RealType t43 = t42*z[0];                const RealType t44 = t35+t41-t43;
01085     const RealType t45 = t44*t44;               const RealType t48 = 1/t25;                   const RealType t51 = t48*t22;
01086     const RealType t54 = 1/t11/t10*t9;          const RealType t55 = _Vm[0]*_Vm[0];           const RealType t56 = 1/t4;
01087     const RealType t57 = t55*t56;               const RealType t60 = 0.5*t54*t57*t1;          const RealType t61 = t29*t55;
01088     const RealType t62 = t56*t1;                const RealType t64 = 0.5*t61*t62;             const RealType t66 = 1/t5/t4;
01089     const RealType t67 = _Vm[0]*t66;            const RealType t70 = 0.5*t30*t67*t1;          const RealType t72 = 0.5*t30*t32;
01090     const RealType t75 = 1/t18/t17*t16;         const RealType t76 = _Vp[0]*_Vp[0];           const RealType t77 = t76*t56;
01091     const RealType t80 = 0.5*t75*t77*t1;        const RealType t81 = t36*t76;                 const RealType t83 = 0.5*t81*t62;
01092     const RealType t84 = _Vp[0]*t66;            const RealType t87 = 0.5*t37*t84*t1;          const RealType t89 = 0.5*t37*t38;
01093     const RealType t90 = _V3[0]*t66;            const RealType t91 = t90*t1;                  const RealType t95 = t12+t19+t20+t21;
01094     const RealType t96 = t95*t95;               const RealType t97 = t96+t3;                  const RealType t98 = sqrt(t97);
01095     const RealType t101 = 1/t98/t97*t96;        const RealType t102 = t35+t41+t43;            const RealType t103 = t102*t102;
01096     const RealType t106 = 1/t98;                const RealType t109 = t106*t95;               const RealType t114 = t32*z[1];
01097     const RealType t116 = 0.5*t30*t114;         const RealType t117 = t38*z[1];               const RealType t119 = 0.5*t37*t117;
01098     const RealType t120 = t42*z[1];             const RealType t121 = t116+t119-t120;         const RealType t130 = t56*z[0]*z[1];
01099     const RealType t132 = 0.5*t54*t55*t130;     const RealType t134 = 0.5*t61*t130;           const RealType t137 = t66*z[0]*z[1];
01100     const RealType t139 = 0.5*t30*_Vm[0]*t137;  const RealType t142 = 0.5*t75*t76*t130;       const RealType t144 = 0.5*t81*t130;
01101     const RealType t147 = 0.5*t37*_Vp[0]*t137;  const RealType t149 = t90*z[0]*z[1];          const RealType t153 = t116+t119+t120;
01102     const RealType t163 = -0.5*t28*t44*t121+0.5*t48*t121*t44+0.5*t51*(-t132+t134-t139-t142+t144-t147+t149)-0.5*t101*t102*t153
01103                           +0.5*t106*t153*t102+0.5*t109*(-t132+t134-t139-t142+t144-t147-t149);
01104     const RealType t165 = 0.5*t30*_Vm[1];       const RealType t167 = 0.5*t37*_Vp[1];         const RealType t168 = t165+t167-_V3[1];
01105     const RealType t172 = t48*t168;             const RealType t175 = t54*_Vm[0];             const RealType t176 = t31*z[0];
01106     const RealType t180 = t29*_Vm[1];           const RealType t183 = t75*_Vp[0];             const RealType t187 = t36*_Vp[1];
01107     const RealType t190 = -0.5*t175*t176*_Vm[1]+0.5*t180*t33-0.5*t183*t176*_Vp[1]+0.5*t187*t39;
01108     const RealType t193 = t165+t167+_V3[1];     const RealType t197 = t106*t193;
01109     const RealType t202 = -0.5*t28*t44*t168+0.5*t172*t44+0.5*t51*t190-0.5*t101*t102*t193+0.5*t197*t102+0.5*t109*t190;
01110     const RealType t203 = t121*t121;            const RealType t210 = 0.5*t54*t57*t2;         const RealType t211 = t56*t2;
01111     const RealType t213 = 0.5*t61*t211;         const RealType t216 = 0.5*t30*t67*t2;         const RealType t219 = 0.5*t75*t77*t2;
01112     const RealType t221 = 0.5*t81*t211;         const RealType t224 = 0.5*t37*t84*t2;         const RealType t225 = t90*t2;
01113     const RealType t229 = t153*t153;            const RealType t243 = t31*z[1];       
01114     const RealType t254 = -0.5*t175*t243*_Vm[1]+0.5*t180*t114-0.5*t183*t243*_Vp[1]+0.5*t187*t117;
01115     const RealType t264 = -0.5*t28*t121*t168+0.5*t172*t121+0.5*t51*t254-0.5*t101*t153*t193+0.5*t197*t153+0.5*t109*t254;
01116     const RealType t265 = t168*t168;            const RealType t270 = _Vm[1]*_Vm[1];          const RealType t275 = _Vp[1]*_Vp[1];
01117     const RealType t280 = -0.5*t54*t270+0.5*t29*t270-0.5*t75*t275+0.5*t36*t275;               const RealType t283 = t193*t193;
01118      
01119     mat[0][0] = -0.5*t28*t45+0.5*t48*t45+0.5*t51*(-t60+t64-t70+t72-t80+t83-t87+t89+t91-t42)-0.5*t101*t103
01120                 +0.5*t106*t103+0.5*t109*(-t60+t64-t70+t72-t80+t83-t87+t89-t91+t42);
01121     mat[0][1] = t163;
01122     mat[0][2] = t202;
01123     mat[1][0] = t163;
01124     mat[1][1] = -0.5*t28*t203+0.5*t48*t203+0.5*t51*(-t210+t213-t216+t72-t219+t221-t224+t89+t225-t42)-0.5*t101*t229
01125                 +0.5*t106*t229+0.5*t109*(-t210+t213-t216+t72-t219+t221-t224+t89-t225+t42);
01126     mat[1][2] = t264;
01127     mat[2][0] = t202;
01128     mat[2][1] = t264;
01129     mat[2][2] = -0.5*t28*t265+0.5*t48*t265+0.5*t51*t280-0.5*t101*t283+0.5*t106*t283+0.5*t109*t280;                                     
01130     
01131     for (int i=0; i<3; i++)                                          
01132       for (int j=0; j<3; j++)                                        
01133         if ( aol::isNaN( mat[i][j] ) ) cerr << "\n******* Found NAN in HessMat!! ******* ";
01134                                                                      
01135   }                                                                  
01136   
01137 };
01138 
01139 /** class for computing several terms for the anisotropy
01140     \f$\tilde{\gamma}(z) = \left( \gamma^p(z) + \epsilon^p \right)^{\frac{1}{p}}\f$
01141     with \f$\gamma(z) = \left( z_1^q + z_2^q \right)^{\frac{1}{q}}\f$.
01142     */
01143 template <typename RealType>
01144 class LpAnisotropy {
01145 private:
01146   RealType _p, _q, _epsilon;
01147 public:
01148   LpAnisotropy ( RealType P, RealType Q, RealType Epsilon ) : _p ( P ), _q ( Q ), _epsilon ( Epsilon ) { }
01149 
01150   RealType p() const { return _p; }
01151 
01152   void implicitPart ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z,
01153                       aol::Matrix22<RealType> &mat ) const {
01154     implicitPart ( z, mat );
01155   }
01156 
01157   void explicitPart ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z,
01158                       aol::Matrix22<RealType> &mat ) const {
01159     explicitPart ( z, mat );
01160   }
01161 
01162   /** factor of the qc::WillmoreStiffImpOp:
01163     *  \f$\frac{\left[ \gamma^p(z) \right]_{zz} }{p \|z\|^{p-1}_{\gamma,\epsilon,p}}\f$
01164     */
01165   void implicitPart ( const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
01166     aol::Vec2<RealType> Z ( aol::Abs ( z[0] ), aol::Abs ( z[1] ) );
01167     RealType a = pow ( Z[0], _q ) + pow ( Z[1], _q );
01168 
01169     if ( a != 0 ) {
01170       RealType factor1 = 0.;
01171       RealType factor2 = 0.;
01172 
01173       if ( _p != _q ) {
01174         factor1 = ( _p - _q ) * pow ( a, _p / _q - 2. );
01175         factor2 = pow ( a, _p / _q - 1. ) * ( _q - 1. );
01176       }
01177 
01178       mat[0][0] = factor1 * pow ( Z[0], 2 * ( _q - 1. ) ) + factor2 * pow ( Z[0], _q - 2. );
01179       mat[1][1] = factor1 * pow ( Z[1], 2 * ( _q - 1. ) ) + factor2 * pow ( Z[1], _q - 2. );
01180 
01181       mat[1][0] = mat[0][1] = factor1 * pow ( Z[0] * Z[1], _q - 1. ) * aol::signum1at0 ( z[0] ) * aol::signum1at0 ( z[1] );
01182 
01183       // Now mat = 0.5 * gamma^2_zz
01184       mat /= pow ( pow ( a, _p / _q ) + pow ( _epsilon, _p ), ( _p - 1. ) / _p ) ;
01185     } else {
01186       mat[0][0] = mat[1][1] = mat[1][0] = mat[1][0] = 0.;
01187     }
01188   }
01189 
01190   /** factor of the qc::WillmoreStiffExpOp:
01191     * \f$\frac{p-1}{p^2} \frac{\gamma^p(z)_z \otimes \gamma^p(z)_z}{\|z\|^{2p-1}_{\gamma,\epsilon,p}}\f$ 
01192     */
01193   // I've removed a little bug from Marc (at least I think so) and optimized the code a little bit.
01194   // For the latest version of Marc's computation see revision 3738.
01195   void explicitPart ( const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
01196     aol::Vec2<RealType> Z ( aol::Abs ( z[0] ), aol::Abs ( z[1] ) );
01197     RealType a = pow ( Z[0], _q ) + pow ( Z[1], _q );
01198 
01199     if ( a != 0 ) {
01200       mat[0][0] = pow ( Z[0], 2. * ( _q - 1. ) );
01201       mat[1][1] = pow ( Z[1], 2. * ( _q - 1. ) );
01202 
01203       mat[1][0] = mat[0][1] = pow ( Z[0], _q - 1. ) * pow ( Z[1], _q - 1. ) * aol::signum1at0 ( z[0] ) * aol::signum1at0 ( z[1] );
01204 
01205       RealType gammaFactor = pow ( a, 2. * ( _p / _q - 1. ) );
01206       RealType normSqr = pow ( pow ( a, _p / _q ) + pow ( _epsilon, _p ), ( 2 * _p - 1. ) / _p );
01207 
01208       mat *= gammaFactor * ( _p - 1. ) / normSqr;
01209     } else {
01210       mat[0][0] = mat[1][1] = mat[1][0] = mat[1][0] = 0.;
01211     }
01212   }
01213 
01214   //! evaluating the function itself \f$\gamma(z) = \left( z_1^q + z_2^q \right)^{\frac{1}{q}}\f$
01215   void gammaFirstDerivative ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
01216     gammaFirstDerivative ( z, v );
01217   }
01218 
01219   /** method computes:
01220    * \f$ \frac{ \gamma^{p-q}(z) \left( |z_i|^{q-1} sgn(z_i) \right)_i }{ \| z \|^{p-1}_{\gamma,\epsilon,p} }  \f$ 
01221    */
01222   void gammaFirstDerivative ( const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
01223     aol::Vec2<RealType> Z ( aol::Abs ( z[0] ), aol::Abs ( z[1] ) );
01224     RealType a = pow ( Z[0], _q ) + pow ( Z[1], _q );
01225     if ( a != 0 ) {
01226       v[0] = pow ( Z[0], _q - 1. ) * aol::signum1at0 ( z[0] );
01227       v[1] = pow ( Z[1], _q - 1. ) * aol::signum1at0 ( z[1] );
01228 
01229       v *= pow ( a, _p / _q - 1. ) *  pow ( pow ( a, _p / _q ) + pow ( _epsilon, _p ), ( 1. - _p ) / _p );
01230     } else {
01231       v[0] = v[1] = 0.;
01232     }
01233   }
01234 
01235   //! evaluating the function itself \f$\gamma(z) = \left( z_1^q + z_2^q \right)^{\frac{1}{q}}\f$
01236   RealType gammaNorm ( const aol::Vec2<RealType> &z ) const {
01237     return pow ( pow ( aol::Abs ( z[0] ), _q ) + pow ( aol::Abs ( z[1] ), _q ) + pow ( _epsilon, _q ), 1. / _q );
01238   }
01239 };
01240 
01241 
01242 
01243 /* *******************************************************************************
01244  * gamma = Lp-Norm in 3d
01245  * ******************************************************************************* */
01246 /*! class for computing several terms for the anisotropy
01247     \f$\tilde{\gamma}(z) = \left( \gamma^p(z) + \epsilon^p \right)^{\frac{1}{p}}\f$
01248     with \f$\gamma(z) = \left( z_1^q + z_2^q \right)^{\frac{1}{q}}\f$.
01249     */
01250 template <typename RealType>
01251 class LpAnisotropy3d {
01252 private:
01253   RealType _p, _q, _epsilon;
01254 public:
01255   LpAnisotropy3d ( RealType P, RealType Q, RealType Epsilon ) : _p ( P ), _q ( Q ), _epsilon ( Epsilon ) {
01256     if ( Q < 1 ) cerr << aol::color::red << "\n\nWARNING: Q is smaller 1, there might be singularities!!\n";
01257   }
01258 
01259 RealType p() const { return _p; }
01260 
01261   void implicitPart ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z,
01262                       aol::Matrix33<RealType> &mat ) const {
01263     implicitPart ( z, mat );
01264   }
01265 
01266   void explicitPart ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z,
01267                       aol::Matrix33<RealType> &mat ) const {
01268     explicitPart ( z, mat );
01269   }
01270 
01271   /** factor of the qc::WillmoreStiffImpOp:
01272     *  \f$\frac{\left[ \gamma^p(z) \right]_{zz} }{p \|z\|^{p-1}_{\gamma,\epsilon,p}}\f$
01273     */
01274   void implicitPart ( const aol::Vec3<RealType> &z, aol::Matrix33<RealType> &mat ) const {
01275     aol::Vec3<RealType> Z ( aol::Abs ( z[0] ), aol::Abs ( z[1] ), aol::Abs ( z[2] ) );
01276     RealType a = pow ( Z[0], _q ) + pow ( Z[1], _q ) + pow ( Z[2], _q );
01277 
01278     if ( a != 0 ) {
01279       RealType factor2 = 0.;
01280       RealType factor1 = 0.;
01281       if ( _p != _q ) {
01282         factor1 = ( _p - _q ) * pow ( a, _p / _q - 2. );
01283         factor2 = ( _q - 1. ) * pow ( a, _p / _q - 1. );
01284       }
01285 
01286       mat[0][0] = factor1 * pow ( Z[0], 2 * ( _q - 1. ) ) + factor2 * pow ( Z[0], _q - 2. );
01287       mat[1][1] = factor1 * pow ( Z[1], 2 * ( _q - 1. ) ) + factor2 * pow ( Z[1], _q - 2. );
01288       mat[2][2] = factor1 * pow ( Z[2], 2 * ( _q - 1. ) ) + factor2 * pow ( Z[2], _q - 2. );
01289 
01290       mat[1][0] = mat[0][1] = factor1 * pow ( Z[0] * Z[1], _q - 1. ) * aol::signum1at0 ( z[0] ) * aol::signum1at0 ( z[1] );
01291       mat[2][0] = mat[0][2] = factor1 * pow ( Z[0] * Z[2], _q - 1. ) * aol::signum1at0 ( z[0] ) * aol::signum1at0 ( z[2] );
01292       mat[2][1] = mat[1][2] = factor1 * pow ( Z[2] * Z[1], _q - 1. ) * aol::signum1at0 ( z[2] ) * aol::signum1at0 ( z[1] );
01293 
01294       // Now mat = 0.5 * gamma^2_zz
01295       mat /= pow ( pow ( a, _p / _q ) + pow ( _epsilon, _p ), ( _p - 1. ) / _p ) ;
01296     } else {
01297       mat[0][0] = mat[0][1] = mat[0][2] = 0.;
01298       mat[1][0] = mat[1][1] = mat[1][2] = 0.;
01299       mat[2][0] = mat[2][1] = mat[2][2] = 0.;
01300     }
01301   }
01302 
01303   /** factor of the qc::WillmoreStiffExpOp:
01304     * \f$\frac{p-1}{p^2} \frac{\gamma^p(z)_z \otimes \gamma^p(z)_z}{\|z\|^{2p-1}_{\gamma,\epsilon,p}}\f$ 
01305     */
01306   void explicitPart ( const aol::Vec3<RealType> &z, aol::Matrix33<RealType> &mat ) const {
01307     aol::Vec3<RealType> Z ( aol::Abs ( z[0] ), aol::Abs ( z[1] ), aol::Abs ( z[2] ) );
01308     RealType a = pow ( Z[0], _q ) + pow ( Z[1], _q ) + pow ( Z[2], _q );
01309 
01310     if ( a != 0 ) {
01311       mat[0][0] = pow ( Z[0], 2. * ( _q - 1. ) );
01312       mat[1][1] = pow ( Z[1], 2. * ( _q - 1. ) );
01313       mat[2][2] = pow ( Z[2], 2. * ( _q - 1. ) );
01314 
01315       mat[1][0] = mat[0][1] = pow ( Z[0], _q - 1. ) * pow ( Z[1], _q - 1. ) * aol::signum1at0 ( z[0] ) * aol::signum1at0 ( z[1] );
01316       mat[2][0] = mat[0][2] = pow ( Z[0], _q - 1. ) * pow ( Z[2], _q - 1. ) * aol::signum1at0 ( z[0] ) * aol::signum1at0 ( z[2] );
01317       mat[1][2] = mat[2][1] = pow ( Z[2], _q - 1. ) * pow ( Z[1], _q - 1. ) * aol::signum1at0 ( z[2] ) * aol::signum1at0 ( z[1] );
01318 
01319       RealType gammaFactor = pow ( a, 2. * ( _p / _q - 1. ) );
01320       RealType normSqr = pow ( pow ( a, _p / _q ) + pow ( _epsilon, _p ), ( 2 * _p - 1. ) / _p );
01321 
01322       mat *= gammaFactor * ( _p - 1. ) / normSqr;
01323     } else {
01324       mat[0][0] = mat[0][1] = mat[0][2] = 0.;
01325       mat[1][0] = mat[1][1] = mat[1][2] = 0.;
01326       mat[2][0] = mat[2][1] = mat[2][2] = 0.;
01327     }
01328   }
01329 
01330 
01331   //! evaluating the function itself \f$\gamma(z) = \left( z_1^q + z_2^q \right)^{\frac{1}{q}}\f$
01332   RealType gamma ( const aol::Vec3<RealType> &z ) const {
01333     return pow ( pow ( aol::Abs ( z[0] ), _q ) + pow ( aol::Abs ( z[1] ), _q )
01334                  + pow ( aol::Abs ( z[2] ), _q ) + pow ( _epsilon, _q ),  1 / _q );
01335   }
01336 
01337 
01338   void gammaFirstDerivative ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
01339     gammaFirstDerivative ( z, v );
01340   }
01341 
01342   /** method computes:
01343    * \f$ \frac{ \gamma^{p-q}(z) \left( |z_i|^{q-1} sgn(z_i) \right)_i }{ \| z \|^{p-1}_{\gamma,\epsilon,p} }  \f$ 
01344      */
01345   void gammaFirstDerivative ( const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
01346     aol::Vec3<RealType> Z ( aol::Abs ( z[0] ), aol::Abs ( z[1] ), aol::Abs ( z[2] ) );
01347     RealType a = pow ( Z[0], _q ) + pow ( Z[1], _q ) + pow ( Z[2], _q );
01348 
01349     if ( a != 0 ) {
01350       v[0] = pow ( Z[0], _q - 1. ) * aol::signum1at0 ( z[0] );
01351       v[1] = pow ( Z[1], _q - 1. ) * aol::signum1at0 ( z[1] );
01352       v[2] = pow ( Z[2], _q - 1. ) * aol::signum1at0 ( z[2] );
01353 
01354       v *= pow ( a, _p / _q - 1. ) *  pow ( pow ( a, _p / _q ) + pow ( _epsilon, _p ), ( 1. - _p ) / _p );
01355     } else {
01356       v[0] = v[1] = v[2] = 0.;
01357     }
01358   }
01359 
01360   //! evaluating the function itself \f$\gamma(z) = \left( z_1^q + z_2^q \right)^{\frac{1}{q}}\f$
01361   RealType gammaNorm ( const aol::Vec3<RealType> &z ) const {
01362     return pow ( pow ( aol::Abs ( z[0] ), _q ) + pow ( aol::Abs ( z[1] ), _q )
01363                  + pow ( aol::Abs ( z[2] ), _q ) + pow ( _epsilon, _q ), 1. / _q );
01364   }
01365 };
01366 
01367 
01368 
01369 //! Class for a pedestal-like Wulffshape in 2D.
01370 template <typename RealType>
01371 class Pedestal2d {
01372 private:
01373   RealType _alpha;                  //! angle of the walls
01374   RealType _alphaHalf;              //! 0.5 * alpha (often needed)
01375   RealType _deltaAbs;               //! Regularization-parameter for the absolute-value
01376   RealType _deltaAbsSqr;            //! Square of Delta (very complicated...)
01377   RealType _deltaRad;               //! Regularization-parameter for the radius
01378   RealType _deltaRadSqr;
01379   RealType _c1, _c2;                //! two constants for the computation (Franck-Diagram-Ansatz)
01380   RealType _sinAlpha, _cosAlpha;    //! two more constants for the computation (WulffShape-Ansatz)
01381 
01382   RealType _phi;                    //! the angle of z in polar coordinates
01383   RealType _r;                      //! |z|   [ z = r (cos phi, sin phi) ]
01384 
01385 public:
01386   Pedestal2d ( RealType alpha, RealType deltaAbs, RealType deltaRad ) :
01387       _alpha ( alpha ), _deltaAbs ( deltaAbs ), _deltaRad ( deltaRad ) {
01388     _alphaHalf = 0.5 * _alpha;
01389     _deltaAbsSqr = _deltaAbs * _deltaAbs;
01390     _deltaRadSqr = _deltaRad * _deltaRad;
01391     _c1 = cos ( _alphaHalf );
01392     _c2 = sin ( _alphaHalf );
01393     _sinAlpha = sin ( _alpha );
01394     _cosAlpha = cos ( _alpha );
01395   }
01396 
01397   //! Evaluating the function itself, i.e. the regularization of
01398   //! \f$ \gamma(z) = \frac12 |z| max\left\{ \frac{\sin(\frac{\alpha}{2}+\varphi)}{\cos \frac{\alpha}{2}}, \frac{\cos(\frac{\alpha}{2}+\varphi)}{\sin \frac{\alpha}{2}} \right\}  \f$.
01399   RealType gammaNormFromDistFct ( const aol::Vec2<RealType> &z ) const {
01400     double r = sqrt ( z.normSqr() + _deltaRadSqr );
01401     double phi = asin ( z[1] / r );    // in [-pi/2, pi/2]
01402 
01403     double CSin = sin ( _alphaHalf + phi );
01404     double CCos = cos ( _alphaHalf + phi );
01405 
01406     double term1 = sqrt ( aol::Sqr ( CSin / _c1 + CCos / _c2 ) + _deltaAbsSqr );
01407     double term2 = sqrt ( aol::Sqr ( CSin / _c1 - CCos / _c2 ) + _deltaAbsSqr );
01408     return 0.5 * r * ( term1 + term2 );
01409   }
01410 
01411   RealType gammaNormFromWulffShape ( const aol::Vec2<RealType> &z ) const {
01412     double r = sqrt ( z.normSqr() + _deltaRadSqr );
01413     double phi = asin ( z[1] / r );    // in [-pi/2, pi/2]
01414 
01415     double CSin = sin ( phi );
01416     double CCos = cos ( phi );
01417 
01418     double term1 = sqrt ( aol::Sqr ( CCos / _sinAlpha ) + _deltaAbsSqr );
01419     double term2 = sqrt ( aol::Sqr ( CSin - CCos * _cosAlpha / _sinAlpha ) + _deltaAbsSqr );
01420 
01421     return r * ( term1 + term2 );
01422   }
01423 
01424   RealType gammaNorm ( const aol::Vec2<RealType> &z ) const {
01425     return gammaNormFromWulffShape ( z );
01426   }
01427 
01428   RealType evaluate ( const aol::Vec2<RealType> &z ) const {
01429     return gammaNorm ( z );
01430   }
01431 
01432   //! evaluating the first derivative: \f$ ... \f$
01433   void gammaFirstDerivative ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
01434     gammaFirstDerivative ( z, v );
01435   }
01436 
01437   //! evaluating the first derivative: \f$ ... \f$
01438   void gammaFirstDerivative ( const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
01439     double r       = z.normSqr();
01440     double r_delta = sqrt ( r + _deltaRadSqr );
01441     double phi     = asin ( z[1] / r_delta );         //! in [-pi/2, pi/2]
01442 
01443     double CSin = sin ( phi );
01444     double CCos = cos ( phi );
01445 
01446     //! \f$ \frac{d}{dr} \gamma \f$
01447     double dGammaR = gammaNorm ( z );
01448     dGammaR *= r / ( r_delta * r_delta * r_delta );   // here is already 1/r_delta from the whole formula incorporated
01449 
01450     double C = _cosAlpha / _sinAlpha;
01451     double S = sqrt ( r_delta * r_delta - z[1] * z[1] );
01452 
01453     //! \f$ \frac{d}{d \varphi} \gamma \f$
01454     double dGammaPhi = ( CSin - C * CCos ) * ( CCos + C * CSin ) / sqrt ( aol::Sqr ( CSin - C * CCos ) + _deltaAbsSqr );
01455     dGammaPhi -= CSin * CCos / ( _sinAlpha * _sinAlpha * sqrt ( C * C + _deltaAbsSqr ) );
01456     //! Some further computing
01457     dGammaPhi /= r_delta;
01458 
01459     v[0] = z[0] * ( dGammaR - dGammaPhi * z[1] / S );
01460     v[1] = z[1] * dGammaR + dGammaPhi * S;
01461 
01462   }
01463 };
01464 
01465 //! Class for a pedestal-like Wulffshape.
01466 template <typename RealType>
01467 class Pedestal3d : public Anisotropy3dGraphInterface<RealType, Pedestal3d<RealType> > {
01468 
01469 public:  
01470   using Anisotropy3dGraphInterface<RealType, Pedestal3d<RealType> >::gammaNorm;
01471   using Anisotropy3dGraphInterface<RealType, Pedestal3d<RealType> >::gammaFirstDerivative;
01472 
01473 private:
01474   RealType _alpha;        //! angle of the walls
01475   RealType _alphaHalf;    //! 0.5 * alpha (often needed)
01476   RealType _deltaAbs;     //! Regularization-parameter for the absolute-value
01477   RealType _deltaAbsSqr;  //! Square of Delta (very complicated...)
01478   RealType _deltaRad;     //! Regularization-parameter for the radius
01479   RealType _deltaRadSqr;
01480   RealType _c1, _c2;      //! two constants for the computation
01481   RealType _sinAlpha, _cosAlpha;    //! two more constants for the computation (WulffShape-Ansatz)
01482 
01483 public:
01484   Pedestal3d ( RealType alpha, RealType deltaAbs, RealType deltaRad ) :
01485       Anisotropy3dGraphInterface<RealType, Pedestal3d<RealType> > ( ),
01486       _alpha ( alpha ), _deltaAbs ( deltaAbs ), _deltaRad ( deltaRad ) {
01487     _alphaHalf = 0.5 * _alpha;
01488     _deltaAbsSqr = _deltaAbs * _deltaAbs;
01489     _deltaRadSqr = _deltaRad * _deltaRad;
01490     _c1 = cos ( _alphaHalf );
01491     _c2 = sin ( _alphaHalf );
01492     _sinAlpha = sin ( _alpha );
01493     _cosAlpha = cos ( _alpha );
01494   }
01495 
01496   //! Evaluating the function itself, i.e. the regularization of
01497   //! \f$ \gamma(z) = \frac12 |z| max\left\{ \frac{\sin(\frac{\alpha}{2}+\varphi)}{\cos \frac{\alpha}{2}}, \frac{\cos(\frac{\alpha}{2}+\varphi)}{\sin \frac{\alpha}{2}} \right\}  \f$.
01498   RealType gammaNormFromFranckDiagram ( const aol::Vec3<RealType> &z ) const {
01499     double r   = sqrt ( z.normSqr() + _deltaRadSqr );
01500     double phi = asin ( z[2] / r );    //! in [-pi/2, pi/2]
01501 
01502     double CSin = sin ( _alphaHalf + phi );
01503     double CCos = cos ( _alphaHalf + phi );
01504 
01505     double term1 = sqrt ( aol::Sqr ( CSin / _c1 + CCos / _c2 ) + _deltaAbsSqr );
01506     double term2 = sqrt ( aol::Sqr ( CSin / _c1 - CCos / _c2 ) + _deltaAbsSqr );
01507     return 0.5 * r * ( term1 + term2 );
01508   }
01509 
01510   RealType gammaNormFromWulffShape ( const aol::Vec3<RealType> &z ) const {
01511     double r   = sqrt ( z.normSqr() + _deltaRadSqr );
01512     double phi = asin ( z[2] / r );    //! in [-pi/2, pi/2]
01513 
01514     double CSin = sin ( phi );
01515     double CCos = cos ( phi );
01516 
01517     double term1 = sqrt ( aol::Sqr ( CCos / _sinAlpha ) + _deltaAbsSqr );
01518     double term2 = sqrt ( aol::Sqr ( CSin - CCos * _cosAlpha / _sinAlpha ) + _deltaAbsSqr );
01519     return r * ( term1 + term2 );
01520   }
01521 
01522   RealType gammaNorm ( const aol::Vec3<RealType> &z ) const {
01523     return gammaNormFromWulffShape ( z );
01524   }
01525 
01526   //! evaluating the first derivative: (too long for TeXing)
01527   void gammaFirstDerivativeFromFranckDiagram ( const aol::Vec3<RealType> &z,
01528                                                aol::Vec3<RealType> &v ) const {
01529     double r   = sqrt ( z.normSqr() + _deltaRadSqr );
01530     double phi = 0.;
01531     if ( r != 0. ) phi = asin ( z[2] / r );  //! in [-pi/2, pi/2]
01532 
01533     double CSin = sin ( _alphaHalf + phi );
01534     double CCos = cos ( _alphaHalf + phi );
01535 
01536     //! \f$ \frac{d}{dr} \gamma \f$
01537     double dGammaR = gammaNorm ( z );
01538     if ( r == 0 )  dGammaR = 0.;
01539     else         dGammaR /= r * r;    // one 1/r is already from the derivative of r
01540 
01541     //! \f$ \frac{d}{d \varphi} \gamma \f$
01542     double dGammaPhi = ( CSin / _c1 + CCos / _c2 ) * ( CCos / _c1 - CSin / _c2 ) / sqrt ( aol::Sqr ( CSin / _c1 + CCos / _c2 ) + _deltaAbsSqr );
01543     dGammaPhi += ( CSin / _c1 - CCos / _c2 ) * ( CCos / _c1 + CSin / _c2 ) / sqrt ( aol::Sqr ( CSin / _c1 - CCos / _c2 ) + _deltaAbsSqr );
01544     dGammaPhi *= 0.5 * r;
01545     //! further summing-up
01546     dGammaPhi /= r * r * r * sqrt ( 1. - z[1] * z[1] / ( r * r ) );
01547 
01548 
01549     v[0] = z[0] * ( dGammaR - z[1] * dGammaPhi );
01550     v[1] = z[1] * dGammaR + ( r * r - z[1] * z[1] ) * dGammaPhi;
01551     v[2] = z[2] * ( dGammaR - z[1] * dGammaPhi );
01552   }
01553 
01554 
01555   //! evaluating the first derivative: (too long for TeXing)
01556   void gammaFirstDerivativeFromWulffShape ( const aol::Vec3<RealType> &z,
01557                                             aol::Vec3<RealType> &v ) const {
01558     double r   = sqrt ( z.normSqr() + _deltaRadSqr );
01559     double phi = asin ( z[2] / r );    //! in [-pi/2, pi/2]
01560 
01561     double CSin = sin ( phi );
01562     double CCos = cos ( phi );
01563 
01564     //! \f$ \frac{d}{dr} \gamma \f$
01565     double dGammaR = gammaNorm ( z );
01566     dGammaR /= r * r;
01567 
01568     double C = _cosAlpha / _sinAlpha;
01569     double S = sqrt ( r * r - z[2] * z[2] );
01570 
01571     //! \f$ \frac{d}{d \varphi} \gamma \f$
01572     double dGammaPhi = ( CSin - C * CCos ) * ( CCos + C * CSin ) / sqrt ( aol::Sqr ( CSin - C * CCos ) + _deltaAbsSqr );
01573     dGammaPhi -= CSin * CCos / ( _sinAlpha * _sinAlpha * sqrt ( C * C + _deltaAbsSqr ) );
01574     //! Some further computing
01575     dGammaPhi /= r;
01576 
01577     v[0] = z[0] * ( dGammaR - dGammaPhi * z[2] / S );
01578     v[1] = z[1] * ( dGammaR - dGammaPhi * z[2] / S );
01579     v[2] = z[2] * dGammaR + dGammaPhi * S;
01580   }
01581 
01582   //! evaluating the first derivative: (too long for TeXing)
01583   void gammaFirstDerivative ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
01584     gammaFirstDerivative ( z, v );
01585   }
01586 
01587 
01588   void gammaFirstDerivative ( const aol::Vec3<RealType> &z,
01589                               aol::Vec3<RealType> &v ) const {
01590     gammaFirstDerivativeFromWulffShape ( z, v );
01591   }
01592 
01593 };
01594 
01595 
01596 
01597 //! Class for a pentagon-like Wulffshape in 2D.
01598 template <typename RealType>
01599 class Pentagon2d {
01600 private:
01601   RealType _deltaAbs;     //! Regularization-parameter for the absolute-value
01602   RealType _deltaAbsSqr;  //! Square of Delta (very complicated...)
01603   RealType _deltaRad;     //! Regularization-parameter for the radius
01604   RealType _deltaRadSqr;
01605 
01606 public:
01607   Pentagon2d ( RealType deltaAbs, RealType deltaRad ) :
01608       _deltaAbs ( deltaAbs ), _deltaRad ( deltaRad ) {
01609     _deltaAbsSqr = _deltaAbs * _deltaAbs;
01610     _deltaRadSqr = _deltaRad * _deltaRad;
01611   }
01612 
01613 
01614   //! Evaluating the function itself, i.e. the regularization of
01615   //! \f$ \gamma(z) = r_{\delta} \left\{ 1 + \sqrt{\sin^2(\frac52(\phi + \frac{\pi}{2})) + \delta^2 } \right\}  \f$.
01616   RealType gammaNorm ( const aol::Vec2<RealType> &z ) const {
01617     double r   = sqrt ( z.normSqr() + _deltaRadSqr );
01618     double phi = asin ( z[1] / r );    //! in [-pi/2, pi/2]
01619 
01620     double CSin = aol::Sqr ( sin ( 2.5 * ( phi + 1.570796326499999999979674536465523715378 ) ) );
01621     return r * ( 1. + sqrt ( CSin + _deltaAbsSqr ) );
01622   }
01623 
01624   //! evaluating the first derivative: \f$ ... \f$
01625   void gammaFirstDerivative ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
01626     gammaFirstDerivative ( z, v );
01627   }
01628 
01629   //! evaluating the first derivative: \f$ ... \f$
01630   void gammaFirstDerivative ( const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
01631     double r   = sqrt ( z.normSqr() + _deltaRadSqr );
01632     double phi = asin ( z[1] / r );    //! in [-pi/2, pi/2]
01633 
01634     double CSin = sin ( 2.5 * ( phi + 1.570796326499999999979674536465523715378 ) );
01635     double CCos = cos ( 2.5 * ( phi + 1.570796326499999999979674536465523715378 ) );
01636 
01637     double root  = sqrt ( CSin * CSin + _deltaAbsSqr );
01638     double term1 = ( 1. + root ) / r;
01639     double term2 = CSin * CCos / ( root * r );
01640 
01641     v[0] = z[0] * term1 - 2.5 * term2 * z[0] * z[1] / sqrt ( r * r - z[1] * z[1] );
01642     v[1] = z[1] * term1 + 2.5 * term2 * sqrt ( r * r - z[1] * z[1] );
01643   }
01644 };
01645 
01646 
01647 //! Class for a rotated pentagon-like Wulffshape in 3D.
01648 template <typename RealType>
01649 class Pentagon3d : public Anisotropy3dGraphInterface<RealType, Pentagon3d<RealType> > {
01650 public:  
01651   using Anisotropy3dGraphInterface<RealType, Pentagon3d<RealType> >::gammaNorm;
01652   using Anisotropy3dGraphInterface<RealType, Pentagon3d<RealType> >::gammaFirstDerivative;
01653 
01654 private:
01655   RealType _deltaAbs;     //! Regularization-parameter for the absolute-value
01656   RealType _deltaAbsSqr;  //! Square of Delta (very complicated...)
01657   RealType _deltaRad;     //! Regularization-parameter for the radius
01658   RealType _deltaRadSqr;
01659 
01660 public:
01661   Pentagon3d ( RealType deltaAbs, RealType deltaRad ) :
01662       Anisotropy3dGraphInterface<RealType, Pentagon3d<RealType> > ( ),
01663       _deltaAbs ( deltaAbs ), _deltaRad ( deltaRad ) {
01664     _deltaAbsSqr = _deltaAbs * _deltaAbs;
01665     _deltaRadSqr = _deltaRad * _deltaRad;
01666   }
01667 
01668   //! Evaluating the function itself, i.e. the regularization of
01669   //! \f$ \gamma(z) = r_{\delta} \left\{ 1 + \sqrt{\sin^2(\frac52(\phi + \frac{\pi}{2})) + \delta^2 } \right\}  \f$.
01670   RealType gammaNorm ( const aol::Vec3<RealType> &z ) const {
01671     RealType r   = sqrt ( z.normSqr() + _deltaRadSqr );
01672     RealType phi = asin ( z[2] / r );    //! in [-pi/2, pi/2]
01673 
01674     RealType CSin = aol::Sqr ( sin ( 2.5 * ( phi + M_PI / 2. ) ) );
01675     return r * ( 1. + sqrt ( CSin + _deltaAbsSqr ) );
01676   }
01677 
01678   //! evaluating the first derivative: \f$ ... \f$
01679   void gammaFirstDerivative ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
01680     gammaFirstDerivative ( z, v );
01681   }
01682 
01683   //! evaluating the first derivative: \f$ ... \f$
01684   void gammaFirstDerivative ( const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
01685     RealType r   = sqrt ( z.normSqr() + _deltaRadSqr );
01686     RealType phi = asin ( z[2] / r );    //! in [-pi/2, pi/2]
01687 
01688     RealType CSin = sin ( 2.5 * ( phi + M_PI / 2. ) );
01689     RealType CCos = cos ( 2.5 * ( phi + M_PI / 2. ) );
01690 
01691     RealType root1 = sqrt ( CSin * CSin + _deltaAbsSqr );
01692     RealType root2 = sqrt ( r * r - z[2] * z[2] );
01693     RealType term1 = ( 1. + root1 ) / r;
01694     RealType term2 = CSin * CCos / ( root1 * r );
01695 
01696     v[0] = z[0] * term1 - 2.5 * term2 * z[0] * z[2] / root2;
01697     v[1] = z[1] * term1 - 2.5 * term2 * z[1] * z[2] / root2;
01698     v[2] = z[2] * term1 + 2.5 * term2 * root2;
01699   }
01700 };
01701 
01702 
01703 //! Class for a rotated polygonal Wulffshape in 3D generated by the function
01704 //! \f$ \gamma(r,\varphi,\vartheta) = r ( a + b| \sin(c(\varphi+d)) | ) \f$.
01705 //! Regularized version: \f$ \gamma(r,\varphi,\vartheta) = r ( a + b \sqrt( \sin^2(c(\varphi+d)) + \delta^2) ) \f$
01706 //! Special cases:
01707 //! \f$ a=1, b=3, c=\frac{3}{2}, d=\frac{\pi}{2}: a rotated triangle                 \f$
01708 //! \f$ a=1, b=3, c=2,           d=0:             a rotated square (diamond like)    \f$
01709 //! \f$ a=1, b=1, c=\frac{5}{2}, d=\frac{\pi}{2}: a rotated pentagon                 \f$
01710 //! \f$ a=1, b=1, c=3,           d=\frac{\pi}{2}: a rotated hexagon                  \f$
01711 //! \f$ a=1, b=1, c=\frac{7}{2}, d=\frac{\pi}{2}: a rotated heptagon                 \f$
01712 //! \f$ a=1, b=1, c=4,           d=0:             a rotated octagon                  \f$
01713 template <typename RealType>
01714 class SinePolygonRotated3d : public Anisotropy3dGraphInterface<RealType, SinePolygonRotated3d<RealType> > {
01715 public:
01716   using Anisotropy3dGraphInterface<RealType, SinePolygonRotated3d<RealType> >::gammaNorm;
01717   using Anisotropy3dGraphInterface<RealType, SinePolygonRotated3d<RealType> >::gammaFirstDerivative;
01718 
01719 private:
01720   const RealType _a;            //! the parameters for specifying the polygon
01721   const RealType _b;
01722   const RealType _c;
01723   const RealType _d;
01724 
01725   const RealType _deltaAbs;     //! Regularization-parameter for the absolute-value
01726   const RealType _deltaRad;     //! Regularization-parameter for the radius
01727   RealType _deltaAbsSqr;        //! Square of DeltaAbs (very complicated...)
01728   RealType _deltaRadSqr;        //! uh, don't know ...
01729 
01730 public:
01731   SinePolygonRotated3d ( RealType a, RealType b, RealType c, RealType d, RealType deltaAbs, RealType deltaRad ) :
01732       Anisotropy3dGraphInterface<RealType, SinePolygonRotated3d<RealType> > ( ),
01733       _a ( a ), _b ( b ), _c ( c ), _d ( d ),
01734       _deltaAbs ( deltaAbs ), _deltaRad ( deltaRad ) {
01735     _deltaAbsSqr = _deltaAbs * _deltaAbs;
01736     _deltaRadSqr = _deltaRad * _deltaRad;
01737   }
01738 
01739   //! Evaluating the function itself, i.e. the regularization of
01740   //! \f$ \gamma(z) = r_{\delta} \left\{ a + b \sqrt{\sin^2( c (\phi + d)) + \delta^2 } \right\}  \f$.
01741   RealType gammaNorm ( const aol::Vec3<RealType> &z ) const {
01742     RealType r   = sqrt ( z.normSqr() + _deltaRadSqr );
01743     RealType phi = asin ( z[2] / r );    //! in [-pi/2, pi/2]
01744 
01745     RealType CSin = aol::Sqr ( sin ( _c * ( phi + _d ) ) );
01746     return r * ( _a + _b*sqrt ( CSin + _deltaAbsSqr ) );
01747   }
01748 
01749   //! evaluating the first derivative: \f$ ... \f$
01750   void gammaFirstDerivative ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
01751     gammaFirstDerivative ( z, v );
01752   }
01753 
01754   //! evaluating the first derivative: \f$ ... \f$
01755   void gammaFirstDerivative ( const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
01756     RealType r   = sqrt ( z.normSqr() + _deltaRadSqr );
01757     RealType phi = asin ( z[2] / r );    //! in [-pi/2, pi/2]
01758 
01759     RealType CSin = sin ( _c * ( phi + _d ) );
01760     RealType CCos = cos ( _c * ( phi + _d ) );
01761 
01762     RealType root1 = sqrt ( CSin * CSin + _deltaAbsSqr );
01763     RealType root2 = sqrt ( r * r - z[2] * z[2] );
01764     RealType term1 = ( _a + _b * root1 ) / r;
01765     RealType term2 = CSin * CCos / ( root1 * r );
01766 
01767     v[0] = z[0] * term1 - _b * _c * term2 * z[0] * z[2] / root2;
01768     v[1] = z[1] * term1 - _b * _c * term2 * z[1] * z[2] / root2;
01769     v[2] = z[2] * term1 + _b * _c * term2 * root2;
01770   }
01771 };
01772 
01773 
01774 
01775 
01776 
01777 /* *******************************************************************************
01778  * gamma = Lp-Norm in 3d (rotated)
01779  * ******************************************************************************* */
01780 
01781 template <typename RealType>
01782 class RotatedLpAnisotropy3d {
01783 private:
01784   RealType _p, _q, _epsilon;
01785   aol::Vec3<RealType> _v;        // direction, on which the x-axxis lies
01786   aol::Matrix33<RealType> _B;    // matrix to perform a basis-transformation
01787   aol::Matrix33<RealType> _BT;   // the transposed transformation-matrix
01788 public:
01789   RotatedLpAnisotropy3d ( RealType P, RealType Q, RealType Epsilon ) : _p ( P ), _q ( Q ), _epsilon ( Epsilon ) {
01790     if ( Q < 1 ) cerr << aol::color::red << "\n\nWARNING: Q is smaller 1, there might be singularities!!\n";
01791     _v[0] = 1.; _v[1] = 0.; _v[2] = 0.;
01792     calcB();
01793   }
01794 
01795   RotatedLpAnisotropy3d ( RealType P, RealType Q, RealType Epsilon, aol::Vec3<RealType> V ) :
01796       _p ( P ), _q ( Q ), _epsilon ( Epsilon ),  _v ( V ) {
01797     if ( Q < 1 ) cerr << aol::color::red << "\n\nWARNING: Q is smaller 1, there might be singularities!!\n";
01798     calcB();
01799   }
01800 
01801 
01802   // calc the matrix B for transforming the coordinate-system for the x-axxis
01803   // lying on v
01804   void calcB() {
01805     // the first row is the direction itself
01806     _B.setRow ( 0, _v );
01807     // the other directions must be orthogonal to _v, but then, they can be
01808     // chosen arbitrary, first choose y orthogonal to _v
01809     // TODO: arranging by size to avoid taking the smallest values (=> stability)
01810     aol::Vec3<RealType> y ( 0, 0, 0 );
01811     if ( _v[1] != 0 ) {
01812       if ( _v[2] != 0 ) {
01813         y[1] = -_v[2];
01814         y[2] = _v[1];
01815         y /= y.norm();
01816       } else y[2] = 1;
01817     } else y[1] = 1;
01818 
01819     _B.setRow ( 1, y );
01820 
01821     // now z = vector-product of v and y, both normed => z is normed too
01822     aol::Vec3<RealType> z ( _v.crossProduct ( y ) );
01823     _B.setRow ( 2, z );
01824 
01825     // finally save the transposed matrix
01826     _BT = _B;
01827     _BT.transpose();
01828   }
01829 
01830 
01831   void printB() {
01832     cerr << "Drehmatrix: \n";
01833     cerr << _BT;
01834   }
01835 
01836   // set the rotate-direction v
01837   void setRotateDirection ( const aol::Vec3<RealType> &v ) {
01838     // only anisotropic, if v is not 0
01839     // if v is not 0, it is very improbable, that v is equal to _v
01840     // (that means = v from the timestep before)
01841     // => it's not worthwile to check this
01842 
01843     if ( v[0] != 0 || v[1] != 0 || v[2] != 0 ) {
01844       _v = v / v.norm();
01845       calcB();
01846     } else cerr << "\n\n'anisotropy:************* VECTOR FOR ROTATING THE ANISOTROPY IS 0 !!!! ************** \n\n";
01847   }
01848 
01849   // set the rotate-direction v
01850   void setRotateDirection ( RealType v0, RealType v1, RealType v2 ) {
01851     if ( v0 != 0 || v1 != 0 || v2 != 0 ) {
01852       _v[0] = v0; _v[1] = v1; _v[2] = v2;
01853       _v /= _v.norm();
01854       calcB();
01855     } else cerr << "\n\n'anisotropy:************* VECTOR FOR ROTATING THE ANISOTROPY IS 0 !!!! ************** \n\n";
01856   }
01857 
01858 
01859 
01860 RealType p() const { return _p; }
01861 
01862 
01863   // evaluating the function itself
01864   RealType gamma ( const aol::Vec3<RealType> &arg ) const {
01865     // first apply B on the argument => rotation
01866     aol::Vec3<RealType> z;
01867     _B.mult ( arg, z );
01868 
01869     return pow ( pow ( aol::Abs ( z[0] ), _q ) + pow ( aol::Abs ( z[1] ), _q )
01870                  + pow ( aol::Abs ( z[2] ), _q ) + pow ( _epsilon, _q ),  1 / _q );
01871   }
01872 
01873 
01874   void gammaFirstDerivative ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
01875     gammaFirstDerivative ( z, v );
01876   }
01877 
01878   void gammaFirstDerivative ( const aol::Vec3<RealType> &arg, aol::Vec3<RealType> &v ) const {
01879     // first apply B on the argument => rotation
01880     aol::Vec3<RealType> z;
01881     aol::Vec3<RealType> temp;
01882     _B.mult ( arg, z );
01883 
01884     aol::Vec3<RealType> Z ( aol::Abs ( z[0] ), aol::Abs ( z[1] ), aol::Abs ( z[2] ) );
01885     RealType a = pow ( Z[0], _q ) + pow ( Z[1], _q ) + pow ( Z[2], _q ) + pow ( _epsilon, _q );
01886 
01887     // It has to be q>1, because for z[i]=0 there appears a singularity
01888     temp[0] = pow ( Z[0], _q - 1. ) * aol::signum1at0 ( z[0] );
01889     temp[1] = pow ( Z[1], _q - 1. ) * aol::signum1at0 ( z[1] );
01890     temp[2] = pow ( Z[2], _q - 1. ) * aol::signum1at0 ( z[2] );
01891     temp *= pow ( a, _p / _q - 1. ) *  pow ( pow ( a, _p / _q ) + pow ( _epsilon, _p ), ( 1. - _p ) / _p );
01892 
01893     // due to the chain rule we have to multiply the rotation-matrix once again
01894 //     cerr << "gamma_z( " << z << ") = " << temp << endl;
01895     _BT.mult ( temp, v );
01896   }
01897 
01898   RealType gammaNorm ( const aol::Vec3<RealType> &arg ) const {
01899     // first apply B on the argument => rotation
01900     aol::Vec3<RealType> z;
01901     _B.mult ( arg, z );
01902 
01903     return pow ( pow ( aol::Abs ( z[0] ), _q ) + pow ( aol::Abs ( z[1] ), _q )
01904                  + pow ( aol::Abs ( z[2] ), _q ) + pow ( _epsilon, _q ), 1. / _q );
01905   }
01906 };
01907 
01908 
01909 
01910 
01911 template <typename RealType>
01912 class RegMaxAnisotropy {
01913 private:
01914   RealType _epsilon;
01915 public:
01916   RegMaxAnisotropy ( RealType eps ) : _epsilon ( eps ) { }
01917 
01918   RealType p() const { return this->_p; }
01919 
01920   void implicitPart ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
01921     implicitPart ( z, mat );
01922   }
01923 
01924   void explicitPart ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
01925     explicitPart ( z, mat );
01926   }
01927 
01928   void implicitPart ( const aol::Vec2<RealType> &/*z*/, aol::Matrix22<RealType> &/*mat*/ ) const {
01929     cerr << "implicit part";
01930   }
01931 
01932   void explicitPart ( const aol::Vec2<RealType> &/*z*/, aol::Matrix22<RealType> &/*mat*/ ) const {
01933     cerr << "explicit part";
01934   }
01935 
01936 
01937   void gammaFirstDerivative ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
01938     gammaFirstDerivative ( z, v );
01939   }
01940 
01941   void gammaFirstDerivative ( const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
01942     RealType denom = ( 1. + aol::Sqr ( _epsilon ) ) * ( aol::Sqr ( z[0] ) + aol::Sqr ( z[1] ) );
01943     v[0] = ( ( 1. + aol::Sqr ( _epsilon ) ) * z[0] - z[1] ) / sqrt ( denom - 2. * z[0] * z[1] )
01944            + ( ( 1. + aol::Sqr ( _epsilon ) ) * z[0] + z[1] ) / sqrt ( denom + 2. * z[0] * z[1] );
01945 
01946     v[1] = ( ( 1. + aol::Sqr ( _epsilon ) ) * z[1] - z[0] ) / sqrt ( denom - 2. * z[0] * z[1] )
01947            + ( ( 1. + aol::Sqr ( _epsilon ) ) * z[1] + z[0] ) / sqrt ( denom + 2. * z[0] * z[1] );
01948     //cerr << "gamma_z( " << z << ") = " << v << endl;
01949   }
01950 
01951   RealType gammaNorm ( const aol::Vec2<RealType> &z ) const {
01952     return sqrt ( aol::Sqr ( z[0] - z[1] ) + aol::Sqr ( _epsilon ) * z.normSqr() ) +
01953            sqrt ( aol::Sqr ( z[0] + z[1] ) + aol::Sqr ( _epsilon ) * z.normSqr() );
01954 
01955   }
01956 };
01957 
01958 
01959 template <typename RealType>
01960 class Isotropy: public Anisotropy3dGraphInterface<RealType, Isotropy<RealType> > {
01961 public:
01962   using Anisotropy3dGraphInterface<RealType, Isotropy<RealType> >::gammaNorm;
01963   using Anisotropy3dGraphInterface<RealType, Isotropy<RealType> >::gammaFirstDerivative; 
01964   using Anisotropy3dGraphInterface<RealType, Isotropy<RealType> >::implicitPart;  
01965   
01966 private:
01967   RealType _epsilon;
01968   
01969 public:
01970   Isotropy ( RealType eps ) : Anisotropy3dGraphInterface<RealType, Isotropy<RealType> > ( ), _epsilon ( eps ) { }
01971 
01972   void implicitPart ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
01973     implicitPart ( z, mat );
01974   }
01975 
01976   void explicitPart ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
01977     explicitPart ( z, mat );
01978   }
01979 
01980   void implicitPart ( const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
01981     mat[0][0] = mat[1][1] = 1.;
01982     mat[1][0] = mat[0][1] = 0.;
01983     mat /= sqrt ( aol::Sqr ( z[0] ) + aol::Sqr ( z[1] ) + aol::Sqr ( _epsilon ) ) ;
01984   }
01985 
01986   void explicitPart ( const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
01987     for ( int i = 0; i < 2; i++ ) for ( int j = 0; j < 2; j++ ) mat[i][j] = z[i] * z[j];
01988     mat /= pow ( sqrt ( aol::Sqr ( z[0] ) + aol::Sqr ( z[1] ) + aol::Sqr ( _epsilon ) ), 3. ) ;
01989   }
01990 
01991 
01992   void gammaFirstDerivative ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
01993     gammaFirstDerivative ( z, v );
01994   }
01995 
01996   void gammaFirstDerivative ( const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
01997     v = z;
01998     v /= gammaNorm ( z );
01999   }
02000 
02001   RealType gammaNorm ( const aol::Vec2<RealType> &z ) const {
02002     return sqrt ( aol::Sqr ( z[0] ) + aol::Sqr ( z[1] ) + aol::Sqr ( _epsilon ) ) ;
02003   }
02004 };
02005 
02006 
02007 template <typename RealType>
02008 class Isotropy3d : public Anisotropy3dGraphInterface<RealType, Isotropy3d<RealType> >{
02009 public:  
02010   using Anisotropy3dGraphInterface<RealType, Isotropy3d<RealType> >::gammaNorm;
02011   using Anisotropy3dGraphInterface<RealType, Isotropy3d<RealType> >::gammaFirstDerivative;
02012   using Anisotropy3dGraphInterface<RealType, Isotropy3d<RealType> >::implicitPart;
02013   using Anisotropy3dGraphInterface<RealType, Isotropy3d<RealType> >::explicitPart;
02014   
02015 private:
02016   RealType _epsilon;
02017 public:
02018   Isotropy3d ( RealType eps ) : 
02019     Anisotropy3dGraphInterface<RealType, Isotropy3d<RealType> > ( ), _epsilon ( eps ) { }
02020 
02021   void implicitPart ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Matrix33<RealType> &mat ) const {
02022     implicitPart ( z, mat );
02023   }
02024 
02025   void explicitPart ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Matrix33<RealType> &mat ) const {
02026     explicitPart ( z, mat );
02027   }
02028 
02029   void implicitPart ( const aol::Vec3<RealType> &z, aol::Matrix33<RealType> &mat ) const {
02030     mat.setIdentity();
02031     mat /= sqrt ( aol::Sqr ( z[0] ) + aol::Sqr ( z[1] ) + aol::Sqr ( z[2] ) + aol::Sqr ( _epsilon ) ) ;
02032   }
02033 
02034   void explicitPart ( const aol::Vec3<RealType> &z, aol::Matrix33<RealType> &mat ) const {
02035     for ( int i = 0; i < 3; i++ ) for ( int j = 0; j < 3; j++ ) mat[i][j] = z[i] * z[j];
02036     mat /= pow ( sqrt ( aol::Sqr ( z[0] ) + aol::Sqr ( z[1] ) + aol::Sqr ( z[2] ) + aol::Sqr ( _epsilon ) ), 3. ) ;
02037   }
02038 
02039 
02040   void gammaFirstDerivative ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
02041     gammaFirstDerivative ( z, v );
02042   }
02043 
02044   void gammaFirstDerivative ( const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
02045     v = z;
02046     v /= gammaNorm ( z );
02047   }
02048 
02049   RealType gammaNorm ( const aol::Vec3<RealType> &z ) const {
02050     return sqrt ( aol::Sqr ( z[0] ) + aol::Sqr ( z[1] ) + aol::Sqr ( z[2] ) + aol::Sqr ( _epsilon ) ) ;
02051   }
02052 };
02053 
02054 
02055 template <typename RealType>
02056 class EllipseAnisotropy {
02057 private:
02058   RealType _a, _b, _epsilon;
02059 public:
02060   EllipseAnisotropy ( RealType a, RealType b, RealType eps ) : _a ( a ), _b ( b ), _epsilon ( eps ) { }
02061 
02062 
02063   void implicitPart ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
02064     implicitPart ( z, mat );
02065   }
02066 
02067   void explicitPart ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
02068     explicitPart ( z, mat );
02069   }
02070 
02071   void gammaFirstDerivative ( const aol::Vec2<RealType> &, const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
02072     gammaFirstDerivative ( z, v );
02073   }
02074 
02075 
02076   void implicitPart ( const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
02077     mat[0][0] = aol::Sqr ( _a );
02078     mat[1][1] = aol::Sqr ( _b );
02079     mat[0][1] = mat[1][0] = 0.;
02080     RealType normSqr = aol::Sqr ( z[0] * _a ) + aol::Sqr ( z[1] * _b ) + _epsilon * _epsilon;
02081     mat /= sqrt ( normSqr );
02082   }
02083 
02084   void explicitPart ( const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
02085     aol::Vec2<RealType> v;
02086     gammaFirstDerivative ( z, v );
02087     mat[0][0] = v[0] * v[0];
02088     mat[1][1] = v[1] * v[1];
02089     mat[0][1] =  mat[1][0] = v[1] * v[0];
02090     RealType normSqr = aol::Sqr ( z[0] * _a ) + aol::Sqr ( z[1] * _b ) + _epsilon * _epsilon;
02091     mat /= sqrt ( normSqr );
02092   }
02093 
02094   void gammaFirstDerivative ( const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
02095     v[0] = aol::Sqr ( _a ) * z[0];
02096     v[1] = aol::Sqr ( _b ) * z[1];
02097     RealType normSqr = aol::Sqr ( z[0] * _a ) + aol::Sqr ( z[1] * _b ) + _epsilon * _epsilon;
02098     v /= sqrt ( normSqr );
02099   }
02100 
02101   RealType gammaNorm ( const aol::Vec2<RealType> &z ) const {
02102     return sqrt ( aol::Sqr ( z[0]*_a ) + aol::Sqr ( z[1] * _b ) + _epsilon*_epsilon );
02103   }
02104 
02105 };
02106 
02107 
02108 // ------------------------ ELLIPSOID ----------------------------------------------
02109 
02110 template <typename RealType>
02111 class EllipsoidAnisotropy : public Anisotropy3dGraphInterface<RealType, EllipsoidAnisotropy<RealType> > {
02112 public:
02113   using Anisotropy3dGraphInterface<RealType, EllipsoidAnisotropy<RealType> >::gammaNorm;
02114   using Anisotropy3dGraphInterface<RealType, EllipsoidAnisotropy<RealType> >::gammaFirstDerivative;
02115   
02116 private:
02117   RealType _a, _b, _c, _epsilon;
02118 public:
02119   EllipsoidAnisotropy ( RealType a, RealType b, RealType c, RealType eps ) :
02120     Anisotropy3dGraphInterface<RealType, EllipsoidAnisotropy<RealType> > ( ),
02121       _a ( a ), _b ( b ), _c ( c ), _epsilon ( eps ) { }
02122 
02123   // evaluating the function itself
02124   RealType gamma ( const aol::Vec3<RealType> &z ) const {
02125     return sqrt ( aol::Sqr ( z[0]*_a ) + aol::Sqr ( z[1]*_b )
02126                   + aol::Sqr ( z[2]*_c ) + _epsilon*_epsilon );
02127   }
02128 
02129   void gammaFirstDerivative ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
02130     gammaFirstDerivative ( z, v );
02131   }
02132 
02133   // evaluating the first derivative
02134   void gammaFirstDerivative ( const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
02135     v[0] = aol::Sqr ( _a ) * z[0];
02136     v[1] = aol::Sqr ( _b ) * z[1];
02137     v[2] = aol::Sqr ( _c ) * z[2];
02138     RealType normSqr = aol::Sqr ( z[0] * _a ) + aol::Sqr ( z[1] * _b )
02139                        + aol::Sqr ( z[2] * _c ) + _epsilon * _epsilon;
02140     v /= sqrt ( normSqr );
02141   }
02142 
02143   RealType gammaNorm ( const aol::Vec3<RealType> &z ) const {
02144     return sqrt ( aol::Sqr ( z[0]*_a ) + aol::Sqr ( z[1] * _b )
02145                   + aol::Sqr ( z[2] * _c ) + _epsilon*_epsilon );
02146   }
02147 
02148 };
02149 
02150 // -----------------------------------------------------------------------------------------------
02151 
02152 
02153 // ------------------------ ROTIERTER ELLIPSOID ----------------------------------------------
02154 // diese Klasse stellt einen Ellipsoid zur Verfuegung, der so rotiert ist,
02155 // dass die x-Achse auf einem gegebenen Vektor v liegt. Die anderen beiden
02156 // Koordinatenachsen sind aufgrund der Symmetrie des Ellipsoids willkuerlich
02157 // gewaehlt.
02158 // -------------------------------------------------------------------------------------------
02159 
02160 template <typename RealType>
02161 class RotatedEllipsoidAnisotropy {
02162 private:
02163   RealType _a, _b, _c, _epsilon;
02164   aol::Vec3<RealType> _v;        // direction, on which the x-axxis lies
02165   aol::Matrix33<RealType> _B;    // matrix to perform a basis-transformation
02166   aol::Matrix33<RealType> _BT;   // the transposed transformation-matrix
02167   bool _isotrop;                // if v is 0, then calc like isotropic mcm => flag
02168 public:
02169   RotatedEllipsoidAnisotropy ( RealType a, RealType b, RealType c, RealType eps ) :
02170       _a ( a ), _b ( b ), _c ( c ), _epsilon ( eps ), _v ( 1., 0., 0. ) {
02171     calcB();
02172     _isotrop = false;
02173   }
02174 
02175   RotatedEllipsoidAnisotropy ( RealType a, RealType b, RealType c, RealType eps, aol::Vec3<RealType> v ) :
02176       _a ( a ), _b ( b ), _c ( c ), _epsilon ( eps ),  _v ( v ) {
02177     calcB();
02178     _isotrop = false;
02179   }
02180 
02181 
02182   // calc the matrix B for transforming the coordinate-system for the x-axxis
02183   // lying on v
02184   void calcB() {
02185     // the first row is the direction itself
02186     _B.setRow ( 0, _v );
02187     // the other directions must be orthogonal to _v, but then, they can be
02188     // chosen arbitrary, first choose y orthogonal to _v
02189     // TODO: arranging by size to avoid taking the smallest values (=> stability)
02190     aol::Vec3<RealType> y ( 0, 0, 0 );
02191     if ( _v[1] != 0 ) {
02192       if ( _v[2] != 0 ) {
02193         y[1] = -_v[2];
02194         y[2] = _v[1];
02195         y /= y.norm();
02196       } else y[2] = 1;
02197     } else y[1] = 1;
02198 
02199     _B.setRow ( 1, y );
02200 
02201     // now z = vector-product of v and y, both normed => z is normed too
02202     aol::Vec3<RealType> z ( _v.crossProduct ( y ) );
02203     _B.setRow ( 2, z );
02204 
02205     // finally save the transposed matrix
02206     _BT = _B;
02207     _BT.transpose();
02208   }
02209 
02210 
02211   // set the rotate-direction v
02212   void setRotateDirection ( const aol::Vec3<RealType> &v ) {
02213     // only anisotropic, if v is not 0
02214     // if v is not 0, it is very improbable, that v is equal to _v
02215     // (that means = v from the timestep before)
02216     // => it's not worthwile to check this
02217     if ( v[0] != 0 || v[1] != 0 || v[2] != 0 ) {
02218       _v = v / v.norm();
02219       calcB();
02220       _isotrop = false;
02221     } else _isotrop = true;
02222   }
02223 
02224   // set the rotate-direction v
02225   void setRotateDirection ( RealType v0, RealType v1, RealType v2 ) {
02226     if ( v0 != 0 || v1 != 0 || v2 != 0 ) {
02227       _v[0] = v0; _v[1] = v1; _v[2] = v2;
02228       _v /= _v.norm();
02229       calcB();
02230       _isotrop = false;
02231     } else _isotrop = true;
02232   }
02233 
02234 
02235 
02236   // evaluating the function itself
02237   RealType gamma ( const aol::Vec3<RealType> &arg ) const {
02238     if ( !_isotrop ) {
02239       // first apply B on the argument => rotation
02240       aol::Vec3<RealType> _z;
02241       _B.mult ( arg, _z );
02242 
02243       return sqrt ( aol::Sqr ( _z[0]*_a ) + aol::Sqr ( _z[1]*_b )
02244                     + aol::Sqr ( _z[2]*_c ) + _epsilon*_epsilon );
02245     } else return arg.norm();
02246 //     1.;    // |z|
02247   }
02248 
02249   void gammaFirstDerivative ( const aol::Vec3<RealType> &, const aol::Vec3<RealType> &z, aol::Vec3<RealType> &v ) const {
02250     gammaFirstDerivative ( z, v );
02251   }
02252 
02253   // evaluating the first derivative
02254   void gammaFirstDerivative ( const aol::Vec3<RealType> &arg, aol::Vec3<RealType> &v ) const {
02255     if ( !_isotrop ) {
02256       aol::Vec3<RealType> _z, temp;
02257       _B.mult ( arg, _z );
02258       temp[0] = aol::Sqr ( _a ) * _z[0];
02259       temp[1] = aol::Sqr ( _b ) * _z[1];
02260       temp[2] = aol::Sqr ( _c ) * _z[2];
02261       RealType normSqr = aol::Sqr ( _z[0] * _a ) + aol::Sqr ( _z[1] * _b )
02262                          + aol::Sqr ( _z[2] * _c ) + _epsilon * _epsilon;
02263       temp /= sqrt ( normSqr );
02264 
02265       // due to the chain rule we have to multiply the rotation-matrix once again
02266       _BT.mult ( temp, v );
02267     } else {
02268       v = arg;
02269       v /= arg.norm();
02270     }
02271   }
02272 
02273   RealType gammaNorm ( const aol::Vec3<RealType> &z ) const {
02274     aol::Vec3<RealType> _z;
02275     _B.mult ( z, _z );
02276     return sqrt ( aol::Sqr ( _z[0]*_a ) + aol::Sqr ( _z[1] * _b )
02277                   + aol::Sqr ( _z[2] * _c ) + _epsilon*_epsilon );
02278   }
02279 
02280 };
02281 
02282 // -----------------------------------------------------------------------------------------------
02283 
02284 
02285 
02286 
02287 template <typename RealType>
02288 class MixedAnisotropy {
02289   LpAnisotropy<RealType>      _a11, _a12;
02290   EllipseAnisotropy<RealType> _a21, _a22;
02291 public:
02292   MixedAnisotropy ( RealType p, RealType a, RealType b, RealType epsilon ) :
02293       _a11 ( 3., p, epsilon ), _a12 ( p, p, epsilon ),
02294       _a21 ( a, b, epsilon ), _a22 ( 1.0, 1.0, epsilon )  {}
02295 
02296   void implicitPart ( const aol::Vec2<RealType> &pt, const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
02297     _a12.implicitPart ( z, mat );
02298     return;
02299 
02300     if ( pt[0] < 0.5 ) {
02301       if ( pt[1] < 0.5 ) {
02302         _a11.implicitPart ( z, mat );
02303       } else {
02304         _a12.implicitPart ( z, mat );
02305       }
02306     } else {
02307       if ( pt[1] < 0.5 ) {
02308         _a21.implicitPart ( z, mat );
02309       } else {
02310         _a22.implicitPart ( z, mat );
02311       }
02312     }
02313     return;
02314 
02315     aol::Matrix22<RealType> mtmp, munten, moben;
02316 
02317     if ( pt[0] < 0.45 ) {
02318       _a11.implicitPart ( z, munten );
02319     } else if ( pt[0] > 0.55 ) {
02320       _a12.implicitPart ( z, munten );
02321     } else {
02322       RealType lambda = ( pt[0] - 0.45 ) / 0.1;
02323       _a11.implicitPart ( z, mtmp );
02324       mtmp *= ( 1. - lambda );
02325       _a12.implicitPart ( z, munten );
02326       munten *= lambda;
02327       munten += mtmp;
02328     }
02329     if ( pt[0] < 0.45 ) {
02330       _a21.implicitPart ( z, moben );
02331     } else if ( pt[0] > 0.55 ) {
02332       _a22.implicitPart ( z, moben );
02333     } else {
02334       RealType lambda = ( pt[0] - 0.45 ) / 0.1;
02335       _a21.implicitPart ( z, mtmp );
02336       mtmp *= ( 1. - lambda );
02337       _a22.implicitPart ( z, moben );
02338       moben *= lambda;
02339       moben += mtmp;
02340     }
02341     if ( pt[1] < 0.45 ) {
02342       mat = munten;
02343       return;
02344     }
02345     if ( pt[1] > 0.55 ) {
02346       mat = moben;
02347       return;
02348     }
02349     RealType lambda = ( pt[1] - 0.45 ) / 0.1;
02350     munten *= ( 1. - lambda );
02351     moben *= lambda;
02352     mat = munten;
02353     mat += moben;
02354   }
02355 
02356   void explicitPart ( const aol::Vec2<RealType> &pt, const aol::Vec2<RealType> &z, aol::Matrix22<RealType> &mat ) const {
02357     _a12.explicitPart ( z, mat );
02358     return;
02359 
02360     if ( pt[0] < 0.5 ) {
02361       if ( pt[1] < 0.5 ) {
02362         _a11.explicitPart ( z, mat );
02363       } else {
02364         _a12.explicitPart ( z, mat );
02365       }
02366     } else {
02367       if ( pt[1] < 0.5 ) {
02368         _a21.explicitPart ( z, mat );
02369       } else {
02370         _a22.explicitPart ( z, mat );
02371       }
02372     }
02373     return;
02374 
02375     aol::Matrix22<RealType> mtmp, munten, moben;
02376 
02377     if ( pt[0] < 0.45 ) {
02378       _a11.explicitPart ( z, munten );
02379     } else if ( pt[0] > 0.55 ) {
02380       _a12.explicitPart ( z, munten );
02381     } else {
02382       RealType lambda = ( pt[0] - 0.45 ) / 0.1;
02383       _a11.explicitPart ( z, mtmp );
02384       mtmp *= ( 1. - lambda );
02385       _a12.explicitPart ( z, munten );
02386       munten *= lambda;
02387       munten += mtmp;
02388     }
02389     if ( pt[0] < 0.45 ) {
02390       _a21.explicitPart ( z, moben );
02391     } else if ( pt[0] > 0.55 ) {
02392       _a22.explicitPart ( z, moben );
02393     } else {
02394       RealType lambda = ( pt[0] - 0.45 ) / 0.1;
02395       _a21.explicitPart ( z, mtmp );
02396       mtmp *= ( 1. - lambda );
02397       _a22.explicitPart ( z, moben );
02398       moben *= lambda;
02399       moben += mtmp;
02400     }
02401     if ( pt[1] < 0.45 ) {
02402       mat = munten;
02403       return;
02404     }
02405     if ( pt[1] > 0.55 ) {
02406       mat = moben;
02407       return;
02408     }
02409     RealType lambda = ( pt[1] - 0.45 ) / 0.1;
02410     munten *= ( 1. - lambda );
02411     moben *= lambda;
02412     mat = munten;
02413     mat += moben;
02414   }
02415 
02416   void gammaFirstDerivative ( const aol::Vec2<RealType> &pt, const aol::Vec2<RealType> &z, aol::Vec2<RealType> &v ) const {
02417     _a12.gammaFirstDerivative ( z, v );
02418     return;
02419 
02420     if ( pt[0] < 0.5 ) {
02421       if ( pt[1] < 0.5 ) {
02422         _a11.gammaFirstDerivative ( z, v );
02423       } else {
02424         _a12.gammaFirstDerivative ( z, v );
02425       }
02426     } else {
02427       if ( pt[1] < 0.5 ) {
02428         _a21.gammaFirstDerivative ( z, v );
02429       } else {
02430         _a22.gammaFirstDerivative ( z, v );
02431       }
02432     }
02433     return;
02434 
02435     aol::Vec2<RealType> mtmp, munten, moben;
02436 
02437     if ( pt[0] < 0.45 ) {
02438       _a11.gammaFirstDerivative ( z, munten );
02439     } else if ( pt[0] > 0.55 ) {
02440       _a12.gammaFirstDerivative ( z, munten );
02441     } else {
02442       RealType lambda = ( pt[0] - 0.45 ) / 0.1;
02443       _a11.gammaFirstDerivative ( z, mtmp );
02444       mtmp *= ( 1. - lambda );
02445       _a12.gammaFirstDerivative ( z, munten );
02446       munten *= lambda;
02447       munten += mtmp;
02448     }
02449     if ( pt[0] < 0.45 ) {
02450       _a21.gammaFirstDerivative ( z, moben );
02451     } else if ( pt[0] > 0.55 ) {
02452       _a22.gammaFirstDerivative ( z, moben );
02453     } else {
02454       RealType lambda = ( pt[0] - 0.45 ) / 0.1;
02455       _a21.gammaFirstDerivative ( z, mtmp );
02456       mtmp *= ( 1. - lambda );
02457       _a22.gammaFirstDerivative ( z, moben );
02458       moben *= lambda;
02459       moben += mtmp;
02460     }
02461     if ( pt[1] < 0.45 ) {
02462       v = munten;
02463       return;
02464     }
02465     if ( pt[1] > 0.55 ) {
02466       v = moben;
02467       return;
02468     }
02469     RealType lambda = ( pt[1] - 0.45 ) / 0.1;
02470     munten *= ( 1. - lambda );
02471     moben *= lambda;
02472     v = munten;
02473     v += moben;
02474   }
02475 
02476   RealType gammaNorm ( const aol::Vec2<RealType> &z ) const {
02477     return _a12.gammaNorm ( z );
02478   }
02479 
02480 };
02481 
02482 } // namespace qc
02483 
02484 #endif
02485 

Generated on Fri Sep 9 2011 21:08:09 for QuocMesh by doxygen 1.7.1