feat!: remove cpp code

1 files changed, 0 insertions, 501 deletions

diff --git a/src/source.cpp b/src/source.cpp
deleted file mode 100755
index 3c35588..0000000
--- a/src/source.cpp
+++ /dev/null
@@ -1,501 +0,0 @@
-// This is not my code.

-// I have no idea how it works.

-// It probably requires paid software.

-// I'm planning my own implementation in Python.

-#include <stdlib.h>

-#include <stdio.h>

-#include <time.h>

-#include <cmath>

-#include "TECIO.h"

-

-#ifndef NULL

-#define NULL 0

-#endif

-

-enum FileType { FULL = 0, GRID = 1, SOLUTION = 2 };

-

-#define BD 20

-#define BX 50

-#define BY 10

-#define NX 2000

-#define NY 200

-#define KT 10

-#define AI -2./3.

-#define DT 0.00001

-#define RX 100

-#define RY 40

-

-#define IY 40

-#define IX 80

-

-#define Re 400

-#define GR 0.0001

-#define T 6000

-

-double dxx[NX], dyy[NY], xi[NX][NY], yi[NX][NY], a[NX*NY], b[NX*NY], c[NX*NY];

-

-double dx[NX][NY], dy[NX][NY];

-

-

-

-

-void PRINT_PLT(void) {

-	double SolTime;

-	INTEGER4 Debug, I, J, III, DIsDouble, VIsDouble, IMax, JMax, KMax, ZoneType, StrandID, ParentZn, IsBlock;

-	INTEGER4 ICellMax, JCellMax, KCellMax, NFConns, FNMode, ShrConn, FileType;

-

-	Debug = 1;

-	VIsDouble = 1;

-	DIsDouble = 1;

-	IMax = NY;

-	JMax = NX;

-	KMax = 1;

-	ZoneType = 0;      /* Ordered */

-	SolTime = 360.0;

-	StrandID = 0;     /* StaticZone */

-	ParentZn = 0;      /* No Parent */

-	IsBlock = 1;      /* Block */

-	ICellMax = 0;

-	JCellMax = 0;

-	KCellMax = 0;

-	NFConns = 0;

-	FNMode = 0;

-	ShrConn = 0;

-	FileType = FULL;

-

-	/*

-	* Open the file and write the tecplot datafile

-	* header information

-	*/

-	I = TECINI112((char*)"SIMPLE DATASET",

-		(char*)"X Y U V",

-		(char*)"UV.plt",

-		(char*)".",

-		&FileType,

-		&Debug,

-		&VIsDouble);

-

-	/*

-	* Write the zone header information.

-	*/

-	I = TECZNE112((char*)"Simple Zone",

-		&ZoneType,

-		&IMax,

-		&JMax,

-		&KMax,

-		&ICellMax,

-		&JCellMax,

-		&KCellMax,

-		&SolTime,

-		&StrandID,

-		&ParentZn,

-		&IsBlock,

-		&NFConns,

-		&FNMode,

-		0,              /* TotalNumFaceNodes */

-		0,              /* NumConnectedBoundaryFaces */

-		0,              /* TotalNumBoundaryConnections */

-		NULL,           /* PassiveVarList */

-		NULL,           /* ValueLocation = Nodal */

-		NULL,           /* SharVarFromZone */

-		&ShrConn);

-

-	/*

-	* Write out the field data.

-	*/

-	III = IMax * JMax;

-	I = TECDAT112(&III, &dx[0][0], &DIsDouble);

-	I = TECDAT112(&III, &dy[0][0], &DIsDouble);

-	I = TECDAT112(&III, &xi[0][0], &DIsDouble);

-	I = TECDAT112(&III, &yi[0][0], &DIsDouble);

-

-	I = TECEND112();

-

-

-}

-

-

-void PRINT_VECT_VTK(double U[NX][NY], double V[NX][NY], double *x, double *y) {

-	int i, j;

-

-	FILE *vect_field;

-

-	/* File grid.vtk*/

-	if ((vect_field = fopen("a.vtk", "w")) == NULL)

-		printf("The grid file could not be open\n");

-	else{

-		fprintf(vect_field, "%d\n", NX);

-		fprintf(vect_field, "%d\n", NY);

-

-		for (i = 0; i <= NX - 1; i++)

-			fprintf(vect_field, "%f\n", x[i]);

-

-		for (j = 0; j <= NY - 1; j++)

-			fprintf(vect_field, "%f\n", y[j]);

-

-		for (i = 0; i <= NY - 1; i++)

-			for (j = 0; j <= NX - 1; j++){

-			fprintf(vect_field, "%f\t %f\t\n", U[j][i], V[j][i]);

-			}

-	}

-	fclose(vect_field);

-

-}

