Skip to content
Snippets Groups Projects
Commit e271ba0c authored by Klaus Böhnlein's avatar Klaus Böhnlein
Browse files

Update

parent c81e9a36
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/microstructure/CorrectorComputer.hh
+ 46
− 11
View file @ e271ba0c
...@@ -159,6 +159,7 @@ public: ...@@ -159,6 +159,7 @@ public:
shared_ptr<FuncScalar> getMu(){return make_shared<FuncScalar>(mu_);} shared_ptr<FuncScalar> getMu(){return make_shared<FuncScalar>(mu_);}
shared_ptr<FuncScalar> getLambda(){return make_shared<FuncScalar>(lambda_);} shared_ptr<FuncScalar> getLambda(){return make_shared<FuncScalar>(lambda_);}
shared_ptr<Material> getMaterial(){return make_shared<Material>(material_);}
// --- Get Correctors // --- Get Correctors
// shared_ptr<VectorCT> getMcontainer(){return make_shared<VectorCT>(mContainer);} // shared_ptr<VectorCT> getMcontainer(){return make_shared<VectorCT>(mContainer);}
...@@ -265,6 +266,11 @@ public: ...@@ -265,6 +266,11 @@ public:
auto elasticityTensor = material_.getElasticityTensor(); auto elasticityTensor = material_.getElasticityTensor();
auto indicatorFunction = material_.getIndicatorFunction();
// auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
// auto indicatorFunction = localFunction(indicatorFunctionGVF);
// indicatorFunction.bind(element);
// LocalBasis-Offset // LocalBasis-Offset
const int localPhiOffset = localView.size(); const int localPhiOffset = localView.size();
...@@ -342,11 +348,16 @@ public: ...@@ -342,11 +348,16 @@ public:
// auto etmp = material_.applyElasticityTensorLocal(defGradientU,quadPos); // auto etmp = material_.applyElasticityTensorLocal(defGradientU,quadPos);
// printmatrix(std::cout, etmp, "etmp", "--"); // printmatrix(std::cout, etmp, "etmp", "--");
// double energyDensity= scalarProduct(etmp,defGradientV); // double energyDensity= scalarProduct(etmp,defGradientV);
double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),defGradientV); // double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),defGradientV);
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,quadPos),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensorLocal(defGradientU,quadPos),defGradientV); // double energyDensity= scalarProduct(material_.applyElasticityTensorLocal(defGradientU,quadPos),defGradientV);
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), defGradientU, defGradientV); //TEST // double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), defGradientU, defGradientV); //TEST
// double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); // also works.. // double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); // also works..
...@@ -358,10 +369,15 @@ public: ...@@ -358,10 +369,15 @@ public:
// "m*phi" & "phi*m" - part // "m*phi" & "phi*m" - part
for( size_t m=0; m<3; m++) for( size_t m=0; m<3; m++)
{ {
double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
// double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
// double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(quadPos)),sym(defGradientV));
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),sym(defGradientV));
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],quadPos),sym(defGradientV));
// double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV); // double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),defGradientV); // double energyDensityGphi= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),defGradientV);
// double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV); double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
// double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], defGradientV); //TEST // double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], defGradientV); //TEST
auto value = energyDensityGphi * quadPoint.weight() * integrationElement; auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
elementMatrix[row][localPhiOffset+m] += value; elementMatrix[row][localPhiOffset+m] += value;
...@@ -383,11 +399,15 @@ public: ...@@ -383,11 +399,15 @@ public:
// std::cout << "mu(quadPos): " << mu(quadPos) << std::endl; // std::cout << "mu(quadPos): " << mu(quadPos) << std::endl;
// std::cout << "lambda(quadPos): " << lambda(quadPos) << std::endl; // std::cout << "lambda(quadPos): " << lambda(quadPos) << std::endl;
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]); // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]);
double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]); // double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]);
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(quadPos)),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],quadPos),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),basisContainer[n]); // double energyDensityGG= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),basisContainer[n]);
// double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]); double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
// double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], basisContainer[n]); //TEST // double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], basisContainer[n]); //TEST
elementMatrix[localPhiOffset+m][localPhiOffset+n] += energyDensityGG * quadPoint.weight() * integrationElement; // += !!!!! (Fixed-Bug) elementMatrix[localPhiOffset+m][localPhiOffset+n] += energyDensityGG * quadPoint.weight() * integrationElement; // += !!!!! (Fixed-Bug)
...@@ -423,6 +443,13 @@ public: ...@@ -423,6 +443,13 @@ public:
// using MatrixRT = FieldMatrix< double, dimworld, dimworld>; // using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
auto elasticityTensor = material_.getElasticityTensor();
auto indicatorFunction = material_.getIndicatorFunction();
// // auto indicatorFunction = material_.getLocalIndicatorFunction();
// auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
// auto indicatorFunction = localFunction(indicatorFunctionGVF);
// indicatorFunction.bind(element);
// Set of shape functions for a single element // Set of shape functions for a single element
const auto& localFiniteElement= localView.tree().child(0).finiteElement(); const auto& localFiniteElement= localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size(); const auto nSf = localFiniteElement.localBasis().size();
...@@ -493,10 +520,14 @@ public: ...@@ -493,10 +520,14 @@ public:
defGradientV = crossSectionDirectionScaling((1.0/gamma_),defGradientV); defGradientV = crossSectionDirectionScaling((1.0/gamma_),defGradientV);
double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),defGradientV); // double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),defGradientV);
// double energyDensity= scalarProduct(material_.applyElasticityTensorLocal((-1.0)*forceTerm(quadPos),quadPos),defGradientV); // double energyDensity= scalarProduct(material_.applyElasticityTensorLocal((-1.0)*forceTerm(quadPos),quadPos),defGradientV);
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV ); double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV );
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)),forceTerm(quadPos), defGradientV ); //TEST // double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)),forceTerm(quadPos), defGradientV ); //TEST
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV ); //TEST // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV ); //TEST
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(element.geometry().global(quadPos)), defGradientV ); //TEST // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(element.geometry().global(quadPos)), defGradientV ); //TEST
...@@ -508,9 +539,13 @@ public: ...@@ -508,9 +539,13 @@ public:
// "f*m"-part // "f*m"-part
for (size_t m=0; m<3; m++) for (size_t m=0; m<3; m++)
{ {
double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),basisContainer[m]); // double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[m]));
// double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(basisContainer[m]));
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),sym(basisContainer[m]));
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),basisContainer[m]);
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),basisContainer[m]); // double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),basisContainer[m]);
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] ); double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] );
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), forceTerm(quadPos),basisContainer[m] ); //TEST // double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), forceTerm(quadPos),basisContainer[m] ); //TEST
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] ); //TEST // double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] ); //TEST
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), forceTerm(element.geometry().global(quadPos)),basisContainer[m] );//TEST // double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), forceTerm(element.geometry().global(quadPos)),basisContainer[m] );//TEST
...@@ -1277,7 +1312,7 @@ public: ...@@ -1277,7 +1312,7 @@ public:
std::cout << "check_Orthogonality:" << "("<< a <<"," << b << "): " << energy << std::endl; std::cout << "check_Orthogonality:" << "("<< a <<"," << b << "): " << energy << std::endl;
if(energy > 1e-1) if(energy > 1e-1)
std::cout << "WARNING: check Orthogonality! apparently (75) not satisfied on discrete level" << std::endl; std::cout << "WARNING: check_Orthogonality failed! apparently orthogonality (75) not satisfied on discrete level" << std::endl;
} }
return 0; return 0;
......
dune/microstructure/prestrainedMaterial.hh
+ 82
− 69
View file @ e271ba0c
...@@ -29,22 +29,33 @@ using std::make_shared; ...@@ -29,22 +29,33 @@ using std::make_shared;
// return !str[h] ? 5381 : (str2int(str, h+1) * 33) ^ str[h]; // return !str[h] ? 5381 : (str2int(str, h+1) * 33) ^ str[h];
// } // }
template<class Domain>
int indicatorFunction_material_1(const Domain& x)
{
double theta=0.25;
if (x[0] < (-1/2+theta) && x[2]<(-1/2+theta))
return 1; //#Phase1
else if (x[1]< (-1/2+theta) && x[2]>(1/2-theta))
return 2; //#Phase2
else
return 0; //#Phase3
}
MatrixRT material_1(const MatrixRT& G, const int& phase) MatrixRT material_1(const MatrixRT& G, const int& phase)
{ {
FieldVector<double,3> mu = {1.0, 2.0, 3.0}; const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
FieldVector<double,3> lambda = {1.0 ,2.0 , 5.0}; const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
if (phase == 1)
if (phase == 0) return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id(); if (phase == 2)
else return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id();
return 2.0 * mu[1] * sym(G) + lambda[2] * trace(sym(G)) * Id(); else
return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id();
} }
template <class GridView> // needed for GridViewFunctions template <class GridView> // needed for GridViewFunctions
class prestrainedMaterial class prestrainedMaterial
{ {
...@@ -102,7 +113,8 @@ protected: ...@@ -102,7 +113,8 @@ protected:
// FuncScalar indicatorFunction_; // FuncScalar indicatorFunction_;
// GridViewFunction<double(const Domain&), GridView> indicatorFunction_; // GridViewFunction<double(const Domain&), GridView> indicatorFunction_;
GridViewFunction<int(const Domain&), GridView> indicatorFunction_; GridViewFunction<int(const Domain&), GridView> indicatorFunctionGVF_;
Func2int indicatorFunction_;
// static const auto indicatorFunction_; // static const auto indicatorFunction_;
// VectorCT x_1_, x_2_, x_3_; // (all) Corrector coefficient vectors // VectorCT x_1_, x_2_, x_3_; // (all) Corrector coefficient vectors
...@@ -148,7 +160,7 @@ public: ...@@ -148,7 +160,7 @@ public:
std::cout << "materialFunctionName_ : " << materialFunctionName_ << std::endl; std::cout << "materialFunctionName_ : " << materialFunctionName_ << std::endl;
Python::Module module = Python::import(materialFunctionName_); // Python::Module module = Python::import(materialFunctionName_);
// elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L")); // elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
...@@ -159,40 +171,45 @@ public: ...@@ -159,40 +171,45 @@ public:
// module.get("Phases").toC<int>(Phases_); // module.get("Phases").toC<int>(Phases_);
auto indicatorFunction = Python::make_function<int>(module.get("indicatorFunction")); // auto indicatorFunction = Python::make_function<int>(module.get("indicatorFunction"));
indicatorFunction_ = Dune::Functions::makeGridViewFunction(indicatorFunction , gridView_); // indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
// indicatorFunction_ = Dune::Functions::makeGridViewFunction(indicatorFunction , gridView_);
L1_ = Python::make_function<MatrixRT>(module.get("L1"));
L2_ = Python::make_function<MatrixRT>(module.get("L2"));
L3_ = Python::make_function<MatrixRT>(module.get("L3"));
// indicatorFunction_ = localFunction(indicatorFunctionGVF);
// L1_ = Python::make_function<MatrixRT>(module.get("L1"));
// L2_ = Python::make_function<MatrixRT>(module.get("L2"));
// L3_ = Python::make_function<MatrixRT>(module.get("L3"));
// Func2TensorParam elasticityTensor_ = Python::make_function<double>(module.get("L")); // Func2TensorParam elasticityTensor_ = Python::make_function<double>(module.get("L"));
// Func2Tensor materialFunction_ = Python::make_function<double>(module.get("f")); // Func2Tensor materialFunction_ = Python::make_function<double>(module.get("f"));
// bool isotropic_ = true; // read from module File TODO // bool isotropic_ = true; // read from module File TODO
// Func2Tensor elasticityTensor_ = Python::make_function<double>(module.get("L"));
// Func2Tensor elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
} }
/////////////////////////////////////////////////////////////////////////////////////////////////////////
// static MatrixRT material_1(const MatrixRT& G, const int& phase) // static int indicatorFunction_material_1(const Domain& x)
// { // {
// FieldVector<double,3> mu = {1.0, 2.0, 3.0}; // double theta=0.25;
// FieldVector<double,3> lambda = {1.0 ,2.0 , 5.0}; // if (x[0] < (-1/2+theta) && x[2]<(-1/2+theta))
// return 1; //#Phase1
// if (phase == 0) // else if (x[1]< (-1/2+theta) && x[2]>(1/2-theta))
// return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id(); // return 2; //#Phase2
// else // else
// return 2.0 * mu[1] * sym(G) + lambda[2] * trace(sym(G)) * Id(); // return 0; //#Phase3
// }
// } // static MatrixRT material_1(const MatrixRT& G, const int& phase)
// {
// FieldVector<double,3> mu = {80.0, 80.0, 60.0};
// FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
// if (phase == 1)
// return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
// if (phase == 2)
// return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id();
// else
// return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id();
// }
//---function that determines elasticity Tensor //---function that determines elasticity Tensor
...@@ -213,43 +230,37 @@ public: ...@@ -213,43 +230,37 @@ public:
{ {
if(name == "material") if(name == "material")
{ {
// indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
// indicatorFunction_ = indicatorFunction_material_1;
indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1<Domain>, gridView_);
indicatorFunction_ = indicatorFunction_material_1<Domain>;
elasticityTensor_ = material_1; elasticityTensor_ = material_1;
} }
else else
DUNE_THROW(Exception, "There exists no material in materialDefinitions.hh with this name "); DUNE_THROW(Exception, "There exists no material in materialDefinitions.hh with this name ");
} }
// }
// void setupElasticityTensor()
// {
// if (materialFunctionName_=="material")
// {
// elasticityTensor_ = &material_1;
// std::cout << "typeid(elasticityTensor_).name() :" << typeid(elasticityTensor_).name() << '\n';
// // std::cout << " type_name<decltype(elasticityTensor_)>() " << type_name<decltype(elasticityTensor_)>() << '\n';
// }
// return;
// }
//--- apply elasticityTensor_ to input Matrix G at position x //--- apply elasticityTensor_ to input Matrix G at position x
MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const
{ {
// Python::Module module = Python::import("material"); // Python::Module module = Python::import("material");
// auto indicatorFunctionTest = Python::make_function<double>(module.get("indicatorFunction")); // auto indicatorFunctionTest = Python::make_function<double>(module.get("indicatorFunction"));
//
// return elasticityTensor_(G,indicatorFunctionTest(x)); // return elasticityTensor_(G,indicatorFunctionTest(x));
int phase = indicatorFunction_(x); // auto IFunction = localFunction(indicatorFunction_);
// auto IFunction = this->getIndicatorFunction();
// auto tmp = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
// auto tmpLocal = LocalF
// return elasticityTensor_(G,IFunction(x));
// return elasticityTensor_(G,indicatorFunction_(x));
return elasticityTensor_(G,indicatorFunctionGVF_(x));
// int phase = indicatorFunction_(x);
// auto tmp = elasticityTensor_(G,phase); // auto tmp = elasticityTensor_(G,phase);
auto tmp = material_1(G,indicatorFunction_(x)); // auto tmp = material_1(G,indicatorFunction_(x));
// printmatrix(std::cout, material_1(G,indicatorFunction_(x)), "material_1(G,indicatorFunction_(x))", "--"); // printmatrix(std::cout, material_1(G,indicatorFunction_(x)), "material_1(G,indicatorFunction_(x))", "--");
// printmatrix(std::cout, elasticityTensor_(G,phase), "elasticityTensor_(G,phase)", "--"); // printmatrix(std::cout, elasticityTensor_(G,phase), "elasticityTensor_(G,phase)", "--");
// return tmp;
return tmp;
// return material_1(G,indicatorFunction_(x)); // return material_1(G,indicatorFunction_(x));
} }
...@@ -267,19 +278,19 @@ public: ...@@ -267,19 +278,19 @@ public:
// } // }
MatrixRT applyElasticityTensorLocal(const MatrixRT& G, const Domain& x) const // MatrixRT applyElasticityTensorLocal(const MatrixRT& G, const Domain& x) const
{ // {
//--- apply elasticityTensor_ to input Matrix G at position x (local coordinates) // //--- apply elasticityTensor_ to input Matrix G at position x (local coordinates)
// MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}}; // // MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
if (indicatorFunction_(x) == 1) // if (indicatorFunction_(x) == 1)
return L1_(G); // return L1_(G);
else if (indicatorFunction_(x) == 2) // else if (indicatorFunction_(x) == 2)
return L2_(G); // return L2_(G);
else // else
return L3_(G); // return L3_(G);
} // }
...@@ -308,7 +319,9 @@ public: ...@@ -308,7 +319,9 @@ public:
// auto getIndicatorFunction() const {return indicatorFunction_;} // auto getIndicatorFunction() const {return indicatorFunction_;}
auto getIndicatorFunction() const {return localFunction(indicatorFunction_);} auto getLocalIndicatorFunction() const {return localFunction(indicatorFunctionGVF_);} //get as localFunction
auto getIndicatorFunctionGVF() const {return indicatorFunctionGVF_;} //get as GridViewFunction
auto getIndicatorFunction() const {return indicatorFunction_;}
// shared_ptr<Func2TensorParam> getElasticityTensor(){return make_shared<Func2TensorParam>(elasticityTensor_);} // shared_ptr<Func2TensorParam> getElasticityTensor(){return make_shared<Func2TensorParam>(elasticityTensor_);}
......
inputs/cellsolver.parset
+ 1
− 1
View file @ e271ba0c
...@@ -25,7 +25,7 @@ print_debug = true #(default=false) ...@@ -25,7 +25,7 @@ print_debug = true #(default=false)
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement ## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
#---------------------------------------------------- #----------------------------------------------------
numLevels= 2 2 numLevels= 2 3
#numLevels = 0 0 # computes all levels from first to second entry #numLevels = 0 0 # computes all levels from first to second entry
#numLevels = 2 2 # 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 #numLevels = 1 3 # computes all levels from first to second entry
......
src/Cell-Problem-New.cc
+ 23
− 24
View file @ e271ba0c
...@@ -346,11 +346,11 @@ int main(int argc, char *argv[]) ...@@ -346,11 +346,11 @@ int main(int argc, char *argv[])
// std::cout << "decltype(elasticityTensor_) " << decltype(elasticityTensor_) << std::endl; // std::cout << "decltype(elasticityTensor_) " << decltype(elasticityTensor_) << std::endl;
std::cout <<"typeid(elasticityTensor).name() :" << typeid(elasticityTensor_).name() << '\n'; std::cout << "typeid(elasticityTensor).name() :" << typeid(elasticityTensor_).name() << '\n';
std::cout << "typeid(TestTensor).name() :" << typeid(TestTensor).name() << '\n'; std::cout << "typeid(TestTensor).name() :" << typeid(TestTensor).name() << '\n';
using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >; using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
std::cout << "Import NOW:" << std::endl; // std::cout << "Import NOW:" << std::endl;
// MatrixFunc symTest = Python::make_function<MatrixRT>(module.get("sym")); // MatrixFunc symTest = Python::make_function<MatrixRT>(module.get("sym"));
// auto indicatorFunction = Python::make_function<double>(module.get("indicatorFunction")); // auto indicatorFunction = Python::make_function<double>(module.get("indicatorFunction"));
...@@ -368,7 +368,7 @@ int main(int argc, char *argv[]) ...@@ -368,7 +368,7 @@ int main(int argc, char *argv[])
// GridView::Element // GridView::Element
// auto localindicatorFunction = localFunction(indicatorFunction); // auto localindicatorFunction = localFunction(indicatorFunction);
std::cout << "typeid(localindicatorFunction).name() :" << typeid(localindicatorFunction).name() << '\n'; // std::cout << "typeid(localindicatorFunction).name() :" << typeid(localindicatorFunction).name() << '\n';
// using MatrixDomainFunc = std::function< MatrixRT(const MatrixRT&,const Domain&)>; // using MatrixDomainFunc = std::function< MatrixRT(const MatrixRT&,const Domain&)>;
// // MatrixFunc elasticityTensor = Python::make_function<MatrixRT>(module.get("L")); // // MatrixFunc elasticityTensor = Python::make_function<MatrixRT>(module.get("L"));
...@@ -397,28 +397,28 @@ int main(int argc, char *argv[]) ...@@ -397,28 +397,28 @@ int main(int argc, char *argv[])
for (const auto& element : elements(Basis_CE.gridView())) // for (const auto& element : elements(Basis_CE.gridView()))
{ // {
localindicatorFunction.bind(element); // localindicatorFunction.bind(element);
int orderQR = 2; // int orderQR = 2;
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR); // const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
for (const auto& quadPoint : quad) // for (const auto& quadPoint : quad)
{ // {
const auto& quadPos = quadPoint.position(); // const auto& quadPos = quadPoint.position();
// std::cout << "localindicatorFunction(quadPos): " << localindicatorFunction(quadPos) << std::endl; // // std::cout << "localindicatorFunction(quadPos): " << localindicatorFunction(quadPos) << std::endl;
// std::cout << "quadPos : " << quadPos << std::endl; // // std::cout << "quadPos : " << quadPos << std::endl;
// auto temp = TestTensor(G1_, element.geometry().global(quadPos)); // // auto temp = TestTensor(G1_, element.geometry().global(quadPos));
// auto temp2 = elasticityTensor_(G1_, element.geometry().global(quadPos)); // // auto temp2 = elasticityTensor_(G1_, element.geometry().global(quadPos));
// std::cout << "material_.applyElasticityTensor:" << std::endl; // // std::cout << "material_.applyElasticityTensor:" << std::endl;
auto tmp3 = material_.applyElasticityTensor(G1_, element.geometry().global(quadPos)); // // auto tmp3 = material_.applyElasticityTensor(G1_, element.geometry().global(quadPos));
// auto tmp3 = material_.applyElasticityTensor(G1_, quadPos); // // auto tmp3 = material_.applyElasticityTensor(G1_, quadPos);
printmatrix(std::cout, tmp3, "tmp3", "--"); // // printmatrix(std::cout, tmp3, "tmp3", "--");
} // }
} // }
...@@ -557,8 +557,7 @@ int main(int argc, char *argv[]) ...@@ -557,8 +557,7 @@ int main(int argc, char *argv[])
// }; // };
}
/*
//------------------------------------------------------------------------------------------------ //------------------------------------------------------------------------------------------------
//--- compute Correctors //--- compute Correctors
...@@ -787,6 +786,6 @@ int main(int argc, char *argv[]) ...@@ -787,6 +786,6 @@ int main(int argc, char *argv[])
std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl; std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
*/
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment