diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 71bb0538e8260d496e97a67f3fa0a3ff18c72de4..c3b7f72cd4ecd6e2baca3b52d51850f5585696fb 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -51,6 +51,7 @@
// #include <dune/fufem/discretizationerror.hh>
// #include <boost/multiprecision/cpp_dec_float.hpp>
// #include <boost/any.hpp>
+
#include <any>
#include <variant>
#include <string>
@@ -59,6 +60,10 @@
using namespace Dune;
+//////////////////////////////////////////////////////////////////////
+// Helper functions for Table-Output
+//////////////////////////////////////////////////////////////////////
+
/*! Center-aligns string within a field of width w. Pads with blank spaces
to enforce alignment. */
std::string center(const std::string s, const int w) {
@@ -77,7 +82,8 @@ std::string center(const std::string s, const int w) {
template<class type>
std::string prd(const type x, const int decDigits, const int width) {
std::stringstream ss;
- ss << std::fixed << std::right;
+// ss << std::fixed << std::right;
+ ss << std::scientific << std::right; // Use scientific Output!
ss.fill(' '); // fill space around displayed #
ss.width(width); // set width around displayed #
ss.precision(decDigits); // set # places after decimal
@@ -89,8 +95,6 @@ std::string prd(const type x, const int decDigits, const int width) {
-
-
template<class Basis>
auto arbitraryComponentwiseIndices(const Basis& basis,
const int elementNumber,
@@ -113,7 +117,7 @@ auto arbitraryComponentwiseIndices(const Basis& basis,
for (int k = 0; k < 3; k++)
{
- auto rowLocal = localView.tree().child(k).localIndex(leafIdx); //Input: LeafIndex! TODO braucht hier (Inverse ) Local-to-Leaf Idx Map
+ auto rowLocal = localView.tree().child(k).localIndex(leafIdx); //Input: LeafIndex! TODO bräuchte hier (Inverse ) Local-to-Leaf Idx Map
// std::cout << "rowLocal:" << rowLocal << std::endl;
row[k] = localView.index(rowLocal);
// std::cout << "row[k]:" << row[k] << std::endl;
@@ -135,7 +139,6 @@ void computeElementStiffnessMatrix(const LocalView& localView,
const double gamma
)
{
-
// Local StiffnessMatrix of the form:
// | phi*phi m*phi |
// | phi *m m*m |
@@ -143,26 +146,18 @@ void computeElementStiffnessMatrix(const LocalView& localView,
using Element = typename LocalView::Element;
const Element element = localView.element();
auto geometry = element.geometry();
-
constexpr int dim = Element::dimension;
- constexpr int nCompo = dim;
-
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
- // using Domain = typename LocalView::GridView::Codim<0>::Geometry::GlobalCoordinate;
- // using FuncScalar = std::function< ScalarRT(const Domain&) >;
- // using Func2Tensor = std::function< Matload_alpha1 ,rixRT(const Domain&) >;
-
+ constexpr int dimWorld = dim;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
elementMatrix.setSize(localView.size()+3, localView.size()+3); //extend by dim ´R_sym^{2x2}
elementMatrix = 0;
-
// LocalBasis-Offset
const int localPhiOffset = localView.size();
- const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
-
// std::cout << "localView.size(): " << localView.size() << std::endl;
// std::cout << "nSf : " << nSf << std::endl;
@@ -173,66 +168,53 @@ void computeElementStiffnessMatrix(const LocalView& localView,
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, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 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", "--");
- // printmatrix(std::cout, basisContainer[1] , "G_2", "--");
-// printmatrix(std::cout, basisContainer[2] , "G_3", "--");
- ////////////////////////////////////////////////////
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-
for (const auto& quadPoint : quad)
{
-
const auto& quadPos = quadPoint.position();
- const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
- const auto integrationElement = geometry.integrationElement(quadPoint.position());
+ const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
+ const auto integrationElement = geometry.integrationElement(quadPos);
std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(quadPoint.position(),
+ localFiniteElement.localBasis().evaluateJacobian(quadPos,
referenceGradients);
-
// Compute the shape function gradients on the grid element
std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
- // std::vector< VectorRT> gradients(referenceGradients.size());
for (size_t i=0; i<gradients.size(); i++)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
-
// "phi*phi"-part
- for (size_t k=0; k < nCompo; k++)
- for (size_t l=0; l< nCompo; l++)
+ for (size_t k=0; k < dimWorld; k++)
+ for (size_t l=0; l< dimWorld; l++)
{
for (size_t i=0; i < nSf; i++)
for (size_t j=0; j < nSf; j++ )
{
-
// (scaled) Deformation gradient of the ansatz basis function
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
-// printmatrix(std::cout, defGradientU , "defGradientU", "--");
+ // printmatrix(std::cout, defGradientU , "defGradientU", "--");
// (scaled) Deformation gradient of the test basis function
MatrixRT defGradientV(0);
defGradientV[l][0] = gradients[j][0]; // Y
defGradientV[l][1] = gradients[j][1]; //X2
defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
-
// symmetric Gradient (Elastic Strains)
MatrixRT strainU, strainV;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
-// strainV[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // same ? genügt strainU
-// printmatrix(std::cout, strainU , "strainU", "--");
+ //printmatrix(std::cout, strainU , "strainU", "--");
}
// St.Venant-Kirchhoff stress
@@ -241,37 +223,29 @@ void computeElementStiffnessMatrix(const LocalView& localView,
MatrixRT stressU(0);
stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
- double trace = 0;
- for (int ii=0; ii<nCompo; ii++)
+ double trace = 0.0;
+ for (int ii=0; ii<dimWorld; ii++)
trace += strainU[ii][ii];
- for (int ii=0; ii<nCompo; ii++)
+ for (int ii=0; ii<dimWorld; ii++)
stressU[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] * 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 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?!
- size_t row = localView.tree().child(l).localIndex(j);
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
+ energyDensity += stressU[ii][jj] * strainV[ii][jj];
+ size_t col = localView.tree().child(k).localIndex(i);
+ size_t row = localView.tree().child(l).localIndex(j);
elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
}
-
}
-
// "m*phi" & "phi*m" -part
- for (size_t l=0; l< nCompo; l++)
+ for (size_t l=0; l< dimWorld; l++)
for (size_t j=0; j < nSf; j++ )
{
-
// (scaled) Deformation gradient of the test basis function
MatrixRT defGradientV(0);
defGradientV[l][0] = gradients[j][0]; // Y
@@ -280,15 +254,14 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// symmetric Gradient (Elastic Strains)
MatrixRT strainV;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
{
strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
}
for( size_t m=0; m<3; m++)
{
-
// < L G_i, sym[D_gamma*nabla phi_j] >
// L G_i* strainV
@@ -297,61 +270,54 @@ void computeElementStiffnessMatrix(const LocalView& localView,
stressG.axpy(2*mu(quadPos), basisContainer[m]); //this += 2mu *strainU
double traceG = 0.0;
- for (int ii=0; ii<nCompo; ii++)
+ for (int ii=0; ii<dimWorld; ii++)
{
-// std::cout << basisContainer[m][ii][ii] << std::endl;
traceG += basisContainer[m][ii][ii];
}
- for (int ii=0; ii<nCompo; ii++)
+ for (int ii=0; ii<dimWorld; ii++)
stressG[ii][ii] += lambda(quadPos) * traceG;
double energyDensityGphi = 0.0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
- energyDensityGphi += stressG[ii][jj] * strainV[ii][jj]; // "m*phi"-part
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
+ energyDensityGphi += stressG[ii][jj] * strainV[ii][jj];
size_t row = localView.tree().child(l).localIndex(j);
-
+
auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
-
elementMatrix[row][localPhiOffset+m] += value;
elementMatrix[localPhiOffset+m][row] += value;
}
-
}
// "m*m"-part
for(size_t m=0; m<3; m++)
for(size_t n=0; n<3; n++)
{
-
// St.Venant-Kirchhoff stress
MatrixRT stressG(0);
stressG.axpy(2*mu(quadPos), basisContainer[m]); //this += 2mu *strainU
double traceG = 0.0;
- for (int ii=0; ii<nCompo; ii++)
+ for (int ii=0; ii<dimWorld; ii++)
traceG += basisContainer[m][ii][ii];
- for (int ii=0; ii<nCompo; ii++)
+ for (int ii=0; ii<dimWorld; ii++)
stressG[ii][ii] += lambda(quadPos) * traceG;
double energyDensityGG = 0.0;
-
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
- energyDensityGG += stressG[ii][jj] * basisContainer[n][ii][jj]; // "m*m"-part
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
+ energyDensityGG += stressG[ii][jj] * basisContainer[n][ii][jj];
elementMatrix[localPhiOffset+m][localPhiOffset+n]= energyDensityGG * quadPoint.weight() * integrationElement;
-
}
}
}
-
// Get the occupation pattern of the stiffness matrix
template<class Basis, class ParameterSet>
void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet& parameterSet)
@@ -361,7 +327,6 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet&
// | phi *m m*m |
auto localView = basis.localView();
-
const int phiOffset = basis.dimension();
nb.resize(basis.size()+3, basis.size()+3);
@@ -380,7 +345,6 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet&
nb.add(row[0],col[0]); // nun würde auch nb.add(row,col) gehen..
}
}
-
for (size_t i=0; i<localView.size(); i++)
{
auto row = localView.index(i);
@@ -400,24 +364,19 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet&
//////////////////////////////////////////////////////////////////
// setOneBaseFunctionToZero
//////////////////////////////////////////////////////////////////
- if(parameterSet.template get<bool>("set_IntegralZero", true)){
+ if(parameterSet.template get<bool>("set_oneBasisFunction_Zero ", true)){
FieldVector<int,3> row;
unsigned int arbitraryLeafIndex = parameterSet.template get<unsigned int>("arbitraryLeafIndex", 0);
unsigned int arbitraryElementNumber = parameterSet.template get<unsigned int>("arbitraryElementNumber", 0);
- const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // macht keinen Unterschied ob hier k oder 0..
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // macht keinen Unterschied ob hier k oder 0..
const auto nSf = localFiniteElement.localBasis().size();
assert(arbitraryLeafIndex < nSf );
-
- std::cout << "arbitraryElementNumber:" << arbitraryElementNumber << std::endl;
- std::cout << "basis.gridView().size(0:" << basis.gridView().size(0) << std::endl;
assert(arbitraryElementNumber < basis.gridView().size(0)); // "arbitraryElementNumber larger than total Number of Elements"
//Determine 3 global indices (one for each component of an arbitrary local FE-function)
row = arbitraryComponentwiseIndices(basis,arbitraryElementNumber,arbitraryLeafIndex);
- printvector(std::cout, row, "row" , "--");
-
for (int k = 0; k<3; k++)
nb.add(row[k],row[k]);
}
@@ -425,7 +384,6 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet&
-
// Compute the source term for a single element
// < L (sym[D_gamma*nabla phi_i] + M_i ), x_3G_alpha >
template<class LocalView, class LocalFunction1, class LocalFunction2, class Vector, class Force>
@@ -445,16 +403,13 @@ void computeElementLoadVector( const LocalView& localView,
const auto element = localView.element();
const auto geometry = element.geometry();
constexpr int dim = Element::dimension;
- constexpr int nCompo = dim;
+ constexpr int dimWorld = dim;
-// using VectorRT = FieldVector< double, nCompo>;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
// Set of shape functions for a single element
const auto& localFiniteElement= localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
- // const auto& localFiniteElement = localView.tree().finiteElement();
-
elementRhs.resize(localView.size() +3);
elementRhs = 0;
@@ -474,10 +429,8 @@ void computeElementLoadVector( const LocalView& localView,
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1); // für simplex
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-
for (const auto& quadPoint : quad)
{
-
const FieldVector<double,dim>& quadPos = quadPoint.position();
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
const double integrationElement = geometry.integrationElement(quadPos);
@@ -486,14 +439,13 @@ void computeElementLoadVector( const LocalView& localView,
localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
- std::vector< FieldVector< double, nCompo> > gradients(referenceGradients.size()); //nComo right? TODO
+ std::vector< FieldVector< double, dim> > gradients(referenceGradients.size());
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++)
+ for (size_t k=0; k < dimWorld; k++)
{
// Deformation gradient of the ansatz basis function
MatrixRT defGradientV(0);
@@ -502,64 +454,58 @@ void computeElementLoadVector( const LocalView& localView,
defGradientV[k][2] = (1.0/gamma)*gradients[i][2]; // X3
// Elastic Strain
- MatrixRT strainV; //strainV never used?
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ MatrixRT strainV;
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
// St.Venant-Kirchhoff stress
MatrixRT stressV(0);
- stressV.axpy(2*mu(quadPos), strainV); //this += 2mu *strainU
+ stressV.axpy(2*mu(quadPos), strainV); //this += 2mu *strainU
double trace = 0.0;
- for (int ii=0; ii<nCompo; ii++)
+ for (int ii=0; ii<dimWorld; ii++)
trace += strainV[ii][ii];
- for (int ii=0; ii<nCompo; ii++)
+ for (int ii=0; ii<dimWorld; ii++)
stressV[ii][ii] += lambda(quadPos) * trace;
-
// Local energy density: stress times strain
double energyDensity = 0.0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
energyDensity += stressV[ii][jj] *forceTerm(quadPos)[ii][jj];
-
size_t row = localView.tree().child(k).localIndex(i);
elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
}
-
// "f*m"-part
for (size_t m=0; m<3; m++)
{
-
// 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++)
+ for (int ii=0; ii<dimWorld; ii++)
traceG += basisContainer[m][ii][ii];
- for (int ii=0; ii<nCompo; ii++)
+ for (int ii=0; ii<dimWorld; ii++)
stressG[ii][ii] += lambda(quadPos) * traceG;
double energyDensityfG = 0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
energyDensityfG += stressG[ii][jj] * forceTerm(quadPos)[ii][jj];
elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
}
-
}
}
-
template<class Basis, class Matrix, class LocalFunction1, class LocalFunction2, class ParameterSet>
void assembleCellStiffness(const Basis& basis,
LocalFunction1& muLocal,
@@ -586,18 +532,12 @@ void assembleCellStiffness(const Basis& basis,
lambdaLocal.bind(element);
const int localPhiOffset = localView.size();
- // std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
-
- // Dune::Matrix<double> elementMatrix;
+ //std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
-// printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
-// localdebugLog << elementMatrix << std::endl;
-
-
- // std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
- // std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
-
+ //printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
+ //std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
+ //std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
//////////////////////////////////////////////////////////////////////////////
// GLOBAL STIFFNES ASSEMBLY
@@ -621,15 +561,12 @@ void assembleCellStiffness(const Basis& basis,
matrix[row][phiOffset+m] += elementMatrix[i][localPhiOffset+m];
matrix[phiOffset+m][row] += elementMatrix[localPhiOffset+m][i];
}
-
for (size_t m=0; m<3; m++ )
for (size_t n=0; n<3; n++ )
{
matrix[phiOffset+m][phiOffset+n] += elementMatrix[localPhiOffset+m][localPhiOffset+n];
}
-
// printmatrix(std::cout, matrix, "StiffnessMatrix", "--");
-
}
}
@@ -645,21 +582,18 @@ void assembleCellLoad(const Basis& basis,
)
{
// std::cout << "assemble load vector." << std::endl;
-
b.resize(basis.size()+3);
b = 0;
auto localView = basis.localView();
-
const int phiOffset = basis.dimension();
- // Transform G_alpha's to GridViewFunctions/LocalFunctions -- why exactly?
+ // Transform G_alpha's to GridViewFunctions/LocalFunctions
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);
@@ -668,12 +602,10 @@ void assembleCellLoad(const Basis& basis,
const int localPhiOffset = localView.size();
// std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
-
// BlockVector<double> elementRhs;
BlockVector<FieldVector<double,1> > elementRhs;
computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, loadFunctional);
// printvector(std::cout, elementRhs, "elementRhs", "--");
-// computeElementLoadTEST(localView, muLocal, lambdaLocal, gamma, elementRhs, loadFunctional);
// printvector(std::cout, elementRhs, "elementRhs", "--");
// LOAD ASSEMBLY
@@ -681,26 +613,18 @@ void assembleCellLoad(const Basis& basis,
{
// The global index of the p-th vertex of the element
auto row = localView.index(p);
-
b[row] += elementRhs[p];
}
-
for (size_t m=0; m<3; m++ )
{
b[phiOffset+m] += elementRhs[localPhiOffset+m];
}
-
-
// printvector(std::cout, b, "b", "--");
}
-
}
-
-
-
template<class Basis, class LocalFunction1, class LocalFunction2, class MatrixFunction>
auto energy(const Basis& basis,
LocalFunction1& muLocal,
@@ -709,13 +633,13 @@ auto energy(const Basis& basis,
const MatrixFunction& matrixFieldFuncB)
{
- // TEST HIGHER PRECISION
+// TEST HIGHER PRECISION
// using float_50 = boost::multiprecision::cpp_dec_float_50;
// float_50 higherPrecEnergy = 0.0;
double energy = 0.0;
constexpr int dim = 3;
- constexpr int nCompo = 3;
+ constexpr int dimWorld = 3;
auto localView = basis.localView();
@@ -726,16 +650,12 @@ auto energy(const Basis& basis,
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
-
-
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
- // TEST
+// TEST
// FieldVector<double,3> testvector = {1.0 , 1.0 , 1.0};
// printmatrix(std::cout, matrixFieldFuncB(testvector) , "matrixFieldB(testvector) ", "--");
-
-
for (const auto& e : elements(basis.gridView()))
{
localView.bind(e);
@@ -744,9 +664,6 @@ auto energy(const Basis& basis,
muLocal.bind(e);
lambdaLocal.bind(e);
-
-// FieldVector<double,1> elementEnergy(0); //!!!
-// printvector(std::cout, elementEnergy, "elementEnergy" , "--");
double elementEnergy = 0.0;
// double elementEnergy_HP = 0.0;
@@ -757,53 +674,45 @@ auto energy(const Basis& basis,
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
- for (const auto& quadPoint : quad) {
-// for (size_t pt=0; pt < quad.size(); pt++) {
-
+ for (const auto& quadPoint : quad)
+ {
-// const FieldVector<double,dim>& quadPos = quadPoint.position();
const auto& quadPos = quadPoint.position();
-// const Domain& quadPos = quad[pt].position();
const double integrationElement = geometry.integrationElement(quadPos);
-
auto strain1 = matrixFieldA(quadPos);
+ auto strain2 = matrixFieldB(quadPos);
// St.Venant-Kirchhoff stress
MatrixRT stressU(0.0);
- stressU.axpy(2.0*muLocal(quadPos), strain1); //this += 2mu *strainU // eigentlich this += 2mu *strain1 ?
+ stressU.axpy(2.0*muLocal(quadPos), strain1); //this += 2mu *strainU
double trace = 0.0;
- for (int ii=0; ii<nCompo; ii++)
+ for (int ii=0; ii<dimWorld; ii++)
trace += strain1[ii][ii];
- for (int ii=0; ii<nCompo; ii++)
+ for (int ii=0; ii<dimWorld; ii++)
stressU[ii][ii] += lambdaLocal(quadPos) * trace;
- auto strain2 = matrixFieldB(quadPos);
-
// Local energy density: stress times strain
double energyDensity = 0.0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
energyDensity += stressU[ii][jj] * strain2[ii][jj];
-// elementEnergy += energyDensity * quad[pt].weight() * integrationElement;
- elementEnergy += energyDensity * quadPoint.weight() * integrationElement; //TODO integratonElement needed here?
-// elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
+ elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
+ //elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
}
energy += elementEnergy;
-// higherPrecEnergy += elementEnergy_HP;
+ //higherPrecEnergy += elementEnergy_HP;
}
- //TEST
+// TEST
// std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
return energy;
}
-
-
template<class Basis, class Matrix, class Vector, class ParameterSet>
void setOneBaseFunctionToZero(const Basis& basis,
Matrix& stiffnessMatrix,
@@ -820,7 +729,6 @@ void setOneBaseFunctionToZero(const Basis& basis,
FieldVector<int,3> row;
unsigned int arbitraryLeafIndex = parameterSet.template get<unsigned int>("arbitraryLeafIndex", 0);
unsigned int arbitraryElementNumber = parameterSet.template get<unsigned int>("arbitraryElementNumber", 0);
-
//Determine 3 global indices (one for each component of an arbitrary local FE-function)
row = arbitraryComponentwiseIndices(basis,arbitraryElementNumber,arbitraryLeafIndex);
@@ -838,17 +746,13 @@ 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
+ // Output : determine all Indices that belong to that component
auto localView = basis.localView();
@@ -861,8 +765,8 @@ auto childToIndexMap(const Basis& basis,
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(); //
+ const auto& localFiniteElement = localView.tree().child(k).finiteElement();
+ const auto nSf = localFiniteElement.localBasis().size();
for(size_t j=0; j<nSf; j++)
{
@@ -883,11 +787,6 @@ auto childToIndexMap(const Basis& basis,
}
-
-
-
-
-
template<class Basis, class Vector>
auto integralMean(const Basis& basis,
Vector& coeffVector
@@ -916,10 +815,10 @@ auto integralMean(const Basis& basis,
for(const auto& quadPoint : quad)
{
- const double integrationElement = element.geometry().integrationElement(quadPoint.position());
- area += quadPoint.weight() * integrationElement;
-
const auto& quadPos = quadPoint.position();
+ const double integrationElement = element.geometry().integrationElement(quadPos);
+ area += quadPoint.weight() * integrationElement;
+
r += localfun(quadPos) * quadPoint.weight() * integrationElement;
}
}
@@ -928,26 +827,18 @@ auto integralMean(const Basis& basis,
}
-
-
-
-
-
-
template<class Basis, class Vector>
auto subtractIntegralMean(const Basis& basis,
Vector& coeffVector
)
{
// Substract correct Integral mean from each associated component function
-
-// auto IM = integralMean(basis, fun);
auto IM = integralMean(basis, coeffVector);
constexpr int dim = 3;
for(size_t k=0; k<dim; k++)
{
- // std::cout << "Integral-Mean: " << IM[k] << std::endl;
+ //std::cout << "Integral-Mean: " << IM[k] << std::endl;
auto idx = childToIndexMap(basis,k);
for ( int i : idx)
@@ -957,25 +848,23 @@ auto subtractIntegralMean(const Basis& basis,
- //////////////////////////////////////////////////
- // Infrastructure for handling periodicity
- //////////////////////////////////////////////////
+//////////////////////////////////////////////////
+// Infrastructure for handling periodicity
+//////////////////////////////////////////////////
// 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)
- {
+ {
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]))
- and (FloatCmp::eq(x[2],y[2]))
- );
- };
-
-
-
-
+ and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1]))
+ and (FloatCmp::eq(x[2],y[2]))
+ );
+ };
+
+////////////////////////////////////////////////////////////// L2-ERROR /////////////////////////////////////////////////////////////////
template<class Basis, class Vector, class MatrixFunction>
@@ -987,16 +876,14 @@ double computeL2SymErrorNew(const Basis& basis,
{
double error = 0.0;
constexpr int dim = 3;
- constexpr int nCompo = 3;
+ constexpr int dimWorld = 3;
auto localView = basis.localView();
auto matrixFieldGVF = Dune::Functions::makeGridViewFunction(matrixFieldFunc, basis.gridView());
auto matrixField = localFunction(matrixFieldGVF);
-// using GridView = typename Basis::GridView;
-// using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
for (const auto& element : elements(basis.gridView()))
{
@@ -1004,37 +891,32 @@ double computeL2SymErrorNew(const Basis& basis,
matrixField.bind(element);
auto geometry = element.geometry();
- const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
-// std::cout << "print nSf:" << nSf << std::endl;
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); //TEST
+
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 ); //TEST
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-
for (const auto& quadPoint : quad)
{
-
const auto& quadPos = quadPoint.position();
- const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
+ const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
const auto integrationElement = geometry.integrationElement(quadPoint.position());
std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(quadPoint.position(),
+ localFiniteElement.localBasis().evaluateJacobian(quadPos,
referenceGradients);
// Compute the shape function gradients on the grid element
std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
- // std::vector< VectorRT> gradients(referenceGradients.size());
for (size_t i=0; i<gradients.size(); i++)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
-
- MatrixRT tmp(0);
+ MatrixRT tmp(0);
double sum = 0.0;
- for (size_t k=0; k < nCompo; k++)
+ for (size_t k=0; k < dimWorld; k++)
for (size_t i=0; i < nSf; i++)
{
@@ -1047,31 +929,26 @@ double computeL2SymErrorNew(const Basis& basis,
defGradientU[k][1] = coeffVector[globalIdx1]*gradients[i][1]; //X2
defGradientU[k][2] = coeffVector[globalIdx1]*(1.0/gamma)*gradients[i][2]; //X3
-// printvector(std::cout, gradients[i], "gradients[i]", "--");
-// std::cout << "coeffVector[globalIdx1]" << coeffVector[globalIdx1] << std::endl;
-
// symmetric Gradient (Elastic Strains)
MatrixRT strainU(0);
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
}
-
// printmatrix(std::cout, strainU, "strainU", "--");
// printmatrix(std::cout, tmp, "tmp", "--");
tmp += strainU;
}
// printmatrix(std::cout, matrixField(quadPos), "matrixField(quadPos)", "--");
-// printmatrix(std::cout, tmp, "tmp", "--"); // TEST Symphi_1 hat ganz andere Struktur als analytic?
+// printmatrix(std::cout, tmp, "tmp", "--"); // TEST Symphi_1 hat ganz andere Struktur als analytic?
// tmp = tmp - matrixField(quadPos);
// printmatrix(std::cout, tmp - matrixField(quadPos), "Difference", "--");
-
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
{
sum += std::pow(tmp[ii][jj]-matrixField(quadPos)[ii][jj],2);
}
@@ -1099,7 +976,7 @@ double computeL2SymError(const Basis& basis,
// This Version computes the SCALAR PRODUCT of the difference ..
double error = 0.0;
constexpr int dim = 3;
- constexpr int nCompo = 3;
+ constexpr int dimWorld = 3;
auto localView = basis.localView();
@@ -1108,7 +985,7 @@ double computeL2SymError(const Basis& basis,
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
for (const auto& element : elements(basis.gridView()))
{
@@ -1116,7 +993,7 @@ double computeL2SymError(const Basis& basis,
matrixField.bind(element);
auto geometry = element.geometry();
- const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
// std::cout << "print nSf:" << nSf << std::endl;
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
@@ -1138,16 +1015,16 @@ double computeL2SymError(const Basis& basis,
// std::vector< VectorRT> gradients(referenceGradients.size());
for (size_t i=0; i<gradients.size(); i++)
- jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]); // TRANSFORMED
+ jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
// MatrixRT defGradientU(0), defGradientV(0);
- for (size_t k=0; k < nCompo; k++)
+ for (size_t k=0; k < dimWorld; k++)
for (size_t i=0; i < nSf; i++)
{
- for (size_t l=0; l < nCompo; l++)
+ for (size_t l=0; l < dimWorld; l++)
for (size_t j=0; j < nSf; j++ )
{
@@ -1169,8 +1046,8 @@ double computeL2SymError(const Basis& basis,
// symmetric Gradient (Elastic Strains)
MatrixRT strainU(0), strainV(0);
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
@@ -1178,47 +1055,46 @@ double computeL2SymError(const Basis& basis,
// Local energy density: stress times strain
double tmp = 0.0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
tmp += coeffVector[globalIdx1]*coeffVector[globalIdx2]* strainU[ii][jj] * strainV[ii][jj];
-
error += tmp * quadPoint.weight() * integrationElement;
}
}
- for (size_t k=0; k < nCompo; k++)
+ for (size_t k=0; k < dimWorld; k++)
for (size_t i=0; i < nSf; i++)
{
size_t localIdx1 = localView.tree().child(k).localIndex(i);
size_t globalIdx1 = localView.index(localIdx1);
// (scaled) Deformation gradient of the ansatz basis function
- MatrixRT defGradientU(0); //TEST (doesnt change anything..)
+ 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
// symmetric Gradient (Elastic Strains)
MatrixRT strainU(0);
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
}
double tmp = 0.0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
tmp += (-2.0)*coeffVector[globalIdx1]*strainU[ii][jj] * matrixField(quadPos)[ii][jj];
error += tmp * quadPoint.weight() * integrationElement;
}
double tmp = 0.0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
tmp += matrixField(quadPos)[ii][jj] * matrixField(quadPos)[ii][jj];
error += tmp * quadPoint.weight() * integrationElement;
@@ -1230,9 +1106,9 @@ double computeL2SymError(const Basis& basis,
+ ////////////////////////////////////////////////////////////// L2-NORM of symmetric Gradient /////////////////////////////////////////////////////////////////
-// TEST
template<class Basis, class Vector>
double computeL2SymNormCoeff(const Basis& basis,
Vector& coeffVector,
@@ -1241,20 +1117,20 @@ double computeL2SymNormCoeff(const Basis& basis,
{
double error = 0.0;
constexpr int dim = 3;
- constexpr int nCompo = 3;
+ constexpr int dimWorld = 3;
auto localView = basis.localView();
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
for (const auto& element : elements(basis.gridView()))
{
localView.bind(element);
auto geometry = element.geometry();
- const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
// std::cout << "print nSf:" << nSf << std::endl;
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
@@ -1273,18 +1149,16 @@ double computeL2SymNormCoeff(const Basis& basis,
// Compute the shape function gradients on the grid element
std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
- // std::vector< VectorRT> gradients(referenceGradients.size());
for (size_t i=0; i<gradients.size(); i++)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
MatrixRT tmp(0);
-
double sum = 0.0;
for (size_t i=0; i < nSf; i++)
- for (size_t k=0; k < nCompo; k++)
+ for (size_t k=0; k < dimWorld; k++)
{
size_t localIdx1 = localView.tree().child(k).localIndex(i); // hier i:leafIdx
@@ -1298,8 +1172,8 @@ double computeL2SymNormCoeff(const Basis& basis,
// symmetric Gradient (Elastic Strains)
MatrixRT strainU(0);
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
}
@@ -1312,8 +1186,8 @@ double computeL2SymNormCoeff(const Basis& basis,
// Local energy density: stress times strain
// double tmp = 0.0;
-// for (int ii=0; ii<nCompo; ii++)
-// for (int jj=0; jj<nCompo; jj++)
+// for (int ii=0; ii<dimWorld; ii++)
+// for (int jj=0; jj<dimWorld; jj++)
// tmp[ii][jj] += coeffVector[globalIdx1]*coeffVector[globalIdx2]* strainU[ii][jj] * strainV[ii][jj];
//
//
@@ -1321,8 +1195,8 @@ double computeL2SymNormCoeff(const Basis& basis,
// }
}
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
sum += std::pow(tmp[ii][jj],2);
error += sum * quadPoint.weight() * integrationElement;
@@ -1350,13 +1224,13 @@ double computeL2SymNormCoeffV2(const Basis& basis,
{
double error = 0.0;
constexpr int dim = 3;
- constexpr int nCompo = 3;
+ constexpr int dimWorld = 3;
auto localView = basis.localView();
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
for (const auto& element : elements(basis.gridView()))
{
@@ -1391,13 +1265,10 @@ double computeL2SymNormCoeffV2(const Basis& basis,
// MatrixRT defGradientU(0), defGradientV(0);
- for (size_t k=0; k < nCompo; k++)
-
-
-
+ for (size_t k=0; k < dimWorld; k++)
for (size_t i=0; i < nSf; i++)
{
- for (size_t l=0; l < nCompo; l++)
+ for (size_t l=0; l < dimWorld; l++)
for (size_t j=0; j < nSf; j++ )
{
@@ -1419,8 +1290,8 @@ double computeL2SymNormCoeffV2(const Basis& basis,
// symmetric Gradient (Elastic Strains)
MatrixRT strainU(0), strainV(0);
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
@@ -1428,8 +1299,8 @@ double computeL2SymNormCoeffV2(const Basis& basis,
// Local energy density: stress times strain
double tmp = 0.0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
tmp += coeffVector[globalIdx1]*coeffVector[globalIdx2]* strainU[ii][jj] * strainV[ii][jj];
@@ -1445,108 +1316,10 @@ double computeL2SymNormCoeffV2(const Basis& basis,
}
-/*
-
-
-
-
-template<class Basis, class Vector>
-double computeL2SymNormCoeffV3(const Basis& basis,
- Vector& coeffVector,
- const double gamma
- )
-{
- double error = 0.0;
- constexpr int dim = 3;
- constexpr int nCompo = 3;
-
- auto localView = basis.localView();
- using GridView = typename Basis::GridView;
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
- for (const auto& element : elements(basis.gridView()))
- {
- localView.bind(element);
- auto geometry = element.geometry();
- const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
- const auto nSf = localFiniteElement.localBasis().size();
-// std::cout << "print nSf:" << nSf << std::endl;
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
- const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
- for (const auto& quadPoint : quad)
- {
-
- const auto& quadPos = quadPoint.position();
- const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
- const auto integrationElement = geometry.integrationElement(quadPoint.position());
-
- std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(quadPoint.position(),
- referenceGradients);
-
- // Compute the shape function gradients on the grid element
- std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
- // std::vector< VectorRT> gradients(referenceGradients.size());
-
- for (size_t i=0; i<gradients.size(); i++)
- jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
-
-
- MatrixRT defGradientU(0), defGradientV(0);
-
-
- for (size_t k=0; k < nCompo; k++)
- for (size_t i=0; i < nSf; i++)
- {
- for (size_t l=0; l < nCompo; l++)
- for (size_t j=0; j < nSf; j++ )
- {
-
- size_t localIdx1 = localView.tree().child(k).localIndex(i); // hier i:leafIdx
- size_t globalIdx1 = localView.index(localIdx1);
- size_t localIdx2 = localView.tree().child(l).localIndex(j);
- size_t globalIdx2 = localView.index(localIdx2);
-
- // (scaled) Deformation gradient of the ansatz basis function
- defGradientU[k][0] = gradients[i][0]; // Y //hier i:leaf oder localIdx?
- defGradientU[k][1] = gradients[i][1]; //X2
- defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
-
- defGradientV[l][0] = gradients[j][0]; // Y
- defGradientV[l][1] = gradients[j][1]; //X2
- defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
-
- // symmetric Gradient (Elastic Strains)
- MatrixRT strainU(0), strainV(0);
- 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]);
- }
-
- // Local energy density: stress times strain
- double tmp = 0.0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
- tmp += coeffVector[globalIdx1]*coeffVector[globalIdx2]* strainU[ii][jj] * strainV[ii][jj];
-
-
- error += tmp * quadPoint.weight() * integrationElement;
- }
- }
-
-
- }
- }
-
- return sqrt(error);
-}
-*/
////////////////////////////////////////////////////////////// L2-NORM /////////////////////////////////////////////////////////////////
@@ -1558,7 +1331,7 @@ double computeL2NormCoeff(const Basis& basis,
{
double error = 0.0;
constexpr int dim = 3;
- constexpr int nCompo = 3;
+ constexpr int dimWorld = 3;
auto localView = basis.localView();
for (const auto& element : elements(basis.gridView()))
@@ -1594,12 +1367,12 @@ double computeL2NormCoeff(const Basis& basis,
std::vector<ValueType> values;*/
double tmp = 0.0;
- for (size_t k=0; k < nCompo; k++) // ANALOG STOKES Beitrag nur wenn k=l!!!
+ for (size_t k=0; k < dimWorld; k++) // ANALOG STOKES Beitrag nur wenn k=l!!!
{
double comp = 0.0;
for (size_t i=0; i < nSf; i++)
{
- // for (size_t l=0; l < nCompo; l++)
+ // for (size_t l=0; l < dimWorld; l++)
// for (size_t j=0; j < nSf; j++ )
// {
size_t localIdx1 = localView.tree().child(k).localIndex(i); // hier i:leafIdx
@@ -1634,12 +1407,12 @@ double computeL2NormCoeffV2(const Basis& basis,
{
double error = 0.0;
constexpr int dim = 3;
- constexpr int nCompo = 3;
+ constexpr int dimWorld = 3;
auto localView = basis.localView();
// using GridView = typename Basis::GridView;
// using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
-// using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+// using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
for (const auto& element : elements(basis.gridView()))
{
@@ -1673,10 +1446,10 @@ double computeL2NormCoeffV2(const Basis& basis,
using ValueType = LocalFEType::Traits::LocalBasisType::Traits::RangeType;
std::vector<ValueType> values;*/
- for (size_t k=0; k < nCompo; k++) // ANALOG STOKES Beitrag nur wenn k=l!!!
+ for (size_t k=0; k < dimWorld; k++) // ANALOG STOKES Beitrag nur wenn k=l!!!
for (size_t i=0; i < nSf; i++)
{
-// for (size_t l=0; l < nCompo; l++)
+// for (size_t l=0; l < dimWorld; l++)
for (size_t j=0; j < nSf; j++ )
{
size_t localIdx1 = localView.tree().child(k).localIndex(i); // hier i:leafIdx
@@ -1709,12 +1482,12 @@ double computeL2NormCoeffV3(const Basis& basis, // Falsch: betrachtet k
{
double error = 0.0;
constexpr int dim = 3;
- constexpr int nCompo = 3;
+ constexpr int dimWorld = 3;
auto localView = basis.localView();
// using GridView = typename Basis::GridView;
// using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
-// using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+// using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
for (const auto& element : elements(basis.gridView()))
{
@@ -1750,7 +1523,7 @@ double computeL2NormCoeffV3(const Basis& basis, // Falsch: betrachtet k
using ValueType = LocalFEType::Traits::LocalBasisType::Traits::RangeType;
std::vector<ValueType> values;*/
- for (size_t k=0; k < nCompo; k++)
+ for (size_t k=0; k < dimWorld; k++)
for (size_t i=0; i < nSf; i++)
{
size_t localIdx1 = localView.tree().child(k).localIndex(i); // hier i:leafIdx
@@ -1778,7 +1551,7 @@ double computeL2NormFunc(const Basis& basis,
// IMPLEMENTATION with makeDiscreteGlobalBasisFunction
double error = 0.0;
constexpr int dim = 3;
- constexpr int nCompo = 3;
+ constexpr int dimWorld = 3;
auto localView = basis.localView();
@@ -1830,7 +1603,7 @@ double computeL2ErrorOld(const Basis& basis,
auto error = 0.0;
constexpr int dim = 3;
- constexpr int nCompo = 3;
+ constexpr int dimWorld = 3;
auto localView = basis.localView();
@@ -1839,7 +1612,7 @@ double computeL2ErrorOld(const Basis& basis,
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
for (const auto& element : elements(basis.gridView()))
@@ -1878,7 +1651,7 @@ double computeL2ErrorOld(const Basis& basis,
MatrixRT defGradientU(0);
- for (size_t k=0; k < nCompo; k++)
+ for (size_t k=0; k < dimWorld; k++)
for (size_t i=0; i < nSf; i++)
{
@@ -1892,8 +1665,8 @@ double computeL2ErrorOld(const Basis& basis,
// symmetric Gradient (Elastic Strains)
MatrixRT strainU(0);
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
{
strainU[ii][jj] = coeffVector[globalIdx] * 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
// strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient //TEST
@@ -1901,8 +1674,8 @@ double computeL2ErrorOld(const Basis& basis,
// Local energy density: stress times strain
double tmp = 0.0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
tmp+= std::pow(strainU[ii][jj] - matrixField(quadPos)[ii][jj],2);
// tmp+= std::pow((coeffVector[globalIdx]*strainU[ii][jj]) - matrixField(quadPos)[ii][jj] ,2);
@@ -1993,7 +1766,7 @@ double computeL2SymErrorSPTest(const Basis& basis,
// This Version computes the SCALAR PRODUCT of the difference ..
double error = 0.0;
constexpr int dim = 3;
- constexpr int nCompo = 3;
+ constexpr int dimWorld = 3;
auto localView = basis.localView();
@@ -2002,7 +1775,7 @@ double computeL2SymErrorSPTest(const Basis& basis,
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
Vector zeroVec;
zeroVec.resize(basis.size());
@@ -2027,8 +1800,8 @@ double computeL2SymErrorSPTest(const Basis& basis,
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-
-
+
+
for (const auto& quadPoint : quad)
{
const auto& quadPos = quadPoint.position();
@@ -2047,7 +1820,7 @@ double computeL2SymErrorSPTest(const Basis& basis,
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]); // TRANSFORMED
- for (size_t k=0; k < nCompo; k++)
+ for (size_t k=0; k < dimWorld; k++)
for (size_t i=0; i < nSf; i++)
{
size_t localIdx1 = localView.tree().child(k).localIndex(i);
@@ -2061,15 +1834,15 @@ double computeL2SymErrorSPTest(const Basis& basis,
// symmetric Gradient (Elastic Strains)
MatrixRT strainU(0);
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
}
double tmp = 0.0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
+ for (int ii=0; ii<dimWorld; ii++)
+ for (int jj=0; jj<dimWorld; jj++)
tmp += (-2.0)*coeffVector[globalIdx1]*strainU[ii][jj] * matrixField(quadPos)[ii][jj];
error += tmp * quadPoint.weight() * integrationElement;
@@ -2104,42 +1877,41 @@ int main(int argc, char *argv[])
ParameterTreeParser::readOptions(argc, argv, parameterSet);
}
- // output setter
+ // Output setter
std::string outputPath = parameterSet.get("outputPath", "../../outputs/output.txt");
std::fstream log;
- std::fstream debugLog;
-// log.open("output2.txt", std::ios::out);
log.open(outputPath ,std::ios::out);
- debugLog.open("../../outputs/debugLog.txt",std::ios::out);
-// log << "Writing this to a file. \n";
-// log.close();
-
constexpr int dim = 3;
- constexpr int nCompo = 3;
- constexpr int order_CE = 1;
+ constexpr int dimWorld = 3;
///////////////////////////////////
// Get Parameters/Data
///////////////////////////////////
double gamma = parameterSet.get<double>("gamma",1.0); // ratio dimension reduction to homogenization
-
+ double alpha = parameterSet.get<double>("alpha", 2.0);
+ double theta = parameterSet.get<double>("theta",1.0/4.0);
+ ///////////////////////////////////
+ // Get Material Parameters
+ ///////////////////////////////////
+ double beta = parameterSet.get<double>("beta",2.0);
+ double mu1 = parameterSet.get<double>("mu1",10.0);;
+ double mu2 = beta*mu1;
+ double lambda1 = parameterSet.get<double>("lambda1",0.0);;
+ double lambda2 = beta*lambda1;
+
// 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
// double eps = parameterSet.get<double>("eps", 0.1); //size of perioticity cell
-
///////////////////////////////////
- // Get Prestrain Parameters
+ // Get Prestrain/Parameters
///////////////////////////////////
unsigned int prestraintype = parameterSet.get<unsigned int>("prestrainType", 2);
double rho1 = parameterSet.get<double>("rho1", 1.0);
- double alpha = parameterSet.get<double>("alpha", 2.0);
- double theta = parameterSet.get<double>("theta",0.3); //TEST theta = 1.0/4.0
- double rho2 = alpha*1.0;
-
+ double rho2 = alpha*rho1;
auto prestrainImp = PrestrainImp(rho1, rho2, theta, width);
auto B_Term = prestrainImp.getPrestrain(prestraintype);
@@ -2147,7 +1919,10 @@ int main(int argc, char *argv[])
log << "alpha: " << alpha << std::endl;
log << "gamma: " << gamma << std::endl;
log << "theta: " << theta << std::endl;
-
+ log << "beta: " << beta << std::endl;
+ log << "material parameters: " << std::endl;
+ log << "mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
+ log << "lambda1: " << lambda1 <<"\nlambda2:" << lambda2 << std::endl;
///////////////////////////////////
// Generate the grid
@@ -2158,7 +1933,6 @@ int main(int argc, char *argv[])
// (shifted) Domain (-1/2,1/2)^3
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});
-
if (cellDomain == 2)
{
// Domain : [0,1)^2 x (-1/2, 1/2)
@@ -2166,66 +1940,12 @@ int main(int argc, char *argv[])
FieldVector<double,dim> upper({1.0, 1.0, 1.0/2.0});
}
-
-
-
-// std::array<int, dim> nElements = parameterSet.get<std::array<int,dim>>("nElements_Cell", {10, 10, 10});
-
-
- ///// LEVEL - APPROACH // TODO
-// int numLevels = parameterSet.get<int>("numLevels", 1);
- std::array<int,2> numLevels = parameterSet.get<std::array<int,2>>("numLevels", {1,3}); // Alternativ : Bereich z.B. [2,4]
+ std::array<int,2> numLevels = parameterSet.get<std::array<int,2>>("numLevels", {1,3});
int levelCounter = 0;
- for(const int &lev : numLevels)
- std::cout << "value of numLevels: " << lev << std::endl;
-
-
-
-
-// // any type
-// std::any a = 1;
-// std::cout << a.type().name() << ": " << std::any_cast<int>(a) << '\n';
-// std::vector<std::any> Storage;
-// std::any Test=1;
-// Storage[0]= 1;
-// std::cout << std::any_cast<int>(Test) << '\n';
-// std::cout << std::any_cast<int>(Storage[0]) << '\n';
-//
-// // Test = "Test";
-//
-// Test = std::make_any<std::string>("Hello World");
-//
-// std::cout << "Test: " << std::any_cast<std::string>(Test) << std::endl;
-//
-// // Test = 5;
-// // std::cout << "Test: " << std::any_cast<int>(Test) << std::endl;
-//
-// Test = 1e-03;
-//
-//
-// std::cout << "Test: " << std::any_cast<double>(Test) << std::endl;
-
- // std::variant
- std::variant<int,double, std::string> a;
- std::variant<std::string> x("abc");
-
-
- std::variant<int, std::string> variant = "Hello";
-
- std::string string_1 = std::get<std::string>( variant ); // get value by type
- std::string string_2 = std::get<1>( variant ); // get value by index
- std::cout << string_1 << std::endl;
- std::cout << string_2 << std::endl;
-
- std::cout << "output variant :" << std::get<std::string>(x) << std::endl;
-
-
-
-
///////////////////////////////////
// Create Data Storage
///////////////////////////////////
@@ -2234,7 +1954,7 @@ int main(int argc, char *argv[])
std::vector<std::variant<std::string, size_t , double>> Storage_Error;
// Storage:: #1 level #2 |q1_a-q1_c| #3 |q2_a-q2_c| #4 |q3_a-q3_c| #5 |b1_a-b1_c| #6 |b2_a-b2_c| #7 |b3_a-b3_c|
- std::vector<std::variant<std::string, size_t , double>> Storage_Quantities; // TODO :better error ?
+ std::vector<std::variant<std::string, size_t , double>> Storage_Quantities;
@@ -2251,35 +1971,15 @@ int main(int argc, char *argv[])
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-// for(int level = 0; level < numLevels; level++)
// for(const size_t &level : numLevels) // explixite Angabe der levels.. {2,4}
for(size_t level = numLevels[0] ; level <= numLevels[1]; level++) // levels von bis.. [2,4]
{
-
std::cout << " ----------------------------------" << std::endl;
std::cout << "Level: " << level << std::endl;
std::cout << " ----------------------------------" << std::endl;
Storage_Error.push_back(level);
Storage_Quantities.push_back(level);
-
std::array<int, dim> nElements = { (int)std::pow(2,level) , (int)std::pow(2,level) , (int)std::pow(2,level) };
std::cout << "Number of Elements in each direction: " << nElements << std::endl;
@@ -2291,57 +1991,37 @@ int main(int argc, char *argv[])
const GridView gridView_CE = grid_CE.leafGridView();
// using ScalarRT = FieldVector< double, 1>;
- // using VectorRT = FieldVector< double, nCompo>;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+ // using VectorRT = FieldVector< double, dimWorld>;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
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> >;
using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >;
- ///////////////////////////////////
- // Get Material Parameters
- ///////////////////////////////////
- double beta = parameterSet.get<double>("beta",2.0);
- double mu1 = parameterSet.get<double>("mu1",10.0);;
- double mu2 = beta*mu1;
- double lambda1 = parameterSet.get<double>("lambda1",0.0);;
- double lambda2 = beta*lambda1;
-
///////////////////////////////////
// Create Lambda-Functions for material Parameters depending on microstructure
///////////////////////////////////
auto muTerm = [mu1, mu2, theta] (const Domain& z) {
- if (abs(z[0]) >= (theta/2.0)) //TODO check ..nullset/boundary
+ if (abs(z[0]) >= (theta/2.0))
return mu1;
- // if (abs(z[0]) > (theta/2.0)) //TEST include boundary (indicatorFunction)
- // return mu1;
- // if (abs(z[0]) < (theta/2) && z[2] < 0) // oder das?
- // return mu2;
else
return mu2;
- // return mu1;
};
auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) {
- if (abs(z[0]) >= (theta/2.0))
+ if (abs(z[0]) >= (theta/2.0))
return lambda1;
- else
+ else
return lambda2;
- };
+ };
auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
auto muLocal = localFunction(muGridF);
auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
auto lambdaLocal = localFunction(lambdaGridF);
- log << "beta: " << beta << std::endl;
- log << "material parameters: " << std::endl;
- log << "mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
- log << "lambda1: " << lambda1 <<"\nlambda2:" << lambda2 << std::endl;
-
///////////////////////////////////
// --- Choose Solver ---
@@ -2350,6 +2030,10 @@ int main(int argc, char *argv[])
// 3 : QR
///////////////////////////////////
unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 1);
+
+
+
+ bool print_debug = parameterSet.get<bool>("print_debug", false);
bool write_corrector_phi1 = parameterSet.get<bool>("write_corrector_phi1", false);
bool write_corrector_phi2 = parameterSet.get<bool>("write_corrector_phi2", false);
bool write_corrector_phi3 = parameterSet.get<bool>("write_corrector_phi3", false);
@@ -2364,43 +2048,28 @@ int main(int argc, char *argv[])
Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
- 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))
- {
- // std::cout << " ---------------------------------------" << std::endl;
for (const auto& v2 : vertices(gridView_CE))
- if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
- {
- periodicIndices.unifyIndexPair({gridView_CE.indexSet().index(v1)}, {gridView_CE.indexSet().index(v2)});
- // equivCounter++;
- }
- }
- // std::cout << "equivCounter:" << equivCounter << std::endl;
- // 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(); // TODO remove
+ if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
+ {
+ periodicIndices.unifyIndexPair({gridView_CE.indexSet().index(v1)}, {gridView_CE.indexSet().index(v2)});
+ }
+ // First order periodic Lagrange-Basis
auto Basis_CE = makeBasis(
gridView_CE,
- power<dim>( //lagrange<1>(),
+ power<dim>(
Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices),
- flatLexicographic())); //
- // blockedInterleaved())); // siehe Periodicbasistest.. funktioniert auch?
- // flatInterleaved()));
+ flatLexicographic()));
- 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();
- debugLog << "phiOffset: "<< phiOffset << std::endl;
/////////////////////////////////////////////////////////
- // Data structures: Stiffness matrix and right hand side vector
+ // Data structures: Stiffness matrix and right hand side vectors
/////////////////////////////////////////////////////////
VectorCT load_alpha1, load_alpha2, load_alpha3;
MatrixCT stiffnessMatrix_CE;
@@ -2415,50 +2084,32 @@ int main(int argc, char *argv[])
substract_integralMean = true;
}
-
/////////////////////////////////////////////////////////
// Define "rhs"-functions
/////////////////////////////////////////////////////////
- Func2Tensor x3G_1 = [] (const Domain& x) {
- return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ Func2Tensor x3G_1 = [] (const Domain& x, double sign=1.0) {
+ return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; //TODO könnte hier sign übergeben?
};
Func2Tensor x3G_2 = [] (const Domain& x) {
return MatrixRT{{0.0, 0.0, 0.0}, {0.0, 1.0*x[2], 0.0}, {0.0, 0.0, 0.0}};
};
- Func2Tensor x3G_3 = [] (const Domain& x) {
+ Func2Tensor x3G_3 = [] (const Domain& x) {
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);
- };
-
+
+ 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);};
// TEST - energy method ///
// different indicatorFunction in muTerm has impact here !!
-
// double GGterm = 0.0;
// GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3G_1, x3G_1 ); // <L i(x3G_alpha) , i(x3G_beta) >
// std::cout << "GGTerm:" << GGterm << std::endl;
// std::cout << " analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
- //
-
-
- //////////////////////////////////////////////////
-
-
-
@@ -2478,12 +2129,13 @@ int main(int argc, char *argv[])
//////////////////////////////////////////////////////////////////////
// Determine global indices of components for arbitrary (local) index
//////////////////////////////////////////////////////////////////////
- int arbitraryLeafIndex = parameterSet.get<unsigned int>("arbitraryLeafIndex", 0); // localIdx..assert Number < #ShapeFcts .. also input Element?
+ int arbitraryLeafIndex = parameterSet.get<unsigned int>("arbitraryLeafIndex", 0); // localIdx..assert Number < #ShapeFcts
int arbitraryElementNumber = parameterSet.get<unsigned int>("arbitraryElementNumber", 0);
FieldVector<int,3> row;
- row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLeafIndex);
- printvector(std::cout, row, "row" , "--");
+ row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLeafIndex);
+ if(print_debug)
+ printvector(std::cout, row, "row" , "--");
/////////////////////////////////////////////////////////
// Assemble the system
@@ -2493,6 +2145,7 @@ int main(int argc, char *argv[])
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha2 ,x3G_2neg);
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3neg);
+
// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix", "--");
// printvector(std::cout, load_alpha1, "load_alpha1", "--");
@@ -2760,57 +2413,53 @@ int main(int argc, char *argv[])
for(size_t a = 0; a < 3; a++)
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) >
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp1 ,x3MatrixBasis[b]); // <L i(M_alpha) + sym(grad phi_alpha), i(x3G_beta) >
- double GGterm = 0.0;
- GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
- // TEST
- //TEST
- std::cout << " analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
+ double GGterm = 0.0;
+ GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
- std::setprecision(std::numeric_limits<float>::digits10);
+ // TEST
+ // std::setprecision(std::numeric_limits<float>::digits10);
- 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;
+ Q[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric...check positiv definitness?
- // std::cout << "Test : " << coeffContainer[a]*loadContainer[b] << std::endl; //same...
+ if (print_debug)
+ {
+ std::cout << "analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
+ std::cout << "GGTerm:" << GGterm << std::endl;
+ std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
+ }
}
printmatrix(std::cout, Q, "Matrix Q", "--");
log << "Computed Matrix Q: " << std::endl;
log << Q << std::endl;
// compute B_hat
- for(size_t a = 0; a < 3; a++)
+ for(size_t a = 0; a<3; a++)
{
- assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp2 ,B_Term); // <L i(M_alpha) + sym(grad phi_alpha), B >
-
- auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , B_Term); // <L i(x3G_alpha) , B>
-
- B_hat[a] = (coeffContainer[a]*tmp2) + GGterm;
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp2 ,B_Term); // <L i(M_alpha) + sym(grad phi_alpha), B >
+ auto GBterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , B_Term); // <L i(x3G_alpha) , B>
+ B_hat[a] = (coeffContainer[a]*tmp2) + GBterm;
-// std::cout << "check this Contribution: " << (coeffContainer[a]*tmp2) << std::endl; //see orthotropic.tex
-// std::cout <<"B_hat[i]: " << B_hat[a] << std::endl;
+ if (print_debug)
+ {
+ std::cout << "check this Contribution: " << (coeffContainer[a]*tmp2) << std::endl; //see orthotropic.tex
+ }
+
}
log << "Computed prestrain B_hat: " << std::endl;
log << B_hat << std::endl;
std::cout << "check this Contribution: " << (coeffContainer[2]*tmp2) << std::endl; //see orthotropic.tex
- // Compute effective Prestrain
+ ///////////////////////////////
+ // Compute effective Prestrain B_eff
+ //////////////////////////////
FieldVector<double, 3> Beff;
-
Q.solve(Beff,B_hat);
log << "Computed prestrain B_eff: " << std::endl;
log << Beff << std::endl;
-
-// std::cout << "Beff : " << std::endl;
-// for(size_t i=0; i<3; i++)
-// {
-// std::cout <<"Beff[i]: " << Beff[i] << std::endl;
-// }
-
+// printvector(std::cout, Beff, "Beff", "--");
@@ -2821,13 +2470,14 @@ int main(int argc, char *argv[])
std::cout<< "q1_c: " << q1_c << std::endl;
std::cout<< "q2_c: " << q2_c << std::endl;
std::cout<< "q3_c: " << q3_c << std::endl;
- std::cout<< "b1hat_c: " << B_hat[0] << std::endl;
- std::cout<< "b2hat_c: " << B_hat[1] << std::endl;
- std::cout<< "b3hat_c: " << B_hat[2] << std::endl;
- std::cout<< "b1eff_c: " << Beff[0] << std::endl;
- std::cout<< "b2eff_c: " << Beff[1] << std::endl;
- std::cout<< "b3eff_c: " << Beff[2] << std::endl;
-
+// std::cout<< "b1hat_c: " << B_hat[0] << std::endl;
+// std::cout<< "b2hat_c: " << B_hat[1] << std::endl;
+// std::cout<< "b3hat_c: " << B_hat[2] << std::endl;
+// std::cout<< "b1eff_c: " << Beff[0] << std::endl;
+// std::cout<< "b2eff_c: " << Beff[1] << std::endl;
+// std::cout<< "b3eff_c: " << Beff[2] << std::endl;
+ printvector(std::cout, B_hat, "B_hat", "--");
+ printvector(std::cout, Beff, "Beff", "--");
std::cout << "Theta: " << theta << std::endl;
std::cout << "Gamma: " << gamma << std::endl;
@@ -2842,15 +2492,13 @@ int main(int argc, char *argv[])
log << "computed b2_hat: " << B_hat[1] << std::endl;
log << "computed b3_hat: " << B_hat[2] << std::endl;
-
-
+
//////////////////////////////////////////////////////////////
// Define Analytic Solutions
//////////////////////////////////////////////////////////////
// double b1 = (mu1*p1/4)*(beta/(theta+(1-theta)*beta))*(1-theta*(1+alpha));
// double b2 = (mu1*p1/8)*(1-theta*(1+beta*alpha));
-
double b1_hat = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
double b2_hat = -(theta/4.0)*mu2;
double b3_hat = 0.0;
@@ -2859,9 +2507,6 @@ int main(int argc, char *argv[])
double b2_eff = (-3.0/2.0)*(theta*mu2)/(mu1*(1-theta)+mu2*theta);
double b3_eff = 0.0;
-
-
-
double q1 = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
@@ -2902,46 +2547,51 @@ int main(int argc, char *argv[])
return MatrixRT{{0.0 , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
-
if(write_L2Error)
{
-
- std::cout << " ----- TEST ----" << std::endl;
-
+ std::cout << " ----- TEST L2ErrorSym ----" << std::endl;
std::cout << "computeL2SymErrorNew: " << computeL2SymErrorNew(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
+ auto L2SymError = computeL2SymError(Basis_CE,phi_1,gamma,symPhi_1_analytic);
+ std::cout << "L2-SymError: " << L2SymError << std::endl;
std::cout << "computeL2SymErrorSPTest: " << computeL2SymErrorSPTest(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
std::cout << " -----------------" << std::endl;
- auto L2SymError = computeL2SymError(Basis_CE,phi_1,gamma,symPhi_1_analytic);
- std::cout << "L2-SymError: " << L2SymError << std::endl;
- // auto L2errorTest = computeL2ErrorTest(Basis_CE,phi_1,gamma,symPhi_1_analytic);
- // std::cout << "L2-ErrorTEST: " << L2errorTest << std::endl;
- auto L2Norm_Symphi = computeL2SymError(Basis_CE,phi_1,gamma,zeroFunction); // Achtung Norm SymPhi_1
+ std::cout << " ----- TEST L2NORMSym ----" << std::endl;
+ auto L2Norm_Symphi = computeL2SymError(Basis_CE,phi_1,gamma,zeroFunction);
std::cout << "L2-Norm(Symphi_1) w zerofunction: " << L2Norm_Symphi << std::endl; // TODO Not Needed
-
-
std::cout << "L2-Norm(SymPhi_1): " << computeL2SymNormCoeff(Basis_CE,phi_1,gamma) << std::endl; // does not make a difference
-
+
VectorCT zeroVec;
zeroVec.resize(Basis_CE.size());
zeroVec = 0;
-
auto L2Norm_SymAnalytic = computeL2SymError(Basis_CE,zeroVec ,gamma,symPhi_1_analytic);
std::cout << "L2-Norm(Symanalytic): " << L2Norm_SymAnalytic << std::endl;
+
+ std::cout << " -----------------" << std::endl;
- log << "L2-Error (symmetric Gradient phi_1):" << L2SymError << std::endl;
- //TEST
+ std::cout << " ----- TEST L2NORM ----" << std::endl;
auto L2Norm = computeL2NormCoeff(Basis_CE,phi_1);
std::cout << "L2Norm(phi_1) - coeff: " << L2Norm << std::endl;
// std::cout << "L2Norm(phi_1) - coeff: " << computeL2NormCoeff(Basis_CE,phi_1) << std::endl;
std::cout << "L2Norm(phi_1) - GVfunc: " << computeL2NormFunc(Basis_CE,phi_1) << std::endl;
std::cout << "L2Norm(phi_1) - coeffV2: " << computeL2NormCoeffV2(Basis_CE,phi_1) << std::endl;
std::cout << "L2Norm(phi_1) - coeffV3: " << computeL2NormCoeffV3(Basis_CE,phi_1) << std::endl;
-
std::cout << "L2ErrorOld:" << computeL2ErrorOld(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
+ std::cout << " -----------------" << std::endl;
+
+
+
+
+
+
+
+
+ log << "L2-Error (symmetric Gradient phi_1):" << L2SymError << std::endl;
+
+
double EOC = 0.0;