Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 
46 using namespace std;
47 using namespace Nektar;
48 
50 
52 {
57 };
58 
59 
60 struct MoveVerts
61 {
62  vector<NekDouble> kspring;
63  vector<int> springVid;
65 };
66 
67 void GetInterfaceVerts(const int compositeID,
69  vector<int> &InterfaceVerts);
70 
72  SpatialDomains::MeshGraphSharedPtr &mesh, string &fieldfile,
74 
75 void GetNewVertexLocation(TiXmlElement *doc,
77  vector<int> &InterfaceVerts,
78  Array<OneD, NekDouble> &xstreak,
79  Array<OneD, NekDouble> &ystreak,
80  Array<OneD,NekDouble> &vertx,
81  Array<OneD,NekDouble> &verty,
82  int maxiter);
83 
84 void TurnOffEdges(TiXmlElement *doc,
85  SpatialDomains::SegGeomMap &meshedges,
86  Array<OneD,MoveVerts> &verts);
87 
88 void RedefineVertices(TiXmlElement *doc,
90 
92  Array<OneD, NekDouble> &dvertx,
93  Array<OneD, NekDouble> &dverty);
94 
95 int main(int argc, char *argv[])
96 {
97  if(argc !=4)
98  {
99  fprintf(stderr,"Usage: ./MoveMeshToCriticalLayer meshfile streakfile outfile\n");
100  exit(1);
101  }
102 
103  //------------------------------------------------------------
104  // Create Session file by reading only first meshfile
105  //-----------------------------------------------------------
107  = LibUtilities::SessionReader::CreateInstance(2, argv);
108 
109  //-------------------------------------------------------------
110  // Read in mesh from input file
111  //------------------------------------------------------------
112  string meshfile(argv[argc-3]);
113  SpatialDomains::MeshGraphSharedPtr mesh = SpatialDomains::MeshGraph::Read(vSession);
114  TiXmlDocument& meshdoc = vSession->GetDocument();
115  TiXmlHandle docHandle(&meshdoc);
116  TiXmlElement* doc = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
117 
118  //------------------------------------------------------------
119  // Get list of vertices along interface
120  //------------------------------------------------------------
121  vector<int> InterfaceVerts;
122  GetInterfaceVerts(300,mesh,InterfaceVerts);
123 
124  //------------------------------------------------------------
125  // Determine location of new interface vertices
126  //------------------------------------------------------------
127  Array<OneD,NekDouble> xstreak(50),ystreak(50);
128  string fieldfile(argv[argc-2]);
129  GetStreakLocation(vSession,mesh,fieldfile,xstreak,ystreak);
130 
131  //------------------------------------------------------------
132  // Move internal mesh using critical layer info and under string analogy
133  //------------------------------------------------------------
134  Array<OneD,NekDouble> dvertx(mesh->GetNvertices(),0.0), dverty(mesh->GetNvertices(),0.0);
135  int maxiter;
136  vSession->LoadParameter("MoveMeshMaxIterations",maxiter,100);
137 
138  GetNewVertexLocation(doc, mesh,InterfaceVerts,xstreak,ystreak,dvertx,dverty,maxiter);
139 
140  //------------------------------------------------------------
141  // Enforce rotational symmetry on mesh
142  //------------------------------------------------------------
143  if(vSession->DefinesSolverInfo("EnforceRotationalSymmetry"))
144  {
145  EnforceRotationalSymmetry(mesh,dvertx,dverty);
146  }
147 
148  //------------------------------------------------------------
149  // Redfine vertices in doc
150  //------------------------------------------------------------
151  RedefineVertices(doc,dvertx,dverty);
152 
153  //------------------------------------------------------------
154  // Write out moved mesh file
155  //------------------------------------------------------------
156  std::string outfile(argv[argc-1]);
157  meshdoc.SaveFile(outfile);
158 }
159 
160 
161 void GetInterfaceVerts(const int compositeID, SpatialDomains::MeshGraphSharedPtr &mesh, vector<int> &InterfaceVerts)
162 {
163  SpatialDomains::Composite composite;
164  composite = mesh->GetComposite(compositeID);
165  int compsize = composite->size();
166  map<int,int> vertmap;
167 
168  for(int i = 0; i < compsize; ++i)
169  {
170  if(vertmap.count((*composite)[i]->GetVid(0)) == 0)
171  {
172  InterfaceVerts.push_back((*composite)[i]->GetVid(0));
173  vertmap[(*composite)[i]->GetVid(0)] = 1;
174  }
175 
176  if(vertmap.count((*composite)[i]->GetVid(1)) == 0)
177  {
178  InterfaceVerts.push_back((*composite)[i]->GetVid(1));
179  vertmap[(*composite)[i]->GetVid(1)] = 1;
180  }
181  }
182 }
183 
185  SpatialDomains::MeshGraphSharedPtr &mesh, string &fieldfile,
187 {
188  //-------------------------------------------------------------
189  // Define Streak Expansion
191 
193  ::AllocateSharedPtr(vSession,mesh);
194  //---------------------------------------------------------------
195 
196  //----------------------------------------------
197  // Import field file.
198  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
199  vector<vector<NekDouble> > fielddata;
200  LibUtilities::Import(fieldfile,fielddef,fielddata);
201  //----------------------------------------------
202 
203  //----------------------------------------------
204  // Copy data from field file
205  string streak_field("w");
206  for(unsigned int i = 0; i < fielddata.size(); ++i)
207  {
208  streak->ExtractDataToCoeffs(fielddef [i],
209  fielddata[i],
210  streak_field,
211  streak->UpdateCoeffs());
212  }
213  //----------------------------------------------
214 
215  NekDouble cr = 0.0;
216 
217  cerr << "Extracting Critical Layer ";
218  Computestreakpositions(streak,xc, yc,cr);
219 }
220 
221 
222 void GetNewVertexLocation(TiXmlElement *doc,
224  vector<int> &InterfaceVerts,
225  Array<OneD, NekDouble> &xstreak,
226  Array<OneD, NekDouble> &ystreak,
227  Array<OneD,NekDouble> &dvertx,
228  Array<OneD,NekDouble> &dverty,
229  int maxiter)
230 {
231  int i,j,k;
232  int nverts = mesh->GetNvertices();
233 
234  SpatialDomains::SegGeomMap meshedges = mesh->GetAllSegGeoms();
235 
236  Array<OneD,MoveVerts> Verts(nverts);
237 
238  // loop mesh edges and fill in verts info
240 
243 
244  int vid0,vid1;
245  NekDouble kspring;
246  NekDouble x,y,x1,y1,z1,x2,y2,z2;
247 
248  // Setup intiial spring and verts
249  for(segIter = meshedges.begin(); segIter != meshedges.end(); ++segIter)
250  {
251  vid0 = (segIter->second)->GetVid(0);
252  vid1 = (segIter->second)->GetVid(1);
253 
254  v0 = (segIter->second)->GetVertex(0);
255  v1 = (segIter->second)->GetVertex(1);
256 
257  kspring = 1.0/v0->dist(*v1);
258 
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);
263 
264  }
265 
266  // Scale spring by total around vertex and turn on all vertices.
267  for(i = 0; i < nverts; ++i)
268  {
269  NekDouble invktot = 0.0;
270  for(j = 0; j < Verts[i].kspring.size(); ++j)
271  {
272  invktot += Verts[i].kspring[j];
273  }
274  invktot = 1.0/invktot;
275  for(j = 0; j < Verts[i].kspring.size(); ++j)
276  {
277  Verts[i].kspring[j] *= invktot;
278  }
279 
280  Verts[i].solve = eSolveXY;
281  }
282 
283 
284  // Turn off all edges defined by composite lists of correct dimension
285  TurnOffEdges(doc,meshedges,Verts);
286 
287  NekDouble z,h0,h1,h2;
288  // Set interface vertices to lie on critical layer
289  for(i = 0; i < InterfaceVerts.size(); ++i)
290  {
291  Verts[InterfaceVerts[i]].solve = eNoSolve;
292  mesh->GetVertex(InterfaceVerts[i])->GetCoords(x,y,z);
293 
294  for(j = 0; j < xstreak.num_elements()-1; ++j)
295  {
296  if((x >= xstreak[j])&&(x <= xstreak[j+1]))
297  {
298  break;
299  }
300  }
301 
302  ASSERTL0(j != xstreak.num_elements(),"Did not find x location along critical layer");
303 
304  k = (j==0)?1: j; // offset index at beginning
305 
306  // quadraticalling interpolate points
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]));
313 
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;
316  }
317 
318  // shift quads in critical layer to move more or less rigidly
319  SpatialDomains::QuadGeomMap quadgeom = mesh->GetAllQuadGeoms();
321  for(quadIter = quadgeom.begin(); quadIter != quadgeom.end(); ++quadIter)
322  {
323  for(i = 0; i < 4; ++i)
324  {
325  vid0 = (quadIter->second)->GetVid(i);
326 
327  switch(Verts[vid0].solve)
328  {
329  case eSolveXY:
330  {
331  mesh->GetVertex(vid0)->GetCoords(x,y,z);
332 
333  // find nearest interface vert
334  mesh->GetVertex(InterfaceVerts[0])->GetCoords(x1,y1,z1);
335  for(j = 0; j < InterfaceVerts.size()-1; ++j)
336  {
337  mesh->GetVertex(InterfaceVerts[j+1])->GetCoords(x2,y2,z2);
338  if((x >= x1)&&(x < x2))
339  {
340  break;
341  }
342  x1 = x2;
343  y1 = y2;
344  }
345 
346  // currently just shift vert as average of two sides
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]];
351  }
352  break;
353  case eSolveY:
354  {
355  mesh->GetVertex(vid0)->GetCoords(x,y,z);
356  mesh->GetVertex(InterfaceVerts[0])->GetCoords(x1,y1,z1);
357 
358  if(fabs(x-x1) < 1e-6)
359  {
360  dverty[vid0] = dverty[InterfaceVerts[0]];
361  }
362  else
363  {
364  dverty[vid0] = dverty[InterfaceVerts[InterfaceVerts.size()-1]];
365  }
366  }
367  break;
368  default:
369  break;
370  }
371  Verts[vid0].solve = eNoSolve;
372  }
373  }
374 
375 
376 
377 
378 
379  // Iterate internal vertices
380  bool ContinueToIterate = true;
381  int cnt = 0;
382  int nsum = 0;
383  NekDouble dsum,dx,dy,sum,prev_sum = 0.0;
384  NekDouble tol = 1e-3,fac;
385  int blend = 50;
386 
387  while (ContinueToIterate)
388  {
389 
390  sum = 0.0;
391  nsum = 0;
392 
393  // use a ramping function to help move interior slowly
394  fac = (cnt < blend)? 1.0/(blend+1.0)*(cnt+1): 1.0;
395 
396  for(i = 0; i < nverts; ++i)
397  {
398  if(Verts[i].solve != eNoSolve)
399  {
400  dx = dy = 0.0;
401 
402  if((Verts[i].solve == eSolveX)||(Verts[i].solve == eSolveXY))
403  {
404  for(j = 0; j < Verts[i].kspring.size(); ++j)
405  {
406  dx += fac*Verts[i].kspring[j]*dvertx[Verts[i].springVid[j]];
407  }
408  }
409 
410  if((Verts[i].solve == eSolveY)||(Verts[i].solve == eSolveXY))
411  {
412  for(j = 0; j < Verts[i].kspring.size(); ++j)
413  {
414  dy += fac*Verts[i].kspring[j]*dverty[Verts[i].springVid[j]];
415  }
416  }
417 
418  dsum = (dx*dx + dy*dy);
419 
420  dvertx[i] = dx;
421  dverty[i] = dy;
422 
423  if(dsum > 1e-16)
424  {
425  //sum += dsum/(dvertx[i]*dvertx[i] + dverty[i]*dverty[i]);
426  sum += dsum;
427  nsum += 1;
428  }
429  }
430  }
431 
432  if(nsum)
433  {
434  sum = sqrt(sum/(NekDouble)nsum);
435 
436  NekDouble chg = sum-prev_sum;
437  prev_sum = sum;
438 
439  cerr << "Iteration " << cnt << " : " << chg << endl;
440 
441  if((chg < tol)&&(cnt > blend))
442  {
443  ContinueToIterate = false;
444  }
445 
446  }
447  else if(cnt > blend)
448  {
449  ContinueToIterate = false;
450 
451  }
452 
453  if(cnt++ > maxiter)
454  {
455  ContinueToIterate = false;
456  }
457  }
458 }
459 
460 
461 // Read Composites from xml document and turn off verts that are along edge composites.
462 void TurnOffEdges(TiXmlElement *doc,
463  SpatialDomains::SegGeomMap &meshedges,
464  Array<OneD,MoveVerts> &Verts)
465 {
466  TiXmlElement* field = doc->FirstChildElement("COMPOSITE");
467  ASSERTL0(field, "Unable to find COMPOSITE tag in file.");
468 
469  int nextCompositeNumber = -1;
470 
471  /// All elements are of the form: "<C ID = "N"> ... </C>".
472  /// Read the ID field first.
473  TiXmlElement *composite = field->FirstChildElement("C");
474 
475  while (composite)
476  {
477  nextCompositeNumber++;
478 
479  int indx;
480  int err = composite->QueryIntAttribute("ID", &indx);
481  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
482 
483  TiXmlNode* compositeChild = composite->FirstChild();
484  // This is primarily to skip comments that may be present.
485  // Comments appear as nodes just like elements.
486  // We are specifically looking for text in the body
487  // of the definition.
488  while(compositeChild && compositeChild->Type() != TiXmlNode::TINYXML_TEXT)
489  {
490  compositeChild = compositeChild->NextSibling();
491  }
492 
493  ASSERTL0(compositeChild, "Unable to read composite definition body.");
494  std::string compositeStr = compositeChild->ToText()->ValueStr();
495 
496  /// Parse out the element components corresponding to type of element.
497  std::istringstream compositeDataStrm(compositeStr.c_str());
498 
499  try
500  {
501  std::string compositeElementStr;
502  compositeDataStrm >> compositeElementStr;
503 
504  std::istringstream tokenStream(compositeElementStr);
505  char type;
506 
507  tokenStream >> type;
508 
509  // in what follows we are assuming there is only one block of data
510  std::string::size_type indxBeg = compositeElementStr.find_first_of('[') + 1;
511  std::string::size_type indxEnd = compositeElementStr.find_last_of(']') - 1;
512 
513  ASSERTL0(indxBeg <= indxEnd, (std::string("Error reading index definition:") + compositeElementStr).c_str());
514 
515  std::string indxStr = compositeElementStr.substr(indxBeg, indxEnd - indxBeg + 1);
516  std::vector<unsigned int> seqVector;
518 
519  bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
520 
521  ASSERTL0(err, (std::string("Error reading composite elements: ") + indxStr).c_str());
522 
523  switch(type)
524  {
525  case 'E': // Turn off vertices along composite edges
526  {
527  int seqlen = seqVector.size();
528  NekDouble x0,y0,z0,x1,y1,z1;
529  int vid0,vid1;
530 
531  for(int i = 0; i < seqlen; ++i)
532  {
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);
537 
538  if(fabs(x0-x1) < 1e-8)
539  {
540  //identify corners by double visit
541  if(Verts[vid0].solve == eSolveX)
542  {
543  Verts[vid0].solve = eNoSolve;
544  }
545  else
546  {
547  Verts[vid0].solve = eSolveY;
548  }
549 
550  if(Verts[vid1].solve == eSolveX)
551  {
552  Verts[vid1].solve = eNoSolve;
553  }
554  else
555  {
556  Verts[vid1].solve = eSolveY;
557  }
558  }
559 
560  if(fabs(y0-y1) < 1e-8)
561  {
562  //identify corners by double visit
563  if(Verts[vid0].solve == eSolveY)
564  {
565  Verts[vid0].solve = eNoSolve;
566  }
567  else
568  {
569  Verts[vid0].solve = eSolveX;
570  }
571 
572  if(Verts[vid1].solve == eSolveY)
573  {
574  Verts[vid1].solve = eNoSolve;
575  }
576  else
577  {
578  Verts[vid1].solve = eSolveX;
579  }
580  }
581  }
582 
583  }
584  break;
585 
586  case 'T': case 'Q': // do nothing
587  {
588  break;
589  }
590  default:
591  NEKERROR(ErrorUtil::efatal, (std::string("Unrecognized composite token: ") + compositeElementStr).c_str());
592  }
593 
594  }
595  catch(...)
596  {
598  (std::string("Unable to read COMPOSITE data for composite: ") + compositeStr).c_str());
599  }
600 
601  /// Keep looking
602  composite = composite->NextSiblingElement("C");
603  }
604 }
605 
606 void RedefineVertices(TiXmlElement *doc,
607  Array<OneD,NekDouble> &dvertx,
608  Array<OneD,NekDouble> &dverty)
609 {
610 
611 
612  TiXmlElement* element = doc->FirstChildElement("VERTEX");
613  ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
614 
615  TiXmlElement *vertex = element->FirstChildElement("V");
616 
617  int indx;
618  int nextVertexNumber = -1;
619  int err; /// Error value returned by TinyXML.
620 
621  vector<NekDouble> xpts,ypts,zpts;
622  NekDouble xval, yval, zval;
623 
624  while (vertex)
625  {
626  nextVertexNumber++;
627 
628  TiXmlAttribute *vertexAttr = vertex->FirstAttribute();
629  std::string attrName(vertexAttr->Name());
630 
631  ASSERTL0(attrName == "ID", (std::string("Unknown attribute name: ") + attrName).c_str());
632 
633  err = vertexAttr->QueryIntValue(&indx);
634  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
635 
636  // Now read body of vertex
637  std::string vertexBodyStr;
638 
639  TiXmlNode *vertexBody = vertex->FirstChild();
640 
641  while (vertexBody)
642  {
643  // Accumulate all non-comment body data.
644  if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
645  {
646  vertexBodyStr += vertexBody->ToText()->Value();
647  vertexBodyStr += " ";
648  }
649 
650  vertexBody = vertexBody->NextSibling();
651  }
652 
653  ASSERTL0(!vertexBodyStr.empty(), "Vertex definitions must contain vertex data.");
654 
655  // Get vertex data from the data string.
656  std::istringstream vertexDataStrm(vertexBodyStr.c_str());
657 
658  try
659  {
660  while(!vertexDataStrm.fail())
661  {
662  vertexDataStrm >> xval >> yval >> zval;
663  }
664 
665  xval += dvertx[indx];
666  yval += dverty[indx];
667 
668  stringstream s;
669  s << scientific << setprecision(8)
670  << xval << " " << yval << " " << zval;
671 
672  vertex->ReplaceChild(vertex->FirstChild(), TiXmlText(s.str()));
673  }
674  catch(...)
675  {
676  ASSERTL0(false, "Unable to read VERTEX data.");
677  }
678  vertex = vertex->NextSiblingElement("V");
679  }
680 
681 }
682 
684  Array<OneD, NekDouble> &dvertx,
685  Array<OneD, NekDouble> &dverty)
686 {
687  int i,j;
688  int nverts = mesh->GetNvertices();
689  Array<OneD, NekDouble> x(nverts),y(nverts);
690  NekDouble xval,yval,zval;
691 
692  for(i = 0; i < nverts; ++i)
693  {
694  mesh->GetVertex(i)->GetCoords(xval,yval,zval);
695  x[i] = xval + dvertx[i];
696  y[i] = yval + dverty[i];
697  }
698 
699  NekDouble xmax = Vmath::Vmax(nverts,x,1);
700  NekDouble tol = 1e-5, dist2,xrot,yrot;
701  Array<OneD,int> index(nverts);
702  // find nearest
703  for(i = 0; i < nverts; ++i)
704  {
705  xrot = -x[i] + xmax;
706  yrot = -y[i];
707  tol = 1.0;
708 
709  for(j = 0; j < nverts; ++j)
710  {
711  dist2 = (x[j]-xrot)*(x[j]-xrot) + (y[j]-yrot)*(y[j]-yrot);
712  if(dist2 < tol)
713  {
714  index[i] = j;
715  tol = dist2;
716  }
717  }
718  }
719 
720  //average points and recalcualte dvertx, dverty
721  for(i = 0; i < nverts; ++i)
722  {
723  mesh->GetVertex(i)->GetCoords(xval,yval,zval);
724 
725  xrot = 0.5*(-x[index[i]] + xmax + x[i]);
726  yrot = 0.5*(-y[index[i]] + y[i]);
727 
728  dvertx[i] = xrot - xval;
729  dverty[i] = yrot - yval;
730  }
731 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:191
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > &ElementIDs)
This function allows for data to be imported from an FLD file when a session and/or communicator is n...
Definition: FieldIO.cpp:279
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:779
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)
STL namespace.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
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)
void GetInterfaceVerts(const int compositeID, SpatialDomains::MeshGraphSharedPtr &mesh, vector< int > &InterfaceVerts)
std::map< int, SegGeomSharedPtr > SegGeomMap
Definition: SegGeom.h:54
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
double NekDouble
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:114
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:442
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60