ASMs3D.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef _ASM_S3D_H
#define _ASM_S3D_H
 
#include "ASMstruct.h"
#include "ASM3D.h"
#include "BasisFunctionCache.h"
#include "Interface.h"
#include "ThreadGroups.h"
 
namespace utl {
  class Point;
}
 
namespace Go {
  class SplineSurface;
  class SplineVolume;
}
 
 
class ASMs3D : public ASMstruct, public ASM3D
{
  struct IJK
  {
    int I = 0; 
    int J = 0; 
    int K = 0; 
  };
 
  using IndexVec = std::vector<IJK> ; 
 
  struct Edge
  {
    int icnod = 0; 
    int incr  = 0; 
    int next();    
  };
 
  struct Face
  {
    int isnod = 0; 
    int incrI = 0; 
    int incrJ = 0; 
    int nnodI = 0; 
    int indxI = 1; 
    int next();    
  };
 
protected:
  class BasisFunctionCache : public ::BasisFunctionCache<3>
  {
  public:
    BasisFunctionCache(const ASMs3D& pch);
 
    BasisFunctionCache(const BasisFunctionCache& cache, int b);
 
    const std::array<size_t,3>& noElms() const { return nel; }
 
  protected:
    bool internalInit() override;
 
    BasisFunctionVals calculatePt(size_t el, size_t gp,
                                  bool reduced) const override;
 
    void calculateAll() override;
 
    size_t index(size_t el, size_t gp, bool reduced) const override;
 
    virtual void setupParameters();
 
    bool setupQuadrature();
 
    const ASMs3D& patch; 
 
    std::array<size_t,3> nel{}; 
 
  private:
    std::array<size_t,3> elmIndex(size_t el) const;
  };
 
public:
  struct BlockNodes
  {
    int  ibnod[8];  
    Edge edges[12]; 
    Face faces[6];  
    int  iinod;     
    int  inc[3];    
    int  nnodI;     
    int  nnodJ;     
    int  indxI;     
    int  indxJ;     
 
    BlockNodes()
    {
      memset(ibnod,0,8*sizeof(int));
      iinod = inc[0] = inc[1] = inc[2] = nnodI = nnodJ = 0;
      indxI = indxJ = 1;
    }
    int next();
  };
 
private:
  struct DirichletFace
  {
    Go::SplineSurface* surf;  
    int                dof;   
    int                code;  
    std::vector<Ipair> nodes; 
 
    DirichletFace(Go::SplineSurface* ss = nullptr, int d = 0, int c = 0)
    : surf(ss), dof(d), code(c) {}
  };
 
public:
  class InterfaceChecker : public ASM::InterfaceChecker
  {
  protected:
    const ASMs3D& myPatch; 
  public:
    explicit InterfaceChecker(const ASMs3D& pch) : myPatch(pch) {}
    virtual short int hasContribution(int, int I, int J, int K) const;
  };
 
  explicit ASMs3D(unsigned char n_f = 3);
  ASMs3D(const ASMs3D& patch, unsigned char n_f);
  ASMs3D(const ASMs3D& patch);
 
  virtual Go::SplineSurface* getBoundary(int dir, int = 1);
  virtual Go::SplineVolume* getBasis(int basis = 1);
  virtual const Go::SplineVolume* getBasis(int basis = 1) const;
  virtual void copyParameterDomain(const ASMbase* other);
 
 
  // Methods for model generation
  // ============================
 
  virtual bool read(std::istream&);
  virtual bool write(std::ostream&, int) const;
 
  virtual bool generateFEMTopology();
 
  virtual void clear(bool retainGeometry = false);
 
  virtual bool addXElms(short int dim, short int item,
                        size_t nXn, IntVec& nodes);
 
  virtual size_t getNodeIndex(int globalNum, bool noAddedNodes = false) const;
  virtual int getNodeID(size_t inod, bool noAddedNodes = false) const;
 
  virtual bool getElementCoordinates(Matrix& X, int iel,
                                     bool forceItg = false) const;
 
  virtual void getNodalCoordinates(Matrix& X, bool geo = false) const;
 
  virtual Vec3 getCoord(size_t inod) const;
 
  virtual Vec3 getElementCenter(int iel) const;
 
  virtual bool updateCoords(const Vector& displ);
 
  virtual void getBoundaryNodes(int lIndex, IntVec& nodes,
                                int basis, int thick = 1,
                                int = 0, bool local = false) const;
  virtual void getBoundary1Nodes(int lEdge, IntVec& nodes, int basis, int = 0,
                                 bool local = false, bool open = false) const;
 
  virtual int getCorner(int I, int J, int K, int basis) const;
 
