diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c67423e6369ea1a3dc19de60ec612ccdf74a284a
--- /dev/null
+++ b/dune/microstructure/matrix_operations.hh
@@ -0,0 +1,228 @@
+#ifndef DUNE_MICROSTRUCTURE_MATRIXOPERATIONS_HH
+#define DUNE_MICROSTRUCTURE_MATRIXOPERATIONS_HH
+
+namespace MatrixOperations {
+
+ using namespace Dune;
+
+ using MatrixRT = FieldMatrix< double, 3, 3>;
+ using VectorRT = FieldVector< double, 3>;
+
+ using std::sin;
+ using std::cos;
+
+ static MatrixRT Id (){
+ MatrixRT Id(0);
+ for (int i=0;i<3;i++)
+ Id[i][i]=1.0;
+ return Id;
+ }
+
+ static MatrixRT sym (MatrixRT M) { // 1/2 (M^T + M)
+ MatrixRT ret(0);
+ for (int i = 0; i< 3; i++)
+ for (int j = 0; j< 3; j++)
+ ret[i][j] = 0.5 * (M[i][j] + M[j][i]);
+ return ret;
+ }
+
+ static VectorRT crossProduct (VectorRT v, VectorRT w) { // v otimes w
+ return {v[1]*w[2] - v[2]*w[1], -1*(v[0]*w[2] - v[2]*w[0]), v[0]*w[1] - v[1]*w[0]};
+ }
+
+ static MatrixRT rankoneTensorproduct (VectorRT v, VectorRT w) { // v otimes w
+ return
+ {{v[0]*w[0], v[1]*w[0], v[2]*w[0]},
+ {v[0]*w[1], v[1]*w[1], v[2]*w[1]},
+ {v[0]*w[2], v[1]*w[2], v[2]*w[2]}};
+ }
+
+ static MatrixRT nematicLiquidCrystal (double p, VectorRT n){ //B = 1/6*p*Id + 1/2*p*(n otimes n)
+ MatrixRT B(0);
+ for (int i=0;i<3;i++)
+ B[i][i]=p/6.0;
+ MatrixRT n_ot_n = rankoneTensorproduct(n,n);
+ n_ot_n*=p/2.0;
+ B += n_ot_n;
+ return B;
+ }
+
+ static MatrixRT biotStrainApprox (VectorRT U, VectorRT k, VectorRT e_cs){ //E_h = (U + k x e_cs, 0, 0)
+ VectorRT k_x_ecs = crossProduct(k, e_cs);
+ VectorRT U_plus_k_x_ecs = U + k_x_ecs;
+ VectorRT e_1 = {1, 0, 0};
+ return rankoneTensorproduct(U_plus_k_x_ecs, e_1);
+ }
+
+ static MatrixRT crossSectionDirectionScaling(double w, MatrixRT M){
+ return {M[0], w*M[1], w*M[2]};
+ }
+
+ static double trace (MatrixRT M){
+ return M[0][0]+ M[1][1] + M[2][2];
+ }
+
+ static double scalarProduct (MatrixRT M1, MatrixRT M2){
+ double sum = 0;
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ sum += M1[i][j] * M2[i][j];
+ return sum;
+ }
+
+ static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2){
+ E1= sym(E1);
+ E2 = sym(E2);
+ double t1 = scalarProduct(E1,E2);
+ t1 *= 2* mu;
+ double t2 = trace(E1)*trace(E2);
+ t2 *= lambda;
+ return t1 + t2;
+ }
+
+ static MatrixRT matrixSqrt(MatrixRT M){
+ std::cout << "matrixSqrt not implemented!!!" << std::endl;//implement this
+ return M;
+ }
+
+ static double lameMu(double E, double nu){
+ return 0.5 * 1.0/(1.0 + nu) * E;
+ }
+
+ static double lameLambda(double E, double nu){
+ return nu/(1.0-2.0*nu) * 1.0/(1.0+nu) * E;
+ }
+
+ static bool isInRotatedPlane(double phi, double x1, double x2){
+ return cos(phi)*x1 + sin(phi)*x2 > 0;
+ }
+
+ static double norm(VectorRT v){
+ return sqrt(pow(v[0],2) + pow(v[1],2) + pow(v[2],2));
+ }
+
+ static double norm(MatrixRT M){
+ return sqrt(
+ pow(M[0][0],2) + pow(M[0][1],2) + pow(M[0][2],2)
+ + pow(M[1][0],2) + pow(M[1][1],2) + pow(M[1][2],2)
+ + pow(M[2][0],2) + pow(M[2][1],2) + pow(M[2][2],2));
+ }
+
+ /*
+ template<double phi>
+ static bool isInRotatedPlane(double x1, double x2){
+ return cos(phi)*x1 + sin(phi)*x2 > 0;
+ }*/
+
+}
+
+
+
+// OPTIONAL H1-Seminorm ...
+
+/*
+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;
+}*/
+
+
+
+
+
+
+/*
+ class BiotStrainApprox
+ {
+
+ // E_h = h^{-1} (R^T F_h - Id)
+ //
+ //E_h = (U + K e_cs, 0, 0)
+ //
+
+ Func2Tensor E_a = [] (const Domain& z) //extra Klasse
+ {
+ return biotStrainApprox({1, 0, 0}, {0, 0, 0}, {0, z[1], z[2]}); //MatrixRT{{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ };
+
+ // twist in x2-x3 plane E_K1 = sym (e1 x (0,x2,x3)^T ot e1)
+ Func2Tensor E_K1 = [] (const Domain& z)
+ {
+ return biotStrainApprox({0, 0, 0}, {1, 0, 0}, {0, z[1], z[2]}); //MatrixRT{{0.0,-0.5*z[2], 0.5*z[1]}, {-0.5*z[2], 0.0 , 0.0 }, {0.5*z[1],0.0,0.0}};
+ };
+
+ // bend strain in x1-x2 plane E_K2 = sym (e2 x (0,x2,x3)^T ot e1)
+ Func2Tensor E_K2 = [] (const Domain& z)
+ {
+ return biotStrainApprox({0, 0, 0}, {0, 1, 0}, {0, z[1], z[2]}); //MatrixRT{{z[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0}};
+ };
+
+ // bend strain in x1-x3 plane E_K3 = sym (e3 x (0,x2,x3)^T ot e1)
+ Func2Tensor E_K3 = [] (const Domain& z)
+ {
+ return biotStrainApprox({0, 0, 0}, {0, 0, 1}, {0, z[1], z[2]}); //MatrixRT{{-z[1], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ };
+ };
+*/
+
+#endif
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index c055cd404df3ba28f18d4fc9a1e4662990530415..4d9563159b37e74326569d99605c3fd7c9706d4e 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -46,6 +46,7 @@
#include <dune/microstructure/prestrainimp.hh>
+#include <dune/microstructure/matrix_operations.hh>
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
// #include <dune/fufem/discretizationerror.hh>
@@ -58,6 +59,7 @@
#include <iomanip>
using namespace Dune;
+using namespace MatrixOperations;
//////////////////////////////////////////////////////////////////////
@@ -118,9 +120,9 @@ auto arbitraryComponentwiseIndices(const Basis& basis,
for (int k = 0; k < 3; k++)
{
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;
+ // std::cout << "rowLocal:" << rowLocal << std::endl;
+ // std::cout << "row[k]:" << row[k] << std::endl;
}
}
elementCounter++;
@@ -142,7 +144,6 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// Local StiffnessMatrix of the form:
// | phi*phi m*phi |
// | phi *m m*m |
-
using Element = typename LocalView::Element;
const Element element = localView.element();
auto geometry = element.geometry();
@@ -188,11 +189,11 @@ void computeElementStiffnessMatrix(const LocalView& localView,
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 < dimWorld; k++)
- for (size_t l=0; l< dimWorld; l++)
+ for (size_t i=0; i < nSf; i++)
{
- for (size_t i=0; i < nSf; i++)
+ // "phi*phi"-part
+ for (size_t l=0; l< dimWorld; l++)
for (size_t j=0; j < nSf; j++ )
{
// (scaled) Deformation gradient of the ansatz basis function
@@ -206,111 +207,29 @@ void computeElementStiffnessMatrix(const LocalView& localView,
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<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]);
- //printmatrix(std::cout, strainU , "strainU", "--");
- }
-
- // St.Venant-Kirchhoff stress
- // < L sym[D_gamma*nabla phi_i], sym[D_gamma*nabla phi_j] >
- // stressU*strainV
- MatrixRT stressU(0);
- stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
-
- double trace = 0.0;
- for (int ii=0; ii<dimWorld; ii++)
- trace += strainU[ii][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<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- energyDensity += stressU[ii][jj] * strainV[ii][jj];
+
+ double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
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< 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
- defGradientV[l][1] = gradients[j][1]; //X2
- defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
-
- // symmetric Gradient (Elastic Strains)
- 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]);
- }
-
- for( size_t m=0; m<3; m++)
- {
- // < L G_i, sym[D_gamma*nabla phi_j] >
- // L G_i* strainV
-
- // St.Venant-Kirchhoff stress
- MatrixRT stressG(0);
- stressG.axpy(2*mu(quadPos), basisContainer[m]); //this += 2mu *strainU
+
+ // "m*phi" & "phi*m" - part
+ for( size_t m=0; m<3; m++)
+ {
+ double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
- double traceG = 0.0;
- for (int ii=0; ii<dimWorld; ii++)
- {
- traceG += basisContainer[m][ii][ii];
+ auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
+ elementMatrix[row][localPhiOffset+m] += value;
+ elementMatrix[localPhiOffset+m][row] += value;
+ }
}
-
- for (int ii=0; ii<dimWorld; ii++)
- stressG[ii][ii] += lambda(quadPos) * traceG;
-
- double energyDensityGphi = 0.0;
- 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<dimWorld; ii++)
- traceG += basisContainer[m][ii][ii];
-
- for (int ii=0; ii<dimWorld; ii++)
- stressG[ii][ii] += lambda(quadPos) * traceG;
-
- double energyDensityGG = 0.0;
- for (int ii=0; ii<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- energyDensityGG += stressG[ii][jj] * basisContainer[n][ii][jj];
-
+ double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
elementMatrix[localPhiOffset+m][localPhiOffset+n]= energyDensityGG * quadPoint.weight() * integrationElement;
}
}
@@ -398,7 +317,6 @@ void computeElementLoadVector( const LocalView& localView,
// | f*phi|
// | --- |
// | f*m |
-
using Element = typename LocalView::Element;
const auto element = localView.element();
const auto geometry = element.geometry();
@@ -453,29 +371,7 @@ void computeElementLoadVector( const LocalView& localView,
defGradientV[k][1] = gradients[i][1]; // X2
defGradientV[k][2] = (1.0/gamma)*gradients[i][2]; // X3
- // Elastic Strain
- 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
-
- double trace = 0.0;
- for (int ii=0; ii<dimWorld; ii++)
- trace += strainV[ii][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<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- energyDensity += stressV[ii][jj] *forceTerm(quadPos)[ii][jj];
-
+ double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientV, forceTerm(quadPos));
size_t row = localView.tree().child(k).localIndex(i);
elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
}
@@ -483,22 +379,7 @@ void computeElementLoadVector( const LocalView& localView,
// "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<dimWorld; ii++)
- traceG += basisContainer[m][ii][ii];
-
- for (int ii=0; ii<dimWorld; ii++)
- stressG[ii][ii] += lambda(quadPos) * traceG;
-
- double energyDensityfG = 0;
- for (int ii=0; ii<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- energyDensityfG += stressG[ii][jj] * forceTerm(quadPos)[ii][jj];
-
+ double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], forceTerm(quadPos));
elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
}
}
@@ -515,7 +396,7 @@ void assembleCellStiffness(const Basis& basis,
ParameterSet& parameterSet
)
{
- std::cout << "assemble stiffnessmatrix begins." << std::endl;
+ std::cout << "assemble Stiffness-Matrix begins." << std::endl;
MatrixIndexSet occupationPattern;
getOccupationPattern(basis, occupationPattern, parameterSet);
@@ -530,8 +411,8 @@ void assembleCellStiffness(const Basis& basis,
localView.bind(element);
muLocal.bind(element);
lambdaLocal.bind(element);
-
const int localPhiOffset = localView.size();
+
//std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
@@ -543,32 +424,25 @@ void assembleCellStiffness(const Basis& basis,
// GLOBAL STIFFNES ASSEMBLY
//////////////////////////////////////////////////////////////////////////////
for (size_t i=0; i<localPhiOffset; i++)
+ for (size_t j=0; j<localPhiOffset; j++ )
{
- auto row = localView.index(i);
-
- for (size_t j=0; j<localPhiOffset; j++ )
- {
+ auto row = localView.index(i);
auto col = localView.index(j);
-
matrix[row][col] += elementMatrix[i][j];
- }
}
for (size_t i=0; i<localPhiOffset; i++)
- for(size_t m=0; m<3; m++)
- {
+ for(size_t m=0; m<3; m++)
+ {
auto row = localView.index(i);
-
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++ )
- {
+ for (size_t n=0; n<3; n++ )
matrix[phiOffset+m][phiOffset+n] += elementMatrix[localPhiOffset+m][localPhiOffset+n];
- }
-// printmatrix(std::cout, matrix, "StiffnessMatrix", "--");
- }
+ // printmatrix(std::cout, matrix, "StiffnessMatrix", "--");
+ }
}
@@ -590,7 +464,7 @@ void assembleCellLoad(const Basis& basis,
// Transform G_alpha's to GridViewFunctions/LocalFunctions
auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
- auto loadFunctional = localFunction(loadGVF); // Dune::Functions:: notwendig?
+ auto loadFunctional = localFunction(loadGVF);
for (const auto& element : elements(basis.gridView()))
{
@@ -602,24 +476,21 @@ 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", "--");
// printvector(std::cout, elementRhs, "elementRhs", "--");
-
- // LOAD ASSEMBLY
+ //////////////////////////////////////////////////////////////////////////////
+ // GLOBAL LOAD ASSEMBLY
+ //////////////////////////////////////////////////////////////////////////////
for (size_t p=0; p<localPhiOffset; p++)
{
- // 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", "--");
+ //printvector(std::cout, b, "b", "--");
}
}
@@ -627,8 +498,8 @@ void assembleCellLoad(const Basis& basis,
template<class Basis, class LocalFunction1, class LocalFunction2, class MatrixFunction>
auto energy(const Basis& basis,
- LocalFunction1& muLocal,
- LocalFunction2& lambdaLocal,
+ LocalFunction1& mu,
+ LocalFunction2& lambda,
const MatrixFunction& matrixFieldFuncA,
const MatrixFunction& matrixFieldFuncB)
{
@@ -636,7 +507,6 @@ auto energy(const Basis& basis,
// 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 dimWorld = 3;
@@ -651,7 +521,6 @@ auto energy(const Basis& basis,
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
-
// TEST
// FieldVector<double,3> testvector = {1.0 , 1.0 , 1.0};
// printmatrix(std::cout, matrixFieldFuncB(testvector) , "matrixFieldB(testvector) ", "--");
@@ -661,44 +530,28 @@ auto energy(const Basis& basis,
localView.bind(e);
matrixFieldA.bind(e);
matrixFieldB.bind(e);
- muLocal.bind(e);
- lambdaLocal.bind(e);
+ mu.bind(e);
+ lambda.bind(e);
double elementEnergy = 0.0;
-// double elementEnergy_HP = 0.0;
+ //double elementEnergy_HP = 0.0;
auto geometry = e.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
+ //int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
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)
{
-
const auto& quadPos = quadPoint.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
-
- double trace = 0.0;
- for (int ii=0; ii<dimWorld; ii++)
- trace += strain1[ii][ii];
-
- for (int ii=0; ii<dimWorld; ii++)
- stressU[ii][ii] += lambdaLocal(quadPos) * trace;
-
- // Local energy density: stress times strain
- double energyDensity = 0.0;
- for (int ii=0; ii<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- energyDensity += stressU[ii][jj] * strain2[ii][jj];
+ double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), strain1, strain2);
elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
//elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
@@ -723,14 +576,12 @@ void setOneBaseFunctionToZero(const Basis& basis,
)
{
std::cout << "Setting one Basis-function to zero" << std::endl;
-
constexpr int dim = 3;
- 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);
+ FieldVector<int,3> row = arbitraryComponentwiseIndices(basis,arbitraryElementNumber,arbitraryLeafIndex);
for (int k = 0; k < dim; k++)
{
@@ -753,7 +604,6 @@ auto childToIndexMap(const Basis& basis,
{
// Input : child/component
// Output : determine all Indices that belong to that component
-
auto localView = basis.localView();
std::vector<int> r = { };
@@ -862,17 +712,13 @@ auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>
};
-
////////////////////////////////////////////////////////////// L2-ERROR /////////////////////////////////////////////////////////////////
-
-
template<class Basis, class Vector, class MatrixFunction>
double computeL2SymError(const Basis& basis,
- Vector& coeffVector,
- const double gamma,
- const MatrixFunction& matrixFieldFunc
- )
+ Vector& coeffVector,
+ const double gamma,
+ const MatrixFunction& matrixFieldFunc)
{
double error = 0.0;
constexpr int dim = 3;
@@ -921,22 +767,11 @@ double computeL2SymError(const Basis& basis,
// (scaled) Deformation gradient of the ansatz basis function
MatrixRT defGradientU(0);
- defGradientU[k][0] = coeffVector[globalIdx1]*gradients[i][0]; // Y //hier i:leaf oder localIdx?
+ defGradientU[k][0] = coeffVector[globalIdx1]*gradients[i][0]; // Y //hier i:leafIdx
defGradientU[k][1] = coeffVector[globalIdx1]*gradients[i][1]; //X2
defGradientU[k][2] = coeffVector[globalIdx1]*(1.0/gamma)*gradients[i][2]; //X3
- // symmetric Gradient (Elastic Strains)
-// MatrixRT strainU(0);
- 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]); // not needed (old version)
- tmp[ii][jj] += 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
- }
-// tmp += strainU;
-// printmatrix(std::cout, strainU, "strainU", "--");
-// printmatrix(std::cout, tmp, "tmp", "--");
-//
+ tmp += sym(defGradientU);
}
// printmatrix(std::cout, matrixField(quadPos), "matrixField(quadPos)", "--");
// printmatrix(std::cout, tmp, "tmp", "--"); // TEST Symphi_1 hat andere Struktur als analytic?
@@ -954,14 +789,10 @@ double computeL2SymError(const Basis& basis,
}
return sqrt(error);
}
-
-
-
////////////////////////////////////////////////////////////// L2-NORM /////////////////////////////////////////////////////////////////
template<class Basis, class Vector>
double computeL2Norm(const Basis& basis,
- Vector& coeffVector
- )
+ Vector& coeffVector)
{
// IMPLEMENTATION with makeDiscreteGlobalBasisFunction
double error = 0.0;
@@ -984,87 +815,11 @@ double computeL2Norm(const Basis& basis,
{
const auto& quadPos = quadPoint.position();
const double integrationElement = element.geometry().integrationElement(quadPos);
- double tmp = localfun(quadPos) * localfun(quadPos);
- error += tmp * quadPoint.weight() * integrationElement;
+ error += localfun(quadPos)*localfun(quadPos) * quadPoint.weight() * integrationElement;
}
}
return sqrt(error);
}
-
-
-
-// OPTIONAL H1-Seminorm ...
-
-/*
-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;
-}*/
-
-
-
-
-
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -1072,10 +827,8 @@ for (const auto& e : elements(basis_.gridView()))
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
int main(int argc, char *argv[])
{
-
MPIHelper::instance(argc, argv);
ParameterTree parameterSet;
@@ -1212,7 +965,6 @@ int main(int argc, char *argv[])
auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
auto lambdaLocal = localFunction(lambdaGridF);
-
///////////////////////////////////
// --- Choose Solver ---
// 1 : CG-Solver
@@ -1221,7 +973,6 @@ int main(int argc, char *argv[])
///////////////////////////////////
unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 1);
-
// Print Options
bool print_debug = parameterSet.get<bool>("print_debug", false);
bool write_corrector_phi1 = parameterSet.get<bool>("write_corrector_phi1", false);
@@ -1230,7 +981,6 @@ int main(int argc, char *argv[])
bool write_L2Error = parameterSet.get<bool>("write_L2Error", false);
bool write_IntegralMean = parameterSet.get<bool>("write_IntegralMean", false);
-
/////////////////////////////////////////////////////////
// Choose a finite element space for Cell Problem
/////////////////////////////////////////////////////////
@@ -1342,7 +1092,6 @@ int main(int argc, char *argv[])
// printvector(std::cout, load_alpha1, "load_alpha1 after setOneBaseFunctionToZero", "--");
}
-
////////////////////////////////////////////////////
// Compute solution
////////////////////////////////////////////////////
@@ -1551,7 +1300,6 @@ int main(int argc, char *argv[])
Basis_CE,
phi_3);
-
/////////////////////////////////////////////////////////
// Create Containers for Basis-Matrices, Correctors and Coefficients
/////////////////////////////////////////////////////////
@@ -1708,7 +1456,6 @@ int main(int argc, char *argv[])
std::cout << "L2SymError: " << L2SymError << std::endl;
std::cout << " -----------------" << std::endl;
-
std::cout << " ----- L2NORMSym ----" << std::endl;
auto L2Norm_Symphi = computeL2SymError(Basis_CE,phi_1,gamma,zeroFunction);
std::cout << "L2-Norm(Symphi_1): " << L2Norm_Symphi << std::endl;
@@ -1735,10 +1482,8 @@ int main(int argc, char *argv[])
if(levelCounter > 0)
{
// Storage_Error:: #1 level #2 L2SymError #3 L2SymErrorOrder #4 L2Norm(sym) #5 L2Norm(sym-analytic) #6 L2Norm(phi_1)
-
// std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter-1)*6)+1 << std::endl; // Besser std::map ???
// std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter)*6)+1 << std::endl; // für Storage_Error[idx] muss idx zur compile time feststehen?!
-
auto ErrorOld = std::get<double>(Storage_Error[((levelCounter-1)*6)+1]);
auto ErrorNew = std::get<double>(Storage_Error[(levelCounter*6)+1]);
//
@@ -1753,7 +1498,6 @@ int main(int argc, char *argv[])
Storage_Error.push_back(L2Norm_SymAnalytic);
Storage_Error.push_back(L2Norm);
}
-
//////////////////////////////////////////////////////////////////////////////////////////////
// Write Data to Matlab / Optimization-Code
//////////////////////////////////////////////////////////////////////////////////////////////
@@ -1769,8 +1513,6 @@ int main(int argc, char *argv[])
BeffMat[0] = Beff;
writeMatrixToMatlab(BeffMat, "../../Matlab-Programs/BMatrix.txt");
-
-
//////////////////////////////////////////////////////////////////////////////////////////////
// Write result to VTK file
//////////////////////////////////////////////////////////////////////////////////////////////