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