From 0b684c69f1daabf7bd10ca55f8dcca2895c4d618 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Mon, 9 Aug 2021 20:52:23 +0200
Subject: [PATCH] backup
---
inputs/cellsolver.parset | 4 +-
src/dune-microstructure.cc | 145 ++++++++++++++++---------------------
2 files changed, 64 insertions(+), 85 deletions(-)
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 06b3cff7..45881378 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -21,7 +21,7 @@ cellDomain = 1
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
########################################################################
-numLevels = 1 3 # computes all levels from first to second entry
+numLevels = 1 2 # computes all levels from first to second entry
#Elements_Cell = 20 20 20 # number elements in each direction (y1 y2 x3)
@@ -33,7 +33,7 @@ numLevels = 1 3 # computes all levels from first to second entry
#nElements_Cell = 10 10 10
#nElements_Cell = 2 2 2
-nElements_Cell = 4 4 4
+#nElements_Cell = 4 4 4
#nElements_Cell = 8 8 8
#nElements_Cell = 16 16 16
#nElements_Cell = 32 32 32
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 0ed65a79..de77f882 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -91,9 +91,6 @@ std::string prd(const type x, const int decDigits, const int width) {
-
-
-
template<class Basis>
auto arbitraryComponentwiseIndices(const Basis& basis,
const int elementNumber,
@@ -117,9 +114,9 @@ auto arbitraryComponentwiseIndices(const Basis& basis,
for (int k = 0; k < 3; k++)
{
auto rowLocal = localView.tree().child(k).localIndex(leafIdx); //Input: LeafIndex! TODO braucht hier (Inverse ) Local-to-Leaf Idx Map
- std::cout << "rowLocal:" << rowLocal << std::endl;
+// std::cout << "rowLocal:" << rowLocal << std::endl;
row[k] = localView.index(rowLocal);
- std::cout << "row[k]:" << row[k] << std::endl;
+// std::cout << "row[k]:" << row[k] << std::endl;
}
}
elementCounter++;
@@ -234,7 +231,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
- // strainV[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // same ? genügt strainU
+// strainV[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // same ? genügt strainU
// printmatrix(std::cout, strainU , "strainU", "--");
}
@@ -265,7 +262,6 @@ void computeElementStiffnessMatrix(const LocalView& localView,
size_t row = localView.tree().child(l).localIndex(j);
elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
-
}
}
@@ -275,8 +271,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
for (size_t l=0; l< nCompo; l++)
for (size_t j=0; j < nSf; j++ )
{
-
-
+
// (scaled) Deformation gradient of the test basis function
MatrixRT defGradientV(0);
defGradientV[l][0] = gradients[j][0]; // Y
@@ -291,7 +286,6 @@ void computeElementStiffnessMatrix(const LocalView& localView,
strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
}
-
for( size_t m=0; m<3; m++)
{
@@ -323,63 +317,10 @@ void computeElementStiffnessMatrix(const LocalView& localView,
elementMatrix[row][localPhiOffset+m] += value;
elementMatrix[localPhiOffset+m][row] += value;
-
}
}
-
-
- // "phi*m"-part
-// for (size_t k=0; k < nCompo; k++)
-// for (size_t i=0; i < nSf; i++)
-// {
-//
-// // (scaled) Deformation gradient of the ansatz basis function
-// MatrixRT defGradientU(0);
-// defGradientU[k][0] = gradients[i][0]; // Y
-// defGradientU[k][1] = gradients[i][1]; //X2
-// defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
-//
-//
-// // symmetric Gradient (Elastic Strains)
-// MatrixRT strainU;
-// for (int ii=0; ii<nCompo; ii++)
-// for (int jj=0; jj<nCompo; jj++)
-// {
-// strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
-// }
-//
-// // St.Venant-Kirchhoff stress
-// MatrixRT stressU(0);
-// stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
-//
-// double trace = 0;
-// for (int ii=0; ii<nCompo; ii++)
-// trace += strainU[ii][ii];
-//
-// for (int ii=0; ii<nCompo; ii++)
-// stressU[ii][ii] += lambda(quadPos) * trace;
-//
-// for( size_t n=0; n<3; n++)
-// {
-//
-// // < L sym[D_gamma*nabla phi_i], G_j >
-// double energyDensityPhiG = 0;
-// for (int ii=0; ii<nCompo; ii++)
-// for (int jj=0; jj<nCompo; jj++)
-// energyDensityPhiG += stressU[ii][jj] * basisContainer[n][ii][jj]; // "phi*m"-part
-//
-// size_t col = localView.tree().child(k).localIndex(i);
-//
-// elementMatrix[localPhiOffset+n][col] += energyDensityPhiG * quadPoint.weight() * integrationElement;
-//
-// }
-// }
-//
-
-
-
// "m*m"-part
for(size_t m=0; m<3; m++)
for(size_t n=0; n<3; n++)
@@ -405,10 +346,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
elementMatrix[localPhiOffset+m][localPhiOffset+n]= energyDensityGG * quadPoint.weight() * integrationElement;
}
-
}
-
-
}
@@ -456,7 +394,7 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet&
for (size_t i=0; i<3; i++ )
for (size_t j=0; j<3; j++ )
{
- nb.add(phiOffset+i,phiOffset+j); // m*m part of StiffnessMatrix
+ nb.add(phiOffset+i,phiOffset+j); // m*m part of StiffnessMatrix
}
}
//////////////////////////////////////////////////////////////////
@@ -467,11 +405,6 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet&
unsigned int arbitraryLeafIndex = parameterSet.template get<unsigned int>("arbitraryLeafIndex", 0);
unsigned int arbitraryElementNumber = parameterSet.template get<unsigned int>("arbitraryElementNumber", 0);
-
-
-// assert(arbitraryLeafIndex < localView.size() ); // "localIndex is larger than #shapeFunctions" TODO ERROR
-
-
const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // macht keinen Unterschied ob hier k oder 0..
const auto nSf = localFiniteElement.localBasis().size();
assert(arbitraryLeafIndex < nSf );
@@ -526,7 +459,6 @@ void computeElementLoadVector( const LocalView& localView,
elementRhs.resize(localView.size() +3);
elementRhs = 0;
-
// LocalBasis-Offset
const int localPhiOffset = localView.size();
@@ -546,7 +478,6 @@ void computeElementLoadVector( const LocalView& localView,
for (const auto& quadPoint : quad)
{
-
const FieldVector<double,dim>& quadPos = quadPoint.position();
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
const double integrationElement = geometry.integrationElement(quadPos);
@@ -2885,6 +2816,7 @@ int main(int argc, char *argv[])
auto q1_c = Q[0][0];
auto q2_c = Q[1][1];
auto q3_c = Q[2][2];
+
std::cout<< "q1_c: " << q1_c << std::endl;
std::cout<< "q2_c: " << q2_c << std::endl;
std::cout<< "q3_c: " << q3_c << std::endl;
@@ -2923,7 +2855,7 @@ int main(int argc, char *argv[])
double b3_hat = 0.0;
double b1_eff = (-3.0/2.0)*theta;
- double b2_eff = (-3.0*theta*mu2)/(mu1*(1-theta)+mu2*theta);
+ double b2_eff = (-3.0/2.0)*(theta*mu2)/(mu1*(1-theta)+mu2*theta);
double b3_eff = 0.0;
@@ -2956,10 +2888,10 @@ int main(int argc, char *argv[])
Storage_Quantities.push_back(std::abs(q1 - q1_c));
Storage_Quantities.push_back(std::abs(q2 - q2_c));
- Storage_Quantities.push_back(std::abs(q3 - q3_c));
- Storage_Quantities.push_back(std::abs(b1_eff - b1eff_c));
- Storage_Quantities.push_back(std::abs(b2_eff - b2eff_c));
- Storage_Quantities.push_back(std::abs(b3_eff - b3eff_c));
+ Storage_Quantities.push_back( q3_c );
+ Storage_Quantities.push_back(std::abs(b1_eff - Beff[0]));
+ Storage_Quantities.push_back(std::abs(b2_eff - Beff[1]));
+ Storage_Quantities.push_back(std::abs(b3_eff - Beff[2]));
auto symPhi_1_analytic = [mu1, mu2, theta, muTerm] (const Domain& x) {
return MatrixRT{{ (((mu1*mu2)/((theta*mu1 +(1-theta)*mu2)*muTerm(x))) - 1)*x[2] , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
@@ -3076,12 +3008,12 @@ int main(int argc, char *argv[])
vtkWriter.write("CellProblem-result");
std::cout << "wrote data to file: CellProblem-result" << std::endl; // better with string for output name..
- log.close();
+
levelCounter++; //TODO
- }
+ } // Level-Loop End
/////////////////////////////////////////
// Print Storage
@@ -3135,6 +3067,13 @@ int main(int argc, char *argv[])
<< center("L2SNorm_ana",tableWidth) << " | "
<< center("L2Norm",tableWidth) << " | " << "\n";
std::cout << std::string(tableWidth*6 + 3*6, '-') << "\n";
+ log << center("Levels",tableWidth) << " | "
+ << center("L2SymError",tableWidth) << " | "
+ << center("Order",tableWidth) << " | "
+ << center("L2SymNorm",tableWidth) << " | "
+ << center("L2SNorm_ana",tableWidth) << " | "
+ << center("L2Norm",tableWidth) << " | " << "\n";
+ log << std::string(tableWidth*6 + 3*6, '-') << "\n";
// std::vector<double> tmpVec; // oder std::variant
@@ -3150,8 +3089,9 @@ int main(int argc, char *argv[])
// std::visit([](auto&& arg){std::cout << std::left << std::setw(12) << std::setfill(' ') << arg;}, v);
// std::visit([tableWidth](auto&& arg){std::cout << center(prd(arg,5,12),tableWidth) << " | ";}, v);
- std::visit([tableWidth](auto&& arg){std::cout << center(prd(arg,5,1),tableWidth) << " | ";}, v);
-
+ std::visit([tableWidth](auto&& arg){std::cout << center(prd(arg,5,1),tableWidth) << " | ";}, v); // Anzahl-Nachkommastellen
+ std::visit([tableWidth, &log](auto&& arg){log << center(prd(arg,5,1),tableWidth) << " & ";}, v);
+
@@ -3163,6 +3103,7 @@ int main(int argc, char *argv[])
// sprintf(buf, pattern2, "2x2x2", std::get<1>(tmpVec[1]), std::get<1>(tmpVec[2]), std::get<1>(tmpVec[2]), std::get<1>(tmpVec[2]), std::get<1>(tmpVec[5]));
std::cout << std::endl;
+ log << std::endl;
}
}
@@ -3189,5 +3130,43 @@ int main(int argc, char *argv[])
// std::cout << buf << std::endl;
+ //////////////// OUTPUT QUANTITIES TABLE //////////////
+
+ std::cout << center("Levels ",tableWidth) << " | "
+ << center("|q1-q1_c|",tableWidth) << " | "
+ << center("|q2-q2_c|",tableWidth) << " | "
+ << center(" q3_c",tableWidth) << " | "
+ << center("|b1-b1_c|",tableWidth) << " | "
+ << center("|b2-b2_c|",tableWidth) << " | "
+ << center("|b3-b3_c|",tableWidth) << " | " << "\n";
+ std::cout << std::string(tableWidth*7 + 3*7, '-') << "\n";
+ log << center("Levels ",tableWidth) << " | "
+ << center("|q1-q1_c|",tableWidth) << " | "
+ << center("|q2-q2_c|",tableWidth) << " | "
+ << center(" q3_c",tableWidth) << " | "
+ << center("|b1-b1_c|",tableWidth) << " | "
+ << center("|b2-b2_c|",tableWidth) << " | "
+ << center("|b3-b3_c|",tableWidth) << " | " << "\n";
+ log << std::string(tableWidth*7 + 3*7, '-') << "\n";
+
+ int StorageCount2 = 0;
+ for(auto& v: Storage_Quantities) {
+
+ std::visit([tableWidth](auto&& arg){std::cout << center(prd(arg,5,1),tableWidth) << " | ";}, v);
+ std::visit([tableWidth, &log](auto&& arg){log << center(prd(arg,5,1),tableWidth) << " & ";}, v);
+
+ StorageCount2++;
+ if(StorageCount2 % 7 == 0 )
+ {
+ std::cout << std::endl;
+ log << std::endl;
+ }
+ }
+
+
+
+
+
+ log.close();
}
--
GitLab