36 #include <boost/iostreams/filter/gzip.hpp>
41 namespace io = boost::iostreams;
45 namespace NekMeshUtils
56 Loki::SingleThreaded> Type;
57 return Type::Instance();
83 string fname =
m_config[
"infile"].as<
string>();
85 if (fname.size() > 3 && fname.substr(fname.size() - 3, 3) ==
".gz")
99 cerr <<
"Error opening file: " << fname << endl;
109 string fname =
m_config[
"outfile"].as<
string>();
111 if (fname.size() > 3 && fname.substr(fname.size() - 3, 3) ==
".gz")
125 cerr <<
"Error opening file: " << fname << endl;
140 vector<ElementSharedPtr> &elmt =
m_mesh->m_element[
m_mesh->m_expDim];
142 m_mesh->m_vertexSet.clear();
144 for (
int i = 0, vid = 0; i < elmt.size(); ++i)
146 for (
int j = 0; j < elmt[i]->GetVertexCount(); ++j)
148 pair<NodeSet::iterator,bool> testIns =
149 m_mesh->m_vertexSet.insert(elmt[i]->GetVertex(j));
153 (*(testIns.first))->m_id = vid++;
157 elmt[i]->SetVertex(j,*testIns.first);
180 if (
m_mesh->m_expDim < 2)
return;
184 vector<ElementSharedPtr> &elmt =
m_mesh->m_element[
m_mesh->m_expDim];
186 m_mesh->m_edgeSet.clear();
191 for (
int i = 0, eid = 0; i < elmt.size(); ++i)
193 for (
int j = 0; j < elmt[i]->GetEdgeCount(); ++j)
195 pair<EdgeSet::iterator,bool> testIns;
197 testIns =
m_mesh->m_edgeSet.insert(ed);
203 ed2->m_elLink.push_back(
204 pair<ElementSharedPtr,int>(elmt[i],j));
209 elmt[i]->SetEdge(j, e2);
210 if (e2->m_edgeNodes.size() == 0 &&
211 ed->m_edgeNodes.size() > 0)
213 e2->m_curveType = ed->m_curveType;
214 e2->m_edgeNodes = ed->m_edgeNodes;
217 if (e2->m_n1->m_id != ed->m_n1->m_id)
219 reverse(e2->m_edgeNodes.begin(),
220 e2->m_edgeNodes.end());
224 e2->m_parentCAD = ed->m_parentCAD;
227 e2->m_elLink.push_back(
228 pair<ElementSharedPtr,int>(elmt[i],j));
235 for (
int i = 0; i <
m_mesh->m_element[1].size(); ++i)
240 vector<NodeSharedPtr> edgeNodes;
242 new Edge(v0, v1, edgeNodes, elmt->GetConf().m_edgeCurveType));
245 if (it ==
m_mesh->m_edgeSet.end())
247 cerr <<
"Cannot find corresponding element edge for "
248 <<
"1D element " << i << endl;
251 elmt->SetEdgeLink(*it);
254 pair<ElementSharedPtr, int> eMap = (*it)->m_elLink.at(0);
255 eMap.first->SetBoundaryLink(eMap.second, i);
258 elmt->SetVertex(0, (*it)->m_n1,
false);
259 elmt->SetVertex(1, (*it)->m_n2,
false);
262 if ((*it)->m_edgeNodes.size() > 0)
265 if (edge->GetVertex(0) == (*it)->m_n1)
267 edge->SetVolumeNodes((*it)->m_edgeNodes);
269 elmt->SetCurveType((*it)->m_curveType);
292 if (
m_mesh->m_expDim < 3)
return;
296 vector<ElementSharedPtr> &elmt =
m_mesh->m_element[
m_mesh->m_expDim];
298 m_mesh->m_faceSet.clear();
301 for (
int i = 0, fid = 0; i < elmt.size(); ++i)
303 for (
int j = 0; j < elmt[i]->GetFaceCount(); ++j)
305 pair<FaceSet::iterator,bool> testIns;
306 testIns =
m_mesh->m_faceSet.insert(elmt[i]->GetFace(j));
310 (*(testIns.first))->m_id = fid++;
311 (*(testIns.first))->m_elLink.push_back(
312 pair<ElementSharedPtr,int>(elmt[i],j));
316 elmt[i]->SetFace(j,*testIns.first);
318 (*(testIns.first))->m_elLink.push_back(
319 pair<ElementSharedPtr,int>(elmt[i],j));
326 for (
int i = 0; i <
m_mesh->m_element[2].size(); ++i)
329 vector<NodeSharedPtr> vertices = elmt->GetVertexList();
330 vector<NodeSharedPtr> faceNodes;
331 vector<EdgeSharedPtr> edgeList = elmt->GetEdgeList();
333 new Face(vertices, faceNodes, edgeList,
334 elmt->GetConf().m_faceCurveType));
337 if (it ==
m_mesh->m_faceSet.end())
339 cout <<
"Cannot find corresponding element face for 2D "
340 <<
"element " << i << endl;
344 elmt->SetFaceLink(*it);
347 for (
int j = 0; j < elmt->GetVertexCount(); ++j)
349 elmt->SetVertex(j, (*it)->m_vertexList[j],
false);
350 elmt->SetEdge(j, (*it)->m_edgeList[j],
false);
353 EdgeSet tmp(edgeList.begin(),edgeList.end());
355 for (
int j = 0; j < elmt->GetEdgeCount(); ++j)
361 e->m_parentCAD = (*f)->m_parentCAD;
366 pair<ElementSharedPtr, int> eMap = (*it)->m_elLink.at(0);
367 eMap.first->SetBoundaryLink(eMap.second, i);
370 if ((*it)->m_faceNodes.size() > 0)
372 elmt->SetVolumeNodes((*it)->m_faceNodes);
373 elmt->SetCurveType((*it)->m_curveType);
389 for (
int i = 0; i <
m_mesh->m_element[
m_mesh->m_expDim].size(); ++i)
406 m_mesh->m_composite.clear();
411 for (
int d = 0; d <=
m_mesh->m_expDim; ++d)
413 vector<ElementSharedPtr> &elmt =
m_mesh->m_element[d];
415 for (
int i = 0; i < elmt.size(); ++i)
418 unsigned int tagid = elmt[i]->GetTagList()[0];
420 it =
m_mesh->m_composite.find(tagid);
422 if (it ==
m_mesh->m_composite.end())
426 pair<CompositeMap::iterator, bool> testIns;
428 tmp->m_tag = elmt[i]->GetTag();
429 if(
m_mesh->m_faceLabels.count(tmp->m_id) != 0)
431 tmp->m_label =
m_mesh->m_faceLabels[tmp->m_id];
434 testIns =
m_mesh->m_composite.insert(
435 pair<unsigned int, CompositeSharedPtr>(tagid,tmp));
439 if (elmt[i]->GetTag() != it->second->m_tag)
441 cout <<
"Different types of elements in same composite!" << endl;
442 cout <<
" -> Composite uses " << it->second->m_tag << endl;
443 cout <<
" -> Element uses " << elmt[i]->GetTag() << endl;
444 cout <<
"Have you specified physical volumes and surfaces?" << endl;
446 it->second->m_items.push_back(elmt[i]);
458 for(eit =
m_mesh->m_edgeSet.begin(); eit !=
m_mesh->m_edgeSet.end(); eit++)
460 (*eit)->m_elLink.clear();
465 for(fit =
m_mesh->m_faceSet.begin(); fit !=
m_mesh->m_faceSet.end(); fit++)
467 (*fit)->m_elLink.clear();
508 set<int> prismsDone, tetsDone;
512 for (i = 0; i <
m_mesh->m_element[3].size(); ++i)
518 prismsDone.insert(i);
528 for (it =
m_mesh->m_vertexSet.begin(); it !=
m_mesh->m_vertexSet.end(); ++it)
535 int prismTris[2][3] = {{0,1,4}, {3,2,5}};
538 bool warnCurvature =
false;
542 boost::unordered_set<int> facesDone;
546 while (prismsDone.size() > 0)
548 vector<ElementSharedPtr> line;
552 PrismLines(*prismsDone.begin(), perFaces, prismsDone, line);
556 for (i = 0; i < line.size(); ++i)
559 vector<int> tags = line[i]->GetTagList();
560 vector<NodeSharedPtr> nodes = line[i]->GetVertexList();
565 line[i]->GetFace(1), line[i]->GetFace(3)
568 fIt[0] = facesDone.find(f[0]->m_id);
569 fIt[1] = facesDone.find(f[1]->m_id);
573 for (j = 0; j < 2; ++j)
575 pIt = perFaces.find(f[j]->m_id);
577 if (pIt == perFaces.end())
582 fIt2 = facesDone.find(pIt->second.first->m_id);
584 if (fIt[j] == facesDone.end() &&
585 fIt2 != facesDone.end())
591 if (fIt[0] != facesDone.end() &&
592 fIt[1] != facesDone.end())
596 ASSERTL0(
false,
"Renumbering error!");
598 else if (fIt[0] == facesDone.end() &&
599 fIt[1] == facesDone.end())
602 for (j = 0; j < 2; ++j)
604 for (k = 0; k < 3; ++k)
614 facesDone.insert(f[0]->m_id);
615 facesDone.insert(f[1]->m_id);
621 int t = fIt[0] == facesDone.end() ? 0 : 1;
623 ASSERTL1(fIt[o] != facesDone.end(),
"Renumbering error");
629 nodes[prismTris[o][0]]->m_id,
630 nodes[prismTris[o][1]]->m_id,
631 nodes[prismTris[o][2]]->m_id
633 int tmp2[3] = {0,1,2};
635 if (tmp1[0] > tmp1[1])
637 swap(tmp1[0], tmp1[1]);
638 swap(tmp2[0], tmp2[1]);
641 if (tmp1[1] > tmp1[2])
643 swap(tmp1[1], tmp1[2]);
644 swap(tmp2[1], tmp2[2]);
647 if (tmp1[0] > tmp1[2])
649 swap(tmp1[0], tmp1[2]);
650 swap(tmp2[0], tmp2[2]);
654 for (j = 0; j < 3; ++j)
663 facesDone.insert(f[t]->m_id);
666 for (j = 0; j < 6; ++j)
668 ASSERTL1(nodes[j]->m_id != -1,
"Renumbering error");
678 for (j = 0; j < 9; ++j)
681 for (k = 0; k < 9; ++k)
684 if (e1->m_n1 == e2->m_n1 && e1->m_n2 == e2->m_n2)
686 e2->m_edgeNodes = e1->m_edgeNodes;
688 else if (e1->m_n1 == e2->m_n1 && e1->m_n2 == e2->m_n2)
690 e2->m_edgeNodes = e1->m_edgeNodes;
691 std::reverse(e2->m_edgeNodes.begin(),
692 e2->m_edgeNodes.end());
700 for (j = 0; j < 5; ++j)
702 if (line[i]->GetFace(j)->m_faceNodes.size() > 0)
704 warnCurvature =
true;
711 m_mesh->m_element[3][line[i]->GetId()] = el;
717 cerr <<
"[ReorderPrisms] WARNING: Face curvature detected in "
718 <<
"some prisms; this will be ignored in further module "
724 for (pIt = perFaces.begin(); pIt != perFaces.end(); ++pIt)
728 vector<int> perVerts = pIt->second.second;
729 int nVerts = perVerts.size();
732 for (j = 0; j < nVerts; ++j)
737 if (n1->m_id == -1 && n2->m_id == -1)
742 else if (n1->m_id != -1 && n2->m_id != -1)
748 ASSERTL0(
false,
"Periodic face renumbering error");
755 for (it2 = tetsDone.begin(); it2 != tetsDone.end(); ++it2)
758 vector<NodeSharedPtr> nodes = el->GetVertexList();
759 vector<int> tags = el->GetTagList();
761 for (i = 0; i < 4; ++i)
763 if (nodes[i]->m_id == -1)
765 nodes[i]->m_id = nodeId++;
776 for (it =
m_mesh->m_vertexSet.begin(); it !=
m_mesh->m_vertexSet.end(); ++it)
778 if ((*it)->m_id == -1)
780 (*it)->m_id = nodeId++;
791 set<int> &prismsDone,
792 vector<ElementSharedPtr> &line)
798 if (it == prismsDone.end())
804 prismsDone.erase(it);
805 line.push_back(
m_mesh->m_element[3][prism]);
808 for (i = 1; i <= 3; i += 2)
814 it2 = perFaces.find(f->m_id);
816 if (it2 != perFaces.end())
818 int id2 = it2->second.first->m_id;
819 nextId = it2->second.first->m_elLink[0].first->GetId();
822 PrismLines(nextId, perFaces, prismsDone, line);
826 if (f->m_elLink.size() == 1)
831 nextId = f->m_elLink[0].first->GetId();
832 if (nextId ==
m_mesh->m_element[3][prism]->GetId())
834 nextId = f->m_elLink[1].first->GetId();
837 PrismLines(nextId, perFaces, prismsDone, line);
849 cerr <<
"WARNING: Unrecognised config option " << key
850 <<
", proceeding anyway." << endl;
853 it->second.beenSet =
true;
855 if (it->second.isBool)
857 it->second.value =
"1";
861 it->second.value = val;
874 cerr <<
"No configuration options for this module." << endl;
880 cerr << setw(10) << it->first <<
": " << it->second.desc
895 if (!it->second.beenSet)
897 it->second.value = it->second.defValue;
909 cerr <<
"Expansion dimension is " <<
m_mesh->m_expDim << endl;
910 cerr <<
"Space dimension is " <<
m_mesh->m_spaceDim << endl;
911 cerr <<
"Read " <<
m_mesh->m_node.size() <<
" nodes" << endl;
912 cerr <<
"Read " <<
m_mesh->GetNumElements() <<
" "
913 <<
m_mesh->m_expDim <<
"-D elements" << endl;
914 cerr <<
"Read " <<
m_mesh->GetNumBndryElements()
915 <<
" boundary elements" << endl;
NEKMESHUTILS_EXPORT void SetDefaults()
Sets default configuration options for those which have not been set.
#define ASSERTL0(condition, msg)
Basic information about an element.
io::filtering_ostream m_mshFile
Output 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.
Represents an edge which joins two points.
NEKMESHUTILS_EXPORT void OpenStream()
Open a file for output.
Represents a face comprised of three or more edges.
std::map< int, std::pair< FaceSharedPtr, std::vector< int > > > PerMap
LibUtilities::NekFactory< ModuleKey, Module, FieldSharedPtr > ModuleFactory
NEKMESHUTILS_EXPORT void RegisterConfig(std::string key, std::string value)
Register a configuration option with a module.
MeshSharedPtr m_mesh
Mesh object.
NEKMESHUTILS_EXPORT void PrintConfig()
Print out all configuration options for a module.
ElementFactory & GetElementFactory()
std::ofstream m_mshFileStream
Input stream.
ModuleFactory & GetModuleFactory()
NEKMESHUTILS_EXPORT void ReorderPrisms(PerMap &perFaces)
Reorder node IDs so that prisms and tetrahedra are aligned correctly.
const char *const ModuleTypeMap[]
virtual NEKMESHUTILS_EXPORT void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
virtual NEKMESHUTILS_EXPORT void ProcessElements()
Generate element IDs.
boost::shared_ptr< Composite > CompositeSharedPtr
Shared pointer to a composite.
Represents a command-line configuration option.
boost::shared_ptr< Node > NodeSharedPtr
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< GeometryVector > Composite
virtual NEKMESHUTILS_EXPORT void ClearElementLinks()
NEKMESHUTILS_EXPORT OutputModule(MeshSharedPtr p_m)
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
NEKMESHUTILS_EXPORT void PrismLines(int prism, PerMap &perFaces, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)
std::ostream & operator<<(std::ostream &os, const ModuleKey &rhs)
virtual NEKMESHUTILS_EXPORT void ProcessVertices()
Extract element vertices.
virtual NEKMESHUTILS_EXPORT void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
boost::shared_ptr< Element > ElementSharedPtr
boost::unordered_set< EdgeSharedPtr, EdgeHash > EdgeSet
boost::shared_ptr< Face > FaceSharedPtr
std::pair< ModuleType, std::string > ModuleKey
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
virtual NEKMESHUTILS_EXPORT void ProcessComposites()
Generate composites.
Provides a generic Factory class.