Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NekMesh/Module.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Module.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: Abstract input/output modules.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include "Module.h"
37 
38 using namespace std;
39 
40 namespace Nektar
41 {
42 namespace Utilities
43 {
44 
45 /**
46  * Returns an instance of the module factory, held as a singleton.
47  */
49 {
50  typedef Loki::SingletonHolder<ModuleFactory,
51  Loki::CreateUsingNew,
52  Loki::NoDestroy,
53  Loki::SingleThreaded> Type;
54  return Type::Instance();
55 }
56 
57 /**
58  * Prints a given module key to a stream.
59  */
60 std::ostream& operator<<(std::ostream& os, const ModuleKey& rhs)
61 {
62  return os << ModuleTypeMap[rhs.first] << ": " << rhs.second;
63 }
64 
65 InputModule::InputModule(MeshSharedPtr m) : Module(m)
66 {
67  m_config["infile"] = ConfigOption(false, "", "Input filename.");
68 }
69 
71 {
72  m_config["outfile"] = ConfigOption(false, "", "Output filename.");
73 }
74 
75 /**
76  * @brief Open a file for input.
77  */
79 {
80  string fname = m_config["infile"].as<string>();
81  m_mshFile.open(fname.c_str());
82  if (!m_mshFile.good())
83  {
84  cerr << "Error opening file: " << fname << endl;
85  abort();
86  }
87 }
88 
89 /**
90  * @brief Open a file for output.
91  */
93 {
94  string fname = m_config["outfile"].as<string>();
95  m_mshFile.open(fname.c_str());
96  if (!m_mshFile.good())
97  {
98  cerr << "Error opening file: " << fname << endl;
99  abort();
100  }
101 }
102 
103 /**
104  * @brief Create a unique set of mesh vertices from elements stored in
105  * Mesh::element.
106  *
107  * Each element is processed in turn and the vertices extracted and
108  * inserted into #m_vertexSet, which at the end of the routine
109  * contains all unique vertices in the mesh.
110  */
112 {
113  vector<ElementSharedPtr> &elmt = m_mesh->m_element[m_mesh->m_expDim];
114 
115  m_mesh->m_vertexSet.clear();
116 
117  for (int i = 0, vid = 0; i < elmt.size(); ++i)
118  {
119  for (int j = 0; j < elmt[i]->GetVertexCount(); ++j)
120  {
121  pair<NodeSet::iterator,bool> testIns =
122  m_mesh->m_vertexSet.insert(elmt[i]->GetVertex(j));
123 
124  if (testIns.second)
125  {
126  (*(testIns.first))->m_id = vid++;
127  }
128  else
129  {
130  elmt[i]->SetVertex(j,*testIns.first);
131  }
132  }
133  }
134 }
135 
136 /**
137  * @brief Create a unique set of mesh edges from elements stored in
138  * Mesh::element.
139  *
140  * All elements are first scanned and a list of unique, enumerated
141  * edges produced in #m_edgeSet. Since each element generated its
142  * edges independently, we must now ensure that each element only uses
143  * edge objects from the #m_edgeSet set This ensures there are no
144  * duplicate edge objects. Finally, we scan the list of elements for
145  * 1-D boundary elements which correspond to an edge in
146  * #m_edgeSet. For such elements, we set its edgeLink to reference the
147  * corresponding edge in #m_edgeSet.
148  *
149  * This routine only proceeds if the expansion dimension is 2 or 3.
150  */
151 void Module::ProcessEdges(bool ReprocessEdges)
152 {
153  if (m_mesh->m_expDim < 2) return;
154 
155  if(ReprocessEdges)
156  {
157  vector<ElementSharedPtr> &elmt = m_mesh->m_element[m_mesh->m_expDim];
158 
159  m_mesh->m_edgeSet.clear();
160 
161  // Clear all edge links
162 
163  // Scan all elements and generate list of unique edges
164  for (int i = 0, eid = 0; i < elmt.size(); ++i)
165  {
166  for (int j = 0; j < elmt[i]->GetEdgeCount(); ++j)
167  {
168  pair<EdgeSet::iterator,bool> testIns;
169  EdgeSharedPtr ed = elmt[i]->GetEdge(j);
170  testIns = m_mesh->m_edgeSet.insert(ed);
171 
172  if (testIns.second)
173  {
174  EdgeSharedPtr ed2 = *testIns.first;
175  ed2->m_id = eid++;
176  ed2->m_elLink.push_back(
177  pair<ElementSharedPtr,int>(elmt[i],j));
178  }
179  else
180  {
181  EdgeSharedPtr e2 = *(testIns.first);
182  elmt[i]->SetEdge(j, e2);
183  if (e2->m_edgeNodes.size() == 0 &&
184  ed->m_edgeNodes.size() > 0)
185  {
186  e2->m_curveType = ed->m_curveType;
187  e2->m_edgeNodes = ed->m_edgeNodes;
188 
189  // Reverse nodes if appropriate.
190  if (e2->m_n1->m_id != ed->m_n1->m_id)
191  {
192  reverse(e2->m_edgeNodes.begin(),
193  e2->m_edgeNodes.end());
194  }
195  }
196 
197 #ifdef NEKTAR_USE_MESHGEN
198  if(ed->onCurve)
199  {
200  e2->onCurve = ed->onCurve;
201  e2->CADCurveId = ed->CADCurveId;
202  e2->CADCurve = ed->CADCurve;
203  }
204 #endif
205 
206  // Update edge to element map.
207  e2->m_elLink.push_back(
208  pair<ElementSharedPtr,int>(elmt[i],j));
209  }
210  }
211  }
212  }
213 
214  // Create links for 1D elements
215  for (int i = 0; i < m_mesh->m_element[1].size(); ++i)
216  {
217  NodeSharedPtr v0 = m_mesh->m_element[1][i]->GetVertex(0);
218  NodeSharedPtr v1 = m_mesh->m_element[1][i]->GetVertex(1);
219  vector<NodeSharedPtr> edgeNodes;
220  EdgeSharedPtr E = boost::shared_ptr<Edge>(
221  new Edge(v0, v1, edgeNodes,
222  m_mesh->m_element[1][i]->GetConf().m_edgeCurveType));
223  EdgeSet::iterator it = m_mesh->m_edgeSet.find(E);
224  if (it == m_mesh->m_edgeSet.end())
225  {
226  cerr << "Cannot find corresponding element edge for "
227  << "1D element " << i << endl;
228  abort();
229  }
230  m_mesh->m_element[1][i]->SetEdgeLink(*it);
231 
232  // Update 2D element boundary map.
233  pair<ElementSharedPtr, int> eMap = (*it)->m_elLink.at(0);
234  eMap.first->SetBoundaryLink(eMap.second, i);
235 
236  // Copy curvature to edge.
237  if ((*it)->m_edgeNodes.size() > 0)
238  {
239  ElementSharedPtr edge = m_mesh->m_element[1][i];
240  if (edge->GetVertex(0) == (*it)->m_n1)
241  {
242  edge->SetVolumeNodes((*it)->m_edgeNodes);
243  }
244  }
245  }
246 }
247 
248 
249 /**
250  * @brief Create a unique set of mesh faces from elements stored in
251  * Mesh::element.
252  *
253  * All elements are scanned and a unique list of enumerated faces is
254  * produced in #m_faceSet. Since elements created their own faces
255  * independently, we examine each element only uses face objects from
256  * #m_faceSet. Duplicate faces of those in #m_face are replaced with
257  * the corresponding entry in #m_faceSet. Finally, we scan the list of
258  * elements for 2-D boundary faces which correspond to faces in
259  * #m_faceSet. For such elements, we set its faceLink to reference the
260  * corresponding face in #m_faceSet.
261  *
262  * This routine only proceeds if the expansion dimension is 3.
263  */
264 void Module::ProcessFaces(bool ReprocessFaces)
265 {
266  if (m_mesh->m_expDim < 3) return;
267 
268  if(ReprocessFaces)
269  {
270  vector<ElementSharedPtr> &elmt = m_mesh->m_element[m_mesh->m_expDim];
271 
272  m_mesh->m_faceSet.clear();
273 
274  // Scan all elements and generate list of unique faces
275  for (int i = 0, fid = 0; i < elmt.size(); ++i)
276  {
277  for (int j = 0; j < elmt[i]->GetFaceCount(); ++j)
278  {
279  pair<FaceSet::iterator,bool> testIns;
280  testIns = m_mesh->m_faceSet.insert(elmt[i]->GetFace(j));
281 
282  if (testIns.second)
283  {
284  (*(testIns.first))->m_id = fid++;
285  (*(testIns.first))->m_elLink.push_back(
286  pair<ElementSharedPtr,int>(elmt[i],j));
287  }
288  else
289  {
290  elmt[i]->SetFace(j,*testIns.first);
291  // Update face to element map.
292  (*(testIns.first))->m_elLink.push_back(
293  pair<ElementSharedPtr,int>(elmt[i],j));
294  }
295  }
296  }
297  }
298 
299  // Create links for 2D elements
300  for (int i = 0; i < m_mesh->m_element[2].size(); ++i)
301  {
302  vector<NodeSharedPtr> vertices = m_mesh->m_element[2][i]->GetVertexList();
303  vector<NodeSharedPtr> faceNodes;
304  vector<EdgeSharedPtr> edgeList = m_mesh->m_element[2][i]->GetEdgeList();
305  FaceSharedPtr F = boost::shared_ptr<Face>(
306  new Face(vertices, faceNodes, edgeList,
307  m_mesh->m_element[2][i]->GetConf().m_faceCurveType));
308  FaceSet::iterator it = m_mesh->m_faceSet.find(F);
309  if (it == m_mesh->m_faceSet.end())
310  {
311  cout << "Cannot find corresponding element face for 2D "
312  << "element " << i << endl;
313  abort();
314  }
315  m_mesh->m_element[2][i]->SetFaceLink(*it);
316 
317  // Update 3D element boundary map.
318  pair<ElementSharedPtr, int> eMap = (*it)->m_elLink.at(0);
319  eMap.first->SetBoundaryLink(eMap.second, i);
320  }
321 }
322 
323 /**
324  * @brief Enumerate elements stored in Mesh::element.
325  *
326  * For all elements of equal dimension to the mesh dimension, we
327  * enumerate sequentially. All other elements in the list should be of
328  * lower dimension and have ID set by a corresponding edgeLink or
329  * faceLink (as set in #ProcessEdges or #ProcessFaces).
330  */
332 {
333  int cnt = 0;
334  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
335  {
336  m_mesh->m_element[m_mesh->m_expDim][i]->SetId(cnt++);
337  }
338 }
339 
340 /**
341  * @brief Generate a list of composites (groups of elements) from tag
342  * IDs stored in mesh vertices/edges/faces/elements.
343  *
344  * Each element is assigned to a composite ID by an input module. First
345  * we scan the element list and generate a list of composite IDs. We
346  * then generate the composite objects and populate them with a second
347  * scan through the element list.
348  */
350 {
351  m_mesh->m_composite.clear();
352 
353  // For each element, check to see if a composite has been
354  // created. If not, create a new composite. Otherwise, add the
355  // element to the composite.
356  for (int d = 0; d <= m_mesh->m_expDim; ++d)
357  {
358  vector<ElementSharedPtr> &elmt = m_mesh->m_element[d];
359 
360  for (int i = 0; i < elmt.size(); ++i)
361  {
363  unsigned int tagid = elmt[i]->GetTagList()[0];
364 
365  it = m_mesh->m_composite.find(tagid);
366 
367  if (it == m_mesh->m_composite.end())
368  {
369  CompositeSharedPtr tmp = boost::shared_ptr<Composite>(
370  new Composite());
371  pair<CompositeMap::iterator, bool> testIns;
372  tmp->m_id = tagid;
373  tmp->m_tag = elmt[i]->GetTag();
374  if(m_mesh->m_faceLabels.count(tmp->m_id) != 0)
375  {
376  tmp->m_label = m_mesh->m_faceLabels[tmp->m_id];
377  }
378 
379  testIns = m_mesh->m_composite.insert(
380  pair<unsigned int, CompositeSharedPtr>(tagid,tmp));
381  it = testIns.first;
382  }
383 
384  if (elmt[i]->GetTag() != it->second->m_tag)
385  {
386  cout << "Different types of elements in same composite!" << endl;
387  cout << " -> Composite uses " << it->second->m_tag << endl;
388  cout << " -> Element uses " << elmt[i]->GetTag() << endl;
389  cout << "Have you specified physical volumes and surfaces?" << endl;
390  }
391  it->second->m_items.push_back(elmt[i]);
392  }
393  }
394 }
395 
396 /**
397  * clear all element link information from mesh entities to be able to reprocess new mesh
398  */
400 {
401  EdgeSet::iterator eit;
402 
403  for(eit = m_mesh->m_edgeSet.begin(); eit != m_mesh->m_edgeSet.end(); eit++)
404  {
405  (*eit)->m_elLink.clear();
406  }
407 
408  FaceSet::iterator fit;
409 
410  for(fit = m_mesh->m_faceSet.begin(); fit != m_mesh->m_faceSet.end(); fit++)
411  {
412  (*fit)->m_elLink.clear();
413  }
414 }
415 
416 /**
417  * @brief Reorder node IDs so that prisms and tetrahedra are aligned
418  * correctly.
419  *
420  * Orientation of prism lines (i.e. a large prism which has been split
421  * into subprisms) cannot be guaranteed when elements are created
422  * one-by-one, or when periodic boundary conditions are used. This
423  * routine uses the following strategy:
424  *
425  * - Destroy the existing node numbering.
426  * - Detect a line of prisms using the PrismLines routine.
427  * - For each line, renumber node IDs consistently so that highest ID
428  * per-element corresponds to the line of collapsed coordinate
429  * points.
430  * - Recreate each prism in the line using the new ordering, and apply
431  * the existing OrientPrism routine to orient nodes accordingly.
432  * - When all prism lines are processed, recreate all tetrahedra using
433  * the existing orientation code.
434  * - Finally renumber any other nodes (i.e. those belonging to
435  * hexahedra).
436  *
437  * The last step is to eliminate duplicate edges/faces and reenumerate.
438  *
439  * NOTE: This routine does not copy face-interior high-order information
440  * yet!
441  */
443 {
444  // Loop over elements and extract any that are prisms.
445  int i, j, k;
446 
447  if (m_mesh->m_expDim < 3)
448  {
449  return;
450  }
451 
452  map<int, int> lines;
453  set<int> prismsDone, tetsDone;
454  PerMap::iterator pIt;
455 
456  // Compile list of prisms and tets.
457  for (i = 0; i < m_mesh->m_element[3].size(); ++i)
458  {
459  ElementSharedPtr el = m_mesh->m_element[3][i];
460 
461  if (el->GetConf().m_e == LibUtilities::ePrism)
462  {
463  prismsDone.insert(i);
464  }
465  else if (el->GetConf().m_e == LibUtilities::eTetrahedron)
466  {
467  tetsDone.insert(i);
468  }
469  }
470 
471  // Destroy existing node numbering.
473  for (it = m_mesh->m_vertexSet.begin(); it != m_mesh->m_vertexSet.end(); ++it)
474  {
475  (*it)->m_id = -1;
476  }
477 
478  // Counter for new node IDs.
479  int nodeId = 0;
480  int prismTris[2][3] = {{0,1,4}, {3,2,5}};
481 
482  // Warning flag for high-order curvature information.
483  bool warnCurvature = false;
484 
485  // facesDone tracks face IDs inside prisms which have already been
486  // aligned.
487  boost::unordered_set<int> facesDone;
489 
490  // Loop over prisms until we've found all lines of prisms.
491  while (prismsDone.size() > 0)
492  {
493  vector<ElementSharedPtr> line;
494 
495  // Call PrismLines to identify all prisms connected to
496  // prismDone.begin() and place them in line[].
497  PrismLines(*prismsDone.begin(), perFaces, prismsDone, line);
498 
499  // Loop over each prism, figure out which line of vertices
500  // contains the vertex with highest ID.
501  for (i = 0; i < line.size(); ++i)
502  {
503  // Copy tags and nodes from existing element.
504  vector<int> tags = line[i]->GetTagList();
505  vector<NodeSharedPtr> nodes = line[i]->GetVertexList();
506 
507  // See if either face of this prism has been renumbered
508  // already.
509  FaceSharedPtr f[2] = {
510  line[i]->GetFace(1), line[i]->GetFace(3)
511  };
512 
513  fIt[0] = facesDone.find(f[0]->m_id);
514  fIt[1] = facesDone.find(f[1]->m_id);
515 
516  // See if either of these faces is periodic. If it is, then
517  // assign ids accordingly.
518  for (j = 0; j < 2; ++j)
519  {
520  pIt = perFaces.find(f[j]->m_id);
521 
522  if (pIt == perFaces.end())
523  {
524  continue;
525  }
526 
527  fIt2 = facesDone.find(pIt->second.first->m_id);
528 
529  if (fIt[j] == facesDone.end() &&
530  fIt2 != facesDone.end())
531  {
532  fIt[j] = fIt2;
533  }
534  }
535 
536  if (fIt[0] != facesDone.end() &&
537  fIt[1] != facesDone.end())
538  {
539  // Should not be the case that both faces have already
540  // been renumbered.
541  ASSERTL0(false, "Renumbering error!");
542  }
543  else if (fIt[0] == facesDone.end() &&
544  fIt[1] == facesDone.end())
545  {
546  // Renumber both faces.
547  for (j = 0; j < 2; ++j)
548  {
549  for (k = 0; k < 3; ++k)
550  {
551  NodeSharedPtr n = nodes[prismTris[j][k]];
552  if (n->m_id == -1)
553  {
554  n->m_id = nodeId++;
555  }
556  }
557  }
558 
559  facesDone.insert(f[0]->m_id);
560  facesDone.insert(f[1]->m_id);
561  }
562  else
563  {
564  // Renumber face. t identifies the face not yet
565  // numbered, o identifies the other face.
566  int t = fIt[0] == facesDone.end() ? 0 : 1;
567  int o = (t+1) % 2;
568  ASSERTL1(fIt[o] != facesDone.end(),"Renumbering error");
569 
570  // Determine which of the three vertices on the 'other'
571  // face corresponds to the highest ID - this signifies
572  // the singular point of the line of prisms.
573  int tmp1[3] = {
574  nodes[prismTris[o][0]]->m_id,
575  nodes[prismTris[o][1]]->m_id,
576  nodes[prismTris[o][2]]->m_id
577  };
578  int tmp2[3] = {0,1,2};
579 
580  if (tmp1[0] > tmp1[1])
581  {
582  swap(tmp1[0], tmp1[1]);
583  swap(tmp2[0], tmp2[1]);
584  }
585 
586  if (tmp1[1] > tmp1[2])
587  {
588  swap(tmp1[1], tmp1[2]);
589  swap(tmp2[1], tmp2[2]);
590  }
591 
592  if (tmp1[0] > tmp1[2])
593  {
594  swap(tmp1[0], tmp1[2]);
595  swap(tmp2[0], tmp2[2]);
596  }
597 
598  // Renumber this face so that highest ID matches.
599  for (j = 0; j < 3; ++j)
600  {
601  NodeSharedPtr n = nodes[prismTris[t][tmp2[j]]];
602  if (n->m_id == -1)
603  {
604  n->m_id = nodeId++;
605  }
606  }
607 
608  facesDone.insert(f[t]->m_id);
609  }
610 
611  for (j = 0; j < 6; ++j)
612  {
613  ASSERTL1(nodes[j]->m_id != -1, "Renumbering error");
614  }
615 
616  // Recreate prism with the new ordering.
617  ElmtConfig conf(LibUtilities::ePrism, 1, false, false, true);
619  LibUtilities::ePrism, conf, nodes, tags);
620 
621  // Now transfer high-order information back into
622  // place. TODO: Face curvature.
623  for (j = 0; j < 9; ++j)
624  {
625  EdgeSharedPtr e1 = line[i]->GetEdge(j);
626  for (k = 0; k < 9; ++k)
627  {
628  EdgeSharedPtr e2 = el->GetEdge(k);
629  if (e1->m_n1 == e2->m_n1 && e1->m_n2 == e2->m_n2)
630  {
631  e2->m_edgeNodes = e1->m_edgeNodes;
632  }
633  else if (e1->m_n1 == e2->m_n1 && e1->m_n2 == e2->m_n2)
634  {
635  e2->m_edgeNodes = e1->m_edgeNodes;
636  std::reverse(e2->m_edgeNodes.begin(),
637  e2->m_edgeNodes.end());
638  }
639  }
640  }
641 
642  // Warn users that we're throwing away face curvature
643  if (!warnCurvature)
644  {
645  for (j = 0; j < 5; ++j)
646  {
647  if (line[i]->GetFace(j)->m_faceNodes.size() > 0)
648  {
649  warnCurvature = true;
650  break;
651  }
652  }
653  }
654 
655  // Replace old prism.
656  m_mesh->m_element[3][line[i]->GetId()] = el;
657  }
658  }
659 
660  if (warnCurvature)
661  {
662  cerr << "[ReorderPrisms] WARNING: Face curvature detected in "
663  << "some prisms; this will be ignored in further module "
664  << "evaluations."
665  << endl;
666  }
667 
668  // Loop over periodic faces, enumerate vertices.
669  for (pIt = perFaces.begin(); pIt != perFaces.end(); ++pIt)
670  {
671  FaceSharedPtr f2 = pIt->second.first;
672  FaceSharedPtr f1 = perFaces[f2->m_id].first;
673  vector<int> perVerts = pIt->second.second;
674  int nVerts = perVerts.size();
675 
676  // Number periodic vertices first.
677  for (j = 0; j < nVerts; ++j)
678  {
679  NodeSharedPtr n1 = f1->m_vertexList[j];
680  NodeSharedPtr n2 = f2->m_vertexList[perVerts[j]];
681 
682  if (n1->m_id == -1 && n2->m_id == -1)
683  {
684  n1->m_id = nodeId++;
685  n2->m_id = nodeId++;
686  }
687  else if (n1->m_id != -1 && n2->m_id != -1)
688  {
689  continue;
690  }
691  else
692  {
693  ASSERTL0(false, "Periodic face renumbering error");
694  }
695  }
696  }
697 
698  // Recreate tets.
699  set<int>::iterator it2;
700  for (it2 = tetsDone.begin(); it2 != tetsDone.end(); ++it2)
701  {
702  ElementSharedPtr el = m_mesh->m_element[3][*it2];
703  vector<NodeSharedPtr> nodes = el->GetVertexList();
704  vector<int> tags = el->GetTagList();
705 
706  for (i = 0; i < 4; ++i)
707  {
708  if (nodes[i]->m_id == -1)
709  {
710  nodes[i]->m_id = nodeId++;
711  }
712  }
713 
714  // Recreate tet.
715  ElmtConfig conf(LibUtilities::eTetrahedron, 1, false, false, true);
716  m_mesh->m_element[3][*it2] = GetElementFactory().CreateInstance(
717  LibUtilities::eTetrahedron, conf, nodes, tags);
718  }
719 
720  // Enumerate rest of vertices.
721  for (it = m_mesh->m_vertexSet.begin(); it != m_mesh->m_vertexSet.end(); ++it)
722  {
723  if ((*it)->m_id == -1)
724  {
725  (*it)->m_id = nodeId++;
726  }
727  }
728 
729  ProcessEdges ();
730  ProcessFaces ();
731  ProcessElements();
732 }
733 
734 void Module::PrismLines(int prism,
735  PerMap &perFaces,
736  set<int> &prismsDone,
737  vector<ElementSharedPtr> &line)
738 {
739  int i;
740  set<int>::iterator it = prismsDone.find(prism);
741  PerMap::iterator it2;
742 
743  if (it == prismsDone.end())
744  {
745  return;
746  }
747 
748  // Remove this prism from the list.
749  prismsDone.erase(it);
750  line.push_back(m_mesh->m_element[3][prism]);
751 
752  // Now find prisms connected to this one through a triangular face.
753  for (i = 1; i <= 3; i += 2)
754  {
755  FaceSharedPtr f = m_mesh->m_element[3][prism]->GetFace(i);
756  int nextId;
757 
758  // See if this face is periodic.
759  it2 = perFaces.find(f->m_id);
760 
761  if (it2 != perFaces.end())
762  {
763  int id2 = it2->second.first->m_id;
764  nextId = it2->second.first->m_elLink[0].first->GetId();
765  perFaces.erase(it2);
766  perFaces.erase(id2);
767  PrismLines(nextId, perFaces, prismsDone, line);
768  }
769 
770  // Nothing else connected to this face.
771  if (f->m_elLink.size() == 1)
772  {
773  continue;
774  }
775 
776  nextId = f->m_elLink[0].first->GetId();
777  if (nextId == m_mesh->m_element[3][prism]->GetId())
778  {
779  nextId = f->m_elLink[1].first->GetId();
780  }
781 
782  PrismLines(nextId, perFaces, prismsDone, line);
783  }
784 }
785 
786 /**
787  * @brief Register a configuration option with a module.
788  */
789 void Module::RegisterConfig(string key, string val)
790 {
792  if (it == m_config.end())
793  {
794  cerr << "WARNING: Unrecognised config option " << key
795  << ", proceeding anyway." << endl;
796  }
797 
798  it->second.beenSet = true;
799 
800  if (it->second.isBool)
801  {
802  it->second.value = "1";
803  }
804  else
805  {
806  it->second.value = val;
807  }
808 }
809 
810 /**
811  * @brief Print out all configuration options for a module.
812  */
813 void Module::PrintConfig()
814 {
816 
817  if (m_config.size() == 0)
818  {
819  cerr << "No configuration options for this module." << endl;
820  return;
821  }
822 
823  for (it = m_config.begin(); it != m_config.end(); ++it)
824  {
825  cerr << setw(10) << it->first << ": " << it->second.desc
826  << endl;
827  }
828 }
829 
830 /**
831  * @brief Sets default configuration options for those which have not
832  * been set.
833  */
834 void Module::SetDefaults()
835 {
837 
838  for (it = m_config.begin(); it != m_config.end(); ++it)
839  {
840  if (!it->second.beenSet)
841  {
842  it->second.value = it->second.defValue;
843  }
844  }
845 }
846 
847 /**
848  * @brief Print a brief summary of information.
849  */
851 {
852  // Compute the number of full-dimensional elements and boundary
853  // elements.
854  cerr << "Expansion dimension is " << m_mesh->m_expDim << endl;
855  cerr << "Space dimension is " << m_mesh->m_spaceDim << endl;
856  cerr << "Read " << m_mesh->m_node.size() << " nodes" << endl;
857  cerr << "Read " << m_mesh->GetNumElements() << " "
858  << m_mesh->m_expDim << "-D elements" << endl;
859  cerr << "Read " << m_mesh->GetNumBndryElements()
860  << " boundary elements" << endl;
861 }
862 
863 }
864 }
LibUtilities::NekFactory< ModuleKey, Module, FieldSharedPtr > ModuleFactory
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Basic information about an element.
Definition: Element.h:58
std::ifstream m_mshFile
Input stream.
pair< ModuleType, string > ModuleKey
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
Represents an edge which joins two points.
Definition: Edge.h:61
Represents a face comprised of three or more edges.
Definition: Face.h:60
map< string, ConfigOption > m_config
List of configuration values.
STL namespace.
MeshSharedPtr m_mesh
Mesh object.
void RegisterConfig(string key, string value)
Register a configuration option with a module.
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
virtual void ClearElementLinks()
virtual void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
std::ostream & operator<<(std::ostream &os, const BasisKey &rhs)
Definition: Basis.cpp:85
virtual void ProcessVertices()
Extract element vertices.
virtual void ProcessElements()
Generate element IDs.
boost::shared_ptr< Composite > CompositeSharedPtr
Shared pointer to a composite.
Definition: Composite.h:121
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
virtual void ProcessComposites()
Generate composites.
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:114
void PrintSummary()
Print summary of elements.
void PrismLines(int prism, PerMap &perFaces, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)
void PrintConfig()
Print out all configuration options for a module.
Represents a command-line configuration option.
std::map< int, std::pair< FaceSharedPtr, std::vector< int > > > PerMap
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
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:137
virtual void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
void ReorderPrisms(PerMap &perFaces)
Reorder node IDs so that prisms and tetrahedra are aligned correctly.
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:52
std::ofstream m_mshFile
Output stream.
boost::shared_ptr< Face > FaceSharedPtr
Shared pointer to a face.
Definition: Face.h:378
void OpenStream()
Open a file for output.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:218
void SetDefaults()
Sets default configuration options for those which have not been set.
void OpenStream()
Open a file for input.
ModuleFactory & GetModuleFactory()
const char *const ModuleTypeMap[]