-

-void PRINT_VTK(double U[NX][NY], double V[NX][NY], double *x, double *y) {

-	int i, j;

-

-	FILE *vect_field;

-

-

-	/* File grid.vtk*/

-	if ((vect_field = fopen("b.vtk", "w")) == NULL)

-		printf("The grid file could not be open\n");

-	else{

-		fprintf(vect_field, "# vtk DataFile Version 2.0\n");

-		fprintf(vect_field, "Sample rectilinear grid\n");

-		fprintf(vect_field, "ASCII\n");

-		fprintf(vect_field, "DATASET RECTILINEAR_GRID\n");

-		fprintf(vect_field, "DIMENSIONS %d %d \n", NX, NY);

-		fprintf(vect_field, "X_COORDINATES %d double\n", NX);

-		for (i = 0; i <= NX - 1; i++)

-			fprintf(vect_field, "%f\n", x[i]);

-		fprintf(vect_field, "Y_COORDINATES %d double\n", NY);

-		for (j = 0; j <= NY - 1; j++)

-			fprintf(vect_field, "%f\n", y[j]);

-//		fprintf(vect_field, "Z_COORDINATES %d double\n", nz);

-//		fprintf(vect_field, "%f\n", z);

-		fprintf(vect_field, "POINT_DATA %d\n", (NX)*(NY));

-		fprintf(vect_field, "VECTORS vectors double\n");

-		for (i = 0; i <= NY - 1; i++)

-			for (j = 0; j <= NX - 1; j++){

-				fprintf(vect_field, "%f\t %f\t %f\n", U[j][i], V[j][i], 0.0);

-			}

-	}

-	fclose(vect_field);

-

-}

-

-

-void inixh(double *a)

-{

-	int i, j, k;

-	for (i = 0; i < NX; i++) {

-		for (j = 0; j < NY; j++) {

-

-			k = i *NY + j;

-			if (i<RX && j>RY)

-				a[k] = 1.;

-			else

-				a[k] = 0.;

-		}

-	}

-}

-void inih(double *a)

-{

-	int i, j, k;

-

-	for (i = 0; i < NX; i++) {

-		for (j = 0; j < NY; j++) {

-			k = i *NY + j;

-			a[k] = 0.;

-		}

-	}

-}

-void USh(double *a, double *b, double *c, double *d)

-{

-	double u1;

-	int i, j, k;

-	for (i = 0; i < NX; i++) {

-		for (j = 0; j < NY; j++) {

-			k = i *NY + j;

-

-			if (i < NX - 2 && i>0 && j < NY - 1 && j>0 && ((i <= RX && j <= RY) != 1)){

-				u1 = -((pow(a[k + NY], 2) - pow(a[k - NY], 2))*NX / 2.0f + (a[k + 1] * 0.5*(b[k] + b[k + NY]) - a[k - 1] * 0.5*(b[k - 1] + b[k + NY - 1]))*NY*KT / 2.0);

-				u1 = u1 + (1.0 / Re)*((a[k + NY] - 2.0*a[k] + a[k - NY])*NX*NX + (a[k + 1] - 2.0*a[k] + a[k - 1])*NY*NY*KT*KT);

-				u1 = a[k] + DT*(u1 - 1.0*NX*(c[k + NY] - c[k]));

-				d[k] = u1;

-			}

-			else

-				d[k] = a[k];

-

-		}

-	}

-}

-void VSh(double *a, double *b, double *c, double *d)

-{

-

-	double v1;

-	int i, j, k;

-	for (i = 0; i < NX; i++) {

-		for (j = 0; j < NY; j++) {

-			k = i *NY + j;

-

-			if (i < NX - 1 && i>0 && j < NY - 2 && j>0 && ((i <= RX && j <= RY) != 1)){

-				v1 = -((b[k + NY] * 0.5*(a[k] + a[k + 1]) - b[k - NY] * 0.5*(a[k - NY] + a[k - NY + 1]))*NX / 2.0 + (b[k + 1] * b[k + 1] - b[k - 1] * b[k - 1])*KT*NY / 2.0);

-				v1 = v1 + (1.0f / Re)*((b[k + NY] - 2.0f*b[k] + b[k - NY])*NX*NX + (b[k + 1] - 2.0f*b[k] + b[k - 1])*NY*NY*KT*KT);

-				v1 = b[k] + DT*(v1 - (1.0f*NY*KT)*(c[k + 1] - c[k]));

-				d[k] = v1;

-			}

-			else

-				d[k] = b[k];

-		}

-	}

-}

