diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 5669c14ae37e35102276852b7cc59e9dfff464f7..c9bb6f88610569cf19654a4ae323481422aa98d2 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -3,7 +3,6 @@
#include <vector>
#include <fstream>
-
#include <iostream>
#include <dune/common/indices.hh>
#include <dune/common/bitsetvector.hh>
@@ -27,17 +26,13 @@
#include <dune/istl/preconditioners.hh>
#include <dune/istl/io.hh>
-
#include <dune/functions/functionspacebases/interpolate.hh>
-
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
-
-
#include <dune/functions/functionspacebases/periodicbasis.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
@@ -52,16 +47,12 @@
#include <dune/microstructure/prestrainimp.hh>
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
-
-
// #include <dune/fufem/discretizationerror.hh>
-
using namespace Dune;
-
template<class LocalView, class Matrix, class localFunction1, class localFunction2>
void computeElementStiffnessMatrix(const LocalView& localView,
Matrix& elementMatrix,
@@ -95,11 +86,11 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// LocalBasis-Offset
const int localPhiOffset = localView.size();
- const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
const auto nSf = localFiniteElement.localBasis().size();
-
-
+// std::cout << "localView.size(): " << localView.size() << std::endl;
+// std::cout << "nSf : " << nSf << std::endl;
///////////////////////////////////////////////
// Basis for R_sym^{2x2} // wird nicht als Funktion benötigt da konstant...
@@ -125,11 +116,9 @@ void computeElementStiffnessMatrix(const LocalView& localView,
const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
const auto integrationElement = geometry.integrationElement(quadPoint.position());
- // std::vector<FieldMatrix<double,1,dim> > referenceGradients; // old
std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(
- quadPoint.position(),
- referenceGradients);
+ localFiniteElement.localBasis().evaluateJacobian(quadPoint.position(),
+ referenceGradients);
// Compute the shape function gradients on the grid element
std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
@@ -152,7 +141,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
defGradientU[k][0] = gradients[i][0]; // Y
defGradientU[k][1] = gradients[i][1]; //X2
defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
-
+// printmatrix(std::cout, defGradientU , "defGradientU", "--");
// (scaled) Deformation gradient of the test basis function
MatrixRT defGradientV(0);
defGradientV[l][0] = gradients[j][0]; // Y
@@ -168,7 +157,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
// strainV[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // same ? genügt strainU
-
+// printmatrix(std::cout, strainU , "strainU", "--");
}
// St.Venant-Kirchhoff stress
@@ -235,14 +224,17 @@ void computeElementStiffnessMatrix(const LocalView& localView,
MatrixRT stressG(0);
stressG.axpy(2*mu(quadPos), basisContainer[m]); //this += 2mu *strainU
- double traceG = 0;
+ double traceG = 0.0;
for (int ii=0; ii<nCompo; ii++)
+ {
+// std::cout << basisContainer[m][ii][ii] << std::endl;
traceG += basisContainer[m][ii][ii];
+ }
for (int ii=0; ii<nCompo; ii++)
stressG[ii][ii] += lambda(quadPos) * traceG;
- double energyDensityGphi = 0;
+ double energyDensityGphi = 0.0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
energyDensityGphi += stressG[ii][jj] * strainV[ii][jj]; // "m*phi"-part
@@ -261,52 +253,52 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// "phi*m"-part
- // for (size_t k=0; k < nCompo; k++)
- // for (size_t i=0; i < nSf; i++)
- // {
- //
- // // (scaled) Deformation gradient of the ansatz basis function
- // MatrixRT defGradientU(0);
- // defGradientU[k][0] = gradients[i][0]; // Y
- // defGradientU[k][1] = gradients[i][1]; //X2
- // defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
- //
- //
- // // symmetric Gradient (Elastic Strains)
- // MatrixRT strainU;
- // 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
- // }
- //
- // // St.Venant-Kirchhoff stress
- // MatrixRT stressU(0);
- // stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
- //
- // double trace = 0;
- // for (int ii=0; ii<nCompo; ii++)
- // trace += strainU[ii][ii];
- //
- // for (int ii=0; ii<nCompo; ii++)
- // stressU[ii][ii] += lambda(quadPos) * trace;
- //
- // for( size_t n=0; n<3; n++)
- // {
- //
- // // < L sym[D_gamma*nabla phi_i], G_j >
- // double energyDensityPhiG = 0;
- // for (int ii=0; ii<nCompo; ii++)
- // for (int jj=0; jj<nCompo; jj++)
- // energyDensityPhiG += stressU[ii][jj] * basisContainer[n][ii][jj]; // "phi*m"-part
- //
- // size_t col = localView.tree().child(k).localIndex(i);
- //
- // elementMatrix[localPhiOffset+n][col] += energyDensityPhiG * quadPoint.weight() * integrationElement;
- //
- // }
- // }
-
+// for (size_t k=0; k < nCompo; k++)
+// for (size_t i=0; i < nSf; i++)
+// {
+//
+// // (scaled) Deformation gradient of the ansatz basis function
+// MatrixRT defGradientU(0);
+// defGradientU[k][0] = gradients[i][0]; // Y
+// defGradientU[k][1] = gradients[i][1]; //X2
+// defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
+//
+//
+// // symmetric Gradient (Elastic Strains)
+// MatrixRT strainU;
+// 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
+// }
+//
+// // St.Venant-Kirchhoff stress
+// MatrixRT stressU(0);
+// stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
+//
+// double trace = 0;
+// for (int ii=0; ii<nCompo; ii++)
+// trace += strainU[ii][ii];
+//
+// for (int ii=0; ii<nCompo; ii++)
+// stressU[ii][ii] += lambda(quadPos) * trace;
+//
+// for( size_t n=0; n<3; n++)
+// {
+//
+// // < L sym[D_gamma*nabla phi_i], G_j >
+// double energyDensityPhiG = 0;
+// for (int ii=0; ii<nCompo; ii++)
+// for (int jj=0; jj<nCompo; jj++)
+// energyDensityPhiG += stressU[ii][jj] * basisContainer[n][ii][jj]; // "phi*m"-part
+//
+// size_t col = localView.tree().child(k).localIndex(i);
+//
+// elementMatrix[localPhiOffset+n][col] += energyDensityPhiG * quadPoint.weight() * integrationElement;
+//
+// }
+// }
+//
@@ -319,14 +311,14 @@ void computeElementStiffnessMatrix(const LocalView& localView,
MatrixRT stressG(0);
stressG.axpy(2*mu(quadPos), basisContainer[m]); //this += 2mu *strainU
- double traceG = 0;
+ double traceG = 0.0;
for (int ii=0; ii<nCompo; ii++)
traceG += basisContainer[m][ii][ii];
for (int ii=0; ii<nCompo; ii++)
stressG[ii][ii] += lambda(quadPos) * traceG;
- double energyDensityGG = 0;
+ double energyDensityGG = 0.0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
@@ -460,7 +452,7 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb)
// Compute the source term for a single element
-// < L sym[D_gamma*nabla phi_i], x_3G_alpha >
+// < L (sym[D_gamma*nabla phi_i] + M_i ), x_3G_alpha >
template<class LocalView, class LocalFunction1, class LocalFunction2, class Vector, class Force>
void computeElementLoadVector( const LocalView& localView,
LocalFunction1& mu,
@@ -489,7 +481,7 @@ void computeElementLoadVector( const LocalView& localView,
// const auto& localFiniteElement = localView.tree().finiteElement();
- elementRhs.resize(localView.size() +3 );
+ elementRhs.resize(localView.size() +3);
elementRhs = 0;
@@ -509,7 +501,8 @@ void computeElementLoadVector( const LocalView& localView,
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
- for (const auto& quadPoint : quad) {
+ for (const auto& quadPoint : quad)
+ {
const FieldVector<double,dim>& quadPos = quadPoint.position();
@@ -587,11 +580,7 @@ void computeElementLoadVector( const LocalView& localView,
elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
}
-
-
-
}
-
}
@@ -622,6 +611,10 @@ void assembleCellStiffness(const Basis& basis,
occupationPattern.exportIdx(matrix);
matrix = 0;
+
+
+ std::fstream localdebugLog;
+ localdebugLog.open("../../outputs/debugLog.txt",std::ios::out);
auto localView = basis.localView();
@@ -641,8 +634,8 @@ void assembleCellStiffness(const Basis& basis,
// Dune::Matrix<double> elementMatrix;
Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
- // printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
-
+// printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
+// localdebugLog << elementMatrix << std::endl;
// std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
@@ -683,7 +676,7 @@ void assembleCellStiffness(const Basis& basis,
matrix[phiOffset+m][phiOffset+n] += elementMatrix[localPhiOffset+m][localPhiOffset+n];
}
-
+// printmatrix(std::cout, matrix, "StiffnessMatrix", "--");
}
@@ -729,6 +722,7 @@ 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", "--");
// LOAD ASSEMBLY
for (size_t p=0; p<localPhiOffset; p++)
@@ -743,6 +737,9 @@ void assembleCellLoad(const Basis& basis,
{
b[phiOffset+m] += elementRhs[localPhiOffset+m];
}
+
+
+// printvector(std::cout, b, "b", "--");
}
}
@@ -799,6 +796,7 @@ auto energy(const Basis& basis,
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), order);
for (size_t pt=0; pt < quad.size(); pt++) {
+
const Domain& quadPos = quad[pt].position();
const double integrationElement = geometry.integrationElement(quadPos);
@@ -1062,8 +1060,8 @@ auto childToIndexMap(const Basis& basis,
for(size_t j=0; j<nSf; j++)
{
- auto Localidx = localView.tree().child(k).localIndex(j);
- r.push_back(localView.index(Localidx));
+ auto Localidx = localView.tree().child(k).localIndex(j); // local indices
+ r.push_back(localView.index(Localidx)); // global indices
}
@@ -1200,12 +1198,21 @@ auto subtractIntegralMean(const Basis& basis,
+ //////////////////////////////////////////////////
+ // Infrastructure for handling periodicity
+ //////////////////////////////////////////////////
+
// Check whether two points are equal on R/Z x R/Z x R
auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>& y)
{
+// 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])));
+ and (FloatCmp::eq(x[2],y[2]))
+ );
};
@@ -1292,27 +1299,30 @@ double computeL2Error(const Basis& basis,
defGradientU[k][0] = gradients[i][0]; // Y
defGradientU[k][1] = gradients[i][1]; //X2
defGradientU[k][2] = (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] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
- }
-
-
- // Local energy density: stress times strain
- double tmp = 0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
- tmp+= std::pow(strainU[ii][jj] - matrixField(quadPos)[ii][jj] ,2);
+ // symmetric Gradient (Elastic Strains)
+ MatrixRT strainU(0);
+ 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
+ }
+
+ 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);
- error += tmp * quadPoint.weight() * integrationElement;
+
+ error += tmp * quadPoint.weight() * integrationElement;
+ }
}
}
@@ -1343,8 +1353,10 @@ int main(int argc, char *argv[])
// output setter
// -- not set up --
std::fstream log;
+ std::fstream debugLog;
// log.open("output2.txt", std::ios::out);
log.open("../../outputs/output.txt",std::ios::out);
+ debugLog.open("../../outputs/debugLog.txt",std::ios::out);
// log << "Writing this to a file. \n";
// log.close();
@@ -1357,7 +1369,7 @@ int main(int argc, char *argv[])
// Get Parameters/Data
///////////////////////////////////
- double gamma = parameterSet.get<double>("gamma",2); // ratio dimension reduction to homogenization
+ double gamma = parameterSet.get<double>("gamma",2.0); // ratio dimension reduction to homogenization
// Material parameter laminat
// double E1 = parameterSet.get<double>("E1", 17e6); //material1
@@ -1412,6 +1424,7 @@ int main(int argc, char *argv[])
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}; //#Elements in each direction
@@ -1500,6 +1513,7 @@ int main(int argc, char *argv[])
log << "\n lambda: " << lambda1 << std::endl;
+
/////////////////////////////////////////////////////////
@@ -1509,18 +1523,33 @@ 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)))
+ {
periodicIndices.unifyIndexPair({gridView_CE.indexSet().index(v1)}, {gridView_CE.indexSet().index(v2)});
-
+ equivCounter++;
+ }
+ }
+ std::cout << "equivCounter:" << equivCounter << std::endl;
+ std::cout <<" Number ofNodes: " << nNodes << std::endl;
+
+ auto perTest = periodicIndices.indexPairSet();
+
+
auto Basis_CE = makeBasis(
gridView_CE,
- power<dim>(
- // lagrange<1>(),
+ power<dim>( //lagrange<1>(),
Functions::Experimental::BasisFactory::periodic(lagrange<1>(), periodicIndices),
flatLexicographic())); //
+// blockedInterleaved())); // siehe Periodicbasistest.. funktioniert auch?
// flatInterleaved()));
@@ -1528,6 +1557,10 @@ 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;
+
+
// std::cout << Basis_CE.preBasis_.children << std::endl; // does not work
// std::cout << Basis_CE::preBasis_::children << std::endl;
@@ -1563,6 +1596,8 @@ int main(int argc, char *argv[])
bool substract_integralMean = true;
// bool set_oneBasisFunction_Zero = false;
// bool substract_integralMean = false;
+ debugLog << " set_oneBasisFunction_Zero: " << set_oneBasisFunction_Zero << std::endl;
+ debugLog << " substract_integralMean: " << substract_integralMean << std::endl;
@@ -1620,10 +1655,10 @@ int main(int argc, char *argv[])
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
- // printmatrix(std::cout, basisContainer[0] , "G_1", "--");
+// 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;
@@ -1635,8 +1670,9 @@ int main(int argc, char *argv[])
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha2 ,x3G_2);
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3);
-
- // printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix before B.C", "--");
+
+// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix before B.C", "--");
+// printvector(std::cout, load_alpha1, "load_alpha1", "--");
// set Integral to zero
@@ -1677,7 +1713,7 @@ int main(int argc, char *argv[])
if(set_oneBasisFunction_Zero)
{
setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3,row);
- // printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
+// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
}
@@ -1830,26 +1866,10 @@ int main(int argc, char *argv[])
// // PRINT Corrector-Basis-Coefficients
- // std::cout << "print coefficient-vector phi_1" << std::endl;
- // for(size_t i=0; i<Basis_CE.size();i++)
- // {
- // std::cout << phi_1[i] << std::endl;
- // }
- // std::cout << "print coefficient-vector phi_2" << std::endl;
- // for(size_t i=0; i<Basis_CE.size();i++)
- // {
- // std::cout << phi_2[i] << std::endl;
- // }
- // std::cout << "print coefficient-vector phi_3" << std::endl;
- // for(size_t i=0; i<Basis_CE.size();i++)
- // {
- // std::cout << phi_3[i] << std::endl;
- // }
- //
-
- // auto df = derivative(correctorFunction_1); // does not work :(
-
-
+
+// printvector(std::cout, phi_1, "Corrector-Phi_1", "--");
+// printvector(std::cout, phi_2, "Corrector-Phi_2", "--");
+// printvector(std::cout, phi_3, "Corrector-Phi_3", "--");
@@ -2010,11 +2030,86 @@ int main(int argc, char *argv[])
//////////////////////////////////////////////////////////////////////////////////////////////
// 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" , "--");
+
+
+ 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);
+ //////////////////////////////////////////////////////////////////////////////////////////////
+
+
// 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_1,
VTK::FieldInfo("corrector phi_1", VTK::FieldInfo::Type::vector, dim));
@@ -2024,6 +2119,7 @@ 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));