diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc index 08911b48fe0b24019942ef0e05d8bb21d5e13587..c055cd404df3ba28f18d4fc9a1e4662990530415 100644 --- a/src/dune-microstructure.cc +++ b/src/dune-microstructure.cc @@ -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(); - }