Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 
39 
40 #include "OutputGmsh.h"
41 #include "../InputModules/InputGmsh.h"
42 
43 using namespace std;
44 using namespace Nektar::NekMeshUtils;
45 
46 namespace Nektar
47 {
48 namespace Utilities
49 {
50 
51 ModuleKey OutputGmsh::className =
53  OutputGmsh::create,
54  "Writes Gmsh msh file.");
55 
56 OutputGmsh::OutputGmsh(MeshSharedPtr m) : OutputModule(m)
57 {
58  map<unsigned int, ElmtConfig> igelmap = InputGmsh::GenElmMap();
60 
61  // Populate #InputGmsh::elmMap and use this to construct an
62  // inverse mapping from %ElmtConfig to Gmsh ID.
63  for (it = igelmap.begin(); it != igelmap.end(); ++it)
64  {
65  elmMap[it->second] = it->first;
66  }
67 
68  m_config["order"] = ConfigOption(false, "-1", "Enforce a polynomial order");
69 }
70 
72 {
73 }
74 
75 /**
76  * @brief Process a mesh to output to Gmsh MSH format.
77  *
78  * Gmsh output is fairly straightforward. The file first contains a
79  * list of nodes, followed by a list of elements. Since
80  * Mesh::vertexSet only contains vertices of the linear elements, we
81  * first loop over the elements so that any high-order vertices can be
82  * enumerated and then added to the node list. We then print out the
83  * list of nodes and finally print the element list.
84  */
86 {
87  if (m_mesh->m_verbose)
88  {
89  cout << "OutputGmsh: Writing file..." << endl;
90  }
91 
92  boost::unordered_map<int, vector<int> > orderingMap;
93  boost::unordered_map<int, vector<int> >::iterator oIt;
94 
95  // Open the file stream.
96  OpenStream();
97 
98  // Write MSH header
99  m_mshFile << "$MeshFormat" << endl
100  << "2.2 0 8" << endl
101  << "$EndMeshFormat" << endl;
102 
103  int id = m_mesh->m_vertexSet.size();
104  vector<ElementSharedPtr> toComplete;
105 
106  int order = m_config["order"].as<int>();
107 
108  if (order != -1)
109  {
110  if (m_mesh->m_verbose)
111  {
112  cout << "Making mesh of order " << order << endl;
113  }
114  }
115  else
116  {
117  // Do first pass over elements of expansion dimension to determine
118  // which elements need completion.
119  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
120  {
121  ElementSharedPtr e = m_mesh->m_element[m_mesh->m_expDim][i];
122  if (e->GetMaxOrder() > order)
123  {
124  order = e->GetMaxOrder();
125  }
126  }
127  }
128 
129  // Convert this mesh into a high-order mesh of uniform order.
130  if (m_mesh->m_verbose)
131  {
132  cout << "Mesh order of " << order << " detected" << endl;
133  }
134 
135  m_mesh->MakeOrder(order, LibUtilities::ePolyEvenlySpaced);
136 
137  // Add edge- and face-interior nodes to vertex set.
138  EdgeSet::iterator eIt;
139  FaceSet::iterator fIt;
140 
141  for (eIt = m_mesh->m_edgeSet.begin(); eIt != m_mesh->m_edgeSet.end(); ++eIt)
142  {
143  m_mesh->m_vertexSet.insert((*eIt)->m_edgeNodes.begin(),
144  (*eIt)->m_edgeNodes.end());
145  }
146 
147  for (fIt = m_mesh->m_faceSet.begin(); fIt != m_mesh->m_faceSet.end(); ++fIt)
148  {
149  m_mesh->m_vertexSet.insert((*fIt)->m_faceNodes.begin(),
150  (*fIt)->m_faceNodes.end());
151  }
152 
153  // Do second pass over elements for volume nodes.
154  for (int d = 1; d <= 3; ++d)
155  {
156  for (int i = 0; i < m_mesh->m_element[d].size(); ++i)
157  {
158  ElementSharedPtr e = m_mesh->m_element[d][i];
159  vector<NodeSharedPtr> volList = e->GetVolumeNodes();
160  m_mesh->m_vertexSet.insert(volList.begin(), volList.end());
161  }
162  }
163 
164  // Create ordered set of nodes - not required but looks nicer.
166  std::set<NodeSharedPtr> tmp(m_mesh->m_vertexSet.begin(),
167  m_mesh->m_vertexSet.end());
168 
169  // Write out nodes section.
170  m_mshFile << "$Nodes" << endl << m_mesh->m_vertexSet.size() << endl;
171 
172  for (it = tmp.begin(); it != tmp.end(); ++it)
173  {
174  m_mshFile << (*it)->m_id+1 << " " << scientific << setprecision(10)
175  << (*it)->m_x << " " << (*it)->m_y << " " << (*it)->m_z
176  << endl;
177  }
178 
179  m_mshFile << "$EndNodes" << endl;
180 
181  // Write elements section. All other sections are not currently
182  // supported (physical names etc).
183  m_mshFile << "$Elements" << endl;
184  m_mshFile << m_mesh->GetNumEntities() << endl;
185 
186  id = 1;
187 
188  for (int d = 1; d <= 3; ++d)
189  {
190  for (int i = 0; i < m_mesh->m_element[d].size(); ++i, ++id)
191  {
192  ElementSharedPtr e = m_mesh->m_element[d][i];
193 
194  // First output element ID and type.
195  int elmtType = elmMap[e->GetConf()];
196  m_mshFile << id << " " << elmtType << " ";
197 
198  // Write out number of element tags and then the tags
199  // themselves.
200  vector<int> tags = e->GetTagList();
201 
202  if (tags.size() == 1)
203  {
204  tags.push_back(tags[0]);
205  tags.push_back(0);
206  }
207 
208  m_mshFile << tags.size() << " ";
209 
210  for (int j = 0; j < tags.size(); ++j)
211  {
212  m_mshFile << tags[j] << " ";
213  }
214 
215  // Finally write out node list. First write vertices, then
216  // internal edge nodes, then face nodes.
217  vector<NodeSharedPtr> nodeList = e->GetVertexList();
218  vector<EdgeSharedPtr> edgeList = e->GetEdgeList();
219  vector<FaceSharedPtr> faceList = e->GetFaceList();
220  vector<NodeSharedPtr> volList = e->GetVolumeNodes();
221 
222  tags.clear();
223 
224  for (int j = 0; j < nodeList.size(); ++j)
225  {
226  tags.push_back(nodeList[j]->m_id);
227  }
228 
229  // Process edge-interior points
230  for (int j = 0; j < edgeList.size(); ++j)
231  {
232  nodeList = edgeList[j]->m_edgeNodes;
233 
234  if (e->GetEdgeOrient(j, edgeList[j]) == StdRegions::eForwards)
235  {
236  for (int k = 0; k < nodeList.size(); ++k)
237  {
238  tags.push_back(nodeList[k]->m_id);
239  }
240  }
241  else
242  {
243  for (int k = nodeList.size() - 1; k >= 0; --k)
244  {
245  tags.push_back(nodeList[k]->m_id);
246  }
247  }
248  }
249 
250  // Process face-interior points
251  for (int j = 0; j < faceList.size(); ++j)
252  {
253  nodeList = faceList[j]->m_faceNodes;
254  int nFaceVerts = faceList[j]->m_vertexList.size();
255  vector<int> faceIds(nFaceVerts), volFaceIds(nFaceVerts);
256 
257  for (int k = 0; k < nFaceVerts; ++k)
258  {
259  faceIds [k] = faceList[j]->m_vertexList[k]->m_id;
260  volFaceIds[k] =
261  e->GetVertexList()[e->GetFaceVertex(j, k)]->m_id;
262  }
263 
264  if (nFaceVerts == 3)
265  {
266  HOTriangle<NodeSharedPtr> hoTri(faceIds, nodeList);
267  hoTri.Align(volFaceIds);
268  for (int k = 0; k < hoTri.surfVerts.size(); ++k)
269  {
270  tags.push_back(hoTri.surfVerts[k]->m_id);
271  }
272  }
273  else
274  {
275  HOQuadrilateral<NodeSharedPtr> hoQuad(faceIds, nodeList);
276  hoQuad.Align(volFaceIds);
277 
278  for (int k = 0; k < hoQuad.surfVerts.size(); ++k)
279  {
280  tags.push_back(hoQuad.surfVerts[k]->m_id);
281  }
282  }
283  }
284 
285  // Process volume nodes
286  for (int j = 0; j < volList.size(); ++j)
287  {
288  tags.push_back(volList[j]->m_id);
289  }
290 
291  // Construct inverse of input reordering. First try to find it in
292  // our cache.
293  oIt = orderingMap.find(elmtType);
294 
295  // If it's not created, then create it.
296  if (oIt == orderingMap.end())
297  {
298  vector<int> reordering = InputGmsh::CreateReordering(elmtType);
299  vector<int> inv(tags.size());
300 
301  ASSERTL1(tags.size() == reordering.size(),
302  "Reordering map size not equal to element tags 1.");
303 
304  for (int j = 0; j < tags.size(); ++j)
305  {
306  inv[reordering[j]] = j;
307  }
308 
309  oIt = orderingMap.insert(make_pair(elmtType, inv)).first;
310  }
311 
312  ASSERTL1(tags.size() == oIt->second.size(),
313  "Reordering map size not equal to element tags 2.");
314 
315  // Finally write element nodes.
316  for (int j = 0; j < tags.size(); ++j)
317  {
318  m_mshFile << tags[oIt->second[j]] + 1 << " ";
319  }
320 
321  m_mshFile << endl;
322  }
323  }
324  m_mshFile << "$EndElements" << endl;
325 }
326 }
327 }
io::filtering_ostream m_mshFile
Output stream.
void Align(std::vector< int > vertId)
Align this surface to a given vertex ID.
Definition: HOAlignment.h:135
boost::unordered_map< NekMeshUtils::ElmtConfig, unsigned int, ElmtConfigHash > elmMap
Definition: OutputGmsh.h:86
NEKMESHUTILS_EXPORT void OpenStream()
Open a file for output.
std::vector< T > surfVerts
The quadrilateral surface vertices – templated so that this can either be nodes or IDs...
Definition: HOAlignment.h:231
Abstract base class for output modules.
STL namespace.
std::vector< T > surfVerts
The triangle surface vertices – templated so that this can either be nodes or IDs.
Definition: HOAlignment.h:65
pair< ModuleType, string > ModuleKey
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:65
static std::vector< int > CreateReordering(unsigned int InputGmshEntity)
Create a reordering map for a given element.
Definition: InputGmsh.cpp:1011
Represents a command-line configuration option.
std::map< std::string, ConfigOption > m_config
List of configuration values.
void Align(std::vector< int > vertId)
Align this surface to a given vertex ID.
Definition: HOAlignment.h:277
A lightweight struct for dealing with high-order triangle alignment.
Definition: HOAlignment.h:50
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
A lightweight struct for dealing with high-order quadrilateral alignment.
Definition: HOAlignment.h:215
virtual void Process()
Write mesh to output file.
Definition: OutputGmsh.cpp:85
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
std::pair< ModuleType, std::string > ModuleKey
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
static std::map< unsigned int, NekMeshUtils::ElmtConfig > GenElmMap()
Definition: InputGmsh.cpp:1685
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215