ASMu2D.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef _ASM_U2D_H
#define _ASM_U2D_H
 
#include "ASMLRSpline.h"
#include "ASM2D.h"
#include "BasisFunctionCache.h"
#include "Interface.h"
#include "LRSpline/LRSpline.h"
#include "ThreadGroups.h"
#include <memory>
 
class FiniteElement;
 
namespace utl {
  class Point;
}
 
namespace Go {
  class SplineSurface;
  struct BasisPtsSf;
  struct BasisDerivsSf;
  struct BasisDerivsSf2;
  struct BasisDerivsSf3;
}
 
namespace LR {
  class LRSplineSurface;
}
 
 
class ASMu2D : public ASMLRSpline, public ASM2D
{
protected:
  class BasisFunctionCache : public ::BasisFunctionCache<2>
  {
  public:
    BasisFunctionCache(const ASMu2D& pch);
 
    BasisFunctionCache(const BasisFunctionCache& cache, int b);
 
    virtual ~BasisFunctionCache() = default;
 
  protected:
    bool internalInit() override;
 
    BasisFunctionVals calculatePt(size_t el, size_t gp, bool reduced) const override;
 
    void calculateAll() override;
 
  protected:
    const ASMu2D& patch; 
 
  private:
    bool setupQuadrature();
  };
 
public:
  class InterfaceChecker : public ASM::InterfaceChecker
  {
  public:
    explicit InterfaceChecker(const ASMu2D& pch);
    virtual ~InterfaceChecker() {}
    virtual short int hasContribution(int iel, int = -1,
                                      int = -1, int = -1) const;
    const RealArray& getIntersections(int iel, int edge,
                                      int* cont = nullptr) const;
 
  protected:
    const ASMu2D& myPatch; 
 
    struct Intersection {
      int continuity = 0; 
      RealArray pts;      
    };
 
    std::map<int,Intersection> intersections;
  };
 
  ASMu2D(unsigned char n_s = 2, unsigned char n_f = 2);
  ASMu2D(const ASMu2D& patch, unsigned char n_f = 0);
  virtual ~ASMu2D() { geomB = nullptr; }
 
  virtual const LR::LRSplineSurface* getBasis(int basis = 1) const;
  virtual LR::LRSplineSurface* getBasis(int basis = 1);
 
 
  // Methods for model generation and refinement
  // ===========================================
 
  virtual bool read(std::istream&);
  virtual bool write(std::ostream&, int) const;
 
  virtual bool generateFEMTopology();
 
  virtual void clear(bool retainGeometry = false);
 
  virtual bool getElementCoordinates(Matrix& X, int iel,
                                     bool forceItg = false) const;
 
  virtual void getElmConnectivities(IntMat& neighs,
                                    int basis = ASM::INTEGRATION_BASIS) const;
 
  virtual IntMat getElmNodes(int basis) const;
 
  virtual void getNodalCoordinates(Matrix& X, bool = false) const
  { this->getCoordinates(X, nsd, *lrspline); }
 
  virtual Vec3 getCoord(size_t inod) const;
 
  virtual bool updateCoords(const Vector& displ);
 
  virtual void getBoundaryNodes(int lIndex, IntVec& nodes, int basis, int = 1,
                                int orient = 0, bool local = false) const;
 
  virtual bool getOrder(int& p1, int& p2, int& p3) const;
 
  virtual size_t getNoNodes(int basis = 0) const;
 
  virtual size_t getNoProjectionNodes() const;
 
  bool diagonalRefine(int minBasisfunctions);
  bool cornerRefine(int minBasisfunctions);
  bool uniformRefine(int minBasisfunctions);
  virtual bool uniformRefine(int dir, int nInsert);
  using ASMLRSpline::refine;
  virtual bool refine(int dir, const RealArray& xi, double scale);
  virtual bool refine(const LR::RefineData& prm, Vectors& sol);
  virtual bool raiseOrder(int ru, int rv);
 
  virtual bool createProjectionBasis(bool init);
 
  virtual void setMinimumSize(double size) { aMin = size; }
  virtual double getMinimumSize(int nrefinements) const;
  virtual bool checkElementSize(int elmId, bool globalNum = true) const;
 
  virtual void extractElmRes(const Matrix& globRes, Matrix& elmRes,
                             size_t internalFirst) const;
 
  void copyRefinement(LR::LRSplineSurface* basis, int multiplicity = 1) const;
 
 
  // Various methods for preprocessing of boundary conditions and patch topology
  // ===========================================================================
 
  virtual void constrainEdge(int dir, bool open, int dof, int code, char basis);
  virtual size_t constrainEdgeLocal(int dir, bool open, int dof, int code,
                                    bool project = false);
 
  virtual void constrainCorner(int I, int J, int dof,
                               int code = 0, char basis = 1);
  virtual void constrainNode(double xi, double eta, int dof, int code = 0);
 
  virtual bool connectPatch(int edge, ASM2D& neighbor, int nedge, bool revers,
                            int = 0, bool coordCheck = true, int thick = 1);
 
 
  // Methods for integration of finite element quantities.
  // These are the main computational methods of the ASM class hierarchy.
  // ====================================================================
 
  using ASMbase::integrate;
  virtual bool integrate(Integrand& integrand,
                         GlobalIntegral& glbInt, const TimeDomain& time);
 
  virtual bool integrate(Integrand& integrand, int lIndex,
                         GlobalIntegral& glbInt, const TimeDomain& time);
 
  virtual bool integrate(Integrand& integrand, GlobalIntegral& glbInt,
                         const TimeDomain& time,
                         const ASM::InterfaceChecker& iChk);
 
  virtual bool diracPoint(Integrand& integrand, GlobalIntegral& glbInt,
                          const double* param, const Vec3& pval);
 
  virtual bool updateDirichlet(const std::map<int,RealFunc*>& func,
                               const std::map<int,VecFunc*>& vfunc, double time,
                               const std::map<int,int>* g2l = nullptr);
 
  virtual int getCorner(int I, int J, int basis) const;
 
protected:
  bool integrate(Integrand& integrand, GlobalIntegral& glbInt,
                 const TimeDomain& time, const Real3DMat& itgPts);
 
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 = false,
                            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 FunctionBase* func, RealArray& vec,
                        int, double time) const;
 
  virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
                            const int* npe, char project = '\0') const;
 
  virtual bool separateProjectionBasis() const;
 
  virtual Field* getProjectedField(const Vector& coefs) const;
 
  virtual Fields* getProjectedFields(const Vector& coefs, size_t = 0) const;
 
  virtual void storeMesh(const std::string& fName, int fType) const;
 
  void setOutputMaster(const ASMu2D* pch)
  { outputMaster = pch; }
 
protected:
  struct DirichletEdge
  {
    LR::LRSplineSurface* lr;       
    LR::parameterEdge   edg;       
    IntVec              MLGE;      
    IntVec              MLGN;      
    IntMat              MNPC;      
    int                 dof;       
    int                 code;      
    int                 corners[2];
 
    DirichletEdge(LR::LRSplineSurface* sf, int dir,
                  int d = 0, int c = 0, int offset = 1);
    bool isCorner(int b) const { return b == corners[0] || b == corners[1]; }
  };
 
  LR::LRSplineSurface* projectSolution(const IntegrandBase& integrand) const;
  LR::LRSplineSurface* scRecovery(const IntegrandBase& integrand) const;
 
  LR::LRSplineSurface* regularInterpolation(const RealArray& upar,
                                            const RealArray& vpar,
                                            const Matrix& points, int basis) const;
 
public:
  virtual LR::LRSpline* evalSolution(const IntegrandBase& integrand) const;
 
  virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
                            const RealArray* gpar, bool = false) const;
 
  bool edgeL2projection(const DirichletEdge& edge, const FunctionBase& values,
                        Real2DMat& result, double time) const;
 
  virtual bool transferGaussPtVars(const LR::LRSpline* old_basis,
                                   const RealArray& oldVar, RealArray& newVar,
                                   int nGauss) const;
  virtual bool transferGaussPtVarsN(const LR::LRSpline* old_basis,
                                    const RealArray& oldVar, RealArray& newVar,
                                    int nGauss) const;
 
  using ASMLRSpline::transferCntrlPtVars;
  virtual bool transferCntrlPtVars(const LR::LRSpline* old_basis,
                                   RealArray& newVar, int nGauss) const;
 
