Nektar++
Functions
MeshMove.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
#include <LibUtilities/LinearAlgebra/Lapack.hpp>
#include <LocalRegions/QuadExp.h>
#include <LocalRegions/SegExp.h>
#include <LocalRegions/TriExp.h>
#include <MultiRegions/ContField.h>
#include <MultiRegions/ExpList.h>
#include <boost/lexical_cast.hpp>
#include <tinyxml.h>

Go to the source code of this file.

Functions

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 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 GenerateAddPointsNewtonIt (NekDouble xi, NekDouble yi, NekDouble &xout, NekDouble &yout, MultiRegions::ExpListSharedPtr function, Array< OneD, NekDouble > derfunction, NekDouble cr)
 
void GenerateMapEidsv1v2 (MultiRegions::ExpListSharedPtr field, Array< OneD, int > &V1, Array< OneD, int > &V2)
 
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)
 
bool checkcommonvert (Array< OneD, int > Vids_laybefore, Array< OneD, int > Vids_c, int Vid)
 
void Cutrepetitions (int nedges, Array< OneD, NekDouble > inarray, Array< OneD, NekDouble > &outarray)
 
int DetermineclosePointxindex (NekDouble x, Array< OneD, NekDouble > xArray)
 
void GenerateNeighbourArrays (int index, int neighpoints, Array< OneD, NekDouble > xArray, Array< OneD, NekDouble > yArray, Array< OneD, NekDouble > &Neighbour_x, Array< OneD, NekDouble > &Neighbour_y)
 
NekDouble LagrangeInterpolant (NekDouble x, int npts, Array< OneD, NekDouble > xpts, Array< OneD, NekDouble > funcvals)
 
void EvaluateTangent (int npoints, Array< OneD, NekDouble > xcQedge, Array< OneD, NekDouble > coeffsinterp, Array< OneD, NekDouble > &txQedge, Array< OneD, NekDouble > &tyQedge)
 
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 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 Orderfunctionx (Array< OneD, NekDouble > inarray_x, Array< OneD, NekDouble > inarray_y, Array< OneD, NekDouble > &outarray_x, Array< OneD, NekDouble > &outarray_y)
 
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 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 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 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 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 CheckSingularQuads (MultiRegions::ExpListSharedPtr Exp, Array< OneD, int > V1, Array< OneD, int > V2, 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)
 
int main (int argc, char *argv[])
 

Function Documentation

◆ checkcommonvert()

bool checkcommonvert ( Array< OneD, int >  Vids_laybefore,
Array< OneD, int >  Vids_c,
int  Vid 
)

Definition at line 2753 of file MeshMove.cpp.

2755{
2756 bool check = false;
2757 for (int u = 0; u < Vids_laybefore.size(); u++)
2758 {
2759 if (Vids_laybefore[u] == Vid || Vids_c[u] == Vid)
2760 {
2761 check = true;
2762 }
2763 cout << Vid << " Vert test=" << Vids_laybefore[u] << endl;
2764 }
2765 return check;
2766}

Referenced by MappingEVids().

◆ CheckSingularQuads()

void CheckSingularQuads ( MultiRegions::ExpListSharedPtr  Exp,
Array< OneD, int >  V1,
Array< OneD, int >  V2,
Array< OneD, NekDouble > &  xnew,
Array< OneD, NekDouble > &  ynew 
)

Definition at line 3713 of file MeshMove.cpp.

3716{
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;
3722 int idbef, idnext;
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++)
3727 {
3728 if ((locQuadExp = (*exp2D)[i]->as<LocalRegions::QuadExp>()))
3729 {
3730 SegGeom = (locQuadExp->GetGeom2D())->GetEdge(0);
3731 idbef = SegGeom->GetGlobalID();
3732 if (xnew[V1[idbef]] <= xnew[V2[idbef]])
3733 {
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);
3739 }
3740 else
3741 {
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);
3747 }
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)
3752 {
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] ]<<"
3758 // y="<<ynew[ V2[idnext] ]<<endl;
3759 if (xV1 == xnew[V1[idnext]] && yV1 == ynew[V1[idnext]])
3760 {
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);
3766 if (i == 23)
3767 {
3768 cout << "case1 x0=" << xV1 << " x1=" << xV2 << endl;
3769 cout << idnext << " 11slope bef =" << slopebef
3770 << " slopenext=" << slopenext << endl;
3771 }
3772 // compare with slope before
3773 if (slopebef / slopenext > 0.84 &&
3774 slopebef / slopenext < 1.18)
3775 {
3776 xnew[V1[idnext]] = xnew[V1[idnext]] - 0.01;
3777 slopenew = (yV2 - yV1) / (xV2 - xnew[V1[idnext]]);
3778
3779 if (abs(slopebef - slopenew) <
3780 abs(slopebef - slopenext))
3781 {
3782 xnew[V1[idnext]] = xnew[V1[idnext]] + 0.02;
3783 slopenew = (yV2 - yV1) / (xV2 - xnew[V1[idnext]]);
3784 }
3785 slopenext = slopenew;
3786 cout << "slopenew=" << slopenew << endl;
3787 cout << "moved x=" << xnew[V1[idnext]] << endl;
3788 }
3789 }
3790 else if (xV2 == xnew[V2[idnext]] && yV2 == ynew[V2[idnext]])
3791 {
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);
3797 if (i == 23)
3798 {
3799 cout << "case2 x0=" << xV1 << " x1=" << xV2 << endl;
3800 cout << idnext << " 22slope bef =" << slopebef
3801 << " slopenext=" << slopenext << endl;
3802 }
3803 // compare with slope before
3804 if (slopebef / slopenext > 0.84 &&
3805 slopebef / slopenext < 1.18)
3806 {
3807 xnew[V2[idnext]] = xnew[V2[idnext]] - 0.01;
3808 slopenew = (yV2 - yV1) / (xV2 - xnew[V2[idnext]]);
3809
3810 if (abs(slopebef - slopenew) <
3811 abs(slopebef - slopenext))
3812 {
3813 xnew[V2[idnext]] = xnew[V2[idnext]] + 0.02;
3814 slopenew = (yV2 - yV1) / (xV2 - xnew[V2[idnext]]);
3815 }
3816 slopenext = slopenew;
3817 cout << "slopenew=" << slopenew << endl;
3818 cout << "moved x=" << xnew[V2[idnext]] << endl;
3819 }
3820 }
3821 else if (xV1 == xnew[V2[idnext]] && yV1 == ynew[V2[idnext]])
3822 {
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);
3828 if (i == 23)
3829 {
3830 cout << "case3 x0=" << xV1 << " x1=" << xV2 << endl;
3831 cout << idnext << " 22slope bef =" << slopebef
3832 << " slopenext=" << slopenext << endl;
3833 }
3834 // compare with slope before
3835 if (slopebef / slopenext > 0.84 &&
3836 slopebef / slopenext < 1.18)
3837 {
3838 xnew[V2[idnext]] = xnew[V2[idnext]] - 0.01;
3839 slopenew = (yV2 - yV1) / (xV2 - xnew[V2[idnext]]);
3840
3841 if (abs(slopebef - slopenew) <
3842 abs(slopebef - slopenext))
3843 {
3844 xnew[V2[idnext]] = xnew[V2[idnext]] + 0.02;
3845 slopenew = (yV2 - yV1) / (xV2 - xnew[V2[idnext]]);
3846 }
3847 slopenext = slopenew;
3848 cout << "slopenew=" << slopenew << endl;
3849 cout << "moved x=" << xnew[V2[idnext]] << endl;
3850 }
3851 }
3852 else if (xV2 == xnew[V1[idnext]] && yV2 == ynew[V1[idnext]])
3853 {
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);
3859 if (i == 23)
3860 {
3861 cout << "case4 x0=" << xV1 << " x1=" << xV2 << endl;
3862 cout << idnext << " 22slope bef =" << slopebef
3863 << " slopenext=" << slopenext << endl;
3864 }
3865 // compare with slope before
3866 if (slopebef / slopenext > 0.84 &&
3867 slopebef / slopenext < 1.18)
3868 {
3869 xnew[V1[idnext]] = xnew[V1[idnext]] - 0.01;
3870 slopenew = (yV2 - yV1) / (xV2 - xnew[V1[idnext]]);
3871
3872 if (abs(slopebef - slopenew) <
3873 abs(slopebef - slopenext))
3874 {
3875 xnew[V1[idnext]] = xnew[V1[idnext]] + 0.02;
3876 slopenew = (yV2 - yV1) / (xV2 - xnew[V1[idnext]]);
3877 }
3878 slopenext = slopenew;
3879 cout << "slopenew=" << slopenew << endl;
3880 cout << "moved x=" << xnew[V1[idnext]] << endl;
3881 }
3882 }
3883 else
3884 {
3885 ASSERTL0(false, "edge not connected");
3886 }
3887 slopebef = slopenext;
3888 }
3889 }
3890 }
3891}
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
Nektar::LocalRegions::QuadExpSharedPtr
std::shared_ptr< QuadExp > QuadExpSharedPtr
Definition: QuadExp.h:255
Nektar::LocalRegions::TriExpSharedPtr
std::shared_ptr< TriExp > TriExpSharedPtr
Definition: TriExp.h:250
Nektar::SpatialDomains::Geometry1DSharedPtr
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry.h:61
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:298

References tinysimd::abs().

Referenced by main().

◆ Computestreakpositions()

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 
)

Definition at line 2102 of file MeshMove.cpp.

2112{
2113 cout << "Computestreakpositions" << endl;
2114 int nq = streak->GetTotPoints();
2115 Array<OneD, NekDouble> coord(2);
2116 // Array<OneD, NekDouble> stvalues(nvertl,-10);
2117 Array<OneD, NekDouble> derstreak(nq);
2118 streak->PhysDeriv(MultiRegions::eY, streak->GetPhys(), derstreak);
2119 int elmtid, offset;
2120 NekDouble U = 0.0, dU = 0.0;
2121 NekDouble F = 1000;
2122
2123 if (verts == true) // only for verts makes sense to init the coord values..
2124 {
2125
2126 // start guess
2127 // yc= (yup+ydown)/2
2128 Vmath::Vadd(xc.size(), yold_up, 1, yold_low, 1, yc, 1);
2129 Vmath::Smul(xc.size(), 0.5, yc, 1, yc, 1);
2130 Vmath::Vcopy(xc.size(), xold_c, 1, xc, 1);
2131 }
2132 else // case of xQ,yQ
2133 {
2134 Vmath::Vcopy(xc.size(), xold_c, 1, xc, 1);
2135 Vmath::Vcopy(xc.size(), yold_c, 1, yc, 1);
2136 }
2137
2138 int its;
2139 int attempt;
2140 NekDouble tol = 1e-3;
2141 NekDouble ytmp;
2142 for (int e = 0; e < npoints; e++)
2143 {
2144 coord[0] = xc[e];
2145 coord[1] = yc[e];
2146
2147 elmtid = streak->GetExpIndex(coord, 0.00001);
2148 offset = streak->GetPhys_Offset(elmtid);
2149 F = 1000;
2150 its = 0;
2151 attempt = 0;
2152 ytmp = coord[1];
2153 // cout<<"start guess: x="<<xc[e]<<" y="<<yc[e]<<endl;
2154 while (abs(F) > 0.000000001)
2155 {
2156
2157 elmtid = streak->GetExpIndex(coord, 0.00001);
2158
2159 // stvalues[e] = streak->GetExp(elmtid)->PhysEvaluate(coord,
2160 // streak->GetPhys() +offset );
2161 // cout<<"elmtid="<<elmtid<<" x="<<coord[0]<<" y="<<coord[1]<<"
2162 // stvalue="<<U<<" dU="<<dU<<endl;
2163
2164 if ((abs(coord[1]) > 1 || elmtid == -1) && attempt == 0 &&
2165 verts == true)
2166 {
2167 // try the old crit lay position:
2168 coord[1] = yold_c[e];
2169 attempt++;
2170 }
2171 else if ((abs(coord[1]) > 1 || elmtid == -1))
2172 {
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))
2182 {
2183 coord[1] = ytmp - 0.01;
2184 }
2185
2186 attempt++;
2187 }
2188 else
2189 {
2190 ASSERTL0(abs(coord[1]) <= 1, " y value out of bound +/-1");
2191 }
2192 offset = streak->GetPhys_Offset(elmtid);
2193 U = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() +
2194 offset);
2195 dU =
2196 streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
2197 coord[1] = coord[1] - (U - cr) / dU;
2198 F = U - cr;
2199 ASSERTL0(coord[0] == xc[e], " x coordinate must remain the same");
2200
2201 its++;
2202 if (its > 200 && abs(F) < 0.00001)
2203 {
2204 cout << "warning streak position obtained with precision:" << F
2205 << endl;
2206 break;
2207 }
2208 else if (its > 1000 && abs(F) < 0.0001)
2209 {
2210 cout << "warning streak position obtained with precision:" << F
2211 << endl;
2212 break;
2213 }
2214 else if (its > 1000)
2215 {
2216 ASSERTL0(false, "no convergence after 1000 iterations");
2217 }
2218 }
2219 yc[e] = coord[1] - (U - cr) / dU;
2220 ASSERTL0(U <= cr + tol, "streak wrong+");
2221 ASSERTL0(U >= cr - tol, "streak wrong-");
2222 // Utilities::Zerofunction(coord[0], coord[1], xtest, ytest, streak,
2223 // derstreak);
2224 cout << "result streakvert x=" << xc[e] << " y=" << yc[e]
2225 << " streak=" << U << endl;
2226 }
2227}
Vmath::Vadd
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.
Definition: Vmath.hpp:180
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825

References tinysimd::abs(), ASSERTL0, Nektar::MultiRegions::eY, Vmath::Smul(), Vmath::Vadd(), and Vmath::Vcopy().

Referenced by main().

◆ Cutrepetitions()

void Cutrepetitions ( int  nedges,
Array< OneD, NekDouble >  inarray,
Array< OneD, NekDouble > &  outarray 
)

Definition at line 2768 of file MeshMove.cpp.

2770{
2771
2772 // determine npedge:
2773 int np_lay = inarray.size();
2774 ASSERTL0(inarray.size() % nedges == 0, " something on number npedge");
2775 // cut out the repetitions:
2776
2777 int cnt = 0;
2778 for (int w = 0; w < np_lay; w++)
2779 {
2780
2781 if (w < np_lay - 1)
2782 {
2783 if (inarray[w] == inarray[w + 1])
2784 {
2785 continue;
2786 }
2787 }
2788 outarray[cnt] = inarray[w];
2789 cnt++;
2790 }
2791
2792 ASSERTL0(cnt == np_lay - (nedges - 1), "wrong cut");
2793}
Nektar::UnitTests::w
std::vector< double > w(NPUPPER)

References ASSERTL0, and Nektar::UnitTests::w().

Referenced by main(), MoveLayerNfixedxpos(), and MoveLayerNnormpos().

◆ DetermineclosePointxindex()

int DetermineclosePointxindex ( NekDouble  x,
Array< OneD, NekDouble >  xArray 
)

Definition at line 2795 of file MeshMove.cpp.

2796{
2797 int npts = xArray.size();
2798 Array<OneD, NekDouble> xcopy(npts);
2799 Vmath::Vcopy(npts, xArray, 1, xcopy, 1);
2800 // subtract xpoint and abs
2801
2802 Vmath::Sadd(npts, -x, xcopy, 1, xcopy, 1);
2803 Vmath::Vabs(npts, xcopy, 1, xcopy, 1);
2804
2805 int index = Vmath::Imin(npts, xcopy, 1);
2806 if (xArray[index] > x) // assume always x[index]< x(HYPHOTHESIS)
2807 {
2808 index = index - 1;
2809 }
2810 return index;
2811}
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.hpp:352
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.hpp:704
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.hpp:194

References Vmath::Imin(), Vmath::Sadd(), Vmath::Vabs(), and Vmath::Vcopy().

Referenced by main(), MoveLayerNfixedxpos(), and MoveLayerNnormpos().

◆ EvaluateTangent()

void EvaluateTangent ( int  npoints,
Array< OneD, NekDouble >  xcQedge,
Array< OneD, NekDouble >  coeffsinterp,
Array< OneD, NekDouble > &  txQedge,
Array< OneD, NekDouble > &  tyQedge 
)

Definition at line 2893 of file MeshMove.cpp.

2897{
2898 Array<OneD, NekDouble> yprime(npoints, 0.0);
2899 int np_pol = coeffsinterp.size();
2900 cout << "evaluatetan with " << np_pol << endl;
2901
2902 // calc derivative poly (cut last entry on b array that is the const)
2903 int derorder;
2904 Array<OneD, NekDouble> yinterp(npoints, 0.0);
2905 int polorder;
2906 for (int q = 0; q < npoints; q++)
2907 {
2908 polorder = np_pol - 1;
2909 derorder = np_pol - 2;
2910 yprime[q] = 0;
2911 for (int d = 0; d < np_pol - 1; d++)
2912 {
2913 yprime[q] += (derorder + 1) * coeffsinterp[d] *
2914 std::pow(xcQedge[q], derorder);
2915 derorder--;
2916 }
2917 // test
2918 for (int a = 0; a < np_pol; a++)
2919 {
2920 yinterp[q] += coeffsinterp[a] * std::pow(xcQedge[q], polorder);
2921 // cout<<"coeff*x^n="<<b[a]*std::pow(xcQedge[q],polorder)<<"
2922 // sum="<<yinterp[q]<<endl;
2923 polorder--;
2924 }
2925 }
2926
2927 // transf yprime into tx,ty:
2928 for (int n = 0; n < npoints; n++)
2929 {
2930 // ABS???!!
2931 txQedge[n] = 0;
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;
2936 }
2937}
Nektar::UnitTests::d
std::vector< double > d(NPUPPER *NPUPPER)
Nektar::UnitTests::q
std::vector< double > q(NPUPPER *NPUPPER)

References Nektar::UnitTests::d(), and Nektar::UnitTests::q().

Referenced by main().

◆ GenerateAddPointsNewtonIt()

void GenerateAddPointsNewtonIt ( NekDouble  xi,
NekDouble  yi,
NekDouble &  xout,
NekDouble &  yout,
MultiRegions::ExpListSharedPtr  function,
Array< OneD, NekDouble >  derfunction,
NekDouble  cr 
)

Definition at line 2229 of file MeshMove.cpp.

2233{
2234 int elmtid, offset;
2235 NekDouble F = 1000;
2236 NekDouble U = 0.0;
2237 NekDouble dU = 0.0;
2238
2239 Array<OneD, NekDouble> coords(2);
2240 coords[0] = xi;
2241 coords[1] = yi;
2242
2243 int attempt = 0;
2244 int its = 0;
2245 NekDouble ytmp = coords[1];
2246 while (abs(F) > 0.00000001)
2247 {
2248
2249 // cout<<"generate newton it xi="<<xi<<" yi="<<yi<<endl;
2250 elmtid = function->GetExpIndex(coords, 0.01);
2251 //@to do if GetType(elmtid)==triangular WRONG!!!
2252 cout << "gen newton xi=" << xi << " yi=" << coords[1]
2253 << " elmtid=" << elmtid << " F=" << F << endl;
2254
2255 if ((abs(coords[1]) > 1 || elmtid == -1))
2256 {
2257
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))
2268 {
2269 coords[1] = ytmp - 0.01;
2270 }
2271
2272 attempt++;
2273 }
2274 else if (abs(coords[1]) < 1.01 && attempt == 0)
2275 {
2276 coords[1] = 1.0;
2277 attempt++;
2278 }
2279 else
2280 {
2281 ASSERTL0(abs(coords[1]) <= 1.00, " y value out of bound +/-1");
2282 }
2283 offset = function->GetPhys_Offset(elmtid);
2284 U = function->GetExp(elmtid)->PhysEvaluate(coords, function->GetPhys() +
2285 offset);
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;
2291 F = U - cr;
2292
2293 its++;
2294 if (its > 200 && abs(F) < 0.00001)
2295 {
2296 cout << "warning streak position obtained with precision:" << F
2297 << endl;
2298 break;
2299 }
2300 else if (its > 1000 && abs(F) < 0.0001)
2301 {
2302 cout << "warning streak position obtained with precision:" << F
2303 << endl;
2304 break;
2305 }
2306 else if (its > 1000)
2307 {
2308 ASSERTL0(false, "no convergence after 1000 iterations");
2309 }
2310
2311 ASSERTL0(coords[0] == xi, " x coordinate must remain the same");
2312 }
2313 xout = xi;
2314 yout = coords[1] - (U - cr) / dU;
2315 cout << "NewtonIt result x=" << xout << " y=" << coords[1] << " U=" << U
2316 << endl;
2317}

References tinysimd::abs(), and ASSERTL0.

Referenced by main().

◆ GenerateMapEidsv1v2()

void GenerateMapEidsv1v2 ( MultiRegions::ExpListSharedPtr  field,
Array< OneD, int > &  V1,
Array< OneD, int > &  V2 
)

Definition at line 2319 of file MeshMove.cpp.

