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