Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
InputSem.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: InputSem.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: Semtex session converter.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
37 
38 #include "InputSem.h"
39 
40 using namespace std;
41 using namespace Nektar::NekMeshUtils;
42 
43 namespace Nektar
44 {
45 namespace Utilities
46 {
47 
48 ModuleKey InputSem::className =
50  InputSem::create,
51  "Reads Semtex session files.");
52 
53 /**
54  * @brief Initialises the InputSem class.
55  */
56 InputSem::InputSem(MeshSharedPtr m) : InputModule(m)
57 {
58 }
59 
61 {
62 }
63 
64 /**
65  * @brief Process a Semtex session file.
66  *
67  * Semtex files are defined by a tokenized markup format. We first
68  * populate #sectionMap which stores the location of the various tags in
69  * the session file so that they can be jumped to, since no ordering is
70  * defined. The converter only requires the NODES and ELEMENTS sections
71  * to exist, but can also read CURVES and SURFACES. High-order curves
72  * rely on the meshfile session.msh to be created with the Semtex
73  * utility meshpr first.
74  *
75  * @param pFilename Filename of Semtex session to read.
76  */
78 {
79  // Open the file stream.
80  OpenStream();
81 
82  if (m_mesh->m_verbose)
83  {
84  cout << "InputSem: Start reading file..." << endl;
85  }
86 
87  // Read through input file and populate the section map.
89  string line, word;
90  stringstream ss;
91  streampos linePos;
92 
93  sectionMap["NODES"] = -1;
94  sectionMap["ELEMENTS"] = -1;
95  sectionMap["CURVES"] = -1;
96  sectionMap["SURFACES"] = -1;
97  sectionMap["GROUPS"] = -1;
98  sectionMap["BCS"] = -1;
99  sectionMap["FIELDS"] = -1;
100 
101  while (!m_mshFile.eof())
102  {
103  linePos = m_mshFile.tellg();
104  getline(m_mshFile, line);
105  ss.clear();
106  ss.str(line);
107  ss >> word;
108 
109  // Iterate over all tokens and see if section exists on this
110  // line.
111  for (it = sectionMap.begin(); it != sectionMap.end(); ++it)
112  {
113  if (word == "<" + it->first || word == "<" + it->first + ">")
114  {
115  sectionMap[it->first] = linePos;
116  }
117  }
118  }
119 
120  // Clear eofbit and go back to the beginning of the file.
121  m_mshFile.clear();
122  m_mshFile.seekg(0);
123 
124  // Check that required sections exist in the file.
125  if (sectionMap["NODES"] == std::streampos(-1))
126  {
127  cerr << "Unable to locate NODES section in session file." << endl;
128  abort();
129  }
130 
131  if (sectionMap["ELEMENTS"] == std::streampos(-1))
132  {
133  cerr << "Unable to locate ELEMENTS section in session file." << endl;
134  abort();
135  }
136 
137  if (sectionMap["SURFACES"] != std::streampos(-1))
138  {
139  if (sectionMap["BCS"] == std::streampos(-1))
140  {
141  cerr << "SURFACES section defined but BCS section not "
142  << "found." << endl;
143  abort();
144  }
145 
146  if (sectionMap["GROUPS"] == std::streampos(-1))
147  {
148  cerr << "SURFACES section defined but GROUPS section not "
149  << "found." << endl;
150  abort();
151  }
152 
153  if (sectionMap["FIELDS"] == std::streampos(-1))
154  {
155  cerr << "SURFACES section defined but FIELDS section not "
156  << "found." << endl;
157  abort();
158  }
159  }
160 
161  m_mesh->m_expDim = 0;
162  string tag;
163  int start, end, nVertices, nEntities, nCurves, nSurf, nGroups, nBCs;
164  int id, i, j, k;
165  vector<double> hoXData, hoYData;
167  ifstream homeshFile;
168 
169  // Begin by reading in list of nodes which define the linear
170  // elements.
171  m_mshFile.seekg(sectionMap["NODES"]);
172  getline(m_mshFile, line);
173  ss.clear();
174  ss.str(line);
175  ss >> word;
176 
177  tag = ss.str();
178  start = tag.find_first_of('=');
179  end = tag.find_first_of('>');
180  nVertices = atoi(tag.substr(start + 1, end).c_str());
181 
182  i = id = 0;
183  while (i < nVertices)
184  {
185  getline(m_mshFile, line);
186  if (line.length() < 7)
187  continue;
188  ss.clear();
189  ss.str(line);
190  double x = 0, y = 0, z = 0;
191  ss >> id >> x >> y >> z;
192 
193  if ((y * y) > 0.000001 && m_mesh->m_spaceDim != 3)
194  {
195  m_mesh->m_spaceDim = 2;
196  }
197  if ((z * z) > 0.000001)
198  {
199  m_mesh->m_spaceDim = 3;
200  }
201  id -= 1; // counter starts at 0
202  m_mesh->m_node.push_back(
203  boost::shared_ptr<Node>(new Node(id, x, y, z)));
204  ++i;
205  }
206 
207  // Now read in elements
208  m_mshFile.seekg(sectionMap["ELEMENTS"]);
209  getline(m_mshFile, line);
210  ss.clear();
211  ss.str(line);
212  ss >> word;
213 
214  tag = ss.str();
215  start = tag.find_first_of('=');
216  end = tag.find_first_of('>');
217  nEntities = atoi(tag.substr(start + 1, end).c_str());
218 
219  i = id = 0;
220  while (i < nEntities)
221  {
222  getline(m_mshFile, line);
223  if (line.length() < 18)
224  {
225  continue;
226  }
227 
228  // Create element tags
229  vector<int> tags;
230  tags.push_back(0); // composite
231 
232  // Read element node list
233  ss.clear();
234  ss.str(line);
235  ss >> id >> word;
236  vector<NodeSharedPtr> nodeList;
237  for (j = 0; j < 4; ++j)
238  {
239  int node = 0;
240  ss >> node;
241  nodeList.push_back(m_mesh->m_node[node - 1]);
242  }
243 
244  // Create element
245  ElmtConfig conf(elType, 1, false, false);
246  ElementSharedPtr E =
247  GetElementFactory().CreateInstance(elType, conf, nodeList, tags);
248 
249  // Determine mesh expansion dimension
250  if (E->GetDim() > m_mesh->m_expDim)
251  {
252  m_mesh->m_expDim = E->GetDim();
253  }
254  m_mesh->m_element[E->GetDim()].push_back(E);
255  ++i;
256  }
257 
258  // Finally, process curves.
259  if (sectionMap["CURVES"] != std::streampos(-1))
260  {
261  int np, nel, nodeId = m_mesh->m_node.size();
262 
263  m_mshFile.seekg(sectionMap["CURVES"]);
264  getline(m_mshFile, line);
265  ss.clear();
266  ss.str(line);
267  ss >> word;
268 
269  tag = ss.str();
270  start = tag.find_first_of('=');
271  end = tag.find_first_of('>');
272  nCurves = atoi(tag.substr(start + 1, end).c_str());
273 
274  // Some session files have empty curves sections; if nCurves
275  // is 0, no nead to load high order mesh file.
276  if (nCurves > 0)
277  {
278  string fname = m_config["infile"].as<string>();
279  int ext = fname.find_last_of('.');
280  string meshfile = fname.substr(0, ext) + ".msh";
281 
282  homeshFile.open(meshfile.c_str());
283  if (!homeshFile.is_open())
284  {
285  cerr << "Cannot open or find mesh file: " << meshfile << endl
286  << "Make sure to run meshpr on your session "
287  << "file first." << endl;
288  abort();
289  }
290 
291  // Make sure we have matching header.
292  getline(homeshFile, line);
293  ss.clear();
294  ss.str(line);
295  ss >> np >> nel >> nel >> nel;
296 
297  if (nel != m_mesh->m_element[m_mesh->m_expDim].size())
298  {
299  cerr << "Number of elements mismatch in mesh file." << endl;
300  abort();
301  }
302 
303  // Now read in all mesh data. This is horribly inefficient
304  // since not all elements are curved, but it is the
305  // easiest way of finding element data.
306  hoXData.resize(nel * np * np);
307  hoYData.resize(nel * np * np);
308 
309  for (j = 0; j < nel * np * np; ++j)
310  {
311  getline(homeshFile, line);
312  ss.clear();
313  ss.str(line);
314  ss >> hoXData[j] >> hoYData[j];
315  }
316 
317  homeshFile.close();
318  }
319 
320  i = id = 0;
321  while (i < nCurves)
322  {
323  getline(m_mshFile, line);
324  if (line.length() < 18)
325  {
326  continue;
327  }
328  int elmt = 0, side = 0;
329  ss.clear();
330  ss.str(line);
331  ss >> id >> elmt >> side >> word;
332  id--;
333  elmt--;
334 
335  vector<NodeSharedPtr> edgeNodes;
336 
337  if (word != "<SPLINE>" && word != "<ARC>")
338  {
339  cerr << "Unknown curve tag: " << word << endl;
340  abort();
341  }
342 
343  // See if we have already retrieved high-order data
344  // for this elements; prevents unnecessary computation
345  // for elements with multiple curves.
346  if (m_mesh->m_element[2][elmt]->GetConf().m_order > 1)
347  {
348  ++i;
349  continue;
350  }
351 
352  // Now set high order data for requested element.
353  for (int side = 0; side < 4; ++side)
354  {
355  int offset = elmt * np * np;
356  int stride = 0;
357 
358  switch (side)
359  {
360  case 0: // Bottom edge
361  offset += 0;
362  stride = 1;
363  break;
364  case 1: // Right edge
365  offset += np - 1;
366  stride = np;
367  break;
368  case 2: // Top edge
369  offset += np * np - 1;
370  stride = -1;
371  break;
372  case 3: // Left edge
373  offset += np * (np - 1);
374  stride = -np;
375  break;
376  default:
377  cerr << "Unknown side for curve id " << id << endl;
378  abort();
379  }
380 
381  for (j = 1; j < np - 1; ++j, ++nodeId)
382  {
383  double x = hoXData[offset + j * stride];
384  double y = hoYData[offset + j * stride];
385  NodeSharedPtr n =
386  boost::shared_ptr<Node>(new Node(nodeId, x, y, 0.0));
387  edgeNodes.push_back(n);
388  }
389  }
390 
391  // Add internal points.
392  for (j = 1; j < np - 1; ++j)
393  {
394  int offset = j * np + 1;
395  for (k = 1; k < np - 1; ++k, ++nodeId)
396  {
397  double x = hoXData[offset + k];
398  double y = hoYData[offset + k];
399  NodeSharedPtr n =
400  boost::shared_ptr<Node>(new Node(nodeId, x, y, 0.0));
401  edgeNodes.push_back(n);
402  }
403  }
404 
405  // Grab existing element from list and retrieve tags and
406  // vertices; insert these into existing edge nodes.
407  ElementSharedPtr e = m_mesh->m_element[2][elmt];
408  vector<NodeSharedPtr> elvert = e->GetVertexList();
409  vector<int> tags = e->GetTagList();
410  edgeNodes.insert(edgeNodes.begin(), elvert.begin(), elvert.end());
411 
412  // Create new element and replace with an incomplete
413  // quadrilateral of the correct order.
414  ElmtConfig conf(elType,
415  np - 1,
416  true,
417  false,
418  true,
420  m_mesh->m_element[2][elmt] = GetElementFactory().CreateInstance(
421  elType, conf, edgeNodes, tags);
422 
423  ++i;
424  }
425  }
426 
427  // Process field names
428  if (sectionMap["FIELDS"] != std::streampos(-1))
429  {
430  m_mshFile.seekg(sectionMap["FIELDS"]);
431  getline(m_mshFile, line);
432  getline(m_mshFile, line);
433  ss.clear();
434  ss.str(line);
435 
436  while (ss >> tag)
437  {
438  m_mesh->m_fields.push_back(tag);
439  }
440  }
441 
442  // Process surfaces if they exist. This is deliberately done after
443  // curves to ensure high-order points are preserved.
444  if (sectionMap["SURFACES"] != std::streampos(-1))
445  {
446  map<string, int> conditionMap;
447  int maxTag = -1;
448 
449  // First read in list of groups, which defines each condition tag.
450  m_mshFile.seekg(sectionMap["GROUPS"]);
451  getline(m_mshFile, line);
452  ss.clear();
453  ss.str(line);
454  ss >> word;
455 
456  tag = ss.str();
457  start = tag.find_first_of('=');
458  end = tag.find_first_of('>');
459  nGroups = atoi(tag.substr(start + 1, end).c_str());
460 
461  i = id = 0;
462  while (i < nGroups)
463  {
464  getline(m_mshFile, line);
465  ss.clear();
466  ss.str(line);
467  ss >> id >> tag;
468  conditionMap[tag] = i++;
469  }
470 
471  maxTag = i;
472 
473  // Now read in actual values for boundary conditions from BCS
474  // section.
475  m_mshFile.seekg(sectionMap["BCS"]);
476  getline(m_mshFile, line);
477  ss.clear();
478  ss.str(line);
479  ss >> word;
480 
481  tag = ss.str();
482  start = tag.find_first_of('=');
483  end = tag.find_first_of('>');
484  nBCs = atoi(tag.substr(start + 1, end).c_str());
485 
486  i = id = 0;
487  while (i < nBCs)
488  {
489  int nF;
490  string tmp;
492  getline(m_mshFile, line);
493  ss.clear();
494  ss.str(line);
495  ss >> id >> tag >> nF;
496 
497  p = ConditionSharedPtr(new Condition());
498  m_mesh->m_condition[conditionMap[tag]] = p;
499 
500  // Read boundary condition.
501  j = 0;
502  while (j < nF)
503  {
504  getline(m_mshFile, line);
505  ss.clear();
506  ss.str(line);
507  ss >> tmp;
508 
509  // First string should be condition type.
510  if (tmp == "<D>")
511  {
512  p->type.push_back(eDirichlet);
513  }
514  else if (tmp == "<N>")
515  {
516  p->type.push_back(eNeumann);
517  }
518  else if (tmp == "<H>")
519  {
520  p->type.push_back(eHOPCondition);
521  p->value.push_back("0");
522  p->field.push_back("p");
523  ++j;
524  continue;
525  }
526  else
527  {
528  cerr << "Unsupported boundary condition type " << tmp
529  << endl;
530  abort();
531  }
532 
533  // Second string should be field.
534  ss >> tmp;
535  p->field.push_back(tmp);
536 
537  // Third string should be equals sign.
538  ss >> tmp;
539  if (tmp != "=")
540  {
541  cerr << "Couldn't read boundary condition type " << tag
542  << endl;
543  abort();
544  }
545 
546  // Fourth string should be value. CAUTION: Assumes
547  // expression is defined without any spaces in it!
548  ss >> tmp;
549  p->value.push_back(tmp);
550 
551  ++j;
552  }
553 
554  // Finally set composite for condition. In this case, all
555  // composites will be lines so there is one set per
556  // composite.
557  p->m_composite.push_back(conditionMap[tag] + 1);
558 
559  ++i;
560  }
561 
562  // Finally read surface information.
563  m_mshFile.seekg(sectionMap["SURFACES"]);
564  getline(m_mshFile, line);
565  ss.clear();
566  ss.str(line);
567  ss >> word;
568 
569  tag = ss.str();
570  start = tag.find_first_of('=');
571  end = tag.find_first_of('>');
572  nSurf = atoi(tag.substr(start + 1, end).c_str());
573 
574  i = id = 0;
575  int elmt, side;
576  int periodicTagId = -1;
577 
578  set<pair<int, int> > visitedPeriodic;
579 
580  while (i < nSurf)
581  {
582  getline(m_mshFile, line);
583  ss.clear();
584  ss.str(line);
585  ss >> id >> elmt >> side >> word;
586  elmt--;
587  side--;
588 
589  if (word == "<P>")
590  {
591  // If this is the first periodic boundary condition
592  // encountered, then set up m_mesh->m_condition with two
593  // periodic conditions.
594  if (periodicTagId == -1)
595  {
596  periodicTagId = maxTag;
598  ConditionSharedPtr out =
600  for (j = 0; j < m_mesh->m_fields.size(); ++j)
601  {
602  in->type.push_back(ePeriodic);
603  out->type.push_back(ePeriodic);
604  in->field.push_back(m_mesh->m_fields[j]);
605  out->field.push_back(m_mesh->m_fields[j]);
606  in->value.push_back("[" + boost::lexical_cast<string>(
607  periodicTagId + 1) +
608  "]");
609  out->value.push_back(
610  "[" + boost::lexical_cast<string>(periodicTagId) +
611  "]");
612  }
613  in->m_composite.push_back(periodicTagId + 1);
614  out->m_composite.push_back(periodicTagId + 2);
615  m_mesh->m_condition[periodicTagId] = in;
616  m_mesh->m_condition[periodicTagId + 1] = out;
617  }
618 
619  int elmtB, sideB;
620 
621  ss >> elmtB >> sideB;
622  elmtB--;
623  sideB--;
624 
625  pair<int, int> c1(elmt, side);
626  pair<int, int> c2(elmtB, sideB);
627 
628  if (visitedPeriodic.count(c1) == 0 &&
629  visitedPeriodic.count(c2) == 0)
630  {
631  visitedPeriodic.insert(make_pair(elmtB, sideB));
632  visitedPeriodic.insert(make_pair(elmt, side));
633  insertEdge(elmt, side, periodicTagId + 1);
634  insertEdge(elmtB, sideB, periodicTagId + 2);
635  }
636  }
637  else if (word == "<B>")
638  {
639  ss >> tag;
640  insertEdge(elmt, side, conditionMap[tag] + 1);
641  }
642  else
643  {
644  cerr << "Unrecognised or unsupported tag " << word << endl;
645  abort();
646  }
647  ++i;
648  }
649  }
650 
651  PrintSummary();
652  m_mshFile.reset();
653 
654  // Process rest of mesh.
655  ProcessVertices();
656  ProcessEdges();
657  ProcessFaces();
658  ProcessElements();
660 }
661 
662 void InputSem::insertEdge(int elmt, int side, int tagId)
663 {
664  EdgeSharedPtr edge = m_mesh->m_element[2][elmt]->GetEdge(side);
665  vector<NodeSharedPtr> edgeNodes = edge->m_edgeNodes;
666  edgeNodes.insert(edgeNodes.begin(), edge->m_n2);
667  edgeNodes.insert(edgeNodes.begin(), edge->m_n1);
668  int order = edgeNodes.size() - 1;
669 
670  vector<int> tags;
671  tags.push_back(tagId);
672 
674  order,
675  order > 1,
676  false,
677  true,
680  LibUtilities::eSegment, conf, edgeNodes, tags);
681  m_mesh->m_element[1].push_back(E);
682 }
683 }
684 }
Defines a boundary condition.
Definition: Mesh.h:71
Basic information about an element.
Definition: ElementConfig.h:50
io::filtering_istream m_mshFile
Input stream.
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
virtual void Process()
Process a Semtex session file.
Definition: InputSem.cpp:77
STL namespace.
pair< ModuleType, string > ModuleKey
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
NEKMESHUTILS_EXPORT void PrintSummary()
Print summary of elements.
NEKMESHUTILS_EXPORT void OpenStream()
Open a file for input.
virtual NEKMESHUTILS_EXPORT void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
virtual NEKMESHUTILS_EXPORT void ProcessElements()
Generate element IDs.
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
std::map< std::string, std::streampos > sectionMap
Maps Semtex sections to positions inside the input file.
Definition: InputSem.h:68
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< Condition > ConditionSharedPtr
Definition: Mesh.h:82
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:136
Abstract base class for input modules.
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:147
virtual NEKMESHUTILS_EXPORT void ProcessVertices()
Extract element vertices.
void insertEdge(int elmt, int side, int tagId)
Definition: InputSem.cpp:662
virtual NEKMESHUTILS_EXPORT void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
std::pair< ModuleType, std::string > ModuleKey
virtual NEKMESHUTILS_EXPORT void ProcessComposites()
Generate composites.
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215