Convert ADCIRC ascii grid files and accompanying simulation output into F5. The levee information contained in the grid file is put into a specific triangular mesh for visualization purposes.
#include <stdio.h> #include <string.h> #include <stdlib.h> #include <assert.h> #include <hdf5.h> #include <F5/F5X.h> #include <F5/F5F.h> #include <F5/F5B.h> #include <F5/F5R.h> #include <F5/F5coordinates.h> #include <F5/F5surface.h> #include <F5/F5defs.h> #include <vector> #include <string> #include <iostream> using namespace std; /* TODO: Should really implement command line argument parsing, instead of reading hard-coded file names in the current directory. Data format can be found here: http://www.adcirc.org/documentv46/Output_file_descriptions.html .14 and .63 are interesting... .14 is given on cmd line it contains vertices and connectivity and more for the types see fort.14 description. Added number of time steps to be read as second a command line argument. */ void AppendSurfaceData( hid_t fileID, const string&gridname, const vector<F5_vec3_point_t>&Coords, const vector<F5_triangle_t>&Triangles, int time_steps_to_load) { char buf[2048]; double time = 0; //FILE*kat63 = fopen("fortkat22.63", "rt"), // *kat74 = fopen("fortkat22.74", "rt"); FILE*kat63 = fopen("fort.63", "rt"), *kat64 = fopen("fort.64", "rt"), *kat74 = fopen("fort.74", "rt"); FILE*elev = 0; #ifdef R_OK if (access("fort.63.gz", R_OK) == 0) { elev = popen("zcat fort.63.gz", "r"); if (elev==0) perror("Error in reading fort.63"); else fputs("Fine, could open fort63 file.\n", stderr); } #endif if (!(kat63 || kat74 || kat64|| elev)) { fputs("No additional field information available.\n", stderr); return; } unsigned int Timesteps = 0, Nodes = 0; // Skip the first 3 lines, ctj if (kat63) { fgets(buf, sizeof(buf), kat63); fprintf(stderr, buf ); fgets(buf, sizeof(buf), kat63); sscanf(buf, "%d %d", &Timesteps, &Nodes); fgets(buf, sizeof(buf), kat63); fprintf(stderr, buf ); } if (kat64) { fgets(buf, sizeof(buf), kat64); fprintf(stderr, buf ); fgets(buf, sizeof(buf), kat64); sscanf(buf, "%d %d", &Timesteps, &Nodes); fgets(buf, sizeof(buf), kat64); fprintf(stderr, buf ); } if (kat74) { fgets(buf, sizeof(buf), kat74); fprintf(stderr, buf ); fgets(buf, sizeof(buf), kat74); sscanf(buf, "%d %d", &Timesteps, &Nodes); fgets(buf, sizeof(buf), kat74); fprintf(stderr, buf ); } if (elev) { fgets(buf, sizeof(buf), elev); fprintf(stderr, buf ); fgets(buf, sizeof(buf), elev); sscanf(buf, " %d %d", &Timesteps, &Nodes); // fprintf(stderr, "----> %s (%d: t=%d, n=%d)\n", // buf,i, Timesteps, Nodes ); fgets(buf, sizeof(buf), elev); fprintf(stderr, buf ); } fprintf(stderr,"ADCIRC Surface Data: %d time steps, %d vertices\n", Timesteps, Nodes); if (Nodes != Coords.size() ) { fputs("Vertices and Coordinates mismatch.", stderr); return; } if(time_steps_to_load > 0 && (unsigned)time_steps_to_load < Timesteps ) Timesteps = time_steps_to_load; vector<float> Elevation; vector<F5_vec3f_t> Wind; vector<F5_vec3f_t> DepthAverageVelocity; if (kat63 || elev) Elevation.resize( Coords.size() ); if (kat74) Wind.resize( Coords.size() ); if (kat64) DepthAverageVelocity.resize( Coords.size() ); F5Path*FirstSurface = 0; for(unsigned tstep=0; tstep < Timesteps; tstep++) { // Time unit is minutes since... driving MM5 data start. // one time step is 30 mins double Time = time + 30*tstep; fprintf(stderr, "Timestep %d -> t=%lg\n", tstep, Time); F5Path *myGrid; if (!FirstSurface) { FirstSurface = F5write_triangular_surface(fileID, Time, gridname.c_str(), &Coords[0], Coords.size(), &Triangles[0], Triangles.size() ); myGrid = FirstSurface; } else myGrid = F5Flink_triangular_surface(FirstSurface, Time, gridname.c_str(), 0, Coords.size(), // share vertices 0, Triangles.size() // share triangles ); if (kat63) { for(unsigned int i=0; i<Coords.size(); i++) { int num; fgets(buf, sizeof(buf), kat63); if ( sscanf(buf, "%d %g\n", &num, &Elevation[i] ) != 2 ) { fprintf(stderr, "Error can't find scalar value for vertex %d\n", i ); return; } // ??? No idea what this is good for ??? if ( Elevation[i] < -3333 ) Elevation[i] = 0.0; } fgets(buf, sizeof(buf), kat63); F5Fwrite_1D(myGrid, "height", Coords.size(), H5T_NATIVE_FLOAT, H5T_NATIVE_FLOAT, &Elevation[0], F5P_DEFAULT); } if (elev) { for(unsigned i=0; i<Coords.size(); i++) {unsigned int num; fgets(buf, sizeof(buf), elev); if ( sscanf(buf, "%d %g\n", &num, &Elevation[i] ) != 2 ) { fprintf(stderr, "Error can't find scalar value for vertex %d\n", i ); return; } // ??? No idea what this is good for ??? //if ( Elevation[i] < -3333 ) Elevation[i] = 0.0; if (num != i+1) fprintf(stderr, "Index mismatch: %s\n", buf); } fgets(buf, sizeof(buf), elev); F5Fwrite_1D(myGrid, "elevation", Coords.size(), H5T_NATIVE_FLOAT, H5T_NATIVE_FLOAT, &Elevation[0], F5P_DEFAULT); } if (kat74) { for(unsigned int i=0; i<Coords.size(); i++) { int num; fgets(buf, sizeof(buf), kat74); if ( sscanf(buf, "%d %g %g\n", &num, &Wind[i].x, &Wind[i].y ) != 3 ) { fprintf(stderr, "Error can't find vector value for vertex %d\n", i ); return; } Wind[i].z = 0.0; } fgets(buf, sizeof(buf), kat74); F5Fwrite_1D(myGrid, "wind", Coords.size(), F5T_VEC3_FLOAT, F5T_VEC3_FLOAT, &Wind[0], F5P_DEFAULT); } if (kat64) { for(unsigned int i=0; i<Coords.size(); i++) { int num; fgets(buf, sizeof(buf), kat64); if ( sscanf(buf, "%d %g %g\n", &num, &DepthAverageVelocity[i].x, &DepthAverageVelocity[i].y ) != 3 ) { fprintf(stderr, "Error can't find vector value for vertex %d\n", i ); return; } DepthAverageVelocity[i].z = 0.0; } fgets(buf, sizeof(buf), kat64); F5Fwrite_1D(myGrid, "depth average velocity", Coords.size(), F5T_VEC3_FLOAT, F5T_VEC3_FLOAT, &DepthAverageVelocity[0], F5P_DEFAULT); } if (myGrid != FirstSurface) F5Fclose( myGrid ); } F5Fclose(FirstSurface); if (kat63) fclose(kat63); if (kat74) fclose(kat74); if (kat64) fclose(kat64); if (elev) pclose(elev); } void flip(const vector<F5_vec3_point_t>&Coords, F5_triangle_t &T) { F5_vec3_point_t A = Coords[T.i], B = Coords[T.j], C = Coords[T.k], BA = { B.x - A.x, B.y - A.y, B.z - A.z }, CA = { C.x - A.x, C.y - A.y, C.z - A.z }; double zDir = BA.x*CA.y - CA.x*BA.y; if (zDir < 0) { int tmp = T.j; T.j = T.k; T.k = tmp; } } double zDir(const vector<F5_vec3_point_t>&Coords, F5_triangle_t &T) { F5_vec3_point_t A = Coords[T.i], B = Coords[T.j], C = Coords[T.k], BA = { B.x - A.x, B.y - A.y, B.z - A.z }, CA = { C.x - A.x, C.y - A.y, C.z - A.z }; double zDir = BA.x*CA.y - CA.x*BA.y; return zDir; } void flip(F5_triangle_t &T) { int tmp = T.j; T.j = T.k; T.k = tmp; } int main(int argc, char*argv[]) { double time = 0.0; hid_t fileID; string gridname; char buf[1024]; int nCoords = 0, nTriangles = 0; FILE*what; int time_steps_to_load = -1; if (argc<2) { puts("Usage: ADCIRCtoF5 <inputfile.grd> <number of time steps>"); puts("Time step argument may be ommitted to convert all time steps"); return 1; } if (argc > 2) { time_steps_to_load = atoi(argv[2]); } /* putenv("F5_VERBOSITY=30");*/ strcpy(buf, argv[1]); { char*c = strrchr(buf, '.'); if (c) { *c = 0; strcat(c, ".f5"); fileID = H5Fcreate(buf, /* The FileName */ H5F_ACC_TRUNC, /* Creation Flags */ H5P_DEFAULT, /* File creation property list identifier */ H5P_DEFAULT); /* File access property list identifier */ } else { puts("invalid filename."); return 3; } *c = 0; gridname = buf; } what = fopen(argv[1], "rt"); if (!what) { printf("Cannot open '%s'. That's not good.\n", argv[1] ); return 2; } fgets(buf, sizeof(buf), what); /* get no of points and no of triangles here*/ if (fscanf(what, "%d %d\n",&nTriangles, &nCoords) !=2) { fprintf(stderr,"Error cant find no of points and triangles in header\n"); return 5; } fprintf(stderr,"Read no of coords %d and triangles %d\n",nCoords, nTriangles); fflush(stderr); vector<F5_vec3_point_t> Coords; Coords.resize( nCoords ); /* VERTEX Information */ {for (int i=0; i<nCoords; i++) { int index; double x,y,z; if (fscanf(what, "%d %lg %lg %lg\n",&index, &x, &y, &z) !=4) { fprintf(stderr,"Error cant find 3 coords for index %d\n",i); return 5; } assert(index == i+ 1); Coords[i].x = x; Coords[i].y = y; Coords[i].z = z; //need to scale down to see something .. it's difficult with //meters/degrees ratio :( //some scaling... which was done before. removed by marcel //Coords[i].z /= 15000.0; //Coords[i].z *= -1; //fprintf(stderr,"Reading coords idx:%d (%f, %f, %f)\n",i,Coords[i].x,Coords[i].y,Coords[i].z); }} /* allocate for the triangle buffer */ vector<F5_triangle_t> Triangles; Triangles.resize( nTriangles ); vector<F5_triangle_t> Levees; vector<vector<int> > LeveeComplex; {for (int i=0; i<nTriangles; i++) { int index; int nodes_nr; int idx[3]; if (fscanf(what, "%d %d %d %d %d\n",&index, &nodes_nr, &idx[0], &idx[1], &idx[2]) !=5) { fprintf(stderr,"Error cant find 3 vertices for index %d\n",i); return 5; } assert(index == i + 1); assert(nodes_nr == 3); Triangles[i].i = idx[0]-1; Triangles[i].j = idx[1]-1; Triangles[i].k = idx[2]-1; //fprintf(stderr,"Reading triangles idx:%d (%d, %d, %d)\n",i,Triangles[i].i,Triangles[i].j, Triangles[i].k); }} /* see http://www.adcirc.org/documentv46/fort_14.html */ vector<vector<int> > WeirComplex; int DoID = -1; char *ADCIRC_ID = getenv("ADCIRC_ID"); if (ADCIRC_ID) DoID = atoi( ADCIRC_ID ); for(int N=0; N<2; N++) { int Number; fgets(buf, sizeof(buf), what); if (sscanf(buf, "%d", &Number) < 1) { puts("ERROR in decoding object number"); break; } printf("Got Object number line: %s", buf); int TotalNodeNumber; fgets(buf, sizeof(buf), what); if (sscanf(buf, "%d", &TotalNodeNumber) < 1) { puts("ERROR in getting the total number of nodes"); break; } printf("Got Object total node number: %s", buf); for(int ID=0; ID<Number; ID++) { int nodes, type; fgets(buf, sizeof(buf), what); int got = sscanf(buf, "%d %d", &nodes, &type); // printf("Got Object line: %s", buf); if (got < 2) { if (got<1) { puts("ERROR in decoding boundary type"); printf("LINE: %s\n", buf); return 1; } // open boundary for(int i=0; i<nodes; i++) { int I; fgets(buf, sizeof(buf), what); sscanf(buf, "%d\n", &I); } puts("open boundary done."); } // NBVV(k,j) include this line only if IBTYPE(k) = 0, 1, 2, 10, 11, 12, 20, 21, 22, 30 else if ( type==0 || type==1 || type==2 || type==10 || type==11 ||type==11 || type==12 || type==20 || type==21 || type==22 || type == 52 // new flux boundary?? ) // have a boundary { int I; for(int i=0; i<nodes; i++) { fgets(buf, sizeof(buf), what); if (sscanf(buf, "%d\n", &I)<1) { puts("ERROR: cannot decode boundary"); break; } } printf("Boundary No. %d done\n", ID); fflush(stdout); } // NBVV(k,j), BARLANHT(k,j), BARLANCFSP(k,j ) - include this line only if IBTYPE(k) = 3, 13, 23 else if ( type==3 || type==13 || type==23 ) { int I; double height, CoeffA; WeirComplex.push_back( vector<int>() ); LeveeComplex.push_back( vector<int>() ); fgets(buf, sizeof(buf), what); if (sscanf(buf, "%d %lg %lg\n", &I, &height, &CoeffA )<3) { puts("ERROR: cannot decode weirs"); break; } I--; // Convert from FORTRAN to C indexing double HeightScale = 0.001; height *= HeightScale; int VI; VI = Coords.size(); F5_vec3_point_t VIc = Coords[I]; VIc.z = height; Coords.push_back( VIc ); // int VJ; // VJ = Coords.size(); // F5_vec3_point_t VJc = Coords[J]; // VJc.z = height; // Coords.push_back( VJc ); std::cerr << "Weir " << ID << "found - adding " << 2*nodes << "triangles (now at " << Triangles.size() << ")" << std::endl; //fprintf(stderr, "Weir %d found - adding %d triangles (now at %d)\n", ID, 2*nodes, (int)Triangles.size() ); fflush(stderr); //bool flipit = false; for(int i=0; i<nodes-1; i++) { int Inew; fgets(buf, sizeof(buf), what); if (sscanf(buf, "%d %lg %lg\n", &Inew, &height, &CoeffA)<3) { puts("ERROR: cannot decode weirs"); break; } Inew--; // Convert from FORTRAN to C indexing /* if (DoID>0 && ID != DoID) continue; // Close weirs { F5_triangle_t left, right; left.i = I; left.j = J; left.k = Inew; right.i=Inew; right.j = J; right.k = Jnew; if (i==0) { flipit = (zDir( Coords, left) <= 0.0) //|| //(zDir( Coords, right) <= 0.0) ; } //if (i>1) continue; // printf("LeftZDir: %lg\n", zDir( Coords, left) ); // printf("RightZDir: %lg\n", zDir( Coords, right) ); WeirComplex.back().push_back( Triangles.size() ); if (flipit) flip(left); Triangles.push_back( left ); WeirComplex.back().push_back( Triangles.size() ); if (flipit) flip(right); Triangles.push_back( right ); } */ /* // Construct Levee { height *= HeightScale; int VInew; VInew = Coords.size(); F5_vec3_point_t VIc = Coords[Inew]; VIc.z = height; Coords.push_back( VIc ); int VJnew; VJnew = Coords.size(); F5_vec3_point_t VJc = Coords[Jnew]; VJc.z = height; Coords.push_back( VJc ); F5_triangle_t Trngl; Trngl.i = I; Trngl.j = VI; Trngl.k = Inew; if (flipit) flip(Trngl); LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl ); Trngl.i = Inew; Trngl.j = VI; Trngl.k = VInew; if (flipit) flip(Trngl); LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl ); Trngl.i = VI; Trngl.j = VJ; Trngl.k = VInew; if (flipit) flip(Trngl); LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl ); Trngl.i = VJ; Trngl.j = VJnew; Trngl.k = VInew; if (flipit) flip(Trngl); LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl ); Trngl.i = VJnew; Trngl.j = J; Trngl.k = Jnew; if (flipit) flip(Trngl); LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl ); Trngl.i = VJnew; Trngl.j = VJ; Trngl.k = J; if (flipit) flip(Trngl); LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl ); VI = VInew; VJ = VJnew; } */ // Reuse indices I = Inew; } } // NBVV(k,j), IBCONN(k,j), BARINHT(k,j), BARINCFSB(k,j), BARINCFSP(k,j) - include this line only if IBTYPE(k) = 4, 24 else if (type==24 || type== 4) // have a weir! { int I, J; double height, CoeffA, CoeffB; WeirComplex.push_back( vector<int>() ); LeveeComplex.push_back( vector<int>() ); fgets(buf, sizeof(buf), what); if (sscanf(buf, "%d %d %lg %lg %lg\n", &I, &J, &height, &CoeffA, &CoeffB)<5) { puts("ERROR: cannot decode weirs"); break; } I--; J--; // Convert from FORTRAN to C indexing double HeightScale = 0.001; height *= HeightScale; int VI; VI = Coords.size(); F5_vec3_point_t VIc = Coords[I]; VIc.z = height; Coords.push_back( VIc ); int VJ; VJ = Coords.size(); F5_vec3_point_t VJc = Coords[J]; VJc.z = height; Coords.push_back( VJc ); std::cerr << "Weir " << ID << " found - adding " << 2*nodes << " triangles (now at " << Triangles.size() << ")" << std::endl; //fprintf(stderr, "Weir %d found - adding %d triangles (now at %d)\n", ID, 2*nodes, (int)Triangles.size() ); fflush(stderr); bool flipit = false; for(int i=0; i<nodes-1; i++) { int Inew, Jnew; fgets(buf, sizeof(buf), what); if (sscanf(buf, "%d %d %lg %lg %lg\n", &Inew, &Jnew, &height, &CoeffA, &CoeffB)<5) { puts("ERROR: cannot decode weirs"); break; } Inew--; Jnew--; // Convert from FORTRAN to C indexing if (DoID>0 && ID != DoID) continue; // Close weirs { F5_triangle_t left, right; left.i = I; left.j = J; left.k = Inew; right.i=Inew; right.j = J; right.k = Jnew; if (i==0) { flipit = (zDir( Coords, left) <= 0.0) //|| //(zDir( Coords, right) <= 0.0) ; } //if (i>1) continue; // printf("LeftZDir: %lg\n", zDir( Coords, left) ); // printf("RightZDir: %lg\n", zDir( Coords, right) ); WeirComplex.back().push_back( Triangles.size() ); if (flipit) flip(left); Triangles.push_back( left ); WeirComplex.back().push_back( Triangles.size() ); if (flipit) flip(right); Triangles.push_back( right ); } // Construct Levee { height *= HeightScale; int VInew; VInew = Coords.size(); F5_vec3_point_t VIc = Coords[Inew]; VIc.z = height; Coords.push_back( VIc ); int VJnew; VJnew = Coords.size(); F5_vec3_point_t VJc = Coords[Jnew]; VJc.z = height; Coords.push_back( VJc ); F5_triangle_t Trngl; Trngl.i = I; Trngl.j = VI; Trngl.k = Inew; if (flipit) flip(Trngl); LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl ); Trngl.i = Inew; Trngl.j = VI; Trngl.k = VInew; if (flipit) flip(Trngl); LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl ); Trngl.i = VI; Trngl.j = VJ; Trngl.k = VInew; if (flipit) flip(Trngl); LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl ); Trngl.i = VJ; Trngl.j = VJnew; Trngl.k = VInew; if (flipit) flip(Trngl); LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl ); Trngl.i = VJnew; Trngl.j = J; Trngl.k = Jnew; if (flipit) flip(Trngl); LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl ); Trngl.i = VJnew; Trngl.j = VJ; Trngl.k = J; if (flipit) flip(Trngl); LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl ); VI = VInew; VJ = VJnew; } // Reuse indices I = Inew; J = Jnew; } } else { puts("ERROR in decoding ADCIRC file"); printf("ERROR LINE: %s\n", buf); printf("Unsupported type = %d", type); return 1; break; } } printf("Boundary Complex %d done.\n", N); } std::cerr << "Triangular surface: " << nCoords << " vertices, " << Triangles.size() << " triangles - " << Triangles.size() - nTriangles << " for weirs" << std::endl; // fprintf(stderr, "Triangular surface: %d vertices, %d triangles - %d for weirs\n", // nCoords, (int)Triangles.size(), (int)Triangles.size() - nTriangles ); fflush(stderr); F5Path *MainGrid = F5write_triangular_surface(fileID, time, gridname.c_str(), &Coords[0], nCoords, &Triangles[0], Triangles.size() ); if (!MainGrid) puts("ERROR in writing grid..."); F5Fclose(MainGrid); std::cerr << "Levee Triangular surface: " << Coords.size() << " vertices, " << Levees.size() << " triangles" << std::endl; // fprintf(stderr, "Levee Triangular surface: %d vertices, %d triangles\n", // (int)Coords.size(), (int)Levees.size() ); fflush(stderr); // No levee data found, don't write if( Levees.size()>0 ) { F5Path *LeveeGrid = F5write_triangular_surface(fileID, time, (gridname + ".levee").c_str(), &Coords[0], Coords.size(), &Levees[0], Levees.size() ); if (!LeveeGrid) puts("ERROR in writing grid..."); F5Fclose(LeveeGrid); } // Due to the levees, we have additional vertices now. // In order to write surface data, we need to crop them // to their original count. Coords.resize( nCoords ); AppendSurfaceData( fileID, gridname, Coords, Triangles, time_steps_to_load ); H5Fclose(fileID); return 0; }
1.6.2