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