2321{
2322
2323 const std::shared_ptr<LocalRegions::ExpansionVector> exp2D =
2324 field->GetExp();
2325 int nel = exp2D->size();
2326 LocalRegions::QuadExpSharedPtr locQuadExp;
2327 LocalRegions::TriExpSharedPtr locTriExp;
2328 SpatialDomains::Geometry1DSharedPtr SegGeom;
2329 int id;
2330 int cnt = 0;
2331 Array<OneD, int> V1tmp(4 * nel, 10000);
2332 Array<OneD, int> V2tmp(4 * nel, 10000);
2333 for (int i = 0; i < nel; i++)
2334 {
2335 if ((locQuadExp = (*exp2D)[i]->as<LocalRegions::QuadExp>()))
2336 {
2337 for (int j = 0; j < locQuadExp->GetNtraces(); ++j)
2338 {
2339 SegGeom = (locQuadExp->GetGeom2D())->GetEdge(j);
2340 id = SegGeom->GetGlobalID();
2341
2342 if (V1tmp[id] == 10000)
2343 {
2344 V1tmp[id] = SegGeom->GetVertex(0)->GetVid();
2345 V2tmp[id] = SegGeom->GetVertex(1)->GetVid();
2346 cnt++;
2347 }
2348 }
2349 }
2350 // in the future the tri edges may be not necessary (if the nedges is
2351 // known)
2352
2353 else if ((locTriExp = (*exp2D)[i]->as<LocalRegions::TriExp>()))
2354 {
2355 for (int j = 0; j < locTriExp->GetNtraces(); ++j)
2356 {
2357 SegGeom = (locTriExp->GetGeom2D())->GetEdge(j);
2358 id = SegGeom->GetGlobalID();
2359
2360 if (V1tmp[id] == 10000)
2361 {
2362 V1tmp[id] = SegGeom->GetVertex(0)->GetVid();
2363 V2tmp[id] = SegGeom->GetVertex(1)->GetVid();
2364 cnt++;
2365 }
2366 }
2367 }
2368 }
2369
2370 V1 = Array<OneD, int>(cnt);
2371 V2 = Array<OneD, int>(cnt);
2372 // populate the output vectors V1,V2
2373 for (int g = 0; g < cnt; g++)
2374 {
2375 V1[g] = V1tmp[g];
2376 ASSERTL0(V1tmp[g] != 10000, "V1 wrong");
2377 V2[g] = V2tmp[g];
2378 ASSERTL0(V2tmp[g] != 10000, "V2 wrong");
2379 }
2380}

References ASSERTL0.

Referenced by main().

◆ GenerateNeighbourArrays()

void GenerateNeighbourArrays ( int  index,
int  neighpoints,
Array< OneD, NekDouble >  xArray,
Array< OneD, NekDouble >  yArray,
Array< OneD, NekDouble > &  Neighbour_x,
Array< OneD, NekDouble > &  Neighbour_y 
)

Definition at line 2813 of file MeshMove.cpp.

2818{
2819 ASSERTL0(neighpoints % 2 == 0, "number of neighbour points should be even");
2820 int leftpoints = (neighpoints / 2) - 1; //-1 because xArray[index]< x
2821 int rightpoints = neighpoints / 2;
2822 int diff;
2823 int start;
2824 // cout<<"index="<<index<<" left="<<leftpoints<<"
2825 // right="<<rightpoints<<endl;
2826 if (index - leftpoints < 0)
2827 {
2828 // cout<"case0"<<endl;
2829 diff = index - leftpoints;
2830 start = 0;
2831 Vmath::Vcopy(neighpoints, &yArray[0], 1, &Neighbour_y[0], 1);
2832 Vmath::Vcopy(neighpoints, &xArray[0], 1, &Neighbour_x[0], 1);
2833 }
2834 else if ((yArray.size() - 1) - index < rightpoints)
2835 {
2836 // cout"case1 closest="<<xArray[index]<<endl;
2837 int rpoints = (yArray.size() - 1) - index; //
2838 diff = rightpoints - rpoints;
2839 // cout<"start index="<<index-leftpoints-diff<<endl;
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);
2843 }
2844 else
2845 {
2846 // cout<<"caseaa"<<endl;
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);
2850 }
2851 /*
2852 for(int t= start; t<start+neighpoints; t++)
2853 {
2854 cout<<"Px="<<xArray[t]<<" "<<yArray[t]<<endl;
2855 }
2856 */
2857 // check if there is any repetition
2858 for (int f = 1; f < neighpoints; f++)
2859 {
2860 ASSERTL0(Neighbour_x[f] != Neighbour_x[f - 1],
2861 " repetition on NeighbourArrays");
2862 }
2863}

References ASSERTL0, and Vmath::Vcopy().

Referenced by main(), MoveLayerNfixedxpos(), and MoveLayerNnormpos().

◆ LagrangeInterpolant()

NekDouble LagrangeInterpolant ( NekDouble  x,
int  npts,
Array< OneD, NekDouble >  xpts,
Array< OneD, NekDouble >  funcvals 
)

Definition at line 2865 of file MeshMove.cpp.

2868{
2869 NekDouble sum = 0.0;
2870 NekDouble LagrangePoly;
2871 // cout<<"lagrange"<<endl;
2872 for (int pt = 0; pt < npts; ++pt)
2873 {
2874 NekDouble h = 1.0;
2875
2876 for (int j = 0; j < pt; ++j)
2877 {
2878 h = h * (x - xpts[j]) / (xpts[pt] - xpts[j]);
2879 }
2880
2881 for (int k = pt + 1; k < npts; ++k)
2882 {
2883 h = h * (x - xpts[k]) / (xpts[pt] - xpts[k]);
2884 }
2885 LagrangePoly = h;
2886
2887 sum += funcvals[pt] * LagrangePoly;
2888 }
2889 // cout<<"result :"<<sum<<endl;
2890 return sum;
2891}

Referenced by main(), MoveLayerNfixedxpos(), and MoveLayerNnormpos().

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 175 of file MeshMove.cpp.

