fei_FEI_Impl.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.                    */
/*--------------------------------------------------------------------*/

#ifndef _fei_FEI_Impl_hpp_
#define _fei_FEI_Impl_hpp_

#include <fei_macros.hpp>
#include <fei_mpi.h>
#include <FEI.hpp>
#include <fei_Factory.hpp>

#include <set>
#include <vector>
#include <map>

namespace fei {

  class FEI_Impl : public FEI, private fei::Logger {
  public:
    FEI_Impl(fei::SharedPtr<LibraryWrapper> wrapper,
          MPI_Comm comm,
          int masterRank=0);

    FEI_Impl(const fei::Factory* factory,
          MPI_Comm comm,
          int masterRank=0);

    virtual ~FEI_Impl();

   fei::SharedPtr<fei::LinearSystem> getLinearSystem();

    int parameters(int numParams, 
                   const char *const* paramStrings);

   int setIDLists(int numMatrices,
                  const int* matrixIDs,
                  int numRHSs,
                  const int* rhsIDs);

   int setSolveType(int solveType);

   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);

   int initCRMult(int numCRNodes,
                  const GlobalID* CRNodeIDs,
                  const int *CRFieldIDs,
                  int& CRID); 

   int initCRPen(int numCRNodes,
                 const GlobalID* CRNodeIDs, 
                 const int *CRFieldIDs,
                 int& CRID); 

   int initComplete();

   int setCurrentMatrix(int matrixID);

   int setCurrentRHS(int rhsID);

   int resetSystem(double s=0.0);

   int resetMatrix(double s=0.0);

   int resetRHSVector(double s=0.0);

  int resetInitialGuess(double s=0.0);

    int loadNodeBCs(int numNodes,
                    const GlobalID *nodeIDs,
                    int fieldID,
                    const int* offsetsIntoField,
                    const double* prescribedValues);

    int loadElemBCs(int numElems,
                    const GlobalID* elemIDs,  
                    int fieldID,
                    const double *const *alpha,  
                    const double *const *beta,  
                    const double *const *gamma);

   int sumInElem(GlobalID elemBlockID,
                 GlobalID elemID,
                 const GlobalID* elemConn,
                 const double* const* elemStiffness,
                 const double* elemLoad,
                 int elemFormat);

   int sumInElemMatrix(GlobalID elemBlockID,
                       GlobalID elemID,
                       const GlobalID* elemConn,
                       const double* const* elemStiffness,
                       int elemFormat);

   int sumInElemRHS(GlobalID elemBlockID,
                    GlobalID elemID,
                    const GlobalID* elemConn,
                    const double* elemLoad);

   int loadCRMult(int CRMultID,
                  int numCRNodes,
                  const GlobalID* CRNodeIDs,
                  const int* CRFieldIDs,
                  const double* CRWeights,
                  double CRValue);

   int loadCRPen(int CRPenID,
                 int numCRNodes,
                 const GlobalID* CRNodeIDs,
                 const int* CRFieldIDs,
                 const double* CRWeights,
                 double CRValue,
                 double penValue);

   int putIntoRHS(int IDType,
                  int fieldID,
                  int numIDs,
                  const GlobalID* IDs,
                  const double* coefficients);

   int sumIntoRHS(int IDType,
                  int fieldID,
                  int numIDs,
                  const GlobalID* IDs,
                  const double* coefficients);

   int setMatScalars(int numScalars,
                     const int* IDs, 
                     const double* scalars);

   int setRHSScalars(int numScalars,
                     const int* IDs,
                     const double* scalars);

   int loadComplete(bool applyBCs=true,
                    bool globalAssemble=true);

   int residualNorm(int whichNorm,
                    int numFields,
                    int* fieldIDs,
                    double* norms);

   int solve(int& status);

   int iterations(int& itersTaken) const;

   int version(const char*& versionString);

   int cumulative_cpu_times(double& initTime,
                            double& loadTime,
                            double& solveTime,
                            double& solnReturnTime);

    int getBlockNodeSolution(GlobalID elemBlockID,  
                             int numNodes, 
                             const GlobalID *nodeIDs, 
                             int *offsets,
                             double *results);

    int getNodalSolution(int numNodes,
                         const GlobalID* nodeIDs,
                         int* offsets,
                         double* results);

    int getBlockFieldNodeSolution(GlobalID elemBlockID,
                                  int fieldID,
                                  int numNodes, 
                                  const GlobalID *nodeIDs, 
                                  double *results);
         
    int getBlockElemSolution(GlobalID elemBlockID,  
                             int numElems, 
                             const GlobalID *elemIDs,
                             int& numElemDOFPerElement,
                             double *results);

   int getNumCRMultipliers(int& numMultCRs);

   int getCRMultIDList(int numMultCRs, int* multIDs);

   int getCRMultipliers(int numCRs,
                        const int* CRIDs,
                        double *multipliers);

// associated "puts" paralleling the solution return services.
// 
// the int sizing parameters are passed for error-checking purposes, so
// that the interface implementation can tell if the passed estimate
// vectors make sense -before- an attempt is made to utilize them as
// initial guesses by unpacking them into the solver's native solution
// vector format (these parameters include lenNodeIDList, lenElemIDList,
// numElemDOF, and numMultCRs -- all other passed params are either 
// vectors or block/constraint-set IDs)

    int putBlockNodeSolution(GlobalID elemBlockID, 
                             int numNodes, 
                             const GlobalID *nodeIDs, 
                             const int *offsets,
                             const double *estimates);

    int putBlockFieldNodeSolution(GlobalID elemBlockID, 
                                  int fieldID, 
                                  int numNodes, 
                                  const GlobalID *nodeIDs, 
                                  const double *estimates);
         
    int putBlockElemSolution(GlobalID elemBlockID,  
                             int numElems, 
                             const GlobalID *elemIDs, 
                             int dofPerElem,
                             const double *estimates);

    int putCRMultipliers(int numMultCRs, 
                         const int* CRIDs,
                         const double* multEstimates);

// utility functions that aid in integrating the FEI calls..............

// support methods for the "gets" and "puts" of the soln services.


    int getBlockNodeIDList(GlobalID elemBlockID,
                           int numNodes,
                           GlobalID *nodeIDs);

   int getBlockElemIDList(GlobalID elemBlockID, 
                          int numElems, 
                          GlobalID* elemIDs);
 
// miscellaneous self-explanatory "read-only" query functions............ 
 
    int getNumSolnParams(GlobalID nodeID, int& numSolnParams) const;

    int getNumElemBlocks(int& numElemBlocks) const;

    int getNumBlockActNodes(GlobalID blockID, int& numNodes) const;

    int getNumBlockActEqns(GlobalID blockID, int& numEqns) const;

    int getNumNodesPerElement(GlobalID blockID, int& nodesPerElem) const;
    
    int getNumEqnsPerElement(GlobalID blockID, int& numEqns) const;

    int getNumBlockElements(GlobalID blockID, int& numElems) const;

    int getNumBlockElemDOF(GlobalID blockID, int& DOFPerElem) const;


    int getParameters(int& numParams, char**& paramStrings);

    //And now a couple of non-FEI query functions that Sandia applications
    //need to augment the matrix-access functions. I (Alan Williams) will
    //argue to have these included in the FEI 2.1 specification update.

    int getFieldSize(int fieldID, int& numScalars);

    int getEqnNumbers(GlobalID ID,
                      int idType, 
                      int fieldID,
                      int& numEqns,
                      int* eqnNumbers);

    int getNodalFieldSolution(int fieldID,
                              int numNodes,
                              const GlobalID* nodeIDs,
                              double* results);

    int getNumLocalNodes(int& numNodes);

    int getLocalNodeIDList(int& numNodes,
                           GlobalID* nodeIDs,
                           int lenNodeIDs);

    int putNodalFieldData(int fieldID,
                          int numNodes,
                          const GlobalID* nodeIDs,
                          const double* nodeData);

  private: //methods
    void basic_initializations();

    int inputRHS(int IDType,
                 int fieldID,
                 int numIDs,
                 const GlobalID* IDs,
                 const double* rhsEntries,
                 bool sumInto);

    int fillNodeset(int blockID) const;

    int aggregateSystem();

    int aggregateSystem_LinSysCore();

  private: //attributes

    std::vector<fei::SharedPtr<LibraryWrapper> > wrapper_;

    int nodeIDType_;
    int elemIDType_;
    int constraintIDType_;

    std::vector<fei::SharedPtr<fei::Factory> > factory_;
    bool createdFactory_;
    fei::SharedPtr<fei::VectorSpace> rowSpace_;
    fei::SharedPtr<fei::MatrixGraph> matGraph_;

    fei::SharedPtr<fei::Vector> x_;
    std::vector<fei::SharedPtr<fei::Vector> > b_;
    std::vector<fei::SharedPtr<fei::Matrix> > A_;
    fei::SharedPtr<fei::LinearSystem> linSys_;

    bool newData_;

    Data *soln_fei_matrix_;
    Data *soln_fei_vector_;

    MPI_Comm comm_;
    int masterRank_;
    int localProc_;
    int numProcs_;

    int numParams_;
    char** paramStrings_;

    std::vector<int> matrixIDs_;
    std::vector<int> rhsIDs_;

    std::vector<double> matScalars_;
    bool matScalarsSet_;
    std::vector<double> rhsScalars_;
    bool rhsScalarsSet_;

    int constraintID_;

    int index_soln_;
    int index_current_;
    int index_current_rhs_row_;

    int solveType_;
    int iterations_;

    bool setSolveTypeCalled_;
    bool initPhaseIsComplete_;

    bool aggregateSystemFormed_;
    int newMatrixDataLoaded_;

    int solveCounter_;

    double initTime_, loadTime_, solveTime_, solnReturnTime_;

    std::vector<int> iwork_;

    mutable std::set<int> nodeset_;
    mutable bool nodeset_filled_;
    mutable int nodeset_blockid_;

    mutable std::map<int,int> block_dof_per_elem_;
    bool any_blocks_have_elem_dof_;
  };//class FEI_Impl

}//namespace fei

#endif