Skip to content
Snippets Groups Projects
Commit 123e6a0c authored by Cristian Weiland's avatar Cristian Weiland
Browse files

Structure u becomes a double * and add header.h


Signed-off-by: default avatarCristian <cw14@inf.ufpr.br>
parent 63f0b79b
Branches
No related tags found
1 merge request!1Develop
header.h 0 → 100644
#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);
#include <stdio.h> #include "header.h"
#include <stdlib.h>
#include <string.h>
#include <math.h>
/* Final version of simplified equation: /* 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 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² + (- ...@@ -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. //#define MAX_SIZE 100000000 // 100 MB.
//int inMemory = 0; //int inMemory = 0;
int Nx, Ny, MaxI;
double Hx, Hy, W; void print_errno(void) {
double UDivisor; printf("%s",strerror(errno));
double Pipi = M_PI * M_PI; // Pipi = Pi * Pi }
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? 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) )); 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; 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 += 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[i-1][j] - u[i+1][j]) / 2 * Hx + (u[i][j-1] - u[i][j+1]) / 2 * 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; res = res / UDivisor;
u[i][j] = res; u[ in(i,j) ] = res;
} }
double residueNorm(double *r) { double residueNorm(double *r) {
...@@ -43,20 +67,7 @@ int main(int argc, char *argv[]) { ...@@ -43,20 +67,7 @@ int main(int argc, char *argv[]) {
int i, j, k; int i, j, k;
FILE *fpExit; FILE *fpExit;
fpExit = fopen(argv[5],"w"); getParams(argc,argv,fpExit);
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]);
Nx = (1/Hx); Nx = (1/Hx);
Ny = (1/Hy); Ny = (1/Hy);
...@@ -76,37 +87,40 @@ int main(int argc, char *argv[]) { ...@@ -76,37 +87,40 @@ int main(int argc, char *argv[]) {
}*/ }*/
double sigma; 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++) { for(j=0; i<Ny; i++) {
u[i][j] = calloc(sizeof(double)); u[i][j] = calloc(sizeof(double));
f[i][j] = malloc(sizeof(double)); f[i][j] = malloc(sizeof(double));
f[i][j] = f(i*Hx,j*Hy); f[i][j] = f(i*Hx,j*Hy);
} }
} }*/
double *u = calloc(sizeof(double), Nx * Ny);
sigma = sinh(M_PI * M_PI); sigma = sinh(M_PI * M_PI);
for(i=0; i<Nx; i++) { // Calculate side extremities (top and bottom are always 0). 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[ in(i,0) ] = sin(2 * M_PI * (M_PI - (i * Hx))) * sigma;
u[i][Nx] = sin(2 * M_PI * (i * Hx)); u[ in(i,Nx) ] = sin(2 * M_PI * (i * Hx));
} }
// Nx columns and Ny rows. // Nx columns and Ny rows.
for(k=0; k<MaxI; k++) { /*for(k=0; k<MaxI; k++) {
for(i=1; i<Ny; i++) { for(i=1; i<Ny; i++) {
sigma = 0; sigma = 0;
for(j=1; j<Nx; j++) { for(j=1; j<Nx; j++) {
/*if(j != i) { /*if(j != i) {
sigma += a[i][j] * fi[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]; sigma -= a[i][i] * fi[i];
//fi[i] = (1 - W) * fi[i] + W/a[i][i] * (b[i] - sigma); //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]) fi[i] = fi[i] + W * ((b[i] - sigma) / a[i][i] - fi[i])
} }
} }*/
return 0; return 0;
} }
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment