49#include <boost/lexical_cast.hpp>
110 int edge,
int npedge);
170 int Npoints,
string s_alp,
175int main(
int argc,
char *argv[])
181 if (argc > 6 || argc < 5)
183 fprintf(stderr,
"Usage: ./MoveMesh meshfile fieldfile changefile "
184 "alpha cr(optional)\n");
191 LibUtilities::SessionReader::CreateInstance(2, argv);
193 SpatialDomains::MeshGraph::Read(vSession);
196 if (argc == 6 && vSession->DefinesSolverInfo(
"INTERFACE") &&
197 vSession->GetSolverInfo(
"INTERFACE") ==
"phase")
199 cr = boost::lexical_cast<NekDouble>(argv[argc - 1]);
207 vSession, graphShPt);
217 string changefile(argv[argc - 2]);
221 string charalp(argv[argc - 1]);
223 cout <<
"read alpha=" << charalp << endl;
227 string fieldfile(argv[argc - 3]);
228 vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
229 vector<vector<NekDouble>> fielddata;
234 vSession, graphShPt,
"w",
true);
238 for (
int i = 0; i < fielddata.size(); i++)
240 streak->ExtractDataToCoeffs(fielddef[i], fielddata[i],
241 fielddef[i]->m_fields[0],
242 streak->UpdateCoeffs());
244 streak->BwdTrans(streak->GetCoeffs(), streak->UpdatePhys());
250 int nIregions, lastIregion = 0;
252 streak->GetBndConditions();
256 int nbnd = bndConditions.size();
257 for (
int r = 0; r < nbnd; r++)
259 if (bndConditions[r]->GetUserDefined() ==
"CalcBC")
268 "there is any boundary region with the tag USERDEFINEDTYPE="
271 cout <<
"nIregions=" << nIregions << endl;
274 streak->GetBndCondExpansions();
277 int nedges = bndfieldx[lastIregion]->GetExpSize();
278 int nvertl = nedges + 1;
292 ((bndfieldx[lastIregion]->GetExp(0)->as<LocalRegions::SegExp>())
295 vertex0->GetCoords(x0, y0, z0);
298 cout <<
"WARNING x0=" << x0 << endl;
303 OrderVertices(nedges, graphShPt, bndfieldx[lastIregion - 1], Vids_low, v1,
304 v2, x_connect, lastedge, xold_low, yold_low);
305 ASSERTL0(Vids_low[v2] != -10,
"Vids_low[v2] is wrong");
307 graphShPt->GetVertex(Vids_low[v2]);
310 cout <<
"x_conn=" << x_connect <<
" yt=" << yt <<
" zt=" << zt
311 <<
" vid=" << Vids_low[v2] << endl;
312 vertex->GetCoords(x_connect, yt, zt);
318 OrderVertices(nedges, graphShPt, bndfieldx[lastIregion - 1], Vids_low,
319 v1, v2, x_connect, lastedge, xold_low, yold_low);
321 graphShPt->GetVertex(Vids_low[v1]);
323 vertex->GetCoords(x_connect, yt, zt);
335 vertex0 = graphShPt->GetVertex(
336 ((bndfieldx[lastIregion]->GetExp(0)->as<LocalRegions::SegExp>())
339 vertex0->GetCoords(x0, y0, z0);
342 cout <<
"WARNING x0=" << x0 << endl;
349 OrderVertices(nedges, graphShPt, bndfieldx[lastIregion - 2], Vids_up, v1,
350 v2, x_connect, lastedge, xold_up, yold_up);
352 graphShPt->GetVertex(Vids_up[v2]);
353 vertexU->GetCoords(x_connect, yt, zt);
359 OrderVertices(nedges, graphShPt, bndfieldx[lastIregion - 2], Vids_up,
360 v1, v2, x_connect, lastedge, xold_up, yold_up);
362 graphShPt->GetVertex(Vids_up[v1]);
365 vertex->GetCoords(x_connect, yt, zt);
377 vertex0 = graphShPt->GetVertex(
378 ((bndfieldx[lastIregion]->GetExp(0)->as<LocalRegions::SegExp>())
381 vertex0->GetCoords(x0, y0, z0);
384 cout <<
"WARNING x0=" << x0 << endl;
392 OrderVertices(nedges, graphShPt, bndfieldx[lastIregion], Vids_c, v1, v2,
393 x_connect, lastedge, xold_c, yold_c);
395 graphShPt->GetVertex(Vids_c[v2]);
398 vertexc->GetCoords(x_connect, yt, zt);
404 OrderVertices(nedges, graphShPt, bndfieldx[lastIregion], Vids_c, v1, v2,
405 x_connect, lastedge, xold_c, yold_c);
407 graphShPt->GetVertex(Vids_c[v1]);
410 vertex->GetCoords(x_connect, yt, zt);
419 for (
int r = 0; r < nvertl; r++)
424 Deltaup[r] = yold_up[r] - yold_c[r];
425 Deltalow[r] = yold_c[r] - yold_low[r];
427 "distance between upper and layer curve is not positive");
429 "distance between lower and layer curve is not positive");
440 vSession, *(bregions.find(lastIregion)->second), graphShPt,
true);
451 if (vSession->DefinesParameter(
"npedge"))
453 npedge = (int)vSession->GetParameter(
"npedge");
462 int nq = streak->GetTotPoints();
465 streak->GetCoords(x, y);
473 xold_c, yold_c, x_c, y_c, cr,
true);
476 for (
int q = 0;
q < nvertl;
q++)
478 if (y_c[
q] < yold_c[
q])
483 Delta_c[
q] =
abs(yold_c[
q] - y_c[
q]);
486 cout << x_c[
q] <<
" " << y_c[
q] << endl;
491 cout <<
"Warning: the critical layer is stationary" << endl;
512 for (
int r = 0; r < nedges; r++)
517 Eid = (bndSegExp->GetGeom1D())->GetGlobalID();
518 id1 = (bndSegExp->GetGeom1D())->GetVid(0);
519 id2 = (bndSegExp->GetGeom1D())->GetVid(1);
520 vertex1 = graphShPt->GetVertex(id1);
521 vertex2 = graphShPt->GetVertex(id2);
523 vertex2->GetCoords(x2, y2, z2);
526 cout <<
"edge=" << r <<
" x1=" << x1 <<
" y1=" << y1 <<
" x2=" << x2
527 <<
" y2=" << y2 << endl;
530 Cpointsx[r] = x1 + (x2 - x1) / 2;
534 if (Cpointsx[r] > x2 || Cpointsx[r] < x1)
536 Cpointsx[r] = -Cpointsx[r];
538 for (
int w = 0;
w < npedge - 2;
w++)
541 Addpointsx[r * (npedge - 2) +
w] =
542 x1 + ((x2 - x1) / (npedge - 1)) * (
w + 1);
543 if (Addpointsx[r * (npedge - 2) +
w] > x2 ||
544 Addpointsx[r * (npedge - 2) +
w] < x1)
546 Addpointsx[r * (npedge - 2) +
w] =
547 -Addpointsx[r * (npedge - 2) +
w];
550 Addpointsy[r * (npedge - 2) +
w] =
551 y_c[r] + ((y_c[r + 1] - y_c[r]) / (x_c[r + 1] - x_c[r])) *
552 (Addpointsx[r * (npedge - 2) +
w] - x1);
556 Addpointsy[r * (npedge - 2) +
w],
557 Addpointsx[r * (npedge - 2) +
w],
558 Addpointsy[r * (npedge - 2) +
w],
572 Cpointsx[r] = x2 + (x1 - x2) / 2;
574 if (Cpointsx[r] > x1 || Cpointsx[r] < x2)
576 Cpointsx[r] = -Cpointsx[r];
578 for (
int w = 0;
w < npedge - 2;
w++)
580 Addpointsx[r * (npedge - 2) +
w] =
581 x2 + ((x1 - x2) / (npedge - 1)) * (
w + 1);
582 if (Addpointsx[r * (npedge - 2) +
w] > x1 ||
583 Addpointsx[r * (npedge - 2) +
w] < x2)
585 Addpointsx[r * (npedge - 2) +
w] =
586 -Addpointsx[r * (npedge - 2) +
w];
590 Addpointsy[r * (npedge - 2) +
w] =
592 ((y_c[r] - y_c[r + 1]) / (x_c[r] - x_c[r + 1])) *
593 (Addpointsx[r * (npedge - 2) +
w] - x2);
598 Addpointsy[r * (npedge - 2) +
w],
599 Addpointsx[r * (npedge - 2) +
w],
600 Addpointsy[r * (npedge - 2) +
w],
613 ASSERTL0(
false,
"point not generated");
631 for (
int a = 0; a < nedges; a++)
634 xcPhys[a * npedge + 0] = x_c[a];
635 ycPhys[a * npedge + 0] = y_c[a];
637 xcPhys[a * npedge + npedge - 1] = x_c[a + 1];
638 ycPhys[a * npedge + npedge - 1] = y_c[a + 1];
640 for (
int b = 0; b < npedge - 2; b++)
642 xcPhys[a * npedge + b + 1] = Addpointsx[a * (npedge - 2) + b];
643 ycPhys[a * npedge + b + 1] = Addpointsy[a * (npedge - 2) + b];
647 cout <<
"xc,yc before tanevaluate" << endl;
648 for (
int v = 0; v < xcPhys.size(); v++)
650 cout << xcPhys[v] <<
" " << ycPhys[v] << endl;
663 MappingEVids(xold_up, yold_up, xold_low, yold_low, xold_c, yold_c, Vids_c,
664 graphShPt, streak, V1, V2, nlays, lay_Eids, lay_Vids);
666 cout <<
"nlays=" << nlays << endl;
670 for (
int g = 0; g < nlays; g++)
681 if (vSession->DefinesParameter(
"Delta"))
683 Delta0 = vSession->GetParameter(
"Delta");
695 int nVertTot = graphShPt->GetNvertices();
711 for (
int i = 0; i < nVertTot; i++)
713 bool mvpoint =
false;
716 vertex->GetCoords(x, y,
z);
725 if (x == 0 && y < yold_low[0] && y > bleft)
730 if (x == xold_c[nvertl - 1] && y > yold_up[nvertl - 1] && y < tright)
735 if (x == xold_c[nvertl - 1] && y < yold_low[nvertl - 1] && y > bright)
740 if (x == 0 && y > yold_up[0] && y < tleft)
745 for (
int j = 0; j < nvertl; j++)
747 if ((xold_up[j] == x) && (yold_up[j] == y))
751 ynew[i] = y_c[j] + Delta0;
755 if ((xold_low[j] == x) && (yold_low[j] == y))
759 ynew[i] = y_c[j] - Delta0;
763 if ((xold_c[j] == x) && (yold_c[j] == y))
769 if (mvpoint ==
false)
774 for (
int k = 0; k < nvertl; k++)
776 if (
abs(x - xold_up[k]) < diff)
778 diff =
abs(x - xold_up[k]);
782 if (y > yold_up[qp_closer] && y < 1)
792 ynew[i] = y_c[qp_closer] + (y - yold_c[qp_closer]) *
793 (1 - y_c[qp_closer]) /
794 (1 - yold_c[qp_closer]);
797 else if (y < yold_low[qp_closer] && y > -1)
807 ynew[i] = y_c[qp_closer] + (y - yold_c[qp_closer]) *
808 (-1 - y_c[qp_closer]) /
809 (-1 - yold_c[qp_closer]);
812 else if (y > yold_c[qp_closer] && y < yold_up[qp_closer])
822 else if (y < yold_c[qp_closer] && y > yold_low[qp_closer])
832 else if (y == 1 || y == -1)
840 if ((ynew[i] > 1 || ynew[i] < -1) &&
841 (y > yold_up[qp_closer] || y < yold_low[qp_closer]))
843 cout <<
"point x=" << xnew[i] <<
" y=" << y
844 <<
" closer x=" << xold_up[qp_closer]
845 <<
" ynew=" << ynew[i] << endl;
846 ASSERTL0(
false,
"shifting out of range");
851 int nqedge = streak->GetExp(0)->GetNumPoints(0);
852 int nquad_lay = (nvertl - 1) * nqedge;
855 bool curv_lay =
true;
856 bool move_norm =
true;
857 int np_lay = (nvertl - 1) * npedge;
867 if (move_norm ==
true)
879 cout <<
"nquad per edge=" << nqedge << endl;
880 for (
int l = 0; l < 2; l++)
882 Edge_newcoords[l] = bndfieldx[lastIregion]
906 bndfieldx[lastIregion]->GetCoords(xcQ, ycQ, zcQ);
915 for (
int l = 0; l < xcQ.size(); l++)
929 Vmath::Vcopy(nquad_lay, ycQ, 1, Cont_y->UpdatePhys(), 1);
932 Cont_y->FwdTransLocalElmt(Cont_y->GetPhys(), coeffsy);
933 Cont_y->BwdTrans(coeffsy, Cont_y->UpdatePhys());
936 cout <<
"xcQ, ycQ" << endl;
937 for (
int s = 0; s < xcQ.size(); s++)
939 cout << xcQ[s] <<
" " << ycQ[s] << endl;
942 bool evaluatetan =
false;
946 for (
int k = 0; k < nedges; k++)
948 Vmath::Vcopy(nqedge, &xcQ[k * nqedge], 1, &xcedgeQ[0], 1);
949 Vmath::Vcopy(nqedge, &ycQ[k * nqedge], 1, &ycedgeQ[0], 1);
951 Edge_newcoords[0]->StdPhysDeriv(xcedgeQ, txedgeQ);
952 Edge_newcoords[1]->StdPhysDeriv(ycedgeQ, tyedgeQ);
954 Vmath::Vmul(nqedge, txedgeQ, 1, txedgeQ, 1, normsQ, 1);
956 Vmath::Vvtvp(nqedge, tyedgeQ, 1, tyedgeQ, 1, normsQ, 1, normsQ, 1);
961 Vmath::Vmul(nqedge, txedgeQ, 1, normsQ, 1, txedgeQ, 1);
962 Vmath::Vmul(nqedge, tyedgeQ, 1, normsQ, 1, tyedgeQ, 1);
967 for (
int u = 0; u < nqedge - 1; u++)
969 incratio = (ycedgeQ[u + 1] - ycedgeQ[u]) /
970 (xcedgeQ[u + 1] - xcedgeQ[u]);
971 cout <<
"incratio=" << incratio << endl;
972 if (
abs(incratio) > 4.0 && evaluatetan ==
false)
974 cout <<
"wrong=" << wrong << endl;
976 ASSERTL0(wrong < 2,
"number edges to change is too high!!");
977 edgeinterp[wrong] = k;
984 cout <<
"tan bef" << endl;
985 for (
int e = 0; e < nqedge; e++)
987 cout << xcedgeQ[e] <<
" " << ycedgeQ[e] <<
" "
988 << txedgeQ[e] << endl;
996 Vmath::Vcopy(npedge, &xcPhysMOD[k * npedge + 0], 1, &xPedges[0],
998 Vmath::Vcopy(npedge, &ycPhysMOD[k * npedge + 0], 1, &yPedges[0],
1001 PolyFit(polyorder, nqedge, xcedgeQ, ycedgeQ, coeffsinterp,
1002 xPedges, yPedges, npedge);
1004 Vmath::Vcopy(npedge, &xPedges[0], 1, &xcPhysMOD[k * npedge + 0],
1006 Vmath::Vcopy(npedge, &yPedges[0], 1, &ycPhysMOD[k * npedge + 0],
1013 Vmath::Vcopy(nqedge, &(txedgeQ[0]), 1, &(txQ[nqedge * k]), 1);
1014 Vmath::Vcopy(nqedge, &(tyedgeQ[0]), 1, &(tyQ[nqedge * k]), 1);
1019 for (
int w = 0;
w < fz.size();
w++)
1021 txQ[
w] = cos(atan(fz[
w]));
1022 tyQ[
w] = sin(atan(fz[
w]));
1023 cout << xcQ[
w] <<
" " << ycQ[
w] <<
" " << fz[
w] << endl;
1028 Vmath::Vcopy(nquad_lay, txQ, 1, Cont_y->UpdatePhys(), 1);
1029 Cont_y->FwdTransLocalElmt(Cont_y->GetPhys(), coeffsy);
1030 Cont_y->BwdTrans(coeffsy, Cont_y->UpdatePhys());
1033 Vmath::Vcopy(nquad_lay, tyQ, 1, Cont_y->UpdatePhys(), 1);
1034 Cont_y->FwdTransLocalElmt(Cont_y->GetPhys(), coeffsy);
1035 Cont_y->BwdTrans(coeffsy, Cont_y->UpdatePhys());
1046 for (
int q = 0;
q < 2;
q++)
1048 edgebef = edgeinterp[
q] - 1;
1050 (txQ[edgebef * nqedge + nqedge - 1] - txQ[edgebef * nqedge]) /
1051 (xcQ[edgebef * nqedge + nqedge - 1] - xcQ[edgebef * nqedge]);
1052 inc = (txQ[edgeinterp[
q] * nqedge + nqedge - 1] -
1053 txQ[edgeinterp[
q] * nqedge]) /
1054 (xcQ[edgeinterp[
q] * nqedge + nqedge - 1] -
1055 xcQ[edgeinterp[
q] * nqedge]);
1056 int npoints = 2 * nqedge;
1062 cout <<
"inc=" << inc <<
" incbef=" << incbefore << endl;
1063 if ((inc / incbefore) > 0.)
1065 cout <<
"before!!" << edgebef << endl;
1077 PolyFit(polyorder, npoints, xQedges, txQedges, coeffsinterp,
1078 xQedges, txQedges, npoints);
1082 &txQ[edgebef * nqedge + 0], 1);
1084 PolyFit(polyorder, npoints, xQedges, tyQedges, coeffsinterp,
1085 xQedges, tyQedges, npoints);
1089 &tyQ[edgebef * nqedge + 0], 1);
1093 cout <<
"after!!" << endl;
1105 PolyFit(polyorder, npoints, xQedges, txQedges, coeffsinterp,
1106 xQedges, txQedges, npoints);
1110 &txQ[edgeinterp[
q] * nqedge + 0], 1);
1112 PolyFit(polyorder, npoints, xQedges, tyQedges, coeffsinterp,
1113 xQedges, tyQedges, npoints);
1117 &tyQ[edgeinterp[
q] * nqedge + 0], 1);
1122 Vmath::Vcopy(nquad_lay, tyQ, 1, Cont_y->UpdatePhys(), 1);
1123 Cont_y->FwdTransLocalElmt(Cont_y->GetPhys(), coeffstmp);
1124 Cont_y->BwdTrans(coeffstmp, Cont_y->UpdatePhys());
1127 Vmath::Vcopy(nquad_lay, txQ, 1, Cont_y->UpdatePhys(), 1);
1128 Cont_y->FwdTransLocalElmt(Cont_y->GetPhys(), coeffstmp);
1129 Cont_y->BwdTrans(coeffstmp, Cont_y->UpdatePhys());
1132 for (
int k = 0; k < nedges; k++)
1139 Vmath::Vcopy(nqedge, &(txQ[nqedge * k]), 1, &(txedgeQ[0]), 1);
1140 Vmath::Vcopy(nqedge, &(tyQ[nqedge * k]), 1, &(tyedgeQ[0]), 1);
1142 Vmath::Vdiv(nqedge, txedgeQ, 1, tyedgeQ, 1, tx_tyedgeQ, 1);
1143 Vmath::Vmul(nqedge, tx_tyedgeQ, 1, tx_tyedgeQ, 1, tx_tyedgeQ, 1);
1144 Vmath::Sadd(nqedge, 1.0, tx_tyedgeQ, 1, nxedgeQ, 1);
1148 Vmath::Smul(nqedge, -1.0, nxedgeQ, 1, nxedgeQ, 1);
1149 Vmath::Vcopy(nqedge, &(nxedgeQ[0]), 1, &(nxQ[nqedge * k]), 1);
1151 Vmath::Vmul(nqedge, nxedgeQ, 1, nxedgeQ, 1, nyedgeQ, 1);
1152 Vmath::Smul(nqedge, -1.0, nyedgeQ, 1, nyedgeQ, 1);
1156 Vmath::Smul(nqedge, -1.0, nyedgeQ, 1, nyedgeQ, 1);
1157 Vmath::Vcopy(nqedge, &(nyedgeQ[0]), 1, &(nyQ[nqedge * k]), 1);
1159 cout <<
"edge:" << k << endl;
1160 cout <<
"tan/normal" << endl;
1161 for (
int r = 0; r < nqedge; r++)
1163 cout << xcQ[k * nqedge + r] <<
" " << txedgeQ[r] <<
" "
1164 << tyedgeQ[r] <<
" " << nxedgeQ[r] <<
" "
1165 << nyedgeQ[r] << endl;
1171 Vmath::Vcopy(nquad_lay, nyQ, 1, Cont_y->UpdatePhys(), 1);
1173 Cont_y->FwdTransLocalElmt(Cont_y->GetPhys(), coeffstmp);
1174 Cont_y->BwdTrans(coeffstmp, Cont_y->UpdatePhys());
1178 Vmath::Zero(Cont_y->GetNcoeffs(), Cont_y->UpdateCoeffs(), 1);
1179 Vmath::Vcopy(nquad_lay, nxQ, 1, Cont_y->UpdatePhys(), 1);
1180 Cont_y->FwdTransLocalElmt(Cont_y->GetPhys(), coeffstmp);
1181 Cont_y->BwdTrans(coeffstmp, Cont_y->UpdatePhys());
1185 for (
int k = 0; k < nedges; k++)
1191 nyQ[(k - 1) * nqedge + nqedge - 1] = nyQ[k * nqedge + 0];
1195 nxQ[(k - 1) * nqedge + nqedge - 1] = nxQ[k * nqedge + 0];
1203 cout <<
"nx,yQbefore" << endl;
1204 for (
int u = 0; u < xcQ.size(); u++)
1206 cout << xcQ[u] <<
" " << nyQ[u] <<
" " << txQ[u] << endl;
1212 cout <<
"nx,yQ" << endl;
1213 for (
int u = 0; u < x_tmpQ.size(); u++)
1215 cout << x_tmpQ[u] <<
" " << tmpnyQ[u] << endl;
1219 for (
int k = 0; k < nedges; k++)
1222 for (
int a = 0; a < npedge; a++)
1226 nxPhys[k * npedge + a] = nxQ[k * nqedge + 0];
1227 nyPhys[k * npedge + a] = nyQ[k * nqedge + 0];
1229 else if (a == npedge - 1)
1231 nxPhys[k * npedge + a] = nxQ[k * nqedge + nqedge - 1];
1232 nyPhys[k * npedge + a] = nyQ[k * nqedge + nqedge - 1];
1247 Addpointsx[k * (npedge - 2) + a - 1], x_tmpQ);
1257 Addpointsx[k * (npedge - 2) + a - 1], 4, Pxinterp,
1268 nxPhys[k * npedge + a] = -
sqrt(
abs(
1269 1 - nyPhys[k * npedge + a] * nyPhys[k * npedge + a]));
1284 nyPhys[(k - 1) * npedge + npedge - 1] =
1285 nyPhys[k * npedge + 0];
1289 nxPhys[(k - 1) * npedge + npedge - 1] =
1290 nxPhys[k * npedge + 0];
1294 cout <<
"xcPhys,," << endl;
1295 for (
int s = 0; s < np_lay; s++)
1298 cout << xcPhysMOD[s] <<
" " << ycPhysMOD[s] <<
" "
1299 << nxPhys[s] <<
" " << nyPhys[s] << endl;
1311 for (
int m = 0; m < nlays; m++)
1318 delta[m] = -(cntlow + 1 - m) * Delta0 / (cntlow + 1);
1322 delta[m] = (m - (cntlow)) * Delta0 / (cntlow + 1);
1328 for (
int h = 0; h < nvertl; h++)
1336 if (move_norm ==
false)
1338 ynew[lay_Vids[m][h]] = y_c[h] + delta[m];
1339 xnew[lay_Vids[m][h]] = x_c[h];
1343 if (h == 0 || h == nvertl - 1)
1345 ynew[lay_Vids[m][h]] = y_c[h] + delta[m];
1346 xnew[lay_Vids[m][h]] = x_c[h];
1350 ynew[lay_Vids[m][h]] =
1351 y_c[h] + delta[m] *
abs(nyPhys[h * npedge + 0]);
1352 xnew[lay_Vids[m][h]] =
1353 x_c[h] + delta[m] *
abs(nxPhys[h * npedge + 0]);
1356 cout <<
"Vid x=" << xnew[lay_Vids[m][h]]
1357 <<
" y=" << ynew[lay_Vids[m][h]] << endl;
1362 cout <<
"edge==" << h << endl;
1366 nyPhys[(h - 1) * npedge + npedge - 1],
1369 nxPhys[(h - 1) * npedge + npedge - 1],
1372 if (move_norm ==
false)
1375 layers_y[m][h * npedge + 0] = y_c[h] + delta[m];
1376 layers_x[m][h * npedge + 0] = xnew[lay_Vids[m][h]];
1378 layers_y[m][h * npedge + npedge - 1] =
1379 y_c[h + 1] + delta[m];
1380 layers_x[m][h * npedge + npedge - 1] =
1381 xnew[lay_Vids[m][h + 1]];
1383 for (
int d = 0;
d < npedge - 2;
d++)
1385 layers_y[m][h * npedge +
d + 1] =
1386 ycPhysMOD[h * npedge +
d + 1] + delta[m];
1388 layers_x[m][h * npedge +
d + 1] =
1389 xcPhysMOD[h * npedge +
d + 1];
1398 tmpy_lay[h * npedge + 0] = y_c[h] + delta[m];
1399 tmpx_lay[h * npedge + 0] = xnew[lay_Vids[m][h]];
1401 tmpy_lay[h * npedge + npedge - 1] =
1403 delta[m] *
abs(nyPhys[h * npedge + npedge - 1]);
1404 tmpx_lay[h * npedge + npedge - 1] =
1406 delta[m] *
abs(nxPhys[h * npedge + npedge - 1]);
1408 else if (h == nedges - 1)
1411 tmpy_lay[h * npedge + 0] =
1412 y_c[h] + delta[m] *
abs(nyPhys[h * npedge + 0]);
1413 tmpx_lay[h * npedge + 0] =
1414 x_c[h] + delta[m] *
abs(nxPhys[h * npedge + 0]);
1416 tmpy_lay[h * npedge + npedge - 1] =
1417 y_c[h + 1] + delta[m];
1418 tmpx_lay[h * npedge + npedge - 1] =
1419 xnew[lay_Vids[m][h + 1]];
1424 tmpy_lay[h * npedge + 0] =
1425 y_c[h] + delta[m] *
abs(nyPhys[h * npedge + 0]);
1426 tmpx_lay[h * npedge + 0] =
1427 x_c[h] + delta[m] *
abs(nxPhys[h * npedge + 0]);
1429 tmpy_lay[h * npedge + npedge - 1] =
1431 delta[m] *
abs(nyPhys[h * npedge + npedge - 1]);
1432 tmpx_lay[h * npedge + npedge - 1] =
1434 delta[m] *
abs(nxPhys[h * npedge + npedge - 1]);
1438 for (
int d = 0;
d < npedge - 2;
d++)
1441 tmpy_lay[h * npedge +
d + 1] =
1442 ycPhysMOD[h * npedge +
d + 1] +
1443 delta[m] *
abs(nyPhys[h * npedge +
d + 1]);
1446 tmpx_lay[h * npedge +
d + 1] =
1447 xcPhysMOD[h * npedge +
d + 1] +
1448 delta[m] *
abs(nxPhys[h * npedge +
d + 1]);
1466 for (
int s = 0; s < np_lay; s++)
1468 cout << tmpx_lay[s] <<
" " << tmpy_lay[s] << endl;
1471 cout <<
"fisrt interp" << endl;
1472 for (
int s = 0; s < np_lay; s++)
1474 cout << tmpx_lay[s] <<
" " << tmpy_lay[s] << endl;
1486 NekDouble boundright = xcPhysMOD[np_lay - 1];
1487 bool outboundleft =
false;
1488 bool outboundright =
false;
1489 if (tmpx_lay[1] < boundleft)
1491 outboundleft =
true;
1493 if (tmpx_lay[np_lay - 2] > boundright)
1495 outboundright =
true;
1503 for (
int r = 0; r < nedges; r++)
1506 if (tmpx_lay[r * npedge + npedge - 1] < boundleft &&
1507 outboundleft ==
true)
1509 vertout[outvert] = r;
1515 if (tmpx_lay[(r + 1) * npedge + npedge - 1] > boundleft)
1517 for (
int s = 0; s < npedge - 2; s++)
1519 if (tmpx_lay[(r + 1) * npedge + s + 1] <
1529 if (tmpx_lay[r * npedge + 0] > boundright &&
1530 outboundright ==
true)
1532 vertout[outvert] = r;
1538 if (tmpx_lay[(r - 1) * npedge + 0] < boundright)
1540 for (
int s = 0; s < npedge - 2; s++)
1542 if (tmpx_lay[(r - 1) * npedge + s + 1] >
1554 outcount = outvert * npedge + 1 + outmiddle;
1556 int replacepointsfromindex = 0;
1557 for (
int c = 0; c < nedges; c++)
1560 if (xcPhysMOD[c * npedge + npedge - 1] <=
1561 tmpx_lay[c * (npedge - (npedge - 2)) + 2] &&
1562 outboundright ==
true)
1564 replacepointsfromindex = c * (npedge - (npedge - 2)) + 2;
1569 if (xcPhysMOD[(nedges - 1 - c) * npedge + 0] >=
1570 tmpx_lay[np_lay - 1 -
1571 (c * (npedge - (npedge - 2)) + 2)] &&
1572 outboundleft ==
true)
1574 replacepointsfromindex =
1575 np_lay - 1 - (c * (npedge - (npedge - 2)) + 2);
1591 if (outboundright ==
true)
1593 pstart = replacepointsfromindex;
1594 shift = np_lay - outcount;
1596 (xcPhysMOD[np_lay - outcount] - xcPhysMOD[pstart]) /
1598 outcount = outcount - 1;
1599 ASSERTL0(tmpx_lay[np_lay - outcount] >
1600 xcPhysMOD[(nedges - 1) * npedge + 0],
1601 "no middle points in the last edge");
1606 pstart = outcount - 1;
1607 increment = (xcPhysMOD[replacepointsfromindex] -
1608 xcPhysMOD[pstart]) /
1611 xcPhysMOD[0 * npedge + npedge - 1],
1612 "no middle points in the first edge");
1629 NekDouble xctmp, ycinterp, nxinterp, nyinterp;
1631 for (
int v = 0; v < outcount; v++)
1633 xctmp = xcPhysMOD[pstart] + (v + 1) * increment;
1647 nxinterp =
sqrt(
abs(1 - nyinterp * nyinterp));
1655 replace_x[v] = xctmp + delta[m] *
abs(nxinterp);
1656 replace_y[v] = ycinterp + delta[m] *
abs(nyinterp);
1657 tmpx_lay[v + shift] = replace_x[v];
1658 tmpy_lay[v + shift] = replace_y[v];
1678 int closepoints = 4;
1684 int pointscount = 0;
1685 for (
int q = 0;
q < np_lay;
q++)
1687 for (
int e = 0; e < nedges; e++)
1689 if (tmpx_lay[
q] <= x_c[e + 1] && tmpx_lay[
q] >= x_c[e])
1693 if (
q == e * npedge + npedge - 1 && pointscount != npedge)
1698 else if (
q == e * npedge + npedge - 1)
1719 lay_Vids[m], layers_x[m], layers_y[m], xnew,
1799 if (curv_lay ==
true)
1808 int npoints = npedge;
1811 for (
int f = 0; f < nedges; f++)
1822 PolyFit(polyorder, npoints, xPedges, yPedges, coeffsinterp,
1823 xPedges, yPedges, npoints);
1827 &layers_y[m][(f)*npedge + 0], 1);
1830 layers_y[m][f * npedge + 0] = ynew[lay_Vids[m][f]];
1831 layers_y[m][f * npedge + npedge - 1] = ynew[lay_Vids[m][f + 1]];
1834 cout <<
" xlay ylay lay:" << m << endl;
1835 for (
int l = 0; l < np_lay; l++)
1838 cout << std::setprecision(8) << layers_x[m][l] <<
" "
1839 << layers_y[m][l] << endl;
1875 cout <<
"lay=" << m << endl;
1878 " different layer ymin val");
1881 " different layer ymax val");
1884 " different layer xmin val");
1887 " different layer xmax val");
1897 npedge, graphShPt, xold_c, yold_c, xold_low, yold_low, xold_up,
1898 yold_up, layers_x[0], layers_y[0], layers_x[nlays - 1],
1899 layers_y[nlays - 1], nxPhys, nyPhys, xnew, ynew);
1980 Down, Up, xnew, ynew, layers_x, layers_y);
1990 cout << std::setprecision(8) <<
"xmin=" <<
Vmath::Vmin(nVertTot, xnew, 1)
1993 " different xmin val");
1995 " different ymin val");
1997 " different xmax val");
1999 " different ymax val");
2004 charalp, layers_x, layers_y, lay_Eids, curv_lay);
2013 int nvertl = nedges + 1;
2017 for (
int j = 0; j < nedges; j++)
2021 edge = (bndSegExplow->GetGeom1D())->GetGlobalID();
2023 for (
int k = 0; k < 2; k++)
2025 Vids_temp[j + k] = (bndSegExplow->GetGeom1D())->GetVid(k);
2027 graphShPt->GetVertex(Vids_temp[j + k]);
2029 vertex->GetCoords(x1, y1, z1);
2030 if (x1 == x_connect && edge != lastedge)
2033 if (x_connect == x0layer)
2035 Vids[v1] = Vids_temp[j + k];
2042 (bndSegExplow->GetGeom1D())->GetVid(1);
2043 Vids[v2] = Vids_temp[j + 1];
2045 graphShPt->GetVertex(Vids[v2]);
2047 vertex->GetCoords(x2, y2, z2);
2054 (bndSegExplow->GetGeom1D())->GetVid(0);
2055 Vids[v2] = Vids_temp[j + 0];
2057 graphShPt->GetVertex(Vids[v2]);
2059 vertex->GetCoords(x2, y2, z2);
2069 (bndSegExplow->GetGeom1D())->GetVid(1);
2070 Vids[v1] = Vids_temp[j + 1];
2072 graphShPt->GetVertex(Vids[v1]);
2074 vertex->GetCoords(x1, y1, z1);
2081 (bndSegExplow->GetGeom1D())->GetVid(0);
2082 Vids[v1] = Vids_temp[j + 0];
2084 graphShPt->GetVertex(Vids[v1]);
2086 vertex->GetCoords(x1, y1, z1);
2094 if (Vids[v1] != -10)
2113 cout <<
"Computestreakpositions" << endl;
2114 int nq = streak->GetTotPoints();
2128 Vmath::Vadd(xc.size(), yold_up, 1, yold_low, 1, yc, 1);
2142 for (
int e = 0; e < npoints; e++)
2147 elmtid = streak->GetExpIndex(coord, 0.00001);
2148 offset = streak->GetPhys_Offset(elmtid);
2154 while (
abs(F) > 0.000000001)
2157 elmtid = streak->GetExpIndex(coord, 0.00001);
2164 if ((
abs(coord[1]) > 1 || elmtid == -1) && attempt == 0 &&
2168 coord[1] = yold_c[e];
2171 else if ((
abs(coord[1]) > 1 || elmtid == -1))
2173 coord[1] = ytmp + 0.01;
2174 elmtid = streak->GetExpIndex(coord, 0.001);
2175 offset = streak->GetPhys_Offset(elmtid);
2176 NekDouble Utmp = streak->GetExp(elmtid)->PhysEvaluate(
2177 coord, streak->GetPhys() + offset);
2178 NekDouble dUtmp = streak->GetExp(elmtid)->PhysEvaluate(
2179 coord, derstreak + offset);
2180 coord[1] = coord[1] - (Utmp - cr) / dUtmp;
2181 if ((
abs(Utmp - cr) >
abs(F)) || (
abs(coord[1]) > 1))
2183 coord[1] = ytmp - 0.01;
2190 ASSERTL0(
abs(coord[1]) <= 1,
" y value out of bound +/-1");
2192 offset = streak->GetPhys_Offset(elmtid);
2193 U = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() +
2196 streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
2197 coord[1] = coord[1] - (U - cr) / dU;
2199 ASSERTL0(coord[0] == xc[e],
" x coordinate must remain the same");
2202 if (its > 200 &&
abs(F) < 0.00001)
2204 cout <<
"warning streak position obtained with precision:" << F
2208 else if (its > 1000 &&
abs(F) < 0.0001)
2210 cout <<
"warning streak position obtained with precision:" << F
2214 else if (its > 1000)
2216 ASSERTL0(
false,
"no convergence after 1000 iterations");
2219 yc[e] = coord[1] - (U - cr) / dU;
2220 ASSERTL0(U <= cr + tol,
"streak wrong+");
2221 ASSERTL0(U >= cr - tol,
"streak wrong-");
2224 cout <<
"result streakvert x=" << xc[e] <<
" y=" << yc[e]
2225 <<
" streak=" << U << endl;
2246 while (
abs(F) > 0.00000001)
2250 elmtid = function->GetExpIndex(coords, 0.01);
2252 cout <<
"gen newton xi=" << xi <<
" yi=" << coords[1]
2253 <<
" elmtid=" << elmtid <<
" F=" << F << endl;
2255 if ((
abs(coords[1]) > 1 || elmtid == -1))
2258 coords[1] = ytmp + 0.01;
2259 elmtid = function->GetExpIndex(coords, 0.01);
2260 offset = function->GetPhys_Offset(elmtid);
2261 NekDouble Utmp = function->GetExp(elmtid)->PhysEvaluate(
2262 coords, function->GetPhys() + offset);
2263 NekDouble dUtmp = function->GetExp(elmtid)->PhysEvaluate(
2264 coords, derfunction + offset);
2265 coords[1] = coords[1] - (Utmp - cr) / dUtmp;
2266 cout <<
"attempt:" << coords[1] << endl;
2267 if ((
abs(Utmp - cr) >
abs(F)) || (
abs(coords[1]) > 1.01))
2269 coords[1] = ytmp - 0.01;
2274 else if (
abs(coords[1]) < 1.01 && attempt == 0)
2281 ASSERTL0(
abs(coords[1]) <= 1.00,
" y value out of bound +/-1");
2283 offset = function->GetPhys_Offset(elmtid);
2284 U = function->GetExp(elmtid)->PhysEvaluate(coords, function->GetPhys() +
2286 dU = function->GetExp(elmtid)->PhysEvaluate(coords,
2287 derfunction + offset);
2288 coords[1] = coords[1] - (U - cr) / dU;
2289 cout << cr <<
"U-cr=" << U - cr <<
" tmp result y:" << coords[1]
2290 <<
" dU=" << dU << endl;
2294 if (its > 200 &&
abs(F) < 0.00001)
2296 cout <<
"warning streak position obtained with precision:" << F
2300 else if (its > 1000 &&
abs(F) < 0.0001)
2302 cout <<
"warning streak position obtained with precision:" << F
2306 else if (its > 1000)
2308 ASSERTL0(
false,
"no convergence after 1000 iterations");
2311 ASSERTL0(coords[0] == xi,
" x coordinate must remain the same");
2314 yout = coords[1] - (U - cr) / dU;
2315 cout <<
"NewtonIt result x=" << xout <<
" y=" << coords[1] <<
" U=" << U
2323 const std::shared_ptr<LocalRegions::ExpansionVector> exp2D =
2325 int nel = exp2D->size();
2333 for (
int i = 0; i < nel; i++)
2335 if ((locQuadExp = (*exp2D)[i]->as<LocalRegions::QuadExp>()))
2337 for (
int j = 0; j < locQuadExp->GetNtraces(); ++j)
2339 SegGeom = (locQuadExp->GetGeom2D())->GetEdge(j);
2340 id = SegGeom->GetGlobalID();
2342 if (V1tmp[
id] == 10000)
2344 V1tmp[id] = SegGeom->GetVertex(0)->GetVid();
2345 V2tmp[id] = SegGeom->GetVertex(1)->GetVid();
2353 else if ((locTriExp = (*exp2D)[i]->as<LocalRegions::TriExp>()))
2355 for (
int j = 0; j < locTriExp->GetNtraces(); ++j)
2357 SegGeom = (locTriExp->GetGeom2D())->GetEdge(j);
2358 id = SegGeom->GetGlobalID();
2360 if (V1tmp[
id] == 10000)
2362 V1tmp[id] = SegGeom->GetVertex(0)->GetVid();
2363 V2tmp[id] = SegGeom->GetVertex(1)->GetVid();
2373 for (
int g = 0; g < cnt; g++)
2376 ASSERTL0(V1tmp[g] != 10000,
"V1 wrong");
2378 ASSERTL0(V2tmp[g] != 10000,
"V2 wrong");
2393 int nlay_Eids = xcold.size() - 1;
2394 int nlay_Vids = xcold.size();
2396 int nVertsTot = mesh->GetNvertices();
2397 cout <<
"nverttot=" << nVertsTot << endl;
2401 cout <<
"init nlays=" << nlays << endl;
2408 cout <<
"yoldup=" << yoldup[0] << endl;
2409 cout <<
"yolddown=" << yolddown[0] << endl;
2411 for (
int r = 0; r < nVertsTot; r++)
2416 vertex->GetCoords(x, y,
z);
2422 if (y <= yoldup[0] && y >= yolddown[0] && y != ycold[0])
2425 tmpVids0[nlays] = r;
2432 cout <<
"nlays=" << nlays << endl;
2444 for (
int w = 0;
w < nlays;
w++)
2447 tmpx0[
w] = tmpx[index];
2448 tmpy0[
w] = tmpf[index];
2449 tmpVids0[
w] = tmpV[index];
2450 tmpf[index] = max + 1000;
2456 for (
int m = 0; m < nlays; m++)
2467 NekDouble xtmp, ytmp, normnext = 0.0, xnext = 0.0, ynext = 0.0, diff;
2468 NekDouble Ubef = 0.0, Utmp = 0.0, Unext = 0.0;
2471 int nTotEdges = V1.size();
2473 for (
int m = 0; m < nlays; m++)
2475 for (
int g = 0; g < nlay_Eids; g++)
2479 for (
int h = 0; h < nTotEdges; h++)
2482 if (tmpVids0[m] == V1[h])
2485 mesh->GetVertex(V2[h]);
2487 vertex->GetCoords(x, y,
z);
2491 Vids_lay[m][0] = V1[h];
2492 Vids_lay[m][1] = V2[h];
2494 mesh->GetVertex(V1[h]);
2496 vertex1->GetCoords(x1, y1, z1);
2498 sqrt((y - y1) * (y - y1) + (x - x1) * (x - x1));
2503 elmtid = streak->GetExpIndex(coord, 0.00001);
2504 offset = streak->GetPhys_Offset(elmtid);
2505 Ubef = streak->GetExp(elmtid)->PhysEvaluate(
2506 coord, streak->GetPhys() + offset);
2509 if (tmpVids0[m] == V2[h])
2512 mesh->GetVertex(V1[h]);
2514 vertex->GetCoords(x, y,
z);
2518 Vids_lay[m][0] = V2[h];
2519 Vids_lay[m][1] = V1[h];
2521 mesh->GetVertex(V2[h]);
2524 sqrt((y - y2) * (y - y2) + (x - x2) * (x - x2));
2529 elmtid = streak->GetExpIndex(coord, 0.00001);
2530 offset = streak->GetPhys_Offset(elmtid);
2531 Ubef = streak->GetExp(elmtid)->PhysEvaluate(
2532 coord, streak->GetPhys() + offset);
2537 cout <<
"Eid=" << Eids_lay[m][0]
2538 <<
" Vids_lay0=" << Vids_lay[m][0]
2539 <<
" Vidslay1=" << Vids_lay[m][1] << endl;
2545 for (
int h = 0; h < nTotEdges; h++)
2549 if ((Vids_lay[m][g] == V1[h] || Vids_lay[m][g] == V2[h]) &&
2550 h != Eids_lay[m][g - 1])
2552 cout <<
"edgetmp=" << h << endl;
2553 ASSERTL0(cnt <= 6,
"wrong number of candidates");
2562 cout <<
"normbef=" << normbef << endl;
2563 cout <<
"Ubefcc=" << Ubef << endl;
2565 for (
int e = 0; e < cnt; e++)
2568 mesh->GetVertex(V1[edgestmp[e]]);
2570 vertex1->GetCoords(x1, y1, z1);
2572 mesh->GetVertex(V2[edgestmp[e]]);
2574 vertex2->GetCoords(x2, y2, z2);
2577 sqrt((y2 - y1) * (y2 - y1) + (x2 - x1) * (x2 - x1));
2579 cout <<
"edgetmp1=" << edgestmp[e] << endl;
2580 cout <<
"V1 x1=" << x1 <<
" y1=" << y1 << endl;
2581 cout <<
"V2 x2=" << x2 <<
" y2=" << y2 << endl;
2582 if (Vids_lay[m][g] == V1[edgestmp[e]])
2589 elmtid = streak->GetExpIndex(coord, 0.00001);
2590 offset = streak->GetPhys_Offset(elmtid);
2591 Utmp = streak->GetExp(elmtid)->PhysEvaluate(
2592 coord, streak->GetPhys() + offset);
2593 diffarray[e] =
abs((xtmp * xbef + ytmp * ybef) /
2594 (normtmp * normbef) -
2596 diffUarray[e] =
abs(Ubef - Utmp);
2597 cout <<
" normtmp=" << normtmp << endl;
2598 cout <<
" Utmpcc=" << Utmp << endl;
2599 cout << xtmp <<
" ytmp=" << ytmp <<
" diff="
2600 <<
abs(((xtmp * xbef + ytmp * ybef) /
2601 (normtmp * normbef)) -
2604 if (
abs((xtmp * xbef + ytmp * ybef) /
2605 (normtmp * normbef) -
2607 y2 <= yoldup[g + 1] && y2 >= yolddown[g + 1] &&
2608 y1 <= yoldup[g] && y1 >= yolddown[g])
2611 Eids_lay[m][g] = edgestmp[e];
2612 Vids_lay[m][g + 1] = V2[edgestmp[e]];
2613 diff =
abs((xtmp * xbef + ytmp * ybef) /
2614 (normtmp * normbef) -
2622 else if (Vids_lay[m][g] == V2[edgestmp[e]])
2629 elmtid = streak->GetExpIndex(coord, 0.00001);
2630 offset = streak->GetPhys_Offset(elmtid);
2631 Utmp = streak->GetExp(elmtid)->PhysEvaluate(
2632 coord, streak->GetPhys() + offset);
2633 diffarray[e] =
abs((xtmp * xbef + ytmp * ybef) /
2634 (normtmp * normbef) -
2636 diffUarray[e] =
abs(Ubef - Utmp);
2637 cout <<
" normtmp=" << normtmp << endl;
2638 cout <<
" Utmpcc=" << Utmp << endl;
2639 cout << xtmp <<
" ytmp=" << ytmp <<
" diff="
2640 <<
abs(((xtmp * xbef + ytmp * ybef) /
2641 (normtmp * normbef)) -
2644 if (
abs((xtmp * xbef + ytmp * ybef) /
2645 (normtmp * normbef) -
2647 y2 <= yoldup[g] && y2 >= yolddown[g] &&
2648 y1 <= yoldup[g + 1] && y1 >= yolddown[g + 1])
2650 Eids_lay[m][g] = edgestmp[e];
2651 Vids_lay[m][g + 1] = V1[edgestmp[e]];
2652 diff =
abs((xtmp * xbef + ytmp * ybef) /
2653 (normtmp * normbef) -
2666 cout <<
"Eid before check=" << Eids_lay[m][g] << endl;
2667 for (
int q = 0;
q < cnt;
q++)
2669 cout <<
q <<
" diff" << diffarray[
q] << endl;
2673 if (m > 0 && m < nlays)
2676 Vids_lay[m][g + 1]);
2680 cout <<
"COMMON VERT" << endl;
2682 diffarray[eid] = 1000;
2686 mesh->GetVertex(V1[edgestmp[eid]]);
2688 vertex1->GetCoords(x1, y1, z1);
2690 mesh->GetVertex(V2[edgestmp[eid]]);
2692 vertex2->GetCoords(x2, y2, z2);
2695 sqrt((y2 - y1) * (y2 - y1) + (x2 - x1) * (x2 - x1));
2697 Eids_lay[m][g] = edgestmp[eid];
2698 if (Vids_lay[m][g] == V1[edgestmp[eid]])
2702 elmtid = streak->GetExpIndex(coord, 0.00001);
2703 offset = streak->GetPhys_Offset(elmtid);
2704 Utmp = streak->GetExp(elmtid)->PhysEvaluate(
2705 coord, streak->GetPhys() + offset);
2706 Vids_lay[m][g + 1] = V2[edgestmp[eid]];
2713 if (Vids_lay[m][g] == V2[edgestmp[eid]])
2717 elmtid = streak->GetExpIndex(coord, 0.00001);
2718 offset = streak->GetPhys_Offset(elmtid);
2719 Utmp = streak->GetExp(elmtid)->PhysEvaluate(
2720 coord, streak->GetPhys() + offset);
2721 Vids_lay[m][g + 1] = V1[edgestmp[eid]];
2730 cout << m <<
"edge aft:" << Eids_lay[m][g]
2731 <<
" Vid=" << Vids_lay[m][g + 1] << endl;
2737 cout <<
"endelse" << normtmp << endl;
2744 for (
int w = 0;
w < nlays;
w++)
2746 for (
int f = 0; f < nlay_Eids; f++)
2748 cout <<
"check=" <<
w <<
" Eid:" << Eids_lay[
w][f] << endl;
2757 for (
int u = 0; u < Vids_laybefore.size(); u++)
2759 if (Vids_laybefore[u] == Vid || Vids_c[u] == Vid)
2763 cout << Vid <<
" Vert test=" << Vids_laybefore[u] << endl;
2773 int np_lay = inarray.size();
2774 ASSERTL0(inarray.size() % nedges == 0,
" something on number npedge");
2778 for (
int w = 0;
w < np_lay;
w++)
2783 if (inarray[
w] == inarray[
w + 1])
2788 outarray[cnt] = inarray[
w];
2792 ASSERTL0(cnt == np_lay - (nedges - 1),
"wrong cut");
2797 int npts = xArray.size();
2806 if (xArray[index] > x)
2819 ASSERTL0(neighpoints % 2 == 0,
"number of neighbour points should be even");
2820 int leftpoints = (neighpoints / 2) - 1;
2821 int rightpoints = neighpoints / 2;
2826 if (index - leftpoints < 0)
2829 diff = index - leftpoints;
2831 Vmath::Vcopy(neighpoints, &yArray[0], 1, &Neighbour_y[0], 1);
2832 Vmath::Vcopy(neighpoints, &xArray[0], 1, &Neighbour_x[0], 1);
2834 else if ((yArray.size() - 1) - index < rightpoints)
2837 int rpoints = (yArray.size() - 1) - index;
2838 diff = rightpoints - rpoints;
2840 start = index - leftpoints - diff;
2841 Vmath::Vcopy(neighpoints, &yArray[start], 1, &Neighbour_y[0], 1);
2842 Vmath::Vcopy(neighpoints, &xArray[start], 1, &Neighbour_x[0], 1);
2847 start = index - leftpoints;
2848 Vmath::Vcopy(neighpoints, &yArray[start], 1, &Neighbour_y[0], 1);
2849 Vmath::Vcopy(neighpoints, &xArray[start], 1, &Neighbour_x[0], 1);
2858 for (
int f = 1; f < neighpoints; f++)
2860 ASSERTL0(Neighbour_x[f] != Neighbour_x[f - 1],
2861 " repetition on NeighbourArrays");
2872 for (
int pt = 0; pt < npts; ++pt)
2876 for (
int j = 0; j < pt; ++j)
2878 h = h * (x - xpts[j]) / (xpts[pt] - xpts[j]);
2881 for (
int k = pt + 1; k < npts; ++k)
2883 h = h * (x - xpts[k]) / (xpts[pt] - xpts[k]);
2887 sum += funcvals[pt] * LagrangePoly;
2899 int np_pol = coeffsinterp.size();
2900 cout <<
"evaluatetan with " << np_pol << endl;
2906 for (
int q = 0;
q < npoints;
q++)
2908 polorder = np_pol - 1;
2909 derorder = np_pol - 2;
2911 for (
int d = 0;
d < np_pol - 1;
d++)
2913 yprime[
q] += (derorder + 1) * coeffsinterp[
d] *
2914 std::pow(xcQedge[
q], derorder);
2918 for (
int a = 0; a < np_pol; a++)
2920 yinterp[
q] += coeffsinterp[a] * std::pow(xcQedge[
q], polorder);
2928 for (
int n = 0; n < npoints; n++)
2932 txQedge[n] = cos((atan((yprime[n]))));
2933 tyQedge[n] = sin((atan((yprime[n]))));
2934 cout << xcQedge[n] <<
" " << yinterp[n] <<
" " << yprime[n]
2935 <<
" " << txQedge[n] <<
" " << tyQedge[n] << endl;
2942 int edge,
int npedge)
2944 int np_pol = xpol.size();
2950 for (
int e = 0; e < N; e++)
2953 for (
int w = 0;
w < N;
w++)
2955 A[N * e + row] = std::pow(xpol[
w], N - 1 - e);
2960 for (
int r = 0; r < np_pol; r++)
2974 std::string message =
2975 "ERROR: The " + std::to_string(-info) +
2976 "th parameter had an illegal parameter for dgetrf";
2981 std::string message =
"ERROR: Element u_" + std::to_string(info) +
2982 std::to_string(info) +
" is 0 from dgetrf";
2988 Lapack::Dgetrs(
'N', N, ncolumns_b,
A.get(), N, ipivot.get(), b.get(), N,
2992 std::string message =
2993 "ERROR: The " + std::to_string(-info) +
2994 "th parameter had an illegal parameter for dgetrf";
2999 std::string message =
"ERROR: Element u_" + std::to_string(info) +
3000 std::to_string(info) +
" is 0 from dgetrf";
3012 for (
int c = 0; c < npedge; c++)
3014 polorder = np_pol - 1;
3016 ycout[edge * (npedge) + c + 1] = 0;
3017 for (
int d = 0;
d < np_pol;
d++)
3019 ycout[edge * (npedge) + c + 1] +=
3020 b[
d] * std::pow(xcout[edge * (npedge) + c + 1], polorder);
3035 int N = polyorder + 1;
3038 cout << npoints << endl;
3039 for (
int u = 0; u < npoints; u++)
3041 cout <<
"c=" << xin[u] <<
" " << fin[u] << endl;
3045 for (
int e = 0; e < N; e++)
3048 for (
int row = 0; row < N; row++)
3050 for (
int w = 0;
w < npoints;
w++)
3052 A[N * e + row] += std::pow(xin[
w], e + row);
3057 for (
int row = 0; row < N; row++)
3059 for (
int h = 0; h < npoints; h++)
3061 b[row] += fin[h] * std::pow(xin[h], row);
3074 std::string message =
3075 "ERROR: The " + std::to_string(-info) +
3076 "th parameter had an illegal parameter for dgetrf";
3081 std::string message =
"ERROR: Element u_" + std::to_string(info) +
3082 std::to_string(info) +
" is 0 from dgetrf";
3087 Lapack::Dgetrs(
'N', N, ncolumns_b,
A.get(), N, ipivot.get(), b.get(), N,
3091 std::string message =
3092 "ERROR: The " + std::to_string(-info) +
3093 "th parameter had an illegal parameter for dgetrf";
3098 std::string message =
"ERROR: Element u_" + std::to_string(info) +
3099 std::to_string(info) +
" is 0 from dgetrf";
3108 for (
int j = 0; j < N; j++)
3110 b[j] = tmp[cnt - 1];
3114 for (
int h = 0; h < N; h++)
3116 cout <<
"coeff:" << b[h] << endl;
3121 for (
int c = 0; c < npout; c++)
3123 polorder = polyorder;
3126 for (
int d = 0;
d < N;
d++)
3129 fout[c] += b[
d] * std::pow(xout[c], polorder);
3152 for (
int w = 0;
w < tmpx.size();
w++)
3155 outarray_x[
w] = tmpx[index];
3156 outarray_y[
w] = tmpy[index];
3157 if (
w < tmpx.size() - 1)
3159 if (tmpx[index] == tmpx[index + 1])
3161 outarray_x[
w + 1] = tmpx[index + 1];
3162 outarray_y[
w + 1] = tmpy[index + 1];
3163 tmpx[index + 1] = max + 1000;
3179 tmpx[index] = max + 1000;
3192 int np_lay = layers_y[0].size();
3194 for (
int h = 1; h < nlays - 1; h++)
3197 for (
int s = 0; s < nvertl; s++)
3200 ASSERTL0(ynew[lay_Vids[h][s]] == -20,
"ynew layers not empty");
3204 ynew[lay_Vids[h][s]] =
3206 h *
abs(ynew[Down[s]] - yc[s]) / (cntlow + 1);
3208 xnew[lay_Vids[h][s]] = xc[s];
3212 layers_y[h][s] = ynew[lay_Vids[h][s]];
3213 layers_x[h][s] = xnew[lay_Vids[h][s]];
3218 ynew[lay_Vids[h][s]] = yc[s] + (h - cntlow) *
3219 abs(ynew[Up[s]] - yc[s]) /
3222 xnew[lay_Vids[h][s]] = xc[s];
3224 layers_y[h][s] = ynew[lay_Vids[h][s]];
3225 layers_x[h][s] = xnew[lay_Vids[h][s]];
3239 int np_lay = xcPhys.size();
3240 int nedges = nvertl - 1;
3247 int closepoints = 4;
3252 for (
int g = 0; g < nedges; g++)
3261 Pxinterp, Pyinterp);
3262 xnew[Vids[g]] = xcPhys[g * npedge + 0];
3263 ylay[g * npedge + 0] = ynew[Vids[g]];
3264 xlay[g * npedge + 0] = xnew[Vids[g]];
3274 xcPhys[g * npedge + npedge - 1], closepoints, Pxinterp, Pyinterp);
3275 xnew[Vids[g + 1]] = xcPhys[g * npedge + npedge - 1];
3276 ylay[g * npedge + npedge - 1] = ynew[Vids[g + 1]];
3277 xlay[g * npedge + npedge - 1] = xnew[Vids[g + 1]];
3280 for (
int r = 0; r < npedge - 2; r++)
3287 ASSERTL0(index <= tmpy.size() - 1,
" index wrong");
3293 xcPhys[g * npedge + r + 1], closepoints, Pxinterp, Pyinterp);
3295 xlay[g * npedge + r + 1] = xcPhys[g * npedge + r + 1];
3316 int np_lay = xcPhys.size();
3317 int nedges = nvertl - 1;
3325 int closepoints = 4;
3330 for (
int g = 0; g < nedges; g++)
3334 ynew[Vids[g]] = tmpy_lay[g * npedge + 0];
3335 xnew[Vids[g]] = tmpx_lay[g * npedge + 0];
3337 ylay[g * npedge + 0] = ynew[Vids[g]];
3338 xlay[g * npedge + 0] = xnew[Vids[g]];
3341 ynew[Vids[g + 1]] = tmpy_lay[g * npedge + npedge - 1];
3342 xnew[Vids[g + 1]] = tmpx_lay[g * npedge + npedge - 1];
3343 ylay[g * npedge + npedge - 1] = ynew[Vids[g + 1]];
3344 xlay[g * npedge + npedge - 1] = xnew[Vids[g + 1]];
3347 for (
int r = 0; r < npedge - 2; r++)
3349 x0 = xlay[g * npedge + 0];
3350 x1 = xlay[g * npedge + npedge - 1];
3351 xtmp = x0 + r * (x1 - x0) / (npedge - 1);
3356 ASSERTL0(index <= tmpy.size() - 1,
" index wrong");
3361 ylay[g * npedge + r + 1] =
3364 xlay[g * npedge + r + 1] = xtmp;
3381 int nvertl = ycold.size();
3382 int nVertTot = mesh->GetNvertices();
3383 for (
int n = 0; n < nVertTot; n++)
3388 vertex->GetCoords(x, y,
z);
3393 for (
int k = 0; k < nvertl; k++)
3395 if (
abs(x - xcold[k]) < tmp)
3397 tmp =
abs(x - xcold[k]);
3409 nplay_closer = (qp_closer - 1) * npedge + npedge - 1;
3412 if (y > yoldup[qp_closer] && y < 1)
3418 ratio = (1 - ylayup[nplay_closer]) / ((1 - yoldup[qp_closer]));
3420 ynew[n] = ylayup[nplay_closer] + (y - yoldup[qp_closer]) * ratio;
3423 else if (y < yolddown[qp_closer] && y > -1)
3426 ratio = (1 + ylaydown[nplay_closer]) / ((1 + yolddown[qp_closer]));
3429 ylaydown[nplay_closer] + (y - yolddown[qp_closer]) * ratio;
3447 int nvertl = xoldup.size();
3448 int nedges = nvertl - 1;
3457 for (
int a = 0; a < nedges; a++)
3462 xnew_down[a] = xlaydown[a * npedge + 0];
3463 ynew_down[a] = ylaydown[a * npedge + 0];
3464 xnew_up[a] = xlayup[a * npedge + 0];
3465 ynew_up[a] = ylayup[a * npedge + 0];
3466 nxvert[a] = nxPhys[a * npedge + 0];
3467 nyvert[a] = nyPhys[a * npedge + 0];
3469 xnew_down[a + 1] = xlaydown[a * npedge + npedge - 1];
3470 ynew_down[a + 1] = ylaydown[a * npedge + npedge - 1];
3471 xnew_up[a + 1] = xlayup[a * npedge + npedge - 1];
3472 ynew_up[a + 1] = ylayup[a * npedge + npedge - 1];
3473 nxvert[a + 1] = nxPhys[a * npedge + npedge - 1];
3474 nyvert[a + 1] = nyPhys[a * npedge + npedge - 1];
3479 xnew_down[a + 1] = xlaydown[a * npedge + npedge - 1];
3480 ynew_down[a + 1] = ylaydown[a * npedge + npedge - 1];
3481 xnew_up[a + 1] = xlayup[a * npedge + npedge - 1];
3482 ynew_up[a + 1] = ylayup[a * npedge + npedge - 1];
3483 nxvert[a + 1] = nxPhys[a * npedge + npedge - 1];
3484 nyvert[a + 1] = nyPhys[a * npedge + npedge - 1];
3495 int nVertTot = mesh->GetNvertices();
3496 for (
int n = 0; n < nVertTot; n++)
3501 vertex->GetCoords(x, y,
z);
3502 int qp_closeroldup = 0, qp_closerolddown = 0;
3508 for (
int k = 0; k < nvertl; k++)
3510 if (
abs(x - xolddown[k]) < diffdown)
3512 diffdown =
abs(x - xolddown[k]);
3513 qp_closerolddown = k;
3515 if (
abs(x - xoldup[k]) < diffup)
3517 diffup =
abs(x - xoldup[k]);
3526 int qp_closerup = 0, qp_closerdown = 0;
3528 for (
int f = 0; f < nvertl; f++)
3530 if (
abs(x - xnew_down[f]) < diffdown)
3532 diffdown =
abs(x - xnew_down[f]);
3535 if (
abs(x - xnew_up[f]) < diffup)
3537 diffup =
abs(x - xnew_up[f]);
3564 int qp_closernormup;
3575 int qp_closernormdown;
3583 if (y > yoldup[qp_closeroldup] && y < 1)
3589 ratio = (1 - ynew_up[qp_closerup]) / ((1 - yoldup[qp_closeroldup]));
3595 ynew_up[qp_closerup] + (y - yoldup[qp_closeroldup]) * ratio;
3603 if (x > (xmax - xmin) / 2. && x < xmax)
3605 ratiox = (xmax - xnew_up[qp_closernormup]) /
3606 (xmax - xoldup[qp_closernormup]);
3607 if ((xmax - xoldup[qp_closernormup]) == 0)
3615 x +
abs(nxvert[qp_closernormup]) *
3616 (xnew_up[qp_closeroldup] - xoldup[qp_closeroldup]) *
3618 ASSERTL0(x > xmin,
" x value <xmin up second half");
3619 ASSERTL0(x < xmax, " x value >xmax up second half
");
3621 else if (x > xmin && x <= (xmax - xmin) / 2.)
3623 // cout<<"up close normold=
"<<qp_closernormoldup<<"
3625 ratiox = (xnew_up[qp_closernormup] - xmin) /
3626 ((xoldup[qp_closernormup] - xmin));
3627 if ((xoldup[qp_closernormup] - xmin) == 0)
3634 x +
abs(nxvert[qp_closernormup]) *
3635 (xnew_up[qp_closeroldup] - xoldup[qp_closeroldup]) *
3638 ASSERTL0(x > xmin,
" x value <xmin up first half");
3639 ASSERTL0(x < xmax, " x value >xmax up first half
");
3642 else if (y < yolddown[qp_closerolddown] && y > -1)
3645 ratio = (1 + ynew_down[qp_closerdown]) /
3646 ((1 + yolddown[qp_closerolddown]));
3648 // ratioy = (1-ynew_down[qp_closernormdown])/
3649 // ( (1-yolddown[qp_closernormolddown]) );
3651 // distance prop to layerlow
3652 ynew[n] = ynew_down[qp_closerdown] +
3653 (y - yolddown[qp_closerolddown]) * ratio;
3654 // ynew[n] = y +abs(nyvert[qp_closernormdown])*
3655 // (ynew_down[qp_closerolddown]-yolddown[qp_closerolddown])*ratioy;
3657 // 0.3*(ynew_down[qp_closerdown]-yolddown[qp_closerdown]); xnew[n] =
3659 // abs(nxvert[qp_closerolddown])*(xnew_down[qp_closerolddown]-xolddown[qp_closerolddown]);
3663 cout<<qp_closerolddown<<" nplaydown=
"<<qp_closerdown<<endl;
3664 cout<<"xolddown=
"<<xolddown[qp_closerolddown]<<"
3665 xnewdown=
"<<xnew_down[qp_closerdown]<<endl; cout<<"xold+
"<<x<<"
3666 xnew+
"<<xnew[n]<<endl;
3670 if (x > (xmax - xmin) / 2. && x < xmax)
3672 ratiox = (xmax - xnew_down[qp_closernormdown]) /
3673 ((xmax - xolddown[qp_closernormdown]));
3674 if ((xmax - xolddown[qp_closernormdown]) == 0)
3678 // xnew[n] = xnew_down[qp_closerdown]
3679 // + (x-xolddown[qp_closerolddown])*ratiox;
3680 xnew[n] = x + abs(nxvert[qp_closernormdown]) *
3681 (xnew_down[qp_closerolddown] -
3682 xolddown[qp_closerolddown]) *
3684 ASSERTL0(x > xmin, " x value <xmin down second half
");
3685 ASSERTL0(x < xmax, " x value >xmax down second half
");
3687 else if (x > xmin && x <= (xmax - xmin) / 2.)
3689 ratiox = (xnew_down[qp_closernormdown] - xmin) /
3690 ((xolddown[qp_closernormdown] - xmin));
3691 if ((xolddown[qp_closernormdown] - xmin) == 0)
3695 // xnew[n] = xnew_down[qp_closerdown]
3696 // + (x-xolddown[qp_closerolddown])*ratiox;
3697 xnew[n] = x + abs(nxvert[qp_closernormdown]) *
3698 (xnew_down[qp_closerolddown] -
3699 xolddown[qp_closerolddown]) *
3701 ASSERTL0(x > xmin, " x value <xmin down first half
");
3702 ASSERTL0(x < xmax, " x value >xmax down first half
");
3706 cout << "xold
" << x << " xnew=
" << xnew[n] << endl;
3707 ASSERTL0(xnew[n] >= xmin, "newx < xmin
");
3708 ASSERTL0(xnew[n] <= xmax, "newx > xmax
");
3713void CheckSingularQuads(MultiRegions::ExpListSharedPtr Exp, Array<OneD, int> V1,
3714 Array<OneD, int> V2, Array<OneD, NekDouble> &xnew,
3715 Array<OneD, NekDouble> &ynew)
3717 const std::shared_ptr<LocalRegions::ExpansionVector> exp2D = Exp->GetExp();
3718 int nel = exp2D->size();
3719 LocalRegions::QuadExpSharedPtr locQuadExp;
3720 LocalRegions::TriExpSharedPtr locTriExp;
3721 SpatialDomains::Geometry1DSharedPtr SegGeom;
3723 NekDouble xV1, yV1, xV2, yV2;
3724 NekDouble slopebef = 0.0, slopenext = 0.0, slopenew = 0.0;
3725 Array<OneD, int> locEids(4);
3726 for (int i = 0; i < nel; i++)
3728 if ((locQuadExp = (*exp2D)[i]->as<LocalRegions::QuadExp>()))
3730 SegGeom = (locQuadExp->GetGeom2D())->GetEdge(0);
3731 idbef = SegGeom->GetGlobalID();
3732 if (xnew[V1[idbef]] <= xnew[V2[idbef]])
3734 xV1 = xnew[V1[idbef]];
3735 yV1 = ynew[V1[idbef]];
3736 xV2 = xnew[V2[idbef]];
3737 yV2 = ynew[V2[idbef]];
3738 slopebef = (yV2 - yV1) / (xV2 - xV1);
3742 xV1 = xnew[V2[idbef]];
3743 yV1 = ynew[V2[idbef]];
3744 xV2 = xnew[V1[idbef]];
3745 yV2 = ynew[V1[idbef]];
3746 slopebef = (yV2 - yV1) / (xV2 - xV1);
3748 // cout<<"00 V1 x=
"<<xnew[ V1[idbef] ]<<" y=
"<<ynew[ V1[idbef]
3749 // ]<<endl; cout<<"00 V2 x=
"<<xnew[ V2[idbef] ]<<" y=
"<<ynew[
3750 // V2[idbef] ]<<endl;
3751 for (int j = 1; j < locQuadExp->GetNtraces(); ++j)
3753 SegGeom = (locQuadExp->GetGeom2D())->GetEdge(j);
3754 idnext = SegGeom->GetGlobalID();
3755 // cout<<"id=
"<<idnext<<" locid=
"<<j<<endl;
3756 // cout<<" V1 x=
"<<xnew[ V1[idnext] ]<<" y=
"<<ynew[
3757 // V1[idnext] ]<<endl; cout<<" V2 x=
"<<xnew[ V2[idnext] ]<<"
3759 if (xV1 == xnew[V1[idnext]] && yV1 == ynew[V1[idnext]])
3761 xV1 = xnew[V1[idnext]];
3762 yV1 = ynew[V1[idnext]];
3763 xV2 = xnew[V2[idnext]];
3764 yV2 = ynew[V2[idnext]];
3765 slopenext = (yV2 - yV1) / (xV2 - xV1);
3768 cout <<
"case1 x0=" << xV1 <<
" x1=" << xV2 << endl;
3769 cout << idnext <<
" 11slope bef =" << slopebef
3770 <<
" slopenext=" << slopenext << endl;
3773 if (slopebef / slopenext > 0.84 &&
3774 slopebef / slopenext < 1.18)
3776 xnew[V1[idnext]] = xnew[V1[idnext]] - 0.01;
3777 slopenew = (yV2 - yV1) / (xV2 - xnew[V1[idnext]]);
3779 if (
abs(slopebef - slopenew) <
3780 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;
3790 else if (xV2 == xnew[V2[idnext]] && yV2 == ynew[V2[idnext]])
3792 xV1 = xnew[V2[idnext]];
3793 yV1 = ynew[V2[idnext]];
3794 xV2 = xnew[V1[idnext]];
3795 yV2 = ynew[V1[idnext]];
3796 slopenext = (yV2 - yV1) / (xV2 - xV1);
3799 cout <<
"case2 x0=" << xV1 <<
" x1=" << xV2 << endl;
3800 cout << idnext <<
" 22slope bef =" << slopebef
3801 <<
" slopenext=" << slopenext << endl;
3804 if (slopebef / slopenext > 0.84 &&
3805 slopebef / slopenext < 1.18)
3807 xnew[V2[idnext]] = xnew[V2[idnext]] - 0.01;
3808 slopenew = (yV2 - yV1) / (xV2 - xnew[V2[idnext]]);
3810 if (
abs(slopebef - slopenew) <
3811 abs(slopebef - slopenext))
3813 xnew[V2[idnext]] = xnew[V2[idnext]] + 0.02;
3814 slopenew = (yV2 - yV1) / (xV2 - xnew[V2[idnext]]);
3816 slopenext = slopenew;
3817 cout <<
"slopenew=" << slopenew << endl;
3818 cout <<
"moved x=" << xnew[V2[idnext]] << endl;
3821 else if (xV1 == xnew[V2[idnext]] && yV1 == ynew[V2[idnext]])
3823 xV1 = xnew[V2[idnext]];
3824 yV1 = ynew[V2[idnext]];
3825 xV2 = xnew[V1[idnext]];
3826 yV2 = ynew[V1[idnext]];
3827 slopenext = (yV2 - yV1) / (xV2 - xV1);
3830 cout <<
"case3 x0=" << xV1 <<
" x1=" << xV2 << endl;
3831 cout << idnext <<
" 22slope bef =" << slopebef
3832 <<
" slopenext=" << slopenext << endl;
3835 if (slopebef / slopenext > 0.84 &&
3836 slopebef / slopenext < 1.18)
3838 xnew[V2[idnext]] = xnew[V2[idnext]] - 0.01;
3839 slopenew = (yV2 - yV1) / (xV2 - xnew[V2[idnext]]);
3841 if (
abs(slopebef - slopenew) <
3842 abs(slopebef - slopenext))
3844 xnew[V2[idnext]] = xnew[V2[idnext]] + 0.02;
3845 slopenew = (yV2 - yV1) / (xV2 - xnew[V2[idnext]]);
3847 slopenext = slopenew;
3848 cout <<
"slopenew=" << slopenew << endl;
3849 cout <<
"moved x=" << xnew[V2[idnext]] << endl;
3852 else if (xV2 == xnew[V1[idnext]] && yV2 == ynew[V1[idnext]])
3854 xV1 = xnew[V1[idnext]];
3855 yV1 = ynew[V1[idnext]];
3856 xV2 = xnew[V2[idnext]];
3857 yV2 = ynew[V2[idnext]];
3858 slopenext = (yV2 - yV1) / (xV2 - xV1);
3861 cout <<
"case4 x0=" << xV1 <<
" x1=" << xV2 << endl;
3862 cout << idnext <<
" 22slope bef =" << slopebef
3863 <<
" slopenext=" << slopenext << endl;
3866 if (slopebef / slopenext > 0.84 &&
3867 slopebef / slopenext < 1.18)
3869 xnew[V1[idnext]] = xnew[V1[idnext]] - 0.01;
3870 slopenew = (yV2 - yV1) / (xV2 - xnew[V1[idnext]]);
3872 if (
abs(slopebef - slopenew) <
3873 abs(slopebef - slopenext))
3875 xnew[V1[idnext]] = xnew[V1[idnext]] + 0.02;
3876 slopenew = (yV2 - yV1) / (xV2 - xnew[V1[idnext]]);
3878 slopenext = slopenew;
3879 cout <<
"slopenew=" << slopenew << endl;
3880 cout <<
"moved x=" << xnew[V1[idnext]] << endl;
3885 ASSERTL0(
false,
"edge not connected");
3887 slopebef = slopenext;
3896 int Npoints,
string s_alp,
3903 TiXmlDocument doc(filename);
3905 TiXmlElement *master = doc.FirstChildElement(
"NEKTAR");
3906 TiXmlElement *mesh = master->FirstChildElement(
"GEOMETRY");
3907 TiXmlElement *element = mesh->FirstChildElement(
"VERTEX");
3910 const char *xscal = element->Attribute(
"XSCALE");
3913 std::string xscalstr = xscal;
3915 xscale = expEvaluator.
Evaluate(expr_id);
3919 newfile = filename.substr(0, filename.find_last_of(
".")) +
"_moved.xml";
3920 doc.SaveFile(newfile);
3923 TiXmlDocument docnew(newfile);
3924 bool loadOkaynew = docnew.LoadFile();
3926 std::string errstr =
"Unable to load file: ";
3928 ASSERTL0(loadOkaynew, errstr.c_str());
3930 TiXmlHandle docHandlenew(&docnew);
3931 TiXmlElement *meshnew =
nullptr;
3932 TiXmlElement *masternew =
nullptr;
3933 TiXmlElement *condnew =
nullptr;
3934 TiXmlElement *Parsnew =
nullptr;
3935 TiXmlElement *parnew =
nullptr;
3939 masternew = docnew.FirstChildElement(
"NEKTAR");
3940 ASSERTL0(masternew,
"Unable to find NEKTAR tag in file.");
3944 condnew = masternew->FirstChildElement(
"CONDITIONS");
3945 Parsnew = condnew->FirstChildElement(
"PARAMETERS");
3946 cout <<
"alpha=" << s_alp << endl;
3947 parnew = Parsnew->FirstChildElement(
"P");
3950 TiXmlNode *node = parnew->FirstChild();
3954 std::string line = node->ToText()->Value();
3958 int beg = line.find_first_not_of(
" ");
3959 int end = line.find_first_of(
"=");
3966 if (end != line.find_last_of(
"="))
3971 if (end == std::string::npos)
3975 lhs = line.substr(line.find_first_not_of(
" "), end - beg);
3976 lhs = lhs.substr(0, lhs.find_last_not_of(
" ") + 1);
3982 boost::to_upper(lhs);
3985 alphastring =
"Alpha = " + s_alp;
3986 parnew->RemoveChild(node);
3987 parnew->LinkEndChild(
new TiXmlText(alphastring));
3991 parnew = parnew->NextSiblingElement(
"P");
3995 meshnew = masternew->FirstChildElement(
"GEOMETRY");
3997 ASSERTL0(meshnew,
"Unable to find GEOMETRY tag in file.");
3999 TiXmlElement *elementnew = meshnew->FirstChildElement(
"VERTEX");
4000 ASSERTL0(elementnew,
"Unable to find mesh VERTEX tag in file.");
4004 elementnew->SetAttribute(
"XSCALE", 1.0);
4006 TiXmlElement *vertexnew = elementnew->FirstChildElement(
"V");
4010 int nextVertexNumber = -1;
4016 TiXmlAttribute *vertexAttr = vertexnew->FirstAttribute();
4017 std::string attrName(vertexAttr->Name());
4019 (std::string(
"Unknown attribute name: ") + attrName).c_str());
4021 err = vertexAttr->QueryIntValue(&indx);
4022 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read attribute ID.");
4024 "Vertex IDs must begin with zero and be sequential.");
4026 std::string vertexBodyStr;
4028 TiXmlNode *vertexBody = vertexnew->FirstChild();
4030 if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
4032 vertexBodyStr += vertexBody->ToText()->Value();
4033 vertexBodyStr +=
" ";
4036 "Vertex definitions must contain vertex data.");
4038 vertexnew->RemoveChild(vertexBody);
4045 s << std::scientific << std::setprecision(8) << newx[nextVertexNumber]
4046 <<
" " << newy[nextVertexNumber] <<
" " << 0.0;
4047 vertexnew->LinkEndChild(
new TiXmlText(s.str()));
4052 vertexnew = vertexnew->NextSiblingElement(
"V");
4056 TiXmlElement *curvednew = meshnew->FirstChildElement(
"CURVED");
4057 ASSERTL0(curvednew,
"Unable to find mesh CURVED tag in file.");
4058 TiXmlElement *edgenew = curvednew->FirstChildElement(
"E");
4061 std::string charindex;
4065 int neids_lay = lay_eids[0].size();
4072 TiXmlAttribute *edgeAttr = edgenew->FirstAttribute();
4073 std::string attrName(edgeAttr->Name());
4074 charindex = edgeAttr->Value();
4075 std::istringstream iss(charindex);
4076 iss >> std::dec >> index;
4078 edgenew->QueryIntAttribute(
"EDGEID", &eid);
4081 for (
int u = 0; u < Eids.size(); u++)
4091 curvednew->RemoveChild(edgenew);
4097 std::string edgeBodyStr;
4099 TiXmlNode *edgeBody = edgenew->FirstChild();
4100 if (edgeBody->Type() == TiXmlNode::TINYXML_TEXT)
4102 edgeBodyStr += edgeBody->ToText()->Value();
4106 "Edge definitions must contain edge data");
4108 edgenew->RemoveChild(edgeBody);
4115 err = edgenew->QueryIntAttribute(
"NUMPOINTS", &numPts);
4117 "Unable to read curve attribute NUMPOINTS.");
4120 edgenew->SetAttribute(
"NUMPOINTS", Npoints);
4121 for (
int u = 0; u < Npoints; u++)
4123 st << std::scientific << std::setprecision(8)
4124 << xcPhys[cnt * Npoints + u] <<
" "
4125 << ycPhys[cnt * Npoints + u] <<
" " << 0.000 <<
" ";
4128 edgenew->LinkEndChild(
new TiXmlText(st.str()));
4149 edgenew = edgenew->NextSiblingElement(
"E");
4153 if (curv_lay ==
true)
4155 cout <<
"write other curved edges" << endl;
4156 TiXmlElement *curved = meshnew->FirstChildElement(
"CURVED");
4158 int nlays = lay_eids.size();
4163 for (
int g = 0; g < nlays; ++g)
4165 for (
int p = 0;
p < neids_lay;
p++)
4168 TiXmlElement *e =
new TiXmlElement(
"E");
4169 e->SetAttribute(
"ID", idcnt++);
4170 e->SetAttribute(
"EDGEID", lay_eids[g][
p]);
4171 e->SetAttribute(
"NUMPOINTS", Npoints);
4172 e->SetAttribute(
"TYPE",
"PolyEvenlySpaced");
4173 for (
int c = 0; c < Npoints; c++)
4175 st << std::scientific << std::setprecision(8)
4176 << x_lay[g][
p * Npoints + c] <<
" "
4177 << y_lay[g][
p * Npoints + c] <<
" " << 0.000 <<
" ";
4180 TiXmlText *t0 =
new TiXmlText(st.str());
4181 e->LinkEndChild(t0);
4182 curved->LinkEndChild(e);
4187 docnew.SaveFile(newfile);
4189 cout <<
"new file: " << newfile << endl;
#define ASSERTL0(condition, msg)
void GenerateMapEidsv1v2(MultiRegions::ExpListSharedPtr field, Array< OneD, int > &V1, Array< OneD, int > &V2)
int main(int argc, char *argv[])
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 EvaluateTangent(int npoints, Array< OneD, NekDouble > xcQedge, Array< OneD, NekDouble > coeffsinterp, Array< OneD, NekDouble > &txQedge, Array< OneD, NekDouble > &tyQedge)
NekDouble LagrangeInterpolant(NekDouble x, int npts, Array< OneD, NekDouble > xpts, Array< OneD, NekDouble > funcvals)
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 Cutrepetitions(int nedges, Array< OneD, NekDouble > inarray, Array< OneD, NekDouble > &outarray)
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 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 CheckSingularQuads(MultiRegions::ExpListSharedPtr Exp, Array< OneD, int > V1, Array< OneD, int > V2, 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)
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 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)
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 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)
void GenerateAddPointsNewtonIt(NekDouble xi, NekDouble yi, NekDouble &xout, NekDouble &yout, MultiRegions::ExpListSharedPtr function, Array< OneD, NekDouble > derfunction, NekDouble cr)
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 GenerateNeighbourArrays(int index, int neighpoints, Array< OneD, NekDouble > xArray, Array< OneD, NekDouble > yArray, Array< OneD, NekDouble > &Neighbour_x, Array< OneD, NekDouble > &Neighbour_y)
bool checkcommonvert(Array< OneD, int > Vids_laybefore, Array< OneD, int > Vids_c, int Vid)
int DetermineclosePointxindex(NekDouble x, Array< OneD, NekDouble > xArray)
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)
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)
Interpreter class for the evaluation of mathematical expressions.
int DefineFunction(const std::string &vlist, const std::string &expr)
Defines a function for the purposes of evaluation.
NekDouble Evaluate(const int id)
Evaluate a function which depends only on constants and/or parameters.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
const BoundaryRegionCollection & GetBoundaryRegions(void) const
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
static void Dgetrf(const int &m, const int &n, double *a, const int &lda, int *ipiv, int &info)
General matrix LU factorisation.
static void Dgetrs(const char &trans, const int &n, const int &nrhs, const double *a, const int &lda, int *ipiv, double *b, const int &ldb, int &info)
General matrix LU backsolve.
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...
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< QuadExp > QuadExpSharedPtr
std::shared_ptr< SegExp > SegExpSharedPtr
std::shared_ptr< TriExp > TriExpSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< ContField > ContFieldSharedPtr
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
std::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::shared_ptr< PointGeom > PointGeomSharedPtr
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
std::vector< double > w(NPUPPER)
std::vector< double > z(NPUPPER)
std::vector< double > d(NPUPPER *NPUPPER)
std::vector< double > q(NPUPPER *NPUPPER)
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
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 Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
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 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 Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/x.
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.
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
void Zero(int n, T *x, const int incx)
Zero vector.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > abs(scalarT< T > in)
scalarT< T > sqrt(scalarT< T > in)