Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <stdio.h>
- #include <stdlib.h>
- #include <math.h>
- #define PARAMETER1 100
- #define nx1 5
- #define nx2 10
- #define nx3 10
- #define ny1 4
- #define ny2 5
- #define ny3 10
- #define Nx nx1+nx2+nx3+1
- #define Ny ny1+ny2+ny3+1
- #define X1_LENGTH 5.0
- #define X2_LENGTH 10.0
- #define X3_LENGTH 50.0
- #define Y1_LENGTH 0.04
- #define Y2_LENGTH 5.0
- #define Y3_LENGTH 20.0
- double f(int j,int i);
- double g(int j,int i);
- double dx(int i);
- double calc_value_pos(int i);
- double calc_value_neg(int i);
- double **matrix;
- double f1(int j,int i);
- double calc_value2_pos(int i);
- double calc_value2_neg(int i);
- double f2(int j,int i);
- int main(void){
- matrix=(double **)malloc((size_t) ((Ny+1)*(Nx+1)+1)*sizeof(double * ));
- for (int j=1; j<=(Ny+1)*(Nx+1); j++) {
- matrix[j] = (double *)malloc((size_t) ((Ny+1)*(Nx+1)+1)* sizeof(double ));
- }
- for (int j=1; j<=Ny+1; j++){
- for (int i=1; i<=Nx+1; i++){
- matrix[j][i]=0.0;
- }
- }
- for (int i=0; i<=Nx; i++){
- for (int j=0; j<=Ny; j++){
- int k=i+(j)*(Nx+1)+1;
- matrix[k][k] = f(j,i); /*including this line causes the problem*/
- //matrix[k][k] = f1(j,i); /* including this line does not cause the problem*/
- }
- }
- return 1;
- }
- double f(int j, int i){
- double output=f1(j,i)+f2(j,i);
- return output;
- }
- double f1(int j, int i){
- double output;
- if (i==0 || j==0 || i==Nx || j==Ny) output=0.0;
- if (i!=0 && j!= 0 && i!= Nx && j!=Ny) {
- output = (g(j,i)*g(j,i+1))/(g(j,i)*calc_value_pos(i)+g(j,i+1)*calc_value_neg(i));
- double output2;
- output2=1.0/( (calc_value_neg(i)/g(j,i)) + (calc_value_pos(i)/g(j,i+1)));
- /*if (output2!=output2) printf("output2: %lf\n",output2);*/
- if (output != output2 && output==output) printf("(j%d,i%d) output:%lf \toutput2:%lf\n",j,i,output,output2);
- }
- if (output!=output) output=0.0;
- return output;
- }
- double f2(int j,int i){
- /* calculates f2 for (j,i) */
- double output;
- if (j==0 || i==0) output=0.0;
- else output = g(j,i)*g(j,i-1)/(g(j,i)*calc_value2_neg(i)+g(j,i-1)*calc_value2_pos(i));
- if (output!=output) output=0.0;
- return output;
- }
- double calc_value2_pos(int i){
- /* calculates calc_value2_pos for X=i */
- double output= dx(i)/2.0;
- return output;
- }
- double calc_value2_neg(int i){
- /* calculates calc_value2_neg for X=i */
- if (i==0) printf("Warning off grid\n");
- double output= dx(i-1)/2.0;
- return output;
- }
- double g(int j,int i){
- double output=PARAMETER1;
- if (j==0 || i==0) output=0.0;
- if (j>=ny1+1 && j<=ny1+ny2 && i>=1 && i<=nx1) output=0.0;
- if (output!=output) printf("calc_D(%d,%d) error: output==nan\n",j,i);
- return output;
- }
- double calc_value_pos(int i){
- if (i==Nx) printf("Warning east pos\n");
- double output;
- output = dx(i+1)/2.0;
- if (output!=output) printf("calc_delta_e_pos(%d) error: output==nan\n",i);
- return output;
- }
- double calc_value_neg(int i){
- double output=dx(i)/2.0;
- if (output!=output) printf("calc_delta_e_neg(%d) error: output==nan\n",i);
- return output;
- }
- double dx(int i){
- double output;
- if (i<=nx1) output=1.0*X1_LENGTH/nx1;
- if (i<=nx1+nx2 && i>=nx1+1) output=1.0*X2_LENGTH/nx2;
- if (i<=nx1+nx2+nx3 && i>=nx1+nx2+1) output=1.0*X3_LENGTH/nx3;
- if (output!=output) printf("delta_x(%d) error: output==nan\n",i);
- return output;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement