diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index e8cbf42b2d07bbdb2d388c80c87eb75d035aed68..dde5d76f8718c72696906f8461f0165fbcaf5c2e 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -188,7 +188,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
energyDensity += stressU[ii][jj] * strainV[ii][jj]; // "phi*phi"-part
// energyDensity += stressU[ii][jj] * strainU[ii][jj];
-// size_t row = localView.tree().child(k).localIndex(i); // kann auf Unterscheidung zwischen k & l verzichten?!
+// size_t row = localView.tree().child(k).localIndex(i); // kann auf Unterscheidung zwischen k & l verzichten?
// size_t col = localView.tree().child(l).localIndex(j); // siehe DUNE-book p.394
// size_t col = localView.tree().child(k).localIndex(j); // Indizes mit k=l genügen .. Kroenecker-Delta_kl NEIN???
size_t col = localView.tree().child(k).localIndex(i); // kann auf Unterscheidung zwischen k & l verzichten?!
@@ -249,7 +249,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
elementMatrix[row][localPhiOffset+m] += value;
- elementMatrix[localPhiOffset+m][row] += value; // ---- reicht das ? --- TODO
+ elementMatrix[localPhiOffset+m][row] += value;
}
@@ -302,9 +302,6 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// elementMatrix[localPhiOffset+n][col] += energyDensityPhiG * quadPoint.weight() * integrationElement;
//
// }
-//
-//
-//
// }
@@ -357,7 +354,7 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb)
const int phiOffset = basis.dimension();
- nb.resize(basis.size()+3, basis.size()+3); // ???
+ nb.resize(basis.size()+3, basis.size()+3); //
for (const auto& element : elements(gridView))
{
@@ -465,11 +462,8 @@ void computeElementLoadVector( const LocalView& localView,
LocalFunction1& mu,
LocalFunction2& lambda,
const double gamma,
-// BlockVector<double>& elementRhs,
Vector& elementRhs,
const Force& forceTerm
-// const std::function<double(FieldVector<double,
-// LocalView::Element::dimension>)> volumeTerm
)
{
// | f*phi|
@@ -477,7 +471,6 @@ void computeElementLoadVector( const LocalView& localView,
// | f*m |
using Element = typename LocalView::Element;
-// auto element = localView.element();
const auto element = localView.element();
const auto geometry = element.geometry();
constexpr int dim = Element::dimension;
@@ -500,7 +493,7 @@ void computeElementLoadVector( const LocalView& localView,
const int localPhiOffset = localView.size();
///////////////////////////////////////////////
- // Basis for R_sym^{2x2} // wird nicht als Funktion benötigt da konstant...
+ // 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}};
@@ -509,13 +502,12 @@ void computeElementLoadVector( const LocalView& localView,
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1); // für simplex
-// const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
for (const auto& quadPoint : quad){
-// const Domain& quadPos = quadPoint.position();
+
const FieldVector<double,dim>& quadPos = quadPoint.position();
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
const double integrationElement = geometry.integrationElement(quadPos);
@@ -617,12 +609,6 @@ void assembleCellStiffness(const Basis& basis,
LocalFunction2& lambdaLocal,
const double gamma,
Matrix& matrix
-// Vector& b,
-// const Force& forceTerm
-// const std::function<
-// double(FieldVector<double,
-// Basis::GridView::dimension>)
-// > volumeTerm
)
{
std::cout << "assemble stiffnessmatrix begins." << std::endl;
@@ -648,33 +634,29 @@ void assembleCellStiffness(const Basis& basis,
const int localPhiOffset = localView.size();
// std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
- Dune::Matrix<double> elementMatrix; //eher das ??
-// Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
+// Dune::Matrix<double> elementMatrix;
+ Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
-
-
// std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
// std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
-// Dune::Matrix<double> TestelementMatrix;
-// computeElementStiffnessMatrixOLD(localView, TestelementMatrix, muLocal, lambdaLocal, gamma);
-// printmatrix(std::cout, TestelementMatrix, "TESTElementMatrix", "--");
+
//////////////////////////////////////////////////////////////////////////////
- // GLOBAL STIFFNES ASSEMBLY ..(hier dreht sich i, j..)
+ // GLOBAL STIFFNES ASSEMBLY
//////////////////////////////////////////////////////////////////////////////
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);
for (size_t j=0; j<localPhiOffset; j++ )
{
- // The global index of the j-th local degree of freedom of the current element
+
auto col = localView.index(j);
- matrix[row][col] += elementMatrix[i][j]; // keine periodicMap nötig..wird durch Basis behandelt..
+ matrix[row][col] += elementMatrix[i][j];
}
}
@@ -766,7 +748,8 @@ void setOneBaseFunctionToZero(const Basis& basis,
Matrix& stiffnessMatrix,
Vector& load_alpha1,
Vector& load_alpha2,
- Vector& load_alpha3
+ Vector& load_alpha3,
+ FieldVector<int,3>& indexVector // oder input : arbitraryLocalIndex und dann hier function computeArbitraryIndex aufrufen
)
{
@@ -789,7 +772,7 @@ void setOneBaseFunctionToZero(const Basis& basis,
auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
row[k] = localView.index(rowLocal);
std::cout << "row[k]:" << row[k] << std::endl;
- load_alpha1[row[k]] = 0.0;
+ load_alpha1[row[k]] = 0.0; //verwende hier indexVector
load_alpha2[row[k]] = 0.0;
load_alpha3[row[k]] = 0.0;
auto cIt = stiffnessMatrix[row[k]].begin();
@@ -925,12 +908,76 @@ void setIntegralZero (const Basis& basis,
+template<class Basis>
+auto childToIndexMap(const Basis& basis,
+ const int k
+ )
+{
+ constexpr int dim = 3;
+
+
+
+ using GridView = typename Basis::GridView;
+ using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+
+
+ using FuncScalar = std::function< FieldVector< double, dim>(const Domain&) >;
+
+ auto localView = basis.localView();
+
+
+
+ std::vector<int> r = { };
+// for (int n: r)
+// std::cout << n << ","<< std::endl;
+
+ // Determine global indizes for each component k = 1,2,3.. in order to subtract correct (component of) integral Mean
+ // (global) Indices that correspond to component k = 1,2,3
+ for(const auto& element : elements(basis.gridView()))
+ {
+
+ localView.bind(element);
+ const auto& localFiniteElement = localView.tree().child(k).finiteElement(); // macht keinen Unterschied ob hier k oder 0..
+ const auto nSf = localFiniteElement.localBasis().size();
+
+
+ for(size_t j=0; j<nSf; j++)
+ {
+ auto Localidx = localView.tree().child(k).localIndex(j);
+ r.push_back(localView.index(Localidx));
+ }
+
+
+
+ }
+
+ // Delete duplicate elements
+ // first remove consecutive (adjacent) duplicates
+ auto last = std::unique(r.begin(), r.end());
+ r.erase(last, r.end());
+ // sort followed by unique, to remove all duplicates
+ std::sort(r.begin(), r.end());
+ last = std::unique(r.begin(), r.end());
+ r.erase(last, r.end());
+
+
+ return r;
+
+
+}
+
+
+
+
+
+
template<class Basis, class LocalFunction>
-double integralMean(const Basis& basis,
+auto integralMean(const Basis& basis,
LocalFunction& fun
)
{
+ // computes integral mean of given LocalFunction
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
@@ -940,7 +987,8 @@ double integralMean(const Basis& basis,
using FuncScalar = std::function< FieldVector< double, dim>(const Domain&) >;
auto localView = basis.localView();
- double r = 0.0;
+// double r = 0.0;
+ FieldVector<double,3> r = {0.0, 0.0, 0.0};
double area = 0.0;
// Compute area integral
@@ -988,8 +1036,9 @@ double integralMean(const Basis& basis,
// auto value = localPhi_1(quadPos);
// auto value = FieldVector<double,dim>{1.0, 1.0, 1.0};
- for(size_t k=0; k<3; k++)
- r += value[k] * quadPoint.weight() * integrationElement;
+// for(size_t k=0; k<3; k++)
+// r += value[k] * quadPoint.weight() * integrationElement;
+ r += value * quadPoint.weight() * integrationElement;
// r += tmp * quadPoint.weight() * integrationElement;
@@ -997,8 +1046,8 @@ double integralMean(const Basis& basis,
}
}
-
- return r/area;
+// return r/area;
+ return (1/area) * r ;
}
@@ -1007,30 +1056,37 @@ double integralMean(const Basis& basis,
-// ---- TODO - MatrixEmbedding function (iota) ---- ???
+template<class Basis, class LocalFunction, class Vector>
+auto subtractIntegralMean(const Basis& basis,
+ LocalFunction& fun,
+ Vector& coeffVector
+ )
+{
+ // subtract correct Integral mean from each associated component function
+ auto IM = integralMean(basis, fun);
+ constexpr int dim = 3;
+ for(size_t k=0; k<dim; k++)
+ {
+ std::cout << "Integral-Mean: " << IM[k] << std::endl;
+ auto idx = childToIndexMap(basis,k);
+
+ for ( int i : idx)
+ coeffVector[i] -= IM[k];
+ }
+}
- // ---------------------------------------------------------------------------------------------
- // -- OPTIONS: use power(periodicbasis) and equivalent only for component 1 & 2 ???
- // use composite-Basis with 2D-Periodic Basis and 1D-Lagrange
+// ---- TODO - MatrixEmbedding function (iota) ---- ?
+
- // Check whether two points are equal on R/Z x R/Z
-// auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>& y)
-// {
-// 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])) );
-// // return ((FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1])));
-// // and (x[0]<1e-8 or x[0] > 0.999));
-// // return (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]));
-// };
@@ -1052,14 +1108,12 @@ int main(int argc, char *argv[])
MPIHelper::instance(argc, argv);
- // input getter
+
ParameterTree parameterSet;
// if (argc < 2)
// ParameterTreeParser::readINITree("../src/cellsolver.parset", parameterSet);
-
// ParameterTreeParser::readINITree(argv[1], parameterSet);
// ParameterTreeParser::readOptions(argc, argv, parameterSet);
-
// output setter
// -- not set up --
@@ -1083,7 +1137,7 @@ int main(int argc, char *argv[])
- // Plate geometry settings
+// // Plate geometry settings
// double width = parameterSet.get<double>("width", 1.0); //geometry cell, cross section
// double len = parameterSet.get<double>("len", 10.0); //length
// double height = parameterSet.get<double>("h", 0.1); //rod thickness
@@ -1111,7 +1165,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 = 1;
+ int nE = 2;
std::array<int,dim> nElements={nE,nE,nE}; //#Elements in each direction
@@ -1121,7 +1175,7 @@ int main(int argc, char *argv[])
const GridView gridView_CE = grid_CE.leafGridView();
- double areaDomain = 1.0;
+
@@ -1129,25 +1183,17 @@ int main(int argc, char *argv[])
using ScalarRT = FieldVector< double, 1>;
using VectorRT = FieldVector< double, nCompo>;
using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
-
-// using Domain = typename gridView_CE::template Codim<0>::Geometry::GlobalCoordinate;
-// using Domain = CellGridType::LeafGridView::Codim<0>::Geometry::GlobalCoordinate;
using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
-
using FuncScalar = std::function< ScalarRT(const Domain&) >;
using Func2Tensor = std::function< MatrixRT(const Domain&) >;
- using VectorCT = BlockVector<FieldVector<double,1> >; // Rob
- // using VectorCT = BlockVector<FieldVector<double,dim> >; // oder das ?
- using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >; // Rob
+ using VectorCT = BlockVector<FieldVector<double,1> >;
+ using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >;
+ // using VectorCT = BlockVector<FieldVector<double,dim> >;
// using MatrixCT = BCRSMatrix<FieldMatrix<double,dim,dim> >;
using ElementMatrixCT = Matrix<FieldMatrix<double,1,1> >;
- // const FuncScalar mu_;
- // const FuncScalar lambda_;
- // double gamma_;
-
double mu1 = 0.5 * 1.0/(1.0 + nu1) * E1;
@@ -1185,9 +1231,10 @@ int main(int argc, char *argv[])
auto Basis_CE = makeBasis(
gridView_CE,
power<dim>(
+// lagrange<1>(),
periodic(lagrange<1>(), periodicIndices),
-// lagrange<order_CE>(),
- flatLexicographic()));//flatInterleaved()
+ flatLexicographic()));//
+// flatInterleaved()));
std::cout << "size feBasis: " << Basis_CE.size() << std::endl;
@@ -1195,29 +1242,24 @@ int main(int argc, char *argv[])
const int phiOffset = Basis_CE.size();
-/////////////////////////////////////////////////////////
-// Stiffness matrix and right hand side vector
-/////////////////////////////////////////////////////////
-
+ // TEST
+ auto X1Basis = subspaceBasis(Basis_CE, 0);
+ auto X2Basis = subspaceBasis(Basis_CE, 1);
+ auto X3Basis = subspaceBasis(Basis_CE, 2);
-// using VelocityVector = BlockVector<FieldVector<double,dim>>;
-// using PressureVector = BlockVector<double>;
-// using Vector = MultiTypeBlockVector<VelocityVector, PressureVector>;
-// using Matrix00 = BCRSMatrix<FieldMatrix<double,dim,dim>>;
-// using Matrix01 = BCRSMatrix<FieldMatrix<double,dim,1>>;
-// using Matrix10 = BCRSMatrix<FieldMatrix<double,1,dim>>;
-// using Matrix11 = BCRSMatrix<double>;
-// using MatrixRow0 = MultiTypeBlockVector<Matrix00, Matrix01>;
-// using MatrixRow1 = MultiTypeBlockVector<Matrix10, Matrix11>;
-// using Matrix = MultiTypeBlockMatrix<MatrixRow0,MatrixRow1>;
+ std::cout << "size X1Basis: " << X1Basis.size() << std::endl;
+ std::cout << "X1Basis.dimension()" << X1Basis.dimension() <<std::endl;
+ std::cout << "size X2Basis: " << X2Basis.size() << std::endl;
+ std::cout << "size X3Basis: " << X3Basis.size() << std::endl;
/////////////////////////////////////////////////////////
-// Assemble the system
+// Stiffness matrix and right hand side vector
/////////////////////////////////////////////////////////
+
VectorCT load_alpha1, load_alpha2, load_alpha3;
// auto rhsBackend = Functions::istlVectorBackend(b);
@@ -1262,13 +1304,19 @@ 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", "--");
-//
-//
-//
+
+
+
+
+
+/////////////////////////////////////////////////////////
+// Assemble the system
+/////////////////////////////////////////////////////////
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);
@@ -1285,27 +1333,39 @@ printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix before B.C", "--");
//
-// set one basis-function to zero
-if(set_oneBasisFunction_Zero)
-{
- setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3);
- printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
-}
+// Determine global indixes of components for arbitrary (local) index TODO :separate Function!
+
+auto localView = Basis_CE.localView();
+
+int arbitraryIndex = 0;
+FieldVector<int,3> row; //fill with Indices..
+ // use first element:
+ auto it = Basis_CE.gridView().template begin<0>();
+ localView.bind(*it);
+ 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;
+ }
+// 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", "--");
+}
-// setzt gesamte Zeile = 0:
-// stiffnessMatrix_CE[2] = 0;
-// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix MOD", "--");
@@ -1385,10 +1445,6 @@ if(set_oneBasisFunction_Zero)
////////////////////////////
// Compute solution
////////////////////////////
-
-// { stokes_solve_begin }
-// Initial iterate: Start from the rhs vector,
-// that way the Dirichlet entries are already correct.
VectorCT x_1 = load_alpha1;
VectorCT x_2 = load_alpha2;
VectorCT x_3 = load_alpha3;
@@ -1414,7 +1470,6 @@ RestartedGMResSolver<VectorCT> solver(
InverseOperatorResult statistics;
-
// solve for different Correctors (alpha = 1,2,3)
solver.apply(x_1, load_alpha1, statistics);
solver.apply(x_2, load_alpha2, statistics);
@@ -1468,7 +1523,7 @@ using SolutionRange = FieldVector<double,dim>;
// auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
// Basis_CE,
-// Dune::Functions::HierarchicVectorWrapper<VectorCT, double>(phi_1)); // ??
+// Dune::Functions::HierarchicVectorWrapper<VectorCT, double>(phi_1));
@@ -1509,14 +1564,11 @@ auto localPhi_1 = localFunction(correctorFunction_1);
-double IM = 0.0;
-IM = integralMean(Basis_CE, localPhi_1);
-std::cout << "Integral-Mean: " << IM << std::endl;
-for(size_t i=0; i<Basis_CE.size();i++)
-{
- phi_1[i] -= IM; // Müsste integralMean von jeder Komponente einzeln berechnen ???
-}
+
+
+subtractIntegralMean(Basis_CE, localPhi_1 , phi_1);
+
// INTEGRAL-MEAN TEST:
@@ -1525,10 +1577,11 @@ auto AVFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
phi_1);
auto AVPhi_1 = localFunction(AVFunction_1);
-double A = 0.0;
-A = integralMean(Basis_CE, AVPhi_1);
-std::cout << "Integral-Mean (TEST) : " << A << std::endl;
-
+auto A = integralMean(Basis_CE, AVPhi_1);
+for(size_t i=0; i<3; i++)
+{
+std::cout << "Integral-Mean (TEST) : " << A[i] << std::endl;
+}
@@ -1549,6 +1602,9 @@ std::cout << "Integral-Mean (TEST) : " << A << std::endl;
+
+
+
//////////////////////////////////////////////////////////////////////////////////////////////
// Write result to VTK file
// We need to subsample, because the dune-grid VTKWriter cannot natively display
@@ -1569,7 +1625,4 @@ std::cout << "Integral-Mean (TEST) : " << A << std::endl;
std::cout << "wrote data to file: CellProblem-result" << std::endl; // better with string for output name..
-
-
-
}