/* Last edited on 2006-05-24 21:33:17 by anamaria */ /* Functions to compute Ex, Ey, and Hz with UMPL conditions */ /* We need to define _GNU_SOURCE because we use the constant NAN. */ #define _GNU_SOURCE #include #include #include #include /* #include */ #include "GradFunc.h" grade_t *cria_grade(int N, int L, double c, double **fundo) { grade_t *M = malloc(sizeof(grade_t)); M->N = N; M->L = L; M->Xinf = matrix_alloc(L, N+1); M->Xsup = matrix_alloc(L, N+1); M->Yinf = matrix_alloc(N + 1 - 2*L, L); M->Ysup = matrix_alloc(N + 1 - 2*L, L); M->c = c; M->fundo = fundo; return M; } double mget(grade_t *M, int ix, int jy) { int N = M->N, L = M->L; if (ix < L) { return M->Xinf[ix][jy]; } if (ix > N-L) { return M->Xsup[ix-(N-L+1)][jy]; } if (jy < L) { return M->Yinf[ix - L][jy]; } if (jy > N-L) { return M->Ysup[ix - L][jy - (N-L+1)]; } return M->fundo[ix][jy] / M->c; } void mset(grade_t *M, int ix, int jy, double valor) { int N = M->N, L = M->L; if (ix < L) { M->Xinf[ix][jy] = valor; return; } if (ix > N-L) { M->Xsup[ix-(N-L+1)][jy] = valor; return; } if (jy < L) { M->Yinf[ix - L][jy] = valor; return; } if (jy > N-L) { M->Ysup[ix - L][jy - (N-L+1)] = valor; return; } fprintf(stderr, "tentativa de alterar a matriz de fundo\n"); exit(1); } double** matrix_alloc(int nl, int nc) { int i; double** dv = (double**) malloc(nl*sizeof(double*)); if (dv == NULL) { fprintf(stderr, "ERROR: allocation failure in dvector(%d, %d)\n", nl, nc); exit(-1); } for( i=0; i < nl; i++) { dv[i] = (double*)malloc(nc*sizeof(double)); if (dv[i] == NULL) { fprintf(stderr, "ERROR: allocation failure in dvector(%d, %d)\n", nl, nc); exit(-1); } } return dv; } /*==============================================================================================*/ double DerivAntSimet(grade_t *M, double dxy, int ix, int jy, int flag) { int N = M->N; int P = 4; /* Number of interpolation points. */ int sk[P]; /* Indices of the samples needed to compute the derivative */ double h[P]; /* Weights needed to compute the derivative */ double v[P]; /* Function values needed to compute the derivative */ int i, k; double du; if (flag==1) k = ix; else k = jy; if (k==0) { sk[0] = 2; sk[1] = 1; sk[2] = 1; sk[3] = 2; h[0] = -1./12.; h[1] = 2./3.; h[2] = 2./3.; h[3] = -1./12.; } else if (k==1) { sk[0] = 1; sk[1] = 0; sk[2] = 2; sk[3] = 3; h[0] = -1./12.; h[1] =-2./3.; h[2] = 2./3.; h[3] = -1./12.; } else if (k==N-1) { sk[0] = N-3; sk[1] = N-2; sk[2] = N; sk[3] = N-1; h[0] = 1./12.; h[1] =-2./3.; h[2] = 2./3.; h[3] = 1./12.; } else if (k==N) { sk[0] = N-2; sk[1] = N-1; sk[2] = N-1; sk[3] = N-2; h[0] = 1./12.; h[1] =-2./3.; h[2] =-2./3.; h[3] = 1./12.; } else { sk[0] = k-2; sk[1] = k-1; sk[2] = k+1; sk[3] = k+2; h[0] = 1./12.; h[1] =-2./3.; h[2] = 2./3.; h[3] =-1./12.; } for (i=0; i N) { fprintf(stderr, "ERROR: deriv(): sjy[%d]=%d\n", i, sk[i]); printf("Estou na derivada Anti-Simétrica \n"); exit(-1); } if (flag==1) v[i] = mget(M, sk[i], jy); else v[i] = mget(M, ix, sk[i]); } /* compute the derivative */ du = 0; for (i=0; i < P; i++) du += v[i] * h[i]; du /= (1*dxy); return du; } /*===========================================================================================*/ double DerivSimet(grade_t *M, double dxy, int ix, int jy, int flag) { int N = M->N; int P = 4; /* Number of interpolation points. */ int sk[P]; /* the positions of the samples needed to compute the derivative */ double h[P]; /* the weights needed to compute the derivative */ double v[P]; /* the function values needed to compute the derivative */ int i, k; double du; if (flag==1) k = ix; else k = jy; if (k==0) { sk[0] = 2; sk[1] = 1; sk[2] = 1; sk[3] = 2; h[0] = 1./12.; h[1] = -2./3.; h[2] = 2./3.; h[3] = -1./12.; } else if (k==1) { sk[0] = 1; sk[1] = 0; sk[2] = 2; sk[3] = 3; h[0] = 1./12.; h[1] =-2./3.; h[2] = 2./3.; h[3] = -1./12.; } else if (k==N-1) { sk[0] = N-3; sk[1] = N-2; sk[2] = N; sk[3] = N-1; h[0] = 1./12.; h[1] =-2./3.; h[2] = 2./3.; h[3] = -1./12.; } else if (k==N) { sk[0] = N-2; sk[1] = N-1; sk[2] = N-1; sk[3] = N-2; h[0] = 1./12.; h[1] =-2./3.; h[2] = 2./3.; h[3] = -1./12.; } else { sk[0] = k-2; sk[1] = k-1; sk[2] = k+1; sk[3] = k+2; h[0] = 1./12.; h[1] =-2./3.; h[2] = 2./3.; h[3] =-1./12.; } for (i=0; i N) { fprintf(stderr, "ERROR: deriv(): six[%d]=%d\n", i, sk[i]); printf("Estou na derivada de uma função simétrica \n"); exit(-1); } if (flag==1) v[i] = mget(M, sk[i], jy); else v[i] = mget(M, ix, sk[i]); } /* compute the derivative */ du = 0; for (i=0; i < P; i++) du += v[i] * h[i]; du /= (1*dxy); return du; } FILE* abre_grafico(int n, double dx, double dy) { FILE* pp = popen("gnuplot >& /dev/null", "w"); if (pp == NULL) { fatal("cannot open gnuplot"); } fputs("set data style lines\n", pp); fputs("set grid\n", pp); fprintf(pp, "set xrange [%g:%g]\n", -0.5*dx, (n+0.5)*dx); fprintf(pp, "set yrange [%g:%g]\n", -0.5*dy, (n+0.5)*dy); fprintf(pp, "set nokey\n"); return pp; } void mostra_matriz ( FILE *pp, char *tit, double time, int n, double dx, double dy, double dz, double **val, int espera ) { char buf[10]; int i, j; fprintf(pp, "set zrange [%g:%g]\n", -1.05*dz, +1.05*dz); fprintf(pp, "set title '%s time = %g ns'\n", tit, time); fputs("splot '-' with lines\n", pp); for (i=0; i <=n; i++) { for (j=0; j <=n; j++) fprintf(pp, "%g %g %g\n", i*dx, j*dy, val[i][j]); fprintf(pp, "\n"); } fputs("e\n", pp); fflush(pp); if (espera) { fputs("Hit return to continue", stderr); fgets(buf, sizeof(buf), stdin); } } void mostra_grade ( FILE *pp, char *tit, double time, int n, double dx, double dy, double dz, grade_t *G, double **val, int espera ) { int i, j; for (i=0; i <=n; i++) for (j=0; j <=n; j++) { val[i][j]=mget(G, i, j); } mostra_matriz(pp, tit, time, n, dx, dy, dz, val, espera); } void fatal(char* s) { fprintf(stderr, "ERROR: %s\n", s); perror(s); exit(-1); } /*===========================================================================================*/