diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index e56f012592c9004b36386aab3d66eacbd0af5f91..b3a9cef8e7af4eedf45c268c146a8b9b3a8efcdb 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -43,16 +43,21 @@ public:
using MatrixCT = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >;
using ElementMatrixCT = Dune::Matrix<double>;
+
+
+ std::shared_ptr<Material> material_;
+
+ // double gamma_;
+
+
protected:
const Basis& basis_;
// const Material& material_;
- Material& material_;
+ // Material& material_;
// Material material_;
- double gamma_;
-
- fstream& log_; // Output-log
+ // fstream& log_; // Output-log
const Dune::ParameterTree& parameterSet_;
MatrixCT stiffnessMatrix_;
@@ -92,14 +97,29 @@ public:
///////////////////////////////
// constructor
///////////////////////////////
+ // CorrectorComputer(const Basis& basis,
+ // Material& material,
+ // std::fstream& log,
+ // const Dune::ParameterTree& parameterSet)
+ // : basis_(basis),
+ // material_(material),
+ // gamma_(material_.getGamma()),
+ // log_(log),
+ // parameterSet_(parameterSet),
+ // matrixBasis_(std::array<VoigtVector<double,3>,3>{matrixToVoigt(Dune::FieldMatrix<double,3,3>({{1, 0, 0}, {0, 0, 0}, {0, 0, 0}})),
+ // matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 0, 0}, {0, 1, 0}, {0, 0, 0}})),
+ // matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 1/std::sqrt(2.0), 0}, {1/std::sqrt(2.0), 0, 0}, {0, 0, 0}}))}),
+ // phiOffset_(basis.size())
+ // {}
+
CorrectorComputer(const Basis& basis,
- Material& material,
- std::fstream& log,
+ std::shared_ptr<Material> material,
+ // std::fstream& log,
const Dune::ParameterTree& parameterSet)
: basis_(basis),
material_(material),
- gamma_(material_.getGamma()),
- log_(log),
+ // gamma_(material_->getGamma()),
+ // log_(log),
parameterSet_(parameterSet),
matrixBasis_(std::array<VoigtVector<double,3>,3>{matrixToVoigt(Dune::FieldMatrix<double,3,3>({{1, 0, 0}, {0, 0, 0}, {0, 0, 0}})),
matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 0, 0}, {0, 1, 0}, {0, 0, 0}})),
@@ -124,11 +144,18 @@ public:
};
- void updateMaterial(Material mat)
+ // void updateMaterial(Material mat)
+ // {
+ // this->material_ = mat;
+ // }
+
+ void updateMaterial(std::shared_ptr<Material> mat)
{
- this->material_ = mat;
+ std::cout << "updateMaterial of CorrectorComputer" << std::endl;
+ material_ = mat;
}
+
///////////////////////////////
// Getter
///////////////////////////////
@@ -136,15 +163,18 @@ public:
Dune::ParameterTree getParameterSet() const {return parameterSet_;}
- fstream* getLog(){return &log_;}
+ // fstream* getLog(){return &log_;}
- double getGamma(){return gamma_;}
+ // double getGamma(){return gamma_;}
+ double getGamma(){return material_->gamma_;}
shared_ptr<MatrixCT> getStiffnessMatrix(){return make_shared<MatrixCT>(stiffnessMatrix_);}
shared_ptr<VectorCT> getLoad_alpha1(){return make_shared<VectorCT>(load_alpha1_);}
shared_ptr<VectorCT> getLoad_alpha2(){return make_shared<VectorCT>(load_alpha2_);}
shared_ptr<VectorCT> getLoad_alpha3(){return make_shared<VectorCT>(load_alpha3_);}
- shared_ptr<Material> getMaterial(){return make_shared<Material>(material_);}
+ // shared_ptr<Material> getMaterial(){return make_shared<Material>(material_);}
+ shared_ptr<Material> getMaterial(){return material_;}
+
// --- Get Correctors
auto getMcontainer(){return mContainer;}
@@ -278,7 +308,7 @@ public:
MatrixRT defGradientV(0);
defGradientV[k] = jacobians[i][0];
- deformationGradient[i][k] = symVoigt(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
+ deformationGradient[i][k] = symVoigt(crossSectionDirectionScaling((1.0/(material_->gamma_)),defGradientV));
}
}
@@ -294,7 +324,7 @@ public:
for (size_t k=0; k < dimworld; k++)
for (size_t i=0; i < nSf; i++)
{
- double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(deformationGradient[i][k],phase),deformationGradient[j][l]);
+ double energyDensity= voigtScalarProduct(material_->applyElasticityTensor(deformationGradient[i][k],phase),deformationGradient[j][l]);
size_t col = localView.tree().child(k).localIndex(i);
@@ -304,7 +334,7 @@ public:
// "m*phi" & "phi*m" - part
for( size_t m=0; m<3; m++)
{
- double energyDensityGphi = voigtScalarProduct(material_.applyElasticityTensor(matrixBasis_[m],phase),deformationGradient[j][l]);
+ double energyDensityGphi = voigtScalarProduct(material_->applyElasticityTensor(matrixBasis_[m],phase),deformationGradient[j][l]);
auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
elementMatrix[row][localPhiOffset+m] += value;
@@ -315,7 +345,7 @@ public:
for(size_t m=0; m<3; m++) //TODO ist symmetric.. reicht die hälfte zu berechnen!!!
for(size_t n=0; n<3; n++)
{
- double energyDensityGG = voigtScalarProduct(material_.applyElasticityTensor(matrixBasis_[m],phase),matrixBasis_[n]);
+ double energyDensityGG = voigtScalarProduct(material_->applyElasticityTensor(matrixBasis_[m],phase),matrixBasis_[n]);
elementMatrix[localPhiOffset+m][localPhiOffset+n] += energyDensityGG * quadPoint.weight() * integrationElement; // += !!!!! (Fixed-Bug)
@@ -346,7 +376,7 @@ public:
const auto element = localView.element();
const auto geometry = element.geometry();
- auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+ auto localIndicatorFunction = material_->getLocalIndicatorFunction();
localIndicatorFunction.bind(element);
// Set of shape functions for a single element
@@ -385,9 +415,9 @@ public:
tmpDefGradientV[k][1] = jacobians[i][0][1]; // X2
tmpDefGradientV[k][2] = jacobians[i][0][2]; // X3
- VoigtVector<double,3> defGradientV = symVoigt(crossSectionDirectionScaling((1.0/gamma_),tmpDefGradientV));
+ VoigtVector<double,3> defGradientV = symVoigt(crossSectionDirectionScaling((1.0/(material_->gamma_)),tmpDefGradientV));
- double energyDensity= voigtScalarProduct(material_.applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),defGradientV);
+ double energyDensity= voigtScalarProduct(material_->applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),defGradientV);
size_t row = localView.tree().child(k).localIndex(i);
elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
@@ -396,7 +426,7 @@ public:
// "f*m"-part
for (size_t m=0; m<3; m++)
{
- double energyDensityfG= voigtScalarProduct(material_.applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),matrixBasis_[m]);
+ double energyDensityfG= voigtScalarProduct(material_->applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),matrixBasis_[m]);
elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
}
@@ -440,7 +470,7 @@ public:
const auto& quadRule = Dune::QuadratureRules<double,dim>::rule(element.type(), orderQR);
// Determine the material phase at each quadrature point
- auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+ auto localIndicatorFunction = material_->getLocalIndicatorFunction();
localIndicatorFunction.bind(element);
std::vector<int> phaseAtQuadPoint(quadRule.size());
@@ -762,7 +792,7 @@ public:
solver.apply(x_2_, load_alpha2_, statistics);
std::cout << "solve linear system for x_3.\n";
solver.apply(x_3_, load_alpha3_, statistics);
- log_ << "Solver-type used: " <<" CG-Solver" << std::endl;
+ // log_ << "Solver-type used: " <<" CG-Solver" << std::endl;
std::cout << "statistics.converged " << statistics.converged << std::endl;
std::cout << "statistics.condition_estimate: " << statistics.condition_estimate << std::endl;
@@ -795,7 +825,7 @@ public:
solver.apply(x_1_, load_alpha1_, statistics); //load_alpha1 now contains the corresponding residual!!
solver.apply(x_2_, load_alpha2_, statistics);
solver.apply(x_3_, load_alpha3_, statistics);
- log_ << "Solver-type used: " <<" GMRES-Solver" << std::endl;
+ // 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;
@@ -805,7 +835,7 @@ public:
else if ( Solvertype==3)// QR - SOLVER
{
std::cout << "------------ QR - Solver ------------" << std::endl;
- log_ << "solveLinearSystems: We use QR solver.\n";
+ // log_ << "solveLinearSystems: We use QR solver.\n";
//TODO install suitesparse
Dune::SPQR<MatrixCT> sPQR(stiffnessMatrix_);
sPQR.setVerbosity(1);
@@ -826,14 +856,14 @@ public:
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;
+ // 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";
+ // log_ << "solveLinearSystems: We use UMFPACK solver.\n";
if(basis_.dimension() > 1e5)
std::cout << "WARNING: Using a direct solver with " << basis_.dimension() << " degrees of freedom may be not feasible or slow." << std::endl;
@@ -860,13 +890,13 @@ public:
// 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;
+ // log_ << "Solver-type used: " <<" UMFPACK-Solver" << std::endl;
}
////////////////////////////////////////////////////////////////////////////////////
else if (Solvertype==5)// CHOLDMOD - SOLVER
{
std::cout << "------------ CHOLMOD - Solver ------------" << std::endl;
- log_ << "solveLinearSystems: We use CHOLMOD solver.\n";
+ // log_ << "solveLinearSystems: We use CHOLMOD solver.\n";
if(basis_.dimension() > 1e5)
std::cout << "WARNING: Using a direct solver with " << basis_.dimension() << " degrees of freedom may be not feasible or slow." << std::endl;
@@ -893,7 +923,7 @@ public:
// 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: " <<" CHOLMOD-Solver" << std::endl;
+ // log_ << "Solver-type used: " <<" CHOLMOD-Solver" << std::endl;
}
std::cout << "--- Corrector computation finished ---" << std::endl;
std::cout << "Time required for solving:" << SolverTimer.elapsed() << std::endl;
@@ -936,14 +966,14 @@ public:
printmatrix(std::cout, mContainer[0], "Corrector-Matrix M_1", "--");
printmatrix(std::cout, mContainer[1], "Corrector-Matrix M_2", "--");
printmatrix(std::cout, mContainer[2], "Corrector-Matrix M_3", "--");
- log_ << "---------- OUTPUT ----------" << std::endl;
- log_ << " --------------------" << std::endl;
- log_ << "Corrector-Matrix M_1: \n" << mContainer[0] << std::endl;
- log_ << " --------------------" << std::endl;
- log_ << "Corrector-Matrix M_2: \n" << mContainer[1] << std::endl;
- log_ << " --------------------" << std::endl;
- log_ << "Corrector-Matrix M_3: \n" << mContainer[2] << std::endl;
- log_ << " --------------------" << std::endl;
+ // log_ << "---------- OUTPUT ----------" << std::endl;
+ // log_ << " --------------------" << std::endl;
+ // log_ << "Corrector-Matrix M_1: \n" << mContainer[0] << std::endl;
+ // log_ << " --------------------" << std::endl;
+ // log_ << "Corrector-Matrix M_2: \n" << mContainer[1] << std::endl;
+ // log_ << " --------------------" << std::endl;
+ // log_ << "Corrector-Matrix M_3: \n" << mContainer[2] << std::endl;
+ // log_ << " --------------------" << std::endl;
if(parameterSet_.get<bool>("write_IntegralMean", false))
@@ -981,24 +1011,24 @@ public:
/////////////////////////////////////////////////////////
// Write Solution (Corrector Coefficients) in Logs
/////////////////////////////////////////////////////////
- if(parameterSet_.get<bool>("write_corrector_phi1", false))
- {
- log_ << "\nSolution of Corrector problems:\n";
- log_ << "\n Corrector_phi1:\n";
- log_ << x_1_ << std::endl;
- }
- if(parameterSet_.get<bool>("write_corrector_phi2", false))
- {
- log_ << "-----------------------------------------------------";
- log_ << "\n Corrector_phi2:\n";
- log_ << x_2_ << std::endl;
- }
- if(parameterSet_.get<bool>("write_corrector_phi3", false))
- {
- log_ << "-----------------------------------------------------";
- log_ << "\n Corrector_phi3:\n";
- log_ << x_3_ << std::endl;
- }
+ // if(parameterSet_.get<bool>("write_corrector_phi1", false))
+ // {
+ // log_ << "\nSolution of Corrector problems:\n";
+ // log_ << "\n Corrector_phi1:\n";
+ // log_ << x_1_ << std::endl;
+ // }
+ // if(parameterSet_.get<bool>("write_corrector_phi2", false))
+ // {
+ // log_ << "-----------------------------------------------------";
+ // log_ << "\n Corrector_phi2:\n";
+ // log_ << x_2_ << std::endl;
+ // }
+ // if(parameterSet_.get<bool>("write_corrector_phi3", false))
+ // {
+ // log_ << "-----------------------------------------------------";
+ // log_ << "\n Corrector_phi3:\n";
+ // log_ << x_3_ << std::endl;
+ // }
}
@@ -1117,9 +1147,9 @@ public:
const auto& quadPos = quadPoint.position();
const auto integrationElement = geometry.integrationElement(quadPos);
- auto scaledSymGrad1 = sym(crossSectionDirectionScaling(1.0/gamma_, localfun_1(quadPos)));
- auto scaledSymGrad2 = sym(crossSectionDirectionScaling(1.0/gamma_, localfun_2(quadPos)));
- auto scaledSymGrad3 = sym(crossSectionDirectionScaling(1.0/gamma_, localfun_3(quadPos)));
+ auto scaledSymGrad1 = sym(crossSectionDirectionScaling(1.0/(material_->gamma_), localfun_1(quadPos)));
+ auto scaledSymGrad2 = sym(crossSectionDirectionScaling(1.0/(material_->gamma_), localfun_2(quadPos)));
+ auto scaledSymGrad3 = sym(crossSectionDirectionScaling(1.0/(material_->gamma_), localfun_3(quadPos)));
error_1 += scalarProduct(scaledSymGrad1,scaledSymGrad1) * quadPoint.weight() * integrationElement;
error_2 += scalarProduct(scaledSymGrad2,scaledSymGrad2) * quadPoint.weight() * integrationElement;
@@ -1153,7 +1183,7 @@ public:
auto localfun_3 = localFunction(GVFunc_3);
const std::array<decltype(localfun_1)*,3> phiDerContainer = {&localfun_1 , &localfun_2 , &localfun_3 };
- auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+ auto localIndicatorFunction = material_->getLocalIndicatorFunction();
for(int a=0; a<3; a++)
for(int b=0; b<3; b++)
@@ -1188,15 +1218,15 @@ public:
const auto& quadPos = quadPoint.position();
const double integrationElement = geometry.integrationElement(quadPos);
- auto Chi = sym(crossSectionDirectionScaling(1.0/gamma_, DerPhi2(quadPos))) + mContainer[b];
+ auto Chi = sym(crossSectionDirectionScaling(1.0/(material_->gamma_), DerPhi2(quadPos))) + mContainer[b];
- const auto strain1 = symVoigt(crossSectionDirectionScaling(1.0/gamma_, DerPhi1(quadPos)));
+ const auto strain1 = symVoigt(crossSectionDirectionScaling(1.0/(material_->gamma_), DerPhi1(quadPos)));
auto G = matrixFieldG(quadPos);
auto tmp = matrixToVoigt(G + mContainer[a]) + strain1;
- double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(tmp,localIndicatorFunction(quadPos)),symVoigt(Chi));
+ double energyDensity= voigtScalarProduct(material_->applyElasticityTensor(tmp,localIndicatorFunction(quadPos)),symVoigt(Chi));
elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
index f856d26c012d88c29adcd2e5fca6b5e085dfe59f..33b5c52669ae3a13dd40adc74b1a32caa25589cd 100644
--- a/dune/microstructure/EffectiveQuantitiesComputer.hh
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -32,8 +32,9 @@ public:
using MatrixPhase = std::function< MatrixRT(const int&) >;
protected:
- CorrectorComputer<Basis,Material>& correctorComputer_;
- const Material& material_;
+ // CorrectorComputer<Basis,Material>& correctorComputer_;
+ std::shared_ptr<CorrectorComputer<Basis,Material>> correctorComputer_;
+ // const Material& material_; //Get this from correctorComputer_ this is now an std::shared_ptr.
public:
@@ -45,24 +46,32 @@ public:
///////////////////////////////
// constructor
///////////////////////////////
- EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer,
- const Material& material)
- : correctorComputer_(correctorComputer),
- material_(material)
+ // EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer,
+ // const Material& material)
+ // : correctorComputer_(correctorComputer),
+ // material_(material)
+ // {}
+ // EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer)
+ // : correctorComputer_(correctorComputer)
+ // {}
+ EffectiveQuantitiesComputer(std::shared_ptr<CorrectorComputer<Basis,Material>> correctorComputer)
+ : correctorComputer_(correctorComputer)
{}
-
///////////////////////////////
// Getter
///////////////////////////////
- CorrectorComputer<Basis,Material> getCorrectorComputer(){return correctorComputer_;}
- const shared_ptr<Basis> getBasis() {return correctorComputer_.getBasis();}
+ // CorrectorComputer<Basis,Material> getCorrectorComputer(){return correctorComputer_;}
+ CorrectorComputer<Basis,Material> getCorrectorComputer(){return *correctorComputer_;}
+
+
+ const shared_ptr<Basis> getBasis() {return *correctorComputer_.getBasis();}
auto getQeff(){return Qeff_;}
auto getBeff(){return Beff_;}
- void updateCorrectorComputer(CorrectorComputer<Basis,Material> correctorComputer)
+ void updateCorrectorComputer(std::shared_ptr<CorrectorComputer<Basis,Material>> correctorComputer)
{
correctorComputer_ = correctorComputer;
}
@@ -74,20 +83,20 @@ public:
std::cout << "starting effective quantities computation..." << std::endl;
Dune::Timer effectiveQuantitiesTimer;
- auto MContainer = correctorComputer_.getMcontainer();
- auto MatrixBasisContainer = correctorComputer_.getMatrixBasiscontainer();
- auto x3MatrixBasisContainer = correctorComputer_.getx3MatrixBasiscontainer();
+ auto MContainer = correctorComputer_->getMcontainer();
+ auto MatrixBasisContainer = correctorComputer_->getMatrixBasiscontainer();
+ auto x3MatrixBasisContainer = correctorComputer_->getx3MatrixBasiscontainer();
- auto gamma = correctorComputer_.getGamma();
- auto basis = *correctorComputer_.getBasis();
- Dune::ParameterTree parameterSet = correctorComputer_.getParameterSet();
+ auto gamma = correctorComputer_->getGamma();
+ auto basis = *correctorComputer_->getBasis();
+ Dune::ParameterTree parameterSet = correctorComputer_->getParameterSet();
- shared_ptr<VectorCT> phiBasis[3] = {correctorComputer_.getCorr_phi1(),
- correctorComputer_.getCorr_phi2(),
- correctorComputer_.getCorr_phi3()
+ shared_ptr<VectorCT> phiBasis[3] = {correctorComputer_->getCorr_phi1(),
+ correctorComputer_->getCorr_phi2(),
+ correctorComputer_->getCorr_phi3()
};
- auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+ auto localIndicatorFunction = (correctorComputer_->material_)->getLocalIndicatorFunction();
Qeff_ = 0 ;
Bhat_ = 0;
@@ -138,13 +147,13 @@ public:
auto X1 = matrixToVoigt(G1 + Chi1);
- double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)), matrixToVoigt(sym(G2)));
+ double energyDensity= voigtScalarProduct((correctorComputer_->material_)->applyElasticityTensor(X1,localIndicatorFunction(quadPos)), matrixToVoigt(sym(G2)));
elementEnergy += energyDensity * quadPoint.weight() * integrationElement; // quad[quadPoint].weight() ???
if (b==0)
{
- auto prestrainPhaseValue = material_.prestrain(localIndicatorFunction(quadPos),element.geometry().global(quadPos));
- auto value = voigtScalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)),matrixToVoigt(sym(prestrainPhaseValue)));
+ auto prestrainPhaseValue = (correctorComputer_->material_)->prestrain(localIndicatorFunction(quadPos),element.geometry().global(quadPos));
+ auto value = voigtScalarProduct((correctorComputer_->material_)->applyElasticityTensor(X1,localIndicatorFunction(quadPos)),matrixToVoigt(sym(prestrainPhaseValue)));
elementPrestrain += value * quadPoint.weight() * integrationElement;
}
@@ -175,14 +184,14 @@ public:
printvector(std::cout, Beff_, "Beff_", "--");
//LOG-Output
- auto& log = *(correctorComputer_.getLog());
- log << "--- Effective moduli --- " << std::endl;
- log << "Qeff_: \n" << Qeff_ << std::endl;
- log << "------------------------ " << std::endl;
- log << "--- Prestrain Output --- " << std::endl;
- log << "Bhat_: " << Bhat_ << std::endl;
- log << "Beff_: " << Beff_ << " (Effective Prestrain)" << std::endl;
- log << "------------------------ " << std::endl;
+ // auto& log = *(correctorComputer_.getLog());
+ // log << "--- Effective moduli --- " << std::endl;
+ // log << "Qeff_: \n" << Qeff_ << std::endl;
+ // log << "------------------------ " << std::endl;
+ // log << "--- Prestrain Output --- " << std::endl;
+ // log << "Bhat_: " << Bhat_ << std::endl;
+ // log << "Beff_: " << Beff_ << " (Effective Prestrain)" << std::endl;
+ // log << "------------------------ " << std::endl;
// std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
std::cout << "Time required to solve for effective quantities: " << effectiveQuantitiesTimer.elapsed() << std::endl;
@@ -208,13 +217,13 @@ public:
const MatrixFunction& matrixFieldFuncB)
{
double energy = 0.0;
- auto mu_ = *correctorComputer_.getMu();
- auto lambda_ = *correctorComputer_.getLambda();
- auto gamma = correctorComputer_.getGamma();
- auto basis = *correctorComputer_.getBasis();
+ auto mu_ = *correctorComputer_->getMu();
+ auto lambda_ = *correctorComputer_->getLambda();
+ auto gamma = correctorComputer_->getGamma();
+ auto basis = *correctorComputer_->getBasis();
auto localView = basis.localView();
- auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+ auto localIndicatorFunction = (correctorComputer_->material_)->getLocalIndicatorFunction();
auto matrixFieldAGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncA, basis.gridView());
auto matrixFieldA = localFunction(matrixFieldAGVF);
@@ -240,7 +249,7 @@ public:
const auto& quadPos = quadPoint.position();
const double integrationElement = geometry.integrationElement(quadPos);
- double energyDensity= scalarProduct(material_.applyElasticityTensor(matrixFieldA(quadPos),localIndicatorFunction(quadPos)),sym(matrixFieldB(quadPos)));
+ double energyDensity= scalarProduct((correctorComputer_->material_)->applyElasticityTensor(matrixFieldA(quadPos),localIndicatorFunction(quadPos)),sym(matrixFieldB(quadPos)));
elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
}
diff --git a/dune/microstructure/microproblem.hh b/dune/microstructure/microproblem.hh
index 887303454e0a7dd808e71d7c0737ecfcc923392b..bfeee8a1d19e03cb8031664f1cf9c28f32427d67 100644
--- a/dune/microstructure/microproblem.hh
+++ b/dune/microstructure/microproblem.hh
@@ -51,6 +51,8 @@
// using namespace Dune;
// using namespace MatrixOperations;
+
+ // method to print types of objects:
template <class T>
constexpr std::string_view type_name()
{
@@ -98,8 +100,11 @@ class MicroProblem
std::fstream log_; // deprecated: only used to define CorrectorComputer object.
- MaterialType material_;
- CorrectorComputer<BasisType, MaterialType> correctorComputer_;
+ // MaterialType material_;
+ std::shared_ptr<MaterialType> material_;
+
+ // CorrectorComputer<BasisType, MaterialType> correctorComputer_;
+ std::shared_ptr<CorrectorComputer<BasisType, MaterialType> > correctorComputer_;
EffectiveQuantitiesComputer<BasisType,MaterialType> effectiveQuantitiesComputer_;
@@ -116,9 +121,12 @@ class MicroProblem
gridView_(grid_.leafGridView()),
basis_(createPeriodicBasis()),
// material_(setupMaterial(gridView_))
- material_(gridView_,microstructure_,parameterSet_,module_),
- correctorComputer_(basis_, material_, log_, parameterSet_),
- effectiveQuantitiesComputer_(correctorComputer_,material_) // Remove material dependency
+ // material_(gridView_,microstructure_,parameterSet_,module_),
+ material_(std::make_shared<MaterialType>(gridView_,microstructure_,parameterSet_,module_)),
+ correctorComputer_(std::make_shared<CorrectorComputer<BasisType, MaterialType>>(basis_, material_, parameterSet_)),
+ // correctorComputer_(basis_, material_, parameterSet_),
+ effectiveQuantitiesComputer_(correctorComputer_)
+ // effectiveQuantitiesComputer_(correctorComputer_,material_) // Remove material dependency
{};
// //Constructor
// MicroProblem(const Dune::ParameterTree& parameterSet,
@@ -198,9 +206,10 @@ class MicroProblem
{
microstructure_ = microstructure;
- material_.updateMicrostructure(microstructure_);
- correctorComputer_.updateMaterial(material_);
- // effectiveQuantitiesComputer_.updateCorrectorComputer(correctorComputer_);
+ // material_.updateMicrostructure(microstructure_);
+ material_->updateMicrostructure(microstructure_);
+ correctorComputer_->updateMaterial(material_);
+ effectiveQuantitiesComputer_.updateCorrectorComputer(correctorComputer_);
}
@@ -246,16 +255,16 @@ class MicroProblem
// std::cout << "Part1 works." << std::endl;
if (parameterSet_.get<bool>("write_materialFunctions", false))
- material_.writeVTKMaterialFunctions(parameterSet_.get<int>("gridLevel", 0));
+ material_->writeVTKMaterialFunctions(parameterSet_.get<int>("gridLevel", 0));
std::cout << "Material setup done" << std::endl;
// // Compute Correctors
// auto correctorComputer = CorrectorComputer(basis_, material, log, parameterSet_);
std::cout << "Starting Corrector assembly" << std::endl;
- correctorComputer_.assemble();
+ correctorComputer_->assemble();
std::cout << "Starting Corrector solve..." << std::endl;
- correctorComputer_.solve();
+ correctorComputer_->solve();
// std::cout << "Corrector done" << std::endl;
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index ba486153a45b38c9e9737a4fa3517435d7d6d8c1..edeaa835add91b6059b5b97eef78e629bf0c4602 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -119,6 +119,7 @@ public:
void updateMicrostructure(const Python::Reference microstructure)
{
+ std::cout << "updateMicrostrucrue of prestrainedMaterial" << std::endl;
microstructure_ = microstructure;
}
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index e8f4ea3f42218a84bb17689a415f204426d5f7a8..f619cbc8cfd0d59d42e5af64ae35b09593c30961 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -59,6 +59,9 @@
using namespace Dune;
using namespace MatrixOperations;
+
+using GridView = Dune::YaspGrid<3, Dune::EquidistantOffsetCoordinates<double, 3>>::LeafGridView;
+
//////////////////////////////////////////////////////////////////////
// Helper functions for Table-Output
//////////////////////////////////////////////////////////////////////
@@ -102,6 +105,7 @@ auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>
};
+
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int main(int argc, char *argv[])
@@ -222,44 +226,50 @@ int main(int argc, char *argv[])
//Setup for a constant microstructure.
Python::Reference microstructure = MicrostructureClass();
-
+
///////////////////////////////////
// Create prestrained material object
///////////////////////////////////
Dune::Timer materialSetupTimer;
+ using MaterialType = prestrainedMaterial<GridView>;
+
// auto material = prestrainedMaterial(gridView_CE,parameterSet,pyModule);
- auto material = prestrainedMaterial(gridView_CE,microstructure,parameterSet,pyModule);
+ // auto material = prestrainedMaterial(gridView_CE,microstructure,parameterSet,pyModule);
+ std::shared_ptr<MaterialType> material = std::make_shared<MaterialType>(gridView_CE,microstructure,parameterSet,pyModule);
std::cout << "Material setup took " << materialSetupTimer.elapsed() << " seconds " << std::endl;
// --- Get scale ratio
// double gamma = parameterSet.get<double>("gamma",1.0);
- std::cout << "scale ratio (gamma) set to : " << gamma << std::endl;
+ std::cout << "scale ratio (gamma) set to : " << material->gamma_ << std::endl;
//--- Compute Correctors
- auto correctorComputer = CorrectorComputer(Basis_CE, material, log, parameterSet);
- correctorComputer.assemble();
- correctorComputer.solve();
+ // auto correctorComputer = CorrectorComputer(Basis_CE, material, log, parameterSet);
+ // auto correctorComputer = CorrectorComputer(Basis_CE, materialPointer, log, parameterSet);
+ std::shared_ptr<CorrectorComputer<decltype(Basis_CE),MaterialType> > correctorComputer = std::make_shared<CorrectorComputer<decltype(Basis_CE),MaterialType>>(Basis_CE, material, parameterSet);
+ correctorComputer->assemble();
+ correctorComputer->solve();
//--- Check Correctors (options):
if(parameterSet.get<bool>("write_L2Error", false))
- correctorComputer.computeNorms();
+ correctorComputer->computeNorms();
if(parameterSet.get<bool>("write_VTK", false))
- correctorComputer.writeCorrectorsVTK(level);
+ correctorComputer->writeCorrectorsVTK(level);
//--- Additional Test: check orthogonality (75) from paper:
if(parameterSet.get<bool>("write_checkOrthogonality", false))
- correctorComputer.check_Orthogonality();
+ correctorComputer->check_Orthogonality();
//--- Check symmetry of stiffness matrix
if(print_debug)
- correctorComputer.checkSymmetry();
+ correctorComputer->checkSymmetry();
//--- Compute effective quantities
- auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,material);
+ // auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,material);
+ auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer);
effectiveQuantitiesComputer.computeEffectiveQuantities();
//--- Write material indicator function to VTK
if (parameterSet.get<bool>("write_materialFunctions", false))
- material.writeVTKMaterialFunctions(level);
+ material->writeVTKMaterialFunctions(level);
//--- Get effective quantities
auto Qeff = effectiveQuantitiesComputer.getQeff();
diff --git a/src/macro-problem.cc b/src/macro-problem.cc
index 32f6cb7b6216c4d4b82e7aaffeb97b31df57e4f7..fcb17c2cf69c00544beef5bc4677f27614c22945 100644
--- a/src/macro-problem.cc
+++ b/src/macro-problem.cc
@@ -67,25 +67,25 @@ const int dim = 2;
using namespace Dune;
- template <class T>
- constexpr std::string_view type_name()
- {
- using namespace std;
- #ifdef __clang__
- string_view p = __PRETTY_FUNCTION__;
- return string_view(p.data() + 34, p.size() - 34 - 1);
- #elif defined(__GNUC__)
- string_view p = __PRETTY_FUNCTION__;
- # if __cplusplus < 201402
- return string_view(p.data() + 36, p.size() - 36 - 1);
- # else
- return string_view(p.data() + 49, p.find(';', 49) - 49);
- # endif
- #elif defined(_MSC_VER)
- string_view p = __FUNCSIG__;
- return string_view(p.data() + 84, p.size() - 84 - 7);
- #endif
- }
+ // template <class T>
+ // constexpr std::string_view type_name()
+ // {
+ // using namespace std;
+ // #ifdef __clang__
+ // string_view p = __PRETTY_FUNCTION__;
+ // return string_view(p.data() + 34, p.size() - 34 - 1);
+ // #elif defined(__GNUC__)
+ // string_view p = __PRETTY_FUNCTION__;
+ // # if __cplusplus < 201402
+ // return string_view(p.data() + 36, p.size() - 36 - 1);
+ // # else
+ // return string_view(p.data() + 49, p.find(';', 49) - 49);
+ // # endif
+ // #elif defined(_MSC_VER)
+ // string_view p = __FUNCSIG__;
+ // return string_view(p.data() + 84, p.size() - 84 - 7);
+ // #endif
+ // }
int main(int argc, char *argv[])
diff --git a/test/readmicrostructuretest.cc b/test/readmicrostructuretest.cc
index 67e4da4155d82b4bf2b89417214ad5c1f46af3bf..0e728264ed1c943b9e2180e5fe12bedacadd7235 100644
--- a/test/readmicrostructuretest.cc
+++ b/test/readmicrostructuretest.cc
@@ -264,17 +264,25 @@ int main(int argc, char *argv[])
}
+
+ using MaterialType = prestrainedMaterial<GridView>;
+ std::shared_ptr<MaterialType> material = std::make_shared<MaterialType>(gridView_CE,microstructure,parameterSet,pyModule);
+
//Setup a material object
- auto material = prestrainedMaterial(gridView_CE,microstructure,parameterSet,pyModule);
+ // auto material = prestrainedMaterial(gridView_CE,microstructure,parameterSet,pyModule);
//Setup a corrector object & try to assemble/solve
- auto correctorComputer = CorrectorComputer(Basis_CE, material, log, parameterSet);
- correctorComputer.assemble();
- correctorComputer.solve();
+ // auto correctorComputer = CorrectorComputer(Basis_CE, material, log, parameterSet);
+ // correctorComputer.assemble();
+ // correctorComputer.solve();
+ std::shared_ptr<CorrectorComputer<decltype(Basis_CE),MaterialType> > correctorComputer = std::make_shared<CorrectorComputer<decltype(Basis_CE),MaterialType>>(Basis_CE, material, parameterSet);
+ correctorComputer->assemble();
+ correctorComputer->solve();
//Setup EffectiveQuantitiesComputer object & get Quantities.
- auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,material);
+ // auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,material);
+ auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer);
effectiveQuantitiesComputer.computeEffectiveQuantities();
auto Qeff = effectiveQuantitiesComputer.getQeff();
auto Beff = effectiveQuantitiesComputer.getBeff();