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