diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 38f9d535e7ca4d8a1f7ac4fac431664f5a06a085..0359ad4c7d753fe4ebfa9b314d6f98a4fda491d1 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -551,7 +551,7 @@ void computeElementLoadVector( const LocalView& localView,
double energyDensity = 0.0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
- energyDensity += stressV[ii][jj] *forceTerm(quadPos)[ii][jj]; //TEST VERWENDE NEGATIVE DER LAST
+ energyDensity += stressV[ii][jj] *forceTerm(quadPos)[ii][jj];
size_t row = localView.tree().child(k).localIndex(i);
@@ -577,7 +577,7 @@ void computeElementLoadVector( const LocalView& localView,
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]; //TEST VERWENDE NEGATIVE DER LAST
+ energyDensityfG += stressG[ii][jj] * forceTerm(quadPos)[ii][jj];
elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
}
@@ -1078,126 +1078,6 @@ void setOneBaseFunctionToZero(const Basis& basis,
-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;
- // std::cout << "element-number: " << elementCounter << std::endl;
- localView.bind(element);
-
-
- if(elementCounter == 1) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
- {
- 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]); // 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] += 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; // TODO Müsste über eine LOOP ?
- }
-
- elementCounter++;
- } //end of quad
-
- } //end for each element
-}
-
-
-
template<class Basis>
@@ -1264,19 +1144,27 @@ auto childToIndexMap(const Basis& basis,
-template<class Basis, class LocalFunction>
+template<class Basis, class Vector>
auto integralMean(const Basis& basis,
- LocalFunction& fun
+ Vector& coeffVector
)
{
// computes integral mean of given LocalFunction
+
+
+
+
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 GVFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis,coeffVector);
+ auto localfun = localFunction(GVFunction);
+
+
auto localView = basis.localView();
// double r = 0.0;
@@ -1307,7 +1195,7 @@ auto integralMean(const Basis& basis,
{
localView.bind(element);
- fun.bind(element);
+ localfun.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);
@@ -1326,7 +1214,7 @@ auto integralMean(const Basis& basis,
// for (size_t i=0; i< gradients.size(); i++)
// jacobian.mv(referenceGradients[i][0], gradients[i]);
- auto value = fun(quadPos);
+ auto value = localfun(quadPos);
// auto value = localPhi_1(quadPos);
// auto value = FieldVector<double,dim>{1.0, 1.0, 1.0};
@@ -1341,7 +1229,7 @@ auto integralMean(const Basis& basis,
}
// return r/area;
- return (1/area) * r ;
+ return (1.0/area) * r ;
}
@@ -1351,15 +1239,16 @@ auto integralMean(const Basis& basis,
-template<class Basis, class LocalFunction, class Vector>
+
+template<class Basis, 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);
+// auto IM = integralMean(basis, fun);
+ auto IM = integralMean(basis, coeffVector);
constexpr int dim = 3;
for(size_t k=0; k<dim; k++)
@@ -1370,7 +1259,6 @@ auto subtractIntegralMean(const Basis& basis,
for ( int i : idx)
coeffVector[i] -= IM[k];
}
-
}
@@ -1698,7 +1586,7 @@ int main(int argc, char *argv[])
FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0, -1.0/2.0});
FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0, 1.0/2.0});
- int nE = 20;
+ int nE = 10;
log << "Number of Elements: " << nE << std::endl;
debugLog << "Number of Elements: " << nE << std::endl;
@@ -1817,7 +1705,7 @@ int main(int argc, char *argv[])
Functions::Experimental::BasisFactory::periodic(lagrange<1>(), periodicIndices),
flatLexicographic())); //
// blockedInterleaved())); // siehe Periodicbasistest.. funktioniert auch?
- // flatInterleaved()));
+// flatInterleaved()));
std::cout << "size feBasis: " << Basis_CE.size() << std::endl;
@@ -1853,15 +1741,15 @@ int main(int argc, char *argv[])
+ bool set_IntegralZero = parameterSet.get<bool>("set_IntegralZero", true);
+ bool set_oneBasisFunction_Zero = false;
+ bool substract_integralMean = false;
+ if(set_IntegralZero)
+ {
+ set_oneBasisFunction_Zero = true;
+ substract_integralMean = true;
+ }
-
-
- // bool set_integral_zero = true;
-
- bool set_oneBasisFunction_Zero = true;
- 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;
@@ -1882,6 +1770,19 @@ int main(int argc, char *argv[])
return MatrixRT{{0.0, 1.0/sqrt(2)*x[2], 0.0}, {1.0/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
+
+
+ Func2Tensor x3G_1neg = [x3G_1] (const Domain& x) {
+ return -1.0*x3G_1(x);
+ };
+
+ Func2Tensor x3G_2neg = [x3G_2] (const Domain& x) {
+ return -1.0*x3G_2(x);;
+ };
+
+ 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", "--");
@@ -1928,36 +1829,24 @@ int main(int argc, char *argv[])
// 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);
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha1 ,x3G_1neg);
+ 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];
- }
+// 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", "--");
// printvector(std::cout, load_alpha1, "load_alpha1", "--");
-
- // 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", "--");
- // }
- //
-
/////////////////////////////////////////////////////////
@@ -2010,8 +1899,6 @@ int main(int argc, char *argv[])
if (Solvertype == 1) // CG - SOLVER
{
-
-
std::cout << "------------ CG - Solver ------------" << std::endl;
MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_CE);
@@ -2101,59 +1988,21 @@ int main(int argc, char *argv[])
phi_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];
}
-
+
+ FieldVector<double,3> m_1, m_2, m_3;
+
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
- //////////////////////////////////////////////////////////////////////////// // ERROR : HIER NOCH KEIN Integralmittelabgezogen!?!?!?!?!
- 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,
- 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) )
@@ -2172,41 +2021,37 @@ int main(int argc, char *argv[])
printmatrix(std::cout, M_3, "Corrector-Matrix M_3", "--");
+
+
+
- auto localPhi_1 = localFunction(correctorFunction_1);
- auto localPhi_2 = localFunction(correctorFunction_2);
- auto localPhi_3 = localFunction(correctorFunction_3);
+ ////////////////////////////////////////////////////////////////////////////
+ // Substract Integral-mean
+ //////////////////////////////////////////////////////////////////////////// // ERROR :
+ using SolutionRange = FieldVector<double,dim>;
std::cout << "check integralMean phi_1: " << std::endl;
- auto A = integralMean(Basis_CE,localPhi_1);
+ 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)
{
std::cout << " --- substracting integralMean --- " << std::endl;
- subtractIntegralMean(Basis_CE, localPhi_1 , phi_1);
- subtractIntegralMean(Basis_CE, localPhi_2 , phi_2);
- subtractIntegralMean(Basis_CE, localPhi_3 , phi_3);
+ subtractIntegralMean(Basis_CE, phi_1);
+ subtractIntegralMean(Basis_CE, phi_2);
+ subtractIntegralMean(Basis_CE, phi_3);
- // AUCH VON X-Vectoren!!!
// printvector(std::cout, x_1 , "x_1", "--" );
- subtractIntegralMean(Basis_CE, localPhi_1 , x_1); //TEST
- subtractIntegralMean(Basis_CE, localPhi_2 , x_2);
- subtractIntegralMean(Basis_CE, localPhi_3 , x_3);
-
+ 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 AVFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
- Basis_CE,
- phi_1);
-
- auto AVPhi_1 = localFunction(AVFunction_1);
- auto A = integralMean(Basis_CE, AVPhi_1);
+ auto A = integralMean(Basis_CE, phi_1);
for(size_t i=0; i<3; i++)
{
std::cout << "Integral-Mean (CHECK) : " << A[i] << std::endl;
@@ -2214,13 +2059,13 @@ int main(int argc, char *argv[])
////////////////////////////////////////////////////
}
-
- // // PRINT Corrector-Basis-Coefficients
-
+// // 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 in Logs
/////////////////////////////////////////////////////////
@@ -2243,19 +2088,19 @@ int main(int argc, char *argv[])
log << x_3 << std::endl;
}
- ////////////////////////////////////////////////////////////////////////////////
- // REMAKE Corrector-Functions after substract_integralMean
- ////////////////////////////////////////////////////////////////////////////////
- auto correctorFunction_1New = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ ////////////////////////////////////////////////////////////////////////////
+ // Make a discrete function from the FE basis and the coefficient vector
+ //////////////////////////////////////////////////////////////////////////// // ERROR :
+ using SolutionRange = FieldVector<double,dim>;
+ auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
Basis_CE,
phi_1);
-
- auto correctorFunction_2New = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ auto correctorFunction_2 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
Basis_CE,
phi_2);
- auto correctorFunction_3New = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ auto correctorFunction_3 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
Basis_CE,
phi_3);
@@ -2298,12 +2143,10 @@ int main(int argc, char *argv[])
//
// }
-
// printvector(std::cout, coeffT_1 , "coeffT_1", "--" );
// printvector(std::cout, x_1 , "x_1", "--" );
-
/////////////////////////////////////////////////////////
// Create Containers for Basis-Matrices, Correctors and Coefficients
/////////////////////////////////////////////////////////
@@ -2311,7 +2154,7 @@ int main(int argc, char *argv[])
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<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};
@@ -2331,21 +2174,20 @@ int main(int argc, char *argv[])
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..
-
+ //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) > verändert man hier x3MatrixBasis????? TODO
+ 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 << "Test : " << coeffContainer[a]*loadContainer[b] << std::endl; //same...
+// std::cout << "Test : " << coeffContainer[a]*loadContainer[b] << std::endl; //same...
}
printmatrix(std::cout, Q, "Matrix Q", "--");
@@ -2565,20 +2407,7 @@ int main(int argc, char *argv[])
// correctorFunction_Test,
// VTK::FieldInfo("corrector Test", VTK::FieldInfo::Type::vector, dim));
-
-
-
- vtkWriter.addVertexData( //TEST
- correctorFunction_1New,
- VTK::FieldInfo("corrector phi_1New", VTK::FieldInfo::Type::vector, dim));
- vtkWriter.addVertexData(
- correctorFunction_2New,
- VTK::FieldInfo("corrector phi_2New", VTK::FieldInfo::Type::vector, dim));
- vtkWriter.addVertexData(
- correctorFunction_3New,
- VTK::FieldInfo("corrector phi_3New", VTK::FieldInfo::Type::vector, dim));
-
-
+
vtkWriter.addVertexData(
correctorFunction_1,
VTK::FieldInfo("corrector phi_1", VTK::FieldInfo::Type::vector, dim));