protected:
 
  // Internal utility methods
  // ========================
 
  virtual void findBoundaryElms(IntVec& elms, int lIndex, int orient) const;
 
  virtual bool assembleL2matrices(SystemMatrix& A, SystemVector& B,
                                  const L2Integrand& integrand,
                                  bool continuous) const;
 
  bool connectBasis(int edge, ASMu2D& neighbor, int nedge, bool revers,
                    int basis = 1, int slave = 0, int master = 0,
                    bool coordCheck = true, int thick = 1);
 
  void getGaussPointParameters(RealArray& uGP, int dir, int nGauss,
                               int iel, const double* xi,
                               const LR::LRSplineSurface* spline = nullptr) const;
 
  bool getGrevilleParameters(RealArray& prm, int dir, int basisNum = 1) const;
 
  double getParametricLength(int iel, int dir) const;
 
  double getElementCorners(int iel, std::vector<Vec3>& XC,
                           RealArray* uC = nullptr) const;
  void getCornerPoints(int iel, std::vector<utl::Point>& XC) const;
 
  bool evaluateBasis(int iel, FiniteElement& fe, int derivs = 0) const;
 
  void computeBasis(double u, double v,
                    Go::BasisPtsSf& bas, int iel,
                    const LR::LRSplineSurface* spline = nullptr) const;
  void computeBasis(double u, double v,
                    Go::BasisDerivsSf& bas, int iel,
                    const LR::LRSplineSurface* spline = nullptr) const;
  void computeBasis(double u, double v,
                    Go::BasisDerivsSf2& bas, int iel,
                    const LR::LRSplineSurface* spline = nullptr) const;
  void computeBasis(double u, double v,
                    Go::BasisDerivsSf3& bas, int iel,
                    const LR::LRSplineSurface* spline = nullptr) const;
 
  virtual int evalPoint(int iel, const double* param, Vec3& X) const;
 
  using ASMLRSpline::generateThreadGroups;
  void generateThreadGroups(const Integrand& integrand, bool silence,
                            bool ignoreGlobalLM);
 
  virtual void generateThreadGroupsFromElms(const IntVec& elms);
 
  virtual void generateProjThreadGroupsFromElms(const IntVec& elms);
 
  virtual void changeNumThreads();
 
  virtual void remapErrors(RealArray& errors,
                           const RealArray& origErr, bool elemErrors) const;
 
  virtual void extendRefinementDomain(IntSet& refineIndices,
                                      const IntSet& neighborIndices) const;
 
  std::shared_ptr<LR::LRSplineSurface> createLRfromTensor();
 
  static std::shared_ptr<LR::LRSplineSurface> createLRNurbs(const Go::SplineSurface& srf);
 
  void generateBezierBasis();
 
  void generateBezierExtraction();
 
  void writePostscriptElementsNurbs (std::shared_ptr<LR::LRSplineSurface> mesh,
                                     std::ostream& out, bool close = true,
                                     int nu = 2, int nv = 2);
  void writePostscriptMeshWithControlPointsNurbs (std::shared_ptr<LR::LRSplineSurface> mesh,
                                                  std::ostream& out, int nu = 2, int nv = 2);
 
public:
  virtual size_t getNoBoundaryElms(char lIndex, char ldim) const;
 
  static void computeBasisNurbs(double u, double v,
                                Go::BasisPtsSf& bas, int iel,
                                const LR::LRSplineSurface& spline);
  static void computeBasisNurbs(double u, double v,
                                Go::BasisDerivsSf& bas, int iel,
                                const LR::LRSplineSurface& spline);
  static void computeBasisNurbs(double u, double v,
                                Go::BasisDerivsSf2& bas, int iel,
                                const LR::LRSplineSurface& spline);
  static void computeBasisNurbs(double u, double v,
                                Go::BasisDerivsSf3& bas, int iel,
                                const LR::LRSplineSurface& spline);
 
  bool rational() const { return is_rational; }
 
  static bool getCoordinates(Matrix& X,
                             unsigned char nsd,
                             const LR::LRSplineSurface& spline,
                             int iel = -1);
 
protected:
  std::shared_ptr<LR::LRSplineSurface> lrspline; 
 
  bool is_rational = false; 
 
  const ASMu2D* outputMaster = nullptr; 
 
  Go::SplineSurface* tensorspline; 
  Go::SplineSurface* tensorPrjBas; 
  // The tensor spline object is kept for backward compatability with the REFINE
  // and RAISEORDER key-words, although we take note that there is a possibility
  // of optimization since all mapping values and Jacobians may be performed on
  // this object for increased efficiency.
 
  std::vector<DirichletEdge> dirich;
 
  const Matrices& bezierExtract; 
  Matrices      myBezierExtract; 
 
  Go::BsplineBasis bezier_u; 
  Go::BsplineBasis bezier_v; 
 
  std::vector<std::unique_ptr<BasisFunctionCache>> myCache;
 
private:
  mutable double aMin; 
 
  bool evaluateBasisNurbs(int iel, FiniteElement& fe,
                          int derivs) const;
};
 
#endif