diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 2510c3a29d504adbc31bfb758c079971c818b667..19e2c7f63b985a1b0e33ff133c334538aa9a68ef 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -585,11 +585,10 @@ void computeElementLoadVector( const LocalView& localView,
}
}
-///////////////////////////////////////////////// TEST ///////////////////////////////////////
-
+///////////////////////////////////////////////// TEST ///////////////////////////////////////
/*
@@ -1012,27 +1011,6 @@ auto energy(const Basis& basis,
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
template<class Basis, class Matrix, class Vector>
void setOneBaseFunctionToZero(const Basis& basis,
Matrix& stiffnessMatrix,
@@ -1270,10 +1248,6 @@ auto subtractIntegralMean(const Basis& basis,
// Check whether two points are equal on R/Z x R/Z x R
auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>& y)
{
-// std::cout << ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
-// and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1]))
-// and (FloatCmp::eq(x[2],y[2])) ) << std::endl;
-
return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1]))
and (FloatCmp::eq(x[2],y[2]))
@@ -1287,13 +1261,12 @@ 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,
@@ -1373,14 +1346,14 @@ double evalSymGrad(const Basis& basis,
}
}
-}
+}*/
template<class Basis, class Vector, class MatrixFunction>
-double computeL2ErrorTest(const Basis& basis,
+double computeL2Error(const Basis& basis,
Vector& coeffVector,
const double gamma,
const MatrixFunction& matrixFieldFunc
@@ -1522,9 +1495,9 @@ double computeL2ErrorTest(const Basis& basis,
-
+/*
template<class Basis, class Vector, class MatrixFunction>
-double computeL2Error(const Basis& basis,
+double computeL2ErrorOld(const Basis& basis,
Vector& coeffVector,
const double gamma,
const MatrixFunction& matrixFieldFunc
@@ -1614,7 +1587,7 @@ double computeL2Error(const Basis& basis,
}
}
return sqrt(error);
-}
+}*/
@@ -1699,13 +1672,8 @@ int main(int argc, char *argv[])
FieldVector<double,dim> upper({1.0, 1.0, 1.0/2.0});
}
-// int nE = 10;
-// log << "Number of Elements: " << nE << std::endl;
-// debugLog << "Number of Elements: " << nE << std::endl;
-
-// std::array<int,dim> nElements={nE,nE,nE};
std::array<int, dim> nElements = parameterSet.get<std::array<int,dim>>("nElements_Cell", {10, 10, 10});
-// std::array<int,dim> nElements={100,100,1}; // TEST
+
log << "Number of Elements in each direction: " << nElements << std::endl;
using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
@@ -1766,7 +1734,7 @@ int main(int argc, char *argv[])
log << "beta: " << beta << std::endl;
log << "material parameters: " << std::endl;
log << "mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
- log << " lambda1: " << lambda1 <<"\nlambda2:" << lambda2 << std::endl;
+ log << "lambda1: " << lambda1 <<"\nlambda2:" << lambda2 << std::endl;
@@ -1782,8 +1750,8 @@ int main(int argc, char *argv[])
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);
/////////////////////////////////////////////////////////
@@ -1804,7 +1772,7 @@ int main(int argc, char *argv[])
if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
{
periodicIndices.unifyIndexPair({gridView_CE.indexSet().index(v1)}, {gridView_CE.indexSet().index(v2)});
- equivCounter++;
+// equivCounter++;
}
}
// std::cout << "equivCounter:" << equivCounter << std::endl;
@@ -1824,7 +1792,7 @@ int main(int argc, char *argv[])
// std::cout << "size feBasis: " << Basis_CE.size() << std::endl;
// std::cout << "basis.dimension()" << Basis_CE.dimension() <<std::endl;
- log << "size feBasis: " << Basis_CE.size() << std::endl;
+ log << "size of FiniteElementBasis: " << Basis_CE.size() << std::endl;
const int phiOffset = Basis_CE.size();
debugLog << "phiOffset: "<< phiOffset << std::endl;
@@ -1849,10 +1817,6 @@ int main(int argc, char *argv[])
/////////////////////////////////////////////////////////
VectorCT load_alpha1, load_alpha2, load_alpha3;
MatrixCT stiffnessMatrix_CE;
- // auto rhsBackend = Functions::istlVectorBackend(b);
- // rhsBackend.resize(Basis_CE);
- // b = 0;
-
bool set_IntegralZero = parameterSet.get<bool>("set_IntegralZero", true);
@@ -1893,31 +1857,6 @@ int main(int argc, char *argv[])
Func2Tensor x3G_3neg = [x3G_3] (const Domain& x) {
return -1.0*x3G_3(x);
};
-// FieldVector<double,3> test= {1.0/4.0 , 0.0 , 0.25};
-// auto Tast = x3G_3(test);
-// printmatrix(std::cout, Tast, "x3G_3", "--");
-
- ///////////////////////////////////////////////////////////////////////
- // double theta = 0.5;
- // double p1 = 1;
- // double p2 = 2;
-
-// using std::abs;
-
-
-// Func2Tensor B = [] (const Domain& x) {
-// double theta = 0.5;
-// double p1 = 1;
-// double p2 = 2;
-// if (abs(x[0])>= (theta/2) && x[2] > 0)
-// return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
-// if (abs(x[0]) < (theta/2) && x[2] < 0)
-// return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
-// else
-// return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
-//
-// };
- //////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////
@@ -1941,18 +1880,7 @@ int main(int argc, char *argv[])
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha2 ,x3G_2neg);
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3neg);
-
- // TEST multiply RHS with -1 :
-// for(int i=0; i< Basis_CE.size()+3 ; i++)
-// {
-// load_alpha1[i] = (-1.0)*load_alpha1[i];
-// load_alpha2[i] = (-1.0)*load_alpha2[i];
-// load_alpha3[i] = (-1.0)*load_alpha3[i];
-// }
-
-
-
-// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix before B.C", "--");
+// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix", "--");
// printvector(std::cout, load_alpha1, "load_alpha1", "--");
@@ -2074,13 +2002,10 @@ int main(int argc, char *argv[])
sPQR.apply(x_3, load_alpha3, statisticsQR);
}
-// printvector(std::cout, load_alpha1BS, "load_alpha1 before SOLVER", "--" );
+// printvector(std::cout, load_alpha1BS, "load_alpha1 before SOLVER", "--" );
// printvector(std::cout, load_alpha1, "load_alpha1 AFTER SOLVER", "--" );
- std::cout << "---------load_alpha2 AFTER SOLVER: -------------" << std::endl;
// printvector(std::cout, load_alpha2, "load_alpha2 AFTER SOLVER", "--" );
-
-
////////////////////////////////////////////////////////////////////////////////////
// Extract phi_alpha & M_alpha coefficients
@@ -2131,9 +2056,6 @@ int main(int argc, char *argv[])
log << "Corrector-Matrix M_2: \n" << M_2 << std::endl;
log << " --------------------" << std::endl;
log << "Corrector-Matrix M_3: \n" << M_3 << std::endl;
-
-
-
////////////////////////////////////////////////////////////////////////////
@@ -2141,11 +2063,14 @@ int main(int argc, char *argv[])
//////////////////////////////////////////////////////////////////////////// // ERROR
using SolutionRange = FieldVector<double,dim>;
- std::cout << "check integralMean phi_1: " << std::endl;
- auto A = integralMean(Basis_CE,phi_1);
- for(size_t i=0; i<3; i++)
+ if(write_IntegralMean)
{
- std::cout << "Integral-Mean : " << A[i] << std::endl;
+ 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 : " << A[i] << std::endl;
+ }
}
if(substract_integralMean)
{
@@ -2153,29 +2078,22 @@ int main(int argc, char *argv[])
subtractIntegralMean(Basis_CE, phi_1);
subtractIntegralMean(Basis_CE, phi_2);
subtractIntegralMean(Basis_CE, phi_3);
-
-// printvector(std::cout, x_1 , "x_1", "--" );
subtractIntegralMean(Basis_CE, x_1);
subtractIntegralMean(Basis_CE, x_2);
subtractIntegralMean(Basis_CE, x_3);
-// printvector(std::cout, x_1 , "x_1 after substract_integralMean", "--" );
- ////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////
// CHECK INTEGRAL-MEAN:
- auto A = integralMean(Basis_CE, phi_1);
- for(size_t i=0; i<3; i++)
+ //////////////////////////////////////////
+ if(write_IntegralMean)
{
- std::cout << "Integral-Mean (CHECK) : " << A[i] << 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;
+ }
}
- ////////////////////////////////////////////////////
}
-// // PRINT Corrector-Basis-Coefficients
-// 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", "--");
-
-
-
/////////////////////////////////////////////////////////
// Write Solution (Corrector Coefficients) in Logs
/////////////////////////////////////////////////////////
@@ -2214,57 +2132,13 @@ int main(int argc, char *argv[])
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-10)
-// x_1[i] = 0.0;
-// if(x_2[i] < 1e-10)
-// x_2[i] = 0.0;
-// if(x_3[i] < 1e-10)
-// 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<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};
@@ -2273,30 +2147,22 @@ int main(int argc, char *argv[])
// 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++)
{
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) >
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 << "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", "--");
@@ -2306,7 +2172,6 @@ int main(int argc, char *argv[])
// 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>
@@ -2322,7 +2187,6 @@ int main(int argc, char *argv[])
// compute effective Prestrain
FieldVector<double, 3> Beff;
-
Q.solve(Beff,B_hat);
log << "Computed prestrain B_eff: " << std::endl;
log << Beff << std::endl;
@@ -2333,32 +2197,16 @@ int main(int argc, char *argv[])
std::cout <<"Beff[i]: " << Beff[i] << std::endl;
}
- // TODO
- // compute q1,q2,q3
- FieldVector<double ,3> g1 = {1.0 , 0 , 0};
- FieldVector<double ,3> g2 = {0 , 1.0 , 0};
- FieldVector<double ,3> g3 = {0 , 0, 1.0};
- FieldVector<double ,3> tmp_1;
- FieldVector<double ,3> tmp_2;
- FieldVector<double ,3> tmp_3;
- // auto tmp = g1- B_hat;
- // std::cout<< tmp << std::endl;
-
- Q.mv(g1,tmp_1);
- auto q1_c = tmp_1 * g1;
- std::cout<< "q1_c: " << q1_c << std::endl;
- Q.mv(g2,tmp_2);
- auto q2_c = tmp_2 * g2;
- std::cout<< "q2_c: " << q2_c << std::endl;
- Q.mv(g3,tmp_3);
- auto q3_c = tmp_3 * g3;
+
+ 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;
-
log << "computed q1: " << q1_c << std::endl;
log << "computed q2: " << q2_c << std::endl;
log << "computed q3: " << q3_c << std::endl;
@@ -2371,8 +2219,6 @@ int main(int argc, char *argv[])
-
- // TODO
//////////////////////////////////////////////////////////////
// Define Analytic Solutions
//////////////////////////////////////////////////////////////
@@ -2402,15 +2248,6 @@ int main(int argc, char *argv[])
log << "q3 should be between q1 and q2" << std::endl;
-
-
- // TODO Define sym grad phi_1
- // - how to compute <mu>_h ?
-
-// Func2Tensor symPhi_1_analytic = [] (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 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}};
};
@@ -2419,114 +2256,31 @@ int main(int argc, char *argv[])
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)", "--");
+ if(write_L2Error)
+ {
+ auto L2error = computeL2Error(Basis_CE,phi_1,gamma,symPhi_1_analytic);
+ std::cout << "L2-Error: " << L2error << std::endl;
+ // auto L2errorTest = computeL2ErrorTest(Basis_CE,phi_1,gamma,symPhi_1_analytic);
+ // std::cout << "L2-ErrorTEST: " << L2errorTest << 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
- auto L2error = computeL2Error(Basis_CE,phi_1,gamma,symPhi_1_analytic);
- std::cout << "L2-Error: " << L2error << std::endl;
-
- auto L2errorTest = computeL2ErrorTest(Basis_CE,phi_1,gamma,symPhi_1_analytic);
- std::cout << "L2-ErrorTEST: " << L2errorTest << 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;
-
- log << "L2-Error (symmetric Gradient phi_1):" << L2errorTest << std::endl;
+ 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;
+
+ log << "L2-Error (symmetric Gradient phi_1):" << L2error << 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();
-
-// 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 (size_t i = 0; i < TestBasis.size()/3; i++)
-// {
-// // auto c = Basis_CE.indices
-// // auto c = Basis_CE.preBasis_.subPreBasis_.transformation_.mappedIdx_;
-// // // Basis_CE.transformIndex(i);
-//
-// std::cout << "i: "<< i << std::endl;
-// Dune::Functions::FlatMultiIndex<unsigned long> idx = {i};
-// SubPreBasis.transformIndex(idx);
-// std::cout << "idx: " << idx << std::endl;
-//
-// fullcoeff[i] = phi_1[idx];
-// fullcoeff[i] = phi_1[idx];
-//
-// }
-//
-
-
-
-// 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(