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 CARTESIAN_H
7 #define CARTESIAN_H
8 
9 #include <set>
10 #include <map>
11 #include <vector>
12 #include <stdio.h>
13 #include "SVector3.h"
14 #include "SPoint3.h"
15 #include "GmshMessage.h"
16 #include "MHexahedron.h"
17 #include "MVertex.h"
18 
19 // A cartesian grid that encompasses an oriented 3-D box, with values
20 // stored at vertices:
21 //
22 //               j
23 //   +---+---+---+---+---+---+
24 //   |   |   |   |   |   |   |
25 // i +---+---+-(i,j)-+---+---+
26 //   |   |   |   |   |   |   |
27 //   +---+---+---+---+---+---+
28 //
29 //   Nodal values and active hexahedral cells are stored and
30 //   referenced by a linear index (i + N * j)
31 //
32 //   The (i,j) cell has nodes (i,j), (i+1,j), (i+1,j+1) and (i,j+1)
33 template <class scalar> class cartesianBox {
34 private:
35   // number of subdivisions along the xi-, eta- and zeta-axis
36   int _nxi, _neta, _nzeta;
37   // origin of the grid and spacing along xi, eta and zeta
38   double _x0, _y0, _z0, _dxi, _deta, _dzeta;
39   // xi-, eta- and zeta-axis directions
40   SVector3 _xiAxis, _etaAxis, _zetaAxis;
41   // set of active cells; the value stored for cell (i,j,k) is the
42   // linear index (i + _nxi * j + _nxi *_neta * k)
43   std::set<int> _activeCells;
44   // map of stored nodal values, indexed by the linear index (i +
45   // (_nxi+1) * j + (_nxi+1) * (_neta+1) * k). Along with the value is
46   // stored a node tag (used for global numbering of the nodes across
47   // the grid levels)
48   typename std::map<int, std::pair<scalar, int> > _nodalValues;
49   // level of the box (coarsest box has highest level; finest box has
50   // level==1)
51   int _level;
52   // pointer to a finer (refined by 2) level box, if any
53   cartesianBox<scalar> *_childBox;
_getNumNodes()54   int _getNumNodes()
55   {
56     int n = 0;
57     for(valIter it = _nodalValues.begin(); it != _nodalValues.end(); it++)
58       if(it->second.second > 0) n++;
59     if(_childBox) n += _childBox->_getNumNodes();
60     return n;
61   }
_printNodes(FILE * f)62   void _printNodes(FILE *f)
63   {
64     for(valIter it = _nodalValues.begin(); it != _nodalValues.end(); ++it) {
65       if(it->second.second > 0) {
66         SPoint3 p = getNodeCoordinates(it->first);
67         fprintf(f, "%d %g %g %g\n", it->second.second, p.x(), p.y(), p.z());
68       }
69     }
70     if(_childBox) _childBox->_printNodes(f);
71   }
_getNumElements(bool simplex)72   int _getNumElements(bool simplex)
73   {
74     int coeff = simplex ? 6 : 1;
75     int n = _activeCells.size() * coeff;
76     if(_childBox) n += _childBox->_getNumElements(simplex);
77     return n;
78   }
79   void _printElements(FILE *f, bool simplex, int startingNum = 1)
80   {
81     int num = startingNum;
82     for(auto it = _activeCells.begin(); it != _activeCells.end(); ++it) {
83       int i, j, k;
84       getCellIJK(*it, i, j, k);
85       if(!simplex) {
86         fprintf(f, "%d 5 3 1 1 1 %d %d %d %d %d %d %d %d\n", num++,
87                 std::abs(getNodeTag(getNodeIndex(i, j, k))),
88                 std::abs(getNodeTag(getNodeIndex(i + 1, j, k))),
89                 std::abs(getNodeTag(getNodeIndex(i + 1, j + 1, k))),
90                 std::abs(getNodeTag(getNodeIndex(i, j + 1, k))),
91                 std::abs(getNodeTag(getNodeIndex(i, j, k + 1))),
92                 std::abs(getNodeTag(getNodeIndex(i + 1, j, k + 1))),
93                 std::abs(getNodeTag(getNodeIndex(i + 1, j + 1, k + 1))),
94                 std::abs(getNodeTag(getNodeIndex(i, j + 1, k + 1))));
95       }
96       else {
97         int idx[6][4] = {
98           {getNodeIndex(i, j + 1, k), getNodeIndex(i, j + 1, k + 1),
99            getNodeIndex(i + 1, j, k + 1), getNodeIndex(i + 1, j + 1, k + 1)},
100           {getNodeIndex(i, j + 1, k), getNodeIndex(i + 1, j + 1, k + 1),
101            getNodeIndex(i + 1, j, k + 1), getNodeIndex(i + 1, j + 1, k)},
102           {getNodeIndex(i, j + 1, k), getNodeIndex(i, j, k + 1),
103            getNodeIndex(i + 1, j, k + 1), getNodeIndex(i, j + 1, k + 1)},
104           {getNodeIndex(i, j + 1, k), getNodeIndex(i + 1, j + 1, k),
105            getNodeIndex(i + 1, j, k + 1), getNodeIndex(i + 1, j, k)},
106           {getNodeIndex(i, j + 1, k), getNodeIndex(i + 1, j, k),
107            getNodeIndex(i + 1, j, k + 1), getNodeIndex(i, j, k)},
108           {getNodeIndex(i, j + 1, k), getNodeIndex(i, j, k),
109            getNodeIndex(i + 1, j, k + 1), getNodeIndex(i, j, k + 1)}};
110         for(int ii = 0; ii < 6; ii++) {
111           fprintf(
112             f, "%d 4 3 1 1 1 %d %d %d %d\n", num++,
113             std::abs(getNodeTag(idx[ii][0])), std::abs(getNodeTag(idx[ii][1])),
114             std::abs(getNodeTag(idx[ii][2])), std::abs(getNodeTag(idx[ii][3])));
115         }
116       }
117     }
118     if(_childBox) _childBox->_printElements(f, simplex, num);
119   }
_printValues(FILE * f)120   void _printValues(FILE *f)
121   {
122     for(valIter it = _nodalValues.begin(); it != _nodalValues.end(); ++it) {
123       if(it->second.second > 0)
124         fprintf(f, "%d %g\n", it->second.second, it->second.first);
125     }
126     if(_childBox) _childBox->_printValues(f);
127   }
128 
129 public:
130   cartesianBox(double X0, double Y0, double Z0, const SVector3 &dxi,
131                const SVector3 &deta, const SVector3 &dzeta, int Nxi, int Neta,
132                int Nzeta, int level = 1)
_nxi(Nxi)133     : _nxi(Nxi), _neta(Neta), _nzeta(Nzeta), _x0(X0), _y0(Y0), _z0(Z0),
134       _dxi(norm(dxi)), _deta(norm(deta)), _dzeta(norm(dzeta)), _xiAxis(dxi),
135       _etaAxis(deta), _zetaAxis(dzeta), _level(level), _childBox(0)
136   {
137     _xiAxis.normalize();
138     _etaAxis.normalize();
139     _zetaAxis.normalize();
140     if(level > 1)
141       _childBox = new cartesianBox<scalar>(
142         X0, Y0, Z0, dxi, deta, dzeta, 2 * Nxi, 2 * Neta, 2 * Nzeta, level - 1);
143   }
getLC()144   double getLC() { return sqrt(_dxi * _dxi + _deta * _deta + _dzeta * _dzeta); }
getNxi()145   int getNxi() { return _nxi; }
getNeta()146   int getNeta() { return _neta; }
getNzeta()147   int getNzeta() { return _nzeta; }
getChildBox()148   cartesianBox<scalar> *getChildBox() { return _childBox; }
getLevel()149   int getLevel() { return _level; }
150   typedef std::set<int>::const_iterator cellIter;
activeCellsBegin()151   cellIter activeCellsBegin() { return _activeCells.begin(); }
activeCellsEnd()152   cellIter activeCellsEnd() { return _activeCells.end(); }
153   typedef typename std::map<int, std::pair<scalar, int> >::iterator valIter;
nodalValuesBegin()154   valIter nodalValuesBegin() { return _nodalValues.begin(); }
nodalValuesEnd()155   valIter nodalValuesEnd() { return _nodalValues.end(); }
setNodalValue(int i,scalar s)156   void setNodalValue(int i, scalar s) { _nodalValues[i].first = s; }
getNodalValuesForCell(int t,std::vector<scalar> & values)157   void getNodalValuesForCell(int t, std::vector<scalar> &values)
158   {
159     int i, j, k;
160     getCellIJK(t, i, j, k);
161     for(int I = 0; I < 2; I++)
162       for(int J = 0; J < 2; J++)
163         for(int K = 0; K < 2; K++) {
164           valIter it = _nodalValues.find(getNodeIndex(i + I, j + J, k + K));
165           if(it != _nodalValues.end())
166             values.push_back(it->second.first);
167           else {
168             Msg::Error("Could not find value i,j,k=%d,%d,%d for cell %d", i + I,
169                        j + J, k + K, t);
170             values.push_back(0.);
171           }
172         }
173   }
getValueContainingPoint(double x,double y,double z)174   double getValueContainingPoint(double x, double y, double z)
175   {
176     SVector3 DP(x - _x0, y - _y0, z - _z0);
177 
178     int t = getCellContainingPoint(x, y, z);
179     int i, j, k;
180     getCellIJK(t, i, j, k);
181 
182     valIter it1 = _nodalValues.find(getNodeIndex(i, j, k));
183     valIter it2 = _nodalValues.find(getNodeIndex(i + 1, j, k));
184     valIter it3 = _nodalValues.find(getNodeIndex(i + 1, j + 1, k));
185     valIter it4 = _nodalValues.find(getNodeIndex(i, j + 1, k));
186     valIter it5 = _nodalValues.find(getNodeIndex(i, j, k + 1));
187     valIter it6 = _nodalValues.find(getNodeIndex(i + 1, j, k + 1));
188     valIter it7 = _nodalValues.find(getNodeIndex(i + 1, j + 1, k + 1));
189     valIter it8 = _nodalValues.find(getNodeIndex(i, j + 1, k + 1));
190 
191     if(it1 == _nodalValues.end())
192       return _childBox->getValueContainingPoint(x, y, z);
193     if(it2 == _nodalValues.end())
194       return _childBox->getValueContainingPoint(x, y, z);
195     if(it3 == _nodalValues.end())
196       return _childBox->getValueContainingPoint(x, y, z);
197     if(it4 == _nodalValues.end())
198       return _childBox->getValueContainingPoint(x, y, z);
199     if(it5 == _nodalValues.end())
200       return _childBox->getValueContainingPoint(x, y, z);
201     if(it6 == _nodalValues.end())
202       return _childBox->getValueContainingPoint(x, y, z);
203     if(it7 == _nodalValues.end())
204       return _childBox->getValueContainingPoint(x, y, z);
205     if(it8 == _nodalValues.end())
206       return _childBox->getValueContainingPoint(x, y, z);
207 
208     double vals[8];
209     vals[0] = it1->second.first;
210     vals[1] = it2->second.first;
211     vals[2] = it3->second.first;
212     vals[3] = it4->second.first;
213     vals[4] = it5->second.first;
214     vals[5] = it6->second.first;
215     vals[6] = it7->second.first;
216     vals[7] = it8->second.first;
217     // for (int i= 0; i< 8; i++) printf("vals %d = %g \n", i, vals[i]);
218 
219     SPoint3 p1 = getNodeCoordinates(it1->first);
220     SPoint3 p2 = getNodeCoordinates(it2->first);
221     SPoint3 p3 = getNodeCoordinates(it3->first);
222     SPoint3 p4 = getNodeCoordinates(it4->first);
223     SPoint3 p5 = getNodeCoordinates(it5->first);
224     SPoint3 p6 = getNodeCoordinates(it6->first);
225     SPoint3 p7 = getNodeCoordinates(it7->first);
226     SPoint3 p8 = getNodeCoordinates(it8->first);
227 
228     MVertex *v1 = new MVertex(p1.x(), p1.y(), p1.z());
229     MVertex *v2 = new MVertex(p2.x(), p2.y(), p2.z());
230     MVertex *v3 = new MVertex(p3.x(), p3.y(), p3.z());
231     MVertex *v4 = new MVertex(p4.x(), p4.y(), p4.z());
232     MVertex *v5 = new MVertex(p5.x(), p5.y(), p5.z());
233     MVertex *v6 = new MVertex(p6.x(), p6.y(), p6.z());
234     MVertex *v7 = new MVertex(p7.x(), p7.y(), p7.z());
235     MVertex *v8 = new MVertex(p8.x(), p8.y(), p8.z());
236 
237     MHexahedron *newElem = new MHexahedron(v1, v2, v3, v4, v5, v6, v7, v8);
238     double uvw[3];
239     double xyz[3] = {x, y, z};
240     newElem->xyz2uvw(xyz, uvw);
241     // printf("uvw =%g %g %g \n", uvw[0],uvw[1],uvw[2]);
242     double val = newElem->interpolate(vals, uvw[0], uvw[1], uvw[2]);
243 
244     delete newElem;
245     delete v1;
246     delete v2;
247     delete v3;
248     delete v4;
249     delete v5;
250     delete v6;
251     delete v7;
252     delete v8;
253 
254     return val;
255   }
getCellContainingPoint(double x,double y,double z)256   int getCellContainingPoint(double x, double y, double z) const
257   {
258     // P = P_0 + xi * _vdx + eta * _vdy + zeta *vdz
259     // DP = P-P_0 * _vdx = xi
260     SVector3 DP(x - _x0, y - _y0, z - _z0);
261     double xi = dot(DP, _xiAxis);
262     double eta = dot(DP, _etaAxis);
263     double zeta = dot(DP, _zetaAxis);
264     int i = xi / _dxi * _nxi;
265     int j = eta / _deta * _neta;
266     int k = zeta / _dzeta * _nzeta;
267     if(i < 0) i = 0;
268     if(i >= _nxi) i = _nxi - 1;
269     if(j < 0) j = 0;
270     if(j >= _neta) j = _neta - 1;
271     if(k < 0) k = 0;
272     if(k >= _nzeta) k = _nzeta - 1;
273     return getCellIndex(i, j, k);
274   }
getNodeCoordinates(int t)275   SPoint3 getNodeCoordinates(int t) const
276   {
277     int i, j, k;
278     getNodeIJK(t, i, j, k);
279     const double xi = i * _dxi / _nxi;
280     const double eta = j * _deta / _neta;
281     const double zeta = k * _dzeta / _nzeta;
282     SVector3 D = xi * _xiAxis + eta * _etaAxis + zeta * _zetaAxis;
283     return SPoint3(_x0 + D.x(), _y0 + D.y(), _z0 + D.z());
284   }
insertActiveCell(int t)285   void insertActiveCell(int t) { _activeCells.insert(t); }
eraseActiveCell(int t)286   void eraseActiveCell(int t) { _activeCells.erase(t); }
activeCellExists(int t)287   bool activeCellExists(int t)
288   {
289     return (_activeCells.find(t) != _activeCells.end());
290   }
getCellIndex(int i,int j,int k)291   int getCellIndex(int i, int j, int k) const
292   {
293     return i + _nxi * j + _nxi * _neta * k;
294   }
getNodeIndex(int i,int j,int k)295   int getNodeIndex(int i, int j, int k) const
296   {
297     return i + (_nxi + 1) * j + (_nxi + 1) * (_neta + 1) * k;
298   }
getNodeTag(int index)299   int getNodeTag(int index)
300   {
301     valIter it = _nodalValues.find(index);
302     if(it != _nodalValues.end())
303       return it->second.second;
304     else
305       return 0;
306   }
getCellIJK(int index,int & i,int & j,int & k)307   void getCellIJK(int index, int &i, int &j, int &k) const
308   {
309     k = index / (_nxi * _neta);
310     j = (index - k * (_nxi * _neta)) / _nxi;
311     i = (index - k * (_nxi * _neta) - j * _nxi);
312   }
getNodeIJK(int index,int & i,int & j,int & k)313   void getNodeIJK(int index, int &i, int &j, int &k) const
314   {
315     k = index / ((_nxi + 1) * (_neta + 1));
316     j = (index - k * ((_nxi + 1) * (_neta + 1))) / (_nxi + 1);
317     i = (index - k * ((_nxi + 1) * (_neta + 1)) - j * (_nxi + 1));
318   }
createNodalValues()319   void createNodalValues()
320   {
321     for(auto it = _activeCells.begin(); it != _activeCells.end(); ++it) {
322       const int &t = *it;
323       int i, j, k;
324       getCellIJK(t, i, j, k);
325       for(int I = 0; I < 2; I++)
326         for(int J = 0; J < 2; J++)
327           for(int K = 0; K < 2; K++)
328             _nodalValues[getNodeIndex(i + I, j + J, k + K)] =
329               std::pair<scalar, int>(0., 0);
330     }
331     if(_childBox) _childBox->createNodalValues();
332   }
333   void renumberNodes(int startingNum = 1, cartesianBox<scalar> *parent = 0)
334   {
335     int num = startingNum;
336     for(valIter it = _nodalValues.begin(); it != _nodalValues.end(); it++) {
337       int i, j, k;
338       getNodeIJK(it->first, i, j, k);
339       if(!parent || i % 2 || j % 2 || k % 2)
340         it->second.second = num++;
341       else {
342         int tag = parent->getNodeTag(parent->getNodeIndex(i / 2, j / 2, k / 2));
343         if(!tag) // FIXME! not sure why this can happen, but it does (bug?)
344           it->second.second = num++;
345         else // the node exists in the coarset grid: store it with negative sign
346           it->second.second = -std::abs(tag);
347       }
348     }
349     if(_childBox) _childBox->renumberNodes(num, this);
350   }
351   void writeMSH(const std::string &fileName, bool simplex = false,
352                 bool writeNodalValues = true)
353   {
354     FILE *f = fopen(fileName.c_str(), "w");
355     if(!f) {
356       Msg::Error("Could not open file '%s'", fileName.c_str());
357       return;
358     }
359     int numNodes = _getNumNodes(), numElements = _getNumElements(simplex);
360     Msg::Info("Writing '%s' (%d nodes, %d elements)", fileName.c_str(),
361               numNodes, numElements);
362     fprintf(f, "$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
363     fprintf(f, "$Nodes\n%d\n", numNodes);
364     _printNodes(f);
365     fprintf(f, "$EndNodes\n");
366     fprintf(f, "$Elements\n%d\n", numElements);
367     _printElements(f, simplex);
368     fprintf(f, "$EndElements\n");
369     if(writeNodalValues) {
370       fprintf(f, "$NodeData\n1\n\"distance\"\n1\n0.0\n3\n0\n1\n%d\n", numNodes);
371       _printValues(f);
372       fprintf(f, "$EndNodeData\n");
373     }
374     fclose(f);
375   }
376 };
377 
378 #endif
379