41 #include <boost/algorithm/string.hpp>
83 int nParam, nElements, nCurves;
84 int i, j, k, nodeCounter = 0;
88 map<LibUtilities::ShapeType, int> domainComposite;
89 map<LibUtilities::ShapeType, vector<vector<NodeSharedPtr> > > elNodes;
90 map<LibUtilities::ShapeType, vector<int> > elIds;
91 boost::unordered_map<int, int> elMap;
92 vector<LibUtilities::ShapeType> elmOrder;
108 cout <<
"InputNek: Start reading file..." << endl;
114 for (i = 0; i < 4; ++i)
119 stringstream s(line);
122 for (i = 0; i < nParam; ++i)
135 for (i = 0; i < j; ++i)
145 for (i = 0; i < j; ++i)
153 bool foundMesh =
false;
157 if (line.find(
"MESH") != string::npos)
166 cerr <<
"Couldn't find MESH tag inside file." << endl;
174 s >> nElements >>
m_mesh->m_expDim;
178 m_mesh->m_fields.push_back(
"u");
179 m_mesh->m_fields.push_back(
"v");
180 if (
m_mesh->m_spaceDim > 2)
182 m_mesh->m_fields.push_back(
"w");
184 m_mesh->m_fields.push_back(
"p");
187 for (i = 0; i < nElements; ++i)
191 if (
m_mesh->m_expDim == 2)
193 if (line.find(
"Qua") != string::npos ||
194 line.find(
"qua") != string::npos)
206 if (line.find(
"Tet") != string::npos ||
207 line.find(
"tet") != string::npos)
211 else if (line.find(
"Hex") != string::npos ||
212 line.find(
"hex") != string::npos)
216 else if (line.find(
"Prism") != string::npos ||
217 line.find(
"prism") != string::npos)
221 else if (line.find(
"Pyr") != string::npos ||
222 line.find(
"pyr") != string::npos)
226 else if (line.find(
"Qua") != string::npos ||
227 line.find(
"qua") != string::npos)
241 for (j = 0; j <
m_mesh->m_expDim; ++j)
246 for (k = 0; k < nNodes; ++k)
253 for (j =
m_mesh->m_expDim; j < 3; ++j)
255 for (k = 0; k < nNodes; ++k)
264 vector<NodeSharedPtr> nodeList;
265 for (k = 0; k < nNodes; ++k)
268 nodeCounter++, vertex[0][k], vertex[1][k], vertex[2][k]));
269 nodeList.push_back(n);
272 elNodes[elType].push_back(nodeList);
273 elIds[elType].push_back(i);
279 for (i = 0; i < elmOrder.size(); ++i)
282 vector<vector<NodeSharedPtr> > &tmp = elNodes[elType];
284 for (j = 0; j < tmp.size(); ++j)
288 domainComposite.find(elType);
289 if (compIt == domainComposite.end())
291 tags.push_back(nComposite);
292 domainComposite[elType] = nComposite;
297 tags.push_back(compIt->second);
300 elMap[elIds[elType][j]] = reorderedId++;
302 vector<NodeSharedPtr> nodeList = tmp[j];
304 for (k = 0; k < nodeList.size(); ++k)
306 pair<NodeSet::iterator, bool> testIns =
307 m_mesh->m_vertexSet.insert(nodeList[k]);
311 nodeList[k] = *(testIns.first);
315 nodeList[k]->m_id = nodeCounter++;
322 elType, conf, nodeList, tags);
323 m_mesh->m_element[E->GetDim()].push_back(E);
329 if (line.find(
"CURVE") == string::npos)
331 cerr <<
"Cannot find curved side data." << endl;
345 for (i = 0; i < nCurves; ++i)
358 s >> word >> curveTag;
361 else if (word ==
"Recon")
367 s >> word >> curveTag;
372 cerr <<
"Unsupported curve type " << word << endl;
384 if (line.find(
"side") == string::npos)
386 cerr <<
"Unable to read number of curved sides" << endl;
392 map<string, pair<NekCurve, string> >
::iterator it;
401 for (i = 0; i < nCurvedSides; ++i)
406 s >> faceId >> elId >> word;
408 elId = elMap[elId - 1];
411 int origFaceId = faceId;
415 boost::shared_ptr<Prism> p =
416 boost::dynamic_pointer_cast<
Prism>(el);
417 if (p->m_orientation == 1)
419 faceId = (faceId + 2) % 6;
421 else if (p->m_orientation == 2)
423 faceId = (faceId + 4) % 6;
428 boost::shared_ptr<Tetrahedron> t =
436 cerr <<
"Unrecognised curve tag " << word <<
" in curved lines"
441 if (it->second.first ==
eRecon)
444 vector<NodeSharedPtr> &tmp = el->GetFace(faceId)->m_vertexList;
445 vector<Node> n(tmp.size());
452 boost::dynamic_pointer_cast<
Prism>(el)->m_orientation;
459 for (j = 0; j < tmp.size(); ++j)
467 for (j = 0; j < tmp.size(); ++j)
475 for (j = 0; j < tmp.size(); ++j)
480 for (j = 0; j < tmp.size(); ++j)
482 int id = tmp[(j + offset) % tmp.size()]->m_id;
484 m_mesh->m_vertexNormals.find(
id);
486 if (vIt ==
m_mesh->m_vertexNormals.end())
488 m_mesh->m_vertexNormals[id] = n[j];
494 m_mesh->m_spherigonSurfs.insert(make_pair(elId, faceId));
496 else if (it->second.first ==
eFile)
499 static int tetFaceVerts[4][3] = {
500 {0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};
502 vector<int> vertId(3);
503 s >> vertId[0] >> vertId[1] >> vertId[2];
510 boost::shared_ptr<Tetrahedron> tet =
512 vector<int> tmpVertId = vertId;
514 for (j = 0; j < 3; ++j)
518 [tetFaceVerts[origFaceId][j]])
521 for (k = 0; k < 3; ++k)
523 int w = f->m_vertexList[k]->m_id;
526 vertId[k] = tmpVertId[j];
534 boost::shared_ptr<Prism> pr =
535 boost::static_pointer_cast<
Prism>(el);
536 if (pr->m_orientation == 1)
538 swap(vertId[2], vertId[1]);
539 swap(vertId[0], vertId[1]);
541 else if (pr->m_orientation == 2)
543 swap(vertId[0], vertId[1]);
544 swap(vertId[2], vertId[1]);
549 boost::shared_ptr<HOSurf>(
new HOSurf(vertId));
551 hoIt =
hoData[word].find(hs);
553 if (hoIt ==
hoData[word].end())
555 cerr <<
"Unable to find high-order surface data "
556 <<
"for element id " << elId + 1 << endl;
566 int Ntot = (*hoIt)->surfVerts.size();
567 int N = ((int)sqrt(8.0 * Ntot + 1.0) - 1) / 2;
572 vector<NodeSharedPtr> tmpVerts = (*hoIt)->surfVerts;
573 for (j = 0; j < tmpVerts.size(); ++j)
575 (*hoIt)->surfVerts[
hoMap[j]] = tmpVerts[j];
578 for (j = 0; j < tmpVerts.size(); ++j)
583 vector<int> faceVertIds(3);
584 faceVertIds[0] = f->m_vertexList[0]->m_id;
585 faceVertIds[1] = f->m_vertexList[1]->m_id;
586 faceVertIds[2] = f->m_vertexList[2]->m_id;
588 for (j = 0; j < f->m_edgeList.size(); ++j)
590 edge = f->m_edgeList[j];
594 if (edge->m_edgeNodes.size() > 0)
602 for (
int k = 0; k < N - 2; ++k)
604 edge->m_edgeNodes.push_back(
605 (*hoIt)->surfVerts[3 + j * (N - 2) + k]);
611 if (edge->m_n1->m_id != faceVertIds[j])
613 reverse(edge->m_edgeNodes.begin(),
614 edge->m_edgeNodes.end());
620 for (
int j = 3 + 3 * (N - 2); j < Ntot; ++j)
622 f->m_faceNodes.push_back((*hoIt)->surfVerts[j]);
634 map<int, vector<pair<int, LibUtilities::ShapeType> > > surfaceCompMap;
647 if (line.find(
"*") != string::npos ||
m_mshFile.eof() ||
658 s >> bcType >> elId >> faceId;
660 elId = elMap[elId - 1];
663 vector<ConditionType> type;
674 for (i = 0; i <
m_mesh->m_fields.size() - 1; ++i)
690 for (i = 0; i <
m_mesh->m_fields.size() - 1; ++i)
693 size_t p = line.find_first_of(
'=');
695 boost::algorithm::trim_copy(line.substr(p + 1)));
708 for (i = 0; i <
m_mesh->m_fields.size(); ++i)
728 cerr <<
"Unknown boundary condition type " << line[0] << endl;
733 c->field =
m_mesh->m_fields;
742 for (it =
m_mesh->m_condition.begin(); it !=
m_mesh->m_condition.end();
752 int compTag, conditionId;
760 if (elm->GetDim() == 3)
766 boost::shared_ptr<Prism> p =
767 boost::dynamic_pointer_cast<
Prism>(elm);
768 if (p->m_orientation == 1)
770 faceId = (faceId + 2) % 6;
772 else if (p->m_orientation == 2)
774 faceId = (faceId + 4) % 6;
779 boost::shared_ptr<Tetrahedron> t =
785 bool tri = f->m_vertexList.size() == 3;
787 vector<NodeSharedPtr> nodeList;
788 nodeList.insert(nodeList.begin(),
789 f->m_vertexList.begin(),
790 f->m_vertexList.end());
802 for (
int i = 0; i < f->m_vertexList.size(); ++i)
804 surfEl->GetEdge(i)->m_edgeNodes = f->m_edgeList[i]->m_edgeNodes;
805 surfEl->GetEdge(i)->m_curveType = f->m_edgeList[i]->m_curveType;
812 vector<NodeSharedPtr> nodeList;
813 nodeList.push_back(f->m_n1);
814 nodeList.push_back(f->m_n2);
835 conditionId =
m_mesh->m_condition.size();
836 compTag = nComposite;
837 c->m_composite.push_back(compTag);
838 m_mesh->m_condition[conditionId] = c;
840 surfaceCompMap[conditionId].push_back(
841 pair<int, LibUtilities::ShapeType>(nComposite, surfElType));
848 map<int, vector<pair<int, LibUtilities::ShapeType> > >
::iterator
850 it2 = surfaceCompMap.find(it->first);
853 if (it2 == surfaceCompMap.end())
856 cerr <<
"Unable to find condition!" << endl;
860 for (j = 0; j < it2->second.size(); ++j)
862 pair<int, LibUtilities::ShapeType> tmp = it2->second[j];
863 if (tmp.second == surfElType)
877 it2->second.push_back(
878 pair<int, LibUtilities::ShapeType>(nComposite, surfElType));
879 compTag = nComposite;
880 m_mesh->m_condition[it->first]->m_composite.push_back(compTag);
884 conditionId = it->first;
889 vector<int> existingTags = surfEl->GetTagList();
891 existingTags.insert(existingTags.begin(), compTag);
892 surfEl->SetTagList(existingTags);
893 surfEl->SetId(nSurfaces);
895 m_mesh->m_element[surfEl->GetDim()].push_back(surfEl);
913 map<string, pair<NekCurve, string> >
::iterator it;
914 int nodeId =
m_mesh->GetNumEntities();
919 string line, fileName = it->second.second;
923 if (it->second.first !=
eFile)
929 dot = fileName.find_last_of(
'.');
930 fileName = fileName.substr(0, dot);
934 hsf.open(fileName.c_str());
937 cerr <<
"Could not open surface file " << fileName << endl;
944 pos = line.find(
"=");
945 if (pos == string::npos)
947 cerr <<
"hsf header error: cannot read number of "
948 <<
"nodal points." << endl;
951 line = line.substr(pos + 1);
952 stringstream ss(line);
955 pos = line.find(
"=");
956 if (pos == string::npos)
958 cerr <<
"hsf header error: cannot read number of "
962 line = line.substr(pos + 1);
967 int Ntot = N * (N + 1) / 2;
974 for (
int i = 0; i < 2; ++i)
985 cerr <<
"hsf header error: cannot read in "
986 <<
"r/s points" << endl;
990 for (
int j = 0; j < Ntot; ++j)
992 ss >> (i == 0 ? r[j] : s[j]);
1006 for (
int i = 0; i < Ntot; ++i)
1008 for (
int j = 0; j < Ntot; ++j)
1010 if (fabs(r[i] - rp[j]) < 1e-5 && fabs(s[i] - sp[j]) < 1e-5)
1022 map<int, vector<NodeSharedPtr> > faceMap;
1023 for (
int i = 0; i < Nface; ++i)
1026 vector<NodeSharedPtr> faceNodes(Ntot);
1027 for (
int j = 0; j < Ntot; ++j, ++nodeId)
1037 for (
int j = 0; j < (N - 1) * (N - 1); ++j)
1041 faceMap[i] = faceNodes;
1047 for (
int i = 0; i < Nface; ++i)
1051 vector<int> nodeIds(3);
1056 ss >> tmp >> fid >> nodeIds[0] >> nodeIds[1] >> nodeIds[2];
1060 cerr <<
"Unable to read hsf connectivity information." << endl;
1064 hoData[it->first].insert(
1080 switch (InputNekEntity)
1107 cerr <<
"unknown Nektar element type" << endl;
A 3-dimensional four-faced element.
Basic information about an element.
pair< ModuleType, string > ModuleKey
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
MeshSharedPtr m_mesh
Mesh object.
ElementFactory & GetElementFactory()
Represents a point in the domain.
virtual void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
NEKMESHUTILS_EXPORT NodeSharedPtr GetVertex(unsigned int i) const
Access a vertex node.
virtual void ProcessElements()
Generate element IDs.
boost::shared_ptr< Node > NodeSharedPtr
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
boost::shared_ptr< Condition > ConditionSharedPtr
virtual void ProcessComposites()
Generate composites.
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 void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
boost::shared_ptr< Element > ElementSharedPtr
boost::shared_ptr< Face > FaceSharedPtr
Shared pointer to a face.
HOTriangle< NodeSharedPtr > HOSurf
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.