Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | List of all members
Nektar::Utilities::ProcessBL Class Reference

This processing module calculates the Jacobian of elements using SpatialDomains::GeomFactors and the Element::GetGeom method. For now it simply prints a list of elements which have negative Jacobian. More...

#include <ProcessBL.h>

Inheritance diagram for Nektar::Utilities::ProcessBL:
Inheritance graph
[legend]
Collaboration diagram for Nektar::Utilities::ProcessBL:
Collaboration graph
[legend]

Public Member Functions

 ProcessBL (MeshSharedPtr m)
 
virtual ~ProcessBL ()
 
virtual void Process ()
 Write mesh to output file. More...
 
- Public Member Functions inherited from Nektar::Utilities::ProcessModule
 ProcessModule ()
 
 ProcessModule (FieldSharedPtr p_f)
 
 ProcessModule (MeshSharedPtr p_m)
 
- Public Member Functions inherited from Nektar::Utilities::Module
 Module (FieldSharedPtr p_f)
 
virtual void Process (po::variables_map &vm)=0
 
void RegisterConfig (string key, string value)
 Register a configuration option with a module. More...
 
void PrintConfig ()
 Print out all configuration options for a module. More...
 
void SetDefaults ()
 Sets default configuration options for those which have not been set. More...
 
bool GetRequireEquiSpaced (void)
 
void SetRequireEquiSpaced (bool pVal)
 
void EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)
 
 Module (MeshSharedPtr p_m)
 
void RegisterConfig (string key, string value)
 
void PrintConfig ()
 
void SetDefaults ()
 
MeshSharedPtr GetMesh ()
 
virtual void ProcessVertices ()
 Extract element vertices. More...
 
virtual void ProcessEdges (bool ReprocessEdges=true)
 Extract element edges. More...
 
virtual void ProcessFaces (bool ReprocessFaces=true)
 Extract element faces. More...
 
virtual void ProcessElements ()
 Generate element IDs. More...
 
virtual void ProcessComposites ()
 Generate composites. More...
 
virtual void ClearElementLinks ()
 

Static Public Member Functions

static boost::shared_ptr< Modulecreate (MeshSharedPtr m)
 Creates an instance of this class. More...
 

Static Public Attributes

static ModuleKey className
 

Additional Inherited Members

- Protected Member Functions inherited from Nektar::Utilities::Module
 Module ()
 
void ReorderPrisms (PerMap &perFaces)
 Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...
 
void PrismLines (int prism, PerMap &perFaces, set< int > &prismsDone, vector< ElementSharedPtr > &line)
 
- Protected Attributes inherited from Nektar::Utilities::Module
FieldSharedPtr m_f
 Field object. More...
 
map< string, ConfigOptionm_config
 List of configuration values. More...
 
bool m_requireEquiSpaced
 
MeshSharedPtr m_mesh
 Mesh object. More...
 

Detailed Description

This processing module calculates the Jacobian of elements using SpatialDomains::GeomFactors and the Element::GetGeom method. For now it simply prints a list of elements which have negative Jacobian.

Definition at line 52 of file ProcessBL.h.

Constructor & Destructor Documentation

Nektar::Utilities::ProcessBL::ProcessBL ( MeshSharedPtr  m)

Definition at line 108 of file ProcessBL.cpp.

References Nektar::Utilities::Module::m_config.

108  : ProcessModule(m)
109 {
110  // BL mesh configuration.
111  m_config["layers"] =
112  ConfigOption(false, "2", "Number of layers to refine.");
113  m_config["nq"] =
114  ConfigOption(false, "5", "Number of points in high order elements.");
115  m_config["surf"] =
116  ConfigOption(false, "", "Tag identifying surface connected to prism.");
117  m_config["r"] =
118  ConfigOption(false, "2.0", "Ratio to use in geometry progression.");
119 }
map< string, ConfigOption > m_config
List of configuration values.
Nektar::Utilities::ProcessBL::~ProcessBL ( )
virtual

Definition at line 121 of file ProcessBL.cpp.

122 {
123 }

Member Function Documentation

static boost::shared_ptr<Module> Nektar::Utilities::ProcessBL::create ( MeshSharedPtr  m)
inlinestatic

Creates an instance of this class.

Definition at line 56 of file ProcessBL.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

57  {
59  }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Nektar::Utilities::ProcessBL::Process ( )
virtual

Write mesh to output file.

Implements Nektar::Utilities::Module.