  bool assignNodeNumbers(BlockNodes& nodes, int basis = 0);
 
  virtual bool checkRightHandSystem();
 
  virtual bool refine(int dir, const RealArray& xi);
  virtual bool uniformRefine(int dir, int nInsert);
  virtual bool raiseOrder(int ru, int rv, int rw, bool setOrder = false);
 
  virtual bool createProjectionBasis(bool init);
 
  virtual bool separateProjectionBasis() const;
 
 
  // Various methods for preprocessing of boundary conditions and patch topology
  // ===========================================================================
 
  virtual void constrainFace(int dir, bool open, int dof,
                             int code = 0, char basis = 1);
  virtual size_t constrainFaceLocal(int dir, bool open, int dof, int code = 0,
                                    bool project = false, char T1 = '\0');
 
  virtual void constrainEdge(int lEdge, bool open, int dof,
                             int code = 0, char basis = 1);
 
  virtual void constrainLine(int fdir, int ldir, double xi, int dof,
                             int code = 0, char basis = 1);
 
  virtual void constrainCorner(int I, int J, int K, int dof,
                               int code = 0, char basis = 1);
  virtual void constrainNode(double xi, double eta, double zeta, int dof,
                             int code = 0);
 
  virtual bool connectPatch(int face, ASM3D& neighbor, int nface, int norient,
                            int = 0, bool coordCheck = true, int thick = 1);
 
  virtual void closeBoundaries(int dir, int basis, int master);
 
  virtual bool collapseFace(int face, int edge = 0, int basis = 1);
 
  virtual bool addRigidCpl(int lindx, int ldim, int basis,
                           int& gMaster, const Vec3& Xmaster, bool extraPt);
 
  virtual void setNodeNumbers(const IntVec& nodes);
 
  virtual bool updateDirichlet(const std::map<int,RealFunc*>& func,
                               const std::map<int,VecFunc*>& vfunc, double time,
                               const std::map<int,int>* g2l, bool tangent);
 
  virtual bool selfInterconnect(const std::vector<Ipair>& nodes, double xtol);
 
 
  // Methods for integration of finite element quantities.
  // These are the main computational methods of the ASM class hierarchy.
  // ====================================================================
 
  virtual bool integrate(Integrand& integrand,
                         GlobalIntegral& glbInt, const TimeDomain& time);
 
  virtual bool integrate(Integrand& integrand, int lIndex,
                         GlobalIntegral& glbInt, const TimeDomain& time);
 
  virtual bool integrateEdge(Integrand& integrand, int lEdge,
                             GlobalIntegral& glbInt, const TimeDomain& time);
 
protected:
  bool integrate(Integrand& integrand, GlobalIntegral& glbInt,
                 const TimeDomain& time, const Real3DMat& itgPts);
 
  bool integrate(Integrand& integrand, GlobalIntegral& glbInt,
                 const TimeDomain& time, const ASM::InterfaceChecker& iChk);
 
public:
 
  // Post-processing methods
  // =======================
 
  virtual int evalPoint(const double* xi, double* param, Vec3& X) const;
 
  virtual int findElementContaining(const double* param) const;
 
  virtual double findPoint(Vec3& X, double* param) const;
 
  virtual bool getGridParameters(RealArray& prm, int dir, int nSegSpan) const;
 
  virtual bool tesselate(ElementBlock& grid, const int* npe) const;
 
  virtual bool evalSolution(Matrix& sField, const Vector& locSol,
                            const int* npe, int n_f, bool piola) const;
 
  virtual bool evalSolution(Matrix& sField, const Vector& locSol,
                            const RealArray* gpar, bool regular = true,
                            int deriv = 0, int = 0) const;
 
  virtual bool evalProjSolution(Matrix& sField, const Vector& locSol,
                                const int* npe, int n_f) const;
 
  virtual bool evaluate(const ASMbase* basis, const Vector& locVec,
                        RealArray& vec, int basisNum) const;
 
  virtual bool evaluate(const Field* field, RealArray& vec, int basisNum) const;
 
  virtual bool evaluate(const FunctionBase* func, RealArray& vec,
                        int basisNum, double time) const;
 
  virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
                            const int* npe = nullptr, char project = 0) const;
 
private:
  Go::SplineVolume* projectSolution(const IntegrandBase& integrand) const;
 
  Go::SplineVolume* projectSolutionLocal(const IntegrandBase&) const;
 
  Go::SplineVolume* projectSolutionLocalApprox(const IntegrandBase&) const;
 
  Go::SplineVolume* projectSolutionLeastSquare(const IntegrandBase&) const;
 
public:
  virtual Go::GeomObject* evalSolution(const IntegrandBase& integrand) const;
 
  virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
                            const RealArray* gpar, bool regular = true) const;
 
