diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index c49ad86676e9aa659e6daa42627ca4c427c1d03e..2dc0af1c0a6048ed3686268669d09cb5c27e491d 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -1,7 +1,7 @@
#include <config.h>
#include <array>
#include <vector>
-
+#include <fstream>
#include <iostream>
@@ -23,8 +23,8 @@
#include <dune/istl/multitypeblockmatrix.hh>
#include <dune/istl/multitypeblockvector.hh>
#include <dune/istl/matrixindexset.hh>
-#include<dune/istl/solvers.hh>
-#include<dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
#include <dune/istl/io.hh>
@@ -64,279 +64,279 @@ using namespace Dune;
template<class LocalView, class Matrix, class localFunction1, class localFunction2>
void computeElementStiffnessMatrix(const LocalView& localView,
- Matrix& elementMatrix,
- const localFunction1& mu,
- const localFunction2& lambda,
- const double gamma
- )
+ Matrix& elementMatrix,
+ const localFunction1& mu,
+ const localFunction2& lambda,
+ const double gamma
+ )
{
- // Local StiffnessMatrix of the form:
- // | phi*phi m*phi |
- // | phi *m m*m |
+ // 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();
+ using Element = typename LocalView::Element;
+ const Element element = localView.element();
+ auto geometry = element.geometry();
- constexpr int dim = Element::dimension;
- constexpr int nCompo = dim;
+ 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&) >;
+ 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&) >;
- elementMatrix.setSize(localView.size()+3, localView.size()+3); //extend by dim ´R_sym^{2x2}
- elementMatrix = 0;
+ elementMatrix.setSize(localView.size()+3, localView.size()+3); //extend by dim ´R_sym^{2x2}
+ elementMatrix = 0;
- // LocalBasis-Offset
- const int localPhiOffset = localView.size();
+ // LocalBasis-Offset
+ const int localPhiOffset = localView.size();
- const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
- const auto nSf = localFiniteElement.localBasis().size();
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
+ const auto nSf = localFiniteElement.localBasis().size();
- ///////////////////////////////////////////////
- // Basis for R_sym^{2x2} // wird nicht als Funktion benötigt da konstant...
- //////////////////////////////////////////////
- MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
- MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
- MatrixRT G_3 {{0.0, 1/sqrt(2), 0.0}, {1/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", "--");
- ////////////////////////////////////////////////////
+ ///////////////////////////////////////////////
+ // Basis for R_sym^{2x2} // wird nicht als Funktion benötigt da konstant...
+ //////////////////////////////////////////////
+ MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+ MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+ MatrixRT G_3 {{0.0, 1/sqrt(2), 0.0}, {1/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);
- const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
- for (const auto& quadPoint : quad)
- {
+ 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; // old
+ 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());
- const auto& quadPos = quadPoint.position();
- const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
- const auto integrationElement = geometry.integrationElement(quadPoint.position());
+ for (size_t i=0; i<gradients.size(); i++)
+ jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
-// std::vector<FieldMatrix<double,1,dim> > referenceGradients; // old
- 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());
+ // "phi*phi"-part
+ for (size_t k=0; k < nCompo; k++)
+ for (size_t l=0; l< nCompo; 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
- for (size_t i=0; i<gradients.size(); i++)
- jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+ // (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
- // "phi*phi"-part
- for (size_t k=0; k < nCompo; k++)
- for (size_t l=0; l< nCompo; l++)
- {
- for (size_t i=0; i < nSf; i++)
- for (size_t j=0; j < nSf; j++ )
- {
+ // symmetric Gradient (Elastic Strains)
+ MatrixRT strainU, strainV;
+ 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]); // 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
- // (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
+ }
+
+ // 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
- // (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
+ double trace = 0;
+ for (int ii=0; ii<nCompo; ii++)
+ trace += strainU[ii][ii];
+ for (int ii=0; ii<nCompo; ii++)
+ stressU[ii][ii] += lambda(quadPos) * trace;
- // symmetric Gradient (Elastic Strains)
- MatrixRT strainU, strainV;
- 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]); // 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
+ // 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);
- // 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
+ elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
+
+ }
+
+ }
+
+
+ // "m*phi" & "phi*m" -part
+ for (size_t l=0; l< nCompo; 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<nCompo; ii++)
+ for (int jj=0; jj<nCompo; 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
+
+ double traceG = 0;
+ for (int ii=0; ii<nCompo; ii++)
+ traceG += basisContainer[m][ii][ii];
+
+ for (int ii=0; ii<nCompo; ii++)
+ stressG[ii][ii] += lambda(quadPos) * traceG;
+
+ double energyDensityGphi = 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
+
+ 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;
+
+ }
+
+ }
+
+
+
+ // "phi*m"-part
+ // for (size_t k=0; k < nCompo; k++)
+ // for (size_t i=0; i < nSf; i++)
+ // {
+ //
+ // // (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
+ //
+ //
+ // // symmetric Gradient (Elastic Strains)
+ // MatrixRT strainU;
+ // 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]); // symmetric gradient
+ // }
+ //
+ // // St.Venant-Kirchhoff stress
+ // MatrixRT stressU(0);
+ // stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
+ //
+ // double trace = 0;
+ // for (int ii=0; ii<nCompo; ii++)
+ // trace += strainU[ii][ii];
+ //
+ // for (int ii=0; ii<nCompo; ii++)
+ // stressU[ii][ii] += lambda(quadPos) * trace;
+ //
+ // for( size_t n=0; n<3; n++)
+ // {
+ //
+ // // < L sym[D_gamma*nabla phi_i], G_j >
+ // double energyDensityPhiG = 0;
+ // for (int ii=0; ii<nCompo; ii++)
+ // for (int jj=0; jj<nCompo; jj++)
+ // energyDensityPhiG += stressU[ii][jj] * basisContainer[n][ii][jj]; // "phi*m"-part
+ //
+ // size_t col = localView.tree().child(k).localIndex(i);
+ //
+ // elementMatrix[localPhiOffset+n][col] += energyDensityPhiG * quadPoint.weight() * integrationElement;
+ //
+ // }
+ // }
- double trace = 0;
- for (int ii=0; ii<nCompo; ii++)
- trace += strainU[ii][ii];
-
- for (int ii=0; ii<nCompo; 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);
-
- elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
-
- }
-
- }
-
-
- // "m*phi" & "phi*m" -part
- for (size_t l=0; l< nCompo; 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<nCompo; ii++)
- for (int jj=0; jj<nCompo; 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
-
- double traceG = 0;
- for (int ii=0; ii<nCompo; ii++)
- traceG += basisContainer[m][ii][ii];
-
- for (int ii=0; ii<nCompo; ii++)
- stressG[ii][ii] += lambda(quadPos) * traceG;
-
- double energyDensityGphi = 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
-
- 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;
-
- }
-
- }
-
-
-
- // "phi*m"-part
-// for (size_t k=0; k < nCompo; k++)
-// for (size_t i=0; i < nSf; i++)
-// {
-//
-// // (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
-//
-//
-// // symmetric Gradient (Elastic Strains)
-// MatrixRT strainU;
-// 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]); // symmetric gradient
-// }
-//
-// // St.Venant-Kirchhoff stress
-// MatrixRT stressU(0);
-// stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
-//
-// double trace = 0;
-// for (int ii=0; ii<nCompo; ii++)
-// trace += strainU[ii][ii];
-//
-// for (int ii=0; ii<nCompo; ii++)
-// stressU[ii][ii] += lambda(quadPos) * trace;
-//
-// for( size_t n=0; n<3; n++)
-// {
-//
-// // < L sym[D_gamma*nabla phi_i], G_j >
-// double energyDensityPhiG = 0;
-// for (int ii=0; ii<nCompo; ii++)
-// for (int jj=0; jj<nCompo; jj++)
-// energyDensityPhiG += stressU[ii][jj] * basisContainer[n][ii][jj]; // "phi*m"-part
-//
-// size_t col = localView.tree().child(k).localIndex(i);
-//
-// elementMatrix[localPhiOffset+n][col] += energyDensityPhiG * quadPoint.weight() * integrationElement;
-//
-// }
-// }
-
-
-
-
- // "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;
- for (int ii=0; ii<nCompo; ii++)
- traceG += basisContainer[m][ii][ii];
-
- for (int ii=0; ii<nCompo; ii++)
- stressG[ii][ii] += lambda(quadPos) * traceG;
-
- double energyDensityGG = 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
-
- elementMatrix[localPhiOffset+m][localPhiOffset+n]= energyDensityGG * quadPoint.weight() * integrationElement;
-
- }
- }
+
+
+ // "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;
+ for (int ii=0; ii<nCompo; ii++)
+ traceG += basisContainer[m][ii][ii];
+
+ for (int ii=0; ii<nCompo; ii++)
+ stressG[ii][ii] += lambda(quadPos) * traceG;
+
+ double energyDensityGG = 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
+
+ elementMatrix[localPhiOffset+m][localPhiOffset+n]= energyDensityGG * quadPoint.weight() * integrationElement;
+
+ }
+
+ }
}
@@ -348,106 +348,106 @@ void computeElementStiffnessMatrix(const LocalView& localView,
template<class Basis>
void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb)
{
- // OccupationPattern:
- // | phi*phi m*phi |
- // | phi *m m*m |
- auto gridView = basis.gridView();
+ // OccupationPattern:
+ // | phi*phi m*phi |
+ // | phi *m m*m |
+ auto gridView = basis.gridView();
- auto localView = basis.localView();
+ auto localView = basis.localView();
- const int phiOffset = basis.dimension();
+ const int phiOffset = basis.dimension();
- nb.resize(basis.size()+3, basis.size()+3); //
+ nb.resize(basis.size()+3, basis.size()+3); //
- for (const auto& element : elements(gridView))
+ for (const auto& element : elements(gridView))
+ {
+ localView.bind(element);
+ for (size_t i=0; i<localView.size(); i++)
{
- localView.bind(element);
- for (size_t i=0; i<localView.size(); i++)
- {
- // The global index of the i-th vertex of the element
- auto row = localView.index(i);
- for (size_t j=0; j<localView.size(); j++ )
- {
- // The global index of the j-th vertex of the element
- auto col = localView.index(j);
-
- nb.add(row[0],col[0]); // nun würde auch nb.add(row,col) gehen..
- }
- }
-
+ // The global index of the i-th vertex of the element
+ auto row = localView.index(i);
+ for (size_t j=0; j<localView.size(); j++ )
+ {
+ // The global index of the j-th vertex of the element
+ auto col = localView.index(j);
- for (size_t i=0; i<localView.size(); i++)
- {
- auto row = localView.index(i);
+ nb.add(row[0],col[0]); // nun würde auch nb.add(row,col) gehen..
+ }
+ }
- for (size_t j=0; j<3; j++ )
- {
- nb.add(row,phiOffset+j); // m*phi part of StiffnessMatrix
- nb.add(phiOffset+j,row); // phi*m part of StiffnessMatrix
- }
- }
+ for (size_t i=0; i<localView.size(); i++)
+ {
+ auto row = localView.index(i);
- for (size_t i=0; i<3; i++ )
- for (size_t j=0; j<3; j++ )
- {
- nb.add(phiOffset+i,phiOffset+j); // m*m part of StiffnessMatrix
- }
+ for (size_t j=0; j<3; j++ )
+ {
+ nb.add(row,phiOffset+j); // m*phi part of StiffnessMatrix
+ nb.add(phiOffset+j,row); // phi*m part of StiffnessMatrix
+ }
}
- //////////////////////////////////////////////////////////////////
- // setIntegralZero:
- int arbitraryIndex = 0;
+ for (size_t i=0; i<3; i++ )
+ for (size_t j=0; j<3; j++ )
+ {
+ nb.add(phiOffset+i,phiOffset+j); // m*m part of StiffnessMatrix
+ }
- FieldVector<int,3> row;
- unsigned int elementCounter = 1;
+ }
- for (const auto& element : elements(basis.gridView()))
- {
- localView.bind(element);
+ //////////////////////////////////////////////////////////////////
+ // setIntegralZero:
+ int arbitraryIndex = 0;
- if(elementCounter == 1) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
- {
- for (int k = 0; k < 3; k++)
- {
- auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
- row[k] = localView.index(rowLocal);
- }
- }
+ FieldVector<int,3> row;
+ unsigned int elementCounter = 1;
- // "global assembly"
- for (int k = 0; k<3; k++)
- for (size_t j=0; j<localView.size(); j++)
- {
- auto col = localView.index(j);
- nb.add(row[k],col);
- }
- elementCounter++;
+ for (const auto& element : elements(basis.gridView()))
+ {
+ localView.bind(element);
+
+ if(elementCounter == 1) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
+ {
+ for (int k = 0; k < 3; k++)
+ {
+ auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
+ row[k] = localView.index(rowLocal);
+ }
}
- /////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////
- // setOneBaseFunctionToZero necessary?
+ // "global assembly"
+ for (int k = 0; k<3; k++)
+ for (size_t j=0; j<localView.size(); j++)
+ {
+ auto col = localView.index(j);
+ nb.add(row[k],col);
+ }
+ elementCounter++;
+ }
+ /////////////////////////////////////////////////////////////////
+
+ ////////////////////////////////////////////////////////////////
+ // setOneBaseFunctionToZero necessary?
- //Determine 3 global indices (one for each component of an arbitrary FE-function).. oder mittels Offset?
- // use first element:
-// auto it = basis.gridView().template begin<0>();
-// localView.bind(*it);
-//
-// int arbitraryIndex = 0;
-// FieldVector<int,3> row; //fill with Indices..
-//
-//
-// for (int k = 0; k < 3; k++)
-// {
-// auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
-// row[k] = localView.index(rowLocal);
-// auto cIt = stiffnessMatrix[row[k]].begin();
-// auto cEndIt = stiffnessMatrix[row[k]].end();
-// for (; cIt!=cEndIt; ++cIt)
-// *cIt = (cIt.index()==row[k]) ? 1.0 : 0.0;
-// }
+ //Determine 3 global indices (one for each component of an arbitrary FE-function).. oder mittels Offset?
+ // use first element:
+ // auto it = basis.gridView().template begin<0>();
+ // localView.bind(*it);
+ //
+ // int arbitraryIndex = 0;
+ // FieldVector<int,3> row; //fill with Indices..
+ //
+ //
+ // for (int k = 0; k < 3; k++)
+ // {
+ // auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
+ // row[k] = localView.index(rowLocal);
+ // auto cIt = stiffnessMatrix[row[k]].begin();
+ // auto cEndIt = stiffnessMatrix[row[k]].end();
+ // for (; cIt!=cEndIt; ++cIt)
+ // *cIt = (cIt.index()==row[k]) ? 1.0 : 0.0;
+ // }
@@ -462,129 +462,129 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb)
// < L sym[D_gamma*nabla phi_i], x_3G_alpha >
template<class LocalView, class LocalFunction1, class LocalFunction2, class Vector, class Force>
void computeElementLoadVector( const LocalView& localView,
- LocalFunction1& mu,
- LocalFunction2& lambda,
+ LocalFunction1& mu,
+ LocalFunction2& lambda,
const double gamma,
- Vector& elementRhs,
- const Force& forceTerm
+ Vector& elementRhs,
+ const Force& forceTerm
)
{
- // | f*phi|
- // | --- |
- // | f*m |
+ // | f*phi|
+ // | --- |
+ // | f*m |
- using Element = typename LocalView::Element;
- const auto element = localView.element();
- const auto geometry = element.geometry();
- constexpr int dim = Element::dimension;
- constexpr int nCompo = dim;
+ using Element = typename LocalView::Element;
+ const auto element = localView.element();
+ const auto geometry = element.geometry();
+ constexpr int dim = Element::dimension;
+ constexpr int nCompo = dim;
- using VectorRT = FieldVector< double, nCompo>;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+ using VectorRT = FieldVector< double, nCompo>;
+ using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
- // 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();
+ // 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;
+ elementRhs.resize(localView.size() +3 );
+ elementRhs = 0;
- // LocalBasis-Offset
- const int localPhiOffset = localView.size();
+ // LocalBasis-Offset
+ const int localPhiOffset = localView.size();
- ///////////////////////////////////////////////
- // Basis for R_sym^{2x2}
- //////////////////////////////////////////////
- MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
- MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
- MatrixRT G_3 {{0.0, 1/sqrt(2), 0.0}, {1/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
- std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
+ ///////////////////////////////////////////////
+ // Basis for R_sym^{2x2}
+ //////////////////////////////////////////////
+ MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+ MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+ MatrixRT G_3 {{0.0, 1/sqrt(2), 0.0}, {1/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1); // für simplex
- const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+ 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){
+ for (const auto& quadPoint : quad) {
- const FieldVector<double,dim>& quadPos = quadPoint.position();
- const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
- const double integrationElement = geometry.integrationElement(quadPos);
+ const FieldVector<double,dim>& quadPos = quadPoint.position();
+ const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+ const double integrationElement = geometry.integrationElement(quadPos);
- std::vector<FieldMatrix<double,1,nCompo> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+ std::vector<FieldMatrix<double,1,nCompo> > referenceGradients;
+ localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
- std::vector< VectorRT > gradients(referenceGradients.size());
- for (size_t i=0; i< gradients.size(); i++)
- jacobian.mv(referenceGradients[i][0], gradients[i]);
+ std::vector< VectorRT > 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++)
- {
- // Deformation gradient of the ansatz basis function
- MatrixRT defGradientV(0);
- defGradientV[k][0] = gradients[i][0]; // Y
- defGradientV[k][1] = gradients[i][1]; // X2
- defGradientV[k][2] = (1.0/gamma)*gradients[i][2]; // X3
+ // "f*phi"-part
+ for (size_t i=0; i < nSf; i++)
+ for (size_t k=0; k < nCompo; k++)
+ {
+ // Deformation gradient of the ansatz basis function
+ MatrixRT defGradientV(0);
+ defGradientV[k][0] = gradients[i][0]; // Y
+ defGradientV[k][1] = gradients[i][1]; // X2
+ 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++)
- strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
+ // Elastic Strain
+ MatrixRT strainV; //strainV never used?
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; 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
+ // St.Venant-Kirchhoff stress
+ MatrixRT stressV(0);
+ stressV.axpy(2*mu(quadPos), strainV); //this += 2mu *strainU
- double trace = 0;
- for (int ii=0; ii<nCompo; ii++)
- trace += strainV[ii][ii];
+ double trace = 0;
+ for (int ii=0; ii<nCompo; ii++)
+ trace += strainV[ii][ii];
- for (int ii=0; ii<nCompo; ii++)
- stressV[ii][ii] += lambda(quadPos) * trace;
+ for (int ii=0; ii<nCompo; ii++)
+ stressV[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 += stressV[ii][jj] * forceTerm(quadPos)[ii][jj];
+ // Local energy density: stress times strain
+ double energyDensity = 0;
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ energyDensity += stressV[ii][jj] * forceTerm(quadPos)[ii][jj];
- size_t row = localView.tree().child(k).localIndex(i);
- elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
+ 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++)
- {
+ // "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
+ // 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++)
- traceG += basisContainer[m][ii][ii];
+ double traceG = 0;
+ for (int ii=0; ii<nCompo; ii++)
+ traceG += basisContainer[m][ii][ii];
- for (int ii=0; ii<nCompo; ii++)
- stressG[ii][ii] += lambda(quadPos) * traceG;
+ for (int ii=0; ii<nCompo; ii++)
+ stressG[ii][ii] += lambda(quadPos) * traceG;
- double energyDensityfG = 0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
- energyDensityfG += stressG[ii][jj] * forceTerm(quadPos)[ii][jj];
+ double energyDensityfG = 0;
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ energyDensityfG += stressG[ii][jj] * forceTerm(quadPos)[ii][jj];
- elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
- }
+ elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
+ }
@@ -608,83 +608,83 @@ void computeElementLoadVector( const LocalView& localView,
template<class Basis, class Matrix, class LocalFunction1, class LocalFunction2>
void assembleCellStiffness(const Basis& basis,
- LocalFunction1& muLocal,
- LocalFunction2& lambdaLocal,
- const double gamma,
- Matrix& matrix
+ LocalFunction1& muLocal,
+ LocalFunction2& lambdaLocal,
+ const double gamma,
+ Matrix& matrix
)
{
- std::cout << "assemble stiffnessmatrix begins." << std::endl;
+ std::cout << "assemble stiffnessmatrix begins." << std::endl;
- MatrixIndexSet occupationPattern;
- getOccupationPattern(basis, occupationPattern);
- occupationPattern.exportIdx(matrix);
- matrix = 0;
+ MatrixIndexSet occupationPattern;
+ getOccupationPattern(basis, occupationPattern);
+ occupationPattern.exportIdx(matrix);
+ matrix = 0;
- auto localView = basis.localView();
- const int phiOffset = basis.dimension();
+ auto localView = basis.localView();
+ const int phiOffset = basis.dimension();
- for (const auto& element : elements(basis.gridView()))
- {
- localView.bind(element);
- muLocal.bind(element);
- lambdaLocal.bind(element);
+ for (const auto& element : elements(basis.gridView()))
+ {
+ localView.bind(element);
+ muLocal.bind(element);
+ lambdaLocal.bind(element);
- const int localPhiOffset = localView.size();
-// std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
+ const int localPhiOffset = localView.size();
+ // std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
-// Dune::Matrix<double> elementMatrix;
- Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
- computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
-// printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
+ // Dune::Matrix<double> elementMatrix;
+ Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
+ computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
+ // printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
-// std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
-// std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
+ // std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
+ // std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
- //////////////////////////////////////////////////////////////////////////////
- // GLOBAL STIFFNES ASSEMBLY
- //////////////////////////////////////////////////////////////////////////////
- for (size_t i=0; i<localPhiOffset; i++)
- {
+ //////////////////////////////////////////////////////////////////////////////
+ // GLOBAL STIFFNES ASSEMBLY
+ //////////////////////////////////////////////////////////////////////////////
+ for (size_t i=0; i<localPhiOffset; i++)
+ {
- auto row = localView.index(i);
+ auto row = localView.index(i);
- for (size_t j=0; j<localPhiOffset; j++ )
- {
+ for (size_t j=0; j<localPhiOffset; j++ )
+ {
- auto col = localView.index(j);
+ auto col = localView.index(j);
- matrix[row][col] += elementMatrix[i][j];
- }
- }
+ matrix[row][col] += elementMatrix[i][j];
+ }
+ }
- for (size_t i=0; i<localPhiOffset; i++)
- for(size_t m=0; m<3; m++)
- {
- auto row = localView.index(i);
+ for (size_t i=0; i<localPhiOffset; i++)
+ 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];
- }
+ 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];
- }
+ 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];
+ }
- }
+ }
}
@@ -692,57 +692,57 @@ void assembleCellStiffness(const Basis& basis,
template<class Basis, class LocalFunction1, class LocalFunction2, class Vector, class Force>
void assembleCellLoad(const Basis& basis,
- LocalFunction1& muLocal,
- LocalFunction2& lambdaLocal,
- const double gamma,
- Vector& b,
- const Force& forceTerm
+ LocalFunction1& muLocal,
+ LocalFunction2& lambdaLocal,
+ const double gamma,
+ Vector& b,
+ const Force& forceTerm
)
{
-// std::cout << "assemble load vector." << std::endl;
+ // std::cout << "assemble load vector." << std::endl;
- b.resize(basis.size()+3);
- b = 0;
+ b.resize(basis.size()+3);
+ b = 0;
- auto localView = basis.localView();
+ auto localView = basis.localView();
- const int phiOffset = basis.dimension();
+ const int phiOffset = basis.dimension();
- // Transform G_alpha's to GridViewFunctions/LocalFunctions -- why exactly?
- auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
- auto loadFunctional = localFunction(loadGVF); // Dune::Functions:: notwendig?
+ // Transform G_alpha's to GridViewFunctions/LocalFunctions -- why exactly?
+ auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
+ auto loadFunctional = localFunction(loadGVF); // Dune::Functions:: notwendig?
- for (const auto& element : elements(basis.gridView()))
- {
+ for (const auto& element : elements(basis.gridView()))
+ {
- localView.bind(element);
- muLocal.bind(element);
- lambdaLocal.bind(element);
- loadFunctional.bind(element);
+ localView.bind(element);
+ muLocal.bind(element);
+ lambdaLocal.bind(element);
+ loadFunctional.bind(element);
- const int localPhiOffset = localView.size();
-// std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
+ 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);
+ // BlockVector<double> elementRhs;
+ BlockVector<FieldVector<double,1> > elementRhs;
+ computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, loadFunctional);
- // 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);
+ // 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];
- }
+ b[row] += elementRhs[p];
+ }
- for (size_t m=0; m<3; m++ )
- {
- b[phiOffset+m] += elementRhs[localPhiOffset+m];
- }
+ for (size_t m=0; m<3; m++ )
+ {
+ b[phiOffset+m] += elementRhs[localPhiOffset+m];
}
+ }
}
@@ -762,73 +762,73 @@ auto energy(const Basis& basis,
LocalFunction2& lambdaLocal,
const MatrixFunction& matrixFieldFuncA,
const MatrixFunction& matrixFieldFuncB)
- {
+{
- auto energy = 0.0;
- constexpr int dim = 3;
- constexpr int nCompo = 3;
+ auto energy = 0.0;
+ constexpr int dim = 3;
+ constexpr int nCompo = 3;
- auto localView = basis.localView();
+ auto localView = basis.localView();
- auto matrixFieldAGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncA, basis.gridView());
- auto matrixFieldA = localFunction(matrixFieldAGVF);
- auto matrixFieldBGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
- auto matrixFieldB = localFunction(matrixFieldBGVF);
+ auto matrixFieldAGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncA, basis.gridView());
+ auto matrixFieldA = localFunction(matrixFieldAGVF);
+ auto matrixFieldBGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
+ auto matrixFieldB = localFunction(matrixFieldBGVF);
- using GridView = typename Basis::GridView;
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+ using GridView = typename Basis::GridView;
+ using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+ using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
- for (const auto& e : elements(basis.gridView()))
- {
- localView.bind(e);
- matrixFieldA.bind(e);
- matrixFieldB.bind(e);
- muLocal.bind(e);
- lambdaLocal.bind(e);
+ for (const auto& e : elements(basis.gridView()))
+ {
+ localView.bind(e);
+ matrixFieldA.bind(e);
+ matrixFieldB.bind(e);
+ muLocal.bind(e);
+ lambdaLocal.bind(e);
- FieldVector<double,1> elementEnergy(0);
+ FieldVector<double,1> elementEnergy(0);
- auto geometry = e.geometry();
- const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+ auto geometry = e.geometry();
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- int order = 2*(dim*localFiniteElement.localBasis().order()-1);
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), order);
+ int order = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), order);
- for (size_t pt=0; pt < quad.size(); pt++) {
- const Domain& quadPos = quad[pt].position();
- const double integrationElement = geometry.integrationElement(quadPos);
+ for (size_t pt=0; pt < quad.size(); pt++) {
+ const Domain& quadPos = quad[pt].position();
+ const double integrationElement = geometry.integrationElement(quadPos);
- auto strain1 = matrixFieldA(quadPos);
+ auto strain1 = matrixFieldA(quadPos);
- // St.Venant-Kirchhoff stress
- MatrixRT stressU(0);
- stressU.axpy(2*muLocal(quadPos), strain1); //this += 2mu *strainU // eigentlich this += 2mu *strain1 ?
+ // St.Venant-Kirchhoff stress
+ MatrixRT stressU(0);
+ stressU.axpy(2*muLocal(quadPos), strain1); //this += 2mu *strainU // eigentlich this += 2mu *strain1 ?
- double trace = 0;
- for (int ii=0; ii<nCompo; ii++)
- trace += strain1[ii][ii];
+ double trace = 0;
+ for (int ii=0; ii<nCompo; ii++)
+ trace += strain1[ii][ii];
- for (int ii=0; ii<nCompo; ii++)
- stressU[ii][ii] += lambdaLocal(quadPos) * trace;
+ for (int ii=0; ii<nCompo; ii++)
+ stressU[ii][ii] += lambdaLocal(quadPos) * trace;
- auto strain2 = matrixFieldB(quadPos);
+ auto strain2 = matrixFieldB(quadPos);
- // 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] * strain2[ii][jj];
+ // 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] * strain2[ii][jj];
- elementEnergy += energyDensity * quad[pt].weight() * integrationElement;
- }
- energy += elementEnergy;
- }
- return energy;
- }
+ elementEnergy += energyDensity * quad[pt].weight() * integrationElement;
+ }
+ energy += elementEnergy;
+ }
+ return energy;
+}
@@ -861,43 +861,43 @@ auto energy(const Basis& basis,
template<class Basis, class Matrix, class Vector>
void setOneBaseFunctionToZero(const Basis& basis,
- Matrix& stiffnessMatrix,
- Vector& load_alpha1,
- Vector& load_alpha2,
- Vector& load_alpha3,
- FieldVector<int,3>& indexVector // oder input : arbitraryLocalIndex und dann hier function computeArbitraryIndex aufrufen
- )
+ Matrix& stiffnessMatrix,
+ Vector& load_alpha1,
+ Vector& load_alpha2,
+ Vector& load_alpha3,
+ FieldVector<int,3>& indexVector // oder input : arbitraryLocalIndex und dann hier function computeArbitraryIndex aufrufen
+ )
{
- std::cout << "setting one Basis-function to zero" << std::endl;
+ std::cout << "setting one Basis-function to zero" << std::endl;
- constexpr int dim = 3;
- auto localView = basis.localView();
+ constexpr int dim = 3;
+ auto localView = basis.localView();
- //Determine 3 global indices (one for each component of an arbitrary FE-function).. oder mittels Offset?
+ //Determine 3 global indices (one for each component of an arbitrary FE-function).. oder mittels Offset?
- // use first element:
- auto it = basis.gridView().template begin<0>();
- localView.bind(*it);
+ // use first element:
+ auto it = basis.gridView().template begin<0>();
+ localView.bind(*it);
- int arbitraryIndex = 0;
- FieldVector<int,3> row; //fill with Indices..
+ int arbitraryIndex = 0;
+ FieldVector<int,3> row; //fill with Indices..
- for (int k = 0; k < dim; k++)
- {
- auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
- row[k] = localView.index(rowLocal);
-// std::cout << "row[k]:" << row[k] << std::endl;
- load_alpha1[row[k]] = 0.0; //verwende hier indexVector
- load_alpha2[row[k]] = 0.0;
- load_alpha3[row[k]] = 0.0;
- auto cIt = stiffnessMatrix[row[k]].begin();
- auto cEndIt = stiffnessMatrix[row[k]].end();
- for (; cIt!=cEndIt; ++cIt)
- *cIt = (cIt.index()==row[k]) ? 1.0 : 0.0;
- }
+ for (int k = 0; k < dim; k++)
+ {
+ auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
+ row[k] = localView.index(rowLocal);
+ // std::cout << "row[k]:" << row[k] << std::endl;
+ load_alpha1[row[k]] = 0.0; //verwende hier indexVector
+ load_alpha2[row[k]] = 0.0;
+ load_alpha3[row[k]] = 0.0;
+ auto cIt = stiffnessMatrix[row[k]].begin();
+ auto cEndIt = stiffnessMatrix[row[k]].end();
+ for (; cIt!=cEndIt; ++cIt)
+ *cIt = (cIt.index()==row[k]) ? 1.0 : 0.0;
+ }
}
@@ -912,114 +912,114 @@ void setIntegralZero (const Basis& basis,
Vector& load_alpha3
)
{
- //
- // 0. moment invariance, 3 dofs
- //
+ //
+ // 0. moment invariance, 3 dofs
+ //
+
+ // if (parameterSet_.get<bool>("set_integral_0")){
+ std::cout << "\nassembling: integral to zero\n";
+
+
+ auto n = basis.size();
+ // auto nNodeBasis = n/nCompo;
+
+ constexpr int dim = 3; // irgendwo herziehen?
+ constexpr int nCompo = dim;
+
+ auto localView = basis.localView();
+
+ int arbitraryIndex = 0; // sollte GlOBALER INDEX sein
+
+ // unsigned long row1;
+ FieldVector<int,3> row; // fill with Indices..
-// if (parameterSet_.get<bool>("set_integral_0")){
- std::cout << "\nassembling: integral to zero\n";
+ //Determine 3 global indices (one for each component of an arbitrary FE-function).. oder mittels Offset?
- auto n = basis.size();
-// auto nNodeBasis = n/nCompo;
+ // for (int k = 0; k < dim; k++)
+ // {
+ //
+ // auto row = localView.tree().child(k).localIndex(arbitraryIndex);
+ // // auto row = perIndexMap_[arbitraryIndex + d*nNodeBasis];
+ // stiffnessMatrix[row] = 0; // setzt gesamte Zeile auf 0 !!!
+ // load_alpha1[row] = 0.0;
+ // load_alpha2[row] = 0.0;
+ // load_alpha3[row] = 0.0;
+ // }
- constexpr int dim = 3; // irgendwo herziehen?
- constexpr int nCompo = dim;
+ unsigned int elementCounter = 1;
+ for (const auto& element : elements(basis.gridView()))
+ {
+ // if (elementCounter % 50 == 1)
+ // std::cout << "integral = zero - treatment - element: " << elementCounter << std::endl;
+ // std::cout << "element-number: " << elementCounter << std::endl;
+ localView.bind(element);
- auto localView = basis.localView();
- int arbitraryIndex = 0; // sollte GlOBALER INDEX sein
+ if(elementCounter == 1) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
+ {
+ for (int k = 0; k < dim; k++)
+ {
+ auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
+ row[k] = localView.index(rowLocal);
+ // std::cout << "row[k]:" << row[k] << std::endl;
+ stiffnessMatrix[row[k]] = 0;
+ load_alpha1[row[k]] = 0.0;
+ load_alpha2[row[k]] = 0.0;
+ load_alpha3[row[k]] = 0.0;
+ }
+ }
-// unsigned long row1;
- FieldVector<int,3> row; // fill with Indices..
+ // "local assembly"
+ BlockVector<FieldVector<double,1> > elementColumns;
+ elementColumns.resize(localView.size()+3);
+ elementColumns = 0;
- //Determine 3 global indices (one for each component of an arbitrary FE-function).. oder mittels Offset?
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+ const auto nSf = localFiniteElement.localBasis().size();
-// for (int k = 0; k < dim; k++)
-// {
-//
-// auto row = localView.tree().child(k).localIndex(arbitraryIndex);
-// // auto row = perIndexMap_[arbitraryIndex + d*nNodeBasis];
-// stiffnessMatrix[row] = 0; // setzt gesamte Zeile auf 0 !!!
-// load_alpha1[row] = 0.0;
-// load_alpha2[row] = 0.0;
-// load_alpha3[row] = 0.0;
-// }
+ int order = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), order);
- unsigned int elementCounter = 1;
- for (const auto& element : elements(basis.gridView()))
+ for (const auto& quadPoint : quad)
+ {
+ // for (size_t pt=0; pt < quad.size(); pt++) {
+
+ const auto& quadPos = quadPoint.position();
+ // const FieldVector<double,dim>& quadPos = quad[pt].position();
+ const double integrationElement = element.geometry().integrationElement(quadPos);
+ const auto jacobian = element.geometry().jacobianInverseTransposed(quadPos);
+
+ std::vector<FieldVector<double,1> > shapeFunctionValues;
+ localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
+
+ // Compute the shape function on the grid element
+ // std::vector<FieldVector<double,1> > BasisFunctionValues(shapeFunctionValues.size());
+ // for (size_t i=0; i<BasisFunctionValues.size(); i++)
+ // jacobian.mv(shapeFunctionValues[i], BasisFunctionValues[i]); // no gradient! not needed!
+
+ for (size_t i=0; i<localView.size(); i++)
+ for (size_t k=0; k<nCompo; k++)
{
-// if (elementCounter % 50 == 1)
-// std::cout << "integral = zero - treatment - element: " << elementCounter << std::endl;
-// std::cout << "element-number: " << elementCounter << std::endl;
- localView.bind(element);
-
-
- if(elementCounter == 1) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
- {
- for (int k = 0; k < dim; k++)
- {
- auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
- row[k] = localView.index(rowLocal);
-// std::cout << "row[k]:" << row[k] << std::endl;
- stiffnessMatrix[row[k]] = 0;
- load_alpha1[row[k]] = 0.0;
- load_alpha2[row[k]] = 0.0;
- load_alpha3[row[k]] = 0.0;
- }
- }
-
-
- // "local assembly"
- BlockVector<FieldVector<double,1> > elementColumns;
- elementColumns.resize(localView.size()+3);
- elementColumns = 0;
-
- const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- const auto nSf = localFiniteElement.localBasis().size();
-
- int order = 2*(dim*localFiniteElement.localBasis().order()-1);
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), order);
-
- for (const auto& quadPoint : quad)
- {
- // for (size_t pt=0; pt < quad.size(); pt++) {
-
- const auto& quadPos = quadPoint.position();
- // const FieldVector<double,dim>& quadPos = quad[pt].position();
- const double integrationElement = element.geometry().integrationElement(quadPos);
- const auto jacobian = element.geometry().jacobianInverseTransposed(quadPos);
-
- std::vector<FieldVector<double,1> > shapeFunctionValues;
- localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
-
- // Compute the shape function on the grid element
-// std::vector<FieldVector<double,1> > BasisFunctionValues(shapeFunctionValues.size());
-// for (size_t i=0; i<BasisFunctionValues.size(); i++)
-// jacobian.mv(shapeFunctionValues[i], BasisFunctionValues[i]); // no gradient! not needed!
-
- for (size_t i=0; i<localView.size(); i++)
- for (size_t k=0; k<nCompo; k++)
- {
- size_t idx = localView.tree().child(k).localIndex(i);
- elementColumns[idx] += shapeFunctionValues[i] * quadPoint.weight() * integrationElement; // bei Periodicbasis einfach mit child(k) arbeiten...
-// elementColumns[idx] += BasisFunctionValues[i] * quadPoint.weight() * integrationElement; // bei Periodicbasis einfach mit child(k) arbeiten...
- }
-
- // "global assembly"
- for (int k = 0; k<dim; k++)
- for (size_t j=0; j<localView.size(); j++)
- {
-// auto colLocal = localView.tree().child(d).localIndex(j);
- auto col = localView.index(j);
- stiffnessMatrix[row[k]][col] += elementColumns[j] * quadPoint.weight() * integrationElement; // TODO Müsste über eine LOOP ?
- }
-
- elementCounter++;
- }//end of quad
-
- }//end for each element
+ size_t idx = localView.tree().child(k).localIndex(i);
+ elementColumns[idx] += shapeFunctionValues[i] * quadPoint.weight() * integrationElement; // bei Periodicbasis einfach mit child(k) arbeiten...
+ // elementColumns[idx] += BasisFunctionValues[i] * quadPoint.weight() * integrationElement; // bei Periodicbasis einfach mit child(k) arbeiten...
+ }
+
+ // "global assembly"
+ for (int k = 0; k<dim; k++)
+ for (size_t j=0; j<localView.size(); j++)
+ {
+ // auto colLocal = localView.tree().child(d).localIndex(j);
+ auto col = localView.index(j);
+ stiffnessMatrix[row[k]][col] += elementColumns[j] * quadPoint.weight() * integrationElement; // TODO Müsste über eine LOOP ?
+ }
+
+ elementCounter++;
+ } //end of quad
+
+ } //end for each element
}
@@ -1028,58 +1028,58 @@ void setIntegralZero (const Basis& basis,
template<class Basis>
auto childToIndexMap(const Basis& basis,
- const int k
- )
+ const int k
+ )
{
- constexpr int dim = 3;
+ constexpr int dim = 3;
- using GridView = typename Basis::GridView;
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+ using GridView = typename Basis::GridView;
+ using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using FuncScalar = std::function< FieldVector< double, dim>(const Domain&) >;
+ using FuncScalar = std::function< FieldVector< double, dim>(const Domain&) >;
- auto localView = basis.localView();
+ auto localView = basis.localView();
- std::vector<int> r = { };
-// for (int n: r)
-// std::cout << n << ","<< std::endl;
+ std::vector<int> r = { };
+ // for (int n: r)
+ // std::cout << n << ","<< std::endl;
- // Determine global indizes for each component k = 1,2,3.. in order to subtract correct (component of) integral Mean
- // (global) Indices that correspond to component k = 1,2,3
- for(const auto& element : elements(basis.gridView()))
- {
+ // Determine global indizes for each component k = 1,2,3.. in order to subtract correct (component of) integral Mean
+ // (global) Indices that correspond to component k = 1,2,3
+ 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();
+ 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();
- for(size_t j=0; j<nSf; j++)
- {
- auto Localidx = localView.tree().child(k).localIndex(j);
- r.push_back(localView.index(Localidx));
- }
+ for(size_t j=0; j<nSf; j++)
+ {
+ auto Localidx = localView.tree().child(k).localIndex(j);
+ r.push_back(localView.index(Localidx));
+ }
- }
+ }
- // Delete duplicate elements
- // first remove consecutive (adjacent) duplicates
- auto last = std::unique(r.begin(), r.end());
- r.erase(last, r.end());
- // sort followed by unique, to remove all duplicates
- std::sort(r.begin(), r.end());
- last = std::unique(r.begin(), r.end());
- r.erase(last, r.end());
+ // Delete duplicate elements
+ // first remove consecutive (adjacent) duplicates
+ auto last = std::unique(r.begin(), r.end());
+ r.erase(last, r.end());
+ // sort followed by unique, to remove all duplicates
+ std::sort(r.begin(), r.end());
+ last = std::unique(r.begin(), r.end());
+ r.erase(last, r.end());
- return r;
+ return r;
}
@@ -1095,77 +1095,77 @@ auto integralMean(const Basis& basis,
LocalFunction& fun
)
{
- // computes integral mean of given LocalFunction
+ // computes integral mean of given LocalFunction
- using GridView = typename Basis::GridView;
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+ using GridView = typename Basis::GridView;
+ using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- constexpr int dim = 3;
- using FuncScalar = std::function< FieldVector< double, dim>(const Domain&) >;
+ constexpr int dim = 3;
+ using FuncScalar = std::function< FieldVector< double, dim>(const Domain&) >;
- auto localView = basis.localView();
-// double r = 0.0;
- FieldVector<double,3> r = {0.0, 0.0, 0.0};
- double area = 0.0;
+ auto localView = basis.localView();
+ // double r = 0.0;
+ FieldVector<double,3> r = {0.0, 0.0, 0.0};
+ double area = 0.0;
- // Compute area integral
- for(const auto& element : elements(basis.gridView()))
- {
+ // Compute area integral
+ for(const auto& element : elements(basis.gridView()))
+ {
- localView.bind(element);
- const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- int order = 2*(dim*localFiniteElement.localBasis().order()-1);
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), order);
+ localView.bind(element);
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+ int order = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), order);
- for (const auto& quadPoint : quad)
- {
- const double integrationElement = element.geometry().integrationElement(quadPoint.position());
- area += quadPoint.weight() * integrationElement;
- }
+ for (const auto& quadPoint : quad)
+ {
+ const double integrationElement = element.geometry().integrationElement(quadPoint.position());
+ area += quadPoint.weight() * integrationElement;
}
-// std::cout << "Domain-Area: " << area << std::endl;
+ }
+ // std::cout << "Domain-Area: " << area << std::endl;
- for(const auto& element : elements(basis.gridView()))
- {
+ for(const auto& element : elements(basis.gridView()))
+ {
- localView.bind(element);
- fun.bind(element);
- const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- int order = 2*(dim*localFiniteElement.localBasis().order()-1);
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), order);
+ localView.bind(element);
+ fun.bind(element);
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+ int order = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), order);
- for (const auto& quadPoint : quad)
- {
+ for (const auto& quadPoint : quad)
+ {
- const auto& quadPos = quadPoint.position();
- const double integrationElement = element.geometry().integrationElement(quadPos);
+ const auto& quadPos = quadPoint.position();
+ const double integrationElement = element.geometry().integrationElement(quadPos);
-// std::vector<FieldMatrix<double,1,nCompo> > referenceGradients;
-// localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
-//
-// std::vector< VectorRT > gradients(referenceGradients.size());
-// for (size_t i=0; i< gradients.size(); i++)
-// jacobian.mv(referenceGradients[i][0], gradients[i]);
+ // std::vector<FieldMatrix<double,1,nCompo> > referenceGradients;
+ // localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+ //
+ // std::vector< VectorRT > gradients(referenceGradients.size());
+ // for (size_t i=0; i< gradients.size(); i++)
+ // jacobian.mv(referenceGradients[i][0], gradients[i]);
- auto value = fun(quadPos);
-// auto value = localPhi_1(quadPos);
-// auto value = FieldVector<double,dim>{1.0, 1.0, 1.0};
+ auto value = fun(quadPos);
+ // auto value = localPhi_1(quadPos);
+ // auto value = FieldVector<double,dim>{1.0, 1.0, 1.0};
-// for(size_t k=0; k<3; k++)
-// r += value[k] * quadPoint.weight() * integrationElement;
- r += value * quadPoint.weight() * integrationElement;
+ // for(size_t k=0; k<3; k++)
+ // r += value[k] * quadPoint.weight() * integrationElement;
+ r += value * quadPoint.weight() * integrationElement;
-// r += tmp * quadPoint.weight() * integrationElement;
+ // r += tmp * quadPoint.weight() * integrationElement;
- }
}
+ }
-// return r/area;
- return (1/area) * r ;
+ // return r/area;
+ return (1/area) * r ;
}
@@ -1181,19 +1181,19 @@ auto subtractIntegralMean(const Basis& basis,
Vector& coeffVector
)
{
- // subtract correct Integral mean from each associated component function
+ // subtract correct Integral mean from each associated component function
- auto IM = integralMean(basis, fun);
- constexpr int dim = 3;
+ auto IM = integralMean(basis, fun);
+ constexpr int dim = 3;
- for(size_t k=0; k<dim; k++)
- {
-// std::cout << "Integral-Mean: " << IM[k] << std::endl;
- auto idx = childToIndexMap(basis,k);
+ for(size_t k=0; k<dim; k++)
+ {
+ // std::cout << "Integral-Mean: " << IM[k] << std::endl;
+ auto idx = childToIndexMap(basis,k);
- for ( int i : idx)
- coeffVector[i] -= IM[k];
- }
+ for ( int i : idx)
+ coeffVector[i] -= IM[k];
+ }
}
@@ -1207,13 +1207,13 @@ auto subtractIntegralMean(const Basis& basis,
- // 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])));
- };
+// 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])));
+ };
@@ -1235,424 +1235,369 @@ auto subtractIntegralMean(const Basis& basis,
int main(int argc, char *argv[])
{
- MPIHelper::instance(argc, argv);
-
-
-
- ParameterTree parameterSet;
-// if (argc < 2)
-// ParameterTreeParser::readINITree("../src/cellsolver.parset", parameterSet);
-// ParameterTreeParser::readINITree(argv[1], parameterSet);
-// ParameterTreeParser::readOptions(argc, argv, parameterSet);
- // output setter
- // -- not set up --
-
-
-
- constexpr int dim = 3;
- constexpr int nCompo = 3;
- constexpr int order_CE = 1;
-
- ///////////////////////////////////
- // Get Parameters/Data
- ///////////////////////////////////
-
- double gamma = parameterSet.get<double>("gamma",10000); // ratio dimension reduction to homogenization
-
- // Material parameter laminat
-// double E1 = parameterSet.get<double>("E1", 17e6); //material1
-// double nu1 = parameterSet.get<double>("nu1", 0.0);
-// double nu1 = parameterSet.get<double>("nu1", 0.3);
-// double E2 = parameterSet.get<double>("E2", 35e6); //material2
-// double nu2 = parameterSet.get<double>("nu2", 0.5);
+ MPIHelper::instance(argc, argv);
-// // 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
+ ParameterTree parameterSet;
+ // if (argc < 2)
+ // ParameterTreeParser::readINITree("../src/cellsolver.parset", parameterSet);
+ // ParameterTreeParser::readINITree(argv[1], parameterSet);
+ // ParameterTreeParser::readOptions(argc, argv, parameterSet);
+
+ // output setter
+ // -- not set up --
+ std::fstream log;
+// log.open("output2.txt", std::ios::out);
+ log.open("../../outputs/output.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;
- ///////////////////////////////////
- // ---Get Prestrain ---
- ///////////////////////////////////
- unsigned int prestraintype = parameterSet.get<unsigned int>("imp", 1);
- double p1 = parameterSet.get<double>("prestress_pressure1", 1.0);
- double p2 = parameterSet.get<double>("prestress_pressure2", 2.0);
- double theta = parameterSet.get<double>("theta",0.3);
-// double theta = 0.5;
- p1 = 1.0;
- double alpha = 2;
- p2 = alpha*1.0;
- prestraintype = 2;
-// prestraintype = 1;
+ ///////////////////////////////////
+ // Get Parameters/Data
+ ///////////////////////////////////
- auto prestrainImp = PrestrainImp(p1, p2, theta, width);
- auto B_Term = prestrainImp.getPrestrain(prestraintype);
+ double gamma = parameterSet.get<double>("gamma",5); // ratio dimension reduction to homogenization
+ // Material parameter laminat
+ // double E1 = parameterSet.get<double>("E1", 17e6); //material1
+ // double nu1 = parameterSet.get<double>("nu1", 0.0);
+ // double nu1 = parameterSet.get<double>("nu1", 0.3);
+ // double E2 = parameterSet.get<double>("E2", 35e6); //material2
+ // double nu2 = parameterSet.get<double>("nu2", 0.5);
- ///////////////////////////////////
- // Generate the grid
- ///////////////////////////////////
- using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
- // Domain : [0,1)^2 x (-1/2, 1/2)
-// FieldVector<double,dim> lower({0.0, 0.0, -1.0/2.0});
-// FieldVector<double,dim> upper({1.0, 1.0, 1.0/2.0});
+ // // 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
- // (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});
- int nE = 10;
+ ///////////////////////////////////
+ // ---Get Prestrain ---
+ ///////////////////////////////////
+ unsigned int prestraintype = parameterSet.get<unsigned int>("imp", 1);
+ double p1 = parameterSet.get<double>("prestress_pressure1", 1.0);
+ double p2 = parameterSet.get<double>("prestress_pressure2", 2.0);
+ double theta = parameterSet.get<double>("theta",0.3);
+ // double theta = 0.5;
+ p1 = 1.0;
+ double alpha = 2;
+ p2 = alpha*1.0;
+ prestraintype = 2;
+ // prestraintype = 1;
- std::array<int,dim> nElements={nE,nE,nE}; //#Elements in each direction
+ auto prestrainImp = PrestrainImp(p1, p2, theta, width);
+ auto B_Term = prestrainImp.getPrestrain(prestraintype);
- CellGridType grid_CE(lower,upper,nElements); //Corrector problem Domain:
+ log << "prestrain imp: " << prestraintype << "\np1 = " << p1 << "\np2 = " << p2 << std::endl;
+ log << "alpha: " << alpha << std::endl;
+ log << "gamma: " << gamma << std::endl;
+
- using GridView = CellGridType::LeafGridView;
- const GridView gridView_CE = grid_CE.leafGridView();
+ ///////////////////////////////////
+ // Generate the grid
+ ///////////////////////////////////
+ using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim> >;
+ // Domain : [0,1)^2 x (-1/2, 1/2)
+ // FieldVector<double,dim> lower({0.0, 0.0, -1.0/2.0});
+ // FieldVector<double,dim> upper({1.0, 1.0, 1.0/2.0});
+ // (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});
+ int nE = 10;
+ log << "Number of Elements: " << nE << std::endl;
+
+ std::array<int,dim> nElements={nE,nE,nE}; //#Elements in each direction
+ CellGridType grid_CE(lower,upper,nElements); //Corrector problem Domain:
+ using GridView = CellGridType::LeafGridView;
+ const GridView gridView_CE = grid_CE.leafGridView();
- using ScalarRT = FieldVector< double, 1>;
- using VectorRT = FieldVector< double, nCompo>;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
- 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> >;
- // using VectorCT = BlockVector<FieldVector<double,dim> >;
- // using MatrixCT = BCRSMatrix<FieldMatrix<double,dim,dim> >;
- using ElementMatrixCT = Matrix<FieldMatrix<double,1,1> >;
+ using ScalarRT = FieldVector< double, 1>;
+ using VectorRT = FieldVector< double, nCompo>;
+ using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+ using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+ using FuncScalar = std::function< ScalarRT(const Domain&) >;
+ using Func2Tensor = std::function< MatrixRT(const Domain&) >;
- double beta = 2;
+ using VectorCT = BlockVector<FieldVector<double,1> >;
+ using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >;
+ // using VectorCT = BlockVector<FieldVector<double,dim> >;
+ // using MatrixCT = BCRSMatrix<FieldMatrix<double,dim,dim> >;
+ using ElementMatrixCT = Matrix<FieldMatrix<double,1,1> >;
-// double mu1 = 10;
- double mu1 = 0.5*17e6;
-// double mu1 = 1000;
- double mu2 = beta*mu1;
-// //
-// double mu1 = 0.5 * 1.0/(1.0 + nu1) * E1;
-// double lambda1 = nu1/(1.0-2.0*nu1) * 1.0/(1.0+nu1) * E1;
- double lambda1 = 0;
- // Create Lambda-Functions for material Parameters
-// auto muTerm = [mu1] (const Domain& z) {
-// return mu1;
-// };
- auto muTerm = [mu1, mu2, theta] (const Domain& z) {
- if (abs(z[0]) >= (theta/2))
- return mu1;
-// if (abs(z[0]) < (theta/2) && z[2] < 0) // oder das?
-// return mu2;
- else
- return mu2;
-// return mu1;
+ double beta = 2;
+ double mu1 = 10;
+// double mu1 = 0.5*17e6;
+
+ // double mu1 = 1000;
+ double mu2 = beta*mu1;
+ // //
- };
+ // double mu1 = 0.5 * 1.0/(1.0 + nu1) * E1;
+ // double lambda1 = nu1/(1.0-2.0*nu1) * 1.0/(1.0+nu1) * E1;
+ double lambda1 = 0;
+ // Create Lambda-Functions for material Parameters
+ // auto muTerm = [mu1] (const Domain& z) {
+ // return mu1;
+ // };
+ auto muTerm = [mu1, mu2, theta] (const Domain& z) {
+ if (abs(z[0]) >= (theta/2))
+ return mu1;
+ // if (abs(z[0]) < (theta/2) && z[2] < 0) // oder das?
+ // return mu2;
+ else
+ return mu2;
+ // return mu1;
- auto lambdaTerm = [lambda1] (const Domain& z) {
- return lambda1;
- };
+ };
- auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
- auto muLocal = localFunction(muGridF);
- auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
- auto lambdaLocal = localFunction(lambdaGridF);
+ auto lambdaTerm = [lambda1] (const Domain& z) {
+ return lambda1;
+ };
- /////////////////////////////////////////////////////////
- // Choose a finite element space for Cell Problem
- /////////////////////////////////////////////////////////
- using namespace Functions::BasisFactory;
- Functions::Experimental::BasisFactory::PeriodicIndexSet periodicIndices;
- // 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))
- 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)});
+ 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 << "\n lambda: " << lambda1 << std::endl;
+
+
+
- auto Basis_CE = makeBasis(
- gridView_CE,
- power<dim>(
-// lagrange<1>(),
- Functions::Experimental::BasisFactory::periodic(lagrange<1>(), periodicIndices),
- flatLexicographic()));//
-// flatInterleaved()));
+ /////////////////////////////////////////////////////////
+ // Choose a finite element space for Cell Problem
+ /////////////////////////////////////////////////////////
+ using namespace Functions::BasisFactory;
+ Functions::Experimental::BasisFactory::PeriodicIndexSet periodicIndices;
- std::cout << "size feBasis: " << Basis_CE.size() << std::endl;
- std::cout << "basis.dimension()" << Basis_CE.dimension() <<std::endl;
+ // 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))
+ 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)});
- const int phiOffset = Basis_CE.size();
+ auto Basis_CE = makeBasis(
+ gridView_CE,
+ power<dim>(
+ // lagrange<1>(),
+ Functions::Experimental::BasisFactory::periodic(lagrange<1>(), periodicIndices),
+ flatLexicographic())); //
+ // flatInterleaved()));
- // TEST
-// auto X1Basis = subspaceBasis(Basis_CE, 0);
-// auto X2Basis = subspaceBasis(Basis_CE, 1);
-// auto X3Basis = subspaceBasis(Basis_CE, 2);
-//
-// std::cout << "size X1Basis: " << X1Basis.size() << std::endl;
-// std::cout << "X1Basis.dimension()" << X1Basis.dimension() <<std::endl;
-// std::cout << "size X2Basis: " << X2Basis.size() << std::endl;
-// std::cout << "size X3Basis: " << X3Basis.size() << std::endl;
+ std::cout << "size feBasis: " << Basis_CE.size() << std::endl;
+ std::cout << "basis.dimension()" << Basis_CE.dimension() <<std::endl;
+ const int phiOffset = Basis_CE.size();
-/////////////////////////////////////////////////////////
-// Stiffness matrix and right hand side vector
-/////////////////////////////////////////////////////////
-VectorCT load_alpha1, load_alpha2, load_alpha3;
-MatrixCT stiffnessMatrix_CE;
-// auto rhsBackend = Functions::istlVectorBackend(b);
-// rhsBackend.resize(Basis_CE);
-// b = 0;
+// std::cout << Basis_CE.preBasis_.children << std::endl; // does not work
+// std::cout << Basis_CE::preBasis_::children << std::endl;
+ // TEST
+ // auto X1Basis = subspaceBasis(Basis_CE, 0);
+ // auto X2Basis = subsppaceBasis(Basis_CE, 1);
+ // auto X3Basis = subspaceBasis(Basis_CE, 2);
+ //
+ // std::cout << "size X1Basis: " << X1Basis.size() << std::endl;
+ // std::cout << "X1Basis.dimension()" << X1Basis.dimension() <<std::endl;
+ // std::cout << "size X2Basis: " << X2Basis.size() << std::endl;
+ // std::cout << "size X3Basis: " << X3Basis.size() << std::endl;
+ /////////////////////////////////////////////////////////
+ // Stiffness matrix and right hand side vector
+ /////////////////////////////////////////////////////////
+ VectorCT load_alpha1, load_alpha2, load_alpha3;
+ MatrixCT stiffnessMatrix_CE;
+ // auto rhsBackend = Functions::istlVectorBackend(b);
+ // rhsBackend.resize(Basis_CE);
+ // b = 0;
-// bool set_integral_zero = true;
-bool set_oneBasisFunction_Zero = true;
-bool substract_integralMean = true;
-// bool set_oneBasisFunction_Zero = false;
-// bool substract_integralMean = false;
-// using Range = FieldVector<double,dim>;
-// auto sourceTerm = [](const FieldVector<double,dim>& x){return Range{0.0, -1.0};};
+ // bool set_integral_zero = true;
+ bool set_oneBasisFunction_Zero = true;
+ bool substract_integralMean = true;
+ // bool set_oneBasisFunction_Zero = false;
+ // bool substract_integralMean = false;
-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_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) {
- return MatrixRT{{0.0, 1/sqrt(2)*x[2], 0.0}, {1/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}}; };
+ // using Range = FieldVector<double,dim>;
+ // auto sourceTerm = [](const FieldVector<double,dim>& x){return Range{0.0, -1.0};};
+ 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_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}};
+ };
-/////////////////////////////////////////////////////////////////////// TODO
-// TODO : PrestrainImp.hh
-// double theta = 0.5;
-// double p1 = 1;
-// double p2 = 2;
+ Func2Tensor x3G_3 = [] (const Domain& x) {
+ return MatrixRT{{0.0, 1/sqrt(2)*x[2], 0.0}, {1/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ };
-using std::abs;
-Func2Tensor B = [] (const Domain& x) {
- double theta = 0.5;
- double p1 = 1;
- double p2 = 2;
- if (abs(x[0])>= (theta/2) && x[2] > 0)
- return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
- if (abs(x[0]) < (theta/2) && x[2] < 0)
- return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
- else
- return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ /////////////////////////////////////////////////////////////////////// TODO
+ // TODO : PrestrainImp.hh
+ // double theta = 0.5;
+ // double p1 = 1;
+ // double p2 = 2;
-};
-//////////////////////////////////////////////////////////////////////////////
+// using std::abs;
+// Func2Tensor B = [] (const Domain& x) {
+// double theta = 0.5;
+// double p1 = 1;
+// double p2 = 2;
+// if (abs(x[0])>= (theta/2) && x[2] > 0)
+// return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+// if (abs(x[0]) < (theta/2) && x[2] < 0)
+// return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
+// else
+// return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+//
+// };
+ //////////////////////////////////////////////////////////////////////////////
-///////////////////////////////////////////////
-// Basis for R_sym^{2x2}
-//////////////////////////////////////////////
-MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
-MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
-MatrixRT G_3 {{0.0, 1/sqrt(2), 0.0}, {1/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
-std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
-// printmatrix(std::cout, basisContainer[0] , "G_1", "--");
-// printmatrix(std::cout, basisContainer[1] , "G_2", "--");
-// printmatrix(std::cout, basisContainer[2] , "´G_3", "--");
+ ///////////////////////////////////////////////
+ // Basis for R_sym^{2x2}
+ //////////////////////////////////////////////
+ MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+ MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+ MatrixRT G_3 {{0.0, 1/sqrt(2), 0.0}, {1/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
+ // printmatrix(std::cout, basisContainer[0] , "G_1", "--");
+ // printmatrix(std::cout, basisContainer[1] , "G_2", "--");
+ // printmatrix(std::cout, basisContainer[2] , "´G_3", "--");
-/////////////////////////////////////////////////////////
-// Assemble the system
-/////////////////////////////////////////////////////////
-assembleCellStiffness(Basis_CE, muLocal, lambdaLocal, gamma, stiffnessMatrix_CE);
-assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha1 ,x3G_1);
-assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha2 ,x3G_2);
-assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3);
+ /////////////////////////////////////////////////////////
+ // Assemble the system
+ /////////////////////////////////////////////////////////
+ assembleCellStiffness(Basis_CE, muLocal, lambdaLocal, gamma, stiffnessMatrix_CE);
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha1 ,x3G_1);
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha2 ,x3G_2);
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3);
-// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix before B.C", "--");
+ // printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix before B.C", "--");
-// set Integral to zero
-// if(set_integral_zero)
-// {
-// setIntegralZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3);
-// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setIntegralZero", "--");
-// }
-//
+ // set Integral to zero
+ // if(set_integral_zero)
+ // {
+ // setIntegralZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3);
+ // printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setIntegralZero", "--");
+ // }
+ //
-// Determine global indixes of components for arbitrary (local) index TODO :separate Function!
+ /////////////////////////////////////////////////////////
+ // Determine global indices of components for arbitrary (local) index TODO :separate Function! arbitraryComponentwiseIndices
+ /////////////////////////////////////////////////////////
+ auto localView = Basis_CE.localView();
-auto localView = Basis_CE.localView();
+ int arbitraryIndex = 0;
-int arbitraryIndex = 0;
+ FieldVector<int,3> row; //fill with Indices..
-FieldVector<int,3> row; //fill with Indices..
-
-
- // use first element:
- auto it = Basis_CE.gridView().template begin<0>();
- localView.bind(*it);
-
- for (int k = 0; k < dim; k++)
- {
- auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
- row[k] = localView.index(rowLocal);
- std::cout << "row[k]:" << row[k] << std::endl;
-
- }
+ // use first element:
+ auto it = Basis_CE.gridView().template begin<0>();
+ localView.bind(*it);
+ for (int k = 0; k < dim; k++)
+ {
+ auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
+ row[k] = localView.index(rowLocal);
+ std::cout << "row[k]:" << row[k] << std::endl;
-// set one basis-function to zero
-if(set_oneBasisFunction_Zero)
-{
- setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3,row);
-// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
-}
+ }
+ ///////////////////////////////////////////////////////////
+
+
+ // set one basis-function to zero
+ if(set_oneBasisFunction_Zero)
+ {
+ setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3,row);
+ // printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
+ }
+ ////////////////////////////
+ // Compute solution
+ ////////////////////////////
+ VectorCT x_1 = load_alpha1;
+ VectorCT x_2 = load_alpha2;
+ VectorCT x_3 = load_alpha3;
+ // Turn the matrix into a linear operator
+ MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_CE);
+ // Fancy (but only) way to not have a preconditioner at all
+ Richardson<VectorCT,VectorCT> preconditioner(1.0);
-/////////////////////////////////////////////////////////
-// Set Dirichlet values.
-// Only velocity components have Dirichlet boundary values
-/////////////////////////////////////////////////////////
-
-
-
-// using BitVector = std::vector<std::array<char,dim> >;
-// BitVector isBoundary;
-//
-// auto isBoundaryBackend = Functions::istlVectorBackend(isBoundary);
-// isBoundaryBackend.resize(Basis);
-//
-//
-//
-// // clamped Boundary Conditions
-// forEachBoundaryDOF(
-// Basis,
-// [&] (auto&& index) {
-// isBoundaryBackend[index] = true;
-// });
-//
-//
-//
-//
-//
-//
-// using Coordinate = GridView::Codim<0> ::Geometry::GlobalCoordinate;
-// // using Range = FieldVector<double,dim+1>;
-// // auto&& g = [](Coordinate x)
-// // {
-// // return Range{0.0, 0.0 , 0.0 };
-// // };
-//
-// auto&& g = [](Coordinate x)
-// {
-// return Range{0.0, 0.0 };
-// };
-//
-// Functions::interpolate(Basis,
-// b,
-// g,
-// isBoundary);
-//
-
-
-////////////////////////////////////////////
-// Modify Dirichlet rows
-////////////////////////////////////////////
-
-// auto localView = Basis.localView();
-// for(const auto& element : elements(gridView))
-// {
-// localView.bind(element);
-// for (size_t i=0; i<localView.size(); ++i)
-// {
-// auto row = localView.index(i);
-// // If row corresponds to a boundary entry,
-// // modify it to be an identity matrix row.
-// if (isBoundaryBackend[row])
-// for (size_t j=0; j<localView.size(); ++j)
-// {
-// auto col = localView.index(j);
-// stiffnessMatrix[row[0]][col[0]][row[1]][col[1]] = (i==j) ? 1 : 0;
-// // matrixEntry(stiffnessMatrix, row, col) = (i==j) ? 1 : 0;
-// }
-// }
-// }
-//
-// printmatrix(std::cout, stiffnessMatrix, "StiffnessMatrix after Dirichlet ", "--");
-
-////////////////////////////
-// Compute solution
-////////////////////////////
-VectorCT x_1 = load_alpha1;
-VectorCT x_2 = load_alpha2;
-VectorCT x_3 = load_alpha3;
-
-
-// Turn the matrix into a linear operator
-MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_CE);
-
-// Fancy (but only) way to not have a preconditioner at all
-Richardson<VectorCT,VectorCT> preconditioner(1.0);
-
-// Construct the iterative solver
-RestartedGMResSolver<VectorCT> solver(
+ // Construct the iterative solver
+ RestartedGMResSolver<VectorCT> solver(
stiffnessOperator, // Operator to invert
preconditioner, // Preconditioner
1e-10, // Desired residual reduction factor
@@ -1661,301 +1606,297 @@ RestartedGMResSolver<VectorCT> solver(
500, // Maximum number of iterations
2); // Verbosity of the solver
-// Object storing some statistics about the solving process
-InverseOperatorResult statistics;
+ // Object storing some statistics about the solving process
+ InverseOperatorResult statistics;
+ // solve for different Correctors (alpha = 1,2,3)
+ solver.apply(x_1, load_alpha1, statistics);
+ solver.apply(x_2, load_alpha2, statistics);
+ solver.apply(x_3, load_alpha3, statistics);
-// solve for different Correctors (alpha = 1,2,3)
-solver.apply(x_1, load_alpha1, statistics);
-solver.apply(x_2, load_alpha2, statistics);
-solver.apply(x_3, load_alpha3, statistics);
-// Extract phi_alpha & M_alpha coefficients
-VectorCT phi_1, phi_2, phi_3;
-phi_1.resize(Basis_CE.size());
-phi_1 = 0;
-phi_2.resize(Basis_CE.size());
-phi_2 = 0;
-phi_3.resize(Basis_CE.size());
-phi_3 = 0;
+ // Extract phi_alpha & M_alpha coefficients
+ VectorCT phi_1, phi_2, phi_3;
+ phi_1.resize(Basis_CE.size());
+ phi_1 = 0;
+ phi_2.resize(Basis_CE.size());
+ phi_2 = 0;
+ phi_3.resize(Basis_CE.size());
+ phi_3 = 0;
-// VectorCT m_1, m_2, m_3;
-// m_1.resize(3);
-// m_1 = 0;
-// m_2.resize(3);
-// m_2 = 0;
-// m_3.resize(3);
-// m_3 = 0;
+ // VectorCT m_1, m_2, m_3;
+ // m_1.resize(3);
+ // m_1 = 0;
+ // m_2.resize(3);
+ // m_2 = 0;
+ // m_3.resize(3);
+ // m_3 = 0;
-FieldVector<double,3> m_1, m_2, m_3;
+ FieldVector<double,3> m_1, m_2, m_3;
-for(size_t i=0; i<Basis_CE.size();i++)
-{
+ for(size_t i=0; i<Basis_CE.size(); i++)
+ {
phi_1[i] = x_1[i];
phi_2[i] = x_2[i];
phi_3[i] = x_3[i];
-}
+ }
-for(size_t i=0; i<3; i++)
-{
+ for(size_t i=0; i<3; i++)
+ {
m_1[i] = x_1[phiOffset+i];
m_2[i] = x_2[phiOffset+i];
m_3[i] = x_3[phiOffset+i];
-}
+ }
-////////////////////////////////////////////////////////////////////////////
-// Make a discrete function from the FE basis and the coefficient vector
-////////////////////////////////////////////////////////////////////////////
-using SolutionRange = FieldVector<double,dim>;
-// using SolutionRange = double;
+ ////////////////////////////////////////////////////////////////////////////
+ // Make a discrete function from the FE basis and the coefficient vector
+ ////////////////////////////////////////////////////////////////////////////
+ using SolutionRange = FieldVector<double,dim>;
+ // using SolutionRange = double;
-// auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
-// Basis_CE,
-// Dune::Functions::HierarchicVectorWrapper<VectorCT, double>(phi_1));
+ // auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ // Basis_CE,
+ // Dune::Functions::HierarchicVectorWrapper<VectorCT, double>(phi_1));
-auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
- Basis_CE,
- phi_1);
+ auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ Basis_CE,
+ phi_1);
-auto correctorFunction_2 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
- Basis_CE,
- phi_2);
+ auto correctorFunction_2 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ Basis_CE,
+ phi_2);
-auto correctorFunction_3 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
- Basis_CE,
- phi_3);
+ auto correctorFunction_3 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ Basis_CE,
+ phi_3);
-// assemble M_alpha's (actually iota(M_alpha) )
+ // assemble M_alpha's (actually iota(M_alpha) )
-MatrixRT M_1(0), M_2(0), M_3(0);
+ MatrixRT M_1(0), M_2(0), M_3(0);
-for(size_t i=0; i<3; i++)
-{
+ for(size_t i=0; i<3; i++)
+ {
M_1 += m_1[i]*basisContainer[i];
M_2 += m_2[i]*basisContainer[i];
M_3 += m_3[i]*basisContainer[i];
-}
-
-std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
-printmatrix(std::cout, M_1, "Corrector-Matrix M_1", "--");
-printmatrix(std::cout, M_2, "Corrector-Matrix M_2", "--");
-printmatrix(std::cout, M_3, "Corrector-Matrix M_3", "--");
-
-
-
-
-auto localPhi_1 = localFunction(correctorFunction_1);
-auto localPhi_2 = localFunction(correctorFunction_2);
-auto localPhi_3 = localFunction(correctorFunction_3);
-
-std::cout << "check integralMean: " << std::endl;
-auto A = integralMean(Basis_CE,localPhi_1);
-for(size_t i=0; i<3; i++)
-{
-std::cout << "Integral-Mean : " << A[i] << std::endl;
-}
-
-if(substract_integralMean)
-{
-std::cout << " --- substracting integralMean --- " << std::endl;
-subtractIntegralMean(Basis_CE, localPhi_1 , phi_1);
-subtractIntegralMean(Basis_CE, localPhi_2 , phi_2);
-subtractIntegralMean(Basis_CE, localPhi_3 , phi_3);
-////////////////////////////////////////////////////////////////////////
-// CHECK INTEGRAL-MEAN:
-auto AVFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
- Basis_CE,
- phi_1);
-
-auto AVPhi_1 = localFunction(AVFunction_1);
-auto A = integralMean(Basis_CE, AVPhi_1);
-for(size_t i=0; i<3; i++)
-{
-std::cout << "Integral-Mean (CHECK) : " << A[i] << std::endl;
-}
-////////////////////////////////////////////////////
-}
-
-
-
-
-// // PRINT COEFFICIENTS
-// std::cout << "print coefficient-vector phi_1" << std::endl;
-// for(size_t i=0; i<Basis_CE.size();i++)
-// {
-// std::cout << phi_1[i] << std::endl;
-// }
-// std::cout << "print coefficient-vector phi_2" << std::endl;
-// for(size_t i=0; i<Basis_CE.size();i++)
-// {
-// std::cout << phi_2[i] << std::endl;
-// }
-// std::cout << "print coefficient-vector phi_3" << std::endl;
-// for(size_t i=0; i<Basis_CE.size();i++)
-// {
-// std::cout << phi_3[i] << std::endl;
-// }
-//
-
-
-
+ }
+ std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
+ printmatrix(std::cout, M_1, "Corrector-Matrix M_1", "--");
+ printmatrix(std::cout, M_2, "Corrector-Matrix M_2", "--");
+ printmatrix(std::cout, M_3, "Corrector-Matrix M_3", "--");
-// auto df = derivative(correctorFunction_1); // does not work :(
+ auto localPhi_1 = localFunction(correctorFunction_1);
+ auto localPhi_2 = localFunction(correctorFunction_2);
+ auto localPhi_3 = localFunction(correctorFunction_3);
+ std::cout << "check integralMean phi_1: " << std::endl;
+ auto A = integralMean(Basis_CE,localPhi_1);
+ for(size_t i=0; i<3; i++)
+ {
+ std::cout << "Integral-Mean : " << A[i] << std::endl;
+ }
+ if(substract_integralMean)
+ {
+ std::cout << " --- substracting integralMean --- " << std::endl;
+ subtractIntegralMean(Basis_CE, localPhi_1 , phi_1);
+ subtractIntegralMean(Basis_CE, localPhi_2 , phi_2);
+ subtractIntegralMean(Basis_CE, localPhi_3 , phi_3);
+ ////////////////////////////////////////////////////////////////////////
+ // CHECK INTEGRAL-MEAN:
+ auto AVFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ Basis_CE,
+ phi_1);
+
+ auto AVPhi_1 = localFunction(AVFunction_1);
+ auto A = integralMean(Basis_CE, AVPhi_1);
+ for(size_t i=0; i<3; i++)
+ {
+ std::cout << "Integral-Mean (CHECK) : " << A[i] << std::endl;
+ }
+ ////////////////////////////////////////////////////
+ }
+ // // PRINT Corrector-Basis-Coefficients
+ // std::cout << "print coefficient-vector phi_1" << std::endl;
+ // for(size_t i=0; i<Basis_CE.size();i++)
+ // {
+ // std::cout << phi_1[i] << std::endl;
+ // }
+ // std::cout << "print coefficient-vector phi_2" << std::endl;
+ // for(size_t i=0; i<Basis_CE.size();i++)
+ // {
+ // std::cout << phi_2[i] << std::endl;
+ // }
+ // std::cout << "print coefficient-vector phi_3" << std::endl;
+ // for(size_t i=0; i<Basis_CE.size();i++)
+ // {
+ // std::cout << phi_3[i] << std::endl;
+ // }
+ //
-const std::array<Func2Tensor, 3> x3MatrixBasis = {x3G_1, x3G_2, x3G_3};
+ // auto df = derivative(correctorFunction_1); // does not work :(
-const std::array<VectorCT, 3> coeffContainer = {x_1, x_2, x_3};
-const std::array<decltype(localPhi_1)*, 3> phiContainer= {&localPhi_1, &localPhi_2, &localPhi_3};
-const std::array<MatrixRT*, 3 > mContainer = {&M_1 , &M_2, & M_3};
+ /////////////////////////////////////////////////////////
+ // Create Containers for Basis-Matrices, Correctors and Coefficients
+ /////////////////////////////////////////////////////////
+ const std::array<Func2Tensor, 3> x3MatrixBasis = {x3G_1, x3G_2, x3G_3};
+ const std::array<VectorCT, 3> coeffContainer = {x_1, x_2, x_3};
+ const std::array<decltype(localPhi_1)*, 3> phiContainer= {&localPhi_1, &localPhi_2, &localPhi_3};
+ const std::array<MatrixRT*, 3 > mContainer = {&M_1 , &M_2, & M_3};
-// compute Q
-MatrixRT Q(0);
-VectorCT tmp1,tmp2;
+ /////////////////////////////////////////////////////////
+ // Compute effective quantities
+ /////////////////////////////////////////////////////////
+ MatrixRT Q(0);
-FieldVector<double,3> B_hat ;
+ VectorCT tmp1,tmp2;
+ FieldVector<double,3> B_hat ;
-for(size_t a = 0; a < 3; a++)
-{
+ // compute Q
+ 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) >
- auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
+ auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
- Q[a][b] = coeffContainer[a]*tmp1 + GGterm; // seems symmetric... check positiv definitness?
+ Q[a][b] = coeffContainer[a]*tmp1 + GGterm; // seems symmetric... check positiv definitness?
}
-}
-printmatrix(std::cout, Q, "Matrix Q", "--");
+ }
+ printmatrix(std::cout, Q, "Matrix Q", "--");
-for(size_t a = 0; a < 3; a++)
-{
+ // compute B_hat
+ 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;
+ std::cout <<"B_hat[i]: " << B_hat[a] << std::endl;
+ }
-}
-for(size_t i=0; i<3; i++)
-{
- std::cout <<"B_hat[i]: " << B_hat[i] << std::endl;
-}
+ // compute effective Prestrain
+ FieldVector<double, 3> Beff;
-// compute effective Prestrain
-FieldVector<double, 3> Beff;
+ Q.solve(Beff,B_hat);
-
-Q.solve(Beff,B_hat);
-
-std::cout << "Beff : " << std::endl;
-for(size_t i=0; i<3; i++)
-{
+ std::cout << "Beff : " << std::endl;
+ for(size_t i=0; i<3; i++)
+ {
std::cout <<"Beff[i]: " << Beff[i] << std::endl;
-}
-
-
-
-// compute q1,q2,q3
-
-FieldVector<double ,3> g1 = {1.0 , 0 , 0};
-FieldVector<double ,3> g2 = {0 , 1.0 , 0};
-FieldVector<double ,3> g3 = {0 , 0, 1.0};
-FieldVector<double ,3> tmp_1;
-FieldVector<double ,3> tmp_2;
-FieldVector<double ,3> tmp_3;
-// auto tmp = g1- B_hat;
-// std::cout<< tmp << std::endl;
-
-Q.mv(g1,tmp_1);
-auto q1_c = tmp_1 * g1;
-std::cout<< "q1_c: " << q1_c << std::endl;
-
-Q.mv(g2,tmp_2);
-auto q2_c = tmp_2 * g2;
-std::cout<< "q2_c: " << q2_c << std::endl;
-
-Q.mv(g3,tmp_3);
-auto q3_c = tmp_3 * g3;
-std::cout<< "q3_c: " << q3_c << std::endl;
-
-
-std::cout << "Gamma: " << gamma << 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 = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
-double b2 = -(theta/4.0)*mu2;
-double b3 = 0.0;
-double q1 = (mu1/6.0)*mu2/(theta*mu1+ (1.0- theta)*mu2);
-double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
-
-std::cout << " --- print analytic solutions --- " << std::endl;
-std::cout << "b1 : " << b1 << std::endl;
-std::cout << "b2 : " << b2 << std::endl;
-std::cout << "b3 : " << b3 << std::endl;
-std::cout << "q1 : " << q1 << std::endl;
-std::cout << "q2 : " << q2 << std::endl;
-std::cout << "q3 should be between q1 and q2" << std::endl;
-
-
-
-//////////////////////////////////////////////////////////////////////////////////////////////
-// Write result to VTK file
-// We need to subsample, because the dune-grid VTKWriter cannot natively display
-// second-order functions
-//////////////////////////////////////////////////////////////////////////////////////////////
+ }
-// SubsamplingVTKWriter<GridView> vtkWriter(
-// gridView,
-// refinementLevels(2));
- VTKWriter<GridView> vtkWriter(gridView_CE);
- vtkWriter.addVertexData(
- correctorFunction_1,
- VTK::FieldInfo("corrector phi_1", VTK::FieldInfo::Type::vector, dim));
- vtkWriter.addVertexData(
- correctorFunction_2,
- VTK::FieldInfo("corrector phi_2", VTK::FieldInfo::Type::vector, dim));
- vtkWriter.addVertexData(
- correctorFunction_3,
- VTK::FieldInfo("corrector phi_3", VTK::FieldInfo::Type::vector, dim));
-// vtkWriter.addVertexData(
-// solutionFunction,
-// VTK::FieldInfo("solution", VTK::FieldInfo::Type::scalar, 1));
- vtkWriter.write("CellProblem-result");
- std::cout << "wrote data to file: CellProblem-result" << std::endl; // better with string for output name..
+ // compute q1,q2,q3
+ FieldVector<double ,3> g1 = {1.0 , 0 , 0};
+ FieldVector<double ,3> g2 = {0 , 1.0 , 0};
+ FieldVector<double ,3> g3 = {0 , 0, 1.0};
+ FieldVector<double ,3> tmp_1;
+ FieldVector<double ,3> tmp_2;
+ FieldVector<double ,3> tmp_3;
+ // auto tmp = g1- B_hat;
+ // std::cout<< tmp << std::endl;
+
+ Q.mv(g1,tmp_1);
+ auto q1_c = tmp_1 * g1;
+ std::cout<< "q1_c: " << q1_c << std::endl;
+
+ Q.mv(g2,tmp_2);
+ auto q2_c = tmp_2 * g2;
+ std::cout<< "q2_c: " << q2_c << std::endl;
+
+ Q.mv(g3,tmp_3);
+ auto q3_c = tmp_3 * g3;
+ std::cout<< "q3_c: " << q3_c << std::endl;
+
+
+ std::cout << "Gamma: " << gamma << std::endl;
+
+ log << "q1: " << q1_c << std::endl;
+ log << "q2: " << q2_c << std::endl;
+ log << "q3: " << q3_c << std::endl;
+ log << "b1: " << Beff[1] << std::endl;
+ log << "b2: " << Beff[2] << std::endl;
+ log << "b3: " << Beff[3] << std::endl;
+ log << "b1_hat: " << B_hat[1] << std::endl;
+ log << "b2_hat: " << B_hat[2] << std::endl;
+ log << "b3_hat: " << B_hat[3] << 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 = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
+ double b2 = -(theta/4.0)*mu2;
+ double b3 = 0.0;
+ double q1 = (mu1/6.0)*mu2/(theta*mu1+ (1.0- theta)*mu2);
+ double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
+
+ std::cout << " --- print analytic solutions --- " << std::endl;
+ std::cout << "b1 : " << b1 << std::endl;
+ std::cout << "b2 : " << b2 << std::endl;
+ std::cout << "b3 : " << b3 << std::endl;
+ std::cout << "q1 : " << q1 << std::endl;
+ std::cout << "q2 : " << q2 << std::endl;
+ std::cout << "q3 should be between q1 and q2" << std::endl;
+
+
+
+ //////////////////////////////////////////////////////////////////////////////////////////////
+ // Write result to VTK file
+ //////////////////////////////////////////////////////////////////////////////////////////////
+
+ // SubsamplingVTKWriter<GridView> vtkWriter(
+ // gridView,
+ // refinementLevels(2));
+ VTKWriter<GridView> vtkWriter(gridView_CE);
+ vtkWriter.addVertexData(
+ correctorFunction_1,
+ VTK::FieldInfo("corrector phi_1", VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.addVertexData(
+ correctorFunction_2,
+ VTK::FieldInfo("corrector phi_2", VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.addVertexData(
+ correctorFunction_3,
+ VTK::FieldInfo("corrector phi_3", VTK::FieldInfo::Type::vector, dim));
+ // vtkWriter.addVertexData(
+ // solutionFunction,
+ // VTK::FieldInfo("solution", VTK::FieldInfo::Type::scalar, 1));
+ vtkWriter.write("CellProblem-result");
+ std::cout << "wrote data to file: CellProblem-result" << std::endl; // better with string for output name..
+
+ log.close();
}