Definition at line 125 of file ProcessBL.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::Utilities::SplitMapHelper::bfaces, Nektar::Utilities::SplitMapHelper::bfacesSize, Nektar::Utilities::SplitMapHelper::conn, Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::LibUtilities::eBoundaryLayerPoints, Nektar::LibUtilities::eBoundaryLayerPointsRev, Nektar::Utilities::SplitMapHelper::edge, Nektar::Utilities::SplitEdgeHelper::edge, Nektar::Utilities::SplitEdgeHelper::edgeVert, Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::eQuadrilateral, Nektar::ParseUtils::GenerateSeqVector(), Nektar::NekMeshUtils::GetElementFactory(), Nektar::Utilities::helper2d(), Nektar::Utilities::SplitMapHelper::inc, Nektar::Utilities::SplitEdgeHelper::inc, Nektar::iterator, Nektar::Utilities::SplitMapHelper::layerOff, Nektar::Utilities::Module::m_config, Nektar::Utilities::Module::m_mesh, Nektar::Utilities::SplitMapHelper::offset, Nektar::Utilities::SplitEdgeHelper::offset, Nektar::LibUtilities::PointsManager(), Nektar::Utilities::Module::ProcessComposites(), Nektar::Utilities::Module::ProcessEdges(), Nektar::Utilities::Module::ProcessElements(), Nektar::Utilities::Module::ProcessFaces(), Nektar::Utilities::Module::ProcessVertices(), Nektar::Utilities::SplitMapHelper::size, Nektar::Utilities::SplitEdgeHelper::size, and Nektar::NekMeshUtils::surf.

