Nektar++
OutputGmsh.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: OutputGmsh.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: Gmsh file format output.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include "../MeshElements.h"
37 #include "OutputGmsh.h"
38 #include "../InputModules/InputGmsh.h"
39 
40 #include <map>
41 #include <vector>
42 #include <string>
43 using namespace std;
44 
45 namespace Nektar
46 {
47  namespace Utilities
48  {
49  ModuleKey OutputGmsh::className =
51  ModuleKey(eOutputModule, "msh"), OutputGmsh::create,
52  "Writes Gmsh msh file.");
53 
54  OutputGmsh::OutputGmsh(MeshSharedPtr m) : OutputModule(m)
55  {
57 
58  // Populate #InputGmsh::elmMap and use this to construct an
59  // inverse mapping from %ElmtConfig to Gmsh ID.
60  for (it = InputGmsh::elmMap.begin(); it != InputGmsh::elmMap.end(); ++it)
61  {
62  elmMap[it->second] = it->first;
63  }
64  }
65 
67  {
68 
69  }
70 
71  /**
72  * @brief Process a mesh to output to Gmsh MSH format.
73  *
74  * Gmsh output is fairly straightforward. The file first contains a
75  * list of nodes, followed by a list of elements. Since
76  * Mesh::vertexSet only contains vertices of the linear elements, we
77  * first loop over the elements so that any high-order vertices can be
78  * enumerated and then added to the node list. We then print out the
79  * list of nodes and finally print the element list.
80  */
82  {
83  if (m_mesh->m_verbose)
84  {
85  cout << "OutputGmsh: Writing file..." << endl;
86  }
87 
88  // Open the file stream.
89  OpenStream();
90 
91  // Write MSH header
92  m_mshFile << "$MeshFormat" << endl
93  << "2.2 0 8" << endl
94  << "$EndMeshFormat" << endl;
95 
96  int id = m_mesh->m_vertexSet.size();
97  vector<ElementSharedPtr> toComplete;
98 
99  int maxOrder = -1;
100 
101  // Do first pass over elements of expansion dimension to determine
102  // which elements need completion.
103  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
104  {
105  ElementSharedPtr e = m_mesh->m_element[m_mesh->m_expDim][i];
106  if (e->GetMaxOrder() > maxOrder)
107  {
108  maxOrder = e->GetMaxOrder();
109  }
110  }
111 
112  //maxOrder = 2;
113  for (int d = 1; d <= 3; ++d)
114  {
115  for (int i = 0; i < m_mesh->m_element[d].size(); ++i)
116  {
117  ElementSharedPtr e = m_mesh->m_element[d][i];
118  if ((e->GetConf().m_order <= 1 && maxOrder > 1) ||
119  (e->GetConf().m_order == maxOrder && e->GetConf().m_faceNodes == false))
120  {
121  toComplete.push_back(e);
122  }
123  // Generate geometry information for this element. This will
124  // be stored locally inside each element.
126  m_mesh->m_element[d][i]->GetGeom(m_mesh->m_spaceDim);
127  }
128  }
129 
130  // Complete these elements.
131  for (int i = 0; i < toComplete.size(); ++i)
132  {
133  toComplete[i]->Complete(maxOrder);
134  }
135 
136  // Do second pass over elements to enumerate high-order vertices.
137  for (int d = 1; d <= 3; ++d)
138  {
139  for (int i = 0; i < m_mesh->m_element[d].size(); ++i)
140  {
141  // Keep track of faces and edges to ensure that high-order
142  // nodes are only added once on common faces/edges.
143  boost::unordered_set<int> edgesDone;
144  boost::unordered_set<int> facesDone;
145  ElementSharedPtr e = m_mesh->m_element[d][i];
146 
147  if (e->GetConf().m_order > 1)
148  {
149  vector<NodeSharedPtr> tmp;
150  vector<EdgeSharedPtr> edgeList = e->GetEdgeList();
151  vector<FaceSharedPtr> faceList = e->GetFaceList();
152  vector<NodeSharedPtr> volList = e->GetVolumeNodes();
153 
154  for (int j = 0; j < edgeList.size(); ++j)
155  {
157  edgesDone.find(edgeList[j]->m_id);
158  if (it == edgesDone.end() || d != 3)
159  {
160  tmp.insert(tmp.end(),
161  edgeList[j]->m_edgeNodes.begin(),
162  edgeList[j]->m_edgeNodes.end());
163  edgesDone.insert(edgeList[j]->m_id);
164  }
165  }
166 
167  for (int j = 0; j < faceList.size(); ++j)
168  {
170  facesDone.find(faceList[j]->m_id);
171  if (it == facesDone.end() || d != 3)
172  {
173  tmp.insert(tmp.end(),
174  faceList[j]->m_faceNodes.begin(),
175  faceList[j]->m_faceNodes.end());
176  facesDone.insert(faceList[j]->m_id);
177  }
178  }
179 
180  tmp.insert(tmp.end(), volList.begin(), volList.end());
181 
182  // Even though faces/edges are at this point unique
183  // across the mesh, still need to test inserts since
184  // high-order nodes may already have been inserted into
185  // the list from an adjoining element or a boundary
186  // element.
187  for (int j = 0; j < tmp.size(); ++j)
188  {
189  pair<NodeSet::iterator, bool> testIns =
190  m_mesh->m_vertexSet.insert(tmp[j]);
191 
192  if (testIns.second)
193  {
194  (*(testIns.first))->m_id = id++;
195  }
196  else
197  {
198  tmp[j]->m_id = (*(testIns.first))->m_id;
199  }
200  }
201  }
202  }
203  }
204 
205  // Create ordered set of nodes - not required but looks nicer.
207  std::set<NodeSharedPtr> tmp(m_mesh->m_vertexSet.begin(),
208  m_mesh->m_vertexSet.end());
209 
210  // Write out nodes section.
211  m_mshFile << "$Nodes" << endl
212  << m_mesh->m_vertexSet.size() << endl;
213 
214  for (it = tmp.begin(); it != tmp.end(); ++it)
215  {
216  m_mshFile << (*it)->m_id << " " << scientific << setprecision(10)
217  << (*it)->m_x << " "
218  << (*it)->m_y << " " << (*it)->m_z
219  << endl;
220  }
221 
222  m_mshFile << "$EndNodes" << endl;
223 
224  // Write elements section. All other sections are not currently
225  // supported (physical names etc).
226  m_mshFile << "$Elements" << endl;
227  m_mshFile << m_mesh->GetNumEntities() << endl;
228 
229  id = 0;
230 
231  for (int d = 1; d <= 3; ++d)
232  {
233  for (int i = 0; i < m_mesh->m_element[d].size(); ++i, ++id)
234  {
235  ElementSharedPtr e = m_mesh->m_element[d][i];
236 
237  // First output element ID and type.
238  m_mshFile << id << " "
239  << elmMap[e->GetConf()] << " ";
240 
241  // Write out number of element tags and then the tags
242  // themselves.
243  vector<int> tags = e->GetTagList();
244 
245  if (tags.size() == 1)
246  {
247  tags.push_back(tags[0]);
248  tags.push_back(0);
249  }
250 
251  m_mshFile << tags.size() << " ";
252 
253  for (int j = 0; j < tags.size(); ++j)
254  {
255  m_mshFile << tags[j] << " ";
256  }
257 
258  // Finally write out node list. First write vertices, then
259  // internal edge nodes, then face nodes.
260  vector<NodeSharedPtr> nodeList = e->GetVertexList ();
261  vector<EdgeSharedPtr> edgeList = e->GetEdgeList ();
262  vector<FaceSharedPtr> faceList = e->GetFaceList ();
263  vector<NodeSharedPtr> volList = e->GetVolumeNodes();
264 
265  tags.clear();
266 
267  for (int j = 0; j < nodeList.size(); ++j)
268  {
269  tags.push_back(nodeList[j]->m_id);
270  }
271 
272  if (e->GetConf().m_order > 1)
273  {
274  for (int j = 0; j < edgeList.size(); ++j)
275  {
276  nodeList = edgeList[j]->m_edgeNodes;
277  for (int k = 0; k < nodeList.size(); ++k)
278  {
279  tags.push_back(nodeList[k]->m_id);
280  //cout << "EDGENODE" << endl;
281  }
282  }
283 
284  for (int j = 0; j < faceList.size(); ++j)
285  {
286  nodeList = faceList[j]->m_faceNodes;
287  for (int k = 0; k < nodeList.size(); ++k)
288  {
289  //cout << "FACENODE" << endl;
290  tags.push_back(nodeList[k]->m_id);
291  }
292  }
293 
294  for (int j = 0; j < volList.size(); ++j)
295  {
296  //cout << "VOLNODE" << endl;
297  tags.push_back(volList[j]->m_id);
298  }
299  }
300 
301  // Re-order tetrahedral vertices.
302  if (e->GetConf().m_e == LibUtilities::eTetrahedron)
303  {
304  int order = e->GetConf().m_order;
305  if (order > 4)
306  {
307  cerr << "Temporary error: Gmsh tets only supported "
308  << "up to 4th order - will fix soon!" << endl;
309  abort();
310  }
311  int pos = 4;
312  // Swap edge 1->3 nodes with edge 2->3 nodes.
313  pos = 4 + 4*(order-1);
314  for (int j = 0; j < order-1; ++j)
315  {
316  swap(tags[j+pos], tags[j+pos+order-1]);
317  }
318  // Reverse ordering of other vertical edge-interior
319  // nodes.
320  reverse(tags.begin()+4+3*(order-1), tags.begin()+4+4*(order-1));
321  reverse(tags.begin()+4+4*(order-1), tags.begin()+4+5*(order-1));
322  reverse(tags.begin()+4+5*(order-1), tags.begin()+4+6*(order-1));
323 
324  // Swap face 2 nodes with face 3.
325  pos = 4 + 6*(order-1) + 2*(order-2)*(order-1)/2;
326  for (int j = 0; j < (order-2)*(order-1)/2; ++j)
327  {
328  swap(tags[j+pos], tags[j+pos+(order-2)*(order-1)/2]);
329  }
330 
331  // Re-order face points. Gmsh ordering (node->face) is:
332  //
333  // Face 0: 0->2->1
334  // Face 1: 0->1->3
335  // Face 2: 0->3->2
336  // Face 3: 3->1->2
337  //
338  // Therefore need to reorder nodes for faces 0, 2 and
339  // 3 to match nodal ordering.
340 
341  // Re-order face 0: transpose
342  vector<int> tmp((order-2)*(order-1)/2);
343  int a = 0;
344  pos = 4 + 6*(order-1);
345  for (int j = 0; j < order-2; ++j)
346  {
347  for (int k = 0; k < order-2-j; ++k, ++a)
348  {
349  tmp[a] = tags[pos+j+k*(2*(order-2)+1-k)/2];
350  }
351  }
352  for (int j = 0; j < (order-1)*(order-2)/2; ++j)
353  {
354  tags[pos+j] = tmp[j];
355  }
356 
357  // Re-order face 2: transpose
358  pos = 4 + 6*(order-1) + 2*(order-2)*(order-1)/2;
359  a = 0;
360  for (int j = 0; j < order-2; ++j)
361  {
362  for (int k = 0; k < order-2-j; ++k, ++a)
363  {
364  tmp[a] = tags[pos+j+k*(2*(order-2)+1-k)/2];
365  }
366  }
367  for (int j = 0; j < (order-1)*(order-2)/2; ++j)
368  {
369  tags[pos+j] = tmp[j];
370  }
371 
372  // Re-order face 3: reflect in y direction
373  pos = 4 + 6*(order-1)+3*(order-2)*(order-1)/2;
374  a = 0;
375  for (int j = 0; j < order-2; ++j)
376  {
377  for (int k = order-3-j; k >= 0; --k, ++a)
378  {
379  tmp[a] = tags[pos+j+k*(2*(order-2)+1-k)/2];
380  }
381  }
382 
383  for (int j = 0; j < (order-1)*(order-2)/2; ++j)
384  {
385  tags[pos+j] = tmp[j];
386  }
387  }
388  // Re-order prism vertices.
389  else if (e->GetConf().m_e == LibUtilities::ePrism)
390  {
391  int order = e->GetConf().m_order;
392  if (order > 2)
393  {
394  cerr << "Temporary error: Gmsh prisms only "
395  << "supported up to 2nd order!" << endl;
396  abort();
397  }
398 
399  // Swap nodes.
400  vector<int> temp (18);
401  temp[0] = tags[0];
402  temp[1] = tags[1];
403  temp[2] = tags[4];
404  temp[3] = tags[3];
405  temp[4] = tags[2];
406  temp[5] = tags[5];
407  temp[6] = tags[6];
408  temp[7] = tags[10];
409  temp[8] = tags[9];
410  temp[9] = tags[11];
411  temp[10] = tags[7];
412  temp[11] = tags[14];
413  temp[12] = tags[8];
414  temp[13] = tags[13];
415  temp[14] = tags[12];
416  temp[15] = tags[15];
417  temp[16] = tags[17];
418  temp[17] = tags[16];
419  for (int k = 0; k < 18; ++k)
420  {
421  tags[k] = temp[k];
422  }
423  }
424 
425  // Finally write element nodes.
426  for (int j = 0; j < tags.size(); ++j)
427  {
428  m_mshFile << tags[j] << " ";
429  }
430 
431  m_mshFile << endl;
432  }
433  }
434  m_mshFile << "$EndElements" << endl;
435  }
436  }
437 }
pair< ModuleType, string > ModuleKey
Abstract base class for output modules.
STL namespace.
MeshSharedPtr m_mesh
Mesh object.
boost::shared_ptr< Element > ElementSharedPtr
Shared pointer to an element.
Definition: MeshElements.h:63
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
virtual void Process()
Write mesh to output file.
Definition: OutputGmsh.cpp:81
std::ofstream m_mshFile
Output stream.
void OpenStream()
Open a file for output.
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
boost::unordered_map< ElmtConfig, unsigned int, ElmtConfigHash > elmMap
Definition: OutputGmsh.h:85
ModuleFactory & GetModuleFactory()
static std::map< unsigned int, ElmtConfig > elmMap
Definition: InputGmsh.h:66
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215