Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 boundary layer elements.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <string>
37 using namespace std;
38 
39 #include "MeshElements.h"
40 #include "ProcessBL.h"
41 
47 #include <LocalRegions/PrismExp.h>
48 #include <LocalRegions/HexExp.h>
49 
50 namespace Nektar
51 {
52  namespace Utilities
53  {
54  ModuleKey ProcessBL::className =
56  ModuleKey(eProcessModule, "bl"), ProcessBL::create,
57  "Refines a prismatic boundary layer.");
58 
59  int **helper2d(int lda, int arr[][2])
60  {
61  int **ret = new int*[lda];
62  for (int i = 0; i < lda; ++i)
63  {
64  ret[i] = new int[2];
65  ret[i][0] = arr[i][0];
66  ret[i][1] = arr[i][1];
67  }
68  return ret;
69  }
70 
71  int **helper2d(int lda, int arr[][4])
72  {
73  int **ret = new int*[lda];
74  for (int i = 0; i < lda; ++i)
75  {
76  ret[i] = new int[4];
77  ret[i][0] = arr[i][0];
78  ret[i][1] = arr[i][1];
79  ret[i][2] = arr[i][2];
80  ret[i][3] = arr[i][3];
81  }
82  return ret;
83  }
84 
86  {
87  int size;
88  int layerOff;
89  int *edge;
90  int *offset;
91  int *inc;
92  int **conn;
94  int *bfaces;
95  };
96 
98  {
99  int size;
100  int *edge;
101  int **edgeVert;
102  int *offset;
103  int *inc;
104  };
105 
106  ProcessBL::ProcessBL(MeshSharedPtr m) : ProcessModule(m)
107  {
108  // BL mesh configuration.
109  m_config["layers"] = ConfigOption(false, "2",
110  "Number of layers to refine.");
111  m_config["nq"] = ConfigOption(false, "5",
112  "Number of points in high order elements.");
113  m_config["surf"] = ConfigOption(false, "",
114  "Tag identifying surface connected to prism.");
115  m_config["r"] = ConfigOption(false, "2.0",
116  "Ratio to use in geometry progression.");
117  }
118 
120  {
121 
122  }
123 
125  {
126  if (m_mesh->m_verbose)
127  {
128  cout << "ProcessBL: Refining prismatic boundary layer..."
129  << endl;
130  }
131 
132  // A set containing all element types which are valid.
133  set<LibUtilities::ShapeType> validElTypes;
134  validElTypes.insert(LibUtilities::ePrism);
135  validElTypes.insert(LibUtilities::eHexahedron);
136 
137  int nodeId = m_mesh->m_vertexSet.size();
138  int nl = m_config["layers"].as<int>();
139  int nq = m_config["nq"]. as<int>();
140 
141  // determine if geometric ratio is string or a constant.
143  NekDouble r = 1;
144  int rExprId = -1;
145  bool ratioIsString = false;
146 
147  if (m_config["r"].isType<NekDouble>())
148  {
149  r = m_config["r"].as<NekDouble>();
150  }
151  else
152  {
153  std::string rstr = m_config["r"].as<string>();
154  rExprId = rEval.DefineFunction("x y z", rstr);
155  ratioIsString = true;
156  }
157 
158  // Prismatic node -> face map.
159  int prismFaceNodes[5][4] = {
160  {0,1,2,3},{0,1,4,-1},{1,2,5,4},{3,2,5,-1},{0,3,5,4}};
161  int hexFaceNodes [6][4] = {
162  {0,1,2,3},{0,1,5,4},{1,2,6,5},{3,2,6,7},{0,3,7,4},{4,5,6,7}};
163  map<LibUtilities::ShapeType, int **> faceNodeMap;
164  faceNodeMap[LibUtilities::ePrism] = helper2d(5, prismFaceNodes);
165  faceNodeMap[LibUtilities::eHexahedron] = helper2d(6, hexFaceNodes);
166 
167  // Default PointsType.
169 
170  // Map which takes element ID to face on surface. This enables
171  // splitting to occur in either y-direction of the prism.
172  boost::unordered_map<int, int> splitEls;
174 
175  // Set up maps which takes an edge (in nektar++ ordering) and return
176  // their offset and stride in the 3d array of collapsed quadrature
177  // points. Note that this map includes only the edges that are on
178  // the triangular faces as the edges in the normal direction are
179  // linear.
180  map<LibUtilities::ShapeType, map<int, SplitMapHelper> > splitMap;
181  int po = nq*(nl+1);
182 
183  SplitMapHelper splitPrism;
184  int splitMapEdgePrism [6] = {0, 2, 4, 5, 6, 7};
185  int splitMapOffsetPrism[6] = {0, nq, 0, nq-1, nq+nq-1, nq};
186  int splitMapIncPrism [6] = {1, 1, po, po, po, po};
187  int splitMapBFacesPrism[3] = {0, 2, 4};
188  int splitMapConnPrism [6][2] = {{0,0}, {1,0}, {1,1},
189  {0,1}, {2,0}, {2,1}};
190  splitPrism.size = 6;
191  splitPrism.layerOff = nq;
192  splitPrism.edge = splitMapEdgePrism;
193  splitPrism.offset = splitMapOffsetPrism;
194  splitPrism.inc = splitMapIncPrism;
195  splitPrism.conn = helper2d(6, splitMapConnPrism);
196  splitPrism.bfacesSize = 3;
197  splitPrism.bfaces = splitMapBFacesPrism;
198  splitMap[LibUtilities::ePrism][1] = splitPrism;
199  splitMap[LibUtilities::ePrism][3] = splitPrism;
200 
201  int ho = nq*(nq-1);
202  int tl = nq*nq;
203  SplitMapHelper splitHex0;
204  int splitMapEdgeHex0 [8] = {0, 1, 2, 3, 8, 9, 10, 11};
205  int splitMapOffsetHex0[8] = {0, nq-1, tl-1, ho, tl, tl+nq-1, 2*tl-1, tl+ho};
206  int splitMapIncHex0 [8] = {1, nq, -1, -nq, 1, nq, -1, -nq};
207  int splitMapBFacesHex0[4] = {1, 2, 3, 4};
208  int splitMapConnHex0 [8][2] = {{0,0}, {1,0}, {2,0}, {3,0},
209  {0,1}, {1,1}, {2,1}, {3,1}};
210  splitHex0.size = 8;
211  splitHex0.layerOff = nq*nq;
212  splitHex0.edge = splitMapEdgeHex0;
213  splitHex0.offset = splitMapOffsetHex0;
214  splitHex0.inc = splitMapIncHex0;
215  splitHex0.conn = helper2d(8, splitMapConnHex0);
216  splitHex0.bfacesSize = 4;
217  splitHex0.bfaces = splitMapBFacesHex0;
218  splitMap[LibUtilities::eHexahedron][0] = splitHex0;
219  splitMap[LibUtilities::eHexahedron][5] = splitHex0;
220 
221  // splitEdge enumerates the edges in the standard prism along which
222  // new nodes should be generated. These edges are the three between
223  // the two triangular faces.
224  //
225  // edgeVertMap specifies the vertices which comprise those edges in
226  // splitEdge; for example splitEdge[0] = 3 which connects vertices 0
227  // and 3.
228  //
229  // edgeOffset holds the offset of each of edges 3, 1 and 8
230  // respectively inside the collapsed coordinate system.
231  map<LibUtilities::ShapeType, map<int, SplitEdgeHelper> > splitEdge;
232 
233  int splitPrismEdges [3] = {3, 1, 8};
234  int splitPrismEdgeVert[3][2] = {{0,3}, {1,2}, {4,5}};
235  int splitPrismOffset [3] = {0, nq-1, nq*(nl+1)*(nq-1)};
236  int splitPrismInc [3] = {nq, nq, nq};
237  SplitEdgeHelper splitPrismEdge;
238  splitPrismEdge.size = 3;
239  splitPrismEdge.edge = splitPrismEdges;
240  splitPrismEdge.edgeVert = helper2d(3, splitPrismEdgeVert);
241  splitPrismEdge.offset = splitPrismOffset;
242  splitPrismEdge.inc = splitPrismInc;
243  splitEdge[LibUtilities::ePrism][1] = splitPrismEdge;
244  splitEdge[LibUtilities::ePrism][3] = splitPrismEdge;
245 
246  int splitHex0Edges [4] = {4, 5, 6, 7};
247  int splitHex0EdgeVert[4][2] = {{0,4}, {1,5}, {2,6}, {3,7}};
248  int splitHex0Offset [4] = {0, nq-1, nq*nq-1, nq*(nq-1) };
249  int splitHex0Inc [4] = {nq*nq, nq*nq, nq*nq, nq*nq};
250  SplitEdgeHelper splitHex0Edge;
251  splitHex0Edge.size = 4;
252  splitHex0Edge.edge = splitHex0Edges;
253  splitHex0Edge.edgeVert = helper2d(4, splitHex0EdgeVert);
254  splitHex0Edge.offset = splitHex0Offset;
255  splitHex0Edge.inc = splitHex0Inc;
256  splitEdge[LibUtilities::eHexahedron][0] = splitHex0Edge;
257  splitEdge[LibUtilities::eHexahedron][5] = splitHex0Edge;
258 
259  map<LibUtilities::ShapeType, map<int, bool> > revPoints;
260  revPoints[LibUtilities::ePrism][1] = false;
261  revPoints[LibUtilities::ePrism][3] = true;
262 
263  revPoints[LibUtilities::eHexahedron][0] = true;
264  revPoints[LibUtilities::eHexahedron][5] = false;
265 
266  // edgeMap associates geometry edge IDs to the (nl+1) vertices which
267  // are generated along that edge when a prism is split, and is used
268  // to avoid generation of duplicate vertices. It is stored as an
269  // unordered map for speed.
270  boost::unordered_map<int, vector<NodeSharedPtr> > edgeMap;
271  boost::unordered_map<int, vector<NodeSharedPtr> >::iterator eIt;
272 
273  string surf = m_config["surf"].as<string>();
274  if (surf.size() > 0)
275  {
276  vector<unsigned int> surfs;
277  ParseUtils::GenerateSeqVector(surf.c_str(), surfs);
278  sort(surfs.begin(), surfs.end());
279 
280  // If surface is defined, process list of elements to find those
281  // that are connected to it.
282  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
283  {
284  ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
285  int nSurf = el->GetFaceCount();
286 
287  for (int j = 0; j < nSurf; ++j)
288  {
289  int bl = el->GetBoundaryLink(j);
290  if (bl == -1)
291  {
292  continue;
293  }
294 
295  ElementSharedPtr bEl = m_mesh->m_element[m_mesh->m_expDim-1][bl];
296  vector<int> tags = bEl->GetTagList();
297  vector<int> inter;
298 
299  sort(tags.begin(), tags.end());
300  set_intersection(surfs.begin(), surfs.end(),
301  tags .begin(), tags .end(),
302  back_inserter(inter));
303  ASSERTL0(inter.size() <= 1,
304  "Intersection of surfaces wrong");
305 
306  if (inter.size() == 1)
307  {
308  if (el->GetConf().m_e == LibUtilities::ePrism)
309  {
310  if (j % 2 == 0)
311  {
312  cerr << "WARNING: Found quadrilateral face "
313  << j << " on surface " << surf
314  << " connected to prism; ignoring."
315  << endl;
316  continue;
317  }
318 
319  if (splitEls.count(el->GetId()) > 0)
320  {
321  cerr << "WARNING: prism already found; "
322  << "ignoring" << endl;
323  }
324 
325  splitEls[el->GetId()] = j;
326  }
327  else if (validElTypes.count(el->GetConf().m_e) == 0)
328  {
329  cerr << "WARNING: Unsupported element type "
330  << "found in surface " << j << "; "
331  << "ignoring" << endl;
332  continue;
333  }
334  }
335  }
336  }
337  }
338  else
339  {
340  // Otherwise, add all prismatic elements and assume face 1 of
341  // the prism lies on the surface.
342  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
343  {
344  ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
345 
346  if (el->GetConf().m_e == LibUtilities::ePrism)
347  {
348  splitEls[el->GetId()] = 1;
349  }
350  else if (validElTypes.count(el->GetConf().m_e) > 0)
351  {
352  splitEls[el->GetId()] = 0;
353  }
354  else
355  {
356  continue;
357  }
358  }
359  }
360 
361  if (splitEls.size() == 0)
362  {
363  cerr << "WARNING: No elements detected to split." << endl;
364  return;
365  }
366 
367  // Erase all elements from the element list. Elements will be
368  // re-added as they are split.
369  vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim];
370  m_mesh->m_element[m_mesh->m_expDim].clear();
371 
372  // Iterate over list of elements of expansion dimension.
373  for (int i = 0; i < el.size(); ++i)
374  {
375  const int elId = el[i]->GetId();
376  sIt = splitEls.find(elId);
377 
378  if (sIt == splitEls.end())
379  {
380  m_mesh->m_element[m_mesh->m_expDim].push_back(el[i]);
381  continue;
382  }
383 
384  const int faceNum = sIt->second;
385  LibUtilities::ShapeType elType = el[i]->GetConf().m_e;
386 
387  SplitMapHelper &sMap = splitMap [elType][faceNum];
388  SplitEdgeHelper &sEdge = splitEdge[elType][faceNum];
389 
390  // Find quadrilateral boundary faces if any
391  std::map<int, int> bLink;
392  for (int j = 0; j < sMap.bfacesSize; ++j)
393  {
394  int bl = el[i]->GetBoundaryLink(sMap.bfaces[j]);
395  if (bl != -1)
396  {
397  bLink[sMap.bfaces[j]] = bl;
398  }
399  }
400 
401  // Get elemental geometry object.
403  boost::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
404  el[i]->GetGeom(m_mesh->m_spaceDim));
405 
406  // Determine whether to use reverse points.
408  revPoints[elType][faceNum] ?
411 
412  // Determine value of r based on geometry.
413  if(ratioIsString)
414  {
415  NekDouble x, y, z;
416  NekDouble x1, y1, z1;
417  int nverts = geom->GetNumVerts();
418 
419  x = y = z = 0.0;
420 
421  for(int i = 0; i < nverts; ++i)
422  {
423  geom->GetVertex(i)->GetCoords(x1,y1,z1);
424  x += x1; y += y1; z += z1;
425  }
426  x /= (NekDouble) nverts;
427  y /= (NekDouble) nverts;
428  z /= (NekDouble) nverts;
429  r = rEval.Evaluate(rExprId,x,y,z,0.0);
430  }
431 
433 
434  if (elType == LibUtilities::ePrism)
435  {
436  // Create basis.
439  LibUtilities::PointsKey(nq,pt));
442  LibUtilities::PointsKey(nl+1, t, r));
445  LibUtilities::PointsKey(nq,pt));
446 
447  // Create local region.
449  boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(
450  geom);
452  B0, B1, B2, g);
453  }
454  else if (elType == LibUtilities::eHexahedron)
455  {
456  // Create basis.
459  LibUtilities::PointsKey(nq,pt));
462  LibUtilities::PointsKey(nl+1, t, r));
463 
464  // Create local region.
466  boost::dynamic_pointer_cast<SpatialDomains::HexGeom>(
467  geom);
469  B0, B0, B1, g);
470  }
471 
472  // Grab co-ordinates.
473  Array<OneD, NekDouble> x(nq*nq*(nl+1));
474  Array<OneD, NekDouble> y(nq*nq*(nl+1));
475  Array<OneD, NekDouble> z(nq*nq*(nl+1));
476  q->GetCoords(x,y,z);
477 
478  int nSplitEdge = sEdge.size;
479  vector<vector<NodeSharedPtr> > edgeNodes(nSplitEdge);
480 
481  // Loop over edges to be split.
482  for (int j = 0; j < nSplitEdge; ++j)
483  {
484  int locEdge = sEdge.edge[j];
485  int edgeId = el[i]->GetEdge(locEdge)->m_id;
486 
487  // Determine whether we have already generated vertices
488  // along this edge.
489  eIt = edgeMap.find(edgeId);
490 
491  if (eIt == edgeMap.end())
492  {
493  // If not then resize storage to hold new points.
494  edgeNodes[j].resize(nl+1);
495 
496  // Re-use existing vertices at endpoints of edge to
497  // avoid duplicating the existing vertices.
498  edgeNodes[j][0] = el[i]->GetVertex(sEdge.edgeVert[j][0]);
499  edgeNodes[j][nl] = el[i]->GetVertex(sEdge.edgeVert[j][1]);
500 
501  // Variable geometric ratio
502  if(ratioIsString)
503  {
504  NekDouble x0,y0,z0;
505  NekDouble x1,y1,z1;
506  NekDouble xm,ym,zm;
507 
508  // -> Find edge end and mid points
509  x0 = x[sEdge.offset[j]];
510  y0 = y[sEdge.offset[j]];
511  z0 = z[sEdge.offset[j]];
512 
513  x1 = x[sEdge.offset[j]+nl*nq];
514  y1 = y[sEdge.offset[j]+nl*nq];
515  z1 = z[sEdge.offset[j]+nl*nq];
516 
517  xm = 0.5*(x0+x1);
518  ym = 0.5*(y0+y1);
519  zm = 0.5*(z0+z1);
520 
521  // evaluate r factor based on mid point value
522  NekDouble rnew;
523  rnew = rEval.Evaluate(rExprId,xm,ym,zm,0.0);
524 
525  // Get basis with new r;
526  LibUtilities::PointsKey Pkey(nl+1, t, rnew);
528  = LibUtilities::PointsManager()[Pkey];
529 
530  const Array<OneD, const NekDouble> z = newP->GetZ();
531 
532  // Create new interior nodes based on this new blend
533  for (int k = 1; k < nl; ++k)
534  {
535  xm = 0.5*(1+z[k])*(x1-x0) + x0;
536  ym = 0.5*(1+z[k])*(y1-y0) + y0;
537  zm = 0.5*(1+z[k])*(z1-z0) + z0;
538  edgeNodes[j][k] = NodeSharedPtr(
539  new Node(nodeId++, xm,ym,zm));
540  }
541  }
542  else
543  {
544  // Create new interior nodes.
545  for (int k = 1; k < nl; ++k)
546  {
547  int pos = sEdge.offset[j] + k*sEdge.inc[j];
548  edgeNodes[j][k] = NodeSharedPtr(
549  new Node(nodeId++, x[pos], y[pos], z[pos]));
550  }
551  }
552 
553  // Store these edges in edgeMap.
554  edgeMap[edgeId] = edgeNodes[j];
555  }
556  else
557  {
558  edgeNodes[j] = eIt->second;
559  }
560  }
561 
562  // Create element layers.
563  for (int j = 0; j < nl; ++j)
564  {
565  // Offset of this layer within the collapsed coordinate
566  // system.
567  int offset = j * sMap.layerOff;
568 
569  // Get corner vertices.
570  vector<NodeSharedPtr> nodeList(sMap.size);
571  for (int k = 0; k < sMap.size; ++k)
572  {
573  nodeList[k] =
574  edgeNodes[sMap.conn[k][0]][j + sMap.conn[k][1]];
575  }
576 
577  // Create the element.
578  ElmtConfig conf(elType, 1, true, true, false);
580  CreateInstance(
581  elType, conf, nodeList, el[i]->GetTagList());
582 
583  // Add high order nodes to split prismatic edges.
584  for (int l = 0; l < sMap.size; ++l)
585  {
586  EdgeSharedPtr HOedge = elmt->GetEdge(
587  sMap.edge[l]);
588  for (int k = 1; k < nq-1; ++k)
589  {
590  int pos = offset + sMap.offset[l] + k*sMap.inc[l];
591  HOedge->m_edgeNodes.push_back(
593  new Node(nodeId++,x[pos],y[pos],z[pos])));
594  }
595  HOedge->m_curveType = pt;
596  }
597 
598  // Change the surface elements to match the layers of
599  // elements on the boundary of the domain.
601  for (it = bLink.begin(); it != bLink.end(); ++it)
602  {
603  int fid = it->first;
604  int bl = it->second;
605 
606  if (j == 0)
607  {
608  // For first layer reuse existing 2D element.
609  ElementSharedPtr e = m_mesh->m_element[m_mesh->m_expDim-1][bl];
610  for (int k = 0; k < 4; ++k)
611  {
612  e->SetVertex(
613  k, nodeList[faceNodeMap[elType][fid][k]]);
614  }
615  }
616  else
617  {
618  // For all other layers create new element.
619  vector<NodeSharedPtr> qNodeList(4);
620  for (int k = 0; k < 4; ++k)
621  {
622  qNodeList[k] = nodeList[faceNodeMap[elType][fid][k]];
623  }
624  vector<int> tagBE;
625  tagBE = m_mesh->m_element[m_mesh->m_expDim-1][bl]->GetTagList();
626  ElmtConfig bconf(LibUtilities::eQuadrilateral,1,true,true,false);
627  ElementSharedPtr boundaryElmt = GetElementFactory().
628  CreateInstance(LibUtilities::eQuadrilateral,bconf,
629  qNodeList,tagBE);
630  m_mesh->m_element[m_mesh->m_expDim-1].push_back(boundaryElmt);
631  }
632  }
633 
634  m_mesh->m_element[m_mesh->m_expDim].push_back(elmt);
635  }
636  }
637 
638  // Re-process mesh to eliminate duplicate vertices and edges.
639  ProcessVertices();
640  ProcessEdges();
641  ProcessFaces();
642  ProcessElements();
644  }
645  }
646 }