diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index f4ee1a89401ebc8ca9e9a2470b3b33f2b30bc045..187ae5a4a74d01ec3fbea38172f2eb87e83f3ecf 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -30,8 +30,7 @@
#include <dune/functions/functionspacebases/interpolate.hh>
-// never used:
-// #include <dune/functions/functionspacebases/taylorhoodbasis.hh>
+
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
@@ -83,9 +82,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// using Func2Tensor = std::function< MatrixRT(const Domain&) >;
-
-
- elementMatrix.setSize(localView.size()+3, localView.size()+3);
+ elementMatrix.setSize(localView.size()+3, localView.size()+3); //extend by dim ´R_sym^{2x2}
elementMatrix = 0;
@@ -104,8 +101,6 @@ void computeElementStiffnessMatrix(const LocalView& localView,
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, 0.5, 0.0}, {0.5, 0.0, 0.0}, {0.0, 0.0, 0.0}};
-
-
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
//print:
printmatrix(std::cout, basisContainer[0] , "G_1", "--");
@@ -120,7 +115,6 @@ void computeElementStiffnessMatrix(const LocalView& localView,
for (const auto& quadPoint : quad)
{
-
const auto& quadPos = quadPoint.position();
const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
const auto integrationElement = geometry.integrationElement(quadPoint.position());
@@ -347,58 +341,6 @@ void computeElementStiffnessMatrix(const LocalView& localView,
-
-
-
-
-
-
-
-
-
-// Set the occupation pattern of the stiffness matrix
-// template<class Basis, class Matrix>
-// void setOccupationPattern(const Basis& basis, Matrix& matrix) // OLD
-// {
-// enum {dim = Basis::GridView::dimension};
-// // MatrixIndexSets store the occupation pattern of a sparse matrix.
-// // They are not particularly efficient, but simple to use.
-// std::array<std::array<MatrixIndexSet, 2>, 2> nb;
-// // Set sizes of the 2x2 submatrices
-// for (size_t i=0; i<2; i++)
-// for (size_t j=0; j<2; j++)
-// nb[i][j].resize(basis.size({i}), basis.size({j}));
-//
-// // A view on the FE basis on a single element
-// auto localView = basis.localView();
-//
-// // Loop over all leaf elements
-// for(const auto& element : elements(basis.gridView()))
-// {
-// // Bind the local view to the current element
-// localView.bind(element);
-// // Add element stiffness matrix onto the global stiffness matrix
-// for (size_t i=0; i<localView.size(); i++)
-// {
-// // Global index of the i-th local degree of freedom of the current element
-// auto row = localView.index(i);
-// for (size_t j=0; j<localView.size(); j++ )
-// {
-// // Global index of the j-th local degree of freedom of the current element
-// auto col = localView.index(j);
-// nb[row[0]][col[0]].add(row[1],col[1]);
-// }
-// }
-// }
-//
-// // Give the matrix the occupation pattern we want.
-// using namespace Indices;
-// nb[0][0].exportIdx(matrix[_0][_0]);
-// nb[0][1].exportIdx(matrix[_0][_1]);
-// nb[1][0].exportIdx(matrix[_1][_0]);
-// nb[1][1].exportIdx(matrix[_1][_1]);
-// }
-
// Get the occupation pattern of the stiffness matrix
template<class Basis>
void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb)
@@ -455,7 +397,6 @@ 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 >
template<class LocalView, class LocalFunction1, class LocalFunction2, class Vector, class Force>
@@ -470,6 +411,10 @@ void computeElementLoadVector( const LocalView& localView,
// LocalView::Element::dimension>)> volumeTerm
)
{
+ // | f*phi|
+ // | --- |
+ // | f*m |
+
using Element = typename LocalView::Element;
// auto element = localView.element();
const auto element = localView.element();
@@ -486,19 +431,21 @@ void computeElementLoadVector( const LocalView& localView,
// const auto& localFiniteElement = localView.tree().finiteElement();
- elementRhs.resize(localView.size());
+ elementRhs.resize(localView.size() +3 );
elementRhs = 0;
-// double q = 1.0; // pressure load
+ // LocalBasis-Offset
+ const int localPhiOffset = localView.size();
+
+ ///////////////////////////////////////////////
+ // Basis for R_sym^{2x2} // wird nicht als Funktion benötigt da konstant...
+ //////////////////////////////////////////////
+ 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, 0.5, 0.0}, {0.5, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
-// for(size_t i=0; i<localFiniteElement.size(); i++)
-// {
-// size_t row = localView.tree().child(0).localIndex(i); // only 1.child ~ function-value dofs play a role
-//
-// elementRhs[row] = (element.geometry().volume()/(dim+1))*q;
-//
-// }
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1); // für simplex
// const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
@@ -506,6 +453,7 @@ void computeElementLoadVector( const LocalView& localView,
for (const auto& quadPoint : quad){
+
// const Domain& quadPos = quadPoint.position();
const FieldVector<double,dim>& quadPos = quadPoint.position();
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
@@ -518,46 +466,74 @@ void computeElementLoadVector( const LocalView& localView,
for (size_t i=0; i< gradients.size(); i++)
jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+ // "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, no: Domain and Range components are swaped
- 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
+ // 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 strainU, strainV; //strainV never used?
+ MatrixRT strainV; //strainV never used?
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]);
+ strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
// St.Venant-Kirchhoff stress
- MatrixRT stressU(0);
- stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
+ MatrixRT stressV(0);
+ stressV.axpy(2*mu(quadPos), strainV); //this += 2mu *strainU
double trace = 0;
for (int ii=0; ii<nCompo; ii++)
- trace += strainU[ii][ii];
+ trace += strainV[ii][ii];
for (int ii=0; ii<nCompo; ii++)
- stressU[ii][ii] += lambda(quadPos) * trace;
- //log_ << "stressU:\n" << stressU << std::endl;
+ stressV[ii][ii] += lambda(quadPos) * trace;
+
+
// Local energy density: stress times strain
double energyDensity = 0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
- energyDensity += stressU[ii][jj] * forceTerm(quadPos)[ii][jj];
- //log_ << "strainE:\n" << strainE(quadPoint.position()) << std::endl;
+ energyDensity += stressV[ii][jj] * forceTerm(quadPos)[ii][jj];
+
size_t row = localView.tree().child(k).localIndex(i);
elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
- //log_ << "\nk: " << k << "\ni: " << i << "\nrow: " << row << std::endl;
+ }
+
+
+ // "f*m"-part
+ for (size_t m=0; m<3; m++)
+ {
- //log << "energyDensity quadPoint.weight() integrationElement: " << energyDensity << " " << quadPoint.weight() << " " << integrationElement << std::endl;
- } //end for k
- } // end for quadPoint
+ // St.Venant-Kirchhoff stress
+ MatrixRT stressG(0);
+ stressG.axpy(2*mu(quadPos), basisContainer[m]); //this += 2mu *strainU
+
+ double traceG = 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 energyDensityfG = 0;
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ energyDensityfG += stressG[ii][jj] * forceTerm(quadPos)[ii][jj];
+
+ elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
+ }
+
+
+
+
+ }
}
@@ -574,14 +550,14 @@ void computeElementLoadVector( const LocalView& localView,
-template<class Basis, class Matrix, class LocalFunction1, class LocalFunction2, class Vector, class Force>
-void assembleCellProblem(const Basis& basis,
+template<class Basis, class Matrix, class LocalFunction1, class LocalFunction2>
+void assembleCellStiffness(const Basis& basis,
LocalFunction1& muLocal,
LocalFunction2& lambdaLocal,
const double gamma,
- Matrix& matrix,
- Vector& b,
- const Force& forceTerm
+ Matrix& matrix
+// Vector& b,
+// const Force& forceTerm
// const std::function<
// double(FieldVector<double,
// Basis::GridView::dimension>)
@@ -593,34 +569,21 @@ void assembleCellProblem(const Basis& basis,
MatrixIndexSet occupationPattern;
getOccupationPattern(basis, occupationPattern);
occupationPattern.exportIdx(matrix);
-// setOccupationPattern(basis, matrix);
-
matrix = 0;
- printmatrix(std::cout, matrix, "matrix built with a MatrixIndexSet", "--");
- b.resize(basis.size());
- b = 0;
auto localView = basis.localView();
-// std::cout << localView.maxSize() << std::endl;
-
const int phiOffset = basis.dimension();
-
- // Transform G_alpha's to GridViewFunctions/LocalFunctions -- why exactly?
- auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
- auto loadFunctional = localFunction(loadGVF); // Dune::Functions:: notwendig?
-
-
for (const auto& element : elements(basis.gridView()))
{
-
localView.bind(element);
muLocal.bind(element);
lambdaLocal.bind(element);
+
const int localPhiOffset = localView.size();
std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
@@ -638,9 +601,9 @@ void assembleCellProblem(const Basis& basis,
// printmatrix(std::cout, TestelementMatrix, "TESTElementMatrix", "--");
//////////////////////////////////////////////////////////////////////////////
- // GLOBAL ASSEMBLY ..(hier dreht sich i, j..)
+ // GLOBAL STIFFNES ASSEMBLY ..(hier dreht sich i, j..)
//////////////////////////////////////////////////////////////////////////////
- for (size_t i=0; i<localPhiOffset i++)
+ for (size_t i=0; i<localPhiOffset; i++)
{
// The global index of the i-th local degree of freedom of the current element
auto row = localView.index(i);
@@ -672,39 +635,211 @@ void assembleCellProblem(const Basis& basis,
+ }
+
+}
+template<class Basis, class LocalFunction1, class LocalFunction2, class Vector, class Force>
+void assembleCellLoad(const Basis& basis,
+ LocalFunction1& muLocal,
+ LocalFunction2& lambdaLocal,
+ const double gamma,
+ Vector& b,
+ const Force& forceTerm
+ )
+{
+
+ std::cout << "assemble load vector." << std::endl;
+ b.resize(basis.size());
+ b = 0;
+ auto localView = basis.localView();
+ const int phiOffset = basis.dimension();
- /*
+ // Transform G_alpha's to GridViewFunctions/LocalFunctions -- why exactly?
+ auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
+ auto loadFunctional = localFunction(loadGVF); // Dune::Functions:: notwendig?
+
+ for (const auto& element : elements(basis.gridView()))
+ {
+
+ localView.bind(element);
+ muLocal.bind(element);
+ lambdaLocal.bind(element);
+ loadFunctional.bind(element);
+
+ const int localPhiOffset = localView.size();
+ std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
- /////////////////////////////////////////////////////////////////////////////////////
- // Now get the local contribution to the right-hand side vector
-// BlockVector<double> localB;
+// BlockVector<double> elementRhs;
BlockVector<FieldVector<double,1>> elementRhs;
computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, loadFunctional);
- for (size_t p=0; p<elementRhs.size(); p++)
+
+ // LOAD ASSEMBLY
+ for (size_t p=0; p<localPhiOffset; p++)
{
// The global index of the p-th vertex of the element
auto row = localView.index(p);
-// b[row[0]][row[1]] += elementRhs[p];
+
b[row] += elementRhs[p];
}
- */
-
+ for (size_t m=0; m<3; m++ )
+ {
+ b[phiOffset+m] += elementRhs[localPhiOffset+m];
+ }
}
+}
+
+
+
+
+template<class Basis, class Matrix, class Vector>
+void setOneBaseFunctionToZero(const Basis& basis,
+ Matrix& stiffnessMatrix,
+ Vector& load_alpha1,
+ Vector& load_alpha2,
+ Vector& load_alpha3
+ )
+{
+
+}
+
+
+
+
+template<class Basis, class Matrix, class Vector>
+void setIntegralZero (const Basis& basis,
+ Matrix& stiffnessMatrix,
+ Vector& load_alpha1,
+ Vector& load_alpha2,
+ Vector& load_alpha3
+ )
+{
+ //
+ // 0. moment invariance, 3 dofs
+ //
+
+// if (parameterSet_.get<bool>("set_integral_0")){
+ std::cout << "\nassembling: integral to zero\n";
+
+
+ auto n = basis.size();
+// auto nNodeBasis = n/nCompo;
+
+ constexpr int dim = 3; // irgendwo herziehen?
+ constexpr int nCompo = dim;
+
+ auto localView = basis.localView();
+
+ int arbitraryIndex = 0; // sollte GlOBALER INDEX sein
+
+// unsigned long row1;
+ FieldVector<int,3> row; // fill with Indices..
+
+
+ //Determine 3 global indices (one for each component of an arbitrary FE-function).. oder mittels Offset?
+
+// for (int k = 0; k < dim; k++)
+// {
+//
+// auto row = localView.tree().child(k).localIndex(arbitraryIndex);
+// // auto row = perIndexMap_[arbitraryIndex + d*nNodeBasis];
+// stiffnessMatrix[row] = 0; // setzt gesamte Zeile auf 0 !!!
+// load_alpha1[row] = 0.0;
+// load_alpha2[row] = 0.0;
+// load_alpha3[row] = 0.0;
+// }
+
+ unsigned int elementCounter = 1;
+ for (const auto& element : elements(basis.gridView()))
+ {
+// if (elementCounter % 50 == 1)
+// std::cout << "integral = zero - treatment - element: " << elementCounter << std::endl;
+
+ localView.bind(element);
+
+
+ if(elementCounter == 1) // get arbitraryIndex(global) for each Component
+ {
+ for (int k = 0; k < dim; k++)
+ {
+ auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
+ row[k] = localView.index(rowLocal);
+ std::cout << "row[k]:" << row[k] << std::endl;
+ stiffnessMatrix[row[k]] = 0;
+ load_alpha1[row[k]] = 0.0;
+ load_alpha2[row[k]] = 0.0;
+ load_alpha3[row[k]] = 0.0;
+ }
+ }
+
+
+ // "local assembly"
+ BlockVector<FieldVector<double,1> > elementColumns;
+ elementColumns.resize(localView.size()+3);
+ elementColumns = 0;
+
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+ const auto nSf = localFiniteElement.localBasis().size();
+
+ int order = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), order);
+
+ for (const auto& quadPoint : quad)
+ {
+ // for (size_t pt=0; pt < quad.size(); pt++) {
+
+ const auto& quadPos = quadPoint.position();
+ // const FieldVector<double,dim>& quadPos = quad[pt].position();
+ const double integrationElement = element.geometry().integrationElement(quadPos);
+ const auto jacobian = element.geometry().jacobianInverseTransposed(quadPos);
+
+ std::vector<FieldVector<double,1> > shapeFunctionValues;
+ localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
+ // Compute the shape function on the grid element
+ std::vector<FieldVector<double,1> > BasisFunctionValues(shapeFunctionValues.size());
+ for (size_t i=0; i<BasisFunctionValues.size(); i++)
+ jacobian.mv(shapeFunctionValues[i], BasisFunctionValues[i]);
+
+ for (size_t i=0; i<localView.size(); i++)
+ for (size_t k=0; k<nCompo; k++)
+ {
+ size_t idx = localView.tree().child(k).localIndex(i);
+ elementColumns[idx] += BasisFunctionValues[i] * quadPoint.weight() * integrationElement; // bei Periodicbasis einfach mit child(k) arbeiten...
+ }
+
+
+
+ // "global assembly"
+ for (int k = 0; k<dim; k++)
+ for (size_t j=0; j<localView.size(); j++)
+ {
+// auto colLocal = localView.tree().child(d).localIndex(j);
+ auto col = localView.index(j);
+ stiffnessMatrix[row[k]][col] += elementColumns[j] * quadPoint.weight() * integrationElement;
+ }
+
+ elementCounter ++;
+ }//end of quad
+
+ }//end for each element
}
+
+
+
+
// ---- TODO - MatrixEmbedding function (iota) ---- ???
@@ -731,6 +866,9 @@ void assembleCellProblem(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)
{
@@ -806,7 +944,7 @@ 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});
- int nE = 2;
+ int nE = 1;
std::array<int,dim> nElements={nE,nE,nE}; //#Elements in each direction
@@ -922,6 +1060,10 @@ MatrixCT stiffnessMatrix_CE;
+bool set_integral_zero = true;
+bool set_oneBasisFunction_Zero = true;
+
+
// using Range = FieldVector<double,dim>;
// auto sourceTerm = [](const FieldVector<double,dim>& x){return Range{0.0, -1.0};};
@@ -929,20 +1071,17 @@ MatrixCT stiffnessMatrix_CE;
Func2Tensor x3G_1 = [] (const Domain& z) {
return MatrixRT{{1.0*z[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; };
-// assembleStiffnesMatrix(Basis_CE, muTerm, lambdaTerm, stiffnessMatrix_CE ,muLocal,lambdaLocal,b,sourceTerm);
+Func2Tensor x3G_2 = [] (const Domain& z) {
+ return MatrixRT{{0.0, 0.0, 0.0}, {0.0, 1.0*z[2], 0.0}, {0.0, 0.0, 0.0}}; };
+
+Func2Tensor x3G_3 = [] (const Domain& z) {
+ return MatrixRT{{0.0, 0.5*z[2], 0.0}, {0.5*z[2], 0.0, 0.0}, {0.0, 0.0, 0.0}}; };
+
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, 0.5, 0.0}, {0.5, 0.0, 0.0}, {0.0, 0.0, 0.0}};
-// printmatrix(std::cout, G_1, "G_1", "--");
-
-
-// std::vector<MatrixRT> basisContainer;
-//
-// basisContainer[0] = {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
-// basisContainer[1] = {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
-// basisContainer[2] = {{0.0, 1/2, 0.0}, {1/2, 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
@@ -954,13 +1093,41 @@ Func2Tensor x3G_1 = [] (const Domain& z) {
-assembleCellProblem(Basis_CE, muLocal, lambdaLocal,gamma, stiffnessMatrix_CE ,load_alpha1 ,x3G_1);
-
+assembleCellStiffness(Basis_CE, muLocal, lambdaLocal, gamma, stiffnessMatrix_CE);
+assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha1 ,x3G_1);
+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", "--");
+// set Integral to zero
+if(set_integral_zero)
+{
+ setIntegralZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3);
+ printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setIntegralZero", "--");
+}
+
+
+
+// set Integral to zero
+if(set_oneBasisFunction_Zero)
+{
+ setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3);
+ printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setIntegralZero", "--");
+}
+
+
+
+
+
+
+// setzt gesamte Zeile = 0:
+// stiffnessMatrix_CE[2] = 0;
+// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix MOD", "--");
+
+
@@ -1042,7 +1209,10 @@ printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix before B.C", "--");
// { stokes_solve_begin }
// Initial iterate: Start from the rhs vector,
// that way the Dirichlet entries are already correct.
-VectorCT x = load_alpha1;
+VectorCT x_1 = load_alpha1;
+VectorCT x_2 = load_alpha2;
+VectorCT x_3 = load_alpha3;
+
// Turn the matrix into a linear operator
MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_CE);
@@ -1063,9 +1233,14 @@ RestartedGMResSolver<VectorCT> solver(
// Object storing some statistics about the solving process
InverseOperatorResult statistics;
-// Solve!
-solver.apply(x, load_alpha1, statistics);
-// { stokes_solve_end }
+
+
+// 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);
+
+
////////////////////////////////////////////////////////////////////////////
// Make a discrete function from the FE basis and the coefficient vector
@@ -1074,9 +1249,14 @@ solver.apply(x, load_alpha1, statistics);
using SolutionRange = FieldVector<double,dim>;
// using SolutionRange = double;
-auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
Basis_CE,
- x);
+ x_1);
+
+
+
+
+
// use only first child of W_h to plot solution
// auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
// Functions::subspaceBasis(dofMapWh, 0),
@@ -1094,7 +1274,7 @@ auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<SolutionRange
// refinementLevels(2));
VTKWriter<GridView> vtkWriter(gridView_CE);
vtkWriter.addVertexData(
- solutionFunction,
+ correctorFunction_1,
VTK::FieldInfo("solution", VTK::FieldInfo::Type::vector, dim));
// vtkWriter.addVertexData(
// solutionFunction,