SNL_FEI_Structure.hpp

Go to the documentation of this file.

#ifndef _SNL_FEI_Structure_hpp_
#define _SNL_FEI_Structure_hpp_
 
/*--------------------------------------------------------------------*/
/*    Copyright 2005 Sandia Corporation.                              */
/*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
/*    non-exclusive license for use of this work by or on behalf      */
/*    of the U.S. Government.  Export of this program may require     */
/*    a license from the United States Government.                    */
/*--------------------------------------------------------------------*/
 
#include <fei_fwd.hpp>
#include <fei_defs.h>
#include <fei_constants.hpp>
#include "fei_TemplateUtils.hpp"
#include <snl_fei_PointBlockMap.hpp>
#include <fei_EqnBuffer.hpp>
#include <fei_FieldDofMap.hpp>
#include <fei_CSRMat.hpp>
#include <fei_CSVec.hpp>
 
#include <fei_NodeCommMgr.hpp>
#include <fei_NodeDatabase.hpp>
#include <fei_NodeDescriptor.hpp>
 
#include <fei_Lookup.hpp>
 
class SNL_FEI_Structure : public Lookup {
 public:
   SNL_FEI_Structure(MPI_Comm comm);
 
   virtual ~SNL_FEI_Structure();
 
   int parameters(int numParams, const char*const* paramStrings);
 
   int initFields(int numFields, const int* fieldSizes, const int* fieldIDs,
                  const int* fieldTypes = NULL);
 
   int initElemBlock(GlobalID elemBlockID,
         int numElements,
         int numNodesPerElement,
         const int* numFieldsPerNode,
         const int* const* nodalFieldIDs,
         int numElemDofFieldsPerElement,
         const int* elemDofFieldIDs,
         int interleaveStrategy);
 
   int initElem(GlobalID elemBlockID,
    GlobalID elemID,
    const GlobalID* elemConn);
 
   int initSlaveVariable(GlobalID slaveNodeID, 
       int slaveFieldID,
       int offsetIntoSlaveField,
       int numMasterNodes,
       const GlobalID* masterNodeIDs,
       const int* masterFieldIDs,
       const double* weights,
       double rhsValue);
 
   int deleteMultCRs();
 
   int initSharedNodes(int numSharedNodes,
           const GlobalID *sharedNodeIDs,  
           const int* numProcsPerNode, 
           const int *const *sharingProcIDs);
 
   // constraint relation initialization
   //- lagrange multiplier formulation
   int initCRMult(int numCRNodes,
      const GlobalID* CRNodes,
      const int *CRFields,
      int& CRID);
 
   // - penalty function formulation
   int initCRPen(int numCRNodes,
     const GlobalID* CRNodes, 
     const int *CRFields,
     int& CRID); 
 
   int initComplete(bool generateGraph = true);
 
   const std::vector<int>& getFieldIDs() const
     { return fieldIDs_; }
 
   const int* getFieldIDsPtr()
     {
       int len = fieldDatabase_->size();
       workarray_.resize(len*2);
       fei::copyToArrays<std::map<int,int> >(*fieldDatabase_, len,
             &workarray_[0],
             &workarray_[0]+len);
       return( &workarray_[0] );
     }
 
   const int* getFieldSizesPtr()
     {
       int len = fieldDatabase_->size();
       workarray_.resize(len*2);
       fei::copyToArrays<std::map<int,int> >(*fieldDatabase_, len,
             &workarray_[0],
             &workarray_[0]+len);
       return( &workarray_[0]+len );
     }
 
   int getNumFields() { return( fieldDatabase_->size() ); };
 
   int getFieldSize(int fieldID)
     {
       std::map<int,int>::const_iterator
   f_iter = fieldDatabase_->find(fieldID);
 
       return(f_iter != fieldDatabase_->end() ? (*f_iter).second : -1);
     }
 
   fei::FieldDofMap<int>& getFieldDofMap()
     { return fieldDofMap_; }
 
   bool isInLocalElement(int nodeNumber);
 
   const int* getNumFieldsPerNode(GlobalID blockID);
 
   const int* const* getFieldIDsTable(GlobalID blockID);
 
   int getEqnNumber(int nodeNumber, int fieldID);
 
   int getOwnerProcForEqn(int eqn);
 
 
   //now the element-block functions
 
   int getNumElemBlocks() {return(blockIDs_.size());};
   const GlobalID* getElemBlockIDs() {return(&blockIDs_[0]);};
 
   void getElemBlockInfo(GlobalID blockID,
                         int& interleaveStrategy, int& lumpingStrategy,
                         int& numElemDOF, int& numElements,
                         int& numNodesPerElem, int& numEqnsPerElem);
 
