fei_LinSysCoreFilter.hpp

00001 #ifndef _LinSysCoreFilter_hpp_
00002 #define _LinSysCoreFilter_hpp_
00003 
00004 /*--------------------------------------------------------------------*/
00005 /*    Copyright 2005 Sandia Corporation.                              */
00006 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00007 /*    non-exclusive license for use of this work by or on behalf      */
00008 /*    of the U.S. Government.  Export of this program may require     */
00009 /*    a license from the United States Government.                    */
00010 /*--------------------------------------------------------------------*/
00011 
00012 #include "fei_macros.hpp"
00013 #include "fei_defs.h"
00014 #include "fei_fwd.hpp"
00015 #include "fei_Filter.hpp"
00016 #include <fei_CSRMat.hpp>
00017 #include <fei_CSVec.hpp>
00018 
00019 namespace fei {
00020   class DirichletBCManager;
00021 }
00022 
00033 class LinSysCoreFilter : public Filter {
00034 
00035  public:
00036    // Constructor.
00037    LinSysCoreFilter(FEI_Implementation* owner, MPI_Comm comm,
00038         SNL_FEI_Structure* probStruct,
00039         LinearSystemCore* lsc,
00040         int masterRank=0);
00041 
00042    //Destructor
00043    virtual ~LinSysCoreFilter();
00044 
00045 
00046    // set a value (usually zeros) throughout the linear system
00047    virtual int resetSystem(double s);
00048    virtual int resetMatrix(double s);
00049    virtual int resetRHSVector(double s);
00050    virtual int resetInitialGuess(double s);
00051 
00052    virtual int deleteMultCRs();
00053 
00054    virtual int loadNodeBCs(int numNodes,
00055                    const GlobalID *nodeIDs,
00056                    int fieldID,
00057                    const int* offsetsIntoField,
00058                    const double* prescribedValues);
00059 
00060    virtual int loadElemBCs(int numElems,
00061                    const GlobalID *elemIDs,
00062                    int fieldID,
00063                    const double *const *alpha,  
00064                    const double *const *beta,
00065                    const double *const *gamma);
00066 
00067    virtual int sumInElem(GlobalID elemBlockID,
00068                  GlobalID elemID,
00069                  const GlobalID* elemConn,
00070                  const double* const* elemStiffness,
00071                  const double* elemLoad,
00072                  int elemFormat);
00073 
00074    virtual int sumInElemMatrix(GlobalID elemBlockID,
00075                        GlobalID elemID,
00076                        const GlobalID* elemConn,
00077                        const double* const* elemStiffness,
00078                        int elemFormat);
00079 
00080    virtual int sumInElemRHS(GlobalID elemBlockID,
00081                     GlobalID elemID,
00082                     const GlobalID* elemConn,
00083                     const double* elemLoad);
00084 
00085     virtual int loadCRMult(int CRMultID, 
00086                    int numCRNodes,
00087                    const GlobalID* CRNodes, 
00088                    const int* CRFields,
00089                    const double* CRWeights,
00090                    double CRValue);
00091 
00092     virtual int loadCRPen(int CRPenID, 
00093                   int numCRNodes, 
00094                   const GlobalID* CRNodes,
00095                   const int *CRFields,
00096                   const double* CRWeights,
00097                   double CRValue,
00098                   double penValue);
00099 
00100    virtual int putIntoRHS(int IDType,
00101       int fieldID,
00102         int numIDs,
00103       const GlobalID* IDs,
00104       const double* rhsEntries);
00105 
00106    virtual int sumIntoRHS(int IDType,
00107       int fieldID,
00108         int numIDs,
00109       const GlobalID* IDs,
00110       const double* rhsEntries);
00111 
00112    virtual int loadComplete();
00113 
00114     // set parameters associated with solver choice, etc.
00115     virtual int parameters(int numParams, const char *const* paramStrings);
00116 
00117     //get residual norms
00118     virtual int residualNorm(int whichNorm, int numFields,
00119                      int* fieldIDs, double* norms, double& residTime);
00120 
00121     // start iterative solution
00122     virtual int solve(int& status, double& sTime);
00123 
00124     // query function iterations performed.
00125     virtual int iterations() const {return(iterations_);};
00126 
00127 // Solution return services.......................................
00128  
00129     // return all nodal solution params on a block-by-block basis 
00130     virtual int getBlockNodeSolution(GlobalID elemBlockID,  
00131                              int numNodes, 
00132                              const GlobalID *nodeIDs, 
00133                              int *offsets,
00134                              double *results);
00135  
00136     virtual int getNodalSolution(int numNodes, 
00137          const GlobalID *nodeIDs, 
00138          int *offsets,
00139          double *results);
00140 
00141     // return nodal solution for one field on a block-by-block basis 
00142     virtual int getBlockFieldNodeSolution(GlobalID elemBlockID,
00143                                   int fieldID,
00144                                   int numNodes, 
00145                                   const GlobalID *nodeIDs, 
00146                                   double *results);
00147          
00148     // return element solution params on a block-by-block basis 
00149     virtual int getBlockElemSolution(GlobalID elemBlockID,  
00150                              int numElems, 
00151                              const GlobalID *elemIDs,
00152                              int& numElemDOFPerElement,
00153                              double *results);
00154 
00155    virtual int getCRMultipliers(int numCRs, const int* CRIDs, double* multipliers);
00156 
00157 // associated "puts" paralleling the solution return services.
00158 // 
00159 // the int sizing parameters are passed for error-checking purposes, so
00160 // that the interface implementation can tell if the passed estimate
00161 // vectors make sense -before- an attempt is made to utilize them as
00162 // initial guesses by unpacking them into the solver's native solution
00163 // vector format (these parameters include lenNodeIDList, lenElemIDList,
00164 // numElemDOF, and numMultCRs -- all other passed params are either 
00165 // vectors or block/constraint-set IDs)
00166 
00167    // put nodal-based solution guess on a block-by-block basis 
00168    virtual int putBlockNodeSolution(GlobalID elemBlockID,
00169                             int numNodes,
00170                             const GlobalID *nodeIDs, 
00171                             const int *offsets,
00172                             const double *estimates);
00173 
00174     // put nodal-based guess for one field on a block-by-block basis 
00175     virtual int putBlockFieldNodeSolution(GlobalID elemBlockID, 
00176                                   int fieldID, 
00177                                   int numNodes, 
00178                                   const GlobalID *nodeIDs, 
00179                                   const double *estimates);
00180          
00181     // put element-based solution guess on a block-by-block basis
00182     virtual int putBlockElemSolution(GlobalID elemBlockID,  
00183                              int numElems, 
00184                              const GlobalID *elemIDs, 
00185                              int dofPerElem,
00186                              const double *estimates);
00187   
00188     virtual int putCRMultipliers(int numMultCRs, 
00189                          const int* CRIDs,
00190                          const double *multEstimates);
00191 
00192 //===== a couple of public non-FEI functions... ================================
00193 //These are intended to be used by an 'outer-layer' class like 
00194 //FEI_Implementation.
00195 //
00196   public:
00197     virtual int getNodalFieldSolution(int fieldID,
00198             int numNodes,
00199             const GlobalID* nodeIDs,
00200             double* results);
00201 
00202     virtual int putNodalFieldData(int fieldID,
00203         int numNodes,
00204         const GlobalID* nodeIDs,
00205         const double* nodeData);
00206 
00207     virtual int putNodalFieldSolution(int fieldID,
00208             int numNodes,
00209             const GlobalID* nodeIDs,
00210             const double* nodeData);
00211 
00212     virtual int unpackSolution();
00213 
00214     void setEqnCommMgr(EqnCommMgr* eqnCommMgr);
00215 
00216    EqnCommMgr* getEqnCommMgr() {return(eqnCommMgr_);};
00217 
00218    virtual int setNumRHSVectors(int numRHSs, int* rhsIDs);
00219    virtual int setCurrentRHS(int rhsID);
00220 
00221    virtual int exchangeRemoteEquations();
00222    virtual int exchangeRemoteBCs(std::vector<int>& essEqns,
00223                                  std::vector<double>& essAlpha,
00224                                  std::vector<double>& essGamma);
00225 
00226    virtual int implementAllBCs();
00227 
00228    virtual int enforceEssentialBCs(const int* eqns, const double* alpha,
00229                                   const double* gamma, int numEqns);
00230 
00231    virtual int enforceRemoteEssBCs(int numEqns, const int* eqns, 
00232          const int* const* colIndices, const int* colIndLens,
00233          const double* const* coefs);
00234 
00235    virtual int initialize();
00236 
00237 //==============================================================================
00238 //private functions for internal implementation of LinSysCoreFilter.
00239 //==============================================================================
00240   private:
00241    int initLinSysCore();
00242    void setLinSysCoreCREqns();
00243 
00244    int unpackRemoteContributions(EqnCommMgr& eqnCommMgr,
00245          int assemblyMode);
00246 
00247    int loadFEDataMultCR(int CRID,
00248       int numCRNodes,
00249       const GlobalID* CRNodes, 
00250       const int* CRFields,
00251       const double* CRWeights,
00252       double CRValue);
00253 
00254    int loadFEDataPenCR(int CRID,
00255            int numCRNodes,
00256            const GlobalID* CRNodes, 
00257            const int* CRFields,
00258            const double* CRWeights,
00259            double CRValue,
00260            double penValue);
00261 
00262    int storeNodalColumnCoefs(int eqn, const NodeDescriptor& node,
00263            int fieldID, int fieldSize,
00264            double* coefs);
00265 
00266    int storeNodalRowCoefs(const NodeDescriptor& node,
00267         int fieldID, int fieldSize,
00268         double* coefs, int eqn);
00269 
00270    int generalElemInput(GlobalID elemBlockID,
00271                         GlobalID elemID,
00272                         const double* const* elemStiffness,
00273                         const double* elemLoad,
00274                         int elemFormat);
00275 
00276    int generalElemInput(GlobalID elemBlockID,
00277                         GlobalID elemID,
00278                         const GlobalID* elemConn,
00279                         const double* const* elemStiffness,
00280                         const double* elemLoad,
00281                         int elemFormat);
00282 
00283    void storeNodalSendIndex(const NodeDescriptor& node, int fieldID, int col);
00284    void storeNodalSendEqn(const NodeDescriptor& node, int fieldID, int col,
00285         double* coefs);
00286    void storeNodalSendIndices(const NodeDescriptor& iNode, int iField,
00287             const NodeDescriptor& jNode, int jField);
00288 
00289    void storePenNodeSendData(const NodeDescriptor& iNode,
00290            int iField, int iFieldSize,
00291            double* iCoefs,
00292            const NodeDescriptor& jNode,
00293            int jField, int jFieldSize,
00294            double* jCoefs,
00295            double penValue, double CRValue);
00296 
00297    int storePenNodeData(const NodeDescriptor& iNode,
00298       int iField, int iFieldSize,
00299       double* iCoefs,
00300       const NodeDescriptor& jNode,
00301       int jField, int jFieldSize,
00302       double* jCoefs,
00303       double penValue, double CRValue);
00304 
00305    void allocElemStuff();
00306 
00307    int resolveConflictingCRs(EqnBuffer& bcEqns);
00308 
00309    int giveToMatrix_symm_noSlaves(int numPtRows,
00310           const int* ptRowNumbers,
00311           const double* const* coefs,
00312           int mode);
00313 
00314    int giveToBlkMatrix_symm_noSlaves(int numPtRows, const int* ptRows,
00315              int numBlkRows, const int* blkRowNumbers,
00316              const int* blkRowSizes,
00317              const double* const* coefs,
00318              int mode);
00319 
00320    int giveToMatrix(int numPtRows, const int* ptRows,
00321                     int numPtCols, const int* ptCols,
00322         const double* const* values,
00323         int mode);
00324 
00325    int giveToLocalReducedMatrix(int numPtRows, const int* ptRows,
00326         int numPtCols, const int* ptCols,
00327         const double* const* values,
00328         int mode);
00329 
00330    int getFromMatrix(int numPtRows, const int* ptRows,
00331                      const int* rowColOffsets, const int* ptCols,
00332                      int numColsPerRow, double** values);
00333 
00334    int sumIntoMatrix(fei::CSRMat& mat);
00335 
00336    int getEqnsFromMatrix(ProcEqns& procEqns, EqnBuffer& eqnData);
00337 
00338    int getEqnsFromRHS(ProcEqns& procEqns, EqnBuffer& eqnData);
00339 
00340    int giveToRHS(int num, const double* values,
00341      const int* indices, int mode);
00342 
00343    int giveToLocalReducedRHS(int num, const double* values,
00344            const int* indices, int mode);
00345 
00346    int getFromRHS(int num, double* values, const int* indices);
00347 
00348    int sumIntoRHS(fei::CSVec& vec);
00349 
00350    int getEqnSolnEntry(int eqnNumber, double& solnValue);
00351 
00352    int getSharedRemoteSolnEntry(int eqnNumber, double& solnValue);
00353 
00354    int getReducedSolnEntry(int eqnNumber, double& solnValue);
00355 
00356    int formResidual(double* residValues, int numLocalEqns);
00357 
00358    int getRemoteSharedEqns(int numPtRows, const int* ptRows,
00359          ProcEqns& remoteProcEqns);
00360 
00361    int resetTheMatrix(double s);
00362    int resetTheRHSVector(double s);
00363 
00364    int assembleEqns(int numPtRows, 
00365                     int numPtCols,
00366                     const int* rowNumbers,
00367                     const int* colIndices,
00368                     const double* const* coefs,
00369                     bool structurallySymmetric,
00370         int numBlkEqns, int* blkEqns, int* blkSizes,
00371         bool useBlkEqns, int mode);
00372 
00373    int assembleReducedEqns();
00374 
00375    int assembleRHS(int numValues, const int* indices, const double* coefs, int mode);
00376 
00377    int assembleReducedRHS();
00378  
00379    void debugOutput(const char* mesg);
00380 
00381    int createEqnCommMgr_put();
00382 
00383 //==============================================================================
00384 //private LinSysCoreFilter variables
00385 //==============================================================================
00386   private:
00387 
00388    int timesInitializeCalled_;
00389 
00390     LinearSystemCore* lsc_;
00391     bool useLookup_;
00392 
00393     int internalFei_;
00394 
00395     bool newMatrixData_;
00396     bool newVectorData_;
00397     bool newConstraintData_;
00398     bool newBCData_;
00399     bool connectivitiesInitialized_;
00400     bool firstRemEqnExchange_;
00401     bool needToCallMatrixLoadComplete_;
00402 
00403     bool resolveConflictRequested_;
00404 
00405     int localStartRow_, localEndRow_, numLocalEqns_, numGlobalEqns_;
00406     int reducedStartRow_, reducedEndRow_, numReducedRows_;
00407     int numLocallyOwnedNodes_, numGlobalNodes_, firstLocalNodeNumber_;
00408 
00409     bool blockMatrix_;
00410     bool tooLateToChooseBlock_;
00411 
00412     int numLocalEqnBlks_;
00413     int localReducedBlkOffset_, numLocalReducedEqnBlks_;
00414 
00415     int iterations_;
00416     int numRHSs_;
00417     int currentRHS_;
00418     std::vector<int> rhsIDs_;
00419 
00420     int outputLevel_;
00421 
00422     MPI_Comm comm_;
00423     int masterRank_;
00424 
00425     SNL_FEI_Structure* problemStructure_;
00426     bool matrixAllocated_;
00427 
00428     std::vector<int> rowIndices_;
00429     std::vector<int> rowColOffsets_, colIndices_;
00430 
00431     fei::FillableMat *Kid_, *Kdi_, *Kdd_;
00432     fei::CSRMat csrD, csrKid, csrKdi, csrKdd, tmpMat1_, tmpMat2_;
00433     fei::CSVec fd_, tmpVec1_;
00434     int reducedEqnCounter_, reducedRHSCounter_;
00435     std::vector<int> rSlave_, cSlave_;
00436 
00437     int nodeIDType_;
00438     fei::DirichletBCManager* bcManager_; //Boundary condition manager
00439 
00440     EqnCommMgr* eqnCommMgr_; //equation communication manager
00441     EqnCommMgr* eqnCommMgr_put_; //only created if users call the 'put'
00442                                  // functions
00443 
00444     int maxElemRows_;
00445     std::vector<int> scatterIndices_;
00446     std::vector<int> blkScatterIndices_;
00447     std::vector<int> iworkSpace_, iworkSpace2_;
00448     std::vector<double> dworkSpace_;
00449     std::vector<const double*> dworkSpace2_;
00450 
00451     double** eStiff_;
00452     double* eStiff1D_;
00453     double* eLoad_;
00454 };
00455 
00456 #endif
00457 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends
Generated on Wed Apr 13 10:08:24 2011 for FEI by  doxygen 1.6.3