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