diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 2d9dd9d71129b1801e35057c8df66aff606773fc..c6fef610267878aac698b61010eff5fc604a1be0 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -50,6 +50,10 @@
// #include <dune/fufem/discretizationerror.hh>
#include <boost/multiprecision/cpp_dec_float.hpp>
+// #include <boost/any.hpp>
+#include <any>
+#include <variant>
+#include <string>
using namespace Dune;
@@ -743,7 +747,7 @@ auto energy(const Basis& basis,
// TEST HIGHER PRECISION
// using float_50 = boost::multiprecision::cpp_dec_float_50;
// float_50 higherPrecEnergy = 0.0;
-
+
double energy = 0.0;
constexpr int dim = 3;
constexpr int nCompo = 3;
@@ -776,11 +780,11 @@ auto energy(const Basis& basis,
lambdaLocal.bind(e);
-// FieldVector<double,1> elementEnergy(0); //!!!
+// FieldVector<double,1> elementEnergy(0); //!!!
// printvector(std::cout, elementEnergy, "elementEnergy" , "--");
double elementEnergy = 0.0;
// double elementEnergy_HP = 0.0;
-
+
auto geometry = e.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
@@ -1038,8 +1042,8 @@ double computeL2SymErrorNew(const Basis& basis,
const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
const auto nSf = localFiniteElement.localBasis().size();
// std::cout << "print nSf:" << nSf << std::endl;
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); //TEST
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 ); //TEST
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); //TEST
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 ); //TEST
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
@@ -1064,11 +1068,11 @@ double computeL2SymErrorNew(const Basis& basis,
MatrixRT tmp(0);
double sum = 0.0;
-
+
for (size_t k=0; k < nCompo; k++)
- for (size_t i=0; i < nSf; i++)
+ for (size_t i=0; i < nSf; i++)
{
-
+
size_t localIdx1 = localView.tree().child(k).localIndex(i); // hier i:leafIdx
size_t globalIdx1 = localView.index(localIdx1);
@@ -1088,22 +1092,22 @@ double computeL2SymErrorNew(const Basis& basis,
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
}
-
+
// printmatrix(std::cout, strainU, "strainU", "--");
// printmatrix(std::cout, tmp, "tmp", "--");
tmp += strainU;
}
// printmatrix(std::cout, matrixField(quadPos), "matrixField(quadPos)", "--");
-// printmatrix(std::cout, tmp, "tmp", "--"); // TEST Symphi_1 hat ganz andere Struktur als analytic?
+// printmatrix(std::cout, tmp, "tmp", "--"); // TEST Symphi_1 hat ganz andere Struktur als analytic?
// tmp = tmp - matrixField(quadPos);
-// printmatrix(std::cout, tmp - matrixField(quadPos), "Difference", "--");
-
-
+// printmatrix(std::cout, tmp - matrixField(quadPos), "Difference", "--");
+
+
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
- {
+ {
sum += std::pow(tmp[ii][jj]-matrixField(quadPos)[ii][jj],2);
}
// std::cout << "sum:" << sum << std::endl;
@@ -1127,7 +1131,7 @@ double computeL2SymError(const Basis& basis,
const MatrixFunction& matrixFieldFunc
)
{
- // This Version computes the SCALAR PRODUCT of the difference ..
+ // This Version computes the SCALAR PRODUCT of the difference ..
double error = 0.0;
constexpr int dim = 3;
constexpr int nCompo = 3;
@@ -1169,7 +1173,7 @@ double computeL2SymError(const Basis& basis,
// std::vector< VectorRT> gradients(referenceGradients.size());
for (size_t i=0; i<gradients.size(); i++)
- jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]); // TRANSFORMED
+ jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]); // TRANSFORMED
// MatrixRT defGradientU(0), defGradientV(0);
@@ -1186,7 +1190,7 @@ double computeL2SymError(const Basis& basis,
size_t globalIdx1 = localView.index(localIdx1);
size_t localIdx2 = localView.tree().child(l).localIndex(j);
size_t globalIdx2 = localView.index(localIdx2);
-
+
MatrixRT defGradientU(0);
// (scaled) Deformation gradient of the ansatz basis function
defGradientU[k][0] = gradients[i][0]; // Y //hier i:leaf oder localIdx?
@@ -1218,7 +1222,7 @@ double computeL2SymError(const Basis& basis,
}
}
-
+
for (size_t k=0; k < nCompo; k++)
for (size_t i=0; i < nSf; i++)
{
@@ -1313,11 +1317,11 @@ double computeL2SymNormCoeff(const Basis& basis,
MatrixRT tmp(0);
double sum = 0.0;
-
- for (size_t i=0; i < nSf; i++)
+
+ for (size_t i=0; i < nSf; i++)
for (size_t k=0; k < nCompo; k++)
{
-
+
size_t localIdx1 = localView.tree().child(k).localIndex(i); // hier i:leafIdx
size_t globalIdx1 = localView.index(localIdx1);
@@ -1334,33 +1338,33 @@ double computeL2SymNormCoeff(const Basis& basis,
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
}
-
+
// printmatrix(std::cout, strainU, "strainU", "--");
// printmatrix(std::cout, tmp, "tmp", "--");
tmp += strainU;
// printmatrix(std::cout, tmp, "tmp", "--");
-
-
+
+
// Local energy density: stress times strain
// double tmp = 0.0;
// for (int ii=0; ii<nCompo; ii++)
// for (int jj=0; jj<nCompo; jj++)
// tmp[ii][jj] += coeffVector[globalIdx1]*coeffVector[globalIdx2]* strainU[ii][jj] * strainV[ii][jj];
-//
-//
+//
+//
// error += tmp * quadPoint.weight() * integrationElement;
// }
}
-
+
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
sum += std::pow(tmp[ii][jj],2);
-
+
error += sum * quadPoint.weight() * integrationElement;
-
+
}
}
-
+
return sqrt(error);
}
@@ -1423,9 +1427,9 @@ double computeL2SymNormCoeffV2(const Basis& basis,
for (size_t k=0; k < nCompo; k++)
-
-
-
+
+
+
for (size_t i=0; i < nSf; i++)
{
for (size_t l=0; l < nCompo; l++)
@@ -1471,7 +1475,7 @@ double computeL2SymNormCoeffV2(const Basis& basis,
}
}
-
+
return sqrt(error);
}
@@ -1574,7 +1578,7 @@ double computeL2SymNormCoeffV3(const Basis& basis,
}
}
-
+
return sqrt(error);
}
*/
@@ -1602,7 +1606,7 @@ double computeL2NormCoeff(const Basis& basis,
// std::cout << "print nSf:" << nSf << std::endl;
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-
+
for (const auto& quadPoint : quad)
{
@@ -1613,19 +1617,19 @@ double computeL2NormCoeff(const Basis& basis,
std::vector< FieldVector<double, 1>> functionValues;
localFiniteElement.localBasis().evaluateFunction(quadPoint.position(),
functionValues);
-
+
// for (auto&& value : functionValues)
// std::cout << value << std::endl;
-
- /*
+
+ /*
using LocalFEType = LagrangeSimplexLocalFiniteElement<double,double,dim,2>;
LocalFEType localFE;
// Get shape function values
using ValueType = LocalFEType::Traits::LocalBasisType::Traits::RangeType;
std::vector<ValueType> values;*/
-
+
double tmp = 0.0;
- for (size_t k=0; k < nCompo; k++) // ANALOG STOKES Beitrag nur wenn k=l!!!
+ for (size_t k=0; k < nCompo; k++) // ANALOG STOKES Beitrag nur wenn k=l!!!
{
double comp = 0.0;
for (size_t i=0; i < nSf; i++)
@@ -1639,13 +1643,13 @@ double computeL2NormCoeff(const Basis& basis,
// size_t localIdx2 = localView.tree().child(l).localIndex(j);
// size_t globalIdx2 = localView.index(localIdx2);
-
+
comp += coeffVector[globalIdx1]*functionValues[i];
-
- // tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
- // tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
+
+ // tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
+ // tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
// std::cout << "tmp:" << tmp << std::endl;
-
+
// error += tmp * quadPoint.weight() * integrationElement;
}
tmp += std::pow(comp,2);
@@ -1682,7 +1686,7 @@ double computeL2NormCoeffV2(const Basis& basis,
// std::cout << "print nSf:" << nSf << std::endl;
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-
+
for (const auto& quadPoint : quad)
{
@@ -1693,18 +1697,18 @@ double computeL2NormCoeffV2(const Basis& basis,
std::vector< FieldVector<double, 1>> functionValues;
localFiniteElement.localBasis().evaluateFunction(quadPoint.position(),
functionValues);
-
+
// for (auto&& value : functionValues)
// std::cout << value << std::endl;
-
- /*
+
+ /*
using LocalFEType = LagrangeSimplexLocalFiniteElement<double,double,dim,2>;
LocalFEType localFE;
// Get shape function values
using ValueType = LocalFEType::Traits::LocalBasisType::Traits::RangeType;
std::vector<ValueType> values;*/
-
- for (size_t k=0; k < nCompo; k++) // ANALOG STOKES Beitrag nur wenn k=l!!!
+
+ for (size_t k=0; k < nCompo; k++) // ANALOG STOKES Beitrag nur wenn k=l!!!
for (size_t i=0; i < nSf; i++)
{
// for (size_t l=0; l < nCompo; l++)
@@ -1718,11 +1722,11 @@ double computeL2NormCoeffV2(const Basis& basis,
double tmp = 0.0;
tmp += coeffVector[globalIdx1]*functionValues[i]*coeffVector[globalIdx2]*functionValues[j];
-
-// tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
-// tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
+
+// tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
+// tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
// std::cout << "tmp:" << tmp << std::endl;
-
+
error += tmp * quadPoint.weight() * integrationElement;
}
}
@@ -1757,30 +1761,30 @@ double computeL2NormCoeffV3(const Basis& basis, // Falsch: betrachtet k
// std::cout << "print nSf:" << nSf << std::endl;
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-
+
for (const auto& quadPoint : quad)
{
const auto& quadPos = quadPoint.position();
// const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
-
+
const auto integrationElement = geometry.integrationElement(quadPoint.position());
std::vector< FieldVector<double, 1>> functionValues;
localFiniteElement.localBasis().evaluateFunction(quadPoint.position(),
functionValues);
-
+
// for (auto&& value : functionValues)
// std::cout << value << std::endl;
-
- /*
+
+ /*
using LocalFEType = LagrangeSimplexLocalFiniteElement<double,double,dim,2>;
LocalFEType localFE;
// Get shape function values
using ValueType = LocalFEType::Traits::LocalBasisType::Traits::RangeType;
std::vector<ValueType> values;*/
-
+
for (size_t k=0; k < nCompo; k++)
for (size_t i=0; i < nSf; i++)
{
@@ -1788,10 +1792,10 @@ double computeL2NormCoeffV3(const Basis& basis, // Falsch: betrachtet k
size_t globalIdx1 = localView.index(localIdx1);
double tmp = 0.0;
-
- tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
+
+ tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
// std::cout << "tmp:" << tmp << std::endl;
-
+
error += tmp * quadPoint.weight() * integrationElement;
}
}
@@ -1811,7 +1815,7 @@ double computeL2NormFunc(const Basis& basis,
constexpr int dim = 3;
constexpr int nCompo = 3;
auto localView = basis.localView();
-
+
using SolutionRange = FieldVector<double,dim>;
auto GVFunc = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(basis,coeffVector);
@@ -1838,10 +1842,10 @@ double computeL2NormFunc(const Basis& basis,
const auto& quadPos = quadPoint.position();
double tmp = localfun(quadPos) * localfun(quadPos);
error += tmp * quadPoint.weight() * integrationElement;
-
+
}
}
-
+
// std::cout << "Domain-Area: " << area << std::endl;
return sqrt(error);
}
@@ -2021,7 +2025,7 @@ double computeL2SymErrorSPTest(const Basis& basis,
const MatrixFunction& matrixFieldFunc
)
{
- // This Version computes the SCALAR PRODUCT of the difference ..
+ // This Version computes the SCALAR PRODUCT of the difference ..
double error = 0.0;
constexpr int dim = 3;
constexpr int nCompo = 3;
@@ -2034,11 +2038,11 @@ double computeL2SymErrorSPTest(const Basis& basis,
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
-
+
Vector zeroVec;
zeroVec.resize(basis.size());
zeroVec = 0;
-
+
double Term1 = std::pow(computeL2SymNormCoeff(basis,coeffVector,gamma),2);
// std::cout << "Term1:" << Term1 << std::endl;
double Term2 = std::pow(computeL2SymError(basis,zeroVec ,gamma,matrixFieldFunc) ,2);
@@ -2056,9 +2060,9 @@ double computeL2SymErrorSPTest(const Basis& basis,
// std::cout << "print nSf:" << nSf << std::endl;
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-
-
-
+
+
+
for (const auto& quadPoint : quad)
{
@@ -2075,7 +2079,7 @@ double computeL2SymErrorSPTest(const Basis& basis,
// std::vector< VectorRT> gradients(referenceGradients.size());
for (size_t i=0; i<gradients.size(); i++)
- jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]); // TRANSFORMED
+ jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]); // TRANSFORMED
for (size_t k=0; k < nCompo; k++)
@@ -2109,7 +2113,7 @@ double computeL2SymErrorSPTest(const Basis& basis,
}
}
-
+
return sqrt(error+Term1+Term2);
}
@@ -2195,651 +2199,712 @@ int main(int argc, char *argv[])
FieldVector<double,dim> lower({0.0, 0.0, -1.0/2.0});
FieldVector<double,dim> upper({1.0, 1.0, 1.0/2.0});
}
-
-
-
+
+
+
std::array<int, dim> nElements = parameterSet.get<std::array<int,dim>>("nElements_Cell", {10, 10, 10});
-
-
- ///// LEVEL - APPROACH // TODO
-/* int levels = parameterSet.get<int>("levels", 1); //besser : numLevels
- std::array<int, dim> nElements = { std::pow(2,levels) , std::pow(2,levels) , std::pow(2,levels) }; */
-
- std::cout << "Number of Elements in each direction: " << nElements << std::endl;
- log << "Number of 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();
-
-// using ScalarRT = FieldVector< double, 1>;
-// using VectorRT = FieldVector< double, nCompo>;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
- 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> >;
+ ///// LEVEL - APPROACH // TODO
+ int numLevels = parameterSet.get<int>("numLevels", 1); //besser : numLevels
- ///////////////////////////////////
- // Get Material Parameters
- ///////////////////////////////////
- double beta = parameterSet.get<double>("beta",2.0);
- double mu1 = parameterSet.get<double>("mu1",10.0);;
- double mu2 = beta*mu1;
- double lambda1 = parameterSet.get<double>("lambda1",0.0);;
- double lambda2 = beta*lambda1;
+// // any type
+// std::any a = 1;
+// std::cout << a.type().name() << ": " << std::any_cast<int>(a) << '\n';
+// std::vector<std::any> Storage;
+// std::any Test=1;
+// Storage[0]= 1;
+// std::cout << std::any_cast<int>(Test) << '\n';
+// std::cout << std::any_cast<int>(Storage[0]) << '\n';
+//
+// // Test = "Test";
+//
+// Test = std::make_any<std::string>("Hello World");
+//
+// std::cout << "Test: " << std::any_cast<std::string>(Test) << std::endl;
+//
+// // Test = 5;
+// // std::cout << "Test: " << std::any_cast<int>(Test) << std::endl;
+//
+// Test = 1e-03;
+//
+//
+// std::cout << "Test: " << std::any_cast<double>(Test) << std::endl;
- ///////////////////////////////////
- // Create Lambda-Functions for material Parameters depending on microstructure
- ///////////////////////////////////
- auto muTerm = [mu1, mu2, theta] (const Domain& z) {
- if (abs(z[0]) >= (theta/2.0)) //TODO check ..nullset/boundary
- return mu1;
-// if (abs(z[0]) > (theta/2.0)) //TEST include boundary (indicatorFunction)
-// return mu1;
- // if (abs(z[0]) < (theta/2) && z[2] < 0) // oder das?
- // return mu2;
- else
- return mu2;
- // return mu1;
- };
-
- auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) {
- if (abs(z[0]) >= (theta/2.0))
- return lambda1;
+
+ // std::variant
+
+ std::variant<int,double, std::string> a;
+ std::variant<std::string> x("abc");
+
+
+ std::variant<int, std::string> variant = "Hello";
+
+ std::string string_1 = std::get<std::string>( variant ); // get value by type
+ std::string string_2 = std::get<1>( variant ); // get value by index
+ std::cout << string_1 << std::endl;
+ std::cout << string_2 << std::endl;
+
+ std::cout << "output variant :" << std::get<std::string>(x) << std::endl;
+
+ std::vector<std::variant<std::string, int , double >> Storage;
+ Storage.push_back(5);
+ Storage.push_back("Second Entry");
+// Storage[1] = "Second Entry";
+
+ std::cout << "Storage[0] :" << std::get<int>(Storage[0]) << std::endl;
+ std::cout << "Storage[0] :" << std::get<1>(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(size_t level = 0; level < numLevels; level++)
+ {
+
+
+ std::array<int, dim> nElements = { (int)std::pow(2,level) , (int)std::pow(2,level) , (int)std::pow(2,level) };
+
+ std::cout << "Number of Elements in each direction: " << nElements << std::endl;
+ log << "Number of 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();
+
+ // using ScalarRT = FieldVector< double, 1>;
+ // using VectorRT = FieldVector< double, nCompo>;
+ using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+ 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> >;
+
+
+ ///////////////////////////////////
+ // Get Material Parameters
+ ///////////////////////////////////
+ double beta = parameterSet.get<double>("beta",2.0);
+ double mu1 = parameterSet.get<double>("mu1",10.0);;
+ double mu2 = beta*mu1;
+ double lambda1 = parameterSet.get<double>("lambda1",0.0);;
+ double lambda2 = beta*lambda1;
+
+ ///////////////////////////////////
+ // Create Lambda-Functions for material Parameters depending on microstructure
+ ///////////////////////////////////
+ auto muTerm = [mu1, mu2, theta] (const Domain& z) {
+ if (abs(z[0]) >= (theta/2.0)) //TODO check ..nullset/boundary
+ return mu1;
+ // if (abs(z[0]) > (theta/2.0)) //TEST include boundary (indicatorFunction)
+ // return mu1;
+ // if (abs(z[0]) < (theta/2) && z[2] < 0) // oder das?
+ // return mu2;
else
- return lambda2;
+ return mu2;
+ // return mu1;
};
- auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
- auto muLocal = localFunction(muGridF);
- auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
- auto lambdaLocal = localFunction(lambdaGridF);
+ auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) {
+ if (abs(z[0]) >= (theta/2.0))
+ return lambda1;
+ else
+ return lambda2;
+ };
+
+ auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
+ auto muLocal = localFunction(muGridF);
+ auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
+ auto lambdaLocal = localFunction(lambdaGridF);
+
+ log << "beta: " << beta << std::endl;
+ log << "material parameters: " << std::endl;
+ log << "mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
+ log << "lambda1: " << lambda1 <<"\nlambda2:" << lambda2 << std::endl;
+
+
+ ///////////////////////////////////
+ // --- Choose Solver ---
+ // 1 : CG-Solver
+ // 2 : GMRES
+ // 3 : QR
+ ///////////////////////////////////
+ unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 1);
+ bool write_corrector_phi1 = parameterSet.get<bool>("write_corrector_phi1", false);
+ bool write_corrector_phi2 = parameterSet.get<bool>("write_corrector_phi2", false);
+ bool write_corrector_phi3 = parameterSet.get<bool>("write_corrector_phi3", false);
+ bool write_L2Error = parameterSet.get<bool>("write_L2Error", false);
+ bool write_IntegralMean = parameterSet.get<bool>("write_IntegralMean", false);
+
+
+ /////////////////////////////////////////////////////////
+ // Choose a finite element space for Cell Problem
+ /////////////////////////////////////////////////////////
+ using namespace Functions::BasisFactory;
+
+ Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
+
+ int equivCounter = 0; // TODO remove
+
+ // 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))
+ {
+ // std::cout << " ---------------------------------------" << std::endl;
+ 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)});
+ // equivCounter++;
+ }
+ }
+ // std::cout << "equivCounter:" << equivCounter << std::endl;
+ // std::cout << "Number of Elements in each direction: " << nE << std::endl;
+ // std::cout << "Number ofNodes: " << gridView_CE.size(dim) << std::endl;
- log << "beta: " << beta << std::endl;
- log << "material parameters: " << std::endl;
- log << "mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
- log << "lambda1: " << lambda1 <<"\nlambda2:" << lambda2 << std::endl;
+ auto perTest = periodicIndices.indexPairSet(); // TODO remove
- ///////////////////////////////////
- // --- Choose Solver ---
- // 1 : CG-Solver
- // 2 : GMRES
- // 3 : QR
- ///////////////////////////////////
- unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 1);
- bool write_corrector_phi1 = parameterSet.get<bool>("write_corrector_phi1", false);
- bool write_corrector_phi2 = parameterSet.get<bool>("write_corrector_phi2", false);
- bool write_corrector_phi3 = parameterSet.get<bool>("write_corrector_phi3", false);
- bool write_L2Error = parameterSet.get<bool>("write_L2Error", false);
- bool write_IntegralMean = parameterSet.get<bool>("write_IntegralMean", false);
+ auto Basis_CE = makeBasis(
+ gridView_CE,
+ power<dim>( //lagrange<1>(),
+ Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices),
+ flatLexicographic())); //
+ // blockedInterleaved())); // siehe Periodicbasistest.. funktioniert auch?
+ // flatInterleaved()));
-
- /////////////////////////////////////////////////////////
- // Choose a finite element space for Cell Problem
- /////////////////////////////////////////////////////////
- using namespace Functions::BasisFactory;
+ std::cout << "size feBasis: " << Basis_CE.size() << std::endl;
+ std::cout << "basis.dimension()" << Basis_CE.dimension() <<std::endl;
+ log << "size of FiniteElementBasis: " << Basis_CE.size() << std::endl;
- Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
+ const int phiOffset = Basis_CE.size();
+ debugLog << "phiOffset: "<< phiOffset << std::endl;
- int equivCounter = 0; // TODO remove
+ /////////////////////////////////////////////////////////
+ // Data structures: Stiffness matrix and right hand side vector
+ /////////////////////////////////////////////////////////
+ VectorCT load_alpha1, load_alpha2, load_alpha3;
+ MatrixCT stiffnessMatrix_CE;
- // 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))
- {
-// std::cout << " ---------------------------------------" << std::endl;
- 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)});
-// equivCounter++;
- }
- }
-// std::cout << "equivCounter:" << equivCounter << std::endl;
-// std::cout << "Number of Elements in each direction: " << nE << std::endl;
-// std::cout << "Number ofNodes: " << gridView_CE.size(dim) << std::endl;
+ bool set_IntegralZero = parameterSet.get<bool>("set_IntegralZero", true);
+ bool set_oneBasisFunction_Zero = false;
+ bool substract_integralMean = false;
+ if(set_IntegralZero)
+ {
+ set_oneBasisFunction_Zero = true;
+ substract_integralMean = true;
+ }
- auto perTest = periodicIndices.indexPairSet(); // TODO remove
- auto Basis_CE = makeBasis(
- gridView_CE,
- power<dim>( //lagrange<1>(),
- Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices),
- flatLexicographic())); //
-// blockedInterleaved())); // siehe Periodicbasistest.. funktioniert auch?
-// flatInterleaved()));
+ /////////////////////////////////////////////////////////
+ // Define "rhs"-functions
+ /////////////////////////////////////////////////////////
+ Func2Tensor x3G_1 = [] (const Domain& x) {
+ return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ };
- std::cout << "size feBasis: " << Basis_CE.size() << std::endl;
- std::cout << "basis.dimension()" << Basis_CE.dimension() <<std::endl;
- log << "size of FiniteElementBasis: " << Basis_CE.size() << std::endl;
+ Func2Tensor x3G_2 = [] (const Domain& x) {
+ return MatrixRT{{0.0, 0.0, 0.0}, {0.0, 1.0*x[2], 0.0}, {0.0, 0.0, 0.0}};
+ };
- const int phiOffset = Basis_CE.size();
- debugLog << "phiOffset: "<< phiOffset << std::endl;
+ Func2Tensor x3G_3 = [] (const Domain& x) {
+ return MatrixRT{{0.0, 1.0/sqrt(2)*x[2], 0.0}, {1.0/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ };
- /////////////////////////////////////////////////////////
- // Data structures: Stiffness matrix and right hand side vector
- /////////////////////////////////////////////////////////
- VectorCT load_alpha1, load_alpha2, load_alpha3;
- MatrixCT stiffnessMatrix_CE;
+ Func2Tensor x3G_1neg = [x3G_1] (const Domain& x) {
+ return -1.0*x3G_1(x);
+ };
+ Func2Tensor x3G_2neg = [x3G_2] (const Domain& x) {
+ return -1.0*x3G_2(x);;
+ };
- bool set_IntegralZero = parameterSet.get<bool>("set_IntegralZero", true);
- bool set_oneBasisFunction_Zero = false;
- bool substract_integralMean = false;
- if(set_IntegralZero)
- {
- set_oneBasisFunction_Zero = true;
- substract_integralMean = true;
- }
+ Func2Tensor x3G_3neg = [x3G_3] (const Domain& x) {
+ return -1.0*x3G_3(x);
+ };
- /////////////////////////////////////////////////////////
- // Define "rhs"-functions
- /////////////////////////////////////////////////////////
- Func2Tensor x3G_1 = [] (const Domain& x) {
- return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
- };
-
- Func2Tensor x3G_2 = [] (const Domain& x) {
- return MatrixRT{{0.0, 0.0, 0.0}, {0.0, 1.0*x[2], 0.0}, {0.0, 0.0, 0.0}};
- };
-
- Func2Tensor x3G_3 = [] (const Domain& x) {
- return MatrixRT{{0.0, 1.0/sqrt(2)*x[2], 0.0}, {1.0/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
- };
-
- Func2Tensor x3G_1neg = [x3G_1] (const Domain& x) {
- return -1.0*x3G_1(x);
- };
-
- Func2Tensor x3G_2neg = [x3G_2] (const Domain& x) {
- return -1.0*x3G_2(x);;
- };
-
- Func2Tensor x3G_3neg = [x3G_3] (const Domain& x) {
- return -1.0*x3G_3(x);
- };
-
-
-
- // TEST - energy method ///
- // different indicatorFunction in muTerm has impact here !!
-
-// double GGterm = 0.0;
-// GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3G_1, x3G_1 ); // <L i(x3G_alpha) , i(x3G_beta) >
-// std::cout << "GGTerm:" << GGterm << std::endl;
-// std::cout << " analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
-//
-
-
- //////////////////////////////////////////////////
-
-
-
-
+ // TEST - energy method ///
+ // different indicatorFunction in muTerm has impact here !!
- ///////////////////////////////////////////////
- // Basis for R_sym^{2x2}
- //////////////////////////////////////////////
- MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
- MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
- MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
- std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
-// printmatrix(std::cout, basisContainer[0] , "G_1", "--");
-// printmatrix(std::cout, basisContainer[1] , "G_2", "--");
-// printmatrix(std::cout, basisContainer[2] , "´G_3", "--");
-// log << basisContainer[0] << std::endl;
+ // double GGterm = 0.0;
+ // GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3G_1, x3G_1 ); // <L i(x3G_alpha) , i(x3G_beta) >
+ // std::cout << "GGTerm:" << GGterm << std::endl;
+ // std::cout << " analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
+ //
- //////////////////////////////////////////////////////////////////////
- // 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" , "--");
- /////////////////////////////////////////////////////////
- // Assemble the system
- /////////////////////////////////////////////////////////
- assembleCellStiffness(Basis_CE, muLocal, lambdaLocal, gamma, stiffnessMatrix_CE, parameterSet);
- 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" , "--");
-//
+ ///////////////////////////////////////////////
+ // Basis for R_sym^{2x2}
+ //////////////////////////////////////////////
+ MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+ MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+ MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
+ // printmatrix(std::cout, basisContainer[0] , "G_1", "--");
+ // printmatrix(std::cout, basisContainer[1] , "G_2", "--");
+ // printmatrix(std::cout, basisContainer[2] , "´G_3", "--");
+ // log << basisContainer[0] << std::endl;
-
+ //////////////////////////////////////////////////////////////////////
+ // 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)
- {
- setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3, parameterSet);
-// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
-// printvector(std::cout, load_alpha1, "load_alpha1 after setOneBaseFunctionToZero", "--");
- }
+ /////////////////////////////////////////////////////////
+ // Assemble the system
+ /////////////////////////////////////////////////////////
+ assembleCellStiffness(Basis_CE, muLocal, lambdaLocal, gamma, stiffnessMatrix_CE, parameterSet);
+ 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", "--");
- ////////////////////////////////////////////////////
- // Compute solution
- ////////////////////////////////////////////////////
- VectorCT x_1 = load_alpha1;
- VectorCT x_2 = load_alpha2;
- VectorCT x_3 = load_alpha3;
- auto load_alpha1BS = 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" , "--");
+ //
-// printvector(std::cout, load_alpha1, "load_alpha1 before SOLVER", "--" );
-// printvector(std::cout, load_alpha2, "load_alpha2 before SOLVER", "--" );
- if (Solvertype == 1) // CG - SOLVER
- {
- std::cout << "------------ CG - Solver ------------" << std::endl;
- MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_CE);
-
-
-
- // Sequential incomplete LU decomposition as the preconditioner
- SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_CE,1.0);
- int iter = parameterSet.get<double>("iterations_CG", 999);
- // Preconditioned conjugate-gradient solver
- CGSolver<VectorCT> solver(op,
- ilu0, //NULL,
- 1e-8, // desired residual reduction factorlack
- iter, // maximum number of iterations
- 2); // verbosity of the solver
- InverseOperatorResult statistics;
- std::cout << "solve linear system for x_1.\n";
- solver.apply(x_1, load_alpha1, statistics);
- std::cout << "solve linear system for x_2.\n";
- solver.apply(x_2, load_alpha2, statistics);
- std::cout << "solve linear system for x_3.\n";
- solver.apply(x_3, load_alpha3, statistics);
- log << "Solver-type used: " <<"\n CG-Solver" << std::endl;
- }
- ////////////////////////////////////////////////////////////////////////////////////
- else if (Solvertype ==2) // GMRES - SOLVER
- {
- std::cout << "------------ GMRES - Solver ------------" << std::endl;
- // Turn the matrix into a linear operator
- MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_CE);
-
- // Fancy (but only) way to not have a preconditioner at all
- Richardson<VectorCT,VectorCT> preconditioner(1.0);
-
- // Construct the iterative solver
- RestartedGMResSolver<VectorCT> solver(
- stiffnessOperator, // Operator to invert
- preconditioner, // Preconditioner
- 1e-10, // Desired residual reduction factor
- 500, // Number of iterations between restarts,
- // here: no restarting
- 500, // Maximum number of iterations
- 2); // Verbosity of the solver
-
- // Object storing some statistics about the solving process
- InverseOperatorResult statistics;
-
- // solve for different Correctors (alpha = 1,2,3)
- solver.apply(x_1, load_alpha1, statistics);
- solver.apply(x_2, load_alpha2, statistics);
- solver.apply(x_3, load_alpha3, statistics);
- log << "Solver-type used: " <<"\n GMRES-Solver" << std::endl;
- }
- ////////////////////////////////////////////////////////////////////////////////////
- else if ( Solvertype == 3)// QR - SOLVER
- {
- std::cout << "------------ QR - Solver ------------" << std::endl;
- log << "solveLinearSystems: We use QR solver.\n";
- //TODO install suitesparse
- SPQR<MatrixCT> sPQR(stiffnessMatrix_CE);
- sPQR.setVerbosity(1);
- InverseOperatorResult statisticsQR;
-
- sPQR.apply(x_1, load_alpha1, statisticsQR);
- sPQR.apply(x_2, load_alpha2, statisticsQR);
- sPQR.apply(x_3, load_alpha3, statisticsQR);
- log << "Solver-type used: " <<"\n QR-Solver" << std::endl;
- }
-// printvector(std::cout, load_alpha1BS, "load_alpha1 before SOLVER", "--" );
-// printvector(std::cout, load_alpha1, "load_alpha1 AFTER SOLVER", "--" );
-// printvector(std::cout, load_alpha2, "load_alpha2 AFTER SOLVER", "--" );
+ // set one basis-function to zero
+ if(set_oneBasisFunction_Zero)
+ {
+ setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3, parameterSet);
+ // printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
+ // printvector(std::cout, load_alpha1, "load_alpha1 after setOneBaseFunctionToZero", "--");
+ }
- ////////////////////////////////////////////////////////////////////////////////////
- // Extract phi_alpha & M_alpha coefficients
- ////////////////////////////////////////////////////////////////////////////////////
- VectorCT phi_1, phi_2, phi_3;
- phi_1.resize(Basis_CE.size());
- phi_1 = 0;
- phi_2.resize(Basis_CE.size());
- phi_2 = 0;
- phi_3.resize(Basis_CE.size());
- phi_3 = 0;
+ ////////////////////////////////////////////////////
+ // Compute solution
+ ////////////////////////////////////////////////////
+ VectorCT x_1 = load_alpha1;
+ VectorCT x_2 = load_alpha2;
+ VectorCT x_3 = load_alpha3;
+ auto load_alpha1BS = load_alpha1;
- for(size_t i=0; i<Basis_CE.size(); i++)
- {
- phi_1[i] = x_1[i];
- phi_2[i] = x_2[i];
- phi_3[i] = x_3[i];
- }
+ // printvector(std::cout, load_alpha1, "load_alpha1 before SOLVER", "--" );
+ // printvector(std::cout, load_alpha2, "load_alpha2 before SOLVER", "--" );
- FieldVector<double,3> m_1, m_2, m_3;
+ if (Solvertype == 1) // CG - SOLVER
+ {
+ std::cout << "------------ CG - Solver ------------" << std::endl;
+ MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_CE);
+
+
+
+ // Sequential incomplete LU decomposition as the preconditioner
+ SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_CE,1.0);
+ int iter = parameterSet.get<double>("iterations_CG", 999);
+ // Preconditioned conjugate-gradient solver
+ CGSolver<VectorCT> solver(op,
+ ilu0, //NULL,
+ 1e-8, // desired residual reduction factorlack
+ iter, // maximum number of iterations
+ 2); // verbosity of the solver
+ InverseOperatorResult statistics;
+ std::cout << "solve linear system for x_1.\n";
+ solver.apply(x_1, load_alpha1, statistics);
+ std::cout << "solve linear system for x_2.\n";
+ solver.apply(x_2, load_alpha2, statistics);
+ std::cout << "solve linear system for x_3.\n";
+ solver.apply(x_3, load_alpha3, statistics);
+ log << "Solver-type used: " <<"\n CG-Solver" << std::endl;
+ }
+ ////////////////////////////////////////////////////////////////////////////////////
- for(size_t i=0; i<3; i++)
- {
- m_1[i] = x_1[phiOffset+i];
- m_2[i] = x_2[phiOffset+i];
- m_3[i] = x_3[phiOffset+i];
- }
+ else if (Solvertype ==2) // GMRES - SOLVER
+ {
- // assemble M_alpha's (actually iota(M_alpha) )
+ std::cout << "------------ GMRES - Solver ------------" << std::endl;
+ // Turn the matrix into a linear operator
+ MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_CE);
+
+ // Fancy (but only) way to not have a preconditioner at all
+ Richardson<VectorCT,VectorCT> preconditioner(1.0);
+
+ // Construct the iterative solver
+ RestartedGMResSolver<VectorCT> solver(
+ stiffnessOperator, // Operator to invert
+ preconditioner, // Preconditioner
+ 1e-10, // Desired residual reduction factor
+ 500, // Number of iterations between restarts,
+ // here: no restarting
+ 500, // Maximum number of iterations
+ 2); // Verbosity of the solver
+
+ // Object storing some statistics about the solving process
+ InverseOperatorResult statistics;
+
+ // solve for different Correctors (alpha = 1,2,3)
+ solver.apply(x_1, load_alpha1, statistics);
+ solver.apply(x_2, load_alpha2, statistics);
+ solver.apply(x_3, load_alpha3, statistics);
+ log << "Solver-type used: " <<"\n GMRES-Solver" << std::endl;
+ }
+ ////////////////////////////////////////////////////////////////////////////////////
+ else if ( Solvertype == 3)// QR - SOLVER
+ {
- MatrixRT M_1(0), M_2(0), M_3(0);
+ std::cout << "------------ QR - Solver ------------" << std::endl;
+ log << "solveLinearSystems: We use QR solver.\n";
+ //TODO install suitesparse
+ SPQR<MatrixCT> sPQR(stiffnessMatrix_CE);
+ sPQR.setVerbosity(1);
+ InverseOperatorResult statisticsQR;
+
+ sPQR.apply(x_1, load_alpha1, statisticsQR);
+ sPQR.apply(x_2, load_alpha2, statisticsQR);
+ sPQR.apply(x_3, load_alpha3, statisticsQR);
+ log << "Solver-type used: " <<"\n QR-Solver" << std::endl;
+ }
+ // printvector(std::cout, load_alpha1BS, "load_alpha1 before SOLVER", "--" );
+ // printvector(std::cout, load_alpha1, "load_alpha1 AFTER SOLVER", "--" );
+ // printvector(std::cout, load_alpha2, "load_alpha2 AFTER SOLVER", "--" );
- for(size_t i=0; i<3; i++)
- {
- M_1 += m_1[i]*basisContainer[i];
- M_2 += m_2[i]*basisContainer[i];
- M_3 += m_3[i]*basisContainer[i];
- }
- std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
- printmatrix(std::cout, M_1, "Corrector-Matrix M_1", "--");
- printmatrix(std::cout, M_2, "Corrector-Matrix M_2", "--");
- printmatrix(std::cout, M_3, "Corrector-Matrix M_3", "--");
- log << "Corrector-Matrix M_1: \n" << M_1 << std::endl;
- log << " --------------------" << std::endl;
- log << "Corrector-Matrix M_2: \n" << M_2 << std::endl;
- log << " --------------------" << std::endl;
- log << "Corrector-Matrix M_3: \n" << M_3 << std::endl;
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Extract phi_alpha & M_alpha coefficients
+ ////////////////////////////////////////////////////////////////////////////////////
+ VectorCT phi_1, phi_2, phi_3;
+ phi_1.resize(Basis_CE.size());
+ phi_1 = 0;
+ phi_2.resize(Basis_CE.size());
+ phi_2 = 0;
+ phi_3.resize(Basis_CE.size());
+ phi_3 = 0;
- ////////////////////////////////////////////////////////////////////////////
- // Substract Integral-mean
- //////////////////////////////////////////////////////////////////////////// // ERROR
- using SolutionRange = FieldVector<double,dim>;
- if(write_IntegralMean)
- {
- std::cout << "check integralMean phi_1: " << std::endl;
- auto A = integralMean(Basis_CE,phi_1);
+ for(size_t i=0; i<Basis_CE.size(); i++)
+ {
+ phi_1[i] = x_1[i];
+ phi_2[i] = x_2[i];
+ phi_3[i] = x_3[i];
+ }
+
+ FieldVector<double,3> m_1, m_2, m_3;
+
for(size_t i=0; i<3; i++)
{
- std::cout << "Integral-Mean : " << A[i] << std::endl;
+ m_1[i] = x_1[phiOffset+i];
+ m_2[i] = x_2[phiOffset+i];
+ m_3[i] = x_3[phiOffset+i];
}
- }
- if(substract_integralMean)
- {
- std::cout << " --- substracting integralMean --- " << std::endl;
- subtractIntegralMean(Basis_CE, phi_1);
- subtractIntegralMean(Basis_CE, phi_2);
- subtractIntegralMean(Basis_CE, phi_3);
- subtractIntegralMean(Basis_CE, x_1);
- subtractIntegralMean(Basis_CE, x_2);
- subtractIntegralMean(Basis_CE, x_3);
- //////////////////////////////////////////
- // CHECK INTEGRAL-MEAN:
- //////////////////////////////////////////
+
+ // assemble M_alpha's (actually iota(M_alpha) )
+
+ MatrixRT M_1(0), M_2(0), M_3(0);
+
+ for(size_t i=0; i<3; i++)
+ {
+ M_1 += m_1[i]*basisContainer[i];
+ M_2 += m_2[i]*basisContainer[i];
+ M_3 += m_3[i]*basisContainer[i];
+ }
+
+ std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
+ printmatrix(std::cout, M_1, "Corrector-Matrix M_1", "--");
+ printmatrix(std::cout, M_2, "Corrector-Matrix M_2", "--");
+ printmatrix(std::cout, M_3, "Corrector-Matrix M_3", "--");
+ log << "Corrector-Matrix M_1: \n" << M_1 << std::endl;
+ log << " --------------------" << std::endl;
+ log << "Corrector-Matrix M_2: \n" << M_2 << std::endl;
+ log << " --------------------" << std::endl;
+ log << "Corrector-Matrix M_3: \n" << M_3 << std::endl;
+
+
+ ////////////////////////////////////////////////////////////////////////////
+ // Substract Integral-mean
+ //////////////////////////////////////////////////////////////////////////// // ERROR
+ using SolutionRange = FieldVector<double,dim>;
+
if(write_IntegralMean)
{
- auto A = integralMean(Basis_CE, phi_1);
+ std::cout << "check integralMean phi_1: " << std::endl;
+ auto A = integralMean(Basis_CE,phi_1);
for(size_t i=0; i<3; i++)
{
- std::cout << "Integral-Mean (CHECK) : " << A[i] << std::endl;
+ std::cout << "Integral-Mean : " << A[i] << std::endl;
+ }
+ }
+ if(substract_integralMean)
+ {
+ std::cout << " --- substracting integralMean --- " << std::endl;
+ subtractIntegralMean(Basis_CE, phi_1);
+ subtractIntegralMean(Basis_CE, phi_2);
+ subtractIntegralMean(Basis_CE, phi_3);
+ subtractIntegralMean(Basis_CE, x_1);
+ subtractIntegralMean(Basis_CE, x_2);
+ subtractIntegralMean(Basis_CE, x_3);
+ //////////////////////////////////////////
+ // CHECK INTEGRAL-MEAN:
+ //////////////////////////////////////////
+ if(write_IntegralMean)
+ {
+ auto A = integralMean(Basis_CE, phi_1);
+ for(size_t i=0; i<3; i++)
+ {
+ std::cout << "Integral-Mean (CHECK) : " << A[i] << std::endl;
+ }
}
}
- }
- /////////////////////////////////////////////////////////
- // Write Solution (Corrector Coefficients) in Logs
- /////////////////////////////////////////////////////////
- log << "\nSolution of Corrector problems:\n";
- if(write_corrector_phi1)
- {
- log << "\n Corrector_phi1:\n";
- log << x_1 << std::endl;
- }
- if(write_corrector_phi2)
- {
- log << "-----------------------------------------------------";
- log << "\n Corrector_phi2:\n";
- log << x_2 << std::endl;
- }
- if(write_corrector_phi3)
- {
- log << "-----------------------------------------------------";
- log << "\n Corrector_phi3:\n";
- log << x_3 << std::endl;
- }
+ /////////////////////////////////////////////////////////
+ // Write Solution (Corrector Coefficients) in Logs
+ /////////////////////////////////////////////////////////
+ log << "\nSolution of Corrector problems:\n";
+ if(write_corrector_phi1)
+ {
+ log << "\n Corrector_phi1:\n";
+ log << x_1 << std::endl;
+ }
+ if(write_corrector_phi2)
+ {
+ log << "-----------------------------------------------------";
+ log << "\n Corrector_phi2:\n";
+ log << x_2 << std::endl;
+ }
+ if(write_corrector_phi3)
+ {
+ log << "-----------------------------------------------------";
+ log << "\n Corrector_phi3:\n";
+ log << x_3 << std::endl;
+ }
- ////////////////////////////////////////////////////////////////////////////
- // 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);
+ ////////////////////////////////////////////////////////////////////////////
+ // 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);
- auto correctorFunction_2 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
- Basis_CE,
- phi_2);
+ auto correctorFunction_2 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ Basis_CE,
+ phi_2);
- auto correctorFunction_3 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
- Basis_CE,
- phi_3);
+ auto correctorFunction_3 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ Basis_CE,
+ phi_3);
- /////////////////////////////////////////////////////////
- // Create Containers for Basis-Matrices, Correctors and Coefficients
- /////////////////////////////////////////////////////////
- const std::array<Func2Tensor, 3> x3MatrixBasis = {x3G_1, x3G_2, x3G_3};
- const std::array<VectorCT, 3> coeffContainer = {x_1, x_2, x_3};
- const std::array<VectorCT, 3> loadContainer = {load_alpha1, load_alpha2, load_alpha3};
- const std::array<MatrixRT*, 3 > mContainer = {&M_1 , &M_2, & M_3};
+ /////////////////////////////////////////////////////////
+ // Create Containers for Basis-Matrices, Correctors and Coefficients
+ /////////////////////////////////////////////////////////
+ const std::array<Func2Tensor, 3> x3MatrixBasis = {x3G_1, x3G_2, x3G_3};
+ const std::array<VectorCT, 3> coeffContainer = {x_1, x_2, x_3};
+ const std::array<VectorCT, 3> loadContainer = {load_alpha1, load_alpha2, load_alpha3};
+ const std::array<MatrixRT*, 3 > mContainer = {&M_1 , &M_2, & M_3};
- /////////////////////////////////////////////////////////
- // Compute effective quantities: Elastic law & Prestrain
- /////////////////////////////////////////////////////////
- MatrixRT Q(0);
- VectorCT tmp1,tmp2;
- FieldVector<double,3> B_hat ;
+ /////////////////////////////////////////////////////////
+ // Compute effective quantities: Elastic law & Prestrain
+ /////////////////////////////////////////////////////////
+ MatrixRT Q(0);
+ VectorCT tmp1,tmp2;
+ FieldVector<double,3> B_hat ;
- // compute Q
- for(size_t a = 0; a < 3; a++)
- for(size_t b=0; b < 3; b++)
+ // compute Q
+ for(size_t a = 0; a < 3; a++)
+ for(size_t b=0; b < 3; b++)
+ {
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp1 ,x3MatrixBasis[b]); // <L i(M_alpha) + sym(grad phi_alpha), i(x3G_beta) >
+
+ double GGterm = 0.0;
+ GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
+ // TEST
+ //TEST
+ std::cout << " analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
+
+ std::setprecision(std::numeric_limits<float>::digits10);
+
+ Q[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric... check positiv definitness?
+ // Q[a][b] = (coeffContainer[a]*loadContainer[b]) + GGterm; // TEST
+ std::cout << "GGTerm:" << GGterm << std::endl;
+ std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
+ // std::cout << "Test : " << coeffContainer[a]*loadContainer[b] << std::endl; //same...
+ }
+ printmatrix(std::cout, Q, "Matrix Q", "--");
+ log << "Computed Matrix Q: " << std::endl;
+ log << Q << std::endl;
+
+ // compute B_hat
+ for(size_t a = 0; a < 3; a++)
{
- assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp1 ,x3MatrixBasis[b]); // <L i(M_alpha) + sym(grad phi_alpha), i(x3G_beta) >
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp2 ,B_Term); // <L i(M_alpha) + sym(grad phi_alpha), B >
- double GGterm = 0.0;
- GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
- // TEST
- //TEST
- std::cout << " analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
-
- std::setprecision(std::numeric_limits<float>::digits10);
-
- Q[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric... check positiv definitness?
-// Q[a][b] = (coeffContainer[a]*loadContainer[b]) + GGterm; // TEST
- std::cout << "GGTerm:" << GGterm << std::endl;
- std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
-// std::cout << "Test : " << coeffContainer[a]*loadContainer[b] << std::endl; //same...
+ auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , B_Term); // <L i(x3G_alpha) , B>
+
+ B_hat[a] = (coeffContainer[a]*tmp2) + GGterm;
+ std::cout << "check this Contribution: " << (coeffContainer[a]*tmp2) << std::endl; //see orthotropic.tex
+ std::cout <<"B_hat[i]: " << B_hat[a] << std::endl;
}
- printmatrix(std::cout, Q, "Matrix Q", "--");
- log << "Computed Matrix Q: " << std::endl;
- log << Q << std::endl;
+ log << "Computed prestrain B_hat: " << std::endl;
+ log << B_hat << std::endl;
- // compute B_hat
- for(size_t a = 0; a < 3; a++)
- {
- assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp2 ,B_Term); // <L i(M_alpha) + sym(grad phi_alpha), B >
- auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , B_Term); // <L i(x3G_alpha) , B>
- B_hat[a] = (coeffContainer[a]*tmp2) + GGterm;
- std::cout << "check this Contribution: " << (coeffContainer[a]*tmp2) << std::endl; //see orthotropic.tex
- std::cout <<"B_hat[i]: " << B_hat[a] << std::endl;
- }
- log << "Computed prestrain B_hat: " << std::endl;
- log << B_hat << std::endl;
+ // Compute effective Prestrain
+ FieldVector<double, 3> Beff;
+ Q.solve(Beff,B_hat);
+ log << "Computed prestrain B_eff: " << std::endl;
+ log << Beff << std::endl;
+ std::cout << "Beff : " << std::endl;
+ for(size_t i=0; i<3; i++)
+ {
+ std::cout <<"Beff[i]: " << Beff[i] << std::endl;
+ }
- // Compute effective Prestrain
- FieldVector<double, 3> Beff;
- Q.solve(Beff,B_hat);
- log << "Computed prestrain B_eff: " << std::endl;
- log << Beff << std::endl;
- std::cout << "Beff : " << std::endl;
- for(size_t i=0; i<3; i++)
- {
- std::cout <<"Beff[i]: " << Beff[i] << std::endl;
- }
+ 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;
+ std::cout << "Gamma: " << gamma << std::endl;
+ std::cout << "Number of Elements in each direction: " << nElements << std::endl;
+ log << "computed q1: " << q1_c << std::endl;
+ log << "computed q2: " << q2_c << std::endl;
+ log << "computed q3: " << q3_c << std::endl;
+ log << "computed b1: " << Beff[0] << std::endl;
+ log << "computed b2: " << Beff[1] << std::endl;
+ log << "computed b3: " << Beff[2] << std::endl;
+ log << "computed b1_hat: " << B_hat[0] << std::endl;
+ log << "computed b2_hat: " << B_hat[1] << std::endl;
+ log << "computed b3_hat: " << B_hat[2] << std::endl;
- 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;
- std::cout << "Gamma: " << gamma << std::endl;
- std::cout << "Number of Elements in each direction: " << nElements << std::endl;
- log << "computed q1: " << q1_c << std::endl;
- log << "computed q2: " << q2_c << std::endl;
- log << "computed q3: " << q3_c << std::endl;
- log << "computed b1: " << Beff[0] << std::endl;
- log << "computed b2: " << Beff[1] << std::endl;
- log << "computed b3: " << Beff[2] << std::endl;
- log << "computed b1_hat: " << B_hat[0] << std::endl;
- 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 = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
- double b2 = -(theta/4.0)*mu2;
- double b3 = 0.0;
- double q1 = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
- double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
-
- std::cout << " --- print analytic solutions --- " << std::endl;
- std::cout << "b1 : " << b1 << std::endl;
- std::cout << "b2 : " << b2 << std::endl;
- std::cout << "b3 : " << b3 << std::endl;
- std::cout << "q1 : " << q1 << std::endl;
- std::cout << "q2 : " << q2 << std::endl;
- std::cout << "q3 should be between q1 and q2" << std::endl;
- log << " --- analytic solutions: --- " << std::endl;
- log << "b1 : " << b1 << std::endl;
- log << "b2 : " << b2 << std::endl;
- log << "b3 : " << b3 << std::endl;
- log << "q1 : " << q1 << std::endl;
- log << "q2 : " << q2 << std::endl;
- log << "q3 should be between q1 and q2" << std::endl;
-
-
- 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}};
- };
-
- auto zeroFunction = [] (const Domain& x) {
- return MatrixRT{{0.0 , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
- };
+ //////////////////////////////////////////////////////////////
+ // 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 = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
+ double b2 = -(theta/4.0)*mu2;
+ double b3 = 0.0;
+ double q1 = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
+ double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
+
+ std::cout << " --- print analytic solutions --- " << std::endl;
+ std::cout << "b1 : " << b1 << std::endl;
+ std::cout << "b2 : " << b2 << std::endl;
+ std::cout << "b3 : " << b3 << std::endl;
+ std::cout << "q1 : " << q1 << std::endl;
+ std::cout << "q2 : " << q2 << std::endl;
+ std::cout << "q3 should be between q1 and q2" << std::endl;
+ log << " --- analytic solutions: --- " << std::endl;
+ log << "b1 : " << b1 << std::endl;
+ log << "b2 : " << b2 << std::endl;
+ log << "b3 : " << b3 << std::endl;
+ log << "q1 : " << q1 << std::endl;
+ log << "q2 : " << q2 << std::endl;
+ log << "q3 should be between q1 and q2" << std::endl;
+
+
+ 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}};
+ };
+
+ auto zeroFunction = [] (const Domain& x) {
+ return MatrixRT{{0.0 , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ };
+
+
+ if(write_L2Error)
+ {
+
+ std::cout << " ----- TEST ----" << std::endl;
+
+ std::cout << "computeL2SymErrorNew: " << computeL2SymErrorNew(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
+ std::cout << "computeL2SymErrorSPTest: " << computeL2SymErrorSPTest(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
+ std::cout << " -----------------" << std::endl;
+
+ auto L2SymError = computeL2SymError(Basis_CE,phi_1,gamma,symPhi_1_analytic);
+ std::cout << "L2-SymError: " << L2SymError << std::endl;
+
+ // auto L2errorTest = computeL2ErrorTest(Basis_CE,phi_1,gamma,symPhi_1_analytic);
+ // std::cout << "L2-ErrorTEST: " << L2errorTest << std::endl;
+
+ auto L2Norm_Symphi = computeL2SymError(Basis_CE,phi_1,gamma,zeroFunction); // Achtung Norm SymPhi_1
+ std::cout << "L2-Norm(Symphi_1) w zerofunction: " << L2Norm_Symphi << std::endl; // TODO Not Needed
+
+
+ std::cout << "L2-Norm(SymPhi_1): " << computeL2SymNormCoeff(Basis_CE,phi_1,gamma) << std::endl; // does not make a difference
+
+ VectorCT zeroVec;
+ zeroVec.resize(Basis_CE.size());
+ zeroVec = 0;
+
+ auto L2Norm_SymAnalytic = computeL2SymError(Basis_CE,zeroVec ,gamma,symPhi_1_analytic);
+ std::cout << "L2-Norm(Symanalytic): " << L2Norm_SymAnalytic << std::endl;
+
+ log << "L2-Error (symmetric Gradient phi_1):" << L2SymError << std::endl;
+
+ //TEST
+ std::cout << "L2Norm(phi_1) - GVfunc: " << computeL2NormFunc(Basis_CE,phi_1) << std::endl;
+ std::cout << "L2Norm(phi_1) - coeff: " << computeL2NormCoeff(Basis_CE,phi_1) << std::endl;
+ std::cout << "L2Norm(phi_1) - coeffV2: " << computeL2NormCoeffV2(Basis_CE,phi_1) << std::endl;
+ std::cout << "L2Norm(phi_1) - coeffV3: " << computeL2NormCoeffV3(Basis_CE,phi_1) << std::endl;
+
+ std::cout << "L2ErrorOld:" << computeL2ErrorOld(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
+ }
+
- if(write_L2Error)
- {
-
- std::cout << " ----- TEST ----" << std::endl;
-
- std::cout << "computeL2SymErrorNew: " << computeL2SymErrorNew(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
- std::cout << "computeL2SymErrorSPTest: " << computeL2SymErrorSPTest(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
- std::cout << " -----------------" << std::endl;
-
- auto L2SymError = computeL2SymError(Basis_CE,phi_1,gamma,symPhi_1_analytic);
- std::cout << "L2-SymError: " << L2SymError << std::endl;
-
- // auto L2errorTest = computeL2ErrorTest(Basis_CE,phi_1,gamma,symPhi_1_analytic);
- // std::cout << "L2-ErrorTEST: " << L2errorTest << std::endl;
-
- auto L2Norm_Symphi = computeL2SymError(Basis_CE,phi_1,gamma,zeroFunction); // Achtung Norm SymPhi_1
- std::cout << "L2-Norm(Symphi_1) w zerofunction: " << L2Norm_Symphi << std::endl; // TODO Not Needed
-
-
- std::cout << "L2-Norm(SymPhi_1): " << computeL2SymNormCoeff(Basis_CE,phi_1,gamma) << std::endl; // does not make a difference
-
- VectorCT zeroVec;
- zeroVec.resize(Basis_CE.size());
- zeroVec = 0;
- auto L2Norm_SymAnalytic = computeL2SymError(Basis_CE,zeroVec ,gamma,symPhi_1_analytic);
- std::cout << "L2-Norm(Symanalytic): " << L2Norm_SymAnalytic << std::endl;
- log << "L2-Error (symmetric Gradient phi_1):" << L2SymError << std::endl;
- //TEST
- std::cout << "L2Norm(phi_1) - GVfunc: " << computeL2NormFunc(Basis_CE,phi_1) << std::endl;
- std::cout << "L2Norm(phi_1) - coeff: " << computeL2NormCoeff(Basis_CE,phi_1) << std::endl;
- std::cout << "L2Norm(phi_1) - coeffV2: " << computeL2NormCoeffV2(Basis_CE,phi_1) << std::endl;
- std::cout << "L2Norm(phi_1) - coeffV3: " << computeL2NormCoeffV3(Basis_CE,phi_1) << std::endl;
-
- std::cout << "L2ErrorOld:" << computeL2ErrorOld(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
- }
-
//////////////////////////////////////////////////////////////////////////////////////////////
- // Write Data for Matlab / Optimization
+ // Write Data to Matlab / Optimization
//////////////////////////////////////////////////////////////////////////////////////////////
-
-
+
+
writeMatrixToMatlab(Q, "matlab.txt");
-
+
//////////////////////////////////////////////////////////////////////////////////////////////
// Write result to VTK file
//////////////////////////////////////////////////////////////////////////////////////////////
@@ -2859,4 +2924,7 @@ int main(int argc, char *argv[])
std::cout << "wrote data to file: CellProblem-result" << std::endl; // better with string for output name..
log.close();
+
+ }
+
}