69 vector<int> &InterfaceVerts);
77 vector<int> &InterfaceVerts,
95 int main(
int argc,
char *argv[])
99 fprintf(stderr,
"Usage: ./MoveMeshToCriticalLayer meshfile streakfile outfile\n");
107 = LibUtilities::SessionReader::CreateInstance(2, argv);
112 string meshfile(argv[argc-3]);
114 TiXmlDocument& meshdoc = vSession->GetDocument();
115 TiXmlHandle docHandle(&meshdoc);
116 TiXmlElement* doc = docHandle.FirstChildElement(
"NEKTAR").FirstChildElement(
"GEOMETRY").Element();
121 vector<int> InterfaceVerts;
128 string fieldfile(argv[argc-2]);
136 vSession->LoadParameter(
"MoveMeshMaxIterations",maxiter,100);
143 if(vSession->DefinesSolverInfo(
"EnforceRotationalSymmetry"))
156 std::string outfile(argv[argc-1]);
157 meshdoc.SaveFile(outfile);
164 composite = mesh->GetComposite(compositeID);
165 int compsize = composite->size();
166 map<int,int> vertmap;
168 for(
int i = 0; i < compsize; ++i)
170 if(vertmap.count((*composite)[i]->GetVid(0)) == 0)
172 InterfaceVerts.push_back((*composite)[i]->GetVid(0));
173 vertmap[(*composite)[i]->GetVid(0)] = 1;
176 if(vertmap.count((*composite)[i]->GetVid(1)) == 0)
178 InterfaceVerts.push_back((*composite)[i]->GetVid(1));
179 vertmap[(*composite)[i]->GetVid(1)] = 1;
198 vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
199 vector<vector<NekDouble> > fielddata;
205 string streak_field(
"w");
206 for(
unsigned int i = 0; i < fielddata.size(); ++i)
208 streak->ExtractDataToCoeffs(fielddef [i],
211 streak->UpdateCoeffs());
217 cerr <<
"Extracting Critical Layer ";
224 vector<int> &InterfaceVerts,
232 int nverts = mesh->GetNvertices();
249 for(segIter = meshedges.begin(); segIter != meshedges.end(); ++segIter)
251 vid0 = (segIter->second)->GetVid(0);
252 vid1 = (segIter->second)->GetVid(1);
254 v0 = (segIter->second)->GetVertex(0);
255 v1 = (segIter->second)->GetVertex(1);
257 kspring = 1.0/v0->
dist(*v1);
259 Verts[vid0].kspring.push_back(kspring);
260 Verts[vid0].springVid.push_back(vid1);
261 Verts[vid1].kspring.push_back(kspring);
262 Verts[vid1].springVid.push_back(vid0);
267 for(i = 0; i < nverts; ++i)
270 for(j = 0; j < Verts[i].kspring.size(); ++j)
272 invktot += Verts[i].kspring[j];
274 invktot = 1.0/invktot;
275 for(j = 0; j < Verts[i].kspring.size(); ++j)
277 Verts[i].kspring[j] *= invktot;
289 for(i = 0; i < InterfaceVerts.size(); ++i)
291 Verts[InterfaceVerts[i]].solve =
eNoSolve;
292 mesh->GetVertex(InterfaceVerts[i])->GetCoords(x,y,z);
294 for(j = 0; j < xstreak.num_elements()-1; ++j)
296 if((x >= xstreak[j])&&(x <= xstreak[j+1]))
302 ASSERTL0(j != xstreak.num_elements(),
"Did not find x location along critical layer");
307 h0 = (x-xstreak[k])*(x-xstreak[k+1])/
308 ((xstreak[k-1]-xstreak[k])*(xstreak[k-1]-xstreak[k+1]));
309 h1 = (x-xstreak[k-1])*(x-xstreak[k+1])/
310 ((xstreak[k]-xstreak[k-1])*(xstreak[k]-xstreak[k+1]));
311 h2 = (x-xstreak[k-1])*(x-xstreak[k])/
312 ((xstreak[k+1]-xstreak[k-1])*(xstreak[k+1]-xstreak[k]));
314 dvertx[InterfaceVerts[i]] = (xstreak[k-1]*h0 + xstreak[k]*h1 + xstreak[k+1]*h2) - x;
315 dverty[InterfaceVerts[i]] = (ystreak[k-1]*h0 + ystreak[k]*h1 + ystreak[k+1]*h2) - y;
321 for(quadIter = quadgeom.begin(); quadIter != quadgeom.end(); ++quadIter)
323 for(i = 0; i < 4; ++i)
325 vid0 = (quadIter->second)->GetVid(i);
327 switch(Verts[vid0].solve)
331 mesh->GetVertex(vid0)->GetCoords(x,y,z);
334 mesh->GetVertex(InterfaceVerts[0])->GetCoords(x1,y1,z1);
335 for(j = 0; j < InterfaceVerts.size()-1; ++j)
337 mesh->GetVertex(InterfaceVerts[j+1])->GetCoords(x2,y2,z2);
338 if((x >= x1)&&(x < x2))
347 dvertx[vid0] = (x2-x)/(x2-x1)*dvertx[InterfaceVerts[j]]+
348 (x-x1)/(x2-x1)*dvertx[InterfaceVerts[j+1]];
349 dverty[vid0] = (x2-x)/(x2-x1)*dverty[InterfaceVerts[j]]+
350 (x-x1)/(x2-x1)*dverty[InterfaceVerts[j+1]];
355 mesh->GetVertex(vid0)->GetCoords(x,y,z);
356 mesh->GetVertex(InterfaceVerts[0])->GetCoords(x1,y1,z1);
358 if(fabs(x-x1) < 1e-6)
360 dverty[vid0] = dverty[InterfaceVerts[0]];
364 dverty[vid0] = dverty[InterfaceVerts[InterfaceVerts.size()-1]];
380 bool ContinueToIterate =
true;
387 while (ContinueToIterate)
394 fac = (cnt < blend)? 1.0/(blend+1.0)*(cnt+1): 1.0;
396 for(i = 0; i < nverts; ++i)
404 for(j = 0; j < Verts[i].kspring.size(); ++j)
406 dx += fac*Verts[i].kspring[j]*dvertx[Verts[i].springVid[j]];
412 for(j = 0; j < Verts[i].kspring.size(); ++j)
414 dy += fac*Verts[i].kspring[j]*dverty[Verts[i].springVid[j]];
418 dsum = (dx*dx + dy*dy);
439 cerr <<
"Iteration " << cnt <<
" : " << chg << endl;
441 if((chg < tol)&&(cnt > blend))
443 ContinueToIterate =
false;
449 ContinueToIterate =
false;
455 ContinueToIterate =
false;
466 TiXmlElement* field = doc->FirstChildElement(
"COMPOSITE");
467 ASSERTL0(field,
"Unable to find COMPOSITE tag in file.");
469 int nextCompositeNumber = -1;
473 TiXmlElement *composite = field->FirstChildElement(
"C");
477 nextCompositeNumber++;
480 int err = composite->QueryIntAttribute(
"ID", &indx);
481 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read attribute ID.");
483 TiXmlNode* compositeChild = composite->FirstChild();
488 while(compositeChild && compositeChild->Type() != TiXmlNode::TINYXML_TEXT)
490 compositeChild = compositeChild->NextSibling();
493 ASSERTL0(compositeChild,
"Unable to read composite definition body.");
494 std::string compositeStr = compositeChild->ToText()->ValueStr();
497 std::istringstream compositeDataStrm(compositeStr.c_str());
501 std::string compositeElementStr;
502 compositeDataStrm >> compositeElementStr;
504 std::istringstream tokenStream(compositeElementStr);
510 std::string::size_type indxBeg = compositeElementStr.find_first_of(
'[') + 1;
511 std::string::size_type indxEnd = compositeElementStr.find_last_of(
']') - 1;
513 ASSERTL0(indxBeg <= indxEnd, (std::string(
"Error reading index definition:") + compositeElementStr).c_str());
515 std::string indxStr = compositeElementStr.substr(indxBeg, indxEnd - indxBeg + 1);
516 std::vector<unsigned int> seqVector;
519 bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
521 ASSERTL0(err, (std::string(
"Error reading composite elements: ") + indxStr).c_str());
527 int seqlen = seqVector.size();
531 for(
int i = 0; i < seqlen; ++i)
533 meshedges[seqVector[i]]->GetVertex(0)->GetCoords(x0,y0,z0);
534 meshedges[seqVector[i]]->GetVertex(1)->GetCoords(x1,y1,z1);
535 vid0 = meshedges[seqVector[i]]->GetVid(0);
536 vid1 = meshedges[seqVector[i]]->GetVid(1);
538 if(fabs(x0-x1) < 1e-8)
541 if(Verts[vid0].solve ==
eSolveX)
550 if(Verts[vid1].solve ==
eSolveX)
560 if(fabs(y0-y1) < 1e-8)
563 if(Verts[vid0].solve ==
eSolveY)
572 if(Verts[vid1].solve ==
eSolveY)
598 (std::string(
"Unable to read COMPOSITE data for composite: ") + compositeStr).c_str());
602 composite = composite->NextSiblingElement(
"C");
612 TiXmlElement* element = doc->FirstChildElement(
"VERTEX");
613 ASSERTL0(element,
"Unable to find mesh VERTEX tag in file.");
615 TiXmlElement *vertex = element->FirstChildElement(
"V");
618 int nextVertexNumber = -1;
621 vector<NekDouble> xpts,ypts,zpts;
628 TiXmlAttribute *vertexAttr = vertex->FirstAttribute();
629 std::string attrName(vertexAttr->Name());
631 ASSERTL0(attrName ==
"ID", (std::string(
"Unknown attribute name: ") + attrName).c_str());
633 err = vertexAttr->QueryIntValue(&indx);
634 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read attribute ID.");
637 std::string vertexBodyStr;
639 TiXmlNode *vertexBody = vertex->FirstChild();
644 if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
646 vertexBodyStr += vertexBody->ToText()->Value();
647 vertexBodyStr +=
" ";
650 vertexBody = vertexBody->NextSibling();
653 ASSERTL0(!vertexBodyStr.empty(),
"Vertex definitions must contain vertex data.");
656 std::istringstream vertexDataStrm(vertexBodyStr.c_str());
660 while(!vertexDataStrm.fail())
662 vertexDataStrm >> xval >> yval >> zval;
665 xval += dvertx[indx];
666 yval += dverty[indx];
669 s << scientific << setprecision(8)
670 << xval <<
" " << yval <<
" " << zval;
672 vertex->ReplaceChild(vertex->FirstChild(), TiXmlText(s.str()));
676 ASSERTL0(
false,
"Unable to read VERTEX data.");
678 vertex = vertex->NextSiblingElement(
"V");
688 int nverts = mesh->GetNvertices();
692 for(i = 0; i < nverts; ++i)
694 mesh->GetVertex(i)->GetCoords(xval,yval,zval);
695 x[i] = xval + dvertx[i];
696 y[i] = yval + dverty[i];
703 for(i = 0; i < nverts; ++i)
709 for(j = 0; j < nverts; ++j)
711 dist2 = (x[j]-xrot)*(x[j]-xrot) + (y[j]-yrot)*(y[j]-yrot);
721 for(i = 0; i < nverts; ++i)
723 mesh->GetVertex(i)->GetCoords(xval,yval,zval);
725 xrot = 0.5*(-x[index[i]] + xmax + x[i]);
726 yrot = 0.5*(-y[index[i]] + y[i]);
728 dvertx[i] = xrot - xval;
729 dverty[i] = yrot - yval;
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
vector< NekDouble > kspring
void GetStreakLocation(LibUtilities::SessionReaderSharedPtr &vSession, SpatialDomains::MeshGraphSharedPtr &mesh, string &fieldfile, Array< OneD, NekDouble > &xc, Array< OneD, NekDouble > &yc)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > ElementiDs)
Imports an FLD file.
void EnforceRotationalSymmetry(SpatialDomains::MeshGraphSharedPtr &mesh, Array< OneD, NekDouble > &dvertx, Array< OneD, NekDouble > &dverty)
void GetNewVertexLocation(TiXmlElement *doc, SpatialDomains::MeshGraphSharedPtr &mesh, vector< int > &InterfaceVerts, Array< OneD, NekDouble > &xstreak, Array< OneD, NekDouble > &ystreak, Array< OneD, NekDouble > &vertx, Array< OneD, NekDouble > &verty, int maxiter)
void TurnOffEdges(TiXmlElement *doc, SpatialDomains::SegGeomMap &meshedges, Array< OneD, MoveVerts > &verts)
void GetInterfaceVerts(const int compositeID, SpatialDomains::MeshGraphSharedPtr &mesh, vector< int > &InterfaceVerts)
std::map< int, SegGeomSharedPtr > SegGeomMap
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
boost::shared_ptr< GeometryVector > Composite
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::map< int, QuadGeomSharedPtr > QuadGeomMap
void RedefineVertices(TiXmlElement *doc, Array< OneD, NekDouble > &dvertx, Array< OneD, NekDouble > &dverty)
NekDouble dist(PointGeom &a)
int main(int argc, char *argv[])
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
boost::shared_ptr< PointGeom > PointGeomSharedPtr