diff --git a/src/dg.cc b/src/dg.cc
index 2739fcea57ec62c4d42528cce238df45f291fa27..d78f919e8e8cfd963f8f135f5da4efe5079b38fd 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -12,9 +12,8 @@
#include <dune/common/parametertreeparser.hh>
#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
-#include <dune/grid/io/file/dgfparser/dgfug.hh>
-#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/onedgrid.hh>
+#include <dune/grid/io/file/dgfparser/dgfoned.hh>
#include <dune/mmesh/mmesh.hh>
@@ -38,8 +37,8 @@ int main(int argc, char** argv)
static constexpr int dim = GRIDDIM;
- //use YaspGrid if dim = 1 and MMesh otherwise
- using Grid = std::conditional<dim == 1, Dune::YaspGrid<dim>,
+ //use OneDGrid if dim = 1 and MMesh otherwise
+ using Grid = std::conditional<dim == 1, Dune::OneDGrid,
Dune::MovingMesh<dim>>::type;
using GridFactory = Dune::DGFGridFactory<Grid>;
using GridView = typename Grid::LeafGridView;
@@ -54,12 +53,21 @@ int main(int argc, char** argv)
Dune::ParameterTree pt;
Dune::ParameterTreeParser::readINITree("parameterDG.ini", pt);
- //assume constant penalty parameter
- const double mu0 = std::stod(pt["mu0"]);
+ double mu0; //penalty parameter (prefactor)
+
+ try
+ {
+ mu0 = std::stod(pt["mu0"]);
+ }
+ catch (std::exception& e)
+ {
+ DUNE_THROW(Dune::Exception,
+ "Invalid penalty parameter or parameter 'mu0' not specified in the "
+ "file parameterDG.ini");
+ }
//create a grid from .dgf file
GridFactory gridFactory( "grids/cube" + std::to_string(dim) + "d.dgf" );
-
const Grid& grid = *gridFactory.grid();
const GridView& gridView = grid.leafGridView();
@@ -87,31 +95,34 @@ int main(int argc, char** argv)
"an infinite series; increase the cut-off if necessary!\n\n";
}
}
- else if (pt["problem"] == "dg4")
+ else
{
- if ( pt.hasKey("dmin") && pt.hasKey("dmax") )
+ double dmin; //minimum aperture
+ double dmax; //maximum aperture
+
+ try
{
- const double dmin = std::stod(pt["dmin"]);
- const double dmax = std::stod(pt["dmax"]);
- problem = new DGProblem4<Coordinate>(dmin, dmax);
+ dmin = std::stod(pt["dmin"]);
+ dmax = std::stod(pt["dmax"]);
}
- else
+ catch (std::exception& e)
{
DUNE_THROW(Dune::Exception,
- "Problem dg4 requires the parameters 'dmin' and 'dmax'. Check the "
- "file parameterDG.ini and make sure that they are specified.");
+ "The problems dg4 and dg5 require the parameters 'dmin' and 'dmax'. "
+ "Check the file parameterDG.ini and make sure that they are "
+ "specified correctly.");
}
- }
- else if (pt["problem"] == "dg5")
- {
- if ( pt.hasKey("dmin") && pt.hasKey("dmax") )
+
+ if (pt["problem"] == "dg4")
+ {
+ problem = new DGProblem4<Coordinate>(dmin, dmax);
+ }
+ else if (pt["problem"] == "dg5")
{
- const double dmin = std::stod(pt["dmin"]);
- const double dmax = std::stod(pt["dmax"]);
problem = new DGProblem5<Coordinate>(dmin, dmax);
- std::cout << "NOTE: You will not see convergence for problem dg5 as it "
- "only fulfills continuity of the normal component of the Darcy "
+ std::cout << "NOTE: You will not see convergence for problem dg5 as it"
+ " only fulfills continuity of the normal component of the Darcy "
"velocity with respect to the normal of the hyperplane Gamma = "
"{ x_1 = 0.5 } but not with respect to the normal of the actual "
"interfaces between the bulk and fracture domains. \n\n";
@@ -119,16 +130,10 @@ int main(int argc, char** argv)
else
{
DUNE_THROW(Dune::Exception,
- "Problem dg4 requires the parameters 'dmin' and 'dmax'. Check the "
- "file parameterDG.ini and make sure that they are specified.");
+ "Invalid problem type or parameter 'problem' not specified in file "
+ "parameterDG.ini (use 'dg1' to 'dg5')");
}
}
- else
- {
- DUNE_THROW(Dune::Exception,
- "Invalid problem type or parameter 'problem' not specified in file "
- "parameterDG.ini (use 'dg1' to 'dg5')");
- }
//determine solver type
SolverType solver = SolverType::Cholmod;
@@ -147,8 +152,8 @@ int main(int argc, char** argv)
}
}
+ //apply DG scheme
DG dg(gridView, mapper, *problem, solver);
-
dg(mu0);
//error and total run time
@@ -162,7 +167,7 @@ int main(int argc, char** argv)
largestEdgeLength = std::max(largestEdgeLength, edge.geometry().volume());
}
- //write error, runtime and maximum edge length
+ //write error, runtime and maximum edge length to file
std::ofstream errorFile("convergenceData", std::ios::out | std::ios::trunc);
if (errorFile.is_open())
diff --git a/src/mmdg.cc b/src/mmdg.cc
index aa8d811651f652318e63668dd327d99a93e0aae5..6897c0d2304293067cb869aaa76407f06a54450a 100644
--- a/src/mmdg.cc
+++ b/src/mmdg.cc
@@ -57,13 +57,35 @@ int main(int argc, char** argv)
Dune::ParameterTree pt;
Dune::ParameterTreeParser::readINITree("parameterMMDG.ini", pt);
- const double mu0 = std::stod(pt["mu0"]); //penalty parameter (prefactor)
- const double dmax = std::stod(pt["dmax"]); //maximum aperture
+ double mu0; //penalty parameter (prefactor)
+ double dmax; //maximum aperture
double xi; //coupling parameter
+ try
+ {
+ mu0 = std::stod(pt["mu0"]);
+ }
+ catch (std::exception& e)
+ {
+ DUNE_THROW(Dune::Exception,
+ "Invalid penalty parameter or parameter 'mu0' not specified in the "
+ "file parameterDG.ini");
+ }
+
+ try
+ {
+ dmax = std::stod(pt["dmax"]);
+ }
+ catch (std::exception& e)
+ {
+ DUNE_THROW(Dune::Exception,
+ "Invalid maximum aperture or parameter 'dmax' not specified in the "
+ "file parameterDG.ini");
+ }
+
Problem* problem; //problem to be solved
- std::string gridType;
+ std::string gridType; //prefix of the grid file name
//determine problem type
if (pt["problem"] == "mmdg1")
@@ -104,18 +126,22 @@ int main(int argc, char** argv)
}
else if (pt["problem"] == "mmdg7")
{
- if (pt.hasKey("dmin"))
+ double dmin; //minimum aperture
+
+ try
{
- problem = new MMDGProblem7<Coordinate>(std::stod(pt["dmin"]), dmax);
- gridType = "mmdg2_C";
- xi = 0.75;
+ dmin = std::stod(pt["dmin"]);
}
- else
+ catch (std::exception& e)
{
DUNE_THROW(Dune::Exception,
"Problem dg4 requires the parameter 'dmin'. Check the "
"file parameterMMDG.ini and make sure that it is specified.");
}
+
+ problem = new MMDGProblem7<Coordinate>(dmin, dmax);
+ gridType = "mmdg2_C";
+ xi = 0.75;
}
else
{
@@ -151,7 +177,7 @@ int main(int argc, char** argv)
}
}
- //create a grid from .dgf file
+ //create the grid from the .msh file
GridFactory gridFactory( "grids/" + gridType + "_" + std::to_string(dim)
+ "d.msh" );
const Grid& grid = *gridFactory.grid();
@@ -163,6 +189,7 @@ int main(int argc, char** argv)
const Mapper mapper(gridView, Dune::mcmgElementLayout());
const IMapper iMapper(iGridView, Dune::mcmgElementLayout());
+ //apply the coupled DG scheme
MMDG mmdg(gridView, mapper, iGridView, iMapper, *problem, solver);
mmdg(mu0, xi);