QuOc

 

hyperelastic.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 __HYPERELASTIC_H
00012 #define __HYPERELASTIC_H
00013 
00014 #include <op.h>
00015 #include <vec.h>
00016 #include <multiVector.h>
00017 #include <linearSmoothOp.h>
00018 #include <FEOpInterface.h>
00019 #include <deformations.h>
00020 
00021 namespace qc {
00022 
00023 template <typename ConfiguratorType, qc::Dimension Dim>
00024 class VolumeGradient
00025       : public aol::FENonlinVectorDiffOpInterface<ConfiguratorType, Dim, Dim, VolumeGradient<ConfiguratorType, Dim> > {
00026 public:
00027   typedef typename ConfiguratorType::RealType RealType;
00028   const qc::GridDefinition &_grid;
00029 protected:
00030 public:
00031   VolumeGradient ( const qc::GridDefinition &Grid )
00032       : aol::FENonlinVectorDiffOpInterface<ConfiguratorType, Dim, Dim, VolumeGradient<ConfiguratorType, Dim> > ( Grid ),
00033       _grid ( Grid ) {};
00034 
00035   RealType gammaPrime ( RealType det ) const {
00036     return ( 2. * det - 2. );
00037   }
00038 
00039   void getNonlinearity ( const aol::auto_container<Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &DiscFuncs,
00040                          const typename ConfiguratorType::ElementType &El,
00041                          int QuadPoint, const typename ConfiguratorType::VecType &/*RefCoord*/,
00042                          aol::Mat<Dim,ConfiguratorType::Dim,RealType> &NL ) const {
00043     typename ConfiguratorType::VecType grad;
00044     typename ConfiguratorType::MatType dphi;
00045 
00046     NL.setZero();
00047     dphi.setZero();
00048 
00049     for ( int i = 0; i < Dim; i++ ) {
00050       DiscFuncs[i].evaluateGradientAtQuadPoint ( El, QuadPoint, grad );
00051       for ( int j = 0; j < Dim; j++ ) {
00052         dphi[i][j] = grad[j];
00053       }
00054     }
00055 
00056     for ( int i = 0; i < Dim; i++ ) {
00057       dphi[i][i] += 1.;
00058     }
00059 
00060     typename ConfiguratorType::MatType aux;
00061     aux.makeCofactorMatrix ( dphi );
00062     NL = aux;
00063 
00064     RealType det = dphi.det();
00065     NL *= gammaPrime ( det );
00066   }
00067 };
00068 
00069 template <typename ConfiguratorType, qc::Dimension Dim>
00070 class VolumeEnergy
00071       : public aol::FENonlinIntegrationVectorInterface < ConfiguratorType,
00072       VolumeEnergy<ConfiguratorType, Dim> > {
00073 public:
00074   typedef typename ConfiguratorType::RealType RealType;
00075   const qc::GridDefinition &_grid;
00076 protected:
00077 public:
00078   VolumeEnergy ( const qc::GridDefinition &Grid )
00079       : aol::FENonlinIntegrationVectorInterface < ConfiguratorType,
00080       VolumeEnergy<ConfiguratorType, Dim> > ( Grid ),
00081       _grid ( Grid ) {};
00082 
00083   RealType evaluateIntegrand ( const aol::auto_container<ConfiguratorType::Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &discrFuncs,
00084                                const typename ConfiguratorType::ElementType &El,
00085                                int QuadPoint, const typename ConfiguratorType::VecType &/*RefCoord*/ ) const {
00086     typename ConfiguratorType::VecType grad;
00087     typename ConfiguratorType::MatType dphi;
00088 
00089     dphi.setZero();
00090 
00091     for ( int i = 0; i < Dim; i++ ) {
00092       discrFuncs[i].evaluateGradientAtQuadPoint ( El, QuadPoint, grad );
00093       for ( int j = 0; j < Dim; j++ ) {
00094         dphi[i][j] = grad[j];
00095       }
00096     }
00097 
00098     for ( int j = 0; j < Dim; j++ ) {
00099       dphi[j][j] += 1.;
00100     }
00101 
00102     RealType det = dphi.det();
00103     return gamma ( det );
00104   }
00105 protected:
00106   RealType gamma ( RealType det ) const {
00107     return aol::Sqr ( det - 1. );
00108   }
00109 };
00110 
00111 
00112 template <typename ConfiguratorType, qc::Dimension Dim>
00113 class LengthGradient
00114       : public aol::FENonlinVectorDiffOpInterface<ConfiguratorType, Dim, Dim, LengthGradient<ConfiguratorType, Dim> > {
00115 public:
00116   typedef typename ConfiguratorType::RealType RealType;
00117   const qc::GridDefinition &_grid;
00118 protected:
00119 public:
00120   LengthGradient ( const qc::GridDefinition &Grid )
00121       : aol::FENonlinVectorDiffOpInterface<ConfiguratorType, Dim, Dim, LengthGradient<ConfiguratorType, Dim> > ( Grid ),
00122       _grid ( Grid ) {};
00123 
00124   void getNonlinearity ( const aol::auto_container<Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &DiscFuncs,
00125                          const typename ConfiguratorType::ElementType &El,
00126                          int QuadPoint, const typename ConfiguratorType::VecType &/*RefCoord*/,
00127                          aol::Mat<Dim,ConfiguratorType::Dim,RealType> &NL ) const {
00128     typename ConfiguratorType::VecType grad;
00129 
00130     for ( int i = 0; i < Dim; i++ ) {
00131       DiscFuncs[i].evaluateGradientAtQuadPoint ( El, QuadPoint, grad );
00132       for ( int j = 0; j < Dim; j++ ) {
00133         NL[i][j] = grad[j];
00134       }
00135     }
00136   }
00137 };
00138 
00139 
00140 template <typename ConfiguratorType, qc::Dimension Dim>
00141 class LengthEnergy
00142       : public aol::FENonlinIntegrationVectorInterface < ConfiguratorType,
00143       LengthEnergy<ConfiguratorType, Dim> > {
00144 public:
00145   typedef typename ConfiguratorType::RealType RealType;
00146   const qc::GridDefinition &_grid;
00147 protected:
00148 public:
00149   LengthEnergy ( const qc::GridDefinition &Grid )
00150       : aol::FENonlinIntegrationVectorInterface < ConfiguratorType,
00151       LengthEnergy<ConfiguratorType, Dim> > ( Grid ),
00152       _grid ( Grid ) {};
00153 
00154   RealType evaluateIntegrand ( const aol::auto_container<ConfiguratorType::Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &discrFuncs,
00155                                const typename ConfiguratorType::ElementType &El,
00156                                int QuadPoint, const typename ConfiguratorType::VecType &/*RefCoord*/ ) const {
00157     typename ConfiguratorType::VecType grad;
00158 
00159     RealType energy = 0.;
00160 
00161     for ( int i = 0; i < Dim; i++ ) {
00162       discrFuncs[i].evaluateGradientAtQuadPoint ( El, QuadPoint, grad );
00163       energy += grad.normSqr();
00164     }
00165     return 0.5 * energy;
00166   }
00167 };
00168 
00169 /**
00170  * Does the same as LengthEnergy, but uses a StiffOp (has to be supplied as argument) to do so.
00171  * instead of using aol::FENonlinIntegrationVectorInterface. If the supplied StiffOp is assembled,
00172  * applying DisplacementLengthEnergy is much faster than applying LengthEnergy.
00173  *
00174  * \author Berkels
00175  */
00176 template <typename ConfiguratorType>
00177 class DisplacementLengthEnergy: public aol::Op<aol::MultiVector<typename ConfiguratorType::RealType>, aol::Scalar<typename ConfiguratorType::RealType> > 
00178 {
00179   typedef typename ConfiguratorType::RealType RealType;
00180   const aol::DiagonalBlockOp<RealType> _blockStiff;
00181 public:
00182 
00183   DisplacementLengthEnergy( const aol::StiffOp<ConfiguratorType> &Stiff )
00184     : _blockStiff ( Stiff ) {
00185   }
00186   virtual void apply( const aol::MultiVector<RealType> &Arg, aol::Scalar<RealType> &Dest ) const {
00187     aol::MultiVector <RealType> mtmp( Arg, aol::STRUCT_COPY );
00188     _blockStiff.apply( Arg, mtmp); // E = 1/2 L_dPhi*Phi
00189     Dest[0] = (0.5)*(Arg * mtmp);
00190   }
00191   virtual void applyAdd( const aol::MultiVector<RealType> &Arg, aol::Scalar<RealType> &Dest) const {
00192     aol::Scalar<RealType> tmp ( Dest );
00193     apply ( Arg, tmp );
00194     Dest += tmp;
00195   }
00196 };
00197 
00198 /**
00199   * \brief This class computes via "apply(...)" or "applyAdd(...)" a hyperelastic energy \f$ \int_\Omega W(I_1,I_2,I_3,x) dx \f$,
00200   * where the hyperelastic invariants are given by \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2 \f$, 
00201   * \f$ I_1=det(\nabla\phi) \f$ and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
00202   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$
00203   * in the different Cartesian directions; the deformation \f$ \phi \f$ is then automatically computed via \f$ \phi=d+\f$ identity.
00204   *
00205   * \author Wirth
00206   */
00207 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
00208 class HyperelasticEnergy
00209       : public aol::FENonlinIntegrationVectorInterface < ConfiguratorType, HyperelasticEnergy<ConfiguratorType,HyperelasticEnergyDensityType> > {
00210 public:
00211   typedef typename ConfiguratorType::RealType RealType;
00212 protected:
00213   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
00214   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
00215 public:
00216   HyperelasticEnergy ( const typename ConfiguratorType::InitType &Grid, const HyperelasticEnergyDensityType &HyperelasticEnergyDensity )
00217       : aol::FENonlinIntegrationVectorInterface < ConfiguratorType, HyperelasticEnergy<ConfiguratorType,HyperelasticEnergyDensityType> > ( Grid ),
00218         _hyperelasticEnergyDensity( HyperelasticEnergyDensity ) {}
00219 
00220   /**
00221     * \brief Returns the hyperelastic energy \f$ W(I_1,I_2,I_3) \f$ for the invariants \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2\f$, \f$ I_3=det(\nabla\phi) \f$.
00222     */
00223   RealType evaluateIntegrand ( const aol::auto_container<ConfiguratorType::Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &DiscFuncs,
00224                                const typename ConfiguratorType::ElementType &El, int QuadPoint,
00225                                const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
00226     // compute the deformation gradient $\nabla\phi$: for each displacment component do...
00227     typename ConfiguratorType::MatType dphi;
00228     for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00229       // evaluate gradient of the displacement $d_i$
00230       DiscFuncs[i].evaluateGradientAtQuadPoint ( El, QuadPoint, dphi[i] );
00231       // compute corresponding deformation gradient of the $i$th component
00232       dphi[i][i] += 1.;
00233     }
00234     return energyDensity( dphi, _hyperelasticEnergyDensity, El, QuadPoint );
00235   }
00236 
00237   /**
00238     * \brief Returns the hyperelastic energy density \f$ W(I_1,I_2,I_3) \f$ for the invariants
00239     * \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2\f$, \f$ I_3=det(\nabla\phi) \f$,
00240     * where \f$ \nabla\phi \f$, \f$ W \f$, and the spatial position are passed as arguments.
00241     */
00242   static inline RealType energyDensity ( const typename ConfiguratorType::MatType &DPhi,
00243                                          const HyperelasticEnergyDensityType &EnergyDensity,
00244                                          const typename ConfiguratorType::ElementType &El, int QuadPoint ) {
00245     // compute the hyperelastic invariants $||\nabla\phi||_F^2,||cof(\nabla\phi)||_F^2,det(\nabla\phi)$
00246     RealType
00247       I1 = DPhi.normSqr(),
00248       I2,
00249       I3 = DPhi.det();
00250     if ( ConfiguratorType::Dim != qc::QC_2D ) {
00251       typename ConfiguratorType::MatType cof;
00252       cof.makeCofactorMatrix( DPhi );
00253       I2 = cof.normSqr();
00254     } else
00255       // in 2d, $||cof(\nabla\phi)||_F$ equals $||\nabla\phi||_F$
00256       I2 = I1;
00257     return EnergyDensity.evaluateAtQuadPoint ( I1, I2, I3, El, QuadPoint );
00258   }
00259 };
00260 
00261 /**
00262   * \brief This class computes via "apply(...)" or "applyAdd(...)" a hyperelastic energy gradient with respect to the displacement,
00263   * i.e. \f$ (<\frac{\partial \int_\Omega W(I_1,I_2,I_3,x) dx}{\partial \phi},\psi_i>)_i \f$, where the \f$ \psi_i \f$ are the FE basis
00264   * functions, the hyperelastic invariants are given by \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2 \f$,
00265   * \f$ I_3=det(\nabla\phi)\f$ and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
00266   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$
00267   * in the different Cartesian directions; the deformation \f$ \phi \f$ is then automatically computed via \f$ \phi=d+\f$ identity.
00268   *
00269   * \author Wirth
00270   */
00271 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
00272 class HyperelasticGradient
00273       : public aol::FENonlinVectorDiffOpInterface<ConfiguratorType, ConfiguratorType::Dim, ConfiguratorType::Dim, HyperelasticGradient<ConfiguratorType,HyperelasticEnergyDensityType> > {
00274 public:
00275   typedef typename ConfiguratorType::RealType RealType;
00276 protected:
00277   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
00278   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
00279 public:
00280   HyperelasticGradient ( const typename ConfiguratorType::InitType &Grid, const HyperelasticEnergyDensityType &HyperelasticEnergyDensity )
00281       : aol::FENonlinVectorDiffOpInterface<ConfiguratorType, ConfiguratorType::Dim, ConfiguratorType::Dim, HyperelasticGradient<ConfiguratorType,HyperelasticEnergyDensityType> > ( Grid ),
00282         _hyperelasticEnergyDensity( HyperelasticEnergyDensity ) {}
00283 
00284   /**
00285     * \brief Returns \f$ 2\frac{\partial W}{\partial I_1}\nabla\phi+2\frac{\partial W}{\partial I_2}\frac{\partial cof(\nabla\phi)}{\partial \phi}+\frac{\partial W}{\partial I_3}cof(\nabla\phi) \f$,
00286     * where \f$ \frac{\partial cof(\nabla\phi)}{\partial \phi} \f$ shall be short for the matrix \f$ M \f$ such that
00287     * \f$ tr(M^T\nabla\psi) = tr[cof(\nabla\phi)cof(\nabla\phi)^T(tr((\nabla\phi)^{-1}\nabla\psi)-\nabla\psi(\nabla\phi)^{-1})] \f$
00288     * for any \f$ \psi \f$.
00289     */
00290   void getNonlinearity ( const aol::auto_container<ConfiguratorType::Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &DiscFuncs,
00291                          const typename ConfiguratorType::ElementType &El,
00292                          int QuadPoint, const typename ConfiguratorType::DomVecType &/*RefCoord*/,
00293                          aol::Mat<ConfiguratorType::Dim,ConfiguratorType::Dim,RealType> &NL ) const {
00294     // compute the deformation gradient $\nabla\phi$: for each displacment component do...
00295     typename ConfiguratorType::MatType dPhi;
00296     for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00297       // evaluate gradient of the displacement $d_i$
00298       DiscFuncs[i].evaluateGradientAtQuadPoint ( El, QuadPoint, dPhi[i] );
00299       // compute corresponding deformation gradient of the $i$th component
00300       dPhi[i][i] += 1.;
00301     }
00302 
00303     firstPiolaKirchhoffStress( dPhi, _hyperelasticEnergyDensity, El, QuadPoint, NL );
00304   }
00305 
00306   /**
00307     * \brief Computes Stress=\f$ \frac{\partial W}{\partial\nabla\phi}(I_1,I_2,I_3) \f$ for the invariants
00308     * \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2\f$, \f$ I_3=det(\nabla\phi) \f$,
00309     * where \f$ \nabla\phi \f$, \f$ W \f$, and the spatial position are passed as arguments.
00310     */
00311   static inline void firstPiolaKirchhoffStress ( const typename ConfiguratorType::MatType &DPhi,
00312                                                  const HyperelasticEnergyDensityType &EnergyDensity,
00313                                                  const typename ConfiguratorType::ElementType &El, int QuadPoint,
00314                                                  aol::Mat<ConfiguratorType::Dim,ConfiguratorType::Dim,RealType> &Stress ) {
00315     Stress.setZero();
00316 
00317     // compute the cofactor matrix of the deformation gradient, $cof(\nabla\phi)$
00318     typename ConfiguratorType::MatType dPhi( DPhi ), cof;
00319     cof.makeCofactorMatrix( DPhi );
00320     // compute the hyperelastic invariants $||\nabla\phi||_F^2,||cof(\nabla\phi)||_F^2,det(\nabla\phi)$
00321     RealType
00322       I1 = DPhi.normSqr(),
00323       I2 = cof.normSqr(),
00324       I3 = DPhi.det();
00325 
00326     // compute the outer derivative of the hyperelastic energy
00327     aol::Vec3<RealType> outerDeriv ( EnergyDensity.evaluateDerivativeAtQuadPoint ( I1, I2, I3, El, QuadPoint ) );
00328 
00329     // compute the first Piola-Kirchhoff stress
00330     // the term belonging to $I_2$...
00331     if ( ConfiguratorType::Dim == qc::QC_3D ) {
00332       Stress[0][0] = DPhi[1][1] * cof[2][2] + DPhi[2][2] * cof[1][1] - DPhi[1][2] * cof[2][1] - DPhi[2][1] * cof[1][2];
00333       Stress[0][1] = DPhi[1][2] * cof[2][0] + DPhi[2][0] * cof[1][2] - DPhi[1][0] * cof[2][2] - DPhi[2][2] * cof[1][0];
00334       Stress[0][2] = DPhi[2][1] * cof[1][0] + DPhi[1][0] * cof[2][1] - DPhi[2][0] * cof[1][1] - DPhi[1][1] * cof[2][0];
00335       Stress[1][0] = DPhi[0][2] * cof[2][1] + DPhi[2][1] * cof[0][2] - DPhi[0][1] * cof[2][2] - DPhi[2][2] * cof[0][1];
00336       Stress[1][1] = DPhi[0][0] * cof[2][2] + DPhi[2][2] * cof[0][0] - DPhi[0][2] * cof[2][0] - DPhi[2][0] * cof[0][2];
00337       Stress[1][2] = DPhi[2][0] * cof[0][1] + DPhi[0][1] * cof[2][0] - DPhi[2][1] * cof[0][0] - DPhi[0][0] * cof[2][1];
00338       Stress[2][0] = DPhi[0][1] * cof[1][2] + DPhi[1][2] * cof[0][1] - DPhi[0][2] * cof[1][1] - DPhi[1][1] * cof[0][2];
00339       Stress[2][1] = DPhi[1][0] * cof[0][2] + DPhi[0][2] * cof[1][0] - DPhi[1][2] * cof[0][0] - DPhi[0][0] * cof[1][2];
00340       Stress[2][2] = DPhi[0][0] * cof[1][1] + DPhi[1][1] * cof[0][0] - DPhi[0][1] * cof[1][0] - DPhi[1][0] * cof[0][1];
00341       Stress *= outerDeriv[1] * 2;
00342       // alternative implementation
00343       /*if ( I3 != 0. ) {
00344         Stress.makeProductABtransposed( cof, cof );
00345         RealType trace = Stress.tr();
00346         for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00347           Stress[i][i] -= trace;
00348         Stress *= cof;
00349         Stress *= - outerDeriv[1] * 2 / I3;
00350       }*/
00351     } else
00352       // in 2d, the term belonging to $I_2$ equals the one belonging to $I_1$
00353       outerDeriv[0] += outerDeriv[1];
00354     // the term belonging to $I_1$...
00355     dPhi *= 2 * outerDeriv[0];
00356     Stress += dPhi;
00357     // the term belonging to $I_3$...
00358     cof *= outerDeriv[2];
00359     Stress += cof;
00360   }
00361 };
00362 
00363 /**
00364   * \brief This class computes via "apply(...)" or "applyAdd(...)" parts of a hyperelastic energy Hessian with respect to the displacement,
00365   * i.e. \f$ (<\frac{\partial^2 \int_\Omega W(I_1,I_2,I_3,x) dx}{\partial \phi^2},\psi_i,\psi_j>)_{i,j} \f$,
00366   * where the \f$ \psi_i \f$ and \f$ \psi_i \f$ are those vector-valued FE basis functions that are only nonzero in the kth and lth component,
00367   * the hyperelastic invariants are given by \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2 \f$,
00368   * \f$ I_3=det(\nabla\phi)\f$, and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
00369   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$
00370   * in the different Cartesian directions; the deformation \f$ \phi \f$ is then automatically computed via \f$ \phi=d+\f$ identity.
00371   *
00372   * \author Wirth
00373   */
00374 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
00375 class HyperelasticSubHessian
00376       : public aol::FELinAsymMatrixWeightedStiffInterface<ConfiguratorType, HyperelasticSubHessian<ConfiguratorType,HyperelasticEnergyDensityType> > {
00377 protected:
00378   typedef typename ConfiguratorType::RealType RealType;
00379 
00380   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
00381   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
00382   // the displacement, at which the sub-Hessian shall be evaluated
00383   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement;
00384   // the indices of the sub-Hessian
00385   const int _k, _l;
00386 public:
00387   HyperelasticSubHessian ( const typename ConfiguratorType::InitType &Grid,
00388                            const HyperelasticEnergyDensityType &HyperelasticEnergyDensity,
00389                            const aol::MultiVector<RealType> &Displacement,
00390                            const int K,
00391                            const int L )
00392       : aol::FELinAsymMatrixWeightedStiffInterface<ConfiguratorType, HyperelasticSubHessian<ConfiguratorType,HyperelasticEnergyDensityType> > ( Grid ),
00393         _hyperelasticEnergyDensity( HyperelasticEnergyDensity ),
00394         _displacement( Grid, Displacement ),
00395         _k( L ), // l and k are exchanged here, since FELinAsymMatrixStiffOp assembles the matrix in the way that \int \nabla\phi^T A \nabla\psi dx becomes \Psi^T L \Phi
00396         _l( K ) {}
00397 
00398   inline void getCoeffMatrix ( const typename ConfiguratorType::ElementType &El, int QuadPoint,
00399                                const typename ConfiguratorType::DomVecType& /*DomRefCoord*/, typename ConfiguratorType::MatType &Matrix ) const {
00400     // compute the deformation gradient $\nabla\phi$
00401     typename ConfiguratorType::MatType dPhi;
00402     _displacement.evaluateGradientAtQuadPoint ( El, QuadPoint, dPhi );
00403     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00404       dPhi[i][i] += 1.;
00405 
00406     // compute the Hessian matrix for the kth component multiplied from the left and the lth from the right
00407     elasticityTensor( dPhi, _hyperelasticEnergyDensity, El, QuadPoint, _k, _l, Matrix );
00408   }
00409 
00410   /**
00411     * \brief Computes ElasticityTensor=\f$ \frac{\partial^2W}{\partial\nabla\phi_K\partial\nabla\phi_L}(I_1,I_2,I_3) \f$ for the invariants
00412     * \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2\f$, \f$ I_3=det(\nabla\phi) \f$,
00413     * where \f$ \nabla\phi \f$, \f$ W \f$, and the spatial position are passed as arguments.
00414     * The result matrix \f$ C= \f$ ElasticityTensor is such that
00415     * \f$ <\frac{\partial^2W}{\partial\nabla\phi_K\partial\nabla\phi_L},\theta_K,\psi_L>=\theta_K\cdot(C\psi_L) \f$
00416     * for a variation \f$ \theta_K \f$ of \f$ \nabla\phi_K \f$ and \f$ \psi_L \f$ of \f$ \nabla\phi_L \f$.
00417     */
00418   static inline void elasticityTensor( const typename ConfiguratorType::MatType &DPhi,
00419                                        const HyperelasticEnergyDensityType &EnergyDensity,
00420                                        const typename ConfiguratorType::ElementType &El, int QuadPoint,
00421                                        const int K, const int L,
00422                                        aol::Mat< ConfiguratorType::Dim, ConfiguratorType::Dim, RealType > &ElasticityTensor ) {
00423     // the cofactor matrix of the deformation gradient
00424     typename ConfiguratorType::MatType cof;
00425     // auxiliary vectors and matrices
00426     typename ConfiguratorType::VecType auxVec1;
00427     typename ConfiguratorType::MatType auxMat1, auxMat2;
00428 
00429     ElasticityTensor.setZero();
00430 
00431     // compute the cofactor matrix of the deformation gradient, $cof(\nabla\phi)$
00432     cof.makeCofactorMatrix( DPhi );
00433     // compute the hyperelastic invariants $||\nabla\phi||_F^2,||cof(\nabla\phi)||_F^2,det(\nabla\phi)$
00434     RealType
00435       I1 = DPhi.normSqr(),
00436       I2 = cof.normSqr(),
00437       I3 = DPhi.det();
00438 
00439     // compute the outer first and second derivative of the hyperelastic energy
00440     aol::Vec3<RealType> outerDeriv ( EnergyDensity.evaluateDerivativeAtQuadPoint ( I1, I2, I3, El, QuadPoint ) );
00441     aol::Matrix33<RealType> outerSecondDeriv ( EnergyDensity.evaluateSecondDerivativeAtQuadPoint ( I1, I2, I3, El, QuadPoint ) );
00442 
00443     // compute auxiliary matrices h, H
00444     typename ConfiguratorType::MatType H( cof.transposed() ), h, hCof;
00445     h.makeProduct( cof, H );
00446     RealType auxScal1 = h.tr();
00447     auxMat1 = h;
00448     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00449       auxMat1[i][i] -= auxScal1;
00450     H.makeProduct( auxMat1, cof );
00451     H /= - I3;
00452     hCof.makeProduct( h, cof );
00453 
00454     // compute the Hessian matrix for the kth component multiplied from the left and the lth from the right
00455     if ( K == L )
00456       for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00457         ElasticityTensor[i][i] = 2 * outerDeriv[0]; // note: \epsilon:\varespilon = \sum_{k,l=1..d}\delta_{kl}\epsilon_{k:}I\varepsilon_{l:}^T
00458     auxVec1.setMultiple( DPhi[K], 4 * outerSecondDeriv[0][0] );
00459     auxVec1.addMultiple( H[K], 4 * outerSecondDeriv[0][1] );
00460     auxVec1.addMultiple( cof[K], 2 * outerSecondDeriv[0][2] );
00461     auxMat1.makeTensorProduct( auxVec1, DPhi[L] );
00462     ElasticityTensor += auxMat1; // note: (A:\epsilon)(B:\varepsilon) = \sum_{k,l=1..d}\epsilon_{k:}A_{k:}^T B_{l:}\varepsilon_{l:}^T
00463 
00464     auxVec1.setMultiple( cof[L], - outerDeriv[2] / I3 );
00465     auxMat1.makeTensorProduct( auxVec1, cof[K] );
00466     ElasticityTensor += auxMat1; // note: tr(A\epsilon^T B\varepsilon^T) = \sum_{k,l=1..d}\epsilon_{k:}A_{l:}^T B_{k:}\varepsilon_{l:}^T
00467     auxVec1.setMultiple( cof[K], outerDeriv[2] / I3 );
00468     auxVec1.addMultiple( DPhi[K], 2 * outerSecondDeriv[2][0] );
00469     auxVec1.addMultiple( H[K], 2 * outerSecondDeriv[2][1] );
00470     auxVec1.addMultiple( cof[K], outerSecondDeriv[2][2] );
00471     auxMat1.makeTensorProduct( auxVec1, cof[L] );
00472     ElasticityTensor += auxMat1;
00473 
00474     auxVec1.setMultiple( cof[K], - 2 * outerDeriv[1] / aol::Sqr( I3 ) );
00475     auxMat1.makeTensorProduct( auxVec1, hCof[L] );
00476     ElasticityTensor += auxMat1;
00477     auxVec1.setMultiple( cof[L], 2 * outerDeriv[1] / aol::Sqr( I3 ) );
00478     auxMat1.makeTensorProduct( auxVec1, hCof[K] );
00479     ElasticityTensor += auxMat1;
00480     auxMat2.transposeFrom( cof );
00481     auxMat1.makeProduct( auxMat2, cof );
00482     auxMat1 *= 2 * outerDeriv[1] / aol::Sqr( I3 ) * h[L][K];
00483     ElasticityTensor += auxMat1; // note: tr(A\epsilon B\varepsilon^T) = \sum_{k,l=1..d}\epsilon_{k:}A_{lk}B\varepsilon_{l:}^T
00484     auxVec1.setMultiple( cof[K],  2 * outerDeriv[1] / aol::Sqr( I3 ) * h.tr() );
00485     auxVec1.addMultiple( hCof[K], - 4 * outerDeriv[1] / aol::Sqr( I3 ) );
00486     auxMat1.makeTensorProduct( auxVec1, cof[L] );
00487     ElasticityTensor += auxMat1;
00488     auxVec1.setMultiple( H[L], - 2 * outerDeriv[1] / I3 );
00489     auxMat1.makeTensorProduct( auxVec1, cof[K] );
00490     ElasticityTensor += auxMat1;
00491     auxVec1.setMultiple( cof[K], 2 * outerDeriv[1] / I3 );
00492     auxVec1.addMultiple( DPhi[K], 4 * outerSecondDeriv[1][0] );
00493     auxVec1.addMultiple( H[K], 4 * outerSecondDeriv[1][1] );
00494     auxVec1.addMultiple( cof[K], 2 * outerSecondDeriv[1][2] );
00495     auxMat1.makeTensorProduct( auxVec1, H[L] );
00496     ElasticityTensor += auxMat1;
00497   }
00498 };
00499 
00500 /**
00501   * \brief This class computes via "apply(...)" or "applyAdd(...)" a hyperelastic energy Hessian with respect to the displacement,
00502   * i.e. \f$ (<\frac{\partial^2 \int_\Omega W(I_1,I_2,I_3,x) dx}{\partial \phi^2},\psi_i,\psi_j>)_{i,j} \f$, where the \f$ \psi_i \f$ are the FE basis
00503   * functions, the hyperelastic invariants are given by \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2 \f$,
00504   * \f$ I_3=det(\nabla\phi)\f$ and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
00505   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$
00506   * in the different Cartesian directions; the deformation \f$ \phi \f$ is then automatically computed via \f$ \phi=d+\f$ identity.
00507   *
00508   * \author Wirth
00509   */
00510 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType, typename SubMatrixType = typename ConfiguratorType::MatrixType>
00511 class HyperelasticHessian
00512       : public aol::Op<aol::MultiVector<typename ConfiguratorType::RealType>, aol::BlockMatrix<SubMatrixType> > {
00513 
00514 protected:
00515   typedef typename ConfiguratorType::RealType RealType;
00516 
00517   // the underlying FE grid
00518   const typename ConfiguratorType::InitType &_grid;
00519   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
00520   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
00521 
00522 public:
00523   HyperelasticHessian ( const typename ConfiguratorType::InitType &Grid,
00524                         const HyperelasticEnergyDensityType &HyperelasticEnergyDensity )
00525       : _grid ( Grid ),
00526         _hyperelasticEnergyDensity( HyperelasticEnergyDensity ) {}
00527 
00528   void applyAdd( const aol::MultiVector<RealType> &Arg, aol::BlockMatrix<SubMatrixType> &Dest ) const {
00529     for ( int k = 0; k < ConfiguratorType::Dim; k++ )
00530       for ( int l = 0; l < ConfiguratorType::Dim; l++ )
00531         HyperelasticSubHessian<ConfiguratorType,HyperelasticEnergyDensityType>( _grid, _hyperelasticEnergyDensity, Arg, k, l ).assembleAddMatrix( Dest.getReference( k, l ) );
00532   }
00533 };
00534 
00535 /**
00536   * \brief This class computes via "apply(...)" or "applyAdd(...)" a hyperelastic energy \f$ \int_\Omega W(I_1,I_2,I_3,x) dx \f$,
00537   * where the hyperelastic invariants are given by \f$ I_1=||\nabla\phi\nabla\Phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi\nabla\Phi)||_F^2 \f$,
00538   * \f$ I_1=det(\nabla\phi\nabla\Phi) \f$ and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
00539   * \f$ \Phi \f$ is a deformation passed to the constructor.
00540   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$
00541   * in the different Cartesian directions; the deformation \f$ \phi \f$ is then automatically computed via \f$ \phi=d+\f$ identity.
00542   *
00543   * \author Wirth
00544   */
00545 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
00546 class HyperelasticDeformEnergy
00547       : public FENonlinDeformIntegrationVectorInterface < ConfiguratorType, HyperelasticDeformEnergy<ConfiguratorType,HyperelasticEnergyDensityType> > {
00548 protected:
00549   typedef typename ConfiguratorType::RealType RealType;
00550 
00551   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
00552   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
00553   // the deformation \Phi
00554   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement;
00555 public:
00556   HyperelasticDeformEnergy ( const typename ConfiguratorType::InitType &Grid,
00557                              const HyperelasticEnergyDensityType &HyperelasticEnergyDensity,
00558                              const aol::MultiVector<RealType> &Displacement ) :
00559     FENonlinDeformIntegrationVectorInterface < ConfiguratorType, HyperelasticDeformEnergy<ConfiguratorType,HyperelasticEnergyDensityType> > ( Grid, Displacement ),
00560     _hyperelasticEnergyDensity( HyperelasticEnergyDensity ),
00561     _displacement( Grid, Displacement ) {}
00562 
00563   RealType evaluateIntegrand ( const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> &DiscFuncs,
00564                                const typename ConfiguratorType::ElementType &El, int QuadPoint, const typename ConfiguratorType::DomVecType &/*LocalCoord*/,
00565                                const typename ConfiguratorType::ElementType &TransformedEl, const typename ConfiguratorType::DomVecType &TransformedLocalCoord ) const {
00566     // compute the deformation gradient $\nabla\Phi$
00567     typename ConfiguratorType::MatType dPhi;
00568     _displacement.evaluateGradientAtQuadPoint( El, QuadPoint, dPhi );
00569     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00570       dPhi[i][i] += 1.;
00571 
00572     // compute the deformation gradient $\nabla\phi$
00573     typename ConfiguratorType::MatType dphi;
00574     DiscFuncs.evaluateGradient( TransformedEl, TransformedLocalCoord, dphi );
00575     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00576       dphi[i][i] += 1.;
00577 
00578     dphi *= dPhi;
00579     return HyperelasticEnergy<ConfiguratorType,HyperelasticEnergyDensityType>::energyDensity( dphi, _hyperelasticEnergyDensity, El, QuadPoint );
00580   }
00581 };
00582 
00583 /**
00584   * \brief This class computes via "apply(...)" or "applyAdd(...)" a hyperelastic energy gradient with respect to the displacement,
00585   * i.e. \f$ (<\frac{\partial \int_\Omega W(I_1,I_2,I_3,x) dx}{\partial \phi},\psi_i>)_i \f$, where the \f$ \psi_i \f$ are the FE basis
00586   * functions, the hyperelastic invariants are given by \f$ I_1=||\nabla\phi\nabla\Phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi\nabla\Phi)||_F^2 \f$,
00587   * \f$ I_3=det(\nabla\phi\nabla\Phi)\f$ and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
00588   * \f$ \Phi \f$ is some prestressing deformation.
00589   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$
00590   * in the different Cartesian directions; the deformation \f$ \phi \f$ is then automatically computed via \f$ \phi=d+\f$ identity.
00591   *
00592   * \author Wirth
00593   */
00594 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
00595 class HyperelasticDeformGradient
00596       : public FENonlinDeformVectorDiffOpInterface<ConfiguratorType, ConfiguratorType::Dim, ConfiguratorType::Dim, HyperelasticDeformGradient<ConfiguratorType,HyperelasticEnergyDensityType> > {
00597 protected:
00598   typedef typename ConfiguratorType::RealType RealType;
00599 
00600   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
00601   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
00602   // the deformation \Phi
00603   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement;
00604 public:
00605   HyperelasticDeformGradient ( const typename ConfiguratorType::InitType &Grid,
00606                                const HyperelasticEnergyDensityType &HyperelasticEnergyDensity,
00607                                const aol::MultiVector<RealType> &Displacement ) :
00608     FENonlinDeformVectorDiffOpInterface<ConfiguratorType, ConfiguratorType::Dim, ConfiguratorType::Dim, HyperelasticDeformGradient<ConfiguratorType,HyperelasticEnergyDensityType> > ( Grid, Displacement ),
00609     _hyperelasticEnergyDensity( HyperelasticEnergyDensity ),
00610     _displacement( Grid, Displacement ) {}
00611 
00612   void getNonlinearity ( const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> &DiscFuncs,
00613                          const typename ConfiguratorType::ElementType &El, int QuadPoint, const typename ConfiguratorType::DomVecType &/*LocalCoord*/,
00614                          const typename ConfiguratorType::ElementType &TransformedEl, const typename ConfiguratorType::DomVecType &TransformedLocalCoord,
00615                          aol::Mat<ConfiguratorType::Dim, ConfiguratorType::Dim, RealType> &NL ) const {
00616     // compute the deformation gradient $\nabla\Phi$
00617     typename ConfiguratorType::MatType dPhi;
00618     _displacement.evaluateGradientAtQuadPoint( El, QuadPoint, dPhi );
00619     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00620       dPhi[i][i] += 1.;
00621 
00622     // compute the deformation gradient $\nabla\phi$
00623     typename ConfiguratorType::MatType dphi;
00624     DiscFuncs.evaluateGradient( TransformedEl, TransformedLocalCoord, dphi );
00625     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00626       dphi[i][i] += 1.;
00627 
00628     dphi *= dPhi;
00629     HyperelasticGradient<ConfiguratorType,HyperelasticEnergyDensityType>::firstPiolaKirchhoffStress( dphi, _hyperelasticEnergyDensity, El, QuadPoint, NL );
00630   }
00631 };
00632 
00633 /**
00634   * \brief This class computes via "apply(...)" or "applyAdd(...)" parts of a hyperelastic energy Hessian with respect to the displacement,
00635   * where the material is prestressed by a deformation \f$ \Phi \f$,
00636   * i.e. \f$ (<\frac{\partial^2 \int_\Omega W(I_1,I_2,I_3,x) dx}{\partial\phi^2},\psi_i\circ\Phi,\psi_j\circ\Phi>)_{i,j} \f$,
00637   * where the \f$ \psi_i \f$ and \f$ \psi_j \f$ are those vector-valued FE basis functions that are only nonzero in the kth and lth component,
00638   * the hyperelastic invariants are given by \f$ I_1=||\nabla\phi\nabla\Phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi\nabla\Phi)||_F^2 \f$,
00639   * \f$ I_3=det(\nabla\phi\nabla\Phi)\f$, and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
00640   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$
00641   * in the different Cartesian directions; the deformation \f$ \phi \f$ is then automatically computed via \f$ \phi=d+\f$ identity.
00642   *
00643   * \author Wirth
00644   */
00645 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
00646 class HyperelasticDeformSubHessian
00647       : public qc::FELinDeformAsymMatrixWeightedStiffInterface<ConfiguratorType, HyperelasticDeformSubHessian<ConfiguratorType,HyperelasticEnergyDensityType> > {
00648 protected:
00649   typedef typename ConfiguratorType::RealType RealType;
00650 
00651   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
00652   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
00653   // the deformations \Phi and \phi
00654   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement, _argDisplacement;
00655   // the indices of the sub-Hessian
00656   const int _k, _l;
00657 
00658 public:
00659   HyperelasticDeformSubHessian ( const typename ConfiguratorType::InitType &Grid,
00660                                  const HyperelasticEnergyDensityType &HyperelasticEnergyDensity,
00661                                  const aol::MultiVector<RealType> &Displacement,
00662                                  const aol::MultiVector<RealType> &ArgDisplacement,
00663                                  const int K,
00664                                  const int L ) :
00665     qc::FELinDeformAsymMatrixWeightedStiffInterface<ConfiguratorType, HyperelasticDeformSubHessian<ConfiguratorType,HyperelasticEnergyDensityType> > ( Grid, Displacement ),
00666     _hyperelasticEnergyDensity( HyperelasticEnergyDensity ),
00667     _displacement( Grid, Displacement ),
00668     _argDisplacement( Grid, ArgDisplacement ),
00669     _k( K ),
00670     _l( L ){}
00671 
00672   inline void getCoeffMatrix ( const typename ConfiguratorType::ElementType &El, int QuadPoint, const typename ConfiguratorType::DomVecType &/*LocalCoord*/,
00673                                const typename ConfiguratorType::ElementType &TransformedEl, const typename ConfiguratorType::DomVecType &TransformedLocalCoord,
00674                                typename ConfiguratorType::MatType &Matrix ) const {
00675     // compute the deformation gradient $\nabla\Phi$
00676     typename ConfiguratorType::MatType dPhi;
00677     _displacement.evaluateGradientAtQuadPoint( El, QuadPoint, dPhi );
00678     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00679       dPhi[i][i] += 1.;
00680 
00681     // compute the deformation gradient $\nabla\phi$
00682     typename ConfiguratorType::MatType dphi;
00683     _argDisplacement.evaluateGradient( TransformedEl, TransformedLocalCoord, dphi );
00684     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00685       dphi[i][i] += 1.;
00686 
00687     // compute the Hessian matrix for the kth component multiplied from the left and the lth from the right
00688     dphi *= dPhi;
00689     HyperelasticSubHessian<ConfiguratorType,HyperelasticEnergyDensityType>::elasticityTensor( dphi, _hyperelasticEnergyDensity, El, QuadPoint, _k, _l, Matrix );
00690   }
00691 };
00692 
00693 /**
00694   * \brief This class computes via "apply(...)" or "applyAdd(...)" a hyperelastic energy Hessian with respect to the displacement,
00695   * where the material is prestressed by a deformation \f$ \Phi \f$,
00696   * i.e. \f$ (<\frac{\partial^2 \int_\Omega W(I_1,I_2,I_3,x) dx}{\partial \phi^2},\psi_i\circ\Phi,\psi_j\circ\Phi>)_{i,j} \f$,
00697   * where the \f$ \psi_i \f$ and \f$ \psi_j \f$ are the vector-valued FE basis functions,
00698   * the hyperelastic invariants are given by \f$ I_1=||\nabla\phi\nabla\Phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi\nabla\Phi)||_F^2 \f$,
00699   * \f$ I_3=det(\nabla\phi\nabla\Phi)\f$, and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
00700   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$
00701   * in the different Cartesian directions; the deformation \f$ \phi \f$ is then automatically computed via \f$ \phi=d+\f$ identity.
00702   *
00703   * \author Wirth
00704   */
00705 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType, typename SubMatrixType = typename ConfiguratorType::MatrixType>
00706 class HyperelasticDeformHessian
00707       : public aol::Op<aol::MultiVector<typename ConfiguratorType::RealType>, aol::BlockMatrix<SubMatrixType> > {
00708 
00709 protected:
00710   typedef typename ConfiguratorType::RealType RealType;
00711 
00712   // the underlying FE grid
00713   const typename ConfiguratorType::InitType &_grid;
00714   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
00715   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
00716   // the prestressing displacement
00717   const aol::MultiVector<RealType> &_displacement;
00718 
00719 public:
00720   HyperelasticDeformHessian ( const typename ConfiguratorType::InitType &Grid,
00721                               const HyperelasticEnergyDensityType &HyperelasticEnergyDensity,
00722                               const aol::MultiVector<RealType> &Displacement ) :
00723     _grid ( Grid ),
00724     _hyperelasticEnergyDensity( HyperelasticEnergyDensity ),
00725     _displacement( Displacement ){}
00726 
00727   void applyAdd( const aol::MultiVector<RealType> &Arg, aol::BlockMatrix<SubMatrixType> &Dest ) const {
00728     for ( int k = 0; k < ConfiguratorType::Dim; k++ )
00729       for ( int l = 0; l < ConfiguratorType::Dim; l++ )
00730         HyperelasticDeformSubHessian<ConfiguratorType,HyperelasticEnergyDensityType>( _grid, _hyperelasticEnergyDensity, _displacement, Arg, k, l ).assembleAddMatrix( Dest.getReference( k, l ) );
00731   }
00732 
00733   void applyAdd( const aol::MultiVector<RealType> &Arg, aol::BlockMatrix<SubMatrixType> &Dest, const aol::BitVector &Mask ) const {
00734     for ( int k = 0; k < ConfiguratorType::Dim; k++ )
00735       for ( int l = 0; l < ConfiguratorType::Dim; l++ )
00736         HyperelasticDeformSubHessian<ConfiguratorType,HyperelasticEnergyDensityType>( _grid, _hyperelasticEnergyDensity, _displacement, Arg, k, l ).assembleAddMatrix( Mask, Dest.getReference( k, l ), aol::NumberTrait<RealType>::one, k == l );
00737   }
00738 };
00739 
00740 /**
00741   * \brief This class computes via "apply(...)" or "applyAdd(...)" a hyperelastic energy \f$ \int_\Omega W(I_1,I_2,I_3,x)det(\nabla\phi_r) dx \f$,
00742   * where the hyperelastic invariants are given by \f$ I_1=||\nabla\phi_l\nabla\phi\nabla\phi_r^{-1}||_F^2 \f$, \f$ I_2=||cof(\nabla\phi_l\nabla\phi\nabla\phi_r^{-1})||_F^2 \f$,
00743   * \f$ I_1=det(\nabla\phi_l\nabla\phi\nabla\phi_r^{-1}) \f$ and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
00744   * \f$ \phi \f$ is a deformation passed to the constructor.
00745   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$ of \f$ \phi_l \f$ and \f$ \phi_r \f$
00746   * in the different Cartesian directions.
00747   *
00748   * \author Wirth
00749   */
00750 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
00751 class HyperelasticPrePostDeformEnergy
00752       : public FENonlinDeformIntegrationVectorInterface < ConfiguratorType, HyperelasticPrePostDeformEnergy<ConfiguratorType,HyperelasticEnergyDensityType>, 2*ConfiguratorType::Dim > {
00753 protected:
00754   typedef typename ConfiguratorType::RealType RealType;
00755 
00756   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
00757   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
00758   // the deformation \phi
00759   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement;
00760 public:
00761   HyperelasticPrePostDeformEnergy ( const typename ConfiguratorType::InitType &Grid,
00762                                     const HyperelasticEnergyDensityType &HyperelasticEnergyDensity,
00763                                     const aol::MultiVector<RealType> &Displacement ) :
00764     FENonlinDeformIntegrationVectorInterface < ConfiguratorType, HyperelasticPrePostDeformEnergy<ConfiguratorType,HyperelasticEnergyDensityType>, 2*ConfiguratorType::Dim > ( Grid, Displacement ),
00765     _hyperelasticEnergyDensity( HyperelasticEnergyDensity ),
00766     _displacement( Grid, Displacement ) {}
00767 
00768   RealType evaluateIntegrand ( const aol::DiscreteVectorFunctionDefault<ConfiguratorType,2*ConfiguratorType::Dim> &DiscFuncs,
00769                                const typename ConfiguratorType::ElementType &El, int QuadPoint, const typename ConfiguratorType::DomVecType &LocalCoord,
00770                                const typename ConfiguratorType::ElementType &TransformedEl, const typename ConfiguratorType::DomVecType &TransformedLocalCoord ) const {
00771     // compute the deformation gradient $\nabla\phi$
00772     typename ConfiguratorType::MatType dPhi;
00773     _displacement.evaluateGradientAtQuadPoint( El, QuadPoint, dPhi );
00774     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00775       dPhi[i][i] += 1.;
00776 
00777     // compute the deformation gradient $\nabla\phi_l$ and $\nabla\phi_r$
00778     typename ConfiguratorType::MatType dPhiL, dPhiR;
00779     for ( int i = 0; i < ConfiguratorType::Dim; i++ ){
00780       DiscFuncs[i].evaluateGradient( TransformedEl, TransformedLocalCoord, dPhiL[i] );
00781       DiscFuncs[ConfiguratorType::Dim+i].evaluateGradient( El, LocalCoord, dPhiR[i] );
00782       dPhiL[i][i] += 1.;
00783       dPhiR[i][i] += 1.;
00784     }
00785 
00786     dPhiL *= dPhi;
00787     dPhiL *= dPhiR.inverse();
00788     return dPhiR.det() * HyperelasticEnergy<ConfiguratorType,HyperelasticEnergyDensityType>::energyDensity( dPhiL, _hyperelasticEnergyDensity, El, QuadPoint );
00789   }
00790 };
00791 
00792 /**
00793   * \brief This class computes via "apply(...)" or "applyAdd(...)" a hyperelastic energy gradient with respect to one of two displacements,
00794   * i.e. \f$ (<\frac{\partial \int_\Omega W(I_1,I_2,I_3,x)det(\nabla\phi_r) dx}{\partial \phi_l},\psi_i>)_i \f$,
00795   * where the \f$ \psi_i \f$ are the vector-valued FE basis functions,
00796   * the hyperelastic invariants are given by \f$ I_1=||\nabla\phi_l\nabla\phi\nabla\phi_r^{-1}||_F^2 \f$, \f$ I_2=||cof(\nabla\phi_l\nabla\phi\nabla\phi_r^{-1})||_F^2 \f$,
00797   * \f$ I_1=det(\nabla\phi_l\nabla\phi\nabla\phi_r^{-1}) \f$ and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
00798   * \f$ \phi \f$ is some prestressing deformation.
00799   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$ of \f$ \phi_l \f$ and \f$ \phi_r \f$
00800   * in the different Cartesian directions.
00801   *
00802   * \author Wirth
00803   */
00804 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
00805 class HyperelasticPrePostDeformGradientPostComp
00806       : public FENonlinDeformVectorDiffOpInterface<ConfiguratorType, 2*ConfiguratorType::Dim, ConfiguratorType::Dim, HyperelasticPrePostDeformGradientPostComp<ConfiguratorType,HyperelasticEnergyDensityType> > {
00807 protected:
00808   typedef typename ConfiguratorType::RealType RealType;
00809 
00810   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
00811   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
00812   // the deformation \phi
00813   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement;
00814 public:
00815   HyperelasticPrePostDeformGradientPostComp ( const typename ConfiguratorType::InitType &Grid,
00816                                               const HyperelasticEnergyDensityType &HyperelasticEnergyDensity,
00817                                               const aol::MultiVector<RealType> &Displacement ) :
00818     FENonlinDeformVectorDiffOpInterface<ConfiguratorType, 2*ConfiguratorType::Dim, ConfiguratorType::Dim, HyperelasticPrePostDeformGradientPostComp<ConfiguratorType,HyperelasticEnergyDensityType> > ( Grid, Displacement ),
00819     _hyperelasticEnergyDensity( HyperelasticEnergyDensity ),
00820     _displacement( Grid, Displacement ) {}
00821 
00822   void getNonlinearity ( const aol::DiscreteVectorFunctionDefault<ConfiguratorType,2*ConfiguratorType::Dim> &DiscFuncs,
00823                          const typename ConfiguratorType::ElementType &El, int QuadPoint, const typename ConfiguratorType::DomVecType &LocalCoord,
00824                          const typename ConfiguratorType::ElementType &TransformedEl, const typename ConfiguratorType::DomVecType &TransformedLocalCoord,
00825                          aol::Mat<ConfiguratorType::Dim, ConfiguratorType::Dim, RealType> &NL ) const {
00826     // compute the deformation gradient $\nabla\phi$
00827     typename ConfiguratorType::MatType dPhi;
00828     _displacement.evaluateGradientAtQuadPoint( El, QuadPoint, dPhi );
00829     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00830       dPhi[i][i] += 1.;
00831 
00832     // compute the deformation gradient $\nabla\phi_l$ and $\nabla\phi_r$
00833     typename ConfiguratorType::MatType dPhiL, dPhiR;
00834     for ( int i = 0; i < ConfiguratorType::Dim; i++ ){
00835       DiscFuncs[i].evaluateGradient( TransformedEl, TransformedLocalCoord, dPhiL[i] );
00836       DiscFuncs[ConfiguratorType::Dim+i].evaluateGradient( El, LocalCoord, dPhiR[i] );
00837       dPhiL[i][i] += 1.;
00838       dPhiR[i][i] += 1.;
00839     }
00840 
00841     // compute \partial_{,A}W(I_1,I_2,I_3)
00842     dPhiL *= dPhi;
00843     dPhiL *= dPhiR.inverse();
00844     HyperelasticGradient<ConfiguratorType,HyperelasticEnergyDensityType>::firstPiolaKirchhoffStress( dPhiL, _hyperelasticEnergyDensity, El, QuadPoint, NL );
00845 
00846     // compute NL = det(\nabla\phi_r) \partial_{,A}W(I_1,I_2,I_3) \nabla\phi_r^{-T}
00847     typename ConfiguratorType::MatType cofDPhiR;
00848     cofDPhiR.makeCofactorMatrix( dPhiR );
00849     NL *= cofDPhiR;
00850   }
00851 };
00852 
00853 /**
00854   * \brief This class computes via "apply(...)" or "applyAdd(...)" a hyperelastic energy gradient with respect to one of two displacements,
00855   * i.e. \f$ (<\frac{\partial \int_\Omega W(I_1,I_2,I_3,x)det(\nabla\phi_r) dx}{\partial \phi_r},\psi_i>)_i \f$,
00856   * where the \f$ \psi_i \f$ are the vector-valued FE basis functions,
00857   * the hyperelastic invariants are given by \f$ I_1=||\nabla\phi_l\nabla\phi\nabla\phi_r^{-1}||_F^2 \f$, \f$ I_2=||cof(\nabla\phi_l\nabla\phi\nabla\phi_r^{-1})||_F^2 \f$,
00858   * \f$ I_1=det(\nabla\phi_l\nabla\phi\nabla\phi_r^{-1}) \f$ and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
00859   * \f$ \phi \f$ is some prestressing deformation.
00860   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$ of \f$ \phi_l \f$ and \f$ \phi_r \f$
00861   * in the different Cartesian directions.
00862   *
00863   * \author Wirth
00864   */
00865 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
00866 class HyperelasticPrePostDeformGradientPreComp
00867       : public aol::FENonlinVectorDiffOpInterface<ConfiguratorType, 2*ConfiguratorType::Dim, ConfiguratorType::Dim, HyperelasticPrePostDeformGradientPreComp<ConfiguratorType,HyperelasticEnergyDensityType> > {
00868 protected:
00869   typedef typename ConfiguratorType::RealType RealType;
00870 
00871   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
00872   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
00873   // the deformation \phi
00874   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement;
00875 public:
00876   HyperelasticPrePostDeformGradientPreComp ( const typename ConfiguratorType::InitType &Grid,
00877                                              const HyperelasticEnergyDensityType &HyperelasticEnergyDensity,
00878                                              const aol::MultiVector<RealType> &Displacement ) :
00879     aol::FENonlinVectorDiffOpInterface<ConfiguratorType, 2*ConfiguratorType::Dim, ConfiguratorType::Dim, HyperelasticPrePostDeformGradientPreComp<ConfiguratorType,HyperelasticEnergyDensityType> > ( Grid ),
00880     _hyperelasticEnergyDensity( HyperelasticEnergyDensity ),
00881     _displacement( Grid, Displacement ) {}
00882 
00883   void getNonlinearity ( const aol::auto_container<2*ConfiguratorType::Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &DiscFuncs,
00884                          const typename ConfiguratorType::ElementType &El, int QuadPoint, const typename ConfiguratorType::DomVecType &LocalCoord,
00885                          aol::Mat<ConfiguratorType::Dim, ConfiguratorType::Dim, RealType> &NL ) const {
00886     // compute the deformation gradient $\nabla\phi$
00887     typename ConfiguratorType::MatType dPhi;
00888     _displacement.evaluateGradientAtQuadPoint( El, QuadPoint, dPhi );
00889     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00890       dPhi[i][i] += 1.;
00891 
00892     // compute $\phi(x)$
00893     typename ConfiguratorType::ElementType transformedEl;
00894     typename ConfiguratorType::DomVecType  transformedCoord;
00895     qc::transformAndClipCoord<ConfiguratorType>( *this->_config, _displacement, El, QuadPoint, LocalCoord, transformedEl, transformedCoord );
00896 
00897     // compute the deformation gradient $\nabla\phi_l$ and $\nabla\phi_r$
00898     typename ConfiguratorType::MatType dPhiL, dPhiR;
00899     for ( int i = 0; i < ConfiguratorType::Dim; i++ ){
00900       DiscFuncs[i].evaluateGradient( transformedEl, transformedCoord, dPhiL[i] );
00901       DiscFuncs[ConfiguratorType::Dim+i].evaluateGradient( El, LocalCoord, dPhiR[i] );
00902       dPhiL[i][i] += 1.;
00903       dPhiR[i][i] += 1.;
00904     }
00905 
00906     // compute W(I_1,I_2,I_3) and \partial_{,A}W(I_1,I_2,I_3)
00907     dPhiL *= dPhi;
00908     dPhiL *= dPhiR.inverse();
00909     const RealType density = HyperelasticEnergy<ConfiguratorType,HyperelasticEnergyDensityType>::energyDensity( dPhiL, _hyperelasticEnergyDensity, El, QuadPoint );
00910     HyperelasticGradient<ConfiguratorType,HyperelasticEnergyDensityType>::firstPiolaKirchhoffStress( dPhiL, _hyperelasticEnergyDensity, El, QuadPoint, NL );
00911 
00912     // compute NL = W(I_1,I_2,I_3) cof(\nabla\phi_r) + det(\nabla\phi_r) \partial_{,A}W(I_1,I_2,I_3) 
00913     typename ConfiguratorType::MatType cofDPhiR, aux;
00914     cofDPhiR.makeCofactorMatrix( dPhiR );
00915     aux.makeProduct( NL, cofDPhiR );
00916     NL.makeProductAtransposedB( dPhiL, aux );
00917     NL *= -1.;
00918     cofDPhiR *= density;
00919     NL += cofDPhiR;
00920   }
00921 };
00922 
00923 /**
00924   * \brief This class computes via "apply(...)" or "applyAdd(...)" part of the hyperelastic energy Hessian with respect to one of two displacements,
00925   * i.e. \f$ (<\frac{\partial^2\int_\Omega W(I_1,I_2,I_3,x)det(\nabla\phi_r) dx}{\partial \phi_l^2},\psi_i,\psi_j>)_{ij} \f$,
00926   * where the \f$ \psi_i \f$ are the vector-valued FE basis functions,
00927   * the hyperelastic invariants are given by \f$ I_1=||\nabla\phi_l\nabla\phi\nabla\phi_r^{-1}||_F^2 \f$, \f$ I_2=||cof(\nabla\phi_l\nabla\phi\nabla\phi_r^{-1})||_F^2 \f$,
00928   * \f$ I_1=det(\nabla\phi_l\nabla\phi\nabla\phi_r^{-1}) \f$ and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
00929   * \f$ \phi \f$ is some prestressing deformation.
00930   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$ of \f$ \phi_l \f$ and \f$ \phi_r \f$
00931   * in the different Cartesian directions.
00932   *
00933   * \author Wirth
00934   */
00935 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
00936 class HyperelasticPrePostDeformSubHessianPostComp
00937       : public qc::FELinDeformAsymMatrixWeightedStiffInterface<ConfiguratorType, HyperelasticPrePostDeformSubHessianPostComp<ConfiguratorType,HyperelasticEnergyDensityType> > {
00938 protected:
00939   typedef typename ConfiguratorType::RealType RealType;
00940 
00941   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
00942   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
00943   // the deformations \phi, \phi_r and \phi_l
00944   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement, _preDisplacement, _postDisplacement;
00945   // the indices of the sub-Hessian
00946   const int _k, _l;
00947 public:
00948   HyperelasticPrePostDeformSubHessianPostComp ( const typename ConfiguratorType::InitType &Grid,
00949                                                 const HyperelasticEnergyDensityType &HyperelasticEnergyDensity,
00950                                                 const aol::MultiVector<RealType> &Displacement,
00951                                                 const aol::MultiVector<RealType> &PreDisplacement,
00952                                                 const aol::MultiVector<RealType> &PostDisplacement,
00953                                                 const int K,
00954                                                 const int L ) :
00955     qc::FELinDeformAsymMatrixWeightedStiffInterface<ConfiguratorType, HyperelasticPrePostDeformSubHessianPostComp<ConfiguratorType,HyperelasticEnergyDensityType> > ( Grid, Displacement ),
00956     _hyperelasticEnergyDensity( HyperelasticEnergyDensity ),
00957     _displacement( Grid, Displacement ),
00958     _preDisplacement( Grid, PreDisplacement ),
00959     _postDisplacement( Grid, PostDisplacement ),
00960     _k( K ),
00961     _l( L ){}
00962 
00963   inline void getCoeffMatrix ( const typename ConfiguratorType::ElementType &El, int QuadPoint, const typename ConfiguratorType::DomVecType &LocalCoord,
00964                                const typename ConfiguratorType::ElementType &TransformedEl, const typename ConfiguratorType::DomVecType &TransformedLocalCoord,
00965                                typename ConfiguratorType::MatType &Matrix ) const {
00966     // compute the deformation gradient $\nabla\Phi$
00967     typename ConfiguratorType::MatType dPhi;
00968     _displacement.evaluateGradientAtQuadPoint( El, QuadPoint, dPhi );
00969     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00970       dPhi[i][i] += 1.;
00971 
00972     // compute the deformation gradient $\nabla\phi_l$ and $\nabla\phi_r$
00973     typename ConfiguratorType::MatType dPhiL, dPhiR;
00974     _preDisplacement.evaluateGradient( El, LocalCoord, dPhiR );
00975     _postDisplacement.evaluateGradient( TransformedEl, TransformedLocalCoord, dPhiL );
00976     for ( int i = 0; i < ConfiguratorType::Dim; i++ ){
00977       dPhiL[i][i] += 1.;
00978       dPhiR[i][i] += 1.;
00979     }
00980 
00981     // compute the Hessian matrix for the kth component multiplied from the left and the lth from the right
00982     typename ConfiguratorType::MatType dPhiRInv, cofDPhiR, aux;
00983     cofDPhiR.makeCofactorMatrix( dPhiR );
00984     dPhiRInv = dPhiR.inverse();
00985     dPhiL *= dPhi;
00986     dPhiL *= dPhiRInv;
00987     HyperelasticSubHessian<ConfiguratorType,HyperelasticEnergyDensityType>::elasticityTensor( dPhiL, _hyperelasticEnergyDensity, El, QuadPoint, _k, _l, Matrix );
00988     aux.makeProduct( Matrix, cofDPhiR );
00989     Matrix.makeProduct( dPhiRInv, aux );
00990   }
00991 };
00992 
00993 /**
00994   * \brief This class computes via "apply(...)" or "applyAdd(...)" part of the hyperelastic energy Hessian with respect to one of two displacements,
00995   * i.e. \f$ (<\frac{\partial^2\int_\Omega W(I_1,I_2,I_3,x)det(\nabla\phi_r) dx}{\partial \phi_r^2},\psi_i,\psi_j>)_{ij} \f$,
00996   * where the \f$ \psi_i \f$ are the vector-valued FE basis functions,
00997   * the hyperelastic invariants are given by \f$ I_1=||\nabla\phi_l\nabla\phi\nabla\phi_r^{-1}||_F^2 \f$, \f$ I_2=||cof(\nabla\phi_l\nabla\phi\nabla\phi_r^{-1})||_F^2 \f$,
00998   * \f$ I_1=det(\nabla\phi_l\nabla\phi\nabla\phi_r^{-1}) \f$ and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
00999   * \f$ \phi \f$ is some prestressing deformation.
01000   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$ of \f$ \phi_l \f$ and \f$ \phi_r \f$
01001   * in the different Cartesian directions.
01002   *
01003   * \author Wirth
01004   */
01005 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
01006 class HyperelasticPrePostDeformSubHessianPreComp
01007       : public aol::FELinAsymMatrixWeightedStiffInterface<ConfiguratorType, HyperelasticPrePostDeformSubHessianPreComp<ConfiguratorType,HyperelasticEnergyDensityType> > {
01008 protected:
01009   typedef typename ConfiguratorType::RealType RealType;
01010 
01011   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
01012   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
01013   // the deformations \phi, \phi_r and \phi_l
01014   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement, _preDisplacement, _postDisplacement;
01015   // the indices of the sub-Hessian
01016   const int _k, _l;
01017 public:
01018   HyperelasticPrePostDeformSubHessianPreComp( const typename ConfiguratorType::InitType &Grid,
01019                                               const HyperelasticEnergyDensityType &HyperelasticEnergyDensity,
01020                                               const aol::MultiVector<RealType> &Displacement,
01021                                               const aol::MultiVector<RealType> &PreDisplacement,
01022                                               const aol::MultiVector<RealType> &PostDisplacement,
01023                                               const int K,
01024                                               const int L ) :
01025     aol::FELinAsymMatrixWeightedStiffInterface<ConfiguratorType, HyperelasticPrePostDeformSubHessianPreComp<ConfiguratorType,HyperelasticEnergyDensityType> > ( Grid ),
01026     _hyperelasticEnergyDensity( HyperelasticEnergyDensity ),
01027     _displacement( Grid, Displacement ),
01028     _preDisplacement( Grid, PreDisplacement ),
01029     _postDisplacement( Grid, PostDisplacement ),
01030     _k( L ), // l and k are exchanged here, since FELinAsymMatrixStiffOp assembles the matrix in the way that \int \nabla\phi^T A \nabla\psi dx becomes \Psi^T L \Phi
01031     _l( K ) {}
01032 
01033   inline void getCoeffMatrix ( const typename ConfiguratorType::ElementType &El, int QuadPoint,
01034                                const typename ConfiguratorType::DomVecType& DomRefCoord, typename ConfiguratorType::MatType &Matrix ) const {
01035     // compute the deformation gradient $\nabla\phi$
01036     typename ConfiguratorType::MatType dPhi;
01037     _displacement.evaluateGradientAtQuadPoint( El, QuadPoint, dPhi );
01038     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
01039       dPhi[i][i] += 1.;
01040 
01041     // compute $\phi(x)$
01042     typename ConfiguratorType::ElementType transformedEl;
01043     typename ConfiguratorType::DomVecType  transformedCoord;
01044     qc::transformAndClipCoord<ConfiguratorType>( *this->_config, _displacement, El, QuadPoint, DomRefCoord, transformedEl, transformedCoord );
01045 
01046     // compute the deformation gradient $\nabla\phi_l$ and $\nabla\phi_r$
01047     typename ConfiguratorType::MatType dPhiL, dPhiR;
01048     _postDisplacement.evaluateGradient( transformedEl, transformedCoord, dPhiL );
01049     _preDisplacement.evaluateGradient( El, DomRefCoord, dPhiR );
01050     for ( int i = 0; i < ConfiguratorType::Dim; i++ ){
01051       dPhiL[i][i] += 1.;
01052       dPhiR[i][i] += 1.;
01053     }
01054 
01055     // compute D\phi_l D\phi D\phi_r^{-1} and \partial_{,A}W
01056     typename ConfiguratorType::MatType defGrad( dPhiL ), cofDPhiR, invDPhiR( dPhiR.inverse() ), dPhiInvDPhiR( dPhi ), dW;
01057     cofDPhiR.makeCofactorMatrix( dPhiR );
01058     dPhiInvDPhiR *= invDPhiR;
01059     defGrad *= dPhiInvDPhiR;
01060     HyperelasticGradient<ConfiguratorType,HyperelasticEnergyDensityType>::firstPiolaKirchhoffStress( defGrad, _hyperelasticEnergyDensity, El, QuadPoint, dW );
01061 
01062     // compute the Hessian matrix; denote variations of \phi_r by \psi and \chi; they only act on the component "_k" and "_l"
01063     // \psi shall belong to vectors multiplied from the left, \chi from the right
01064     // use tr(A D\psi B D\chi) = ( D\psi_{k:} B_{:l} )( D\chi_{l:} A_{:k} ) = D\psi_{k:} ( B_{:l} A_{:k}^T ) D\chi_{l:}
01065     typename ConfiguratorType::MatType aux1, aux2;
01066     typename ConfiguratorType::VecType vec1;
01067     // compute M s.t. (D\psi)_{k}M(D\chi)_{l}^T = [-\partial_{,A}W:(D\phi_l D\phi D\phi_r^{-1} D\psi D\phi_r^{-1})][cof(D\phi_r):D\chi]
01068     aux1.makeProductAtransposedB( defGrad, dW );
01069     aux2.makeProductABtransposed( aux1, invDPhiR );
01070     Matrix.makeTensorProduct( aux2[_k], cofDPhiR[_l] );
01071     Matrix *= -1.;
01072     // compute M s.t. (D\psi)_{k}M(D\chi)_{l}^T = [-\partial_{,A}W:(D\phi_l D\phi D\phi_r^{-1} D\chi D\phi_r^{-1})][cof(D\phi_r):D\psi]
01073     aux1.makeTensorProduct( cofDPhiR[_k], aux2[_l] );
01074     Matrix -= aux1;
01075     // compute M s.t. (D\psi)_{k}M(D\chi)_{l}^T = W<\partial_{,AA}det D\phi_r,D\psi,D\chi>
01076     //   = [(cof D\phi_r;D\psi)(cof D\phi_r;D\chi)-tr(cof D\phi_r D\psi^T cof D\phi_r D\chi^T)] / det D\phi_r
01077     aux1.makeTensorProduct( cofDPhiR[_k], cofDPhiR[_l] );
01078     aux2.makeTensorProduct( cofDPhiR[_l], cofDPhiR[_k] );
01079     aux1 -= aux2;
01080     Matrix.addMultiple( aux1, HyperelasticEnergy<ConfiguratorType,HyperelasticEnergyDensityType>::energyDensity( defGrad, _hyperelasticEnergyDensity, El, QuadPoint ) / dPhiR.det() );
01081     // compute M s.t. (D\psi)_{k}M(D\chi)_{l}^T = det D\phi_r<\partial_{,AA}W,D\phi_l D\phi D\phi_r^{-1} D\psi D\phi_r^{-1},D\phi_l D\phi D\phi_r^{-1} D\chi D\phi_r^{-1}>
01082     aux1.setZero();
01083     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
01084       for ( int j = 0; j < ConfiguratorType::Dim; j++ ) {
01085         HyperelasticSubHessian<ConfiguratorType,HyperelasticEnergyDensityType>::elasticityTensor( defGrad, _hyperelasticEnergyDensity, El, QuadPoint, i, j, aux2 );
01086         aux1.addMultiple( aux2, defGrad[j][_l] * defGrad[i][_k] );
01087       }
01088     aux1 *= cofDPhiR;
01089     aux2.makeProduct( invDPhiR, aux1 );
01090     Matrix += aux2;
01091     // compute M s.t. (D\psi)_{k}M(D\chi)_{l}^T = det D\phi_r[\partial_{,A}W:(D\phi_l D\phi D\phi_r^{-1} [D\psi D\phi_r^{-1} D\chi + D\chi D\phi_r^{-1} D\psi] D\phi_r^{-1})]
01092     aux1.makeProductAtransposedB( defGrad, dW );
01093     aux1 *= cofDPhiR;
01094     invDPhiR.getColumn( _l, vec1 );
01095     aux2.makeTensorProduct( vec1, aux1[_k] );
01096     Matrix += aux2;
01097     invDPhiR.getColumn( _k, vec1 );
01098     aux2.makeTensorProduct( aux1[_l], vec1 );
01099     Matrix += aux2;
01100   }
01101 };
01102 
01103 /**
01104   * \brief This class computes via "apply(...)" or "applyAdd(...)" part of the hyperelastic energy Hessian with respect to two displacements,
01105   * i.e. \f$ (<\frac{\partial^2\int_\Omega W(I_1,I_2,I_3,x)det(\nabla\phi_r) dx}{\partial\phi_l\partial\phi_r},\psi_i,\psi_j>)_{ij} \f$,
01106   * where the \f$ \psi_i \f$ are the vector-valued FE basis functions,
01107   * the hyperelastic invariants are given by \f$ I_1=||\nabla\phi_l\nabla\phi\nabla\phi_r^{-1}||_F^2 \f$, \f$ I_2=||cof(\nabla\phi_l\nabla\phi\nabla\phi_r^{-1})||_F^2 \f$,
01108   * \f$ I_1=det(\nabla\phi_l\nabla\phi\nabla\phi_r^{-1}) \f$ and the function \f$ W(I_1,I_2,I_3,x) \f$ is passed to the constructor (e.g. HyerelasticEnergyDensityDefault).
01109   * \f$ \phi \f$ is some prestressing deformation.
01110   * In its first argument, "apply" and "applyAdd" expect a Multivector which represents the displacements \f$ d_i \f$ of \f$ \phi_l \f$ and \f$ \phi_r \f$
01111   * in the different Cartesian directions.
01112   * If ``PhiRVariationFromRight''=true, then the resulting matrix is the Hessian with respect to a variation of the ``L''th component
01113   * of \f$ \phi_r \f$ from the right and a vatiation of the ``K''th component of \f$ \phi_l \f$ from the left;
01114   * else \f$ \phi_l \f$ and \f$ \phi_r \f$ are exchanged.
01115   *
01116   * \author Wirth
01117   */
01118 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
01119 class HyperelasticPrePostDeformSubHessianMixedPart
01120       : public qc::FELinDeformAsymMatrixWeightedStiffInterface<ConfiguratorType, HyperelasticPrePostDeformSubHessianMixedPart<ConfiguratorType,HyperelasticEnergyDensityType> > {
01121 protected:
01122   typedef typename ConfiguratorType::RealType RealType;
01123   typedef typename ConfiguratorType::MatType    MatType;
01124   typedef typename ConfiguratorType::DomVecType DomVecType;
01125 
01126   // the hyperelastic energy density $W$ as function of the three deformation gradient invariants
01127   const HyperelasticEnergyDensityType &_hyperelasticEnergyDensity;
01128   // the deformations \phi, \phi_r and \phi_l
01129   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement, _preDisplacement, _postDisplacement;
01130   // the indices of the sub-Hessian
01131   const int _k, _l;
01132   // flag indicating whether the test direction for phi_r is multiplied from the right or the left
01133   const bool _phiRVariationFromRight;
01134 public:
01135   HyperelasticPrePostDeformSubHessianMixedPart( const typename ConfiguratorType::InitType &Grid,
01136                                                 const HyperelasticEnergyDensityType &HyperelasticEnergyDensity,
01137                                                 const aol::MultiVector<RealType> &Displacement,
01138                                                 const aol::MultiVector<RealType> &PreDisplacement,
01139                                                 const aol::MultiVector<RealType> &PostDisplacement,
01140                                                 const int K,
01141                                                 const int L,
01142                                                 const bool PhiRVariationFromRight ) :
01143     qc::FELinDeformAsymMatrixWeightedStiffInterface<ConfiguratorType, HyperelasticPrePostDeformSubHessianMixedPart<ConfiguratorType,HyperelasticEnergyDensityType> > ( Grid, Displacement, aol::ASSEMBLED, PhiRVariationFromRight ? qc::LEFT : qc::RIGHT ),
01144     _hyperelasticEnergyDensity( HyperelasticEnergyDensity ),
01145     _displacement( Grid, Displacement ),
01146     _preDisplacement( Grid, PreDisplacement ),
01147     _postDisplacement( Grid, PostDisplacement ),
01148     _k( K ),
01149     _l( L ),
01150     _phiRVariationFromRight( PhiRVariationFromRight ){}
01151 
01152   inline void getCoeffMatrix ( const typename ConfiguratorType::ElementType &El, int QuadPoint, const DomVecType &LocalCoord,
01153                                const typename ConfiguratorType::ElementType &TransformedEl, const DomVecType &TransformedLocalCoord,
01154                                MatType &Matrix ) const {
01155     // compute the deformation gradient $\nabla\phi$
01156     typename ConfiguratorType::MatType dPhi;
01157     _displacement.evaluateGradientAtQuadPoint( El, QuadPoint, dPhi );
01158     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
01159       dPhi[i][i] += 1.;
01160 
01161     // compute the deformation gradient $\nabla\phi_l$ and $\nabla\phi_r$
01162     typename ConfiguratorType::MatType dPhiL, dPhiR;
01163     _postDisplacement.evaluateGradient( TransformedEl, TransformedLocalCoord, dPhiL );
01164     _preDisplacement.evaluateGradient( El, LocalCoord, dPhiR );
01165     for ( int i = 0; i < ConfiguratorType::Dim; i++ ){
01166       dPhiL[i][i] += 1.;
01167       dPhiR[i][i] += 1.;
01168     }
01169 
01170     // compute D\phi_l D\phi D\phi_r^{-1} and \partial_{,A}W
01171     typename ConfiguratorType::MatType defGrad( dPhiL ), cofDPhiR, invDPhiR( dPhiR.inverse() ), dPhiInvDPhiR( dPhi ), dW;
01172     cofDPhiR.makeCofactorMatrix( dPhiR );
01173     dPhiInvDPhiR *= invDPhiR;
01174     defGrad *= dPhiInvDPhiR;
01175     HyperelasticGradient<ConfiguratorType,HyperelasticEnergyDensityType>::firstPiolaKirchhoffStress( defGrad, _hyperelasticEnergyDensity, El, QuadPoint, dW );
01176 
01177     // compute the Hessian matrix, pretending the gradient of the variation of \phi_r is multiplied from the right
01178     // denote variations of \phi_r and \phi_l by \psi_r and \psi_l; these variations only act on the component "compPhiR" and "compPhiL"
01179     int compPhiR = _phiRVariationFromRight ? _l : _k,
01180         compPhiL = _phiRVariationFromRight ? _k : _l;
01181     typename ConfiguratorType::MatType aux1, aux2;
01182     typename ConfiguratorType::VecType vec1;
01183     // compute M s.t. (D\psi_l)_{"compPhiL"}M(D\psi_r)_{"compPhiR"}^T  = [\partial_{,A}W:(D\psi_l D\phi D\phi_r^{-1})][cof(D\phi_r):D\psi_r]
01184     aux1.makeProductABtransposed( dW, dPhiInvDPhiR );
01185     Matrix.makeTensorProduct( aux1[compPhiL], cofDPhiR[compPhiR] );
01186     // compute M s.t. (D\psi_l)_{"compPhiL"}M(D\psi_r)_{"compPhiR"}^T  = -det(D\phi_r)[\partial_{,A}W:(D\psi_l D\phi D\phi_r^{-1} D\psi_r D\phi_r^{-1})]
01187     aux1.makeProduct( dW, cofDPhiR );
01188     dPhiInvDPhiR.getColumn( compPhiR, vec1 );
01189     aux2.makeTensorProduct( vec1, aux1[compPhiL] );
01190     Matrix -= aux2;
01191     // compute M s.t. (D\psi_l)_{"compPhiL"}M(D\psi_r)_{"compPhiR"}^T  = -det(D\phi_r)[\partial_{,AA}W(D\psi_l D\phi D\phi_r^{-1})]:(D\phi_l D\phi D\phi_r^{-1} D\psi_r D\phi_r^{-1})
01192     aux1.setZero();
01193     for ( int j = 0; j < ConfiguratorType::Dim; j++ ) {
01194       HyperelasticSubHessian<ConfiguratorType,HyperelasticEnergyDensityType>::elasticityTensor( defGrad, _hyperelasticEnergyDensity, El, QuadPoint, compPhiL, j, aux2 );
01195       aux1.addMultiple( aux2, defGrad[j][compPhiR] );
01196     }
01197     aux2.makeProduct( dPhiInvDPhiR, aux1 );
01198     aux2 *= cofDPhiR;
01199     Matrix -= aux2;
01200     // adapt local matrix if gradient of the variation of \phi_r is multiplied from the left
01201     if ( !_phiRVariationFromRight )
01202       Matrix.transpose();
01203   }
01204 };
01205 
01206 /**
01207   * \brief This class represents a standard implementation of the hyperelastic energy density \f$ W(I_1,I_2,I_3) \f$
01208   * with respect to the invariants \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2 \f$, \f$ I_3=det(\nabla\phi) \f$
01209   * of the deformation gradient \f$\nabla\phi\f$. In addition to \f$I_1,I_2,I_3\f$, the methods
01210   * "evaluate(...)" and "evaluateDerivative(...)" have to be passed a spatial position, since
01211   * in principle the energy density may be space dependent.
01212   *
01213   * \author Wirth
01214   */
01215 template <typename ConfiguratorType>
01216 class HyperelasticEnergyDensityDefault{
01217 protected:
01218   typedef typename ConfiguratorType::RealType RealType;
01219   // weights of the different energy components
01220   const RealType _weightLengthEnergy;
01221   const RealType _weightSurfaceEnergy;
01222   const RealType _weightVolumeEnergy;
01223   // Maximum compression allowed. If this compression is reached, energy is set to a high value.
01224   // If the threshold I3t is reached, a gradient descent step with minimum step length should still be Armijo-rule-
01225   // compatible, i.e. in this case ln(I3)-ln(I3t)>[1/I3t*(I3-I3t)]/2 with I3=I3t+\tau_{min}*1/I3t, hence I3t\circ1e-5.
01226   const RealType _compressionThreshold;
01227   const RealType _valueAtThreshold;
01228 public:
01229   HyperelasticEnergyDensityDefault( const RealType WeightLengthEnergy, const RealType WeightSurfaceEnergy, const RealType WeightVolumeEnergy ) :
01230     _weightLengthEnergy ( WeightLengthEnergy ),
01231     _weightSurfaceEnergy ( WeightSurfaceEnergy ),
01232     _weightVolumeEnergy ( WeightVolumeEnergy ),
01233     _compressionThreshold ( 1.e-10 ),
01234     _valueAtThreshold ( aol::NumberTrait<RealType>::Inf ) {}
01235 
01236   //! returns \f$\frac\mu2I_1+\frac\lambda4I_3^2-(\mu+\frac\lambda2)\log I_3-\frac\mu2d\frac\lambda4\f$ for parameters \f$\mu,\lambda\f$, dimension \f$d\f$
01237   inline RealType evaluate ( const RealType I1, const RealType /*I2*/, const RealType I3,
01238                              const typename ConfiguratorType::ElementType &/*El*/, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01239     return _weightLengthEnergy * I1 / 2
01240            + _weightVolumeEnergy * aol::Sqr( I3 ) / 4
01241            - ( _weightLengthEnergy + _weightVolumeEnergy / 2 ) * ( I3 <= _compressionThreshold ? log( _compressionThreshold ) + ( I3 - _compressionThreshold ) / _compressionThreshold  -_valueAtThreshold : log( I3 ) )
01242            - ( ConfiguratorType::Dim / 2. ) * _weightLengthEnergy - _weightVolumeEnergy / 4;
01243   }
01244 
01245   //! returns \f$\frac\mu2I_1+\frac\lambda4I_3^2-(\mu+\frac\lambda2)\log I_3-\frac\mu2d\frac\lambda4\f$ for parameters \f$\mu,\lambda\f$, dimension \f$d\f$
01246   inline RealType evaluateAtQuadPoint ( const RealType I1, const RealType /*I2*/, const RealType I3,
01247                                         const typename ConfiguratorType::ElementType &/*El*/, int /*QuadPoint*/ ) const {
01248     return _weightLengthEnergy * I1 / 2
01249            + _weightVolumeEnergy * aol::Sqr( I3 ) / 4
01250            - ( _weightLengthEnergy + _weightVolumeEnergy / 2 ) * ( I3 <= _compressionThreshold ? log( _compressionThreshold ) + ( I3 - _compressionThreshold ) / _compressionThreshold  -_valueAtThreshold : log( I3 ) )
01251            - ( ConfiguratorType::Dim / 2. ) * _weightLengthEnergy - _weightVolumeEnergy / 4;
01252   }
01253 
01254   //! returns \f$(\frac\mu2,0,\frac\lambda2I_3-(\mu+\frac\lambda2)\frac1{I_3})
01255   inline aol::Vec3<RealType> evaluateDerivative ( const RealType /*I1*/, const RealType /*I2*/, const RealType I3,
01256                                                   const typename ConfiguratorType::ElementType &/*El*/, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01257     return aol::Vec3<RealType>( _weightLengthEnergy / 2,
01258                                 0,
01259                                 _weightVolumeEnergy * I3 / 2 - ( _weightLengthEnergy + _weightVolumeEnergy / 2 ) / aol::Max( I3, _compressionThreshold ) );
01260   }
01261 
01262   //! returns \f$(\frac\mu2,0,\frac\lambda2I_3-(\mu+\frac\lambda2)\frac1{I_3})
01263   inline aol::Vec3<RealType> evaluateDerivativeAtQuadPoint ( const RealType /*I1*/, const RealType /*I2*/, const RealType I3,
01264                                                              const typename ConfiguratorType::ElementType &/*El*/, int /*QuadPoint*/ ) const {
01265     return aol::Vec3<RealType>( _weightLengthEnergy / 2,
01266                                 0,
01267                                 _weightVolumeEnergy * I3 / 2 - ( _weightLengthEnergy + _weightVolumeEnergy / 2 ) / aol::Max( I3, _compressionThreshold ) );
01268   }
01269 
01270   //! returns diag\f$(0,0,\frac\lambda2+(\mu+\frac\lambda2)\frac1{I_3^2})
01271   inline aol::Matrix33<RealType> evaluateSecondDerivative ( const RealType /*I1*/, const RealType /*I2*/, const RealType I3,
01272                                                             const typename ConfiguratorType::ElementType &/*El*/, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01273     return aol::Matrix33<RealType>( 0, 0, 0, 0, 0, 0, 0, 0,
01274                                     _weightVolumeEnergy / 2 + ( I3 <= _compressionThreshold ? 0. : ( _weightLengthEnergy + _weightVolumeEnergy / 2 ) / aol::Sqr( I3 ) ) );
01275   }
01276 
01277   //! returns diag\f$(0,0,\frac\lambda2+(\mu+\frac\lambda2)\frac1{I_3^2})
01278   inline aol::Matrix33<RealType> evaluateSecondDerivativeAtQuadPoint ( const RealType /*I1*/, const RealType /*I2*/, const RealType I3,
01279                                                                        const typename ConfiguratorType::ElementType &/*El*/, int /*QuadPoint*/ ) const {
01280     return aol::Matrix33<RealType>( 0, 0, 0, 0, 0, 0, 0, 0,
01281                                     _weightVolumeEnergy / 2 + ( I3 <= _compressionThreshold ? 0. : ( _weightLengthEnergy + _weightVolumeEnergy / 2 ) / aol::Sqr( I3 ) ) );
01282   }
01283 };
01284 
01285 /**
01286   * \brief Implements the St Venant-Kirchhoff energy density \f$ W(I_1,I_2,I_3)=\frac\lambda2(I_1-d)^2+\frac\mu4((I_1-1)^2+d-1-2\zeta)+\nu\Gamma(I_3) \f$
01287   * with \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2 \f$, \f$ I_3=det(\nabla\phi) \f$, \f$ \zeta=I_3^2 \f$ in 2D and \f$ \zeta=I_2 \f$ in 3D,
01288   * \f$ \Gamma(I_3)=\frac1{I_3^2}+I_3^2-2 \f$, and the deformation gradient \f$\nabla\phi\f$, where \f$ \Gamma \f$ in fact does not belong to the law and
01289   * only prohibits complete compression.
01290   *
01291   * \author Wirth
01292   */
01293 template <typename ConfiguratorType>
01294 class HyperelasticEnergyDensityStVenantKirchhoff{
01295 protected:
01296   typedef typename ConfiguratorType::RealType RealType;
01297   // Lame constants and weight of term preventing complete volume compression (which in fact does not belong to St Venant-Kirchhoff)
01298   const RealType _lambdaQuarter, _lambdaEighth;
01299   const RealType _mu, _muQuarter;
01300   const RealType _weightAntiCompression;
01301   const RealType _dimMinusOne;
01302 public:
01303   HyperelasticEnergyDensityStVenantKirchhoff( const RealType Lambda, const RealType Mu, const RealType WeightAntiCompression ) :
01304     _lambdaQuarter( Lambda / 4 ),
01305     _lambdaEighth( Lambda / 8 ),
01306     _mu( Mu ),
01307     _muQuarter( Mu / 4 ),
01308     _weightAntiCompression( WeightAntiCompression ),
01309     _dimMinusOne( ConfiguratorType::Dim - 1 ) {}
01310 
01311   inline RealType evaluate ( const RealType I1, const RealType I2, const RealType I3,
01312                              const typename ConfiguratorType::ElementType &/*El*/, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01313     return _lambdaEighth * aol::Sqr( I1 - ConfiguratorType::Dim )
01314            + _muQuarter * ( aol::Sqr( I1 - 1 ) + _dimMinusOne - 2 * ( ConfiguratorType::Dim == qc::QC_2D ? aol::Sqr( I3 ) : I2 ) )
01315            + _weightAntiCompression * ( 1 / aol::Sqr( I3 ) + aol::Sqr( I3 ) - 2 );
01316   }
01317 
01318   inline RealType evaluateAtQuadPoint ( const RealType I1, const RealType I2, const RealType I3,
01319                                         const typename ConfiguratorType::ElementType &/*El*/, int /*QuadPoint*/ ) const {
01320     return _lambdaEighth * aol::Sqr( I1 - ConfiguratorType::Dim )
01321            + _muQuarter * ( aol::Sqr( I1 - 1 ) + _dimMinusOne - 2 * ( ConfiguratorType::Dim == qc::QC_2D ? aol::Sqr( I3 ) : I2 ) )
01322            + _weightAntiCompression * ( 1 / aol::Sqr( I3 ) + aol::Sqr( I3 ) - 2 );
01323   }
01324 
01325   inline aol::Vec3<RealType> evaluateDerivative ( const RealType I1, const RealType /*I2*/, const RealType I3,
01326                                                   const typename ConfiguratorType::ElementType &/*El*/, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01327     return aol::Vec3<RealType>( _lambdaQuarter * ( I1 - ConfiguratorType::Dim ) + _mu / 2 * ( I1 - 1 ),
01328                                 ( ConfiguratorType::Dim == qc::QC_2D ? 0 : - _mu / 2 ),
01329                                 ( ConfiguratorType::Dim == qc::QC_2D ? - _mu * I3 : 0 ) + _weightAntiCompression * 2 * ( I3 - 1 / aol::Cub( I3 ) ) );
01330   }
01331 
01332   inline aol::Vec3<RealType> evaluateDerivativeAtQuadPoint ( const RealType I1, const RealType /*I2*/, const RealType I3,
01333                                                              const typename ConfiguratorType::ElementType &/*El*/, int /*QuadPoint*/ ) const {
01334     return aol::Vec3<RealType>( _lambdaQuarter * ( I1 - ConfiguratorType::Dim ) + _mu / 2 * ( I1 - 1 ),
01335                                 ( ConfiguratorType::Dim == qc::QC_2D ? 0 : - _mu / 2 ),
01336                                 ( ConfiguratorType::Dim == qc::QC_2D ? - _mu * I3 : 0 ) + _weightAntiCompression * 2 * ( I3 - 1 / aol::Cub( I3 ) ) );
01337   }
01338 
01339   inline aol::Matrix33<RealType> evaluateSecondDerivative ( const RealType /*I1*/, const RealType /*I2*/, const RealType /*I3*/,
01340                                                             const typename ConfiguratorType::ElementType &/*El*/, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01341     return aol::Matrix33<RealType>( _lambdaQuarter + _mu / 2, 0, 0, 0, 0, 0, 0, 0, ( ConfiguratorType::Dim == qc::QC_2D ? - _mu : 0 ) );
01342   }
01343 
01344   inline aol::Matrix33<RealType> evaluateSecondDerivativeAtQuadPoint ( const RealType /*I1*/, const RealType /*I2*/, const RealType /*I3*/,
01345                                                                        const typename ConfiguratorType::ElementType &/*El*/, int /*QuadPoint*/ ) const {
01346     return aol::Matrix33<RealType>( _lambdaQuarter + _mu / 2, 0, 0, 0, 0, 0, 0, 0, ( ConfiguratorType::Dim == qc::QC_2D ? - _mu : 0 ) );
01347   }
01348 };
01349 
01350 /**
01351   * \brief Implements the 2D linear energy density \f$ W(I_1,I_2,I_3)=\alpha(I_1+2-2\sqrt{I_1+2\det I_3})=\alpha dist(\nabla\phi,SO(2))^2 \f$
01352   * with \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2 \f$, \f$ I_3=det(\nabla\phi) \f$
01353   *
01354   * \author Wirth
01355   */
01356 template <typename ConfiguratorType>
01357 class HyperelasticEnergyDensity2DLinear{
01358 protected:
01359   typedef typename ConfiguratorType::RealType RealType;
01360   const RealType _alpha;
01361 
01362 public:
01363   HyperelasticEnergyDensity2DLinear( const RealType Alpha ) :
01364     _alpha( Alpha ) {
01365     if ( ConfiguratorType::Dim != qc::QC_2D )
01366       throw aol::Exception ( "qc::HyperelasticEnergyDensity2DLinear currently only implemented for 2D", __FILE__, __LINE__ );
01367   }
01368 
01369   inline RealType evaluate ( const RealType I1, const RealType /*I2*/, const RealType I3,
01370                              const typename ConfiguratorType::ElementType &/*El*/, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01371     return _alpha * ( I1 + 2 - 2 * sqrt( I1 + 2 * I3 ) );
01372   }
01373 
01374   inline RealType evaluateAtQuadPoint ( const RealType I1, const RealType /*I2*/, const RealType I3,
01375                                         const typename ConfiguratorType::ElementType &/*El*/, int /*QuadPoint*/ ) const {
01376     return _alpha * ( I1 + 2 - 2 * sqrt( I1 + 2 * I3 ) );
01377   }
01378 
01379   inline aol::Vec3<RealType> evaluateDerivative ( const RealType I1, const RealType /*I2*/, const RealType I3,
01380                                                   const typename ConfiguratorType::ElementType &/*El*/, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01381     RealType auxTerm = - 1. / sqrt( I1 + 2 * I3 );
01382     return aol::Vec3<RealType>( _alpha * ( 1 + auxTerm ), 0, _alpha * 2 * auxTerm );
01383   }
01384 
01385   inline aol::Vec3<RealType> evaluateDerivativeAtQuadPoint ( const RealType I1, const RealType /*I2*/, const RealType I3,
01386                                                              const typename ConfiguratorType::ElementType &/*El*/, int /*QuadPoint*/ ) const {
01387     RealType auxTerm = - 1. / sqrt( I1 + 2 * I3 );
01388     return aol::Vec3<RealType>( _alpha * ( 1 + auxTerm ), 0, _alpha * 2 * auxTerm );
01389   }
01390 
01391   inline aol::Matrix33<RealType> evaluateSecondDerivative ( const RealType I1, const RealType /*I2*/, const RealType I3,
01392                                                             const typename ConfiguratorType::ElementType &/*El*/, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01393     RealType auxTerm = aol::Cub( 1. / sqrt( I1 + 2 * I3 ) );
01394     return aol::Matrix33<RealType>( _alpha / 2. * auxTerm, 0, _alpha * auxTerm,
01395                                     0, 0, 0,
01396                                     _alpha * auxTerm, 0, _alpha * 2 * auxTerm );
01397   }
01398 
01399   inline aol::Matrix33<RealType> evaluateSecondDerivativeAtQuadPoint ( const RealType I1, const RealType /*I2*/, const RealType I3,
01400                                                                        const typename ConfiguratorType::ElementType &/*El*/, int /*QuadPoint*/ ) const {
01401     RealType auxTerm = aol::Cub( 1. / sqrt( I1 + 2 * I3 ) );
01402     return aol::Matrix33<RealType>( _alpha / 2. * auxTerm, 0, _alpha * auxTerm,
01403                                     0, 0, 0,
01404                                     _alpha * auxTerm, 0, _alpha * 2 * auxTerm );
01405   }
01406 };
01407 
01408 /**
01409   * \brief Implements a generalized Mooney-Rivlin energy density \f$ W(I_1,I_2,I_3)=aI_1^{\frac{p}2}+bI_2^{\frac{q}2}+c\Gamma(I_3)-ad^{\frac{p}2}-bd^{\frac{q}2} \f$
01410   * with \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2 \f$, \f$ I_3=det(\nabla\phi) \f$,
01411   * \f$ \Gamma(I_3)=\frac1{I_3^s}+\frac{s}rI_3^r-1-\frac{s}r \f$, and the deformation gradient \f$\nabla\phi\f$.
01412   *
01413   * \author Wirth
01414   */
01415 template <typename ConfiguratorType>
01416 class HyperelasticEnergyDensityGeneralizedMooneyRivlin{
01417 protected:
01418   typedef typename ConfiguratorType::RealType RealType;
01419   // material parameters
01420   const RealType _pHalf, _qHalf;
01421   const RealType _lengthFac, _areaFac, _volumeFac;
01422   const RealType _r, _s;
01423 public:
01424   HyperelasticEnergyDensityGeneralizedMooneyRivlin( const RealType P, const RealType Q, const RealType LengthFac, const RealType AreaFac, const RealType VolumeFac, const RealType R, const RealType S ) :
01425     _pHalf( P / 2 ),
01426     _qHalf( Q / 2 ),
01427     _lengthFac( LengthFac ),
01428     _areaFac( AreaFac ),
01429     _volumeFac( VolumeFac ),
01430     _r( R ),
01431     _s( S ) {}
01432 
01433   inline RealType evaluate ( const RealType I1, const RealType I2, const RealType I3,
01434                              const typename ConfiguratorType::ElementType &/*El*/, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01435     return _lengthFac * ( pow( I1, _pHalf ) - pow( ConfiguratorType::Dim, _pHalf ) )
01436            + _areaFac * ( pow( I2, _qHalf ) - pow( ConfiguratorType::Dim, _qHalf ) )
01437            + _volumeFac * ( pow( I3, -_s ) - 1 + ( _s / _r ) * ( pow( I3, _r ) - 1 ) );
01438   }
01439 
01440   inline RealType evaluateAtQuadPoint ( const RealType I1, const RealType I2, const RealType I3,
01441                                         const typename ConfiguratorType::ElementType &/*El*/, int /*QuadPoint*/ ) const {
01442     return _lengthFac * ( pow( I1, _pHalf ) - pow( ConfiguratorType::Dim, _pHalf ) )
01443            + _areaFac * ( pow( I2, _qHalf ) - pow( ConfiguratorType::Dim, _qHalf ) )
01444            + _volumeFac * ( pow( I3, -_s ) - 1 + ( _s / _r ) * ( pow( I3, _r ) - 1 ) );
01445   }
01446 
01447   inline aol::Vec3<RealType> evaluateDerivative ( const RealType I1, const RealType I2, const RealType I3,
01448                                                   const typename ConfiguratorType::ElementType &/*El*/, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01449     return aol::Vec3<RealType>( _lengthFac * _pHalf * pow( I1, _pHalf - 1 ),
01450                                 _areaFac * _qHalf * pow( I2, _qHalf - 1 ),
01451                                 _volumeFac * _s * ( pow( I3, _r - 1 ) - pow( I3, -_s - 1 ) ) );
01452   }
01453 
01454   inline aol::Vec3<RealType> evaluateDerivativeAtQuadPoint ( const RealType I1, const RealType I2, const RealType I3,
01455                                                              const typename ConfiguratorType::ElementType &/*El*/, int /*QuadPoint*/ ) const {
01456     return aol::Vec3<RealType>( _lengthFac * _pHalf * pow( I1, _pHalf - 1 ),
01457                                 _areaFac * _qHalf * pow( I2, _qHalf - 1 ),
01458                                 _volumeFac * _s * ( pow( I3, _r - 1 ) - pow( I3, -_s - 1 ) ) );
01459   }
01460 
01461   inline aol::Matrix33<RealType> evaluateSecondDerivative ( const RealType I1, const RealType I2, const RealType I3,
01462                                                             const typename ConfiguratorType::ElementType &/*El*/, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01463     return aol::Matrix33<RealType>( _lengthFac * _pHalf * ( _pHalf - 1 ) * pow( I1, _pHalf - 2 ), 0, 0,
01464                                     0, _areaFac * _qHalf * ( _qHalf - 1 ) * pow( I2, _qHalf - 2 ), 0,
01465                                     0, 0, _volumeFac * _s * ( ( _r - 1 ) * pow( I3, _r - 2 ) + ( _s + 1 ) * pow( I3, -_s - 2 ) ) );
01466   }
01467 
01468   inline aol::Matrix33<RealType> evaluateSecondDerivativeAtQuadPoint ( const RealType I1, const RealType I2, const RealType I3,
01469                                                                        const typename ConfiguratorType::ElementType &/*El*/, int /*QuadPoint*/ ) const {
01470     return aol::Matrix33<RealType>( _lengthFac * _pHalf * ( _pHalf - 1 ) * pow( I1, _pHalf - 2 ), 0, 0,
01471                                     0, _areaFac * _qHalf * ( _qHalf - 1 ) * pow( I2, _qHalf - 2 ), 0,
01472                                     0, 0, _volumeFac * _s * ( ( _r - 1 ) * pow( I3, _r - 2 ) + ( _s + 1 ) * pow( I3, -_s - 2 ) ) );
01473   }
01474 };
01475 
01476 /**
01477   * \brief This class represents a spatially weighted hyperelastic energy density \f$ w(x)W(I_1,I_2,I_3) \f$
01478   * with respect to the invariants \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2 \f$, \f$ I_3=det(\nabla\phi) \f$
01479   * of the deformation gradient \f$\nabla\phi\f$, where the hyperelastic energy \f$W\f$ and the weight \f$w\f$ have to be
01480   * passed to the constructor. \f$W\f$ can e.g. be "HyperelasticEnergyDensityDefault".
01481   *
01482   * \author Wirth
01483   */
01484 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
01485 class WeightedHyperelasticEnergyDensity {
01486 private:
01487   typedef typename ConfiguratorType::RealType RealType;
01488   const HyperelasticEnergyDensityType &_energy;
01489   const aol::DiscreteFunctionDefault<ConfiguratorType> _weight;
01490 public:
01491   WeightedHyperelasticEnergyDensity( const HyperelasticEnergyDensityType &Energy, const typename ConfiguratorType::InitType &Grid, const aol::Vector<RealType> &Weight ) :
01492     _energy( Energy ),
01493     _weight( Grid, Weight ) {}
01494 
01495   inline RealType evaluate ( const RealType I1, const RealType I2, const RealType I3,
01496                              const typename ConfiguratorType::ElementType &El, const typename ConfiguratorType::DomVecType &RefCoord ) const {
01497     return _weight.evaluate( El, RefCoord ) * _energy.evaluate( I1, I2, I3, El, RefCoord );
01498   }
01499 
01500   inline RealType evaluateAtQuadPoint ( const RealType I1, const RealType I2, const RealType I3,
01501                                         const typename ConfiguratorType::ElementType &El, int QuadPoint ) const {
01502     return _weight.evaluateAtQuadPoint( El, QuadPoint ) * _energy.evaluateAtQuadPoint( I1, I2, I3, El, QuadPoint );
01503   }
01504 
01505   inline aol::Vec3<RealType> evaluateDerivative ( const RealType I1, const RealType I2, const RealType I3,
01506                                                   const typename ConfiguratorType::ElementType &El, const typename ConfiguratorType::DomVecType &RefCoord ) const {
01507     return _weight.evaluate( El, RefCoord ) * _energy.evaluateDerivative( I1, I2, I3, El, RefCoord );
01508   }
01509 
01510   inline aol::Vec3<RealType> evaluateDerivativeAtQuadPoint ( const RealType I1, const RealType I2, const RealType I3,
01511                                                              const typename ConfiguratorType::ElementType &El, int QuadPoint ) const {
01512     return _weight.evaluateAtQuadPoint( El, QuadPoint ) * _energy.evaluateDerivativeAtQuadPoint( I1, I2, I3, El, QuadPoint );
01513   }
01514 
01515   inline aol::Matrix33<RealType> evaluateSecondDerivative ( const RealType I1, const RealType I2, const RealType I3,
01516                                                             const typename ConfiguratorType::ElementType &El, const typename ConfiguratorType::DomVecType &RefCoord ) const {
01517     aol::Matrix33<RealType> result( _energy.evaluateSecondDerivative( I1, I2, I3, El, RefCoord ) );
01518     result *= _weight.evaluateAtQuadPoint( El, RefCoord );
01519     return result;
01520   }
01521 
01522   inline aol::Matrix33<RealType> evaluateSecondDerivativeAtQuadPoint ( const RealType I1, const RealType I2, const RealType I3,
01523                                                                        const typename ConfiguratorType::ElementType &El, int QuadPoint ) const {
01524     aol::Matrix33<RealType> result( _energy.evaluateSecondDerivativeAtQuadPoint( I1, I2, I3, El, QuadPoint ) );
01525     result *= _weight.evaluateAtQuadPoint( El, QuadPoint );
01526     return result;
01527   }
01528 };
01529 
01530 /**
01531   * \brief This class represents the hyperelastic energy density \f$ W(I_1,I_2,I_3) + f(I3) \f$
01532   * with respect to the invariants \f$ I_1=||\nabla\phi||_F^2 \f$, \f$ I_2=||cof(\nabla\phi)||_F^2 \f$, \f$ I_3=det(\nabla\phi) \f$
01533   * of the deformation gradient \f$\nabla\phi\f$, where the hyperelastic energy \f$W\f$ has to be
01534   * passed to the constructor and \f$ f(x)=-\log(3x)+\frac5{12} \f$ for \f$ x<\frac13 \f$,
01535   * \f$ f(x)=-\frac92(\frac92x+1)(x-\frac23)^3 \f$ for \f$ \frac13< x <\frac23 \f$, and
01536   * \f$ f(x)=0 \f$ for \f$ x>\frac23 \f$.
01537   *
01538   * \author Wirth
01539   */
01540 template <typename ConfiguratorType, typename HyperelasticEnergyDensityType>
01541 class AntiCompressionHyperelasticEnergyDensity {
01542 private:
01543   typedef typename ConfiguratorType::RealType RealType;
01544   const HyperelasticEnergyDensityType &_energy;
01545 public:
01546   AntiCompressionHyperelasticEnergyDensity( const HyperelasticEnergyDensityType &Energy ) :
01547     _energy( Energy ) {}
01548 
01549   inline RealType evaluate ( const RealType I1, const RealType I2, const RealType I3,
01550                              const typename ConfiguratorType::ElementType &El, const typename ConfiguratorType::DomVecType &RefCoord ) const {
01551     return _energy.evaluate( I1, I2, I3, El, RefCoord )
01552            + ( I3 < 0 ? aol::NumberTrait<RealType>::Inf
01553            : ( 3 * I3 < 1 ? -log( 3 * I3 ) + 5. / 12.
01554            : ( 3 * I3 < 2 ? -4.5 * ( 4.5 * I3 + 1 ) * aol::Cub( I3 - 2. / 3. ) : 0 ) ) );
01555   }
01556 
01557   inline RealType evaluateAtQuadPoint ( const RealType I1, const RealType I2, const RealType I3,
01558                                         const typename ConfiguratorType::ElementType &El, int QuadPoint ) const {
01559     return _energy.evaluateAtQuadPoint( I1, I2, I3, El, QuadPoint )
01560            + ( I3 < 0 ? aol::NumberTrait<RealType>::Inf
01561            : ( 3 * I3 < 1 ? -log( 3 * I3 ) + 5. / 12.
01562            : ( 3 * I3 < 2 ? -4.5 * ( 4.5 * I3 + 1 ) * aol::Cub( I3 - 2. / 3. ) : 0 ) ) );
01563   }
01564 
01565   inline aol::Vec3<RealType> evaluateDerivative ( const RealType I1, const RealType I2, const RealType I3,
01566                                                   const typename ConfiguratorType::ElementType &El, const typename ConfiguratorType::DomVecType &RefCoord ) const {
01567     aol::Vec3<RealType> result = _energy.evaluateDerivative( I1, I2, I3, El, RefCoord );
01568     result[2] += ( I3 < 0 ? -1
01569                : ( 3 * I3 < 1 ? -1. / I3
01570                : ( 3 * I3 < 2 ? -4.5 * ( 4.5 * aol::Cub( I3 - 2. / 3. ) + 3 * ( 4.5 * I3 + 1 ) * aol::Sqr( I3 - 2. / 3. ) ) : 0 ) ) );
01571     return result;
01572   }
01573 
01574   inline aol::Vec3<RealType> evaluateDerivativeAtQuadPoint ( const RealType I1, const RealType I2, const RealType I3,
01575                                                              const typename ConfiguratorType::ElementType &El, int QuadPoint ) const {
01576     aol::Vec3<RealType> result = _energy.evaluateDerivativeAtQuadPoint( I1, I2, I3, El, QuadPoint );
01577     result[2] += ( I3 < 0 ? -1
01578                : ( 3 * I3 < 1 ? -1. / I3
01579                : ( 3 * I3 < 2 ? -4.5 * ( 4.5 * aol::Cub( I3 - 2. / 3. ) + 3 * ( 4.5 * I3 + 1 ) * aol::Sqr( I3 - 2. / 3. ) ) : 0 ) ) );
01580     return result;
01581   }
01582 
01583   inline aol::Matrix33<RealType> evaluateSecondDerivative ( const RealType I1, const RealType I2, const RealType I3,
01584                                                             const typename ConfiguratorType::ElementType &El, const typename ConfiguratorType::DomVecType &RefCoord ) const {
01585     aol::Matrix33<RealType> result( _energy.evaluateSecondDerivative( I1, I2, I3, El, RefCoord ) );
01586     result[2][2] += ( I3 < 0 ? 0
01587                   : ( 3 * I3 < 1 ? 1. / aol::Sqr( I3 )
01588                   : ( 3 * I3 < 2 ? -4.5 * ( 27 * aol::Sqr( I3 - 2. / 3. ) + 6 * ( 4.5 * I3 + 1 ) * ( I3 - 2. / 3. ) ) : 0 ) ) );
01589     return result;
01590   }
01591 
01592   inline aol::Matrix33<RealType> evaluateSecondDerivativeAtQuadPoint ( const RealType I1, const RealType I2, const RealType I3,
01593                                                                        const typename ConfiguratorType::ElementType &El, int QuadPoint ) const {
01594     aol::Matrix33<RealType> result( _energy.evaluateSecondDerivativeAtQuadPoint( I1, I2, I3, El, QuadPoint ) );
01595     result[2][2] += ( I3 < 0 ? 0
01596                   : ( 3 * I3 < 1 ? 1. / aol::Sqr( I3 )
01597                   : ( 3 * I3 < 2 ? -4.5 * ( 27 * aol::Sqr( I3 - 2. / 3. ) + 6 * ( 4.5 * I3 + 1 ) * ( I3 - 2. / 3. ) ) : 0 ) ) );
01598     return result;
01599   }
01600 };
01601 
01602 
01603 template <typename RealType, typename VecType>
01604 class RegularizedLinCombDescentDirOp : public aol::Op<VecType, VecType> { };
01605 
01606 // specialization
01607 template <class RealType>
01608 class RegularizedLinCombDescentDirOp<RealType, aol::MultiVector<RealType> > : public aol::Op<aol::MultiVector<RealType> > {
01609 protected:
01610   typedef aol::MultiVector<RealType> VecType;
01611   const qc::GridDefinition &_grid;
01612   RealType _smoothSigma;
01613   const aol::Op<VecType, VecType> &_massOpInv;
01614   aol::LinCombOp<VecType> _linCombOpReg;
01615   aol::LinCombOp<VecType> _linCombOpUnreg;
01616 
01617   mutable VecType *_gradient;
01618 
01619 public:
01620   RegularizedLinCombDescentDirOp ( const qc::GridDefinition &grid,
01621                                    RealType smoothSigma,
01622                                    const aol::Op<VecType, VecType> &massOpInv )
01623       : _grid ( grid ), _smoothSigma ( smoothSigma ), _massOpInv ( massOpInv ), _gradient ( NULL ) {}
01624 
01625   virtual ~RegularizedLinCombDescentDirOp() {
01626     delete _gradient;
01627   }
01628 
01629   const VecType& getGradient() const {
01630     if ( !_gradient ) {
01631       throw aol::Exception ( "gradient not computed yet.", __FILE__, __LINE__ );
01632     }
01633     return *_gradient;
01634   }
01635 
01636   void appendReg ( const aol::Op<aol::MultiVector<RealType> > &op, RealType factor ) {
01637     if ( factor != 0. ) {
01638       _linCombOpReg.appendReference ( op, factor );
01639     }
01640   }
01641 
01642   void appendUnreg ( const aol::Op<aol::MultiVector<RealType> > &op, RealType factor ) {
01643     if ( factor != 0. ) {
01644       _linCombOpUnreg.appendReference ( op, factor );
01645     }
01646   }
01647 
01648   virtual void applyAdd ( const aol::MultiVector<RealType> &arg, aol::MultiVector<RealType> &dest ) const {
01649     VecType  tmp ( arg );
01650     if ( !_gradient ) {
01651       _gradient = new VecType ( arg );
01652     }
01653     VecType &gradient = *_gradient;
01654 
01655     // collect all gradients and reverse sign
01656     _linCombOpReg.apply ( arg, tmp );
01657 
01658     cerr << sqrt ( tmp.normSqr() ) << endl;
01659     // multiply with inverse lumped mass matrix
01660     _massOpInv.apply ( tmp, gradient );
01661     cerr << sqrt ( gradient.normSqr() ) << "   max x " << gradient[0].getMaxValue() << "   max y " << gradient[1].getMaxValue() << endl;
01662 
01663 
01664     qc::LinearSmoothOp<RealType> _linSmooth;
01665     _linSmooth.setCurrentGrid ( _grid );
01666     _linSmooth.setSigma ( _smoothSigma );
01667     for ( int c = 0; c < dest.numComponents(); c++ ) {
01668       _linSmooth.applyAdd ( gradient[c], dest[c] );
01669     }
01670 
01671     _linCombOpUnreg.apply ( arg, tmp );
01672     _massOpInv.applyAdd ( tmp, gradient );
01673     _massOpInv.applyAdd ( tmp, dest );
01674 
01675     dest *= -1.;
01676 
01677   }
01678 
01679 };
01680 
01681 //! Computes \f$ \frac{1}{2}\int |Du+Du^T|^2 \f$, where marg=u.
01682 template <typename ConfiguratorType>
01683 class SymmetricLengthEnergy : public aol::FENonlinIntegrationVectorInterface < ConfiguratorType,
01684                                                                                SymmetricLengthEnergy<ConfiguratorType> > {
01685 public:
01686   typedef typename ConfiguratorType::RealType RealType;
01687 
01688   SymmetricLengthEnergy ( const typename ConfiguratorType::InitType &Grid )
01689     : aol::FENonlinIntegrationVectorInterface < ConfiguratorType,
01690                                                 SymmetricLengthEnergy<ConfiguratorType> > ( Grid ) {};
01691 
01692   RealType evaluateIntegrand ( const aol::auto_container<ConfiguratorType::Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &DiscrFuncs,
01693                                const typename ConfiguratorType::ElementType &El,
01694                                int QuadPoint, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01695     typename ConfiguratorType::VecType grad;
01696     typename ConfiguratorType::MatType matrix;
01697     
01698     for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
01699       DiscrFuncs[i].evaluateGradientAtQuadPoint ( El, QuadPoint, grad );
01700       for ( int j = 0; j < ConfiguratorType::Dim; j++ ) {
01701         matrix[i][j] += grad[j];
01702         matrix[j][i] += grad[j];
01703       }
01704     }
01705 
01706     return 0.5 * matrix.normSqr();
01707   }
01708 };
01709 
01710 //! Computes the first variation of SymmetricLengthEnergy, e.q. \f$ 2 int (Du+Du^T):D\zeta \f$, where marg=u.
01711 template <typename ConfiguratorType>
01712 class VariationOfSymmetricLengthEnergy
01713   : public aol::FENonlinVectorDiffOpInterface<ConfiguratorType, ConfiguratorType::Dim, ConfiguratorType::Dim, VariationOfSymmetricLengthEnergy<ConfiguratorType> > {
01714 public:
01715   typedef typename ConfiguratorType::RealType RealType;
01716 
01717   VariationOfSymmetricLengthEnergy ( const typename ConfiguratorType::InitType &Grid )
01718     : aol::FENonlinVectorDiffOpInterface<ConfiguratorType, ConfiguratorType::Dim, ConfiguratorType::Dim, VariationOfSymmetricLengthEnergy<ConfiguratorType> > ( Grid ) {};
01719 
01720   void getNonlinearity ( const aol::auto_container<ConfiguratorType::Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &DiscFuncs,
01721                          const typename ConfiguratorType::ElementType &El,
01722                          int QuadPoint, const typename ConfiguratorType::DomVecType &/*RefCoord*/,
01723                          typename aol::Mat<ConfiguratorType::Dim, ConfiguratorType::Dim, RealType> &NL ) const {
01724     typename ConfiguratorType::VecType grad;
01725     NL.setZero();
01726     for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
01727       DiscFuncs[i].evaluateGradientAtQuadPoint ( El, QuadPoint, grad );
01728       for ( int j = 0; j < ConfiguratorType::Dim; j++ ) {
01729         NL[i][j] += grad[j];
01730         NL[j][i] += grad[j];
01731       }
01732     }
01733     NL *= 2;
01734   }
01735 };
01736 
01737 //! Computes \f$ \frac{1}{2}\int (\partial_iu_j-\partial_ju_i \f$, where marg=u.
01738 template <typename ConfiguratorType>
01739 class SkewSymmetricMean : public aol::FENonlinIntegrationVectorInterface < ConfiguratorType,
01740                                                                            SkewSymmetricMean<ConfiguratorType> > {
01741 public:
01742   typedef typename ConfiguratorType::RealType RealType;
01743 protected:
01744   const int _i, _j;
01745 public:
01746   SkewSymmetricMean ( const typename ConfiguratorType::InitType &Grid,
01747                       const int I,
01748                       const int J )
01749     : aol::FENonlinIntegrationVectorInterface < ConfiguratorType,
01750                                                   SkewSymmetricMean<ConfiguratorType> > ( Grid ),
01751       _i(I),
01752       _j(J) {};
01753 
01754   RealType evaluateIntegrand ( const aol::auto_container<ConfiguratorType::Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &DiscrFuncs,
01755                                const typename ConfiguratorType::ElementType &El,
01756                                int QuadPoint, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
01757     typename ConfiguratorType::VecType grad;
01758     DiscrFuncs[_j].evaluateGradientAtQuadPoint ( El, QuadPoint, grad );
01759     RealType temp = grad[_i];
01760     DiscrFuncs[_i].evaluateGradientAtQuadPoint ( El, QuadPoint, grad );
01761     temp -= grad[_j];
01762     return 0.5 * temp;
01763   }
01764 };
01765 
01766 /**
01767  * \brief Computes via the method "apply(...)" \f$\frac1{h^d}(\int_{\Omega}(\|D\phi\|,\|cof(D\phi)\|,det(D\phi))^T\psi_idx)_i\f$,
01768  * where the domain \f$\Omega\f$ is given by a grid (passed to the constructor), \f$\psi_i\f$ denotes the finite element functions,
01769  * and \f$D\phi\f$ denotes the deformation gradient of a deformation \f$\phi=identity+d\f$, where the displacement \f$d\f$
01770  * is passed to "apply(...)". Expressed differently, the three invariants of the deformation gradient are evaluated on the grid.
01771  *
01772  * \author Wirth
01773  */
01774 template <typename ConfiguratorType>
01775 class DeformationGradientInvariants :
01776   public aol::FENonlinVectorOpInterface<ConfiguratorType, ConfiguratorType::Dim, 3, DeformationGradientInvariants<ConfiguratorType> > {
01777 
01778 private:
01779   typedef typename ConfiguratorType::RealType RealType;
01780 
01781   const RealType _gridCellVolume;
01782 
01783 public:
01784   DeformationGradientInvariants( const qc::GridDefinition &Grid ) :
01785     // load the grid
01786     aol::FENonlinVectorOpInterface<ConfiguratorType, ConfiguratorType::Dim, 3, DeformationGradientInvariants<ConfiguratorType> >( Grid ),
01787     _gridCellVolume( std::pow( Grid.H(), ConfiguratorType::Dim ) ) {
01788   }
01789 
01790   /**
01791    * \brief Computes \f$(\|D\phi\|,\|cof(D\phi)\|,det(D\phi))^T\f$.
01792    */
01793   void getNonlinearity( const aol::auto_container< ConfiguratorType::Dim, aol::DiscreteFunctionDefault< ConfiguratorType > > &DiscFuncs,
01794                         const typename ConfiguratorType::ElementType &El,
01795                         int QuadPoint,
01796                         const typename ConfiguratorType::DomVecType /*&LocalCoord*/,
01797                         aol::Vec< 3, RealType > &NL ) const {
01798     // compute the transposed displacement gradient and put it into "dphi"
01799     typename ConfiguratorType::MatType dphi;
01800     for ( int j=0; j<ConfiguratorType::Dim; j++ )
01801       DiscFuncs[j].evaluateGradientAtQuadPoint ( El, QuadPoint, dphi[j] );
01802     // compute the transposed deformation gradient "dphi"
01803     for ( int j=0; j<ConfiguratorType::Dim; j++ )
01804       dphi[j][j] += 1.;
01805     // return norm, cofactor-norm, and determinant, weighted by the inverse volume/area of a grid cell
01806     typename ConfiguratorType::MatType cof;
01807     cof.makeCofactorMatrix( dphi );
01808     NL[0] = dphi.norm();
01809     NL[1] = cof.norm();
01810     NL[2] = dphi.det();
01811     NL /= _gridCellVolume;
01812   }
01813 };
01814 
01815 } // end namespace qc
01816 
01817 #endif

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