1 // @HEADER 2 // ************************************************************************ 3 // 4 // Rapid Optimization Library (ROL) Package 5 // Copyright (2014) Sandia Corporation 6 // 7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 8 // license for use of this work by or on behalf of the U.S. Government. 9 // 10 // Redistribution and use in source and binary forms, with or without 11 // modification, are permitted provided that the following conditions are 12 // met: 13 // 14 // 1. Redistributions of source code must retain the above copyright 15 // notice, this list of conditions and the following disclaimer. 16 // 17 // 2. Redistributions in binary form must reproduce the above copyright 18 // notice, this list of conditions and the following disclaimer in the 19 // documentation and/or other materials provided with the distribution. 20 // 21 // 3. Neither the name of the Corporation nor the names of the 22 // contributors may be used to endorse or promote products derived from 23 // this software without specific prior written permission. 24 // 25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 36 // 37 // Questions? Contact lead developers: 38 // Drew Kouri (dpkouri@sandia.gov) and 39 // Denis Ridzal (dridzal@sandia.gov) 40 // 41 // ************************************************************************ 42 // @HEADER 43 44 /*! \file DofManager.hpp 45 \brief Sets up data structures for the indexing of the global degrees 46 of freedom (Dof's). 47 */ 48 49 #include "Teuchos_ParameterList.hpp" 50 #include "Intrepid_FieldContainer.hpp" 51 52 template <class Real> 53 class DofManager { 54 55 typedef Intrepid::FieldContainer<Real> FCR; 56 typedef Intrepid::FieldContainer<int> FCI; 57 58 enum FieldType { 59 VELX, VELY, PRES 60 }; 61 62 /* Backward-facing step geometry. 63 64 *************************************************** 65 * * * * 66 * 3 * 4 * 5 * 67 * * * * 68 ********** * * * * * * * * * * * * * * * * * * * ** 69 * 1 * 2 * 70 * * * 71 ****************************************** 72 */ 73 74 private: 75 Real channelH_; // channel height (height of regions 1+4) 76 Real channelW_; // channel width (width of regions 3+4+5) 77 Real stepH_; // step height (height of region 1) 78 Real stepW_; // step width (width of region 3) 79 Real observeW_; // width of observation region (width of region 1) 80 81 int nx1_; 82 int nx2_; 83 int nx3_; 84 int nx4_; 85 int nx5_; 86 int ny1_; 87 int ny2_; 88 int ny3_; 89 int ny4_; 90 int ny5_; 91 92 int ref_; 93 94 int numCells_; 95 int numNodes_; 96 int numEdges_; 97 int numMeshEntities_; 98 99 FCR meshNodes_; 100 FCI meshCellToNodeMap_; 101 FCI meshCellToEdgeMap_; 102 103 int numDofs_; 104 int numPres_; 105 int numVelX_; 106 int numVelY_; 107 int numLocalDofs_; 108 int numLocalPres_; 109 int numLocalVelX_; 110 int numLocalVelY_; 111 112 FCI dofArray_; 113 FCI cellDofs_; 114 FCI cellPres_; 115 FCI cellVelX_; 116 FCI cellVelY_; 117 118 public: 119 DofManager(Teuchos::ParameterList & parlist)120 DofManager(Teuchos::ParameterList &parlist) { 121 // Geometry data. 122 channelH_ = parlist.sublist("Navier Stokes").get( "Channel height", 1.0); 123 channelW_ = parlist.sublist("Navier Stokes").get( "Channel width", 8.0); 124 stepH_ = parlist.sublist("Navier Stokes").get( "Step height", 0.5); 125 stepW_ = parlist.sublist("Navier Stokes").get( "Step width", 1.0); 126 observeW_ = parlist.sublist("Navier Stokes").get("Observation width", 3.0); 127 // Mesh data. 128 ref_ = parlist.sublist("Navier Stokes").get( "Refinement level", 1); 129 nx1_ = parlist.sublist("Navier Stokes").get( "Observation region NX", 4*ref_); 130 ny1_ = parlist.sublist("Navier Stokes").get( "Observation region NY", 5*ref_); 131 nx2_ = parlist.sublist("Navier Stokes").get("Laminar flow region NX", 2*ref_); 132 ny2_ = ny1_; 133 nx3_ = parlist.sublist("Navier Stokes").get( "Inflow region NX", 1*ref_); 134 ny3_ = parlist.sublist("Navier Stokes").get( "Inflow region NY", 2*ref_); 135 nx4_ = nx1_; 136 ny4_ = ny3_; 137 nx5_ = nx2_; 138 ny5_ = ny3_; 139 numCells_ = (nx1_ + nx2_)*ny1_ + (nx3_ + nx1_ + nx2_)*ny3_; 140 numNodes_ = (nx1_ + nx2_ + 1)*ny1_ + (nx3_ + nx1_ + nx2_ + 1)*(ny3_ + 1); 141 numEdges_ = (2*(nx1_ + nx2_)+1)*ny1_ + (2*(nx3_ + nx1_ + nx2_)+1)*ny3_ + (nx3_ + nx1_ + nx2_); 142 numMeshEntities_ = numCells_ + numNodes_ + numEdges_; 143 // Mesh and Dof data structures. 144 numLocalPres_ = 4; 145 numLocalVelX_ = 9; 146 numLocalVelY_ = 9; 147 numLocalDofs_ = numLocalPres_ + numLocalVelX_ + numLocalVelY_; 148 computeMeshNodes(meshNodes_); 149 computeMeshCellToNodeMap(meshCellToNodeMap_); 150 computeMeshCellToEdgeMap(meshCellToEdgeMap_); 151 computeDofArray(dofArray_, numDofs_); 152 computeCellDofs(cellDofs_, dofArray_, meshCellToNodeMap_, meshCellToEdgeMap_); 153 computeCellPres(cellPres_, cellDofs_); 154 computeCellVelX(cellVelX_, cellDofs_); 155 computeCellVelY(cellVelY_, cellDofs_); 156 } 157 158 computeMeshNodes(FCR & nodes) const159 void computeMeshNodes(FCR &nodes) const { 160 161 Real dy1 = stepH_ / ny1_; 162 Real dy3 = (channelH_ - stepH_) / ny3_; 163 Real dx1 = observeW_ / nx1_; 164 Real dx2 = (channelW_ - stepW_ - observeW_) / nx2_; 165 Real dx3 = stepW_ / nx3_; 166 int nodeCt = 0; 167 nodes.resize(numNodes_, 2); 168 169 // bottom region 170 for (int j=0; j<ny1_; ++j) { 171 for (int i=0; i<=nx1_; ++i) { 172 nodes(nodeCt, 0) = stepW_ + i*dx1; 173 nodes(nodeCt, 1) = j*dy1; 174 ++nodeCt; 175 } 176 for (int i=0; i<nx2_; ++i) { 177 nodes(nodeCt, 0) = stepW_ + observeW_ + (i+1)*dx2; 178 nodes(nodeCt, 1) = j*dy1; 179 ++nodeCt; 180 } 181 } 182 183 // top region 184 for (int j=0; j<=ny3_; ++j) { 185 for (int i=0; i<=nx3_; ++i) { 186 nodes(nodeCt, 0) = i*dx3; 187 nodes(nodeCt, 1) = stepH_ + j*dy3; 188 ++nodeCt; 189 } 190 for (int i=0; i<nx1_; ++i) { 191 nodes(nodeCt, 0) = stepW_ + (i+1)*dx1; 192 nodes(nodeCt, 1) = stepH_ + j*dy3; 193 ++nodeCt; 194 } 195 for (int i=0; i<nx2_; ++i) { 196 nodes(nodeCt, 0) = stepW_ + observeW_ + (i+1)*dx2; 197 nodes(nodeCt, 1) = stepH_ + j*dy3; 198 ++nodeCt; 199 } 200 } 201 202 } // computeMeshNodes 203 204 computeMeshCellToNodeMap(FCI & ctn)205 void computeMeshCellToNodeMap(FCI &ctn) { 206 207 int cellCt = 0; 208 ctn.resize(numCells_, 4); 209 210 // bottom region 211 for (int j=0; j<ny1_-1; ++j) { 212 for (int i=0; i<nx1_+nx2_; ++i) { 213 ctn(cellCt, 0) = j*(nx1_+nx2_+1) + i; 214 ctn(cellCt, 1) = j*(nx1_+nx2_+1) + (i+1); 215 ctn(cellCt, 2) = (j+1)*(nx1_+nx2_+1) + (i+1); 216 ctn(cellCt, 3) = (j+1)*(nx1_+nx2_+1) + i; 217 ++cellCt; 218 } 219 } 220 221 // transition region 222 for (int i=0; i<nx1_+nx2_; ++i) { 223 ctn(cellCt, 0) = (ny1_-1)*(nx1_+nx2_+1) + i; 224 ctn(cellCt, 1) = (ny1_-1)*(nx1_+nx2_+1) + (i+1); 225 ctn(cellCt, 2) = ny1_*(nx1_+nx2_+1) + nx3_ + (i+1); 226 ctn(cellCt, 3) = ny1_*(nx1_+nx2_+1) + nx3_ + i; 227 ++cellCt; 228 } 229 230 // top region 231 for (int j=0; j<ny3_; ++j) { 232 for (int i=0; i<nx3_+nx1_+nx2_; ++i) { 233 ctn(cellCt, 0) = ny1_*(nx1_+nx2_+1) + j*(nx3_+nx1_+nx2_+1) + i; 234 ctn(cellCt, 1) = ny1_*(nx1_+nx2_+1) + j*(nx3_+nx1_+nx2_+1) + (i+1); 235 ctn(cellCt, 2) = ny1_*(nx1_+nx2_+1) + (j+1)*(nx3_+nx1_+nx2_+1) + (i+1); 236 ctn(cellCt, 3) = ny1_*(nx1_+nx2_+1) + (j+1)*(nx3_+nx1_+nx2_+1) + i; 237 ++cellCt; 238 } 239 } 240 241 } // computeMeshCellToNodeMap 242 243 computeMeshCellToEdgeMap(FCI & cte)244 void computeMeshCellToEdgeMap(FCI &cte) { 245 246 int cellCt = 0; 247 cte.resize(numCells_, 4); 248 249 // bottom region 250 for (int j=0; j<ny1_-1; ++j) { 251 for (int i=0; i<nx1_+nx2_; ++i) { 252 cte(cellCt, 0) = j*(2*(nx1_+nx2_)+1) + i; 253 cte(cellCt, 1) = j*(2*(nx1_+nx2_)+1) + (nx1_+nx2_) + (i+1); 254 cte(cellCt, 2) = (j+1)*(2*(nx1_+nx2_)+1) + i; 255 cte(cellCt, 3) = j*(2*(nx1_+nx2_)+1) + (nx1_+nx2_) + i; 256 ++cellCt; 257 } 258 } 259 260 // transition region 261 for (int i=0; i<nx1_+nx2_; ++i) { 262 cte(cellCt, 0) = (ny1_-1)*(2*(nx1_+nx2_)+1) + i; 263 cte(cellCt, 1) = (ny1_-1)*(2*(nx1_+nx2_)+1) + (nx1_+nx2_) + (i+1); 264 cte(cellCt, 2) = ny1_*(2*(nx1_+nx2_)+1) + nx3_ + i; 265 cte(cellCt, 3) = (ny1_-1)*(2*(nx1_+nx2_)+1) + (nx1_+nx2_) + i; 266 ++cellCt; 267 } 268 269 // top region 270 for (int j=0; j<ny3_; ++j) { 271 for (int i=0; i<nx3_+nx1_+nx2_; ++i) { 272 cte(cellCt, 0) = ny1_*(2*(nx1_+nx2_)+1) + j*(2*(nx3_+nx1_+nx2_)+1) + i; 273 cte(cellCt, 1) = ny1_*(2*(nx1_+nx2_)+1) + j*(2*(nx3_+nx1_+nx2_)+1) + (nx3_+nx1_+nx2_) + (i+1); 274 cte(cellCt, 2) = ny1_*(2*(nx1_+nx2_)+1) + (j+1)*(2*(nx3_+nx1_+nx2_)+1) + i; 275 cte(cellCt, 3) = ny1_*(2*(nx1_+nx2_)+1) + j*(2*(nx3_+nx1_+nx2_)+1) + (nx3_+nx1_+nx2_) + i; 276 ++cellCt; 277 } 278 } 279 280 } // computeMeshCellToEdgeMap 281 282 computeDofArray(FCI & dof,int & numDofs,const int & globOrder=1)283 void computeDofArray(FCI &dof, int &numDofs, const int &globOrder=1) { 284 285 dof.resize(numMeshEntities_, 3); 286 int dofCt = -1; 287 288 if (globOrder == 1) { 289 // 290 // This is the independent node id --> edge id --> cell id ordering: 291 // 292 // VelX VelY Pres 293 // -------------- 294 // n1 1 1 1 295 // n2 1 1 1 296 // ... 297 // e1 1 1 0 298 // e2 1 1 0 299 // ... 300 // c1 1 1 0 301 // c2 1 1 0 302 // ... 303 // 304 for (int i=0; i<numNodes_; ++i) { 305 // dof(i, 1) --> 1 Velocity X 306 // dof(i, 1) --> 1 Velocity Y 307 // dof(i, 2) --> 1 Pressure 308 dof(i, 0) = ++dofCt; // VelX 309 dof(i, 1) = ++dofCt; // VelY 310 dof(i, 2) = ++dofCt; // Pres 311 } 312 for (int i=numNodes_; i<numNodes_+numEdges_; ++i) { 313 // dof(i, 1) --> 1 Velocity X 314 // dof(i, 1) --> 1 Velocity Y 315 // dof(i, 2) --> 0 Pressure 316 dof(i, 0) = ++dofCt; // VelX 317 dof(i, 1) = ++dofCt; // VelY 318 dof(i, 2) = -1; 319 } 320 for (int i=numNodes_+numEdges_; i<numMeshEntities_; ++i) { 321 // dof(i, 1) --> 1 Velocity X 322 // dof(i, 1) --> 1 Velocity Y 323 // dof(i, 2) --> 0 Pressure 324 dof(i, 0) = ++dofCt; // VelX 325 dof(i, 1) = ++dofCt; // VelY 326 dof(i, 2) = -1; 327 } 328 } 329 else if (globOrder==2) { 330 // 331 // This is the cell-driven cell node id --> cell edge id --> cell id ordering: 332 // 333 // VelX VelY Pres 334 // -------------- 335 // cell 1 336 // n1 1 1 1 337 // n2 1 1 1 338 // ... 339 // e1 1 1 0 340 // e2 1 1 0 341 // ... 342 // c1 1 1 0 343 // c2 1 1 0 344 // ... 345 // cell 2 346 // ... 347 // (unique ids only, of course) 348 // 349 // TO BE IMPLEMENTED!!! 350 } 351 numDofs = dofCt + 1; 352 } // computeDofArray 353 354 computeCellDofs(FCI & cdofs,const FCI & dofs,const FCI & ctn,const FCI & cte)355 void computeCellDofs(FCI &cdofs, const FCI &dofs, const FCI &ctn, const FCI &cte) { 356 357 cdofs.resize(numCells_, numLocalDofs_); 358 359 for (int i=0; i<numCells_; ++i) { 360 int ct = -1; 361 int gid0 = ctn(i,0); 362 int gid1 = ctn(i,1); 363 int gid2 = ctn(i,2); 364 int gid3 = ctn(i,3); 365 int gid4 = numNodes_ + cte(i,0); 366 int gid5 = numNodes_ + cte(i,1); 367 int gid6 = numNodes_ + cte(i,2); 368 int gid7 = numNodes_ + cte(i,3); 369 int gid8 = numNodes_ + numEdges_ + i; 370 cdofs(i,++ct) = dofs(gid0, 0); 371 cdofs(i,++ct) = dofs(gid0, 1); 372 cdofs(i,++ct) = dofs(gid0, 2); 373 cdofs(i,++ct) = dofs(gid1, 0); 374 cdofs(i,++ct) = dofs(gid1, 1); 375 cdofs(i,++ct) = dofs(gid1, 2); 376 cdofs(i,++ct) = dofs(gid2, 0); 377 cdofs(i,++ct) = dofs(gid2, 1); 378 cdofs(i,++ct) = dofs(gid2, 2); 379 cdofs(i,++ct) = dofs(gid3, 0); 380 cdofs(i,++ct) = dofs(gid3, 1); 381 cdofs(i,++ct) = dofs(gid3, 2); 382 cdofs(i,++ct) = dofs(gid4, 0); 383 cdofs(i,++ct) = dofs(gid4, 1); 384 cdofs(i,++ct) = dofs(gid5, 0); 385 cdofs(i,++ct) = dofs(gid5, 1); 386 cdofs(i,++ct) = dofs(gid6, 0); 387 cdofs(i,++ct) = dofs(gid6, 1); 388 cdofs(i,++ct) = dofs(gid7, 0); 389 cdofs(i,++ct) = dofs(gid7, 1); 390 cdofs(i,++ct) = dofs(gid8, 0); 391 cdofs(i,++ct) = dofs(gid8, 1); 392 } 393 } // computeCellDofs 394 395 computeCellPres(FCI & cpres,const FCI & cdofs)396 void computeCellPres(FCI &cpres, const FCI &cdofs) { 397 cpres.resize(numCells_, numLocalPres_); 398 for (int i=0; i<numCells_; ++i) { 399 for (int j=0; j<numLocalPres_; ++j) { 400 cpres(i,j) = cdofs(i, j*3+2); 401 } 402 } 403 } // computeCellPres 404 405 computeCellVelX(FCI & cvelx,const FCI & cdofs)406 void computeCellVelX(FCI &cvelx, const FCI &cdofs) { 407 cvelx.resize(numCells_, numLocalVelX_); 408 for (int i=0; i<numCells_; ++i) { 409 for (int j=0; j<numLocalPres_; ++j) { 410 cvelx(i,j) = cdofs(i, j*3); 411 } 412 for (int k=0; k<numLocalVelX_-numLocalPres_; ++k) { 413 cvelx(i,k+numLocalPres_) = cdofs(i, 4*3 + k*2); 414 } 415 } 416 } // computeCellVelX 417 418 computeCellVelY(FCI & cvely,const FCI & cdofs)419 void computeCellVelY(FCI &cvely, const FCI &cdofs) { 420 cvely.resize(numCells_, numLocalVelY_); 421 for (int i=0; i<numCells_; ++i) { 422 for (int j=0; j<numLocalPres_; ++j) { 423 cvely(i,j) = cdofs(i, j*3+1); 424 } 425 for (int k=0; k<numLocalVelY_-numLocalPres_; ++k) { 426 cvely(i,k+numLocalPres_) = cdofs(i, 4*3 + k*2 + 1); 427 } 428 } 429 } // computeCellVelY 430 431 getCellDofs(const int & cid,const FieldType & ft) const432 std::vector<int>& getCellDofs(const int& cid, const FieldType& ft) const { 433 std::vector<int> dofVec; 434 int nPres = 4; 435 int nVelX = 9; 436 int nVelY = 9; 437 switch(ft) { 438 case VELX: 439 break; 440 441 case VELY: 442 break; 443 444 case PRES: 445 for (int i=0; i<nPres; ++i) { 446 dofVec.push_back(cellDofs_(cid)); 447 } 448 break; 449 } 450 return dofVec; 451 } // getDofArray 452 453 setRefinementLevel(const int & refLevel)454 void setRefinementLevel(const int &refLevel) { 455 ref_ = refLevel; 456 } // setRefinementLevel 457 458 getNumCells() const459 int getNumCells() const { 460 return numCells_; 461 } // getNumCells 462 463 getNumNodes() const464 int getNumNodes() const { 465 return numNodes_; 466 } // getNumNodes 467 468 getNumEdges() const469 int getNumEdges() const { 470 return numEdges_; 471 } // getNumEdges 472 473 getMeshNodes() const474 const FCR& getMeshNodes() const { 475 return meshNodes_; 476 } // getMeshNodes 477 478 getMeshCellToNodeMap() const479 const FCI& getMeshCellToNodeMap() const { 480 return meshCellToNodeMap_; 481 } // getMeshCellToNodeMap 482 483 getMeshCellToEdgeMap() const484 const FCI& getMeshCellToEdgeMap() const { 485 return meshCellToEdgeMap_; 486 } // getMeshCellToEdgeMap 487 488 getNumDofs() const489 int getNumDofs() const { 490 return numDofs_; 491 } // getNumDofs 492 493 getDofArray() const494 const FCI& getDofArray() const { 495 return dofArray_; 496 } // getDofArray 497 498 getCellDofs() const499 const FCI& getCellDofs() const { 500 return cellDofs_; 501 } // getCellDofs 502 503 getCellPres() const504 const FCI& getCellPres() const { 505 return cellPres_; 506 } // getCellPres 507 508 getCellVelX() const509 const FCI& getCellVelX() const { 510 return cellVelX_; 511 } // getCellVelX 512 513 getCellVelY() const514 const FCI& getCellVelY() const { 515 return cellVelY_; 516 } // getCellVelY 517 518 }; 519