   int addBlock(GlobalID blockID);
 
   int getBlockDescriptor(GlobalID blockID, BlockDescriptor*& block);
 
   int getBlockDescriptor_index(int index, BlockDescriptor*& block);
 
   int getIndexOfBlock(GlobalID blockID) const;
   int getBlockID(unsigned index) const
   {
     if (index < blockIDs_.size()) return blockIDs_[index];
     return -1;
   }
 
   int allocateBlockConnectivity(GlobalID blockID);
   void destroyConnectivityTables();
 
   ConnectivityTable& getBlockConnectivity(GlobalID blockID);
 
   void getScatterIndices_ID(GlobalID blockID, GlobalID elemID,
                             int interleaveStrategy,
                             int* scatterIndices);
 
   void getScatterIndices_index(int blockIndex, int elemIndex,
                             int interleaveStrategy,
                             int* scatterIndices);
 
   int getBlkScatterIndices_index(int blockIndex,
          int elemIndex,
          int* scatterIndices);
 
   void getScatterIndices_ID(GlobalID blockID, GlobalID elemID,
                             int interleaveStrategy,
                             int* scatterIndices,
           int* blkScatterIndices,
           int* blkSizes);
 
   void getScatterIndices_index(int blockIndex, int elemIndex,
                             int interleaveStrategy,
                             int* scatterIndices,
           int* blkScatterIndices,
        int* blkSizes);
 
 
   //now the shared-node lookup functions from the Lookup interface.
 
   int getNumSharedNodes() {return(nodeCommMgr_->getNumSharedNodes());};
 
   const int* getSharedNodeNumbers() {
      return(&(nodeCommMgr_->getSharedNodeNumbers())[0]);
   };
 
   const int* getSharedNodeProcs(int nodeNumber) {
      int index = nodeCommMgr_->getSharedNodeIndex_num(nodeNumber);
      if (index < 0) return(NULL);
      return(&(nodeCommMgr_->getSharedNodeProcs(index))[0]);
   };
 
   int getNumSharingProcs(int nodeNumber) {
      int index = nodeCommMgr_->getSharedNodeIndex_num(nodeNumber);
      if (index < 0) return(-1);
      return(nodeCommMgr_->getSharedNodeProcs(index).size());
   };
 
   int getNumSubdomains(int nodeNumber) {
     const NodeDescriptor* node = NULL;
     int err = nodeDatabase_->getNodeWithNumber(nodeNumber, node);
     if (err != 0) return(-1);
     GlobalID nodeID = node->getGlobalNodeID();
     return(nodeCommMgr_->getSharedNodeNumSubdomains(nodeID));
   }
 
   int* getSubdomainList(int nodeNumber) {
     const NodeDescriptor* node = NULL;
     int err = nodeDatabase_->getNodeWithNumber(nodeNumber, node);
     if (err != 0) return(NULL);
     GlobalID nodeID = node->getGlobalNodeID();
     return(&(*(nodeCommMgr_->getSharedNodeSubdomainList(nodeID)))[0]);
 
   }
 
   int translateToReducedNodeNumber(int nodeNumber, int proc);
 
 
   int getAssociatedNodeNumber(int eqnNumber)
     {
       int eqn = translateFromReducedEqn(eqnNumber);
       int nodeNumber = nodeDatabase_->getAssociatedNodeNumber(eqn);
       int reducedNodeNumber = -1;
       if (nodeNumber >= 0) {
   reducedNodeNumber = translateToReducedNodeNumber(nodeNumber, localProc_);
       }
       return( reducedNodeNumber );
     }
 
   int getAssociatedFieldID(int eqnNumber)
     {
       int eqn = translateFromReducedEqn(eqnNumber);
       return( nodeDatabase_->getAssociatedFieldID(eqn) );
     }
 
 
   //now the point-eqn to block-eqn queries...
 
   bool isExactlyBlkEqn(int ptEqn) {
      return(blkEqnMapper_->isExactlyBlkEqn(ptEqn));
   };
 
   int ptEqnToBlkEqn(int ptEqn) {
      return(blkEqnMapper_->eqnToBlkEqn(ptEqn));
   };
 
   int getOffsetIntoBlkEqn(int blkEqn, int ptEqn) {
      return(blkEqnMapper_->getBlkEqnOffset(blkEqn, ptEqn));
   };
 
   int getBlkEqnSize(int blkEqn) {
      return(blkEqnMapper_->getBlkEqnSize(blkEqn));
   }
 
 
   int getNumActiveNodes() {return(nodeDatabase_->getNodeIDs().size());}
 