126 {
127  if (m_mesh->m_verbose)
128  {
129  cout << "ProcessBL: Refining prismatic boundary layer..." << 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.
142  LibUtilities::AnalyticExpressionEvaluator rEval;
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] = {{0, 1, 2, 3},
162  {0, 1, 5, 4},
163  {1, 2, 6, 5},
164  {3, 2, 6, 7},
165  {0, 3, 7, 4},
166  {4, 5, 6, 7}};
167  map<LibUtilities::ShapeType, int **> faceNodeMap;
168  faceNodeMap[LibUtilities::ePrism] = helper2d(5, prismFaceNodes);
169  faceNodeMap[LibUtilities::eHexahedron] = helper2d(6, hexFaceNodes);
170 
171  // Default PointsType.
173 
174  // Map which takes element ID to face on surface. This enables
175  // splitting to occur in either y-direction of the prism.
176  boost::unordered_map<int, int> splitEls;
178 
179  // Set up maps which takes an edge (in nektar++ ordering) and return
180  // their offset and stride in the 3d array of collapsed quadrature
181  // points. Note that this map includes only the edges that are on
182  // the triangular faces as the edges in the normal direction are
183  // linear.
184  map<LibUtilities::ShapeType, map<int, SplitMapHelper> > splitMap;
185  int po = nq * (nl + 1);
186 
187  SplitMapHelper splitPrism;
188  int splitMapEdgePrism[6] = {0, 2, 4, 5, 6, 7};
189  int splitMapOffsetPrism[6] = {0, nq, 0, nq - 1, nq + nq - 1, nq};
190  int splitMapIncPrism[6] = {1, 1, po, po, po, po};
191  int splitMapBFacesPrism[3] = {0, 2, 4};
192  int splitMapConnPrism[6][2] = {
193  {0, 0}, {1, 0}, {1, 1}, {0, 1}, {2, 0}, {2, 1}};
194  splitPrism.size = 6;
195  splitPrism.layerOff = nq;
196  splitPrism.edge = splitMapEdgePrism;
197  splitPrism.offset = splitMapOffsetPrism;
198  splitPrism.inc = splitMapIncPrism;
199  splitPrism.conn = helper2d(6, splitMapConnPrism);
200  splitPrism.bfacesSize = 3;
201  splitPrism.bfaces = splitMapBFacesPrism;
202  splitMap[LibUtilities::ePrism][1] = splitPrism;
203  splitMap[LibUtilities::ePrism][3] = splitPrism;
204 
205  int ho = nq * (nq - 1);
206  int tl = nq * nq;
207  SplitMapHelper splitHex0;
208  int splitMapEdgeHex0[8] = {0, 1, 2, 3, 8, 9, 10, 11};
209  int splitMapOffsetHex0[8] = {
210  0, nq - 1, tl - 1, ho, tl, tl + nq - 1, 2 * tl - 1, tl + ho};
211  int splitMapIncHex0[8] = {1, nq, -1, -nq, 1, nq, -1, -nq};
212  int splitMapBFacesHex0[4] = {1, 2, 3, 4};
213  int splitMapConnHex0[8][2] = {
214  {0, 0}, {1, 0}, {2, 0}, {3, 0}, {0, 1}, {1, 1}, {2, 1}, {3, 1}};
215  splitHex0.size = 8;
216  splitHex0.layerOff = nq * nq;
217  splitHex0.edge = splitMapEdgeHex0;
218  splitHex0.offset = splitMapOffsetHex0;
219  splitHex0.inc = splitMapIncHex0;
220  splitHex0.conn = helper2d(8, splitMapConnHex0);
221  splitHex0.bfacesSize = 4;
222  splitHex0.bfaces = splitMapBFacesHex0;
223  splitMap[LibUtilities::eHexahedron][0] = splitHex0;
224  splitMap[LibUtilities::eHexahedron][5] = splitHex0;
225 
226  // splitEdge enumerates the edges in the standard prism along which
227  // new nodes should be generated. These edges are the three between
228  // the two triangular faces.
229  //
230  // edgeVertMap specifies the vertices which comprise those edges in
231  // splitEdge; for example splitEdge[0] = 3 which connects vertices 0
232  // and 3.
233  //
234  // edgeOffset holds the offset of each of edges 3, 1 and 8
235  // respectively inside the collapsed coordinate system.
236  map<LibUtilities::ShapeType, map<int, SplitEdgeHelper> > splitEdge;
237 
238  int splitPrismEdges[3] = {3, 1, 8};
239  int splitPrismEdgeVert[3][2] = {{0, 3}, {1, 2}, {4, 5}};
240  int splitPrismOffset[3] = {0, nq - 1, nq * (nl + 1) * (nq - 1)};
241  int splitPrismInc[3] = {nq, nq, nq};
242  SplitEdgeHelper splitPrismEdge;
243  splitPrismEdge.size = 3;
244  splitPrismEdge.edge = splitPrismEdges;
245  splitPrismEdge.edgeVert = helper2d(3, splitPrismEdgeVert);
246  splitPrismEdge.offset = splitPrismOffset;
247  splitPrismEdge.inc = splitPrismInc;
248  splitEdge[LibUtilities::ePrism][1] = splitPrismEdge;
249  splitEdge[LibUtilities::ePrism][3] = splitPrismEdge;
250 
251  int splitHex0Edges[4] = {4, 5, 6, 7};
252  int splitHex0EdgeVert[4][2] = {{0, 4}, {1, 5}, {2, 6}, {3, 7}};
253  int splitHex0Offset[4] = {0, nq - 1, nq * nq - 1, nq * (nq - 1)};
254  int splitHex0Inc[4] = {nq * nq, nq * nq, nq * nq, nq * nq};
255  SplitEdgeHelper splitHex0Edge;
256  splitHex0Edge.size = 4;
257  splitHex0Edge.edge = splitHex0Edges;
258  splitHex0Edge.edgeVert = helper2d(4, splitHex0EdgeVert);
259  splitHex0Edge.offset = splitHex0Offset;
260  splitHex0Edge.inc = splitHex0Inc;
261  splitEdge[LibUtilities::eHexahedron][0] = splitHex0Edge;
262  splitEdge[LibUtilities::eHexahedron][5] = splitHex0Edge;
263 
264  map<LibUtilities::ShapeType, map<int, bool> > revPoints;
265  revPoints[LibUtilities::ePrism][1] = true;
266  revPoints[LibUtilities::ePrism][3] = false;
267 
268  revPoints[LibUtilities::eHexahedron][0] = true;
269  revPoints[LibUtilities::eHexahedron][5] = false;
270 
271  // edgeMap associates geometry edge IDs to the (nl+1) vertices which
272  // are generated along that edge when a prism is split, and is used
273  // to avoid generation of duplicate vertices. It is stored as an
274  // unordered map for speed.
275  boost::unordered_map<int, vector<NodeSharedPtr> > edgeMap;
276  boost::unordered_map<int, vector<NodeSharedPtr> >::iterator eIt;
277 
278  string surf = m_config["surf"].as<string>();
279  if (surf.size() > 0)
280  {
281  vector<unsigned int> surfs;
282  ParseUtils::GenerateSeqVector(surf.c_str(), surfs);
283  sort(surfs.begin(), surfs.end());
284 
285  // If surface is defined, process list of elements to find those
286  // that are connected to it.
287  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
288  {
289  ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
290  int nSurf = el->GetFaceCount();
291 
292  for (int j = 0; j < nSurf; ++j)
293  {
294  int bl = el->GetBoundaryLink(j);
295  if (bl == -1)
296  {
297  continue;
298  }
299 
300  ElementSharedPtr bEl =
301  m_mesh->m_element[m_mesh->m_expDim - 1][bl];
302  vector<int> tags = bEl->GetTagList();
303  vector<int> inter;
304 
305  sort(tags.begin(), tags.end());
306  set_intersection(surfs.begin(),
307  surfs.end(),
308  tags.begin(),
309  tags.end(),
310  back_inserter(inter));
311  ASSERTL0(inter.size() <= 1, "Intersection of surfaces wrong");
312 
313  if (inter.size() == 1)
314  {
315  if (el->GetConf().m_e == LibUtilities::ePrism)
316  {
317  if (j % 2 == 0)
318  {
319  cerr << "WARNING: Found quadrilateral face " << j
320  << " on surface " << surf
321  << " connected to prism; ignoring." << endl;
322  continue;
323  }
324 
325  if (splitEls.count(el->GetId()) > 0)
326  {
327  cerr << "WARNING: prism already found; "
328  << "ignoring" << endl;
329  }
330 
331  splitEls[el->GetId()] = j;
332  }
333  else if (validElTypes.count(el->GetConf().m_e) == 0)
334  {
335  cerr << "WARNING: Unsupported element type "
336  << "found in surface " << j << "; "
337  << "ignoring" << endl;
338  continue;
339  }
340  }
341  }
342  }
343  }
344  else
345  {
346  // Otherwise, add all prismatic elements and assume face 1 of
347  // the prism lies on the surface.
348  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
349  {
350  ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
351 
352  if (el->GetConf().m_e == LibUtilities::ePrism)
353  {
354  splitEls[el->GetId()] = 1;
355  }
356  else if (validElTypes.count(el->GetConf().m_e) > 0)
357  {
358  splitEls[el->GetId()] = 0;
359  }
360  else
361  {
362  continue;
363  }
364  }
365  }
366 
367  if (splitEls.size() == 0)
368  {
369  cerr << "WARNING: No elements detected to split." << endl;
370  return;
371  }
372 
373  // Erase all elements from the element list. Elements will be
374  // re-added as they are split.
375  vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim];
376  m_mesh->m_element[m_mesh->m_expDim].clear();
377 
378  map<int, SpatialDomains::Geometry3DSharedPtr> geomMap;
379  for (int i = 0; i < el.size(); ++i)
380  {
381  const int elId = el[i]->GetId();
382  sIt = splitEls.find(elId);
383  if (sIt == splitEls.end())
384  {
385  continue;
386  }
387 
388  // Get elemental geometry object and put into map.
389  geomMap[elId] = boost::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
390  el[i]->GetGeom(m_mesh->m_spaceDim));
391  }
392 
393  // Iterate over list of elements of expansion dimension.
394  for (int i = 0; i < el.size(); ++i)
395  {
396  const int elId = el[i]->GetId();
397  sIt = splitEls.find(elId);
398 
399  if (sIt == splitEls.end())
400  {
401  m_mesh->m_element[m_mesh->m_expDim].push_back(el[i]);
402  continue;
403  }
404 
405  SpatialDomains::Geometry3DSharedPtr geom = geomMap[elId];
406 
407  const int faceNum = sIt->second;
408  LibUtilities::ShapeType elType = el[i]->GetConf().m_e;
409 
410  SplitMapHelper &sMap = splitMap[elType][faceNum];
411  SplitEdgeHelper &sEdge = splitEdge[elType][faceNum];
412 
413  // Find quadrilateral boundary faces if any
414  std::map<int, int> bLink;
415  for (int j = 0; j < sMap.bfacesSize; ++j)
416  {
417  int bl = el[i]->GetBoundaryLink(sMap.bfaces[j]);
418  if (bl != -1)
419  {
420  bLink[sMap.bfaces[j]] = bl;
421  }
422  }
423 
424  // Determine whether to use reverse points.
426  revPoints[elType][faceNum] ? LibUtilities::eBoundaryLayerPoints
428 
429  // Determine value of r based on geometry.
430  if (ratioIsString)
431  {
432  NekDouble x, y, z;
433  NekDouble x1, y1, z1;
434  int nverts = geom->GetNumVerts();
435 
436  x = y = z = 0.0;
437 
438  for (int i = 0; i < nverts; ++i)
439  {
440  geom->GetVertex(i)->GetCoords(x1, y1, z1);
441  x += x1;
442  y += y1;
443  z += z1;
444  }
445  x /= (NekDouble)nverts;
446  y /= (NekDouble)nverts;
447  z /= (NekDouble)nverts;
448  r = rEval.Evaluate(rExprId, x, y, z, 0.0);
449  }
450 
452 
453  if (elType == LibUtilities::ePrism)
454  {
455  // Create basis.
456  LibUtilities::BasisKey B0(
457  LibUtilities::eModified_A, nq, LibUtilities::PointsKey(nq, pt));
458  LibUtilities::BasisKey B1(LibUtilities::eModified_A,
459  2,
460  LibUtilities::PointsKey(nl + 1, t, r));
461  LibUtilities::BasisKey B2(
462  LibUtilities::eModified_B, nq, LibUtilities::PointsKey(nq, pt));
463 
464  // Create local region.
466  boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(geom);
468  B0, B1, B2, g);
469  }
470  else if (elType == LibUtilities::eHexahedron)
471  {
472  // Create basis.
473  LibUtilities::BasisKey B0(
474  LibUtilities::eModified_A, nq, LibUtilities::PointsKey(nq, pt));
475  LibUtilities::BasisKey B1(LibUtilities::eModified_A,
476  2,
477  LibUtilities::PointsKey(nl + 1, t, r));
478 
479  // Create local region.
481  boost::dynamic_pointer_cast<SpatialDomains::HexGeom>(geom);
483  B0, B0, B1, g);
484  }
485 
486  // Grab co-ordinates.
487  Array<OneD, NekDouble> x(nq * nq * (nl + 1));
488  Array<OneD, NekDouble> y(nq * nq * (nl + 1));
489  Array<OneD, NekDouble> z(nq * nq * (nl + 1));
490  q->GetCoords(x, y, z);
491 
492  int nSplitEdge = sEdge.size;
493  vector<vector<NodeSharedPtr> > edgeNodes(nSplitEdge);
494 
495  // Loop over edges to be split.
496  for (int j = 0; j < nSplitEdge; ++j)
497  {
498  int locEdge = sEdge.edge[j];
499  int edgeId = el[i]->GetEdge(locEdge)->m_id;
500 
501  // Determine whether we have already generated vertices
502  // along this edge.
503  eIt = edgeMap.find(edgeId);
504 
505  if (eIt == edgeMap.end())
506  {
507  // If not then resize storage to hold new points.
508  edgeNodes[j].resize(nl + 1);
509 
510  // Re-use existing vertices at endpoints of edge to
511  // avoid duplicating the existing vertices.
512  edgeNodes[j][0] = el[i]->GetVertex(sEdge.edgeVert[j][0]);
513  edgeNodes[j][nl] = el[i]->GetVertex(sEdge.edgeVert[j][1]);
514 
515  // Variable geometric ratio
516  if (ratioIsString)
517  {
518  NekDouble x0, y0, z0;
519  NekDouble x1, y1, z1;
520  NekDouble xm, ym, zm;
521 
522  // -> Find edge end and mid points
523  x0 = x[sEdge.offset[j]];
524  y0 = y[sEdge.offset[j]];
525  z0 = z[sEdge.offset[j]];
526 
527  x1 = x[sEdge.offset[j] + nl * nq];
528  y1 = y[sEdge.offset[j] + nl * nq];
529  z1 = z[sEdge.offset[j] + nl * nq];
530 
531  xm = 0.5 * (x0 + x1);
532  ym = 0.5 * (y0 + y1);
533  zm = 0.5 * (z0 + z1);
534 
535  // evaluate r factor based on mid point value
536  NekDouble rnew;
537  rnew = rEval.Evaluate(rExprId, xm, ym, zm, 0.0);
538 
539  // Get basis with new r;
540  LibUtilities::PointsKey Pkey(nl + 1, t, rnew);
543 
544  const Array<OneD, const NekDouble> z = newP->GetZ();
545 
546  // Create new interior nodes based on this new blend
547  for (int k = 1; k < nl; ++k)
548  {
549  xm = 0.5 * (1 + z[k]) * (x1 - x0) + x0;
550  ym = 0.5 * (1 + z[k]) * (y1 - y0) + y0;
551  zm = 0.5 * (1 + z[k]) * (z1 - z0) + z0;
552  edgeNodes[j][k] =
553  NodeSharedPtr(new Node(nodeId++, xm, ym, zm));
554  }
555  }
556  else
557  {
558  // Create new interior nodes.
559  for (int k = 1; k < nl; ++k)
560  {
561  int pos = sEdge.offset[j] + k * sEdge.inc[j];
562  edgeNodes[j][k] = NodeSharedPtr(
563  new Node(nodeId++, x[pos], y[pos], z[pos]));
564  }
565  }
566 
567  // Store these edges in edgeMap.
568  edgeMap[edgeId] = edgeNodes[j];
569  }
570  else
571  {
572  edgeNodes[j] = eIt->second;
573  }
574  }
575 
576  // Create element layers.
577  for (int j = 0; j < nl; ++j)
578  {
579  // Offset of this layer within the collapsed coordinate
580  // system.
581  int offset = j * sMap.layerOff;
582 
583  // Get corner vertices.
584  vector<NodeSharedPtr> nodeList(sMap.size);
585  for (int k = 0; k < sMap.size; ++k)
586  {
587  nodeList[k] = edgeNodes[sMap.conn[k][0]][j + sMap.conn[k][1]];
588  }
589 
590  // Create the element.
591  ElmtConfig conf(elType, 1, true, true, false);
593  elType, conf, nodeList, el[i]->GetTagList());
594 
595  // Add high order nodes to split prismatic edges.
596  for (int l = 0; l < sMap.size; ++l)
597  {
598  EdgeSharedPtr HOedge = elmt->GetEdge(sMap.edge[l]);
599  for (int k = 1; k < nq - 1; ++k)
600  {
601  int pos = offset + sMap.offset[l] + k * sMap.inc[l];
602  HOedge->m_edgeNodes.push_back(NodeSharedPtr(
603  new Node(nodeId++, x[pos], y[pos], z[pos])));
604  }
605  HOedge->m_curveType = pt;
606  }
607 
608  // Change the surface elements to match the layers of
609  // elements on the boundary of the domain.
611  for (it = bLink.begin(); it != bLink.end(); ++it)
612  {
613  int fid = it->first;
614  int bl = it->second;
615 
616  vector<NodeSharedPtr> qNodeList(4);
617  for (int k = 0; k < 4; ++k)
618  {
619  qNodeList[k] = nodeList[faceNodeMap[elType][fid][k]];
620  }
621  vector<int> tagBE;
622  tagBE =
623  m_mesh->m_element[m_mesh->m_expDim - 1][bl]->GetTagList();
624  ElmtConfig bconf(
625  LibUtilities::eQuadrilateral, 1, true, true, false);
626  ElementSharedPtr boundaryElmt =
628  LibUtilities::eQuadrilateral, bconf, qNodeList, tagBE);
629 
630  // Overwrite first layer boundary element with new
631  // boundary element, otherwise push this back to end of
632  // the boundary list
633  if (j == 0)
634  {
635  m_mesh->m_element[m_mesh->m_expDim - 1][bl] = boundaryElmt;
636  }
637  else
638  {
639  m_mesh->m_element[m_mesh->m_expDim - 1].push_back(
640  boundaryElmt);
641  }
642  }
643 
644  m_mesh->m_element[m_mesh->m_expDim].push_back(elmt);
645  }
646  }
647 
648  // Re-process mesh to eliminate duplicate vertices and edges.
649  ProcessVertices();
650  ProcessEdges();
651  ProcessFaces();
652  ProcessElements();
654 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Basic information about an element.
Definition: Element.h:58
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.
map< string, ConfigOption > m_config
List of configuration values.
Principle Modified Functions .
Definition: BasisType.h:49
MeshSharedPtr m_mesh
Mesh object.
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:61
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:79
Represents a point in the domain.
Definition: Node.h:60
virtual void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
virtual void ProcessVertices()
Extract element vertices.
virtual void ProcessElements()
Generate element IDs.
Principle Modified Functions .
Definition: BasisType.h:50
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
PointsManagerT & PointsManager(void)
double NekDouble
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
virtual void ProcessComposites()
Generate composites.
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:196
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:67
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:109
virtual void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
1D power law distribution for boundary layer points
Definition: PointsType.h:66
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:52
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
boost::shared_ptr< Geometry3D > Geometry3DSharedPtr
Definition: Geometry3D.h:52

Member Data Documentation

ModuleKey Nektar::Utilities::ProcessBL::className
static
Initial value:
"Refines a prismatic boundary layer.")

Definition at line 60 of file ProcessBL.h.