1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle 2 // 3 // See the LICENSE.txt file in the Gmsh root directory for license information. 4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues. 5 6 #ifndef BEZIER_BASIS_H 7 #define BEZIER_BASIS_H 8 9 #include <set> 10 #include <map> 11 #include <vector> 12 #include "fullMatrix.h" 13 #include "FuncSpaceData.h" 14 #include "BasisFactory.h" 15 16 class MElement; 17 class bezierBasisRaiser; 18 class bezierCoeff; 19 20 class bezierBasis { 21 private: 22 // Number of corner coeff which are 'real' values of the expanded function 23 int _numLagCoeff; 24 int _dimSimplex; 25 const FuncSpaceData _funcSpaceData; 26 bezierBasisRaiser *_raiser; 27 fullMatrix<double> _exponents; 28 fullMatrix<double> _matrixLag2Bez; 29 fullVector<double> _ordered1dBezPoints; 30 fullMatrix<double> _samplingPntsLagDomain; 31 32 friend class bezierBasisRaiser; 33 friend class bezierCoeff; 34 35 public: 36 // Constructors 37 bezierBasis(FuncSpaceData data); 38 ~bezierBasis(); 39 40 // get methods getDim()41 inline int getDim() const { return _exponents.size2(); } getType()42 inline int getType() const { return _funcSpaceData.getType(); } getOrder()43 inline int getOrder() const { return _funcSpaceData.getSpaceOrder(); } getDimSimplex()44 inline int getDimSimplex() const { return _dimSimplex; } getNumCoeff()45 inline int getNumCoeff() const { return _exponents.size1(); } getNumLagCoeff()46 inline int getNumLagCoeff() const { return _numLagCoeff; } getFuncSpaceData()47 inline FuncSpaceData getFuncSpaceData() const { return _funcSpaceData; } 48 const bezierBasisRaiser *getRaiser() const; 49 getSamplingPointsToComputeBezierCoeff()50 inline const fullMatrix<double> &getSamplingPointsToComputeBezierCoeff() const 51 { 52 return _samplingPntsLagDomain; 53 } 54 55 private: 56 void _construct(); 57 void _constructPyr(); 58 }; 59 60 class bezierBasisRaiser { 61 // Let f, g, h be three function whose Bezier coefficients are given. 62 // This class allows to compute the Bezier coefficients of f*g and f*g*h. 63 private: 64 class _data { 65 friend class bezierBasisRaiser; 66 67 private: 68 int i, j, k; 69 double val; 70 71 public: 72 _data(double vv, int ii, int jj = -1, int kk = -1) i(ii)73 : i(ii), j(jj), k(kk), val(vv) 74 { 75 } 76 }; 77 std::vector<std::vector<_data> > _raiser2, _raiser3; 78 const bezierBasis *_bfs; 79 80 public: bezierBasisRaiser(const bezierBasis * bezier)81 bezierBasisRaiser(const bezierBasis *bezier) : _bfs(bezier) 82 { 83 _fillRaiserData(); 84 }; 85 86 void computeCoeff(const fullVector<double> &coeffA, 87 const fullVector<double> &coeffB, 88 fullVector<double> &coeffSquare) const; 89 void computeCoeff(const fullVector<double> &coeffA, 90 const fullVector<double> &coeffB, 91 const fullVector<double> &coeffC, 92 fullVector<double> &coeffCubic) const; 93 94 private: 95 void _fillRaiserData(); 96 void _fillRaiserDataPyr(); 97 }; 98 99 class bezierCoeffMemoryPool { 100 // This class is to avoid multiple allocation / deallocation during 101 // the subdivision algorithm. 102 private: 103 std::vector<double> _memory; 104 std::size_t _sizeBlocks; 105 std::size_t _numUsedBlocks; 106 std::size_t _currentIndexOfSearch; 107 std::size_t _endOfSearch; 108 // if a reallocation is performed, the pointers must be updated, we need to 109 // know which bezierCoeff have to be updated: 110 std::vector<bezierCoeff *> _bezierCoeff; 111 112 public: 113 bezierCoeffMemoryPool(); ~bezierCoeffMemoryPool()114 ~bezierCoeffMemoryPool() {} 115 116 // before to be used, the size of the blocks has to be specified 117 void setSizeBlocks(std::size_t size); 118 119 double *giveBlock(bezierCoeff *bez); // gives a block of size _sizeBlocks[num] 120 void releaseBlock(double *block, bezierCoeff *bez); 121 void freeMemory(); 122 123 private: 124 void _checkEnoughMemory(); 125 }; 126 127 class bezierCoeff { 128 private: 129 int _numPool; 130 FuncSpaceData _funcSpaceData; 131 const bezierBasis *_basis; 132 int _r, _c; // size of the matrix 133 double *_data; // pointer on the first element 134 bool _ownData; // to know if data should be freed when object is deleted 135 136 static bezierCoeffMemoryPool *_pool0; 137 static bezierCoeffMemoryPool *_pool1; 138 static fullMatrix<double> _sub; 139 // FIXME: not thread safe. We shoud use one pool and one _sub per thread. 140 // The best would be to give the pool to the constructor. 141 // (the pools should be created and deleted e.g. by the plugin 142 // AnalyseMeshQuality) 143 144 public: bezierCoeff()145 bezierCoeff(){}; 146 bezierCoeff(const bezierCoeff &other, bool swap = false); 147 148 bezierCoeff(const FuncSpaceData fsData, 149 const fullVector<double> &orderedLagCoeff, int numOfPool = -1); 150 bezierCoeff(const FuncSpaceData fsData, 151 const fullMatrix<double> &orderedLagCoeff, int numOfPool = -1); 152 // [orderedLagCoeff] : the Lagrange coefficients in the order given by 153 // function gmshGenerateOrderedPoints(..) 154 // -> use bezierBasis::getSamplingPointsToComputeBezierCoeff(..) method 155 // to get the sampling points at which compute those coefficients. 156 // [numOfPool] : the number of the pool (0 or 1) that should be used. 157 // To activate this functionality, first call usePools(..) function. 158 159 ~bezierCoeff(); 160 161 static void usePools(std::size_t size0, std::size_t size1); 162 static void releasePools(); 163 void updateDataPtr(long diff); 164 getPolynomialOrder()165 inline int getPolynomialOrder() const 166 { 167 return _funcSpaceData.getSpaceOrder(); 168 } getNumCoeff()169 inline int getNumCoeff() const { return _r; } getNumColumns()170 inline int getNumColumns() const { return _c; } getNumCornerCoeff()171 inline int getNumCornerCoeff() const { return _basis->getNumLagCoeff(); } 172 int getIdxCornerCoeff(int i) const; getCornerCoeff(int k)173 inline double getCornerCoeff(int k) const 174 { 175 return _data[getIdxCornerCoeff(k)]; 176 } getCornerCoeff(int k,int j)177 inline double getCornerCoeff(int k, int j) const 178 { 179 return _data[getIdxCornerCoeff(k) + _r * j]; 180 } 181 void getCornerCoeffs(fullVector<double> &) const; 182 void getCornerCoeffs(fullMatrix<double> &) const; getDataPtr()183 inline double *getDataPtr() { return _data; } getBezierBasis()184 inline const bezierBasis *getBezierBasis() const { return _basis; } 185 setMatrixAsProxy(fullMatrix<double> & m)186 inline void setMatrixAsProxy(fullMatrix<double> &m) const 187 { 188 m.setAsProxy(_data, _r, _c); 189 } setVectorAsProxy(fullVector<double> & v)190 inline void setVectorAsProxy(fullVector<double> &v) const 191 { 192 v.setAsProxy(_data, _r); 193 } 194 195 void subdivide(std::vector<bezierCoeff *> &subCoeff) const; 196 operator()197 inline double operator()(int i) const { return _data[i]; } operator()198 inline double operator()(int i, int j) const { return _data[i + _r * j]; } 199 operator()200 inline double &operator()(int i) { return _data[i]; } operator()201 inline double &operator()(int i, int j) { return _data[i + _r * j]; } 202 203 private: 204 enum SubdivisionTet { 205 subdivU, 206 subdivV, 207 subdivW, 208 node0CrossEdge12, 209 node3CrossEdge12, 210 node1CrossEdge03, 211 node2CrossEdge03 212 }; 213 static void _subdivideTet(SubdivisionTet which, int n, bezierCoeff &coeff); 214 215 static void _subdivide(fullMatrix<double> &coeff, int npts, int start); 216 static void _subdivide(fullMatrix<double> &coeff, int npts, int start, 217 int inc); 218 static void _subdivideTriangle(const bezierCoeff &coeff, int start, 219 std::vector<bezierCoeff *> &subCoeff); 220 static void _subdivideTetrahedron(const bezierCoeff &coeff, 221 std::vector<bezierCoeff *> &vSubCoeff); 222 static void _subdivideQuadrangle(const bezierCoeff &coeff, 223 std::vector<bezierCoeff *> &subCoeff); 224 static void _subdivideHexahedron(const bezierCoeff &coeff, 225 std::vector<bezierCoeff *> &subCoeff); 226 static void _subdividePrism(const bezierCoeff &coeff, 227 std::vector<bezierCoeff *> &subCoeff); 228 static void _subdividePyramid(const bezierCoeff &coeff, 229 std::vector<bezierCoeff *> &subCoeff); 230 static void _copy(const bezierCoeff &from, int start, int num, 231 bezierCoeff &to); 232 static void _copyLine(const fullMatrix<double> &allSub, int n, int starti, 233 bezierCoeff &sub); 234 static void _copyQuad(const fullMatrix<double> &allSub, int n, int starti, 235 int startj, bezierCoeff &sub); 236 static void _copyHex(const fullMatrix<double> &allSub, int n, int starti, 237 int startj, int startk, bezierCoeff &sub); 238 static void _copyPyr(const fullMatrix<double> &allSub, int nij, int nk, 239 int starti, int startj, int startk, bezierCoeff &sub); _ij2Index(int i,int j,int n)240 inline static int _ij2Index(int i, int j, int n) 241 { 242 return i + j * n - j * (j - 1) / 2; 243 } _ijk2Index(int i,int j,int k,int n)244 inline static int _ijk2Index(int i, int j, int k, int n) 245 { 246 // the whole tet - the top tet + current triangle 247 return (n + 2) * (n + 1) * n / 6 - (n - k + 2) * (n - k + 1) * (n - k) / 6 + 248 j * (n - k) - j * (j - 1) / 2 + i; 249 } 250 void _computeCoefficients(const double *lagCoeffData); 251 }; 252 253 #endif 254