FEI_Implementation.hpp

#ifndef _FEI_Implementation_hpp_
#define _FEI_Implementation_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_mpi.h>
#include <fei_SharedPtr.hpp>
#include <fei_defs.h>

#include <FEI.hpp>

class FEI_Implementation : public FEI {

 public:
   FEI_Implementation(fei::SharedPtr<LibraryWrapper> libWrapper,
                      MPI_Comm comm,
                      int masterRank=0);

   virtual ~FEI_Implementation();

//public FEI functions:


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

//Structural initialization functions.............................

   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* CRNodes,
                  const int *CRFields,
                  int& CRID); 

   int initCRPen(int numCRNodes,
                 const GlobalID* CRNodes, 
                 const int *CRFields,
                 int& CRID); 

   int initComplete();

// FEI data loading sequence..........................................

   int setCurrentMatrix(int matID);

   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 CRID,
                  int numCRNodes,
                  const GlobalID* CRNodes,
                  const int* CRFields,
                  const double* CRWeights,
                  double CRValue);

   int loadCRPen(int CRID,
                 int numCRNodes,
                 const GlobalID* CRNodes,
                 const int* CRFields,
                 const double* CRWeights,
                 double CRValue,
                 double penValue);

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

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

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

// Equation solution services.....................................

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

// Solution return services.......................................
 
    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);

  const SNL_FEI_Structure& getProblemStructure() const
  { return *problemStructure_; }

  SNL_FEI_Structure& getProblemStructure()
  { return *problemStructure_; }

  //============================================================================
  private: //functions

    FEI_Implementation(const FEI_Implementation& src);
    FEI_Implementation& operator=(const FEI_Implementation& src);

    void deleteIDs();
    void deleteRHSScalars();

    int allocateInternalFEIs();

    void debugOut(const char* msg);
    void debugOut(const char* msg, int whichFEI);

    void buildLinearSystem();
    int aggregateSystem();

    void messageAbort(const char* msg);
    void notAllocatedAbort(const char* name);
    void needParametersAbort(const char* name);
    void badParametersAbort(const char* name);

    void setDebugOutput(const char* path, const char* name);

  //============================================================================
  private: //member variables

    fei::SharedPtr<LibraryWrapper> wrapper_;
    fei::SharedPtr<LinearSystemCore> linSysCore_;
    std::vector<fei::SharedPtr<LinearSystemCore> > lscArray_;
    bool haveLinSysCore_;
    bool haveFEData_;
    SNL_FEI_Structure* problemStructure_;
    Filter** filter_;

    int numInternalFEIs_;
    bool internalFEIsAllocated_;

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

    bool IDsAllocated_;

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

    int index_soln_filter_;
    int index_current_filter_;
    int index_current_rhs_row_;

    int solveType_;

    bool setSolveTypeCalled_;
    bool initPhaseIsComplete_;

    bool aggregateSystemFormed_;
    int newMatrixDataLoaded_;

    Data *soln_fei_matrix_;
    Data *soln_fei_vector_;

    MPI_Comm comm_;

    int masterRank_;
    int localRank_;
    int numProcs_;

    int outputLevel_;

    int solveCounter_;
    int debugOutput_;
#ifdef FEI_HAVE_IOSFWD
    std::ostream* dbgOStreamPtr_;
#else
    ostream* dbgOStreamPtr_;
#endif
    bool dbgFileOpened_;
#ifdef FEI_HAVE_IOSFWD
    std::ofstream* dbgFStreamPtr_;
#else
    ofstream* dbgFStreamPtr_;
#endif

    double initTime_, loadTime_, solveTime_, solnReturnTime_;

    int numParams_;
    char** paramStrings_;
};

#endif