diff --git a/appl/lecture/msm/mcwhorter/dumux-in b/appl/lecture/msm/mcwhorter/dumux-in index e7f31b2f86ebe69d4c325e6b9527fa2916655bcd..c8100147f3457aa2dc95eb38d980022e18e9cc8a 100644 --- a/appl/lecture/msm/mcwhorter/dumux-in +++ b/appl/lecture/msm/mcwhorter/dumux-in @@ -20,5 +20,5 @@ viscosityNW = 1e-3 [problem] -DiscretizationLength = 0.05 +DiscretizationLength = 0.08 tEnd = 1000 diff --git a/appl/lecture/msm/mcwhorter/mcwhorter_ff.cc b/appl/lecture/msm/mcwhorter/mcwhorter_ff.cc index d7ef56268a61a7dca2f7871646a585a19b6c2a35..532aea6100022f9a307a0012814abbb04541b0bd 100644 --- a/appl/lecture/msm/mcwhorter/mcwhorter_ff.cc +++ b/appl/lecture/msm/mcwhorter/mcwhorter_ff.cc @@ -69,7 +69,7 @@ int main(int argc, char** argv) Dune::FieldVector<Scalar, dim> upperRight(2.0); // UpperRight[1] = 1; Dune::FieldVector<int, dim> cellNumbers(cellNumber); - cellNumbers[0] = 26; +// cellNumbers[0] = 26; double tEnd = Params::tree().get<double>("problem.tEnd"); double dt = tEnd; diff --git a/appl/lecture/msm/mcwhorter/mcwhorterproblem.hh b/appl/lecture/msm/mcwhorter/mcwhorterproblem.hh index db39bedc7b19e7783354ab05f8c70c41d75b0012..c0a071074368f950871bfcb6c0e5cd748377d0b6 100644 --- a/appl/lecture/msm/mcwhorter/mcwhorterproblem.hh +++ b/appl/lecture/msm/mcwhorter/mcwhorterproblem.hh @@ -253,6 +253,10 @@ public: if (time < 0) return; + double discretizationLength = Params::tree().template get<double>("problem.DiscretizationLength"); + int cellNumberX = static_cast<int>(2/discretizationLength); + discretizationLength = 2.0/cellNumberX; //Might be slightly different from the input parameter + std::ofstream dataFile; dataFile.open("vipplot.vgf", std::fstream::app); @@ -264,30 +268,43 @@ public: dataFile << "# title Mc Whorther Problem\n"; dataFile << "# x-label x\n"; dataFile << "# y-label Saturation\n"; - dataFile << "# color 0 0 255\n"; dataFile << "# x-range -0.1 2.1\n"; dataFile << "# y-range 0 1\n"; -// dataFile << "# symbol none"; -// dataFile << "# linestyle solid"; - ElementIterator eItEnd = this->gridView().template end<0>(); - for (ElementIterator eIt = this->gridView().template begin<0>(); eIt != eItEnd; ++eIt) - { - int currentIdx = this->variables().index(*eIt); - Scalar sat = this->variables().saturation()[currentIdx]; - Scalar leftPos = eIt->geometry().corner(0)[0]; - Scalar rightPos = eIt->geometry().corner(1)[0]; - dataFile << leftPos <<" "<<sat<<" "<<rightPos<<" "<<sat<<"\n"; - } + dataFile << "# color 0 0 255\n"; + + Scalar curSat = this->variables().saturation()[0]; + + // Draw first piece separately, as no vertical line is needed + dataFile << 0 << " " << curSat << " " + << discretizationLength << " " << curSat <<std::endl; + for (int i=1; i < cellNumberX; i++) + { + // Vertical Line + dataFile << i*discretizationLength << " " + << curSat << " "; + curSat = this->variables().saturation()[i]; + dataFile << i*discretizationLength << " " << curSat << std::endl; + //Horizontal Line + dataFile << i*discretizationLength << " " << curSat << " " + << (i+1)*discretizationLength << " " << curSat << std::endl; + } dataFile << "# color 255 0 0\n"; - for (ElementIterator eIt = this->gridView().template begin<0>(); eIt != eItEnd; ++eIt) - { - int currentIdx = this->variables().index(*eIt); - Scalar sat = analyticSolution_.AnalyticSolution()[currentIdx]; - Scalar leftPos = eIt->geometry().corner(0)[0]; - Scalar rightPos = eIt->geometry().corner(1)[0]; - dataFile << leftPos <<" "<<sat<<" "<<rightPos<<" "<<sat<<"\n"; - } + + curSat = analyticSolution_.AnalyticSolution()[0]; + dataFile << 0 << " " << curSat << " " + << discretizationLength << " " << curSat <<std::endl; + for (int i=1; i < cellNumberX; i++) + { + // Vertical Line + dataFile << i*discretizationLength << " " + << curSat << " "; + curSat = analyticSolution_.AnalyticSolution()[i]; + dataFile << i*discretizationLength << " " << curSat << std::endl; + //Horizontal Line + dataFile << i*discretizationLength << " " << curSat << " " + << (i+1)*discretizationLength << " " << curSat << std::endl; + } dataFile << "# legend Numerical Solution,Analytic Solution\n"; dataFile.close(); }