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

Cleanup Storage

parent b6b34348
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/dune-microstructure.cc
+ 19
− 170
View file @ 2753f77a
......@@ -1154,8 +1154,6 @@ int main(int argc, char *argv[])
int levelCounter = 0;
///////////////////////////////////
// Create Data Storage
///////////////////////////////////
......@@ -1166,20 +1164,6 @@ int main(int argc, char *argv[])
// Storage:: #1 level #2 |q1_a-q1_c| #3 |q2_a-q2_c| #4 |q3_a-q3_c| #5 |b1_a-b1_c| #6 |b2_a-b2_c| #7 |b3_a-b3_c|
std::vector<std::variant<std::string, size_t , double>> Storage_Quantities;
// Storage:: #1 level #2 L2SymError #3 L2SymErrorOrder #4 L2Norm(sym) #5 L2Norm(sym-analytic) #6 L2Norm(phi_1)
// TODO Storage as Vectors? (q1,q2,q3) ,(b1,b2,b3) ....
// Storage.push_back(5.0);
// Storage.push_back("Second Entry");
// std::cout << "Storage[0] :" << std::get<double>(Storage[0]) << std::endl;
// std::cout << "Storage[0] :" << std::get<2>(Storage[0]) << std::endl;
// std::cout << "output variant :" << std::get<std::string>(Storage[1]) << std::endl;
// std::cout << "output variant :" << std::get<0>(Storage[1]) << std::endl;
// for(const size_t &level : numLevels) // explixite Angabe der levels.. {2,4}
for(size_t level = numLevels[0] ; level <= numLevels[1]; level++) // levels von bis.. [2,4]
......@@ -1200,16 +1184,12 @@ int main(int argc, char *argv[])
using GridView = CellGridType::LeafGridView;
const GridView gridView_CE = grid_CE.leafGridView();
// using ScalarRT = FieldVector< double, 1>;
// using VectorRT = FieldVector< double, dimWorld>;
using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
// using FuncScalar = std::function< ScalarRT(const Domain&) >;
using Func2Tensor = std::function< MatrixRT(const Domain&) >;
using VectorCT = BlockVector<FieldVector<double,1> >;
using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >;
///////////////////////////////////
// Create Lambda-Functions for material Parameters depending on microstructure
///////////////////////////////////
......@@ -1242,7 +1222,7 @@ int main(int argc, char *argv[])
unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 1);
// Print Options
bool print_debug = parameterSet.get<bool>("print_debug", false);
bool write_corrector_phi1 = parameterSet.get<bool>("write_corrector_phi1", false);
bool write_corrector_phi2 = parameterSet.get<bool>("write_corrector_phi2", false);
......@@ -1255,11 +1235,8 @@ int main(int argc, char *argv[])
// Choose a finite element space for Cell Problem
/////////////////////////////////////////////////////////
using namespace Functions::BasisFactory;
Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
// Don't do the following in real life: It has quadratic run-time in the number of vertices.
for (const auto& v1 : vertices(gridView_CE))
for (const auto& v2 : vertices(gridView_CE))
......@@ -1284,7 +1261,6 @@ int main(int argc, char *argv[])
VectorCT load_alpha1, load_alpha2, load_alpha3;
MatrixCT stiffnessMatrix_CE;
bool set_IntegralZero = parameterSet.get<bool>("set_IntegralZero", true);
bool set_oneBasisFunction_Zero = false;
bool substract_integralMean = false;
......@@ -1322,7 +1298,6 @@ int main(int argc, char *argv[])
// std::cout << " analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
///////////////////////////////////////////////
// Basis for R_sym^{2x2}
//////////////////////////////////////////////
......@@ -1354,28 +1329,11 @@ int main(int argc, char *argv[])
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha1 ,x3G_1neg);
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha2 ,x3G_2neg);
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3neg);
// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix", "--");
// printvector(std::cout, load_alpha1, "load_alpha1", "--");
//////////////////////////////////////////////////////////////////////
// Determine global indices of components for arbitrary (local) index
//////////////////////////////////////////////////////////////////////
// int arbitraryLeafIndex = parameterSet.get<unsigned int>("arbitraryLeafIndex", 0); // localIdx..assert Number < #ShapeFcts .. also input Element?
// int arbitraryElementNumber = parameterSet.get<unsigned int>("arbitraryElementNumber", 0);
//
// FieldVector<int,3> row;
// row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLeafIndex);
// printvector(std::cout, row, "row" , "--");
//
// set one basis-function to zero
if(set_oneBasisFunction_Zero)
{
......@@ -1385,7 +1343,6 @@ int main(int argc, char *argv[])
}
////////////////////////////////////////////////////
// Compute solution
////////////////////////////////////////////////////
......@@ -1394,7 +1351,6 @@ int main(int argc, char *argv[])
VectorCT x_3 = load_alpha3;
auto load_alpha1BS = load_alpha1;
// printvector(std::cout, load_alpha1, "load_alpha1 before SOLVER", "--" );
// printvector(std::cout, load_alpha2, "load_alpha2 before SOLVER", "--" );
......@@ -1474,11 +1430,9 @@ int main(int argc, char *argv[])
// printvector(std::cout, load_alpha1, "load_alpha1 AFTER SOLVER", "--" );
// printvector(std::cout, load_alpha2, "load_alpha2 AFTER SOLVER", "--" );
////////////////////////////////////////////////////////////////////////////////////
// Extract phi_alpha & M_alpha coefficients
////////////////////////////////////////////////////////////////////////////////////
VectorCT phi_1, phi_2, phi_3;
phi_1.resize(Basis_CE.size());
phi_1 = 0;
......@@ -1487,7 +1441,6 @@ int main(int argc, char *argv[])
phi_3.resize(Basis_CE.size());
phi_3 = 0;
for(size_t i=0; i<Basis_CE.size(); i++)
{
phi_1[i] = x_1[i];
......@@ -1503,7 +1456,6 @@ int main(int argc, char *argv[])
m_2[i] = x_2[phiOffset+i];
m_3[i] = x_3[phiOffset+i];
}
// assemble M_alpha's (actually iota(M_alpha) )
MatrixRT M_1(0), M_2(0), M_3(0);
......@@ -1587,7 +1539,6 @@ int main(int argc, char *argv[])
////////////////////////////////////////////////////////////////////////////
// Make a discrete function from the FE basis and the coefficient vector
//////////////////////////////////////////////////////////////////////////// // ERROR
using SolutionRange = FieldVector<double,dim>;
auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
Basis_CE,
phi_1);
......@@ -1601,7 +1552,6 @@ int main(int argc, char *argv[])
phi_3);
/////////////////////////////////////////////////////////
// Create Containers for Basis-Matrices, Correctors and Coefficients
/////////////////////////////////////////////////////////
......@@ -1611,7 +1561,6 @@ int main(int argc, char *argv[])
const std::array<MatrixRT*, 3 > mContainer = {&M_1 , &M_2, & M_3};
/////////////////////////////////////////////////////////
// Compute effective quantities: Elastic law & Prestrain
/////////////////////////////////////////////////////////
......@@ -1619,7 +1568,7 @@ int main(int argc, char *argv[])
VectorCT tmp1,tmp2;
FieldVector<double,3> B_hat ;
// compute Q
// Compute effective elastic law Q
for(size_t a = 0; a < 3; a++)
for(size_t b=0; b < 3; b++)
{
......@@ -1655,7 +1604,6 @@ int main(int argc, char *argv[])
{
std::cout << "check this Contribution: " << (coeffContainer[a]*tmp2) << std::endl; //see orthotropic.tex
}
}
log << "Computed prestrain B_hat: " << std::endl;
log << B_hat << std::endl;
......@@ -1671,8 +1619,6 @@ int main(int argc, char *argv[])
log << Beff << std::endl;
// printvector(std::cout, Beff, "Beff", "--");
auto q1_c = Q[0][0];
auto q2_c = Q[1][1];
auto q3_c = Q[2][2];
......@@ -1702,11 +1648,9 @@ int main(int argc, char *argv[])
log << "computed b2_hat: " << B_hat[1] << std::endl;
log << "computed b3_hat: " << B_hat[2] << std::endl;
//////////////////////////////////////////////////////////////
// Define Analytic Solutions
//////////////////////////////////////////////////////////////
// double b1 = (mu1*p1/4)*(beta/(theta+(1-theta)*beta))*(1-theta*(1+alpha));
// double b2 = (mu1*p1/8)*(1-theta*(1+beta*alpha));
double b1_hat = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
......@@ -1780,15 +1724,10 @@ int main(int argc, char *argv[])
std::cout << "L2Norm(phi_1): " << L2Norm << std::endl;
std::cout << " -----------------" << std::endl;
log << "L2-Error (symmetric Gradient phi_1):" << L2SymError << std::endl;
log << "L2-Norm(Symphi_1): " << L2Norm_Symphi<< std::endl;
log << "L2-Norm(SymAnalytic): " << L2Norm_SymAnalytic << std::endl;
double EOC = 0.0;
Storage_Error.push_back(L2SymError);
......@@ -1796,12 +1735,9 @@ int main(int argc, char *argv[])
if(levelCounter > 0)
{
// Storage_Error:: #1 level #2 L2SymError #3 L2SymErrorOrder #4 L2Norm(sym) #5 L2Norm(sym-analytic) #6 L2Norm(phi_1)
// TODO Storage_Error as Vectors? (q1,q2,q3) ,(b1,b2,b3) ....
// Storage_Error.push_back(5);
// Storage_Error.push_back("Second Entry");
std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter-1)*6)+1 << std::endl; // Besser std::map ???
std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter)*6)+1 << std::endl; // für Storage_Error[idx] muss idx zur compile time feststehen?!
// std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter-1)*6)+1 << std::endl; // Besser std::map ???
// std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter)*6)+1 << std::endl; // für Storage_Error[idx] muss idx zur compile time feststehen?!
auto ErrorOld = std::get<double>(Storage_Error[((levelCounter-1)*6)+1]);
auto ErrorNew = std::get<double>(Storage_Error[(levelCounter*6)+1]);
......@@ -1819,7 +1755,7 @@ int main(int argc, char *argv[])
}
//////////////////////////////////////////////////////////////////////////////////////////////
// Write Data to Matlab / Optimization
// Write Data to Matlab / Optimization-Code
//////////////////////////////////////////////////////////////////////////////////////////////
// writeMatrixToMatlab(Q, "matlab.txt");
......@@ -1850,59 +1786,15 @@ int main(int argc, char *argv[])
correctorFunction_3,
VTK::FieldInfo("corrector phi_3", VTK::FieldInfo::Type::vector, dim));
vtkWriter.write("CellProblem-result");
std::cout << "wrote data to file: CellProblem-result" << std::endl; // better with string for output name..
std::cout << "wrote data to file: CellProblem-result" << std::endl; // better with string for output name..
levelCounter++; //TODO
levelCounter++;
} // Level-Loop End
/////////////////////////////////////////
// Print Storage
/////////////////////////////////////////
// std::cout << "PRINT STORAGE: " << std::endl;
// int StorageCount = 0;
// for(auto& v: Storage) {
//
// if(StorageCount % 6 == 0 )
// std::cout<< "Level: ";
//
//
// std::visit([](auto&& arg){std::cout << arg;}, v);
// std::cout << "-----" << std::endl;
// StorageCount++;
// }
/////////////////////////////////////////
// Output Table
/////////////////////////////////////////
// NICE TABLE
// char buf[256];
// char pattern[] = "%10s %10s %7s %10s %10s %10s ";
// char pattern2[] = "%10s %10f %7.0f %10f %10f %10f ";
// sprintf(buf, pattern, "Level", "L2SymError", "Order", "L2SymNorm", "L2SNorm_ana", "L2Norm");
// std::cout << buf << std::endl;
// std::cout << std::setw(10) << " " << std::setw(10) << "Level " << "L2SymError " << std::setw(10) << "Order " << std::setw(10) << "L2SymNorm " << std::setw(10) << "L2SNorm_ana " << std::setw(10) << "L2Norm " << std::endl;
// std::cout << std::left << std::setw(10) << std::setfill(' ') << "Level" << "L2SymError " << std::setw(10) << "Order " << std::setw(10) << "L2SymNorm " << std::setw(10) << "L2SNorm_ana " << std::setw(10) << "L2Norm " << std::endl;
//WORKS:
// std::cout << std::left << std::setw(12) << std::setfill(' ') << "Level ";
// std::cout << std::left << std::setw(12) << std::setfill(' ') << "L2SymError ";
// std::cout << std::left << std::setw(12) << std::setfill(' ') << "Order ";
// std::cout << std::left << std::setw(12) << std::setfill(' ') << "L2SymNorm ";
// std::cout << std::left << std::setw(12) << std::setfill(' ') << "L2SNorm_ana ";
// std::cout << std::left << std::setw(12) << std::setfill(' ') << "L2Norm ";
// std::cout << std::endl;
// std::cout << "-----------------------------------------------------------------------------------------"<< std::endl;
//
// WORKS
int tableWidth = 12;
std::cout << center("Levels",tableWidth) << " | "
<< center("L2SymError",tableWidth) << " | "
......@@ -1920,62 +1812,22 @@ int main(int argc, char *argv[])
log << std::string(tableWidth*6 + 3*6, '-') << "\n";
// std::vector<double> tmpVec; // oder std::variant
std::vector<std::variant<const size_t , double>> tmpVec;
int StorageCount = 0;
for(auto& v: Storage_Error) {
// if(StorageCount % 6 == 0 )
// std::cout<< "Level: ";
// 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); // Anzahl-Nachkommastellen
std::visit([tableWidth, &log](auto&& arg){log << center(prd(arg,5,1),tableWidth) << " & ";}, v);
// std::visit([tmpVec](auto&& arg){tmpVec.push_back(arg);}, v);
StorageCount++;
if(StorageCount % 6 == 0 )
for(auto& v: Storage_Error)
{
// 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;
}
}
// for(double x=1.5; x<200; x +=x*2) {
// std::cout << prd(x,1,10) << " | "
// << prd(x*x,2,10) << " | "
// << prd(x*x/8.0,4,10) << "\n";
// }
// NICE TABLE
// char buf[256];
// char pattern[] = "%10s %10s %7s %10s %10s %10s ";
// char pattern2[] = "%10s %10f %7.0f %10f %10f %10f ";
// sprintf(buf, pattern, "Level", "L2SymError", "Order", "L2SymNorm", "L2SNorm_ana", "L2Norm");
// std::cout << buf << std::endl;
// sprintf(buf, pattern2, "2x2x2", 10.44, 7.78, 16.27, 1.99, 48.92);
// std::cout << buf << std::endl;
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);
StorageCount++;
if(StorageCount % 6 == 0 )
{
std::cout << std::endl;
log << std::endl;
}
}
//////////////// OUTPUT QUANTITIES TABLE //////////////
std::cout << center("Levels ",tableWidth) << " | "
<< center("|q1-q1_c|",tableWidth) << " | "
<< center("|q2-q2_c|",tableWidth) << " | "
......@@ -1994,8 +1846,8 @@ int main(int argc, char *argv[])
log << std::string(tableWidth*7 + 3*7, '-') << "\n";
int StorageCount2 = 0;
for(auto& v: Storage_Quantities) {
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);
......@@ -2009,8 +1861,5 @@ int main(int argc, char *argv[])
log.close();
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment