////////////////////////////////////////////////////////////////////////
// THIS FILE IS THE INTERFACE FILE BETWEEN MATLAB AND C //
////////////////////////////////////////////////////////////////////////
// Authors: Matt Ueckermann, Pierre Lermusiaux, Pat Haley for MIT course 2.29
////////////////////////////////////
// DECLARE STANDARD HEADER FILES //
////////////////////////////////////
#include "math.h"
#include "mex.h"
#include <stdlib.h>
#include <string.h>
#define MIN(a, b) ((a)>(b)?(b):(a))
#define MAX(a, b) ((a)<(b)?(b):(a))
int sign(float v) {
return((v<0)-(v>0));
}
/////////////////////////////
// DEFINE LIMITER FUNCTION //
/////////////////////////////
double TVD(double phi) {
//MC limiter
return(MAX(0.0, MIN(MIN(0.5 * (1.0 + phi), 2.0), 2.0 * phi)));
//minmod limiter
// return(((phi) > (0) ? MIN(phi,1) : 0.0));
//superbee
// return(MAX(MAX(0,MIN(1.0, 2.0 * phi)),MIN(2,phi)));
//Upwind
// return(0.0);
// Lax Wendroff (NOT TVD)
// return(1.0);
}
////////////////////////////
// MASTER FUNCTION FILE //
////////////////////////////
//////////////////////////////////////////
// MATLAB INTERFACE FUNCTION DEFINITION //
//////////////////////////////////////////
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
// Declare input variables
double dx, dy, dt, *u, *v, *rho, NBlen;
int Nx, Ny, *Nodeu, *Nodev, *NodeP, Nbcs, NbcsDP, NbcsDu, NbcsDv;
int *NBids;
// Declare output variables
double *Fu;
//Declare working variables
int i, j;
double fl, fr, fu, fd, theta, U, sumRHO, diffRHO;
double lenx, leny, len, lenxxy, lenxyy;
// Some minor input/output error checking
if (!(nrhs==16)) {
mexErrMsgTxt("Need 16 input arguments");
}
if (!(nlhs==2)) {
mexErrMsgTxt("Need 2 outputs arguments");
}
// Get access to input data from MATLAB
// Pointers are used for large array
// Scalar are re-created in memory
/********!!!! WARNING !!!************/
// WARNING !!!! note that C is expecting to be pointing to type int32 from
// Matlab. Therefore, for this code to work, the C argument must be
// int32. This is done in matlab by int32(C);
// The function is then called like [Bnew Cnew] = axB(5,B,int32(C));
// Use double(Cnew) if you want to convert it back (this might slow things).
/********!!!! WARNING !!!************/
//[Fu Fv] = UVadvect_DO_CDS(Nx, Ny, dx, dy, dt, ui, vi, uj, vj,
// Nodeu, Nodev, idsu, idsv)
Nx = (int) mxGetScalar( prhs[0] );
Ny = (int) mxGetScalar( prhs[1] );
dx = (double) mxGetScalar( prhs[2] );
dy = (double) mxGetScalar( prhs[3] );
dt = (double) mxGetScalar( prhs[4] );
rho = (double*) mxGetPr( prhs[5] );
u = (double*) mxGetPr( prhs[6] );
v = (double*) mxGetPr( prhs[7] );
NodeP = (int*) mxGetPr( prhs[8] );
Nodeu = (int*) mxGetPr( prhs[9] );
Nodev = (int*) mxGetPr( prhs[10] );
Nbcs = (int) mxGetScalar( prhs[11] );
NbcsDP = (int) mxGetScalar( prhs[12] );
NbcsDu = (int) mxGetScalar( prhs[13] );
NbcsDv = (int) mxGetScalar( prhs[14] );
NBlen = (double) mxGetScalar( prhs[15] );
// Example of outputting text to MATLAB
// if (time ==0)
// mexPrintf("Welcome message for time %f\n",time);
// These macros allow easy access to values of rho and u at
// UU
// U
// LL L X R RR
// D
// DD
// with respect to the current cell
#define RHOM rho[NodeP[j + i * (Ny + 2)] - 1]
#define RHOUU rho[NodeP[j - 2 + i * (Ny + 2)] - 1]
#define RHOU rho[NodeP[j - 1 + i *(Ny + 2)] - 1]
#define RHOD rho[NodeP[j + 1 + i * (Ny + 2)] - 1]
#define RHODD rho[NodeP[j + 2 + i * (Ny + 2)] - 1]
#define RHOLL rho[NodeP[j + (i - 2) * (Ny + 2)] - 1]
#define RHOL rho[NodeP[j + (i - 1) * (Ny + 2)] - 1]
#define RHOR rho[NodeP[j + (i + 1) * (Ny + 2)] - 1]
#define RHORR rho[NodeP[j + (i + 2) * (Ny + 2)] - 1]
#define RHORD rho[NodeP[j + 1 + (i + 1) * (Ny + 2)] - 1]
#define RHORU rho[NodeP[j - 1 + (i + 1) * (Ny + 2)] - 1]
#define RHOLD rho[NodeP[j + 1 + (i - 1) * (Ny + 2)] - 1]
#define RHOLU rho[NodeP[j - 1 + (i - 1) * (Ny + 2)] - 1]
#define RHORRD rho[NodeP[j + 1 + (i + 2) * (Ny + 2)] - 1]
#define RHORRU rho[NodeP[j - 1 + (i + 2) * (Ny + 2)] - 1]
#define RHOLLD rho[NodeP[j + 1 + (i - 2) * (Ny + 2)] - 1]
#define RHOLLU rho[NodeP[j - 1 + (i - 2) * (Ny + 2)] - 1]
#define RHORDD rho[NodeP[j + 2 + (i + 1) * (Ny + 2)] - 1]
#define RHORUU rho[NodeP[j - 2 + (i + 1) * (Ny + 2)] - 1]
#define RHOLDD rho[NodeP[j + 2 + (i - 1) * (Ny + 2)] - 1]
#define RHOLUU rho[NodeP[j - 2 + (i - 1) * (Ny + 2)] - 1]
#define RHORRDD rho[NodeP[j + 2 + (i + 2) * (Ny + 2)] - 1]
#define RHORRUU rho[NodeP[j - 2 + (i + 2) * (Ny + 2)] - 1]
#define RHOLLDD rho[NodeP[j + 2 + (i - 2) * (Ny + 2)] - 1]
#define RHOLLUU rho[NodeP[j - 2 + (i - 2) * (Ny + 2)] - 1]
#define UU v[Nodev[j - 1 + (i) * (Ny + 1)] - 1]
#define UD v[Nodev[j + (i) * (Ny + 1)] - 1]
#define UL u[Nodeu[j + (i - 1) * (Ny + 2)] - 1]
#define UR u[Nodeu[j + (i)* (Ny + 2)] - 1]
// ALLOCATE AND POINT TO OUTPUTS
plhs[0] = mxCreateDoubleMatrix((Nx)* Ny , 1, mxREAL);
Fu = (double*)mxGetPr(plhs[0]);
plhs[1] = mxCreateNumericMatrix(Nx* Ny , 1, mxINT32_CLASS, mxREAL);
NBids = (int*)mxGetPr(plhs[1]);
// Calculate lengths to determine narrow band region
lenx = dx*NBlen;
leny = dy*NBlen;
len = sqrt(lenx * lenx + leny * leny);
lenxxy = sqrt(4 * lenx * lenx + leny * leny);
lenxyy = sqrt(lenx * lenx + 4 * leny * leny);
//for(i=1;i<Nx+1;i++) {
// for(j=1;j<Ny+1;j++) {
//TO avoid segfaults, only go within 2 CV's of the boundary
for(i=2; i<Nx; i++) {
for(j=2; j<Ny; j++) {
//If this is a boundary node, don't do any work
if (NodeP[j +(i)*(Ny+2)]<=Nbcs) {
Fu[j - 1 + (i - 1) * Ny] = 0.0;
}
else {
//If this is in the narrowband do work
double lenx;
double leny;
double tmpright;
double tmpleft;
double tmptop;
double tmpbottom;
lenx = dx * NBlen;
leny = dy * NBlen;
len = sqrt(lenx * lenx + leny * leny);
if ((fabs(RHOM) < MAX(lenx, leny)) || \
(fabs(RHOU) < leny) || (fabs(RHOUU) < 2 * leny) || \
(fabs(RHOD) < leny) || (fabs(RHODD) < 2 * leny) || \
(fabs(RHOL) < lenx) || (fabs(RHOLL) < 2 * lenx) || \
(fabs(RHOR) < lenx) || (fabs(RHORR) < 2 * lenx) || \
(fabs(RHORD) < len) || (fabs(RHORU) < len) || \
(fabs(RHOLD) < len) || (fabs(RHOLU) < len) || \
(fabs(RHORRD) < lenxxy)||(fabs(RHORRU) < lenxxy)|| \
(fabs(RHOLLD) < lenxxy)||(fabs(RHOLLU) < lenxxy)|| \
(fabs(RHORDD) < lenxyy)||(fabs(RHORUU) < lenxyy)|| \
(fabs(RHOLDD) < lenxyy)||(fabs(RHOLUU) < lenxyy)||\
(fabs(RHORRDD) < 2 * len)||(fabs(RHORRUU) < 2 * len)|| \
(fabs(RHOLLDD) < 2 * len)||(fabs(RHOLLUU) < 2 * len)) {
//Set the narrow band flag to indicate that this cell is
//active
NBids[j-1 + (i-1) * (Ny)] = 1;
// if (RHOR<=RHOM){tmpleft = 1;}
// else {tmpleft = 0;}
// if (RHOR>=RHOM){tmpright = 1;}
// else {tmpright=0;}
// fl = (((1-tmpleft)*fabs(RHOR - RHOM)/dx)+((1-tmpright)*fabs(RHOM-RHOL)/dx))/(MAX(tmpleft+tmpright,1));
// fl = (((1-tmpleft)*fabs(RHOM - RHOL)/dx)+((1-tmpright)*fabs(RHOR-RHOM)/dx))/(MAX(tmpleft+tmpright,1));
// fl = (RHOM - MIN(RHOL,RHOR)) / dx; (Older version)
// fl = (RHOM - MIN(MIN(RHOL,RHOM),MIN(RHOM,RHOR))) / dx; (Works partially)
if (RHOL<=RHOM && RHOM>=RHOR)
{fl = (fabs(RHOM-RHOR) + fabs(RHOM-RHOL))/(2*dx); }
if (RHOL>=RHOM && RHOM>=RHOR)
{fl = (fabs(RHOM-RHOR))/(dx); }
if (RHOL<=RHOM && RHOM<=RHOR)
{fl = (fabs(RHOM-RHOL))/(dx); }
if (RHOL>=RHOM && RHOM<=RHOR)
{fl = 0; }
// Prevent Obstacles from affecting flux at cells adjacent to them//
if (NodeP[j +(i-1)*(Ny+2)]<=Nbcs && NodeP[j +(i+1)*(Ny+2)]<=Nbcs) // If cell to the left and cell to right are within obstacle
{fl = 0;} //Then no flux in x direction in CV
if (NodeP[j +(i-1)*(Ny+2)]<=Nbcs && NodeP[j +(i+1)*(Ny+2)]>Nbcs) // If only cell to the left is within obstacle
{if(RHOR<=RHOM)
{fl = (fabs(RHOM-RHOR))/(dx);}
else
{fl = 0;} //Take flux only if upwind condition satisfied, otherwise flux is zero
} //Then take flux contribution due to the right cell
if (NodeP[j +(i-1)*(Ny+2)]>Nbcs && NodeP[j +(i+1)*(Ny+2)]<=Nbcs)//If only cell to the right is within obstacle
{if(RHOL<=RHOM) //Is information coming towards the obstacle? If yes, it is upwind direction and take flux
{fl = (fabs(RHOM-RHOL))/(dx);}
else
{fl=0;} //Otherwise no flux in this cell
} //Then take flux contribution due to the left cell
// If both sides don't have obstacles, then take flux originally computed.
// if (RHOD<=RHOM){tmptop = 1;}
//else {tmptop = 0;}
// if (RHOD>=RHOM){tmpbottom = 1;}
//else {tmpbottom=0;}
//fd = (((1-tmptop)*fabs(RHOD - RHOM)/dy)+((1-tmpbottom)*fabs(RHOM-RHOU)/dy))/(MAX(tmptop+tmpbottom,1));
//fd = (((1-tmptop)*fabs(RHOM - RHOU)/dy)+((1-tmpbottom)*fabs(RHOD-RHOM)/dy))/(MAX(tmptop+tmpbottom,1));
//fd = (RHOM - MIN(MIN(RHOU,RHOM),MIN(RHOM,RHOD))) / dy; (Works partially)
//fd = (RHOM - MIN(RHOU,RHOD)) / dy; (Older version)
if (RHOU<=RHOM && RHOM>=RHOD)
{fd = (fabs(RHOM-RHOU) + fabs(RHOM-RHOD))/(2*dy); }
if (RHOU>=RHOM && RHOM>=RHOD)
{fd = (fabs(RHOM-RHOD))/(dy); }
if (RHOU<=RHOM && RHOM<=RHOD)
{fd = (fabs(RHOM-RHOU))/(dy); }
if (RHOU>=RHOM && RHOM<=RHOD)
{fd = 0; }
// Prevent Obstacles from affecting flux at cells adjacent to them//
if (NodeP[j+1 +(i)*(Ny+2)]<=Nbcs && NodeP[j-1 +(i)*(Ny+2)]<=Nbcs) // If cell to the top and cell to bottom are within obstacle
{fd = 0;} //Then no flux in y direction in CV
if (NodeP[j+1 +(i)*(Ny+2)]<=Nbcs && NodeP[j-1 +(i)*(Ny+2)]>Nbcs) // If only cell to the bottom is within obstacle
{if(RHOU<=RHOM)
{fd = (fabs(RHOM-RHOU))/(dy);}//Is it upwind? Then take flux contribution due to the upper cell
else
{fd=0;} //Else, no flux
}
if (NodeP[j+1 +(i)*(Ny+2)]>Nbcs && NodeP[j-1 +(i)*(Ny+2)]<=Nbcs) //If only cell to the top is within obstacle
{if(RHOD<=RHOM)
{fd = (fabs(RHOM-RHOD))/(dy);}
else
{fd=0;}
} //Then take flux contribution due to the lower cell
// If both sides don't have obstacles, then take flux originally computed.
//Calculate Total flux change for this cell
Fu[j-1+(i-1)*Ny] = sqrt(fl * fl + fd * fd);
}
else { //If this is not inside the narrowband, set the flux zero
//and set the narrow band flag to indicate that the cell
//is inactive
Fu[j-1+(i-1)*Ny] = 0.0;
NBids[j-1 +(i-1) * Ny] = 0;
}
}
}
}
}