diff --git a/header.h b/header.h
new file mode 100644
index 0000000000000000000000000000000000000000..ce8bc5a66a5634c3897d9f6ae5ef13d402779762
--- /dev/null
+++ b/header.h
@@ -0,0 +1,18 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <errno.h>
+
+#define ARGS_NUM 9
+
+int Nx, Ny, MaxI;
+double Hx, Hy, W;
+double UDivisor;
+double Pipi = M_PI * M_PI; // Pipi = Pi * Pi
+
+void print_errno(void);
+void getParams(int argc, char* argv[], FILE *fp);
+double f(double x, double y);
+void u(int i, int j, double **u, double **f);
+double residueNorm(double *r);
diff --git a/pdeSolver.c b/pdeSolver.c
index ac9382c7752693cf620bb175c056e35236d4f75f..afd011be23a9e5616b07599227fd9bcf4630377f 100644
--- a/pdeSolver.c
+++ b/pdeSolver.c
@@ -1,7 +1,4 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
+#include "header.h"
/* Final version of simplified equation:
u(i,j) = f(x,y) + (u(i+1,j) + u(i-1,j))/Δx² + (u(i,j+1) + u(i,j-1))/Δy² + (-u(i+1,j)+u(i-1,j))/2Δx + (-u(i,j+1)+u(i,j-1))/2Δy
@@ -12,21 +9,48 @@ u(i,j) = f(x,y) + (u(i+1,j) + u(i-1,j))/Δx² + (u(i,j+1) + u(i,j-1))/Δy² + (-
//#define MAX_SIZE 100000000 // 100 MB.
//int inMemory = 0;
-int Nx, Ny, MaxI;
-double Hx, Hy, W;
-double UDivisor;
-double Pipi = M_PI * M_PI; // Pipi = Pi * Pi
+
+void print_errno(void) {
+ printf("%s",strerror(errno));
+}
+
+void getParams(int argc, char* argv[], FILE *fp) {
+ if(argc != ARGS_NUM) {
+ fprintf(stderr,"Wrong number of arguments.\n");
+ exit(-1);
+ }
+ int i;
+ for(i=1; i<ARGS_NUM; i+=2) {
+ if(strcmp(argv[i],"-hx") == 0) {
+ Hx = atof(argv[i+1]);
+ } else if(strcmp(argv[i],"-hy") == 0) {
+ Hy = atof(argv[i+1]);
+ } else if(strcmp(argv[i],"-i") == 0) {
+ MaxI = atoi(argv[i+1]);
+ } else if(strcmp(argv[i],"-o") == 0) {
+ fp = fopen(argv[i+1],"w");
+ } else {
+ fprintf(stderr,"Incorrect parameter.\n");
+ exit(-1);
+ }
+ }
+}
double f(double x, double y) { // If stored in memory, it will need Nx * Ny * 8 bytes of memory. Should we do it?
return 4 * M_PI * ( ( sin(2 * M_PI * x) ) * ( sinh(M_PI * y) ) + ( sin(2 * Pipi - M_PI * x) ) * ( sinh(Pipi - M_PI * y) ));
}
-void u(int i, int j, double **u, double **f) {
+inline int in(int i, int j) {
+// Calcula o indice do vetor, como se fosse uma matriz.
+ return i*Ny + j;
+}
+
+void calcU(int i, int j, double *u) {
double res = 0;
- res += f[i][j] + (u[i+1][j] + u[i-1][j]) / Hx * Hx + (u[i][j+1] + u[i][j-1]) / Hy * Hy;
- res += (u[i-1][j] - u[i+1][j]) / 2 * Hx + (u[i][j-1] - u[i][j+1]) / 2 * Hy;
+ res += f(i,j) + (u[ in(i+1,j) ] + u[ in(i-1,j) ] ) / Hx * Hx + (u[ in(i,j+1) ] + u[ in(i,j-1) ]) / Hy * Hy;
+ res += (u[ in(i-1,j) ] - u[ in(i+1,j) ]) / 2 * Hx + (u[ in(i,j-1) ] - u[ in(i,j+1) ]) / 2 * Hy;
res = res / UDivisor;
- u[i][j] = res;
+ u[ in(i,j) ] = res;
}
double residueNorm(double *r) {
@@ -43,20 +67,7 @@ int main(int argc, char *argv[]) {
int i, j, k;
FILE *fpExit;
- fpExit = fopen(argv[5],"w");
-
- if(argc == 9) {
- if((strcmp(argv[1],"-hx") != 0) || (strcmp(argv[3],"-hy") != 0) || (strcmp(argv[5],"-i") != 0) || (strcmp(argv[7],"-o") != 0)) {
- puts("Please use the following format: './pdeSolver -hx <Hx> -hy <Hy> -i <maxIter> -o arquivo_saida.'");
- exit(-1);
- }
- } else {
- puts("Please use the following format: './pdeSolver -hx <Hx> -hy <Hy> -i <maxIter> -o arquivo_saida.'");
- }
-
- Hx = atof(argv[2]);
- Hy = atof(argv[4]);
- MaxI = atoi(argv[6]);
+ getParams(argc,argv,fpExit);
Nx = (1/Hx);
Ny = (1/Hy);
@@ -76,37 +87,40 @@ int main(int argc, char *argv[]) {
}*/
double sigma;
- double **u,**f;
+ //double **u,**f;
- for(i=0; i<Nx; i++) {
+ /*for(i=0; i<Nx; i++) {
for(j=0; i<Ny; i++) {
u[i][j] = calloc(sizeof(double));
f[i][j] = malloc(sizeof(double));
f[i][j] = f(i*Hx,j*Hy);
}
- }
+ }*/
+
+ double *u = calloc(sizeof(double), Nx * Ny);
+
sigma = sinh(M_PI * M_PI);
for(i=0; i<Nx; i++) { // Calculate side extremities (top and bottom are always 0).
- u[i][0] = sin(2 * M_PI * (M_PI - (i * Hx))) * sigma;
- u[i][Nx] = sin(2 * M_PI * (i * Hx));
+ u[ in(i,0) ] = sin(2 * M_PI * (M_PI - (i * Hx))) * sigma;
+ u[ in(i,Nx) ] = sin(2 * M_PI * (i * Hx));
}
// Nx columns and Ny rows.
- for(k=0; k<MaxI; k++) {
+ /*for(k=0; k<MaxI; k++) {
for(i=1; i<Ny; i++) {
sigma = 0;
for(j=1; j<Nx; j++) {
/*if(j != i) {
sigma += a[i][j] * fi[i];
}*/
- sigma += a[i][j] * fi[i];
+ /*sigma += a[i][j] * fi[i];
}
sigma -= a[i][i] * fi[i];
//fi[i] = (1 - W) * fi[i] + W/a[i][i] * (b[i] - sigma);
fi[i] = fi[i] + W * ((b[i] - sigma) / a[i][i] - fi[i])
}
- }
+ }*/
return 0;
}
\ No newline at end of file