41 #include <boost/algorithm/string.hpp>
57 InputNek5000::create,
"Reads Nektar rea file.");
82 int nParam, nElements, nCurves;
83 int i, j, k, nodeCounter = 0;
93 cout <<
"InputNek5000: Start reading file..." << endl;
99 for (i = 0; i < 4; ++i)
104 stringstream s(line);
107 for (i = 0; i < nParam; ++i)
120 for (i = 0; i < j; ++i)
130 for (i = 0; i < j; ++i)
138 bool foundMesh =
false;
142 if (line.find(
"MESH") != string::npos)
151 cerr <<
"Couldn't find MESH tag inside file." << endl;
159 s >> nElements >>
m_mesh->m_expDim;
163 m_mesh->m_fields.push_back(
"u");
164 m_mesh->m_fields.push_back(
"v");
165 if (
m_mesh->m_spaceDim > 2)
167 m_mesh->m_fields.push_back(
"w");
169 m_mesh->m_fields.push_back(
"p");
172 for (i = 0; i < nElements; ++i)
177 if (
m_mesh->m_expDim == 2)
182 for (j = 0; j < 2; ++j)
187 for (k = 0; k < 4; ++k)
199 for (j = 0; j < 6; ++j)
204 int offset = j > 2 ? 4 : 0;
205 for (k = 0; k < 4; ++k)
207 s >> vertex[offset + k][j % 3];
215 vector<NodeSharedPtr> nodeList(nNodes);
216 for (k = 0; k < nNodes; ++k)
218 nodeList[k] = boost::shared_ptr<Node>(
220 0, vertex[k][0], vertex[k][1], vertex[k][2]));
221 pair<NodeSet::iterator, bool> testIns =
222 m_mesh->m_vertexSet.insert(nodeList[k]);
226 nodeList[k] = *(testIns.first);
230 nodeList[k]->m_id = nodeCounter++;
234 vector<int> tags(1, 0);
237 elType, conf, nodeList, tags);
238 m_mesh->m_element[E->GetDim()].push_back(E);
243 if (line.find(
"CURVE") == string::npos)
245 cerr <<
"Cannot find curved side data." << endl;
263 int nek2nekedge[12] = {
264 0, 1, 2, 3, 8, 9, 10, 11, 4, 5, 6, 7
270 int nek2nekface[6] = {
276 for (i = 0; i < nCurves; ++i)
284 if (nElements < 1000)
287 s.str(line.substr(0, 3));
290 s.str(line.substr(3, 3));
292 line = line.substr(6);
294 else if (nElements < 1000000)
297 s.str(line.substr(0, 2));
300 s.str(line.substr(2, 6));
302 line = line.substr(8);
307 s.str(line.substr(0, 2));
310 s.str(line.substr(2, 12));
312 line = line.substr(14);
318 for (j = 0; j < 5; ++j)
326 side = nek2nekedge[side];
337 int convexity = radius < 0 ? -1 : 1;
338 radius = fabs(radius);
346 Node P1(*(edge->m_n1)), P2(*(edge->m_n2));
348 if (fabs(P1.m_z - P2.
m_z) > 1e-8)
350 cout <<
"warning: detected non x-y edge." << endl;
352 P1.m_z = P2.
m_z = 0.0;
354 Node unitNormal, link, centroid, centre;
355 Node midpoint = (P1 + P2)*0.5, dx = P2 - P1;
358 unitNormal.
m_x = -dx.m_y / l;
359 unitNormal.
m_y = dx.m_x / l;
361 if (2.0 * radius < l)
363 cerr <<
"failure" << endl;
367 semiangle = asin (0.5 * l / radius);
371 vector<NodeSharedPtr> elNodes = el->GetVertexList();
372 int nNodes = elNodes.size();
374 for (
int i = 0; i < nNodes; ++i)
377 Node tmp(*elNodes[i]);
383 link = centroid - midpoint;
386 centre = midpoint + unitNormal * (
sign * cos(semiangle) *
390 theta1 = atan2 (P1.m_y - centre.
m_y, P1.m_x - centre.
m_x);
391 theta2 = atan2 (P2.
m_y - centre.
m_y, P2.
m_x - centre.
m_x);
392 dtheta = theta2 - theta1;
394 if (fabs(dtheta) > 2.0*semiangle + 1e-15)
396 dtheta += (dtheta < 0.0) ? 2.0*M_PI : -2.0*M_PI;
399 edge->m_edgeNodes.clear();
401 for (j = 1; j < nq-1; ++j) {
402 phi = theta1 + dtheta * 0.5 * (rp[j] + 1.0);
405 centre.
m_x + radius * cos(phi),
406 centre.
m_y + radius * sin(phi),
408 edge->m_edgeNodes.push_back(asd);
416 cerr <<
"Curve type '" << curveType <<
"' on side " << side
417 <<
" of element " << elmt <<
" is unsupported;"
418 <<
"will ignore." << endl;
421 cerr <<
"Unknown curve type '" << curveType <<
"' on side "
422 << side <<
" of element " << elmt <<
"; will ignore."
432 if (line.find(
"BOUNDARY") == string::npos)
434 cerr <<
"Cannot find boundary conditions." << endl;
439 boost::unordered_set<pair<int, int> > periodicIn;
440 int periodicInId = -1, periodicOutId = -1;
445 int perIn = 0, perOut = 0;
453 if (line.find(
"*") != string::npos)
465 s.str(line.substr(0, 4));
475 if (nElements < 1000)
479 s.str(line.substr(4, 3));
482 s.str(line.substr(7, 3));
484 line = line.substr(10);
486 else if (nElements < 100000)
490 s.str(line.substr(4, 5));
493 s.str(line.substr(9, 1));
495 line = line.substr(10);
497 else if (nElements < 1000000)
501 s.str(line.substr(4, 6));
503 side = lineCnt % (2 *
m_mesh->m_expDim);
504 line = line.substr(9);
510 s.str(line.substr(4, 12));
512 side = lineCnt % (2 *
m_mesh->m_expDim);
513 line = line.substr(15);
519 for (i = 0; i < 5; ++i)
533 std::string fields[] = {
"u",
"v",
"w",
"p" };
544 for (i = 0; i <
m_mesh->m_fields.size() - 1; ++i)
546 c->field.push_back(fields[i]);
547 c->value.push_back(
"0");
552 c->field.push_back(fields[3]);
553 c->value.push_back(
"0");
561 int perElmt = (int)(data[0] + 0.5) - 1;
562 int perFace = (int)(data[1] + 0.5) - 1;
565 if (periodicInId == -1)
567 periodicInId =
m_mesh->m_condition.size();
568 periodicOutId =
m_mesh->m_condition.size()+1;
572 bool hasIn = periodicIn.find(make_pair(perElmt, perFace)) !=
577 swap(periodicInId, periodicOutId);
582 periodicIn.insert(make_pair(elmt, side));
586 std::string periodicInStr =
"[" +
587 boost::lexical_cast<
string>(periodicInId) +
"]";
588 std::string periodicOutStr =
"[" +
589 boost::lexical_cast<
string>(periodicOutId) +
"]";
591 for (i = 0; i <
m_mesh->m_fields.size() - 1; ++i)
593 c->field.push_back(fields[i]);
594 c->value.push_back(periodicOutStr);
598 c->field.push_back(fields[3]);
599 c->value.push_back(periodicOutStr);
607 c->m_composite.push_back(nComposite++);
608 c2->m_composite.push_back(nComposite++);
610 c2->field = c->field;
612 for (i = 0; i < c->type.size(); ++i)
614 c2->value.push_back(periodicInStr);
617 m_mesh->m_condition[periodicInId] = c;
618 m_mesh->m_condition[periodicOutId] = c2;
623 swap(periodicInId, periodicOutId);
633 int compTag, conditionId;
637 if (el->GetDim() == 3)
640 vector<NodeSharedPtr> nodeList;
641 nodeList.insert(nodeList.begin(),
642 f->m_vertexList.begin(),
643 f->m_vertexList.end());
651 conf, nodeList, tags);
654 for (
int i = 0; i < f->m_vertexList.size(); ++i)
656 surfEl->GetEdge(i)->m_edgeNodes = f->m_edgeList[i]->m_edgeNodes;
657 surfEl->GetEdge(i)->m_curveType = f->m_edgeList[i]->m_curveType;
664 vector<NodeSharedPtr> nodeList;
665 nodeList.push_back(f->m_n1);
666 nodeList.push_back(f->m_n2);
682 for (it =
m_mesh->m_condition.begin(); it !=
m_mesh->m_condition.end();
695 conditionId =
m_mesh->m_condition.size();
696 compTag = nComposite;
697 c->m_composite.push_back(compTag);
698 m_mesh->m_condition[conditionId] = c;
702 compTag = c->m_composite[0];
707 vector<int> existingTags = surfEl->GetTagList();
709 existingTags.insert(existingTags.begin(), compTag);
710 surfEl->SetTagList(existingTags);
711 surfEl->SetId(nSurfaces);
713 m_mesh->m_element[surfEl->GetDim()].push_back(surfEl);
717 if (lineCnt != nElements * (
m_mesh->m_expDim * 2))
719 cerr <<
"Warning: boundary conditions may not have been correctly read "
720 <<
"from Nek5000 input file." << endl;
725 cerr <<
"Warning: number of periodic faces does not match." << endl;
737 if (periodicInId != -1)
740 ->m_composite[0]]->m_reorder =
false;
742 ->m_composite[0]]->m_reorder =
false;
NEKMESHUTILS_EXPORT NekDouble dot(const Node &pSrc) const
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.
#define sign(a, b)
return the sign(b)*a
pair< ModuleType, string > ModuleKey
NekDouble m_y
Y-coordinate.
MeshSharedPtr m_mesh
Mesh object.
ElementFactory & GetElementFactory()
Represents a point in the domain.
virtual NEKMESHUTILS_EXPORT void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
virtual NEKMESHUTILS_EXPORT 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
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
NekDouble m_x
X-coordinate.
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
NekDouble m_z
Z-coordinate.
std::pair< ModuleType, std::string > ModuleKey
virtual NEKMESHUTILS_EXPORT void ProcessComposites()
Generate composites.
1D Gauss-Lobatto-Legendre quadrature points
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.