176{
177 NekDouble cr;
178 // set cr =0
179 cr = 0;
180 // change argc from 6 to 5 allow the loading of cr to be optional
181 if (argc > 6 || argc < 5)
182 {
183 fprintf(stderr, "Usage: ./MoveMesh meshfile fieldfile changefile "
184 "alpha cr(optional)\n");
185 exit(1);
186 }
187
188 // ATTEnTION !!! with argc=2 you impose that vSession refers to is
189 // argv[1]=meshfile!!!!!
190 LibUtilities::SessionReaderSharedPtr vSession =
191 LibUtilities::SessionReader::CreateInstance(2, argv);
192 SpatialDomains::MeshGraphSharedPtr graphShPt =
193 SpatialDomains::MeshGraph::Read(vSession);
194 //----------------------------------------------
195
196 if (argc == 6 && vSession->DefinesSolverInfo("INTERFACE") &&
197 vSession->GetSolverInfo("INTERFACE") == "phase")
198 {
199 cr = boost::lexical_cast<NekDouble>(argv[argc - 1]);
200 argc = 5;
201 }
202
203 // Also read and store the boundary conditions
204 SpatialDomains::BoundaryConditionsSharedPtr boundaryConditions;
205 boundaryConditions =
206 MemoryManager<SpatialDomains::BoundaryConditions>::AllocateSharedPtr(
207 vSession, graphShPt);
208 SpatialDomains::BoundaryConditions bcs(vSession, graphShPt);
209 //----------------------------------------------
210
211 // the mesh file should have 2 component: set output fields
212 // fields has to be of the SAME dimension of the mesh (that's why there is
213 // the changefile as an input)
214 // a contfield is needed to extract boundary conditions!!!
215
216 // store name of the file to change
217 string changefile(argv[argc - 2]);
218 //----------------------------------------------
219
220 // store the value of alpha
221 string charalp(argv[argc - 1]);
222 // NekDouble alpha = boost::lexical_cast<NekDouble>(charalp);
223 cout << "read alpha=" << charalp << endl;
224
225 //---------------------------------------------
226 // Import field file.
227 string fieldfile(argv[argc - 3]);
228 vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
229 vector<vector<NekDouble>> fielddata;
230 //----------------------------------------------
231
232 MultiRegions::ExpListSharedPtr streak;
233 streak = MemoryManager<MultiRegions::ContField>::AllocateSharedPtr(
234 vSession, graphShPt, "w", true);
235
236 LibUtilities::Import(fieldfile, fielddef, fielddata);
237
238 for (int i = 0; i < fielddata.size(); i++)
239 {
240 streak->ExtractDataToCoeffs(fielddef[i], fielddata[i],
241 fielddef[i]->m_fields[0],
242 streak->UpdateCoeffs());
243 }
244 streak->BwdTrans(streak->GetCoeffs(), streak->UpdatePhys());
245
246 //------------------------------------------------
247 // determine the I regions (3 region expected)
248 // hypothesys: the layes have the same number of points
249
250 int nIregions, lastIregion = 0;
251 const Array<OneD, SpatialDomains::BoundaryConditionShPtr> bndConditions =
252 streak->GetBndConditions();
253 Array<OneD, int> Iregions = Array<OneD, int>(bndConditions.size(), -1);
254
255 nIregions = 0;
256 int nbnd = bndConditions.size();
257 for (int r = 0; r < nbnd; r++)
258 {
259 if (bndConditions[r]->GetUserDefined() == "CalcBC")
260 {
261 lastIregion = r;
262 Iregions[r] = r;
263 nIregions++;
264 }
265 }
266
267 ASSERTL0(nIregions > 0,
268 "there is any boundary region with the tag USERDEFINEDTYPE="
269 "CalcBC"
270 " specified");
271 cout << "nIregions=" << nIregions << endl;
272 // set expansion along a layers
273 Array<OneD, MultiRegions::ExpListSharedPtr> bndfieldx =
274 streak->GetBndCondExpansions();
275 //--------------------------------------------------------
276 // determine the points in the lower and upper curve...
277 int nedges = bndfieldx[lastIregion]->GetExpSize();
278 int nvertl = nedges + 1;
279 Array<OneD, int> Vids_low(nvertl, -10);
280 Array<OneD, NekDouble> xold_low(nvertl);
281 Array<OneD, NekDouble> yold_low(nvertl);
282 Array<OneD, NekDouble> zi(nvertl);
283
284 // order the ids on the lower curve lastIregion starting from the id on x=0
285 NekDouble x_connect;
286 NekDouble x0, y0, z0, yt = 0, zt = 0;
287 int lastedge = -1;
288 int v1, v2;
289 // first point for x_connect=0(or-1.6 for the full mesh (-pi,pi) )
290 x_connect = 0;
291 SpatialDomains::PointGeomSharedPtr vertex0 = graphShPt->GetVertex(
292 ((bndfieldx[lastIregion]->GetExp(0)->as<LocalRegions::SegExp>())
293 ->GetGeom1D())
294 ->GetVid(0));
295 vertex0->GetCoords(x0, y0, z0);
296 if (x0 != 0.0)
297 {
298 cout << "WARNING x0=" << x0 << endl;
299 x_connect = x0;
300 }
301 v1 = 0;
302 v2 = 1;
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");
306 SpatialDomains::PointGeomSharedPtr vertex =
307 graphShPt->GetVertex(Vids_low[v2]);
308
309 // update x_connect
310 cout << "x_conn=" << x_connect << " yt=" << yt << " zt=" << zt
311 << " vid=" << Vids_low[v2] << endl;
312 vertex->GetCoords(x_connect, yt, zt);
313
314 int i = 2;
315 while (i < nvertl)
316 {
317 v1 = i;
318 OrderVertices(nedges, graphShPt, bndfieldx[lastIregion - 1], Vids_low,
319 v1, v2, x_connect, lastedge, xold_low, yold_low);
320 SpatialDomains::PointGeomSharedPtr vertex =
321 graphShPt->GetVertex(Vids_low[v1]);
322 // update x_connect (lastedge is updated on the OrderVertices function)
323 vertex->GetCoords(x_connect, yt, zt);
324 i++;
325 }
326
327 //------------------------------------------------------------------------
328 // order in the same way the id of the upper curve lastIregion-1 starting
329 // from x=0:
330 Array<OneD, int> Vids_up(nvertl, -10);
331 Array<OneD, NekDouble> xold_up(nvertl);
332 Array<OneD, NekDouble> yold_up(nvertl);
333 // first point for x_connect=0 (or-1.6)
334 x_connect = 0;
335 vertex0 = graphShPt->GetVertex(
336 ((bndfieldx[lastIregion]->GetExp(0)->as<LocalRegions::SegExp>())
337 ->GetGeom1D())
338 ->GetVid(0));
339 vertex0->GetCoords(x0, y0, z0);
340 if (x0 != 0.0)
341 {
342 cout << "WARNING x0=" << x0 << endl;
343 x_connect = x0;
344 }
345 lastedge = -1;
346
347 v1 = 0;
348 v2 = 1;
349 OrderVertices(nedges, graphShPt, bndfieldx[lastIregion - 2], Vids_up, v1,
350 v2, x_connect, lastedge, xold_up, yold_up);
351 SpatialDomains::PointGeomSharedPtr vertexU =
352 graphShPt->GetVertex(Vids_up[v2]);
353 vertexU->GetCoords(x_connect, yt, zt);
354
355 i = 2;
356 while (i < nvertl)
357 {
358 v1 = i;
359 OrderVertices(nedges, graphShPt, bndfieldx[lastIregion - 2], Vids_up,
360 v1, v2, x_connect, lastedge, xold_up, yold_up);
361 SpatialDomains::PointGeomSharedPtr vertex =
362 graphShPt->GetVertex(Vids_up[v1]);
363 // cout<<"VIdup="<<Vids_up[v1]<<endl;
364 // update x_connect (lastedge is updated on the OrderVertices function)
365 vertex->GetCoords(x_connect, yt, zt);
366 i++;
367 }
368
369 //-----------------------------------------------------------------------------------
370 // order in the same way the id of the layer curve lastIregion starting from
371 // x=0:
372 Array<OneD, int> Vids_c(nvertl, -10);
373 Array<OneD, NekDouble> xold_c(nvertl);
374 Array<OneD, NekDouble> yold_c(nvertl);
375 // first point for x_connect=0(or-1.6)
376 x_connect = 0;
377 vertex0 = graphShPt->GetVertex(
378 ((bndfieldx[lastIregion]->GetExp(0)->as<LocalRegions::SegExp>())
379 ->GetGeom1D())
380 ->GetVid(0));
381 vertex0->GetCoords(x0, y0, z0);
382 if (x0 != 0.0)
383 {
384 cout << "WARNING x0=" << x0 << endl;
385 x_connect = x0;
386 }
387 lastedge = -1;
388
389 v1 = 0;
390 v2 = 1;
391
392 OrderVertices(nedges, graphShPt, bndfieldx[lastIregion], Vids_c, v1, v2,
393 x_connect, lastedge, xold_c, yold_c);
394 SpatialDomains::PointGeomSharedPtr vertexc =
395 graphShPt->GetVertex(Vids_c[v2]);
396
397 // update x_connect
398 vertexc->GetCoords(x_connect, yt, zt);
399
400 i = 2;
401 while (i < nvertl)
402 {
403 v1 = i;
404 OrderVertices(nedges, graphShPt, bndfieldx[lastIregion], Vids_c, v1, v2,
405 x_connect, lastedge, xold_c, yold_c);
406 SpatialDomains::PointGeomSharedPtr vertex =
407 graphShPt->GetVertex(Vids_c[v1]);
408 // cout<<"Vids cl="<<Vids_low[v1]<<endl;
409 // update x_connect (lastedge is updated on the OrderVertices function)
410 vertex->GetCoords(x_connect, yt, zt);
411 i++;
412 }
413
414 // calculate the distances between the layer and the upper/lower curve
415
416 Array<OneD, NekDouble> Deltaup(nvertl, -200);
417 Array<OneD, NekDouble> Deltalow(nvertl, -200);
418
419 for (int r = 0; r < nvertl; r++)
420 {
421 // Always positive!!!
422 // cout<<"i="<<r<<" yup="<<yold_up[r]<<" yc="<<yold_c[r]<<"
423 // ylow="<<yold_low[r]<<endl;
424 Deltaup[r] = yold_up[r] - yold_c[r];
425 Deltalow[r] = yold_c[r] - yold_low[r];
426 ASSERTL0(Deltaup[r] > 0,
427 "distance between upper and layer curve is not positive");
428 ASSERTL0(Deltalow[r] > 0,
429 "distance between lower and layer curve is not positive");
430 }
431 //------------------------------------------------------------------------
432
433 // fieds to force continuity:
434 const SpatialDomains::BoundaryRegionCollection &bregions =
435 bcs.GetBoundaryRegions();
436 MultiRegions::ContFieldSharedPtr Cont_y;
437 MultiRegions::ExpListSharedPtr yexp;
438
439 yexp = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
440 vSession, *(bregions.find(lastIregion)->second), graphShPt, true);
441 Cont_y = MemoryManager<MultiRegions::ContField>::AllocateSharedPtr(vSession,
442 *yexp);
443 //--------------------------------------
444 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
445 // generate additional points using the Newton iteration
446 // determine the xposition for every edge (in the middle even if it
447 // is not necessary
448 // PARAMETER which determines the number of points per edge @todo put as an
449 // input
450 int npedge;
451 if (vSession->DefinesParameter("npedge"))
452 {
453 npedge = (int)vSession->GetParameter("npedge");
454 }
455 else
456 {
457 npedge = 5; // default value
458 }
459 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
460 // find the points where u=0 and determine the sign of the shift and the
461 // delta
462 int nq = streak->GetTotPoints();
463 Array<OneD, NekDouble> x(nq);
464 Array<OneD, NekDouble> y(nq);
465 streak->GetCoords(x, y);
466 Array<OneD, NekDouble> x_c(nvertl);
467 Array<OneD, NekDouble> y_c(nvertl, -200);
468 Array<OneD, NekDouble> tmp_w(nvertl, 200);
469 Array<OneD, int> Sign(nvertl, 1);
470 Array<OneD, NekDouble> Delta_c(nvertl, -200);
471
472 Computestreakpositions(nvertl, streak, xold_up, yold_up, xold_low, yold_low,
473 xold_c, yold_c, x_c, y_c, cr, true);
474 // if the curve is low the old layer point, it has to shift down
475 NekDouble shift = 0;
476 for (int q = 0; q < nvertl; q++)
477 {
478 if (y_c[q] < yold_c[q])
479 {
480 Sign[q] = -1;
481 }
482 // calculate delta
483 Delta_c[q] = abs(yold_c[q] - y_c[q]);
484 // check the shifting of the layer:
485 shift += Delta_c[q];
486 cout << x_c[q] << " " << y_c[q] << endl;
487 }
488 // cout<<"shift="<<shift<<endl;
489 if (shift < 0.001)
490 {
491 cout << "Warning: the critical layer is stationary" << endl;
492 }
493 //------------------------------------------------------------------
494
495 // additional points arrays
496 Array<OneD, NekDouble> Cpointsx(nedges);
497 Array<OneD, NekDouble> Cpointsy(nedges, 0.0);
498 Array<OneD, int> Eids(nedges);
499 Array<OneD, NekDouble> Addpointsx(nedges * (npedge - 2), 0.0);
500 Array<OneD, NekDouble> Addpointsy(nedges * (npedge - 2), 0.0);
501 // calculate the dU_dy
502 Array<OneD, NekDouble> dU(streak->GetTotPoints());
503 streak->PhysDeriv(MultiRegions::eY, streak->GetPhys(), dU);
504
505 LocalRegions::SegExpSharedPtr bndSegExp;
506
507 int Eid, id1, id2;
508 NekDouble x1, y1, z1;
509 NekDouble x2, y2, z2;
510 SpatialDomains::PointGeomSharedPtr vertex1;
511 SpatialDomains::PointGeomSharedPtr vertex2;
512 for (int r = 0; r < nedges; r++)
513 {
514
515 bndSegExp =
516 bndfieldx[lastIregion]->GetExp(r)->as<LocalRegions::SegExp>();
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);
522 vertex1->GetCoords(x1, y1, z1);
523 vertex2->GetCoords(x2, y2, z2);
524 // cout<<"edge="<<r<<" x1="<<x1<<" x2="<<x2<<endl;
525 // cout<<"edge="<<r<<" y1="<<y1<<" y2="<<y2<<endl;
526 cout << "edge=" << r << " x1=" << x1 << " y1=" << y1 << " x2=" << x2
527 << " y2=" << y2 << endl;
528 if (x2 > x1)
529 {
530 Cpointsx[r] = x1 + (x2 - x1) / 2;
531 // cout<<"edge="<<r<<" x1="<<x1<<" x2="<<x2<<"
532 // Cx="<<Cpointsx[r]<<endl; cout<<"edge="<<r<<" x1="<<x1<<"
533 // y1="<<y1<<" x2="<<x2<<" y2="<<y2<<endl;
534 if (Cpointsx[r] > x2 || Cpointsx[r] < x1)
535 {
536 Cpointsx[r] = -Cpointsx[r];
537 }
538 for (int w = 0; w < npedge - 2; w++)
539 {
540
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)
545 {
546 Addpointsx[r * (npedge - 2) + w] =
547 -Addpointsx[r * (npedge - 2) + w];
548 }
549 // initial guess along the line defined by the NEW verts y_c
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);
553 // Addpointsy[r*(npedge-2) +w] = y1 +
554 // ((y2-y1)/(x2-x1))*(Addpointsx[r*(npedge-2) +w]-x1);
555 GenerateAddPointsNewtonIt(Addpointsx[r * (npedge - 2) + w],
556 Addpointsy[r * (npedge - 2) + w],
557 Addpointsx[r * (npedge - 2) + w],
558 Addpointsy[r * (npedge - 2) + w],
559 streak, dU, cr);
560
561 // Lay_x->UpdatePhys()[r*npedge +1 +w]= Addpointsx[r*(npedge-2)
562 // +w]; Lay_y->UpdatePhys()[r*npedge +1 +w]=
563 // Addpointsy[r*(npedge-2) +w];
564 }
565 // Lay_x->UpdatePhys()[r*npedge +0] =x1;
566 // Lay_y->UpdatePhys()[r*npedge +0] =y1;
567 // Lay_x->UpdatePhys()[r*npedge +npedge-1] =x2;
568 // Lay_y->UpdatePhys()[r*npedge +npedge-1] =y2;
569 }
570 else if (x1 > x2)
571 {
572 Cpointsx[r] = x2 + (x1 - x2) / 2;
573 // cout<<"edge="<<r<<" y1="<<y1<<" y2="<<y2<<endl;
574 if (Cpointsx[r] > x1 || Cpointsx[r] < x2)
575 {
576 Cpointsx[r] = -Cpointsx[r];
577 }
578 for (int w = 0; w < npedge - 2; w++)
579 {
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)
584 {
585 Addpointsx[r * (npedge - 2) + w] =
586 -Addpointsx[r * (npedge - 2) + w];
587 }
588
589 // initial guess along the line defined by the NEW verts y_c
590 Addpointsy[r * (npedge - 2) + w] =
591 y_c[r + 1] +
592 ((y_c[r] - y_c[r + 1]) / (x_c[r] - x_c[r + 1])) *
593 (Addpointsx[r * (npedge - 2) + w] - x2);
594
595 // Addpointsy[r*(npedge-2) +w] = y2 +
596 // ((y1-y2)/(x1-x2))*(Addpointsx[r*(npedge-2) +w]-x2);
597 GenerateAddPointsNewtonIt(Addpointsx[r * (npedge - 2) + w],
598 Addpointsy[r * (npedge - 2) + w],
599 Addpointsx[r * (npedge - 2) + w],
600 Addpointsy[r * (npedge - 2) + w],
601 streak, dU, cr);
602 // Lay_x->UpdatePhys()[r*npedge +1]= Addpointsx[r*(npedge-2)
603 // +w]; Lay_y->UpdatePhys()[r*npedge +1]=
604 // Addpointsy[r*(npedge-2) +w];
605 }
606 // Lay_x->UpdatePhys()[r*npedge +0] =x2;
607 // Lay_y->UpdatePhys()[r*npedge +0] =y2;
608 // Lay_x->UpdatePhys()[r*npedge +npedge-1] =x1;
609 // Lay_y->UpdatePhys()[r*npedge +npedge-1] =y1;
610 }
611 else
612 {
613 ASSERTL0(false, "point not generated");
614 }
615 // cout<<"calculate cpoints coords"<<endl;
616 // Cpointsy[r] = y1 + (y2-y1)/2;
617 // cout<<"central point:"<<endl;
618 // GenerateAddPointsNewtonIt( Cpointsx[r], Cpointsy[r],Cpointsx[r],
619 // Cpointsy[r],
620 // streak, dU,cr);
621 // NekDouble diff = Cpointsy[r]-Addpointsy[r*(npedge-2)];
622 // cout<<"diff="<<diff<<endl;
623 Eids[r] = Eid;
624 }
625 //-------------------------------------------------------------
626
627 // fill the xPhys,yPhys array( may necessary after)
628 Array<OneD, NekDouble> xcPhys(nedges * npedge, 0.0);
629 Array<OneD, NekDouble> ycPhys(nedges * npedge, 0.0);
630
631 for (int a = 0; a < nedges; a++)
632 {
633 // v1
634 xcPhys[a * npedge + 0] = x_c[a];
635 ycPhys[a * npedge + 0] = y_c[a];
636 // v2
637 xcPhys[a * npedge + npedge - 1] = x_c[a + 1];
638 ycPhys[a * npedge + npedge - 1] = y_c[a + 1];
639
640 for (int b = 0; b < npedge - 2; b++)
641 {
642 xcPhys[a * npedge + b + 1] = Addpointsx[a * (npedge - 2) + b];
643 ycPhys[a * npedge + b + 1] = Addpointsy[a * (npedge - 2) + b];
644 }
645 }
646
647 cout << "xc,yc before tanevaluate" << endl;
648 for (int v = 0; v < xcPhys.size(); v++)
649 {
650 cout << xcPhys[v] << " " << ycPhys[v] << endl;
651 }
652
653 //-------------------------------------------------
654
655 // V1[eid],V2[eid] vertices associate with the edge Id=eid
656 Array<OneD, int> V1;
657 Array<OneD, int> V2;
658
659 GenerateMapEidsv1v2(streak, V1, V2);
660 Array<OneD, Array<OneD, int>> lay_Eids;
661 Array<OneD, Array<OneD, int>> lay_Vids;
662 int nlays = 0;
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);
665
666 cout << "nlays=" << nlays << endl;
667 Array<OneD, Array<OneD, NekDouble>> layers_y(nlays);
668 Array<OneD, Array<OneD, NekDouble>> layers_x(nlays);
669 // initialise layers_y,lay_eids
670 for (int g = 0; g < nlays; g++)
671 {
672 layers_y[g] = Array<OneD, NekDouble>((nvertl - 1) * npedge);
673 }
674
675 /////////////////////////////////////////////////////////////////////////////
676 ////////////////////////////////////////////////////////////////////////////
677 //££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££
678 //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
679 //@todo set Delta0 from session file
680 NekDouble Delta0;
681 if (vSession->DefinesParameter("Delta"))
682 {
683 Delta0 = vSession->GetParameter("Delta");
684 }
685 else
686 {
687 Delta0 = 0.1; // default value
688 }
689
690 //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
691 //££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££
692 ////////////////////////////////////////////////////////////////////////////
693 /////////////////////////////////////////////////////////////////////////////
694 // save the coords of the old vertices
695 int nVertTot = graphShPt->GetNvertices();
696 // arrays to check with the new values
697 Array<OneD, NekDouble> xold(nVertTot);
698 Array<OneD, NekDouble> yold(nVertTot);
699 // calculate the new cordinates of the vertices
700 Array<OneD, NekDouble> xnew(nVertTot);
701 Array<OneD, NekDouble> ynew(nVertTot, -20);
702 Array<OneD, int> Up(nvertl); // Vids lay Up
703 Array<OneD, int> Down(nvertl); // Vids lay Down
704 int cntup = 0;
705 int cntlow = 0;
706 // Vids needed only if a layers has to be moved
707 NekDouble bleft = -10;
708 NekDouble tright = 10;
709 NekDouble bright = -10;
710 NekDouble tleft = 10;
711 for (int i = 0; i < nVertTot; i++)
712 {
713 bool mvpoint = false;
714 SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(i);
715 NekDouble x, y, z;
716 vertex->GetCoords(x, y, z);
717
718 xold[i] = x;
719 yold[i] = y;
720
721 // x coord doesn't change
722 xnew[i] = x;
723 // cout<<"x="<<x<<" y="<<y<<endl;
724 // bottom, left (x=0, y<ydown)
725 if (x == 0 && y < yold_low[0] && y > bleft)
726 {
727 bleft = y;
728 }
729 // top, right
730 if (x == xold_c[nvertl - 1] && y > yold_up[nvertl - 1] && y < tright)
731 {
732 tright = y;
733 }
734 // bottom, right]
735 if (x == xold_c[nvertl - 1] && y < yold_low[nvertl - 1] && y > bright)
736 {
737 bright = y;
738 }
739 // top, left
740 if (x == 0 && y > yold_up[0] && y < tleft)
741 {
742 tleft = y;
743 }
744 // find the corresponding yold_l and deltaold
745 for (int j = 0; j < nvertl; j++)
746 {
747 if ((xold_up[j] == x) && (yold_up[j] == y))
748 {
749 // ratio = (1-y)*(1+y)/( (1-yold_c[j])*(1+yold_c[j]) );
750 // ynew[i] = y + Sign[j]*Delta_c[j]*ratio;
751 ynew[i] = y_c[j] + Delta0;
752 mvpoint = true;
753 Up[j] = i;
754 }
755 if ((xold_low[j] == x) && (yold_low[j] == y))
756 {
757 // ratio = (1-y)*(1+y)/( (1-yold_c[j])*(1+yold_c[j]) );
758 // ynew[i] = y + Sign[j]*Delta_c[j]*ratio;
759 ynew[i] = y_c[j] - Delta0;
760 mvpoint = true;
761 Down[j] = i;
762 }
763 if ((xold_c[j] == x) && (yold_c[j] == y))
764 {
765 ynew[i] = y_c[j];
766 mvpoint = true;
767 }
768 }
769 if (mvpoint == false)
770 {
771 // determine the closer xold_up
772 NekDouble diff = 1000;
773 int qp_closer = 0;
774 for (int k = 0; k < nvertl; k++)
775 {
776 if (abs(x - xold_up[k]) < diff)
777 {
778 diff = abs(x - xold_up[k]);
779 qp_closer = k;
780 }
781 }
782 if (y > yold_up[qp_closer] && y < 1)
783 {
784
785 // ratio = (1-y)*(1+y)/(
786 // (1-yold_c[qp_closer])*(1+yold_c[qp_closer]) ); ynew[i] = y +
787 // Sign[qp_closer]*Delta_c[qp_closer]*ratio;
788
789 // ratio = (1-y)*(1+y)/(
790 // (1-yold_up[qp_closer])*(1+yold_up[qp_closer]) ); distance
791 // prop to layerup
792 ynew[i] = y_c[qp_closer] + (y - yold_c[qp_closer]) *
793 (1 - y_c[qp_closer]) /
794 (1 - yold_c[qp_closer]);
795 // cout<<"upper zone y="<<y<<" ratio="<<ratio<<endl;
796 }
797 else if (y < yold_low[qp_closer] && y > -1)
798 {
799
800 // ratio = (1-y)*(1+y)/(
801 // (1-yold_c[qp_closer])*(1+yold_c[qp_closer]) ); ynew[i] = y +
802 // Sign[qp_closer]*Delta_c[qp_closer]*ratio;
803
804 // ratio = (1-y)*(1+y)/(
805 // (1-yold_low[qp_closer])*(1+yold_low[qp_closer]) ); distance
806 // prop to layerlow
807 ynew[i] = y_c[qp_closer] + (y - yold_c[qp_closer]) *
808 (-1 - y_c[qp_closer]) /
809 (-1 - yold_c[qp_closer]);
810 }
811
812 else if (y > yold_c[qp_closer] && y < yold_up[qp_closer])
813 {
814 if (x == 0)
815 {
816 cntup++;
817 }
818 // ratio = (1-y)*(1+y)/(
819 // (1-yold_c[qp_closer])*(1+yold_c[qp_closer]) ); ynew[i] = y +
820 // Sign[qp_closer]*Delta_c[qp_closer]*ratio;
821 }
822 else if (y < yold_c[qp_closer] && y > yold_low[qp_closer])
823 {
824 if (x == 0)
825 {
826 cntlow++;
827 }
828 // ratio = (1-y)*(1+y)/(
829 // (1-yold_c[qp_closer])*(1+yold_c[qp_closer]) ); ynew[i] = y +
830 // Sign[qp_closer]*Delta_c[qp_closer]*ratio;
831 }
832 else if (y == 1 || y == -1) // bcs don't move
833 {
834 ynew[i] = y;
835 // cout<<"point x="<<xnew[i]<<" y="<<y<<" closer
836 // x="<<xold_up[qp_closer]<<endl;
837 }
838
839 // internal layers are not moved yet so...
840 if ((ynew[i] > 1 || ynew[i] < -1) &&
841 (y > yold_up[qp_closer] || y < yold_low[qp_closer]))
842 {
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");
847 }
848 }
849 }
850
851 int nqedge = streak->GetExp(0)->GetNumPoints(0);
852 int nquad_lay = (nvertl - 1) * nqedge;
853 Array<OneD, NekDouble> coeffstmp(Cont_y->GetNcoeffs(), 0.0);
854 // curve the edges around the NEW critical layer (bool to turn on/off)
855 bool curv_lay = true;
856 bool move_norm = true;
857 int np_lay = (nvertl - 1) * npedge; // nedges*npedge (Eq. Points!!!)
858
859 Array<OneD, NekDouble> xcQ(nqedge * nedges, 0.0);
860 Array<OneD, NekDouble> ycQ(nqedge * nedges, 0.0);
861 Array<OneD, NekDouble> zcQ(nqedge * nedges, 0.0);
862 Array<OneD, NekDouble> nxPhys(npedge * nedges, 0.0);
863 Array<OneD, NekDouble> nyPhys(npedge * nedges, 0.0);
864 Array<OneD, NekDouble> nxQ(nqedge * nedges, 0.0);
865 Array<OneD, NekDouble> nyQ(nqedge * nedges, 0.0);
866
867 if (move_norm == true)
868 {
869 // np_lay = (nvertl-1)*nqedge;//nedges*nqedge (Q points!!!)
870 // extract crit lay normals (through tangents):
871 Array<OneD, NekDouble> xcPhysMOD(xcPhys.size());
872 Array<OneD, NekDouble> ycPhysMOD(ycPhys.size());
873 // copy(temporary the xcPhys,ycPhys into xcPhysMOD, ycPhysMOD)
874 Vmath::Vcopy(xcPhys.size(), xcPhys, 1, xcPhysMOD, 1);
875 Vmath::Vcopy(xcPhys.size(), ycPhys, 1, ycPhysMOD, 1);
876
877 Array<OneD, StdRegions::StdExpansion1DSharedPtr> Edge_newcoords(2);
878
879 cout << "nquad per edge=" << nqedge << endl;
880 for (int l = 0; l < 2; l++)
881 {
882 Edge_newcoords[l] = bndfieldx[lastIregion]
883 ->GetExp(0)
884 ->as<StdRegions::StdExpansion1D>();
885 }
886 Array<OneD, NekDouble> xnull(nqedge);
887 Array<OneD, NekDouble> ynull(nqedge);
888 Array<OneD, NekDouble> xcedgeQ(nqedge, 0.0);
889 Array<OneD, NekDouble> ycedgeQ(nqedge, 0.0);
890 Array<OneD, NekDouble> txedgeQ(nqedge, 0.0);
891 Array<OneD, NekDouble> tyedgeQ(nqedge, 0.0);
892 Array<OneD, NekDouble> normsQ(nqedge, 0.0);
893
894 Array<OneD, NekDouble> txQ(nqedge * nedges, 0.0);
895 Array<OneD, NekDouble> tyQ(nqedge * nedges, 0.0);
896 Array<OneD, NekDouble> tx_tyedgeQ(nqedge, 0.0);
897 Array<OneD, NekDouble> nxedgeQ(nqedge, 0.0);
898 Array<OneD, NekDouble> nyedgeQ(nqedge, 0.0);
899 Array<OneD, const NekDouble> ntempQ(nqedge);
900
901 Array<OneD, NekDouble> nxedgePhys(npedge, 0.0);
902 Array<OneD, NekDouble> nyedgePhys(npedge, 0.0);
903
904 Array<OneD, NekDouble> CoordsPhys(2);
905
906 bndfieldx[lastIregion]->GetCoords(xcQ, ycQ, zcQ);
907 // determine the NEW crit lay quad points values(lagrangeInterpolant):
908 // interp(from xcPhys, ycPhys).
909 Array<OneD, NekDouble> x_tmp(np_lay - (nedges - 1), 0.0);
910 Array<OneD, NekDouble> y_tmp(np_lay - (nedges - 1), 0.0);
911 Array<OneD, NekDouble> closex(4, 0.0);
912 Array<OneD, NekDouble> closey(4, 0.0);
913 Cutrepetitions(nedges, xcPhysMOD, x_tmp);
914 Cutrepetitions(nedges, ycPhysMOD, y_tmp);
915 for (int l = 0; l < xcQ.size(); l++)
916 {
917 Cutrepetitions(nedges, xcPhysMOD, x_tmp);
918 Cutrepetitions(nedges, ycPhysMOD, y_tmp);
919 int indclose = DetermineclosePointxindex(xcQ[l], x_tmp);
920
921 // generate neighbour arrays
922 GenerateNeighbourArrays(indclose, 4, x_tmp, y_tmp, closex, closey);
923
924 ycQ[l] = LagrangeInterpolant(xcQ[l], 4, closex, closey);
925 }
926
927 // force continuity
928
929 Vmath::Vcopy(nquad_lay, ycQ, 1, Cont_y->UpdatePhys(), 1);
930 Array<OneD, NekDouble> coeffsy(Cont_y->GetNcoeffs(), 0.0);
931
932 Cont_y->FwdTransLocalElmt(Cont_y->GetPhys(), coeffsy);
933 Cont_y->BwdTrans(coeffsy, Cont_y->UpdatePhys());
934 Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(), 1, ycQ, 1);
935
936 cout << "xcQ, ycQ" << endl;
937 for (int s = 0; s < xcQ.size(); s++)
938 {
939 cout << xcQ[s] << " " << ycQ[s] << endl;
940 }
941 // ASSERTL0(false, "dsdfs");
942 bool evaluatetan = false;
943 NekDouble incratio;
944 Array<OneD, int> edgeinterp(2);
945 int wrong = 0;
946 for (int k = 0; k < nedges; k++)
947 {
948 Vmath::Vcopy(nqedge, &xcQ[k * nqedge], 1, &xcedgeQ[0], 1);
949 Vmath::Vcopy(nqedge, &ycQ[k * nqedge], 1, &ycedgeQ[0], 1);
950 // calc the NEW tangent values
951 Edge_newcoords[0]->StdPhysDeriv(xcedgeQ, txedgeQ);
952 Edge_newcoords[1]->StdPhysDeriv(ycedgeQ, tyedgeQ);
953 // norms=tx*tx
954 Vmath::Vmul(nqedge, txedgeQ, 1, txedgeQ, 1, normsQ, 1);
955 // norms=tx*tx+ty*ty
956 Vmath::Vvtvp(nqedge, tyedgeQ, 1, tyedgeQ, 1, normsQ, 1, normsQ, 1);
957
958 Vmath::Vsqrt(nqedge, normsQ, 1, normsQ, 1);
959 Vmath::Sdiv(nqedge, 1.0, normsQ, 1, normsQ, 1);
960
961 Vmath::Vmul(nqedge, txedgeQ, 1, normsQ, 1, txedgeQ, 1);
962 Vmath::Vmul(nqedge, tyedgeQ, 1, normsQ, 1, tyedgeQ, 1);
963
964 // try evaluate tangent if the incremental ratio is high
965
966 evaluatetan = false;
967 for (int u = 0; u < nqedge - 1; u++)
968 {
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)
973 {
974 cout << "wrong=" << wrong << endl;
975 evaluatetan = true;
976 ASSERTL0(wrong < 2, "number edges to change is too high!!");
977 edgeinterp[wrong] = k;
978 wrong++;
979 }
980 }
981
982 if (evaluatetan)
983 {
984 cout << "tan bef" << endl;
985 for (int e = 0; e < nqedge; e++)
986 {
987 cout << xcedgeQ[e] << " " << ycedgeQ[e] << " "
988 << txedgeQ[e] << endl;
989 }
990
991 // OR: fit
992 int polyorder = 3;
993 Array<OneD, NekDouble> coeffsinterp(polyorder + 1);
994 Array<OneD, NekDouble> yPedges(npedge, 0.0);
995 Array<OneD, NekDouble> xPedges(npedge, 0.0);
996 Vmath::Vcopy(npedge, &xcPhysMOD[k * npedge + 0], 1, &xPedges[0],
997 1);
998 Vmath::Vcopy(npedge, &ycPhysMOD[k * npedge + 0], 1, &yPedges[0],
999 1);
1000
1001 PolyFit(polyorder, nqedge, xcedgeQ, ycedgeQ, coeffsinterp,
1002 xPedges, yPedges, npedge);
1003 // update values
1004 Vmath::Vcopy(npedge, &xPedges[0], 1, &xcPhysMOD[k * npedge + 0],
1005 1);
1006 Vmath::Vcopy(npedge, &yPedges[0], 1, &ycPhysMOD[k * npedge + 0],
1007 1);
1008
1009 EvaluateTangent(nqedge, xcedgeQ, coeffsinterp, txedgeQ,
1010 tyedgeQ);
1011 }
1012 // copy also the tx,ty
1013 Vmath::Vcopy(nqedge, &(txedgeQ[0]), 1, &(txQ[nqedge * k]), 1);
1014 Vmath::Vcopy(nqedge, &(tyedgeQ[0]), 1, &(tyQ[nqedge * k]), 1);
1015 }
1016
1017 Array<OneD, NekDouble> fz(nedges * nqedge, 0.0);
1018 bndfieldx[lastIregion]->PhysDeriv(MultiRegions::eX, ycQ, fz);
1019 for (int w = 0; w < fz.size(); w++)
1020 {
1021 txQ[w] = cos(atan(fz[w]));
1022 tyQ[w] = sin(atan(fz[w]));
1023 cout << xcQ[w] << " " << ycQ[w] << " " << fz[w] << endl;
1024 }
1025 // ASSERTL0(false, "bobo");
1026 // force continuity tx,ty
1027 // tx
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());
1031 Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(), 1, txQ, 1);
1032 // ty
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());
1036 Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(), 1, tyQ, 1);
1037
1038 // check if the tan points have the same curvature otherwise interp
1039 NekDouble inc, incbefore;
1040
1041 // build-up the fit for the tan using the edge with
1042 // the same derivative sign (before or after)
1043
1044 int edgebef;
1045
1046 for (int q = 0; q < 2; q++)
1047 {
1048 edgebef = edgeinterp[q] - 1;
1049 incbefore =
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;
1057
1058 Array<OneD, NekDouble> yQedges(npoints, 0.0);
1059 Array<OneD, NekDouble> xQedges(npoints, 0.0);
1060 Array<OneD, NekDouble> tyQedges(npoints, 0.0);
1061 Array<OneD, NekDouble> txQedges(npoints, 0.0);
1062 cout << "inc=" << inc << " incbef=" << incbefore << endl;
1063 if ((inc / incbefore) > 0.)
1064 {
1065 cout << "before!!" << edgebef << endl;
1066 int polyorder = 2;
1067 Array<OneD, NekDouble> coeffsinterp(polyorder + 1);
1068 Vmath::Vcopy(npoints, &xcQ[edgebef * nqedge + 0], 1,
1069 &xQedges[0], 1);
1070 Vmath::Vcopy(npoints, &ycQ[edgebef * nqedge + 0], 1,
1071 &yQedges[0], 1);
1072 Vmath::Vcopy(npoints, &txQ[edgebef * nqedge + 0], 1,
1073 &txQedges[0], 1);
1074 Vmath::Vcopy(npoints, &tyQ[edgebef * nqedge + 0], 1,
1075 &tyQedges[0], 1);
1076
1077 PolyFit(polyorder, npoints, xQedges, txQedges, coeffsinterp,
1078 xQedges, txQedges, npoints);
1079
1080 // copy back the values:
1081 Vmath::Vcopy(npoints, &txQedges[0], 1,
1082 &txQ[edgebef * nqedge + 0], 1);
1083
1084 PolyFit(polyorder, npoints, xQedges, tyQedges, coeffsinterp,
1085 xQedges, tyQedges, npoints);
1086
1087 // copy back the values:
1088 Vmath::Vcopy(npoints, &tyQedges[0], 1,
1089 &tyQ[edgebef * nqedge + 0], 1);
1090 }
1091 else
1092 {
1093 cout << "after!!" << endl;
1094 int polyorder = 2;
1095 Array<OneD, NekDouble> coeffsinterp(polyorder + 1);
1096 Vmath::Vcopy(npoints, &xcQ[edgeinterp[q] * nqedge + 0], 1,
1097 &xQedges[0], 1);
1098 Vmath::Vcopy(npoints, &ycQ[edgeinterp[q] * nqedge + 0], 1,
1099 &yQedges[0], 1);
1100 Vmath::Vcopy(npoints, &txQ[edgeinterp[q] * nqedge + 0], 1,
1101 &txQedges[0], 1);
1102 Vmath::Vcopy(npoints, &tyQ[edgeinterp[q] * nqedge + 0], 1,
1103 &tyQedges[0], 1);
1104
1105 PolyFit(polyorder, npoints, xQedges, txQedges, coeffsinterp,
1106 xQedges, txQedges, npoints);
1107
1108 // copy back the values:
1109 Vmath::Vcopy(npoints, &txQedges[0], 1,
1110 &txQ[edgeinterp[q] * nqedge + 0], 1);
1111
1112 PolyFit(polyorder, npoints, xQedges, tyQedges, coeffsinterp,
1113 xQedges, tyQedges, npoints);
1114
1115 // copy back the values:
1116 Vmath::Vcopy(npoints, &tyQedges[0], 1,
1117 &tyQ[edgeinterp[q] * nqedge + 0], 1);
1118 }
1119 }
1120 // force continuity of the tangent
1121 // tyQ
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());
1125 Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(), 1, tyQ, 1);
1126 // txQ
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());
1130 Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(), 1, txQ, 1);
1131
1132 for (int k = 0; k < nedges; k++)
1133 {
1134 // determine the normal from eqs(APART FROM SIGN):
1135 // tx^2 +ty^2= 1 = nx^2 + ny^2;
1136 // t\cdot n=0= tx*nx +ty*ny
1137 // result: nx = ( 1+(tx/ty)^2 )^(-1/2)
1138
1139 Vmath::Vcopy(nqedge, &(txQ[nqedge * k]), 1, &(txedgeQ[0]), 1);
1140 Vmath::Vcopy(nqedge, &(tyQ[nqedge * k]), 1, &(tyedgeQ[0]), 1);
1141
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);
1145 Vmath::Vsqrt(nqedge, nxedgeQ, 1, nxedgeQ, 1);
1146 Vmath::Sdiv(nqedge, 1.0, nxedgeQ, 1, nxedgeQ, 1);
1147 // normal DOWNWARDS!!! mult by -1
1148 Vmath::Smul(nqedge, -1.0, nxedgeQ, 1, nxedgeQ, 1);
1149 Vmath::Vcopy(nqedge, &(nxedgeQ[0]), 1, &(nxQ[nqedge * k]), 1);
1150 // ny = (1-nx ^2)^(1/2)
1151 Vmath::Vmul(nqedge, nxedgeQ, 1, nxedgeQ, 1, nyedgeQ, 1);
1152 Vmath::Smul(nqedge, -1.0, nyedgeQ, 1, nyedgeQ, 1);
1153 Vmath::Sadd(nqedge, 1.0, nyedgeQ, 1, nyedgeQ, 1);
1154 Vmath::Vsqrt(nqedge, nyedgeQ, 1, nyedgeQ, 1);
1155 // normal DOWNWARDS!!! mult by -1
1156 Vmath::Smul(nqedge, -1.0, nyedgeQ, 1, nyedgeQ, 1);
1157 Vmath::Vcopy(nqedge, &(nyedgeQ[0]), 1, &(nyQ[nqedge * k]), 1);
1158
1159 cout << "edge:" << k << endl;
1160 cout << "tan/normal" << endl;
1161 for (int r = 0; r < nqedge; r++)
1162 {
1163 cout << xcQ[k * nqedge + r] << " " << txedgeQ[r] << " "
1164 << tyedgeQ[r] << " " << nxedgeQ[r] << " "
1165 << nyedgeQ[r] << endl;
1166 }
1167 }
1168
1169 // force continuity:
1170 // REMEMBER: the Fwd/Bwd operation get wrong with the border values!!!
1171 Vmath::Vcopy(nquad_lay, nyQ, 1, Cont_y->UpdatePhys(), 1);
1172
1173 Cont_y->FwdTransLocalElmt(Cont_y->GetPhys(), coeffstmp);
1174 Cont_y->BwdTrans(coeffstmp, Cont_y->UpdatePhys());
1175 Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(), 1, nyQ, 1);
1176
1177 Vmath::Zero(nquad_lay, Cont_y->UpdatePhys(), 1);
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());
1182 Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(), 1, nxQ, 1);
1183
1184 // force the normal at interface point to be equal
1185 for (int k = 0; k < nedges; k++)
1186 {
1187 if (k > 0)
1188 {
1189 // nyPhys[f*npedge +0] =
1190 // (nyPhys[(f-1)*npedge+npedge-1]+nyPhys[f*npedge+0])/2.;
1191 nyQ[(k - 1) * nqedge + nqedge - 1] = nyQ[k * nqedge + 0];
1192 // nx= (1-ny^2)^{1/2}
1193 // nxPhys[f*npedge+0]=
1194 // sqrt(1- nyPhys[f*npedge+0]*nyPhys[f*npedge+0]);
1195 nxQ[(k - 1) * nqedge + nqedge - 1] = nxQ[k * nqedge + 0];
1196 }
1197 }
1198
1199 Array<OneD, NekDouble> x_tmpQ(nquad_lay - (nedges - 1));
1200 // Array<OneD, NekDouble>tmpnxQ(nquad_lay-(nedges-1));
1201 Array<OneD, NekDouble> tmpnyQ(nquad_lay - (nedges - 1));
1202
1203 cout << "nx,yQbefore" << endl;
1204 for (int u = 0; u < xcQ.size(); u++)
1205 {
1206 cout << xcQ[u] << " " << nyQ[u] << " " << txQ[u] << endl;
1207 }
1208
1209 Cutrepetitions(nedges, xcQ, x_tmpQ);
1210 // Cutrepetitions(nedges, nxQ, tmpnxQ);
1211 Cutrepetitions(nedges, nyQ, tmpnyQ);
1212 cout << "nx,yQ" << endl;
1213 for (int u = 0; u < x_tmpQ.size(); u++)
1214 {
1215 cout << x_tmpQ[u] << " " << tmpnyQ[u] << endl;
1216 }
1217
1218 // interpolate the values into phys points(curved points)
1219 for (int k = 0; k < nedges; k++)
1220 {
1221 // cout<<"edge:"<<k<<endl;
1222 for (int a = 0; a < npedge; a++)
1223 {
1224 if (a == 0) // verts pos no interp necessary
1225 {
1226 nxPhys[k * npedge + a] = nxQ[k * nqedge + 0];
1227 nyPhys[k * npedge + a] = nyQ[k * nqedge + 0];
1228 }
1229 else if (a == npedge - 1) // verts pos no interp necessary
1230 {
1231 nxPhys[k * npedge + a] = nxQ[k * nqedge + nqedge - 1];
1232 nyPhys[k * npedge + a] = nyQ[k * nqedge + nqedge - 1];
1233 // cout<<":last"<<nyQ[k*nqedge+a]<<endl;
1234 }
1235 else
1236 {
1237 // use lagrange interpolant to get the
1238 // normal at phys(equispaced points)
1239
1240 // order normal functions(cut out repetitions)
1241 // QUAD POINTS
1242
1243 int index;
1244 // determine closest index:
1245
1246 index = DetermineclosePointxindex(
1247 Addpointsx[k * (npedge - 2) + a - 1], x_tmpQ);
1248
1249 Array<OneD, NekDouble> Pxinterp(4);
1250 Array<OneD, NekDouble> Pyinterp(4);
1251
1252 // generate neighbour arrays (y):
1253 GenerateNeighbourArrays(index, 4, x_tmpQ, tmpnyQ, Pxinterp,
1254 Pyinterp);
1255 // interp the new normal components(y)
1256 nyPhys[k * npedge + a] = LagrangeInterpolant(
1257 Addpointsx[k * (npedge - 2) + a - 1], 4, Pxinterp,
1258 Pyinterp);
1259 /*
1260 //generate neighbour arrays (x):
1261 GenerateNeighbourArrays(index,4,x_tmpQ,tmpnxQ,Pxinterp,Pyinterp);
1262 //interp the new normal components(x)
1263 nxPhys[k*npedge +a]=
1264 LagrangeInterpolant(Addpointsx[k*(npedge-2)
1265 +a],4,Pxinterp,Pyinterp );
1266 */
1267 // nx=-(1-ny*ny){1/2} the normal is DOWNWARDS!!!
1268 nxPhys[k * npedge + a] = -sqrt(abs(
1269 1 - nyPhys[k * npedge + a] * nyPhys[k * npedge + a]));
1270 /*
1271 //(put the middle points as quad points)
1272 //GaussLobattoLegendre points in the middle:
1273 nxPhys[k*npedge +a] = nxedgeQ[a];
1274 nyPhys[k*npedge +a] = nyedgeQ[a];
1275 ASSERTL0(npedge< nqedge," quad points too low");
1276 */
1277 }
1278
1279 // force the normal at interface point to be equal
1280 if (k > 0)
1281 {
1282 // nyPhys[f*npedge +0] =
1283 // (nyPhys[(f-1)*npedge+npedge-1]+nyPhys[f*npedge+0])/2.;
1284 nyPhys[(k - 1) * npedge + npedge - 1] =
1285 nyPhys[k * npedge + 0];
1286 // nx= (1-ny^2)^{1/2}
1287 // nxPhys[f*npedge+0]=
1288 // sqrt(1- nyPhys[f*npedge+0]*nyPhys[f*npedge+0]);
1289 nxPhys[(k - 1) * npedge + npedge - 1] =
1290 nxPhys[k * npedge + 0];
1291 }
1292 }
1293 }
1294 cout << "xcPhys,," << endl;
1295 for (int s = 0; s < np_lay; s++)
1296 {
1297
1298 cout << xcPhysMOD[s] << " " << ycPhysMOD[s] << " "
1299 << nxPhys[s] << " " << nyPhys[s] << endl;
1300 }
1301
1302 // determine the new coords of the vertices and the curve points
1303 // for each edge
1304
1305 // int np_lay = (nvertl-1)*npedge;//nedges*npedge
1306
1307 // NB delta=ynew-y_c DEPENDS ON the coord trnaf ynew= y+Delta_c*ratio!!!
1308 Array<OneD, NekDouble> delta(nlays);
1309 Array<OneD, NekDouble> tmpy_lay(np_lay);
1310 Array<OneD, NekDouble> tmpx_lay(np_lay);
1311 for (int m = 0; m < nlays; m++)
1312 {
1313 // delta[m] = (ynew[lay_Vids[m][0]] - y_c[0])/1.0;
1314
1315 // depends on Delta0
1316 if (m < cntlow + 1)
1317 {
1318 delta[m] = -(cntlow + 1 - m) * Delta0 / (cntlow + 1);
1319 }
1320 else
1321 {
1322 delta[m] = (m - (cntlow)) * Delta0 / (cntlow + 1);
1323 }
1324
1325 layers_x[m] = Array<OneD, NekDouble>(np_lay);
1326 // cout<<"delta="<<delta[m]<<" cntlow="<<cntlow<<endl;
1327
1328 for (int h = 0; h < nvertl; h++)
1329 {
1330 // shift using the dinstance delta from the crit layer AT x=0
1331 // for each layer
1332
1333 // cout<<m<<"Vid:"<<lay_Vids[m][h]<<" mod from
1334 // y="<<ynew[lay_Vids[m][h] ]<<" to y="<<y_c[h]
1335 // +delta[m]<<endl;
1336 if (move_norm == false)
1337 {
1338 ynew[lay_Vids[m][h]] = y_c[h] + delta[m];
1339 xnew[lay_Vids[m][h]] = x_c[h];
1340 }
1341 else
1342 {
1343 if (h == 0 || h == nvertl - 1) // nx=0,ny=1 at the borders
1344 {
1345 ynew[lay_Vids[m][h]] = y_c[h] + delta[m];
1346 xnew[lay_Vids[m][h]] = x_c[h];
1347 }
1348 else
1349 {
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]);
1354 }
1355 }
1356 cout << "Vid x=" << xnew[lay_Vids[m][h]]
1357 << " y=" << ynew[lay_Vids[m][h]] << endl;
1358 if (h < nedges
1359 //&& curv_lay==true
1360 )
1361 {
1362 cout << "edge==" << h << endl;
1363 if (h > 0) // check normal consistency
1364 {
1365 ASSERTL0(nyPhys[h * npedge + 0] ==
1366 nyPhys[(h - 1) * npedge + npedge - 1],
1367 " normaly wrong");
1368 ASSERTL0(nxPhys[h * npedge + 0] ==
1369 nxPhys[(h - 1) * npedge + npedge - 1],
1370 " normalx wrong");
1371 }
1372 if (move_norm == false)
1373 {
1374 // v1
1375 layers_y[m][h * npedge + 0] = y_c[h] + delta[m];
1376 layers_x[m][h * npedge + 0] = xnew[lay_Vids[m][h]];
1377 // v2
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]];
1382 // middle points (shift crit lay points by delta):
1383 for (int d = 0; d < npedge - 2; d++)
1384 {
1385 layers_y[m][h * npedge + d + 1] =
1386 ycPhysMOD[h * npedge + d + 1] + delta[m];
1387 // Addpointsy[h*(npedge-2) +d] +delta[m];
1388 layers_x[m][h * npedge + d + 1] =
1389 xcPhysMOD[h * npedge + d + 1];
1390 // Addpointsx[h*(npedge-2) +d];
1391 }
1392 }
1393 else
1394 {
1395 if (h == 0) // nx=0,ny=1 at the borders
1396 {
1397 // v1
1398 tmpy_lay[h * npedge + 0] = y_c[h] + delta[m];
1399 tmpx_lay[h * npedge + 0] = xnew[lay_Vids[m][h]];
1400 // v2
1401 tmpy_lay[h * npedge + npedge - 1] =
1402 y_c[h + 1] +
1403 delta[m] * abs(nyPhys[h * npedge + npedge - 1]);
1404 tmpx_lay[h * npedge + npedge - 1] =
1405 x_c[h + 1] +
1406 delta[m] * abs(nxPhys[h * npedge + npedge - 1]);
1407 }
1408 else if (h == nedges - 1) // nx=0,ny=1 at the borders
1409 {
1410 // v1
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]);
1415 // v2
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]];
1420 }
1421 else
1422 {
1423 // v1
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]);
1428 // v2
1429 tmpy_lay[h * npedge + npedge - 1] =
1430 y_c[h + 1] +
1431 delta[m] * abs(nyPhys[h * npedge + npedge - 1]);
1432 tmpx_lay[h * npedge + npedge - 1] =
1433 x_c[h + 1] +
1434 delta[m] * abs(nxPhys[h * npedge + npedge - 1]);
1435 }
1436
1437 // middle points
1438 for (int d = 0; d < npedge - 2; d++)
1439 {
1440
1441 tmpy_lay[h * npedge + d + 1] =
1442 ycPhysMOD[h * npedge + d + 1] +
1443 delta[m] * abs(nyPhys[h * npedge + d + 1]);
1444 // Addpointsy[h*(npedge-2) +d] +
1445 // 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]);
1449 // Addpointsx[h*(npedge-2) +d] +
1450 // delta[m]*abs(nxPhys[h*npedge +d+1]);
1451
1452 // NB ycQ,xcQ refers to nqedge NOT npedge!!!
1453 // tmpy_lay[h*npedge +d+1] = ycQ[h*nqedge +d+1] +
1454 // delta[m]*abs(nyPhys[h*npedge +d+1]);
1455 // tmpx_lay[h*npedge +d+1]= xcQ[h*nqedge +d+1] +
1456 // delta[m]*abs(nxPhys[h*npedge +d+1]);
1457 // cout<<"xmoved="<<tmpx_lay[h*npedge +d+1]<<"
1458 // xold="<<xcQ[h*nqedge +d+1]<<endl;
1459 // ASSERTL0(tmpx_lay[h*npedge +d+1]>0," middle point
1460 // with x<0")
1461 }
1462 }
1463 } // close edges
1464 } // close verts h
1465
1466 for (int s = 0; s < np_lay; s++)
1467 {
1468 cout << tmpx_lay[s] << " " << tmpy_lay[s] << endl;
1469 }
1470 Orderfunctionx(tmpx_lay, tmpy_lay, tmpx_lay, tmpy_lay);
1471 cout << "fisrt interp" << endl;
1472 for (int s = 0; s < np_lay; s++)
1473 {
1474 cout << tmpx_lay[s] << " " << tmpy_lay[s] << endl;
1475 }
1476
1477 // ASSERTL0(false, "dasd");
1478 // ASSERTL0(tmpx_lay[h*npedge +0]>=0," vert 0 x<0");
1479 // ASSERTL0(tmpx_lay[h*npedge +npedge-1]>0," vert 1 x<0");
1480
1481 // check if the x coord is 'outofbound' and calculate the
1482 // number of outofbound points
1483
1484 // determine which boudn has been overcome:
1485 NekDouble boundleft = xcPhysMOD[0];
1486 NekDouble boundright = xcPhysMOD[np_lay - 1];
1487 bool outboundleft = false;
1488 bool outboundright = false;
1489 if (tmpx_lay[1] < boundleft)
1490 {
1491 outboundleft = true;
1492 }
1493 if (tmpx_lay[np_lay - 2] > boundright)
1494 {
1495 outboundright = true;
1496 }
1497
1498 int outvert = 0;
1499 int outmiddle = 0;
1500 int outcount = 0;
1501
1502 Array<OneD, int> vertout(nvertl);
1503 for (int r = 0; r < nedges; r++)
1504 {
1505 // check point outofboundleft
1506 if (tmpx_lay[r * npedge + npedge - 1] < boundleft &&
1507 outboundleft == true) // assume the neg coords start from 0
1508 {
1509 vertout[outvert] = r;
1510 outvert++;
1511
1512 if (r < nedges - 1)
1513 {
1514 // check if after the last negvert there are neg points
1515 if (tmpx_lay[(r + 1) * npedge + npedge - 1] > boundleft)
1516 {
1517 for (int s = 0; s < npedge - 2; s++)
1518 {
1519 if (tmpx_lay[(r + 1) * npedge + s + 1] <
1520 boundleft)
1521 {
1522 outmiddle++;
1523 }
1524 }
1525 }
1526 }
1527 }
1528 // check point outofboundright
1529 if (tmpx_lay[r * npedge + 0] > boundright &&
1530 outboundright == true) // assume the neg coords start from 0
1531 {
1532 vertout[outvert] = r;
1533 outvert++;
1534
1535 if (r > 0)
1536 {
1537 // check if after the last outvert there are out points
1538 if (tmpx_lay[(r - 1) * npedge + 0] < boundright)
1539 {
1540 for (int s = 0; s < npedge - 2; s++)
1541 {
1542 if (tmpx_lay[(r - 1) * npedge + s + 1] >
1543 boundright)
1544 {
1545 outmiddle++;
1546 }
1547 }
1548 }
1549 }
1550 }
1551 }
1552
1553 // calc number of point to replace
1554 outcount = outvert * npedge + 1 + outmiddle;
1555 // determine from(until) which index the replacement will start
1556 int replacepointsfromindex = 0;
1557 for (int c = 0; c < nedges; c++)
1558 {
1559 // assume at least 1 middle point per edge
1560 if (xcPhysMOD[c * npedge + npedge - 1] <=
1561 tmpx_lay[c * (npedge - (npedge - 2)) + 2] &&
1562 outboundright == true)
1563 {
1564 replacepointsfromindex = c * (npedge - (npedge - 2)) + 2;
1565 break;
1566 }
1567
1568 // assume at least 1 middle point per edge
1569 if (xcPhysMOD[(nedges - 1 - c) * npedge + 0] >=
1570 tmpx_lay[np_lay - 1 -
1571 (c * (npedge - (npedge - 2)) + 2)] &&
1572 outboundleft == true)
1573 {
1574 replacepointsfromindex =
1575 np_lay - 1 - (c * (npedge - (npedge - 2)) + 2);
1576 break;
1577 }
1578 }
1579
1580 // cout<<"out="<<outcount<<endl;
1581 // cout<<"replacefrom="<<replacepointsfromindex<<endl;
1582 // if xcoord is neg find the first positive xcoord
1583
1584 if (outcount > 1)
1585 {
1586 // determine x new coords:
1587 // distribute the point all over the layer
1588 int pstart, shift;
1589 NekDouble increment;
1590
1591 if (outboundright == true)
1592 {
1593 pstart = replacepointsfromindex;
1594 shift = np_lay - outcount;
1595 increment =
1596 (xcPhysMOD[np_lay - outcount] - xcPhysMOD[pstart]) /
1597 (outcount + 1);
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");
1602 }
1603 else
1604 {
1605 shift = 1;
1606 pstart = outcount - 1;
1607 increment = (xcPhysMOD[replacepointsfromindex] -
1608 xcPhysMOD[pstart]) /
1609 (outcount + 1);
1610 ASSERTL0(tmpx_lay[pstart] <
1611 xcPhysMOD[0 * npedge + npedge - 1],
1612 "no middle points in the first edge");
1613 }
1614
1615 // interp to points between posindex and posindex-1
1616 Array<OneD, NekDouble> replace_x(outcount);
1617 Array<OneD, NekDouble> replace_y(outcount);
1618 // order normal functions(cut out repetitions)
1619 Array<OneD, NekDouble> x_tmp(np_lay - (nedges - 1));
1620 Array<OneD, NekDouble> y_tmp(np_lay - (nedges - 1));
1621 Array<OneD, NekDouble> tmpny(np_lay - (nedges - 1));
1622 Cutrepetitions(nedges, xcPhysMOD, x_tmp);
1623 Cutrepetitions(nedges, ycPhysMOD, y_tmp);
1624 Cutrepetitions(nedges, nyPhys, tmpny);
1625 // init neigh arrays
1626 Array<OneD, NekDouble> closex(4);
1627 Array<OneD, NekDouble> closey(4);
1628 Array<OneD, NekDouble> closeny(4);
1629 NekDouble xctmp, ycinterp, nxinterp, nyinterp;
1630
1631 for (int v = 0; v < outcount; v++)
1632 {
1633 xctmp = xcPhysMOD[pstart] + (v + 1) * increment;
1634
1635 // determine closest point index:
1636 int index = DetermineclosePointxindex(xctmp, x_tmp);
1637 // cout<<" vert="<<index<<endl;
1638
1639 // generate neighbour arrays (ny)
1640 GenerateNeighbourArrays(index, 4, x_tmp, tmpny, closex,
1641 closeny);
1642
1643 // interp:
1644 nyinterp = LagrangeInterpolant(xctmp, 4, closex, closeny);
1645
1646 // calc nxinterp
1647 nxinterp = sqrt(abs(1 - nyinterp * nyinterp));
1648
1649 // generata neighbour arrays (yc)
1650 GenerateNeighbourArrays(index, 4, x_tmp, y_tmp, closex,
1651 closey);
1652 // interp:
1653 ycinterp = LagrangeInterpolant(xctmp, 4, closex, closey);
1654 // calc new coord
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];
1659
1660 // cout<<"xinterp="<<replace_x[v]<<"
1661 // yinterp="<<replace_y[v]<<endl;
1662 }
1663 } // end outcount if
1664
1665 /*
1666 for(int s=0; s<np_lay; s++)
1667 {
1668
1669 cout<<tmpx_lay[s]<<" "<<tmpy_lay[s]<<endl;
1670
1671 }
1672 if(m== 0)
1673 {
1674 //ASSERTL0(false, "ssa");
1675 }
1676 */
1677
1678 int closepoints = 4;
1679
1680 Array<OneD, NekDouble> Pyinterp(closepoints);
1681 Array<OneD, NekDouble> Pxinterp(closepoints);
1682
1683 // check if any edge has less than npedge points
1684 int pointscount = 0;
1685 for (int q = 0; q < np_lay; q++)
1686 {
1687 for (int e = 0; e < nedges; e++)
1688 {
1689 if (tmpx_lay[q] <= x_c[e + 1] && tmpx_lay[q] >= x_c[e])
1690 {
1691 pointscount++;
1692 }
1693 if (q == e * npedge + npedge - 1 && pointscount != npedge)
1694 {
1695 // cout<<"edge with few points :"<<e<<endl;
1696 pointscount = 0;
1697 }
1698 else if (q == e * npedge + npedge - 1)
1699 {
1700 pointscount = 0;
1701 }
1702 }
1703 }
1704 //----------------------------------------------------------
1705 /*
1706 cout<<"notordered"<<endl;
1707 for(int g=0; g<tmpx_lay.size(); g++)
1708 {
1709 cout<<tmpx_lay[g]<<" "<<tmpy_lay[g]<<endl;
1710 }
1711 */
1712
1713 // cout<<nedges<<"nedges"<<npedge<<" np_lay="<<np_lay<<endl;
1714 // calc lay coords
1715 // MoveLayerNfixedxpos(nvertl, npedge, xcPhysMOD, tmpx_lay,
1716 // tmpy_lay, lay_Vids[m], layers_x[m],
1717 // layers_y[m],xnew,ynew);
1718 MoveLayerNnormpos(nvertl, npedge, xcPhysMOD, tmpx_lay, tmpy_lay,
1719 lay_Vids[m], layers_x[m], layers_y[m], xnew,
1720 ynew);
1721
1722 /*
1723 //generate x,y arrays without lastedgepoint
1724 //(needed to interp correctly)
1725 Array<OneD, NekDouble>tmpx(np_lay-(nedges-1));
1726 Array<OneD, NekDouble>tmpy(np_lay-(nedges-1));
1727
1728 Cutrepetitions(nedges, tmpx_lay, tmpx);
1729 Cutrepetitions(nedges, tmpy_lay, tmpy);
1730 //order points in x:
1731 int index;
1732 Array<OneD, NekDouble> copyarray_x(tmpx.size());
1733 Array<OneD, NekDouble> copyarray_y(tmpx.size());
1734 Orderfunctionx(tmpx, tmpy, tmpx, tmpy);
1735
1736 //determine the neighbour points (-3;+3)
1737 for(int g=0; g< nvertl; g++)
1738 {
1739 //verts
1740 //cout<<"determine value for vert x="<<x_c[g]<<endl;
1741 //determine closest index:
1742
1743 index=
1744 DetermineclosePointxindex( x_c[g], tmpx);
1745 //generate neighbour arrays:
1746 GenerateNeighbourArrays(index,
1747 closepoints,tmpx,tmpy,Pxinterp,Pyinterp);
1748 //write vert coords
1749 ynew[lay_Vids[m][g] ]=
1750 LagrangeInterpolant(x_c[g],closepoints,Pxinterp,Pyinterp );
1751 xnew[lay_Vids[m][g] ]= x_c[g];
1752
1753
1754 if(g<nedges)
1755 {
1756
1757 //v1
1758 layers_y[m][g*npedge +0] = ynew[lay_Vids[m][g] ];
1759 layers_x[m][g*npedge +0] = xnew[lay_Vids[m][g] ];
1760 //v2
1761
1762 //determine closest index:
1763 index=
1764 DetermineclosePointxindex( x_c[g+1], tmpx);
1765 //generate neighbour arrays:
1766 GenerateNeighbourArrays(index,
1767 closepoints,tmpx,tmpy,Pxinterp,Pyinterp); layers_y[m][g*npedge
1768 +npedge-1] =
1769 LagrangeInterpolant(x_c[g+1],closepoints,Pxinterp,Pyinterp );
1770 layers_x[m][g*npedge +npedge-1] = x_c[g+1];
1771
1772 //middle points
1773 for(int r=0; r< npedge-2; r++)
1774 {
1775 //determine closest point index:
1776 index =
1777 DetermineclosePointxindex( xcPhysMOD[g*npedge +r+1], tmpx);
1778 //cout<<" vert+"<<index<<endl;
1779
1780 ASSERTL0( index<= tmpy.size()-1, " index wrong");
1781 //generate neighbour arrays Pyinterp,Pxinterp
1782 GenerateNeighbourArrays(index,
1783 closepoints,tmpx,tmpy,Pxinterp,Pyinterp);
1784
1785 layers_y[m][g*npedge +r+1]=
1786 LagrangeInterpolant(
1787 xcPhysMOD[g*npedge +r+1],closepoints,Pxinterp,Pyinterp );
1788 //cout<<"x value="<<xcPhysMOD[g*npedge +r+1]<<endl;
1789 layers_x[m][g*npedge +r+1]= xcPhysMOD[g*npedge +r+1];
1790
1791 }
1792
1793
1794 }//if edge closed g
1795 }// nvertl closed
1796 */
1797 // check if there are points out of range:
1798 // cout<<Vmath::Vmax(np_lay,layers_y[m],1)<<endl;
1799 if (curv_lay == true)
1800 {
1801 // ASSERTL0(Vmath::Vmax(np_lay,layers_y[m],1)<
1802 // Vmath::Vmax(nVertTot,yold,1),"point>ymax");
1803 // ASSERTL0(Vmath::Vmin(np_lay,layers_y[m],1)>
1804 // Vmath::Vmin(nVertTot,yold,1),"point<ymin");
1805 }
1806
1807 // force Polycontinuity of the layer(3 order poly every 2 edges):
1808 int npoints = npedge;
1809 Array<OneD, NekDouble> xPedges(npoints);
1810 Array<OneD, NekDouble> yPedges(npoints);
1811 for (int f = 0; f < nedges; f++)
1812 {
1813 int polyorder = 2;
1814
1815 Array<OneD, NekDouble> coeffsinterp(polyorder + 1);
1816
1817 Vmath::Vcopy(npoints, &layers_x[m][(f)*npedge + 0], 1,
1818 &xPedges[0], 1);
1819 Vmath::Vcopy(npoints, &layers_y[m][(f)*npedge + 0], 1,
1820 &yPedges[0], 1);
1821
1822 PolyFit(polyorder, npoints, xPedges, yPedges, coeffsinterp,
1823 xPedges, yPedges, npoints);
1824
1825 // copy back the values:
1826 Vmath::Vcopy(npoints, &yPedges[0], 1,
1827 &layers_y[m][(f)*npedge + 0], 1);
1828
1829 // verts still the same:
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]];
1832 }
1833
1834 cout << " xlay ylay lay:" << m << endl;
1835 for (int l = 0; l < np_lay; l++)
1836 {
1837 // cout<<tmpx_lay[l]<<" "<<tmpy_lay[l]<<endl;
1838 cout << std::setprecision(8) << layers_x[m][l] << " "
1839 << layers_y[m][l] << endl;
1840 }
1841 /*
1842 cout<<"nverts"<<endl;
1843 for(int l=0; l<nvertl; l++)
1844 {
1845 cout<<std::setprecision(8)<<xnew[lay_Vids[m][l] ]<<"
1846 "<<ynew[lay_Vids[m][l] ]<<endl;
1847 }
1848 */
1849
1850 // ASSERTL0(false, "as");
1851
1852 // if the layers coords are too steep use two edges verts to get an
1853 // poly interp third order
1854 /*
1855 bool polyinterp=false;
1856 for(int b=0; b< nedges; b++)
1857 {
1858 for(int u=0; u<npedge; u++)
1859 {
1860 if(
1861 abs(layers_y[m][b*npedge+u+1]- layers_y[m][b*npedge
1862 +u])/(layers_x[m][b*npedge+u+1]-layers_x[m][b*npedge+u]))>4.0 )
1863 {
1864 polyinterp=true;
1865 break;
1866 }
1867 cout<<"incratio="<<incratio<<endl;
1868 }
1869 //
1870 //Evaluatelay_tan
1871
1872 }
1873 */
1874
1875 cout << "lay=" << m << endl;
1876 ASSERTL0(Vmath::Vmin(nVertTot, yold, 1) <
1877 Vmath::Vmin(np_lay, layers_y[m], 1),
1878 " different layer ymin val");
1879 ASSERTL0(Vmath::Vmax(nVertTot, yold, 1) >
1880 Vmath::Vmax(np_lay, layers_y[m], 1),
1881 " different layer ymax val");
1882 ASSERTL0(Vmath::Vmin(nVertTot, xold, 1) ==
1883 Vmath::Vmin(np_lay, layers_x[m], 1),
1884 " different layer xmin val");
1885 ASSERTL0(Vmath::Vmax(nVertTot, xold, 1) ==
1886 Vmath::Vmax(np_lay, layers_x[m], 1),
1887 " different layer xmax val");
1888 }
1889
1890 // MoveOutsidePointsfixedxpos(npedge, graphShPt,xold_c, yold_c,
1891 // xold_low, yold_low, xold_up, yold_up, layers_y[0],
1892 // layers_y[nlays-1],
1893 // xnew, ynew);
1894 // lastIregion -1 = laydown
1895 // lastIregion -2 = layup
1896 MoveOutsidePointsNnormpos(
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);
1900
1901 /*
1902 //update vertices coords outside layers region
1903 NekDouble ratio;
1904 for(int n=0; n<nVertTot; n++)
1905 {
1906 NekDouble ratio;
1907 SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(n);
1908 NekDouble x,y,z;
1909 vertex->GetCoords(x,y,z);
1910 int qp_closer;
1911 NekDouble diff;
1912 int qp_closerup, qp_closerdown;
1913 NekDouble diffup, diffdown;
1914 //determine the closer xold_up
1915 NekDouble tmp=1000;
1916 diffup =1000;
1917 diffdown = 1000;
1918 for(int k=0; k<nvertl; k++)
1919 {
1920 if(abs(x-xold_c[k]) < tmp)
1921 {
1922 tmp = abs(x-xold_c[k]);
1923 qp_closer=k;
1924 }
1925 }
1926
1927 //find nplay_closer
1928 int nplay_closer;
1929 if(qp_closer==0)
1930 {
1931 nplay_closer=0;//first vert
1932 }
1933 else
1934 {
1935 nplay_closer= (qp_closer-1)*npedge +npedge-1;
1936 }
1937
1938 if( y>yold_up[qp_closer] && y<1 )//nlays-1 is layerup
1939 {
1940
1941 // ratio =
1942 (1-layers_y[nlays-1][qp_closer])*(1-y_c[qp_closer])/
1943 // ( (1-yold_up[n])*(1-yold_c[qp_closer]) );
1944 ratio = (1-layers_y[nlays-1][nplay_closer])/
1945 ( (1-yold_up[qp_closer]) );
1946 //distance prop to layerup
1947 ynew[n] = layers_y[nlays-1][nplay_closer]
1948 + (y-yold_up[qp_closer])*ratio;
1949 xnew[n] = x;
1950
1951 }
1952 else if( y< yold_low[qp_closer] && y>-1 )//0 is layerdown
1953 {
1954
1955 ratio = (1+layers_y[0][nplay_closer])/
1956 ( (1+yold_low[qp_closer]) );
1957 //distance prop to layerlow
1958 ynew[n] = layers_y[0][nplay_closer]
1959 + (y-yold_low[qp_closer])*ratio;
1960 xnew[n] = x;
1961 }
1962
1963 }
1964 */
1965
1966 /*
1967 //update verts coords of critlay(in case EvaluateLayTnaget has been
1968 called)
1969 //it's not necessary !!!
1970 for(int e=0; e<nvertl; e++)
1971 {
1972 ynew[CritLay[e]] = y_c[e];
1973 }
1974 */
1975
1976 } // move_norm bool
1977 else // move vertically
1978 {
1979 MoveLayersvertically(nlays, nvertl, cntlow, cntup, lay_Vids, x_c, y_c,
1980 Down, Up, xnew, ynew, layers_x, layers_y);
1981 }
1982
1983 // check singular quads:
1984 CheckSingularQuads(streak, V1, V2, xnew, ynew);
1985 // check borders of the new mesh verts:
1986 // cout<<std::setprecision(8)<<"yoldmax="<<Vmath::Vmax(nVertTot,
1987 // yold,1)<<endl;
1988 // cout<<std::setprecision(8)<<"ynewmax="<<Vmath::Vmax(nVertTot,ynew,1)<<endl;
1989
1990 cout << std::setprecision(8) << "xmin=" << Vmath::Vmin(nVertTot, xnew, 1)
1991 << endl;
1992 ASSERTL0(Vmath::Vmin(nVertTot, xold, 1) == Vmath::Vmin(nVertTot, xnew, 1),
1993 " different xmin val");
1994 ASSERTL0(Vmath::Vmin(nVertTot, yold, 1) == Vmath::Vmin(nVertTot, ynew, 1),
1995 " different ymin val");
1996 ASSERTL0(Vmath::Vmax(nVertTot, xold, 1) == Vmath::Vmax(nVertTot, xnew, 1),
1997 " different xmax val");
1998 ASSERTL0(Vmath::Vmax(nVertTot, yold, 1) == Vmath::Vmax(nVertTot, ynew, 1),
1999 " different ymax val");
2000
2001 // replace the vertices with the new ones
2002
2003 Replacevertices(changefile, xnew, ynew, xcPhys, ycPhys, Eids, npedge,
2004 charalp, layers_x, layers_y, lay_Eids, curv_lay);
2005}
GenerateMapEidsv1v2
void GenerateMapEidsv1v2(MultiRegions::ExpListSharedPtr field, Array< OneD, int > &V1, Array< OneD, int > &V2)
Definition: MeshMove.cpp:2319
MoveLayerNnormpos
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)
Definition: MeshMove.cpp:3308
EvaluateTangent
void EvaluateTangent(int npoints, Array< OneD, NekDouble > xcQedge, Array< OneD, NekDouble > coeffsinterp, Array< OneD, NekDouble > &txQedge, Array< OneD, NekDouble > &tyQedge)
Definition: MeshMove.cpp:2893
LagrangeInterpolant
NekDouble LagrangeInterpolant(NekDouble x, int npts, Array< OneD, NekDouble > xpts, Array< OneD, NekDouble > funcvals)
Definition: MeshMove.cpp:2865
Computestreakpositions
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)
Definition: MeshMove.cpp:2102
Cutrepetitions
void Cutrepetitions(int nedges, Array< OneD, NekDouble > inarray, Array< OneD, NekDouble > &outarray)
Definition: MeshMove.cpp:2768
MoveOutsidePointsNnormpos
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)
Definition: MeshMove.cpp:3435
MappingEVids
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)
Definition: MeshMove.cpp:2382
CheckSingularQuads
void CheckSingularQuads(MultiRegions::ExpListSharedPtr Exp, Array< OneD, int > V1, Array< OneD, int > V2, Array< OneD, NekDouble > &xnew, Array< OneD, NekDouble > &ynew)
Definition: MeshMove.cpp:3713
Orderfunctionx
void Orderfunctionx(Array< OneD, NekDouble > inarray_x, Array< OneD, NekDouble > inarray_y, Array< OneD, NekDouble > &outarray_x, Array< OneD, NekDouble > &outarray_y)
Definition: MeshMove.cpp:3138
Replacevertices
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)
Definition: MeshMove.cpp:3893
OrderVertices
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)
Definition: MeshMove.cpp:2007
PolyFit
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)
Definition: MeshMove.cpp:3029
GenerateAddPointsNewtonIt
void GenerateAddPointsNewtonIt(NekDouble xi, NekDouble yi, NekDouble &xout, NekDouble &yout, MultiRegions::ExpListSharedPtr function, Array< OneD, NekDouble > derfunction, NekDouble cr)
Definition: MeshMove.cpp:2229
MoveLayersvertically
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)
Definition: MeshMove.cpp:3183
GenerateNeighbourArrays
void GenerateNeighbourArrays(int index, int neighpoints, Array< OneD, NekDouble > xArray, Array< OneD, NekDouble > yArray, Array< OneD, NekDouble > &Neighbour_x, Array< OneD, NekDouble > &Neighbour_y)
Definition: MeshMove.cpp:2813
DetermineclosePointxindex
int DetermineclosePointxindex(NekDouble x, Array< OneD, NekDouble > xArray)
Definition: MeshMove.cpp:2795
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:81
Nektar::StdRegions::StdExpansion::GetCoords
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
Definition: StdExpansion.h:579
Nektar::LibUtilities::Import
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > &ElementIDs)
This function allows for data to be imported from an FLD file when a session and/or communicator is n...
Definition: FieldIO.cpp:288
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:113
Nektar::LocalRegions::SegExpSharedPtr
std::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:248
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
Nektar::MultiRegions::ContFieldSharedPtr
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:268
Nektar::SpatialDomains::BoundaryRegionCollection
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:210
Nektar::SpatialDomains::BoundaryConditionsSharedPtr
std::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
Definition: Conditions.h:289
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
Nektar::SpatialDomains::PointGeomSharedPtr
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
Nektar::UnitTests::z
std::vector< double > z(NPUPPER)
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
Vmath::Vmul
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.
Definition: Vmath.hpp:72
Vmath::Vmin
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.hpp:725
Vmath::Vvtvp
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
Definition: Vmath.hpp:366
Vmath::Sdiv
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/x.
Definition: Vmath.hpp:154
Vmath::Vdiv
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.
Definition: Vmath.hpp:126
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.hpp:644
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

