41 #include <boost/algorithm/string.hpp>
63 true,
"0",
"If defined then assume input rea is for scalar problem");
85 int nParam, nElements, nCurves;
86 int i, j, k, nodeCounter = 0;
90 map<LibUtilities::ShapeType, int> domainComposite;
91 map<LibUtilities::ShapeType, vector<vector<NodeSharedPtr> > > elNodes;
92 map<LibUtilities::ShapeType, vector<int> > elIds;
93 boost::unordered_map<int, int> elMap;
94 vector<LibUtilities::ShapeType> elmOrder;
96 bool scalar =
m_config[
"scalar"].as<
bool>();
112 cout <<
"InputNek: Start reading file..." << endl;
118 for (i = 0; i < 4; ++i)
123 stringstream s(line);
126 for (i = 0; i < nParam; ++i)
139 for (i = 0; i < j; ++i)
149 for (i = 0; i < j; ++i)
157 bool foundMesh =
false;
161 if (line.find(
"MESH") != string::npos)
170 cerr <<
"Couldn't find MESH tag inside file." << endl;
178 s >> nElements >>
m_mesh->m_expDim;
182 m_mesh->m_fields.push_back(
"u");
186 m_mesh->m_fields.push_back(
"v");
187 if (
m_mesh->m_spaceDim > 2)
189 m_mesh->m_fields.push_back(
"w");
191 m_mesh->m_fields.push_back(
"p");
195 for (i = 0; i < nElements; ++i)
199 if (
m_mesh->m_expDim == 2)
201 if (line.find(
"Qua") != string::npos ||
202 line.find(
"qua") != string::npos)
214 if (line.find(
"Tet") != string::npos ||
215 line.find(
"tet") != string::npos)
219 else if (line.find(
"Hex") != string::npos ||
220 line.find(
"hex") != string::npos)
224 else if (line.find(
"Prism") != string::npos ||
225 line.find(
"prism") != string::npos)
229 else if (line.find(
"Pyr") != string::npos ||
230 line.find(
"pyr") != string::npos)
234 else if (line.find(
"Qua") != string::npos ||
235 line.find(
"qua") != string::npos)
249 for (j = 0; j <
m_mesh->m_expDim; ++j)
254 for (k = 0; k < nNodes; ++k)
261 for (j =
m_mesh->m_expDim; j < 3; ++j)
263 for (k = 0; k < nNodes; ++k)
272 vector<NodeSharedPtr> nodeList;
273 for (k = 0; k < nNodes; ++k)
276 nodeCounter++, vertex[0][k], vertex[1][k], vertex[2][k]));
277 nodeList.push_back(n);
280 elNodes[elType].push_back(nodeList);
281 elIds[elType].push_back(i);
287 for (i = 0; i < elmOrder.size(); ++i)
290 vector<vector<NodeSharedPtr> > &tmp = elNodes[elType];
292 for (j = 0; j < tmp.size(); ++j)
296 domainComposite.find(elType);
297 if (compIt == domainComposite.end())
299 tags.push_back(nComposite);
300 domainComposite[elType] = nComposite;
305 tags.push_back(compIt->second);
308 elMap[elIds[elType][j]] = reorderedId++;
310 vector<NodeSharedPtr> nodeList = tmp[j];
312 for (k = 0; k < nodeList.size(); ++k)
314 pair<NodeSet::iterator, bool> testIns =
315 m_mesh->m_vertexSet.insert(nodeList[k]);
319 nodeList[k] = *(testIns.first);
323 nodeList[k]->m_id = nodeCounter++;
330 elType, conf, nodeList, tags);
331 m_mesh->m_element[E->GetDim()].push_back(E);
337 if (line.find(
"CURVE") == string::npos)
339 cerr <<
"Cannot find curved side data." << endl;
353 for (i = 0; i < nCurves; ++i)
366 s >> word >> curveTag;
369 else if (word ==
"Recon")
375 s >> word >> curveTag;
380 cerr <<
"Unsupported curve type " << word << endl;
392 if (line.find(
"side") == string::npos)
394 cerr <<
"Unable to read number of curved sides" << endl;
400 map<string, pair<NekCurve, string> >
::iterator it;
409 for (i = 0; i < nCurvedSides; ++i)
414 s >> faceId >> elId >> word;
416 elId = elMap[elId - 1];
419 int origFaceId = faceId;
423 boost::shared_ptr<Prism>
p =
424 boost::dynamic_pointer_cast<
Prism>(el);
425 if (p->m_orientation == 1)
427 faceId = (faceId + 2) % 6;
429 else if (p->m_orientation == 2)
431 faceId = (faceId + 4) % 6;
436 boost::shared_ptr<Tetrahedron> t =
444 cerr <<
"Unrecognised curve tag " << word <<
" in curved lines"
449 if (it->second.first ==
eRecon)
452 vector<NodeSharedPtr> &tmp = el->GetFace(faceId)->m_vertexList;
453 vector<Node> n(tmp.size());
460 boost::dynamic_pointer_cast<
Prism>(el)->m_orientation;
467 for (j = 0; j < tmp.size(); ++j)
475 for (j = 0; j < tmp.size(); ++j)
483 for (j = 0; j < tmp.size(); ++j)
488 for (j = 0; j < tmp.size(); ++j)
490 int id = tmp[(j + offset) % tmp.size()]->m_id;
492 m_mesh->m_vertexNormals.find(
id);
494 if (vIt ==
m_mesh->m_vertexNormals.end())
496 m_mesh->m_vertexNormals[id] = n[j];
502 m_mesh->m_spherigonSurfs.insert(make_pair(elId, faceId));
504 else if (it->second.first ==
eFile)
507 static int tetFaceVerts[4][3] = {
508 {0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};
510 vector<int> vertId(3);
511 s >> vertId[0] >> vertId[1] >> vertId[2];
518 boost::shared_ptr<Tetrahedron> tet =
520 vector<int> tmpVertId = vertId;
522 for (j = 0; j < 3; ++j)
526 [tetFaceVerts[origFaceId][j]])
529 for (k = 0; k < 3; ++k)
531 int w = f->m_vertexList[k]->m_id;
534 vertId[k] = tmpVertId[j];
542 boost::shared_ptr<Prism> pr =
543 boost::static_pointer_cast<
Prism>(el);
544 if (pr->m_orientation == 1)
546 swap(vertId[2], vertId[1]);
547 swap(vertId[0], vertId[1]);
549 else if (pr->m_orientation == 2)
551 swap(vertId[0], vertId[1]);
552 swap(vertId[2], vertId[1]);
557 boost::shared_ptr<HOSurf>(
new HOSurf(vertId));
559 hoIt =
hoData[word].find(hs);
561 if (hoIt ==
hoData[word].end())
563 cerr <<
"Unable to find high-order surface data "
564 <<
"for element id " << elId + 1 << endl;
574 int Ntot = (*hoIt)->surfVerts.size();
575 int N = ((int)sqrt(8.0 * Ntot + 1.0) - 1) / 2;
580 vector<NodeSharedPtr> tmpVerts = (*hoIt)->surfVerts;
581 for (j = 0; j < tmpVerts.size(); ++j)
583 (*hoIt)->surfVerts[
hoMap[j]] = tmpVerts[j];
586 for (j = 0; j < tmpVerts.size(); ++j)
591 vector<int> faceVertIds(3);
592 faceVertIds[0] = f->m_vertexList[0]->m_id;
593 faceVertIds[1] = f->m_vertexList[1]->m_id;
594 faceVertIds[2] = f->m_vertexList[2]->m_id;
596 for (j = 0; j < f->m_edgeList.size(); ++j)
598 edge = f->m_edgeList[j];
602 if (edge->m_edgeNodes.size() > 0)
610 for (
int k = 0; k < N - 2; ++k)
612 edge->m_edgeNodes.push_back(
613 (*hoIt)->surfVerts[3 + j * (N - 2) + k]);
619 if (edge->m_n1->m_id != faceVertIds[j])
621 reverse(edge->m_edgeNodes.begin(),
622 edge->m_edgeNodes.end());
628 for (
int j = 3 + 3 * (N - 2); j < Ntot; ++j)
630 f->m_faceNodes.push_back((*hoIt)->surfVerts[j]);
642 map<int, vector<pair<int, LibUtilities::ShapeType> > > surfaceCompMap;
655 if (line.find(
"*") != string::npos ||
m_mshFile.eof() ||
666 s >> bcType >> elId >> faceId;
668 elId = elMap[elId - 1];
671 vector<ConditionType> type;
677 if ((elm->GetDim() == 2 && faceId >= elm->GetEdgeCount()) ||
678 (elm->GetDim() == 3 && faceId >= elm->GetFaceCount()))
698 for (i = 0; i <
m_mesh->m_fields.size() - 1; ++i)
718 size_t p = line.find_first_of(
'=');
720 boost::algorithm::trim_copy(line.substr(p + 1)));
725 for (i = 0; i <
m_mesh->m_fields.size() - 1; ++i)
728 size_t p = line.find_first_of(
'=');
730 boost::algorithm::trim_copy(line.substr(p + 1)));
751 for (i = 0; i <
m_mesh->m_fields.size(); ++i)
772 cerr <<
"Unknown boundary condition type " << line[0] << endl;
777 c->field =
m_mesh->m_fields;
786 for (it =
m_mesh->m_condition.begin(); it !=
m_mesh->m_condition.end();
796 int compTag, conditionId;
803 if (elm->GetDim() == 3)
809 boost::shared_ptr<Prism>
p =
810 boost::dynamic_pointer_cast<
Prism>(elm);
811 if (p->m_orientation == 1)
813 faceId = (faceId + 2) % 6;
815 else if (p->m_orientation == 2)
817 faceId = (faceId + 4) % 6;
822 boost::shared_ptr<Tetrahedron> t =
828 bool tri = f->m_vertexList.size() == 3;
830 vector<NodeSharedPtr> nodeList;
831 nodeList.insert(nodeList.begin(),
832 f->m_vertexList.begin(),
833 f->m_vertexList.end());
845 for (
int i = 0; i < f->m_vertexList.size(); ++i)
847 surfEl->GetEdge(i)->m_edgeNodes = f->m_edgeList[i]->m_edgeNodes;
848 surfEl->GetEdge(i)->m_curveType = f->m_edgeList[i]->m_curveType;
851 else if (faceId < elm->GetEdgeCount())
855 vector<NodeSharedPtr> nodeList;
856 nodeList.push_back(f->m_n1);
857 nodeList.push_back(f->m_n2);
883 conditionId =
m_mesh->m_condition.size();
884 compTag = nComposite;
885 c->m_composite.push_back(compTag);
886 m_mesh->m_condition[conditionId] = c;
888 surfaceCompMap[conditionId].push_back(
889 pair<int, LibUtilities::ShapeType>(nComposite, surfElType));
896 map<int, vector<pair<int, LibUtilities::ShapeType> > >
::iterator
898 it2 = surfaceCompMap.find(it->first);
901 if (it2 == surfaceCompMap.end())
904 cerr <<
"Unable to find condition!" << endl;
908 for (j = 0; j < it2->second.size(); ++j)
910 pair<int, LibUtilities::ShapeType> tmp = it2->second[j];
911 if (tmp.second == surfElType)
925 it2->second.push_back(
926 pair<int, LibUtilities::ShapeType>(nComposite, surfElType));
927 compTag = nComposite;
928 m_mesh->m_condition[it->first]->m_composite.push_back(compTag);
932 conditionId = it->first;
937 vector<int> existingTags = surfEl->GetTagList();
939 existingTags.insert(existingTags.begin(), compTag);
940 surfEl->SetTagList(existingTags);
941 surfEl->SetId(nSurfaces);
943 m_mesh->m_element[surfEl->GetDim()].push_back(surfEl);
961 map<string, pair<NekCurve, string> >
::iterator it;
962 int nodeId =
m_mesh->GetNumEntities();
967 string line, fileName = it->second.second;
971 if (it->second.first !=
eFile)
977 dot = fileName.find_last_of(
'.');
978 fileName = fileName.substr(0, dot);
982 hsf.open(fileName.c_str());
985 cerr <<
"Could not open surface file " << fileName << endl;
992 pos = line.find(
"=");
993 if (pos == string::npos)
995 cerr <<
"hsf header error: cannot read number of "
996 <<
"nodal points." << endl;
999 line = line.substr(pos + 1);
1000 stringstream ss(line);
1003 pos = line.find(
"=");
1004 if (pos == string::npos)
1006 cerr <<
"hsf header error: cannot read number of "
1007 <<
"faces." << endl;
1010 line = line.substr(pos + 1);
1015 int Ntot = N * (N + 1) / 2;
1022 for (
int i = 0; i < 2; ++i)
1033 cerr <<
"hsf header error: cannot read in "
1034 <<
"r/s points" << endl;
1038 for (
int j = 0; j < Ntot; ++j)
1040 ss >> (i == 0 ? r[j] : s[j]);
1054 for (
int i = 0; i < Ntot; ++i)
1056 for (
int j = 0; j < Ntot; ++j)
1058 if (fabs(r[i] - rp[j]) < 1e-5 && fabs(s[i] - sp[j]) < 1e-5)
1070 map<int, vector<NodeSharedPtr> > faceMap;
1071 for (
int i = 0; i < Nface; ++i)
1074 vector<NodeSharedPtr> faceNodes(Ntot);
1075 for (
int j = 0; j < Ntot; ++j, ++nodeId)
1085 for (
int j = 0; j < (N - 1) * (N - 1); ++j)
1089 faceMap[i] = faceNodes;
1095 for (
int i = 0; i < Nface; ++i)
1099 vector<int> nodeIds(3);
1104 ss >> tmp >> fid >> nodeIds[0] >> nodeIds[1] >> nodeIds[2];
1108 cerr <<
"Unable to read hsf connectivity information." << endl;
1112 hoData[it->first].insert(
1128 switch (InputNekEntity)
1155 cerr <<
"unknown Nektar element type" << endl;
A 3-dimensional four-faced element.
Basic information about an element.
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.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
A 3-dimensional five-faced element (2 triangles, 3 quadrilaterals).
boost::shared_ptr< HOSurf > HOSurfSharedPtr
pair< ModuleType, string > ModuleKey
MeshSharedPtr m_mesh
Mesh object.
ElementFactory & GetElementFactory()
NEKMESHUTILS_EXPORT NodeSharedPtr GetVertex(unsigned int i) const
Access a vertex node.
virtual NEKMESHUTILS_EXPORT void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
virtual NEKMESHUTILS_EXPORT void ProcessElements()
Generate element IDs.
Represents a command-line configuration option.
boost::shared_ptr< Node > NodeSharedPtr
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< Condition > ConditionSharedPtr
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.
virtual NEKMESHUTILS_EXPORT void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
boost::shared_ptr< Element > ElementSharedPtr
boost::shared_ptr< Face > FaceSharedPtr
std::pair< ModuleType, std::string > ModuleKey
HOTriangle< NodeSharedPtr > HOSurf
virtual NEKMESHUTILS_EXPORT void ProcessComposites()
Generate composites.
1D Gauss-Lobatto-Legendre quadrature points
2D Nodal Electrostatic Points on a Triangle
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.