/*==================================================================================== ! Purpose(s): ! ! Contain(s): ! ! Author(s): Javier Lopez Lara (jav.lopez@unican.es) ! Gabriel Barajas (barajasg@unican.es) !====================================================================================*/ #include "CArrayMasks.h" #include "CVars.h" #include "C2DArray.h" #include "CStress.h" #include "CMisc.h" #include "CCalcFricVel.h" #include "CErrors.h" #include "math.h" #include #include #include "CMathMisc.h" #ifdef __INTEL_COMPILER #define PYMM ckepsilon #else #define PYMM ckepsilon_ #endif extern "C" void _FORTRAN PYMM (_DOUBLE_ARRAY AAC, _DOUBLE_ARRAY AAR, _DOUBLE_ARRAY AAT, _DOUBLE_ARRAY AU, _DOUBLE_ARRAY AV, _DOUBLE_ARRAY AUn, _DOUBLE_ARRAY AVn, _DOUBLE_ARRAY AXNut, _DOUBLE_ARRAY AXNutn, _DOUBLE_ARRAY ATauxx, _DOUBLE_ARRAY ATauyy, _DOUBLE_ARRAY ATauxy, _DOUBLE_ARRAY ATauxxTurb, _DOUBLE_ARRAY ATauyyTurb, _DOUBLE_ARRAY ATauxyTurb, _DOUBLE_ARRAY APorosity, _DOUBLE_ARRAY AK, _DOUBLE_ARRAY AEp, _DOUBLE_ARRAY AKn, _DOUBLE_ARRAY AEpn, _DOUBLE_ARRAY AF, _DOUBLE_ARRAY AFn, _DOUBLE_ARRAY AUtmp, _DOUBLE_ARRAY AVtmp, _DOUBLE_ARRAY AXCoefs, _DOUBLE_ARRAY AYCoefs, _DOUBLE_ARRAY AXi, _DOUBLE_ARRAY AYj, _DOUBLE_ARRAY ADelX, _DOUBLE_ARRAY ADelY, _DOUBLE AMu, _DOUBLE ANu, _DOUBLE ACMu, _DOUBLE AEddyCoef, _DOUBLE ARhoF, _DOUBLE ADelT, _INT AKb, _INT AKt, _INT AKl, _INT AKr, _INT AIMax, _INT AJMax, _INT_ARRAY ANf, _INT_ARRAY ANpc, _INT Aim1, _INT Ajm1, _BOOL_ARRAY AGage, _DOUBLE_ARRAY AD50, _INT APorosityModel, _INT AWavemakerType, _INT AKEModel, _DOUBLE ASigE, _DOUBLE ASigK, _DOUBLE ARatio, _DOUBLE AUTurb, _DOUBLE Asd, _INT Anks, _DOUBLE AC1e, _DOUBLE AC2e, _DOUBLE AC1ea, _DOUBLE AC1eb, _DOUBLE t, _INT nrough, _DOUBLE Acrest, _INT_ARRAY ALowTurbRect, _INT Anloc ) { /*!======================================================================= ! Purpose(s): ! ! Author(s): Javier Lopez Lara (jav.lopez@unican.es) ! Gabriel Barajas (barajasg@unican.es) !=======================================================================*/ // Constants. const double Gamma = 1.0L; double xmaxxnut = ANu*1.0e-6L; // Roughness. double nks=Anks*1.0; double Ks=nks*Asd; double prodmax = 0.0; Acrest = 1.0e16L; int ValNod, ValNod2, ValNod4, ValNod5; double FricV, oldEp, oldK, Prod, TotProd, KMax, EpMax, Eps2KPor, EpsPor, K2EpsPor; double c1ee, c2ee, xmaxxk, xmaxxep, KPor, ttots; double CMu; int ni, nj; // Type of boundary conditions. bool KbRigid = (AKb==bcRigidFreeSlip || AKb==bcRigidNoSlip); bool KtRigid = (AKt==bcRigidFreeSlip || AKt==bcRigidNoSlip); bool KlRigid = (AKl==bcRigidFreeSlip || AKl==bcRigidNoSlip); bool KrRigid = (AKr==bcRigidFreeSlip || AKr==bcRigidNoSlip); bool KbContOutFlow = (AKb==bcContinuativeOutflow); bool KtContOutFlow = (AKt==bcContinuativeOutflow); bool KlContOutFlow = (AKl==bcContinuativeOutflow); bool KrContOutFlow = (AKr==bcContinuativeOutflow); // These arrays are used to store turbulence values at probe locations. They are allocated at // the begining of the program. static T2DArray disp_t(Anloc, AJMax), fkx_t(Anloc, AJMax), fky_t(Anloc, AJMax), prod_t(Anloc, AJMax), visx_t(Anloc, AJMax), visy_t(Anloc, AJMax); // Finds the rectangle where low turbulence condition is enforced. TIntRect LowTurbRect(ALowTurbRect); // Does it have to write AGages at this time step? bool HasToWriteGages=(Anloc>0); // Is body surface rough? bool IsRoughSurface=(nrough==1); { //Just a block for stablishing the scope of the following variables. // Creates the mesh object. TMesh Mesh(AIMax, AJMax, AXCoefs, AYCoefs, AXi, AYj, ADelX, ADelY); // Creates the masks. TCMask wF (Mesh, AF); TCMask wFn (Mesh, AFn); TCMask wAC (Mesh, AAC); TRMask wAR (Mesh, AAR); TTMask wAT (Mesh, AAT); TRMask wU (Mesh, AU); TTMask wV (Mesh, AV); TRMask wUtmp (Mesh, AUtmp); TTMask wVtmp (Mesh, AVtmp); TCMask wK (Mesh, AK); TCMask wEp (Mesh, AEp); TCMask wKn (Mesh, AKn); TCMask wEpn (Mesh, AEpn); TCMask wPorosity (Mesh, APorosity); TCMask wNuT (Mesh, AXNut); TCMask wNuTn (Mesh, AXNutn); TIntMask wNf (Mesh, ANf); TIntMask wNpc (Mesh, ANpc); TCMask wTauxxTurb (Mesh, ATauxxTurb); TTRMask wTauxyTurb (Mesh, ATauxyTurb); TCMask wTauyyTurb (Mesh, ATauyyTurb); TRMask wUAR (Mesh); TTMask wVAT (Mesh); TBLMask wUyVx (Mesh); TCMask wUx (Mesh); TCMask wVy (Mesh); TCMask wEpnPor (Mesh); TLMask wEpnPorDx (Mesh); TBMask wEpnPorDy (Mesh); TLMask wEpnPorDNutSigENux(Mesh); TBMask wEpnPorDNutSigENuy(Mesh); TCMask wKPor (Mesh); TLMask wKPorDx (Mesh); TBMask wKPorDy (Mesh); TLMask wKPorDNutSigKNux (Mesh); TBMask wKPorDNutSigKNuy (Mesh); TBLMask wTauxyTurbUyVx (Mesh); TLMask wKnx (Mesh); TBMask wKny (Mesh); TLMask wEpnx (Mesh); TBMask wEpny (Mesh); Mesh.DefineArea(2, AIMax-1, 2, AJMax-1); Mesh.First2(); do { int i=Mesh.I; int j=Mesh.J; bool IsLSide=(i==2); bool IsRSide=(i==AIMax-1); bool IsBSide=(j==2); bool IsTSide=(j==AJMax-1); Prod=0.0; EpsPor=0.0; KPor=0.0; Eps2KPor=0.0; K2EpsPor=0.0; // Calculates which cells have porous media. int PorousCells; if (APorosityModel!=pmNone) PorousCells=wNpc.NodesNotEqualTo(1); else PorousCells=0; bool IsPorousCell=((PorousCells & C_B)!=0); // Calculate which cells are on the surface. int SurfaceCells=wNf.NodesNotEqualTo(ctCellFullOfFluid); bool IsSurface=((SurfaceCells & LTRB_B)!=0); // Calculte which cells are empty. int EmptyCells=wNf.NodesEqualTo(ctEmptyCell); // Finds the nodes that will be considered in various derivatives. ValNod=wFn.NodesGreaterThan(1e-3) & wAC.NodesGreaterThan(1e-3); // Cells completely full of fluid and with no porous media. ValNod2=~SurfaceCells & ~PorousCells; ValNod4=ValidNodes(!(KlContOutFlow && IsLSide), !(KrContOutFlow && IsRSide) , !(KtContOutFlow && IsTSide), !(KbContOutFlow && IsBSide)); ValNod4=ValNod4 & ~EmptyCells; ValNod5=BLValidNodesFromC(ValNod2); // This code has to be executed at each loop iteration. wUAR=wU*wAR; wVAT=wV*wAT; TVector VelC(wUAR.IC(), wVAT.IC()); TVector SqrVelC(Sqr(VelC.x), Sqr(VelC.y)); wUyVx=BLM_Dy(wUAR)+BLM_Dx(wVAT); wUyVx=TIntBLMask(wUyVx, 0.0, ValNod5); wUx=CM_Dx(wUAR); wVy=CM_Dy(wVAT); // Additional effort to make large scale production zero in porous media wUx=TIntCMask(wUx, 0.0, ~PorousCells); wVy=TIntCMask(wVy, 0.0, ~PorousCells); wEpnPor=wEpn/wPorosity; wEpnx=LM_Dx(wEpn); wEpny=BM_Dy(wEpn); wEpnPorDx=LM_Dx(wEpnPor, ValNod); wEpnPorDy=BM_Dy(wEpnPor, ValNod); double EpnPorDx=C_Dx(wEpnPor, Gamma, VelC, ValNod); double EpnPorDy=C_Dy(wEpnPor, Gamma, VelC, ValNod); wEpnPorDNutSigENux=wEpnPorDx*(wNuTn.MIL(ValNod4)/ASigE+ANu); wEpnPorDNutSigENuy=wEpnPorDy*(wNuTn.MIB(ValNod4)/ASigE+ANu); wKPor=wK/wPorosity; wKnx=LM_Dx(wKn); wKny=BM_Dy(wKn); wKPorDx=LM_Dx(wKPor, ValNod); wKPorDy=BM_Dy(wKPor, ValNod); double KPorDx=C_Dx(wKPor, Gamma, VelC, ValNod); double KPorDy=C_Dy(wKPor, Gamma, VelC, ValNod); wKPorDNutSigKNux=wKPorDx*(wNuTn.MIL(ValNod4)/ASigK+ANu); wKPorDNutSigKNuy=wKPorDy*(wNuTn.MIB(ValNod4)/ASigK+ANu); wTauxyTurbUyVx=wTauxyTurb*wUyVx; if (Mesh.IsInDefinedArea()) // If the cell is inside the area of calculus... { // Although the boundary and initial condition are symmetric, // the pressure solution, which solves the Poisson equation // iteratively and does not enforce symmetric property, can be // asymmetric. However, such asymmetry is very small if no // turbulence is involved. When turbulence is present, the // asymmetry can be blowed up through eddy visdcosity (Reynold // stress) and eventually affect the mean velocity. To have // symmetric shape, both turbulence model (k_epsilon.F) and final // velocity calculation (accel.F) are enforced for symmetry. // Noted that such asymmetry might be true in real world. // Be cautious to interpret the results using imposed symmetry. // omit calculation of turbulence for obstacle or void cell if (wAC.c()<=1e-6 || wNf.c()==ctEmptyCell) // If cell is closed to flow or empty... { *wK.C=0.0; *wNuT.C=0.0; *wEp.C=0.0; } else if (AWavemakerType==wtSourceRegion && LowTurbRect.IsInside(i, j)) // If wavemaker is a source function and cell is inside the low turbulence rectangle // (and cell is open to flow)... { // Enforce low turbulence level around the source region. *wK.C = 0.5*Sqr(AUTurb); *wNuT.C = 0.1*ANu; *wEp.C = ACMu*Sqr(wK.c())/(wNuT.c()*wPorosity.c()); } else // If wavemaker is not a source function or cell is outside the low turbulence rectangle // (and cell is open to flow)... { double SqrV=SqrMod(VelC); double n=wPorosity.c(); CMu=ACMu; if (AKEModel==keNonLinearEddyVisc) { double smax=(wKn.c()/wEpn.c())*max2(fabs(wUx.c()), fabs(wVy.c())); CMu = (2.0/3.0)/(7.4+2.0*smax); } if (!IsSurface && IsPorousCell && SqrV>1e-32) // If it is a porous cell not at the free surface and velocity is not too small... { // Computes small scale turbulence due to the presence of porous media. EpsPor=39.0*pow((1.0-n), 2.5)*pow(SqrV, 1.5)/(n*AD50[wNpc.c()-1]); KPor=3.7*(1.0-n)*SqrV/sqrt(n); Eps2KPor=Sqr(EpsPor)/KPor; K2EpsPor=Sqr(KPor)/(EpsPor/n); } //if // it is found that for plane jet problem, the wall turbulence // BC at specified pressure boundaries, which is physically // incorrect, ensure the low turbulence at boundaries and thus // provide stable solution for long time simulation. Skipping the // following wall calculation results in numerical instability. // The reason is that zero turbulence gradient treatment as used // here has accumulated errors, especially for no free surface // flow. Better boundary condition for turbulence or larger // computational domain is needed to resolve the problem. // Determines if the cell has a solid interface, and its characteristics. TVector ObsSrfce; double CDist, TRDist; bool IsObsIntrfce=IsObstacleInterface(Mesh, wAC, wAR, wAT, KlRigid, KrRigid, KtRigid, KbRigid, ObsSrfce, CDist, TRDist); if (IsObsIntrfce) // If cell is near an obstacle... { // Apply wall-bc for fluid-solid interface if (!IsPorousCell) // If cell is not a porous one... { // Calculate a unitary vector tangent to the outward normal from the solid surface. // We choose the one that is 90ยบ clockwise from the normal one. The direction of this // vector is important because we store the frictional velocity with its sign. TVector ObsSrfceTan=NormalVector(ObsSrfce); // Find the projection of V at the center of the cell along the direction parallel to the surface. double VC=DotProduct(ObsSrfceTan, VelC); // For the frictional velocity at the corner uses the value of the velocity interpolated at // this point. Original COBRAS code used the value of velocity at the R or T points (depending // whether the obstacle was horizontal or vertical), this is more correct, but can't be used // with obstacles not parallel to mesh lines. TVector VelTR(wUAR.ITR(), wVAT.ITR()); // CHECK this: // TVector VelTR(wU.ITR(), wV.ITR()); // Find the projection of V at TR along the direction parallel to the surface. double VTR=DotProduct(ObsSrfceTan, VelTR); // CHECK this: // This will be deleted as soon the rest of the code is updated. *wUtmp.C=VelC.x; *wVtmp.C=VelC.y; if (fabs(VC)<1.0e-6) // If the projection of the velocity along the surface is really small... { *wK.C = 0.5*Sqr(AUTurb); *wNuT.C = 0.1*ANu; *wEp.C = ACMu*Sqr(wK.c())/(wNuT.c()*wPorosity.c()); } else // If the projection of the velocity along the surface is big enough... { FricV=CalcFrictionVelocity(E, VC, CDist, ANu, Kappa, IsRoughSurface, Ks); // Constant CMu is used along boundary to avoid singularity (Rodi, 1980; Pengzhi thesis, pg. 22) *wK.C = Sqr(FricV)/sqrt(ACMu); *wEp.C = Cube(fabs(FricV))/(Kappa*CDist); *wNuT.C = ACMu*Sqr(wK.c())/wEp.c(); } //if } //if else // If cell is a porous one... { // Specify k and eps value near the porous-solid interface. *wK.C = KPor; *wEp.C = EpsPor; // Constant CMu is used along boundary to avoid singularity (Rodi, 1980; Pengzhi thesis, pg. 22) *wNuT.C = ACMu*K2EpsPor; } //if } else // For cells not in the interface fluid/obstacle/mesh boundaries... { // Solve the explicit KEps transport equations. // In the following code, Gamma is a coefficient that controls the relative importance // between the central difference and the upwind scheme (similar to Alpha in the // momentum advection). // Epsilon equation // ---------------- // Compute the advective term double fe=DotProduct(VelC, TVector(EpnPorDx, EpnPorDy)); TVector Epnxxyy=TVector(C_Dx(wEpnx), C_Dy(wEpny)); // Compute viscous term (uses a term to correct the negative diffusion by FTCS scheme). double VisE=n*C_Div(wEpnPorDNutSigENux, wEpnPorDNutSigENuy)+ 0.5*(1.0-Gamma)*ADelT*DotProduct(SqrVelC, Epnxxyy); // Compute the production term (see Lemos p96-97) // Assumes large scale turbulence generation is null in porous media. if (!IsPorousCell) // If it isn't a porous media... { Prod=(wTauxxTurb.c()*wUx.c()+wTauyyTurb.c()*wVy.c()+wTauxyTurbUyVx.IC(ValNod5))/n; } //if // There are some cases where negative production is not allowed. // - Near free surface. // - In cells with very small turbulence where molecular viscous effect is important and // turbulence calculation is inaccurate. if (Prod<0.0 && (IsSurface || wK.c()<0.5*Sqr(AUTurb*5.0))) { Prod=0.0; } //if TotProd=Prod+EpsPor*n; // Compute time scale of turbulence over that of mean stress ttots = (wK.c()/wEp.c())*sqrt(2.0*(Sqr(wUx.c())+Sqr(wVy.c()))+Sqr(wUyVx.IC())); // Compute the new Ep oldEp=wEp.c(); c1ee=AC1e+AC1ea-AC1ea/(1.0+AC1eb*ttots); c2ee=AC2e; // when RNG theory is used, c2e has chance to go negative when // S*k/epsilon is large; Adapt explicit scheme to avoid // numerical instability when this (rarely) happens. double DEpDt=-fe+VisE+c1ee*wEpn.c()/wKn.c()*Prod-c2ee*Sqr(wEpn.c())/wKn.c(); if (c2ee>=0.0) DEpDt=DEpDt+c2ee*Eps2KPor*n; *wEp.C=wEpn.c()+ADelT*DEpDt; *wEp.C=(wEpn.c()+ADelT*(-fe+VisE+c1ee*wEpn.c()/wKn.c()*Prod+c2ee*Eps2KPor*n))/ ((1.0+ADelT*c2ee*wEpn.c()/wKn.c())); // ERROR CONDITIONS. if (wEp.c()<1.0e-16) { printf("Epsilon too small \n "); *wEp.C=oldEp; } //if if (wEp.c()>1.0e16 || wKn.c()<1.0e-16) { printf("Check turbulence parameters! \n"); *wEp.C=oldEp; } // K equation // ---------- double fk=DotProduct(VelC, TVector(KPorDx, KPorDy)); TVector Knxxyy=TVector(C_Dx(wKnx), C_Dy(wKny)); // Compute viscous term (uses a term to correct the negative diffusion by FTCS scheme). double VisK=n*C_Div(wKPorDNutSigKNux, wKPorDNutSigKNuy)+ 0.5*(1.0-Gamma)*ADelT*DotProduct(SqrVelC, Knxxyy); oldK=wK.c(); double DKDt=-fk+VisK+Prod-wEp.c()+n*EpsPor; // if RNG is used, eps can grow very rapidly due to negative c2ee // use semi-implicit method to calculate k to ensure stability if (wK.c()<=0.0 && c2ee<=0.0) { DKDt=-fk+VisK+Prod; printf("rng, negative k using explicit method \n"); } //if *wK.C=wKn.c()+ADelT*DKDt; // ERROR CONDITIONS. if (wK.c()<1.0e-16) { printf("k too small \n"); *wK.C=oldK; } //if if (wK.c()<=1.0e-16 || wEp.c()<=1.0e-16) { printf("*** Check for scheme! \n"); *wK.C=0.5*Sqr(AUTurb); *wNuT.C=5.0*ANu*Sqr(ARatio); *wEp.C=ACMu*Sqr(wK.c())/(wNuT.c()*n); } //if *wNuT.C=CMu*Sqr(wK.c())/(wEp.c()*n); // enforcing damping function so that the diffusion of k, eps // will be calculated in consistent with stress tau calculation // (k_epsilon model only) // Lemos (1992), pp. 27, Harlow & Nakayama, 1967 (Physics of Fluids). // Cd is modified locally to reduce effective eddy // viscosity for low turbulence. *wNuT.C=Sqr(wNuT.c())*(1.0-exp(-AEddyCoef*ANu/max2(1.0e-20L, wNuT.c())))/(ANu*AEddyCoef); // monitoring the stability condition if (!IsPorousCell) { // Saves information about the cell with the greatest NuT. if (wNuT.c()>=xmaxxnut) { xmaxxnut=wNuT.c(); xmaxxk=wK.c(); xmaxxep=wEp.c(); prodmax=Prod; ni=i; nj=j; } //if } //if // Compute the maximum DelT set by turbulence. if (TotProd<0.0) { // for k equation, it is linear decay; for epsilon equation, // it is exponential decay! calculated differently Acrest = min3(Acrest, 0.5*wK.c()/(-TotProd+wEp.c()),wK.c()/(-AC1e*TotProd+AC2e*wEp.c())); } //if } //if } //if } //if } while(Mesh.Next2()); } //Just a block // Set BC for K, Ep and NuT. // ------------------------- { //Just a for stablishing the scope of the following variables block // Creates the mesh object. TMesh Mesh(AIMax, AJMax, AXCoefs, AYCoefs, AXi, AYj, ADelX, ADelY); // Creates the masks. TCMask wAC (Mesh, AAC); TCMask wK (Mesh, AK); TCMask wEp (Mesh, AEp); TCMask wKn (Mesh, AKn); TCMask wEpn (Mesh, AEpn); TCMask wNuT (Mesh, AXNut); TCMask wPorosity (Mesh, APorosity); TIntMask wNf (Mesh, ANf); TCMask wOldK (Mesh); TCMask wOldEp (Mesh); TCMask wOldNuT (Mesh); Mesh.DefineArea(2, AIMax-1, 2, AJMax-1); Mesh.First2(); do { int i=Mesh.I; int j=Mesh.J; // Stores old K,Ep & NuT Values. wOldK=wK; wOldEp=wEp; wOldNuT=wNuT; if (Mesh.IsInDefinedArea()) // If the cell is inside the area of calculus... { // Special treatment on free surface if (wAC.c()>1e-6 && (wNf.c()>=ctRightFreeSurface && wNf.c()<=ctIsolatedCell) && (AKEModel==ke1EqKLModel || AKEModel==keReynoldsTensor)) // If cell is open to flow and it is at the free surface or an isolated cell // (and in the case of the specified KEModels)... { KMax= max4(wK.r(), wK.l(), wK.t(), wK.b()); EpMax=max4(wEp.r(), wEp.l(), wEp.t(), wEp.b()); // Set the value of K-Epsilon at the free surface cells coping it from the // neighbouhooding cells. // // Warning: Notice that there's no "break" instruction after every "case". the "break" is // inside an if clause. This is intentional and simulates the behavior of the // original code. switch(wNf.c()) // Depending on the type of fluid interface at the cell... { case ctRightFreeSurface: if (wAC.l()>EM6 && wNf.l()==ctCellFullOfFluid) { *wK.C= wOldK.l(); *wEp.C= wOldEp.l(); *wNuT.C=wOldNuT.l(); break; } //if case ctLeftFreeSurface: if (wAC.r()>EM6 && wNf.r()==ctCellFullOfFluid) { *wK.C= wOldK.r(); *wEp.C= wOldEp.r(); *wNuT.C=wOldNuT.r(); break; } //if case ctBottomFreeSurface: if (wAC.t()>EM6 && wNf.t()==ctCellFullOfFluid) { *wK.C= wOldK.t(); *wEp.C= wOldEp.t(); *wNuT.C=wOldNuT.t(); break; } //if case ctTopFreeSurface: if (wAC.b()>EM6 && wNf.b()==ctCellFullOfFluid) { *wK.C= wOldK.b(); *wEp.C= wOldEp.b(); *wNuT.C=wOldNuT.b(); break; } //if *wK.C= KMax; *wEp.C= EpMax; *wNuT.C=ACMu*Sqr(wK.c())/(wEp.c()*wPorosity.c()); break; case ctIsolatedCell: *wK.C= wKn.c(); *wEp.C= wEpn.c(); *wNuT.C=ACMu*Sqr(wK.c())/(wEp.c()*wPorosity.c()); break; } //switch if (wEp.c()<1.0e-6 || wK.c()<1.0e-6) // If K or Eps are too small... { printf("Eps or K small \n"); *wK.C= wOldK.c(); *wEp.C= wOldEp.c(); *wNuT.C=wOldNuT.c(); } //if } //if // BC at mesh sides. if ((i==2) && (AKl!=bcSpecInflowOutFlow)) // Left side (and not flow specified at it)... { wK.Set_l(wK.c()); wEp.Set_l(wEp.c()); wNuT.Set_l(wNuT.c()); } //if else if (i==AIMax-1) // Right side... { wK.Set_r(wK.c()); wEp.Set_r(wEp.c()); wNuT.Set_r(wNuT.c()); } //if if (j==2) // Bottom side... { wK.Set_b(wK.c()); wEp.Set_b(wEp.c()); wNuT.Set_b(wNuT.c()); } //if else if (j==AJMax-1) // Top side... { wK.Set_t(wK.c()); wEp.Set_t(wEp.c()); wNuT.Set_t(wNuT.c()); } //if } //if } while(Mesh.Next2()); } //Just a block }