References tinysimd::abs(), ASSERTL0, CheckSingularQuads(), Computestreakpositions(), Cutrepetitions(), Nektar::UnitTests::d(), DetermineclosePointxindex(), EvaluateTangent(), Nektar::MultiRegions::eX, Nektar::MultiRegions::eY, GenerateAddPointsNewtonIt(), GenerateMapEidsv1v2(), GenerateNeighbourArrays(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), Nektar::StdRegions::StdExpansion::GetCoords(), Nektar::LibUtilities::Import(), LagrangeInterpolant(), MappingEVids(), MoveLayerNnormpos(), MoveLayersvertically(), MoveOutsidePointsNnormpos(), Orderfunctionx(), OrderVertices(), PolyFit(), Nektar::UnitTests::q(), Replacevertices(), Vmath::Sadd(), Vmath::Sdiv(), Vmath::Smul(), tinysimd::sqrt(), Vmath::Vcopy(), Vmath::Vdiv(), Vmath::Vmax(), Vmath::Vmin(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vvtvp(), Nektar::UnitTests::w(), Nektar::UnitTests::z(), and Vmath::Zero().

◆ MappingEVids()

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 
)

Definition at line 2382 of file MeshMove.cpp.

2392{
2393 int nlay_Eids = xcold.size() - 1;
2394 int nlay_Vids = xcold.size();
2395
2396 int nVertsTot = mesh->GetNvertices();
2397 cout << "nverttot=" << nVertsTot << endl;
2398 // find the first vert of each layer
2399 // int qp_closerup,qp_closerdown;
2400 // NekDouble diffup,diffdown;
2401 cout << "init nlays=" << nlays << endl;
2402 // tmp first vert array (of course <100!!)
2403 Array<OneD, int> tmpVids0(100);
2404 Array<OneD, NekDouble> tmpx0(100);
2405 Array<OneD, NekDouble> tmpy0(100);
2406 Array<OneD, NekDouble> x2D(nVertsTot);
2407 Array<OneD, NekDouble> y2D(nVertsTot);
2408 cout << "yoldup=" << yoldup[0] << endl;
2409 cout << "yolddown=" << yolddown[0] << endl;
2410
2411 for (int r = 0; r < nVertsTot; r++)
2412 {
2413
2414 SpatialDomains::PointGeomSharedPtr vertex = mesh->GetVertex(r);
2415 NekDouble x, y, z;
2416 vertex->GetCoords(x, y, z);
2417 x2D[r] = x;
2418 y2D[r] = y;
2419 if (xcold[0] == x)
2420 {
2421 // check if the vert is inside layer zone
2422 if (y <= yoldup[0] && y >= yolddown[0] && y != ycold[0])
2423 {
2424 // cout<<"x="<<x<<" y="<<y<<endl;
2425 tmpVids0[nlays] = r;
2426 tmpx0[nlays] = x;
2427 tmpy0[nlays] = y;
2428 nlays++;
2429 }
2430 }
2431 }
2432 cout << "nlays=" << nlays << endl;
2433
2434 // reorder first verts from xlower to xhigher
2435 Array<OneD, NekDouble> tmpx(nlays);
2436 Array<OneD, NekDouble> tmpf(nlays);
2437 Array<OneD, int> tmpV(nlays);
2438 // local copy to prevent overwriting
2439 Vmath::Vcopy(nlays, tmpx0, 1, tmpx, 1);
2440 Vmath::Vcopy(nlays, tmpy0, 1, tmpf, 1);
2441 Vmath::Vcopy(nlays, tmpVids0, 1, tmpV, 1);
2442 NekDouble max = Vmath::Vmax(tmpf.size(), tmpf, 1);
2443 int index = 0;
2444 for (int w = 0; w < nlays; w++)
2445 {
2446 index = Vmath::Imin(tmpf.size(), tmpf, 1);
2447 tmpx0[w] = tmpx[index];
2448 tmpy0[w] = tmpf[index];
2449 tmpVids0[w] = tmpV[index];
2450 tmpf[index] = max + 1000;
2451 }
2452
2453 // initialise layers arrays
2454 Eids_lay = Array<OneD, Array<OneD, int>>(nlays);
2455 Vids_lay = Array<OneD, Array<OneD, int>>(nlays);
2456 for (int m = 0; m < nlays; m++)
2457 {
2458 Eids_lay[m] = Array<OneD, int>(nlay_Eids);
2459 Vids_lay[m] = Array<OneD, int>(nlay_Vids);
2460 }
2461
2462 // determine the others verts and edge for each layer
2463 NekDouble normbef = 0.0;
2464 NekDouble normtmp = 0.0;
2465 NekDouble xbef = 0.0;
2466 NekDouble ybef = 0.0;
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;
2469 Array<OneD, NekDouble> coord(2);
2470 int elmtid, offset;
2471 int nTotEdges = V1.size();
2472 Array<OneD, int> edgestmp(6);
2473 for (int m = 0; m < nlays; m++)
2474 {
2475 for (int g = 0; g < nlay_Eids; g++)
2476 {
2477 if (g == 0)
2478 {
2479 for (int h = 0; h < nTotEdges; h++)
2480 {
2481
2482 if (tmpVids0[m] == V1[h])
2483 {
2484 SpatialDomains::PointGeomSharedPtr vertex =
2485 mesh->GetVertex(V2[h]);
2486 NekDouble x, y, z;
2487 vertex->GetCoords(x, y, z);
2488 if (x != tmpx0[m])
2489 {
2490 Eids_lay[m][0] = h;
2491 Vids_lay[m][0] = V1[h];
2492 Vids_lay[m][1] = V2[h];
2493 SpatialDomains::PointGeomSharedPtr vertex1 =
2494 mesh->GetVertex(V1[h]);
2495 NekDouble x1, y1, z1;
2496 vertex1->GetCoords(x1, y1, z1);
2497 normbef =
2498 sqrt((y - y1) * (y - y1) + (x - x1) * (x - x1));
2499 ybef = (y - y1);
2500 xbef = (x - x1);
2501 coord[0] = x;
2502 coord[1] = y;
2503 elmtid = streak->GetExpIndex(coord, 0.00001);
2504 offset = streak->GetPhys_Offset(elmtid);
2505 Ubef = streak->GetExp(elmtid)->PhysEvaluate(
2506 coord, streak->GetPhys() + offset);
2507 }
2508 }
2509 if (tmpVids0[m] == V2[h])
2510 {
2511 SpatialDomains::PointGeomSharedPtr vertex =
2512 mesh->GetVertex(V1[h]);
2513 NekDouble x, y, z;
2514 vertex->GetCoords(x, y, z);
2515 if (x != tmpx0[m])
2516 {
2517 Eids_lay[m][0] = h;
2518 Vids_lay[m][0] = V2[h];
2519 Vids_lay[m][1] = V1[h];
2520 SpatialDomains::PointGeomSharedPtr vertex2 =
2521 mesh->GetVertex(V2[h]);
2522 NekDouble x2 = 0.0, y2 = 0.0;
2523 normbef =
2524 sqrt((y - y2) * (y - y2) + (x - x2) * (x - x2));
2525 ybef = (y - y2);
2526 xbef = (x - x2);
2527 coord[0] = x;
2528 coord[1] = y;
2529 elmtid = streak->GetExpIndex(coord, 0.00001);
2530 offset = streak->GetPhys_Offset(elmtid);
2531 Ubef = streak->GetExp(elmtid)->PhysEvaluate(
2532 coord, streak->GetPhys() + offset);
2533 }
2534 }
2535 }
2536
2537 cout << "Eid=" << Eids_lay[m][0]
2538 << " Vids_lay0=" << Vids_lay[m][0]
2539 << " Vidslay1=" << Vids_lay[m][1] << endl;
2540 }
2541 else
2542 {
2543 // three verts(edges) candidates
2544 int cnt = 0;
2545 for (int h = 0; h < nTotEdges; h++)
2546 {
2547 // cout<<Vids_lay[m][g]<<" V1="<<V1[h]<<"
2548 // V2[h]="<<V2[h]<<endl;
2549 if ((Vids_lay[m][g] == V1[h] || Vids_lay[m][g] == V2[h]) &&
2550 h != Eids_lay[m][g - 1])
2551 {
2552 cout << "edgetmp=" << h << endl;
2553 ASSERTL0(cnt <= 6, "wrong number of candidates");
2554 edgestmp[cnt] = h;
2555 cnt++;
2556 }
2557 }
2558
2559 diff = 1000;
2560 Array<OneD, NekDouble> diffarray(cnt);
2561 Array<OneD, NekDouble> diffUarray(cnt);
2562 cout << "normbef=" << normbef << endl;
2563 cout << "Ubefcc=" << Ubef << endl;
2564 // choose the right candidate
2565 for (int e = 0; e < cnt; e++)
2566 {
2567 SpatialDomains::PointGeomSharedPtr vertex1 =
2568 mesh->GetVertex(V1[edgestmp[e]]);
2569 NekDouble x1, y1, z1;
2570 vertex1->GetCoords(x1, y1, z1);
2571 SpatialDomains::PointGeomSharedPtr vertex2 =
2572 mesh->GetVertex(V2[edgestmp[e]]);
2573 NekDouble x2, y2, z2;
2574 vertex2->GetCoords(x2, y2, z2);
2575
2576 normtmp =
2577 sqrt((y2 - y1) * (y2 - y1) + (x2 - x1) * (x2 - x1));
2578
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]])
2583 {
2584
2585 ytmp = (y2 - y1);
2586 xtmp = (x2 - x1);
2587 coord[0] = x2;
2588 coord[1] = y2;
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) -
2595 1);
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)) -
2602 1)
2603 << endl;
2604 if (abs((xtmp * xbef + ytmp * ybef) /
2605 (normtmp * normbef) -
2606 1) < diff &&
2607 y2 <= yoldup[g + 1] && y2 >= yolddown[g + 1] &&
2608 y1 <= yoldup[g] && y1 >= yolddown[g])
2609 {
2610
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) -
2615 1);
2616 normnext = normtmp;
2617 ynext = ytmp;
2618 xnext = xtmp;
2619 Unext = Utmp;
2620 }
2621 }
2622 else if (Vids_lay[m][g] == V2[edgestmp[e]])
2623 {
2624
2625 ytmp = (y1 - y2);
2626 xtmp = (x1 - x2);
2627 coord[0] = x1;
2628 coord[1] = y1;
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) -
2635 1);
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)) -
2642 1)
2643 << endl;
2644 if (abs((xtmp * xbef + ytmp * ybef) /
2645 (normtmp * normbef) -
2646 1) < diff &&
2647 y2 <= yoldup[g] && y2 >= yolddown[g] &&
2648 y1 <= yoldup[g + 1] && y1 >= yolddown[g + 1])
2649 {
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) -
2654 1);
2655 normnext = normtmp;
2656 ynext = ytmp;
2657 xnext = xtmp;
2658 Unext = Utmp;
2659 }
2660 }
2661 else
2662 {
2663 ASSERTL0(false, "eid not found");
2664 }
2665 }
2666 cout << "Eid before check=" << Eids_lay[m][g] << endl;
2667 for (int q = 0; q < cnt; q++)
2668 {
2669 cout << q << " diff" << diffarray[q] << endl;
2670 }
2671 // check if the eid has a vert in common with another layer
2672 bool check = false;
2673 if (m > 0 && m < nlays)
2674 {
2675 check = checkcommonvert(Vids_lay[m - 1], Vids_c,
2676 Vids_lay[m][g + 1]);
2677 }
2678 if (check == true)
2679 {
2680 cout << "COMMON VERT" << endl;
2681 int eid = Vmath::Imin(cnt, diffarray, 1);
2682 diffarray[eid] = 1000;
2683 eid = Vmath::Imin(cnt, diffarray, 1);
2684
2685 SpatialDomains::PointGeomSharedPtr vertex1 =
2686 mesh->GetVertex(V1[edgestmp[eid]]);
2687 NekDouble x1, y1, z1;
2688 vertex1->GetCoords(x1, y1, z1);
2689 SpatialDomains::PointGeomSharedPtr vertex2 =
2690 mesh->GetVertex(V2[edgestmp[eid]]);
2691 NekDouble x2, y2, z2;
2692 vertex2->GetCoords(x2, y2, z2);
2693
2694 normtmp =
2695 sqrt((y2 - y1) * (y2 - y1) + (x2 - x1) * (x2 - x1));
2696
2697 Eids_lay[m][g] = edgestmp[eid];
2698 if (Vids_lay[m][g] == V1[edgestmp[eid]])
2699 {
2700 coord[0] = x2;
2701 coord[1] = y2;
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]];
2707 normnext = normtmp;
2708 ynext = (y2 - y1);
2709 xnext = (x2 - x1);
2710 ;
2711 Unext = Utmp;
2712 }
2713 if (Vids_lay[m][g] == V2[edgestmp[eid]])
2714 {
2715 coord[0] = x1;
2716 coord[1] = y1;
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]];
2722 normnext = normtmp;
2723 ynext = (y1 - y2);
2724 xnext = (x1 - x2);
2725 ;
2726 Unext = Utmp;
2727 }
2728 }
2729
2730 cout << m << "edge aft:" << Eids_lay[m][g]
2731 << " Vid=" << Vids_lay[m][g + 1] << endl;
2732 normbef = normnext;
2733 ybef = ynext;
2734 xbef = xnext;
2735 Ubef = Unext;
2736
2737 cout << "endelse" << normtmp << endl;
2738 } // end else
2739
2740 } // end g it
2741
2742 } // end m it
2743
2744 for (int w = 0; w < nlays; w++)
2745 {
2746 for (int f = 0; f < nlay_Eids; f++)
2747 {
2748 cout << "check=" << w << " Eid:" << Eids_lay[w][f] << endl;
2749 }
2750 }
2751}
checkcommonvert
bool checkcommonvert(Array< OneD, int > Vids_laybefore, Array< OneD, int > Vids_c, int Vid)
Definition: MeshMove.cpp:2753

