diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index d5b7c9b7365b91ea82567d57d50301c23b505e04..4d0cf39d6130d96e9ce9fc4c9e4c8a90bcffefbc 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -71,7 +71,7 @@ protected:
MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
MatrixRT G2_ {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
MatrixRT G3_ {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
- std::array<MatrixRT,3 > basisContainer_ = {G1_, G2_, G3_};
+ std::array<MatrixRT,3 > MatrixBasisContainer_ = {G1_, G2_, G3_};
Func2Tensor x3G_1_ = [] (const Domain& x) {
return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; //TODO könnte hier sign übergeben?
@@ -85,7 +85,7 @@ protected:
return MatrixRT{{0.0, (1.0/sqrt(2.0))*x[2], 0.0}, {(1.0/sqrt(2.0))*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
- const std::array<Func2Tensor, 3> x3MatrixBasis_ = {x3G_1_, x3G_2_, x3G_3_};
+ const std::array<Func2Tensor, 3> x3MatrixBasisContainer_ = {x3G_1_, x3G_2_, x3G_3_};
// --- Offset between basis indices
const int phiOffset_;
@@ -186,17 +186,29 @@ public:
shared_ptr<VectorCT> getLoad_alpha2(){return make_shared<VectorCT>(load_alpha2_);}
shared_ptr<VectorCT> getLoad_alpha3(){return make_shared<VectorCT>(load_alpha3_);}
- shared_ptr<VectorCT> getMu(){return make_shared<FuncScalar>(mu_);}
- shared_ptr<VectorCT> getLambda(){return make_shared<FuncScalar>(lambda_);}
+ shared_ptr<FuncScalar> getMu(){return make_shared<FuncScalar>(mu_);}
+ shared_ptr<FuncScalar> getLambda(){return make_shared<FuncScalar>(lambda_);}
// --- Get Correctors
- shared_ptr<VectorCT> getMcontainer(){return make_shared<VectorCT>(mContainer);}
+ // shared_ptr<VectorCT> getMcontainer(){return make_shared<VectorCT>(mContainer);}
+ // auto getMcontainer(){return make_shared<std::array<MatrixRT*, 3 > >(mContainer);}
+ auto getMcontainer(){return mContainer;}
shared_ptr<std::array<VectorCT, 3>> getPhicontainer(){return make_shared<std::array<VectorCT, 3>>(phiContainer);}
+
+
+ // shared_ptr<std::array<VectorRT, 3>> getBasiscontainer(){return make_shared<std::array<VectorRT, 3>>(basisContainer_);}
+ auto getMatrixBasiscontainer(){return make_shared<std::array<MatrixRT,3 >>(MatrixBasisContainer_);}
+ // auto getx3MatrixBasiscontainer(){return make_shared<std::array<Func2Tensor, 3>>(x3MatrixBasisContainer_);}
+ auto getx3MatrixBasiscontainer(){return x3MatrixBasisContainer_;}
+
+
+
+
// shared_ptr<VectorCT> getCorr_a(){return make_shared<VectorCT>(x_a_);}
- // shared_ptr<VectorCT> getCorr_K1(){return make_shared<VectorCT>(x_K1_);}
- // shared_ptr<VectorCT> getCorr_K2(){return make_shared<VectorCT>(x_K2_);}
- // shared_ptr<VectorCT> getCorr_K3(){return make_shared<VectorCT>(x_K3_);}
+ shared_ptr<VectorCT> getCorr_phi1(){return make_shared<VectorCT>(phi_1_);}
+ shared_ptr<VectorCT> getCorr_phi2(){return make_shared<VectorCT>(phi_2_);}
+ shared_ptr<VectorCT> getCorr_phi3(){return make_shared<VectorCT>(phi_3_);}
// Get the occupation pattern of the stiffness matrix
@@ -968,9 +980,9 @@ public:
for(size_t i=0; i<3; i++)
{
- M1_ += m_1_[i]*basisContainer_[i];
- M2_ += m_2_[i]*basisContainer_[i];
- M3_ += m_3_[i]*basisContainer_[i];
+ M1_ += m_1_[i]*MatrixBasisContainer_[i];
+ M2_ += m_2_[i]*MatrixBasisContainer_[i];
+ M3_ += m_3_[i]*MatrixBasisContainer_[i];
}
std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
@@ -1198,7 +1210,7 @@ public:
auto DerPhi1 = *phiDerContainer[a];
auto DerPhi2 = *phiDerContainer[b];
- auto matrixFieldGGVF = Dune::Functions::makeGridViewFunction(x3MatrixBasis_[a], basis_.gridView());
+ auto matrixFieldGGVF = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer_[a], basis_.gridView());
auto matrixFieldG = localFunction(matrixFieldGGVF);
// auto matrixFieldBGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
// auto matrixFieldB = localFunction(matrixFieldBGVF);
diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
index 12151857570b04bc46aeec54cec2a697a8e2387e..c3d6e20655da5991796f8209085cc386d54e8f6f 100644
--- a/dune/microstructure/EffectiveQuantitiesComputer.hh
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -7,6 +7,7 @@
#include <dune/microstructure/matrix_operations.hh>
#include <dune/microstructure/CorrectorComputer.hh>
+using namespace Dune;
using namespace MatrixOperations;
using std::shared_ptr;
using std::make_shared;
@@ -24,7 +25,6 @@ public:
static const int dimworld = 3;
// static const int nCells = 4;
-
static const int dim = Basis::GridView::dimension;
using Domain = typename CorrectorComputer<Basis>::Domain;
@@ -50,6 +50,7 @@ public:
FieldVector<double, dim> Bhat_; //effective loads induced by prestrain <LF_i, B>_L2
FieldVector<double, dim> Beff_; //effective strains Mb = ak
+
// corrector parts
VectorCT phi_E_TorusCV_; //phi_i * (a,K)_i
VectorCT phi_perp_TorusCV_;
@@ -86,13 +87,14 @@ public:
///////////////////////////////
// EffectiveQuantitiesComputer(CorrectorComputer<Basis>& correctorComputer, Func2Tensor prestrain)
// : correctorComputer_(correctorComputer), prestrain_(prestrain)
- EffectiveQuantitiesComputer(CorrectorComputer<Basis>& correctorComputer)
- : correctorComputer_(correctorComputer)
+ EffectiveQuantitiesComputer(CorrectorComputer<Basis>& correctorComputer, Func2Tensor prestrain)
+ : correctorComputer_(correctorComputer), prestrain_(prestrain)
{
// computePrestressLoadCV();
// computeEffectiveStrains();
-
+ // Q_ = 0;
+ // Q_ = {{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}};
// compute_phi_E_TorusCV();
// compute_phi_perp_TorusCV();
// compute_phi_TorusCV();
@@ -105,131 +107,9 @@ public:
}
-private:
-
- void computeQ()
- {
-
- auto test = correctorComputer_.getLoad_alpha1();
- auto phiContainer = correctorComputer_.getPhicontainer();
-
- // std::cout << mu_ << "mu_ can be used this way"
- // auto mu =
- // auto lamba =
- // auto gamma =
-
- return ;
- }
-
-
-
-
-
-
-
-
-
- // void computePrestressLoadCV()
- // {
- // cout << "computePrestressLoadCV ..." << endl;
- // auto cellAssembler = cellSolver_.getCellAssembler();
- // cellAssembler.assembleLoadVector(prestrain_, B_load_TorusCV_);
- // }
-
- // void computeEffectiveStrains()
- // {
- // cout << "Compute effective strains ..." << endl;
- // auto cellAssembler = cellSolver_.getCellAssembler();
-
- // //E_h = (U + K e_cs, 0, 0)
-
- // // lateral stretch strain E_U1 = e1 ot e1 = MatrixRT{{1, 0, 0}, {0, 0, 0}, {0, 0, 0}};
- // Func2Tensor E_a = [&] (const Domain& z) { return biotStrainApprox({1, 0, 0}, {0, 0, 0}, {0, z[dim-2], z[dim-1]}); };
-
- // // twist in x2-x3 plane E_K1 = (e1 x (0,x2,x3)^T ot e1) = MatrixRT{{0,-z[2], z[1]}, {0, 0 , 0 }, {0,0,0}};
- // Func2Tensor E_K1 = [&] (const Domain& z) { return biotStrainApprox({0, 0, 0}, {1, 0, 0}, {0, z[dim-2], z[dim-1]}); };
-
- // // bend strain in x1-x2 plane E_K2 = (e2 x (0,x2,x3)^T ot e1) = MatrixRT{{z[2], 0, 0}, {0, 0, 0}, {0, 0, 0}};
- // Func2Tensor E_K2 = [&] (const Domain& z) { return biotStrainApprox({0, 0, 0}, {0, 1, 0}, {0, z[dim-2], z[dim-1]}); };
-
- // // bend strain in x1-x3 plane E_K3 = (e3 x (0,x2,x3)^T ot e1) = MatrixRT{{-z[1], 0, 0}, {0, 0, 0}, {0, 0, 0}};
- // Func2Tensor E_K3 = [&] (const Domain& z) { return biotStrainApprox({0, 0, 0}, {0, 0, 1}, {0, z[dim-2], z[dim-1]}); };
-
-
- // shared_ptr<VectorCT> phi_Ei_TorusCV[nCells] = {cellSolver_.getCorr_a(),
- // cellSolver_.getCorr_K1(),
- // cellSolver_.getCorr_K1(),
- // cellSolver_.getCorr_K1()};
-
- // shared_ptr<VectorCT> loadsCell[nCells] = {cellAssembler.getLoad_a(),
- // cellAssembler.getLoad_K1(),
- // cellAssembler.getLoad_K2(),
- // cellAssembler.getLoad_K3()};
-
- // const std::array<shared_ptr<Func2Tensor>, nCells> strainsE = {make_shared<Func2Tensor>(E_a),
- // make_shared<Func2Tensor>(E_K1),
- // make_shared<Func2Tensor>(E_K2),
- // make_shared<Func2Tensor>(E_K3)};
-
-
- // for (unsigned int k = 0; k < nCells; k++)
- // for (unsigned int l = 0; l < nCells; l++){
- // std::cout << "<L Fi,Fj> ... " << k << l << std::endl;
- // M_[k][l] = cellAssembler.energyScalarProduct_FiFj(*strainsE[k], *phi_Ei_TorusCV[k], *strainsE[l], *phi_Ei_TorusCV[l]);
- // }
-
-
- // for (unsigned int k = 0; k < nCells; k++){
- // std::cout << "<L Fi,B> ... " << k << std::endl;
- // b_[k] = cellAssembler.energyScalarProduct_FiB(*strainsE[k], *phi_Ei_TorusCV[k], prestrain_);
- // }
-
-
- // M_.solve( aK_ ,b_);
- // std::cout << "\nstretch twist bend1 bend2\naK = "<< aK_ << std::endl;
- // }
-
- // void compute_phi_E_TorusCV()
- // {
- // cout << "compute_phi_E_TorusCV ..." << endl;
- // phi_E_TorusCV_ = *(cellSolver_.getCorr_a());
- // phi_E_TorusCV_ *= aK_[0];
- // auto tmp_xK1 = *(cellSolver_.getCorr_K1());
- // auto tmp_xK2 = *(cellSolver_.getCorr_K2());
- // auto tmp_xK3 = *(cellSolver_.getCorr_K3());
- // tmp_xK1 *= aK_[1];
- // tmp_xK2 *= aK_[2];
- // tmp_xK3 *= aK_[3];
- // phi_E_TorusCV_ += tmp_xK1;
- // phi_E_TorusCV_ += tmp_xK2;
- // phi_E_TorusCV_ += tmp_xK3;
- // }
-
- // void compute_phi_perp_TorusCV()
- // {
- // cout << "compute_phi_perp_TorusCV ..." << endl;
- // cellSolver_.solvePrestrainCorrector(phi_perp_TorusCV_, B_load_TorusCV_);
- // }
-
- // void compute_phi_TorusCV()
- // {
- // cout << "compute_phi_TorusCV ..." << endl;
- // phi_TorusCV_.resize(phi_E_TorusCV_.size());
- // for (int i=0; i < phi_TorusCV_.size(); i++)
- // phi_TorusCV_[i] = phi_E_TorusCV_[i] + phi_perp_TorusCV_[i];
- // }
-
-
-
-
-
-
-
-
- public:
- ///////////////////////////////
- // getter
- ///////////////////////////////
+ ///////////////////////////////
+ // getter
+ ///////////////////////////////
CorrectorComputer<Basis> getCorrectorComputer(){return correctorComputer_;}
const shared_ptr<Basis> getBasis()
@@ -237,17 +117,149 @@ private:
return correctorComputer_.getBasis();
}
+ auto getQeff(){return Q_;}
+ auto getBeff(){return Beff_;}
+ // -----------------------------------------------------------------
+ // --- Compute Effective Quantities
+ void computeEffectiveQuantities()
+ {
-
+ // Get everything.. better TODO: with Inheritance?
+ // auto test = correctorComputer_.getLoad_alpha1();
+ auto phiContainer = correctorComputer_.getPhicontainer();
+ auto MContainer = correctorComputer_.getMcontainer();
+ auto MatrixBasisContainer = correctorComputer_.getMatrixBasiscontainer();
+ auto x3MatrixBasisContainer = correctorComputer_.getx3MatrixBasiscontainer();
+ auto mu_ = *correctorComputer_.getMu();
+ auto lambda_ = *correctorComputer_.getLambda();
+ auto gamma = correctorComputer_.getGamma();
+ auto basis = *correctorComputer_.getBasis();
+
+ auto test = correctorComputer_.getCorr_phi1();
+
+ shared_ptr<VectorCT> phiBasis[3] = {correctorComputer_.getCorr_phi1(),
+ correctorComputer_.getCorr_phi2(),
+ correctorComputer_.getCorr_phi3()
+ };
+
+ auto prestrainGVF = Dune::Functions::makeGridViewFunction(prestrain_, basis.gridView());
+ auto prestrainFunctional = localFunction(prestrainGVF);
+
+ Q_ = 0 ;
+ Bhat_ = 0;
+
+ for(size_t a=0; a < 3; a++)
+ for(size_t b=0; b < 3; b++)
+ {
+ double energy = 0.0;
+ double prestrain = 0.0;
+ auto localView = basis.localView();
+ // auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiContainer[a]));
+ auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiBasis[a]));
+ // auto GVFunc_b = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,phiContainer[b]));
+ auto localfun_a = localFunction(GVFunc_a);
+ // auto localfun_b = localFunction(GVFunc_b);
+
+ ///////////////////////////////////////////////////////////////////////////////
+
+ auto matrixFieldG1GVF = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer[a], basis.gridView());
+ auto matrixFieldG1 = localFunction(matrixFieldG1GVF);
+ auto matrixFieldG2GVF = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer[b], basis.gridView());
+ auto matrixFieldG2 = localFunction(matrixFieldG2GVF);
+
+ auto muGridF = Dune::Functions::makeGridViewFunction(mu_, basis.gridView());
+ auto mu = localFunction(muGridF);
+ auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambda_, basis.gridView());
+ auto lambda= localFunction(lambdaGridF);
+
+ using GridView = typename Basis::GridView;
+
+ for (const auto& e : elements(basis.gridView()))
+ {
+ localView.bind(e);
+ matrixFieldG1.bind(e);
+ matrixFieldG2.bind(e);
+ localfun_a.bind(e);
+ // DerPhi2.bind(e);
+ mu.bind(e);
+ lambda.bind(e);
+ prestrainFunctional.bind(e);
+
+ double elementEnergy = 0.0;
+ double elementPrestrain = 0.0;
+
+ auto geometry = e.geometry();
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+
+ // int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ // int orderQR = 0;
+ // int orderQR = 1;
+ // int orderQR = 2;
+ // int orderQR = 3;
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
+
+ for (const auto& quadPoint : quad)
+ {
+ const auto& quadPos = quadPoint.position();
+ const double integrationElement = geometry.integrationElement(quadPos);
+
+ auto Chi1 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_a(quadPos))) + *MContainer[a];
+
+
+ auto G1 = matrixFieldG1(quadPos);
+ auto G2 = matrixFieldG2(quadPos);
+ // auto G1 = matrixFieldG1(e.geometry().global(quadPos)); //TEST
+ // auto G2 = matrixFieldG2(e.geometry().global(quadPos)); //TEST
+
+ auto X1 = G1 + Chi1;
+ // auto X2 = G2 + Chi2;
+
+
+ double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, G2);
+ elementEnergy += energyDensity * quadPoint.weight() * integrationElement; // quad[quadPoint].weight() ???
+ if (b==0)
+ {
+ elementPrestrain += linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, prestrainFunctional(quadPos));
+ }
+ }
+ energy += elementEnergy;
+ prestrain += elementPrestrain;
+
+
+ }
+ Q_[a][b] = energy;
+ if (b==0)
+ Bhat_[a] = prestrain;
+ }
+ printmatrix(std::cout, Q_, "Matrix Q_", "--");
+ printvector(std::cout, Bhat_, "Bhat_", "--");
+
+
+ ///////////////////////////////
+ // Compute effective Prestrain B_eff (by solving linear system)
+ //////////////////////////////
+ Q_.solve(Beff_,Bhat_);
+
+
+ //LOG-Output
+ auto& log = *(correctorComputer_.getLog());
+ log << "--- Prestrain Output --- " << std::endl;
+ log << "Bhat_: " << Bhat_ << std::endl;
+ log << "Beff_: " << Beff_ << " (Effective Prestrain)" << std::endl;
+ log << "------------------------ " << std::endl;
-
+ // TEST
+ // std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
+ return ;
+ }
-}; // end class
+}; // end class
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index c9a46000116651bf2f3f4ebd1bc041446e6d84f8..fea6d3c9fd58b11746b58aceadff55c96ef34615 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -9,8 +9,12 @@
# Choose Output-path for logfile
#############################################
### Remove/Comment this when running via Python-Script:
-outputPath=/home/stefan/DUNE/dune-microstructure/outputs
-#outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs
+#outputPath=/home/stefan/DUNE/dune-microstructure/outputs
+outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs
+
+
+### DEBUG Option:
+print_debug = true
#############################################
@@ -27,7 +31,7 @@ cellDomain=1
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
#----------------------------------------------------
-numLevels= 2 4
+numLevels= 2 3
#numLevels = 0 0 # computes all levels from first to second entry
#numLevels = 2 2 # computes all levels from first to second entry
#numLevels = 1 3 # computes all levels from first to second entry
@@ -83,8 +87,8 @@ theta=0.25
#--- choose composite-Implementation:
-geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries
-#geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/geometries
+#geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries
+geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/geometries
material_prestrain_imp= "material_neukamm" #(Python-indicator-function with same name determines material phases)
#material_prestrain_imp= "two_phase_material_2" #(Python-indicator-function with same name determines material phases)
#material_prestrain_imp= "two_phase_material_3" #(Python-indicator-function with same name determines material phases)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index f3f3b05f03424e97cdfdaca5742385a080db0069..42b95a4c3836f9c74ea51de3686a23ccbfc99df5 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -3,7 +3,8 @@
#set(CMAKE_BUILD_TYPE Debug)#
-set(programs Cell-Problem)
+set(programs Cell-Problem
+ Cell-Problem-New)
foreach(_program ${programs})
diff --git a/src/Cell-Problem-New.cc b/src/Cell-Problem-New.cc
index e565d5f49da746bdea7a0250b629adf3d226e2a6..4d583a128eaaaf3e7f65de7dfa098f9d1f3b1465 100644
--- a/src/Cell-Problem-New.cc
+++ b/src/Cell-Problem-New.cc
@@ -1718,20 +1718,78 @@ int main(int argc, char *argv[])
// define type of FE-Basis...
+
+ //Read from Parset...
+ int Phases = parameterSet.get<int>("Phases", 1);
+ std::string materialFunction = parameterSet.get<std::string>("materialFunction", "material");
+
+ Python::Module module = Python::import("material");
+ auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+ switch (Phases)
+ {
+ case 1: //homogeneous material
+ {
+ std::array<double,1> mu_ = parameterSet.get<std::array<double,1>>("mu", {1.0});
+ Python::Module module = Python::import(materialFunction);
+ auto indicatorFunction = Python::make_function<double>(module.get("f")); // get indicator function
+ auto muTerm = [mu_] (const Domain& x) {return mu_;};
+ break;
+ }
+ case 2:
+ {
+ std::array<double,2> mu_ = parameterSet.get<std::array<double,2>>("mu", {1.0,3.0});
+ Python::Module module = Python::import(materialFunction);
+ auto indicatorFunction = Python::make_function<double>(module.get("f")); // get indicator function
+ auto muTerm = [mu_,indicatorFunction] (const Domain& x)
+ {
+ if (indicatorFunction(x) == 1)
+ return mu_[0];
+ else
+ return mu_[1];
+ };
+ break;
+ }
+ // case 3:
+ // auto muTerm = [mu,indicatorFunction] (const Domain& x)
+ // {
+ // if (indicatorFunction(x) == 1)
+ // return mu[0];
+ // else if (indicatorFunction(x) == 2)
+ // return mu[1];
+ // else
+ // return mu[2];
+ // };
+ // break;
+ }
+
+
+
+
+
+
+
+
+
// typedef Dune::Functions::LagrangeBasis<GridView, order> FEBasis;
- // auto correctorComputer = CorrectorComputer(Basis_CE, muTerm, lambdaTerm, gamma, log, parameterSet);
+ auto correctorComputer = CorrectorComputer(Basis_CE, muTerm, lambdaTerm, gamma, log, parameterSet);
- // correctorComputer.solve();
- // correctorComputer.computeNorms();
- // correctorComputer.writeCorrectorsVTK(level);
- // correctorComputer.check_Orthogonality();
+ correctorComputer.solve();
+ correctorComputer.computeNorms();
+ correctorComputer.writeCorrectorsVTK(level);
+ correctorComputer.check_Orthogonality();
// // -------------------------------------------
- // auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer);
- // // effectiveModelComputer.getQuantities(Q_eff,B_eff);
- // effectiveQuantitiesComputer.computeQ();
+ auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,B_Term);
+ effectiveQuantitiesComputer.computeEffectiveQuantities();
+
+ //TEST
+ auto QT = effectiveQuantitiesComputer.getQeff();
+ auto Beff_ = effectiveQuantitiesComputer.getBeff();
+ printmatrix(std::cout, QT, "Matrix Q_T", "--");
+ printvector(std::cout, Beff_, "Beff", "--");
std::cout << "\n WORKED \n";
// break;