SIMbase.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef _SIM_BASE_H
#define _SIM_BASE_H
 
#include "SIMadmin.h"
#include "SIMdependency.h"
#include "TimeDomain.h"
#include "Property.h"
#include "MatVec.h"
#include <set>
 
class IntegrandBase;
class NormBase;
class ForceBase;
class AnaSol;
class SAM;
class AlgEqSystem;
class LinSolParams;
class SystemMatrix;
class SystemVector;
class FunctionBase;
class RealFunc;
class VecFunc;
class TractionFunc;
class ScalarFunc;
class Vec4;
class Vec3;
 
 
struct Mode
{
  int    eigNo;   
  double eigVal;  
  double damping; 
  Vector eigVec;  
  Vector eqnVec;  
 
  Mode() : eigNo(0), eigVal(0.0), damping(0.0) {}
  bool orthonormalize(const SystemMatrix& mat);
  bool computeDamping(const SystemMatrix& mat);
};
 
 
class SIMbase : public SIMadmin, public SIMdependency
{
protected:
  explicit SIMbase(IntegrandBase* itg);
 
public:
  virtual ~SIMbase();
 
 
  // Model input and pre-processing methods
  // ======================================
 
  bool readModel(const char* fileName);
 
  virtual bool createFEMmodel(char resetNumb) = 0;
 
  virtual void clearProperties();
 
  virtual void initForSingleStep() {}
  virtual void initForMultiStep() {}
 
  virtual bool preprocessC(const IntVec& ignored, bool fixDup, double time0);
 
  virtual bool merge(SIMbase* that,
                     const std::map<int,int>* old2new = nullptr, int poff = 0);
 
  bool initSystem(LinAlg::MatrixType mType,
                  size_t nMats = 1, size_t nVec = 1, size_t nScl = 0,
                  bool withRF = false);
 
  bool initSystem(const SIMbase* that);
 
  void initLHSbuffers();
 
  bool setAssociatedRHS(size_t iMat, size_t iVec);
 
  bool setMode(int mode, bool needIntegr = true, bool resetSol = false);
 
  void setIntegrationPrm(unsigned short int i, double prm);
 
  void setQuadratureRule(size_t ng, bool redimBuffers = false,
                         bool printQP = false);
 
  virtual bool printProblem() const;
 
  const IntegrandBase* getProblem() const { return myProblem; }
 
  virtual ASM::InterfaceChecker* getInterfaceChecker(size_t) const
  { return nullptr; }
 
  void clearProblem() { myProblem = nullptr; }
 
  virtual std::string getName() const { return "SIMbase"; }
  virtual bool mixedProblem() const { return false; }
 
  const LinSolParams* getSolParams() const { return mySolParams; }
 
  virtual unsigned short int getNoParamDim() const = 0;
  virtual size_t getNoSpaceDim() const { return nsd; }
  size_t getNoFields(int basis = 0) const;
  size_t getNoDOFs(bool subSim = false) const;
  size_t getNoNodes(int basis = 0) const;
  size_t getNoElms(bool includeXelms = false, bool includeZelms = false) const;
  size_t getNoSolutions(bool allocated = false) const;
  int getNoPatches() const { return nGlPatches; }
  size_t getNoEquations() const;
  size_t getNoConstraints() const;
  virtual size_t getNoRHS() const;
  unsigned char getNoBasis() const;
 
  char getNodeType(int inod) const;
  Vec4 getNodeCoord(int inod) const;
  bool isFixed(int inod, int dof = 123) const;
  int getGlobalNode(int node) const;
  int getLocalNode(int node) const;
  bool getElmNodes(IntVec& mnpc, int iel) const;
  virtual std::vector<IntVec> getElmConnectivities() const = 0;
 
  void getBoundaryNodes(int pcode, IntVec& glbNodes,
                        std::vector<Vec3>* XYZ = nullptr) const;
 
  int findClosestNode(const Vec3&) const;
 
  IntVec getNodeSet(const std::string& setName) const;
 
  bool initDirichlet(double time = 0.0);
  bool hasTimeDependentDirichlet() const;
 
  virtual bool updateDirichlet(double time = 0.0,
                               const Vector* prevSol = nullptr,
                               bool tangent = false);
 
