From db8034ca00aad56b98a71c1d3b19e200fe3120a0 Mon Sep 17 00:00:00 2001 From: eug-vs Date: Mon, 23 May 2022 16:37:42 +0400 Subject: feat!: remove cpp code --- src/source.cpp | 501 --------------------------------------------------------- 1 file changed, 501 deletions(-) delete mode 100755 src/source.cpp 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 -#include -#include -#include -#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 (iRY) - 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