diff --git a/dune.module b/dune.module
index cbf9d40aef30499b5e199ea6d34cefd3b675fd8e..a9fb18ed5ba8e920393d55bcf7eaee6f3f11401f 100644
--- a/dune.module
+++ b/dune.module
@@ -7,6 +7,6 @@ Module: dune-microstructure
Version: 1.0
Maintainer: klaus.boehnlein@tu-dresden.de
# Required build dependencies
-Depends: dune-common dune-istl dune-geometry dune-localfunctions dune-uggrid dune-grid dune-typetree dune-functions dune-fufem
+Depends: dune-common dune-istl dune-geometry dune-localfunctions dune-uggrid dune-grid dune-typetree dune-functions dune-fufem dune-solvers
# Optional build dependencies
#Suggests:
diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
index 7efd263b9a8990bf23e5f569521d33894cb3ddc5..a4a563bcdb1373b1c4831832e058b15e801b247f 100644
--- a/dune/microstructure/matrix_operations.hh
+++ b/dune/microstructure/matrix_operations.hh
@@ -100,29 +100,29 @@ namespace MatrixOperations {
// return t1 + t2;
// }
- static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2) // CHANGED
- {
- auto t1 = 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
- return scalarProduct(t1,sym(E2));
- }
-//
// static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2) // CHANGED
// {
// auto t1 = 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
-// auto tmp1 = scalarProduct(t1,sym(E2));
-//
-// auto t2 = 2.0 * mu * sym(E2) + lambda * trace(sym(E2)) * Id();
-// auto tmp2 = scalarProduct(t2,sym(E1));
-//
-// if(abs(tmp1-tmp2) > 1e-8 )
-// {
-// std::cout << "linearizedStVenantKirchhoffDensity NOT Symmetric!" << std::endl;
-// std::cout << "tmp1: " << tmp1 << std::endl;
-// std::cout << "tmp2: " << tmp2 << std::endl;
-// }
-// return tmp1;
-//
+// return scalarProduct(t1,sym(E2));
// }
+// //
+ static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2) // CHANGED
+ {
+ auto t1 = 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
+ auto tmp1 = scalarProduct(t1,sym(E2));
+
+ auto t2 = 2.0 * mu * sym(E2) + lambda * trace(sym(E2)) * Id();
+ auto tmp2 = scalarProduct(t2,sym(E1));
+
+ if(abs(tmp1-tmp2) > 1e-8 )
+ {
+ std::cout << "linearizedStVenantKirchhoffDensity NOT Symmetric!" << std::endl;
+ std::cout << "tmp1: " << tmp1 << std::endl;
+ std::cout << "tmp2: " << tmp2 << std::endl;
+ }
+ return tmp1;
+
+ }
// --- Generalization: Define Quadratic QuadraticForm
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 0ff8aa2ac586a6f92af2586e0410bb1c22ba4ad9..8cb35c60ee4b2aa87e9abb98d234abdb2371ddc1 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -41,7 +41,7 @@ numLevels= 2 4
########################################################################################
# --- Choose scale ratio gamma:
-gamma=0.1
+gamma=1.0
#############################################
@@ -65,10 +65,10 @@ lambda1=5.0
# better: pass material parameters as a vector
# Poisson ratio nu = 1/4 => lambda = 0.4*mu
-mu=1.2 1.2 1
+#mu=1.2 1.2 1
#lambda=0.48 0.48 0.4
#rho=1.0 0
-#mu=80 80 60
+mu=80 80 60
lambda=80 80 25
@@ -124,9 +124,9 @@ set_IntegralZero = true
#############################################
-# Solver Type: #1: CG - SOLVER (default), #2: GMRES - SOLVER, #3: QR - SOLVER
+# Solver Type: #1: CG - SOLVER (default), #2: GMRES - SOLVER, #3: QR - SOLVER #4: UMFPACK - SOLVER
#############################################
-Solvertype = 1
+Solvertype = 3
Solver_verbosity = 0 #(default = 2) degree of information for solver output
@@ -137,6 +137,6 @@ Solver_verbosity = 0 #(default = 2) degree of information for solver output
#write_corrector_phi1 = false
#write_corrector_phi2 = false
#write_corrector_phi3 = false
-#write_corrector_phi1 = true
-#write_corrector_phi2 = true
-#write_corrector_phi3 = true
+write_corrector_phi1 = true
+write_corrector_phi2 = true
+write_corrector_phi3 = true
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index c8c4561a290df25949f7914d10a896175899901d..6fb3b2acb23b1b72b0cbf74fdaa492ad2eff6856 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -52,6 +52,8 @@
#include <dune/microstructure/matrix_operations.hh>
#include <dune/microstructure/vtk_filler.hh> //TEST
+#include <dune/solvers/solvers/umfpacksolver.hh> //TEST
+
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
@@ -418,7 +420,6 @@ void computeElementLoadVector( const LocalView& localView,
const double integrationElement = geometry.integrationElement(quadPos);
std::vector<FieldMatrix<double,1,dim> > referenceGradients;
-
localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
std::vector< FieldVector< double, dim> > gradients(referenceGradients.size());
@@ -1764,6 +1765,10 @@ int main(int argc, char *argv[])
solver.apply(x_2, load_alpha2, statistics);
solver.apply(x_3, load_alpha3, statistics);
log << "Solver-type used: " <<" GMRES-Solver" << std::endl;
+
+ std::cout << "statistics.converged " << statistics.converged << std::endl;
+ std::cout << "statistics.condition_estimate: " << statistics.condition_estimate << std::endl;
+ std::cout << "statistics.iterations: " << statistics.iterations << std::endl;
}
////////////////////////////////////////////////////////////////////////////////////
else if ( Solvertype == 3)// QR - SOLVER
@@ -1777,9 +1782,53 @@ int main(int argc, char *argv[])
InverseOperatorResult statisticsQR;
sPQR.apply(x_1, load_alpha1, statisticsQR);
+ std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
+ std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
+ std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
sPQR.apply(x_2, load_alpha2, statisticsQR);
+ std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
+ std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
+ std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
sPQR.apply(x_3, load_alpha3, statisticsQR);
+ std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
+ std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
+ std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
log << "Solver-type used: " <<" QR-Solver" << std::endl;
+
+
+ }
+ else if ( Solvertype == 4)// UMFPACK - SOLVER
+ {
+
+ std::cout << "------------ UMFPACK - Solver ------------" << std::endl;
+ log << "solveLinearSystems: We use UMFPACK solver.\n";
+
+ Dune::Solvers::UMFPackSolver<MatrixCT,VectorCT> solver;
+ solver.setProblem(stiffnessMatrix_CE,x_1,load_alpha1);
+// solver.preprocess();
+ solver.solve();
+ solver.setProblem(stiffnessMatrix_CE,x_2,load_alpha2);
+// solver.preprocess();
+ solver.solve();
+ solver.setProblem(stiffnessMatrix_CE,x_3,load_alpha3);
+// solver.preprocess();
+ solver.solve();
+
+// sPQR.apply(x_1, load_alpha1, statisticsQR);
+// std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
+// std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
+// std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
+// sPQR.apply(x_2, load_alpha2, statisticsQR);
+// std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
+// std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
+// std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
+// sPQR.apply(x_3, load_alpha3, statisticsQR);
+// std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
+// std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
+// std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
+ log << "Solver-type used: " <<" UMFPACK-Solver" << std::endl;
+
+
}
// printvector(std::cout, load_alpha1BS, "load_alpha1 before SOLVER", "--" );
// printvector(std::cout, load_alpha1, "load_alpha1 AFTER SOLVER", "--" );
@@ -1917,8 +1966,28 @@ int main(int argc, char *argv[])
const std::array<VectorCT, 3> loadContainer = {load_alpha1, load_alpha2, load_alpha3};
const std::array<MatrixRT*, 3 > mContainer = {&M_1 , &M_2, & M_3};
+ //TEST
+
+ auto zeroFunction = [] (const Domain& x) {
+ return MatrixRT{{0.0 , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ };
+ auto L2Norm_1 = computeL2Norm(Basis_CE,phi_1);
+ auto L2Norm_Symphi_1 = computeL2SymError(Basis_CE,phi_1,gamma,zeroFunction);
+ auto L2Norm_2 = computeL2Norm(Basis_CE,phi_2);
+ auto L2Norm_Symphi_2 = computeL2SymError(Basis_CE,phi_2,gamma,zeroFunction);
+ auto L2Norm_3 = computeL2Norm(Basis_CE,phi_3);
+ auto L2Norm_Symphi_3 = computeL2SymError(Basis_CE,phi_3,gamma,zeroFunction);
+ std::cout<< "L2Norm - Corrector 1: " << L2Norm_1 << std::endl;
+ std::cout<< "L2Norm (symgrad) - Corrector 1: " << L2Norm_Symphi_1 << std::endl;
+ std::cout<< "L2Norm - Corrector 2: " << L2Norm_2 << std::endl;
+ std::cout<< "L2Norm (symgrad) - Corrector 2: " << L2Norm_Symphi_2 << std::endl;
+ std::cout<< "L2Norm - Corrector 3: " << L2Norm_3 << std::endl;
+ std::cout<< "L2Norm (symgrad) - Corrector 3: " << L2Norm_Symphi_3 << std::endl;
+ std::cout<< "Frobenius-Norm of M_1: " << M_1.frobenius_norm() << std::endl;
+ std::cout<< "Frobenius-Norm of M_2: " << M_2.frobenius_norm() << std::endl;
+ std::cout<< "Frobenius-Norm of M_3: " << M_3.frobenius_norm() << std::endl;
//TEST
@@ -2096,10 +2165,10 @@ int main(int argc, char *argv[])
auto q2 = Q[1][1];
auto q3 = Q[2][2];
- std::cout<< "q1 : " << q1 << std::endl;
- std::cout<< "q2 : " << q2 << std::endl;
+// std::cout<< "q1 : " << q1 << std::endl;
+// std::cout<< "q2 : " << q2 << std::endl;
std::cout.precision(10);
- std::cout << "q3 : " << std::fixed << q3 << std::endl;
+// std::cout << "q3 : " << std::fixed << q3 << std::endl;
// std::cout<< "q3 : " << q3 << std::endl;
std::cout<< std::fixed << std::setprecision(6) << "q_onetwo=" << Q[0][1] << std::endl;
// std::cout<< std::fixed << std::setprecision(6) << "q_onetwo=" << Q[1][0] << std::endl; //TEST
@@ -2247,12 +2316,11 @@ int main(int argc, char *argv[])
return MatrixRT{{ (((mu1*mu2)/((theta*mu1 +(1-theta)*mu2)*muTerm(x))) - 1)*x[2] , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
- auto zeroFunction = [] (const Domain& x) {
- return MatrixRT{{0.0 , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
- };
+
if(write_L2Error)
{
+
// std::cout << " ----- L2ErrorSym ----" << std::endl;
auto L2SymError = computeL2SymError(Basis_CE,phi_1,gamma,symPhi_1_analytic);
// std::cout << "L2SymError: " << L2SymError << std::endl;
@@ -2269,7 +2337,7 @@ int main(int argc, char *argv[])
// std::cout << " -----------------" << std::endl;
// std::cout << " ----- L2NORM ----" << std::endl;
- auto L2Norm = computeL2Norm(Basis_CE,phi_1);
+ auto L2Norm = computeL2Norm(Basis_CE,phi_1);
// std::cout << "L2Norm(phi_1): " << L2Norm << std::endl;
// std::cout << " -----------------" << std::endl;