3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
gstatrandomfield.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_GSTAT_RANDOM_FIELD_HH
25#define DUMUX_GSTAT_RANDOM_FIELD_HH
26
27#include <iostream>
28
29#include <dune/grid/common/mcmgmapper.hh>
30#include <dune/grid/io/file/vtk.hh>
31
32namespace Dumux {
33
47template<class GridView, class Scalar>
49{
50 enum { dim = GridView::dimension };
51 enum { dimWorld = GridView::dimensionworld };
52
53 using DataVector = std::vector<Scalar>;
54 using Element = typename GridView::Traits::template Codim<0>::Entity;
55 using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
56
57public:
58 // Add field types if you want to implement e.g. tensor permeabilities.
60
67 GstatRandomField(const GridView& gridView, const ElementMapper& elementMapper)
68 : gridView_(gridView)
69 , elementMapper_(elementMapper)
70 , data_(gridView.size(0))
71 {}
72
89 void create(const std::string& gstatControlFile,
90 const std::string& gstatInputFile = "gstatInput.txt",
91 const std::string& gstatOutputFile = "permeab.dat",
92 FieldType fieldType = FieldType::scalar,
93 bool createNew = true)
94 {
95 fieldType_ = fieldType;
96 if (createNew)
97 {
98#if !HAVE_GSTAT
99 DUNE_THROW(Dune::InvalidStateException, "Requested data field generation with gstat"
100 << " but gstat was not found on your system. Set GSTAT_ROOT to the path where gstat "
101 << " is installed and pass it to CMake, e.g. through an opts file.");
102#else
103 std::ofstream gstatInput(gstatInputFile);
104 for (const auto& element : elements(gridView_))
105 // get global coordinates of cell centers
106 gstatInput << element.geometry().center() << std::endl;
107
108 gstatInput.close();
109
110 std::string syscom;
111 syscom = GSTAT_EXECUTABLE;
112 syscom += " ";
113 syscom += gstatControlFile;
114
115 if (!gstatInput.good())
116 {
117 DUNE_THROW(Dune::IOError, "Reading the gstat control file: "
118 << gstatControlFile << " failed." << std::endl);
119 }
120
121 if (system(syscom.c_str()))
122 {
123 DUNE_THROW(Dune::IOError, "Executing gstat failed.");
124 }
125#endif
126 }
127
128 std::ifstream gstatOutput(gstatOutputFile);
129 if (!gstatOutput.good())
130 {
131 DUNE_THROW(Dune::IOError, "Reading from file: "
132 << gstatOutputFile << " failed." << std::endl);
133 }
134
135 std::string line;
136 std::getline(gstatOutput, line);
137
138 Scalar trash, dataValue;
139 for (const auto& element : elements(gridView_))
140 {
141 std::getline(gstatOutput, line);
142 std::istringstream curLine(line);
143 if (dim == 1)
144 curLine >> trash >> dataValue;
145 else if (dim == 2)
146 curLine >> trash >> trash >> dataValue;
147 else if (dim == 3)
148 curLine >> trash >> trash >> trash >> dataValue;
149 else
150 DUNE_THROW(Dune::InvalidStateException, "Invalid dimension " << dim);
151
152 data_[elementMapper_.index(element)] = dataValue;
153 }
154 gstatOutput.close();
155
156 // post processing
157 using std::pow;
158 if (fieldType_ == FieldType::log10)
159 std::for_each(data_.begin(), data_.end(), [](Scalar& s){ s = pow(10.0, s); });
160 }
161
163 Scalar data(const Element& e) const
164 {
165 return data_[elementMapper_.index(e)];
166 }
167
169 const DataVector& data() const
170 {
171 return data_;
172 }
173
175 void writeVtk(const std::string& vtkName,
176 const std::string& dataName = "data") const
177 {
178 Dune::VTKWriter<GridView> vtkwriter(gridView_);
179 vtkwriter.addCellData(data_, dataName);
180
181 DataVector logPerm;
182 if (fieldType_ == FieldType::log10)
183 {
184 logPerm = data_;
185 using std::log10;
186 std::for_each(logPerm.begin(), logPerm.end(), [](Scalar& s){ s = log10(s); });
187 vtkwriter.addCellData(logPerm, "log10 of " + dataName);
188 }
189 vtkwriter.write(vtkName, Dune::VTK::OutputType::ascii);
190 }
191
192private:
193 const GridView gridView_;
194 const ElementMapper& elementMapper_;
195 DataVector data_;
196 FieldType fieldType_;
197};
198
199} // end namespace Dumux
200
201#endif
Definition: adapt.hh:29
Creating random fields using gstat.
Definition: gstatrandomfield.hh:49
void writeVtk(const std::string &vtkName, const std::string &dataName="data") const
Write the data to a vtk file.
Definition: gstatrandomfield.hh:175
Scalar data(const Element &e) const
Return an entry of the data vector.
Definition: gstatrandomfield.hh:163
void create(const std::string &gstatControlFile, const std::string &gstatInputFile="gstatInput.txt", const std::string &gstatOutputFile="permeab.dat", FieldType fieldType=FieldType::scalar, bool createNew=true)
Creates a new field with random variables, if desired. Otherwise creates a data field from already av...
Definition: gstatrandomfield.hh:89
GstatRandomField(const GridView &gridView, const ElementMapper &elementMapper)
Constructor.
Definition: gstatrandomfield.hh:67
const DataVector & data() const
Return the data vector for analysis or external vtk output.
Definition: gstatrandomfield.hh:169
FieldType
Definition: gstatrandomfield.hh:59
@ scalar
Definition: gstatrandomfield.hh:59
@ log10
Definition: gstatrandomfield.hh:59