diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index c9bb6f88610569cf19654a4ae323481422aa98d2..f88f61e0bcb6dd9a71ae5bb552bfb6bf991334d5 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -23,6 +23,7 @@
#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>
@@ -105,6 +106,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// printmatrix(std::cout, basisContainer[2] , "G_3", "--");
////////////////////////////////////////////////////
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
@@ -496,7 +498,7 @@ void computeElementLoadVector( const LocalView& localView,
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};
-
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1); // für simplex
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
@@ -537,7 +539,7 @@ void computeElementLoadVector( const LocalView& localView,
MatrixRT stressV(0);
stressV.axpy(2*mu(quadPos), strainV); //this += 2mu *strainU
- double trace = 0;
+ double trace = 0.0;
for (int ii=0; ii<nCompo; ii++)
trace += strainV[ii][ii];
@@ -546,7 +548,7 @@ void computeElementLoadVector( const LocalView& localView,
// Local energy density: stress times strain
- double energyDensity = 0;
+ double energyDensity = 0.0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
energyDensity += stressV[ii][jj] * forceTerm(quadPos)[ii][jj];
@@ -583,15 +585,174 @@ void computeElementLoadVector( const LocalView& localView,
}
}
+///////////////////////////////////////////////// TEST ///////////////////////////////////////
+
+
+
+
+
+/*
+
+// Compute the source term for a single element
+// < L (sym[D_gamma*nabla phi_i] + M_i ), x_3G_alpha >
+template<class LocalView, class LocalFunction1, class LocalFunction2, class Vector, class Force>
+void computeElementLoadTEST( const LocalView& localView,
+ LocalFunction1& mu,
+ LocalFunction2& lambda,
+ const double gamma,
+ Vector& elementRhs,
+ const Force& forceTerm
+ )
+{
+ // | f*phi|
+ // | --- |
+ // | f*m |
+
+ using Element = typename LocalView::Element;
+ const auto element = localView.element();
+ const auto geometry = element.geometry();
+ constexpr int dim = Element::dimension;
+ constexpr int nCompo = dim;
+
+ using VectorRT = FieldVector< double, nCompo>;
+ using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+
+ // Set of shape functions for a single element
+ const auto& localFiniteElement= localView.tree().child(0).finiteElement();
+ const auto nSf = localFiniteElement.localBasis().size();
+ // const auto& localFiniteElement = localView.tree().finiteElement();
+
+
+ elementRhs.resize(localView.size() +3);
+ elementRhs = 0;
+
+
+ // LocalBasis-Offset
+ const int localPhiOffset = localView.size();
+
+ ///////////////////////////////////////////////
+ // 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};
+
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1); // für simplex
+ const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+ for (const auto& quadPoint : quad)
+ {
+ const FieldVector<double,dim>& quadPos = quadPoint.position();
+ const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+ const double integrationElement = geometry.integrationElement(quadPos);
+ std::vector<FieldMatrix<double,1,nCompo> > referenceGradients;
+ localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+ std::vector< VectorRT > gradients(referenceGradients.size());
+ for (size_t i=0; i< gradients.size(); i++)
+ jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+
+ auto strain1 = forceTerm(quadPos);
+
+ // "f*phi"-part
+ for (size_t i=0; i < nSf; i++)
+ for (size_t k=0; k < nCompo; k++)
+ {
+ // Deformation gradient of the ansatz basis function
+ MatrixRT defGradientV(0);
+ defGradientV[k][0] = gradients[i][0]; // Y
+ defGradientV[k][1] = gradients[i][1]; // X2
+ defGradientV[k][2] = (1.0/gamma)*gradients[i][2]; // X3
+ // Elastic Strain
+ MatrixRT strainV; //strainV never used?
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
+
+
+
+
+
+ // St.Venant-Kirchhoff stress
+ MatrixRT stressV(0);
+ stressV.axpy(2*mu(quadPos), strain1); //this += 2mu *strainU
+ double trace = 0.0;
+ for (int ii=0; ii<nCompo; ii++)
+ trace += strain1[ii][ii];
+
+ for (int ii=0; ii<nCompo; ii++)
+ stressV[ii][ii] += lambda(quadPos) * trace;
+
+
+ // Local energy density: stress times strain
+ double energyDensity = 0.0;
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ energyDensity += stressV[ii][jj] * strainV[ii][jj];
+
+
+ size_t row = localView.tree().child(k).localIndex(i);
+ elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
+ }
+
+
+ // "f*m"-part
+ for (size_t m=0; m<3; m++)
+ {
+
+ // St.Venant-Kirchhoff stress
+ MatrixRT stressG(0);
+ stressG.axpy(2*mu(quadPos), strain1); //this += 2mu *strainU
+
+ double traceG = 0;
+ for (int ii=0; ii<nCompo; ii++)
+ traceG += strain1[ii][ii];
+
+ for (int ii=0; ii<nCompo; ii++)
+ stressG[ii][ii] += lambda(quadPos) * traceG;
+
+ double energyDensityfG = 0;
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ energyDensityfG += stressG[ii][jj] * basisContainer[m][ii][jj];
+
+ elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
+ }
+
+ }
+}
+
+
+
+
+
+*/
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////
@@ -722,6 +883,8 @@ void assembleCellLoad(const Basis& basis,
// BlockVector<double> elementRhs;
BlockVector<FieldVector<double,1> > elementRhs;
computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, loadFunctional);
+// printvector(std::cout, elementRhs, "elementRhs", "--");
+// computeElementLoadTEST(localView, muLocal, lambdaLocal, gamma, elementRhs, loadFunctional);
// printvector(std::cout, elementRhs, "elementRhs", "--");
// LOAD ASSEMBLY
@@ -762,7 +925,7 @@ auto energy(const Basis& basis,
const MatrixFunction& matrixFieldFuncB)
{
- auto energy = 0.0;
+ double energy = 0.0;
constexpr int dim = 3;
constexpr int nCompo = 3;
@@ -776,6 +939,13 @@ auto energy(const Basis& basis,
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+
+
+
+ // TEST
+// FieldVector<double,3> testvector = {1.0 , 1.0 , 1.0};
+// printmatrix(std::cout, matrixFieldFuncB(testvector) , "matrixFieldB(testvector) ", "--");
+
for (const auto& e : elements(basis.gridView()))
@@ -792,12 +962,16 @@ auto energy(const Basis& basis,
auto geometry = e.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- int order = 2*(dim*localFiniteElement.localBasis().order()-1);
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), order);
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+3; // TEST
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
- for (size_t pt=0; pt < quad.size(); pt++) {
+ for (const auto& quadPoint : quad) {
+// for (size_t pt=0; pt < quad.size(); pt++) {
+
- const Domain& quadPos = quad[pt].position();
+ const FieldVector<double,dim>& quadPos = quadPoint.position();
+// const Domain& quadPos = quad[pt].position();
const double integrationElement = geometry.integrationElement(quadPos);
@@ -807,7 +981,7 @@ auto energy(const Basis& basis,
MatrixRT stressU(0);
stressU.axpy(2*muLocal(quadPos), strain1); //this += 2mu *strainU // eigentlich this += 2mu *strain1 ?
- double trace = 0;
+ double trace = 0.0;
for (int ii=0; ii<nCompo; ii++)
trace += strain1[ii][ii];
@@ -817,12 +991,13 @@ auto energy(const Basis& basis,
auto strain2 = matrixFieldB(quadPos);
// Local energy density: stress times strain
- double energyDensity = 0;
+ double energyDensity = 0.0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
energyDensity += stressU[ii][jj] * strain2[ii][jj];
- elementEnergy += energyDensity * quad[pt].weight() * integrationElement;
+// elementEnergy += energyDensity * quad[pt].weight() * integrationElement;
+ elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
}
energy += elementEnergy;
}
@@ -1114,8 +1289,10 @@ auto integralMean(const Basis& basis,
localView.bind(element);
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- int order = 2*(dim*localFiniteElement.localBasis().order()-1);
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), order);
+
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
for (const auto& quadPoint : quad)
{
@@ -1131,7 +1308,7 @@ auto integralMean(const Basis& basis,
localView.bind(element);
fun.bind(element);
- const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
int order = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), order);
@@ -1223,6 +1400,92 @@ auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>
// ---- TODO - MatrixEmbedding function (iota) ---- ?
+
+
+
+
+
+
+
+template<class Basis, class Vector, class MatrixFunction, class Domain> //TODO
+double evalSymGrad(const Basis& basis,
+ Vector& coeffVector,
+ const double gamma,
+ Domain& x // evaluation Point in reference coordinates
+ )
+{
+
+
+ constexpr int dim = 3;
+ constexpr int nCompo = 3;
+
+ auto localView = basis.localView();
+
+
+ using GridView = typename Basis::GridView;
+// using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+ using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+
+
+
+
+
+ for (const auto& element : elements(basis.gridView()))
+ {
+
+ localView.bind(element);
+ auto geometry = element.geometry();
+
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
+ const auto nSf = localFiniteElement.localBasis().size();
+
+
+
+ const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(x);
+
+
+ std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
+ localFiniteElement.localBasis().evaluateJacobian(x,
+ referenceGradients);
+
+ // Compute the shape function gradients on the grid element
+ std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
+ // std::vector< VectorRT> gradients(referenceGradients.size());
+
+ for (size_t i=0; i<gradients.size(); i++)
+ jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+
+ MatrixRT defGradientU(0);
+
+
+ for (size_t k=0; k < nCompo; k++)
+ for (size_t i=0; i < nSf; i++) //Error: these are local fcts!
+ {
+
+ size_t localIdx = localView.tree().child(k).localIndex(i);
+ size_t globalIdx = localView.index(localIdx);
+
+ // (scaled) Deformation gradient of the ansatz basis function
+ defGradientU[k][0] += coeffVector[globalIdx]* gradients[i][0]; // Y
+ defGradientU[k][1] += coeffVector[globalIdx]* gradients[i][1]; //X2
+ defGradientU[k][2] += coeffVector[globalIdx]*(1.0/gamma)*gradients[i][2]; //X3
+
+
+
+ // symmetric Gradient (Elastic Strains)
+ MatrixRT strainU(0);
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ {
+ strainU[ii][jj] = coeffVector[globalIdx] * 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
+ }
+
+
+ }
+ }
+
+}
@@ -1263,6 +1526,7 @@ double computeL2Error(const Basis& basis,
const auto nSf = localFiniteElement.localBasis().size();
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
@@ -1294,6 +1558,9 @@ double computeL2Error(const Basis& basis,
for (size_t k=0; k < nCompo; k++)
for (size_t i=0; i < nSf; i++)
{
+
+ size_t localIdx = localView.tree().child(k).localIndex(i);
+ size_t globalIdx = localView.index(localIdx);
// (scaled) Deformation gradient of the ansatz basis function
defGradientU[k][0] = gradients[i][0]; // Y
@@ -1307,17 +1574,17 @@ double computeL2Error(const Basis& basis,
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
{
- strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
+ strainU[ii][jj] = coeffVector[globalIdx] * 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
}
- size_t localIdx = localView.tree().child(k).localIndex(i);
- size_t globalIdx = localView.index(localIdx);
+
// 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+= std::pow(coeffVector[globalIdx]*strainU[ii][jj] - matrixField(quadPos)[ii][jj] ,2);
+ tmp+= std::pow(strainU[ii][jj] - matrixField(quadPos)[ii][jj],2);
+// tmp+= std::pow((coeffVector[globalIdx]*strainU[ii][jj]) - matrixField(quadPos)[ii][jj] ,2);
@@ -1398,6 +1665,7 @@ int main(int argc, char *argv[])
p1 = 1.0;
double alpha = 2;
p2 = alpha*1.0;
+
prestraintype = 2;
// prestraintype = 1;
@@ -1426,7 +1694,9 @@ int main(int argc, char *argv[])
log << "Number of Elements: " << nE << std::endl;
debugLog << "Number of Elements: " << nE << std::endl;
- std::array<int,dim> nElements={nE,nE,nE}; //#Elements in each direction
+ std::array<int,dim> nElements={nE,nE,nE};
+// printvector(std::cout, coeffT_1 , "coeffT_1", "--" );
+// printvector(std::cout, x_1 , "x_1", "--" );; //#Elements in each direction
CellGridType grid_CE(lower,upper,nElements); //Corrector problem Domain:
@@ -1457,7 +1727,7 @@ int main(int argc, char *argv[])
-// double beta = 1;
+// double beta = 1; //homogeneous case
double beta = 2;
@@ -1524,13 +1794,12 @@ int main(int argc, char *argv[])
Functions::Experimental::BasisFactory::PeriodicIndexSet periodicIndices;
int equivCounter = 0;
- int nNodes = 0;
+
// 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;
- nNodes++;
for (const auto& v2 : vertices(gridView_CE))
if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
{
@@ -1539,7 +1808,8 @@ int main(int argc, char *argv[])
}
}
std::cout << "equivCounter:" << equivCounter << std::endl;
- std::cout <<" Number ofNodes: " << nNodes << std::endl;
+ std::cout << "Number of Elements in each direction: " << nE << std::endl;
+ std::cout << "Number ofNodes: " << gridView_CE.size(dim) << std::endl;
auto perTest = periodicIndices.indexPairSet();
@@ -1557,7 +1827,6 @@ int main(int argc, char *argv[])
std::cout << "basis.dimension()" << Basis_CE.dimension() <<std::endl;
const int phiOffset = Basis_CE.size();
-
debugLog << "phiOffset: "<< phiOffset << std::endl;
@@ -1620,8 +1889,7 @@ int main(int argc, char *argv[])
// auto Tast = x3G_3(test);
// printmatrix(std::cout, Tast, "x3G_3", "--");
- /////////////////////////////////////////////////////////////////////// TODO
- // TODO : PrestrainImp.hh
+ ///////////////////////////////////////////////////////////////////////
// double theta = 0.5;
// double p1 = 1;
// double p2 = 2;
@@ -1644,8 +1912,6 @@ int main(int argc, char *argv[])
//////////////////////////////////////////////////////////////////////////////
-
-
///////////////////////////////////////////////
// Basis for R_sym^{2x2}
//////////////////////////////////////////////
@@ -1654,7 +1920,6 @@ int main(int argc, char *argv[])
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", "--");
@@ -1683,6 +1948,8 @@ int main(int argc, char *argv[])
// }
//
+
+
/////////////////////////////////////////////////////////
// Determine global indices of components for arbitrary (local) index TODO :separate Function! arbitraryComponentwiseIndices
/////////////////////////////////////////////////////////
@@ -1712,20 +1979,54 @@ int main(int argc, char *argv[])
// set one basis-function to zero
if(set_oneBasisFunction_Zero)
{
- setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3,row);
-// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
+ setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3, row);
+// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
+// printvector(std::cout, load_alpha1, "load_alpha1 after setOneBaseFunctionToZero", "--");
}
- ////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Compute solution
////////////////////////////
VectorCT x_1 = load_alpha1;
VectorCT x_2 = load_alpha2;
VectorCT x_3 = load_alpha3;
+
+// printvector(std::cout, load_alpha2, "load_alpha2 before SOLVER", "--" );
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // CG - SOLVER
+
+// 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);
+
+
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ // GMRES - SOLVER
+
+
// Turn the matrix into a linear operator
MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_CE);
@@ -1749,11 +2050,94 @@ int main(int argc, char *argv[])
solver.apply(x_1, load_alpha1, statistics);
solver.apply(x_2, load_alpha2, statistics);
solver.apply(x_3, load_alpha3, statistics);
+
+
+
+
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ // QR - SOLVER
+
+// std::cout << "We use QR solver.\n";
+// 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);
+ std::cout << "---------load_alpha2 AFTER SOLVER: -------------" << std::endl;
+// printvector(std::cout, load_alpha2, "load_alpha2 AFTER SOLVER", "--" );
+
+
+
+ ////////////////////////////////////////////////////////////////////////////////////
+
+
+
+ // Write solutions in logs // TODO --------------------------------------
+// if (parameterSet.get<int>("write_solutions_corrector_problems") == 1){
+// log << "\nSolution of Corrector problems:\n";
+//
+// auto sizeComp = x_a.size()/nCompo;
+// std::vector<VectorRT> x_a_vec(sizeComp);
+// std::vector<VectorRT> x_K1_vec(sizeComp);
+// std::vector<VectorRT> x_K2_vec(sizeComp);
+// std::vector<VectorRT> x_K3_vec(sizeComp);
+//
+// for (int i=0; i < x_a_vec.size(); i++){
+// x_a_vec[i] = VectorRT{x_a[i], x_a[i + sizeComp], x_a[i + 2*sizeComp]};
+// x_K1_vec[i] = VectorRT{x_K1[i], x_K1[i + sizeComp], x_K1[i + 2*sizeComp]};
+// x_K2_vec[i] = VectorRT{x_K2[i], x_K2[i + sizeComp], x_K2[i + 2*sizeComp]};
+// x_K3_vec[i] = VectorRT{x_K3[i], x_K3[i + sizeComp], x_K3[i + 2*sizeComp]};
+// }
+//
+// log << "\nxa:\n";
+// for (int i=0; i < x_a_vec.size(); i++)
+// log << i << ": " << x_a_vec[i] << std::endl;
+//
+// log << "\nxK1:\n";
+// for (int i=0; i < x_K1_vec.size(); i++)
+// log << i << ": " << x_K1_vec[i] << std::endl;
+//
+// log << "\nxK2:\n";
+// for (int i=0; i < x_K2_vec.size(); i++)
+// log << i << ": " << x_K2_vec[i] << std::endl;
+//
+// log << "\nxK3:\n";
+// for (int i=0; i < x_K3_vec.size(); i++)
+// log << i << ": " << x_K3_vec[i] << std::endl;
+//
+// /*
+// log << "\nxa:\n";
+// for (int i=0; i < x_a.size(); i++)
+// log << i << " " << x_a[i] << std::endl;
+//
+// log << "\nxK1:\n";
+// for (int i=0; i < x_K1.size(); i++)
+// log << i << " " << x_K1[i] << std::endl;
+//
+// log << "\nxK2:\n";
+// for (int i=0; i < x_K2.size(); i++)
+// log << i << " " << x_K2[i] << std::endl;
+//
+// log << "\nxK3:\n";
+// for (int i=0; i < x_K3.size(); i++)
+// log << i << " " << x_K3[i] << std::endl;
+// */
+// }
+
+
+ ////////////////////////////////////////////////////////////////////////////////////
// Extract phi_alpha & M_alpha coefficients
+ ////////////////////////////////////////////////////////////////////////////////////
+
VectorCT phi_1, phi_2, phi_3;
phi_1.resize(Basis_CE.size());
phi_1 = 0;
@@ -1762,13 +2146,7 @@ int main(int argc, char *argv[])
phi_3.resize(Basis_CE.size());
phi_3 = 0;
- // VectorCT m_1, m_2, m_3;
- // m_1.resize(3);
- // m_1 = 0;
- // m_2.resize(3);
- // m_2 = 0;
- // m_3.resize(3);
- // m_3 = 0;
+
FieldVector<double,3> m_1, m_2, m_3;
@@ -1791,7 +2169,7 @@ int main(int argc, char *argv[])
////////////////////////////////////////////////////////////////////////////
// Make a discrete function from the FE basis and the coefficient vector
- ////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////// // ERROR : HIER NOCH KEIN Integralmittelabgezogen!?!?!?!?!
using SolutionRange = FieldVector<double,dim>;
// using SolutionRange = double;
@@ -1813,6 +2191,16 @@ int main(int argc, char *argv[])
phi_3);
+
+
+
+
+
+
+
+
+
+
// assemble M_alpha's (actually iota(M_alpha) )
MatrixRT M_1(0), M_2(0), M_3(0);
@@ -1849,6 +2237,14 @@ int main(int argc, char *argv[])
subtractIntegralMean(Basis_CE, localPhi_1 , phi_1);
subtractIntegralMean(Basis_CE, localPhi_2 , phi_2);
subtractIntegralMean(Basis_CE, localPhi_3 , phi_3);
+
+ // AUCH VON X-Vectoren!!!
+// printvector(std::cout, x_1 , "x_1", "--" );
+ subtractIntegralMean(Basis_CE, localPhi_1 , x_1); //TEST
+ subtractIntegralMean(Basis_CE, localPhi_2 , x_2);
+ subtractIntegralMean(Basis_CE, localPhi_3 , x_3);
+
+// printvector(std::cout, x_1 , "x_1 after substract_integralMean", "--" );
////////////////////////////////////////////////////////////////////////
// CHECK INTEGRAL-MEAN:
auto AVFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
@@ -1867,18 +2263,79 @@ int main(int argc, char *argv[])
// // PRINT Corrector-Basis-Coefficients
-// printvector(std::cout, phi_1, "Corrector-Phi_1", "--");
-// printvector(std::cout, phi_2, "Corrector-Phi_2", "--");
-// printvector(std::cout, phi_3, "Corrector-Phi_3", "--");
+// printvector(std::cout, phi_1, "Coefficients Corrector-Phi_1", "--");
+// printvector(std::cout, phi_2, "Coefficients Corrector-Phi_2", "--");
+// printvector(std::cout, phi_3, "Coefficients Corrector-Phi_3", "--");
+
+
+
+
+ ////////////////////////////////////////////////////////////////////////////////
+ // REMAKE Corrector-Functions after substract_integralMean
+
+ auto correctorFunction_1New = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ Basis_CE,
+ phi_1);
+
+
+ auto correctorFunction_2New = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ Basis_CE,
+ phi_2);
+ auto correctorFunction_3New = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ Basis_CE,
+ phi_3);
+ // Extract phi_alpha & M_alpha coefficients
+ VectorCT coeffT_1,coeffT_2,coeffT_3;
+ coeffT_1.resize(Basis_CE.size()+3);
+ coeffT_1 = 0;
+ coeffT_2.resize(Basis_CE.size()+3);
+ coeffT_2 = 0;
+ coeffT_3.resize(Basis_CE.size()+3);
+ coeffT_3 = 0;
+
+ for(int i=0; i< Basis_CE.size() ; i++)
+ {
+ coeffT_1[i] = phi_1[i];
+ coeffT_2[i] = phi_2[i];
+ coeffT_3[i] = phi_3[i];
+ }
+
+ coeffT_1[Basis_CE.size()] = m_1[0];
+ coeffT_1[Basis_CE.size()+1] = m_1[1];
+ coeffT_1[Basis_CE.size()+2] = m_1[2];
+ coeffT_2[Basis_CE.size()] = m_2[0];
+ coeffT_2[Basis_CE.size()+1] = m_2[1];
+ coeffT_2[Basis_CE.size()+2] = m_2[2];
+ coeffT_3[Basis_CE.size()] = m_3[0];
+ coeffT_3[Basis_CE.size()+1] = m_3[1];
+ coeffT_3[Basis_CE.size()+2] = m_3[2];
+
+ // TEST
+// for(int i = 0; i<Basis_CE.size()+3 ; i++)
+// {
+// if(x_1[i] < 1e-14)
+// x_1[i] = 0.0;
+// if(x_2[i] < 1e-14)
+// x_2[i] = 0.0;
+// if(x_3[i] < 1e-14)
+// x_3[i] = 0.0;
+//
+// }
+
+
+// printvector(std::cout, coeffT_1 , "coeffT_1", "--" );
+// printvector(std::cout, x_1 , "x_1", "--" );
/////////////////////////////////////////////////////////
// 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<decltype(localPhi_1)*, 3> phiContainer= {&localPhi_1, &localPhi_2, &localPhi_3};
+ const std::array<VectorCT, 3> coeffContainer = {x_1, x_2, x_3};
+// const std::array<VectorCT, 3> coeffContainer = {coeffT_1, coeffT_2, coeffT_3}; // TEST
+ const std::array<VectorCT, 3> loadContainer = {load_alpha1, load_alpha2, load_alpha3};
+ const std::array<decltype(localPhi_1)*, 3> phiContainer= {&localPhi_1, &localPhi_2, &localPhi_3}; // Achtung hie wurde kein Integralmittelabgezogen!!!!!
const std::array<MatrixRT*, 3 > mContainer = {&M_1 , &M_2, & M_3};
@@ -1895,18 +2352,25 @@ int main(int argc, char *argv[])
// 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) >
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp1 ,x3MatrixBasis[b]); // <L i(M_alpha) + sym(grad phi_alpha), i(x3G_beta) >
+ //entspricht load_alpha1,2,3..
+
+// printvector(std::cout, tmp1, "tmp1 ", "--");
+// printvector(std::cout, load_alpha1, "load_alpha1", "--" );
+// printvector(std::cout, loadContainer[b], "loadContainer[b]", "--" );
+
+ double GGterm = 0.0;
+ GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) > verändert man hier x3MatrixBasis????? TODO
- auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
+ Q[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric... check positiv definitness?
+// Q[a][b] = (coeffContainer[a]*loadContainer[b]) + GGterm; // TEST
- Q[a][b] = coeffContainer[a]*tmp1 + GGterm; // seems symmetric... check positiv definitness?
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", "--");
@@ -1918,7 +2382,7 @@ int main(int argc, char *argv[])
auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , B_Term); // <L i(x3G_alpha) , B>
- B_hat[a] = coeffContainer[a]*tmp2 + GGterm;
+ B_hat[a] = (coeffContainer[a]*tmp2) + GGterm;
std::cout <<"B_hat[i]: " << B_hat[a] << std::endl;
}
@@ -2012,60 +2476,76 @@ int main(int argc, char *argv[])
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}};
+ };
+
FieldVector<double,3> testvector = {1.0/4.0 , 0.0 , 0.25};
- std::cout << "t[2]: " << testvector[2] << std::endl;
- std::cout << "muTerm(t):" << muTerm(testvector) << std::endl;
- auto Teest = symPhi_1_analytic(testvector);
- printmatrix(std::cout, Teest, "symPhi_1_analytic(t)", "--");
+// std::cout << "t[2]: " << testvector[2] << std::endl;
+// std::cout << "muTerm(t):" << muTerm(testvector) << std::endl;
+// auto Teest = symPhi_1_analytic(testvector);
+// printmatrix(std::cout, Teest, "symPhi_1_analytic(t)", "--");
auto L2error = computeL2Error(Basis_CE,phi_1,gamma,symPhi_1_analytic);
std::cout << "L2-Error: " << L2error << std::endl;
-
+
+ auto L2Norm_Symphi = computeL2Error(Basis_CE,phi_1,gamma,zeroFunction); // Achtung Norm SymPhi_1
+ std::cout << "L2-Norm(Symphi_1): " << L2Norm_Symphi << std::endl; // TODO Not Needed
+
+ VectorCT zeroVec;
+ zeroVec.resize(Basis_CE.size());
+ zeroVec = 0;
+
+ auto L2Norm_analytic = computeL2Error(Basis_CE,zeroVec ,gamma,symPhi_1_analytic);
+ std::cout << "L2-Norm(analytic): " << L2Norm_analytic << std::endl;
//////////////////////////////////////////////////////////////////////////////////////////////
// Write result to VTK file
//////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////// TEST //////////////////////////////////////
- auto TestBasis = makeBasis(
- gridView_CE,
- power<dim>(
- lagrange<1>(),
- flatLexicographic()));
-
-
- auto& SubPreBasis = Basis_CE.preBasis().subPreBasis();
-
- VectorCT fullcoeff(TestBasis.size());
-
- std::cout << "TestBasis.size(): "<< TestBasis.size() << std::endl;
- auto testlocalView = TestBasis.localView();
- auto perlocalView = Basis_CE.localView();
+// auto TestBasis = makeBasis(
+// gridView_CE,
+// power<dim>(
+// lagrange<1>(),
+// flatLexicographic()));
+//
+//
+// auto& SubPreBasis = Basis_CE.preBasis().subPreBasis();
+//
+// VectorCT fullcoeff(TestBasis.size());
+//
+// std::cout << "TestBasis.size(): "<< TestBasis.size() << std::endl;
+// auto testlocalView = TestBasis.localView();
+// auto perlocalView = Basis_CE.localView();
// std::cout << "testlocalView.size(): " << testlocalView.size() << std::endl;
// std::cout << "perlocalView.size(): " << perlocalView.size() << std::endl;
- for (const auto& element : elements(gridView_CE))
- {
- testlocalView.bind(element);
- perlocalView.bind(element);
-// std::cout << "testlocalView.size(): " << testlocalView.size() << std::endl;
-// std::cout << "perlocalView.size(): " << perlocalView.size() << std::endl;
-
- for(size_t i=0; i<testlocalView.size(); i++)
- {
- auto testIdx = testlocalView.index(i);
- auto perIdx = perlocalView.index(i);
-// std::cout << "testIdx: " << testIdx << std::endl;
-// std::cout << "perIdx: " << perIdx << std::endl;
- fullcoeff[testIdx] = phi_1[perIdx];
- }
-
- }
+// for (const auto& element : elements(gridView_CE))
+// {
+// testlocalView.bind(element);
+// perlocalView.bind(element);
+// // std::cout << "testlocalView.size(): " << testlocalView.size() << std::endl;
+// // std::cout << "perlocalView.size(): " << perlocalView.size() << std::endl;
+//
+// for(size_t i=0; i<testlocalView.size(); i++)
+// {
+// auto testIdx = testlocalView.index(i);
+// auto perIdx = perlocalView.index(i);
+// // std::cout << "testIdx: " << testIdx << std::endl;
+// // std::cout << "perIdx: " << perIdx << std::endl;
+// fullcoeff[testIdx] = phi_1[perIdx];
+// }
+//
+// }
+//
+//
+//
// for (size_t i = 0; i < TestBasis.size()/3; i++)
// {
@@ -2083,31 +2563,44 @@ int main(int argc, char *argv[])
//
// }
//
- printvector(std::cout, fullcoeff, "fullcoeff" , "--");
- printvector(std::cout, phi_1, "phi_1" , "--");
+
- for (auto itr = perTest.begin(); itr != perTest.end(); itr++)
- {
-// std::cout << (*itr).first << std::endl;
-// std::cout << (*itr).second << std::endl;
-// fullcoeff[(*itr).first] = (*itr).second;
- }
-
-
- auto correctorFunction_Test = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
- TestBasis,
- fullcoeff);
- //////////////////////////////////////////////////////////////////////////////////////////////
+// printvector(std::cout, fullcoeff, "fullcoeff" , "--");
+// printvector(std::cout, phi_1, "phi_1" , "--");
+
+
+// auto correctorFunction_Test = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+// TestBasis,
+// fullcoeff);
+
+ //////////////////////////////////////////////////////////////////////////////////////////////
// SubsamplingVTKWriter<GridView> vtkWriter(
// gridView,
// refinementLevels(2));
+
+
VTKWriter<GridView> vtkWriter(gridView_CE);
- vtkWriter.addVertexData(
- correctorFunction_Test,
- VTK::FieldInfo("corrector Test", VTK::FieldInfo::Type::vector, dim));
+
+
+// vtkWriter.addVertexData(
+// correctorFunction_Test,
+// VTK::FieldInfo("corrector Test", VTK::FieldInfo::Type::vector, dim));
+
+
+
+
+ vtkWriter.addVertexData( //TEST
+ correctorFunction_1New,
+ VTK::FieldInfo("corrector phi_1New", VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.addVertexData(
+ correctorFunction_2New,
+ VTK::FieldInfo("corrector phi_2New", VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.addVertexData(
+ correctorFunction_3New,
+ VTK::FieldInfo("corrector phi_3New", VTK::FieldInfo::Type::vector, dim));
vtkWriter.addVertexData(
@@ -2119,10 +2612,8 @@ int main(int argc, char *argv[])
vtkWriter.addVertexData(
correctorFunction_3,
VTK::FieldInfo("corrector phi_3", VTK::FieldInfo::Type::vector, dim));
-//
- // vtkWriter.addVertexData(
- // solutionFunction,
- // VTK::FieldInfo("solution", VTK::FieldInfo::Type::scalar, 1));
+
+
vtkWriter.write("CellProblem-result");
std::cout << "wrote data to file: CellProblem-result" << std::endl; // better with string for output name..