   NodeDatabase& getNodeDatabase() { return( *nodeDatabase_ ); }
 
   std::map<GlobalID,int>& getActiveNodeIDList()
     { return( nodeDatabase_->getNodeIDs() ); }
 
   std::vector<int>& getGlobalNodeOffsets() {return(globalNodeOffsets_);}
   std::vector<int>& getGlobalEqnOffsets() {return(globalEqnOffsets_);}
   std::vector<int>& getGlobalBlkEqnOffsets() {return(globalBlkEqnOffsets_);}
 
   NodeCommMgr& getNodeCommMgr() {return(*nodeCommMgr_);}
   EqnCommMgr&  getEqnCommMgr()  {return(*eqnCommMgr_ );}
 
   void initializeEqnCommMgr();
 
   void getEqnInfo(int& numGlobalEqns, int& numLocalEqns,
                   int& localStartRow, int& localEndRow);
 
   int getEqnNumbers(GlobalID ID, int idType, int fieldID,
         int& numEqns, int* eqnNumbers);
 
   int getEqnNumbers(int numIDs, const GlobalID* IDs,
         int idType, int fieldID,
         int& numEqns, int* eqnNumbers);
 
   void getEqnBlkInfo(int& numGlobalEqnBlks, int& numLocalEqnBlks,
                      int& localBlkOffset);
 
   snl_fei::PointBlockMap& getBlkEqnMapper() {return(*blkEqnMapper_);}
 
   void destroyMatIndices();
 
   int getNumMultConstRecords() {return(multCRs_.size());};
 
   std::map<GlobalID,snl_fei::Constraint<GlobalID>*>&
     getMultConstRecords()
     {return(multCRs_);};
 
   int getMultConstRecord(int CRID, snl_fei::Constraint<GlobalID>*& multCR)
   {
     std::map<int,snl_fei::Constraint<GlobalID>*>::iterator
       cr_iter = multCRs_.find(CRID);
     int returncode = -1;
     if (cr_iter != multCRs_.end()) {
       multCR = (*cr_iter).second;
       returncode = 0;
     }
 
     return( returncode );
   }
 
   int getNumPenConstRecords() {return(penCRs_.size());}
   std::map<GlobalID,snl_fei::Constraint<GlobalID>*>&
     getPenConstRecords()
     { return(penCRs_); }
 
   int getPenConstRecord(int CRID, snl_fei::Constraint<GlobalID>*& penCR)
   {
     std::map<int,snl_fei::Constraint<GlobalID>*>::iterator
       cr_iter = penCRs_.find(CRID);
     int returncode = -1;
     if (cr_iter != penCRs_.end()) {
       penCR = (*cr_iter).second;
       returncode = 0;
     }
 
     return( returncode );
   }
 
   void addSlaveVariable(SlaveVariable* svar) {slaveVars_->push_back(svar);}
 
   int calculateSlaveEqns(MPI_Comm comm);
 
   fei::FillableMat* getSlaveDependencies() {return(slaveMatrix_);}
 
   EqnBuffer* getSlaveEqns() { return(slaveEqns_); }
 
   int numSlaveEquations() { return(numSlvs_); }
 
   bool isSlaveEqn(int eqn);
 
   bool translateToReducedEqn(int eqn, int& reducedEqn);
 
   int translateToReducedEqns(EqnCommMgr& eqnCommMgr);
 
   int translateToReducedEqns(EqnBuffer& eqnBuf);
 
   int translateToReducedEqns(ProcEqns& procEqns);
 
   int translateMatToReducedEqns(fei::CSRMat& mat);
 
   int translateFromReducedEqn(int reducedEqn);
 
   int getMasterEqnNumbers(int slaveEqn, std::vector<int>*& masterEqns);
 
   int getMasterEqnCoefs(int slaveEqn, std::vector<double>*& masterCoefs);
 
   int getMasterEqnRHS(int slaveEqn, double& rhsValue);
 
   int getNumGlobalEqns() { return( numGlobalEqns_ ); }
   int getNumLocalEqns() { return( numLocalEqns_ ); }
   int getFirstLocalEqn() { return( localStartRow_ ); }
   int getLastLocalEqn() { return( localEndRow_ ); }
 
   int getFirstReducedEqn() { return( reducedStartRow_ ); }
   int getLastReducedEqn()  { return( reducedEndRow_   ); }
 
