OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Using Datamapper.h in writeVtkData

Browse Source
This commit is contained in:
Joakim Hove 2012-06-28 13:20:18 +02:00
parent d782ba15cb
commit 01f43b18d7
2 changed files with 257 additions and 258 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

501
opm/core/utility/writeVtkData.cpp
View File

@ -18,6 +18,7 @@
*/ */
#include <opm/core/utility/writeVtkData.hpp> #include <opm/core/utility/writeVtkData.hpp>
#include <opm/core/utility/DataMap.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <boost/lexical_cast.hpp> #include <boost/lexical_cast.hpp>
@ -34,63 +35,63 @@ namespace Opm
{ {
void writeVtkData(const std::tr1::array<int, 3>& dims, void writeVtkData(const std::tr1::array<int, 3>& dims,
const std::tr1::array<double, 3>& cell_size, const std::tr1::array<double, 3>& cell_size,
const DataMap& data, const DataMap& data,
std::ostream& os) std::ostream& os)
{ {
// Dimension is hardcoded in the prototype and the next two lines, // Dimension is hardcoded in the prototype and the next two lines,
// but the rest is flexible (allows dimension == 2 or 3). // but the rest is flexible (allows dimension == 2 or 3).
int dimension = 3; int dimension = 3;
int num_cells = dims[0]*dims[1]*dims[2]; int num_cells = dims[0]*dims[1]*dims[2];
ASSERT(dimension == 2 || dimension == 3); ASSERT(dimension == 2 || dimension == 3);
ASSERT(num_cells == dims[0]*dims[1]* (dimension == 2 ? 1 : dims[2])); ASSERT(num_cells == dims[0]*dims[1]* (dimension == 2 ? 1 : dims[2]));
os << "# vtk DataFile Version 2.0\n"; os << "# vtk DataFile Version 2.0\n";
os << "Structured Grid\n \n"; os << "Structured Grid\n \n";
os << "ASCII \n"; os << "ASCII \n";
os << "DATASET STRUCTURED_POINTS\n"; os << "DATASET STRUCTURED_POINTS\n";
os << "DIMENSIONS " os << "DIMENSIONS "
<< dims[0] + 1 << " " << dims[0] + 1 << " "
<< dims[1] + 1 << " "; << dims[1] + 1 << " ";
if (dimension == 3) { if (dimension == 3) {
os << dims[2] + 1; os << dims[2] + 1;
} else { } else {
os << 1; os << 1;
} }
os << "\n"; os << "\n";
os << "ORIGIN " << 0.0 << " " << 0.0 << " " << 0.0 << "\n"; os << "ORIGIN " << 0.0 << " " << 0.0 << " " << 0.0 << "\n";
os << "SPACING " << cell_size[0] << " " << cell_size[1]; os << "SPACING " << cell_size[0] << " " << cell_size[1];
if (dimension == 3) { if (dimension == 3) {
os << " " << cell_size[2]; os << " " << cell_size[2];
} else { } else {
os << " " << 0.0; os << " " << 0.0;
} }
os << "\n"; os << "\n";
os << "\nCELL_DATA " << num_cells << '\n'; os << "\nCELL_DATA " << num_cells << '\n';
for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) { for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) {
std::string name = dit->first; std::string name = dit->first;
os << "SCALARS " << name << " float" << '\n'; os << "SCALARS " << name << " float" << '\n';
os << "LOOKUP_TABLE " << name << "_table " << '\n'; os << "LOOKUP_TABLE " << name << "_table " << '\n';
const std::vector<double>& field = *(dit->second); const std::vector<double>& field = *(dit->second);
// We always print only the first data item for every // We always print only the first data item for every
// cell, using 'stride'. // cell, using 'stride'.
// This is a hack to get water saturation nicely. // This is a hack to get water saturation nicely.
// \TODO: Extend to properly printing vector data. // \TODO: Extend to properly printing vector data.
const int stride = field.size()/num_cells; const int stride = field.size()/num_cells;
const int num_per_line = 5; const int num_per_line = 5;
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
os << field[stride*c] << ' '; os << field[stride*c] << ' ';
if (c % num_per_line == num_per_line - 1 if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) { || c == num_cells - 1) {
os << '\n'; os << '\n';
} }
} }
} }
} }
typedef std::map<std::string, std::string> PMap; typedef std::map<std::string, std::string> PMap;
@ -98,219 +99,219 @@ namespace Opm
struct Tag struct Tag
{ {
Tag(const std::string& tag, const PMap& props, std::ostream& os) Tag(const std::string& tag, const PMap& props, std::ostream& os)
: name_(tag), os_(os) : name_(tag), os_(os)
{ {
indent(os); indent(os);
os << "<" << tag; os << "<" << tag;
for (PMap::const_iterator it = props.begin(); it != props.end(); ++it) { for (PMap::const_iterator it = props.begin(); it != props.end(); ++it) {
os << " " << it->first << "=\"" << it->second << "\""; os << " " << it->first << "=\"" << it->second << "\"";
} }
os << ">\n"; os << ">\n";
++indent_; ++indent_;
} }
Tag(const std::string& tag, std::ostream& os) Tag(const std::string& tag, std::ostream& os)
: name_(tag), os_(os) : name_(tag), os_(os)
{ {
indent(os); indent(os);
os << "<" << tag << ">\n"; os << "<" << tag << ">\n";
++indent_; ++indent_;
} }
~Tag() ~Tag()
{ {
--indent_; --indent_;
indent(os_); indent(os_);
os_ << "</" << name_ << ">\n"; os_ << "</" << name_ << ">\n";
} }
static void indent(std::ostream& os) static void indent(std::ostream& os)
{ {
for (int i = 0; i < indent_; ++i) { for (int i = 0; i < indent_; ++i) {
os << " "; os << " ";
} }
} }
private: private:
static int indent_; static int indent_;
std::string name_; std::string name_;
std::ostream& os_; std::ostream& os_;
}; };
int Tag::indent_ = 0; int Tag::indent_ = 0;
void writeVtkData(const UnstructuredGrid& grid, void writeVtkData(const UnstructuredGrid& grid,
const DataMap& data, const DataMap& data,
std::ostream& os) std::ostream& os)
{ {
if (grid.dimensions != 3) { if (grid.dimensions != 3) {
THROW("Vtk output for 3d grids only"); THROW("Vtk output for 3d grids only");
} }
os.precision(12); os.precision(12);
os << "<?xml version=\"1.0\"?>\n"; os << "<?xml version=\"1.0\"?>\n";
PMap pm; PMap pm;
pm["type"] = "UnstructuredGrid"; pm["type"] = "UnstructuredGrid";
Tag vtkfiletag("VTKFile", pm, os); Tag vtkfiletag("VTKFile", pm, os);
Tag ugtag("UnstructuredGrid", os); Tag ugtag("UnstructuredGrid", os);
int num_pts = grid.number_of_nodes; int num_pts = grid.number_of_nodes;
int num_cells = grid.number_of_cells; int num_cells = grid.number_of_cells;
pm.clear(); pm.clear();
pm["NumberOfPoints"] = boost::lexical_cast<std::string>(num_pts); pm["NumberOfPoints"] = boost::lexical_cast<std::string>(num_pts);
pm["NumberOfCells"] = boost::lexical_cast<std::string>(num_cells); pm["NumberOfCells"] = boost::lexical_cast<std::string>(num_cells);
Tag piecetag("Piece", pm, os); Tag piecetag("Piece", pm, os);
{ {
Tag pointstag("Points", os); Tag pointstag("Points", os);
pm.clear(); pm.clear();
pm["type"] = "Float64"; pm["type"] = "Float64";
pm["Name"] = "Coordinates"; pm["Name"] = "Coordinates";
pm["NumberOfComponents"] = "3"; pm["NumberOfComponents"] = "3";
pm["format"] = "ascii"; pm["format"] = "ascii";
Tag datag("DataArray", pm, os); Tag datag("DataArray", pm, os);
for (int i = 0; i < num_pts; ++i) { for (int i = 0; i < num_pts; ++i) {
Tag::indent(os); Tag::indent(os);
os << grid.node_coordinates[3*i + 0] << ' ' os << grid.node_coordinates[3*i + 0] << ' '
<< grid.node_coordinates[3*i + 1] << ' ' << grid.node_coordinates[3*i + 1] << ' '
<< grid.node_coordinates[3*i + 2] << '\n'; << grid.node_coordinates[3*i + 2] << '\n';
} }
} }
{ {
Tag cellstag("Cells", os); Tag cellstag("Cells", os);
pm.clear(); pm.clear();
pm["type"] = "Int32"; pm["type"] = "Int32";
pm["NumberOfComponents"] = "1"; pm["NumberOfComponents"] = "1";
pm["format"] = "ascii"; pm["format"] = "ascii";
std::vector<int> cell_numpts; std::vector<int> cell_numpts;
cell_numpts.reserve(num_cells); cell_numpts.reserve(num_cells);
{ {
pm["Name"] = "connectivity"; pm["Name"] = "connectivity";
Tag t("DataArray", pm, os); Tag t("DataArray", pm, os);
int hf = 0; int hf = 0;
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
std::set<int> cell_pts; std::set<int> cell_pts;
for (; hf < grid.cell_facepos[c+1]; ++hf) { for (; hf < grid.cell_facepos[c+1]; ++hf) {
int f = grid.cell_faces[hf]; int f = grid.cell_faces[hf];
const int* fnbeg = grid.face_nodes + grid.face_nodepos[f]; const int* fnbeg = grid.face_nodes + grid.face_nodepos[f];
const int* fnend = grid.face_nodes + grid.face_nodepos[f+1]; const int* fnend = grid.face_nodes + grid.face_nodepos[f+1];
cell_pts.insert(fnbeg, fnend); cell_pts.insert(fnbeg, fnend);
} }
cell_numpts.push_back(cell_pts.size()); cell_numpts.push_back(cell_pts.size());
Tag::indent(os); Tag::indent(os);
std::copy(cell_pts.begin(), cell_pts.end(), std::copy(cell_pts.begin(), cell_pts.end(),
std::ostream_iterator<int>(os, " ")); std::ostream_iterator<int>(os, " "));
os << '\n'; os << '\n';
} }
} }
{ {
pm["Name"] = "offsets"; pm["Name"] = "offsets";
Tag t("DataArray", pm, os); Tag t("DataArray", pm, os);
int offset = 0; int offset = 0;
const int num_per_line = 10; const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) { if (c % num_per_line == 0) {
Tag::indent(os); Tag::indent(os);
} }
offset += cell_numpts[c]; offset += cell_numpts[c];
os << offset << ' '; os << offset << ' ';
if (c % num_per_line == num_per_line - 1 if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) { || c == num_cells - 1) {
os << '\n'; os << '\n';
} }
} }
} }
std::vector<int> cell_foffsets; std::vector<int> cell_foffsets;
cell_foffsets.reserve(num_cells); cell_foffsets.reserve(num_cells);
{ {
pm["Name"] = "faces"; pm["Name"] = "faces";
Tag t("DataArray", pm, os); Tag t("DataArray", pm, os);
const int* fp = grid.cell_facepos; const int* fp = grid.cell_facepos;
int offset = 0; int offset = 0;
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
Tag::indent(os); Tag::indent(os);
os << fp[c+1] - fp[c] << '\n'; os << fp[c+1] - fp[c] << '\n';
++offset; ++offset;
for (int hf = fp[c]; hf < fp[c+1]; ++hf) { for (int hf = fp[c]; hf < fp[c+1]; ++hf) {
int f = grid.cell_faces[hf]; int f = grid.cell_faces[hf];
const int* np = grid.face_nodepos; const int* np = grid.face_nodepos;
int f_num_pts = np[f+1] - np[f]; int f_num_pts = np[f+1] - np[f];
Tag::indent(os); Tag::indent(os);
os << f_num_pts << ' '; os << f_num_pts << ' ';
++offset; ++offset;
std::copy(grid.face_nodes + np[f], std::copy(grid.face_nodes + np[f],
grid.face_nodes + np[f+1], grid.face_nodes + np[f+1],
std::ostream_iterator<int>(os, " ")); std::ostream_iterator<int>(os, " "));
os << '\n'; os << '\n';
offset += f_num_pts; offset += f_num_pts;
} }
cell_foffsets.push_back(offset); cell_foffsets.push_back(offset);
} }
} }
{ {
pm["Name"] = "faceoffsets"; pm["Name"] = "faceoffsets";
Tag t("DataArray", pm, os); Tag t("DataArray", pm, os);
const int num_per_line = 10; const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) { if (c % num_per_line == 0) {
Tag::indent(os); Tag::indent(os);
} }
os << cell_foffsets[c] << ' '; os << cell_foffsets[c] << ' ';
if (c % num_per_line == num_per_line - 1 if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) { || c == num_cells - 1) {
os << '\n'; os << '\n';
} }
} }
} }
{ {
pm["type"] = "UInt8"; pm["type"] = "UInt8";
pm["Name"] = "types"; pm["Name"] = "types";
Tag t("DataArray", pm, os); Tag t("DataArray", pm, os);
const int num_per_line = 10; const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) { if (c % num_per_line == 0) {
Tag::indent(os); Tag::indent(os);
} }
os << "42 "; os << "42 ";
if (c % num_per_line == num_per_line - 1 if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) { || c == num_cells - 1) {
os << '\n'; os << '\n';
} }
} }
} }
} }
{ {
pm.clear(); pm.clear();
if (data.find("saturation") != data.end()) { if (data.find("saturation") != data.end()) {
pm["Scalars"] = "saturation"; pm["Scalars"] = "saturation";
} else if (data.find("pressure") != data.end()) { } else if (data.find("pressure") != data.end()) {
pm["Scalars"] = "pressure"; pm["Scalars"] = "pressure";
} }
Tag celldatatag("CellData", pm, os); Tag celldatatag("CellData", pm, os);
pm.clear(); pm.clear();
pm["NumberOfComponents"] = "1"; pm["NumberOfComponents"] = "1";
pm["format"] = "ascii"; pm["format"] = "ascii";
pm["type"] = "Float64"; pm["type"] = "Float64";
for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) { for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) {
pm["Name"] = dit->first; pm["Name"] = dit->first;
const std::vector<double>& field = *(dit->second); const std::vector<double>& field = *(dit->second);
const int num_comps = field.size()/grid.number_of_cells; const int num_comps = field.size()/grid.number_of_cells;
pm["NumberOfComponents"] = boost::lexical_cast<std::string>(num_comps); pm["NumberOfComponents"] = boost::lexical_cast<std::string>(num_comps);
Tag ptag("DataArray", pm, os); Tag ptag("DataArray", pm, os);
const int num_per_line = num_comps == 1 ? 5 : num_comps; const int num_per_line = num_comps == 1 ? 5 : num_comps;
for (int item = 0; item < num_cells*num_comps; ++item) { for (int item = 0; item < num_cells*num_comps; ++item) {
if (item % num_per_line == 0) { if (item % num_per_line == 0) {
Tag::indent(os); Tag::indent(os);
} }
double value = field[item]; double value = field[item];
if (std::fabs(value) < std::numeric_limits<double>::min()) { if (std::fabs(value) < std::numeric_limits<double>::min()) {
// Avoiding denormal numbers to work around // Avoiding denormal numbers to work around
// bug in Paraview. // bug in Paraview.
value = 0.0; value = 0.0;
} }
os << value << ' '; os << value << ' ';
if (item % num_per_line == num_per_line - 1 if (item % num_per_line == num_per_line - 1
|| item == num_cells - 1) { || item == num_cells - 1) {
os << '\n'; os << '\n';
} }
} }
} }
} }
} }
} // namespace Opm } // namespace Opm