  virtual bool updateConfiguration(const Vector&) { return true; }
  virtual bool updateRotations(const RealArray&, double = 0.0) { return true; }
 
  bool updateGrid(const RealArray& displ);
  bool updateGrid(const std::string& field);
 
  void setRefined(int nref) { isRefined = nref; }
  int getRefined() const { return isRefined; }
 
  bool hasElementActivator(double t1 = -1.0, double t0 = 0.0) const;
 
  void updateForNewElements(Vector& solution, const TimeDomain& time,
                            int verbose = 0) const;
 
 
  // Computational methods
  // =====================
 
  virtual bool assembleSystem(const TimeDomain& time, const Vectors& prevSol,
                              bool newLHSmatrix = true, bool poorConvg = false);
 
  bool assembleSystem(double t0 = 0.0, const Vectors& pSol = Vectors())
  { return this->assembleSystem(TimeDomain(t0),pSol); }
 
  bool extractLoadVec(Vector& loadVec, size_t idx = 0,
                      const char* hd = nullptr) const;
  bool extractScalars(RealArray& values) const;
  double extractScalar(size_t idx = 0) const;
 
  bool applyDirichlet(Vector& glbVec) const;
 
  bool solveEqSystem(Vector& solution, size_t idxRHS, double* rCond,
                     int printSol = 0, bool dumpEqSys = false,
                     const char* compName = "displacement");
 
  virtual bool solveSystem(Vector& solution, int printSol, double* rCond,
                           const char* compName = "displacement",
                           size_t idxRHS = 0)
  {
    return this->solveEqSystem(solution,idxRHS,rCond,printSol,true,compName);
  }
 
  bool solveSystem(Vector& solution, int printSol = 0,
                   const char* compName = "displacement")
  {
    return this->solveSystem(solution,printSol,nullptr,compName);
  }
 
  bool solveSystem(Vectors& solution, int printSol = 0,
                   const char* cmpName = "displacement");
 
  void getWorstDofs(const Vector& x, const Vector& r,
                    size_t nWorst, double eps, int iteNorm,
                    std::map<std::pair<int,int>,RealArray>& worst) const;
 
  virtual void iterationNorms(const Vector& x, const Vector& r, double& eNorm,
                              double& rNorm, double& dNorm) const;
 
  double solutionNorms(const Vector& x, double* inf = nullptr,
                       size_t* ind = nullptr, size_t nf = 0,
                       char type = 'D') const;
 
  bool solutionNorms(const TimeDomain& time,
                     const Vectors& psol, const Vectors& ssol,
                     Vectors& gNorm, Matrix* eNorm = nullptr,
                     const char* name = nullptr);
  bool solutionNorms(const TimeDomain& time, const Vectors& psol,
                     Vectors& gNorm, Matrix* eNorm = nullptr)
  { return this->solutionNorms(time,psol,Vectors(),gNorm,eNorm); }
  bool solutionNorms(const Vector& psol, const Vectors& ssol,
                     Matrix& eNorm, Vectors& gNorm,
                     const char* name = nullptr)
  {
    return this->solutionNorms(TimeDomain(),Vectors(1,psol),ssol,
                               gNorm,&eNorm,name);
  }
  bool solutionNorms(const Vector& psol, Matrix& eNorm, Vectors& gNorm)
  {
    return this->solutionNorms(TimeDomain(),Vectors(1,psol),Vectors(),
                               gNorm,&eNorm);
  }
  bool solutionNorms(const Vector& psol, const Vectors& ssol, Vectors& gNorm,
                     const char* name = nullptr)
  {
    return this->solutionNorms(TimeDomain(),Vectors(1,psol),ssol,gNorm,
                               nullptr,name);
  }
 
  virtual void printStep(int istep, const TimeDomain& time) const;
 
  virtual void printSolutionSummary(const Vector& solution, int printSol = 0,
                                    const char* compName = nullptr,
                                    std::streamsize outPrec = 0);
 
  bool getCurrentReactions(RealArray& RF, const Vector& psol,
                           int pcode = 0) const;
  bool haveReactions(int pcode = 0) const;
  virtual const RealArray* getReactionForces() const;
 
