QuOc

 

deformations.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 __DEFORMATIONS_H
00012 #define __DEFORMATIONS_H
00013 
00014 #include <generator.h>
00015 #include <Newton.h>
00016 #include <multiArray.h>
00017 
00018 namespace qc {
00019 
00020 template <typename ConfiguratorType>
00021 void deformImageWithCoarseDeformation ( const aol::Vector<typename ConfiguratorType::RealType> &Image,
00022                                         const typename ConfiguratorType::InitType &Finegrid,
00023                                         const typename ConfiguratorType::InitType &Coarsegrid,
00024                                         qc::Array<typename ConfiguratorType::RealType>  &DeformedImage,
00025                                         const aol::MultiVector<typename ConfiguratorType::RealType> &Phidofs ) {
00026   typedef typename ConfiguratorType::ArrayType ArrayType;
00027   typedef typename ConfiguratorType::RealType RealType;
00028 
00029   aol::auto_container< ConfiguratorType::Dim, const ArrayType > phi;
00030   ArrayType     aImage ( Image, Finegrid );
00031 
00032   for ( int c = 0; c < Phidofs.numComponents(); c++ ) {
00033     ArrayType temp ( Phidofs[c], Coarsegrid );
00034     phi.set_copy ( c, temp );
00035   }
00036 
00037   //ArrayType templ( Image[Finegrid.getGridDepth()], Finegrid );
00038 
00039   const RealType h = Finegrid.H();
00040   const int num = Finegrid.getWidth();
00041 
00042   typename ConfiguratorType::InitType::OldAllNodeIterator fnit;
00043   for ( fnit = Finegrid._nBeginIt; fnit != Finegrid._nEndIt; ++fnit ) {
00044     RealType scaling =  ( static_cast<RealType> ( Coarsegrid.getWidth() - 1 ) ) / ( static_cast<RealType> ( Finegrid.getWidth() - 1 ) );
00045     typename ConfiguratorType::VecType ds;
00046     typename ConfiguratorType::VecType position;
00047     for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00048       position[i] = ( *fnit ) [i];
00049     }
00050     position *= scaling;
00051     for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00052       ds[i] = ( *fnit ) [i] + phi[i].interpolate ( position ) / h;
00053       ds[i] = aol::Max ( aol::Min ( ds[i], static_cast<RealType> ( num - 1 ) ), aol::ZOTrait<RealType>::zero );
00054     }
00055     RealType value =  aImage.interpolate ( ds );
00056     DeformedImage.set ( *fnit, value );
00057   }
00058   return;
00059 
00060 
00061 }
00062 
00063 /**
00064  * \brief DeformedImage = Image circ Phi
00065  *
00066  * \author Berkels
00067  */
00068 template <typename ConfiguratorType>
00069 void DeformImage ( const aol::Vector<typename ConfiguratorType::RealType> &Image,
00070                    const typename ConfiguratorType::InitType &Grid,
00071                    aol::Vector<typename ConfiguratorType::RealType> &DeformedImage,
00072                    const aol::MultiVector<typename ConfiguratorType::RealType> &Phi,
00073                    const bool ExtendWithConstant = true,
00074                    const typename ConfiguratorType::RealType ExtensionConstant = 0 ) {
00075   typedef typename ConfiguratorType::RealType RealType;
00076   const typename ConfiguratorType::ArrayType imageArray ( Image, Grid, aol::FLAT_COPY );
00077   typename ConfiguratorType::ArrayType deformedImageArray ( DeformedImage, Grid, aol::FLAT_COPY );
00078   const qc::MultiArray<RealType, ConfiguratorType::Dim> phiMArray ( Grid, Phi, aol::FLAT_COPY ); 
00079 
00080   const aol::Vec3<int> gridSize = Grid.getSize();
00081   const RealType h = static_cast<RealType> ( Grid.H() );
00082 
00083   typename ConfiguratorType::InitType::OldAllNodeIterator fnit;
00084   for ( fnit = Grid._nBeginIt; fnit != Grid._nEndIt; ++fnit ) {
00085     typename ConfiguratorType::VecType ds;
00086 
00087     bool transformPositionInDomain = true;
00088 
00089     for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00090       ds[i] = ( *fnit ) [i] + phiMArray[i].get ( *fnit ) / h;
00091       if ( ExtendWithConstant == true ) {
00092         if ( ds[i] < aol::ZOTrait<RealType>::zero || ds[i] > static_cast<RealType> ( gridSize[i] - 1 ) )
00093           transformPositionInDomain = false;
00094       } else
00095         ds[i] = aol::Clamp ( ds[i], aol::ZOTrait<RealType>::zero, static_cast<RealType> ( gridSize[i] - 1 ) );
00096     }
00097 
00098     RealType value = ExtensionConstant;
00099     if ( transformPositionInDomain == true )
00100       value = imageArray.interpolate ( ds );
00101     deformedImageArray.set ( *fnit, value );
00102   }
00103 }
00104 
00105 /**
00106  * @brief DeformedImage = Image circ Phi,
00107  * where deformed image is extended by a second image in regions, where the deformed image is not defined.
00108  *
00109  * @author Wirth
00110  */
00111 template <typename ConfiguratorType>
00112 void DeformImage ( const aol::Vector<typename ConfiguratorType::RealType> &Image,
00113                    const typename ConfiguratorType::InitType &Grid,
00114                    aol::Vector<typename ConfiguratorType::RealType> &DeformedImage,
00115                    const aol::MultiVector<typename ConfiguratorType::RealType> &Phi,
00116                    const aol::Vector<typename ConfiguratorType::RealType> &ExtendImage ) {
00117   typedef typename ConfiguratorType::RealType RealType;
00118   typename ConfiguratorType::ArrayType imageArray ( Image, Grid, aol::FLAT_COPY );
00119   typename ConfiguratorType::ArrayType deformedImageArray ( DeformedImage, Grid, aol::FLAT_COPY );
00120   typename ConfiguratorType::ArrayType extendImageArray ( ExtendImage, Grid, aol::FLAT_COPY );
00121 
00122   const int num = Grid.getNumX();
00123   RealType h = Grid.H();
00124 
00125   typename ConfiguratorType::InitType::OldFullNodeIterator fnit;
00126   qc::FastILexMapper<ConfiguratorType::Dim> mapper ( Grid );
00127   for ( fnit = Grid.begin(); fnit != Grid.end(); ++fnit ) {
00128     typename ConfiguratorType::VecType ds;
00129 
00130     bool transformPositionInDomain = true;
00131 
00132     for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00133       ds[i] = ( *fnit ) [i] + Phi[i].get ( mapper.getGlobalIndex ( *fnit ) ) / h;
00134       if ( ds[i] < aol::ZOTrait<RealType>::zero || ds[i] > static_cast<RealType> ( num - 1 ) )
00135         transformPositionInDomain = false;
00136     }
00137     if ( transformPositionInDomain == true )
00138       deformedImageArray.set ( *fnit, imageArray.interpolate ( ds ) );
00139     else
00140       deformedImageArray.set ( *fnit, extendImageArray.get ( *fnit ) );
00141   }
00142   return;
00143 }
00144 
00145 template <typename ConfiguratorType>
00146 void DeformImageFromMVecPart ( const aol::Vector<typename ConfiguratorType::RealType> &Image,
00147                                const typename ConfiguratorType::InitType &Grid,
00148                                aol::Vector<typename ConfiguratorType::RealType> &DeformedImage,
00149                                const aol::MultiVector<typename ConfiguratorType::RealType> &MultiVec,
00150                                const bool ExtendWithZero,
00151                                const int StartIndex ) {
00152   aol::MultiVector<typename ConfiguratorType::RealType> phi ( 0, 0 );
00153   for ( int i = StartIndex; i < StartIndex + ConfiguratorType::Dim; i++ )
00154     phi.appendReference ( MultiVec[i] );
00155   DeformImage<ConfiguratorType> ( Image, Grid, DeformedImage, phi, ExtendWithZero );
00156 }
00157 
00158 /**
00159  * @brief DeformedImage = Image\f$\circ\phi^{-1}\f$, where \f$\phi\f$=D+identity.
00160  * "VectorType" may be either aol::Vector or aol::MultiVector.
00161  *
00162  * @author Wirth
00163  */
00164 template <typename ConfiguratorType, typename VectorType>
00165 void InvDeformImage ( const VectorType &Image,
00166                       const typename ConfiguratorType::InitType &Grid,
00167                       VectorType &DeformedImage,
00168                       const aol::MultiVector<typename ConfiguratorType::RealType> &D ) {
00169 
00170   // deform the image
00171   qc::TransformFunction<typename ConfiguratorType::RealType,ConfiguratorType::Dim> transformFunction( Grid );
00172   transformFunction.setDeformation( D );
00173   transformFunction.apply( Image, DeformedImage );
00174 }
00175 
00176 /**
00177  * @brief DeformedImage = Image\f$\circ\phi^{-1}\f$, where \f$\phi\f$=D+identity.
00178  * The deformed image is extended by "ExtendImage", where it is not defined.
00179  * "VectorType" may be either aol::Vector or aol::MultiVector.
00180  *
00181  * @author Wirth
00182  */
00183 template <typename ConfiguratorType, typename VectorType>
00184 void InvDeformImage ( const VectorType &Image,
00185                       const typename ConfiguratorType::InitType &Grid,
00186                       VectorType &DeformedImage,
00187                       const aol::MultiVector<typename ConfiguratorType::RealType> &D,
00188                       const VectorType &ExtendImage ) {
00189 
00190   // deform the image
00191   qc::TransformFunction<typename ConfiguratorType::RealType,ConfiguratorType::Dim> transformFunction( Grid );
00192   transformFunction.setDeformation( D );
00193   typename qc::BitArray<ConfiguratorType::Dim> mask;
00194   transformFunction.transform( Image, DeformedImage, mask, ExtendImage );
00195 }
00196 
00197 /**
00198  * @brief DeformedImage = Image\f$\circ\phi^{-1}\f$, where \f$\phi\f$=D+identity.
00199  * The deformed image is extended by "ExtendImage", where it is not defined.
00200  * "VectorType" may be either aol::Vector or aol::MultiVector.
00201  *
00202  * @author Wirth
00203  */
00204 template <typename ConfiguratorType, typename VectorType>
00205 void InvDeformImage ( const VectorType &Image,
00206                       const typename ConfiguratorType::InitType &Grid,
00207                       VectorType &DeformedImage,
00208                       const aol::MultiVector<typename ConfiguratorType::RealType> &D,
00209                       const VectorType &ExtendImage,
00210                       typename qc::BitArray<ConfiguratorType::Dim> &InvDispDefined ) {
00211 
00212   // deform the image
00213   qc::TransformFunction<typename ConfiguratorType::RealType,ConfiguratorType::Dim> transformFunction( Grid );
00214   transformFunction.setDeformation( D );
00215   transformFunction.transform( Image, DeformedImage, InvDispDefined, ExtendImage );
00216 }
00217 
00218 /**
00219  * @brief Copies Image to ExtendedImage, but replaces values at all nodes where ValueDefined is false such that L2-norm of gradient is minimized.
00220  * "VectorType" may be either aol::Vector or aol::MultiVector.
00221  *
00222  * @author Wirth
00223  */
00224 template <typename ConfiguratorType, typename VectorType>
00225 void SmoothlyExtendImage ( const VectorType &Image,
00226                            const typename ConfiguratorType::InitType &Grid,
00227                            VectorType &ExtendedImage,
00228                            // VC++ 2008 somehow doesn't see here that ConfiguratorType::Dim is of type qc::Dimension,
00229                            // thus we have to use qc::BitArrayTrait which can use a simple int.
00230                            const typename qc::BitArrayTrait<ConfiguratorType::Dim>::ArrayType &ValueDefined ) {
00231   typedef typename ConfiguratorType::RealType RealType;
00232 
00233   aol::MultiVector<RealType> image, extendedImage;
00234   image.appendReference( Image );
00235   extendedImage.appendReference( ExtendedImage );
00236 
00237   // smoothly extend Image $u$ to ExtendedImage $v$ by minimizing $\int_\Omega|\nabla v|^2dx$
00238   // (i.e.solving $\Delta v=0$) under boundary conditions $v=u$ where ValueDefined = true
00239 
00240   // compute boundary conditions
00241   aol::MultiVector<RealType> bc( image, aol::DEEP_COPY );
00242   for ( int i = 0; i < bc.numComponents(); i++ )
00243     for ( int j = 0; j < ValueDefined.size(); j++ )
00244       if ( !ValueDefined[j] )
00245         bc[i][j] = 0.;
00246   // compute right hand side
00247   aol::StiffOp<ConfiguratorType> stiffOp( Grid );
00248   aol::MultiVector<RealType> rhs( image, aol::STRUCT_COPY );
00249   for ( int i = 0; i < rhs.numComponents(); i++ ){
00250     stiffOp.apply( bc[i], rhs[i] );
00251     rhs[i] *= -1.;
00252     // build in boundary conditions
00253     for ( int j = 0; j < ValueDefined.size(); j++ )
00254       if ( ValueDefined[j] )
00255         rhs[i][j] = bc[i][j];
00256   }
00257   // find the solution
00258   typename ConfiguratorType::MatrixType systemMatrix( Grid );
00259   stiffOp.assembleAddMatrix( systemMatrix, &ValueDefined, true );
00260   aol::DiagonalPreconditioner<aol::Vector<RealType> > preCond( systemMatrix );
00261   aol::PCGInverse<aol::Vector<RealType> > inv( systemMatrix, preCond, 1.e-16, 1000, aol::STOPPING_ABSOLUTE );
00262   inv.setQuietMode( true );
00263   for ( int j = 0; j < rhs.numComponents(); j++ )
00264     inv.apply( rhs[j], extendedImage[j] );
00265 }
00266 
00267 /**
00268  * @brief DeformedImage = \f$L(\f$Image\f$\circ\phi^{-1})\f$, where \f$\phi\f$=D+identity.
00269  * The deformed image is smoothly extended via the operator \f$L\f$, where the inverse deformation is not defined.
00270  * \f$L\f$ minimizes the L2-norm of the gradient.
00271  * "VectorType" may be either aol::Vector or aol::MultiVector.
00272  *
00273  * @author Wirth
00274  */
00275 template <typename ConfiguratorType, typename VectorType>
00276 void InvDeformAndSmoothlyExtendImage ( const VectorType &Image,
00277                                        const typename ConfiguratorType::InitType &Grid,
00278                                        VectorType &DeformedImage,
00279                                        const aol::MultiVector<typename ConfiguratorType::RealType> &D ) {
00280 
00281   // deform the image
00282   typename qc::BitArray<ConfiguratorType::Dim> invDispDefined;
00283   VectorType deformationResult( Image, aol::STRUCT_COPY );
00284   qc::TransformFunction<typename ConfiguratorType::RealType,ConfiguratorType::Dim> transformFunction( Grid );
00285   transformFunction.setDeformation( D );
00286   transformFunction.transform( Image, deformationResult, invDispDefined );
00287   // smoothly extend the image, where the inverse deformation was not defined
00288   qc::SmoothlyExtendImage<ConfiguratorType,VectorType>( deformationResult, Grid, DeformedImage, invDispDefined );
00289 }
00290 
00291 /**
00292  * @brief identity + ResultDisplacement = \f$L(\phi_1\circ\phi_2^{-1})\f$, where \f$\phi_i\f$=Displacementi+identity.
00293  * The resulting displacement is smoothly extended via the operator \f$L\f$ where the displacement is not defined.
00294  * \f$L\f$ minimizes the L2-norm of the gradient.
00295  *
00296  * @author Wirth
00297  */
00298 template <typename ConfiguratorType>
00299 void InvConcatAndSmoothlyExtendDeformations ( const aol::MultiVector<typename ConfiguratorType::RealType> &Displacement1,
00300                                               const aol::MultiVector<typename ConfiguratorType::RealType> &Displacement2,
00301                                               const typename ConfiguratorType::InitType &Grid,
00302                                               aol::MultiVector<typename ConfiguratorType::RealType> &ResultDisplacement ) {
00303   // compute the deformation \phi_1
00304   aol::MultiVector<typename ConfiguratorType::RealType> identity( Grid ), deformation1( Displacement1 ), resultDisplacement( Grid );
00305   qc::DataGenerator<ConfiguratorType>( Grid ).generateIdentity( identity );
00306   deformation1 += identity;
00307 
00308   // concatenate \phi_1 with \phi_2^{-1}
00309   typename qc::BitArray<ConfiguratorType::Dim> invDispDefined;
00310   qc::TransformFunction<typename ConfiguratorType::RealType,ConfiguratorType::Dim> transformFunction( Grid );
00311   transformFunction.setDeformation( Displacement2 );
00312   transformFunction.transform( deformation1, resultDisplacement, invDispDefined );
00313 
00314   // compute the corresponding displacement
00315   resultDisplacement -= identity;
00316 
00317   // smoothly extend the displacement where it is not yet defined
00318   qc::SmoothlyExtendImage<ConfiguratorType,aol::MultiVector<typename ConfiguratorType::RealType> >( resultDisplacement, Grid, ResultDisplacement, invDispDefined );
00319 }
00320 
00321 /**
00322  * @brief DeformedImage = \f$L(\f$Image\f$\circ\phi)\f$, where \f$\phi\f$=D+identity.
00323  * The deformed image is smoothly extended via the operator \f$L\f$, where the image is not defined.
00324  * \f$L\f$ minimizes the L2-norm of the gradient.
00325  * "VectorType" may be either aol::Vector or aol::MultiVector.
00326  *
00327  * @author Wirth
00328  */
00329 template <typename ConfiguratorType, typename VectorType>
00330 void DeformAndSmoothlyExtendImage ( const VectorType &Image,
00331                                     const typename ConfiguratorType::InitType &Grid,
00332                                     VectorType &DeformedImage,
00333                                     const aol::MultiVector<typename ConfiguratorType::RealType> &D ) {
00334   aol::MultiVector<typename ConfiguratorType::RealType> image, deformedImage;
00335   image.appendReference( Image );
00336   deformedImage.appendReference( DeformedImage );
00337 
00338   // deform the image
00339   aol::MultiVector<typename ConfiguratorType::RealType> deformationResult( image, aol::STRUCT_COPY );
00340   aol::Vector<typename ConfiguratorType::RealType> extendImage( Grid.getNumberOfNodes() );
00341   extendImage.setAll( aol::NumberTrait<typename ConfiguratorType::RealType>::NaN );
00342   for ( int i = 0; i < image.numComponents(); i++ )
00343     DeformImage<ConfiguratorType>( image[i], Grid, deformationResult[i], D, extendImage );
00344   typename qc::BitArray<ConfiguratorType::Dim> imageDefined( Grid );
00345   for ( int i = 0; i < imageDefined.size(); i++ )
00346     imageDefined.set( i, !aol::isNaN( deformationResult[0][i] ) );
00347 
00348   // smoothly extend the image where it is not yet defined
00349   qc::SmoothlyExtendImage<ConfiguratorType,aol::MultiVector<typename ConfiguratorType::RealType> >( deformationResult, Grid, deformedImage, imageDefined );
00350 }
00351 
00352 /**
00353  * @brief identity + ResultDisplacement = \f$L(\phi_1\circ\phi_2)\f$, where \f$\phi_i\f$=Displacementi+identity.
00354  * The resulting displacement is smoothly extended via the operator \f$L\f$ where the displacement is not defined.
00355  * \f$L\f$ minimizes the L2-norm of the gradient.
00356  *
00357  * @author Wirth
00358  */
00359 template <typename ConfiguratorType>
00360 void ConcatAndSmoothlyExtendDeformations ( const aol::MultiVector<typename ConfiguratorType::RealType> &Displacement1,
00361                                            const aol::MultiVector<typename ConfiguratorType::RealType> &Displacement2,
00362                                            const typename ConfiguratorType::InitType &Grid,
00363                                            aol::MultiVector<typename ConfiguratorType::RealType> &ResultDisplacement ) {
00364   aol::MultiVector<typename ConfiguratorType::RealType> resultDisplacement( Displacement1, aol::STRUCT_COPY );
00365 
00366   // concatenate the displacement with the deformation
00367   aol::Vector<typename ConfiguratorType::RealType> extendImage( Grid.getNumberOfNodes() );
00368   extendImage.setAll( aol::NumberTrait<typename ConfiguratorType::RealType>::NaN );
00369   for ( int i = 0; i < Displacement1.numComponents(); i++ )
00370     DeformImage<ConfiguratorType>( Displacement1[i], Grid, resultDisplacement[i], Displacement2, extendImage );
00371 
00372   // compute the correct deformation and then displacement
00373   resultDisplacement += Displacement2;
00374 
00375   // smoothly extend the displacement where it is not yet defined
00376   typename qc::BitArray<ConfiguratorType::Dim> imageDefined( qc::GridSize<ConfiguratorType::Dim>::createFrom( Grid ) );
00377   for ( int i = 0; i < imageDefined.size(); i++ )
00378     imageDefined.set( i, !aol::isNaN( resultDisplacement[0][i] ) );
00379   
00380   cerr<<Grid.getGridDepth()<<endl;
00381   cerr<<resultDisplacement[0].size();
00382   cerr<<ResultDisplacement[0].size();
00383   
00384   qc::SmoothlyExtendImage<ConfiguratorType,aol::MultiVector<typename ConfiguratorType::RealType> >( resultDisplacement, Grid, ResultDisplacement, imageDefined );
00385 }
00386 
00387 /**
00388  * Loads image (possibly colored) named "InputImageFileName"
00389  * deforms it with "Identity + (InputFileNameDefX, InputFileNameDefY)"
00390  * and writes the result to OutputFileName.
00391  *
00392  * \author Berkels
00393  */
00394 template <typename ConfiguratorType>
00395 void deformAndSaveColoredImage ( const char *InputImageFileName, const char *InputFileNameDefX, const char *InputFileNameDefY, const char *OutputFileName, const bool UseInverseDeformation = true ){
00396   typedef typename ConfiguratorType::RealType RealType;
00397   qc::ScalarArray<RealType,qc::QC_2D> defX( InputFileNameDefX );
00398   qc::ScalarArray<RealType,qc::QC_2D> defY( InputFileNameDefY );
00399   aol::MultiVector<RealType> deformation ( 0, 0 );
00400   deformation.appendReference ( defX );
00401   deformation.appendReference ( defY );
00402 
00403   qc::MultiArray< RealType, 2, 3 > deformedImage( defX.getNumX(), defX.getNumY() );
00404   deformedImage.loadPNG( InputImageFileName );
00405 
00406   const int depth = qc::logBaseTwo ( defX.getNumX() );
00407   qc::GridDefinition grid( depth, qc::QC_2D );
00408 
00409   qc::ScalarArray<RealType,qc::QC_2D> tmp( deformedImage[0], aol::STRUCT_COPY );
00410   for ( int i = 0; i < 3; i++ ){
00411     if ( UseInverseDeformation )
00412       qc::DeformImage<ConfiguratorType>( deformedImage[i], grid, tmp, deformation);
00413     else
00414       qc::InvDeformImage<ConfiguratorType>( deformedImage[i], grid, tmp, deformation);
00415     deformedImage[i] = tmp;
00416   }
00417 
00418   deformedImage.setQuietMode( true );
00419   deformedImage.savePNG( OutputFileName );
00420 }
00421 
00422 template <typename ConfiguratorType>
00423 void concatenateTwoDeformations ( const qc::GridDefinition &Grid,
00424                                   const aol::MultiVector<typename ConfiguratorType::RealType> &Deformation1,
00425                                   const aol::MultiVector<typename ConfiguratorType::RealType> &Deformation2,
00426                                   aol::MultiVector<typename ConfiguratorType::RealType> &Dest ) {
00427   typedef typename ConfiguratorType::RealType RealType;
00428 
00429   qc::MultiArray<RealType, ConfiguratorType::Dim> deformation2Array ( Grid, Deformation2 );
00430   qc::MultiArray<RealType, ConfiguratorType::Dim> destArray ( Grid, Dest );
00431 
00432   const int num = Grid.getWidth();
00433   RealType h = Grid.H();
00434 
00435   qc::GridDefinition::OldFullNodeIterator fnit;
00436   qc::FastILexMapper<ConfiguratorType::Dim> mapper ( Grid );
00437   for ( fnit = Grid.begin(); fnit != Grid.end(); ++fnit ) {
00438     typename ConfiguratorType::VecType ds;
00439 
00440     bool transformPositionInDomain = true;
00441 
00442     for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00443       ds[i] = ( *fnit ) [i] + Deformation1[i].get ( mapper.getGlobalIndex ( *fnit ) ) / h;
00444       if ( ds[i] < aol::ZOTrait<RealType>::zero || ds[i] > static_cast<RealType> ( num - 1 ) )
00445         transformPositionInDomain = false;
00446     }
00447 
00448     for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00449       RealType value = 0.;
00450       if ( transformPositionInDomain == true ) {
00451         value = ( deformation2Array[i].interpolate ( ds ) ) + Deformation1[i].get ( mapper.getGlobalIndex ( *fnit ) );
00452       }
00453       destArray[i].set ( *fnit, value );
00454     }
00455   }
00456   return;
00457 }
00458 
00459 /**
00460  * Shifts coordinates by Offset and checks whether the transformed coordinates are still
00461  * in [0,1]^d. Returns true if yes and false if not. Offset contains the offset relative
00462  * to [0,1]^d, i.e. the pixel offset is Offset / H
00463  *
00464  * \author Berkels
00465  */
00466 template <typename ConfiguratorType>
00467 inline bool transformCoord ( const typename ConfiguratorType::InitType &Grid,
00468                              const typename ConfiguratorType::VecType &Coord,
00469                              const typename ConfiguratorType::VecType &Offset,
00470                              qc::Element &TransformedEl,
00471                              typename ConfiguratorType::VecType &TransformedLocalCoord ) {
00472   typename ConfiguratorType::VecType transformedCoord;
00473   for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00474     transformedCoord[i] = Coord[i] + Offset[i] / Grid.H();
00475     TransformedEl[i] = static_cast<short> ( transformedCoord[i] );
00476     TransformedLocalCoord[i] = transformedCoord[i] - TransformedEl[i];
00477 
00478     if ( ( TransformedEl[i] == ( (Grid.getSize())[i] - 1 ) ) && ( TransformedLocalCoord[i] == 0 ) ) {
00479       TransformedEl[i] = (Grid.getSize())[i] - 2;
00480       TransformedLocalCoord[i] = 1;
00481     }
00482 
00483     if ( ( transformedCoord[i] < 0. ) || ( TransformedEl[i] < 0 ) || ( TransformedEl[i] >= ( (Grid.getSize())[i] - 1 ) ) ) {
00484       return false;
00485     }
00486   }
00487   return true;
00488 }
00489 
00490 template <typename ConfiguratorType>
00491 inline bool transformCoord ( const typename ConfiguratorType::InitType &Grid,
00492                              const typename ConfiguratorType::ElementType &El,
00493                              const typename ConfiguratorType::VecType &RefCoord,
00494                              const typename ConfiguratorType::VecType &Offset,
00495                              qc::Element &TransformedEl,
00496                              typename ConfiguratorType::VecType &TransformedLocalCoord ) {
00497   typename ConfiguratorType::VecType coord;
00498   for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00499     coord[i] = El[i] + RefCoord[i];
00500   }
00501   return qc::transformCoord<ConfiguratorType> ( Grid, coord, Offset, TransformedEl, TransformedLocalCoord );
00502 }
00503 
00504 template <typename ConfiguratorType>
00505 inline bool transformCoord ( const typename ConfiguratorType::InitType &Grid,
00506                              const aol::auto_container<ConfiguratorType::Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &DiscrTransformation,
00507                              const typename ConfiguratorType::ElementType &El,
00508                              const int QuadPoint,
00509                              const typename ConfiguratorType::VecType &RefCoord,
00510                              qc::Element &TransformedEl,
00511                              typename ConfiguratorType::VecType &TransformedLocalCoord ) {
00512   typename ConfiguratorType::VecType offset;
00513   for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00514     offset[i] = DiscrTransformation[i].evaluateAtQuadPoint ( El, QuadPoint );
00515   }
00516   return qc::transformCoord<ConfiguratorType> ( Grid, El, RefCoord, offset, TransformedEl, TransformedLocalCoord );
00517 }
00518 
00519 template <typename ConfiguratorType>
00520 inline bool transformCoord ( const typename ConfiguratorType::InitType &Grid,
00521                              const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> &DiscrTransformation,
00522                              const typename ConfiguratorType::ElementType &El,
00523                              const int QuadPoint,
00524                              const typename ConfiguratorType::VecType &RefCoord,
00525                              qc::Element &TransformedEl,
00526                              typename ConfiguratorType::VecType &TransformedLocalCoord ) {
00527   typename ConfiguratorType::VecType offset;
00528   DiscrTransformation.evaluateAtQuadPoint ( El, QuadPoint, offset );
00529   return qc::transformCoord<ConfiguratorType> ( Grid, El, RefCoord, offset, TransformedEl, TransformedLocalCoord );
00530 }
00531 
00532 //! Shifts a point (given by element El and local coordinates RefCoord) by Offset and computes corresponding point (TransformedEl, TransformedLocalCoord).
00533 /** If the point is shifted outside [0,1[^d, its coordinates are clipped to the interval [0,1[, and the method returns false.
00534  */
00535 template <typename ConfiguratorType>
00536 inline bool transformAndClipCoord ( const ConfiguratorType &Config,
00537                                     const typename ConfiguratorType::ElementType &El,
00538                                     const typename ConfiguratorType::DomVecType &RefCoord,
00539                                     const typename ConfiguratorType::VecType &Offset,
00540                                     typename ConfiguratorType::ElementType &TransformedEl,
00541                                     typename ConfiguratorType::DomVecType &TransformedLocalCoord ) {
00542   bool insideFlag = true;
00543   typename ConfiguratorType::VecType coord;
00544   Config.getGlobalCoords( El, RefCoord, coord );
00545   coord += Offset;
00546   for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00547     if ( coord[i] < 0. ) {
00548       coord[i] = 0.;
00549       insideFlag = false;
00550     } else if ( coord[i] >= 1. ) {
00551       coord[i] = 1. - std::numeric_limits<typename ConfiguratorType::RealType>::epsilon(); // largest number < 1.0 (for clipping to [0,1[)
00552       insideFlag = false;
00553     }
00554   }
00555   Config.getLocalCoords( coord, TransformedEl, TransformedLocalCoord );
00556   return insideFlag;
00557 }
00558 
00559 //! Shifts a point (given by element El and local coordinates RefCoord) by Offset and computes corresponding point (TransformedEl, TransformedLocalCoord).
00560 /** If the point is shifted outside [0,1[^d, its coordinates are clipped to the interval [0,1[, and CoordinateWithinLimits is false for the clipped coordinates.
00561  */
00562 template <typename ConfiguratorType>
00563 inline void transformAndClipCoord ( const ConfiguratorType &Config,
00564                                     const typename ConfiguratorType::ElementType &El,
00565                                     const typename ConfiguratorType::DomVecType &RefCoord,
00566                                     const typename ConfiguratorType::VecType &Offset,
00567                                     typename ConfiguratorType::ElementType &TransformedEl,
00568                                     typename ConfiguratorType::DomVecType &TransformedLocalCoord,
00569                                     aol::Vec<ConfiguratorType::Dim,bool> &CoordinateWithinLimits ) {
00570   CoordinateWithinLimits.setAll( true );
00571   typename ConfiguratorType::VecType coord;
00572   Config.getGlobalCoords( El, RefCoord, coord );
00573   coord += Offset;
00574   for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00575     if ( coord[i] < 0. ) {
00576       coord[i] = 0.;
00577       CoordinateWithinLimits[i] = false;
00578     } else if ( coord[i] >= 1. ) {
00579       coord[i] = 1. - std::numeric_limits<typename ConfiguratorType::RealType>::epsilon(); // largest number < 1.0 (for clipping to [0,1[)
00580       CoordinateWithinLimits[i] = false;
00581     }
00582   }
00583   Config.getLocalCoords( coord, TransformedEl, TransformedLocalCoord );
00584 }
00585 
00586 //! Shifts a point (given by element El and local coordinates RefCoord) by the transformation DiscrTransformation and computes corresponding point (TransformedEl, TransformedLocalCoord).
00587 /** If the point is shifted outside [0,1]^d, its coordinates are clipped to the interval [0,1], and the method returns false.
00588  */
00589 template <typename ConfiguratorType>
00590 inline bool transformAndClipCoord ( const ConfiguratorType &Config,
00591                                     const aol::auto_container<ConfiguratorType::Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &DiscrTransformation,
00592                                     const typename ConfiguratorType::ElementType &El,
00593                                     const int QuadPoint,
00594                                     const typename ConfiguratorType::DomVecType &RefCoord,
00595                                     typename ConfiguratorType::ElementType &TransformedEl,
00596                                     typename ConfiguratorType::DomVecType &TransformedLocalCoord ) {
00597   typename ConfiguratorType::VecType offset;
00598   for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00599     offset[i] = DiscrTransformation[i].evaluateAtQuadPoint ( El, QuadPoint );
00600   return qc::transformAndClipCoord<ConfiguratorType> ( Config, El, RefCoord, offset, TransformedEl, TransformedLocalCoord );
00601 }
00602 
00603 //! Shifts a point (given by element El and local coordinates RefCoord) by the transformation DiscrTransformation and computes corresponding point (TransformedEl, TransformedLocalCoord).
00604 /** If the point is shifted outside [0,1]^d, its coordinates are clipped to the interval [0,1], and CoordinateWithinLimits is false for the clipped coordinates.
00605  */
00606 template <typename ConfiguratorType>
00607 inline void transformAndClipCoord ( const ConfiguratorType &Config,
00608                                     const aol::auto_container<ConfiguratorType::Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &DiscrTransformation,
00609                                     const typename ConfiguratorType::ElementType &El,
00610                                     const int QuadPoint,
00611                                     const typename ConfiguratorType::DomVecType &RefCoord,
00612                                     typename ConfiguratorType::ElementType &TransformedEl,
00613                                     typename ConfiguratorType::DomVecType &TransformedLocalCoord,
00614                                     aol::Vec<ConfiguratorType::Dim,bool> &CoordinateWithinLimits ) {
00615   typename ConfiguratorType::VecType offset;
00616   for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00617     offset[i] = DiscrTransformation[i].evaluateAtQuadPoint ( El, QuadPoint );
00618   qc::transformAndClipCoord<ConfiguratorType> ( Config, El, RefCoord, offset, TransformedEl, TransformedLocalCoord, CoordinateWithinLimits );
00619 }
00620 
00621 //! Shifts a point (given by element El and local coordinates RefCoord) by the transformation DiscrTransformation and computes corresponding point (TransformedEl, TransformedLocalCoord).
00622 /** If the point is shifted outside [0,1]^d, its coordinates are clipped to the interval [0,1], and the method returns false.
00623  */
00624 template <typename ConfiguratorType>
00625 inline bool transformAndClipCoord ( const ConfiguratorType &Config,
00626                                     const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> &DiscrTransformation,
00627                                     const typename ConfiguratorType::ElementType &El,
00628                                     const int QuadPoint,
00629                                     const typename ConfiguratorType::DomVecType &RefCoord,
00630                                     typename ConfiguratorType::ElementType &TransformedEl,
00631                                     typename ConfiguratorType::DomVecType &TransformedLocalCoord ) {
00632   typename ConfiguratorType::VecType offset;
00633   DiscrTransformation.evaluateAtQuadPoint ( El, QuadPoint, offset );
00634   return qc::transformAndClipCoord<ConfiguratorType> ( Config, El, RefCoord, offset, TransformedEl, TransformedLocalCoord );
00635 }
00636 
00637 //! Shifts a point (given by element El and local coordinates RefCoord) by the transformation DiscrTransformation and computes corresponding point (TransformedEl, TransformedLocalCoord).
00638 /** If the point is shifted outside [0,1]^d, its coordinates are clipped to the interval [0,1], and CoordinateWithinLimits is false for the clipped coordinates.
00639  */
00640 template <typename ConfiguratorType>
00641 inline void transformAndClipCoord ( const ConfiguratorType &Config,
00642                                     const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> &DiscrTransformation,
00643                                     const typename ConfiguratorType::ElementType &El,
00644                                     const int QuadPoint,
00645                                     const typename ConfiguratorType::DomVecType &RefCoord,
00646                                     typename ConfiguratorType::ElementType &TransformedEl,
00647                                     typename ConfiguratorType::DomVecType &TransformedLocalCoord,
00648                                     aol::Vec<ConfiguratorType::Dim,bool> &CoordinateWithinLimits ) {
00649   typename ConfiguratorType::VecType offset;
00650   DiscrTransformation.evaluateAtQuadPoint ( El, QuadPoint, offset );
00651   qc::transformAndClipCoord<ConfiguratorType> ( Config, El, RefCoord, offset, TransformedEl, TransformedLocalCoord, CoordinateWithinLimits );
00652 }
00653 
00654 //! Shifts coordinates by TranslationOffset and then deforms the shifted coordinates with DiscrTransformation.
00655 //! Returns true if the transformed coordinates are still in [0,1]^d after each step and false if not.
00656 //! TranslationOffset contains the offset relative to [0,1]^d, i.e. the pixel offset is Offset / H
00657 template <typename ConfiguratorType>
00658 inline bool translateAndTransformCoord ( const typename ConfiguratorType::InitType &Grid,
00659                                          const aol::auto_container<ConfiguratorType::Dim, aol::DiscreteFunctionDefault<ConfiguratorType> > &DiscrTransformation,
00660                                          const typename ConfiguratorType::ElementType &El,
00661                                          const typename ConfiguratorType::DomVecType &RefCoord,
00662                                          const typename ConfiguratorType::VecType &TranslationOffset,
00663                                          typename ConfiguratorType::ElementType &TransformedEl,
00664                                          typename ConfiguratorType::DomVecType &TransformedLocalCoord ) {
00665   typename ConfiguratorType::VecType translated_local_coord;
00666   qc::Element translated_el;
00667   // translation by TranslationOffset -> store result in translated_el, translated_local_coord
00668   if ( qc::transformCoord<ConfiguratorType> ( Grid, El, RefCoord, TranslationOffset, translated_el, translated_local_coord ) ) {
00669     // deformation by DiscreteFunctionDefault
00670     typename ConfiguratorType::VecType offset;
00671     for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00672       offset[i] = DiscrTransformation[i].evaluate ( translated_el, translated_local_coord );
00673     }
00674     return qc::transformCoord<ConfiguratorType> ( Grid, translated_el, translated_local_coord, offset, TransformedEl, TransformedLocalCoord );
00675   } else
00676     return false;
00677 }
00678 
00679 /**
00680  *  @brief Save the MultiVector to file
00681  *  @author Jingfeng Han
00682 **/
00683 template<typename ConfiguratorType>
00684 bool SaveMultiVector ( const qc::GridDefinition &Grid,
00685                        const aol::MultiVector<typename ConfiguratorType::RealType> &Phi,
00686                        const char *Filename ) {
00687   typedef typename ConfiguratorType::RealType RealType;
00688   typedef typename ConfiguratorType::ArrayType ArrayType;
00689   ofstream of ( Filename, ios::binary );
00690   if ( !of )
00691     return false;
00692 
00693   int width = Grid.getWidth ();
00694   int depth = Grid.getGridDepth();
00695   int NumOfBytes = sizeof ( RealType );
00696   int NumOfNodes = Grid.getNumberOfNodes();
00697 
00698   if ( NumOfBytes == sizeof ( float ) ) {
00699     of << "F" << "\n" << width << "\n" << depth << "\n" << NumOfNodes << "\n";
00700     cerr << "F" << "\n" << width << "\n" << depth << "\n" << NumOfNodes << "\n";
00701   } else {
00702     of << "D" << "\n" << width << "\n" << depth << "\n" << NumOfNodes << "\n";
00703     cerr << "D" << "\n" << width << "\n" << depth << "\n" << NumOfNodes << "\n";
00704   }
00705   RealType* dummy = NULL;
00706   dummy = new RealType[ConfiguratorType::Dim*NumOfNodes];
00707   if ( !dummy ) {
00708     return false;
00709   }
00710 
00711 
00712   //ofstream tmpOf("P:/local/tmp/save.txt");
00713 
00714   qc::GridDefinition::OldFullNodeIterator fnit;
00715   qc::FastILexMapper<ConfiguratorType::Dim> mapper ( Grid );
00716   int i = 0;
00717   for ( fnit = Grid.begin(); fnit != Grid.end(); ++fnit ) {
00718 
00719     for ( int n = 0;n < ConfiguratorType::Dim;n++ ) {
00720       dummy[i++] = Phi[n][mapper.getGlobalIndex ( *fnit ) ];
00721       //tmpOf<< dummy[i-1]<<" ";
00722     }
00723 
00724   }
00725 
00726   of.write ( reinterpret_cast<char*> ( dummy ), ConfiguratorType::Dim*NumOfNodes*sizeof ( RealType ) );
00727   of.close();
00728   delete[] dummy;
00729 
00730   return true;
00731 }
00732 
00733 /**
00734  *  @brief load the MultiVector from file
00735  *  @author Jingfeng Han
00736 **/
00737 template<typename ConfiguratorType>
00738 bool LoadMultiVector ( int &LevelOfDeformation,
00739                        aol::MultiVector<typename ConfiguratorType::RealType>* &Phi,
00740                        const char *Filename ) {
00741   typedef typename ConfiguratorType::RealType RealType;
00742   typedef typename ConfiguratorType::ArrayType ArrayType;
00743   ifstream in ( Filename, ios::binary );
00744   if ( !in ) {
00745     return false;
00746   }
00747   char          tmp[256];
00748   in.get ( tmp, 255, '\n' );
00749   int NumOfBytes = sizeof ( RealType );
00750   if ( tmp[0] != 'F' && tmp[0] != 'D' ) {
00751     return false;
00752   }
00753   if ( tmp[0] == 'F' && ( NumOfBytes != sizeof ( float ) ) ) {
00754     return false;
00755   }
00756   if ( tmp[0] == 'D' && ( NumOfBytes != sizeof ( double ) ) ) {
00757     return false;
00758   }
00759   cerr << "type: " << tmp << endl;
00760   int inWidth = 0;
00761   int NumOfNodes = 0;
00762 
00763   READ_COMMENTS ( in );
00764   in >> inWidth;
00765   cerr << "width: " << inWidth << endl;
00766 
00767   READ_COMMENTS ( in );
00768   in >> LevelOfDeformation;
00769   cerr << "inDepth: " << LevelOfDeformation << endl;
00770   if ( inWidth != static_cast<int> ( pow ( 2., LevelOfDeformation ) + 1. ) ) {
00771     cerr << "end " << endl;
00772     return false;
00773   }
00774 
00775   READ_COMMENTS ( in );
00776   in >> NumOfNodes;
00777   cerr << "NumOfNodes: " << NumOfNodes << endl;
00778   if ( NumOfNodes != static_cast<int> ( pow ( pow ( 2., LevelOfDeformation ) + 1., ConfiguratorType::Dim ) ) ) {
00779     cerr << "end " << endl;
00780     return false;
00781   }
00782   in.ignore();
00783 
00784   Phi = new aol::MultiVector<RealType> ( ConfiguratorType::Dim, NumOfNodes );
00785 
00786   RealType* dummy = NULL;
00787   dummy = new RealType[ConfiguratorType::Dim*NumOfNodes];
00788   in.read ( reinterpret_cast<char*> ( dummy ),  ConfiguratorType::Dim*NumOfNodes * sizeof ( RealType ) );
00789   in.close();
00790 
00791   qc::GridDefinition::OldFullNodeIterator fnit;
00792   qc::GridDefinition grid ( LevelOfDeformation, ConfiguratorType::Dim );
00793   qc::FastILexMapper<ConfiguratorType::Dim> mapper ( grid );
00794   int i = 0;
00795 
00796   for ( fnit = grid.begin(); fnit != grid.end(); ++fnit ) {
00797 
00798     for ( int n = 0; n < ConfiguratorType::Dim; n++ ) {
00799       ( *Phi ) [n][mapper.getGlobalIndex ( *fnit ) ] = dummy[i++];
00800     }
00801 
00802   }
00803 
00804   delete[] dummy;
00805   return true;
00806 }
00807 
00808 /**
00809  * \brief This class computes via "apply(...)": \f$ \int_\Omega w(x)(d-Ax-b)^2 dx \f$ in "Dest",
00810  * where the domain \f$ \Omega \f$, the weight \f$ w \f$, and the displacement \f$ d \f$ are passed to the constructor, and
00811  * \f$ b \f$ and \f$ A \f$ are passed to "apply(...)".
00812  *
00813  * \author Wirth
00814  */
00815 template <typename ConfiguratorType>
00816 class AffineDisplacementProjectionEnergy :
00817   public aol::Op< aol::Mat<ConfiguratorType::Dim,ConfiguratorType::Dim+1,typename ConfiguratorType::RealType>, aol::Scalar<typename ConfiguratorType::RealType> > {
00818 
00819 private:
00820   typedef typename ConfiguratorType::RealType RealType;
00821   typedef typename ConfiguratorType::ArrayType VectorType;
00822 
00823   // image grid
00824   const typename ConfiguratorType::InitType &_grid;
00825   // auxiliary operator and corresponding matrix
00826   const aol::WeightedMassOp<ConfiguratorType> _massOp;
00827   qc::FastUniformGridMatrix<RealType, ConfiguratorType::Dim> _massMatrix;
00828   // displacement
00829   const aol::MultiVector<RealType> &_d;
00830 
00831 public:
00832   AffineDisplacementProjectionEnergy( const typename ConfiguratorType::InitType &Grid,
00833                                       const aol::MultiVector<RealType> &D,
00834                                       const aol::Vector<RealType> &W ) :
00835     _grid( Grid ),
00836     _massOp( Grid, W, aol::ASSEMBLED ),
00837     _massMatrix( Grid ),
00838     _d( D ) {
00839     _massOp.assembleAddMatrix( _massMatrix, 1.0 );
00840   }
00841 
00842   /**
00843    * \brief Returns \f$ \int_\Omega w(x)(d-Ax-b)^2 dx \f$ in "Dest",
00844    * where the domain \f$ \Omega \f$, the weight \f$ w \f$, and the displacement \f$ d \f$ were passed to the constructor, and
00845    * \f$ b \f$ is the last column of "Arg" and \f$ A \f$ the rest.
00846    */
00847   void applyAdd( const aol::Mat<ConfiguratorType::Dim,ConfiguratorType::Dim+1,RealType> &Arg, aol::Scalar<RealType> &Dest ) const {
00848     // obtain matrix $A+I$ and vector $b$
00849     typename ConfiguratorType::VecType shift;
00850     typename ConfiguratorType::MatType transformMatrix;
00851     Arg.getCol( ConfiguratorType::Dim, shift );
00852     Arg.getSubMatrix( 0, 0, transformMatrix );
00853     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00854       transformMatrix[i][i] += 1;
00855     // compute $-(d-Ax-b)$
00856     aol::MultiVector<RealType> displacement( _d, aol::STRUCT_COPY );
00857     (qc::DataGenerator<ConfiguratorType>( _grid )).generateAffineDisplacement( transformMatrix, shift, displacement );
00858     displacement -= _d;
00859     // return $\int_\Omega w(x)(d-Ax-b)^2 dx$
00860     aol::MultiVector<RealType> tmp( _d, aol::STRUCT_COPY );
00861     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00862       _massMatrix.apply( displacement[i], tmp[i] );
00863     Dest += displacement * tmp;
00864   }
00865 };
00866 
00867 /**
00868  * \brief This class computes via "apply(...)" the gradient of \f$ \int_\Omega w(x)(d-Ax-b)^2 dx \f$ wrt \f$ A \f$ and \f$ b \f$ in "Dest",
00869  * where the domain \f$ \Omega \f$ and the displacement \f$ d \f$ are passed to the constructor, and
00870  * \f$ b \f$ and \f$ A \f$ are passed to "apply(...)".
00871  *
00872  * \author Wirth
00873  */
00874 template <typename ConfiguratorType>
00875 class AffineDisplacementProjectionGradient :
00876   public aol::Op< aol::Mat<ConfiguratorType::Dim,ConfiguratorType::Dim+1,typename ConfiguratorType::RealType> > {
00877 
00878 private:
00879   typedef typename ConfiguratorType::RealType RealType;
00880   typedef typename ConfiguratorType::ArrayType VectorType;
00881 
00882   // image grid
00883   const typename ConfiguratorType::InitType &_grid;
00884   // auxiliary operator and corresponding matrix
00885   const aol::WeightedMassOp<ConfiguratorType> _massOp;
00886   qc::FastUniformGridMatrix<RealType, ConfiguratorType::Dim> _massMatrix;
00887   // displacement
00888   const aol::MultiVector<RealType> &_d;
00889 
00890 public:
00891   AffineDisplacementProjectionGradient( const typename ConfiguratorType::InitType &Grid,
00892                                         const aol::MultiVector<RealType> &D,
00893                                         const aol::Vector<RealType> &W ) :
00894     _grid( Grid ),
00895     _massOp( Grid, W, aol::ASSEMBLED ),
00896     _massMatrix( Grid ),
00897     _d( D ) {
00898     _massOp.assembleAddMatrix( _massMatrix, 1.0 );
00899   }
00900 
00901   /**
00902    * \brief Returns variation of \f$ \int_\Omega w(x)(d-Ax-b)^2 dx \f$ in "Dest",
00903    * where the domain \f$ \Omega \f$ and the displacement \f$ d \f$ were passed to the constructor, and
00904    * \f$ b \f$ is the last column of "Arg" and \f$ A \f$ the rest.
00905    */
00906   void applyAdd( const aol::Mat<ConfiguratorType::Dim,ConfiguratorType::Dim+1,RealType> &Arg, aol::Mat<ConfiguratorType::Dim,ConfiguratorType::Dim+1,RealType> &Dest ) const {
00907     // obtain matrix $A+I$ and vector $b$
00908     typename ConfiguratorType::VecType shift;
00909     aol::Mat<ConfiguratorType::Dim,ConfiguratorType::Dim,RealType> transformMatrix;
00910     Arg.getCol( ConfiguratorType::Dim, shift );
00911     Arg.getSubMatrix( 0, 0, transformMatrix );
00912     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00913       transformMatrix[i][i] += 1;
00914     // compute $-(d-Ax-b)$
00915     aol::MultiVector<RealType> displacement( _d, aol::STRUCT_COPY );
00916     (qc::DataGenerator<ConfiguratorType>( _grid )).generateAffineDisplacement( transformMatrix, shift, displacement );
00917     displacement -= _d;
00918     // return $\int_\Omega -2w(x)(d-Ax-b) dx$ and $\int_\Omega -2w(x)(d-Ax-b) x^T dx$
00919     aol::Vector<RealType> ones( _d[0].size() );
00920     aol::MultiVector<RealType> id( _d, aol::STRUCT_COPY );
00921     ones.setAll( 1 );
00922     (qc::DataGenerator<ConfiguratorType>( _grid )).generateIdentity( id );
00923     aol::Vector<RealType> tmp( _d[0].size() );
00924     _massMatrix.apply( ones, tmp );
00925     for ( int j = 0; j < ConfiguratorType::Dim; j++ )
00926       Dest[j][ConfiguratorType::Dim] += tmp * displacement[j];
00927     for ( int i = 0; i < ConfiguratorType::Dim; i++ ) {
00928       _massMatrix.apply( id[i], tmp );
00929       for ( int j = 0; j < ConfiguratorType::Dim; j++ )
00930         Dest[j][i] += tmp * displacement[j];
00931     }
00932     Dest *= 2;
00933   }
00934 };
00935 
00936 /**
00937  * \brief This class computes via "apply(...)": \f$ \int_\Omega w(x)(\phi-R(\alpha_i)(x-x_0)-b-x_0)^2 dx \f$ in "Dest",
00938  * where the domain \f$ \Omega \f$, the weight \f$ w \f$, and the displacement \f$ \phi-id \f$ are passed to the constructor, and
00939  * \f$ b \f$ and the angles \f$ \alpha_i \f$ for the rotation matrix \f$ R(\alpha_i) \f$ are passed to "apply(...)".
00940  * If the offset \f$ x_0 \f$ was not passed to the constructor it is taken to be zero.
00941  *
00942  * \author Wirth
00943  */
00944 template <typename ConfiguratorType>
00945 class RigidDisplacementProjectionEnergy :
00946   public aol::Op< aol::Vector<typename ConfiguratorType::RealType>, aol::Scalar<typename ConfiguratorType::RealType> > {
00947 
00948 private:
00949   typedef typename ConfiguratorType::RealType RealType;
00950 
00951   const AffineDisplacementProjectionEnergy<ConfiguratorType> _affineProjector;
00952   // the offset $x_0$
00953   const typename ConfiguratorType::VecType _offset;
00954 
00955 public:
00956   RigidDisplacementProjectionEnergy( const typename ConfiguratorType::InitType &Grid,
00957                                      const aol::MultiVector<RealType> &D,
00958                                      const aol::Vector<RealType> &W ) :
00959     _affineProjector( Grid, D, W ) {
00960   }
00961 
00962 public:
00963   RigidDisplacementProjectionEnergy( const typename ConfiguratorType::InitType &Grid,
00964                                      const aol::MultiVector<RealType> &D,
00965                                      const aol::Vector<RealType> &W,
00966                                      const typename ConfiguratorType::VecType &Offset ) :
00967     _affineProjector( Grid, D, W ),
00968     _offset( Offset ) {
00969   }
00970 
00971   /**
00972    * \brief Returns \f$ \int_\Omega w(x)(\phi-R(\alpha_i)(x-x_0)-b-x_0)^2 dx \f$ in "Dest",
00973    * where the domain \f$ \Omega \f$, the weight \f$ w \f$, and the displacement \f$ \phi-id \f$ were passed to the constructor, and
00974    * \f$ b \f$ and the angles \f$ \alpha_i \f$ for the rotation matrix \f$ R(\alpha_i) \f$ are passed as "Arg".
00975    * The first dim entries of "Arg" constitute the vector \f$ b \f$, and the rest are rotation angles for rotation around
00976    * the z-, y- and x- axis (in 3D).
00977    * If the offset \f$ x_0 \f$ was not passed to the constructor it is taken to be zero.
00978    */
00979   void applyAdd( const aol::Vector<RealType> &Arg, aol::Scalar<RealType> &Dest ) const {
00980     // obtain translation b and put it into the last column of a rectangular matrix
00981     aol::Mat<ConfiguratorType::Dim,ConfiguratorType::Dim+1,RealType> totalMatrix;
00982     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
00983       totalMatrix[i][ConfiguratorType::Dim] = Arg[i];
00984     // obtain rotation R
00985     int index1[3] = { 0, 0, 1 };
00986     int index2[3] = { 1, 2, 2 };
00987     int numberAngles = Arg.size() - ConfiguratorType::Dim;
00988     typename ConfiguratorType::MatType transformMatrix, rotationMatrix;
00989     transformMatrix.setIdentity();
00990     for ( int i = 0; i < numberAngles; i++ ) {
00991       RealType alpha = Arg[ConfiguratorType::Dim + i];
00992       rotationMatrix.setIdentity();
00993       rotationMatrix[index1[i]][index1[i]] = cos( alpha );
00994       rotationMatrix[index2[i]][index2[i]] = cos( alpha );
00995       rotationMatrix[index1[i]][index2[i]] = -sin( alpha );
00996       rotationMatrix[index2[i]][index1[i]] = sin( alpha );
00997       transformMatrix *= rotationMatrix;
00998     }
00999     // put A=R-I into the first columns of the rectangular matrix
01000     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
01001       transformMatrix[i][i] -= 1;
01002     totalMatrix.setSubMatrix( 0, 0, transformMatrix );
01003     // add $x_0-Rx_0$ to the translation
01004     typename ConfiguratorType::VecType offset;
01005     transformMatrix.mult( _offset, offset );
01006     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
01007       totalMatrix[i][ConfiguratorType::Dim] -= offset[i];
01008     // compute the energy
01009     _affineProjector.applyAdd( totalMatrix, Dest );
01010   }
01011 };
01012 
01013 /**
01014  * \brief This class computes via "apply(...)" the gradient of \f$ \int_\Omega w(x)(\phi-R(\alpha_i)(x-x_0)-b-x_0)^2 dx \f$ wrt \f$ b \f$ and the \f$ \alpha_i \f$ in "Dest",
01015  * where the domain \f$ \Omega \f$, the weight \f$ w \f$, and the displacement \f$ \phi-id \f$ are passed to the constructor, and
01016  * \f$ b \f$ and the angles \f$ \alpha_i \f$ for the rotation matrix \f$ R(\alpha_i) \f$ are passed to "apply(...)".
01017  * If the offset \f$ x_0 \f$ was not passed to the constructor it is taken to be zero.
01018  *
01019  * \author Wirth
01020  */
01021 template <typename ConfiguratorType>
01022 class RigidDisplacementProjectionGradient :
01023   public aol::Op< aol::Vector<typename ConfiguratorType::RealType>, aol::Vector<typename ConfiguratorType::RealType> > {
01024 
01025 private:
01026   typedef typename ConfiguratorType::RealType RealType;
01027 
01028   const AffineDisplacementProjectionGradient<ConfiguratorType> _affineProjectorGradient;
01029   // the offset $x_0$
01030   const typename ConfiguratorType::VecType _offset;
01031 
01032 public:
01033   RigidDisplacementProjectionGradient( const typename ConfiguratorType::InitType &Grid,
01034                                        const aol::MultiVector<RealType> &D,
01035                                        const aol::Vector<RealType> &W ) :
01036     _affineProjectorGradient( Grid, D, W ) {
01037   }
01038 
01039   RigidDisplacementProjectionGradient( const typename ConfiguratorType::InitType &Grid,
01040                                        const aol::MultiVector<RealType> &D,
01041                                        const aol::Vector<RealType> &W,
01042                                        const typename ConfiguratorType::VecType &Offset ) :
01043     _affineProjectorGradient( Grid, D, W ),
01044     _offset( Offset ) {
01045   }
01046 
01047   /**
01048    * \brief Returns gradient of \f$ \int_\Omega w(x)(\phi-R(\alpha_i)(x-x_0)-b-x_0)^2 dx \f$ wrt \f$ b \f$ and the \f$ \alpha_i \f$ in "Dest",
01049    * where the domain \f$ \Omega \f$, the weight \f$ w \f$, and the displacement \f$ \phi-id \f$ were passed to the constructor, and
01050    * \f$ b \f$ and the angles \f$ \alpha_i \f$ for the rotation matrix \f$ R(\alpha_i) \f$ are passed as "Arg".
01051    * The first dim entries of "Arg" constitute the vector \f$ b \f$, and the rest are rotation angles for rotation around
01052    * the z-, y- and x- axis (in 3D).
01053    * If the offset \f$ x_0 \f$ was not passed to the constructor it is taken to be zero.
01054    */
01055   void applyAdd( const aol::Vector<RealType> &Arg, aol::Vector<RealType> &Dest ) const {
01056     // obtain translation b and put it into the last column of a rectangular matrix
01057     aol::Mat<ConfiguratorType::Dim,ConfiguratorType::Dim+1,RealType> totalMatrix;
01058     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
01059       totalMatrix[i][ConfiguratorType::Dim] = Arg[i];
01060     // obtain rotation R
01061     int index1[3] = { 0, 0, 1 };
01062     int index2[3] = { 1, 2, 2 };
01063     int numberAngles = Arg.size() - ConfiguratorType::Dim;
01064     typename ConfiguratorType::MatType transformMatrix, derivMatrix;
01065     typename ConfiguratorType::VecType derivVector;
01066     std::vector<typename ConfiguratorType::MatType> rotationMatrices;
01067     rotationMatrices.resize( numberAngles );
01068     transformMatrix.setIdentity();
01069     for ( int i = 0; i < numberAngles; i++ ) {
01070       RealType alpha = Arg[ConfiguratorType::Dim + i];
01071       rotationMatrices[i].setIdentity();
01072       rotationMatrices[i][index1[i]][index1[i]] = cos( alpha );
01073       rotationMatrices[i][index2[i]][index2[i]] = cos( alpha );
01074       rotationMatrices[i][index1[i]][index2[i]] = -sin( alpha );
01075       rotationMatrices[i][index2[i]][index1[i]] = sin( alpha );
01076       transformMatrix *= rotationMatrices[i];
01077     }
01078     // put A=R-I into the first columns of the rectangular matrix
01079     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
01080       transformMatrix[i][i] -= 1;
01081     totalMatrix.setSubMatrix( 0, 0, transformMatrix );
01082     // add $x_0-Rx_0$ to the translation
01083     typename ConfiguratorType::VecType offset;
01084     transformMatrix.mult( _offset, offset );
01085     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
01086       totalMatrix[i][ConfiguratorType::Dim] -= offset[i];
01087     // obtain the derivative of the projection energy wrt A and b
01088     aol::Mat<ConfiguratorType::Dim,ConfiguratorType::Dim+1,RealType> deriv;
01089     _affineProjectorGradient.apply( totalMatrix, deriv );
01090     // put the derivative wrt b into the result vector
01091     deriv.getCol( ConfiguratorType::Dim, derivVector );
01092     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
01093       Dest[i] += derivVector[i];
01094     // compute the derivative wrt the rotation matrix R
01095     deriv.getSubMatrix( 0, 0, derivMatrix );
01096     transformMatrix.makeTensorProduct( derivVector, _offset );
01097     derivMatrix -= transformMatrix;
01098     // compute the derivative wrt the rotation angles and put it into the result vector
01099     for ( int i = 0; i < numberAngles; i++ ) {
01100       std::vector<typename ConfiguratorType::MatType> derivRotMatrices( rotationMatrices );
01101       RealType alpha = Arg[ConfiguratorType::Dim + i];
01102       derivRotMatrices[i].setZero();
01103       derivRotMatrices[i][index1[i]][index1[i]] = -sin( alpha );
01104       derivRotMatrices[i][index2[i]][index2[i]] = -sin( alpha );
01105       derivRotMatrices[i][index1[i]][index2[i]] = -cos( alpha );
01106       derivRotMatrices[i][index2[i]][index1[i]] = cos( alpha );
01107       transformMatrix.setIdentity();
01108       for ( int j = 0; j < numberAngles; j++ )
01109         transformMatrix *= derivRotMatrices[j];
01110       Dest[ConfiguratorType::Dim + i] += transformMatrix.dotProduct( derivMatrix );
01111     }
01112   }
01113 };
01114 
01115 /**
01116  * This class represents the weighted mass matrix \f$ \left(\int_\Omega w(\phi(x))\varphi_i(x)\varphi_j(x)dx\right)_{ij} \f$.
01117  * Here, the \f$ \varphi \f$ represent the finite element base functions, and the domain \f$ \Omega \f$ is represented by the grid
01118  * which is passed to the constructor. Furthermore, \f$ w \f$ and \f$ d \f$ are passed to the constructor as aol::Vector and
01119  * aol::MultiVector respectively, where \f$ d \f$ is the displacement such that the deformation \f$ \phi \f$ is given by \f$ d \f$+identity.
01120  *
01121  * Same usage as aol::MassOp.
01122  *
01123  * \author Wirth
01124  */
01125 template <typename ConfiguratorType, aol::GridGlobalIndexMode IndexMode = aol::QUOC_GRID_INDEX_MODE>
01126 class DeformedWeightMassOp :
01127   public aol::FELinScalarWeightedMassInterface<ConfiguratorType, DeformedWeightMassOp<ConfiguratorType, IndexMode>, IndexMode> {
01128 
01129 protected:
01130   typedef typename ConfiguratorType::RealType RealType;
01131   // the displacement $d$ in all Cartesian directions
01132   const aol::DiscreteVectorFunctionDefault< ConfiguratorType, ConfiguratorType::Dim > _d;
01133   // the weight $w$
01134   const aol::DiscreteFunctionDefault<ConfiguratorType> _w;
01135 
01136 public:
01137   DeformedWeightMassOp( const typename ConfiguratorType::InitType &Grid,
01138                         const aol::Vector<RealType> &W,
01139                         const aol::MultiVector<RealType> &D,
01140                         aol::OperatorType OpType = aol::ONTHEFLY ) :
01141     // initialise the grid
01142     aol::FELinScalarWeightedMassInterface<ConfiguratorType, DeformedWeightMassOp<ConfiguratorType, IndexMode>, IndexMode>( Grid, OpType ),
01143     _d( Grid, D ),
01144     _w( Grid, W ) {}
01145 
01146   /**
01147    * Returns \f$ w(\phi(x)) \f$ at the point $x$ specified by element, quadrature point, and local coordinates.
01148    */
01149   inline RealType getCoeff( const qc::Element &El, int QuadPoint, const typename ConfiguratorType::VecType& RefCoord ) const {
01150     // compute $\phi(x)$
01151     typename ConfiguratorType::VecType transformedLocalCoord;
01152     typename ConfiguratorType::ElementType transformedEl;
01153     if( !qc::transformCoord<ConfiguratorType> ( this->_grid, _d, El, QuadPoint, RefCoord, transformedEl, transformedLocalCoord ) )
01154       // if point $x$ was displaced outside the grid, return 0
01155       return 0.;
01156     // return $w(\phi(x))$
01157     return _w.evaluate( transformedEl, transformedLocalCoord );
01158   }
01159 };
01160 
01161 /**
01162  * This class represents the weighted mass matrix \f$ \left(\int_\Omega w(\phi(x))^2\varphi_i(x)\varphi_j(x)dx\right)_{ij} \f$.
01163  * Here, the \f$ \varphi \f$ represent the finite element base functions, and the domain \f$ \Omega \f$ is represented by the grid
01164  * which is passed to the constructor. Furthermore, \f$ w \f$ and \f$ d \f$ are passed to the constructor as aol::Vector and
01165  * aol::MultiVector respectively, where \f$ d \f$ is the displacement such that the deformation \f$ \phi \f$ is given by \f$ d \f$+identity.
01166  *
01167  * Same usage as aol::MassOp.
01168  *
01169  * \author Wirth
01170  */
01171 template <typename ConfiguratorType, aol::GridGlobalIndexMode IndexMode = aol::QUOC_GRID_INDEX_MODE>
01172 class SquaredDeformedWeightMassOp :
01173   public aol::FELinScalarWeightedMassInterface<ConfiguratorType, SquaredDeformedWeightMassOp<ConfiguratorType, IndexMode>, IndexMode> {
01174 
01175 protected:
01176   typedef typename ConfiguratorType::RealType RealType;
01177   // the displacement $d$ in all Cartesian directions
01178   const aol::DiscreteVectorFunctionDefault< ConfiguratorType, ConfiguratorType::Dim > _d;
01179   // the weight $w$ to be squared
01180   const aol::DiscreteFunctionDefault<ConfiguratorType> _w;
01181 
01182 public:
01183   SquaredDeformedWeightMassOp( const typename ConfiguratorType::InitType &Grid,
01184                                const aol::Vector<RealType> &W,
01185                                const aol::MultiVector<RealType> &D,
01186                                aol::OperatorType OpType = aol::ONTHEFLY ) :
01187     // initialise the grid
01188     aol::FELinScalarWeightedMassInterface<ConfiguratorType, SquaredDeformedWeightMassOp<ConfiguratorType, IndexMode>, IndexMode>( Grid, OpType ),
01189     _d( Grid, D ),
01190     _w( Grid, W ) {}
01191 
01192   /**
01193    * Returns \f$ w(\phi(x))^2 \f$ at the point $x$ specified by element, quadrature point, and local coordinates.
01194    */
01195   inline RealType getCoeff( const qc::Element &El, int QuadPoint, const typename ConfiguratorType::VecType& RefCoord ) const {
01196     // compute $\phi(x)$
01197     typename ConfiguratorType::VecType transformedLocalCoord;
01198     typename ConfiguratorType::ElementType transformedEl;
01199     if( !qc::transformCoord<ConfiguratorType> ( this->_grid, _d, El, QuadPoint, RefCoord, transformedEl, transformedLocalCoord ) )
01200       // if point $x$ was displaced outside the grid, return 0
01201       return 0.;
01202     // return $w(\phi(x))^2$
01203     return aol::Sqr( _w.evaluate( transformedEl, transformedLocalCoord ) );
01204   }
01205 };
01206 
01207 /**
01208  * This class represents the weighted mass matrix \f$ \left(\int_\Omega w(\phi^{-1}(x))\varphi_i(x)\varphi_j(x)|det(\nabla\phi^{-1}(x))|dx\right)_{ij} \f$.
01209  * Here, the \f$ \varphi \f$ represent the finite element base functions, and the domain \f$ \Omega \f$ is represented by the grid
01210  * which is passed to the constructor. Furthermore, \f$ w \f$ and \f$ d \f$ are passed to the constructor as aol::Vector and
01211  * aol::MultiVector respectively, where \f$ d \f$ is the displacement such that the deformation \f$ \phi \f$ is given by \f$ d \f$+identity.
01212  *
01213  * Same usage as aol::MassOp.
01214  *
01215  * \author Wirth
01216  */
01217 template <typename ConfiguratorType, aol::GridGlobalIndexMode IndexMode = aol::QUOC_GRID_INDEX_MODE>
01218 class InvDeformedWeightMassOp :
01219   public aol::FELinScalarWeightedMassInterface<ConfiguratorType, InvDeformedWeightMassOp<ConfiguratorType, IndexMode>, IndexMode> {
01220 
01221 protected:
01222   typedef typename ConfiguratorType::RealType RealType;
01223   // the displacement $d$ in all Cartesian directions
01224   const aol::DiscreteVectorFunctionDefault< ConfiguratorType, ConfiguratorType::Dim > _d;
01225   // the weight $w$
01226   const aol::DiscreteFunctionDefault<ConfiguratorType> _w;
01227   // the inverse displacement, $\phi^{-1}$-identity, in all Cartesian directions as vector of the dofs and as multilinear interpolation
01228   aol::MultiVector<RealType> _invDispDOFs;
01229   const aol::DiscreteVectorFunctionDefault< ConfiguratorType, ConfiguratorType::Dim > _invDisp;
01230   // array encoding whether $\phi^{-1}$ is defined at the corresponding node
01231   typename qc::BitArray<ConfiguratorType::Dim> _invPhiDefined;
01232 
01233 public:
01234   InvDeformedWeightMassOp( const typename ConfiguratorType::InitType &Grid,
01235                            const aol::Vector<RealType> &W,
01236                            const aol::MultiVector<RealType> &D,
01237                            aol::OperatorType OpType = aol::ONTHEFLY ) :
01238     // initialise the grid
01239     aol::FELinScalarWeightedMassInterface<ConfiguratorType, InvDeformedWeightMassOp<ConfiguratorType, IndexMode>, IndexMode>( Grid, OpType ),
01240     _d( Grid, D ),
01241     _w( Grid, W ),
01242     _invDispDOFs( D, aol::STRUCT_COPY ),
01243     _invDisp( Grid, _invDispDOFs ),
01244     _invPhiDefined( Grid ) {
01245     // generate the identity function
01246     aol::MultiVector<RealType> identity( D, aol::STRUCT_COPY );
01247     qc::DataGenerator<ConfiguratorType>( Grid ).generateIdentity( identity );
01248     // compute $\phi^{-1}$-identity
01249     qc::TransformFunction<RealType, ConfiguratorType::Dim> transformOp( Grid );
01250     transformOp.setDeformation( D );
01251     transformOp.transform( identity, _invDispDOFs, _invPhiDefined );
01252     _invDispDOFs -= identity;
01253   }
01254 
01255   /**
01256    * Returns \f$ w\circ\phi^{-1}|det(\nabla\phi^{-1})| \f$ evaluated at the point specified by the element "El" and quadrature point "QuadPoint".
01257    */
01258   inline RealType getCoeff( const qc::Element &El, int QuadPoint, const typename ConfiguratorType::VecType& RefCoord ) const {
01259     // if $\phi^{-1}$ is not defined at this position, return 0
01260     if ( !_invPhiDefined.elementTrue( El ) )
01261       return 0.;
01262 
01263     // compute $\phi^{-1}(x)$
01264     typename ConfiguratorType::VecType transformedLocalCoord;
01265     typename ConfiguratorType::ElementType transformedEl;
01266     if ( !qc::transformCoord<ConfiguratorType> ( this->_grid, _invDisp, El, QuadPoint, RefCoord, transformedEl, transformedLocalCoord ) )
01267       return 0.;
01268 
01269     // compute $D\phi(\phi^{-1}(x))$
01270     typename ConfiguratorType::MatType dPhi;
01271     _d.evaluateGradient( transformedEl, transformedLocalCoord, dPhi );
01272     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
01273       dPhi[i][i] += 1.;
01274 
01275     // compute $w(\phi^{-1}(x))|det(D\phi^{-1}(x))|$
01276     return _w.evaluate( transformedEl, transformedLocalCoord ) / aol::Abs( dPhi.det() );
01277   }
01278 };
01279 
01280 /**
01281  * This class represents the weighted mass matrix \f$ \left(\int_\Omega w(\phi^{-1}(x))^2\varphi_i(x)\varphi_j(x)|det(\nabla\phi^{-1}(x))|dx\right)_{ij} \f$.
01282  * Here, the \f$ \varphi \f$ represent the finite element base functions, and the domain \f$ \Omega \f$ is represented by the grid
01283  * which is passed to the constructor. Furthermore, \f$ w \f$ and \f$ d \f$ are passed to the constructor as aol::Vector and
01284  * aol::MultiVector respectively, where \f$ d \f$ is the displacement such that the deformation \f$ \phi \f$ is given by \f$ d \f$+identity.
01285  *
01286  * Same usage as aol::MassOp.
01287  *
01288  * \author Wirth
01289  */
01290 template <typename ConfiguratorType, aol::GridGlobalIndexMode IndexMode = aol::QUOC_GRID_INDEX_MODE>
01291 class SquaredInvDeformedWeightMassOp :
01292   public aol::FELinScalarWeightedMassInterface<ConfiguratorType, SquaredInvDeformedWeightMassOp<ConfiguratorType, IndexMode>, IndexMode> {
01293 
01294 protected:
01295   typedef typename ConfiguratorType::RealType RealType;
01296   // the displacement $d$ in all Cartesian directions
01297   const aol::DiscreteVectorFunctionDefault< ConfiguratorType, ConfiguratorType::Dim > _d;
01298   // the weight $w$
01299   const aol::DiscreteFunctionDefault<ConfiguratorType> _w;
01300   // the inverse displacement, $\phi^{-1}$-identity, in all Cartesian directions as vector of the dofs and as multilinear interpolation
01301   aol::MultiVector<RealType> _invDispDOFs;
01302   const aol::DiscreteVectorFunctionDefault< ConfiguratorType, ConfiguratorType::Dim > _invDisp;
01303   // array encoding whether $\phi^{-1}$ is defined at the corresponding node
01304   typename qc::BitArray<ConfiguratorType::Dim> _invPhiDefined;
01305 
01306 public:
01307   SquaredInvDeformedWeightMassOp( const typename ConfiguratorType::InitType &Grid,
01308                                   const aol::Vector<RealType> &W,
01309                                   const aol::MultiVector<RealType> &D,
01310                                   aol::OperatorType OpType = aol::ONTHEFLY ) :
01311     // initialise the grid
01312     aol::FELinScalarWeightedMassInterface<ConfiguratorType, SquaredInvDeformedWeightMassOp<ConfiguratorType, IndexMode>, IndexMode>( Grid, OpType ),
01313     _d( Grid, D ),
01314     _w( Grid, W ),
01315     _invDispDOFs( D, aol::STRUCT_COPY ),
01316     _invDisp( Grid, _invDispDOFs ),
01317     _invPhiDefined( Grid ) {
01318     // generate the identity function
01319     aol::MultiVector<RealType> identity( D, aol::STRUCT_COPY );
01320     qc::DataGenerator<ConfiguratorType>( Grid ).generateIdentity( identity );
01321     // compute $\phi^{-1}$-identity
01322     qc::TransformFunction<RealType, ConfiguratorType::Dim> transformOp( Grid );
01323     transformOp.setDeformation( D );
01324     transformOp.transform( identity, _invDispDOFs, _invPhiDefined );
01325     _invDispDOFs -= identity;
01326   }
01327 
01328   /**
01329    * Returns \f$ (w\circ\phi^{-1})^2|det(\nabla\phi^{-1})| \f$ evaluated at the point specified by the element "El" and quadrature point "QuadPoint".
01330    */
01331   inline RealType getCoeff( const qc::Element &El, int QuadPoint, const typename ConfiguratorType::VecType& RefCoord ) const {
01332     // if $\phi^{-1}$ is not defined at this position, return 0
01333     if ( !_invPhiDefined.elementTrue( El ) )
01334       return 0.;
01335 
01336     // compute $\phi^{-1}(x)$
01337     typename ConfiguratorType::DomVecType transformedLocalCoord;
01338     typename ConfiguratorType::ElementType transformedEl;
01339     if ( !qc::transformCoord<ConfiguratorType> ( this->_grid, _invDisp, El, QuadPoint, RefCoord, transformedEl, transformedLocalCoord ) )
01340       return 0.;
01341 
01342     // compute $D\phi(\phi^{-1}(x))$
01343     typename ConfiguratorType::MatType dPhi;
01344     _d.evaluateGradient( transformedEl, transformedLocalCoord, dPhi );
01345     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
01346       dPhi[i][i] += 1.;
01347 
01348     // compute $w(\phi^{-1}(x))|det(D\phi^{-1}(x))|$
01349     return aol::Sqr( _w.evaluate( transformedEl, transformedLocalCoord ) / aol::Abs( dPhi.det() ) );
01350   }
01351 };
01352 
01353 //! Interface to compute \f$(\int_\Omega f(\phi,\Phi(x),x) \varphi_i\circ\Phi(x) dx)_i\f$,
01354 /** where \f$\Phi\f$ is a deformation on \f$\Omega\f$, \f$\varphi_i\f$ is the \f$i\f$th FE basis function,
01355  *  and \f$\phi\f$ is the argument of the operator. If \f$\Phi(x)\f$ comes to lie outside \f$\Omega\f$,
01356  *  it is projected onto the nearest point in \f$\partial\Omega\f$.
01357  *
01358  *  \author Wirth
01359  */
01360 template <typename ConfiguratorType, typename Imp>
01361 class FENonlinDeformOpInterface :
01362   public aol::FENonlinOpInterfaceWithMovedTestFunctions<ConfiguratorType, FENonlinDeformOpInterface<ConfiguratorType, Imp> > {
01363 
01364 protected:
01365   typedef typename ConfiguratorType::RealType RealType;
01366 
01367   // the deformation \Phi
01368   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement;
01369 
01370 public:
01371   explicit FENonlinDeformOpInterface( const typename ConfiguratorType::InitType &Grid, const aol::MultiVector<RealType> &Displacement ) :
01372     aol::FENonlinOpInterfaceWithMovedTestFunctions<ConfiguratorType, FENonlinDeformOpInterface<ConfiguratorType, Imp> >( Grid ),
01373     _displacement( Grid, Displacement ) {}
01374 
01375   //! Represents \f$f(\phi,\Phi(x),x)\f$ and has to be implemented in the derived class.
01376   void getNonlinearity ( const aol::DiscreteFunctionDefault<ConfiguratorType> &DiscFunc,
01377                          const typename ConfiguratorType::ElementType &El, int QuadPoint, const typename ConfiguratorType::DomVecType &LocalCoord,
01378                          const typename ConfiguratorType::ElementType &TransformedEl, const typename ConfiguratorType::DomVecType &TransformedLocalCoord,
01379                          typename ConfiguratorType::RealType &NL ) const {
01380     this->asImp().getNonlinearity ( DiscFunc, El, QuadPoint, LocalCoord, TransformedEl, TransformedLocalCoord, NL );
01381   }
01382 
01383   //! Returns \f$\Phi(x)\f$.
01384   bool getTransformedPosition ( const aol::DiscreteFunctionDefault<ConfiguratorType> &/*DiscFunc*/,
01385                                 const typename ConfiguratorType::ElementType &El,
01386                                 int QuadPoint, const typename ConfiguratorType::DomVecType &RefCoord,
01387                                 typename ConfiguratorType::ElementType &TransformedEl,
01388                                 typename ConfiguratorType::DomVecType &TransformedLocalCoord ) const {
01389     qc::transformAndClipCoord<ConfiguratorType>( *this->_config, _displacement, El, QuadPoint, RefCoord, TransformedEl, TransformedLocalCoord );
01390     return true;
01391   }
01392 
01393 protected:
01394   // Barton-Nackman trick
01395   inline Imp& asImp() { return static_cast<Imp&> ( *this ); }
01396   inline const Imp& asImp() const { return static_cast<const Imp&> ( *this ); }
01397 };
01398 
01399 //! Interface to compute \f$\int_\Omega f(\phi,\Phi(x),x)dx\f$,
01400 /** where \f$\Phi\f$ is a deformation on \f$\Omega\f$ and \f$\phi\f$ is the argument of the operator.
01401  *  If \f$\Phi(x)\f$ comes to lie outside \f$\Omega\f$,
01402  *  it is projected onto the nearest point in \f$\partial\Omega\f$.
01403  *
01404  *  \author Wirth
01405  */
01406 template <typename ConfiguratorType, typename Imp, int NumComponents = ConfiguratorType::Dim >
01407 class FENonlinDeformIntegrationVectorInterface
01408       : public aol::Op<aol::MultiVector<typename ConfiguratorType::RealType>, aol::Scalar<typename ConfiguratorType::RealType> > {
01409 
01410 protected:
01411   typedef typename ConfiguratorType::RealType RealType;
01412 
01413   aol::DeleteFlagPointer<const ConfiguratorType> _config;
01414   // the deformation \Phi
01415   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement;
01416 
01417 public:
01418   explicit FENonlinDeformIntegrationVectorInterface ( const typename ConfiguratorType::InitType &Grid, const aol::MultiVector<RealType> &Displacement ) :
01419     _config ( new ConfiguratorType ( Grid ), true ), _displacement( Grid, Displacement ) {}
01420 
01421   explicit FENonlinDeformIntegrationVectorInterface ( const ConfiguratorType &Config, const aol::MultiVector<RealType> &Displacement ) :
01422     _config ( &Config, false ), _displacement( _config->getInitializer(), Displacement ) {}
01423 
01424   void applyAdd ( const aol::MultiVector<RealType> &Arg, aol::Scalar<RealType> &Dest ) const {
01425 
01426     const aol::DiscreteVectorFunctionDefault<ConfiguratorType,NumComponents> arg( _config->getInitializer(), Arg );
01427 
01428     typename ConfiguratorType::ElementType transformedEl;
01429     typename ConfiguratorType::DomVecType  transformedCoord;
01430 
01431     typedef typename ConfiguratorType::ElementIteratorType IteratorType;
01432     const typename IteratorType::EndType end_it = _config->end();
01433 
01434     // traverse the elements of the grid and on each element the quadrature points
01435     for ( IteratorType it = _config->begin(); it != end_it; ++it ) {
01436       const typename ConfiguratorType::BaseFuncSetType &bfs = _config->getBaseFunctionSet ( *it );
01437       const int numQuadPoints = bfs.numQuadPoints( );
01438 
01439       RealType elIntegral = aol::ZTrait<RealType>::zero;
01440       for ( int q = 0; q < numQuadPoints; q++ ) {
01441         // compute \Phi(x) ("transformedEl", "transformedCoord")
01442         typename ConfiguratorType::DomVecType localCoord( bfs.getRefCoord ( q ) );
01443         qc::transformAndClipCoord<ConfiguratorType>( *_config, _displacement, *it, q, localCoord, transformedEl, transformedCoord );
01444 
01445         // compute the integrand at the quadrature points
01446         elIntegral += this->asImp().evaluateIntegrand ( arg, *it, q, localCoord, transformedEl, transformedCoord ) * bfs.getWeight( q );
01447       }
01448 
01449       Dest += elIntegral * _config->vol( *it );
01450     }
01451   }
01452 
01453   //! Represents \f$f(\phi,\Phi(x),x)\f$ and has to be implemented in the derived class.
01454   RealType evaluateIntegrand ( const aol::DiscreteVectorFunctionDefault<ConfiguratorType,NumComponents> &DiscFuncs,
01455                                const typename ConfiguratorType::ElementType &El, int QuadPoint, const typename ConfiguratorType::DomVecType &LocalCoord,
01456                                const typename ConfiguratorType::ElementType &TransformedEl, const typename ConfiguratorType::DomVecType &TransformedLocalCoord ) const {
01457     return this->asImp().evaluateIntegrand ( DiscFuncs, El, QuadPoint, LocalCoord, TransformedEl, TransformedLocalCoord );
01458   }
01459 
01460 protected:
01461   // barton-nackman
01462   inline Imp& asImp() { return static_cast<Imp&> ( *this ); }
01463   inline const Imp& asImp() const { return static_cast<const Imp&> ( *this ); }
01464 };
01465 
01466 //! Interface to compute \f$(\int_\Omega f(\phi,\Phi(x),x):\nabla(\varphi_i\circ\Phi(x)) dx)_i\f$,
01467 /** where \f$\Phi\f$ is a deformation on \f$\Omega\f$, \f$\varphi_i\f$ is the \f$i\f$th vector-valued FE basis function,
01468  *  and \f$\phi\f$ is the argument of the operator. If \f$\Phi(x)\f$ comes to lie outside \f$\Omega\f$,
01469  *  it is projected onto the nearest point in \f$\partial\Omega\f$.
01470  *
01471  *  \author Wirth
01472  */
01473 template <typename ConfiguratorType, int NumCompArg, int NumCompDest, typename Imp>
01474 class FENonlinDeformVectorDiffOpInterface : public aol::Op<aol::MultiVector<typename ConfiguratorType::RealType> > {
01475 
01476 protected:
01477   typedef typename ConfiguratorType::RealType RealType;
01478 
01479   aol::DeleteFlagPointer<const ConfiguratorType> _config;
01480   // the deformation \Phi
01481   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement;
01482 
01483 public:
01484   explicit FENonlinDeformVectorDiffOpInterface ( const typename ConfiguratorType::InitType &Grid, const aol::MultiVector<RealType> &Displacement ) :
01485     _config ( new ConfiguratorType ( Grid ), true ), _displacement( Grid, Displacement ) {}
01486 
01487   explicit FENonlinDeformVectorDiffOpInterface ( const ConfiguratorType &Config, const aol::MultiVector<RealType> &Displacement ) :
01488     _config ( &Config, false ), _displacement( _config->getInitializer(), Displacement ) {}
01489 
01490   void applyAdd ( const aol::MultiVector<RealType> &Arg, aol::MultiVector<RealType> &Dest ) const {
01491 
01492     const aol::DiscreteVectorFunctionDefault<ConfiguratorType,NumCompArg> arg( _config->getInitializer(), Arg );
01493 
01494     typename ConfiguratorType::ElementType transformedEl;
01495     typename ConfiguratorType::DomVecType  transformedCoord;
01496     typename ConfiguratorType::MatType     dPhi;
01497     typename ConfiguratorType::VecType     grad;
01498     aol::Vec<NumCompDest, typename ConfiguratorType::RealType> integrand;
01499 
01500     typedef typename ConfiguratorType::ElementIteratorType IteratorType;
01501     const typename IteratorType::EndType end_it = _config->end();
01502 
01503     // traverse the elements of the grid and on each element the quadrature points
01504     for ( IteratorType it = _config->begin(); it != end_it; ++it ) {
01505       const typename ConfiguratorType::BaseFuncSetType &bfs = _config->getBaseFunctionSet ( *it );
01506       const int numQuadPoints = bfs.numQuadPoints( );
01507 
01508       for ( int q = 0; q < numQuadPoints; q++ ) {
01509         // compute \Phi(x) ("transformedEl", "transformedCoord") and \nabla\Phi
01510         typename ConfiguratorType::DomVecType localCoord( bfs.getRefCoord ( q ) );
01511         qc::transformAndClipCoord<ConfiguratorType>( *_config, _displacement, *it, q, localCoord, transformedEl, transformedCoord );
01512         _displacement.evaluateGradientAtQuadPoint ( *it, q, dPhi );
01513         for ( int i = 0; i < ConfiguratorType::Dim; i++ )
01514           dPhi[i][i] += 1.;
01515 
01516         // compute f(\phi,\Phi(x),x) and its product with \nabla\Phi^T
01517         aol::Mat<NumCompDest, ConfiguratorType::Dim, RealType> f, fdPhi;
01518         this->asImp().getNonlinearity ( arg, *it, q, localCoord, transformedEl, transformedCoord, f );
01519         dPhi.transpose();
01520         fdPhi.makeProduct( f, dPhi );
01521 
01522         // add f\nabla\Phi^T\nabla\phi_i to the ith position in the result
01523         const int numLocalDofs = _config->getNumLocalDofs ( transformedEl );
01524         const typename ConfiguratorType::BaseFuncSetType &tbfs = _config->getBaseFunctionSet ( transformedEl );
01525         for ( int dof = 0; dof < numLocalDofs; ++dof ) {
01526           tbfs.evaluateGradient ( dof, transformedCoord, grad );
01527           fdPhi.mult( grad, integrand );
01528           integrand *= bfs.getWeight( q ) * _config->vol( *it );
01529           for ( int j = 0; j < NumCompDest; ++j )
01530             Dest[j][_config->localToGlobal( transformedEl, dof )] += integrand[j];
01531         }
01532       }
01533     }
01534   }
01535 
01536   //! Represents \f$f(\phi,\Phi(x),x)\f$ and has to be implemented in the derived class.
01537   void getNonlinearity ( const aol::DiscreteFunctionDefault<ConfiguratorType> &DiscFunc,
01538                          const typename ConfiguratorType::ElementType &El, int QuadPoint, const typename ConfiguratorType::DomVecType &LocalCoord,
01539                          const typename ConfiguratorType::ElementType &TransformedEl, const typename ConfiguratorType::DomVecType &TransformedLocalCoord,
01540                          aol::Mat<NumCompDest, ConfiguratorType::Dim, RealType> &NL ) const {
01541     this->asImp().getNonlinearity ( DiscFunc, El, QuadPoint, LocalCoord, TransformedEl, TransformedLocalCoord, NL );
01542   }
01543 
01544 protected:
01545   // barton-nackman
01546   inline Imp& asImp() { return static_cast<Imp&> ( *this ); }
01547   inline const Imp& asImp() const { return static_cast<const Imp&> ( *this ); }
01548 };
01549 
01550 //! type signifying whether in the definition of a linear FE operator
01551 //! all basis functions are deformed or none or only the ones belonging to vectors that are multiplied from the right or left
01552 enum FELinOpBaseFunDeformModeType { BOTH, RIGHT, LEFT, NONE };
01553 
01554 inline bool leftBasisFunDeformed( const FELinOpBaseFunDeformModeType defMode ) { return !( defMode % 2 ); }
01555 
01556 inline bool rightBasisFunDeformed( const FELinOpBaseFunDeformModeType defMode ) { return ( defMode < LEFT ); }
01557 
01558 //! Interface for linear FE operators with basis functions deformed by a deformation \f$\Phi\f$.
01559 /**
01560  *  \author Wirth
01561  */
01562 template <typename ConfiguratorType, typename Imp>
01563 class FELinDeformOpInterface : public aol::Op<aol::Vector<typename ConfiguratorType::RealType> > {
01564 
01565 protected:
01566   typedef typename ConfiguratorType::RealType RealType;
01567 
01568   aol::DeleteFlagPointer<const ConfiguratorType> _config;
01569   mutable typename ConfiguratorType::MatrixType *_mat;
01570   aol::OperatorType _opType;
01571   FELinOpBaseFunDeformModeType _deformationMode;
01572   // the displacement corresponding to the deformation \Phi
01573   const aol::DiscreteVectorFunctionDefault<ConfiguratorType,ConfiguratorType::Dim> _displacement;
01574 
01575 public:
01576   explicit FELinDeformOpInterface ( const typename ConfiguratorType::InitType &Grid, const aol::MultiVector<RealType> &Displacement, aol::OperatorType OpType = aol::ONTHEFLY, FELinOpBaseFunDeformModeType DeformationMode = BOTH ) :
01577     _config ( new ConfiguratorType(Grid), true ),
01578     _mat ( NULL ),
01579     _opType ( OpType ),
01580     _deformationMode( DeformationMode ),
01581     _displacement( Grid, Displacement ) {}
01582 
01583   virtual ~FELinDeformOpInterface( ) {
01584     delete _mat;
01585   }
01586 
01587   //! Clears the assembled matrix.
01588   void reset( ) {
01589     if ( _mat )
01590       delete _mat;
01591   }
01592 
01593   void applyAdd ( const aol::Vector<RealType> &Arg, aol::Vector<RealType> &Dest ) const {
01594     switch ( _opType ) {
01595     case aol::ONTHEFLY:
01596       multiplyOnTheFly ( Arg, Dest );
01597       break;
01598     case aol::ASSEMBLED:
01599       if ( !_mat )
01600         assembleMatrix();
01601       _mat->applyAdd ( Arg, Dest );
01602       break;
01603     };
01604   }
01605 
01606   typename ConfiguratorType::MatrixType& getMatrix( ) {
01607     if ( !_mat ) {
01608       _mat = _config->createNewMatrix( );
01609       assembleMatrix( );
01610     }
01611     //return dynamic_cast<typename ConfiguratorType::MatrixType&>(*_mat);
01612     return *_mat;
01613   }
01614 
01615 protected:
01616   void multiplyOnTheFly ( const aol::Vector<RealType> &Arg, aol::Vector<RealType> &Dest ) const {
01617 
01618     typename ConfiguratorType::ElementType transformedEl;
01619     typename ConfiguratorType::DomVecType  transformedCoord;
01620     aol::Mat<ConfiguratorType::maxNumLocalDofs,ConfiguratorType::maxNumLocalDofs,RealType> localMatrix;
01621     int globalDofsR[ ConfiguratorType::maxNumLocalDofs ],
01622         globalDofsL[ ConfiguratorType::maxNumLocalDofs ];
01623 
01624     typedef typename ConfiguratorType::ElementIteratorType IteratorType;
01625     const typename IteratorType::EndType end_it = _config->end();
01626 
01627     // traverse the elements of the grid and on each element the quadrature points
01628     for ( IteratorType it = _config->begin(); it != end_it; ++it ) {
01629       const typename ConfiguratorType::BaseFuncSetType &bfs = _config->getBaseFunctionSet ( *it );
01630       const int numQuadPoints = bfs.numQuadPoints( );
01631 
01632       for ( int q = 0; q < numQuadPoints; q++ ) {
01633         // compute \Phi(x) ("transformedEl", "transformedCoord")
01634         typename ConfiguratorType::DomVecType localCoord( bfs.getRefCoord ( q ) );
01635         qc::transformAndClipCoord<ConfiguratorType>( *_config, _displacement, *it, q, localCoord, transformedEl, transformedCoord );
01636 
01637         // assemble the local matrix belonging to the current quadrature point
01638         this->asImp().prepareLocalMatrix ( *it, q, localCoord, transformedEl, transformedCoord, localMatrix );
01639 
01640         // get the global indices to the current DoFs belonging to vectors multiplied from the right
01641         const typename ConfiguratorType::ElementType& elR( rightBasisFunDeformed( _deformationMode ) ? transformedEl : *it );
01642         const int numLocalDofsR = _config->getNumLocalDofs ( elR );
01643         for ( int i = 0; i < numLocalDofsR; ++i )
01644           globalDofsR[ i ] = _config->localToGlobal ( elR, i );
01645 
01646         // get the global indices to the current DoFs belonging to vectors multiplied from the left
01647         const typename ConfiguratorType::ElementType& elL( leftBasisFunDeformed( _deformationMode ) ? transformedEl : *it );
01648         const int numLocalDofsL = _config->getNumLocalDofs ( elL );
01649         for ( int i = 0; i < numLocalDofsL; ++i )
01650           globalDofsL[ i ] = _config->localToGlobal ( elL, i );
01651 
01652         // add the locally computed value to the global result
01653         for ( int i = 0; i < numLocalDofsL; ++i )
01654           for ( int j = 0; j < numLocalDofsR; ++j )
01655             Dest[ globalDofsL[ i ] ] += localMatrix[ i ][ j ] * Arg[ globalDofsR[ j ] ] ;
01656       }
01657     }
01658   }
01659 
01660   void assembleMatrix( ) const {
01661     if ( _mat )
01662       delete _mat;
01663     _mat = _config->createNewMatrix( );
01664     assembleAddMatrix ( *_mat );
01665   }
01666 
01667 public:
01668   /** (this assembled matrix * Factor) is added to Mat  */
01669   template <typename MatrixType>
01670   void assembleAddMatrix ( MatrixType &Mat, const RealType Factor = aol::NumberTrait<RealType>::one ) const {
01671 
01672     typename ConfiguratorType::ElementType transformedEl;
01673     typename ConfiguratorType::DomVecType  transformedCoord;
01674     aol::Mat<ConfiguratorType::maxNumLocalDofs,ConfiguratorType::maxNumLocalDofs,RealType> localMatrix;
01675     int globalDofsR[ ConfiguratorType::maxNumLocalDofs ],
01676         globalDofsL[ ConfiguratorType::maxNumLocalDofs ];
01677 
01678     typedef typename ConfiguratorType::ElementIteratorType IteratorType;
01679     const typename IteratorType::EndType end_it = _config->end();
01680 
01681     // traverse the elements of the grid and on each element the quadrature points
01682     for ( IteratorType it = _config->begin(); it != end_it; ++it ) {
01683       const typename ConfiguratorType::BaseFuncSetType &bfs = _config->getBaseFunctionSet ( *it );
01684       const int numQuadPoints = bfs.numQuadPoints( );
01685 
01686       for ( int q = 0; q < numQuadPoints; q++ ) {
01687         // compute \Phi(x) ("transformedEl", "transformedCoord")
01688         typename ConfiguratorType::DomVecType localCoord( bfs.getRefCoord ( q ) );
01689         qc::transformAndClipCoord<ConfiguratorType>( *_config, _displacement, *it, q, localCoord, transformedEl, transformedCoord );
01690 
01691         // assemble the local matrix belonging to the current quadrature point
01692         this->asImp().prepareLocalMatrix ( *it, q, localCoord, transformedEl, transformedCoord, localMatrix );
01693 
01694         // get the global indices to the current DoFs belonging to vectors multiplied from the right
01695         const typename ConfiguratorType::ElementType& elR( rightBasisFunDeformed( _deformationMode ) ? transformedEl : *it );
01696         const int numLocalDofsR = _config->getNumLocalDofs ( elR );
01697         for ( int i = 0; i < numLocalDofsR; ++i )
01698           globalDofsR[ i ] = _config->localToGlobal ( elR, i );
01699 
01700         // get the global indices to the current DoFs belonging to vectors multiplied from the left
01701         const typename ConfiguratorType::ElementType& elL( leftBasisFunDeformed( _deformationMode ) ? transformedEl : *it );
01702         const int numLocalDofsL = _config->getNumLocalDofs ( elL );
01703         for ( int i = 0; i < numLocalDofsL; ++i )
01704           globalDofsL[ i ] = _config->localToGlobal ( elL, i );
01705 
01706         // add the locally assembled matrix to the globally assembled matrix
01707         for ( int i = 0; i < numLocalDofsL; ++i )
01708           for ( int j = 0; j < numLocalDofsR; ++j )
01709             Mat.add ( globalDofsL[ i ], globalDofsR[ j ], Factor * localMatrix [ i ][ j ] );
01710       }
01711     }
01712   }
01713 
01714   template <typename MatrixType>
01715   void assembleAddMatrix ( const aol::BitVector &Mask, MatrixType &Mat, const RealType Factor = aol::NumberTrait<RealType>::one, const bool SetMaskedRowsToIdentity = true ) const {
01716 
01717     typename ConfiguratorType::ElementType transformedEl;
01718     typename ConfiguratorType::DomVecType  transformedCoord;
01719     aol::Mat<ConfiguratorType::maxNumLocalDofs,ConfiguratorType::maxNumLocalDofs,RealType> localMatrix;
01720     int globalDofsR[ ConfiguratorType::maxNumLocalDofs ],
01721         globalDofsL[ ConfiguratorType::maxNumLocalDofs ];
01722 
01723     typedef typename ConfiguratorType::ElementIteratorType IteratorType;
01724     const typename IteratorType::EndType end_it = _config->end();
01725 
01726     // traverse the elements of the grid and on each element the quadrature points
01727     for ( IteratorType it = _config->begin(); it != end_it; ++it ) {
01728       const typename ConfiguratorType::BaseFuncSetType &bfs = _config->getBaseFunctionSet ( *it );
01729       const int numQuadPoints = bfs.numQuadPoints( );
01730 
01731       for ( int q = 0; q < numQuadPoints; q++ ) {
01732         // compute \Phi(x) ("transformedEl", "transformedCoord")
01733         typename ConfiguratorType::DomVecType localCoord( bfs.getRefCoord ( q ) );
01734         qc::transformAndClipCoord<ConfiguratorType>( *_config, _displacement, *it, q, localCoord, transformedEl, transformedCoord );
01735 
01736         // assemble the local matrix belonging to the current quadrature point
01737         this->asImp().prepareLocalMatrix ( *it, q, localCoord, transformedEl, transformedCoord, localMatrix );
01738 
01739         // get the global indices to the current DoFs belonging to vectors multiplied from the right
01740         const typename ConfiguratorType::ElementType& elR( rightBasisFunDeformed( _deformationMode ) ? transformedEl : *it );
01741         const int numLocalDofsR = _config->getNumLocalDofs ( elR );
01742         for ( int i = 0; i < numLocalDofsR; ++i )
01743           globalDofsR[ i ] = _config->localToGlobal ( elR, i );
01744 
01745         // get the global indices to the current DoFs belonging to vectors multiplied from the left
01746         const typename ConfiguratorType::ElementType& elL( leftBasisFunDeformed( _deformationMode ) ? transformedEl : *it );
01747         const int numLocalDofsL = _config->getNumLocalDofs ( elL );
01748         for ( int i = 0; i < numLocalDofsL; ++i )
01749           globalDofsL[ i ] = _config->localToGlobal ( elL, i );
01750 
01751         // add the locally assembled matrix to the globally assembled matrix,
01752         // but write the ith row only if the ith node is not a Dirichlet node and the jth column only if the jth node is not a Dirichlet node
01753         for ( int i = 0; i < numLocalDofsL; ++i )
01754           if ( !Mask[ globalDofsL[ i ] ] )
01755             for ( int j = 0; j < numLocalDofsR; ++j )
01756               if ( !Mask[ globalDofsR[ j ] ] )
01757                 Mat.add ( globalDofsL[ i ], globalDofsR[ j ], Factor * localMatrix [ i ][ j ] );
01758       }
01759     }
01760 
01761     // set ones on the diagonal for Dirichlet nodes
01762     if ( SetMaskedRowsToIdentity )
01763       for ( int i = 0; i < Mask.size(); i++ )
01764         if ( Mask[i] )
01765           Mat.add( i, i, Factor );
01766   }
01767 
01768 protected:
01769   // barton-nackman
01770   inline Imp& asImp() { return static_cast<Imp&> ( *this ); }
01771   inline const Imp& asImp() const { return static_cast<const Imp&> ( *this ); }
01772 };
01773 
01774 //! Interface to compute \f$(\int_\Omega w(x,\Phi(x))(\varphi_i\circ\Phi)(\varphi_j\circ\Phi) dx)_{ij}\f$,
01775 /** where \f$\Phi\f$ is a deformation on \f$\Omega\f$, \f$\varphi_i\f$ is the \f$i\f$th FE basis function,
01776  *  and the scalar weight \f$w\f$ has to be provided in derived classes. If \f$\Phi(x)\f$ comes to lie outside \f$\Omega\f$,
01777  *  it is projected onto the nearest point in \f$\partial\Omega\f$.
01778  *
01779  *  \author Wirth
01780  */
01781 template <typename ConfiguratorType, typename Imp>
01782 class FELinDeformScalarWeightedMassInterface :
01783   public FELinDeformOpInterface< ConfiguratorType, FELinDeformScalarWeightedMassInterface<ConfiguratorType, Imp> > {
01784 
01785 protected:
01786   typedef typename ConfiguratorType::RealType   RealType;
01787   typedef typename ConfiguratorType::DomVecType DomVecType;
01788 
01789 public:
01790   explicit FELinDeformScalarWeightedMassInterface ( const typename ConfiguratorType::InitType &Grid, const aol::MultiVector<RealType> &Displacement, aol::OperatorType OpType = aol::ONTHEFLY ) :
01791     FELinDeformOpInterface<ConfiguratorType, FELinDeformScalarWeightedMassInterface<ConfiguratorType, Imp> > ( Grid, Displacement, OpType ) {}
01792 
01793   //! This function has to be provided in the implementation (derived class) of the interface.
01794   inline RealType getCoeff ( const typename ConfiguratorType::ElementType &El, int QuadPoint, const DomVecType &LocalCoord,
01795                              const typename ConfiguratorType::ElementType &TransformedEl, const DomVecType &TransformedLocalCoord ) const {
01796     return this->asImp().getCoeff ( El, QuadPoint, LocalCoord, TransformedEl, TransformedLocalCoord );
01797   }
01798 
01799   //! Performs the numerical quadrature of the bilinear form and saves the values locally.
01800   inline void prepareLocalMatrix ( const typename ConfiguratorType::ElementType &El, int QuadPoint, const DomVecType &LocalCoord,
01801                                    const typename ConfiguratorType::ElementType &TransformedEl, const DomVecType &TransformedLocalCoord,
01802                                    aol::Mat<ConfiguratorType::maxNumLocalDofs, ConfiguratorType::maxNumLocalDofs, RealType> &LocalMatrix ) const {
01803 
01804     aol::Vec<ConfiguratorType::maxNumLocalDofs,RealType> basisFunctionValues;
01805 
01806     // compute the weight and the basis function values
01807     const RealType coeff = getCoeff ( El, QuadPoint, LocalCoord, TransformedEl, TransformedLocalCoord );
01808     const int numDofs = this->_config->getNumLocalDofs ( TransformedEl );
01809     const typename ConfiguratorType::BaseFuncSetType &tbfs = this->_config->getBaseFunctionSet ( TransformedEl );
01810     for ( int i = 0; i < numDofs; i++ )
01811       basisFunctionValues[i] = tbfs.evaluate( i, TransformedLocalCoord );
01812 
01813     // assemble the local mass matrix
01814     const typename ConfiguratorType::BaseFuncSetType &bfs = this->_config->getBaseFunctionSet ( El );
01815     LocalMatrix.makeTensorProduct( basisFunctionValues, basisFunctionValues );
01816     LocalMatrix *= coeff * bfs.getWeight ( QuadPoint ) * this->_config->vol ( El );
01817   }
01818 
01819 protected:
01820   inline Imp &asImp() { return static_cast<Imp&> ( *this ); }
01821   inline const Imp &asImp() const { return static_cast<const Imp&> ( *this ); }
01822 };
01823 
01824 //! Interface to compute \f$(\int_\Omega A(x,\Phi(x))\nabla(\varphi_j\circ\Phi)\cdot\nabla(\varphi_i\circ\Phi) dx)_{ij}\f$,
01825 /** where \f$\Phi\f$ is a deformation on \f$\Omega\f$, \f$\varphi_i\f$ is the \f$i\f$th FE basis function,
01826  *  and the (asymmetric) matrix \f$A\f$ has to be provided in derived classes. If \f$\Phi(x)\f$ comes to lie outside \f$\Omega\f$,
01827  *  it is projected onto the nearest point in \f$\partial\Omega\f$.
01828  *
01829  *  \author Wirth
01830  */
01831 template <typename ConfiguratorType, typename Imp>
01832 class FELinDeformAsymMatrixWeightedStiffInterface :
01833   public FELinDeformOpInterface< ConfiguratorType,FELinDeformAsymMatrixWeightedStiffInterface<ConfiguratorType, Imp> > {
01834 
01835 protected:
01836   typedef typename ConfiguratorType::RealType   RealType;
01837   typedef typename ConfiguratorType::MatType    MatType;
01838   typedef typename ConfiguratorType::VecType    VecType;
01839   typedef typename ConfiguratorType::DomVecType DomVecType;
01840 
01841 public:
01842   explicit FELinDeformAsymMatrixWeightedStiffInterface ( const typename ConfiguratorType::InitType &Grid, const aol::MultiVector<RealType> &Displacement, aol::OperatorType OpType = aol::ONTHEFLY, FELinOpBaseFunDeformModeType DeformationMode = BOTH ) :
01843     FELinDeformOpInterface<ConfiguratorType, FELinDeformAsymMatrixWeightedStiffInterface<ConfiguratorType, Imp> > ( Grid, Displacement, OpType, DeformationMode ) {}
01844 
01845   //! This function has to be provided in the implementation (derived class) of the interface. It returns the matrix \f$A\f$.
01846   inline void getCoeffMatrix ( const typename ConfiguratorType::ElementType &El, int QuadPoint, const DomVecType &LocalCoord,
01847                                const typename ConfiguratorType::ElementType &TransformedEl, const DomVecType &TransformedLocalCoord,
01848                                MatType &Matrix ) const {
01849     this->asImp().getCoeffMatrix ( El, QuadPoint, LocalCoord, TransformedEl, TransformedLocalCoord, Matrix );
01850   }
01851 
01852   //! Performs the numerical quadrature of the bilinear form and saves the values locally.
01853   inline void prepareLocalMatrix ( const typename ConfiguratorType::ElementType &El, int QuadPoint, const DomVecType &LocalCoord,
01854                                    const typename ConfiguratorType::ElementType &TransformedEl, const DomVecType &TransformedLocalCoord,
01855                                    aol::Mat<ConfiguratorType::maxNumLocalDofs, ConfiguratorType::maxNumLocalDofs, RealType> &LocalMatrix ) const {
01856 
01857     MatType mat, dPhi;
01858     aol::Mat<ConfiguratorType::maxNumLocalDofs,ConfiguratorType::Dim,RealType> gradL, gradR, auxMat2;
01859 
01860     // compute the matrices A and \nabla\Phi and the product \nabla\Phi A \nabla\Phi^T (or corresponding product for _deformationMode!=BOTH)
01861     getCoeffMatrix ( El, QuadPoint, LocalCoord, TransformedEl, TransformedLocalCoord, mat ); // A
01862     this->_displacement.evaluateGradientAtQuadPoint ( El, QuadPoint, dPhi );                 // \nabla\Phi
01863     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
01864       dPhi[i][i] += 1.;
01865     MatType auxMat( mat );
01866     if ( leftBasisFunDeformed( this->_deformationMode ) )
01867       auxMat.makeProduct( dPhi, mat );
01868     mat = auxMat;
01869     if ( rightBasisFunDeformed( this->_deformationMode ) )
01870       mat *= dPhi.transposed();                                                              // the product
01871 
01872     // compute the gradients of the basis functions belonging to vectors multiplied from the right
01873     const typename ConfiguratorType::ElementType& elR( rightBasisFunDeformed( this->_deformationMode ) ? TransformedEl : El );
01874     const int numDofsR = this->_config->getNumLocalDofs ( elR );
01875     const typename ConfiguratorType::BaseFuncSetType &tbfsR = this->_config->getBaseFunctionSet ( elR );
01876     for ( int i = 0; i < numDofsR; i++ )
01877       tbfsR.evaluateGradient( i, rightBasisFunDeformed( this->_deformationMode ) ? TransformedLocalCoord : LocalCoord, static_cast<VecType&>( gradR[i] ) );
01878 
01879     // compute the gradients of the basis functions belonging to vectors multiplied from the left
01880     const typename ConfiguratorType::ElementType& elL( leftBasisFunDeformed( this->_deformationMode ) ? TransformedEl : El );
01881     const int numDofsL = this->_config->getNumLocalDofs ( elL );
01882     const typename ConfiguratorType::BaseFuncSetType &tbfsL = this->_config->getBaseFunctionSet ( elL );
01883     for ( int i = 0; i < numDofsL; i++ )
01884       tbfsL.evaluateGradient( i, rightBasisFunDeformed( this->_deformationMode ) ? TransformedLocalCoord : LocalCoord, static_cast<VecType&>( gradL[i] ) );
01885 
01886     // assemble the local stiffness matrix
01887     const typename ConfiguratorType::BaseFuncSetType &bfs = this->_config->getBaseFunctionSet ( El );
01888     auxMat2.makeProduct( gradL, mat );
01889     LocalMatrix.makeProduct( auxMat2, gradR.transposed() );
01890     LocalMatrix *= bfs.getWeight ( QuadPoint ) * this->_config->vol ( El );
01891   }
01892 
01893 protected:
01894   inline Imp &asImp() { return static_cast<Imp&> ( *this ); }
01895   inline const Imp &asImp() const { return static_cast<const Imp&> ( *this ); }
01896 };
01897 
01898 //! Class to compute \f$(\int_\Omega w(x)^2(\varphi_i\circ\Phi)(\varphi_j\circ\Phi) dx)_{ij}\f$,
01899 /** where \f$\Phi\f$ is a deformation on \f$\Omega\f$, \f$\varphi_i\f$ is the \f$i\f$th FE basis function,
01900  *  and the scalar weight \f$w\f$ is given as a discretized function. If \f$\Phi(x)\f$ comes to lie outside \f$\Omega\f$,
01901  *  it is projected onto the nearest point in \f$\partial\Omega\f$.
01902  *
01903  *  \author Wirth
01904  */
01905 template <typename ConfiguratorType>
01906 class SquaredWeightDeformMassOp :
01907   public FELinDeformScalarWeightedMassInterface<ConfiguratorType, SquaredWeightDeformMassOp<ConfiguratorType> > {
01908 
01909 protected:
01910   typedef typename ConfiguratorType::RealType   RealType;
01911   typedef typename ConfiguratorType::DomVecType DomVecType;
01912 
01913   // the weight to be squared
01914   const aol::DiscreteFunctionDefault<ConfiguratorType> _weight;
01915 
01916 public:
01917   explicit SquaredWeightDeformMassOp ( const typename ConfiguratorType::InitType &Grid, const aol::MultiVector<RealType> &Displacement, const aol::Vector<RealType> &Weight, aol::OperatorType OpType = aol::ONTHEFLY ) :
01918     FELinDeformScalarWeightedMassInterface<ConfiguratorType, SquaredWeightDeformMassOp<ConfiguratorType> > ( Grid, Displacement, OpType ),
01919     _weight( Grid, Weight ) {}
01920 
01921   inline RealType getCoeff ( const typename ConfiguratorType::ElementType &El, int QuadPoint, const DomVecType &/*LocalCoord*/,
01922                              const typename ConfiguratorType::ElementType &/*TransformedEl*/, const DomVecType &/*TransformedLocalCoord*/ ) const {
01923     return aol::Sqr( _weight.evaluateAtQuadPoint( El, QuadPoint ) );
01924   }
01925 };
01926 
01927 /** This class adds a (linearized or nonlinear) rigid body motion onto a given displacement such that mean displacement and moment are zero.
01928  *
01929  *  \author Wirth
01930  */
01931 template <typename ConfiguratorType>
01932 class meanZeroShiftAndRotationProjector :
01933   public aol::Op<aol::MultiVector<typename ConfiguratorType::RealType> > {
01934 
01935 protected:
01936   typedef typename ConfiguratorType::RealType RealType;
01937   const typename ConfiguratorType::InitType &_grid;
01938   const bool _linearized;
01939   aol::Vector<RealType> _integrator;
01940   aol::MultiVector<RealType> _momentIntegrator;
01941   aol::MultiVector<RealType> _identity;
01942   RealType _volume;
01943   typename ConfiguratorType::VecType _center;
01944   typename ConfiguratorType::MatType _moments;
01945 
01946   class Displacement3DMomentNormSqr :
01947     public aol::Op< aol::Vec3<RealType>, aol::Scalar<RealType> > {
01948   private:
01949     aol::Matrix33<RealType> _outerProd;
01950   public:
01951     explicit Displacement3DMomentNormSqr( const aol::Matrix33<RealType> OuterProd ) :
01952       _outerProd( OuterProd ) {}
01953     void applyAdd( const aol::Vec3<RealType> &Arg, aol::Scalar<RealType> &Dest ) const {
01954       // compute the rotation matrix from the given three angles
01955       aol::Matrix33<RealType> rotMat, auxMat;
01956       rotMat.setRotationAboutAxis( Arg[0], 0 );
01957       for ( int i = 1; i < 3; i++ ) {
01958         auxMat.setRotationAboutAxis( Arg[i], i );
01959         rotMat *= auxMat;
01960       }
01961       // compute the corresponding moment
01962       int index1[3] = { 1, 2, 0 };
01963       int index2[3] = { 2, 0, 1 };
01964       aol::Vec3<RealType> moment;
01965       for ( int i = 0; i < qc::QC_3D; i++ )
01966         for ( int j = 0; j < qc::QC_3D; j++ )
01967           moment[i] += _outerProd[index1[i]][j] * rotMat[index2[i]][j] - _outerProd[index2[i]][j] * rotMat[index1[i]][j];
01968       // return the square of the moment
01969       Dest += moment.normSqr();
01970     }
01971   };
01972 
01973   class Displacement3DMomentNormSqrGradient :
01974     public aol::Op< aol::Vec3<RealType>, aol::Vec3<RealType> > {
01975   private:
01976     aol::Matrix33<RealType> _outerProd;
01977   public:
01978     explicit Displacement3DMomentNormSqrGradient( const aol::Matrix33<RealType> OuterProd ) :
01979       _outerProd( OuterProd ) {}
01980     void applyAdd( const aol::Vec3<RealType> &Arg, aol::Vec3<RealType> &Dest ) const {
01981       // compute the rotation matrix from the given three angles
01982       aol::Matrix33<RealType> rotMat;
01983       rotMat.setIdentity();
01984       std::vector< aol::Matrix33<RealType> > axRotMat( 3 ), axRotMatDeriv( 3 );
01985       for ( int i = 0; i < 3; i++ ) {
01986         axRotMat[i].setRotationAboutAxis( Arg[i], i );
01987         rotMat *= axRotMat[i];
01988         axRotMatDeriv[i].setRotationAboutAxis( Arg[i] + ( aol::NumberTrait<long double>::pi / 2. ), i );
01989         for ( int j = 0; j < qc::QC_3D; j++ )
01990           axRotMatDeriv[i][j][j] = 0.;
01991       }
01992       // compute the corresponding moment
01993       int index1[3] = { 1, 2, 0 };
01994       int index2[3] = { 2, 0, 1 };
01995       aol::Vec3<RealType> moment;
01996       for ( int i = 0; i < qc::QC_3D; i++ )
01997         for ( int j = 0; j < qc::QC_3D; j++ )
01998           moment[i] += _outerProd[index1[i]][j] * rotMat[index2[i]][j] - _outerProd[index2[i]][j] * rotMat[index1[i]][j];
01999       // compute and return the moment derivatives
02000       for( int k = 0; k < qc::QC_3D; k++ ) {
02001         // compute derivative of rotation matrix wrt kth angle
02002         aol::Matrix33<RealType> rotMat;
02003         rotMat.setIdentity();
02004         for ( int i = 0; i < qc::QC_3D; i++ )
02005           if ( i != k )
02006             rotMat *= axRotMat[i];
02007           else
02008             rotMat *= axRotMatDeriv[i];
02009         aol::Vec3<RealType> momentDeriv;
02010         for ( int i = 0; i < qc::QC_3D; i++ )
02011           for ( int j = 0; j < qc::QC_3D; j++ )
02012             momentDeriv[i] += _outerProd[index1[i]][j] * rotMat[index2[i]][j] - _outerProd[index2[i]][j] * rotMat[index1[i]][j];
02013         Dest[k] += 2 * ( moment * momentDeriv );
02014       }
02015     }
02016   };
02017 
02018 public:
02019   explicit meanZeroShiftAndRotationProjector( const typename ConfiguratorType::InitType &Grid, const aol::Vector<RealType> &Weight, const bool Linearized = true ) :
02020     _grid( Grid ),
02021     _linearized( Linearized ),
02022     _integrator( Grid.getNumberOfNodes() ),
02023     _momentIntegrator( ConfiguratorType::Dim, Grid.getNumberOfNodes() ),
02024     _identity( ConfiguratorType::Dim, Grid.getNumberOfNodes() ) {
02025 
02026     // compute the integration operator, \int_\Omega w ... dx
02027     aol::MassOp<ConfiguratorType>( Grid ).apply( Weight, _integrator );
02028 
02029     // compute the moment integration operator, \int_O w x ... dx
02030     qc::DataGenerator<ConfiguratorType>( Grid ).generateIdentity( _identity );
02031 #ifdef _OPENMP
02032 #pragma omp parallel for
02033 #endif
02034     for ( int i = 0; i < ConfiguratorType::Dim; i++ )
02035       aol::WeightedMassOp<ConfiguratorType>( Grid, Weight ).apply( _identity[i], _momentIntegrator[i] );
02036 
02037     // compute volume, \int_\Omega w dx
02038     aol::Vector<RealType> ones( Grid.getNumberOfNodes() );
02039     ones.setAll( 1. );
02040     _volume = _integrator * ones;
02041 
02042     // compute the center of mass, x_c = \int_\Omega w x dx / \int_\Omega w dx
02043     for ( int j = 0; j < ConfiguratorType::Dim; j++ )
02044       _center[j] = _momentIntegrator[j] * ones;
02045     _center /= _volume;
02046 
02047     // compute second moments, \int_\Omega w (x-x_c)_i (x-x_c)_j dx
02048     for ( int k = 0; k < ConfiguratorType::Dim; k++ )
02049       for ( int j = 0; j < ConfiguratorType::Dim; j++ )
02050         _moments[k][j] = _momentIntegrator[k] * _identity[j] - _volume * _center[k] * _center[j];
02051   }
02052 
02053   //! computes total shift, \f$ \int_\Omega w u dx \f$
02054   typename ConfiguratorType::VecType totalShift( const aol::MultiVector<RealType> &Arg ) const {
02055     typename ConfiguratorType::VecType shift;
02056     for ( int j = 0; j < ConfiguratorType::Dim; j++ )
02057       shift[j] = _integrator * Arg[j];
02058     return shift;
02059   }
02060 
02061   //! computes total moment, \f$ \int_\Omega w x \times u dx \f$
02062   aol::Vec3<RealType> totalMoment( const aol::MultiVector<RealType> &Arg ) const {
02063     aol::Vec3<RealType> moment;
02064     moment[2] = _momentIntegrator[0] * Arg[1] - _momentIntegrator[1] * Arg[0];
02065     if ( ConfiguratorType::Dim == qc::QC_3D ) {
02066       moment[0] = _momentIntegrator[1] * Arg[2] - _momentIntegrator[2] * Arg[1];
02067       moment[1] = _momentIntegrator[2] * Arg[0] - _momentIntegrator[0] * Arg[2];
02068     }
02069     return moment;
02070   }
02071 
02072   //! computes a skew-symmetric matrix \f$ R \f$ such that \f$ \int_\Omega w x \times [R(x-x_c)] dx = m \f$,
02073   //! where \f$ x_c \f$ is \f$ \int_\Omega w x dx \f$ and \f$ m \f$ the given moment
02074   aol::Matrix33<RealType> momentMatrixLinearized( const aol::Vec3<RealType> &Moment ) const {
02075     aol::Matrix33<RealType> matrix, auxMat;
02076     if ( ConfiguratorType::Dim == qc::QC_2D ) {
02077       matrix[0][1] = Moment[2] / ( _moments[0][0] + _moments[1][1] );
02078       matrix[1][0] = -matrix[0][1];
02079     } else {
02080       auxMat[0][0] =  _moments[0][2];
02081       auxMat[0][1] = -_moments[0][1];
02082       auxMat[0][2] = -_moments[1][1] - _moments[2][2];
02083       auxMat[1][0] =  _moments[1][2];
02084       auxMat[1][1] =  _moments[0][0] + _moments[2][2];
02085       auxMat[1][2] =  _moments[0][1];
02086       auxMat[2][0] = -_moments[0][0] - _moments[1][1];
02087       auxMat[2][1] = -_moments[1][2];
02088       auxMat[2][2] =  _moments[0][2];
02089       aol::Vec3<RealType> abc;
02090       abc = auxMat.inverse() * Moment;
02091       matrix[0][1] =  abc[0];
02092       matrix[0][2] =  abc[1];
02093       matrix[1][0] = -abc[0];
02094       matrix[1][2] =  abc[2];
02095       matrix[2][0] = -abc[1];
02096       matrix[2][1] = -abc[2];
02097     }
02098     return matrix;
02099   }
02100 
02101   //! computes a rotation matrix \f$ R \f$ such that \f$ \int_\Omega w x \times (R d) dx = 0 \f$,
02102   //! where \f$ d = x + u - x_c \f$ is given in ``ShiftedDeform''
02103   typename ConfiguratorType::MatType momentMatrixNonlinear( const aol::MultiVector<RealType> &ShiftedDeform ) const {
02104     typename ConfiguratorType::MatType matrix;
02105     aol::Matrix33<RealType> outerProd;
02106     // create (\int_\Omega w x_k(u_j(x)+x_j-x_{cj}) dx)_{kj}
02107     for ( int k = 0; k < ConfiguratorType::Dim; k++ )
02108       for ( int j = 0; j < ConfiguratorType::Dim; j++ )
02109         outerProd[k][j] = _momentIntegrator[k] * ShiftedDeform[j];
02110     // compute rotation matrix such that \int_\Omega w x \times [R(x+u-x_c)+x_c-x] dx = \int_\Omega w x \times [R(x+u-x_c)] dx = 0
02111     if ( ConfiguratorType::Dim == qc::QC_2D ) {
02112       RealType aux1 = outerProd[0][0] + outerProd[1][1],
02113                aux2 = outerProd[0][1] - outerProd[1][0];
02114       RealType sin = sqrt( 1. / ( 1. + aol::Sqr( aux1 / aux2 ) ) ) * ( aux1 * aux2 < 0 ? 1 : -1 ),
02115                cos = sqrt( 1. - aol::Sqr( sin ) ); // note: \int_\Omega (w x \times u)_3 dx = \int_\Omega (w x \times ( u + x - x_c) ) dx
02116       matrix[0][1] = -sin;
02117       matrix[1][0] = sin;
02118       matrix[0][0] = cos;
02119       matrix[1][1] = cos;
02120     } else {
02121       aol::Vec3<RealType> initialAngles, angles;
02122       Displacement3DMomentNormSqr e( outerProd );
02123       Displacement3DMomentNormSqrGradient de( outerProd );
02124       aol::QuasiNewtonIteration< RealType, aol::Vec3<RealType>, aol::Vec3<RealType> > descent( e, de, 1000, 1.e-20, 10, false, NULL );
02125       descent.getNewtonInfo().setQuietMode( true );
02126       descent.setTimestepController( aol::NewtonInfo<RealType>::WOLFE );
02127       descent.apply( initialAngles, angles );
02128       aol::Matrix33<RealType> auxMat, resultMat;
02129       resultMat.setRotationAboutAxis( angles[0], 0 );
02130       for ( int i = 1; i < 3; i++ ) {
02131         auxMat.setRotationAboutAxis( angles[i], i );
02132         resultMat *= auxMat;
02133       }
02134       matrix.setSubMatrix( 0, 0, resultMat );
02135     }
02136     return matrix;
02137   }
02138 
02139   void apply( const aol::MultiVector<RealType> &Arg, aol::MultiVector<RealType> &Dest ) const {
02140     if ( &Dest != &Arg )
02141       Dest = Arg;
02142 
02143     // compute mean shift and make it zero
02144     typename ConfiguratorType::VecType shift( totalShift( Arg ) );
02145     for ( int j = 0; j < ConfiguratorType::Dim; j++ )
02146       Dest[j].addToAll( - shift[j] / _volume );
02147 
02148     // compute the moment of the displacement (with zero total shift), \int_\Omega w x \times u dx, and make it zero
02149     aol::Vec3<RealType> moment( totalMoment( Dest ) );
02150     aol::MultiVector<RealType> momentDisp( Arg, aol::DEEP_COPY );
02151     if ( _linearized ) {
02152       aol::Matrix33<RealType> matrix( momentMatrixLinearized( moment ) );
02153       typename ConfiguratorType::MatType momentMat;
02154       for ( int j = 0; j < ConfiguratorType::Dim; j++ )
02155         for ( int k = 0; k < ConfiguratorType::Dim; k++ )
02156           momentMat[j][k] = - matrix[j][k];
02157       shift = momentMat * _center;
02158       shift *= -1.;
02159       for ( int j = 0; j < ConfiguratorType::Dim; j++ )
02160         momentMat[j][j] = 1.;
02161       qc::DataGenerator<ConfiguratorType>( _grid ).generateAffineDisplacement( momentMat, shift, momentDisp );
02162     } else {
02163       // use the following: given the center x_c of \Omega and a deformation \phi and displacement u=\phi-id with zero mean,
02164       // find rotation \psi or rotation matrix R such that \phi_0=R(\phi-x_c)+x_c with displacement u_0=u+(R-I)(u+id-x_c) has zero moment
02165       // create u+id-x_c
02166       momentDisp += _identity;
02167       for ( int j = 0; j < ConfiguratorType::Dim; j++ )
02168         momentDisp[j].addToAll( - _center[j] );
02169       // find R and compute R-I
02170       typename ConfiguratorType::MatType rotMat( momentMatrixNonlinear( momentDisp ) );
02171       for ( int j = 0; j < ConfiguratorType::Dim; j++ )
02172         rotMat[j][j] -= 1.;
02173       // compute (R-I)(u+id-x_c)
02174       for ( int k = 0; k < momentDisp[0].size(); k++ ) {
02175         typename ConfiguratorType::VecType vec, newVec;
02176         for ( int j = 0; j < ConfiguratorType::Dim; j++ )
02177           vec[j] = momentDisp[j][k];
02178         newVec = rotMat * vec;
02179         for ( int j = 0; j < ConfiguratorType::Dim; j++ )
02180           momentDisp[j][k] = newVec[j];
02181       }
02182     }
02183     Dest += momentDisp;
02184   }
02185 
02186   void applyAdd( const aol::MultiVector<RealType> &/*Arg*/, aol::MultiVector<RealType> &/*Dest*/ ) const {
02187     throw aol::Exception( "applyAdd makes no sense in meanZeroShiftAndRotationProjector" );
02188   }
02189 };
02190 
02191 
02192 } // end namespace qc
02193 
02194 #endif // __DEFORMATIONS

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