#include <config.h>
#include <array>
#include <vector>
#include <fstream>
#include <iostream>
#include <dune/common/indices.hh>
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/float_cmp.hh>
#include <dune/common/math.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/istl/matrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/multitypeblockmatrix.hh>
#include <dune/istl/multitypeblockvector.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/spqr.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/io.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/periodicbasis.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <dune/functions/functionspacebases/boundarydofs.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/microstructure/matrix_operations.hh>
#include <dune/microstructure/CorrectorComputer.hh>
#include <dune/microstructure/EffectiveQuantitiesComputer.hh>
#include <dune/microstructure/prestrainedMaterial.hh>
#include <dune/solvers/solvers/umfpacksolver.hh>
#include <dune/istl/eigenvalue/test/matrixinfo.hh>
#include <dune/fufem/dunepython.hh>
#include <any>
#include <variant>
#include <string>
#include <iomanip>
#include <filesystem>
using namespace Dune;
using namespace MatrixOperations;
/**
* @brief In the special case of a parametrized laminate (Lemma 4.5) in
* [Böhnlein,Neukamm,Padilla-Garza,Sander - A homogenized bending theory for prestrained plates]
* some of the effective quantities can be computed analytically.
* This test compares the analytical quantities with numerical results.
*/
//////////////////////////////////////////////////
// Infrastructure for handling periodicity
//////////////////////////////////////////////////
// Check whether two points are equal on R/Z x R/Z x R
auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>& y)
{
return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1]))
and (FloatCmp::eq(x[2],y[2]))
);
};
int main(int argc, char *argv[])
{
MPIHelper::instance(argc, argv);
Dune::Timer globalTimer;
//--- setup Log-File
std::fstream log;
std::cout << "Current path is " << std::filesystem::current_path() << '\n';
std::filesystem::path file_path = (__FILE__);
std::cout<< "File path: " << file_path<<std::endl;
std::cout << "dir_path: " << file_path.parent_path() << std::endl;
std::string dir_path = file_path.parent_path();
// ParameterTree parameterSet;
// if (argc < 2)
// ParameterTreeParser::readINITree(dir_path + "/parametrizedlaminate.parset", parameterSet);
// else
// {
// ParameterTreeParser::readINITree(argv[1], parameterSet);
// ParameterTreeParser::readOptions(argc, argv, parameterSet);
// }
// ParameterTree parameterSet;
// if (argc < 2)
// ParameterTreeParser::readINITree("parset/parametrizedlaminate.parset", parameterSet);
// else
// {
// ParameterTreeParser::readINITree(argv[1], parameterSet);
// ParameterTreeParser::readOptions(argc, argv, parameterSet);
// }
// ParameterTree parameterSet;
// if (argc < 2)
// ParameterTreeParser::readINITree("/parset/parametrizedlaminate.parset", parameterSet);
// else
// {
// ParameterTreeParser::readINITree(argv[1], parameterSet);
// ParameterTreeParser::readOptions(argc, argv, parameterSet);
// }
//--- Start Python interpreter
// Python::start();
// Python::Reference main = Python::import("__main__");
// Python::run("import math");
// Python::runStream()
// << std::endl << "import sys"
// << std::endl << "sys.path.append('" << dir_path << "')"
// << std::endl;
// Start Python interpreter
Python::start();
auto pyMain = Python::main();
pyMain.runStream()
<< std::endl << "import math"
<< std::endl << "import sys"
<< std::endl << "sys.path.append('" << dir_path << "')" << std::endl;
auto pyModule = pyMain.import("parametrized_laminate");
ParameterTree parameterSet;
pyModule.get("parameterSet").toC(parameterSet);
// read possible further parameters from the command line
ParameterTreeParser::readOptions(argc, argv, parameterSet);
// Print all parameters, to make them appear in the log file
std::cout << "Input parameters:" << std::endl;
parameterSet.report();
constexpr int dim = 3;
constexpr int dimWorld = 3;
///////////////////////////////////
// Generate the grid
///////////////////////////////////
// --- Corrector Problem Domain (-1/2,1/2)^3:
FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0, -1.0/2.0});
FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0, 1.0/2.0});
int gridLevel = 3;
std::array<int, dim> nElements = {(int)std::pow(2,gridLevel) ,(int)std::pow(2,gridLevel) ,(int)std::pow(2,gridLevel)};
std::cout << "Number of Grid-Elements in each direction: " << nElements << std::endl;
using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
CellGridType grid_CE(lower,upper,nElements);
using GridView = CellGridType::LeafGridView;
const GridView gridView_CE = grid_CE.leafGridView();
//--- Choose a finite element space for Cell Problem
using namespace Functions::BasisFactory;
Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
//--- get PeriodicIndices for periodicBasis (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))
if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
{
periodicIndices.unifyIndexPair({gridView_CE.indexSet().index(v1)}, {gridView_CE.indexSet().index(v2)});
}
//--- setup first order periodic Lagrange-Basis
auto Basis_CE = makeBasis(
gridView_CE,
power<dim>(
Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices),
flatLexicographic()
));
///////////////////////////////////
// Create prestrained material object
///////////////////////////////////
auto material = prestrainedMaterial(gridView_CE,parameterSet,pyModule);
// --- Get scale ratio
double gamma = parameterSet.get<double>("gamma",1.0);
//------------------------------------------------------------------------------------------------
//--- Compute Correctors
auto correctorComputer = CorrectorComputer(Basis_CE, material, gamma, log, parameterSet);
correctorComputer.assemble();
correctorComputer.solve();
//--- Compute effective quantities
auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,material);
effectiveQuantitiesComputer.computeEffectiveQuantities();
//--- get effective quantities
auto Qeff = effectiveQuantitiesComputer.getQeff();
auto Beff = effectiveQuantitiesComputer.getBeff();
printmatrix(std::cout, Qeff, "Matrix Qeff", "--");
printvector(std::cout, Beff, "Beff", "--");
/**
* @brief Compute analytical effective quantitites
*
*/
//Parameters
double theta = 0.25;
double mu_1 = 1.0;
double theta_mu = 2.0;
double mu_2 = theta_mu * mu_1;
double rho_1 = 1.0;
double theta_rho = 2.0;
double rho_2 = rho_1*theta_rho;
//Effective quantities
double q_1 = mu_1 * (theta_mu / (6*(theta + (1 - theta)*theta_mu)));
double q_2 = (mu_1/6.0) * ((1-theta) + theta*theta_mu);
double Beff_1 = (3.0*rho_1/2.0) * (1-(theta*(1+theta_rho)));
double Beff_2 = (3.0*rho_1/2.0) * ((1-theta*(1+theta_mu*theta_rho))/(1-theta+theta*theta_mu));
double Beff_3 = 0.0;
// std::cout << "q_1: " << q_1 << std::endl;
// std::cout << "q_2: " << q_2 << std::endl;
// std::cout << "Beff_1: " << Beff_1 << std::endl;
// std::cout << "Beff_2: " << Beff_2 << std::endl;
// std::cout << "Beff_3: " << Beff_3 << std::endl;
/**
* @brief Compare numerical and analytical results.
*
*/
if( std::abs(q_1 - Qeff[0][0]) > 1e-2)
{
std::cout << "std::abs(q_1 - Qeff[0][0]):" << std::abs(q_1 - Qeff[0][0]) << std::endl;
std::cerr << std::setprecision(9);
std::cerr << "Computed Qeff[0][0] is " << Qeff[0][0] << " but " << q_1 << " was expected!" << std::endl;
return 1;
}
if( std::abs(q_2 - Qeff[1][1]) > 1e-2)
{
std::cerr << std::setprecision(9);
std::cerr << "Computed Qeff[1][1] is " << Qeff[1][1] << " but " << q_2 << " was expected!" << std::endl;
return 1;
}
if( std::abs(Beff_1 - Beff[0]) > 1e-2)
{
std::cerr << std::setprecision(9);
std::cerr << "Computed Beff[0] is " << Beff[0] << " but " << Beff_1 << " was expected!" << std::endl;
return 1;
}
if( std::abs(Beff_2 - Beff[1]) > 1e-2)
{
std::cerr << std::setprecision(9);
std::cerr << "Computed Beff[1] is " << Beff[1] << " but " << Beff_2 << " was expected!" << std::endl;
return 1;
}
if( std::abs(Beff_3 - Beff[2]) > 1e-2)
{
std::cerr << std::setprecision(9);
std::cerr << "Computed Beff[2] is " << Beff[2] << " but " << Beff_3 << " was expected!" << std::endl;
return 1;
}
if( (Qeff[2][2] < Qeff[0][0]) || (Qeff[2][2] > Qeff[1][1]) )
{
std::cerr << std::setprecision(9);
std::cerr << "Computed Qeff[2][2] is " << Qeff[2][2] << " does not lie between Qeff[0][0]:"
<< Qeff[0][0] << " and Qeff[1][1]:" << Qeff[1][1] << " as was expected!" << std::endl;
return 1;
}
std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
std::cout << "Test parametrizedlaminatetest passed." << std::endl;
return 0;
}