63 #include <boost/lexical_cast.hpp>
116 int edge,
int npedge);
117 void PolyFit(
int polyorder,
int npoints,
173 int main(
int argc,
char *argv[])
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=0;
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);
1941 int nvertl = nedges+1;
1945 for(
int j=0; j<nedges; j++)
1949 edge = (bndSegExplow->GetGeom1D())->GetEid();
1951 for(
int k=0; k<2; k++)
1953 Vids_temp[j+k]=(bndSegExplow->GetGeom1D())->GetVid(k);
1956 vertex->GetCoords(x1,y1,z1);
1957 if(x1==x_connect && edge!=lastedge)
1960 if(x_connect==x0layer)
1962 Vids[v1]=Vids_temp[j+k];
1968 Vids_temp[j+1]=(bndSegExplow->GetGeom1D())->GetVid(1);
1969 Vids[v2]=Vids_temp[j+1];
1972 vertex->GetCoords(x2,y2,z2);
1978 Vids_temp[j+0]=(bndSegExplow->GetGeom1D())->GetVid(0);
1979 Vids[v2]=Vids_temp[j+0];
1982 vertex->GetCoords(x2,y2,z2);
1991 Vids_temp[j+1]=(bndSegExplow->GetGeom1D())->GetVid(1);
1992 Vids[v1]=Vids_temp[j+1];
1995 vertex->GetCoords(x1,y1,z1);
2001 Vids_temp[j+0]=(bndSegExplow->GetGeom1D())->GetVid(0);
2002 Vids[v1]=Vids_temp[j+0];
2005 vertex->GetCoords(x1,y1,z1);
2030 cout<<
"Computestreakpositions"<<endl;
2031 int nq = streak->GetTotPoints();
2045 Vmath::Vadd(xc.num_elements(), yold_up,1,yold_low,1, yc,1);
2059 for(
int e=0; e<npoints; e++)
2064 elmtid = streak->GetExpIndex(coord,0.00001);
2065 offset = streak->GetPhys_Offset(elmtid);
2071 while( abs(F)> 0.000000001)
2074 elmtid = streak->GetExpIndex(coord,0.00001);
2079 if( (abs(coord[1])>1 || elmtid==-1)
2080 && attempt==0 && verts==
true
2084 coord[1] = yold_c[e];
2087 else if( (abs(coord[1])>1 || elmtid==-1) )
2089 coord[1] = ytmp +0.01;
2090 elmtid = streak->GetExpIndex(coord,0.001);
2091 offset = streak->GetPhys_Offset(elmtid);
2092 NekDouble Utmp = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2093 NekDouble dUtmp = streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
2094 coord[1] = coord[1] - (Utmp-cr)/dUtmp;
2095 if( (abs(Utmp-cr)>abs(F))||(abs(coord[1])>1) )
2097 coord[1] = ytmp -0.01;
2104 ASSERTL0(abs(coord[1])<= 1,
" y value out of bound +/-1");
2106 offset = streak->GetPhys_Offset(elmtid);
2107 U = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2108 dU = streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
2109 coord[1] = coord[1] - (U-cr)/dU;
2111 ASSERTL0( coord[0]==xc[e],
" x coordinate must remain the same");
2114 if(its>200 && abs(F)<0.00001)
2116 cout<<
"warning streak position obtained with precision:"<<F<<endl;
2119 else if(its>1000 && abs(F)< 0.0001)
2121 cout<<
"warning streak position obtained with precision:"<<F<<endl;
2126 ASSERTL0(
false,
"no convergence after 1000 iterations");
2129 yc[e] = coord[1] - (U-cr)/dU;
2130 ASSERTL0( U<= cr + tol,
"streak wrong+");
2131 ASSERTL0( U>= cr -tol,
"streak wrong-");
2133 cout<<
"result streakvert x="<<xc[e]<<
" y="<<yc[e]<<
" streak="<<U<<endl;
2154 while( abs(F)> 0.00000001)
2158 elmtid =
function->GetExpIndex(coords, 0.01);
2160 cout<<
"gen newton xi="<<xi<<
" yi="<<coords[1]<<
" elmtid="<<elmtid<<
" F="<<F<<endl;
2162 if( (abs(coords[1])>1 || elmtid==-1) )
2165 coords[1] = ytmp +0.01;
2166 elmtid =
function->GetExpIndex(coords,0.01);
2167 offset =
function->GetPhys_Offset(elmtid);
2168 NekDouble Utmp =
function->GetExp(elmtid)->PhysEvaluate(coords, function->GetPhys() + offset);
2169 NekDouble dUtmp =
function->GetExp(elmtid)->PhysEvaluate(coords, derfunction + offset);
2170 coords[1] = coords[1] - (Utmp-cr)/dUtmp;
2171 cout<<
"attempt:"<<coords[1]<<endl;
2172 if( (abs(Utmp-cr)>abs(F))||(abs(coords[1])>1.01) )
2174 coords[1] = ytmp -0.01;
2179 else if( abs(coords[1])<1.01 &&attempt==0)
2186 ASSERTL0(abs(coords[1])<= 1.00,
" y value out of bound +/-1");
2188 offset =
function->GetPhys_Offset(elmtid);
2189 U =
function->GetExp(elmtid)->PhysEvaluate(coords, function->GetPhys() + offset);
2190 dU =
function->GetExp(elmtid)->PhysEvaluate(coords, derfunction + offset);
2191 coords[1] = coords[1] - (U-cr)/dU;
2192 cout<<cr<<
"U-cr="<<U-cr<<
" tmp result y:"<<coords[1]<<
" dU="<<dU<<endl;
2196 if(its>200 && abs(F)<0.00001)
2198 cout<<
"warning streak position obtained with precision:"<<F<<endl;
2201 else if(its>1000 && abs(F)< 0.0001)
2203 cout<<
"warning streak position obtained with precision:"<<F<<endl;
2208 ASSERTL0(
false,
"no convergence after 1000 iterations");
2212 ASSERTL0( coords[0]==xi,
" x coordinate must remain the same");
2215 yout = coords[1] - (U-cr)/dU;
2216 cout<<
"NewtonIt result x="<<xout<<
" y="<<coords[1]<<
" U="<<U<<endl;
2223 const boost::shared_ptr<LocalRegions::ExpansionVector> exp2D = field->GetExp();
2224 int nel = exp2D->size();
2232 for(
int i=0; i<nel; i++)
2234 if((locQuadExp = (*exp2D)[i]->as<LocalRegions::QuadExp>()))
2236 for(
int j = 0; j < locQuadExp->GetNedges(); ++j)
2238 SegGeom = (locQuadExp->GetGeom2D())->GetEdge(j);
2239 id = SegGeom->GetEid();
2240 if( V1tmp[
id] == 10000)
2242 V1tmp[id]= SegGeom->GetVertex(0)->GetVid();
2243 V2tmp[id]= SegGeom->GetVertex(1)->GetVid();
2250 else if((locTriExp = (*exp2D)[i]->as<LocalRegions::TriExp>()))
2252 for(
int j = 0; j < locTriExp->GetNedges(); ++j)
2254 SegGeom = (locTriExp->GetGeom2D())->GetEdge(j);
2255 id = SegGeom->GetEid();
2257 if( V1tmp[
id] == 10000)
2259 V1tmp[id]= SegGeom->GetVertex(0)->GetVid();
2260 V2tmp[id]= SegGeom->GetVertex(1)->GetVid();
2272 for(
int g=0; g<cnt; g++)
2275 ASSERTL0(V1tmp[g]!=10000,
"V1 wrong");
2277 ASSERTL0(V2tmp[g]!=10000,
"V2 wrong");
2300 int nlay_Eids = xcold.num_elements()-1;
2301 int nlay_Vids = xcold.num_elements();
2303 int nVertsTot = mesh->GetNvertices();
2304 cout<<
"nverttot="<<nVertsTot<<endl;
2308 cout<<
"init nlays="<<nlays<<endl;
2315 cout<<
"yoldup="<<yoldup[0]<<endl;
2316 cout<<
"yolddown="<<yolddown[0]<<endl;
2318 for(
int r=0; r< nVertsTot; r++)
2323 vertex->GetCoords(x,y,z);
2330 y<= yoldup[0] && y>= yolddown[0]
2343 cout<<
"nlays="<<nlays<<endl;
2355 for(
int w=0; w< nlays; w++)
2358 tmpx0[w]= tmpx[index];
2359 tmpy0[w]= tmpf[index];
2360 tmpVids0[w] = tmpV[index];
2361 tmpf[index] = max+1000;
2372 for(
int m=0; m<nlays; m++)
2383 NekDouble xtmp,ytmp,normnext=0.0,xnext=0.0,ynext=0.0,diff;
2384 NekDouble Ubef = 0.0, Utmp = 0.0, Unext = 0.0;
2387 int nTotEdges = V1.num_elements();
2389 for(
int m=0; m<nlays; m++)
2391 for(
int g=0; g<nlay_Eids; g++)
2395 for(
int h=0; h< nTotEdges; h++)
2398 if( tmpVids0[m]== V1[h] )
2402 vertex->GetCoords(x,y,z);
2406 Vids_lay[m][0] = V1[h];
2407 Vids_lay[m][1] = V2[h];
2409 = mesh->GetVertex(V1[h]);
2411 vertex1->GetCoords(x1,y1,z1);
2412 normbef= sqrt( (y-y1)*(y-y1)+(x-x1)*(x-x1) );
2417 elmtid = streak->GetExpIndex(coord,0.00001);
2418 offset = streak->GetPhys_Offset(elmtid);
2419 Ubef = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2424 if( tmpVids0[m]== V2[h] )
2428 vertex->GetCoords(x,y,z);
2432 Vids_lay[m][0] = V2[h];
2433 Vids_lay[m][1] = V1[h];
2435 = mesh->GetVertex(V2[h]);
2437 normbef= sqrt( (y-y2)*(y-y2)+(x-x2)*(x-x2) );
2442 elmtid = streak->GetExpIndex(coord,0.00001);
2443 offset = streak->GetPhys_Offset(elmtid);
2444 Ubef = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2451 cout<<
"Eid="<<Eids_lay[m][0]<<
" Vids_lay0="<<Vids_lay[m][0]<<
" Vidslay1="<<Vids_lay[m][1]<<endl;
2458 for(
int h=0; h< nTotEdges; h++)
2461 if( (Vids_lay[m][g]==V1[h] || Vids_lay[m][g]==V2[h]) && h!= Eids_lay[m][g-1])
2463 cout<<
"edgetmp="<<h<<endl;
2464 ASSERTL0(cnt<=6,
"wrong number of candidates");
2473 cout<<
"normbef="<<normbef<<endl;
2474 cout<<
"Ubefcc="<<Ubef<<endl;
2476 for(
int e=0; e< cnt; e++)
2480 vertex1->GetCoords(x1,y1,z1);
2483 vertex2->GetCoords(x2,y2,z2);
2485 normtmp= sqrt( (y2-y1)*(y2-y1)+(x2-x1)*(x2-x1) );
2487 cout<<
"edgetmp1="<<edgestmp[e]<<endl;
2488 cout<<
"V1 x1="<<x1<<
" y1="<<y1<<endl;
2489 cout<<
"V2 x2="<<x2<<
" y2="<<y2<<endl;
2490 if( Vids_lay[m][g]==V1[edgestmp[e]] )
2498 elmtid = streak->GetExpIndex(coord,0.00001);
2499 offset = streak->GetPhys_Offset(elmtid);
2500 Utmp = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2501 diffarray[e] = abs((xtmp*xbef+ytmp*ybef)/(normtmp*normbef)-1);
2502 diffUarray[e] = abs(Ubef-Utmp);
2503 cout<<
" normtmp="<<normtmp<<endl;
2504 cout<<
" Utmpcc="<<Utmp<<endl;
2505 cout<<xtmp<<
" ytmp="<<ytmp<<
" diff="<<abs(((xtmp*xbef+ytmp*ybef)/(normtmp*normbef))-1)<<endl;
2507 abs( (xtmp*xbef+ytmp*ybef)/(normtmp*normbef)-1)<diff
2508 && y2<= yoldup[g+1] && y2>= yolddown[g+1]
2509 && y1<= yoldup[g] && y1>= yolddown[g]
2513 Eids_lay[m][g] = edgestmp[e];
2514 Vids_lay[m][g+1] = V2[edgestmp[e]];
2515 diff = abs((xtmp*xbef+ytmp*ybef)/(normtmp*normbef)-1);
2522 else if( Vids_lay[m][g]==V2[edgestmp[e]] )
2530 elmtid = streak->GetExpIndex(coord,0.00001);
2531 offset = streak->GetPhys_Offset(elmtid);
2532 Utmp = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2533 diffarray[e] = abs((xtmp*xbef+ytmp*ybef)/(normtmp*normbef)-1);
2534 diffUarray[e] = abs(Ubef-Utmp);
2535 cout<<
" normtmp="<<normtmp<<endl;
2536 cout<<
" Utmpcc="<<Utmp<<endl;
2537 cout<<xtmp<<
" ytmp="<<ytmp<<
" diff="<<abs(((xtmp*xbef+ytmp*ybef)/(normtmp*normbef))-1)<<endl;
2539 abs((xtmp*xbef+ytmp*ybef)/(normtmp*normbef)-1)<diff
2540 && y2<= yoldup[g] && y2>= yolddown[g]
2541 && y1<= yoldup[g+1] && y1>= yolddown[g+1]
2544 Eids_lay[m][g] = edgestmp[e];
2545 Vids_lay[m][g+1] = V1[edgestmp[e]];
2546 diff = abs((xtmp*xbef+ytmp*ybef)/(normtmp*normbef)-1);
2562 cout<<
"Eid before check="<<Eids_lay[m][g]<<endl;
2563 for(
int q=0; q<cnt; q++)
2565 cout<<q<<
" diff"<<diffarray[q]<<endl;
2575 cout<<
"COMMON VERT"<<endl;
2577 diffarray[eid]=1000;
2583 vertex1->GetCoords(x1,y1,z1);
2586 vertex2->GetCoords(x2,y2,z2);
2588 normtmp= sqrt( (y2-y1)*(y2-y1)+(x2-x1)*(x2-x1) );
2590 Eids_lay[m][g] = edgestmp[eid];
2591 if(Vids_lay[m][g] == V1[edgestmp[eid]])
2595 elmtid = streak->GetExpIndex(coord,0.00001);
2596 offset = streak->GetPhys_Offset(elmtid);
2597 Utmp = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2598 Vids_lay[m][g+1] = V2[edgestmp[eid]];
2606 if(Vids_lay[m][g] == V2[edgestmp[eid]])
2610 elmtid = streak->GetExpIndex(coord,0.00001);
2611 offset = streak->GetPhys_Offset(elmtid);
2612 Utmp = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2613 Vids_lay[m][g+1] = V1[edgestmp[eid]];
2624 cout<<m<<
"edge aft:"<<Eids_lay[m][g]<<
" Vid="<<Vids_lay[m][g+1]<<endl;
2630 cout<<
"endelse"<<normtmp<<endl;
2642 for(
int w=0; w< nlays; w++)
2644 for(
int f=0; f< nlay_Eids; f++)
2646 cout<<
"check="<<w<<
" Eid:"<<Eids_lay[w][f]<<endl;
2656 for(
int u=0; u< Vids_laybefore.num_elements(); u++)
2658 if( Vids_laybefore[u]==Vid || Vids_c[u]==Vid)
2662 cout<<Vid<<
" Vert test="<<Vids_laybefore[u]<<endl;
2674 int np_lay = inarray.num_elements();
2675 ASSERTL0(inarray.num_elements()%nedges==0,
" something on number npedge");
2679 for(
int w=0; w< np_lay; w++)
2684 if(inarray[w] ==inarray[w+1])
2689 outarray[cnt]= inarray[w];
2694 ASSERTL0( cnt== np_lay-(nedges-1),
"wrong cut");
2700 int npts = xArray.num_elements();
2709 if(xArray[index]> x)
2720 ASSERTL0( neighpoints%2==0,
"number of neighbour points should be even");
2721 int leftpoints = (neighpoints/2)-1;
2722 int rightpoints = neighpoints/2;
2726 if(index-leftpoints<0)
2729 diff = index-leftpoints;
2731 Vmath::Vcopy(neighpoints, &yArray[0],1,&Neighbour_y[0],1);
2732 Vmath::Vcopy(neighpoints, &xArray[0],1,&Neighbour_x[0],1);
2734 else if( (yArray.num_elements()-1)-index < rightpoints)
2737 int rpoints = (yArray.num_elements()-1)-index;
2738 diff = rightpoints-rpoints;
2740 start = index-leftpoints-diff;
2741 Vmath::Vcopy(neighpoints, &yArray[start],1,&Neighbour_y[0],1);
2742 Vmath::Vcopy(neighpoints, &xArray[start],1,&Neighbour_x[0],1);
2747 start = index-leftpoints;
2748 Vmath::Vcopy(neighpoints, &yArray[start],1,&Neighbour_y[0],1);
2749 Vmath::Vcopy(neighpoints, &xArray[start],1,&Neighbour_x[0],1);
2758 for(
int f=1; f< neighpoints; f++)
2760 ASSERTL0(Neighbour_x[f]!=Neighbour_x[f-1],
" repetition on NeighbourArrays");
2771 for(
int pt=0;pt<
npts;++pt)
2775 for(
int j=0;j<pt; ++j)
2777 h = h * (x - xpts[j])/(xpts[pt]-xpts[j]);
2780 for(
int k=pt+1;k<
npts;++k)
2782 h = h * (x - xpts[k])/(xpts[pt]-xpts[k]);
2786 sum += funcvals[pt]*LagrangePoly;
2798 int np_pol= coeffsinterp.num_elements();
2799 cout<<
"evaluatetan with "<<np_pol<<endl;
2805 for(
int q=0; q< npoints; q++)
2810 for(
int d=0; d< np_pol-1; d++)
2812 yprime[q] += (derorder +1)*coeffsinterp[d]*std::pow(xcQedge[q],derorder);
2816 for(
int a=0; a< np_pol; a++)
2818 yinterp[q] += coeffsinterp[a]*std::pow(xcQedge[q],polorder);
2826 for(
int n=0; n< npoints; n++)
2830 txQedge[n] = cos((atan((yprime[n]))));
2831 tyQedge[n] = sin((atan((yprime[n]))));
2832 cout<<xcQedge[n]<<
" "<<yinterp[n]<<
" "<<yprime[n]<<
" "<<txQedge[n]<<
" "<<tyQedge[n]<<endl;
2839 int edge,
int npedge)
2841 int np_pol = xpol.num_elements();
2847 for(
int e=0; e<N; e++)
2850 for(
int w=0; w < N; w++)
2852 A[N*e+row] = std::pow( xpol[w], N-1-e);
2857 for(
int r= 0; r< np_pol; r++)
2868 Lapack::Dgetrf( N, N, A.get(), N, ipivot.get(), info);
2871 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
2872 "th parameter had an illegal parameter for dgetrf";
2877 std::string message =
"ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
2878 boost::lexical_cast<std::string>(info) +
" is 0 from dgetrf";
2884 Lapack::Dgetrs(
'N', N, ncolumns_b , A.get() , N, ipivot.get(), b.get(), N, info);
2887 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
2888 "th parameter had an illegal parameter for dgetrf";
2893 std::string message =
"ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
2894 boost::lexical_cast<std::string>(info) +
" is 0 from dgetrf";
2908 for(
int c=0; c< npedge; c++)
2912 ycout[edge*(npedge)+c+1]=0;
2913 for(
int d=0; d< np_pol; d++)
2915 ycout[edge*(npedge)+c+1] += b[d]
2916 *std::pow(xcout[edge*(npedge)+c+1],polorder);
2934 int N = polyorder+1;
2937 cout<<npoints<<endl;
2938 for(
int u=0; u<npoints; u++)
2940 cout<<
"c="<<xin[u]<<
" "<<
2945 for(
int e=0; e<N; e++)
2948 for(
int row=0; row<N; row++)
2950 for(
int w=0; w < npoints; w++)
2952 A[N*e+row] += std::pow( xin[w], e+row);
2957 for(
int row= 0; row< N; row++)
2959 for(
int h=0; h< npoints; h++)
2961 b[row] += fin[h]*std::pow(xin[h],row);
2975 Lapack::Dgetrf( N, N, A.get(), N, ipivot.get(), info);
2979 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
2980 "th parameter had an illegal parameter for dgetrf";
2985 std::string message =
"ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
2986 boost::lexical_cast<std::string>(info) +
" is 0 from dgetrf";
2991 Lapack::Dgetrs(
'N', N, ncolumns_b , A.get() , N, ipivot.get(), b.get(), N, info);
2994 std::string message =
"ERROR: The " + boost::lexical_cast<std::string>(-info) +
2995 "th parameter had an illegal parameter for dgetrf";
3000 std::string message =
"ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
3001 boost::lexical_cast<std::string>(info) +
" is 0 from dgetrf";
3010 for(
int j=0; j<N; j++)
3016 for(
int h=0; h<N; h++)
3018 cout<<
"coeff:"<<b[h]<<endl;
3023 for(
int c=0; c< npout; c++)
3028 for(
int d=0; d< N; d++)
3032 *std::pow(xout[c],polorder);
3052 Vmath::Vcopy(inarray_x.num_elements() , inarray_x,1,tmpx,1);
3053 Vmath::Vcopy(inarray_x.num_elements() , inarray_y,1,tmpy,1);
3058 for(
int w=0; w<tmpx.num_elements(); w++)
3061 outarray_x[w]= tmpx[index];
3062 outarray_y[w]= tmpy[index];
3063 if(w< tmpx.num_elements()-1)
3065 if(tmpx[index] == tmpx[index+1])
3067 outarray_x[w+1]= tmpx[index+1];
3068 outarray_y[w+1]= tmpy[index+1];
3069 tmpx[index+1] = max+1000;
3085 tmpx[index] = max+1000;
3097 int np_lay = layers_y[0].num_elements();
3099 for(
int h=1; h<nlays-1; h++)
3102 for(
int s=0; s<nvertl; s++)
3105 ASSERTL0(ynew[ lay_Vids[h][s] ]==-20,
"ynew layers not empty");
3109 ynew[ lay_Vids[h][s] ] = ynew[Down[s]]+ h*abs(ynew[Down[s]] - yc[s])/(cntlow+1);
3111 xnew[lay_Vids[h][s] ] = xc[s];
3115 layers_y[h][s] = ynew[ lay_Vids[h][s] ];
3116 layers_x[h][s] = xnew[ lay_Vids[h][s] ];
3121 ynew[ lay_Vids[h][s] ] = yc[s] + (h-cntlow)*abs(ynew[Up[s]] - yc[s])/(cntup+1);
3123 xnew[lay_Vids[h][s] ] = xc[s];
3125 layers_y[h][s] = ynew[ lay_Vids[h][s] ];
3126 layers_x[h][s] = xnew[ lay_Vids[h][s] ];
3140 int np_lay = xcPhys.num_elements();
3141 int nedges = nvertl-1;
3148 int closepoints = 4;
3153 for(
int g=0; g< nedges; g++)
3162 xnew[Vids[g] ]= xcPhys[g*npedge+0];
3163 ylay[g*npedge +0] = ynew[ Vids[g] ];
3164 xlay[g*npedge +0] = xnew[ Vids[g] ];
3172 ynew[Vids[g+1] ]=
LagrangeInterpolant(xcPhys[g*npedge +npedge-1],closepoints,Pxinterp,Pyinterp );
3173 xnew[Vids[g+1] ]= xcPhys[g*npedge +npedge-1];
3174 ylay[g*npedge +npedge-1] = ynew[Vids[g+1] ];
3175 xlay[g*npedge +npedge-1] = xnew[Vids[g+1] ];
3180 for(
int r=0; r< npedge-2; r++)
3188 ASSERTL0( index<= tmpy.num_elements()-1,
" index wrong");
3192 ylay[g*npedge +r+1]=
3194 xcPhys[g*npedge +r+1],closepoints,Pxinterp,Pyinterp );
3196 xlay[g*npedge +r+1]= xcPhys[g*npedge +r+1];
3219 int np_lay = xcPhys.num_elements();
3220 int nedges = nvertl-1;
3228 int closepoints = 4;
3233 for(
int g=0; g< nedges; g++)
3237 ynew[Vids[g] ]= tmpy_lay[g*npedge+0];
3238 xnew[Vids[g] ]= tmpx_lay[g*npedge+0];
3241 ylay[g*npedge +0] = ynew[ Vids[g] ];
3242 xlay[g*npedge +0] = xnew[ Vids[g] ];
3245 ynew[Vids[g+1] ]= tmpy_lay[g*npedge+npedge-1];
3246 xnew[Vids[g+1] ]= tmpx_lay[g*npedge+npedge-1];
3247 ylay[g*npedge +npedge-1] = ynew[Vids[g+1] ];
3248 xlay[g*npedge +npedge-1] = xnew[Vids[g+1] ];
3253 for(
int r=0; r< npedge-2; r++)
3255 x0 = xlay[g*npedge +0];
3256 x1 = xlay[g*npedge +npedge-1];
3257 xtmp = x0 + r*(x1-x0)/(npedge-1);
3263 ASSERTL0( index<= tmpy.num_elements()-1,
" index wrong");
3267 ylay[g*npedge +r+1]=
3269 xtmp,closepoints,Pxinterp,Pyinterp );
3271 xlay[g*npedge +r+1]= xtmp;
3287 int nvertl = ycold.num_elements();
3288 int nVertTot = mesh->GetNvertices();
3289 for(
int n=0; n<nVertTot; n++)
3294 vertex->GetCoords(x,y,z);
3299 for(
int k=0; k<nvertl; k++)
3301 if(abs(x-xcold[k]) < tmp)
3303 tmp = abs(x-xcold[k]);
3316 nplay_closer= (qp_closer-1)*npedge +npedge-1;
3320 if( y>yoldup[qp_closer] && y<1 )
3325 ratio = (1-ylayup[nplay_closer])/
3326 ( (1-yoldup[qp_closer]) );
3328 ynew[n] = ylayup[nplay_closer]
3329 + (y-yoldup[qp_closer])*ratio;
3333 else if( y< yolddown[qp_closer] && y>-1 )
3336 ratio = (1+ylaydown[nplay_closer])/
3337 ( (1+yolddown[qp_closer]) );
3339 ynew[n] = ylaydown[nplay_closer]
3340 + (y-yolddown[qp_closer])*ratio;
3371 int nvertl = xoldup.num_elements();
3372 int nedges = nvertl-1;
3381 for(
int a=0; a< nedges;a++)
3386 xnew_down[a] = xlaydown[a*npedge+0];
3387 ynew_down[a] = ylaydown[a*npedge+0];
3388 xnew_up[a] = xlayup[a*npedge+0];
3389 ynew_up[a] = ylayup[a*npedge+0];
3390 nxvert[a] = nxPhys[a*npedge+0];
3391 nyvert[a] = nyPhys[a*npedge+0];
3393 xnew_down[a+1] = xlaydown[a*npedge+npedge-1];
3394 ynew_down[a+1] = ylaydown[a*npedge+npedge-1];
3395 xnew_up[a+1] = xlayup[a*npedge+npedge-1];
3396 ynew_up[a+1] = ylayup[a*npedge+npedge-1];
3397 nxvert[a+1] = nxPhys[a*npedge+npedge-1];
3398 nyvert[a+1] = nyPhys[a*npedge+npedge-1];
3403 xnew_down[a+1] = xlaydown[a*npedge+npedge-1];
3404 ynew_down[a+1] = ylaydown[a*npedge+npedge-1];
3405 xnew_up[a+1] = xlayup[a*npedge+npedge-1];
3406 ynew_up[a+1] = ylayup[a*npedge+npedge-1];
3407 nxvert[a+1] = nxPhys[a*npedge+npedge-1];
3408 nyvert[a+1] = nyPhys[a*npedge+npedge-1];
3420 int nVertTot = mesh->GetNvertices();
3421 for(
int n=0; n<nVertTot; n++)
3426 vertex->GetCoords(x,y,z);
3427 int qp_closeroldup = 0, qp_closerolddown = 0;
3436 for(
int k=0; k<nvertl; k++)
3438 if(abs(x-xolddown[k]) < diffdown)
3440 diffdown = abs(x-xolddown[k]);
3443 if(abs(x-xoldup[k]) < diffup)
3445 diffup = abs(x-xoldup[k]);
3455 int qp_closerup = 0, qp_closerdown = 0;
3457 for(
int f=0; f< nvertl; f++)
3459 if(abs(x-xnew_down[f]) < diffdown)
3461 diffdown = abs(x-xnew_down[f]);
3464 if(abs(x-xnew_up[f]) < diffup)
3466 diffup = abs(x-xnew_up[f]);
3493 int qp_closernormup;
3504 int qp_closernormdown;
3515 if( y>yoldup[qp_closeroldup] && y<1 )
3520 ratio = (1-ynew_up[qp_closerup])/
3521 ( (1-yoldup[qp_closeroldup]) );
3526 ynew[n] = ynew_up[qp_closerup]
3527 + (y-yoldup[qp_closeroldup])*ratio;
3533 if(x> (xmax-xmin)/2. && x< xmax)
3535 ratiox = (xmax-xnew_up[qp_closernormup])/
3536 (xmax-xoldup[qp_closernormup]) ;
3537 if( (xmax-xoldup[qp_closernormup])==0)
3544 xnew[n] = x + abs(nxvert[qp_closernormup])*(xnew_up[qp_closeroldup]-xoldup[qp_closeroldup])*ratiox;
3545 ASSERTL0(x>xmin,
" x value <xmin up second half");
3546 ASSERTL0(x<xmax," x value >xmax up second half
");
3548 else if( x> xmin && x<= (xmax-xmin)/2.)
3550 //cout<<"up close normold=
"<<qp_closernormoldup<<" closenorm=
"<<qp_closernormup<<endl;
3551 ratiox = (xnew_up[qp_closernormup]-xmin)/
3552 ( (xoldup[qp_closernormup]-xmin) );
3553 if( (xoldup[qp_closernormup]-xmin)==0)
3557 //xnew[n] = xnew_up[qp_closerup]
3558 // + (x-xoldup[qp_closeroldup])*ratiox;
3559 xnew[n] = x + abs(nxvert[qp_closernormup])*(xnew_up[qp_closeroldup]-xoldup[qp_closeroldup])*ratiox;
3560 //cout<<"up xold=
"<<x<<" xnew=
"<<xnew[n]<<endl;
3561 ASSERTL0(x>xmin," x value <xmin up first half
");
3562 ASSERTL0(x<xmax," x value >xmax up first half
");
3567 else if( y< yolddown[qp_closerolddown] && y>-1 )
3570 ratio = (1+ynew_down[qp_closerdown])/
3571 ( (1+yolddown[qp_closerolddown]) );
3573 // ratioy = (1-ynew_down[qp_closernormdown])/
3574 // ( (1-yolddown[qp_closernormolddown]) );
3576 //distance prop to layerlow
3577 ynew[n] = ynew_down[qp_closerdown]
3578 + (y-yolddown[qp_closerolddown])*ratio;
3579 //ynew[n] = y +abs(nyvert[qp_closernormdown])*
3580 // (ynew_down[qp_closerolddown]-yolddown[qp_closerolddown])*ratioy;
3581 //ynew[n] = y + 0.3*(ynew_down[qp_closerdown]-yolddown[qp_closerdown]);
3582 //xnew[n] = x + abs(nxvert[qp_closerolddown])*(xnew_down[qp_closerolddown]-xolddown[qp_closerolddown]);
3586 cout<<qp_closerolddown<<" nplaydown=
"<<qp_closerdown<<endl;
3587 cout<<"xolddown=
"<<xolddown[qp_closerolddown]<<" xnewdown=
"<<xnew_down[qp_closerdown]<<endl;
3588 cout<<"xold+
"<<x<<" xnew+
"<<xnew[n]<<endl;
3592 if(x> (xmax-xmin)/2. && x <xmax)
3594 ratiox = (xmax-xnew_down[qp_closernormdown])/
3595 ( (xmax-xolddown[qp_closernormdown]) );
3596 if( (xmax-xolddown[qp_closernormdown])==0)
3600 //xnew[n] = xnew_down[qp_closerdown]
3601 // + (x-xolddown[qp_closerolddown])*ratiox;
3603 abs(nxvert[qp_closernormdown])*(xnew_down[qp_closerolddown]-xolddown[qp_closerolddown])*ratiox;
3604 ASSERTL0(x>xmin," x value <xmin down second half
");
3605 ASSERTL0(x<xmax," x value >xmax down second half
");
3607 else if( x>xmin && x<= (xmax-xmin)/2.)
3609 ratiox = (xnew_down[qp_closernormdown]-xmin)/
3610 ( (xolddown[qp_closernormdown]-xmin) );
3611 if( (xolddown[qp_closernormdown]-xmin)==0)
3615 //xnew[n] = xnew_down[qp_closerdown]
3616 // + (x-xolddown[qp_closerolddown])*ratiox;
3618 abs(nxvert[qp_closernormdown])*(xnew_down[qp_closerolddown]-xolddown[qp_closerolddown])*ratiox;
3619 ASSERTL0(x>xmin," x value <xmin down first half
");
3620 ASSERTL0(x<xmax," x value >xmax down first half
");
3625 cout<<"xold
"<<x<<" xnew=
"<<xnew[n]<<endl;
3626 ASSERTL0(xnew[n] >= xmin, "newx < xmin
");
3627 ASSERTL0(xnew[n]<= xmax, "newx > xmax
");
3632 void CheckSingularQuads( MultiRegions::ExpListSharedPtr Exp,
3633 Array<OneD, int> V1, Array<OneD, int> V2,
3634 Array<OneD, NekDouble>& xnew,Array<OneD, NekDouble>& ynew)
3636 const boost::shared_ptr<LocalRegions::ExpansionVector> exp2D = Exp->GetExp();
3637 int nel = exp2D->size();
3638 LocalRegions::QuadExpSharedPtr locQuadExp;
3639 LocalRegions::TriExpSharedPtr locTriExp;
3640 SpatialDomains::Geometry1DSharedPtr SegGeom;
3642 NekDouble xV1, yV1, xV2,yV2;
3643 NekDouble slopebef,slopenext,slopenew;
3644 Array<OneD, int> locEids(4);
3645 for(int i=0; i<nel; i++)
3647 if((locQuadExp = (*exp2D)[i]->as<LocalRegions::QuadExp>()))
3649 SegGeom = (locQuadExp->GetGeom2D())->GetEdge(0);
3650 idbef = SegGeom->GetEid();
3651 if(xnew[ V1[idbef] ]<= xnew[ V2[idbef] ])
3653 xV1 = xnew[ V1[idbef] ];
3654 yV1 = ynew[ V1[idbef] ];
3655 xV2 = xnew[ V2[idbef] ];
3656 yV2 = ynew[ V2[idbef] ];
3657 slopebef = (yV2 -yV1)/(xV2 -xV1);
3661 xV1 = xnew[ V2[idbef] ];
3662 yV1 = ynew[ V2[idbef] ];
3663 xV2 = xnew[ V1[idbef] ];
3664 yV2 = ynew[ V1[idbef] ];
3665 slopebef = (yV2 -yV1)/(xV2 -xV1);
3667 //cout<<"00 V1 x=
"<<xnew[ V1[idbef] ]<<" y=
"<<ynew[ V1[idbef] ]<<endl;
3668 //cout<<"00 V2 x=
"<<xnew[ V2[idbef] ]<<" y=
"<<ynew[ V2[idbef] ]<<endl;
3669 for(int j = 1; j < locQuadExp->GetNedges(); ++j)
3671 SegGeom = (locQuadExp->GetGeom2D())->GetEdge(j);
3672 idnext = SegGeom->GetEid();
3673 //cout<<"id=
"<<idnext<<" locid=
"<<j<<endl;
3674 //cout<<" V1 x=
"<<xnew[ V1[idnext] ]<<" y=
"<<ynew[ V1[idnext] ]<<endl;
3675 //cout<<" V2 x=
"<<xnew[ V2[idnext] ]<<" y=
"<<ynew[ V2[idnext] ]<<endl;
3676 if(xV1 == xnew[ V1[idnext] ] && yV1 == ynew[ V1[idnext] ] )
3678 xV1 = xnew[ V1[idnext] ];
3679 yV1 = ynew[ V1[idnext] ];
3680 xV2 = xnew[ V2[idnext] ];
3681 yV2 = ynew[ V2[idnext] ];
3682 slopenext = (yV2 -yV1)/(xV2 -xV1);
3685 cout<<"case1 x0=
"<<xV1<<" x1=
"<<xV2<<endl;
3686 cout<<idnext<<" 11slope bef =
"<<slopebef<<" slopenext=
"<<slopenext<<endl;
3688 //compare with slope before
3689 if( slopebef/slopenext>0.84 && slopebef/slopenext <1.18)
3691 xnew[ V1[idnext] ] = xnew[ V1[idnext] ] -0.01;
3692 slopenew = (yV2-yV1)/(xV2- xnew[ V1[idnext] ]);
3694 if( abs(slopebef-slopenew) < abs(slopebef-slopenext) )
3696 xnew[ V1[idnext] ] = xnew[ V1[idnext] ] +0.02;
3697 slopenew = (yV2-yV1)/(xV2- xnew[ V1[idnext] ]);
3699 slopenext = slopenew;
3700 cout<<"slopenew=
"<<slopenew<<endl;
3701 cout<<"moved x=
"<<xnew[ V1[idnext] ]<<endl;
3704 else if(xV2 == xnew[ V2[idnext] ] && yV2 == ynew[ V2[idnext] ] )
3706 xV1 = xnew[ V2[idnext] ];
3707 yV1 = ynew[ V2[idnext] ];
3708 xV2 = xnew[ V1[idnext] ];
3709 yV2 = ynew[ V1[idnext] ];
3710 slopenext = (yV2 -yV1)/(xV2 -xV1);
3713 cout<<"case2 x0=
"<<xV1<<" x1=
"<<xV2<<endl;
3714 cout<<idnext<<" 22slope bef =
"<<slopebef<<" slopenext=
"<<slopenext<<endl;
3716 //compare with slope before
3717 if( slopebef/slopenext>0.84 && slopebef/slopenext <1.18)
3719 xnew[ V2[idnext] ] = xnew[ V2[idnext] ] -0.01;
3720 slopenew = (yV2-yV1)/(xV2- xnew[ V2[idnext] ]);
3722 if( abs(slopebef-slopenew) < abs(slopebef-slopenext) )
3724 xnew[ V2[idnext] ] = xnew[ V2[idnext] ] +0.02;
3725 slopenew = (yV2-yV1)/(xV2- xnew[ V2[idnext] ]);
3728 slopenext = slopenew;
3729 cout<<"slopenew=
"<<slopenew<<endl;
3730 cout<<"moved x=
"<<xnew[ V2[idnext] ]<<endl;
3733 else if(xV1 == xnew[ V2[idnext] ] && yV1 == ynew[ V2[idnext] ] )
3735 xV1 = xnew[ V2[idnext] ];
3736 yV1 = ynew[ V2[idnext] ];
3737 xV2 = xnew[ V1[idnext] ];
3738 yV2 = ynew[ V1[idnext] ];
3739 slopenext = (yV2 -yV1)/(xV2 -xV1);
3742 cout<<"case3 x0=
"<<xV1<<" x1=
"<<xV2<<endl;
3743 cout<<idnext<<" 22slope bef =
"<<slopebef<<" slopenext=
"<<slopenext<<endl;
3745 //compare with slope before
3746 if( slopebef/slopenext>0.84 && slopebef/slopenext <1.18)
3748 xnew[ V2[idnext] ] = xnew[ V2[idnext] ] -0.01;
3749 slopenew = (yV2-yV1)/(xV2- xnew[ V2[idnext] ]);
3751 if( abs(slopebef-slopenew) < abs(slopebef-slopenext) )
3753 xnew[ V2[idnext] ] = xnew[ V2[idnext] ] +0.02;
3754 slopenew = (yV2-yV1)/(xV2- xnew[ V2[idnext] ]);
3756 slopenext = slopenew;
3757 cout<<"slopenew=
"<<slopenew<<endl;
3758 cout<<"moved x=
"<<xnew[ V2[idnext] ]<<endl;
3762 else if(xV2 == xnew[ V1[idnext] ] && yV2 == ynew[ V1[idnext] ] )
3764 xV1 = xnew[ V1[idnext] ];
3765 yV1 = ynew[ V1[idnext] ];
3766 xV2 = xnew[ V2[idnext] ];
3767 yV2 = ynew[ V2[idnext] ];
3768 slopenext = (yV2 -yV1)/(xV2 -xV1);
3771 cout<<"case4 x0=
"<<xV1<<" x1=
"<<xV2<<endl;
3772 cout<<idnext<<" 22slope bef =
"<<slopebef<<" slopenext=
"<<slopenext<<endl;
3774 //compare with slope before
3775 if( slopebef/slopenext>0.84 && slopebef/slopenext <1.18)
3777 xnew[ V1[idnext] ] = xnew[ V1[idnext] ] -0.01;
3778 slopenew = (yV2-yV1)/(xV2- xnew[ V1[idnext] ]);
3780 if( abs(slopebef-slopenew) < abs(slopebef-slopenext) )
3782 xnew[ V1[idnext] ] = xnew[ V1[idnext] ] +0.02;
3783 slopenew = (yV2-yV1)/(xV2- xnew[ V1[idnext] ]);
3785 slopenext = slopenew;
3786 cout<<"slopenew=
"<<slopenew<<endl;
3787 cout<<"moved x=
"<<xnew[ V1[idnext] ]<<endl;
3793 ASSERTL0(false, "edge not connected
");
3795 slopebef = slopenext;
3804 void Replacevertices(string filename, Array<OneD, NekDouble> newx,
3805 Array<OneD, NekDouble> newy,
3806 Array<OneD, NekDouble> xcPhys, Array<OneD, NekDouble> ycPhys,
3807 Array<OneD, int>Eids, int Npoints, string s_alp,
3808 Array<OneD, Array<OneD, NekDouble> > x_lay,
3809 Array<OneD, Array<OneD, NekDouble> > y_lay,
3810 Array<OneD, Array<OneD, int > >lay_eids, bool curv_lay)
3812 //load existing file
3814 TiXmlDocument doc(filename);
3815 //load xscale parameter (if exists)
3816 TiXmlElement* master = doc.FirstChildElement("NEKTAR
");
3817 TiXmlElement* mesh = master->FirstChildElement("GEOMETRY
");
3818 TiXmlElement* element = mesh->FirstChildElement("VERTEX
");
3819 NekDouble xscale = 1.0;
3820 LibUtilities::AnalyticExpressionEvaluator expEvaluator;
3821 const char *xscal = element->Attribute("XSCALE
");
3824 std::string xscalstr = xscal;
3825 int expr_id = expEvaluator.DefineFunction("",xscalstr);
3826 xscale = expEvaluator.Evaluate(expr_id);
3830 // Save a new XML file.
3831 newfile = filename.substr(0, filename.find_last_of(".
"))+"_moved.xml
";
3832 doc.SaveFile( newfile );
3834 //write the new vertices
3835 TiXmlDocument docnew(newfile);
3836 bool loadOkaynew = docnew.LoadFile();
3838 std::string errstr = "Unable to load file:
";
3840 ASSERTL0(loadOkaynew, errstr.c_str());
3842 TiXmlHandle docHandlenew(&docnew);
3843 TiXmlElement* meshnew = NULL;
3844 TiXmlElement* masternew = NULL;
3845 TiXmlElement* condnew = NULL;
3846 TiXmlElement* Parsnew = NULL;
3847 TiXmlElement* parnew = NULL;
3849 // Master tag within which all data is contained.
3852 masternew = docnew.FirstChildElement("NEKTAR
");
3853 ASSERTL0(masternew, "Unable to
find NEKTAR tag in file.
");
3855 //set the alpha value
3857 condnew = masternew->FirstChildElement("CONDITIONS
");
3858 Parsnew = condnew->FirstChildElement("PARAMETERS
");
3859 cout<<"alpha=
"<<s_alp<<endl;
3860 parnew = Parsnew->FirstChildElement("P");
3863 TiXmlNode *node = parnew->FirstChild();
3866 // Format is "paramName = value
"
3867 std::string line = node->ToText()->Value();
3871 int beg = line.find_first_not_of("
");
3872 int end = line.find_first_of("=
");
3873 // Check for no parameter name
3874 if (beg == end) throw 1;
3875 // Check for no parameter value
3876 if (end != line.find_last_of("=
")) throw 1;
3877 // Check for no equals sign
3878 if (end == std::string::npos) throw 1;
3879 lhs = line.substr(line.find_first_not_of(" "), end-beg);
3880 lhs = lhs.substr(0, lhs.find_last_not_of(" ")+1);
3882 //rhs = line.substr(line.find_last_of("=
")+1);
3883 //rhs = rhs.substr(rhs.find_first_not_of(" "));
3884 //rhs = rhs.substr(0, rhs.find_last_not_of(" ")+1);
3886 boost::to_upper(lhs);
3889 alphastring = "Alpha =
"+ s_alp;
3890 parnew->RemoveChild(node);
3891 parnew->LinkEndChild(new TiXmlText(alphastring) );
3895 parnew = parnew->NextSiblingElement("P");
3899 // Find the Mesh tag and same the dim and space attributes
3900 meshnew = masternew->FirstChildElement("GEOMETRY
");
3902 ASSERTL0(meshnew, "Unable to
find GEOMETRY tag in file.
");
3903 // Now read the vertices
3904 TiXmlElement* elementnew = meshnew->FirstChildElement("VERTEX
");
3905 ASSERTL0(elementnew, "Unable to
find mesh VERTEX tag in file.
");
3909 elementnew->SetAttribute("XSCALE
",1.0);
3911 TiXmlElement *vertexnew = elementnew->FirstChildElement("V
");
3917 int nextVertexNumber = -1;
3922 //delete the old one
3923 TiXmlAttribute *vertexAttr = vertexnew->FirstAttribute();
3924 std::string attrName(vertexAttr->Name());
3925 ASSERTL0(attrName == "ID
", (std::string("Unknown attribute name:
") + attrName).c_str());
3927 err = vertexAttr->QueryIntValue(&indx);
3928 ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.
");
3929 ASSERTL0(indx == nextVertexNumber, "Vertex IDs must begin with zero and be sequential.
");
3931 std::string vertexBodyStr;
3932 // Now read body of vertex
3933 TiXmlNode *vertexBody = vertexnew->FirstChild();
3934 // Accumulate all non-comment body data.
3935 if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
3937 vertexBodyStr += vertexBody->ToText()->Value();
3938 vertexBodyStr += " ";
3940 ASSERTL0(!vertexBodyStr.empty(), "Vertex definitions must contain vertex data.
");
3941 //remove the old coordinates
3942 vertexnew->RemoveChild(vertexBody);
3944 //cout<<"writing.. v:
"<<nextVertexNumber<<endl;
3946 //we need at least 5 digits (setprecision 5) to get the streak position with
3948 s << std::scientific << std::setprecision(8) << newx[nextVertexNumber] << " "
3949 << newy[nextVertexNumber] << " " << 0.0;
3950 vertexnew->LinkEndChild(new TiXmlText(s.str()));
3951 //TiXmlNode *newvertexBody = vertexnew->FirstChild();
3952 //string newvertexbodystr= newvertexBody->SetValue(s.str());
3953 //vertexnew->ReplaceChild(vertexBody,new TiXmlText(newvertexbodystr));
3955 vertexnew = vertexnew->NextSiblingElement("V
");
3960 //read the curved tag
3961 TiXmlElement* curvednew = meshnew->FirstChildElement("CURVED
");
3962 ASSERTL0(curvednew, "Unable to
find mesh CURVED tag in file.
");
3963 TiXmlElement *edgenew = curvednew->FirstChildElement("E
");
3965 //ID is different from index...
3966 std::string charindex;
3970 int neids_lay = lay_eids[0].num_elements();
3971 //if edgenew belongs to the crit lay replace it, else delete it.
3977 TiXmlAttribute *edgeAttr = edgenew->FirstAttribute();
3978 std::string attrName(edgeAttr->Name());
3979 charindex = edgeAttr->Value();
3980 std::istringstream iss(charindex);
3981 iss >> std::dec >> index;
3983 edgenew->QueryIntAttribute("EDGEID
", &eid);
3984 //cout<<"eid=
"<<eid<<" neid=
"<<Eids.num_elements()<<endl;
3985 //find the corresponding index curve point
3986 for(int u=0; u<Eids.num_elements(); u++)
3988 //cout<<"Eids=
"<<Eids[u]<<" eid=
"<<eid<<endl;
3996 curvednew->RemoveChild(edgenew);
3997 //ASSERTL0(false, "edge to update not found
");
4002 std::string edgeBodyStr;
4003 //read the body of the edge
4004 TiXmlNode *edgeBody = edgenew->FirstChild();
4005 if(edgeBody->Type() == TiXmlNode::TINYXML_TEXT)
4007 edgeBodyStr += edgeBody->ToText()->Value();
4010 ASSERTL0(!edgeBodyStr.empty(), "Edge definitions must contain edge data
");
4011 //remove the old coordinates
4012 edgenew->RemoveChild(edgeBody);
4013 //write the new points coordinates
4014 //we need at least 5 digits (setprecision 5) to get the streak position with
4017 //Determine the number of points
4018 err = edgenew->QueryIntAttribute("NUMPOINTS
", &numPts);
4019 ASSERTL0(err == TIXML_SUCCESS, "Unable to read curve attribute NUMPOINTS.
");
4022 edgenew->SetAttribute("NUMPOINTS
", Npoints);
4023 for(int u=0; u< Npoints; u++)
4025 st << std::scientific <<
4026 std::setprecision(8) <<xcPhys[cnt*Npoints+u]
4027 << " " << ycPhys[cnt*Npoints+u] << " " << 0.000<<" ";
4030 edgenew->LinkEndChild(new TiXmlText(st.str()));
4034 st << std::scientific << std::setprecision(8) << x_crit[v1] << " "
4035 << y_crit[v1] << " " << 0.000<<" ";
4036 for(int a=0; a< Npoints-2; a++)
4038 st << std::scientific << std::setprecision(8) <<
4039 " "<<Pcurvx[indexeid*(Npoints-2) +a]<<" "<<Pcurvy[indexeid*(Npoints-2) +a]
4043 st << std::scientific << std::setprecision(8) <<
4044 " "<<x_crit[v2]<<" "<< y_crit[v2] <<" "<< 0.000;
4045 edgenew->LinkEndChild(new TiXmlText(st.str()));
4052 edgenew = edgenew->NextSiblingElement("E
");
4056 //write also the others layers curve points
4057 if(curv_lay == true)
4059 cout<<"write other curved edges
"<<endl;
4060 TiXmlElement * curved = meshnew->FirstChildElement("CURVED
");
4062 int nlays = lay_eids.num_elements();
4064 //TiXmlComment * comment = new TiXmlComment();
4065 //comment->SetValue(" new edges
");
4066 //curved->LinkEndChild(comment);
4067 for (int g=0; g< nlays; ++g)
4069 for(int p=0; p< neids_lay; p++)
4072 TiXmlElement * e = new TiXmlElement( "E
" );
4073 e->SetAttribute("ID
", idcnt++);
4074 e->SetAttribute("EDGEID
", lay_eids[g][p]);
4075 e->SetAttribute("NUMPOINTS
", Npoints);
4076 e->SetAttribute("TYPE
", "PolyEvenlySpaced
");
4077 for(int c=0; c< Npoints; c++)
4079 st << std::scientific << std::setprecision(8) <<x_lay[g][p*Npoints +c]
4080 << " " << y_lay[g][p*Npoints +c] << " " << 0.000<<" ";
4084 TiXmlText * t0 = new TiXmlText(st.str());
4085 e->LinkEndChild(t0);
4086 curved->LinkEndChild(e);
4094 docnew.SaveFile( newfile );
4096 cout<<"new file:
"<<newfile<<endl;
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)
bool checkcommonvert(Array< OneD, int > Vids_laybefore, Array< OneD, int > Vids_c, int Vid)
#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)
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...
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)
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
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 Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
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 Cutrepetitions(int nedges, Array< OneD, NekDouble > inarray, Array< OneD, NekDouble > &outarray)
int main(int argc, char *argv[])
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.
void MoveOutsidePointsfixedxpos(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 > ylaydown, Array< OneD, NekDouble > ylayup, Array< OneD, NekDouble > &xnew, Array< OneD, NekDouble > &ynew)
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)
boost::shared_ptr< QuadExp > QuadExpSharedPtr
Represents a vertex in the mesh.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
void MoveLayerNfixedxpos(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 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
boost::shared_ptr< Geometry1D > Geometry1DSharedPtr
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
int DetermineclosePointxindex(NekDouble x, Array< OneD, NekDouble > xArray)
const BoundaryRegionCollection & GetBoundaryRegions(void) const
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 PolyInterp(Array< OneD, NekDouble > xpol, Array< OneD, NekDouble > ypol, Array< OneD, NekDouble > &coeffsinterp, Array< OneD, NekDouble > &xcout, Array< OneD, NekDouble > &ycout, int edge, int npedge)
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 Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
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)
boost::shared_ptr< TriExp > TriExpSharedPtr