References tinysimd::abs(), ASSERTL0, checkcommonvert(), Vmath::Imin(), Nektar::UnitTests::q(), tinysimd::sqrt(), Vmath::Vcopy(), Vmath::Vmax(), Nektar::UnitTests::w(), and Nektar::UnitTests::z().

Referenced by main().

◆ MoveLayerNfixedxpos()

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 
)

Definition at line 3231 of file MeshMove.cpp.

3238{
3239 int np_lay = xcPhys.size();
3240 int nedges = nvertl - 1;
3241 Array<OneD, NekDouble> tmpx(np_lay - (nedges - 1));
3242 Array<OneD, NekDouble> tmpy(np_lay - (nedges - 1));
3243 Cutrepetitions(nedges, tmpx_lay, tmpx);
3244 Cutrepetitions(nedges, tmpy_lay, tmpy);
3245 // order points in x:
3246 int index;
3247 int closepoints = 4;
3248 Array<OneD, NekDouble> Pxinterp(closepoints);
3249 Array<OneD, NekDouble> Pyinterp(closepoints);
3250 Orderfunctionx(tmpx, tmpy, tmpx, tmpy);
3251 // determine the neighbour points (-3;+3)
3252 for (int g = 0; g < nedges; g++)
3253 {
3254 // write vert coords
3255 // v1
3256 index = DetermineclosePointxindex(xcPhys[g * npedge + 0], tmpx);
3257 // generate neighbour arrays:
3258 GenerateNeighbourArrays(index, closepoints, tmpx, tmpy, Pxinterp,
3259 Pyinterp);
3260 ynew[Vids[g]] = LagrangeInterpolant(xcPhys[g * npedge + 0], closepoints,
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]];
3265
3266 // v2
3267 // determine closest index:
3268 index =
3269 DetermineclosePointxindex(xcPhys[g * npedge + npedge - 1], tmpx);
3270 // generate neighbour arrays:
3271 GenerateNeighbourArrays(index, closepoints, tmpx, tmpy, Pxinterp,
3272 Pyinterp);
3273 ynew[Vids[g + 1]] = LagrangeInterpolant(
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]];
3278
3279 // middle points
3280 for (int r = 0; r < npedge - 2; r++)
3281 {
3282
3283 // determine closest point index:
3284 index = DetermineclosePointxindex(xcPhys[g * npedge + r + 1], tmpx);
3285 // cout<<" vert+"<<index<<endl;
3286
3287 ASSERTL0(index <= tmpy.size() - 1, " index wrong");
3288 // generate neighbour arrays Pyinterp,Pxinterp
3289 GenerateNeighbourArrays(index, closepoints, tmpx, tmpy, Pxinterp,
3290 Pyinterp);
3291
3292 ylay[g * npedge + r + 1] = LagrangeInterpolant(
3293 xcPhys[g * npedge + r + 1], closepoints, Pxinterp, Pyinterp);
3294 // cout<<"x value="<<xcPhysMOD[g*npedge +r+1]<<endl;
3295 xlay[g * npedge + r + 1] = xcPhys[g * npedge + r + 1];
3296
3297 /*
3298 for(int t=0; t<6; t++)
3299 {
3300 cout<<"Px="<<Pxinterp[t]<<" "<<Pyinterp[t]<<endl;
3301 }
3302 */
3303 }
3304
3305 } // edge closed g
3306}

