gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
cartesian.h
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 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
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
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  }
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  }
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  }
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)
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();
140  if(level > 1)
142  X0, Y0, Z0, dxi, deta, dzeta, 2 * Nxi, 2 * Neta, 2 * Nzeta, level - 1);
143  }
144  double getLC() { return sqrt(_dxi * _dxi + _deta * _deta + _dzeta * _dzeta); }
145  int getNxi() { return _nxi; }
146  int getNeta() { return _neta; }
147  int getNzeta() { return _nzeta; }
149  int getLevel() { return _level; }
150  typedef std::set<int>::const_iterator cellIter;
151  cellIter activeCellsBegin() { return _activeCells.begin(); }
153  typedef typename std::map<int, std::pair<scalar, int> >::iterator valIter;
154  valIter nodalValuesBegin() { return _nodalValues.begin(); }
155  valIter nodalValuesEnd() { return _nodalValues.end(); }
156  void setNodalValue(int i, scalar s) { _nodalValues[i].first = s; }
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  }
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  }
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  }
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  }
285  void insertActiveCell(int t) { _activeCells.insert(t); }
286  void eraseActiveCell(int t) { _activeCells.erase(t); }
287  bool activeCellExists(int t)
288  {
289  return (_activeCells.find(t) != _activeCells.end());
290  }
291  int getCellIndex(int i, int j, int k) const
292  {
293  return i + _nxi * j + _nxi * _neta * k;
294  }
295  int getNodeIndex(int i, int j, int k) const
296  {
297  return i + (_nxi + 1) * j + (_nxi + 1) * (_neta + 1) * k;
298  }
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  }
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  }
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  }
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
cartesianBox::_dxi
double _dxi
Definition: cartesian.h:38
cartesianBox::activeCellsEnd
cellIter activeCellsEnd()
Definition: cartesian.h:152
D
#define D
Definition: DefaultOptions.h:24
MElement::xyz2uvw
virtual void xyz2uvw(double xyz[3], double uvw[3]) const
Definition: MElement.cpp:1128
cartesianBox::_childBox
cartesianBox< scalar > * _childBox
Definition: cartesian.h:53
cartesianBox::insertActiveCell
void insertActiveCell(int t)
Definition: cartesian.h:285
dot
double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b)
Definition: STensor3.h:71
cartesianBox::getLC
double getLC()
Definition: cartesian.h:144
cartesianBox::_dzeta
double _dzeta
Definition: cartesian.h:38
Msg::Info
static void Info(const char *fmt,...)
Definition: GmshMessage.cpp:587
cartesianBox::getValueContainingPoint
double getValueContainingPoint(double x, double y, double z)
Definition: cartesian.h:174
cartesianBox::_etaAxis
SVector3 _etaAxis
Definition: cartesian.h:40
cartesianBox::createNodalValues
void createNodalValues()
Definition: cartesian.h:319
cartesianBox::_activeCells
std::set< int > _activeCells
Definition: cartesian.h:43
MVertex
Definition: MVertex.h:24
cartesianBox::getNodeTag
int getNodeTag(int index)
Definition: cartesian.h:299
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
cartesianBox::writeMSH
void writeMSH(const std::string &fileName, bool simplex=false, bool writeNodalValues=true)
Definition: cartesian.h:351
SPoint3
Definition: SPoint3.h:14
LegendrePolynomials::f
void f(int n, double u, double *val)
Definition: orthogonalBasis.cpp:77
cartesianBox::_xiAxis
SVector3 _xiAxis
Definition: cartesian.h:40
SVector3
Definition: SVector3.h:16
SVector3.h
cartesianBox::getLevel
int getLevel()
Definition: cartesian.h:149
cartesianBox::eraseActiveCell
void eraseActiveCell(int t)
Definition: cartesian.h:286
GmshMessage.h
cartesianBox::renumberNodes
void renumberNodes(int startingNum=1, cartesianBox< scalar > *parent=0)
Definition: cartesian.h:333
cartesianBox::getChildBox
cartesianBox< scalar > * getChildBox()
Definition: cartesian.h:148
cartesianBox::getNeta
int getNeta()
Definition: cartesian.h:146
cartesianBox::valIter
std::map< int, std::pair< scalar, int > >::iterator valIter
Definition: cartesian.h:153
cartesianBox::_printValues
void _printValues(FILE *f)
Definition: cartesian.h:120
cartesianBox::getNzeta
int getNzeta()
Definition: cartesian.h:147
SPoint3::x
double x(void) const
Definition: SPoint3.h:125
cartesianBox::_printElements
void _printElements(FILE *f, bool simplex, int startingNum=1)
Definition: cartesian.h:79
MHexahedron.h
MVertex.h
norm
void norm(const double *vec, double *norm)
Definition: gmshLevelset.cpp:202
cartesianBox::_z0
double _z0
Definition: cartesian.h:38
cartesianBox::setNodalValue
void setNodalValue(int i, scalar s)
Definition: cartesian.h:156
cartesianBox::getNodeCoordinates
SPoint3 getNodeCoordinates(int t) const
Definition: cartesian.h:275
SPoint3::y
double y(void) const
Definition: SPoint3.h:127
cartesianBox::getCellContainingPoint
int getCellContainingPoint(double x, double y, double z) const
Definition: cartesian.h:256
MHexahedron
Definition: MHexahedron.h:28
cartesianBox::_neta
int _neta
Definition: cartesian.h:36
cartesianBox::_deta
double _deta
Definition: cartesian.h:38
cartesianBox::_getNumElements
int _getNumElements(bool simplex)
Definition: cartesian.h:72
cartesianBox::getNodalValuesForCell
void getNodalValuesForCell(int t, std::vector< scalar > &values)
Definition: cartesian.h:157
cartesianBox::_level
int _level
Definition: cartesian.h:51
cartesianBox::getNodeIndex
int getNodeIndex(int i, int j, int k) const
Definition: cartesian.h:295
cartesianBox::_nodalValues
std::map< int, std::pair< scalar, int > > _nodalValues
Definition: cartesian.h:48
cartesianBox::_zetaAxis
SVector3 _zetaAxis
Definition: cartesian.h:40
cartesianBox::cellIter
std::set< int >::const_iterator cellIter
Definition: cartesian.h:150
cartesianBox
Definition: cartesian.h:33
cartesianBox::cartesianBox
cartesianBox(double X0, double Y0, double Z0, const SVector3 &dxi, const SVector3 &deta, const SVector3 &dzeta, int Nxi, int Neta, int Nzeta, int level=1)
Definition: cartesian.h:130
cartesianBox::_getNumNodes
int _getNumNodes()
Definition: cartesian.h:54
cartesianBox::nodalValuesEnd
valIter nodalValuesEnd()
Definition: cartesian.h:155
z
const double z
Definition: GaussQuadratureQuad.cpp:56
cartesianBox::nodalValuesBegin
valIter nodalValuesBegin()
Definition: cartesian.h:154
cartesianBox::getCellIndex
int getCellIndex(int i, int j, int k) const
Definition: cartesian.h:291
cartesianBox::getNodeIJK
void getNodeIJK(int index, int &i, int &j, int &k) const
Definition: cartesian.h:313
cartesianBox::activeCellExists
bool activeCellExists(int t)
Definition: cartesian.h:287
SPoint3::z
double z(void) const
Definition: SPoint3.h:129
cartesianBox::_nzeta
int _nzeta
Definition: cartesian.h:36
cartesianBox::_x0
double _x0
Definition: cartesian.h:38
cartesianBox::_y0
double _y0
Definition: cartesian.h:38
cartesianBox::activeCellsBegin
cellIter activeCellsBegin()
Definition: cartesian.h:151
SPoint3.h
cartesianBox::_printNodes
void _printNodes(FILE *f)
Definition: cartesian.h:62
MElement::interpolate
double interpolate(double val[], double u, double v, double w, int stride=1, int order=-1)
Definition: MElement.cpp:1216
cartesianBox::_nxi
int _nxi
Definition: cartesian.h:36
cartesianBox::getNxi
int getNxi()
Definition: cartesian.h:145
cartesianBox::getCellIJK
void getCellIJK(int index, int &i, int &j, int &k) const
Definition: cartesian.h:307
SVector3::normalize
double normalize()
Definition: SVector3.h:38