  bool systemModes(std::vector<Mode>& solution,
                   int nev, int ncv, int iop, double shift,
                   size_t iA = 0, size_t iB = 1);
  bool systemModes(std::vector<Mode>& solution, size_t iA = 0, size_t iB = 1)
  {
    return this->systemModes(solution,opt.nev,opt.ncv,opt.eig,opt.shift,iA,iB);
  }
 
  virtual bool haveBoundaryReactions(bool = false) const { return false; }
 
  bool assembleForces(const Vector& solution, double t0,
                      RealArray* R, Vector* S = nullptr);
 
 
  // Post-processing methods
  // =======================
 
  virtual bool project(Matrix& ssol, const Vector& psol,
                       SIMoptions::ProjectionMethod method = SIMoptions::GLOBAL,
                       const TimeDomain& time = TimeDomain()) const;
  bool project(Vector& ssol, const Vector& psol,
               SIMoptions::ProjectionMethod method = SIMoptions::GLOBAL,
               size_t iComp = 0) const;
 
  bool projectAnaSol(Vector& ssol, SIMoptions::ProjectionMethod method) const;
 
  bool project(RealArray& values, const FunctionBase* f,
               int basis = 1, int iField = 0, int nFields = 1,
               SIMoptions::ProjectionMethod method = SIMoptions::GLOBAL,
               double time = 0.0) const;
 
  bool evalSecondarySolution(Matrix& field, int pindx) const;
 
  virtual bool fieldProjections() const;
 
  virtual bool haveAnaSol() const { return mySol ? true : false; }
  virtual bool haveDualSol() const { return dualField ? true : false; }
 
  NormBase* getNormIntegrand() const;
  ForceBase* getBoundaryForceIntegrand(const Vec3* X0 = nullptr) const;
  ForceBase* getNodalForceIntegrand() const;
 
  const SAM* getSAM() const { return mySam; }
 
protected:
  VecFunc* getVecFunc(size_t patch, Property::Type ptype) const;
 
  virtual bool addConstraint(int patch, int lndx, int ldim, int dirs, int code,
                             int& ngnod, char basis = 1, bool ovrD = false);
 
  virtual void preprocessA() {}
  virtual bool preprocessBeforeAsmInit(int&) { return true; }
  virtual bool preprocessB() { return true; }
  virtual void preprocessResultPoints() = 0;
 
  int renumberNodes(bool renumMNPC = false);
  virtual bool renumberNodes(const std::map<int,int>&) { return true; }
 
  virtual bool extractPatchSolution(IntegrandBase* problem,
                                    const Vectors& sol, size_t pindx) const;
 
public:
  using SIMdependency::registerDependency;
  void registerDependency(const std::string& name, SIMdependency* sim,
                          short int nvc = 1, unsigned char basis = 1);
 
  bool extractPatchSolution(const Vectors& sol, size_t pindx) const
  { return this->extractPatchSolution(myProblem,sol,pindx); }
 
  size_t extractPatchSolution(const RealArray& sol, RealArray& vec,
                              const ASMbase* pch,
                              unsigned char nndof = 0,
                              unsigned char basis = 0) const;
 
  bool injectPatchSolution(RealArray& sol, const RealArray& vec,
                           const ASMbase* pch,
                           unsigned char nndof = 0,
                           unsigned char basis = 0) const;
 
  bool extractElmRes(const Matrix& glbRes, Vector& elRes,
                     size_t row, int pindx, std::string& name,
                     const char* prefix = nullptr) const;
  bool extractElmRes(const Matrix& glbRes, Vector& elRes,
                     size_t row, int pindx) const;
  bool extractElmRes(const Vector& glbRes, Vector& elRes, int pindx) const;
 
  int getLocalPatchIndex(int patchNo) const;
 
  const PatchVec& getFEModel() const { return myModel; }
  ASMbase* getPatch(int idx, bool glbIndex = false) const;
 
  bool setPatchMaterial(size_t patch) const;
 
  const std::map<int,int>& getGlob2LocMap() const { return myGlb2Loc; }
 
  PropertyVec::const_iterator begin_prop() const { return myProps.begin(); }
  PropertyVec::const_iterator end_prop() const { return myProps.end(); }
 