14
opm/core/utility/writeVtkData.hpp
View File

@ -26,25 +26,23 @@
#include <vector> #include <vector>
#include <tr1/array> #include <tr1/array>
#include <iosfwd> #include <iosfwd>
#include <opm/core/utility/DataMap.hpp>
struct UnstructuredGrid; struct UnstructuredGrid;
namespace Opm namespace Opm
{ {
/// Intended to map strings (giving the output field names) to data.
typedef std::map<std::string, const std::vector<double>*> DataMap;
/// Vtk output for cartesian grids. /// Vtk output for cartesian grids.
void writeVtkData(const std::tr1::array<int, 3>& dims, void writeVtkData(const std::tr1::array<int, 3>& dims,
const std::tr1::array<double, 3>& cell_size, const std::tr1::array<double, 3>& cell_size,
const DataMap& data, const DataMap& data,
std::ostream& os); std::ostream& os);
/// Vtk output for general grids. /// Vtk output for general grids.
void writeVtkData(const UnstructuredGrid& grid, void writeVtkData(const UnstructuredGrid& grid,
const DataMap& data, const DataMap& data,
std::ostream& os); std::ostream& os);
} // namespace Opm } // namespace Opm
#endif // OPM_WRITEVTKDATA_HEADER_INCLUDED #endif // OPM_WRITEVTKDATA_HEADER_INCLUDED