References ASSERTL0, Cutrepetitions(), DetermineclosePointxindex(), GenerateNeighbourArrays(), LagrangeInterpolant(), and Orderfunctionx().

◆ MoveLayerNnormpos()

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 
)

Definition at line 3308 of file MeshMove.cpp.

3315{
3316 int np_lay = xcPhys.size();
3317 int nedges = nvertl - 1;
3318 NekDouble x0, x1, xtmp;
3319 Array<OneD, NekDouble> tmpx(np_lay - (nedges - 1));
3320 Array<OneD, NekDouble> tmpy(np_lay - (nedges - 1));
3321 Cutrepetitions(nedges, tmpx_lay, tmpx);
3322 Cutrepetitions(nedges, tmpy_lay, tmpy);
3323 // order points in x:
3324 int index;
3325 int closepoints = 4;
3326 Array<OneD, NekDouble> Pxinterp(closepoints);
3327 Array<OneD, NekDouble> Pyinterp(closepoints);
3328 Orderfunctionx(tmpx, tmpy, tmpx, tmpy);
3329 // determine the neighbour points (-3;+3)
3330 for (int g = 0; g < nedges; g++)
3331 {
3332 // write vert coords
3333 // v1
3334 ynew[Vids[g]] = tmpy_lay[g * npedge + 0];
3335 xnew[Vids[g]] = tmpx_lay[g * npedge + 0];
3336
3337 ylay[g * npedge + 0] = ynew[Vids[g]];
3338 xlay[g * npedge + 0] = xnew[Vids[g]];
3339
3340 // v2
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]];
3345
3346 // middle points
3347 for (int r = 0; r < npedge - 2; r++)
3348 {
3349 x0 = xlay[g * npedge + 0];
3350 x1 = xlay[g * npedge + npedge - 1];
3351 xtmp = x0 + r * (x1 - x0) / (npedge - 1);
3352 // determine closest point index:
3353 index = DetermineclosePointxindex(xtmp, tmpx);
3354 // cout<<" vert+"<<index<<endl;
3355
3356 ASSERTL0(index <= tmpy.size() - 1, " index wrong");
3357 // generate neighbour arrays Pyinterp,Pxinterp
3358 GenerateNeighbourArrays(index, closepoints, tmpx, tmpy, Pxinterp,
3359 Pyinterp);
3360
3361 ylay[g * npedge + r + 1] =
3362 LagrangeInterpolant(xtmp, closepoints, Pxinterp, Pyinterp);
3363 // cout<<"x value="<<xcPhysMOD[g*npedge +r+1]<<endl;
3364 xlay[g * npedge + r + 1] = xtmp;
3365 }
3366
3367 } // edge closed g
3368}

