ASMu3D.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef _ASM_U3D_H
#define _ASM_U3D_H
 
#include "ASMLRSpline.h"
#include "ASM3D.h"
#include "BasisFunctionCache.h"
#include "LRSpline/LRSpline.h"
#include "ThreadGroups.h"
#include <memory>
 
class FiniteElement;
 
namespace utl {
  class Point;
}
 
namespace Go {
  class SplineVolume;
}
 
namespace LR {
  class LRSplineVolume;
}
 
 
class ASMu3D : public ASMLRSpline, public ASM3D
{
protected:
  class BasisFunctionCache : public ::BasisFunctionCache<3>
  {
  public:
    BasisFunctionCache(const ASMu3D& pch, bool useBezier = true);
 
    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;
 
    BasisFunctionVals calculatePrm(FiniteElement& fe,
                                   const std::array<double,3>& du,
                                   size_t el, size_t gp, bool reduced) const;
 
    void calculateAll() override;
 
    struct BezierExtract {
      Matrix N;     
      Matrix dNdu; 
      Matrix dNdv; 
      Matrix dNdw; 
    };
 
    bool bezierEnabled; 
    BezierExtract mainB; 
    BezierExtract reducedB; 
 
    const ASMu3D& patch; 
 
  private:
    bool setupQuadrature();
  };
 
public:
  explicit ASMu3D(unsigned char n_f = 3);
  ASMu3D(const ASMu3D& patch, unsigned char n_f = 0);
  virtual ~ASMu3D() {}
 
  virtual const LR::LRSplineVolume* getBasis(int basis = 1) const;
  virtual LR::LRSplineVolume* 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;
 
  virtual Vec3 getCoord(size_t inod) const;
 
  virtual bool updateCoords(const Vector& displ);
 
  virtual void getBoundaryNodes(int lIndex, IntVec& nodes, int basis,
                                int, int orient, bool local) const;
 
  virtual void getBoundary1Nodes(int lEdge, IntVec& nodes, int basis,
                                 int orient, bool local,
                                 bool open = false) const;
 
  virtual int getCorner(int I, int J, int K, int basis) const;
 
  virtual bool getOrder(int& p1, int& p2, int& p3) const;
 
  virtual size_t getNoNodes(int basis = 0) const;
 
  virtual size_t getNoProjectionNodes() const;
 
  virtual bool uniformRefine(int dir, int nInsert);
  using ASMLRSpline::refine;
  virtual bool refine(int dir, const RealArray& xi);
  virtual bool refine(const LR::RefineData& prm, Vectors& sol);
  virtual bool raiseOrder(int ru, int rv, int rw, bool setOrder);
 
  virtual bool createProjectionBasis(bool init);
 
  virtual void setMinimumSize(double size) { vMin = size; }
  virtual double getMinimumSize(int nrefinements) const;
  virtual bool checkElementSize(int elmId, bool globalNum = true) const;
 
  void copyRefinement(LR::LRSplineVolume* basis, int multiplicity = 1) 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);
 
 
  // Methods for integration of finite element quantities.
  // These are the main computational methods of the ASM class hierarchy.
  // ====================================================================
 
  virtual void getNoBouPoints(size_t& nPt, char ldim, char lindx);
 
  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);
 
  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);
 
  // 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;
 
  void setOutputMaster(const ASMu3D* pch)
  { outputMaster = pch; }
 
  virtual void extractElmRes(const Matrix& globRes, Matrix& elmRes,
                             size_t internalFirst) const;
 
protected:
  struct DirichletFace
  {
    LR::LRSplineVolume* lr;        
    LR::parameterEdge   edg;       
    IntVec              MLGE;      
    IntVec              MLGN;      
    IntMat              MNPC;      
    int                 dof;       
    int                 code;      
    int                 corners[4];
 
    DirichletFace(LR::LRSplineVolume* sv, int dir,
                  int d = 0, int c = 0, int offset = 1);
    bool isCorner(int b) const;
  };
 
  LR::LRSplineVolume* projectSolution(const IntegrandBase& integrand) const;
  LR::LRSplineVolume* scRecovery(const IntegrandBase& integrand) const;
 
  LR::LRSplineVolume* regularInterpolation(const RealArray& upar,
                                           const RealArray& vpar,
                                           const RealArray& wpar,
                                           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 faceL2projection(const DirichletFace& face, 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;
 
  using ASMLRSpline::generateThreadGroups;
  virtual void generateThreadGroups(const Integrand& integrand, bool silence,
                                    bool ignoreGlobalLM);
 
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 face, ASMu3D& 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,
                               int iel, const double* xi,
                               const LR::LRSplineVolume* spline = nullptr) const;
 
  bool getGrevilleParameters(RealArray& prm, int dir, int basisNum = 1) const;
 
  bool getQuasiInterplParameters(RealArray& prm, int dir) const;
 
  double getParametricArea(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;
 
  void evaluateBasis(int iel, double u, double v, double w,
                     Vector& N, Matrix& dNdu, int basis = 1) const;
 
  void evaluateBasis(Vector& N, Matrix& dNdu,
                     const Matrix& C, const Matrix& B) const;
 
  void evaluateBasis(int iel, FiniteElement& fe, Matrix& dNdu,
                     int basis = 1) const;
 
  void evaluateBasis(int iel, double u, double v, double w,
                     Vector& N, Matrix& dNdu, Matrix3D& d2Ndu2, int basis = 1) const;
 
  virtual int evalPoint(int iel, const double* param, Vec3& X) const;
 
  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::LRSplineVolume> createLRfromTensor();
 
public:
  virtual size_t getNoBoundaryElms(char lIndex, char ldim) const;
 
protected:
  std::shared_ptr<LR::LRSplineVolume> lrspline; 
 
  const ASMu3D* outputMaster = nullptr; 
 
  Go::SplineVolume* tensorspline; 
  Go::SplineVolume* 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<DirichletFace> dirich;
 
  const Matrices& bezierExtract; 
  Matrices      myBezierExtract; 
 
  std::vector<std::unique_ptr<BasisFunctionCache>> myCache;
 
private:
  mutable double vMin; 
};
 
#endif