Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ProcessBL.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessBL.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Refine prismatic or quadrilateral boundary layer elements.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <string>
37 
43 #include <LocalRegions/HexExp.h>
44 #include <LocalRegions/PrismExp.h>
45 #include <LocalRegions/QuadExp.h>
46 
48 #include "ProcessBL.h"
49 
50 using namespace std;
51 using namespace Nektar::NekMeshUtils;
52 
53 namespace Nektar
54 {
55 namespace Utilities
56 {
57 ModuleKey ProcessBL::className = GetModuleFactory().RegisterCreatorFunction(
59  ProcessBL::create,
60  "Refines a prismatic or quadrilateral boundary layer.");
61 
62 int **helper2d(int lda, int arr[][2])
63 {
64  int **ret = new int *[lda];
65  for (int i = 0; i < lda; ++i)
66  {
67  ret[i] = new int[2];
68  ret[i][0] = arr[i][0];
69  ret[i][1] = arr[i][1];
70  }
71  return ret;
72 }
73 
74 int **helper2d(int lda, int arr[][4])
75 {
76  int **ret = new int *[lda];
77  for (int i = 0; i < lda; ++i)
78  {
79  ret[i] = new int[4];
80  ret[i][0] = arr[i][0];
81  ret[i][1] = arr[i][1];
82  ret[i][2] = arr[i][2];
83  ret[i][3] = arr[i][3];
84  }
85  return ret;
86 }
87 
89 {
90  int size;
91  int layerOff;
92  int *edge;
93  int *offset;
94  int *inc;
95  int **conn;
97  int *bfaces;
98 };
99 
101 {
102  int size;
103  int *edge;
104  int **edgeVert;
105  int *offset;
106  int *inc;
107 };
108 
109 ProcessBL::ProcessBL(MeshSharedPtr m) : ProcessModule(m)
110 {
111  // BL mesh configuration.
112  m_config["layers"] =
113  ConfigOption(false, "2", "Number of layers to refine.");
114  m_config["nq"] =
115  ConfigOption(false, "5", "Number of points in high order elements.");
116  m_config["surf"] =
117  ConfigOption(false, "", "Tag identifying surface connected to prism.");
118  m_config["r"] =
119  ConfigOption(false, "2.0", "Ratio to use in geometry progression.");
120 }
121 
123 {
124 }
125 
127 {
128  if (m_mesh->m_verbose)
129  {
130  cout << "ProcessBL: Refining boundary layer..." << endl;
131  }
132  int dim = m_mesh->m_expDim;
133  switch (dim)
134  {
135  case 2:
136  {
137  BoundaryLayer2D();
138  }
139  break;
140 
141  case 3:
142  {
143  BoundaryLayer3D();
144  }
145  break;
146 
147  default:
148  ASSERTL0(0,"Dimension not supported");
149  break;
150  }
151 
152  // Re-process mesh to eliminate duplicate vertices and edges.
153  ProcessVertices();
154  ProcessEdges();
155  ProcessFaces();
156  ProcessElements();
158 }
159 
161 {
162  int nodeId = m_mesh->m_vertexSet.size();
163  int nl = m_config["layers"].as<int>();
164  int nq = m_config["nq"]. as<int>();
165 
166  // determine if geometric ratio is string or a constant.
168  NekDouble r = 1;
169  int rExprId = -1;
170  bool ratioIsString = false;
171 
172  if (m_config["r"].isType<NekDouble>())
173  {
174  r = m_config["r"].as<NekDouble>();
175  }
176  else
177  {
178  std::string rstr = m_config["r"].as<string>();
179  rExprId = rEval.DefineFunction("x y z", rstr);
180  ratioIsString = true;
181  }
182 
183  // Default PointsType.
185 
186  // Map which takes element ID to edge on surface. This enables
187  // splitting to occur in either y-direction of the prism.
188  map<int, int> splitEls;
189 
190  // edgeMap associates geometry edge IDs to the (nl+1) vertices which
191  // are generated along that edge when a prism is split, and is used
192  // to avoid generation of duplicate vertices. It is stored as an
193  // unordered map for speed.
194  boost::unordered_map<int, vector<NodeSharedPtr> > edgeMap;
195  boost::unordered_map<int, vector<NodeSharedPtr> >::iterator eIt;
196 
197  string surf = m_config["surf"].as<string>();
198  if (surf.size() > 0)
199  {
200  vector<unsigned int> surfs;
201  ParseUtils::GenerateSeqVector(surf.c_str(), surfs);
202  sort(surfs.begin(), surfs.end());
203 
204  // If surface is defined, process list of elements to find those
205  // that are connected to it.
206  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
207  {
208  ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
209  int nSurf = el->GetEdgeCount();
210 
211  for (int j = 0; j < nSurf; ++j)
212  {
213  int bl = el->GetBoundaryLink(j);
214  if (bl == -1)
215  {
216  continue;
217  }
218 
219  ElementSharedPtr bEl = m_mesh->m_element[m_mesh->m_expDim-1][bl];
220  vector<int> tags = bEl->GetTagList();
221  vector<int> inter;
222 
223  sort(tags.begin(), tags.end());
224  set_intersection(surfs.begin(), surfs.end(),
225  tags .begin(), tags .end(),
226  back_inserter(inter));
227  ASSERTL0(inter.size() <= 1,
228  "Intersection of surfaces wrong");
229 
230  if (inter.size() == 1)
231  {
232  if (el->GetConf().m_e != LibUtilities::eQuadrilateral)
233  {
234  cerr << "WARNING: Found non-quad element "
235  << "to split in surface " << surf
236  << "; ignoring" << endl;
237  continue;
238  }
239 
240  if (splitEls.count(el->GetId()) > 0)
241  {
242  cerr << "WARNING: quad already found; "
243  << "ignoring" << endl;
244  continue;
245  }
246 
247  splitEls[el->GetId()] = j;
248  }
249  }
250  }
251  }
252  else
253  {
254  ASSERTL0(false, "Surface must be specified.");
255  }
256 
257  if (splitEls.size() == 0)
258  {
259  cerr << "WARNING: No elements detected to split." << endl;
260  return;
261  }
262 
263  // Erase all elements from the element list. Elements will be
264  // re-added as they are split.
265  vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim];
266  m_mesh->m_element[m_mesh->m_expDim].clear();
267 
268  // Iterate over list of elements of expansion dimension.
269  for (int i = 0; i < el.size(); ++i)
270  {
271 
272 
273  if (splitEls.count(el[i]->GetId()) == 0)
274  {
275  m_mesh->m_element[m_mesh->m_expDim].push_back(el[i]);
276  continue;
277  }
278 
279  // Find other boundary faces if any
280  std::map<int, int> bLink;
281  for (int j = 0; j < 4; j++)
282  {
283  int bl = el[i]->GetBoundaryLink(j);
284  if ( (bl != -1) && ( j != splitEls[el[i]->GetId()]) )
285  {
286  bLink[j] = bl;
287  }
288  }
289 
290  // Get elemental geometry object.
292  boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
293  el[i]->GetGeom(m_mesh->m_spaceDim));
294 
295  // Determine whether to use reverse points.
296  // (if edges 1 or 2 are on the surface)
297  LibUtilities::PointsType t = ( (splitEls[el[i]->GetId()]+1) %4) < 2 ?
300 
301 
302  if(ratioIsString) // determine value of r base on geom
303  {
304  NekDouble x,y,z;
305  NekDouble x1,y1,z1;
306  int nverts = geom->GetNumVerts();
307 
308  x = y = z = 0.0;
309 
310  for(int i = 0; i < nverts; ++i)
311  {
312  geom->GetVertex(i)->GetCoords(x1,y1,z1);
313  x += x1; y += y1; z += z1;
314  }
315  x /= (NekDouble) nverts;
316  y /= (NekDouble) nverts;
317  z /= (NekDouble) nverts;
318  r = rEval.Evaluate(rExprId,x,y,z,0.0);
319  }
320 
321  // Create basis.
324  LibUtilities::PointsKey(nq,pt));
327  LibUtilities::PointsKey(nl+1, t, r));
328 
329  // Create local region.
331  if (splitEls[el[i]->GetId()] % 2)
332  {
334  B1,B0,geom);
335  }
336  else
337  {
339  B0,B1,geom);
340  }
341 
342  // Grab co-ordinates.
343  Array<OneD, NekDouble> x(nq*(nl+1));
344  Array<OneD, NekDouble> y(nq*(nl+1));
345  Array<OneD, NekDouble> z(nq*(nl+1));
346  q->GetCoords(x,y,z);
347 
348  vector<vector<NodeSharedPtr> > edgeNodes(2);
349 
350  // Loop over edges to be split.
351  for (int j = 0; j < 2; ++j)
352  {
353  int locEdge = (splitEls[el[i]->GetId()]+1+2*j)%4;
354  int edgeId = el[i]->GetEdge(locEdge)->m_id;
355 
356  // Determine whether we have already generated vertices
357  // along this edge.
358  eIt = edgeMap.find(edgeId);
359 
360  if (eIt == edgeMap.end())
361  {
362  // If not then resize storage to hold new points.
363  edgeNodes[j].resize(nl+1);
364 
365  // Re-use existing vertices at endpoints of edge to
366  // avoid duplicating the existing vertices.
367  edgeNodes[j][0] = el[i]->GetVertex(locEdge);
368  edgeNodes[j][nl] = el[i]->GetVertex((locEdge+1)%4);
369 
370  // Variable geometric ratio
371  if(ratioIsString)
372  {
373  NekDouble x0,y0;
374  NekDouble x1,y1;
375  NekDouble xm,ym,zm=0.0;
376 
377  // -> Find edge end and mid points
378  x0 = edgeNodes[j][0]->m_x;
379  y0 = edgeNodes[j][0]->m_y;
380 
381  x1 = edgeNodes[j][nl]->m_x;
382  y1 = edgeNodes[j][nl]->m_y;
383 
384  xm = 0.5*(x0+x1);
385  ym = 0.5*(y0+y1);
386 
387  // evaluate r factor based on mid point value
388  NekDouble rnew;
389  rnew = rEval.Evaluate(rExprId,xm,ym,zm,0.0);
390 
391  // Get basis with new r;
394  LibUtilities::PointsKey Pkey(nl+1, t, rnew);
396  = LibUtilities::PointsManager()[Pkey];
397 
398  const Array<OneD, const NekDouble> z = newP->GetZ();
399 
400  // Create new interior nodes based on this new blend
401  for (int k = 1; k < nl; ++k)
402  {
403  xm = 0.5*(1+z[k])*(x1-x0) + x0;
404  ym = 0.5*(1+z[k])*(y1-y0) + y0;
405  zm = 0.0;
406  edgeNodes[j][k] = NodeSharedPtr(
407  new Node(nodeId++, xm,ym,zm));
408  }
409  }
410  else
411  {
412  // Create new interior nodes.
413  int pos;
414  for (int k = 1; k < nl; ++k)
415  {
416  switch (locEdge)
417  {
418  case 0:
419  pos = k;
420  break;
421  case 1:
422  pos = nq -1 + k*nq;
423  break;
424  case 2:
425  pos = nq*(nl+1) -1 - k;
426  break;
427  case 3:
428  pos = nq*nl - k*nq;
429  break;
430  default:
431  ASSERTL0(0,"Quad edge should be < 4.");
432  break;
433  }
434  edgeNodes[j][k] = NodeSharedPtr(
435  new Node(nodeId++, x[pos], y[pos], z[pos]));
436  }
437  }
438 
439  // Store these edges in edgeMap.
440  edgeMap[edgeId] = edgeNodes[j];
441  }
442  else
443  {
444  // Check orientation
445  if (eIt->second[0] == el[i]->GetVertex(locEdge))
446  {
447  // Same orientation: copy nodes
448  edgeNodes[j] = eIt->second;
449  }
450  else
451  {
452  // Reversed orientation: copy in reversed order
453  edgeNodes[j].resize(nl+1);
454  for (int k = 0; k < nl+1; ++k)
455  {
456  edgeNodes[j][k] = eIt->second[nl-k];
457  }
458  }
459 
460  }
461  }
462 
463  // Create element layers.
464  for (int j = 0; j < nl; ++j)
465  {
466  // Get corner vertices.
467  vector<NodeSharedPtr> nodeList(4);
468  switch (splitEls[el[i]->GetId()])
469  {
470  case 0:
471  {
472  nodeList[0] = edgeNodes[1][nl-j ];
473  nodeList[1] = edgeNodes[0][j ];
474  nodeList[2] = edgeNodes[0][j+1];
475  nodeList[3] = edgeNodes[1][nl-j-1];
476  break;
477  }
478  case 1:
479  {
480  nodeList[0] = edgeNodes[1][j];
481  nodeList[1] = edgeNodes[1][j+1];
482  nodeList[2] = edgeNodes[0][nl-j-1];
483  nodeList[3] = edgeNodes[0][nl-j];
484  break;
485  }
486  case 2:
487  {
488  nodeList[0] = edgeNodes[0][nl-j ];
489  nodeList[1] = edgeNodes[1][j ];
490  nodeList[2] = edgeNodes[1][j+1];
491  nodeList[3] = edgeNodes[0][nl-j-1];
492  break;
493  }
494  case 3:
495  {
496  nodeList[0] = edgeNodes[0][j];
497  nodeList[1] = edgeNodes[0][j+1];
498  nodeList[2] = edgeNodes[1][nl-j-1];
499  nodeList[3] = edgeNodes[1][nl-j];
500  break;
501  }
502  }
503  // Create the element.
504  ElmtConfig conf(LibUtilities::eQuadrilateral, 1, true, false, true);
506  CreateInstance(
507  LibUtilities::eQuadrilateral,conf,nodeList,el[i]->GetTagList());
508 
509  // Add high order nodes to split edges.
510  for (int l = 0; l < 2; ++l)
511  {
512  int locEdge = (splitEls[el[i]->GetId()]+2*l)%4;
513  EdgeSharedPtr HOedge = elmt->GetEdge(
514  locEdge);
515  int pos;
516  for (int k = 1; k < nq-1; ++k)
517  {
518  switch (locEdge)
519  {
520  case 0:
521  pos = j*nq+ k;
522  break;
523  case 1:
524  pos = j+1 + k*(nl+1);
525  break;
526  case 2:
527  pos = (j+1)*nq + (nq-1) - k;
528  break;
529  case 3:
530  pos = (nl+1)*(nq-1) + j - k*(nl+1);
531  break;
532  default:
533  ASSERTL0(0,"Quad edge should be < 4.");
534  break;
535  }
536  HOedge->m_edgeNodes.push_back(
538  new Node(nodeId++,x[pos],y[pos],0.0)));
539  }
540  HOedge->m_curveType = pt;
541  }
542 
543  // Change the elements on the boundary
544  // to match the layers
546  for (it = bLink.begin(); it != bLink.end(); ++it)
547  {
548  int eid = it->first;
549  int bl = it->second;
550 
551  if (j == 0)
552  {
553  // For first layer reuse existing 2D element.
554  ElementSharedPtr e = m_mesh->m_element[m_mesh->m_expDim-1][bl];
555  for (int k = 0; k < 2; ++k)
556  {
557  e->SetVertex(
558  k, nodeList[(eid+k)%4]);
559  }
560  }
561  else
562  {
563  // For all other layers create new element.
564  vector<NodeSharedPtr> qNodeList(2);
565  for (int k = 0; k < 2; ++k)
566  {
567  qNodeList[k] = nodeList[(eid+k)%4];
568  }
569  vector<int> tagBE;
570  tagBE = m_mesh->m_element[m_mesh->m_expDim-1][bl]->GetTagList();
571  ElmtConfig bconf(LibUtilities::eSegment,1,true,true,false);
572  ElementSharedPtr boundaryElmt = GetElementFactory().
573  CreateInstance(LibUtilities::eSegment,bconf,
574  qNodeList,tagBE);
575  m_mesh->m_element[m_mesh->m_expDim-1].push_back(boundaryElmt);
576  }
577  }
578 
579  m_mesh->m_element[m_mesh->m_expDim].push_back(elmt);
580  }
581  }
582 }
583 
585 {
586  // A set containing all element types which are valid.
587  set<LibUtilities::ShapeType> validElTypes;
588  validElTypes.insert(LibUtilities::ePrism);
589  validElTypes.insert(LibUtilities::eHexahedron);
590 
591  int nodeId = m_mesh->m_vertexSet.size();
592  int nl = m_config["layers"].as<int>();
593  int nq = m_config["nq"].as<int>();
594 
595  // determine if geometric ratio is string or a constant.
597  NekDouble r = 1;
598  int rExprId = -1;
599  bool ratioIsString = false;
600 
601  if (m_config["r"].isType<NekDouble>())
602  {
603  r = m_config["r"].as<NekDouble>();
604  }
605  else
606  {
607  std::string rstr = m_config["r"].as<string>();
608  rExprId = rEval.DefineFunction("x y z", rstr);
609  ratioIsString = true;
610  }
611 
612  // Prismatic node -> face map.
613  int prismFaceNodes[5][4] = {
614  {0, 1, 2, 3}, {0, 1, 4, -1}, {1, 2, 5, 4}, {3, 2, 5, -1}, {0, 3, 5, 4}};
615  int hexFaceNodes[6][4] = {{0, 1, 2, 3},
616  {0, 1, 5, 4},
617  {1, 2, 6, 5},
618  {3, 2, 6, 7},
619  {0, 3, 7, 4},
620  {4, 5, 6, 7}};
621  map<LibUtilities::ShapeType, int **> faceNodeMap;
622  faceNodeMap[LibUtilities::ePrism] = helper2d(5, prismFaceNodes);
623  faceNodeMap[LibUtilities::eHexahedron] = helper2d(6, hexFaceNodes);
624 
625  // Default PointsType.
627 
628  // Map which takes element ID to face on surface. This enables
629  // splitting to occur in either y-direction of the prism.
630  boost::unordered_map<int, int> splitEls;
632 
633  // Set up maps which takes an edge (in nektar++ ordering) and return
634  // their offset and stride in the 3d array of collapsed quadrature
635  // points. Note that this map includes only the edges that are on
636  // the triangular faces as the edges in the normal direction are
637  // linear.
638  map<LibUtilities::ShapeType, map<int, SplitMapHelper> > splitMap;
639  int po = nq * (nl + 1);
640 
641  SplitMapHelper splitPrism;
642  int splitMapEdgePrism[6] = {0, 2, 4, 5, 6, 7};
643  int splitMapOffsetPrism[6] = {0, nq, 0, nq - 1, nq + nq - 1, nq};
644  int splitMapIncPrism[6] = {1, 1, po, po, po, po};
645  int splitMapBFacesPrism[3] = {0, 2, 4};
646  int splitMapConnPrism[6][2] = {
647  {0, 0}, {1, 0}, {1, 1}, {0, 1}, {2, 0}, {2, 1}};
648  splitPrism.size = 6;
649  splitPrism.layerOff = nq;
650  splitPrism.edge = splitMapEdgePrism;
651  splitPrism.offset = splitMapOffsetPrism;
652  splitPrism.inc = splitMapIncPrism;
653  splitPrism.conn = helper2d(6, splitMapConnPrism);
654  splitPrism.bfacesSize = 3;
655  splitPrism.bfaces = splitMapBFacesPrism;
656  splitMap[LibUtilities::ePrism][1] = splitPrism;
657  splitMap[LibUtilities::ePrism][3] = splitPrism;
658 
659  int ho = nq * (nq - 1);
660  int tl = nq * nq;
661  SplitMapHelper splitHex0;
662  int splitMapEdgeHex0[8] = {0, 1, 2, 3, 8, 9, 10, 11};
663  int splitMapOffsetHex0[8] = {
664  0, nq - 1, tl - 1, ho, tl, tl + nq - 1, 2 * tl - 1, tl + ho};
665  int splitMapIncHex0[8] = {1, nq, -1, -nq, 1, nq, -1, -nq};
666  int splitMapBFacesHex0[4] = {1, 2, 3, 4};
667  int splitMapConnHex0[8][2] = {
668  {0, 0}, {1, 0}, {2, 0}, {3, 0}, {0, 1}, {1, 1}, {2, 1}, {3, 1}};
669  splitHex0.size = 8;
670  splitHex0.layerOff = nq * nq;
671  splitHex0.edge = splitMapEdgeHex0;
672  splitHex0.offset = splitMapOffsetHex0;
673  splitHex0.inc = splitMapIncHex0;
674  splitHex0.conn = helper2d(8, splitMapConnHex0);
675  splitHex0.bfacesSize = 4;
676  splitHex0.bfaces = splitMapBFacesHex0;
677  splitMap[LibUtilities::eHexahedron][0] = splitHex0;
678  splitMap[LibUtilities::eHexahedron][5] = splitHex0;
679 
680  // splitEdge enumerates the edges in the standard prism along which
681  // new nodes should be generated. These edges are the three between
682  // the two triangular faces.
683  //
684  // edgeVertMap specifies the vertices which comprise those edges in
685  // splitEdge; for example splitEdge[0] = 3 which connects vertices 0
686  // and 3.
687  //
688  // edgeOffset holds the offset of each of edges 3, 1 and 8
689  // respectively inside the collapsed coordinate system.
690  map<LibUtilities::ShapeType, map<int, SplitEdgeHelper> > splitEdge;
691 
692  int splitPrismEdges[3] = {3, 1, 8};
693  int splitPrismEdgeVert[3][2] = {{0, 3}, {1, 2}, {4, 5}};
694  int splitPrismOffset[3] = {0, nq - 1, nq * (nl + 1) * (nq - 1)};
695  int splitPrismInc[3] = {nq, nq, nq};
696  SplitEdgeHelper splitPrismEdge;
697  splitPrismEdge.size = 3;
698  splitPrismEdge.edge = splitPrismEdges;
699  splitPrismEdge.edgeVert = helper2d(3, splitPrismEdgeVert);
700  splitPrismEdge.offset = splitPrismOffset;
701  splitPrismEdge.inc = splitPrismInc;
702  splitEdge[LibUtilities::ePrism][1] = splitPrismEdge;
703  splitEdge[LibUtilities::ePrism][3] = splitPrismEdge;
704 
705  int splitHex0Edges[4] = {4, 5, 6, 7};
706  int splitHex0EdgeVert[4][2] = {{0, 4}, {1, 5}, {2, 6}, {3, 7}};
707  int splitHex0Offset[4] = {0, nq - 1, nq * nq - 1, nq * (nq - 1)};
708  int splitHex0Inc[4] = {nq * nq, nq * nq, nq * nq, nq * nq};
709  SplitEdgeHelper splitHex0Edge;
710  splitHex0Edge.size = 4;
711  splitHex0Edge.edge = splitHex0Edges;
712  splitHex0Edge.edgeVert = helper2d(4, splitHex0EdgeVert);
713  splitHex0Edge.offset = splitHex0Offset;
714  splitHex0Edge.inc = splitHex0Inc;
715  splitEdge[LibUtilities::eHexahedron][0] = splitHex0Edge;
716  splitEdge[LibUtilities::eHexahedron][5] = splitHex0Edge;
717 
718  map<LibUtilities::ShapeType, map<int, bool> > revPoints;
719  revPoints[LibUtilities::ePrism][1] = true;
720  revPoints[LibUtilities::ePrism][3] = false;
721 
722  revPoints[LibUtilities::eHexahedron][0] = true;
723  revPoints[LibUtilities::eHexahedron][5] = false;
724 
725  // edgeMap associates geometry edge IDs to the (nl+1) vertices which
726  // are generated along that edge when a prism is split, and is used
727  // to avoid generation of duplicate vertices. It is stored as an
728  // unordered map for speed.
729  boost::unordered_map<int, vector<NodeSharedPtr> > edgeMap;
730  boost::unordered_map<int, vector<NodeSharedPtr> >::iterator eIt;
731 
732  string surf = m_config["surf"].as<string>();
733  if (surf.size() > 0)
734  {
735  vector<unsigned int> surfs;
736  ParseUtils::GenerateSeqVector(surf.c_str(), surfs);
737  sort(surfs.begin(), surfs.end());
738 
739  // If surface is defined, process list of elements to find those
740  // that are connected to it.
741  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
742  {
743  ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
744  int nSurf = el->GetFaceCount();
745 
746  for (int j = 0; j < nSurf; ++j)
747  {
748  int bl = el->GetBoundaryLink(j);
749  if (bl == -1)
750  {
751  continue;
752  }
753 
754  ElementSharedPtr bEl =
755  m_mesh->m_element[m_mesh->m_expDim - 1][bl];
756  vector<int> tags = bEl->GetTagList();
757  vector<int> inter;
758 
759  sort(tags.begin(), tags.end());
760  set_intersection(surfs.begin(),
761  surfs.end(),
762  tags.begin(),
763  tags.end(),
764  back_inserter(inter));
765  ASSERTL0(inter.size() <= 1, "Intersection of surfaces wrong");
766 
767  if (inter.size() == 1)
768  {
769  if (el->GetConf().m_e == LibUtilities::ePrism)
770  {
771  if (j % 2 == 0)
772  {
773  cerr << "WARNING: Found quadrilateral face " << j
774  << " on surface " << surf
775  << " connected to prism; ignoring." << endl;
776  continue;
777  }
778 
779  if (splitEls.count(el->GetId()) > 0)
780  {
781  cerr << "WARNING: prism already found; "
782  << "ignoring" << endl;
783  }
784 
785  splitEls[el->GetId()] = j;
786  }
787  else if (validElTypes.count(el->GetConf().m_e) == 0)
788  {
789  cerr << "WARNING: Unsupported element type "
790  << "found in surface " << j << "; "
791  << "ignoring" << endl;
792  continue;
793  }
794  }
795  }
796  }
797  }
798  else
799  {
800  // Otherwise, add all prismatic elements and assume face 1 of
801  // the prism lies on the surface.
802  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
803  {
804  ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
805 
806  if (el->GetConf().m_e == LibUtilities::ePrism)
807  {
808  splitEls[el->GetId()] = 1;
809  }
810  else if (validElTypes.count(el->GetConf().m_e) > 0)
811  {
812  splitEls[el->GetId()] = 0;
813  }
814  else
815  {
816  continue;
817  }
818  }
819  }
820 
821  if (splitEls.size() == 0)
822  {
823  cerr << "WARNING: No elements detected to split." << endl;
824  return;
825  }
826 
827  // Erase all elements from the element list. Elements will be
828  // re-added as they are split.
829  vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim];
830  m_mesh->m_element[m_mesh->m_expDim].clear();
831 
832  map<int, SpatialDomains::Geometry3DSharedPtr> geomMap;
833  for (int i = 0; i < el.size(); ++i)
834  {
835  const int elId = el[i]->GetId();
836  sIt = splitEls.find(elId);
837  if (sIt == splitEls.end())
838  {
839  continue;
840  }
841 
842  // Get elemental geometry object and put into map.
843  geomMap[elId] = boost::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
844  el[i]->GetGeom(m_mesh->m_spaceDim));
845  }
846 
847  // Iterate over list of elements of expansion dimension.
848  for (int i = 0; i < el.size(); ++i)
849  {
850  const int elId = el[i]->GetId();
851  sIt = splitEls.find(elId);
852 
853  if (sIt == splitEls.end())
854  {
855  m_mesh->m_element[m_mesh->m_expDim].push_back(el[i]);
856  continue;
857  }
858 
859  SpatialDomains::Geometry3DSharedPtr geom = geomMap[elId];
860 
861  const int faceNum = sIt->second;
862  LibUtilities::ShapeType elType = el[i]->GetConf().m_e;
863 
864  SplitMapHelper &sMap = splitMap[elType][faceNum];
865  SplitEdgeHelper &sEdge = splitEdge[elType][faceNum];
866 
867  // Find quadrilateral boundary faces if any
868  std::map<int, int> bLink;
869  for (int j = 0; j < sMap.bfacesSize; ++j)
870  {
871  int bl = el[i]->GetBoundaryLink(sMap.bfaces[j]);
872  if (bl != -1)
873  {
874  bLink[sMap.bfaces[j]] = bl;
875  }
876  }
877 
878  // Determine whether to use reverse points.
880  revPoints[elType][faceNum] ? LibUtilities::eBoundaryLayerPoints
882 
883  // Determine value of r based on geometry.
884  if (ratioIsString)
885  {
886  NekDouble x, y, z;
887  NekDouble x1, y1, z1;
888  int nverts = geom->GetNumVerts();
889 
890  x = y = z = 0.0;
891 
892  for (int i = 0; i < nverts; ++i)
893  {
894  geom->GetVertex(i)->GetCoords(x1, y1, z1);
895  x += x1;
896  y += y1;
897  z += z1;
898  }
899  x /= (NekDouble)nverts;
900  y /= (NekDouble)nverts;
901  z /= (NekDouble)nverts;
902  r = rEval.Evaluate(rExprId, x, y, z, 0.0);
903  }
904 
906 
907  if (elType == LibUtilities::ePrism)
908  {
909  // Create basis.
913  2,
914  LibUtilities::PointsKey(nl + 1, t, r));
917 
918  // Create local region.
920  boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(geom);
922  B0, B1, B2, g);
923  }
924  else if (elType == LibUtilities::eHexahedron)
925  {
926  // Create basis.
930  2,
931  LibUtilities::PointsKey(nl + 1, t, r));
932 
933  // Create local region.
935  boost::dynamic_pointer_cast<SpatialDomains::HexGeom>(geom);
937  B0, B0, B1, g);
938  }
939 
940  // Grab co-ordinates.
941  Array<OneD, NekDouble> x(nq * nq * (nl + 1));
942  Array<OneD, NekDouble> y(nq * nq * (nl + 1));
943  Array<OneD, NekDouble> z(nq * nq * (nl + 1));
944  q->GetCoords(x, y, z);
945 
946  int nSplitEdge = sEdge.size;
947  vector<vector<NodeSharedPtr> > edgeNodes(nSplitEdge);
948 
949  // Loop over edges to be split.
950  for (int j = 0; j < nSplitEdge; ++j)
951  {
952  int locEdge = sEdge.edge[j];
953  int edgeId = el[i]->GetEdge(locEdge)->m_id;
954 
955  // Determine whether we have already generated vertices
956  // along this edge.
957  eIt = edgeMap.find(edgeId);
958 
959  if (eIt == edgeMap.end())
960  {
961  // If not then resize storage to hold new points.
962  edgeNodes[j].resize(nl + 1);
963 
964  // Re-use existing vertices at endpoints of edge to
965  // avoid duplicating the existing vertices.
966  edgeNodes[j][0] = el[i]->GetVertex(sEdge.edgeVert[j][0]);
967  edgeNodes[j][nl] = el[i]->GetVertex(sEdge.edgeVert[j][1]);
968 
969  // Variable geometric ratio
970  if (ratioIsString)
971  {
972  NekDouble x0, y0, z0;
973  NekDouble x1, y1, z1;
974  NekDouble xm, ym, zm;
975 
976  // -> Find edge end and mid points
977  x0 = x[sEdge.offset[j]];
978  y0 = y[sEdge.offset[j]];
979  z0 = z[sEdge.offset[j]];
980 
981  x1 = x[sEdge.offset[j] + nl * nq];
982  y1 = y[sEdge.offset[j] + nl * nq];
983  z1 = z[sEdge.offset[j] + nl * nq];
984 
985  xm = 0.5 * (x0 + x1);
986  ym = 0.5 * (y0 + y1);
987  zm = 0.5 * (z0 + z1);
988 
989  // evaluate r factor based on mid point value
990  NekDouble rnew;
991  rnew = rEval.Evaluate(rExprId, xm, ym, zm, 0.0);
992 
993  // Get basis with new r;
994  LibUtilities::PointsKey Pkey(nl + 1, t, rnew);
997 
998  const Array<OneD, const NekDouble> z = newP->GetZ();
999 
1000  // Create new interior nodes based on this new blend
1001  for (int k = 1; k < nl; ++k)
1002  {
1003  xm = 0.5 * (1 + z[k]) * (x1 - x0) + x0;
1004  ym = 0.5 * (1 + z[k]) * (y1 - y0) + y0;
1005  zm = 0.5 * (1 + z[k]) * (z1 - z0) + z0;
1006  edgeNodes[j][k] =
1007  NodeSharedPtr(new Node(nodeId++, xm, ym, zm));
1008  }
1009  }
1010  else
1011  {
1012  // Create new interior nodes.
1013  for (int k = 1; k < nl; ++k)
1014  {
1015  int pos = sEdge.offset[j] + k * sEdge.inc[j];
1016  edgeNodes[j][k] = NodeSharedPtr(
1017  new Node(nodeId++, x[pos], y[pos], z[pos]));
1018  }
1019  }
1020 
1021  // Store these edges in edgeMap.
1022  edgeMap[edgeId] = edgeNodes[j];
1023  }
1024  else
1025  {
1026  edgeNodes[j] = eIt->second;
1027  }
1028  }
1029 
1030  // Create element layers.
1031  for (int j = 0; j < nl; ++j)
1032  {
1033  // Offset of this layer within the collapsed coordinate
1034  // system.
1035  int offset = j * sMap.layerOff;
1036 
1037  // Get corner vertices.
1038  vector<NodeSharedPtr> nodeList(sMap.size);
1039  for (int k = 0; k < sMap.size; ++k)
1040  {
1041  nodeList[k] = edgeNodes[sMap.conn[k][0]][j + sMap.conn[k][1]];
1042  }
1043 
1044  // Create the element.
1045  ElmtConfig conf(elType, 1, true, true, false);
1047  elType, conf, nodeList, el[i]->GetTagList());
1048 
1049  // Add high order nodes to split prismatic edges.
1050  for (int l = 0; l < sMap.size; ++l)
1051  {
1052  EdgeSharedPtr HOedge = elmt->GetEdge(sMap.edge[l]);
1053  for (int k = 1; k < nq - 1; ++k)
1054  {
1055  int pos = offset + sMap.offset[l] + k * sMap.inc[l];
1056  HOedge->m_edgeNodes.push_back(NodeSharedPtr(
1057  new Node(nodeId++, x[pos], y[pos], z[pos])));
1058  }
1059  HOedge->m_curveType = pt;
1060  }
1061 
1062  // Change the surface elements to match the layers of
1063  // elements on the boundary of the domain.
1065  for (it = bLink.begin(); it != bLink.end(); ++it)
1066  {
1067  int fid = it->first;
1068  int bl = it->second;
1069 
1070  vector<NodeSharedPtr> qNodeList(4);
1071  for (int k = 0; k < 4; ++k)
1072  {
1073  qNodeList[k] = nodeList[faceNodeMap[elType][fid][k]];
1074  }
1075  vector<int> tagBE;
1076  tagBE =
1077  m_mesh->m_element[m_mesh->m_expDim - 1][bl]->GetTagList();
1078  ElmtConfig bconf(
1079  LibUtilities::eQuadrilateral, 1, true, true, false);
1080  ElementSharedPtr boundaryElmt =
1082  LibUtilities::eQuadrilateral, bconf, qNodeList, tagBE);
1083 
1084  // Overwrite first layer boundary element with new
1085  // boundary element, otherwise push this back to end of
1086  // the boundary list
1087  if (j == 0)
1088  {
1089  m_mesh->m_element[m_mesh->m_expDim - 1][bl] = boundaryElmt;
1090  }
1091  else
1092  {
1093  m_mesh->m_element[m_mesh->m_expDim - 1].push_back(
1094  boundaryElmt);
1095  }
1096  }
1097 
1098  m_mesh->m_element[m_mesh->m_expDim].push_back(elmt);
1099  }
1100  }
1101 }
1102 }
1103 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Basic information about an element.
Definition: ElementConfig.h:50
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
Principle Modified Functions .
Definition: BasisType.h:49
STL namespace.
pair< ModuleType, string > ModuleKey
boost::shared_ptr< Points< NekDouble > > PointsSharedPtr
boost::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:110
int ** helper2d(int lda, int arr[][2])
Definition: ProcessBL.cpp:62
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:79
virtual NEKMESHUTILS_EXPORT void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
virtual void Process()
Write mesh to output file.
Definition: ProcessBL.cpp:126
Principle Modified Functions .
Definition: BasisType.h:50
virtual NEKMESHUTILS_EXPORT void ProcessElements()
Generate element IDs.
Represents a command-line configuration option.
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
Definition: Points.h:58
boost::shared_ptr< QuadExp > QuadExpSharedPtr
Definition: QuadExp.h:297
double NekDouble
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
This class defines evaluator of analytic (symbolic) mathematical expressions. Expressions are allowed...
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:136
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
1D power law distribution for boundary layer points
Definition: PointsType.h:69
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:147
3D geometry information
Definition: Geometry3D.h:70
Abstract base class for processing modules.
virtual NEKMESHUTILS_EXPORT void ProcessVertices()
Extract element vertices.
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:109
1D power law distribution for boundary layer points
Definition: PointsType.h:68
virtual NEKMESHUTILS_EXPORT void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
std::pair< ModuleType, std::string > ModuleKey
Describes the specification for a Basis.
Definition: Basis.h:50
virtual NEKMESHUTILS_EXPORT void ProcessComposites()
Generate composites.
int ** helper2d(int lda, int arr[][4])
Definition: ProcessBL.cpp:74
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
boost::shared_ptr< Geometry3D > Geometry3DSharedPtr
Definition: Geometry3D.h:52