/* Last edited on 2023-02-25 16:15:46 by stolfi */ /* * spr2D.c * 29/11/2006 * * Demonstrates sparse decompositions based on iterated * interpolation in 2D Gaussian * Available geometries: empty domain, closed box, open horn. */ /* We need to define _GNU_SOURCE to get {asprintf}. */ #define _GNU_SOURCE #include #include #include #include #include #include #include #include #include /* The maximum number of samples that we can have: */ #define NPx 100 #define NPy 100 #define SIGMA_MAX (1e7) /* Nominal conductivity {sigma} of metals. */ /* NON-HOMOGENEOUS MEDIUM Non-homogeneous media are simulated by setting the appropriate elements of the {sigma} field. */ void setup_initial_state ( int geom, int nx, int ny, double dx, double dy, double **sigma, double **hz, double **ex, double **ey ); /* Sets the conductivity map {sigma} according to the code {geom} that describes the domain's geometry. Currently: 0 -> homogeneous non-conductive medium (vacuum). 1 -> closed square conductive box. 2 -> conductive microwave horn. Also sets the initial state of the magnetic field {hz} and of the electric field {ex,ey} to values appropriate to that geometry. */ void setup_geometry_box(int nx, int ny, double **sigma); /* Draws into {sigma[0..nx][0..ny]} the outline of a microwave horn, with open mouth and highly conductive sides. */ void setup_geometry_horn(int nx, int ny, double **sigma); /* Draws into {sigma[0..nx][0..ny]} a closed square box, with empty interior and highly conductive walls. */ /* DATA PLOTTING */ void show_field ( mwplot_t *wp, int nx, int ny, double dx, double dy, double **A, double vMag, char *zLabel, double time ); /* generates on {wp} a 3D graph of the field {A}, assumed to be an array of samples with {nx} rows and {ny} columns, with step sizes {dx} and {dy}, respectively. The vertical scale is set assuming the values of {A} are in {[-vMag _ +vMag]}. The Z axis will be labeled with {zLabel}, and the plot will have a title showing the specified {time}, assumed to be in nanoseconds. If the plotter {wp} is configured to write to disk instead of the screen, the plot file will be called "{zLabel}-{NNNNNN}.{ext}", where {NNNNNN} is the {time} value rounded to integer, and {ext} is the appropriate extension ("eps", "png", etc.) */ void show_grid ( mwplot_t *wp, int nx, int ny, double dx, double dy, double **A, double vMag, char *zLabel, double time ); /* Generates on {wp} a 2D plot of the matrix {A}, assumed to be an array of samples with {nx} rows and {ny} columns, with step sizes {dx} and {dy}, respectively. The color scale is set assuming the values of {A} are in {[-vMag _ +vMag]}. The plot will be a grid of colore squares. It will have a title showing the specified {time}, assumed to be in nanoseconds. If the plotter {wp} is configured to write to disk instead of the screen, the plot file will be called "{zLabel}-{NNNNNN}.{ext}", where {NNNNNN} is the {time} value rounded to integer, and {ext} is the appropriate extension ("eps", "png", etc.) */ void get_grid_structure(int nx, int ny, double **F, double **A); /* Sets {A[0..nx][0..ny]} to a map showing the occupied elements of the sparse arrray {F[0..nx][0..ny]}. Specifically. each {A[i][j]} is set to 0.0 if {F[i][j]} is NAN, and to 1 otherwise. The array {A} can then be displayed with {show_grid}. */ void get_field_range(int nx, int ny, double **F, double *vLo, double *vHi, double *vMag); /* Stores in {*vLo} and {*vHi}, respectively, the minimum and maximum values seen in {F[0..nx][0..ny]} Also stores in {*vHi} the maximum absolute value among all elements -- that is, {max{fabs(vHi),fabs(vLo)}}. */ /* UTILITY ROUTINES */ void fatal(char* s); double** dvector(int nl, int nc); /* Allocates an array with {nl} rows and {nc} columns. */ void wait_for_user(void); /* Writes "Hit return to continue" and waits for user to do so. */ double get_double_from_user(char *prompt, double defval); int get_int_from_user(char *prompt, int defval); /* Writes a message "Enter {promt} [{defval}]:" to {stderr}; then waits for the user to type a number ({double} or {int}, rspectively) + 'ENTER', and returns that number. Returns {defval} if the user just types 'ENTER'. */ void print_grid_statistics(int nx, int ny, double **A); /* Counts the number of valid (not NaN) entries in {A[0..nx][0..ny]}, compares with total number of entries {(nx+1)*(ny+1)}. */ /* MAIN PROGRAM */ int main() { /* Fundamental constants */ double cc = 2.99792458e8; /* speed of light in free space */ double muz = 4.0 * M_PI * 1.0e-7; /* permeability of free space */ double epsz =1.0 / (cc*cc*muz); /* permittivity of free space */ /* double aimp = sqrt(muz/epsz); */ /* wave impedance in free space */ int i, j, k; /* Starts the gnuplot process for 3D field plots: */ mwplot_t *wpf = mwplot_start_gnuplot(640,480); mwplot_set_terminal(wpf, /*x11*/ 0); /* Starts the gnuplot process for 2D grid plots */ mwplot_t *wpg = mwplot_start_gnuplot(640,480); /* 0,0 for debug output */ mwplot_set_terminal(wpg, /*x11*/ 0); mwplot_set_gridmap_palette(wpg, /*yellow--red*/ 0); fprintf(stderr, "Max samples in coarsest level = %d\n", NPx); int n0 = get_int_from_user("the number of samples in the coarsest level", 8); int ell = get_int_from_user("the number of decompositions levels", 4); int ell_atual = ell; /* actual number of decomp levels with significant wavelet coeffs */ int level_max; /* number of decomposition levels */ int interp_points; /* 2 or 4 point interpolation (linear or cubic) */ double eps = get_double_from_user("the error/threshold for sparse representation", 1e-6); double eps1 = eps/10; /* error/threshold for the sparse representation. */ interp_points = get_int_from_user("the number of interpolation points", 4); int n = n0 * (1< (*vHi)) { (*vHi) = Fij; } } } double aLo = fabs(*vLo), aHi = fabs(*vHi); (*vMag) = (aLo > aHi ? aLo : aHi); } void get_grid_structure(int nx, int ny, double **F, double **A) { /* Copy {hz} to {az}, changing NaNs to zeros: */ int i, j; for (i = 0; i <= nx; i++) for (j = 0; j <= ny; j++) { A[i][j] = (isnan(F[i][j]) ? 0.0 : 1.0); } /* DEBUG -- add some {-1} marks to debug plot ranges: */ auto void mark(int im, int jm); void mark(int im, int jm) { if (A[im][jm] == 0) { A[im][jm] = -1; } } for (i = 0; i < 5; i++) { mark(i,0); mark(i,i); mark(0,i); mark(nx-i,ny); mark(nx-i,ny-i); mark(nx,ny-i); } } void print_grid_statistics(int nx, int ny, double **A) { int i, j; int nvalid = 0; for (i = 0; i <= nx; i++) for (j = 0; j <= ny; j++) { if (!isnan(A[i][j])) { nvalid = nvalid+1; } } fprintf(stderr, "sparse grid: %7d points\n", nvalid); int ntotal=(nx+1)*(ny+1); fprintf(stderr, "full grid: %7d points\n", ntotal); fprintf(stderr, "fill ratio: %7.5f\n", ((double)nvalid)/((double)ntotal)); fprintf(stderr, "\n"); } void wait_for_user(void) { char buf[10]; fputs("\nHit return to continue", stderr); fgets(buf, sizeof(buf), stdin); } #define BUFSZ 1000 double get_double_from_user(char *prompt, double defval) { char buf[BUFSZ]; fprintf(stderr,"Enter %s", prompt); if (! isnan(defval)) { fprintf(stderr," [%g]", defval); } fprintf(stderr,": "); fgets(buf, BUFSZ, stdin); char *data = buf; while (((*data) == ' ') || ((*data) == '\011') || ((*data) == '\015')) { data++; } if (((*data) != '\n') && ((*data) != 0)) { /* Parse input as integer: */ char *rest = NULL; double dval = strtod(data, &rest); /* Check number syntax: */ while(((*rest) == ' ') || ((*rest) == '\011') || ((*rest) == '\015')) { rest++; } if ((*rest) != '\n') { fprintf(stderr, "** invalid character '%c' = '\\%03o'\n", (*rest), (*rest)); exit(-1); } return dval; } else { /* Empty input: */ return defval; } } int get_int_from_user(char *prompt, int defval) { /* All 32-bit ints are exactly representable as {double}s, so: */ double dval = get_double_from_user(prompt, (double)defval); /* Check whether value is an exact integer: */ int ival = (int)dval; if (dval != ((double)ival)) { fprintf(stderr, "** value '%g' must be an integer\n", dval); exit(-1); } return ival; } void fatal(char* s) { fprintf(stderr, "ERROR: %s\n", s); perror(s); exit(-1); } double** dvector(int nl, int nc) { int i; double** dv = (double**) malloc(nc*sizeof(double*)); if (dv == NULL) { fprintf(stderr, "ERROR: allocation failure in dvector(%d, %d)\n", nl, nc); exit(-1); } for( i=0; i