From dbc79ba992373cb58e62e65ecd184bc7348d7476 Mon Sep 17 00:00:00 2001
From: Klaus Boehnlein <klaus.boehnlein@tu-dresden.de>
Date: Wed, 16 Jun 2021 00:48:40 +0200
Subject: [PATCH] corrected some errors
---
src/dune-microstructure.cc | 359 +++++++++++++++++++++++++++++++++----
1 file changed, 320 insertions(+), 39 deletions(-)
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 187ae5a4..1c8835b5 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -50,6 +50,9 @@
#include <dune/microstructure/prestrainimp.hh>
+#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
+
+
using namespace Dune;
@@ -392,6 +395,64 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb)
}
}
+
+ //////////////////////////////////////////////////////////////////
+ // setIntegralZero:
+ int arbitraryIndex = 0;
+
+ FieldVector<int,3> row;
+ unsigned int elementCounter = 1;
+
+ for (const auto& element : elements(basis.gridView()))
+ {
+ localView.bind(element);
+
+ if(elementCounter == 1) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
+ {
+ for (int k = 0; k < 3; k++)
+ {
+ auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
+ row[k] = localView.index(rowLocal);
+ }
+ }
+
+ // "global assembly"
+ for (int k = 0; k<3; k++)
+ for (size_t j=0; j<localView.size(); j++)
+ {
+ auto col = localView.index(j);
+ nb.add(row[k],col);
+ }
+ elementCounter++;
+ }
+ /////////////////////////////////////////////////////////////////
+
+ ////////////////////////////////////////////////////////////////
+ // setOneBaseFunctionToZero necessary?
+
+ //Determine 3 global indices (one for each component of an arbitrary FE-function).. oder mittels Offset?
+ // use first element:
+// auto it = basis.gridView().template begin<0>();
+// localView.bind(*it);
+//
+// int arbitraryIndex = 0;
+// FieldVector<int,3> row; //fill with Indices..
+//
+//
+// for (int k = 0; k < 3; k++)
+// {
+// auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
+// row[k] = localView.index(rowLocal);
+// auto cIt = stiffnessMatrix[row[k]].begin();
+// auto cEndIt = stiffnessMatrix[row[k]].end();
+// for (; cIt!=cEndIt; ++cIt)
+// *cIt = (cIt.index()==row[k]) ? 1.0 : 0.0;
+// }
+
+
+
+
+
}
@@ -585,7 +646,7 @@ void assembleCellStiffness(const Basis& basis,
const int localPhiOffset = localView.size();
- std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
+// std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
Dune::Matrix<double> elementMatrix; //eher das ??
// Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
@@ -653,7 +714,7 @@ void assembleCellLoad(const Basis& basis,
std::cout << "assemble load vector." << std::endl;
- b.resize(basis.size());
+ b.resize(basis.size()+3);
b = 0;
auto localView = basis.localView();
@@ -709,6 +770,34 @@ void setOneBaseFunctionToZero(const Basis& basis,
)
{
+ constexpr int dim = 3;
+ auto localView = basis.localView();
+
+
+
+ //Determine 3 global indices (one for each component of an arbitrary FE-function).. oder mittels Offset?
+
+ // use first element:
+ auto it = basis.gridView().template begin<0>();
+ localView.bind(*it);
+
+ int arbitraryIndex = 0;
+ FieldVector<int,3> row; //fill with Indices..
+
+ 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;
+ load_alpha1[row[k]] = 0.0;
+ load_alpha2[row[k]] = 0.0;
+ load_alpha3[row[k]] = 0.0;
+ auto cIt = stiffnessMatrix[row[k]].begin();
+ auto cEndIt = stiffnessMatrix[row[k]].end();
+ for (; cIt!=cEndIt; ++cIt)
+ *cIt = (cIt.index()==row[k]) ? 1.0 : 0.0;
+ }
+
}
@@ -762,11 +851,11 @@ void setIntegralZero (const Basis& basis,
{
// if (elementCounter % 50 == 1)
// std::cout << "integral = zero - treatment - element: " << elementCounter << std::endl;
-
+// std::cout << "element-number: " << elementCounter << std::endl;
localView.bind(element);
- if(elementCounter == 1) // get arbitraryIndex(global) for each Component
+ if(elementCounter == 1) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
{
for (int k = 0; k < dim; k++)
{
@@ -805,29 +894,28 @@ void setIntegralZero (const Basis& basis,
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]);
+// std::vector<FieldVector<double,1> > BasisFunctionValues(shapeFunctionValues.size());
+// for (size_t i=0; i<BasisFunctionValues.size(); i++)
+// jacobian.mv(shapeFunctionValues[i], BasisFunctionValues[i]); // no gradient! not needed!
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...
+ elementColumns[idx] += shapeFunctionValues[i] * quadPoint.weight() * integrationElement; // bei Periodicbasis einfach mit child(k) arbeiten...
+// 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;
+ stiffnessMatrix[row[k]][col] += elementColumns[j] * quadPoint.weight() * integrationElement; // TODO Müsste über eine LOOP ?
}
- elementCounter ++;
+ elementCounter++;
}//end of quad
}//end for each element
@@ -838,6 +926,85 @@ void setIntegralZero (const Basis& basis,
+template<class Basis, class LocalFunction>
+double integralMean(const Basis& basis,
+ LocalFunction& fun
+ )
+{
+
+ using GridView = typename Basis::GridView;
+ using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+
+
+ constexpr int dim = 3;
+ using FuncScalar = std::function< FieldVector< double, dim>(const Domain&) >;
+
+ auto localView = basis.localView();
+ double r = 0.0;
+ double area = 0.0;
+
+ // Compute area integral
+ for(const auto& element : elements(basis.gridView()))
+ {
+
+ 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);
+
+ for (const auto& quadPoint : quad)
+ {
+ const double integrationElement = element.geometry().integrationElement(quadPoint.position());
+ area += quadPoint.weight() * integrationElement;
+ }
+ }
+ std::cout << "Domain-Area: " << area << std::endl;
+
+
+ for(const auto& element : elements(basis.gridView()))
+ {
+
+ localView.bind(element);
+ fun.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);
+
+ for (const auto& quadPoint : quad)
+ {
+
+ const auto& quadPos = quadPoint.position();
+ const double integrationElement = element.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 value = fun(quadPos);
+// 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;
+
+
+// r += tmp * quadPoint.weight() * integrationElement;
+
+ }
+ }
+
+
+ return r/area;
+
+}
+
+
+
+
// ---- TODO - MatrixEmbedding function (iota) ---- ???
@@ -944,7 +1111,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 = 3;
std::array<int,dim> nElements={nE,nE,nE}; //#Elements in each direction
@@ -954,7 +1121,7 @@ int main(int argc, char *argv[])
const GridView gridView_CE = grid_CE.leafGridView();
-
+ double areaDomain = 1.0;
@@ -1026,11 +1193,13 @@ int main(int argc, char *argv[])
std::cout << "size feBasis: " << Basis_CE.size() << std::endl;
std::cout << "basis.dimension()" << Basis_CE.dimension() <<std::endl;
+ const int phiOffset = Basis_CE.size();
+
/////////////////////////////////////////////////////////
// Stiffness matrix and right hand side vector
/////////////////////////////////////////////////////////
-// { linear_algebra_setup_begin }
+
// using VelocityVector = BlockVector<FieldVector<double,dim>>;
// using PressureVector = BlockVector<double>;
// using Vector = MultiTypeBlockVector<VelocityVector, PressureVector>;
@@ -1041,7 +1210,7 @@ int main(int argc, char *argv[])
// using MatrixRow0 = MultiTypeBlockVector<Matrix00, Matrix01>;
// using MatrixRow1 = MultiTypeBlockVector<Matrix10, Matrix11>;
// using Matrix = MultiTypeBlockMatrix<MatrixRow0,MatrixRow1>;
-// { linear_algebra_setup_end }
+
/////////////////////////////////////////////////////////
@@ -1079,18 +1248,18 @@ Func2Tensor x3G_3 = [] (const Domain& z) {
- 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};
-
- printmatrix(std::cout, basisContainer[0] , "G_1", "--");
- printmatrix(std::cout, basisContainer[1] , "G_2", "--");
- printmatrix(std::cout, basisContainer[2] , "´G_3", "--");
-
+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};
+// printmatrix(std::cout, basisContainer[0] , "G_1", "--");
+// printmatrix(std::cout, basisContainer[1] , "G_2", "--");
+// printmatrix(std::cout, basisContainer[2] , "´G_3", "--");
+//
+//
+//
assembleCellStiffness(Basis_CE, muLocal, lambdaLocal, gamma, stiffnessMatrix_CE);
@@ -1103,19 +1272,19 @@ 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", "--");
-}
-
+// 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
+// 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 setIntegralZero", "--");
+ printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
}
@@ -1123,6 +1292,12 @@ if(set_oneBasisFunction_Zero)
+
+
+
+
+
+
// setzt gesamte Zeile = 0:
// stiffnessMatrix_CE[2] = 0;
// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix MOD", "--");
@@ -1242,6 +1417,42 @@ solver.apply(x_3, load_alpha3, statistics);
+// Extract phi_alpha & M_alpha coefficients
+VectorCT phi_1, phi_2, phi_3;
+phi_1.resize(Basis_CE.size());
+phi_1 = 0;
+phi_2.resize(Basis_CE.size());
+phi_2 = 0;
+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;
+
+
+
+for(size_t i=0; i<Basis_CE.size();i++)
+{
+ phi_1[i] = x_1[i];
+ phi_2[i] = x_2[i];
+ phi_3[i] = x_3[i];
+}
+
+for(size_t i=0; i<3; i++)
+{
+ m_1[i] = x_1[phiOffset+i];
+ m_2[i] = x_2[phiOffset+i];
+ m_3[i] = x_3[phiOffset+i];
+}
+
+
////////////////////////////////////////////////////////////////////////////
// Make a discrete function from the FE basis and the coefficient vector
////////////////////////////////////////////////////////////////////////////
@@ -1249,18 +1460,88 @@ solver.apply(x_3, load_alpha3, statistics);
using SolutionRange = FieldVector<double,dim>;
// using SolutionRange = double;
+
+// auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+// Basis_CE,
+// Dune::Functions::HierarchicVectorWrapper<VectorCT, double>(phi_1)); // ??
+
+
+
auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
Basis_CE,
- x_1);
+ phi_1);
+
+
+auto correctorFunction_2 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ Basis_CE,
+ phi_2);
+
+auto correctorFunction_3 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ Basis_CE,
+ phi_3);
+
+
+// assemble M_alpha's (actually iota(M_alpha) )
+
+MatrixRT M_1(0), M_2(0), M_3(0);
+
+for(size_t i=0; i<3; i++)
+{
+ M_1 += m_1[i]*basisContainer[i];
+ M_2 += m_2[i]*basisContainer[i];
+ M_3 += m_3[i]*basisContainer[i];
+}
+
+
+printmatrix(std::cout, M_1, "Corrector-Matrix M_1", "--");
+printmatrix(std::cout, M_2, "Corrector-Matrix M_1", "--");
+printmatrix(std::cout, M_3, "Corrector-Matrix M_1", "--");
+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 ???
+}
+
+
+// INTEGRAL-MEAN TEST:
+auto AVFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ Basis_CE,
+ 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;
+
+
+
+
+// FieldVector<double,3> ones = {1.0, 1.0, 1.0};
+//
+// phi_1 - IM*;
+
+
+// // const-function
+// auto f1 = [](const auto& x)
+// {
+// return IM; // does not work..
+// };
+//
+//
+// VectorCT IM_1;
+// interpolate(Basis_CE, IM_1, f1);
- // use only first child of W_h to plot solution
-// auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
-// Functions::subspaceBasis(dofMapWh, 0),
-// x);
//////////////////////////////////////////////////////////////////////////////////////////////
--
GitLab