diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index edc90c05b307810a5ffe8ed1a2edd4a11833627c..e0f840452d6494c0e7cd04b9aa1e308ba29b1711 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -54,6 +54,50 @@ using namespace Dune;
+
+template<class Basis>
+auto arbitraryComponentwiseIndices(const Basis& basis,
+ const int elementNumber,
+ const int localIdx
+ )
+{
+ // (Local Approach -- works for non Lagrangian-Basis too)
+ // Input : elementNumber & localIdx
+ // Output : determine all Component-Indices that correspond to that basis-function
+
+// constexpr int dim = 3;
+
+// using GridView = typename Basis::GridView;
+// using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+//
+
+ auto localView = basis.localView();
+
+ FieldVector<int,3> row;
+ int elementCounter = 0;
+
+ for (const auto& element : elements(basis.gridView()))
+ {
+ if(elementCounter == elementNumber) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
+ {
+ localView.bind(element);
+
+ for (int k = 0; k < 3; k++)
+ {
+ auto rowLocal = localView.tree().child(k).localIndex(localIdx);
+ row[k] = localView.index(rowLocal);
+// std::cout << "row[k]:" << row[k] << std::endl;
+ }
+ }
+ elementCounter++;
+ }
+
+ return row;
+}
+
+
+
+
template<class LocalView, class Matrix, class localFunction1, class localFunction2>
void computeElementStiffnessMatrix(const LocalView& localView,
Matrix& elementMatrix,
@@ -339,8 +383,8 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// Get the occupation pattern of the stiffness matrix
-template<class Basis>
-void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb)
+template<class Basis, class ParameterSet>
+void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet& parameterSet)
{
// OccupationPattern:
// | phi*phi m*phi |
@@ -389,62 +433,43 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb)
}
}
-
+ //////////////////////////////////////////////////////////////////"localIndex is larger than #shapeFunctions",
+ // setOneBaseFunctionToZero
//////////////////////////////////////////////////////////////////
- // 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
+ if(parameterSet.template get<bool>("set_IntegralZero", true)){
+ FieldVector<int,3> row;
+ unsigned int elementCounter = 1;
+ unsigned int arbitraryLocalIndex = parameterSet.template get<unsigned int>("arbitraryLocalIndex", 0);
+ unsigned int arbitraryElementNumber = parameterSet.template get<unsigned int>("arbitraryElementNumber", 0);
+
+ assert(arbitraryLocalIndex < localView.size() );
+ assert(arbitraryElementNumber < basis.gridView().size(0));
+
+ //Determine 3 global indices (one for each component of an arbitrary local FE-function)
+ row = arbitraryComponentwiseIndices(basis,arbitraryElementNumber,arbitraryLocalIndex);
+
+
+// for (const auto& element : elements(basis.gridView()))
// {
-// for (int k = 0; k < 3; k++)
-// {
-// auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
-// row[k] = localView.index(rowLocal);
-// }
-// }
+// localView.bind(element);
//
-// // "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? //TODO besser: if(parameterSet_.get<bool>("set_one_base_function_to_0" )){
-
- //Determine 3 global indices (one for each component of an arbitrary FE-function).. oder mittels Offset?
- // use first element:
- 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);
- }
- }
-
+// if(elementCounter == arbitraryElementNumber) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
+// {
+// for (int k = 0; k < 3; k++)
+// {
+// auto rowLocal = localView.tree().child(k).localIndex(arbitraryLocalIndex);
+// row[k] = localView.index(rowLocal);
+// }
+// }
+// // "global assembly"
+// for (int k = 0; k<3; k++)
+// nb.add(row[k],row[k]);
+//
+// elementCounter++;
+// }
// "global assembly"
for (int k = 0; k<3; k++)
nb.add(row[k],row[k]);
-
- elementCounter++;
}
@@ -732,9 +757,6 @@ void computeElementLoadTEST( const LocalView& localView,
}
-
-
-
*/
@@ -756,18 +778,19 @@ void computeElementLoadTEST( const LocalView& localView,
-template<class Basis, class Matrix, class LocalFunction1, class LocalFunction2>
+template<class Basis, class Matrix, class LocalFunction1, class LocalFunction2, class ParameterSet>
void assembleCellStiffness(const Basis& basis,
LocalFunction1& muLocal,
LocalFunction2& lambdaLocal,
const double gamma,
- Matrix& matrix
+ Matrix& matrix,
+ ParameterSet& parameterSet
)
{
std::cout << "assemble stiffnessmatrix begins." << std::endl;
MatrixIndexSet occupationPattern;
- getOccupationPattern(basis, occupationPattern);
+ getOccupationPattern(basis, occupationPattern, parameterSet);
occupationPattern.exportIdx(matrix);
matrix = 0;
@@ -1007,40 +1030,40 @@ auto energy(const Basis& basis,
-
-
-
-
-template<class Basis, class Matrix, class Vector>
+template<class Basis, class Matrix, class Vector, class ParameterSet>
void setOneBaseFunctionToZero(const Basis& basis,
Matrix& stiffnessMatrix,
Vector& load_alpha1,
Vector& load_alpha2,
Vector& load_alpha3,
- FieldVector<int,3>& indexVector // oder input : arbitraryLocalIndex und dann hier function computeArbitraryIndex aufrufen
+ ParameterSet& parameterSet
+// FieldVector<int,3>& indexVector // TODO oder input : arbitraryLocalIndex und dann hier function computeArbitraryIndex aufrufen
)
{
-
- std::cout << "setting one Basis-function to zero" << std::endl;
+ std::cout << "Setting one Basis-function to zero" << std::endl;
constexpr int dim = 3;
- auto localView = basis.localView();
+// auto localView = basis.localView();
+ FieldVector<int,3> row;
+ unsigned int arbitraryLocalIndex = parameterSet.template get<unsigned int>("arbitraryLocalIndex", 0);
+ unsigned int arbitraryElementNumber = parameterSet.template get<unsigned int>("arbitraryElementNumber", 0);
- //Determine 3 global indices (one for each component of an arbitrary FE-function).. oder mittels Offset?
+ //Determine 3 global indices (one for each component of an arbitrary local FE-function)
+ row = arbitraryComponentwiseIndices(basis,arbitraryElementNumber,arbitraryLocalIndex);
// use first element:
- auto it = basis.gridView().template begin<0>();
- localView.bind(*it);
-
- int arbitraryIndex = 0;
- FieldVector<int,3> row; //fill with Indices..
+// 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);
+// 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; //verwende hier indexVector
load_alpha2[row[k]] = 0.0;
@@ -1058,11 +1081,14 @@ void setOneBaseFunctionToZero(const Basis& basis,
+
template<class Basis>
auto childToIndexMap(const Basis& basis,
const int k
)
{
+ // Input : child/component
+ // Output : determine all Indices that correspond to that component
constexpr int dim = 3;
@@ -1112,8 +1138,6 @@ auto childToIndexMap(const Basis& basis,
return r;
-
-
}
@@ -1208,7 +1232,6 @@ auto integralMean(const Basis& basis,
// return r/area;
return (1.0/area) * r ;
-
}
@@ -1353,7 +1376,7 @@ double evalSymGrad(const Basis& basis,
template<class Basis, class Vector, class MatrixFunction>
-double computeL2Error(const Basis& basis,
+double computeL2SymError(const Basis& basis,
Vector& coeffVector,
const double gamma,
const MatrixFunction& matrixFieldFunc
@@ -1590,6 +1613,96 @@ double computeL2ErrorOld(const Basis& basis,
}*/
+/*
+template<class VectorCT>
+double semi_H1_vectorScalarProduct(VectorCT& coeffVector1, VectorCT& coeffVector2)
+{
+
+double ret(0);
+auto localView = basis_.localView();
+
+for (const auto& e : elements(basis_.gridView()))
+ {
+ localView.bind(e);
+
+ auto geometry = e.geometry();
+
+ 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(e.type(), order);
+
+ for (const auto& quadPoint : quad)
+ {
+ double elementScp(0);
+
+ auto pos = quadPoint.position();
+ const auto jacobian = geometry.jacobianInverseTransposed(pos);
+ const double integrationElement = geometry.integrationElement(pos);
+
+ std::vector<FieldMatrix<double,1,dim> > referenceGradients;
+ localFiniteElement.localBasis().evaluateJacobian(pos, referenceGradients);
+
+ std::vector< VectorRT > gradients(referenceGradients.size());
+ for (size_t i=0; i< gradients.size(); i++)
+ jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+ for (size_t i=0; i < nSf; i++)
+ for (size_t j=0; j < nSf; j++)
+ for (size_t k=0; k < dimworld; k++)
+ for (size_t l=0; l < dimworld; l++)
+ {
+ MatrixRT defGradient1(0);
+ defGradient1[k] = gradients[i];
+
+ MatrixRT defGradient2(0);
+ defGradient2[l] = gradients[j];
+
+ size_t localIdx1 = localView.tree().child(k).localIndex(i);
+ size_t globalIdx1 = localView.index(localIdx1);
+
+ size_t localIdx2 = localView.tree().child(l).localIndex(j);
+ size_t globalIdx2 = localView.index(localIdx2);
+
+ double scp(0);
+ for (int ii=0; ii<dimworld; ii++)
+ for (int jj=0; jj<dimworld; jj++)
+ scp += coeffVector1[globalIdx1] * defGradient1[ii][jj] * coeffVector2[globalIdx2] * defGradient2[ii][jj];
+
+ elementScp += scp;
+ }
+
+ ret += elementScp * quadPoint.weight() * integrationElement;
+ } //qp
+ } //e
+ return ret;
+}*/
+
+
+
+
+
+
+
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+
@@ -1763,9 +1876,8 @@ int main(int argc, char *argv[])
Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
- int equivCounter = 0;
+ int equivCounter = 0; // TODO remove
-
// 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))
{
@@ -1781,7 +1893,8 @@ int main(int argc, char *argv[])
// std::cout << "Number of Elements in each direction: " << nE << std::endl;
// std::cout << "Number ofNodes: " << gridView_CE.size(dim) << std::endl;
- auto perTest = periodicIndices.indexPairSet();
+
+ auto perTest = periodicIndices.indexPairSet(); // TODO remove
auto Basis_CE = makeBasis(
@@ -1792,8 +1905,8 @@ int main(int argc, char *argv[])
// blockedInterleaved())); // siehe Periodicbasistest.. funktioniert auch?
// flatInterleaved()));
-// std::cout << "size feBasis: " << Basis_CE.size() << std::endl;
-// std::cout << "basis.dimension()" << Basis_CE.dimension() <<std::endl;
+ std::cout << "size feBasis: " << Basis_CE.size() << std::endl;
+ std::cout << "basis.dimension()" << Basis_CE.dimension() <<std::endl;
log << "size of FiniteElementBasis: " << Basis_CE.size() << std::endl;
const int phiOffset = Basis_CE.size();
@@ -1877,7 +1990,7 @@ int main(int argc, char *argv[])
/////////////////////////////////////////////////////////
// Assemble the system
/////////////////////////////////////////////////////////
- assembleCellStiffness(Basis_CE, muLocal, lambdaLocal, gamma, stiffnessMatrix_CE);
+ assembleCellStiffness(Basis_CE, muLocal, lambdaLocal, gamma, stiffnessMatrix_CE, parameterSet);
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);
@@ -1887,35 +2000,40 @@ int main(int argc, char *argv[])
- /////////////////////////////////////////////////////////
- // Determine global indices of components for arbitrary (local) index TODO :separate Function! arbitraryComponentwiseIndices
- /////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////
+ // Determine global indices of components for arbitrary (local) index
+ //////////////////////////////////////////////////////////////////////
auto localView = Basis_CE.localView();
- int arbitraryIndex = 0;
-
- FieldVector<int,3> row; //fill with Indices..
+ int arbitraryLocalIndex = parameterSet.get<unsigned int>("arbitraryLocalIndex", 0); // localIdx..assert Number < #ShapeFcts .. also input Element?
+ int arbitraryElementNumber = parameterSet.get<unsigned int>("arbitraryElementNumber", 0);
+ FieldVector<int,3> row; //fill with Indices..
+ row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLocalIndex);
+// printvector(std::cout, row, "row" , "--");
- // 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);
+// // 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(arbitraryLocalIndex);
+// 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);
+ setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3, parameterSet);
// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
// printvector(std::cout, load_alpha1, "load_alpha1 after setOneBaseFunctionToZero", "--");
}
@@ -2133,6 +2251,7 @@ int main(int argc, char *argv[])
auto correctorFunction_3 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
Basis_CE,
phi_3);
+
/////////////////////////////////////////////////////////
@@ -2163,8 +2282,8 @@ int main(int argc, char *argv[])
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 << "GGTerm:" << GGterm << std::endl;
+ std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
// std::cout << "Test : " << coeffContainer[a]*loadContainer[b] << std::endl; //same...
}
printmatrix(std::cout, Q, "Matrix Q", "--");
@@ -2186,7 +2305,7 @@ int main(int argc, char *argv[])
- // compute effective Prestrain
+ // Compute effective Prestrain
FieldVector<double, 3> Beff;
Q.solve(Beff,B_hat);
@@ -2261,20 +2380,20 @@ int main(int argc, char *argv[])
if(write_L2Error)
{
- auto L2error = computeL2Error(Basis_CE,phi_1,gamma,symPhi_1_analytic);
+ auto L2error = computeL2SymError(Basis_CE,phi_1,gamma,symPhi_1_analytic);
std::cout << "L2-Error: " << L2error << std::endl;
// auto L2errorTest = computeL2ErrorTest(Basis_CE,phi_1,gamma,symPhi_1_analytic);
// std::cout << "L2-ErrorTEST: " << L2errorTest << std::endl;
- auto L2Norm_Symphi = computeL2Error(Basis_CE,phi_1,gamma,zeroFunction); // Achtung Norm SymPhi_1
- std::cout << "L2-Norm(Symphi_1): " << L2Norm_Symphi << std::endl; // TODO Not Needed
+ auto L2Norm_Symphi = computeL2SymError(Basis_CE,phi_1,gamma,zeroFunction); // Achtung Norm SymPhi_1
+ std::cout << "L2-Norm(Symphi_1): " << L2Norm_Symphi << std::endl; // TODO Not Needed
VectorCT zeroVec;
zeroVec.resize(Basis_CE.size());
zeroVec = 0;
- auto L2Norm_analytic = computeL2Error(Basis_CE,zeroVec ,gamma,symPhi_1_analytic);
+ auto L2Norm_analytic = computeL2SymError(Basis_CE,zeroVec ,gamma,symPhi_1_analytic);
std::cout << "L2-Norm(analytic): " << L2Norm_analytic << std::endl;
log << "L2-Error (symmetric Gradient phi_1):" << L2error << std::endl;