References ASSERTL0, Cutrepetitions(), DetermineclosePointxindex(), GenerateNeighbourArrays(), LagrangeInterpolant(), and Orderfunctionx().

Referenced by main().

◆ MoveLayersvertically()

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 
)

Definition at line 3183 of file MeshMove.cpp.

3191{
3192 int np_lay = layers_y[0].size();
3193 // 0<h<nlays-1 to fill only the 'internal' layers(no up,low);
3194 for (int h = 1; h < nlays - 1; h++)
3195 {
3196 layers_x[h] = Array<OneD, NekDouble>(np_lay);
3197 for (int s = 0; s < nvertl; s++)
3198 {
3199 // check if ynew is still empty
3200 ASSERTL0(ynew[lay_Vids[h][s]] == -20, "ynew layers not empty");
3201 if (h < cntlow + 1) // layers under the crit lay
3202 {
3203 // y= ylow+delta
3204 ynew[lay_Vids[h][s]] =
3205 ynew[Down[s]] +
3206 h * abs(ynew[Down[s]] - yc[s]) / (cntlow + 1);
3207 // put the layer vertical
3208 xnew[lay_Vids[h][s]] = xc[s];
3209 // cout<<"ynew="<<ynew[ lay_Vids[h][s] ]<<"
3210 // ydown="<<ynew[Down[s]]<< " delta="<<abs(ynew[Down[s]] -
3211 // y_c[s])/(cntlow+1)<<endl; until now layers_y=yold
3212 layers_y[h][s] = ynew[lay_Vids[h][s]];
3213 layers_x[h][s] = xnew[lay_Vids[h][s]];
3214 }
3215 else
3216 {
3217 // y = yc+delta
3218 ynew[lay_Vids[h][s]] = yc[s] + (h - cntlow) *
3219 abs(ynew[Up[s]] - yc[s]) /
3220 (cntup + 1);
3221 // put the layer vertical
3222 xnew[lay_Vids[h][s]] = xc[s];
3223 // until now layers_y=yold
3224 layers_y[h][s] = ynew[lay_Vids[h][s]];
3225 layers_x[h][s] = xnew[lay_Vids[h][s]];
3226 }
3227 }
3228 }
3229}

References tinysimd::abs(), and ASSERTL0.

Referenced by main().

◆ MoveOutsidePointsfixedxpos()

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 
)

Definition at line 3370 of file MeshMove.cpp.

3379{
3380 // update vertices coords outside layers region
3381 int nvertl = ycold.size();
3382 int nVertTot = mesh->GetNvertices();
3383 for (int n = 0; n < nVertTot; n++)
3384 {
3385 NekDouble ratio;
3386 SpatialDomains::PointGeomSharedPtr vertex = mesh->GetVertex(n);
3387 NekDouble x, y, z;
3388 vertex->GetCoords(x, y, z);
3389 int qp_closer = 0;
3390 // determine the closer xold_up
3391 NekDouble tmp = 1000;
3392
3393 for (int k = 0; k < nvertl; k++)
3394 {
3395 if (abs(x - xcold[k]) < tmp)
3396 {
3397 tmp = abs(x - xcold[k]);
3398 qp_closer = k;
3399 }
3400 }
3401 // find nplay_closer
3402 int nplay_closer;
3403 if (qp_closer == 0)
3404 {
3405 nplay_closer = 0; // first vert
3406 }
3407 else
3408 {
3409 nplay_closer = (qp_closer - 1) * npedge + npedge - 1;
3410 }
3411
3412 if (y > yoldup[qp_closer] && y < 1) // nlays-1 is layerup
3413 {
3414
3415 // ratio =
3416 // (1-layers_y[nlays-1][qp_closer])*(1-y_c[qp_closer])/
3417 // ( (1-yold_up[n])*(1-yold_c[qp_closer]) );
3418 ratio = (1 - ylayup[nplay_closer]) / ((1 - yoldup[qp_closer]));
3419 // distance prop to layerup
3420 ynew[n] = ylayup[nplay_closer] + (y - yoldup[qp_closer]) * ratio;
3421 xnew[n] = x;
3422 }
3423 else if (y < yolddown[qp_closer] && y > -1) // 0 is layerdown
3424 {
3425
3426 ratio = (1 + ylaydown[nplay_closer]) / ((1 + yolddown[qp_closer]));
3427 // distance prop to layerlow
3428 ynew[n] =
3429 ylaydown[nplay_closer] + (y - yolddown[qp_closer]) * ratio;
3430 xnew[n] = x;
3431 }
3432 }
3433}

References tinysimd::abs(), and Nektar::UnitTests::z().

◆ MoveOutsidePointsNnormpos()

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 
)

Definition at line 3435 of file MeshMove.cpp.

3445{
3446 // determine the new verts up/down pos:
3447 int nvertl = xoldup.size();
3448 int nedges = nvertl - 1;
3449 Array<OneD, NekDouble> xnew_down(nvertl);
3450 Array<OneD, NekDouble> ynew_down(nvertl);
3451 Array<OneD, NekDouble> xnew_up(nvertl);
3452 Array<OneD, NekDouble> ynew_up(nvertl);
3453 Array<OneD, NekDouble> nxvert(nvertl);
3454 Array<OneD, NekDouble> nyvert(nvertl);
3455 Array<OneD, NekDouble> norm(nvertl);
3456 Array<OneD, NekDouble> tmp(nvertl);
3457 for (int a = 0; a < nedges; a++)
3458 {
3459 if (a == 0)
3460 {
3461 // v1
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];
3468 // v2
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];
3475 }
3476 else
3477 {
3478 // v2
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];
3485 }
3486 }
3487
3488 NekDouble xmax = Vmath::Vmax(nvertl, xoldup, 1);
3489 NekDouble xmin = Vmath::Vmin(nvertl, xoldup, 1);
3490
3491 // update vertices coords outside layers region
3492 NekDouble ratiox;
3493 // NekDouble ratioy;
3494
3495 int nVertTot = mesh->GetNvertices();
3496 for (int n = 0; n < nVertTot; n++)
3497 {
3498 NekDouble ratio;
3499 SpatialDomains::PointGeomSharedPtr vertex = mesh->GetVertex(n);
3500 NekDouble x, y, z;
3501 vertex->GetCoords(x, y, z);
3502 int qp_closeroldup = 0, qp_closerolddown = 0;
3503 NekDouble diffup, diffdown;
3504 // determine the closer xold_up,down
3505 diffdown = 1000;
3506 diffup = 1000;
3507
3508 for (int k = 0; k < nvertl; k++)
3509 {
3510 if (abs(x - xolddown[k]) < diffdown)
3511 {
3512 diffdown = abs(x - xolddown[k]);
3513 qp_closerolddown = k;
3514 }
3515 if (abs(x - xoldup[k]) < diffup)
3516 {
3517 diffup = abs(x - xoldup[k]);
3518 qp_closeroldup = k;
3519 }
3520 }
3521
3522 // find nplay_closer
3523 diffdown = 1000;
3524 diffup = 1000;
3525
3526 int qp_closerup = 0, qp_closerdown = 0;
3527
3528 for (int f = 0; f < nvertl; f++)
3529 {
3530 if (abs(x - xnew_down[f]) < diffdown)
3531 {
3532 diffdown = abs(x - xnew_down[f]);
3533 qp_closerdown = f;
3534 }
3535 if (abs(x - xnew_up[f]) < diffup)
3536 {
3537 diffup = abs(x - xnew_up[f]);
3538 qp_closerup = f;
3539 }
3540 }
3541
3542 // int qp_closernormoldup;
3543 Vmath::Sadd(nvertl, -x, xoldup, 1, tmp, 1);
3544 Vmath::Vmul(nvertl, tmp, 1, tmp, 1, tmp, 1);
3545 Vmath::Sadd(nvertl, -y, yoldup, 1, norm, 1);
3546 Vmath::Vmul(nvertl, norm, 1, norm, 1, norm, 1);
3547 Vmath::Vadd(nvertl, norm, 1, tmp, 1, norm, 1);
3548 // qp_closernormoldup = Vmath::Imin(nvertl, norm,1);
3549
3550 Vmath::Zero(nvertl, norm, 1);
3551 Vmath::Zero(nvertl, tmp, 1);
3552
3553 // int qp_closernormolddown;
3554 Vmath::Sadd(nvertl, -x, xolddown, 1, tmp, 1);
3555 Vmath::Vmul(nvertl, tmp, 1, tmp, 1, tmp, 1);
3556 Vmath::Sadd(nvertl, -y, yolddown, 1, norm, 1);
3557 Vmath::Vmul(nvertl, norm, 1, norm, 1, norm, 1);
3558 Vmath::Vadd(nvertl, norm, 1, tmp, 1, norm, 1);
3559 // qp_closernormolddown = Vmath::Imin(nvertl, norm,1);
3560
3561 Vmath::Zero(nvertl, norm, 1);
3562 Vmath::Zero(nvertl, tmp, 1);
3563
3564 int qp_closernormup;
3565 Vmath::Sadd(nvertl, -x, xnew_up, 1, tmp, 1);
3566 Vmath::Vmul(nvertl, tmp, 1, tmp, 1, tmp, 1);
3567 Vmath::Sadd(nvertl, -y, ynew_up, 1, norm, 1);
3568 Vmath::Vmul(nvertl, norm, 1, norm, 1, norm, 1);
3569 Vmath::Vadd(nvertl, norm, 1, tmp, 1, norm, 1);
3570 qp_closernormup = Vmath::Imin(nvertl, norm, 1);
3571
3572 Vmath::Zero(nvertl, norm, 1);
3573 Vmath::Zero(nvertl, tmp, 1);
3574
3575 int qp_closernormdown;
3576 Vmath::Sadd(nvertl, -x, xnew_down, 1, tmp, 1);
3577 Vmath::Vmul(nvertl, tmp, 1, tmp, 1, tmp, 1);
3578 Vmath::Sadd(nvertl, -y, ynew_down, 1, norm, 1);
3579 Vmath::Vmul(nvertl, norm, 1, norm, 1, norm, 1);
3580 Vmath::Vadd(nvertl, norm, 1, tmp, 1, norm, 1);
3581 qp_closernormdown = Vmath::Imin(nvertl, norm, 1);
3582
3583 if (y > yoldup[qp_closeroldup] && y < 1)
3584 {
3585
3586 // ratio =
3587 // (1-layers_y[nlays-1][qp_closer])*(1-y_c[qp_closer])/
3588 // ( (1-yold_up[n])*(1-yold_c[qp_closer]) );
3589 ratio = (1 - ynew_up[qp_closerup]) / ((1 - yoldup[qp_closeroldup]));
3590
3591 // ratioy = (1-ynew_up[qp_closernormup])/
3592 // ( (1-yoldup[qp_closernormoldup]) );
3593 // distance prop to layerup
3594 ynew[n] =
3595 ynew_up[qp_closerup] + (y - yoldup[qp_closeroldup]) * ratio;
3596 // ynew[n] = y
3597 // +abs(nyvert[qp_closernormup])*(ynew_up[qp_closeroldup]-yoldup[qp_closeroldup])*ratioy;
3598
3599 // ynew[n] = y + 0.3*(ynew_up[qp_closerup]-yoldup[qp_closerup]);
3600 // xnew[n] = x +
3601 // abs(nxvert[qp_closeroldup])*(xnew_up[qp_closeroldup]-xoldup[qp_closeroldup]);
3602
3603 if (x > (xmax - xmin) / 2. && x < xmax)
3604 {
3605 ratiox = (xmax - xnew_up[qp_closernormup]) /
3606 (xmax - xoldup[qp_closernormup]);
3607 if ((xmax - xoldup[qp_closernormup]) == 0)
3608 {
3609 ratiox = 1.0;
3610 }
3611
3612 // xnew[n] = xnew_up[qp_closerup]
3613 // + (x-xoldup[qp_closerup])*ratiox;
3614 xnew[n] =
3615 x + abs(nxvert[qp_closernormup]) *
3616 (xnew_up[qp_closeroldup] - xoldup[qp_closeroldup]) *
3617 ratiox;
3618 ASSERTL0(x > xmin, " x value <xmin up second half");
3619 ASSERTL0(x < xmax, " x value >xmax up second half");
3620 }
3621 else if (x > xmin && x <= (xmax - xmin) / 2.)
3622 {
3623 // cout<<"up close normold="<<qp_closernormoldup<<"
3624 // closenorm="<<qp_closernormup<<endl;
3625 ratiox = (xnew_up[qp_closernormup] - xmin) /
3626 ((xoldup[qp_closernormup] - xmin));
3627 if ((xoldup[qp_closernormup] - xmin) == 0)
3628 {
3629 ratiox = 1.0;
3630 }
3631 // xnew[n] = xnew_up[qp_closerup]
3632 // + (x-xoldup[qp_closeroldup])*ratiox;
3633 xnew[n] =
3634 x + abs(nxvert[qp_closernormup]) *
3635 (xnew_up[qp_closeroldup] - xoldup[qp_closeroldup]) *
3636 ratiox;
3637 // cout<<"up xold="<<x<<" xnew="<<xnew[n]<<endl;
3638 ASSERTL0(x > xmin, " x value <xmin up first half");
3639 ASSERTL0(x < xmax, " x value >xmax up first half");
3640 }
3641 }
3642 else if (y < yolddown[qp_closerolddown] && y > -1)
3643 {
3644
3645 ratio = (1 + ynew_down[qp_closerdown]) /
3646 ((1 + yolddown[qp_closerolddown]));
3647
3648 // ratioy = (1-ynew_down[qp_closernormdown])/
3649 // ( (1-yolddown[qp_closernormolddown]) );
3650
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;
3656 // ynew[n] = y +
3657 // 0.3*(ynew_down[qp_closerdown]-yolddown[qp_closerdown]); xnew[n] =
3658 // x +
3659 // abs(nxvert[qp_closerolddown])*(xnew_down[qp_closerolddown]-xolddown[qp_closerolddown]);
3660 /*
3661 if(n==74)
3662 {
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;
3667 }
3668 */
3669
3670 if (x > (xmax - xmin) / 2. && x < xmax)
3671 {
3672 ratiox = (xmax - xnew_down[qp_closernormdown]) /
3673 ((xmax - xolddown[qp_closernormdown]));
3674 if ((xmax - xolddown[qp_closernormdown]) == 0)
3675 {
3676 ratiox = 1.0;
3677 }
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]) *
3683 ratiox;
3684 ASSERTL0(x > xmin, " x value <xmin down second half");
3685 ASSERTL0(x < xmax, " x value >xmax down second half");
3686 }
3687 else if (x > xmin && x <= (xmax - xmin) / 2.)
3688 {
3689 ratiox = (xnew_down[qp_closernormdown] - xmin) /
3690 ((xolddown[qp_closernormdown] - xmin));
3691 if ((xolddown[qp_closernormdown] - xmin) == 0)
3692 {
3693 ratiox = 1.0;
3694 }
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]) *
3700 ratiox;
3701 ASSERTL0(x > xmin, " x value <xmin down first half");
3702 ASSERTL0(x < xmax, " x value >xmax down first half");
3703 }
3704 }
3705
3706 cout << "xold" << x << " xnew=" << xnew[n] << endl;
3707 ASSERTL0(xnew[n] >= xmin, "newx < xmin");
3708 ASSERTL0(xnew[n] <= xmax, "newx > xmax");
3709
3710 } // verts closed
3711}

References tinysimd::abs(), ASSERTL0, Vmath::Imin(), Vmath::Sadd(), Vmath::Vadd(), Vmath::Vmax(), Vmath::Vmin(), Vmath::Vmul(), Nektar::UnitTests::z(), and Vmath::Zero().

Referenced by main().

◆ Orderfunctionx()

void Orderfunctionx ( Array< OneD, NekDouble >  inarray_x,
Array< OneD, NekDouble >  inarray_y,
Array< OneD, NekDouble > &  outarray_x,
Array< OneD, NekDouble > &  outarray_y 
)

Definition at line 3138 of file MeshMove.cpp.

3142{
3143 Array<OneD, NekDouble> tmpx(inarray_x.size());
3144 Array<OneD, NekDouble> tmpy(inarray_x.size());
3145 // local copy to prevent overwriting
3146 Vmath::Vcopy(inarray_x.size(), inarray_x, 1, tmpx, 1);
3147 Vmath::Vcopy(inarray_x.size(), inarray_y, 1, tmpy, 1);
3148
3149 // order function with respect to x
3150 int index;
3151 NekDouble max = Vmath::Vmax(tmpx.size(), tmpx, 1);
3152 for (int w = 0; w < tmpx.size(); w++)
3153 {
3154 index = Vmath::Imin(tmpx.size(), tmpx, 1);
3155 outarray_x[w] = tmpx[index];
3156 outarray_y[w] = tmpy[index];
3157 if (w < tmpx.size() - 1) // case of repetitions
3158 {
3159 if (tmpx[index] == tmpx[index + 1])
3160 {
3161 outarray_x[w + 1] = tmpx[index + 1];
3162 outarray_y[w + 1] = tmpy[index + 1];
3163 tmpx[index + 1] = max + 1000;
3164 w++;
3165 }
3166 }
3167 /*
3168 if(w>0)//case of repetitions
3169 {
3170 if(inarray_x[index] == tmpx[index-1])
3171 {
3172 outarray_x[w+1]= tmpx[index-1];
3173 outarray_y[w+1]= tmpy[index-1];
3174 tmpx[index-1] = max+1000;
3175 w++;
3176 }
3177 }
3178 */
3179 tmpx[index] = max + 1000;
3180 }
3181}

References Vmath::Imin(), Vmath::Vcopy(), Vmath::Vmax(), and Nektar::UnitTests::w().

Referenced by main(), MoveLayerNfixedxpos(), and MoveLayerNnormpos().

◆ OrderVertices()

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 
)

Definition at line 2007 of file MeshMove.cpp.

2012{
2013 int nvertl = nedges + 1;
2014 int edge;
2015 Array<OneD, int> Vids_temp(nvertl, -10);
2016 NekDouble x0layer = 0.0;
2017 for (int j = 0; j < nedges; j++)
2018 {
2019 LocalRegions::SegExpSharedPtr bndSegExplow =
2020 bndfield->GetExp(j)->as<LocalRegions::SegExp>();
2021 edge = (bndSegExplow->GetGeom1D())->GetGlobalID();
2022 // cout<<" edge="<<edge<<endl;
2023 for (int k = 0; k < 2; k++)
2024 {
2025 Vids_temp[j + k] = (bndSegExplow->GetGeom1D())->GetVid(k);
2026 SpatialDomains::PointGeomSharedPtr vertex =
2027 graphShPt->GetVertex(Vids_temp[j + k]);
2028 NekDouble x1, y1, z1;
2029 vertex->GetCoords(x1, y1, z1);
2030 if (x1 == x_connect && edge != lastedge)
2031 {
2032 // first 2 points
2033 if (x_connect == x0layer)
2034 {
2035 Vids[v1] = Vids_temp[j + k];
2036 x[v1] = x1;
2037 y[v1] = y1;
2038
2039 if (k == 0)
2040 {
2041 Vids_temp[j + 1] =
2042 (bndSegExplow->GetGeom1D())->GetVid(1);
2043 Vids[v2] = Vids_temp[j + 1];
2044 SpatialDomains::PointGeomSharedPtr vertex =
2045 graphShPt->GetVertex(Vids[v2]);
2046 NekDouble x2, y2, z2;
2047 vertex->GetCoords(x2, y2, z2);
2048 x[v2] = x2;
2049 y[v2] = y2;
2050 }
2051 if (k == 1)
2052 {
2053 Vids_temp[j + 0] =
2054 (bndSegExplow->GetGeom1D())->GetVid(0);
2055 Vids[v2] = Vids_temp[j + 0];
2056 SpatialDomains::PointGeomSharedPtr vertex =
2057 graphShPt->GetVertex(Vids[v2]);
2058 NekDouble x2, y2, z2;
2059 vertex->GetCoords(x2, y2, z2);
2060 x[v2] = x2;
2061 y[v2] = y2;
2062 }
2063 }
2064 else // other vertices
2065 {
2066 if (k == 0)
2067 {
2068 Vids_temp[j + 1] =
2069 (bndSegExplow->GetGeom1D())->GetVid(1);
2070 Vids[v1] = Vids_temp[j + 1];
2071 SpatialDomains::PointGeomSharedPtr vertex =
2072 graphShPt->GetVertex(Vids[v1]);
2073 NekDouble x1, y1, z1;
2074 vertex->GetCoords(x1, y1, z1);
2075 x[v1] = x1;
2076 y[v1] = y1;
2077 }
2078 if (k == 1)
2079 {
2080 Vids_temp[j + 0] =
2081 (bndSegExplow->GetGeom1D())->GetVid(0);
2082 Vids[v1] = Vids_temp[j + 0];
2083 SpatialDomains::PointGeomSharedPtr vertex =
2084 graphShPt->GetVertex(Vids[v1]);
2085 NekDouble x1, y1, z1;
2086 vertex->GetCoords(x1, y1, z1);
2087 x[v1] = x1;
2088 y[v1] = y1;
2089 }
2090 }
2091 break;
2092 }
2093 }
2094 if (Vids[v1] != -10)
2095 {
2096 lastedge = edge;
2097 break;
2098 }
2099 }
2100}

