Nektar++
MoveMeshToCriticalLayer.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File MoveMeshToCriticalLayer.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: MoveMesh to critical layer
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <cstdio>
37 #include <cstdlib>
38 #include <iomanip>
39 
40 #include <tinyxml.h>
42 #include <MultiRegions/ExpList.h>
43 #include <MultiRegions/ExpList2D.h>
45 
47 
49 {
54 };
55 
56 
57 struct MoveVerts
58 {
59  vector<NekDouble> kspring;
60  vector<int> springVid;
62 };
63 
64 void GetInterfaceVerts(const int compositeID,
66  vector<int> &InterfaceVerts);
67 
69  SpatialDomains::MeshGraphSharedPtr &mesh, string &fieldfile,
71 
72 void GetNewVertexLocation(TiXmlElement *doc,
74  vector<int> &InterfaceVerts,
75  Array<OneD, NekDouble> &xstreak,
76  Array<OneD, NekDouble> &ystreak,
77  Array<OneD,NekDouble> &vertx,
78  Array<OneD,NekDouble> &verty,
79  int maxiter);
80 
81 void TurnOffEdges(TiXmlElement *doc,
82  SpatialDomains::SegGeomMap &meshedges,
83  Array<OneD,MoveVerts> &verts);
84 
85 void RedefineVertices(TiXmlElement *doc,
87 
89  Array<OneD, NekDouble> &dvertx,
90  Array<OneD, NekDouble> &dverty);
91 
92 int main(int argc, char *argv[])
93 {
94  if(argc !=4)
95  {
96  fprintf(stderr,"Usage: ./MoveMeshToCriticalLayer meshfile streakfile outfile\n");
97  exit(1);
98  }
99 
100  //------------------------------------------------------------
101  // Create Session file by reading only first meshfile
102  //-----------------------------------------------------------
104  = LibUtilities::SessionReader::CreateInstance(2, argv);
105 
106  //-------------------------------------------------------------
107  // Read in mesh from input file
108  //------------------------------------------------------------
109  string meshfile(argv[argc-3]);
110  SpatialDomains::MeshGraphSharedPtr mesh = SpatialDomains::MeshGraph::Read(vSession);
111  TiXmlDocument& meshdoc = vSession->GetDocument();
112  TiXmlHandle docHandle(&meshdoc);
113  TiXmlElement* doc = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
114 
115  //------------------------------------------------------------
116  // Get list of vertices along interface
117  //------------------------------------------------------------
118  vector<int> InterfaceVerts;
119  GetInterfaceVerts(300,mesh,InterfaceVerts);
120 
121  //------------------------------------------------------------
122  // Determine location of new interface vertices
123  //------------------------------------------------------------
124  Array<OneD,NekDouble> xstreak(50),ystreak(50);
125  string fieldfile(argv[argc-2]);
126  GetStreakLocation(vSession,mesh,fieldfile,xstreak,ystreak);
127 
128  //------------------------------------------------------------
129  // Move internal mesh using critical layer info and under string analogy
130  //------------------------------------------------------------
131  Array<OneD,NekDouble> dvertx(mesh->GetNvertices(),0.0), dverty(mesh->GetNvertices(),0.0);
132  int maxiter;
133  vSession->LoadParameter("MoveMeshMaxIterations",maxiter,100);
134 
135  GetNewVertexLocation(doc, mesh,InterfaceVerts,xstreak,ystreak,dvertx,dverty,maxiter);
136 
137  //------------------------------------------------------------
138  // Enforce rotational symmetry on mesh
139  //------------------------------------------------------------
140  if(vSession->DefinesSolverInfo("EnforceRotationalSymmetry"))
141  {
142  EnforceRotationalSymmetry(mesh,dvertx,dverty);
143  }
144 
145  //------------------------------------------------------------
146  // Redfine vertices in doc
147  //------------------------------------------------------------
148  RedefineVertices(doc,dvertx,dverty);
149 
150  //------------------------------------------------------------
151  // Write out moved mesh file
152  //------------------------------------------------------------
153  std::string outfile(argv[argc-1]);
154  meshdoc.SaveFile(outfile);
155 }
156 
157 
158 void GetInterfaceVerts(const int compositeID, SpatialDomains::MeshGraphSharedPtr &mesh, vector<int> &InterfaceVerts)
159 {
160  SpatialDomains::Composite composite;
161  composite = mesh->GetComposite(compositeID);
162  int compsize = composite->size();
163  map<int,int> vertmap;
164 
165  for(int i = 0; i < compsize; ++i)
166  {
167  if(vertmap.count((*composite)[i]->GetVid(0)) == 0)
168  {
169  InterfaceVerts.push_back((*composite)[i]->GetVid(0));
170  vertmap[(*composite)[i]->GetVid(0)] = 1;
171  }
172 
173  if(vertmap.count((*composite)[i]->GetVid(1)) == 0)
174  {
175  InterfaceVerts.push_back((*composite)[i]->GetVid(1));
176  vertmap[(*composite)[i]->GetVid(1)] = 1;
177  }
178  }
179 }
180 
182  SpatialDomains::MeshGraphSharedPtr &mesh, string &fieldfile,
184 {
185  //-------------------------------------------------------------
186  // Define Streak Expansion
188 
190  ::AllocateSharedPtr(vSession,mesh);
191  //---------------------------------------------------------------
192 
193  //----------------------------------------------
194  // Import field file.
195  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
196  vector<vector<NekDouble> > fielddata;
197  LibUtilities::Import(fieldfile,fielddef,fielddata);
198  //----------------------------------------------
199 
200  //----------------------------------------------
201  // Copy data from field file
202  string streak_field("w");
203  for(unsigned int i = 0; i < fielddata.size(); ++i)
204  {
205  streak->ExtractDataToCoeffs(fielddef [i],
206  fielddata[i],
207  streak_field,
208  streak->UpdateCoeffs());
209  }
210  //----------------------------------------------
211 
212  NekDouble cr = 0.0;
213 
214  cerr << "Extracting Critical Layer ";
215  Computestreakpositions(streak,xc, yc,cr);
216 }
217 
218 
219 void GetNewVertexLocation(TiXmlElement *doc,
221  vector<int> &InterfaceVerts,
222  Array<OneD, NekDouble> &xstreak,
223  Array<OneD, NekDouble> &ystreak,
224  Array<OneD,NekDouble> &dvertx,
225  Array<OneD,NekDouble> &dverty,
226  int maxiter)
227 {
228  int i,j,k;
229  int nverts = mesh->GetNvertices();
230 
231  SpatialDomains::SegGeomMap meshedges = mesh->GetAllSegGeoms();
232 
233  Array<OneD,MoveVerts> Verts(nverts);
234 
235  // loop mesh edges and fill in verts info
237 
240 
241  int vid0,vid1;
242  NekDouble kspring;
243  NekDouble x,y,x1,y1,z1,x2,y2,z2;
244 
245  // Setup intiial spring and verts
246  for(segIter = meshedges.begin(); segIter != meshedges.end(); ++segIter)
247  {
248  vid0 = (segIter->second)->GetVid(0);
249  vid1 = (segIter->second)->GetVid(1);
250 
251  v0 = (segIter->second)->GetVertex(0);
252  v1 = (segIter->second)->GetVertex(1);
253 
254  kspring = 1.0/v0->dist(*v1);
255 
256  Verts[vid0].kspring.push_back(kspring);
257  Verts[vid0].springVid.push_back(vid1);
258  Verts[vid1].kspring.push_back(kspring);
259  Verts[vid1].springVid.push_back(vid0);
260 
261  }
262 
263  // Scale spring by total around vertex and turn on all vertices.
264  for(i = 0; i < nverts; ++i)
265  {
266  NekDouble invktot = 0.0;
267  for(j = 0; j < Verts[i].kspring.size(); ++j)
268  {
269  invktot += Verts[i].kspring[j];
270  }
271  invktot = 1.0/invktot;
272  for(j = 0; j < Verts[i].kspring.size(); ++j)
273  {
274  Verts[i].kspring[j] *= invktot;
275  }
276 
277  Verts[i].solve = eSolveXY;
278  }
279 
280 
281  // Turn off all edges defined by composite lists of correct dimension
282  TurnOffEdges(doc,meshedges,Verts);
283 
284  NekDouble z,h0,h1,h2;
285  // Set interface vertices to lie on critical layer
286  for(i = 0; i < InterfaceVerts.size(); ++i)
287  {
288  Verts[InterfaceVerts[i]].solve = eNoSolve;
289  mesh->GetVertex(InterfaceVerts[i])->GetCoords(x,y,z);
290 
291  for(j = 0; j < xstreak.num_elements()-1; ++j)
292  {
293  if((x >= xstreak[j])&&(x <= xstreak[j+1]))
294  {
295  break;
296  }
297  }
298 
299  ASSERTL0(j != xstreak.num_elements(),"Did not find x location along critical layer");
300 
301  k = (j==0)?1: j; // offset index at beginning
302 
303  // quadraticalling interpolate points
304  h0 = (x-xstreak[k])*(x-xstreak[k+1])/
305  ((xstreak[k-1]-xstreak[k])*(xstreak[k-1]-xstreak[k+1]));
306  h1 = (x-xstreak[k-1])*(x-xstreak[k+1])/
307  ((xstreak[k]-xstreak[k-1])*(xstreak[k]-xstreak[k+1]));
308  h2 = (x-xstreak[k-1])*(x-xstreak[k])/
309  ((xstreak[k+1]-xstreak[k-1])*(xstreak[k+1]-xstreak[k]));
310 
311  dvertx[InterfaceVerts[i]] = (xstreak[k-1]*h0 + xstreak[k]*h1 + xstreak[k+1]*h2) - x;
312  dverty[InterfaceVerts[i]] = (ystreak[k-1]*h0 + ystreak[k]*h1 + ystreak[k+1]*h2) - y;
313  }
314 
315  // shift quads in critical layer to move more or less rigidly
316  SpatialDomains::QuadGeomMap quadgeom = mesh->GetAllQuadGeoms();
318  for(quadIter = quadgeom.begin(); quadIter != quadgeom.end(); ++quadIter)
319  {
320  for(i = 0; i < 4; ++i)
321  {
322  vid0 = (quadIter->second)->GetVid(i);
323 
324  switch(Verts[vid0].solve)
325  {
326  case eSolveXY:
327  {
328  mesh->GetVertex(vid0)->GetCoords(x,y,z);
329 
330  // find nearest interface vert
331  mesh->GetVertex(InterfaceVerts[0])->GetCoords(x1,y1,z1);
332  for(j = 0; j < InterfaceVerts.size()-1; ++j)
333  {
334  mesh->GetVertex(InterfaceVerts[j+1])->GetCoords(x2,y2,z2);
335  if((x >= x1)&&(x < x2))
336  {
337  break;
338  }
339  x1 = x2;
340  y1 = y2;
341  }
342 
343  // currently just shift vert as average of two sides
344  dvertx[vid0] = (x2-x)/(x2-x1)*dvertx[InterfaceVerts[j]]+
345  (x-x1)/(x2-x1)*dvertx[InterfaceVerts[j+1]];
346  dverty[vid0] = (x2-x)/(x2-x1)*dverty[InterfaceVerts[j]]+
347  (x-x1)/(x2-x1)*dverty[InterfaceVerts[j+1]];
348  }
349  break;
350  case eSolveY:
351  {
352  mesh->GetVertex(vid0)->GetCoords(x,y,z);
353  mesh->GetVertex(InterfaceVerts[0])->GetCoords(x1,y1,z1);
354 
355  if(fabs(x-x1) < 1e-6)
356  {
357  dverty[vid0] = dverty[InterfaceVerts[0]];
358  }
359  else
360  {
361  dverty[vid0] = dverty[InterfaceVerts[InterfaceVerts.size()-1]];
362  }
363  }
364  break;
365  default:
366  break;
367  }
368  Verts[vid0].solve = eNoSolve;
369  }
370  }
371 
372 
373 
374 
375 
376  // Iterate internal vertices
377  bool ContinueToIterate = true;
378  int cnt = 0;
379  int nsum = 0;
380  NekDouble dsum,dx,dy,sum,prev_sum = 0.0;
381  NekDouble tol = 1e-3,fac;
382  int blend = 50;
383 
384  while (ContinueToIterate)
385  {
386 
387  sum = 0.0;
388  nsum = 0;
389 
390  // use a ramping function to help move interior slowly
391  fac = (cnt < blend)? 1.0/(blend+1.0)*(cnt+1): 1.0;
392 
393  for(i = 0; i < nverts; ++i)
394  {
395  if(Verts[i].solve != eNoSolve)
396  {
397  dx = dy = 0.0;
398 
399  if((Verts[i].solve == eSolveX)||(Verts[i].solve == eSolveXY))
400  {
401  for(j = 0; j < Verts[i].kspring.size(); ++j)
402  {
403  dx += fac*Verts[i].kspring[j]*dvertx[Verts[i].springVid[j]];
404  }
405  }
406 
407  if((Verts[i].solve == eSolveY)||(Verts[i].solve == eSolveXY))
408  {
409  for(j = 0; j < Verts[i].kspring.size(); ++j)
410  {
411  dy += fac*Verts[i].kspring[j]*dverty[Verts[i].springVid[j]];
412  }
413  }
414 
415  dsum = (dx*dx + dy*dy);
416 
417  dvertx[i] = dx;
418  dverty[i] = dy;
419 
420  if(dsum > 1e-16)
421  {
422  //sum += dsum/(dvertx[i]*dvertx[i] + dverty[i]*dverty[i]);
423  sum += dsum;
424  nsum += 1;
425  }
426  }
427  }
428 
429  if(nsum)
430  {
431  sum = sqrt(sum/(NekDouble)nsum);
432 
433  NekDouble chg = sum-prev_sum;
434  prev_sum = sum;
435 
436  cerr << "Iteration " << cnt << " : " << chg << endl;
437 
438  if((chg < tol)&&(cnt > blend))
439  {
440  ContinueToIterate = false;
441  }
442 
443  }
444  else if(cnt > blend)
445  {
446  ContinueToIterate = false;
447 
448  }
449 
450  if(cnt++ > maxiter)
451  {
452  ContinueToIterate = false;
453  }
454  }
455 }
456 
457 
458 // Read Composites from xml document and turn off verts that are along edge composites.
459 void TurnOffEdges(TiXmlElement *doc,
460  SpatialDomains::SegGeomMap &meshedges,
461  Array<OneD,MoveVerts> &Verts)
462 {
463  TiXmlElement* field = doc->FirstChildElement("COMPOSITE");
464  ASSERTL0(field, "Unable to find COMPOSITE tag in file.");
465 
466  int nextCompositeNumber = -1;
467 
468  /// All elements are of the form: "<C ID = "N"> ... </C>".
469  /// Read the ID field first.
470  TiXmlElement *composite = field->FirstChildElement("C");
471 
472  while (composite)
473  {
474  nextCompositeNumber++;
475 
476  int indx;
477  int err = composite->QueryIntAttribute("ID", &indx);
478  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
479 
480  TiXmlNode* compositeChild = composite->FirstChild();
481  // This is primarily to skip comments that may be present.
482  // Comments appear as nodes just like elements.
483  // We are specifically looking for text in the body
484  // of the definition.
485  while(compositeChild && compositeChild->Type() != TiXmlNode::TINYXML_TEXT)
486  {
487  compositeChild = compositeChild->NextSibling();
488  }
489 
490  ASSERTL0(compositeChild, "Unable to read composite definition body.");
491  std::string compositeStr = compositeChild->ToText()->ValueStr();
492 
493  /// Parse out the element components corresponding to type of element.
494  std::istringstream compositeDataStrm(compositeStr.c_str());
495 
496  try
497  {
498  std::string compositeElementStr;
499  compositeDataStrm >> compositeElementStr;
500 
501  std::istringstream tokenStream(compositeElementStr);
502  char type;
503 
504  tokenStream >> type;
505 
506  // in what follows we are assuming there is only one block of data
507  std::string::size_type indxBeg = compositeElementStr.find_first_of('[') + 1;
508  std::string::size_type indxEnd = compositeElementStr.find_last_of(']') - 1;
509 
510  ASSERTL0(indxBeg <= indxEnd, (std::string("Error reading index definition:") + compositeElementStr).c_str());
511 
512  std::string indxStr = compositeElementStr.substr(indxBeg, indxEnd - indxBeg + 1);
513  std::vector<unsigned int> seqVector;
515 
516  bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
517 
518  ASSERTL0(err, (std::string("Error reading composite elements: ") + indxStr).c_str());
519 
520  switch(type)
521  {
522  case 'E': // Turn off vertices along composite edges
523  {
524  int seqlen = seqVector.size();
525  NekDouble x0,y0,z0,x1,y1,z1;
526  int vid0,vid1;
527 
528  for(int i = 0; i < seqlen; ++i)
529  {
530  meshedges[seqVector[i]]->GetVertex(0)->GetCoords(x0,y0,z0);
531  meshedges[seqVector[i]]->GetVertex(1)->GetCoords(x1,y1,z1);
532  vid0 = meshedges[seqVector[i]]->GetVid(0);
533  vid1 = meshedges[seqVector[i]]->GetVid(1);
534 
535  if(fabs(x0-x1) < 1e-8)
536  {
537  //identify corners by double visit
538  if(Verts[vid0].solve == eSolveX)
539  {
540  Verts[vid0].solve = eNoSolve;
541  }
542  else
543  {
544  Verts[vid0].solve = eSolveY;
545  }
546 
547  if(Verts[vid1].solve == eSolveX)
548  {
549  Verts[vid1].solve = eNoSolve;
550  }
551  else
552  {
553  Verts[vid1].solve = eSolveY;
554  }
555  }
556 
557  if(fabs(y0-y1) < 1e-8)
558  {
559  //identify corners by double visit
560  if(Verts[vid0].solve == eSolveY)
561  {
562  Verts[vid0].solve = eNoSolve;
563  }
564  else
565  {
566  Verts[vid0].solve = eSolveX;
567  }
568 
569  if(Verts[vid1].solve == eSolveY)
570  {
571  Verts[vid1].solve = eNoSolve;
572  }
573  else
574  {
575  Verts[vid1].solve = eSolveX;
576  }
577  }
578  }
579 
580  }
581  break;
582 
583  case 'T': case 'Q': // do nothing
584  {
585  break;
586  }
587  default:
588  NEKERROR(ErrorUtil::efatal, (std::string("Unrecognized composite token: ") + compositeElementStr).c_str());
589  }
590 
591  }
592  catch(...)
593  {
595  (std::string("Unable to read COMPOSITE data for composite: ") + compositeStr).c_str());
596  }
597 
598  /// Keep looking
599  composite = composite->NextSiblingElement("C");
600  }
601 }
602 
603 void RedefineVertices(TiXmlElement *doc,
604  Array<OneD,NekDouble> &dvertx,
605  Array<OneD,NekDouble> &dverty)
606 {
607 
608 
609  TiXmlElement* element = doc->FirstChildElement("VERTEX");
610  ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
611 
612  TiXmlElement *vertex = element->FirstChildElement("V");
613 
614  int indx;
615  int nextVertexNumber = -1;
616  int err; /// Error value returned by TinyXML.
617 
618  vector<NekDouble> xpts,ypts,zpts;
619  NekDouble xval, yval, zval;
620 
621  while (vertex)
622  {
623  nextVertexNumber++;
624 
625  TiXmlAttribute *vertexAttr = vertex->FirstAttribute();
626  std::string attrName(vertexAttr->Name());
627 
628  ASSERTL0(attrName == "ID", (std::string("Unknown attribute name: ") + attrName).c_str());
629 
630  err = vertexAttr->QueryIntValue(&indx);
631  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
632 
633  // Now read body of vertex
634  std::string vertexBodyStr;
635 
636  TiXmlNode *vertexBody = vertex->FirstChild();
637 
638  while (vertexBody)
639  {
640  // Accumulate all non-comment body data.
641  if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
642  {
643  vertexBodyStr += vertexBody->ToText()->Value();
644  vertexBodyStr += " ";
645  }
646 
647  vertexBody = vertexBody->NextSibling();
648  }
649 
650  ASSERTL0(!vertexBodyStr.empty(), "Vertex definitions must contain vertex data.");
651 
652  // Get vertex data from the data string.
653  std::istringstream vertexDataStrm(vertexBodyStr.c_str());
654 
655  try
656  {
657  while(!vertexDataStrm.fail())
658  {
659  vertexDataStrm >> xval >> yval >> zval;
660  }
661 
662  xval += dvertx[indx];
663  yval += dverty[indx];
664 
665  stringstream s;
666  s << scientific << setprecision(8)
667  << xval << " " << yval << " " << zval;
668 
669  vertex->ReplaceChild(vertex->FirstChild(), TiXmlText(s.str()));
670  }
671  catch(...)
672  {
673  ASSERTL0(false, "Unable to read VERTEX data.");
674  }
675  vertex = vertex->NextSiblingElement("V");
676  }
677 
678 }
679 
681  Array<OneD, NekDouble> &dvertx,
682  Array<OneD, NekDouble> &dverty)
683 {
684  int i,j;
685  int nverts = mesh->GetNvertices();
686  Array<OneD, NekDouble> x(nverts),y(nverts);
687  NekDouble xval,yval,zval;
688 
689  for(i = 0; i < nverts; ++i)
690  {
691  mesh->GetVertex(i)->GetCoords(xval,yval,zval);
692  x[i] = xval + dvertx[i];
693  y[i] = yval + dverty[i];
694  }
695 
696  NekDouble xmax = Vmath::Vmax(nverts,x,1);
697  NekDouble tol = 1e-5, dist2,xrot,yrot;
698  Array<OneD,int> index(nverts);
699  // find nearest
700  for(i = 0; i < nverts; ++i)
701  {
702  xrot = -x[i] + xmax;
703  yrot = -y[i];
704  tol = 1.0;
705 
706  for(j = 0; j < nverts; ++j)
707  {
708  dist2 = (x[j]-xrot)*(x[j]-xrot) + (y[j]-yrot)*(y[j]-yrot);
709  if(dist2 < tol)
710  {
711  index[i] = j;
712  tol = dist2;
713  }
714  }
715  }
716 
717  //average points and recalcualte dvertx, dverty
718  for(i = 0; i < nverts; ++i)
719  {
720  mesh->GetVertex(i)->GetCoords(xval,yval,zval);
721 
722  xrot = 0.5*(-x[index[i]] + xmax + x[i]);
723  yrot = 0.5*(-y[index[i]] + y[i]);
724 
725  dvertx[i] = xrot - xval;
726  dverty[i] = yrot - yval;
727  }
728 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:158
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:756
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
Definition: MeshPartition.h:50
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.
Definition: FieldIO.cpp:106
void Computestreakpositions(MultiRegions::ExpListSharedPtr &streak, Array< OneD, NekDouble > &xc, Array< OneD, NekDouble > &yc, NekDouble cr, NekDouble trans)
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)
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: ExpList.h:1340
void GetInterfaceVerts(const int compositeID, SpatialDomains::MeshGraphSharedPtr &mesh, vector< int > &InterfaceVerts)
std::map< int, SegGeomSharedPtr > SegGeomMap
Definition: SegGeom.h:54
double NekDouble
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:111
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
vector< int > springVid
std::map< int, QuadGeomSharedPtr > QuadGeomMap
Definition: QuadGeom.h:57
void RedefineVertices(TiXmlElement *doc, Array< OneD, NekDouble > &dvertx, Array< OneD, NekDouble > &dverty)
NekDouble dist(PointGeom &a)
Definition: PointGeom.cpp:186
int main(int argc, char *argv[])
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60