////////////////////////////////////////////////////////////////////////
// THIS FILE IS THE INTERFACE FILE BETWEEN MATLAB AND C //
////////////////////////////////////////////////////////////////////////
// Authors: Matt Ueckermann, Themis Sapsis,
// 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>
////////////////////////////
// MASTER FUNCTION FILE //
////////////////////////////
//////////////////////////////////////////
// MATLAB INTERFACE FUNCTION DEFINITION //
//////////////////////////////////////////
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
// Declare input variables
double dx, dy, *ui, *vi, *uj, *vj;
int Nx, Ny, *Nodeu, *Nodev, Nbcs, NbcsDu, NbcsDv;
// Declare output variables
double *Fu, *Fv;
//Declare working variables
int i, j;
double fw, fe, fn, fs;
// Some minor input/output error checking
if (!(nrhs==13)) {
mexErrMsgTxt("Need 13 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, 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] );
ui = (double*) mxGetPr( prhs[4] );
vi = (double*) mxGetPr( prhs[5] );
uj = (double*) mxGetPr( prhs[6] );
vj = (double*) mxGetPr( prhs[7] );
Nodeu = (int*) mxGetPr( prhs[8] );
Nodev = (int*) mxGetPr( prhs[9] );
Nbcs = (int) mxGetScalar( prhs[10] );
NbcsDu = (int) mxGetScalar( prhs[11] );
NbcsDv = (int) mxGetScalar( prhs[12] );
// idsu = (int*) mxGetPr( prhs[10] );
// Nidsu = mxGetM ( prhs[10] );
// idsv = (int*) mxGetPr( prhs[11] );
// Nidsv = mxGetM ( prhs[11] );
// Example of outputting text to MATLAB
// if (time ==0)
// mexPrintf("Welcome message for time %f\n",time);
// ALLOCATE AND POINT TO OUTPUTS
plhs[0] = mxCreateDoubleMatrix((Nx-1)* Ny , 1, mxREAL);
plhs[1] = mxCreateDoubleMatrix( Nx *(Ny-1), 1, mxREAL);
Fu = (double*)mxGetPr(plhs[0]);
Fv = (double*)mxGetPr(plhs[1]);
for(i=1;i<Nx;i++) {
for(j=1;j<Ny+1;j++) {
if (Nodeu[j+(i)*(Ny+2)]<=Nbcs) // don't do work
{
Fu[j-1+(i-1)*Ny] = 0.0;
}
else{ //do work
// West Face
if(Nodeu[j +(i-1)*(Ny+2)]>Nbcs) {
fw = (0.25/dx)*( ui[Nodeu[j +(i-1)*(Ny+2)]-1] + ui[Nodeu[j +i*(Ny+2)]-1] )\
*( uj[Nodeu[j +(i-1)*(Ny+2)]-1] + uj[Nodeu[j +i*(Ny+2)]-1] );
}
else //treat boundary condition
{
if (Nodeu[j +(i-1)*(Ny+2)]>NbcsDu) { //Neumann
fw = (1.0/dx)*( 0.5*ui[Nodeu[j +(i-1)*(Ny+2)]-1] + ui[Nodeu[j +i*(Ny+2)]-1] )\
*( 0.5*uj[Nodeu[j +(i-1)*(Ny+2)]-1] + uj[Nodeu[j +i*(Ny+2)]-1] );
}
else { //same as before - Dirichlet
fw = (0.25/dx)*( ui[Nodeu[j +(i-1)*(Ny+2)]-1] + ui[Nodeu[j +i*(Ny+2)]-1] )\
*( uj[Nodeu[j +(i-1)*(Ny+2)]-1] + uj[Nodeu[j +i*(Ny+2)]-1] );
}
}
//East face
if (Nodeu[j +(i+1)*(Ny+2)]>Nbcs) {
fe = (0.25/dx)*( ui[Nodeu[j +(i+1)*(Ny+2)]-1] + ui[Nodeu[j +i*(Ny+2)]-1] )\
*( uj[Nodeu[j +(i+1)*(Ny+2)]-1] + uj[Nodeu[j +i*(Ny+2)]-1] );
}
else //treat boundary condition
{
if (Nodeu[j +(i+1)*(Ny+2)]>NbcsDu) { //Neumann
//mexPrintf("Neumann %d,%d, %f\n",Nodeu[j+(i+1)*(Ny+2)]-1,NbcsDu,uj[Nodeu[j +(i+1)*(Ny+2)]-1] );
fe = (1.0/dx)*( 0.5*ui[Nodeu[j +(i+1)*(Ny+2)]-1] + ui[Nodeu[j +i*(Ny+2)]-1] )\
*( 0.5*uj[Nodeu[j +(i+1)*(Ny+2)]-1] + uj[Nodeu[j +i*(Ny+2)]-1] );
}
else {//same as before - Dirichlet
//mexPrintf("Dirichlet %d,%d, %f\n",Nodeu[j+(i+1)*(Ny+2)]-1,NbcsDu,uj[Nodeu[j +(i+1)*(Ny+2)]-1] );
fe = (0.25/dx)*( ui[Nodeu[j +(i+1)*(Ny+2)]-1] + ui[Nodeu[j +i*(Ny+2)]-1] )\
*( uj[Nodeu[j +(i+1)*(Ny+2)]-1] + uj[Nodeu[j +i*(Ny+2)]-1] );
}
}
//South Face
if (Nodeu[j+1+ i*(Ny+2)]>Nbcs) {
fs = (0.25/dy)*( ui[Nodeu[j+1 + i*(Ny+2) ]-1] + ui[Nodeu[j +i*(Ny+2)]-1] )\
*( vj[Nodev[j +(i+1)*(Ny+1)]-1] + vj[Nodev[j +i*(Ny+1)]-1] );
}
else { //treat boundary conditions
fs = 0.0;
if (Nodeu[j+1+(i)*(Ny+2)]>NbcsDu) { //Neumann on u
fs = 0.5*ui[Nodeu[j+1 + i*(Ny+2) ]-1] + ui[Nodeu[j +i*(Ny+2)]-1];
}
else { //same as before on u - Dirichlet
fs = 0.5*(ui[Nodeu[j+1 + i*(Ny+2) ]-1] + ui[Nodeu[j +i*(Ny+2)]-1]);
}
if (Nodev[j +(i+1)*(Ny+1)]>NbcsDv && Nodev[j +(i+1)*(Ny+1)]<=Nbcs) { //Neumann on v
fs = (0.5/dy)*fs*( vj[Nodev[j +(i+1)*(Ny+1)]-1] + vj[Nodev[j +i*(Ny+1)]-1]\
+ vj[Nodev[j-1 +(i+1)*(Ny+1)]-1] + vj[Nodev[j-1 +i*(Ny+1)]-1]);
}
else { //same as before on v - Dirichlet
fs = (0.5/dy)*fs*( vj[Nodev[j +(i+1)*(Ny+1)]-1] + vj[Nodev[j +i*(Ny+1)]-1]);
}
}
//North Face
if (Nodeu[j-1+ i*(Ny+2)]>Nbcs) {
fn = (0.25/dy)*( ui[Nodeu[j-1 + i*(Ny+2) ]-1] + ui[Nodeu[j +i*(Ny+2)]-1] )\
*( vj[Nodev[j-1 +(i+1)*(Ny+1)]-1] + vj[Nodev[j-1+i*(Ny+1)]-1] );
}
else { //treat boundary conditions
fn = 0.0;
if (Nodeu[j-1+(i)*(Ny+2)]>NbcsDu) { //Neumann on u
fn = ( 0.5*ui[Nodeu[j-1 + i*(Ny+2) ]-1] + ui[Nodeu[j +i*(Ny+2)]-1] );
}
else { //Dirchlet on u -- same as before
fn = 0.5*( ui[Nodeu[j-1 + i*(Ny+2) ]-1] + ui[Nodeu[j +i*(Ny+2)]-1] );
}
if (Nodev[j-1 +(i+1)*(Ny+1)]>NbcsDv && Nodev[j-1 +(i+1)*(Ny+1)]<=Nbcs) { //Neumann on v
fn = (0.5/dy)*fn*( vj[Nodev[j +(i+1)*(Ny+1)]-1] + vj[Nodev[j +i*(Ny+1)]-1]\
+ vj[Nodev[j-1 +(i+1)*(Ny+1)]-1] + vj[Nodev[j-1+i*(Ny+1)]-1]);
}
else { //Dirichlet on v -- same as before
fn = (0.5/dy)*fn*( vj[Nodev[j-1 +(i+1)*(Ny+1)]-1] + vj[Nodev[j-1+i*(Ny+1)]-1]);
}
}
//U advection, central difference scheme
Fu[j-1+(i-1)*Ny] = fw-fe+fs-fn;
}
}
}
for(i=1;i<Nx+1;i++) {
for(j=1;j<Ny;j++) {
if (Nodev[j +(i)*(Ny+1)]<=Nbcs) // don't do work
{
Fv[j-1+(i-1)*(Ny-1)] = 0.0;
}
else{ //do work
//V advection, central difference sceme
if(Nodev[j +(i-1)*(Ny+1)]>Nbcs) {
//West face
fw = (0.25/dx)*( vi[Nodev[j +(i-1)*(Ny+1)]-1] + vi[Nodev[j +i*(Ny+1)]-1] )\
*( uj[Nodeu[j+1 +(i-1)*(Ny+2)]-1] + uj[Nodeu[j +(i-1)*(Ny+2)]-1] );
}
else //treat boundary condition
{
fw=0.0;
//West face
if (Nodev[j +(i-1)*(Ny+1)]-1<=NbcsDv) { //Dirchlet on v --same as before
fw = 0.5*( vi[Nodev[j +(i-1)*(Ny+1)]-1] + vi[Nodev[j +i*(Ny+1)]-1] );
}
else { //Neuman on v
fw = (0.5*vi[Nodev[j +(i-1)*(Ny+1)]-1] + vi[Nodev[j +i*(Ny+1)]-1] );
}
if (Nodeu[j+1 +(i-1)*(Ny+2)]<=NbcsDu || Nodeu[j+1 +(i-1)*(Ny+2)]>Nbcs) { //Dirchlet on u -- same as before
fw = (0.5/dx)*fw*( uj[Nodeu[j+1 +(i-1)*(Ny+2)]-1] + uj[Nodeu[j +(i-1)*(Ny+2)]-1] );
}
else { //Neuman on u
fw = (0.5/dx)*fw\
*( uj[Nodeu[j+1 +(i-1)*(Ny+2)]-1] + uj[Nodeu[j +(i-1)*(Ny+2)]-1]+ \
uj[Nodeu[j+1 +(i )*(Ny+2)]-1] + uj[Nodeu[j +(i )*(Ny+2)]-1]);
}
}
if (Nodev[j +(i+1)*(Ny+1)]>Nbcs) {
//East face
fe = (0.25/dx)*( vi[Nodev[j +(i+1)*(Ny+1)]-1] + vi[Nodev[j +i*(Ny+1)]-1] )\
*( uj[Nodeu[j+1 + i*(Ny+2) ]-1] + uj[Nodeu[j +i*(Ny+2)]-1] );
}
else //treat boundary condition
{
fe=0.0;
//East face
if (Nodev[j +(i+1)*(Ny+1)]-1<=NbcsDv) { //Dirichlet on v --same as before
fe = 0.5*( vi[Nodev[j +(i+1)*(Ny+1)]-1] + vi[Nodev[j +i*(Ny+1)]-1]);
}
else { //Neumann on v
fe = ( 0.5*vi[Nodev[j +(i+1)*(Ny+1)]-1] + vi[Nodev[j +i*(Ny+1)]-1]);
}
if (Nodeu[j+1 +(i)*(Ny+2)]<=NbcsDu || Nodeu[j+1 +(i)*(Ny+2)]>Nbcs) { //Dirchlet on u -- same as before
fe = (0.5/dx)*fe*( uj[Nodeu[j+1 + i*(Ny+2) ]-1] + uj[Nodeu[j +i*(Ny+2)]-1] );
}
else { //Neuman on u
fe = (0.5/dx)*fe\
*( uj[Nodeu[j+1+ i*(Ny+2) ]-1] + uj[Nodeu[j+ i*(Ny+2)]-1] +\
uj[Nodeu[j+1+(i-1)*(Ny+2) ]-1] + uj[Nodeu[j+(i-1)*(Ny+2)]-1]);
}
}
if (Nodev[j+1+ i*(Ny+1)]>Nbcs) {
//South Face
fs = (0.25/dy)*( vi[Nodev[j+1 + i*(Ny+1) ]-1] + vi[Nodev[j +i*(Ny+1)]-1] )\
*( vj[Nodev[j+1 + i*(Ny+1) ]-1] + vj[Nodev[j +i*(Ny+1)]-1] );
}
else { //Treat boundary condition
if (Nodev[j+1+(i)*(Ny+1)]-1<=NbcsDv) { //Dirchlet on v -- same as before
fs = (0.25/dy)*( vi[Nodev[j+1 + i*(Ny+1) ]-1] + vi[Nodev[j +i*(Ny+1)]-1] )\
*( vj[Nodev[j+1 + i*(Ny+1) ]-1] + vj[Nodev[j +i*(Ny+1)]-1] );
}
else { //Neuman on v
fs = (1.0/dy)*( 0.5*vi[Nodev[j+1 + i*(Ny+1) ]-1] + vi[Nodev[j +i*(Ny+1)]-1] )\
*( 0.5*vj[Nodev[j+1 + i*(Ny+1) ]-1] + vj[Nodev[j +i*(Ny+1)]-1] );
}
}
//North Face
if (Nodev[j-1+ i*(Ny+1)]>Nbcs) {
fn = (0.25/dy)*( vi[Nodev[j-1 + i*(Ny+1) ]-1] + vi[Nodev[j +i*(Ny+1)]-1] )\
*( vj[Nodev[j-1 + i*(Ny+1) ]-1] + vj[Nodev[j +i*(Ny+1)]-1] );
}
else { //treat boundary condition
if (Nodev[j-1+(i)*(Ny+1)]-1<=NbcsDv) { //Dirichlet --same as before
fn = (0.25/dy)*( vi[Nodev[j-1 + i*(Ny+1) ]-1] + vi[Nodev[j +i*(Ny+1)]-1] )\
*( vj[Nodev[j-1 + i*(Ny+1) ]-1] + vj[Nodev[j +i*(Ny+1)]-1] );
}
else { //Neuman
fn = (1.0/dy)*( 0.5*vi[Nodev[j-1 + i*(Ny+1) ]-1] + vi[Nodev[j +i*(Ny+1)]-1] )\
*( 0.5*vj[Nodev[j-1 + i*(Ny+1) ]-1] + vj[Nodev[j +i*(Ny+1)]-1] );
}
}
Fv[j-1+(i-1)*(Ny-1)] = fw-fe+fs-fn;
}
}
}
}