Referenced by main().

◆ PolyFit()

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 
)

Definition at line 3029 of file MeshMove.cpp.

3033{
3034
3035 int N = polyorder + 1;
3036 Array<OneD, NekDouble> A(N * N, 0.0);
3037 Array<OneD, NekDouble> b(N, 0.0);
3038 cout << npoints << endl;
3039 for (int u = 0; u < npoints; u++)
3040 {
3041 cout << "c=" << xin[u] << " " << fin[u] << endl;
3042 }
3043 // fill column by column
3044 // e counts cols
3045 for (int e = 0; e < N; e++)
3046 {
3047
3048 for (int row = 0; row < N; row++)
3049 {
3050 for (int w = 0; w < npoints; w++)
3051 {
3052 A[N * e + row] += std::pow(xin[w], e + row);
3053 }
3054 }
3055 }
3056
3057 for (int row = 0; row < N; row++)
3058 {
3059 for (int h = 0; h < npoints; h++)
3060 {
3061 b[row] += fin[h] * std::pow(xin[h], row);
3062 // cout<<b="<<b[row]<<" y_c="<<y_c[r]<<endl;
3063 }
3064 }
3065
3066 // cout<<"A elements="<<A.size()<<endl;
3067 Array<OneD, int> ipivot(N);
3068 int info = 0;
3069 // Lapack::Dgesv( N, 1, A.get(), N, ipivot.get(), b.get(), N, info);
3070 Lapack::Dgetrf(N, N, A.get(), N, ipivot.get(), info);
3071
3072 if (info < 0)
3073 {
3074 std::string message =
3075 "ERROR: The " + std::to_string(-info) +
3076 "th parameter had an illegal parameter for dgetrf";
3077 ASSERTL0(false, message.c_str());
3078 }
3079 else if (info > 0)
3080 {
3081 std::string message = "ERROR: Element u_" + std::to_string(info) +
3082 std::to_string(info) + " is 0 from dgetrf";
3083 ASSERTL0(false, message.c_str());
3084 }
3085 // N means no transponse (direct matrix)
3086 int ncolumns_b = 1;
3087 Lapack::Dgetrs('N', N, ncolumns_b, A.get(), N, ipivot.get(), b.get(), N,
3088 info);
3089 if (info < 0)
3090 {
3091 std::string message =
3092 "ERROR: The " + std::to_string(-info) +
3093 "th parameter had an illegal parameter for dgetrf";
3094 ASSERTL0(false, message.c_str());
3095 }
3096 else if (info > 0)
3097 {
3098 std::string message = "ERROR: Element u_" + std::to_string(info) +
3099 std::to_string(info) + " is 0 from dgetrf";
3100 ASSERTL0(false, message.c_str());
3101 }
3102
3103 // the lower coeff the lower is the grad
3104 // reverse:
3105 Array<OneD, NekDouble> tmp(N);
3106 Vmath::Vcopy(N, b, 1, tmp, 1);
3107 int cnt = N;
3108 for (int j = 0; j < N; j++)
3109 {
3110 b[j] = tmp[cnt - 1];
3111 cnt--;
3112 }
3113
3114 for (int h = 0; h < N; h++)
3115 {
3116 cout << "coeff:" << b[h] << endl;
3117 }
3118
3119 // ovewrite the ycPhysMOD:
3120 int polorder;
3121 for (int c = 0; c < npout; c++)
3122 {
3123 polorder = polyorder;
3124 // put ycPhysMOD to zero
3125 fout[c] = 0;
3126 for (int d = 0; d < N; d++)
3127 {
3128
3129 fout[c] += b[d] * std::pow(xout[c], polorder);
3130 polorder--;
3131 }
3132 }
3133
3134 // write coeffs
3135 Vmath::Vcopy(N, b, 1, coeffsinterp, 1);
3136}
Lapack::Dgetrf
static void Dgetrf(const int &m, const int &n, double *a, const int &lda, int *ipiv, int &info)
General matrix LU factorisation.
Definition: Lapack.hpp:262
Lapack::Dgetrs
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.
Definition: Lapack.hpp:269

References ASSERTL0, Nektar::UnitTests::d(), Lapack::Dgetrf(), Lapack::Dgetrs(), Vmath::Vcopy(), and Nektar::UnitTests::w().

Referenced by main().

◆ PolyInterp()

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 
)

Definition at line 2939 of file MeshMove.cpp.

2943{
2944 int np_pol = xpol.size();
2945 int N = np_pol;
2946 Array<OneD, NekDouble> A(N * N, 1.0);
2947 Array<OneD, NekDouble> b(N);
2948 int row = 0;
2949 // fill column by column
2950 for (int e = 0; e < N; e++)
2951 {
2952 row = 0;
2953 for (int w = 0; w < N; w++)
2954 {
2955 A[N * e + row] = std::pow(xpol[w], N - 1 - e);
2956 row++;
2957 }
2958 }
2959 row = 0;
2960 for (int r = 0; r < np_pol; r++)
2961 {
2962 b[row] = ypol[r];
2963 // cout<<"b="<<b[row]<<" y_c="<<y_c[r]<<endl;
2964 row++;
2965 }
2966
2967 // cout<<"A elements="<<A.size()<<endl;
2968 Array<OneD, int> ipivot(N);
2969 int info = 0;
2970 // Lapack::Dgesv( N, 1, A.get(), N, ipivot.get(), b.get(), N, info);
2971 Lapack::Dgetrf(N, N, A.get(), N, ipivot.get(), info);
2972 if (info < 0)
2973 {
2974 std::string message =
2975 "ERROR: The " + std::to_string(-info) +
2976 "th parameter had an illegal parameter for dgetrf";
2977 ASSERTL0(false, message.c_str());
2978 }
2979 else if (info > 0)
2980 {
2981 std::string message = "ERROR: Element u_" + std::to_string(info) +
2982 std::to_string(info) + " is 0 from dgetrf";
2983 ASSERTL0(false, message.c_str());
2984 }
2985
2986 // N means no transponse (direct matrix)
2987 int ncolumns_b = 1;
2988 Lapack::Dgetrs('N', N, ncolumns_b, A.get(), N, ipivot.get(), b.get(), N,
2989 info);
2990 if (info < 0)
2991 {
2992 std::string message =
2993 "ERROR: The " + std::to_string(-info) +
2994 "th parameter had an illegal parameter for dgetrf";
2995 ASSERTL0(false, message.c_str());
2996 }
2997 else if (info > 0)
2998 {
2999 std::string message = "ERROR: Element u_" + std::to_string(info) +
3000 std::to_string(info) + " is 0 from dgetrf";
3001 ASSERTL0(false, message.c_str());
3002 }
3003 /*
3004 for(int h=0; h<np_pol; h++)
3005 {
3006 cout<<"coeff:"<<b[h]<<endl;
3007 }
3008 */
3009
3010 // ovewrite the ycPhysMOD:
3011 int polorder;
3012 for (int c = 0; c < npedge; c++)
3013 {
3014 polorder = np_pol - 1;
3015 // put ycPhysMOD to zero
3016 ycout[edge * (npedge) + c + 1] = 0;
3017 for (int d = 0; d < np_pol; d++)
3018 {
3019 ycout[edge * (npedge) + c + 1] +=
3020 b[d] * std::pow(xcout[edge * (npedge) + c + 1], polorder);
3021 polorder--;
3022 }
3023 }
3024
3025 // write coeffs
3026 Vmath::Vcopy(np_pol, b, 1, coeffsinterp, 1);
3027}

References ASSERTL0, Nektar::UnitTests::d(), Lapack::Dgetrf(), Lapack::Dgetrs(), Vmath::Vcopy(), and Nektar::UnitTests::w().

◆ Replacevertices()

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 
)

Pull out lhs and rhs and eliminate any spaces.

Definition at line 3893 of file MeshMove.cpp.

3900{
3901 // load existing file
3902 string newfile;
3903 TiXmlDocument doc(filename);
3904 // load xscale parameter (if exists)
3905 TiXmlElement *master = doc.FirstChildElement("NEKTAR");
3906 TiXmlElement *mesh = master->FirstChildElement("GEOMETRY");
3907 TiXmlElement *element = mesh->FirstChildElement("VERTEX");
3908 NekDouble xscale = 1.0;
3909 LibUtilities::Interpreter expEvaluator;
3910 const char *xscal = element->Attribute("XSCALE");
3911 if (xscal)
3912 {
3913 std::string xscalstr = xscal;
3914 int expr_id = expEvaluator.DefineFunction("", xscalstr);
3915 xscale = expEvaluator.Evaluate(expr_id);
3916 }
3917
3918 // Save a new XML file.
3919 newfile = filename.substr(0, filename.find_last_of(".")) + "_moved.xml";
3920 doc.SaveFile(newfile);
3921
3922 // write the new vertices
3923 TiXmlDocument docnew(newfile);
3924 bool loadOkaynew = docnew.LoadFile();
3925
3926 std::string errstr = "Unable to load file: ";
3927 errstr += newfile;
3928 ASSERTL0(loadOkaynew, errstr.c_str());
3929
3930 TiXmlHandle docHandlenew(&docnew);
3931 TiXmlElement *meshnew = nullptr;
3932 TiXmlElement *masternew = nullptr;
3933 TiXmlElement *condnew = nullptr;
3934 TiXmlElement *Parsnew = nullptr;
3935 TiXmlElement *parnew = nullptr;
3936
3937 // Master tag within which all data is contained.
3938
3939 masternew = docnew.FirstChildElement("NEKTAR");
3940 ASSERTL0(masternew, "Unable to find NEKTAR tag in file.");
3941
3942 // set the alpha value
3943 string alphastring;
3944 condnew = masternew->FirstChildElement("CONDITIONS");
3945 Parsnew = condnew->FirstChildElement("PARAMETERS");
3946 cout << "alpha=" << s_alp << endl;
3947 parnew = Parsnew->FirstChildElement("P");
3948 while (parnew)
3949 {
3950 TiXmlNode *node = parnew->FirstChild();
3951 if (node)
3952 {
3953 // Format is "paramName = value"
3954 std::string line = node->ToText()->Value();
3955 std::string lhs;
3956 std::string rhs;
3957 /// Pull out lhs and rhs and eliminate any spaces.
3958 int beg = line.find_first_not_of(" ");
3959 int end = line.find_first_of("=");
3960 // Check for no parameter name
3961 if (beg == end)
3962 {
3963 throw 1;
3964 }
3965 // Check for no parameter value
3966 if (end != line.find_last_of("="))
3967 {
3968 throw 1;
3969 }
3970 // Check for no equals sign
3971 if (end == std::string::npos)
3972 {
3973 throw 1;
3974 }
3975 lhs = line.substr(line.find_first_not_of(" "), end - beg);
3976 lhs = lhs.substr(0, lhs.find_last_not_of(" ") + 1);
3977
3978 // rhs = line.substr(line.find_last_of("=")+1);
3979 // rhs = rhs.substr(rhs.find_first_not_of(" "));
3980 // rhs = rhs.substr(0, rhs.find_last_not_of(" ")+1);
3981
3982 boost::to_upper(lhs);
3983 if (lhs == "ALPHA")
3984 {
3985 alphastring = "Alpha = " + s_alp;
3986 parnew->RemoveChild(node);
3987 parnew->LinkEndChild(new TiXmlText(alphastring));
3988 }
3989 }
3990
3991 parnew = parnew->NextSiblingElement("P");
3992 }
3993
3994 // Find the Mesh tag and same the dim and space attributes
3995 meshnew = masternew->FirstChildElement("GEOMETRY");
3996
3997 ASSERTL0(meshnew, "Unable to find GEOMETRY tag in file.");
3998 // Now read the vertices
3999 TiXmlElement *elementnew = meshnew->FirstChildElement("VERTEX");
4000 ASSERTL0(elementnew, "Unable to find mesh VERTEX tag in file.");
4001 // set xscale 1!!
4002 if (xscale != 1.0)
4003 {
4004 elementnew->SetAttribute("XSCALE", 1.0);
4005 }
4006 TiXmlElement *vertexnew = elementnew->FirstChildElement("V");
4007
4008 int indx;
4009 int err, numPts;
4010 int nextVertexNumber = -1;
4011
4012 while (vertexnew)
4013 {
4014 nextVertexNumber++;
4015 // delete the old one
4016 TiXmlAttribute *vertexAttr = vertexnew->FirstAttribute();
4017 std::string attrName(vertexAttr->Name());
4018 ASSERTL0(attrName == "ID",
4019 (std::string("Unknown attribute name: ") + attrName).c_str());
4020
4021 err = vertexAttr->QueryIntValue(&indx);
4022 ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
4023 ASSERTL0(indx == nextVertexNumber,
4024 "Vertex IDs must begin with zero and be sequential.");
4025
4026 std::string vertexBodyStr;
4027 // Now read body of vertex
4028 TiXmlNode *vertexBody = vertexnew->FirstChild();
4029 // Accumulate all non-comment body data.
4030 if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
4031 {
4032 vertexBodyStr += vertexBody->ToText()->Value();
4033 vertexBodyStr += " ";
4034 }
4035 ASSERTL0(!vertexBodyStr.empty(),
4036 "Vertex definitions must contain vertex data.");
4037 // remove the old coordinates
4038 vertexnew->RemoveChild(vertexBody);
4039 // write the new one
4040 // cout<<"writing.. v:"<<nextVertexNumber<<endl;
4041 stringstream s;
4042 // we need at least 5 digits (setprecision 5) to get the streak position
4043 // with
4044 // precision 10^-10
4045 s << std::scientific << std::setprecision(8) << newx[nextVertexNumber]
4046 << " " << newy[nextVertexNumber] << " " << 0.0;
4047 vertexnew->LinkEndChild(new TiXmlText(s.str()));
4048 // TiXmlNode *newvertexBody = vertexnew->FirstChild();
4049 // string newvertexbodystr= newvertexBody->SetValue(s.str());
4050 // vertexnew->ReplaceChild(vertexBody,new TiXmlText(newvertexbodystr));
4051
4052 vertexnew = vertexnew->NextSiblingElement("V");
4053 }
4054
4055 // read the curved tag
4056 TiXmlElement *curvednew = meshnew->FirstChildElement("CURVED");
4057 ASSERTL0(curvednew, "Unable to find mesh CURVED tag in file.");
4058 TiXmlElement *edgenew = curvednew->FirstChildElement("E");
4059 int cnt = -1;
4060 // ID is different from index...
4061 std::string charindex;
4062 int eid;
4063 int index;
4064 int indexeid;
4065 int neids_lay = lay_eids[0].size();
4066 // if edgenew belongs to the crit lay replace it, else delete it.
4067 while (edgenew)
4068 {
4069 indexeid = -1;
4070 cnt++;
4071 // get the index...
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;
4077 // get the eid
4078 edgenew->QueryIntAttribute("EDGEID", &eid);
4079 // cout<<"eid="<<eid<<" neid="<<Eids.size()<<endl;
4080 // find the corresponding index curve point
4081 for (int u = 0; u < Eids.size(); u++)
4082 {
4083 // cout<<"Eids="<<Eids[u]<<" eid="<<eid<<endl;
4084 if (Eids[u] == eid)
4085 {
4086 indexeid = u;
4087 }
4088 }
4089 if (indexeid == -1)
4090 {
4091 curvednew->RemoveChild(edgenew);
4092 // ASSERTL0(false, "edge to update not found");
4093 }
4094 else
4095 {
4096
4097 std::string edgeBodyStr;
4098 // read the body of the edge
4099 TiXmlNode *edgeBody = edgenew->FirstChild();
4100 if (edgeBody->Type() == TiXmlNode::TINYXML_TEXT)
4101 {
4102 edgeBodyStr += edgeBody->ToText()->Value();
4103 edgeBodyStr += " ";
4104 }
4105 ASSERTL0(!edgeBodyStr.empty(),
4106 "Edge definitions must contain edge data");
4107 // remove the old coordinates
4108 edgenew->RemoveChild(edgeBody);
4109 // write the new points coordinates
4110 // we need at least 5 digits (setprecision 5) to get the streak
4111 // position with
4112 // precision 10^-10
4113
4114 // Determine the number of points
4115 err = edgenew->QueryIntAttribute("NUMPOINTS", &numPts);
4116 ASSERTL0(err == TIXML_SUCCESS,
4117 "Unable to read curve attribute NUMPOINTS.");
4118
4119 stringstream st;
4120 edgenew->SetAttribute("NUMPOINTS", Npoints);
4121 for (int u = 0; u < Npoints; u++)
4122 {
4123 st << std::scientific << std::setprecision(8)
4124 << xcPhys[cnt * Npoints + u] << " "
4125 << ycPhys[cnt * Npoints + u] << " " << 0.000 << " ";
4126 }
4127
4128 edgenew->LinkEndChild(new TiXmlText(st.str()));
4129
4130 /*
4131 st << std::scientific <<
4132 std::setprecision(8) << x_crit[v1] << " "
4133 << y_crit[v1] << " " << 0.000<<" ";
4134 for(int a=0; a< Npoints-2; a++)
4135 {
4136 st << std::scientific <<
4137 std::setprecision(8) << " "<<Pcurvx[indexeid*(Npoints-2)
4138 +a]<<" "<<Pcurvy[indexeid*(Npoints-2) +a]
4139 <<" "<<0.000<<" ";
4140
4141 }
4142 st << std::scientific <<
4143 std::setprecision(8) << " "<<x_crit[v2]<<" "<< y_crit[v2]
4144 <<" "<< 0.000; edgenew->LinkEndChild(new TiXmlText(st.str()));
4145
4146 */
4147 }
4148
4149 edgenew = edgenew->NextSiblingElement("E");
4150 }
4151
4152 // write also the others layers curve points
4153 if (curv_lay == true)
4154 {
4155 cout << "write other curved edges" << endl;
4156 TiXmlElement *curved = meshnew->FirstChildElement("CURVED");
4157 int idcnt = 300;
4158 int nlays = lay_eids.size();
4159
4160 // TiXmlComment * comment = new TiXmlComment();
4161 // comment->SetValue(" new edges ");
4162 // curved->LinkEndChild(comment);
4163 for (int g = 0; g < nlays; ++g)
4164 {
4165 for (int p = 0; p < neids_lay; p++)
4166 {
4167 stringstream st;
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++)
4174 {
4175 st << std::scientific << std::setprecision(8)
4176 << x_lay[g][p * Npoints + c] << " "
4177 << y_lay[g][p * Npoints + c] << " " << 0.000 << " ";
4178 }
4179
4180 TiXmlText *t0 = new TiXmlText(st.str());
4181 e->LinkEndChild(t0);
4182 curved->LinkEndChild(e);
4183 }
4184 }
4185 }
4186
4187 docnew.SaveFile(newfile);
4188
4189 cout << "new file: " << newfile << endl;
4190}
Nektar::LibUtilities::Interpreter
Interpreter class for the evaluation of mathematical expressions.
Definition: Interpreter.h:76
Nektar::LibUtilities::Interpreter::DefineFunction
int DefineFunction(const std::string &vlist, const std::string &expr)
Defines a function for the purposes of evaluation.
Definition: Interpreter/Interpreter.cpp:2242
Nektar::LibUtilities::Interpreter::Evaluate
NekDouble Evaluate(const int id)
Evaluate a function which depends only on constants and/or parameters.
Definition: Interpreter/Interpreter.cpp:2248
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350

References ASSERTL0, Nektar::LibUtilities::Interpreter::DefineFunction(), Nektar::LibUtilities::Interpreter::Evaluate(), and CellMLToNektar.cellml_metadata::p.

Referenced by main().

Generated on Wed Jan 31 2024 09:35:59 for Nektar++ by doxygen 1.9.4