From 5c4926bca4e6479ae49521683abe0b0f4d3cfcda Mon Sep 17 00:00:00 2001 From: Ilhan Özgen Xian Date: Wed, 9 Jun 2021 18:56:14 -0700 Subject: add meshing functionality --- .gitignore | 4 +- Makefile | 27 +++++- header.h | 229 ++++++++++++++++++++++++++------------------ main.c | 319 +++++++++++++++++++++++++++++++++++++++++++------------------ 4 files changed, 386 insertions(+), 193 deletions(-) diff --git a/.gitignore b/.gitignore index ae1fdfa..5ce281c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,6 @@ output/ *~ mammoth -*.o \ No newline at end of file +*.o +.DS_STORE +TAGS \ No newline at end of file diff --git a/Makefile b/Makefile index 7174c85..6840954 100644 --- a/Makefile +++ b/Makefile @@ -1,11 +1,20 @@ # -# Makefile +# Makefile for mammoth # +# Author: +# Dr. Ilhan Ozgen-Xian +# Lawrence Berkeley National Laboratory 2021 +# +# Type "make" to compile mammoth +# + PRG=mammoth SRC=main.c +HDR=header.h -CC=/usr/local/bin/gcc-10 +CC=/usr/local/bin/gcc-11 +CTAGS=uctags TRIANGLE_DIR=/Users/IOzgen/Documents/work/coding/tools/triangle @@ -14,7 +23,7 @@ LFLAGS=-I$(TRIANGLE_DIR) CPU= RUNFLAGS= INCLUDES= -LIB= +LIB=$(TRIANGLE_DIR)/triangle.o # # @@ -23,6 +32,7 @@ LIB= OBJ=$(SRC:.c=.o) all: $(PRG) + make tags run: $(PRG) $(RUNFLAGS) ./$(PRG) @@ -34,9 +44,16 @@ $(PRG): $(OBJ) $(CC) $(LFLAGS) $(CFLAGS) -c $< -o $@ clean: - rm -vf $(PRG) $(OBJ) output/* + rm -vf $(PRG) $(OBJ) rm -vf $(P) -.PHONY: all clean +cleanall: + rm -vf output/* + make clean + +tags: $(SRC) + $(CTAGS) -Re $(SRC) $(HDR) + +.PHONY: all clean cleanall tags # eof diff --git a/header.h b/header.h index 3988abb..1391b5b 100644 --- a/header.h +++ b/header.h @@ -43,6 +43,27 @@ int multiresolution_meshing(int argc, char** argv); * ---------------------------------------------------------------------- */ +int sign(double x) +{ + int flag = (x > 0) ? 1 : -1; + if (x == 0) { + flag = 0; + } + return flag; +} + +double max(double a, double b) +{ + double results = (a > b) ? a : b; + return results; +} + +double min(double a, double b) +{ + double results = (a < b) ? a : b; + return results; +} + int ipow(int base, int exponent) { int result = 1; @@ -116,101 +137,6 @@ int transpose(double* matrix, int nrow, int ncol) { return 0; } -void report(io, markers, reporttriangles, reportneighbors, reportsegments, - reportedges, reportnorms) -struct triangulateio *io; -int markers; -int reporttriangles; -int reportneighbors; -int reportsegments; -int reportedges; -int reportnorms; -{ - int i, j; - - for (i = 0; i < io->numberofpoints; i++) { - printf("Point %4d:", i); - for (j = 0; j < 2; j++) { - printf(" %.6g", io->pointlist[i * 2 + j]); - } - if (io->numberofpointattributes > 0) { - printf(" attributes"); - } - for (j = 0; j < io->numberofpointattributes; j++) { - printf(" %.6g", - io->pointattributelist[i * io->numberofpointattributes + j]); - } - if (markers) { - printf(" marker %d\n", io->pointmarkerlist[i]); - } else { - printf("\n"); - } - } - printf("\n"); - - if (reporttriangles || reportneighbors) { - for (i = 0; i < io->numberoftriangles; i++) { - if (reporttriangles) { - printf("Triangle %4d points:", i); - for (j = 0; j < io->numberofcorners; j++) { - printf(" %4d", io->trianglelist[i * io->numberofcorners + j]); - } - if (io->numberoftriangleattributes > 0) { - printf(" attributes"); - } - for (j = 0; j < io->numberoftriangleattributes; j++) { - printf(" %.6g", io->triangleattributelist[i * - io->numberoftriangleattributes + j]); - } - printf("\n"); - } - if (reportneighbors) { - printf("Triangle %4d neighbors:", i); - for (j = 0; j < 3; j++) { - printf(" %4d", io->neighborlist[i * 3 + j]); - } - printf("\n"); - } - } - printf("\n"); - } - - if (reportsegments) { - for (i = 0; i < io->numberofsegments; i++) { - printf("Segment %4d points:", i); - for (j = 0; j < 2; j++) { - printf(" %4d", io->segmentlist[i * 2 + j]); - } - if (markers) { - printf(" marker %d\n", io->segmentmarkerlist[i]); - } else { - printf("\n"); - } - } - printf("\n"); - } - - if (reportedges) { - for (i = 0; i < io->numberofedges; i++) { - printf("Edge %4d points:", i); - for (j = 0; j < 2; j++) { - printf(" %4d", io->edgelist[i * 2 + j]); - } - if (reportnorms && (io->edgelist[i * 2 + 1] == -1)) { - for (j = 0; j < 2; j++) { - printf(" %.6g", io->normlist[i * 2 + j]); - } - } - if (markers) { - printf(" marker %d\n", io->edgemarkerlist[i]); - } else { - printf("\n"); - } - } - printf("\n"); - } -} - /* * ---------------------------------------------------------------------- * meshing functions @@ -300,3 +226,116 @@ int get_border(double *data, double nodata, int nrow, int ncol) return 0; } + +void write_avs(struct triangulateio *io, int counter) +{ + + char fname[1024]; + sprintf(fname, "output/mesh%d.inp", counter); + + FILE *fp; + fp = fopen(fname, "w"); + + fprintf(fp, "%d %d %d %d %d\n", io->numberofpoints, io->numberoftriangles, 0, 1, 0); + + for (int i = 0; i < io->numberofpoints; ++ i) { + fprintf(fp, "%d %f %f %f\n", i+1, io->pointlist[2*i], io->pointlist[2*i+1], 0.0); + } + + for (int i = 0; i < io->numberoftriangles; ++ i) { + fprintf(fp, "%d %d tri %d %d %d\n", i+1, 1, io->trianglelist[3*i], io->trianglelist[3*i+1], io->trianglelist[3*i+2]); + } + + fclose(fp); + +} + +void report(struct triangulateio *io, int markers, int reporttriangles, + int reportneighbors, int reportsegments, int reportedges, + int reportnorms) +{ + + int i, j; + + for (i = 0; i < io->numberofpoints; i++) { + printf("Point %4d:", i); + for (j = 0; j < 2; j++) { + printf(" %.6g", io->pointlist[i * 2 + j]); + } + if (io->numberofpointattributes > 0) { + printf(" attributes"); + } + for (j = 0; j < io->numberofpointattributes; j++) { + printf(" %.6g", + io->pointattributelist[i * io->numberofpointattributes + j]); + } + if (markers) { + printf(" marker %d\n", io->pointmarkerlist[i]); + } else { + printf("\n"); + } + } + printf("\n"); + + if (reporttriangles || reportneighbors) { + for (i = 0; i < io->numberoftriangles; i++) { + if (reporttriangles) { + printf("Triangle %4d points:", i); + for (j = 0; j < io->numberofcorners; j++) { + printf(" %4d", io->trianglelist[i * io->numberofcorners + j]); + } + if (io->numberoftriangleattributes > 0) { + printf(" attributes"); + } + for (j = 0; j < io->numberoftriangleattributes; j++) { + printf(" %.6g", io->triangleattributelist[i * + io->numberoftriangleattributes + j]); + } + printf("\n"); + } + if (reportneighbors) { + printf("Triangle %4d neighbors:", i); + for (j = 0; j < 3; j++) { + printf(" %4d", io->neighborlist[i * 3 + j]); + } + printf("\n"); + } + } + printf("\n"); + } + + if (reportsegments) { + for (i = 0; i < io->numberofsegments; i++) { + printf("Segment %4d points:", i); + for (j = 0; j < 2; j++) { + printf(" %4d", io->segmentlist[i * 2 + j]); + } + if (markers) { + printf(" marker %d\n", io->segmentmarkerlist[i]); + } else { + printf("\n"); + } + } + printf("\n"); + } + + if (reportedges) { + for (i = 0; i < io->numberofedges; i++) { + printf("Edge %4d points:", i); + for (j = 0; j < 2; j++) { + printf(" %4d", io->edgelist[i * 2 + j]); + } + if (reportnorms && (io->edgelist[i * 2 + 1] == -1)) { + for (j = 0; j < 2; j++) { + printf(" %.6g", io->normlist[i * 2 + j]); + } + } + if (markers) { + printf(" marker %d\n", io->edgemarkerlist[i]); + } else { + printf("\n"); + } + } + printf("\n"); + } +} diff --git a/main.c b/main.c index f59bdad..1a1a9d3 100644 --- a/main.c +++ b/main.c @@ -1,4 +1,4 @@ -/* -*- mode: c; style: linux -*- */ +/* -*- mode: c -*- */ /* Purpose: Given input data in the form of a 2D array, perform a * discrete Haar wavelet transform. @@ -10,7 +10,7 @@ /* Main logic */ -int main(int argc, char**argv) +int main(int argc, char** argv) { printf("\n\n"); @@ -21,23 +21,23 @@ int main(int argc, char**argv) printf(" | | | |/ ___ \\| | | | | | | |_| || | | _ |\n"); printf(" |_| |_/_/ \\_\\_| |_|_| |_|\\___/ |_| |_| |_|\n\n"); - printf(" +----------------------------+\n"); - printf(" | |\n"); - printf(" | . ........... |\n"); - printf(" | ................... |\n"); - printf(" | .... .. ........... .. |\n"); - printf(" | ... . .............. . |\n"); - printf(" | .. .............. |\n"); - printf(" | .. .. .. |\n"); - printf(" | .. .. |\n"); - printf(" | |\n"); - printf(" +----------------------------+\n\n"); - - printf(" I. Ozgen-Xian, G. Kesserwani, D. Caviedes-Voullieme,\n"); - printf(" S. Molins, Z. Xu, D. Dwivedi, J. D. Moulton, C. I. Steefel,\n"); - printf(" Wavelet-based local mesh refinement for rainfall-runoff\n"); - printf(" simulations, Journal of Hydroinformatics, 22, 1059-1077, 2020.\n"); - printf(" doi: 10.2166/hydro.2020.198\n\n"); + /* printf(" +----------------------------+\n"); */ + /* printf(" | |\n"); */ + /* printf(" | . ........... |\n"); */ + /* printf(" | ................... |\n"); */ + /* printf(" | .... .. ........... .. |\n"); */ + /* printf(" | ... . .............. . |\n"); */ + /* printf(" | .. .............. |\n"); */ + /* printf(" | .. .. .. |\n"); */ + /* printf(" | .. .. |\n"); */ + /* printf(" | |\n"); */ + /* printf(" +----------------------------+\n\n"); */ + + printf(" I. Ozgen-Xian, G. Kesserwani, D. Caviedes-Voullieme,\n"); + printf(" S. Molins, Z. Xu, D. Dwivedi, J. D. Moulton, C. I. Steefel,\n"); + printf(" Wavelet-based local mesh refinement for rainfall-runoff\n"); + printf(" simulations, Journal of Hydroinformatics, 22, 1059-1077, 2020.\n"); + printf(" doi: 10.2166/hydro.2020.198\n\n"); @@ -161,9 +161,9 @@ int multiresolution_meshing(int argc, char** argv) boundary_coords[3*counter] = index_x * cellsize + xllcorner; boundary_coords[3*counter+1] = (nrow - index_y) * cellsize + yllcorner; boundary_coords[3*counter+2] = z; - #if VERBOSE +#if VERBOSE printf("[ok] Read point %d/%d\n", counter + 1, n_bounds); - #endif +#endif } fclose(bp); @@ -193,6 +193,7 @@ int multiresolution_meshing(int argc, char** argv) levels = calloc(ncol * nrow, sizeof(int)); + for (int i = 0; i < nrow; i ++) { for (int j = 0; j < ncol; j ++) { levels[j + ncol * i] = (int) nodata; @@ -220,15 +221,19 @@ int multiresolution_meshing(int argc, char** argv) fclose(lp); -#if DEBUG_IO for (int i = 0; i < nrow; i ++) { for (int j = 0; j < ncol; j ++) { - levels[j + ncol * i] = _tmp_levels[j + ncol * i]; + levels[j + ncol * i] = _tmp_levels[j + ncol_padded * i]; +#if DEBUG_IO printf(" %d ", levels[j + ncol * i]); +#endif } +#if DEBUG_IO printf("\n"); - } #endif + } + + printf("[OK] Level data read\n"); /* * ---------------------------------------------------------------------- @@ -236,99 +241,172 @@ int multiresolution_meshing(int argc, char** argv) * ---------------------------------------------------------------------- */ - double a0 = cellsize * cellsize; - + double a0 = cellsize * cellsize * ipow(2, nlevels-1); + + printf("[OK] Building mesh with largest cell area %.2f m2 {\n\n", a0); + struct triangulateio in; - struct triangulateio mid; - struct triangulateio out; - struct triangulateio vorout; in.numberofpoints = n_bounds; in.numberofpointattributes = 1; - + in.numberofsegments = n_bounds; + in.numberofholes = 0; + in.numberofregions = 0; + in.pointlist = (double *) calloc(in.numberofpoints * 2, sizeof(double)); in.pointattributelist = (double *) calloc(in.numberofpoints * in.numberofpointattributes, sizeof(double)); in.pointmarkerlist = (int *) calloc(in.numberofpoints, sizeof(int)); + in.segmentlist = (int *) calloc(in.numberofsegments * 2, sizeof(int)); + in.segmentmarkerlist = (int *) calloc(in.numberofsegments, sizeof(int)); + + /* Each of the columns of this array contains the constraint's x and + y coordinates are at indices [1] and [2], followed by the + regional attribute at index [3], followed by the maximum area at + index [4]. So we have size(regionlist,1)==4. Note that each + regional attribute is used only if you select the 'A' switch, and + each area constraint is used only if you select the 'a' switch + (with no number following), but omitting one of these switches + does not change the memory layout. Input only, although the + pointer is copied to the output structure for convenience. */ + in.regionlist = (double *) calloc(in.numberofregions * 4, sizeof(double)); - in.numberofsegments = 0; - in.numberofholes = 0; - in.numberofregions = 1; - - for (int i = 0; i < n_bounds; i ++) { + for (int i = 0; i < in.numberofpoints; ++ i) { + in.pointmarkerlist[i] = 1; + //in.pointattributelist[i] = boundary_coords[3 * i + 2]; // set z coordinates + for (int j = 0; j < 2; ++ j) { + in.pointlist[2 * i + j] = boundary_coords[3 * i + j]; + } + } - in.pointlist[2*i] = boundary_coords[3*i]; - in.pointlist[2*i+1] = boundary_coords[3*i+1]; - in.pointattributelist[i] = 1.0; - in.pointmarkerlist[i] = 0; - + for (int i = 0; i < in.numberofsegments-1; ++ i) { + in.segmentlist[2*i] = in.numberofsegments-i; + in.segmentlist[2*i+1] = in.numberofsegments-1-i; } - in.regionlist = (double *) calloc(in.numberofregions * 4, sizeof(double)); - in.regionlist[3] = a0; + in.segmentlist[2*in.numberofpoints-2] = 1; + in.segmentlist[2*in.numberofpoints-1] = in.numberofsegments; - /* - * Make necessary initializations so that Triangle can return a - * triangulation in `mid' and a voronoi diagram in `vorout'. - */ + struct triangulateio out[nlevels]; + + out[0].pointlist = (double *) NULL; + out[0].pointattributelist = (double *) NULL; + out[0].pointmarkerlist = (int *) NULL; + out[0].trianglelist = (int *) NULL; + out[0].triangleattributelist = (double *) NULL; + out[0].neighborlist = (int *) NULL; + out[0].segmentlist = (int *) NULL; + out[0].segmentmarkerlist = (int *) NULL; - mid.pointlist = (double *) NULL; - mid.pointattributelist = (double *) NULL; - mid.pointmarkerlist = (int *) NULL; - mid.trianglelist = (int *) NULL; - mid.triangleattributelist = (double *) NULL; - mid.segmentlist = (int *) NULL; - mid.segmentmarkerlist = (int *) NULL; - mid.edgelist = (int *) NULL; - mid.edgemarkerlist = (int *) NULL; + printf("[ok] Initial triangulation..."); + + char triangleargs[1024]; + sprintf(triangleargs, "QDpa%f", a0); + triangulate(triangleargs, &in, &out[0], (struct triangulateio *) NULL); - vorout.pointlist = (double *) NULL; - vorout.pointattributelist = (double *) NULL; - vorout.edgelist = (int *) NULL; - vorout.normlist = (double *) NULL; + printf(" done\n"); - triangulate("pczevn", &in, &mid, &vorout); + write_avs(&out[0], 1); - printf("Initial triangulation:\n\n"); - report(&mid, 1, 1, 1, 1, 1, 0); - printf("Initial Voronoi diagram:\n\n"); - report(&vorout, 0, 0, 0, 0, 1, 1); + for (int current_level = 1; current_level < nlevels; ++ current_level) { - mid.trianglearealist = (double *) calloc(mid.numberoftriangles, sizeof(double)); + out[current_level-1].trianglearealist = (double *) calloc(out[current_level-1].numberoftriangles, sizeof(double)); - for (int i = 0; i < mid.numberoftriangles; i ++) { - mid.trianglearealist[i] = a0; - } + // tag the stuff + for (int i = 0; i < out[current_level-1].numberoftriangles; ++ i) { - out.pointlist = (double *) NULL; - out.pointattributelist = (double *) NULL; - out.trianglelist = (int *) NULL; - out.triangleattributelist = (double *) NULL; + int p1 = out[current_level-1].trianglelist[3*i] - 1; + int p2 = out[current_level-1].trianglelist[3*i+1] - 1; + int p3 = out[current_level-1].trianglelist[3*i+2] - 1; + + double x1 = out[current_level-1].pointlist[2*p1] - xllcorner; + double y1 = out[current_level-1].pointlist[2*p1+1] - yllcorner; + + double x2 = out[current_level-1].pointlist[2*p2] - xllcorner; + double y2 = out[current_level-1].pointlist[2*p2+1] - yllcorner; + + double x3 = out[current_level-1].pointlist[2*p3] - xllcorner; + double y3 = out[current_level-1].pointlist[2*p3+1] - yllcorner; + + double xmin = min(x1, min(x2, x3)); + double ymin = min(y1, min(y2, y3)); + double xmax = max(x1, max(x2, x3)); + double ymax = max(y1, max(y2, y3)); + + int ixmin = (int) floor(xmin / cellsize); + int iymin = (int) floor(ymin / cellsize); + int ixmax = (int) floor(xmax / cellsize); + int iymax = (int) floor(ymax / cellsize); + + /* double xc = (x1 + x2 + x3) / 3.0; */ + /* double yc = (y1 + y2 + y3) / 3.0; */ + /* double a = 0.5 * (- y2 * x3 + y1 * (-x2 + x3) + x1 * (y2 - y3) + x2 * y3); */ + /* double c = 1.0 / (2.0 * a); */ + + /* double f1 = y1 * x2; */ + /* double f2 = x1 * y3; */ + /* double f3 = y3 - y1; */ + /* double f4 = x1 - x3; */ + + /* double g1 = x1 * y2; */ + /* double g2 = y1 * x2; */ + /* double g3 = y1 - y2; */ + /* double g4 = x2 - x1; */ + + int flag = -1; + + /* printf("[ok] Triangle #%d/%d\n", i+1, mid.numberoftriangles); */ + /* printf("\t %d %d %d / %d\n", p1, p2, p3, mid.numberofpoints); */ + /* printf("\t %d %d : %d %d\n", ixmin, ixmax, iymin, iymax); */ + /* printf("\t %f %f : %f %f\n", xmin, xmax, ymin, ymax); */ + + for (int iy = iymin; iy < iymax; ++ iy) { + for (int ix = ixmin; ix < ixmax; ++ ix) { + /* printf("%d < %d?\n", flag, levels[ix + ncol * iy]); */ + flag = (flag < levels[ix + ncol * (nrow - iy)]) ? levels[ix + ncol * (nrow -iy)] : flag; + } + } + + /* printf("\t flag: %d...", flag); */ + + if (flag >= current_level) { + /* printf(" tagged\n"); */ + out[current_level-1].trianglearealist[i] = a0 / ipow(2, current_level); + } + + } + + out[current_level].pointlist = (double *) NULL; + out[current_level].pointattributelist = (double *) NULL; + out[current_level].pointmarkerlist = (int *) NULL; + out[current_level].trianglelist = (int *) NULL; + out[current_level].triangleattributelist = (double *) NULL; + out[current_level].segmentlist = (int *) NULL; + out[current_level].segmentmarkerlist = (int *) NULL; + + printf("[ok] Triangulation #%d ...", current_level + 1); + + triangulate("QDarp", &out[current_level-1], &out[current_level], (struct triangulateio *) NULL); + + printf(" done\n"); + + write_avs(&out[current_level], current_level+1); + + } + + printf("[OK] } done\n"); - triangulate("prazB", &mid, &out, (struct triangulateio *) NULL); free(in.pointlist); free(in.pointattributelist); free(in.pointmarkerlist); free(in.regionlist); - free(mid.pointlist); - free(mid.pointattributelist); - free(mid.pointmarkerlist); - free(mid.trianglelist); - free(mid.triangleattributelist); - free(mid.trianglearealist); - free(mid.neighborlist); - free(mid.segmentlist); - free(mid.segmentmarkerlist); - free(mid.edgelist); - free(mid.edgemarkerlist); - free(vorout.pointlist); - free(vorout.pointattributelist); - free(vorout.edgelist); - free(vorout.normlist); - free(out.pointlist); - free(out.pointattributelist); - free(out.trianglelist); - free(out.triangleattributelist); + + free(out[0].pointlist); + free(out[0].pointattributelist); + free(out[0].pointmarkerlist); + free(out[0].trianglelist); + free(out[0].triangleattributelist); + free(out[0].neighborlist); free(levels); free(boundary_coords); @@ -686,6 +764,63 @@ int multiresolution_analysis(int argc, char** argv) printf("done\n"); +/* FILE *bp; */ +/* FILE *op; */ + +/* bp = fopen("output/boundary.asc", "r"); */ + +/* int index_x; */ +/* int index_y; */ +/* double z; */ + +/* int n_bounds = 0; */ +/* char _tmp_char; */ + +/* while (!feof(bp)) { */ +/* _tmp_char = fgetc(bp); */ +/* if (_tmp_char == '\n') { */ +/* n_bounds ++; */ +/* } */ +/* } */ + +/* printf("[ok] Found %d points\n", n_bounds); */ + +/* rewind(bp); */ + +/* double *boundary_coords = calloc(3 * n_bounds, sizeof(double)); */ + +/* for (int counter = 0; counter < n_bounds; counter ++) { */ +/* fscanf(bp, "%d %d %lf", &index_x, &index_y, &z); */ +/* boundary_coords[3*counter] = index_x * cellsize + xllcorner; */ +/* boundary_coords[3*counter+1] = (nrow - index_y) * cellsize + yllcorner; */ +/* boundary_coords[3*counter+2] = z; */ +/* #if VERBOSE */ +/* printf("[ok] Read point %d/%d\n", counter + 1, n_bounds); */ +/* #endif */ +/* } */ + +/* fclose(bp); */ + +/* op = fopen("output/boundary.poly", "w"); */ + +/* fprintf(op, "%d %d %d %d\n", n_bounds, 2, 1, 0); */ + +/* for (int i = 0; i < n_bounds; ++ i) { */ +/* fprintf(op, "%d %f %f %f\n", i+1, boundary_coords[3*i], boundary_coords[3*i+1], boundary_coords[3*i+2]); */ +/* } */ + +/* fprintf(op, "%d %d\n", n_bounds, 0); */ + +/* for (int i = 0; i < n_bounds-1; ++ i) { */ +/* fprintf(op, "%d %d %d\n", i+1, n_bounds-i, n_bounds-1-i); */ +/* } */ + +/* fprintf(op, "%d %d %d\n0 0\n", n_bounds, 1, n_bounds); */ + +/* fclose(op); */ + +/* free(boundary_coords); */ + //FILE *dp; // //dp = fopen("output/boundary.asc", "r"); -- cgit v1.2.3-13-gbd6f