   int getNumGlobalEqnBlks() { return( numGlobalEqnBlks_ ); }
   int getNumLocalEqnBlks() { return( numLocalEqnBlks_ ); }
   int getNumLocalReducedEqnBlks() { return( numLocalReducedEqnBlks_ ); }
   int getGlobalMaxBlkSize() { return(globalMaxBlkSize_); }
 
   int getNumLocalReducedEqns() { return( numLocalReducedRows_ ); }
 
   int getMatrixRowLengths(std::vector<int>& rowLengths);
   int getMatrixStructure(int** colIndices, std::vector<int>& rowLengths);
 
   int getMatrixStructure(int** ptColIndices, std::vector<int>& ptRowLengths,
        int** blkColIndices, int* blkIndices_1D,
        std::vector<int>& blkRowLengths,
        std::vector<int>& numPtRowsPerBlkRow);
 
   static int gatherSlaveEqns(MPI_Comm comm,
            EqnCommMgr* eqnCommMgr,
            EqnBuffer* slaveEqns);
 
   static int removeCouplings(EqnBuffer& eqnbuf, int& levelsOfCoupling);
 
   int calcTotalNumElemDOF();
   int calcNumMultCREqns();
 
   MPI_Comm getCommunicator() const { return( comm_ ); }
 
#ifdef FEI_HAVE_IOSFWD
   int setDbgOut(std::ostream& ostr, const char* path, const char* feiName);
#else
   int setDbgOut(ostream& ostr, const char* path, const char* feiName);
#endif
 
 private:
 
   NodeDescriptor* findNode(GlobalID nodeID);
   NodeDescriptor& findNodeDescriptor(GlobalID nodeID);
 
   int writeEqn2NodeMap();
 
   int getElemNodeDescriptors(int blockIndex, int elemIndex,
                              NodeDescriptor** nodes);
 
   int getNodeIndices_simple(NodeDescriptor** nodes, int numNodes,
           int fieldID,
           int* scatterIndices, int& offset);
 
   int getNodeIndices_simple(NodeDescriptor** nodes, int numNodes,
           int fieldID,
           int* scatterIndices, int& offset,
           int* blkScatterIndices,
           int* blkSizes, int& blkOffset);
 
   int getNodeMajorIndices(NodeDescriptor** nodes, int numNodes,
                           int** fieldIDs, int* fieldsPerNode,
                           int* scatterIndices, int& offset);
 
   int getNodeBlkIndices(NodeDescriptor** nodes, int numNodes,
       int* scatterIndices, int& offset);
 
   int getNodeMajorIndices(NodeDescriptor** nodes, int numNodes,
                           int** fieldIDs, int* fieldsPerNode,
                           int* scatterIndices, int& offset,
         int* blkScatterIndices,
         int* blkSizes, int& blkOffset);
 
   int getNodeMajorIndices(NodeDescriptor** nodes, int numNodes,
                           std::vector<int>* fieldIDs,
                           std::vector<int>& fieldsPerNode,
                           std::vector<int>& scatterIndices);
 
   int getFieldMajorIndices(NodeDescriptor** nodes, int numNodes,
                            int** fieldIDs, int* fieldsPerNode,
                            int* scatterIndices, int& offset);
 
   int getFieldMajorIndices(NodeDescriptor** nodes, int numNodes,
                                          std::vector<int>* fieldIDs,
                                          std::vector<int>& fieldsPerNode,
                                          std::vector<int>& scatterIndices);
 
   void calcGlobalEqnInfo(int numLocallyOwnedNodes, 
        int numLocalEqns,
        int numLocalEqnBlks);
 
   int finalizeActiveNodes();
   int finalizeNodeCommMgr();
   bool activeNodesInitialized();
 
   int formMatrixStructure();
 
   int initElemBlockStructure();
   int initMultCRStructure();
   int initPenCRStructure();
   int createMatrixPosition(int row, int col, const char* callingFunction);
   int createMatrixPositions(int row, int numCols, int* cols,
           const char* callingFunction);
   int createMatrixPositions(fei::CSRMat& mat);
 
   int createSymmEqnStructure(std::vector<int>& scatterIndices);
   int createBlkSymmEqnStructure(std::vector<int>& scatterIndices);
 
   int storeElementScatterIndices(std::vector<int>& scatterIndices);
   int storeElementScatterIndices_noSlaves(std::vector<int>& scatterIndices);
   int storeElementScatterBlkIndices_noSlaves(std::vector<int>& scatterIndices);
 
   void storeLocalNodeIndices(NodeDescriptor& iNode, int iField,
                                     NodeDescriptor& jNode, int jField);
   void storeNodalColumnIndices(int eqn, NodeDescriptor& node,
                 int fieldID);
   void storeNodalRowIndices(NodeDescriptor& node, int fieldID, int eqn);
   void storeNodalSendIndex(NodeDescriptor& node, int fieldID, int col);
   void storeNodalSendIndices(NodeDescriptor& iNode, int iField,
                                  NodeDescriptor& jNode, int jField);
 
   int assembleReducedStructure();
 
   bool nodalEqnsAllSlaves(const NodeDescriptor* node, std::vector<int>& slaveEqns);
 
   int initializeBlkEqnMapper();
 
   int setNumNodesAndEqnsPerBlock();
 
   void destroyBlockRoster();
 
#ifdef FEI_HAVE_IOSFWD
   std::ostream& dbgOut() { return( *dbgOStreamPtr_ ); }
#else
   ostream& dbgOut() { return( *dbgOStreamPtr_ ); }
#endif
 
   void addCR(int CRID,
       snl_fei::Constraint<GlobalID>*& cr,
       std::map<GlobalID,snl_fei::Constraint<GlobalID>* >& crDB);
 
   int setNodalEqnInfo();
   void setElemDOFEqnInfo();
   int setMultCREqnInfo();
 
   SNL_FEI_Structure(const SNL_FEI_Structure& src);
   SNL_FEI_Structure& operator=(const SNL_FEI_Structure& src);
 
   MPI_Comm comm_;
 
   int localProc_, masterProc_, numProcs_;
 
   std::vector<int> fieldIDs_;
   std::vector<int> fieldSizes_;
   std::map<int,int>* fieldDatabase_;
   fei::FieldDofMap<int> fieldDofMap_;
   std::vector<int> workarray_;
 
   std::vector<GlobalID> blockIDs_;
   std::vector<BlockDescriptor*> blocks_;
   std::vector<ConnectivityTable*> connTables_;
 
   NodeDatabase* nodeDatabase_;
 
   bool activeNodesInitialized_;
 
   std::vector<int> globalNodeOffsets_;
   std::vector<int> globalEqnOffsets_;
   std::vector<int> globalBlkEqnOffsets_;
 
   std::vector<SlaveVariable*>* slaveVars_;
   EqnBuffer* slaveEqns_;
   std::vector<int>* slvEqnNumbers_;
   int numSlvs_, lowestSlv_, highestSlv_;
   fei::FillableMat* slaveMatrix_;
   std::vector<int> globalNumNodesVanished_;
   std::vector<int> localVanishedNodeNumbers_;
 
   NodeCommMgr* nodeCommMgr_;
   EqnCommMgr* eqnCommMgr_;
   EqnCommMgr* slvCommMgr_;
 
   int numGlobalEqns_;
   int numLocalEqns_;
   int localStartRow_;
   int localEndRow_;
 
   int numLocalNodalEqns_;
   int numLocalElemDOF_;
   int numLocalMultCRs_;
 
   int reducedStartRow_, reducedEndRow_, numLocalReducedRows_;
   fei::FillableMat *Kid_, *Kdi_, *Kdd_;
   fei::CSRMat csrD, csrKid, csrKdi, csrKdd, tmpMat1_, tmpMat2_;
   int reducedEqnCounter_, reducedRHSCounter_;
   std::vector<int> rSlave_, cSlave_;
   std::vector<NodeDescriptor*> work_nodePtrs_;
 
   bool structureFinalized_;
   bool generateGraph_;
 
   fei::ctg_set<int>* sysMatIndices_;
 
   bool blockMatrix_;
   int numGlobalEqnBlks_;
   int numLocalEqnBlks_;
   int numLocalReducedEqnBlks_;
   int localBlkOffset_;
   int localReducedBlkOffset_;
   int globalMaxBlkSize_;
 
   int firstLocalNodeNumber_;
   int numGlobalNodes_;
 
   fei::ctg_set<int>* sysBlkMatIndices_;
   bool matIndicesDestroyed_;
 
   std::vector<int> workSpace_;
 
   snl_fei::PointBlockMap* blkEqnMapper_;
 
   std::map<GlobalID, snl_fei::Constraint<GlobalID>* > multCRs_;
 
   std::map<GlobalID, snl_fei::Constraint<GlobalID>* > penCRs_;
 
   bool checkSharedNodes_;
 
   std::string name_;
 
   int outputLevel_;
   bool debugOutput_;
   std::string dbgPath_;
#ifdef FEI_HAVE_IOSFWD
   std::ostream* dbgOStreamPtr_;
#else
   ostream* dbgOStreamPtr_;
#endif
   bool setDbgOutCalled_;
};
 
#endif