diff --git a/src/Compute_MuGamma.cc b/src/Compute_MuGamma.cc
index ad14fab74274f0ad10e19b4be13010f1f2eaf54e..3cef8ce93319663790cb25915e468b82f4450186 100644
--- a/src/Compute_MuGamma.cc
+++ b/src/Compute_MuGamma.cc
@@ -707,28 +707,34 @@ int main(int argc, char *argv[])
ParameterTree parameterSet;
if (argc < 2)
- ParameterTreeParser::readINITree("../../inputs/cellsolver.parset", parameterSet);
+ ParameterTreeParser::readINITree("../../inputs/computeMuGamma.parset", parameterSet);
else
{
ParameterTreeParser::readINITree(argv[1], parameterSet);
ParameterTreeParser::readOptions(argc, argv, parameterSet);
}
-
+ /////////////////////////////////
+ // SET OUTPUT
+ /////////////////////////////////
+ std::string outputPath = parameterSet.get("outputPath", "/home/klaus/Desktop/DUNE/dune-microstructure/outputs");
+ std::fstream log;
+ log.open(outputPath + "/outputMuGamma.txt" ,std::ios::out);
+ std::cout << "outputPath:" << outputPath << std::endl;
+
+
//////////////////////////////////
// Generate the grid
//////////////////////////////////
constexpr int dim = 2;
-// constexpr int dim = 3; //TEST
-
// QUAD-GRID
FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0});
FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0});
-// std::array<int, dim> nElements = {128 ,128};
- std::array<int, dim> nElements = {64,64};
-// std::array<int, dim> nElements = {32 ,32};
// std::array<int, dim> nElements = {16,16};
-// std::array<int, dim> nElements = {8,8};
+
+ std::array<int,2> nElements = parameterSet.get<std::array<int,2>>("nElements", {32,32});
+ std::cout << "Number of Elements in each direction: " << nElements << std::endl;
+ log << "Number of Elements in each direction: " << nElements << std::endl;
using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
@@ -738,45 +744,45 @@ int main(int argc, char *argv[])
using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
// ----------- INPUT PARAMETERS -----------------------------
+ std::string imp = parameterSet.get<std::string>("material_prestrain_imp", "analytical_Example");
+ log << "material_prestrain used: "<< imp << std::endl;
double gamma = parameterSet.get<double>("gamma", 1.0);
double theta = parameterSet.get<double>("theta", 1.0/4.0);
double mu1 = parameterSet.get<double>("mu1", 10.0);
double beta = parameterSet.get<double>("beta", 2.0);
double mu2 = beta*mu1;
-
std::cout << "Gamma:" << gamma << std::endl;
std::cout << "Theta:" << theta << std::endl;
std::cout << "mu1:" << mu1 << std::endl;
std::cout << "mu2:" << mu2 << std::endl;
std::cout << "beta:" << beta << std::endl;
+ log << "----- Input Parameters -----: " << std::endl;
+ log << "gamma: " << gamma << std::endl;
+ log << "theta: " << theta << std::endl;
+ log << "beta: " << beta << std::endl;
+ log << "material parameters: " << std::endl;
+ log << "mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
+
+// auto muTerm = [mu1, mu2, theta] (const Domain& z) {
+//
+// // if (abs(z[0]) > (theta/2.0))
+// if (abs(z[0]) >= (theta/2.0))
+// return mu1;
+// else
+// return mu2;
+// };
- auto muTerm = [mu1, mu2, theta] (const Domain& z) {
-
-// if (abs(z[0]) > (theta/2.0))
- if (abs(z[0]) >= (theta/2.0))
- return mu1;
- else
- return mu2;
- };
-
-
auto materialImp = IsotropicMaterialImp<dim>();
- auto muTermT = materialImp.getMu(parameterSet);
-
-// auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView);
-// auto muLocal = localFunction(muGridF);
- auto muGridF = Dune::Functions::makeGridViewFunction(muTermT, gridView);
+ auto muTerm = materialImp.getMu(parameterSet);
+ auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView);
auto muLocal = localFunction(muGridF);
/////////////////////////////////////////////////////////
// Stiffness matrix and right hand side vector
/////////////////////////////////////////////////////////
-// using Matrix = BCRSMatrix<double>;
-// using Vector = BlockVector<double>;
using Vector = BlockVector<FieldVector<double,1> >;
using Matrix = BCRSMatrix<FieldMatrix<double,1,1> >;
-
Matrix stiffnessMatrix;
Vector b;
@@ -795,14 +801,9 @@ int main(int argc, char *argv[])
auto basis = makeBasis(gridView, Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices)); // flatLexicographic()?
-
-// auto x3Func = [](const FieldVector<double,dim>& x){return x[1];}; //2D-Version
-// auto x3Func = [](const FieldVector<double,dim>& x){return x[2];}; // 3D-Version
-// auto forceTerm = [](const FieldVector<double,dim>& x){return 0.0;}; // gives q3= 2.08333 (close)...
-
auto forceTerm = [](const FieldVector<double,dim>& x){return x[1];}; //2D-Version
-// auto forceTermNEG = [](const FieldVector<double,dim>& x){return -1.0*x[1];}; //2D-Version
-// auto forceTerm = [](const FieldVector<double,dim>& x){return x[2];}; // 3D-Version
+ auto ForceGridF = Dune::Functions::makeGridViewFunction(forceTerm, gridView);
+ auto ForceLocal = localFunction(ForceGridF);
assemblePoissonProblem(basis, stiffnessMatrix, b, muLocal, forceTerm, gamma);
// printmatrix(std::cout, stiffnessMatrix, "StiffnessMatrix", "--");
@@ -810,18 +811,6 @@ int main(int argc, char *argv[])
- /////////////////////////////////////////////////////////////////////////////
- // Write matrix and load vector to files, to be used in later examples
- /////////////////////////////////////////////////////////////////////////////
-
-// // { matrix_rhs_writing_begin }
-// std::string baseName = "Dune-PeriodicPoisson";
-// //+ std::to_string(grid->maxLevel()) + "-refinements";
-// storeMatrixMarket(stiffnessMatrix, baseName + "-matrix.mtx");
-// storeMatrixMarket(b, baseName + "-rhs.mtx");
-// // { matrix_rhs_writing_end }
-//
-
@@ -829,27 +818,8 @@ int main(int argc, char *argv[])
// Compute solution
///////////////////////////
- // { algebraic_solving_begin }
- // Choose an initial iterate that fulfills the Dirichlet conditions
Vector x(basis.size());
x = b;
-
- // Turn the matrix into a linear operator
-// MatrixAdapter<Matrix,Vector,Vector> linearOperator(stiffnessMatrix);
-//
-// // Sequential incomplete LU decomposition as the preconditioner
-// SeqILU<Matrix,Vector,Vector> preconditioner(stiffnessMatrix, 1.0); // Relaxation factor
-// // Preconditioned conjugate gradient solver
-// CGSolver<Vector> cg(linearOperator,
-// preconditioner,
-// 1e-5, // Desired residual reduction factor
-// 50,
-// // Maximum number of iterations
-// 2);
-// // Verbosity of the solver
-// // Object storing some statistics about the solving process
-// InverseOperatorResult statistics;
-
std::cout << "------------ CG - Solver ------------" << std::endl;
MatrixAdapter<Matrix, Vector, Vector> op(stiffnessMatrix);
@@ -867,10 +837,7 @@ int main(int argc, char *argv[])
2); // verbosity of the solver
InverseOperatorResult statistics;
// Solve!
-// cg.apply(x, b, statistics);
solver.apply(x, b, statistics);
- // { algebraic_solving_end }
-
// std::cout << "------------ GMRES - Solver ------------" << std::endl;
// // Turn the matrix into a linear operator
@@ -896,65 +863,42 @@ int main(int argc, char *argv[])
// solver.apply(x, b, statistics);
//
+// -----------------------------------------------------------------------------------------------------
-
-
-
-
-
-
-
-
-// using SolutionRange = FieldVector<double,dim>;
-// using SolutionRange = FieldVector<double,1>;
using SolutionRange = double;
auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(basis, x);
-
+ // -------- PRINT OUTPUT --------
// printvector(std::cout, x, "coeffVector", "--" );
std::cout << "Gamma:" << gamma << std::endl;
-
- auto ForceGridF = Dune::Functions::makeGridViewFunction(forceTerm, gridView);
- auto ForceLocal = localFunction(ForceGridF);
-
-
-
-// double mu_gamma = computeMuGamma(basis, x, gamma, muLocal, x3Func);
-// double mu_gamma = computeMuGamma(basis, x, gamma, muLocal, forceTerm);
double mu_gamma = computeMuGamma(basis, x, gamma, muLocal, ForceLocal);
std::cout << "mu_gamma:" << mu_gamma << std::endl;
+ log << "----- OUTPUT: -----: " << std::endl;
+ log << "mu_gamma=" << mu_gamma << std::endl;
-
- std::cout.precision(10);
- std::cout << "mu_gamma:" << std::fixed << mu_gamma << std::endl;
-
+// std::cout.precision(10);
+// std::cout << "mu_gamma:" << std::fixed << mu_gamma << std::endl;
- // TEST "computeMuGamma"
-// Vector zeroVec;
-// zeroVec.resize(Basis_CE.size());
- Vector zeroVec(basis.size());
- zeroVec = 0;
-
- std::cout << "TEST computeMuGamma:" << computeMuGamma(basis, zeroVec, gamma, muLocal, ForceLocal)<< std::endl;
+
+// Vector zeroVec(basis.size());
+// zeroVec = 0;
+// std::cout << "TEST computeMuGamma:" << computeMuGamma(basis, zeroVec, gamma, muLocal, ForceLocal)<< std::endl;
std::cout << " --- print analytic solutions --- " << std::endl;
double q1 = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
-
-
- double hm = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0);
- double gm = mu1*((1-theta)+theta*beta) *(1.0/6.0);
+// double hm = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0);
+// double am = mu1*((1-theta)+theta*beta) *(1.0/6.0);
std::cout << "q1 : " << q1 << std::endl;
std::cout << "q2 : " << q2 << std::endl;
std::cout << "q3 should be between q1 and q2" << std::endl;
-
// std::cout << "hm : " << hm << std::endl;
-// std::cout << "gm : " << gm << std::endl;
+// std::cout << "am : " << am << std::endl;
if(mu_gamma > q2)