protected:
 
  // Internal utility methods
  // ========================
 
  virtual void findBoundaryElms(IntVec& elms, int lIndex, int = 0) const;
 
  virtual bool assembleL2matrices(SystemMatrix& A, SystemVector& B,
                                  const L2Integrand& integrand,
                                  bool continuous) const;
 
  bool connectBasis(int face, ASMs3D& neighbor, int nface, int norient,
                    int basis = 1, int slave = 0, int master = 0,
                    bool coordCheck = true, int thick = 1);
 
  void getGaussPointParameters(RealArray& uGP, int dir,
                               int nGauss, const double* xi) const;
 
  bool getGrevilleParameters(RealArray& prm, int dir, int basisNum = 1) const;
 
  bool getQuasiInterplParameters(RealArray& prm, int dir) const;
 
  bool getElementCoordinatesPrm(Matrix& X, double u, double v, double w) const;
 
  double getParametricVolume(int iel) const;
  double getParametricArea(int iel, int dir) const;
 
  virtual void getElementBorders(int iel, double* u) const;
 
  void getElementBorders(int i1, int i2, int i3,
                         double* u, double* v, double* w) const;
 
  double getElementCorners(int i1, int i2, int i3, std::vector<Vec3>& XC,
                           RealArray* uC = nullptr) const;
  void getCornerPoints(int i1, int i2, int i3,
                       std::vector<utl::Point>& XC) const;
 
  void generateThreadGroups(size_t strip1, size_t strip2, size_t strip3,
                            bool silence, bool ignoreGlobalLM);
 
  void generateThreadGroups(size_t strip1, size_t strip2, size_t strip3,
                            char lIndex, bool silence, bool);
 
  virtual void changeNumThreads();
 
  virtual int getFirstItgElmNode() const { return 0; }
  virtual int getLastItgElmNode() const;
 
public:
  virtual void generateThreadGroups(const Integrand& integrand, bool silence,
                                    bool ignoreGlobalLM);
 
  virtual void generateThreadGroups(char lIndex, bool silence, bool);
 
  virtual void generateThreadGroupsFromElms(const IntVec& elms);
 
  virtual void generateProjThreadGroupsFromElms(const IntVec& elms);
 
  static void scatterInd(int n1, int n2, int n3, int p1, int p2, int p3,
                         const int* start, IntVec& index);
 
  virtual bool getOrder(int& p1, int& p2, int& p3) const;
 
  virtual bool getSize(int& n1, int& n2, int& n3, int basis = 0) const;
 
  virtual bool getNoStructElms(int& n1, int& n2, int& n3) const;
 
  virtual void getElmConnectivities(IntMat& neigh,
                                    int basis = ASM::INTEGRATION_BASIS) const;
 
  virtual IntMat getElmNodes(int basis) const;
 
  virtual size_t getNoBoundaryElms(char lIndex, char ldim) const;
 
  virtual size_t getNoNodes(int basis = 0) const;
  virtual size_t getNoProjectionNodes() const;
 
  virtual bool getParameterDomain(Real2DMat& u, IntVec* corners) const;
 
  virtual void evaluateBasis(double u, double v, double w, Vector& N) const;
 
  void extractBasis(double u, double v, double w, Vector& N, Matrix& dNdu,
                    bool fromRight = true) const;
  void extractBasis(double u, double v, double w, Vector& N,
                    Matrix& dNdu, Matrix3D& d2Ndu2,
                    bool fromRight = true) const;
  void extractBasis(double u, double v, double w, int dir, int p, Vector& dN,
                    bool fromRight = true) const;
 
  virtual Field* getProjectedField(const Vector& coefs) const;
 
  virtual Fields* getProjectedFields(const Vector& coefs, size_t = 0) const;
 
private:
  int coeffInd(size_t inod) const;
 
  int findStartNode(int& n1, int& n2, int& n3, char basis) const;
 
  bool getFaceSize(int& n1, int& n2, int basis, int face) const;
 
  static void createMNPC(const Go::SplineVolume* svol, IntMat& MNPC);
 
protected:
  std::shared_ptr<Go::SplineVolume> svol;  
  bool                              swapW; 
 
  const IndexVec& nodeInd; 
  IndexVec      myNodeInd; 
 
  std::map<size_t,size_t> xnMap; 
  std::map<size_t,size_t> nxMap; 
 
  std::vector<DirichletFace> dirich;
 
  ThreadGroups                threadGroupsVol;
  std::map<char,ThreadGroups> threadGroupsFace;
 
  std::vector<std::unique_ptr<BasisFunctionCache>> myCache;
};
 
#endif