179 if(argc > 6 || argc < 5)
182 "Usage: ./MoveMesh meshfile fieldfile changefile alpha cr(optional)\n");
188 = LibUtilities::SessionReader::CreateInstance(2, argv);
192 vSession->DefinesSolverInfo(
"INTERFACE")
193 && vSession->GetSolverInfo(
"INTERFACE")==
"phase" )
195 cr = boost::lexical_cast<
NekDouble>(argv[argc-1]);
200 string meshfile(argv[argc-4]);
218 string changefile(argv[argc-2]);
222 string charalp (argv[argc-1]);
224 cout<<
"read alpha="<<charalp<<endl;
228 string fieldfile(argv[argc-3]);
229 vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
230 vector<vector<NekDouble> > fielddata;
239 for(
int i=0; i<fielddata.size(); i++)
241 streak->ExtractDataToCoeffs(fielddef[i], fielddata[i], fielddef[i]->m_fields[0], streak->UpdateCoeffs());
243 streak->BwdTrans_IterPerExp(streak->GetCoeffs(), streak->UpdatePhys());
249 int nIregions, lastIregion;
254 int nbnd= bndConditions.num_elements();
255 for(
int r=0; r<nbnd; r++)
257 if(bndConditions[r]->GetUserDefined()==
"CalcBC")
265 ASSERTL0(nIregions>0,
"there is any boundary region with the tag USERDEFINEDTYPE=""CalcBC"" specified");
266 cout<<
"nIregions="<<nIregions<<endl;
271 int nedges = bndfieldx[lastIregion]->GetExpSize();
272 int nvertl = nedges +1 ;
288 ( (bndfieldx[lastIregion]->GetExp(0)->as<LocalRegions::SegExp>())
293 vertex0->GetCoords(x0,y0,z0);
296 cout<<
"WARNING x0="<<x0<<endl;
302 Vids_low, v1, v2 , x_connect ,lastedge, xold_low,yold_low);
303 ASSERTL0(Vids_low[v2]!=-10,
"Vids_low[v2] is wrong");
307 cout<<
"x_conn="<<x_connect<<
" yt="<<yt<<
" zt="<<zt<<
" vid="<<Vids_low[v2]<<endl;
308 vertex->GetCoords(x_connect,yt,zt);
315 Vids_low, v1, v2 , x_connect, lastedge, xold_low, yold_low );
318 vertex->GetCoords(x_connect,yt,zt);
332 ( (bndfieldx[lastIregion]->GetExp(0)->as<LocalRegions::SegExp>())
337 vertex0->GetCoords(x0,y0,z0);
340 cout<<
"WARNING x0="<<x0<<endl;
348 Vids_up, v1, v2 , x_connect ,lastedge, xold_up, yold_up);
350 vertexU->GetCoords(x_connect,yt,zt);
357 Vids_up, v1, v2 , x_connect, lastedge, xold_up, yold_up );
361 vertex->GetCoords(x_connect,yt,zt);
373 graphShPt->GetVertex(((bndfieldx[lastIregion]->GetExp(0)
374 ->as<LocalRegions::SegExp>())->GetGeom1D())->GetVid(0));
375 vertex0->GetCoords(x0,y0,z0);
378 cout<<
"WARNING x0="<<x0<<endl;
387 Vids_c, v1, v2 , x_connect ,lastedge, xold_c, yold_c);
391 vertexc->GetCoords(x_connect,yt,zt);
398 Vids_c, v1, v2 , x_connect, lastedge, xold_c, yold_c );
402 vertex->GetCoords(x_connect,yt,zt);
411 for(
int r=0; r<nvertl; r++)
415 Deltaup[r] = yold_up[r] - yold_c[r];
416 Deltalow[r] = yold_c[r] - yold_low[r];
417 ASSERTL0(Deltaup[r]>0,
"distance between upper and layer curve is not positive");
418 ASSERTL0(Deltalow[r]>0,
"distance between lower and layer curve is not positive");
439 if(vSession->DefinesParameter(
"npedge"))
441 npedge = (int)vSession->GetParameter(
"npedge");
449 int nq= streak->GetTotPoints();
452 streak->GetCoords(x,y);
461 xold_low, yold_low, xold_c, yold_c, x_c, y_c,cr,
true);
464 for(
int q=0; q<nvertl; q++)
466 if(y_c[q] < yold_c[q])
471 Delta_c[q] = abs(yold_c[q]-y_c[q]);
474 cout<<x_c[q]<<
" "<<y_c[q]<<endl;
479 cout<<
"Warning: the critical layer is stationary"<<endl;
502 for(
int r=0; r<nedges; r++)
505 bndSegExp = bndfieldx[lastIregion]->GetExp(r)
507 Eid = (bndSegExp->GetGeom1D())->GetEid();
508 id1 = (bndSegExp->GetGeom1D())->GetVid(0);
509 id2 = (bndSegExp->GetGeom1D())->GetVid(1);
510 vertex1 = graphShPt->GetVertex(id1);
511 vertex2 = graphShPt->GetVertex(id2);
513 vertex2->GetCoords(x2,y2,z2);
516 cout<<
"edge="<<r<<
" x1="<<x1<<
" y1="<<y1<<
" x2="<<x2<<
" y2="<<y2<<endl;
519 Cpointsx[r] = x1 +(x2-x1)/2;
522 if( Cpointsx[r]>x2 || Cpointsx[r]< x1)
524 Cpointsx[r] = -Cpointsx[r];
526 for(
int w=0; w< npedge-2; w++)
529 Addpointsx[r*(npedge-2) +w] = x1 +((x2-x1)/(npedge - 1))*(w+1);
530 if( Addpointsx[r*(npedge-2) +w] > x2 || Addpointsx[r*(npedge-2) +w] < x1)
532 Addpointsx[r*(npedge-2) +w] = -Addpointsx[r*(npedge-2) +w];
535 Addpointsy[r*(npedge-2) +w] = y_c[r] + ((y_c[r+1]-y_c[r])/(x_c[r+1]-x_c[r]))*(Addpointsx[r*(npedge-2) +w]-x1);
538 Addpointsx[r*(npedge-2) +w], Addpointsy[r*(npedge-2) +w], streak, dU,cr);
551 Cpointsx[r] = x2+ (x1-x2)/2;
553 if( Cpointsx[r] > x1 || Cpointsx[r] < x2)
555 Cpointsx[r] = -Cpointsx[r];
557 for(
int w=0; w< npedge-2; w++)
559 Addpointsx[r*(npedge-2) +w] = x2 +((x1-x2)/(npedge - 1))*(w+1);
560 if( Addpointsx[r*(npedge-2) +w] > x1 || Addpointsx[r*(npedge-2) +w] < x2)
562 Addpointsx[r*(npedge-2) +w] = -Addpointsx[r*(npedge-2) +w];
566 Addpointsy[r*(npedge-2) +w] = y_c[r+1] + ((y_c[r]-y_c[r+1])/(x_c[r]-x_c[r+1]))*(Addpointsx[r*(npedge-2) +w]-x2);
570 Addpointsx[r*(npedge-2) +w], Addpointsy[r*(npedge-2) +w], streak, dU,cr);
581 ASSERTL0(
false,
"point not generated");
600 for(
int a=0; a<nedges; a++)
603 xcPhys[a*npedge+0] = x_c[a];
604 ycPhys[a*npedge+0] = y_c[a];
606 xcPhys[a*npedge+npedge-1] = x_c[a+1];
607 ycPhys[a*npedge+npedge-1] = y_c[a+1];
609 for(
int b=0; b<npedge-2; b++)
611 xcPhys[a*npedge +b+1] = Addpointsx[a*(npedge-2)+b];
612 ycPhys[a*npedge +b+1] = Addpointsy[a*(npedge-2)+b];
616 cout<<
"xc,yc before tanevaluate"<<endl;
617 for(
int v=0; v< xcPhys.num_elements(); v++)
619 cout<<xcPhys[v]<<
" "<<ycPhys[v]<<endl;
632 MappingEVids(xold_up, yold_up, xold_low, yold_low, xold_c, yold_c, Vids_c,
633 graphShPt,streak, V1, V2, nlays, lay_Eids, lay_Vids);
637 cout<<
"nlays="<<nlays<<endl;
641 for(
int g=0; g<nlays; g++)
654 if(vSession->DefinesParameter(
"Delta"))
656 Delta0 = vSession->GetParameter(
"Delta");
668 int nVertTot = graphShPt->GetNvertices();
684 for(
int i=0; i<nVertTot; i++)
689 vertex->GetCoords(x,y,z);
699 if(x==0 && y< yold_low[0]
705 if(x== xold_c[nvertl-1] && y> yold_up[nvertl-1]
711 if(x==xold_c[nvertl-1] && y<yold_low[nvertl-1]
717 if(x== 0 && y> yold_up[0]
723 for(
int j=0; j<nvertl; j++)
725 if((xold_up[j]==x)&&(yold_up[j]==y))
729 ynew[i] = y_c[j] +Delta0;
733 if((xold_low[j]==x)&&(yold_low[j]==y))
737 ynew[i] = y_c[j] -Delta0;
741 if((xold_c[j]==x)&&(yold_c[j]==y))
752 for(
int k=0; k<nvertl; k++)
754 if(abs(x-xold_up[k]) < diff)
756 diff = abs(x-xold_up[k]);
760 if( y>yold_up[qp_closer] && y< 1)
768 ynew[i] = y_c[qp_closer] +(y-yold_c[qp_closer])*
769 (1-y_c[qp_closer])/(1-yold_c[qp_closer]);
775 else if(y<yold_low[qp_closer] && y> -1)
783 ynew[i] = y_c[qp_closer] + (y-yold_c[qp_closer] )*
784 (-1-y_c[qp_closer])/(-1-yold_c[qp_closer]);
788 else if ( y>yold_c[qp_closer] && y < yold_up[qp_closer])
795 else if (y<yold_c[qp_closer] && y > yold_low[qp_closer])
797 if(x==0){ cntlow++; }
802 else if( y==1 || y==-1)
809 if( (ynew[i]>1 || ynew[i]<-1)
810 && ( y>yold_up[qp_closer] || y<yold_low[qp_closer]) )
812 cout<<
"point x="<<xnew[i]<<
" y="<<y<<
" closer x="<<xold_up[qp_closer]<<
" ynew="<<ynew[i]<<endl;
813 ASSERTL0(
false,
"shifting out of range");
823 int nqedge = streak->GetExp(0)->GetNumPoints(0);
824 int nquad_lay = (nvertl-1)*nqedge;
829 int np_lay = (nvertl-1)*npedge;
839 if( move_norm==
true )
846 Vmath::Vcopy(xcPhys.num_elements(),xcPhys,1,xcPhysMOD,1);
847 Vmath::Vcopy(xcPhys.num_elements(),ycPhys,1,ycPhysMOD,1);
851 cout<<
"nquad per edge="<<nqedge<<endl;
852 for(
int l=0; l<2; l++)
854 Edge_newcoords[l] = bndfieldx[lastIregion]->GetExp(0)
879 bndfieldx[lastIregion]->GetCoords(xcQ, ycQ, zcQ);
888 for(
int l=0; l< xcQ.num_elements(); l++)
898 xcQ[l],4,closex,closey );
911 Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffsy);
912 Cont_y->BwdTrans_IterPerExp( coeffsy, Cont_y->UpdatePhys());
915 cout<<
"xcQ, ycQ"<<endl;
916 for(
int s=0; s<xcQ.num_elements(); s++)
918 cout<<xcQ[s]<<
" "<<ycQ[s]<<endl;
921 bool evaluatetan=
false;
925 for(
int k=0; k<nedges; k++)
930 Edge_newcoords[0]->StdPhysDeriv(xcedgeQ,txedgeQ);
931 Edge_newcoords[1]->StdPhysDeriv(ycedgeQ,tyedgeQ);
935 Vmath::Vvtvp(nqedge,tyedgeQ,1,tyedgeQ,1,normsQ,1,normsQ,1);
950 for(
int u=0; u<nqedge-1; u++)
952 incratio = (ycedgeQ[u+1]- ycedgeQ[u])/(xcedgeQ[u+1]- xcedgeQ[u]);
953 cout<<
"incratio="<<incratio<<endl;
954 if(abs(incratio)> 4.0 && evaluatetan==false )
956 cout<<
"wrong="<<wrong<<endl;
958 ASSERTL0(wrong<2,
"number edges to change is too high!!");
966 cout<<
"tan bef"<<endl;
967 for(
int e=0; e< nqedge; e++)
969 cout<<xcedgeQ[e]<<
" "<<ycedgeQ[e]<<
" "<<txedgeQ[e]<<endl;
977 Vmath::Vcopy(npedge, &xcPhysMOD[k*npedge+0],1,&xPedges[0],1);
978 Vmath::Vcopy(npedge, &ycPhysMOD[k*npedge+0],1,&yPedges[0],1);
980 PolyFit(polyorder,nqedge, xcedgeQ,ycedgeQ, coeffsinterp, xPedges,yPedges, npedge);
982 Vmath::Vcopy(npedge, &xPedges[0],1, &xcPhysMOD[k*npedge+0],1);
983 Vmath::Vcopy(npedge, &yPedges[0],1, &ycPhysMOD[k*npedge+0],1);
990 Vmath::Vcopy(nqedge, &(txedgeQ[0]), 1, &(txQ[nqedge*k]),1);
991 Vmath::Vcopy(nqedge, &(tyedgeQ[0]), 1, &(tyQ[nqedge*k]),1);
996 for(
int w=0; w< fz.num_elements(); w++)
998 txQ[w] = cos(atan(fz[w]));
999 tyQ[w] = sin(atan(fz[w]));
1000 cout<<xcQ[w]<<
" "<<ycQ[w]<<
" "<<fz[w]<<endl;
1005 Vmath::Vcopy(nquad_lay, txQ,1, Cont_y->UpdatePhys(),1);
1006 Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffsy);
1007 Cont_y->BwdTrans_IterPerExp( coeffsy, Cont_y->UpdatePhys());
1010 Vmath::Vcopy(nquad_lay, tyQ,1, Cont_y->UpdatePhys(),1);
1011 Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffsy);
1012 Cont_y->BwdTrans_IterPerExp( coeffsy, Cont_y->UpdatePhys());
1024 for(
int q=0; q<2; q++)
1026 edgebef = edgeinterp[q]-1;
1027 incbefore = (txQ[edgebef*nqedge+nqedge-1]-txQ[edgebef*nqedge])/
1028 (xcQ[edgebef*nqedge+nqedge-1]-xcQ[edgebef*nqedge]);
1029 inc = (txQ[edgeinterp[q]*nqedge+nqedge-1]-txQ[edgeinterp[q]*nqedge])/
1030 (xcQ[edgeinterp[q]*nqedge+nqedge-1]-xcQ[edgeinterp[q]*nqedge]);
1031 int npoints = 2*nqedge;
1037 cout<<
"inc="<<inc<<
" incbef="<<incbefore<<endl;
1038 if( (inc/incbefore)>0. )
1040 cout<<
"before!!"<<edgebef<<endl;
1043 Vmath::Vcopy(npoints, &xcQ[edgebef*nqedge+0],1,&xQedges[0],1);
1044 Vmath::Vcopy(npoints, &ycQ[edgebef*nqedge+0],1,&yQedges[0],1);
1045 Vmath::Vcopy(npoints, &txQ[edgebef*nqedge+0],1,&txQedges[0],1);
1046 Vmath::Vcopy(npoints, &tyQ[edgebef*nqedge+0],1,&tyQedges[0],1);
1050 coeffsinterp, xQedges,txQedges, npoints);
1053 Vmath::Vcopy(npoints,&txQedges[0],1, &txQ[edgebef*nqedge+0],1);
1057 coeffsinterp, xQedges,tyQedges, npoints);
1060 Vmath::Vcopy(npoints,&tyQedges[0],1, &tyQ[edgebef*nqedge+0],1);
1065 cout<<
"after!!"<<endl;
1068 Vmath::Vcopy(npoints, &xcQ[edgeinterp[q]*nqedge+0],1,&xQedges[0],1);
1069 Vmath::Vcopy(npoints, &ycQ[edgeinterp[q]*nqedge+0],1,&yQedges[0],1);
1070 Vmath::Vcopy(npoints, &txQ[edgeinterp[q]*nqedge+0],1,&txQedges[0],1);
1071 Vmath::Vcopy(npoints, &tyQ[edgeinterp[q]*nqedge+0],1,&tyQedges[0],1);
1076 coeffsinterp, xQedges,txQedges, npoints);
1079 Vmath::Vcopy(npoints,&txQedges[0],1, &txQ[edgeinterp[q]*nqedge+0],1);
1083 coeffsinterp, xQedges,tyQedges, npoints);
1086 Vmath::Vcopy(npoints,&tyQedges[0],1, &tyQ[edgeinterp[q]*nqedge+0],1);
1095 Vmath::Vcopy(nquad_lay, tyQ,1, Cont_y->UpdatePhys(),1);
1096 Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffstmp);
1097 Cont_y->BwdTrans_IterPerExp( coeffstmp, Cont_y->UpdatePhys());
1100 Vmath::Vcopy(nquad_lay, txQ,1, Cont_y->UpdatePhys(),1);
1101 Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffstmp);
1102 Cont_y->BwdTrans_IterPerExp( coeffstmp, Cont_y->UpdatePhys());
1105 for(
int k=0; k<nedges; k++)
1113 Vmath::Vcopy(nqedge, &(txQ[nqedge*k]),1, &(txedgeQ[0]), 1);
1114 Vmath::Vcopy(nqedge, &(tyQ[nqedge*k]),1, &(tyedgeQ[0]), 1);
1116 Vmath::Vdiv(nqedge, txedgeQ,1,tyedgeQ,1,tx_tyedgeQ,1);
1117 Vmath::Vmul(nqedge, tx_tyedgeQ,1,tx_tyedgeQ,1,tx_tyedgeQ,1);
1123 Vmath::Vcopy(nqedge, &(nxedgeQ[0]),1, &(nxQ[nqedge*k]),1);
1125 Vmath::Vmul(nqedge, nxedgeQ,1,nxedgeQ,1,nyedgeQ,1);
1131 Vmath::Vcopy(nqedge, &(nyedgeQ[0]), 1, &(nyQ[nqedge*k]),1);
1134 cout<<
"edge:"<<k<<endl;
1135 cout<<
"tan/normal"<<endl;
1136 for(
int r=0; r<nqedge; r++)
1138 cout<<xcQ[k*nqedge+r]<<
" "<<txedgeQ[r]<<
" "<<tyedgeQ[r]<<
" "
1139 <<nxedgeQ[r]<<
" "<<nyedgeQ[r]<<endl;
1145 Vmath::Vcopy(nquad_lay, nyQ,1, Cont_y->UpdatePhys(),1);
1147 Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffstmp);
1148 Cont_y->BwdTrans_IterPerExp( coeffstmp, Cont_y->UpdatePhys());
1152 Vmath::Zero(Cont_y->GetNcoeffs(),Cont_y->UpdateCoeffs(),1);
1153 Vmath::Vcopy(nquad_lay, nxQ,1, Cont_y->UpdatePhys(),1);
1154 Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffstmp);
1155 Cont_y->BwdTrans_IterPerExp( coeffstmp, Cont_y->UpdatePhys());
1159 for(
int k=0; k<nedges; k++)
1165 nyQ[(k-1)*nqedge+nqedge-1]=
1170 nxQ[(k-1)*nqedge+nqedge-1]=
1179 cout<<
"nx,yQbefore"<<endl;
1180 for(
int u=0; u<xcQ.num_elements(); u++)
1182 cout<<xcQ[u]<<
" "<<nyQ[u]<<
" "<<txQ[u]<<endl;
1188 cout<<
"nx,yQ"<<endl;
1189 for(
int u=0; u<x_tmpQ.num_elements(); u++)
1191 cout<<x_tmpQ[u]<<
" "<<tmpnyQ[u]<<endl;
1195 for(
int k=0; k<nedges; k++)
1198 for(
int a=0; a<npedge; a++)
1202 nxPhys[k*npedge +a]= nxQ[k*nqedge +0];
1203 nyPhys[k*npedge +a]= nyQ[k*nqedge +0];
1206 else if(a== npedge-1)
1208 nxPhys[k*npedge +a]= nxQ[k*nqedge +nqedge-1];
1209 nyPhys[k*npedge +a]= nyQ[k*nqedge +nqedge-1];
1233 nyPhys[k*npedge +a]=
1243 nxPhys[k*npedge +a]= -sqrt(abs(1- nyPhys[k*npedge +a]*nyPhys[k*npedge +a]));
1259 nyPhys[(k-1)*npedge+npedge-1]=
1264 nxPhys[(k-1)*npedge+npedge-1]=
1269 cout<<
"xcPhys,,"<<endl;
1270 for(
int s=0; s<np_lay; s++)
1273 cout<<xcPhysMOD[s]<<
" "<<ycPhysMOD[s]<<
" "<<nxPhys[s]<<
" "<<nyPhys[s]<<endl;
1286 for(
int m=0; m<nlays; m++)
1293 delta[m] = -(cntlow+1-m)*Delta0/(cntlow+1);
1297 delta[m] = ( m-(cntlow) )*Delta0/(cntlow+1);
1304 for(
int h=0; h< nvertl; h++)
1311 if(move_norm==
false)
1313 ynew[lay_Vids[m][h] ]= y_c[h] +delta[m];
1314 xnew[lay_Vids[m][h] ]= x_c[h];
1318 if(h==0 || h==nvertl-1 )
1320 ynew[lay_Vids[m][h] ]= y_c[h] +delta[m];
1321 xnew[lay_Vids[m][h] ]= x_c[h];
1325 ynew[lay_Vids[m][h] ]= y_c[h] +delta[m]*abs(nyPhys[h*npedge+0]);
1326 xnew[lay_Vids[m][h] ]= x_c[h] +delta[m]*abs(nxPhys[h*npedge+0]);
1329 cout<<
"Vid x="<<xnew[lay_Vids[m][h] ]<<
" y="<<ynew[lay_Vids[m][h] ]<<endl;
1334 cout<<
"edge=="<<h<<endl;
1337 ASSERTL0( nyPhys[h*npedge+0]==nyPhys[(h-1)*npedge+npedge-1],
" normaly wrong");
1338 ASSERTL0( nxPhys[h*npedge+0]==nxPhys[(h-1)*npedge+npedge-1],
" normalx wrong");
1341 if(move_norm==
false)
1344 layers_y[m][h*npedge +0] = y_c[h] +delta[m];
1345 layers_x[m][h*npedge +0] = xnew[lay_Vids[m][h] ];
1347 layers_y[m][h*npedge +npedge-1] = y_c[h+1] +delta[m];
1348 layers_x[m][h*npedge +npedge-1] = xnew[lay_Vids[m][h+1] ];
1350 for(
int d=0; d< npedge-2; d++)
1352 layers_y[m][h*npedge +d+1]= ycPhysMOD[h*npedge +d+1] +delta[m];
1354 layers_x[m][h*npedge +d+1]= xcPhysMOD[h*npedge +d+1];
1363 tmpy_lay[h*npedge +0] = y_c[h] +delta[m];
1364 tmpx_lay[h*npedge +0] = xnew[lay_Vids[m][h] ];
1366 tmpy_lay[h*npedge +npedge-1] =
1367 y_c[h+1] +delta[m]*abs(nyPhys[h*npedge +npedge-1]);
1368 tmpx_lay[h*npedge +npedge-1] =
1369 x_c[h+1] +delta[m]*abs(nxPhys[h*npedge +npedge-1]);
1371 else if(h==nedges-1)
1374 tmpy_lay[h*npedge +0] =
1375 y_c[h] +delta[m]*abs(nyPhys[h*npedge +0]);
1376 tmpx_lay[h*npedge +0] =
1377 x_c[h] +delta[m]*abs(nxPhys[h*npedge +0]);
1379 tmpy_lay[h*npedge +npedge-1] = y_c[h+1] +delta[m];
1380 tmpx_lay[h*npedge +npedge-1] = xnew[lay_Vids[m][h+1] ];
1385 tmpy_lay[h*npedge +0] =
1386 y_c[h] +delta[m]*abs(nyPhys[h*npedge +0]);
1387 tmpx_lay[h*npedge +0] =
1388 x_c[h] +delta[m]*abs(nxPhys[h*npedge +0]);
1390 tmpy_lay[h*npedge +npedge-1] =
1391 y_c[h+1] +delta[m]*abs(nyPhys[h*npedge +npedge-1]);
1392 tmpx_lay[h*npedge +npedge-1] =
1393 x_c[h+1] +delta[m]*abs(nxPhys[h*npedge +npedge-1]);
1397 for(
int d=0; d< npedge-2; d++)
1400 tmpy_lay[h*npedge +d+1] = ycPhysMOD[h*npedge +d+1] +
1401 delta[m]*abs(nyPhys[h*npedge +d+1]);
1404 tmpx_lay[h*npedge +d+1]= xcPhysMOD[h*npedge +d+1] +
1405 delta[m]*abs(nxPhys[h*npedge +d+1]);
1422 for(
int s=0; s<np_lay; s++)
1424 cout<<tmpx_lay[s]<<
" "<<tmpy_lay[s]<<endl;
1427 cout<<
"fisrt interp"<<endl;
1428 for(
int s=0; s<np_lay; s++)
1430 cout<<tmpx_lay[s]<<
" "<<tmpy_lay[s]<<endl;
1442 NekDouble boundright = xcPhysMOD[np_lay-1];
1443 bool outboundleft=
false;
1444 bool outboundright=
false;
1445 if(tmpx_lay[1]< boundleft )
1447 outboundleft =
true;
1449 if(tmpx_lay[np_lay-2] > boundright )
1451 outboundright =
true;
1459 for(
int r=0; r< nedges; r++)
1462 if(tmpx_lay[r*npedge + npedge-1]< boundleft && outboundleft==
true )
1470 if( tmpx_lay[(r+1)*npedge + npedge-1]> boundleft )
1472 for(
int s=0; s<npedge-2; s++)
1474 if(tmpx_lay[(r+1)*npedge + s+1]< boundleft)
1484 if(tmpx_lay[r*npedge + 0]> boundright && outboundright==
true )
1492 if( tmpx_lay[(r-1)*npedge + 0]< boundright )
1494 for(
int s=0; s<npedge-2; s++)
1496 if(tmpx_lay[(r-1)*npedge + s+1]> boundright)
1508 outcount = outvert*npedge+1+ outmiddle;
1510 int replacepointsfromindex=0;
1511 for(
int c=0; c<nedges; c++)
1514 if(xcPhysMOD[c*npedge+npedge-1] <= tmpx_lay[c*(npedge-(npedge-2)) +2] && outboundright==
true)
1516 replacepointsfromindex = c*(npedge-(npedge-2))+2;
1521 if(xcPhysMOD[(nedges-1 -c)*npedge+0] >= tmpx_lay[np_lay-1 -(c*(npedge-(npedge-2)) +2)] && outboundleft==
true)
1523 replacepointsfromindex = np_lay-1 -(c*(npedge-(npedge-2)) +2);
1539 if( outboundright==
true)
1541 pstart = replacepointsfromindex;
1542 shift = np_lay-outcount;
1543 increment = (xcPhysMOD[np_lay-outcount]-xcPhysMOD[pstart])/(outcount+1);
1544 outcount = outcount-1;
1545 ASSERTL0(tmpx_lay[np_lay-outcount]>xcPhysMOD[(nedges-1)*npedge+0],
"no middle points in the last edge");
1551 increment = (xcPhysMOD[replacepointsfromindex]-xcPhysMOD[pstart])/(outcount+1);
1552 ASSERTL0(tmpx_lay[pstart]<xcPhysMOD[0*npedge +npedge-1],
"no middle points in the first edge");
1569 NekDouble xctmp,ycinterp,nxinterp,nyinterp;
1571 for(
int v=0; v<outcount;v++)
1573 xctmp = xcPhysMOD[pstart]+(v+1)*increment;
1586 xctmp,4,closex,closeny );
1589 nxinterp = sqrt(abs(1-nyinterp*nyinterp));
1596 replace_x[v] = xctmp +delta[m]*abs(nxinterp);
1597 replace_y[v] = ycinterp +delta[m]*abs(nyinterp);
1598 tmpx_lay[ v+shift ] = replace_x[v];
1599 tmpy_lay[ v+shift ] = replace_y[v];
1620 int closepoints = 4;
1627 for(
int q=0; q<np_lay; q++)
1629 for(
int e=0; e<nedges; e++)
1631 if(tmpx_lay[q]<= x_c[e+1] && tmpx_lay[q]>= x_c[e])
1635 if(q == e*npedge +npedge-1 && pointscount!=npedge )
1640 else if(q == e*npedge +npedge-1)
1660 lay_Vids[m], layers_x[m], layers_y[m],xnew,ynew);
1743 int npoints = npedge;
1746 for(
int f=0; f<nedges; f++)
1752 Vmath::Vcopy(npoints, &layers_x[m][(f)*npedge+0],1,&xPedges[0],1);
1753 Vmath::Vcopy(npoints, &layers_y[m][(f)*npedge+0],1,&yPedges[0],1);
1757 coeffsinterp, xPedges,yPedges, npoints);
1760 Vmath::Vcopy(npoints,&yPedges[0],1, &layers_y[m][(f)*npedge+0],1);
1763 layers_y[m][f*npedge+0]= ynew[lay_Vids[m][f]];
1764 layers_y[m][f*npedge+npedge-1]= ynew[lay_Vids[m][f+1]];
1767 cout<<
" xlay ylay lay:"<<m<<endl;
1768 for(
int l=0; l<np_lay; l++)
1771 cout<<std::setprecision(8)<<layers_x[m][l]<<
" "<<layers_y[m][l]<<endl;
1805 cout<<
"lay="<<m<<endl;
1807 " different layer ymin val");
1809 " different layer ymax val");
1811 " different layer xmin val");
1813 " different layer xmax val");
1823 layers_x[0], layers_y[0], layers_x[nlays-1], layers_y[nlays-1],nxPhys, nyPhys,xnew, ynew);
1905 lay_Vids, x_c, y_c, Down, Up, xnew, ynew, layers_x, layers_y);
1916 cout<<std::setprecision(8)<<
"xmin="<<
Vmath::Vmin(nVertTot, xnew,1)<<endl;
1918 " different xmin val");
1920 " different ymin val");
1922 " different xmax val");
1924 " different ymax val");
1930 Replacevertices(changefile, xnew , ynew, xcPhys, ycPhys, Eids, npedge, charalp, layers_x,layers_y, lay_Eids, curv_lay);
void PolyFit(int polyorder, int npoints, Array< OneD, NekDouble > xin, Array< OneD, NekDouble > fin, Array< OneD, NekDouble > &coeffsinterp, Array< OneD, NekDouble > &xout, Array< OneD, NekDouble > &fout, int npout)
#define ASSERTL0(condition, msg)
boost::shared_ptr< ContField1D > ContField1DSharedPtr
void OrderVertices(int nedges, SpatialDomains::MeshGraphSharedPtr graphShPt, MultiRegions::ExpListSharedPtr &bndfield, Array< OneD, int > &Vids, int v1, int v2, NekDouble x_connect, int &lastedge, Array< OneD, NekDouble > &x, Array< OneD, NekDouble > &y)
void GenerateMapEidsv1v2(MultiRegions::ExpListSharedPtr field, Array< OneD, int > &V1, Array< OneD, int > &V2)
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
void MoveOutsidePointsNnormpos(int npedge, SpatialDomains::MeshGraphSharedPtr mesh, Array< OneD, NekDouble > xcold, Array< OneD, NekDouble > ycold, Array< OneD, NekDouble > xolddown, Array< OneD, NekDouble > yolddown, Array< OneD, NekDouble > xoldup, Array< OneD, NekDouble > yoldup, Array< OneD, NekDouble > xlaydown, Array< OneD, NekDouble > ylaydown, Array< OneD, NekDouble > xlayup, Array< OneD, NekDouble > ylayup, Array< OneD, NekDouble > nxPhys, Array< OneD, NekDouble > nyPhys, Array< OneD, NekDouble > &xnew, Array< OneD, NekDouble > &ynew)
void GenerateAddPointsNewtonIt(NekDouble xi, NekDouble yi, NekDouble &xout, NekDouble &yout, MultiRegions::ExpListSharedPtr function, Array< OneD, NekDouble > derfunction, NekDouble cr)
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...
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
void MoveLayerNnormpos(int nvertl, int npedge, Array< OneD, NekDouble > xcPhys, Array< OneD, NekDouble > tmpx_lay, Array< OneD, NekDouble > tmpy_lay, Array< OneD, int > Vids, Array< OneD, NekDouble > &xlay, Array< OneD, NekDouble > &ylay, Array< OneD, NekDouble > &xnew, Array< OneD, NekDouble > &ynew)
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
void MappingEVids(Array< OneD, NekDouble > xoldup, Array< OneD, NekDouble > yoldup, Array< OneD, NekDouble > xolddown, Array< OneD, NekDouble > yolddown, Array< OneD, NekDouble > xcold, Array< OneD, NekDouble > ycold, Array< OneD, int > Vids_c, SpatialDomains::MeshGraphSharedPtr mesh, MultiRegions::ExpListSharedPtr streak, Array< OneD, int > V1, Array< OneD, int > V2, int &nlays, Array< OneD, Array< OneD, int > > &Eids_lay, Array< OneD, Array< OneD, int > > &Vids_lay)
void GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2=NullNekDouble1DArray, Array< OneD, NekDouble > &coords_3=NullNekDouble1DArray)
this function returns the physical coordinates of the quadrature points of the expansion ...
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 Cutrepetitions(int nedges, Array< OneD, NekDouble > inarray, Array< OneD, NekDouble > &outarray)
boost::shared_ptr< SegExp > SegExpSharedPtr
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
NekDouble LagrangeInterpolant(NekDouble x, int npts, Array< OneD, NekDouble > xpts, Array< OneD, NekDouble > funcvals)
void CheckSingularQuads(MultiRegions::ExpListSharedPtr Exp, Array< OneD, int > V1, Array< OneD, int > V2, Array< OneD, NekDouble > &xnew, Array< OneD, NekDouble > &ynew)
void EvaluateTangent(int npoints, Array< OneD, NekDouble > xcQedge, Array< OneD, NekDouble > coeffsinterp, Array< OneD, NekDouble > &txQedge, Array< OneD, NekDouble > &tyQedge)
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
void Computestreakpositions(int nvertl, MultiRegions::ExpListSharedPtr streak, Array< OneD, NekDouble > xold_up, Array< OneD, NekDouble > yold_up, Array< OneD, NekDouble > xold_low, Array< OneD, NekDouble > yold_low, Array< OneD, NekDouble > xold_c, Array< OneD, NekDouble > yold_c, Array< OneD, NekDouble > &xc, Array< OneD, NekDouble > &yc, NekDouble cr, bool verts)
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
void Orderfunctionx(Array< OneD, NekDouble > inarray_x, Array< OneD, NekDouble > inarray_y, Array< OneD, NekDouble > &outarray_x, Array< OneD, NekDouble > &outarray_y)
boost::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
int DetermineclosePointxindex(NekDouble x, Array< OneD, NekDouble > xArray)
void MoveLayersvertically(int nlays, int nvertl, int cntlow, int cntup, Array< OneD, Array< OneD, int > > lay_Vids, Array< OneD, NekDouble > xc, Array< OneD, NekDouble > yc, Array< OneD, int > Down, Array< OneD, int > Up, Array< OneD, NekDouble > &xnew, Array< OneD, NekDouble > &ynew, Array< OneD, Array< OneD, NekDouble > > &layers_x, Array< OneD, Array< OneD, NekDouble > > &layers_y)
void Zero(int n, T *x, const int incx)
Zero vector.
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
void Replacevertices(string filename, Array< OneD, NekDouble > newx, Array< OneD, NekDouble > newy, Array< OneD, NekDouble > xcPhys, Array< OneD, NekDouble > ycPhys, Array< OneD, int >Eids, int Npoints, string s_alp, Array< OneD, Array< OneD, NekDouble > > x_lay, Array< OneD, Array< OneD, NekDouble > > y_lay, Array< OneD, Array< OneD, int > >lay_eids, bool curv_lay)
boost::shared_ptr< PointGeom > PointGeomSharedPtr
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void GenerateNeighbourArrays(int index, int neighpoints, Array< OneD, NekDouble > xArray, Array< OneD, NekDouble > yArray, Array< OneD, NekDouble > &Neighbour_x, Array< OneD, NekDouble > &Neighbour_y)