-void sumah(double *a, double *b)

-{

-	int i, j, k;

-	for (i = 0; i < NX; i++) {

-		for (j = 0; j < NY; j++) {

-			k = i *NY + j;

-			b[k] += a[k];

-		}

-	}

-}

-void rubnixh(double *a)

-{

-	int i, j, k;

-	for (i = 0; i < NX; i++) {

-		for (j = 0; j < NY; j++) {

-			k = i *NY + j;

-

-			if (i == 0 && j>RY)

-				a[k] = 1.;

-			if (i == NX - 2)

-				a[k] = a[k - NY];

-			if (j == 0)

-				a[k] = 0.;

-			if (j == NY - 1)

-				a[k] = a[k - 1];

-			if (j == RY && i < RX)

-				a[k] = 0.;

-			if (i == RX && j <= RY)

-				a[k] = a[k + NY] * AI;

-		}

-	}

-}

-void rubniyh(double *a)

-{

-	int i, j, k;

-	for (i = 0; i < NX; i++) {

-		for (j = 0; j < NY; j++) {

-			k = i *NY + j;

-

-			if (j == 0)

-				a[k] = AI*a[k + 1];

-			if (j == NY - 2)

-				a[k] = a[k - 1];

-			if (i == 0)

-				a[k] = 0.;

-			if (i == NX - 1)

-				a[k] = a[k - NY];

-			if (j == RY && i <= RX)

-				a[k] = AI*a[k + 1];

-			if (i == RX && j < RY)

-				a[k] = 0.;

-		}

-	}

-}

-void rubniph(double *a)

-{

-	int i, j, k;

-	for (i = 0; i < NX; i++) {

-		for (j = 0; j < NY; j++) {

-			k = i *NY + j;

-

-

-			if (j == 0 && i != 0 && i != NX - 1)

-				a[k] = a[k + 1];

-			if (j == NY - 1 && i != 0 && i != NX - 1)

-				a[k] = 0.0;

-			if (i == 0 && j != 0 && j != NY - 1)

-				a[k] = 0;

-			if (i == NX - 1 && j != 0 && j != NY - 1)

-				a[k] = 0.;

-			if (i == RX && j <= RY)

-				a[k] = a[k + NY];

-			if (i <= RX && j == RY)

-				a[k] = a[k + 1];

-		}

-	}

-}

-void korekh(double *a, double *b, double *c)

-{

-	int i, j, k;

-	for (i = 0; i < NX; i++) {

-		for (j = 0; j < NY; j++) {

-			k = i *NY + j;

-

-			if (i>0 && i < NX - 2 && ((i <= RX && j <= RY) != 1)){

-				a[k] -= DT*NX*(c[k + NY] - c[k]);

-			}

-

-			if (j>0 && j < NY - 2 && ((i <= RX && j <= RY) != 1)){

-				b[k] -= DT*NY*KT*(c[k + 1] - c[k]);

-			}

-

-		}

-	}

-}

-

-void poissh(double *a, double *b, double *c, double *d)

-{

-	int i, j, k;

-	double corr;

-	for (i = 0; i < NX; i++) {

-		for (j = 0; j < NY; j++) {

-			k = i *NY + j;

-			corr = c[k];

-			if (i > 0 && j > 0 && i < NX - 1 && j < NY - 1 && ((i <= RX && j <= RY) != 1)){

-				c[k] = -(1.0 / (2 * (DT*NX*NX + DT*NY*NY*KT*KT)))*((-DT*NX*NX)*c[k + NY] + (-DT*NX*NX)*c[k - NY] + (-DT*NY*NY*KT*KT)*c[k + 1] + (-DT*NY*NY*KT*KT)*c[k - 1] + (1.0*NX)*(a[k] - a[k - NY]) + (1.0*NY*KT)*(b[k] - b[k - 1]));

-			}

-			corr -= c[k];

-			if (corr < 0)

-				corr *= (-1);

-			d[k] = corr;

-		}

-	}

-}

-

-

-

-int main()

