QuOc

 

tpCFEStandardOp.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 __TPCFESTANDARDOP_H
00012 #define __TPCFESTANDARDOP_H
00013 
00014 #include <op.h>
00015 #include <progressBar.h>
00016 #include <quoc.h>
00017 
00018 #include <tpCFEElement.h>
00019 #include <tpCFEGrid.h>
00020 #include <tpCFEMatrices.h>
00021 
00022 
00023 namespace { // nameless
00024 // the following workaround is necessary because depending on constraint type, the constructor must or must not be given the grid.
00025 
00026 template< typename MatrixType, tpcfe::ConstraintType CT, typename NodalCoeffType >
00027 class CFENewMatrixForConfiguratorCreator {
00028 public:
00029   static MatrixType* create ( const tpcfe::CFEGrid< typename MatrixType::DataType, CT, NodalCoeffType > &Grid ) {
00030     return ( new MatrixType ( Grid ) );
00031   }
00032 };
00033 
00034 template< typename MatrixType, typename NodalCoeffType>
00035 class CFENewMatrixForConfiguratorCreator< MatrixType, tpcfe::CFE_DOMAIN, NodalCoeffType > {
00036 public:
00037   static MatrixType* create ( const tpcfe::CFEGrid< typename MatrixType::DataType, tpcfe::CFE_DOMAIN, NodalCoeffType > &Grid ) {
00038     if ( Grid.hasDomain() ) {
00039       const int numDofs = Grid.largestDofIndex();
00040       return ( new MatrixType ( numDofs, numDofs ) );
00041     } else {
00042       return ( new MatrixType ( Grid ) );
00043     }
00044   }
00045 };
00046 
00047 } // end of nameless namespace
00048 
00049 
00050 namespace tpcfe {
00051 
00052 /** Utility class that is used for on-the-fly multiplication of CFE-operators.
00053  */
00054 template < typename RealType >
00055 class CFEOnTheFlyMultiplicator {
00056 protected:
00057   const aol::Vector < RealType > &_arg;
00058   aol::Vector < RealType > &_dest;
00059 public:
00060   CFEOnTheFlyMultiplicator ( const aol::Vector < RealType > &Arg,
00061                              aol::Vector < RealType > &Dest ) : _arg ( Arg ), _dest ( Dest ) {}
00062 
00063   void processContribution ( const int row, const int col, const RealType value ) {
00064     _dest[row] += value * _arg[col];
00065   }
00066 
00067   const char* getProgressBarText ( ) const {
00068     return ( "Applying" );
00069   }
00070 };
00071 
00072 /** Utility class that is used for the assembly of CFE-matrices
00073  */
00074 template < typename RealType, typename MatrixType, ConstraintType CT >
00075 class CFEMatrixAssembler {
00076 protected:
00077   MatrixType &_mat;
00078 
00079 public:
00080   template< typename GridType >
00081   CFEMatrixAssembler ( MatrixType &Mat, const GridType & ) : _mat ( Mat ) {}
00082 
00083   void processContribution ( const int row, const int col, const RealType value ) {
00084     if ( value != aol::NumberTrait<RealType>::zero )
00085       _mat.add ( row, col, value );
00086   }
00087 
00088   const char* getProgressBarText ( ) const {
00089     return ( "Assembling" );
00090   }
00091 
00092 };
00093 
00094 /** Utility class that is used for the assembly of CFE-matrices in
00095  *  the case of the CFE_DOMAIN mode, where a redistribution of the degrees
00096  *  of freedom has been done.
00097  */
00098 template < typename RealType, typename MatrixType >
00099 class CFEMatrixAssembler< RealType, MatrixType, CFE_DOMAIN > {
00100 protected:
00101   typedef tpcfe::CFEGrid<RealType, CFE_DOMAIN> GridType;
00102 
00103   MatrixType  &_mat;
00104   const GridType    &_grid;
00105 
00106 public:
00107   CFEMatrixAssembler ( MatrixType &Mat, const GridType &Grid ) : _mat ( Mat ), _grid ( Grid ) {}
00108 
00109   void processContribution ( const int row, const int col, const RealType value ) {
00110     int rowIndex, colIndex;
00111 
00112     // if the given a row or col index is negative, the dof is a virtual node
00113     // the global index is obtained by removing the "-" sign
00114     // for all other dofs the global index is stored in _dofMap
00115     if ( row < 0 ) rowIndex = -row;
00116     else rowIndex = _grid.remappedDof ( row );
00117     if ( col < 0 ) colIndex = -col;
00118     else colIndex = _grid.remappedDof ( col );
00119 
00120     if ( value != static_cast<RealType> ( 0.0 ) ) _mat.add ( rowIndex, colIndex, value );
00121   }
00122 
00123   const char* getProgressBarText ( ) const {
00124     return ( "Assembling" );
00125   }
00126 };
00127 
00128 
00129 /** CFE configurator. This class provides access to the grid and the matrices
00130  *  which shall be used by the CFE. Here we have the realization for 3D only
00131  */
00132 template < typename _GridType, typename _MatrixType >
00133 class CFEConfigurator {
00134 public:
00135   typedef typename _GridType::RealType                       RealType;
00136   typedef _GridType                                          GridType;
00137   typedef _MatrixType                                        MatrixType;
00138 
00139   typedef tpcfe::CFEElement < RealType >                     ElementType;
00140 
00141   typedef aol::Mat<4, 4, RealType>                           TetraMatrixType;
00142 
00143   static const qc::Dimension                                 DimOfWorld = qc::QC_3D;
00144   static const short                                         maxNumLocalDofs = 8;
00145   static const tpcfe::ConstraintType                         CT = GridType::CT;
00146 
00147   /** A standard constructor
00148    */
00149   explicit CFEConfigurator ( const GridType & Grid ) : _grid ( Grid ) {}
00150 
00151   /** Give out a new matrix
00152    */
00153   MatrixType *createNewMatrix () const {
00154     return ( CFENewMatrixForConfiguratorCreator< MatrixType, CT, typename GridType::NodalCoeffType >::create ( _grid ) );
00155   }
00156 
00157   /** Return the grid
00158    */
00159   const GridType &grid() const { return _grid; }
00160 
00161 protected:
00162   const GridType & _grid;
00163 };
00164 
00165 /** Standard Interface for CFE operators. It is an extension of aol::Op and
00166  *  thus can be plugged into iterative solvers.
00167  *  New operators should not directly derive from CFELinOpInterface, because for
00168  *  most purposes the class CFEStandardOp already provides more functionality.
00169  */
00170 template < typename T_ConfiguratorType, typename Imp >
00171 class CFELinOpInterface : public aol::Op < aol::Vector < typename T_ConfiguratorType::RealType > > {
00172 public:
00173   typedef T_ConfiguratorType                                   ConfiguratorType;
00174   typedef aol::Vector < typename ConfiguratorType::RealType >  VectorType;
00175 
00176 private:
00177   typedef typename ConfiguratorType::RealType                  RealType;
00178   typedef typename ConfiguratorType::MatrixType                MatrixType;
00179   typedef typename ConfiguratorType::GridType                  GridType;
00180 
00181 protected:
00182   aol::OperatorType      _opType;
00183   const ConfiguratorType _config;
00184 
00185 public:
00186   bool  _quietMode;
00187 
00188 protected:
00189   mutable MatrixType * _mat;
00190 
00191 public:
00192   /** A constructor which takes a configurator
00193    */
00194   explicit CFELinOpInterface ( ConfiguratorType Config, aol::OperatorType OpType = aol::ONTHEFLY ) :
00195       _opType ( OpType ), _config ( Config ), _quietMode ( false ), _mat ( NULL ) { }
00196 
00197   /** This is a constructor which builds its own configurator with the given grid
00198    */
00199   explicit CFELinOpInterface ( const GridType &Grid, aol::OperatorType OpType = aol::ONTHEFLY ) :
00200       _opType ( OpType ), _config ( Grid ), _quietMode ( false ), _mat ( NULL ) { }
00201 
00202   /** Return the matrix which is stored for this configurator
00203    */
00204   MatrixType &getMatrixRef() const {
00205     if ( _opType == aol::ASSEMBLED && !_mat )
00206       assembleMatrix();
00207 
00208     if ( _mat ) {
00209       return *_mat;
00210     } else {
00211       throw aol::Exception ( "tpcfe::CFELinOpInterface::getMatrixRef: _mat not set", __FILE__, __LINE__ );
00212     }
00213   }
00214 
00215   /** Write an image which shows the non-zero pattern of the associated matrix
00216    */
00217   void writeNonZeroPattern ( const char *fileName ) const {
00218     _mat->writeNonZeroPattern ( fileName );
00219   }
00220 
00221   /** The standard applyAdd does an "on-the-fly" traversal or a matrix multiplication.
00222    *  If the matrix has not been assembled yet, it will be assembled now
00223    */
00224   void applyAdd ( const aol::Vector <RealType > &Arg,
00225                   aol::Vector < RealType > &Dest ) const {
00226     switch ( _opType ) {
00227       case aol::ONTHEFLY: {
00228           tpcfe::CFEOnTheFlyMultiplicator<RealType> otfm ( Arg, Dest );
00229           loopOverElements<tpcfe::CFEOnTheFlyMultiplicator<RealType> > ( otfm );
00230           break;
00231         }
00232       case aol::ASSEMBLED : {
00233           if ( !_mat ) assembleMatrix();
00234           _mat->applyAdd ( Arg, Dest );
00235           break;
00236         }
00237     }
00238   }
00239 
00240   /** Compute the contribution of the given element.
00241    *  This is a standard method which calls @see tpcfe::CFEStandardOp::addSubTetraContrib()
00242    *  or @see tpcfe::CFEStandardOp::computeNonInterfacedHexaContrib()
00243    *  on the derived subclasses depending on whether this element is
00244    *  interfaced or not.
00245    *  What is done with the contribution is determined by the given instance of
00246    *  ContribProcessorType. In the standard cases this
00247    *  will be either CFEOnTheFlyMultiplicator or CFEMatrixAssembler
00248    */
00249   template <typename ContribProcessorType>
00250   void computeElementContribution ( tpcfe::CFEElement<RealType> &e, ContribProcessorType &cpt ) const {
00251 
00252     asImp().preProcessLocalContribution ( e, cpt );
00253 
00254     if ( ! ( e.cfeType().representsInterfaced() ) ) {
00255       // Element is not interfaced
00256       asImp().computeNonInterfacedHexaContrib ( e, cpt );
00257     } else {
00258       // Prepare all data needed for this element
00259       e.computeAssembleData ( _config.grid() );
00260 
00261       // Iterate over tetrahedra of this element
00262       for ( tpcfe::CFETetraInElementIterator<RealType> it ( e ); it.notAtEnd(); ++it ) {
00263         switch ( it->getVolType() ) {
00264           case CFETopoTetra::P:
00265           case CFETopoTetra::N:
00266             asImp().addSubTetraContrib ( e, *it, cpt );
00267             break;
00268           case CFETopoTetra::PPPN:
00269           case CFETopoTetra::NNNP:
00270             asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00271             asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00272             asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00273             asImp().addSubTetraContrib ( e, *it, cpt );
00274             break;
00275           case CFETopoTetra::NNNPPP:
00276           case CFETopoTetra::PPPNNN:
00277             asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00278             asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00279             asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00280             asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00281             asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00282             asImp().addSubTetraContrib ( e, *it, cpt );
00283             break;
00284           default:
00285             break;
00286         }
00287       }
00288     }
00289     asImp().postProcessLocalContribution ( e, cpt );
00290   }
00291 
00292   /** Compute the contribution of the given element. Works exactly as
00293    *  computeElementContribution, but considers elements with the correct
00294    *  sign only
00295    */
00296   template <typename ContribProcessorType>
00297   void computeSignElementContribution ( tpcfe::CFEElement<RealType> &e,
00298                                         ContribProcessorType &cpt,
00299                                         const signed char sign = -1             ) const {
00300 
00301     const unsigned char pureCFEType = e.cfeType()._pureType;
00302 
00303     if ( ( sign == -1 && pureCFEType == 0x00 ) || ( sign == + 1 && pureCFEType == 0xFF ) ) return;
00304 
00305     asImp().preProcessLocalContribution ( e, cpt );
00306 
00307     if ( ( sign == + 1 && pureCFEType == 0x00 ) ||
00308          ( sign == -1 && pureCFEType == 0xFF ) ) {
00309       // Element is not interfaced
00310       asImp().computeNonInterfacedHexaContrib ( e, cpt );
00311     } else {
00312       // Prepare all data needed for this element
00313       e.computeAssembleData ( _config.grid() );
00314       // Iterate over tetrahedra of this element
00315       for ( tpcfe::CFETetraInElementIterator<RealType> it ( e ); it.notAtEnd(); ++it ) {
00316         switch ( it->getVolType() ) {
00317           case CFETopoTetra::P:
00318             if ( sign == + 1 ) asImp().addSubTetraContrib ( e, *it, cpt );
00319             break;
00320           case CFETopoTetra::N:
00321             if ( sign == - 1 ) asImp().addSubTetraContrib ( e, *it, cpt );
00322             break;
00323           case CFETopoTetra::PPPN:
00324             if ( sign == + 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00325             if ( sign == + 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00326             if ( sign == + 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00327             if ( sign == - 1 ) asImp().addSubTetraContrib ( e, *it, cpt );
00328             break;
00329           case CFETopoTetra::NNNP:
00330             if ( sign == - 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00331             if ( sign == - 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00332             if ( sign == - 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00333             if ( sign == + 1 ) asImp().addSubTetraContrib ( e, *it, cpt );
00334             break;
00335           case CFETopoTetra::NNNPPP:
00336             if ( sign == - 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00337             if ( sign == - 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00338             if ( sign == - 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00339             if ( sign == + 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00340             if ( sign == + 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00341             if ( sign == + 1 ) asImp().addSubTetraContrib ( e, *it, cpt );
00342             break;
00343           case CFETopoTetra::PPPNNN:
00344             if ( sign == + 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00345             if ( sign == + 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00346             if ( sign == + 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00347             if ( sign == - 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00348             if ( sign == - 1 ) asImp().addSubTetraContrib ( e, *it, cpt ); ++it;
00349             if ( sign == - 1 ) asImp().addSubTetraContrib ( e, *it, cpt );
00350             break;
00351           default:
00352             break;
00353         }
00354       }
00355     }
00356     asImp().postProcessLocalContribution ( e, cpt );
00357   }
00358 
00359 
00360   /** Assemble the matrix into mat
00361    */
00362   template < typename MatrixType >
00363   void assembleAddMatrix ( MatrixType & mat ) const {
00364     tpcfe::CFEMatrixAssembler<RealType, MatrixType, ConfiguratorType::CT>  assembler ( mat, _config.grid() );
00365     loopOverElements<tpcfe::CFEMatrixAssembler<RealType, MatrixType, ConfiguratorType::CT> > ( assembler );
00366   }
00367 
00368 
00369 private:
00370   template <typename ASSEMBLER_TYPE>
00371   void loopOverElements ( ASSEMBLER_TYPE &ma ) const {
00372     aol::ProgressBar<> pb ( ma.getProgressBarText() );
00373     pb.start ( _config.grid().getNumberOfElements() );
00374 
00375     if ( _config.grid().hasDomain() ) { // the grid has a domain set
00376 
00377       const qc::GridSize<qc::QC_3D> gridSize ( _config.grid() );
00378       for ( typename GridType::FullElementIterator it ( _config.grid() ); it.notAtEnd(); ++it ) {
00379         if ( !_quietMode ) pb++;
00380         tpcfe::CFEElement<RealType> el ( *it, gridSize, _config.grid().getElType ( *it ) );
00381         if ( _config.grid().getCT() == CFE_CDWI_TPOS || _config.grid().getCT() == CFE_CDWI_LIEHR ) {
00382           const int structureNo = el.cfeType()._structureNo;
00383 
00384           // integrate not interfaced elements in domain
00385           if ( ! ( el.cfeType().representsInterfaced() ) ) {
00386             asImp().computeSignElementContribution ( el, ma, -1 );
00387           } else {
00388             // integrate correct part of domain-cutted elements
00389             if ( structureNo == MAX_STRUCT_ID ) {
00390               asImp().computeSignElementContribution ( el, ma, -1 );
00391             }
00392             // integrate interface-cutted elements
00393             else {
00394               asImp().computeElementContribution ( el, ma );
00395             }
00396           }
00397         } else {
00398           asImp().computeSignElementContribution ( el, ma, -1 );
00399         }
00400       }
00401 
00402     } else {
00403 
00404       const qc::GridSize<qc::QC_3D> gridSize ( _config.grid() );
00405       for ( typename GridType::FullElementIterator it ( _config.grid() ); it.notAtEnd(); ++it ) {
00406         if ( !_quietMode ) pb++;
00407         tpcfe::CFEElement<RealType> el ( *it, gridSize, _config.grid().getElType ( *it ) );
00408         asImp().computeElementContribution ( el, ma );
00409       }
00410 
00411     }
00412     if ( !_quietMode ) pb.finish();
00413   }
00414 
00415 
00416 public:
00417   /** Make the (possibly already assembled) matrix invalid by deleting it
00418    */
00419   void invalidateMatrix() {
00420     cerr << "Invalidating matrix " << endl;
00421     if ( _mat != NULL ) delete _mat;
00422     _mat = NULL;
00423   }
00424 
00425   /** Print the matrix object
00426    */
00427   void dump ( ostream &out ) {
00428     if ( !_mat ) throw aol::Exception ( "CFELinOpInterface::dump Matrix object does not exist", __FILE__, __LINE__ );
00429     //          _mat->dump(out);
00430     _mat->dump();
00431   }
00432 
00433   /** Print one line of the matrix object
00434    */
00435   void dumpLine ( const int which, ostream &out ) {
00436     if ( !_mat ) throw aol::Exception ( "CFELinOpInterface::dump Matrix object does not exist", __FILE__, __LINE__ );
00437     _mat->dumpLine ( which, out );
00438   }
00439 
00440   /** (Re)allocate matrix for this operator
00441    */
00442   void allocateMatrix () const {
00443     if ( _mat ) delete _mat;
00444     _mat = _config.createNewMatrix ();
00445   }
00446 
00447   /** Ask the configurator to build a new matrix and then assemble it
00448    */
00449   void assembleMatrix () const {
00450     allocateMatrix();
00451     assembleAddMatrix ( *_mat );
00452   }
00453 
00454   /** Checks whether the matrices of the operators are strictly equal
00455    */
00456   bool operator== ( CFELinOpInterface &other ) const {
00457     return ( *_mat ) == other.getMatrixRef();
00458   }
00459 
00460   /** Checks whether the matrices of the operators are approximately equal, i.e. fails if
00461    *  the entries differ by more than epsilon
00462    */
00463   bool approxEqual ( const CFELinOpInterface &other, RealType epsilon ) const {
00464     return _mat->approxEqual ( other.getMatrixRef(), epsilon );
00465   }
00466 
00467   /** Return the configurator
00468    */
00469   const ConfiguratorType &getConfig() const { return _config; }
00470 
00471   /** call makeRowEntries on Matrix, for compatibility with GaussSeidelInverse
00472    */
00473   void makeRowEntries ( std::vector<typename aol::Row<RealType>::RowEntry > &vec, const int RowNum ) const {
00474     if ( _mat ) {
00475       _mat->makeRowEntries ( vec, RowNum );
00476     } else {
00477       throw aol::Exception ( "tpcfe::CFELinOpInterface: Cannot makeRowEntries: _mat not set", __FILE__, __LINE__ );
00478     }
00479   }
00480 
00481 protected:
00482   /** Barton-Nackman-Trick avoids virtual functions because of
00483    *  the static cast
00484    */
00485   inline Imp & asImp () {
00486     return static_cast < Imp & > ( *this );
00487   }
00488 
00489   inline const Imp & asImp () const {
00490     return static_cast < const Imp & > ( *this );
00491   }
00492 
00493   ~CFELinOpInterface() {
00494     if ( _mat )
00495       delete ( _mat );
00496   }
00497 
00498   CFELinOpInterface ( const CFELinOpInterface<T_ConfiguratorType, Imp> &other ) {
00499     throw aol::UnimplementedCodeException ( "tpcfe::CFELinOpInterface copy constructor not implemented", __FILE__, __LINE__ );
00500   }
00501 
00502 };
00503 
00504 
00505 /** A base class for standard operators which use precomputed matrices stored in
00506  *  lookup tables
00507  *  \ingroup tpCFE
00508  */
00509 template < typename ConfiguratorType, typename ImpType >
00510 class CFEStandardOp : public CFELinOpInterface < ConfiguratorType, ImpType > {
00511 protected:
00512   typedef typename ConfiguratorType::RealType                                                          RealType;
00513   typedef typename ConfiguratorType::GridType                                                          GridType;
00514   typedef typename ConfiguratorType::MatrixType                                                        MatrixType;
00515 
00516   typedef typename GridType::VNType                                                                    VNType;
00517   typedef aol::Mat<ConfiguratorType::maxNumLocalDofs, ConfiguratorType::maxNumLocalDofs, RealType>     LocalMatrixType;
00518   typedef typename ConfiguratorType::TetraMatrixType                                                   TetraMatrixType;
00519 
00520 protected:
00521   LocalMatrixType  _defaultLocalHexaMatrix;       //!< contains the hexahedron element matrix for a non-interfaced element. This matrix is obtained by adding up all _defaultLocalTetraMatrices for the 6 standard tetrahedra
00522   TetraMatrixType  _defaultLocalTetraMatrix[6];   //!< stores the matrix for standard tetrahedra which have volume 1./6.
00523 
00524   mutable TetraMatrixType _localTetraMatrix;      //!< stores the local sub-tetra matrix. is filled in prepareTetraMatrix
00525 
00526   const GridType &_grid;
00527 
00528 public:
00529   explicit CFEStandardOp ( const typename ConfiguratorType::GridType & Grid, const aol::OperatorType OpType = aol::ONTHEFLY ) : CFELinOpInterface< ConfiguratorType, ImpType > ( Grid, OpType ), _grid ( Grid ) {}
00530 
00531 protected:
00532   /** Create the default matrices
00533    */
00534   void init() {
00535     this->asImp().createDefaultTetraMatrices();
00536     this->asImp().buildHexaMatrix();
00537   }
00538 
00539   /** Assemble the local matrix for a whole hexahedron
00540    *  from the matrices of the standard tetrahedra
00541    */
00542   void buildHexaMatrix() {
00543 
00544     _defaultLocalHexaMatrix.setZero();
00545 
00546     for ( CFETopoTetraIterator it ( CFEType ( -1, 0 ) ); it.notAtEnd(); ++it ) {
00547       const int parent = it->getParent();
00548 
00549       for ( short i = 0; i < 4; ++i ) {
00550         const int i1 = ( *it ) ( i, 0 );
00551         for ( short j = 0; j < 4; ++j ) {
00552           const int j1 = ( *it ) ( j, 0 );
00553           _defaultLocalHexaMatrix[i1][j1] += _defaultLocalTetraMatrix[parent][i][j];
00554         }
00555       }
00556     }
00557   }
00558 
00559 public:
00560 
00561   /** Print out the local hexa matrix
00562    */
00563   void dumpLocalHexaMatrix ( ostream &out ) const {
00564     for ( short i = 0; i < 8; ++i ) {
00565       for ( short j = 0; j < 8; ++j ) {
00566         out << _defaultLocalHexaMatrix[i][j] << " ";
00567       }
00568       out << endl;
00569     }
00570   }
00571 
00572   /** Preprocess any computations. The term "local" means Hexahedron not Tetrahedron or Subtetrahedron.
00573    */
00574   template <typename ContribProcessorType>
00575   inline void preProcessLocalContribution ( const CFEElement< RealType>& /*e*/,
00576                                             ContribProcessorType& /*cpt*/ ) const {}
00577 
00578   /** Postprocess any computations. The term "local" means Hexahedron not Tetrahedron or Subtetrahedron.
00579    */
00580   template <typename ContribProcessorType>
00581   inline void postProcessLocalContribution ( const CFEElement< RealType>& /*e*/,
00582                                              ContribProcessorType& /*cpt*/ ) const {}
00583 
00584   /** copy the local matrix for an element
00585    */
00586   template <typename ContribProcessorType>
00587   inline void computeNonInterfacedHexaContrib ( const CFEElement< RealType> &e,
00588                                                 ContribProcessorType &cpt ) const {
00589     for ( short i = 0; i < 8; ++i ) {
00590       const int i1 = e.globalIndex ( i );
00591       for ( short j = 0; j < 8; ++j ) {
00592         const int j1 = e.globalIndex ( j );
00593         cpt.processContribution ( i1, j1, _defaultLocalHexaMatrix[i][j] );
00594       }
00595     }
00596   }
00597 
00598   /** Return the matrix contribution the coupling i, j would give for standard
00599    *  basis functions i and j.
00600    *  In this default implementation it is assumed that the sub-tetra-matrix is
00601    *  just a scaled version of the parental std-tetra-matrix. Note that this is not
00602    *  satisfied in all cases. It holds for the standard mass matrices but not for the
00603    *  standard stiffness matrices!
00604    */
00605   inline const RealType getTetraContrib ( const CFEElement< RealType >& /*e*/, const CFETetra<RealType>& /*t*/,
00606                                           const int /*parent*/, const int i, const int j ) const {
00607     return _localTetraMatrix[i][j];
00608   }
00609 
00610   /** Add the constribution of a sub tetrahedron according to the construction of cfe weights
00611    */
00612   template <typename ContribProcessorType>
00613   inline void addSubTetraContrib ( const CFEElement< RealType > &e,
00614                                    const CFETetra<RealType> & t,
00615                                    ContribProcessorType &cpt        ) const {
00616 
00617     const short parent = t.getParent();
00618     this->asImp().preprocessLocalTetraMatrix ( e, t );
00619 
00620     for ( short i = 0; i < 4; ++i ) {
00621       unsigned char sI = t.isVirtualNode ( i );
00622       for ( short j = 0; j < 4; ++j ) {
00623         unsigned char sJ = t.isVirtualNode ( j );
00624         sJ <<= 1;
00625 
00626         const RealType contrib = this->asImp().getTetraContrib ( e, t, parent, i, j );
00627         switch ( sI | sJ ) {
00628           case 0: { // Both nodes are no virtual nodes
00629             const int i1 = e.globalIndex ( t ( i, 0 ) );
00630             const int j1 = e.globalIndex ( t ( j, 0 ) );  // Global index for node j
00631 
00632             cpt.processContribution ( i1, j1, contrib );
00633             break;
00634           }
00635           case 1: { // Node i is a virtual node
00636             const int j1 = e.globalIndex ( t ( j, 0 ) );  // Global index for node j
00637 
00638             // Global indices of nodes defining the edge on which the virtual node lies
00639             const int en1 = e.globalIndex ( t ( i, 0 ) );
00640             const int en2 = e.globalIndex ( t ( i, 1 ) );
00641 
00642             // Get the virtual node object from the grid
00643             const VNType& n = this->_grid.getVirtualNodeRef ( en1, en2 );
00644 
00645             //Ignore basis function if virtual node is a dirichlet node
00646             if ( n._isDirichlet ) { continue; }
00647 
00648             const int numIConstraints = n.numConstraints();
00649             for ( int l = 0; l < numIConstraints; ++l ) {
00650               const RealType coeff = n.weight ( l );
00651               const int i1 = n.constraint ( l );
00652 
00653               cpt.processContribution ( i1, j1, contrib*coeff );
00654             }
00655             //cerr << "1";
00656             break;
00657           }
00658           case 2: { // Node j is a virtual node
00659             const int i1 = e.globalIndex ( t ( i, 0 ) );
00660 
00661             // Global indices of nodes defining the edge on which the virtual node lies
00662             const int en1 = e.globalIndex ( t ( j, 0 ) );
00663             const int en2 = e.globalIndex ( t ( j, 1 ) );
00664 
00665             // Get the virtual node object from the grid
00666             const VNType& n = this->_grid.getVirtualNodeRef ( en1, en2 );
00667 
00668             //Ignore basis function if virtual node is a dirichlet node
00669             if ( n._isDirichlet ) { continue; }
00670 
00671             const int numJConstraints = n.numConstraints();
00672             for ( int l = 0; l < numJConstraints; ++l ) {
00673               const RealType coeff = n.weight ( l );
00674               const int j1 = n.constraint ( l );
00675 
00676               cpt.processContribution ( i1, j1, contrib*coeff );
00677             }
00678             //cerr << "2";
00679             break;
00680           }
00681           case 3: { // Both nodes are virtual nodes
00682             // Global indices of nodes defining the edge on which the virtual node lies
00683             const int eIn1 = e.globalIndex ( t ( i, 0 ) );
00684             const int eIn2 = e.globalIndex ( t ( i, 1 ) );
00685 
00686             const int eJn1 = e.globalIndex ( t ( j, 0 ) );
00687             const int eJn2 = e.globalIndex ( t ( j, 1 ) );
00688 
00689             // Get the virtual node objects from the grid
00690             const VNType& nI = this->_grid.getVirtualNodeRef ( eIn1, eIn2 );
00691             const VNType& nJ = this->_grid.getVirtualNodeRef ( eJn1, eJn2 );
00692 
00693             //Ignore basis function if virtual node is a dirichlet node
00694             if ( nI._isDirichlet ) continue;
00695             if ( nJ._isDirichlet ) continue;
00696 
00697             const int numIConstraints = nI.numConstraints();
00698             const int numJConstraints = nJ.numConstraints();
00699             for ( int k = 0; k < numIConstraints; ++k ) {
00700               const RealType coeffI = nI.weight ( k );
00701               const int i1 = nI.constraint ( k );
00702               for ( int l = 0; l < numJConstraints; ++l ) {
00703                 const RealType coeffJ = nJ.weight ( l );
00704                 const int j1 = nJ.constraint ( l );
00705 
00706                 cpt.processContribution ( i1, j1, contrib*coeffI*coeffJ );
00707               }
00708             }
00709             //cerr << "3";
00710             break;
00711           }
00712           default:
00713             cerr << static_cast< int > ( sI | sJ ) << endl;
00714             throw aol::Exception ( "CFELinOpInterface::addSubTetraContrib: Cannot be!", __FILE__, __LINE__ );
00715         }
00716       }
00717     }
00718   }
00719 
00720   /** Do any computations which have to be performed before the assembly loop.
00721    *  E.g. it is more efficient to assemble the local tetra matrix here and
00722    *  distribute it only in the assembly loop.
00723    *  Not virtual because we are using Barton-Nackman
00724    */
00725   void preprocessLocalTetraMatrix ( const CFEElement< RealType > &e, const CFETetra<RealType> & t ) const {}
00726 
00727   /** Do any computations -- whatever that may be -- which have to be performed after the assembly loop.
00728    *  Not virtual because we are using Barton-Nackman
00729    */
00730   void postprocessLocalTetraMatrix ( const CFEElement< RealType > &e, const CFETetra<RealType> & t ) const {}
00731 
00732   // end class CFEStandardOp
00733 };
00734 
00735 
00736 
00737 /** Composite Finite Element Mass Matrix
00738  */
00739 template < typename _ConfiguratorType >
00740 class CFEMassOp : public CFEStandardOp < _ConfiguratorType, CFEMassOp < _ConfiguratorType > > {
00741 public:
00742   typedef CFEStandardOp < _ConfiguratorType, CFEMassOp < _ConfiguratorType > > BaseType;
00743   typedef typename BaseType::ConfiguratorType ConfiguratorType;
00744   typedef typename BaseType::VectorType       VectorType;
00745 
00746 protected:
00747   typedef typename ConfiguratorType::RealType RealType;
00748 
00749 
00750 public:
00751   explicit CFEMassOp ( const typename ConfiguratorType::GridType & Grid, aol::OperatorType OpType = aol::ONTHEFLY ) :
00752       CFEStandardOp < ConfiguratorType,  CFEMassOp < ConfiguratorType > > ( Grid, OpType ) {
00753     this->init();
00754   }
00755 
00756   /** Create the standard matrices of the non-interfaced standard tetrahedra
00757    */
00758   void createDefaultTetraMatrices() {
00759     const RealType stdTetraVolFactor = aol::Cub ( this->_grid.H() );
00760 
00761     for ( short l = 0; l < 6; ++l )
00762       for ( short i = 0; i < 4; ++i )
00763         for ( short j = 0; j < 4; ++j )
00764           this->_defaultLocalTetraMatrix[l][i][j] = tpcfe::CFELookup<RealType>::_localMassMatrixRefTetra[i][j] * stdTetraVolFactor;
00765   }
00766 
00767   /** fill the matrix for a sub-tetrahedron
00768    */
00769   void preprocessLocalTetraMatrix ( const CFEElement< RealType > &/*e*/, const CFETetra<RealType> & t ) const {
00770     const RealType volumeFactor = t.getVolume() / tpcfe::CFELookup<RealType>::_stdTetraVolume * aol::Cub ( this->_grid.H() );
00771 
00772     for ( short i = 0; i < 4; ++i )
00773       for ( short j = 0; j < 4; ++j )
00774         this->_localTetraMatrix[i][j] = tpcfe::CFELookup<RealType>::_localMassMatrixRefTetra[i][j] * volumeFactor;
00775   }
00776 
00777   // end class CFEMassOp
00778 };
00779 
00780 
00781 
00782 /** Composite Finite Element Stiffness Matrix
00783  */
00784 template < typename _ConfiguratorType >
00785 class CFEStiffOp : public CFEStandardOp < _ConfiguratorType, CFEStiffOp < _ConfiguratorType > > {
00786 public:
00787   typedef CFEStandardOp < _ConfiguratorType, CFEStiffOp < _ConfiguratorType > > BaseType;
00788   typedef typename BaseType::ConfiguratorType ConfiguratorType;
00789   typedef typename BaseType::VectorType       VectorType;
00790 
00791 protected:
00792   typedef typename ConfiguratorType::RealType RealType;
00793 
00794 public:
00795   explicit CFEStiffOp ( const typename ConfiguratorType::GridType & Grid, const aol::OperatorType OpType = aol::ONTHEFLY ) :
00796       CFEStandardOp < ConfiguratorType,  CFEStiffOp< ConfiguratorType > > ( Grid, OpType ) {
00797     this->init();
00798   }
00799 
00800   /** Create the standard matrices of the non-interfaced standard tetrahedra
00801    */
00802   void createDefaultTetraMatrices() {
00803     const RealType scale = this->_grid.H();
00804 
00805     for ( short l = 0; l < 6; ++l ) {
00806       for ( short i = 0; i < 4; ++i ) {
00807         for ( short j = 0; j < 4; ++j ) {
00808           this->_defaultLocalTetraMatrix[l][i][j] =
00809             tpcfe::CFELookup<RealType>::_localStiffnessMatrixStdTetra[l][i][j] * scale;
00810         }
00811       }
00812     }
00813   }
00814 
00815   /** fill the matrix for a sub-tetrahedron
00816    */
00817   void preprocessLocalTetraMatrix ( const CFEElement< RealType > &/*e*/, const CFETetra<RealType> & t ) const {
00818     // Prepare stiffness matrix
00819     t.computeInverseTransformation();
00820     t.computeStiffnessMatrix();
00821 
00822     // Compute scaling factor
00823     const RealType det = t.determinant();
00824     const RealType scale = this->_grid.H() / ( 6 * det ); // note: t._lsm is the local stiffness matrix scaled by det^2, hence this factor needs to be divided out
00825 
00826     for ( short i = 0; i < 4; ++i ) {
00827       for ( short j = 0; j < 4; ++j ) {
00828         this->_localTetraMatrix[i][j] = t._lsm.get ( i, j ) * scale;
00829       }
00830     }
00831 
00832   }
00833 
00834   // end class CFEStiffOp
00835 };
00836 
00837 
00838 /** This class performs the computation of average (e.g. diffusivity or elasticity) coefficients for one (cubic) CFEElement.
00839  *  \author Schwen (based on code by Preusser)
00840  */
00841 template < typename RealType, typename NodalCoeffType >
00842 class CFEWeightProvider {
00843 protected:
00844   const qc::AArray< NodalCoeffType, qc::QC_3D > &_weights;
00845   aol::RandomAccessContainer<NodalCoeffType>     _localWeights;
00846   NodalCoeffType                                 _plusWeight, _minusWeight;
00847   aol::Vector < signed char >                     _sign;
00848 
00849 public:
00850   CFEWeightProvider ( const qc::AArray< NodalCoeffType, qc::QC_3D > &Weights, const CFEElement<RealType> &Element )
00851     : _weights ( Weights ), _localWeights ( 8 ), _plusWeight( aol::ZTrait<NodalCoeffType>::zero ), _minusWeight( aol::ZTrait<NodalCoeffType>::zero ), _sign ( 8 ) {
00852 
00853     updateForElement ( Element );
00854   }
00855 
00856   //! reuse this weightProvider for a different element (of the same grid)
00857   void updateForElement ( const CFEElement<RealType> &Element ) {
00858     const CFEType cfeType = Element.cfeType();
00859     int numMinusNodes = 0, numPlusNodes = 0;
00860     NodalCoeffType pWeight = aol::ZTrait<NodalCoeffType>::zero, mWeight = aol::ZTrait<NodalCoeffType>::zero;
00861 
00862     for ( int i = 0, bitMask = 1; i < 8; ++i, bitMask <<= 1 ) {
00863       const int index = Element.globalIndex ( i );
00864       const NodalCoeffType v = _weights [ index ];
00865       _sign[i] = ( cfeType._pureType & bitMask ) ? -1 : + 1;
00866       if ( _sign[i] == -1 ) {
00867         mWeight += v;
00868         ++numMinusNodes;
00869       } else {
00870         pWeight += v;
00871         ++numPlusNodes;
00872       }
00873       _localWeights[i] = v;
00874     }
00875     if ( numMinusNodes != 0 ) {
00876       _minusWeight = mWeight;
00877       _minusWeight /= static_cast<RealType>(numMinusNodes);
00878     }
00879     if ( numPlusNodes != 0 ) {
00880       _plusWeight = pWeight;
00881       _plusWeight /= static_cast<RealType>(numPlusNodes);
00882     }
00883   }
00884 
00885   // no implementation of destructor, copy constructor necessary
00886   // assignment operator should fail due to const reference member
00887 
00888   //! Return the weight value at the given global index
00889   inline NodalCoeffType nodalWeight ( const int globalIndex ) const {
00890     return ( _weights[ globalIndex ] );
00891   }
00892 
00893   //! Return the weight value at the given local index for the current element
00894   inline NodalCoeffType getLocalWeight ( const int localIndex ) const {
00895     return ( _localWeights[localIndex] );
00896   }
00897 
00898   //! Return the average weight at the given side of the interface
00899   inline NodalCoeffType meanWeight ( const signed char sign ) const {
00900     return ( sign == -1 ) ? _minusWeight : _plusWeight;
00901   }
00902 
00903   //! Return the average weight on the inside
00904   inline NodalCoeffType meanMinusWeight() const {
00905     return ( _minusWeight );
00906   }
00907 
00908   //! Return the average weight on the outside
00909   inline NodalCoeffType meanPlusWeight() const {
00910     return ( _plusWeight );
00911   }
00912 
00913   //! Return the sign of the level set function at a local node
00914   inline signed char getSign ( const int index ) const {
00915     return ( _sign[index] );
00916   }
00917 
00918 protected:
00919   //! standard constructor cannot work due to const reference member
00920   CFEWeightProvider();
00921 };
00922 
00923 
00924 /** Implements a standard operator for bilinear forms with non-constant coefficients
00925  *  \author Preusser
00926  */
00927 template < typename ConfiguratorType, typename ImpType >
00928 class CFEWeightedStandardOp : public CFEStandardOp < ConfiguratorType, ImpType > {
00929 protected:
00930   typedef typename ConfiguratorType::RealType   RealType;
00931 
00932   const qc::AArray<RealType, qc::QC_3D> &_nodalCoeffs;
00933   mutable CFEWeightProvider< RealType, RealType > _weightProvider; //!< store weight data for current element
00934 
00935 public:
00936   CFEWeightedStandardOp ( const qc::AArray<RealType, qc::QC_3D> &NodalCoeffs,
00937                           const typename ConfiguratorType::GridType & Grid,
00938                           const aol::OperatorType OpType = aol::ONTHEFLY ) : CFEStandardOp< ConfiguratorType, ImpType > ( Grid, OpType ), _nodalCoeffs ( NodalCoeffs ), _weightProvider ( NodalCoeffs, CFEElement<RealType>() ) {
00939     if ( qc::GridSize<qc::QC_3D> ( Grid ) != qc::GridSize<qc::QC_3D>( _nodalCoeffs ) )
00940       throw aol::Exception ( "CFEWeightedStandardOp::CFEWeightedStandardOp dimensions of grid and weights do not match", __FILE__, __LINE__ );
00941   }
00942 
00943   /** Clear all matrix entries in the local matrix and gather the local weights
00944    */
00945   template <typename ContribProcessorType>
00946   inline void preProcessLocalContribution ( const CFEElement< RealType> &e,
00947                                             ContribProcessorType &/*cpt*/ ) const {
00948     _weightProvider.updateForElement ( e );
00949   }
00950 
00951   /** Assemble the local matrix for a non-interfaced element
00952    */
00953   template <typename ContribProcessorType>
00954   inline void computeNonInterfacedHexaContrib ( const CFEElement< RealType> &e,
00955                                                 ContribProcessorType &cpt ) const {
00956     for ( short i = 0; i < 8; ++i ) {
00957       const int i1 = e.globalIndex ( i );
00958       const RealType w = _weightProvider.getLocalWeight ( i );
00959       for ( short j = 0; j < 8; ++j ) {
00960         const int j1 = e.globalIndex ( j );
00961         cpt.processContribution ( i1, j1, w*this->_defaultLocalHexaMatrix[i][j] );
00962       }
00963     }
00964   }
00965 
00966   /** Multiply the local tetra matrix by the weight of the parent node of the tetra
00967    */
00968 
00969   inline const RealType getTetraContrib ( const CFEElement< RealType >& /*e*/,
00970                                           const CFETetra<RealType> & t,
00971                                           const int /*parent*/,
00972                                           const int i, const int j ) const {
00973 
00974     const RealType weight = _weightProvider.meanWeight ( t.getSign() );
00975     return weight*this->_localTetraMatrix[i][j];
00976   }
00977 
00978   // end class CFEWeightedStandardOp
00979 };
00980 
00981 
00982 
00983 /** Provides a mass matrix which works with real cfe-basis functions (with jumping gradients)
00984  */
00985 template < typename _ConfiguratorType >
00986 class CFEMassOpWI : public CFEWeightedStandardOp < _ConfiguratorType, CFEMassOpWI < _ConfiguratorType > > {
00987 public:
00988   typedef CFEStandardOp < _ConfiguratorType, CFEMassOpWI < _ConfiguratorType > > BaseType;
00989   typedef typename BaseType::ConfiguratorType ConfiguratorType;
00990   typedef typename BaseType::VectorType       VectorType;
00991 
00992 protected:
00993   typedef typename ConfiguratorType::RealType RealType;
00994 
00995 public:
00996   CFEMassOpWI ( const qc::AArray<RealType, qc::QC_3D> &NodalCoeffs,
00997                 const typename ConfiguratorType::GridType &Grid,
00998                 aol::OperatorType OpType = aol::ONTHEFLY ) :
00999       CFEWeightedStandardOp < ConfiguratorType,  CFEMassOpWI < ConfiguratorType > > ( NodalCoeffs, Grid, OpType ) {
01000     this->init();
01001   }
01002 
01003 
01004   /** Create the standard matrices of the non-interfaced standard tetrahedra
01005    */
01006   void createDefaultTetraMatrices() {
01007     const RealType stdTetraVolFactor = aol::Cub ( this->_grid.H() );
01008 
01009     for ( int l = 0; l < 6; ++l )
01010       for ( int i = 0; i < 4; ++i )
01011         for ( int j = 0; j < 4; ++j )
01012           this->_defaultLocalTetraMatrix[l][i][j] = tpcfe::CFELookup<RealType>::_localMassMatrixRefTetra[i][j] * stdTetraVolFactor;
01013   }
01014 
01015   /** fill the matrix for a sub-tetrahedron
01016    */
01017   void preprocessLocalTetraMatrix ( const CFEElement< RealType >& /*e*/, const CFETetra<RealType> & t ) const {
01018     const RealType volumeFactor = t.getVolume() / tpcfe::CFELookup<RealType>::_stdTetraVolume * aol::Cub ( this->_grid.H() );
01019 
01020     for ( int i = 0; i < 4; ++i )
01021       for ( int j = 0; j < 4; ++j )
01022         this->_localTetraMatrix[i][j] = tpcfe::CFELookup<RealType>::_localMassMatrixRefTetra[i][j] * volumeFactor;
01023   }
01024 
01025   // end class CFEMassOpWI
01026 };
01027 
01028 
01029 
01030 /** Class for a stiffness-matrix with a spatially dependent coefficient
01031  *  Provides assymbly of matrices for bilinear forms like
01032  *  \f[  (a(x) \nabla\phi_i, \nabla\phi_j)_{0,2,E} \f]
01033  *  \author Preusser
01034  */
01035 template < typename _ConfiguratorType >
01036 class CFEStiffOpWI : public CFEWeightedStandardOp < _ConfiguratorType, CFEStiffOpWI < _ConfiguratorType > > {
01037 public:
01038   typedef CFEStandardOp < _ConfiguratorType, CFEStiffOpWI < _ConfiguratorType > > BaseType;
01039   typedef typename BaseType::ConfiguratorType ConfiguratorType;
01040   typedef typename BaseType::VectorType       VectorType;
01041 
01042 protected:
01043   typedef typename ConfiguratorType::RealType RealType;
01044 
01045 public:
01046   CFEStiffOpWI ( const qc::AArray<RealType, qc::QC_3D> &NodalCoeffs,
01047                  const typename ConfiguratorType::GridType & Grid, const aol::OperatorType OpType = aol::ONTHEFLY ) :
01048       CFEWeightedStandardOp < ConfiguratorType,  CFEStiffOpWI< ConfiguratorType > > ( NodalCoeffs, Grid, OpType ) {
01049     this->init();
01050   }
01051 
01052   /** Create the standard matrices of the non-interfaced standard tetrahedra
01053    */
01054   void createDefaultTetraMatrices() {
01055     const RealType scale = this->_grid.H();
01056 
01057     for ( short l = 0; l < 6; ++l )
01058       for ( short i = 0; i < 4; ++i )
01059         for ( short j = 0; j < 4; ++j )
01060           this->_defaultLocalTetraMatrix[l][i][j] = tpcfe::CFELookup<RealType>::_localStiffnessMatrixStdTetra[l][i][j] * scale;
01061   }
01062 
01063   /** fill the matrix for a sub-tetrahedron
01064    */
01065   void preprocessLocalTetraMatrix ( const CFEElement< RealType >& /*e*/, const CFETetra<RealType> & t ) const {
01066     // Prepare stiffness matrix
01067     t.computeInverseTransformation();
01068     t.computeStiffnessMatrix();
01069 
01070     // Compute scaling factor
01071     const RealType det  = t.determinant();
01072     const RealType scale = this->_grid.H() / ( 6 * det ) ;
01073 
01074     for ( short i = 0; i < 4; ++i )
01075       for ( short j = 0; j < 4; ++j )
01076         this->_localTetraMatrix[i][j] = t._lsm.get ( i, j ) * scale;
01077   }
01078 
01079   // end class CFEStiffOpWI
01080 };
01081 
01082 }
01083 
01084 #endif

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