/*==================================================================================== ! Purpose(s): ! ! Contain(s): ! ! Author(s): Javier Lopez Lara (jav.lopez@unican.es) ! Gabriel Barajas (barajasg@unican.es) !====================================================================================*/ #include "CMisc.h" #include "CMathMisc.h" #include "CErrors.h" #include #ifdef __INTEL_COMPILER #define PYM cdefinelowturbrect #else #define PYM cdefinelowturbrect_ #endif int IdxPack(int AL, int AR, int AT, int AB) { if (AL<0 || AR<0 || AT<0 || AB<0 || AL>255 || AR>255 || AT>255 || AB>255) { printf("Invalid porosity index \n"); return 0; } return AL+(AR << 8)+(AT << 16)+(AB << 24); } //---------------------------------------------------------------------------------------- void IdxUnpack(int AValue, int &AL, int &AR, int &AT, int &AB) { AL=AValue & 255; AR=(AValue >> 8) & 255; AT=(AValue >> 16) & 255; AB=(AValue >> 24) & 255; } //---------------------------------------------------------------------------------------- //---------------------------------------------------------------------------------------- //---------------------------------------------------------------------------------------- bool IsObstacleInterface(const TMesh &AMesh, TMask &wAc, TMask &wAr, TMask &wAt, bool AKlRigid, bool AKrRigid, bool AKtRigid, bool AKbRigid, TVector &ANormVec, double &ACDist, double &ATRDist) /*!======================================================================= ! Purpose(s): This function returns true if the cell has an obstacle interface. In this case it also returns: ! - ANormVec: an unitary vector perpendicular to this interface (outward from the obstacle). ! - ACDist: Distance from the central point to the obstacle. ! - ATRDist: Distance from top-right point to the obstacle. ! ! Author(s): Javier Lopez Lara (jav.lopez@unican.es) ! Gabriel Barajas (barajasg@unican.es) !=======================================================================*/ // Desc: { bool Res=false; double TotSideOp=wAt.c()+wAt.b()+wAr.c()+wAr.l(); double TotOp=wAc.r()+wAc.l()+wAc.t()+wAc.b(); double provDistMult=1.0; if (AMesh.I==2 && AKlRigid) // If cell is at rigid left boundary... { ANormVec=TVector(1.0, 0.0); Res=true; } else if (AMesh.I==AMesh.Width-1 && AKrRigid) // If cell is at rigid right boundary... { ANormVec=TVector(-1.0, 0.0); Res=true; } else if (AMesh.J==2 && AKbRigid) // If cell is at rigid bottom boundary... { ANormVec=TVector(0.0, 1.0); Res=true; } else if (AMesh.J==AMesh.Height-1 && AKtRigid) // If cell is at rigid top boundary... { ANormVec=TVector(0.0, -1.0); Res=true; } else if (wAc.c()==1.0 && TotOp<3.01) // If cell is open to flow and near cells close to flow (we search for almost one fully closed, thats why // we do TotOp<3.01 instead of TotOp<4.0, hopefully this will close more the search). { if (wAc.l()<1e-6) // If cell has a closed cell at its left... { ANormVec=TVector(1.0, 0.0); Res=true; } else if (wAc.r()<1e-6) // If cell has a closed cell at its right... { ANormVec=TVector(-1.0, 0.0); Res=true; } else if (wAr.b()<1e-6) // If cell has a closed cell at its bottom... { ANormVec=TVector(0.0, 1.0); Res=true; } else if (wAr.t()<1e-6) // If cell has a closed cell at its top... { ANormVec=TVector(0.0, -1.0); Res=true; } //if..else if... } //---------------------------------------------------------------------------------------- else if (wAc.c()==1.0 && TotOp<4.0) // If cell is open to flow and near cells close to flow (we search for almost one fully closed, thats why // we do TotOp<3.01 instead of TotOp<4.0, hopefully this will close more the search). { if (wAc.b()<1.0) // If cell has a partially closed cell at its bottom... { ANormVec=TVector(0.0, 1.0); Res=true; provDistMult=3.0; } else if (wAc.t()<1.0) // If cell has a partially closed cell at its top... { ANormVec=TVector(0.0, -1.0); Res=true; provDistMult=3.0; } //if..else if... else if (wAc.l()<1.0) // If cell has a partially closed cell at its left... { ANormVec=TVector(1.0, 0.0); Res=true; provDistMult=3.0; } else if (wAc.r()<1.0) // If cell has a partially closed cell at its right... { ANormVec=TVector(-1.0, 0.0); Res=true; provDistMult=3.0; } } //---------------------------------------------------------------------------------------- else if (TotSideOp<4.0) // If the cell is partially opened to flow... { if (TotSideOp<=2.0) // If sides are more closed than opened... { if (max2(wAr.c(), wAr.l())>max2(wAt.c(), wAt.b())) // If the obstacle is more vertical than horizontal... { if (wAr.c()min2(wAt.c(), wAt.b())) // If the obstacle is more horizontal than vertical... { if (wAt.c()>wAt.b()) // If the obstacle is at the bottom... { ANormVec=TVector(0.0, 1.0); Res=true; } else // If the obstacle is at the top... { ANormVec=TVector(0.0, -1.0); Res=true; } //if } else // If the obstacle is more vertical than horizontal... { if (wAr.c()>wAr.l()) // If the obstacle is at the left... { ANormVec=TVector(1.0, 0.0); Res=true; } else // If the obstacle is at the right... { ANormVec=TVector(-1.0, 0.0); Res=true; } //if } //if } //if } //if // The distance is calculated as half the absolute value of the projection of normal vector and vector // (ADelX, ADelY). At this moment, both distances are considered equal. ACDist=provDistMult*0.5*fabs(DotProduct(TVector(AMesh.DelX, AMesh.DelY), ANormVec)); ATRDist=ACDist; return Res; }; //---------------------------------------------------------------------------------------- extern "C" void _FORTRAN PYM (_INT_ARRAY ARect, _BOOL AIsExNonTurbArea, _INT AISourceS, _INT AISourceE, _INT AJSourceS, _INT AJSourceE, _INT ALowTurbLength) { /*!======================================================================= ! Purpose(s): Defines the low turbulence rectangle. ARect is an integer array of size 4. The reason for using ! this array and not a TIntRect is because we need to store the rectangle in a Fortran's array. ! The order of the values in the array is: L, R, B, T ! ! Author(s): Javier Lopez Lara (jav.lopez@unican.es) ! Gabriel Barajas (barajasg@unican.es) !=======================================================================*/ const int LowTurbWidth = 5; if (AIsExNonTurbArea) // If region occupies everything from left side to ALowTurbLength... { ARect[0]=0; ARect[1]=ALowTurbLength; ARect[2]=-1; ARect[3]=100000000; // This number should be bigger than any possible JMax. } else // If region must surround the source... { ARect[0]=AISourceS-LowTurbWidth; ARect[1]=AISourceE+LowTurbWidth; ARect[2]=AJSourceS-LowTurbWidth; ARect[3]=AJSourceE+LowTurbWidth; } //if }; //----------------------------------------------------------------------------------------