Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 (NekMeshUtils::MeshSharedPtr m)
 
virtual ~ProcessBL ()
 
void BoundaryLayer2D ()
 
void BoundaryLayer3D ()
 
virtual void Process ()
 Write mesh to output file. More...
 
- Public Member Functions inherited from Nektar::NekMeshUtils::ProcessModule
NEKMESHUTILS_EXPORT ProcessModule (MeshSharedPtr p_m)
 
- Public Member Functions inherited from Nektar::NekMeshUtils::Module
NEKMESHUTILS_EXPORT Module (MeshSharedPtr p_m)
 
NEKMESHUTILS_EXPORT void RegisterConfig (std::string key, std::string value)
 Register a configuration option with a module. More...
 
NEKMESHUTILS_EXPORT void PrintConfig ()
 Print out all configuration options for a module. More...
 
NEKMESHUTILS_EXPORT void SetDefaults ()
 Sets default configuration options for those which have not been set. More...
 
NEKMESHUTILS_EXPORT MeshSharedPtr GetMesh ()
 
virtual NEKMESHUTILS_EXPORT void ProcessVertices ()
 Extract element vertices. More...
 
virtual NEKMESHUTILS_EXPORT void ProcessEdges (bool ReprocessEdges=true)
 Extract element edges. More...
 
virtual NEKMESHUTILS_EXPORT void ProcessFaces (bool ReprocessFaces=true)
 Extract element faces. More...
 
virtual NEKMESHUTILS_EXPORT void ProcessElements ()
 Generate element IDs. More...
 
virtual NEKMESHUTILS_EXPORT void ProcessComposites ()
 Generate composites. More...
 
virtual NEKMESHUTILS_EXPORT void ClearElementLinks ()
 

Static Public Member Functions

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

Static Public Attributes

static NekMeshUtils::ModuleKey className
 

Additional Inherited Members

- Protected Member Functions inherited from Nektar::NekMeshUtils::Module
NEKMESHUTILS_EXPORT void ReorderPrisms (PerMap &perFaces)
 Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...
 
NEKMESHUTILS_EXPORT void PrismLines (int prism, PerMap &perFaces, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)
 
- Protected Attributes inherited from Nektar::NekMeshUtils::Module
MeshSharedPtr m_mesh
 Mesh object. More...
 
std::map< std::string,
ConfigOption
m_config
 List of configuration values. 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 ( NekMeshUtils::MeshSharedPtr  m)

Definition at line 109 of file ProcessBL.cpp.

References Nektar::NekMeshUtils::Module::m_config.

109  : 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 }
Represents a command-line configuration option.
NEKMESHUTILS_EXPORT ProcessModule(MeshSharedPtr p_m)
std::map< std::string, ConfigOption > m_config
List of configuration values.
Nektar::Utilities::ProcessBL::~ProcessBL ( )
virtual

Definition at line 122 of file ProcessBL.cpp.

123 {
124 }

Member Function Documentation

void Nektar::Utilities::ProcessBL::BoundaryLayer2D ( )

Definition at line 160 of file ProcessBL.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::eBoundaryLayerPoints, Nektar::LibUtilities::eBoundaryLayerPointsRev, Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eSegment, Nektar::ParseUtils::GenerateSeqVector(), Nektar::NekMeshUtils::GetElementFactory(), Nektar::iterator, Nektar::NekMeshUtils::Module::m_config, Nektar::NekMeshUtils::Module::m_mesh, class_topology::Node, and Nektar::LibUtilities::PointsManager().

Referenced by Process().

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.
167  LibUtilities::AnalyticExpressionEvaluator rEval;
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.
322  LibUtilities::BasisKey B0(
324  LibUtilities::PointsKey(nq,pt));
325  LibUtilities::BasisKey B1(
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 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Basic information about an element.
Definition: ElementConfig.h:50
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
boost::shared_ptr< Points< NekDouble > > PointsSharedPtr
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:79
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
PointsManagerT & PointsManager(void)
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< 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
1D power law distribution for boundary layer points
Definition: PointsType.h:68
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
void Nektar::Utilities::ProcessBL::BoundaryLayer3D ( )

Definition at line 584 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::NekMeshUtils::Module::m_config, Nektar::NekMeshUtils::Module::m_mesh, class_topology::Node, Nektar::Utilities::SplitMapHelper::offset, Nektar::Utilities::SplitEdgeHelper::offset, Nektar::LibUtilities::PointsManager(), Nektar::Utilities::SplitMapHelper::size, and Nektar::Utilities::SplitEdgeHelper::size.

Referenced by Process().

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.
596  LibUtilities::AnalyticExpressionEvaluator rEval;
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.
910  LibUtilities::BasisKey B0(
911  LibUtilities::eModified_A, nq, LibUtilities::PointsKey(nq, pt));
912  LibUtilities::BasisKey B1(LibUtilities::eModified_A,
913  2,
914  LibUtilities::PointsKey(nl + 1, t, r));
915  LibUtilities::BasisKey B2(
916  LibUtilities::eModified_B, nq, LibUtilities::PointsKey(nq, pt));
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.
927  LibUtilities::BasisKey B0(
928  LibUtilities::eModified_A, nq, LibUtilities::PointsKey(nq, pt));
929  LibUtilities::BasisKey B1(LibUtilities::eModified_A,
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 }
#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.
Principle Modified Functions .
Definition: BasisType.h:49
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
Principle Modified Functions .
Definition: BasisType.h:50
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
PointsManagerT & PointsManager(void)
double NekDouble
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
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< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:109
1D power law distribution for boundary layer points
Definition: PointsType.h:68
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
boost::shared_ptr< Geometry3D > Geometry3DSharedPtr
Definition: Geometry3D.h:52
static boost::shared_ptr<Module> Nektar::Utilities::ProcessBL::create ( NekMeshUtils::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::NekMeshUtils::Module.

Definition at line 126 of file ProcessBL.cpp.

References ASSERTL0, BoundaryLayer2D(), BoundaryLayer3D(), Nektar::NekMeshUtils::Module::m_mesh, Nektar::NekMeshUtils::Module::ProcessComposites(), Nektar::NekMeshUtils::Module::ProcessEdges(), Nektar::NekMeshUtils::Module::ProcessElements(), Nektar::NekMeshUtils::Module::ProcessFaces(), and Nektar::NekMeshUtils::Module::ProcessVertices().

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 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
virtual NEKMESHUTILS_EXPORT void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
virtual NEKMESHUTILS_EXPORT void ProcessElements()
Generate element IDs.
virtual NEKMESHUTILS_EXPORT void ProcessVertices()
Extract element vertices.
virtual NEKMESHUTILS_EXPORT void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
virtual NEKMESHUTILS_EXPORT void ProcessComposites()
Generate composites.

Member Data Documentation

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

Definition at line 60 of file ProcessBL.h.