  SystemMatrix* getRayleighDampingMatrix(size_t iM = 1, size_t iK = 0) const;
  SystemMatrix* getLHSmatrix(size_t idx = 0, bool copy = false) const;
  SystemVector* getRHSvector(size_t idx = 0, bool copy = false) const;
  void addToRHSvector(size_t idx, const SystemVector& vec, double scale = 1.0);
 
  RealFunc* getSclFunc(int code) const;
 
  void setMDflag(char flag) { mdFlag = flag; }
  bool isFirst() const { return mdFlag <= 1; }
 
  void dumpEqSys(bool initialBlankLine = false);
  void dumpSolVec(const Vector& x,bool isExpanded = true, bool expOnly = false);
 
  virtual int printNRforces(const IntVec& = {}) const { return 0; }
 
protected:
  char getMDflag() const { return mdFlag; }
  virtual void shiftGlobalNums(int, int) {}
 
  virtual bool initMaterial(size_t) { return true; }
  virtual bool initBodyLoad(size_t) { return true; }
  virtual bool initNeumann(size_t) { return true; }
 
  virtual bool assembleDiscreteTerms(const IntegrandBase*,
                                     const TimeDomain&) { return true; }
  virtual bool assembleDiscreteItems(const IntegrandBase* itg,
                                     const TimeDomain& time, const Vectors&)
  { return this->assembleDiscreteTerms(itg,time); }
 
  virtual double externalEnergy(const Vectors& psol, const TimeDomain&) const;
 
  virtual bool postProcessNorms(Vectors&, Matrix*) { return true; }
 
  void generateThreadGroups(const Property& p, bool silence = false);
 
  void changeNumThreads();
 
  bool addMADOF(unsigned char basis, unsigned char nndof, bool other = true);
 
  double* theExtEnerg() { return &extEnergy; }
  const double* getExtEnerg() const { return &extEnergy; }
 
private:
  const IntVec& getMADOF(unsigned char basis, unsigned char nndof) const;
 
public:
  static bool ignoreDirichlet; 
  static bool preserveNOrder;  
 
protected:
  typedef std::map<int,RealFunc*>     SclFuncMap;
  typedef std::map<int,VecFunc*>      VecFuncMap;
  typedef std::map<int,TractionFunc*> TracFuncMap;
 
  typedef std::multimap<int,IntegrandBase*> IntegrandMap;
 
  // Model attributes
  unsigned char  nsd;       
  PatchVec       myModel;   
  PropertyVec    myProps;   
  SclFuncMap     myScalars; 
  VecFuncMap     myVectors; 
  TracFuncMap    myTracs;   
  IntegrandBase* myProblem; 
  IntegrandMap   myInts;    
  AnaSol*        mySol;     
  FunctionBase*  dualField; 
 
  std::vector<FunctionBase*> extrFunc; 
 
  int  nGlbNodes; 
  int  isRefined; 
  bool lagMTOK;   
  bool fixZeros;  
 
  struct DumpData
  {
    std::string   fname;  
    LinAlg::StorageFormat format; 
    std::set<int> step;   
    bool          expand; 
    int           count;  
    double        eps;    
 
    DumpData() : format(LinAlg::FLAT), expand(false), count(0), eps(1.0e-6) {}
 
    bool doDump() { return !fname.empty() && step.find(++count) != step.end(); }
  };
 
  // Post-processing attributes
  std::vector<DumpData> lhsDump; 
  std::vector<DumpData> rhsDump; 
  std::vector<DumpData> solDump; 
 
  // Parallel computing attributes
  int               nGlPatches; 
  IntVec            myPatches;  
  IntVec            myLoc2Glb;  
  std::map<int,int> myGlb2Loc;  
  const std::map<int,int>* g2l; 
  std::map<int,int> myDegenElm; 
  std::set<int>     myDupNodes; 
 
  // Equation solver attributes
  AlgEqSystem*  myEqSys;     
  SAM*          mySam;       
  LinSolParams* mySolParams; 
  LinSolParams* myGl2Params; 
 
private:
  size_t nIntGP; 
  size_t nBouGP; 
  size_t nDofS;  
  char   mdFlag; 
 
  mutable std::map<int,IntVec> extraMADOFs;
 
  mutable double extEnergy;  
  mutable Vector prevForces; 
};
 
#endif