From 2de96c84ab5a38c1dd4ad3d9ebb68b1c6201f816 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Sun, 28 Aug 2022 17:08:18 +0200
Subject: [PATCH] Update
---
dune/microstructure/CorrectorComputer.hh | 57 ++++++--
dune/microstructure/prestrainedMaterial.hh | 151 +++++++++++----------
inputs/cellsolver.parset | 2 +-
src/Cell-Problem-New.cc | 47 ++++---
4 files changed, 152 insertions(+), 105 deletions(-)
diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 17809b83..323ddbc8 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -159,6 +159,7 @@ public:
shared_ptr<FuncScalar> getMu(){return make_shared<FuncScalar>(mu_);}
shared_ptr<FuncScalar> getLambda(){return make_shared<FuncScalar>(lambda_);}
+ shared_ptr<Material> getMaterial(){return make_shared<Material>(material_);}
// --- Get Correctors
// shared_ptr<VectorCT> getMcontainer(){return make_shared<VectorCT>(mContainer);}
@@ -265,6 +266,11 @@ public:
auto elasticityTensor = material_.getElasticityTensor();
+ auto indicatorFunction = material_.getIndicatorFunction();
+ // auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
+ // auto indicatorFunction = localFunction(indicatorFunctionGVF);
+ // indicatorFunction.bind(element);
+
// LocalBasis-Offset
const int localPhiOffset = localView.size();
@@ -342,11 +348,16 @@ public:
// auto etmp = material_.applyElasticityTensorLocal(defGradientU,quadPos);
// printmatrix(std::cout, etmp, "etmp", "--");
// 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 = 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 = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); // also works..
@@ -358,10 +369,15 @@ public:
// "m*phi" & "phi*m" - part
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(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
auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
elementMatrix[row][localPhiOffset+m] += value;
@@ -383,11 +399,15 @@ public:
// std::cout << "mu(quadPos): " << mu(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(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 = 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
elementMatrix[localPhiOffset+m][localPhiOffset+n] += energyDensityGG * quadPoint.weight() * integrationElement; // += !!!!! (Fixed-Bug)
@@ -423,6 +443,13 @@ public:
// 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
const auto& localFiniteElement= localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
@@ -493,10 +520,14 @@ public:
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 = 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(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV ); //TEST
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(element.geometry().global(quadPos)), defGradientV ); //TEST
@@ -508,9 +539,13 @@ public:
// "f*m"-part
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 = 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(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
@@ -1277,7 +1312,7 @@ public:
std::cout << "check_Orthogonality:" << "("<< a <<"," << b << "): " << energy << std::endl;
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;
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 274cf421..ac423af3 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -29,22 +29,33 @@ using std::make_shared;
// 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)
{
- FieldVector<double,3> mu = {1.0, 2.0, 3.0};
- FieldVector<double,3> lambda = {1.0 ,2.0 , 5.0};
-
- if (phase == 0)
- return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
- else
- return 2.0 * mu[1] * sym(G) + lambda[2] * trace(sym(G)) * Id();
-
+ const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
+ const 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();
}
+
template <class GridView> // needed for GridViewFunctions
class prestrainedMaterial
{
@@ -102,7 +113,8 @@ protected:
// FuncScalar 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_;
// VectorCT x_1_, x_2_, x_3_; // (all) Corrector coefficient vectors
@@ -148,7 +160,7 @@ public:
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"));
@@ -159,40 +171,45 @@ public:
// module.get("Phases").toC<int>(Phases_);
- auto indicatorFunction = Python::make_function<int>(module.get("indicatorFunction"));
- indicatorFunction_ = Dune::Functions::makeGridViewFunction(indicatorFunction , gridView_);
-
-
+ // auto indicatorFunction = Python::make_function<int>(module.get("indicatorFunction"));
+ // 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"));
// Func2Tensor materialFunction_ = Python::make_function<double>(module.get("f"));
// 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)
- // {
- // FieldVector<double,3> mu = {1.0, 2.0, 3.0};
- // FieldVector<double,3> lambda = {1.0 ,2.0 , 5.0};
-
- // if (phase == 0)
- // return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
- // else
- // return 2.0 * mu[1] * sym(G) + lambda[2] * trace(sym(G)) * Id();
-
- // }
+/////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+
+// static 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
+// }
+// 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
@@ -213,43 +230,37 @@ public:
{
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;
}
else
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
MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const
{
-
// Python::Module module = Python::import("material");
// auto indicatorFunctionTest = Python::make_function<double>(module.get("indicatorFunction"));
- //
// 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 = 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, elasticityTensor_(G,phase), "elasticityTensor_(G,phase)", "--");
-
- return tmp;
+ // return tmp;
// return material_1(G,indicatorFunction_(x));
-
}
@@ -267,19 +278,19 @@ public:
// }
- MatrixRT applyElasticityTensorLocal(const MatrixRT& G, const Domain& x) const
- {
- //--- 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 applyElasticityTensorLocal(const MatrixRT& G, const Domain& x) const
+ // {
+ // //--- 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}};
- if (indicatorFunction_(x) == 1)
- return L1_(G);
- else if (indicatorFunction_(x) == 2)
- return L2_(G);
- else
- return L3_(G);
+ // if (indicatorFunction_(x) == 1)
+ // return L1_(G);
+ // else if (indicatorFunction_(x) == 2)
+ // return L2_(G);
+ // else
+ // return L3_(G);
- }
+ // }
@@ -308,7 +319,9 @@ public:
// 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_);}
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 9f0192d8..4fac1d02 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -25,7 +25,7 @@ print_debug = true #(default=false)
## {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 = 2 2 # computes all levels from first to second entry
#numLevels = 1 3 # computes all levels from first to second entry
diff --git a/src/Cell-Problem-New.cc b/src/Cell-Problem-New.cc
index a98ecfbc..4b024380 100644
--- a/src/Cell-Problem-New.cc
+++ b/src/Cell-Problem-New.cc
@@ -346,11 +346,11 @@ int main(int argc, char *argv[])
// 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';
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"));
// auto indicatorFunction = Python::make_function<double>(module.get("indicatorFunction"));
@@ -368,7 +368,7 @@ int main(int argc, char *argv[])
// GridView::Element
// 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&)>;
// // MatrixFunc elasticityTensor = Python::make_function<MatrixRT>(module.get("L"));
@@ -397,28 +397,28 @@ int main(int argc, char *argv[])
- for (const auto& element : elements(Basis_CE.gridView()))
- {
- localindicatorFunction.bind(element);
- int orderQR = 2;
- const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
- for (const auto& quadPoint : quad)
- {
- const auto& quadPos = quadPoint.position();
+ // for (const auto& element : elements(Basis_CE.gridView()))
+ // {
+ // localindicatorFunction.bind(element);
+ // int orderQR = 2;
+ // const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+ // for (const auto& quadPoint : quad)
+ // {
+ // 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;
- // auto temp = TestTensor(G1_, element.geometry().global(quadPos));
- // auto temp2 = elasticityTensor_(G1_, element.geometry().global(quadPos));
- // std::cout << "material_.applyElasticityTensor:" << std::endl;
- auto tmp3 = material_.applyElasticityTensor(G1_, element.geometry().global(quadPos));
- // auto tmp3 = material_.applyElasticityTensor(G1_, quadPos);
- printmatrix(std::cout, tmp3, "tmp3", "--");
- }
- }
+ // // std::cout << "quadPos : " << quadPos << std::endl;
+ // // auto temp = TestTensor(G1_, element.geometry().global(quadPos));
+ // // auto temp2 = elasticityTensor_(G1_, element.geometry().global(quadPos));
+ // // std::cout << "material_.applyElasticityTensor:" << std::endl;
+ // // auto tmp3 = material_.applyElasticityTensor(G1_, element.geometry().global(quadPos));
+ // // auto tmp3 = material_.applyElasticityTensor(G1_, quadPos);
+ // // printmatrix(std::cout, tmp3, "tmp3", "--");
+ // }
+ // }
@@ -557,8 +557,7 @@ int main(int argc, char *argv[])
// };
- }
-/*
+
//------------------------------------------------------------------------------------------------
//--- compute Correctors
@@ -787,6 +786,6 @@ int main(int argc, char *argv[])
std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
-*/
+
}
--
GitLab