53 Loki::SingleThreaded> Type;
54 return Type::Instance();
80 string fname =
m_config[
"infile"].as<
string>();
84 cerr <<
"Error opening file: " << fname << endl;
94 string fname =
m_config[
"outfile"].as<
string>();
98 cerr <<
"Error opening file: " << fname << endl;
113 vector<ElementSharedPtr> &elmt =
m_mesh->m_element[
m_mesh->m_expDim];
115 m_mesh->m_vertexSet.clear();
117 for (
int i = 0, vid = 0; i < elmt.size(); ++i)
119 for (
int j = 0; j < elmt[i]->GetVertexCount(); ++j)
121 pair<NodeSet::iterator,bool> testIns =
122 m_mesh->m_vertexSet.insert(elmt[i]->GetVertex(j));
126 (*(testIns.first))->m_id = vid++;
130 elmt[i]->SetVertex(j,*testIns.first);
153 if (
m_mesh->m_expDim < 2)
return;
157 vector<ElementSharedPtr> &elmt =
m_mesh->m_element[
m_mesh->m_expDim];
159 m_mesh->m_edgeSet.clear();
164 for (
int i = 0, eid = 0; i < elmt.size(); ++i)
166 for (
int j = 0; j < elmt[i]->GetEdgeCount(); ++j)
168 pair<EdgeSet::iterator,bool> testIns;
170 testIns =
m_mesh->m_edgeSet.insert(ed);
176 ed2->m_elLink.push_back(
177 pair<ElementSharedPtr,int>(elmt[i],j));
182 elmt[i]->SetEdge(j, e2);
183 if (e2->m_edgeNodes.size() == 0 &&
184 ed->m_edgeNodes.size() > 0)
186 e2->m_curveType = ed->m_curveType;
187 e2->m_edgeNodes = ed->m_edgeNodes;
190 if (e2->m_n1->m_id != ed->m_n1->m_id)
192 reverse(e2->m_edgeNodes.begin(),
193 e2->m_edgeNodes.end());
197 #ifdef NEKTAR_USE_MESHGEN
200 e2->onCurve = ed->onCurve;
201 e2->CADCurveId = ed->CADCurveId;
202 e2->CADCurve = ed->CADCurve;
207 e2->m_elLink.push_back(
208 pair<ElementSharedPtr,int>(elmt[i],j));
215 for (
int i = 0; i <
m_mesh->m_element[1].size(); ++i)
219 vector<NodeSharedPtr> edgeNodes;
221 new Edge(v0, v1, edgeNodes,
222 m_mesh->m_element[1][i]->GetConf().m_edgeCurveType));
224 if (it ==
m_mesh->m_edgeSet.end())
226 cerr <<
"Cannot find corresponding element edge for "
227 <<
"1D element " << i << endl;
230 m_mesh->m_element[1][i]->SetEdgeLink(*it);
233 pair<ElementSharedPtr, int> eMap = (*it)->m_elLink.at(0);
234 eMap.first->SetBoundaryLink(eMap.second, i);
237 if ((*it)->m_edgeNodes.size() > 0)
240 if (edge->GetVertex(0) == (*it)->m_n1)
242 edge->SetVolumeNodes((*it)->m_edgeNodes);
266 if (
m_mesh->m_expDim < 3)
return;
270 vector<ElementSharedPtr> &elmt =
m_mesh->m_element[
m_mesh->m_expDim];
272 m_mesh->m_faceSet.clear();
275 for (
int i = 0, fid = 0; i < elmt.size(); ++i)
277 for (
int j = 0; j < elmt[i]->GetFaceCount(); ++j)
279 pair<FaceSet::iterator,bool> testIns;
280 testIns =
m_mesh->m_faceSet.insert(elmt[i]->GetFace(j));
284 (*(testIns.first))->m_id = fid++;
285 (*(testIns.first))->m_elLink.push_back(
286 pair<ElementSharedPtr,int>(elmt[i],j));
290 elmt[i]->SetFace(j,*testIns.first);
292 (*(testIns.first))->m_elLink.push_back(
293 pair<ElementSharedPtr,int>(elmt[i],j));
300 for (
int i = 0; i <
m_mesh->m_element[2].size(); ++i)
302 vector<NodeSharedPtr> vertices =
m_mesh->m_element[2][i]->GetVertexList();
303 vector<NodeSharedPtr> faceNodes;
304 vector<EdgeSharedPtr> edgeList =
m_mesh->m_element[2][i]->GetEdgeList();
306 new Face(vertices, faceNodes, edgeList,
307 m_mesh->m_element[2][i]->GetConf().m_faceCurveType));
309 if (it ==
m_mesh->m_faceSet.end())
311 cout <<
"Cannot find corresponding element face for 2D "
312 <<
"element " << i << endl;
315 m_mesh->m_element[2][i]->SetFaceLink(*it);
318 pair<ElementSharedPtr, int> eMap = (*it)->m_elLink.at(0);
319 eMap.first->SetBoundaryLink(eMap.second, i);
334 for (
int i = 0; i <
m_mesh->m_element[
m_mesh->m_expDim].size(); ++i)
351 m_mesh->m_composite.clear();
356 for (
int d = 0; d <=
m_mesh->m_expDim; ++d)
358 vector<ElementSharedPtr> &elmt =
m_mesh->m_element[d];
360 for (
int i = 0; i < elmt.size(); ++i)
363 unsigned int tagid = elmt[i]->GetTagList()[0];
365 it =
m_mesh->m_composite.find(tagid);
367 if (it ==
m_mesh->m_composite.end())
371 pair<CompositeMap::iterator, bool> testIns;
373 tmp->m_tag = elmt[i]->GetTag();
374 if(
m_mesh->m_faceLabels.count(tmp->m_id) != 0)
376 tmp->m_label =
m_mesh->m_faceLabels[tmp->m_id];
379 testIns =
m_mesh->m_composite.insert(
380 pair<unsigned int, CompositeSharedPtr>(tagid,tmp));
384 if (elmt[i]->GetTag() != it->second->m_tag)
386 cout <<
"Different types of elements in same composite!" << endl;
387 cout <<
" -> Composite uses " << it->second->m_tag << endl;
388 cout <<
" -> Element uses " << elmt[i]->GetTag() << endl;
389 cout <<
"Have you specified physical volumes and surfaces?" << endl;
391 it->second->m_items.push_back(elmt[i]);
403 for(eit =
m_mesh->m_edgeSet.begin(); eit !=
m_mesh->m_edgeSet.end(); eit++)
405 (*eit)->m_elLink.clear();
410 for(fit =
m_mesh->m_faceSet.begin(); fit !=
m_mesh->m_faceSet.end(); fit++)
412 (*fit)->m_elLink.clear();
453 set<int> prismsDone, tetsDone;
457 for (i = 0; i <
m_mesh->m_element[3].size(); ++i)
463 prismsDone.insert(i);
473 for (it =
m_mesh->m_vertexSet.begin(); it !=
m_mesh->m_vertexSet.end(); ++it)
480 int prismTris[2][3] = {{0,1,4}, {3,2,5}};
483 bool warnCurvature =
false;
487 boost::unordered_set<int> facesDone;
491 while (prismsDone.size() > 0)
493 vector<ElementSharedPtr> line;
497 PrismLines(*prismsDone.begin(), perFaces, prismsDone, line);
501 for (i = 0; i < line.size(); ++i)
504 vector<int> tags = line[i]->GetTagList();
505 vector<NodeSharedPtr> nodes = line[i]->GetVertexList();
510 line[i]->GetFace(1), line[i]->GetFace(3)
513 fIt[0] = facesDone.find(f[0]->m_id);
514 fIt[1] = facesDone.find(f[1]->m_id);
518 for (j = 0; j < 2; ++j)
520 pIt = perFaces.find(f[j]->m_id);
522 if (pIt == perFaces.end())
527 fIt2 = facesDone.find(pIt->second.first->m_id);
529 if (fIt[j] == facesDone.end() &&
530 fIt2 != facesDone.end())
536 if (fIt[0] != facesDone.end() &&
537 fIt[1] != facesDone.end())
541 ASSERTL0(
false,
"Renumbering error!");
543 else if (fIt[0] == facesDone.end() &&
544 fIt[1] == facesDone.end())
547 for (j = 0; j < 2; ++j)
549 for (k = 0; k < 3; ++k)
559 facesDone.insert(f[0]->m_id);
560 facesDone.insert(f[1]->m_id);
566 int t = fIt[0] == facesDone.end() ? 0 : 1;
568 ASSERTL1(fIt[o] != facesDone.end(),
"Renumbering error");
574 nodes[prismTris[o][0]]->m_id,
575 nodes[prismTris[o][1]]->m_id,
576 nodes[prismTris[o][2]]->m_id
578 int tmp2[3] = {0,1,2};
580 if (tmp1[0] > tmp1[1])
582 swap(tmp1[0], tmp1[1]);
583 swap(tmp2[0], tmp2[1]);
586 if (tmp1[1] > tmp1[2])
588 swap(tmp1[1], tmp1[2]);
589 swap(tmp2[1], tmp2[2]);
592 if (tmp1[0] > tmp1[2])
594 swap(tmp1[0], tmp1[2]);
595 swap(tmp2[0], tmp2[2]);
599 for (j = 0; j < 3; ++j)
608 facesDone.insert(f[t]->m_id);
611 for (j = 0; j < 6; ++j)
613 ASSERTL1(nodes[j]->m_id != -1,
"Renumbering error");
623 for (j = 0; j < 9; ++j)
626 for (k = 0; k < 9; ++k)
629 if (e1->m_n1 == e2->m_n1 && e1->m_n2 == e2->m_n2)
631 e2->m_edgeNodes = e1->m_edgeNodes;
633 else if (e1->m_n1 == e2->m_n1 && e1->m_n2 == e2->m_n2)
635 e2->m_edgeNodes = e1->m_edgeNodes;
636 std::reverse(e2->m_edgeNodes.begin(),
637 e2->m_edgeNodes.end());
645 for (j = 0; j < 5; ++j)
647 if (line[i]->GetFace(j)->m_faceNodes.size() > 0)
649 warnCurvature =
true;
656 m_mesh->m_element[3][line[i]->GetId()] = el;
662 cerr <<
"[ReorderPrisms] WARNING: Face curvature detected in "
663 <<
"some prisms; this will be ignored in further module "
669 for (pIt = perFaces.begin(); pIt != perFaces.end(); ++pIt)
673 vector<int> perVerts = pIt->second.second;
674 int nVerts = perVerts.size();
677 for (j = 0; j < nVerts; ++j)
682 if (n1->m_id == -1 && n2->m_id == -1)
687 else if (n1->m_id != -1 && n2->m_id != -1)
693 ASSERTL0(
false,
"Periodic face renumbering error");
700 for (it2 = tetsDone.begin(); it2 != tetsDone.end(); ++it2)
703 vector<NodeSharedPtr> nodes = el->GetVertexList();
704 vector<int> tags = el->GetTagList();
706 for (i = 0; i < 4; ++i)
708 if (nodes[i]->m_id == -1)
710 nodes[i]->m_id = nodeId++;
721 for (it =
m_mesh->m_vertexSet.begin(); it !=
m_mesh->m_vertexSet.end(); ++it)
723 if ((*it)->m_id == -1)
725 (*it)->m_id = nodeId++;
736 set<int> &prismsDone,
737 vector<ElementSharedPtr> &line)
743 if (it == prismsDone.end())
749 prismsDone.erase(it);
750 line.push_back(
m_mesh->m_element[3][prism]);
753 for (i = 1; i <= 3; i += 2)
759 it2 = perFaces.find(f->m_id);
761 if (it2 != perFaces.end())
763 int id2 = it2->second.first->m_id;
764 nextId = it2->second.first->m_elLink[0].first->GetId();
767 PrismLines(nextId, perFaces, prismsDone, line);
771 if (f->m_elLink.size() == 1)
776 nextId = f->m_elLink[0].first->GetId();
777 if (nextId ==
m_mesh->m_element[3][prism]->GetId())
779 nextId = f->m_elLink[1].first->GetId();
782 PrismLines(nextId, perFaces, prismsDone, line);
794 cerr <<
"WARNING: Unrecognised config option " << key
795 <<
", proceeding anyway." << endl;
798 it->second.beenSet =
true;
800 if (it->second.isBool)
802 it->second.value =
"1";
806 it->second.value = val;
819 cerr <<
"No configuration options for this module." << endl;
825 cerr << setw(10) << it->first <<
": " << it->second.desc
840 if (!it->second.beenSet)
842 it->second.value = it->second.defValue;
854 cerr <<
"Expansion dimension is " <<
m_mesh->m_expDim << endl;
855 cerr <<
"Space dimension is " <<
m_mesh->m_spaceDim << endl;
856 cerr <<
"Read " <<
m_mesh->m_node.size() <<
" nodes" << endl;
857 cerr <<
"Read " <<
m_mesh->GetNumElements() <<
" "
858 <<
m_mesh->m_expDim <<
"-D elements" << endl;
859 cerr <<
"Read " <<
m_mesh->GetNumBndryElements()
860 <<
" boundary elements" << endl;
LibUtilities::NekFactory< ModuleKey, Module, FieldSharedPtr > ModuleFactory
#define ASSERTL0(condition, msg)
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.
Represents an edge which joins two points.
Represents a face comprised of three or more edges.
map< string, ConfigOption > m_config
List of configuration values.
MeshSharedPtr m_mesh
Mesh object.
void RegisterConfig(string key, string value)
Register a configuration option with a module.
ElementFactory & GetElementFactory()
virtual void ClearElementLinks()
virtual void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
std::ostream & operator<<(std::ostream &os, const BasisKey &rhs)
virtual void ProcessVertices()
Extract element vertices.
virtual void ProcessElements()
Generate element IDs.
OutputModule(FieldSharedPtr p_f)
boost::shared_ptr< Composite > CompositeSharedPtr
Shared pointer to a composite.
boost::shared_ptr< Node > NodeSharedPtr
virtual void ProcessComposites()
Generate composites.
boost::shared_ptr< GeometryVector > Composite
void PrismLines(int prism, PerMap &perFaces, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)
void PrintConfig()
Print out all configuration options for a module.
Represents a command-line configuration option.
std::map< int, std::pair< FaceSharedPtr, std::vector< int > > > PerMap
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.
void ReorderPrisms(PerMap &perFaces)
Reorder node IDs so that prisms and tetrahedra are aligned correctly.
boost::shared_ptr< Element > ElementSharedPtr
std::ofstream m_mshFile
Output stream.
boost::shared_ptr< Face > FaceSharedPtr
Shared pointer to a face.
void OpenStream()
Open a file for output.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
void SetDefaults()
Sets default configuration options for those which have not been set.
ModuleFactory & GetModuleFactory()
const char *const ModuleTypeMap[]