-{

-

-	clock_t t2, t1;

-	FILE *fp, *fp1, *fp2;

-

-	double *a_h, *a_h1, *b_h, *b_h1, *c_h, *c_h1, *corr_h;

-

-	double  dx1, dy1;

-	int j, i, t, k;

-

-	fp = fopen("x.dat", "w");

-	fp1 = fopen("y.dat", "w");

-	fp2 = fopen("tlak.dat", "w");

-	t1 = clock();

-

-	a_h = (double*)malloc(sizeof(double)*NX*NY);

-	a_h1 = (double*)malloc(sizeof(double)*NX*NY);

-	b_h = (double*)malloc(sizeof(double)*NX*NY);

-	b_h1 = (double*)malloc(sizeof(double)*NX*NY);

-	c_h = (double*)malloc(sizeof(double)*NX*NY);

-	c_h1 = (double*)malloc(sizeof(double)*NX*NY);

-	corr_h = (double*)malloc(sizeof(double)*NX*NY);

-

-	inixh(a_h);

-	inixh(a_h1);

-	inih(b_h);

-	inih(b_h1);

-	inih(c_h);

-	inih(c_h1);

-

-

-	printf("Start");

-	for (t = 0; t<T; t++)

-	{

-

-		USh(a_h, b_h, c_h, a_h1);

-		VSh(a_h, b_h, c_h, b_h1);

-		rubnixh(a_h1);

-		rubniyh(b_h1);

-

-		for (i = 0; i<30; i++)

-			poissh(a_h1, b_h1, c_h1, corr_h);

-

-		rubniph(c_h1);

-		korekh(a_h1, b_h1, c_h1);

-		rubnixh(a_h1);

-		rubniyh(b_h1);

-		sumah(c_h1, c_h);

-		inih(c_h1);

-		rubniph(c_h);

-		USh(a_h1, b_h, c_h, a_h);

-		VSh(a_h, b_h1, c_h, b_h);

-		rubnixh(a_h);

-		rubniyh(b_h);

-

-		for (i = 0; i<30; i++)

-	    poissh(a_h, b_h, c_h1, corr_h);

-

-

-		rubniph(c_h1);

-		korekh(a_h, b_h, c_h1);

-		rubnixh(a_h);

-		rubniyh(b_h);

-		sumah(c_h1, c_h);

-		inih(c_h1);

-		rubniph(c_h);

-

-

-	}

-	for (i = 0; i<NX; i++)

-		for (j = 0; j<NY; j++) if (i != NX - 1)

-			fprintf(fp, "\n%d %d %f", i, j, a_h[i*NY + j]);

-	for (i = 0; i<NX; i++)

-		for (j = 0; j<NY; j++) if (j != NY - 1)

-			fprintf(fp1, "\n%d %d %f", i, j, b_h[i*NY + j]);

-	for (i = 0; i<NX; i++)

-		for (j = 0; j<NY; j++)

-			fprintf(fp2, "\n%d %d %f", i, j, c_h[i*NY + j]);

-	t2 = clock();

-	printf(" %f\n", ((double)(t2 - t1)) / CLOCKS_PER_SEC);

-

-	dx1 = 1.0*NX - 1;

-	dy1 = KT*NY - 1;

-	for (i = 0; i<NX; i++)

-		dxx[i] = 1.0*i / dx1;

-	for (i = 0; i<NY; i++)

-		dyy[i] = 1.0*i / dy1;

-	for (i = 0; i<NX; i++)

-	{

-		xi[i][0] = 0.0;

-		yi[i][0] = 0.0;

-		xi[i][NY - 1] = 0.0;

-		yi[i][NY - 1] = 0.0;

-	}

-	for (i = 0; i<NY; i++)

-	{

-		xi[0][i] = 1.0;

-		yi[0][i] = 0.0;

-		xi[NX - 1][i] = 0.0;

-		yi[NX - 1][i] = 0.0;

-	}

-	for (i = 1; i <= NX - 2; i++)

-		for (j = 1; j <= NY - 1; j++)

-		{

-		k = (i - 1)*NY + j;

-

-		xi[i][j] = 0.5*(a_h[k] + a_h[k + NY]);

-		}

-	for (i = 1; i <= NX - 1; i++)

-		for (j = 1; j <= NY - 2; j++)

-		{

-		k = i*NY + j - 1;

-

-		yi[i][j] = 0.5*(b_h[k] + b_h[k + 1]);

-		}

-

-	PRINT_VECT_VTK(xi, yi, dxx, dyy);

-	PRINT_VTK(xi, yi, dxx, dyy);

-

-	for (j = 0; j < NX; j++)

-	for (i = 0; i < NY ; i++)

-		 {

-			dx[j][i] = (1.0*j / dx1);

-			dy[j][i]=  (1.0*i / dy1);

-		}

-

-

-

-	PRINT_PLT();

-

-

-	printf("END");

-

-

-    return 0;

-}

-

-

-