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