From 48cd13228f916aee1bb2f3931b64aea3260941d5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 31 Jan 2012 22:33:48 +0100 Subject: [PATCH] Vtk output for general grids. Grid read from deck can use DIMENS or SPECGRID. --- examples/spu_2p.cpp | 252 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 241 insertions(+), 11 deletions(-) diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index c904a9d6..7137171d 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -62,6 +62,7 @@ #include #include +#include #include #include @@ -102,16 +103,23 @@ namespace Opm const std::vector& zcorn = deck.getFloatingPointValue("ZCORN"); const std::vector& coord = deck.getFloatingPointValue("COORD"); const std::vector& actnum = deck.getIntegerValue("ACTNUM"); - const std::vector& dimens = deck.getIntegerValue("DIMENS"); + std::vector dims; + if (deck.hasField("DIMENS")) { + dims = deck.getIntegerValue("DIMENS"); + } else if (deck.hasField("SPECGRID")) { + dims = deck.getSPECGRID().dimensions; + } else { + THROW("Deck must have either DIMENS or SPECGRID."); + } // Collect in input struct for preprocessing. struct grdecl grdecl; grdecl.zcorn = &zcorn[0]; grdecl.coord = &coord[0]; grdecl.actnum = &actnum[0]; - grdecl.dims[0] = dimens[0]; - grdecl.dims[1] = dimens[1]; - grdecl.dims[2] = dimens[2]; + grdecl.dims[0] = dims[0]; + grdecl.dims[1] = dims[1]; + grdecl.dims[2] = dims[2]; // Process and compute. preprocess(&grdecl, 0.0, &cpg_); @@ -281,19 +289,18 @@ compute_porevolume(const UnstructuredGrid* g, template -void outputState(const std::tr1::array& grid_dims, - const std::tr1::array& cell_size, +void outputState(const UnstructuredGrid* grid, const State& state, const int step, const std::string& output_dir) { std::ostringstream vtkfilename; - vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtk"; + vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; std::ofstream vtkfile(vtkfilename.str().c_str()); if (!vtkfile) { THROW("Failed to open " << vtkfilename.str()); } - writeVtkDataAllCartesian(grid_dims, cell_size, state, vtkfile); + writeVtkDataGeneralGrid(grid, state, vtkfile); } @@ -358,6 +365,230 @@ void writeVtkDataAllCartesian(const std::tr1::array& dims, } } +typedef std::map PMap; + + +struct Tag +{ + Tag(const std::string& tag, const PMap& props, std::ostream& os) + : name_(tag), os_(os) + { + indent(os); + os << "<" << tag; + for (PMap::const_iterator it = props.begin(); it != props.end(); ++it) { + os << " " << it->first << "=\"" << it->second << "\""; + } + os << ">\n"; + ++indent_; + } + Tag(const std::string& tag, std::ostream& os) + : name_(tag), os_(os) + { + indent(os); + os << "<" << tag << ">\n"; + ++indent_; + } + ~Tag() + { + --indent_; + indent(os_); + os_ << "\n"; + } + static void indent(std::ostream& os) + { + for (int i = 0; i < indent_; ++i) { + os << " "; + } + } +private: + static int indent_; + std::string name_; + std::ostream& os_; +}; + +int Tag::indent_ = 0; + + +template +void writeVtkDataGeneralGrid(const UnstructuredGrid* grid, + const State& state, + std::ostream& os) +{ + if (grid->dimensions != 3) { + THROW("Vtk output for 3d grids only"); + } + os.precision(12); + os << "\n"; + PMap pm; + pm["type"] = "UnstructuredGrid"; + Tag vtkfiletag("VTKFile", pm, os); + Tag ugtag("UnstructuredGrid", os); + int num_pts = grid->number_of_nodes; + int num_cells = grid->number_of_cells; + pm.clear(); + pm["NumberOfPoints"] = boost::lexical_cast(num_pts); + pm["NumberOfCells"] = boost::lexical_cast(num_cells); + Tag piecetag("Piece", pm, os); + { + Tag pointstag("Points", os); + pm.clear(); + pm["type"] = "Float64"; + pm["Name"] = "Coordinates"; + pm["NumberOfComponents"] = "3"; + pm["format"] = "ascii"; + Tag datag("DataArray", pm, os); + for (int i = 0; i < num_pts; ++i) { + Tag::indent(os); + os << grid->node_coordinates[3*i + 0] << ' ' + << grid->node_coordinates[3*i + 1] << ' ' + << grid->node_coordinates[3*i + 2] << '\n'; + } + } + { + Tag cellstag("Cells", os); + pm.clear(); + pm["type"] = "Int32"; + pm["NumberOfComponents"] = "1"; + pm["format"] = "ascii"; + std::vector cell_numpts; + cell_numpts.reserve(num_cells); + { + pm["Name"] = "connectivity"; + Tag t("DataArray", pm, os); + int hf = 0; + for (int c = 0; c < num_cells; ++c) { + std::set cell_pts; + for (; hf < grid->cell_facepos[c+1]; ++hf) { + int f = grid->cell_faces[hf]; + const int* fnbeg = grid->face_nodes + grid->face_nodepos[f]; + const int* fnend = grid->face_nodes + grid->face_nodepos[f+1]; + cell_pts.insert(fnbeg, fnend); + } + cell_numpts.push_back(cell_pts.size()); + Tag::indent(os); + std::copy(cell_pts.begin(), cell_pts.end(), + std::ostream_iterator(os, " ")); + os << '\n'; + } + } + { + pm["Name"] = "offsets"; + Tag t("DataArray", pm, os); + int offset = 0; + const int num_per_line = 10; + for (int c = 0; c < num_cells; ++c) { + if (c % num_per_line == 0) { + Tag::indent(os); + } + offset += cell_numpts[c]; + os << offset << ' '; + if (c % num_per_line == num_per_line - 1 + || c == num_cells - 1) { + os << '\n'; + } + } + } + std::vector cell_foffsets; + cell_foffsets.reserve(num_cells); + { + pm["Name"] = "faces"; + Tag t("DataArray", pm, os); + const int* fp = grid->cell_facepos; + int offset = 0; + for (int c = 0; c < num_cells; ++c) { + Tag::indent(os); + os << fp[c+1] - fp[c] << '\n'; + ++offset; + for (int hf = fp[c]; hf < fp[c+1]; ++hf) { + int f = grid->cell_faces[hf]; + const int* np = grid->face_nodepos; + int f_num_pts = np[f+1] - np[f]; + Tag::indent(os); + os << f_num_pts << ' '; + ++offset; + std::copy(grid->face_nodes + np[f], + grid->face_nodes + np[f+1], + std::ostream_iterator(os, " ")); + os << '\n'; + offset += f_num_pts; + } + cell_foffsets.push_back(offset); + } + } + { + pm["Name"] = "faceoffsets"; + Tag t("DataArray", pm, os); + const int num_per_line = 10; + for (int c = 0; c < num_cells; ++c) { + if (c % num_per_line == 0) { + Tag::indent(os); + } + os << cell_foffsets[c] << ' '; + if (c % num_per_line == num_per_line - 1 + || c == num_cells - 1) { + os << '\n'; + } + } + } + { + pm["type"] = "UInt8"; + pm["Name"] = "types"; + Tag t("DataArray", pm, os); + const int num_per_line = 10; + for (int c = 0; c < num_cells; ++c) { + if (c % num_per_line == 0) { + Tag::indent(os); + } + os << "42 "; + if (c % num_per_line == num_per_line - 1 + || c == num_cells - 1) { + os << '\n'; + } + } + } + } + { + pm.clear(); + pm["Scalars"] = "saturation"; + Tag celldatatag("CellData", pm, os); + pm.clear(); + pm["type"] = "Int32"; + pm["NumberOfComponents"] = "1"; + pm["format"] = "ascii"; + pm["type"] = "Float64"; + { + pm["Name"] = "pressure"; + Tag ptag("DataArray", pm, os); + const int num_per_line = 5; + for (int c = 0; c < num_cells; ++c) { + if (c % num_per_line == 0) { + Tag::indent(os); + } + os << state.pressure()[c] << ' '; + if (c % num_per_line == num_per_line - 1 + || c == num_cells - 1) { + os << '\n'; + } + } + } + { + pm["Name"] = "saturation"; + Tag ptag("DataArray", pm, os); + const int num_per_line = 5; + for (int c = 0; c < num_cells; ++c) { + if (c % num_per_line == 0) { + Tag::indent(os); + } + os << state.saturation()[2*c] << ' '; + if (c % num_per_line == num_per_line - 1 + || c == num_cells - 1) { + os << '\n'; + } + } + } + } +} + @@ -499,7 +730,6 @@ main(int argc, char** argv) // If we have a "deck_filename", grid and props will be read from that. bool use_deck = param.has("deck_filename"); - // UnstructuredGrid* grid = 0; boost::scoped_ptr grid; boost::scoped_ptr props; if (use_deck) { @@ -590,7 +820,7 @@ main(int argc, char** argv) << "\n" << std::endl; if (output) { - // outputState(grid->c_grid()_dims, cell_size, state, pstep, output_dir); + outputState(grid->c_grid(), state, pstep, output_dir); } psolver.solve(grid->c_grid(), totmob, src, state); @@ -622,7 +852,7 @@ main(int argc, char** argv) } if (output) { - // outputState(grid->c_grid()_dims, cell_size, state, num_psteps, output_dir); + outputState(grid->c_grid(), state, num_psteps, output_dir); } destroy_transport_source(tsrc);