diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 3e4f5f20bd47c6008534aa430fee9ac7f6e89274..d9b1d490bb1173b15f5b9338caefe833e7cb05af 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -15,7 +15,8 @@ cellDomain = 1
#############################################
#nElements_Cell = 20 20 20 # number elements in each direction (y1 y2 x3)
-nElements_Cell = 15 15 15
+#nElements_Cell = 100 100 2
+nElements_Cell = 4 4 4
width = 1.0
@@ -38,7 +39,7 @@ lambda1 = 0.0
theta = 0.3 # volume fraction
-prestraintype = 2
+prestrainType = 2
rho1 = 1.0
@@ -82,6 +83,9 @@ Func2Tensor B2_ = [this] (const Domain& x) { // ISOTROPIC PRESSURE
set_IntegralZero = true
+arbitraryLocalIndex = 7
+arbitraryElementNumber = 3
+
#############################################
# Solver Type
diff --git a/src/cellsolver.parset b/src/cellsolver.parset
deleted file mode 100644
index e743ae4a9ebf789a1fd01004c6e2a269369cf191..0000000000000000000000000000000000000000
--- a/src/cellsolver.parset
+++ /dev/null
@@ -1,134 +0,0 @@
-#1=yes 0=no
-
-#path for logfile
-#outputPath = "../clamped_cosserat_rod/Be1_t4s1"
-
-#############################################
-# EOC settings
-#############################################
-
-number_computations = 5
-nElements_Cell_Coarse = 1 1 1
-
-printSol = 1
-
-##
-refinement_level = 10
-order = 1
-
-
-#############################################
-# Rod grid parameters
-#############################################
-
-nElements_Cell = 100 1 1 # number elements (y x2 x3) 1:2:2 recommented, min 2 for y axis
-width = 2.0 # [0,1]xwidth*S
-len = 2.0 # Length of rod
-eps = 0.25 # Periodicity
-h = 0.05 # 3d rod thickness
-
-
-#############################################
-# Rod problem specifications
-#############################################
-
-#homotopy reference configuration:
-# A: u_ref satisfy bc
-# B: u_ref do not bc
-# C: u_ref satisfy bc but other homotopy
-#boundary condition:
-# a: isometric
-# b: extended
-# c: compressed
-# d: ring
-# e: init from strains
-#GFE order
-# 1: first order
-# 2: second order ..
-#Material no
-# A: stiff isotrop
-# B: soft isotrop
-# C: stiff lateral isotrop (?)
-# D: soft lateral isotrop (?)
-
-#
-# es extend shear
-# tb twist bend
-#
-# e1 not explicit in name
-#
-#
-
-# boundary types:
-# isometric (1)
-# extended (2)
-# compressed (3)
-# ring (4)
-
-# centerline implementations:
-# ellipse (3)
-
-testcase = Be
-
-U_ref = 1 0 0
-K_ref = 4 0 0
-
-U_init = 1 1 0
-K_init = 0 0 0
-
-#Qelas = 755e+7 755 1963
-#Qcurv = 1.23 1.23 0.94
-
-#Qelas = 4.4e+6 4.4e+6 4.4e+6 #stretch shear shear
-#Qcurv = 2.7e+6 10e+8 2.7e+6 #twist bend bend
-Qelas = 1963 755 755
-Qcurv = 0.94 1.23 1.23
-
-
-
-
-
-
-
-
-#############################################
-# Cosserat rod BVB solver parameters
-#############################################
-
-# Initial load increment
-loadIncrement = 0.01
-
-# Tolerance of the trust region solver
-tolerance = 1e-6
-
-# Max number of steps of the trust region solver
-maxTrustRegionSteps = 60
-
-# Initial trust-region radius
-initialTrustRegionRadius = 1
-
-# Number of multigrid iterations per trust-region step
-numIt = 200
-
-# Number of presmoothing steps (S=Solver to distingues between material parameter)
-nu1_S = 3
-
-# Number of postsmoothing steps
-nu2_S = 3
-
-# Number of coarse grid corrections
-mu_S = 1
-
-# Number of base solver iterations
-baseIt = 100
-
-# Tolerance of the multigrid solver
-mgTolerance = 1e-10
-
-# Tolerance of the base grid solver
-baseTolerance = 1e-8
-
-instrumented = no
-
-
-
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index e0f840452d6494c0e7cd04b9aa1e308ba29b1711..c14ea1cb7d8e8dbefe5465955f2ff40b987b313e 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -64,13 +64,6 @@ auto arbitraryComponentwiseIndices(const Basis& basis,
// (Local Approach -- works for non Lagrangian-Basis too)
// Input : elementNumber & localIdx
// Output : determine all Component-Indices that correspond to that basis-function
-
-// constexpr int dim = 3;
-
-// using GridView = typename Basis::GridView;
-// using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
-//
-
auto localView = basis.localView();
FieldVector<int,3> row;
@@ -84,14 +77,14 @@ auto arbitraryComponentwiseIndices(const Basis& basis,
for (int k = 0; k < 3; k++)
{
- auto rowLocal = localView.tree().child(k).localIndex(localIdx);
+ auto rowLocal = localView.tree().child(k).localIndex(localIdx); // TODO braucht hier Local to Leaf Idx Map
+ std::cout << "rowLocal:" << rowLocal << std::endl;
row[k] = localView.index(rowLocal);
-// std::cout << "row[k]:" << row[k] << std::endl;
+ std::cout << "row[k]:" << row[k] << std::endl;
}
}
elementCounter++;
}
-
return row;
}
@@ -389,15 +382,14 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet&
// OccupationPattern:
// | phi*phi m*phi |
// | phi *m m*m |
- auto gridView = basis.gridView();
-
+
auto localView = basis.localView();
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(basis.gridView()))
{
localView.bind(element);
for (size_t i=0; i<localView.size(); i++)
@@ -408,12 +400,10 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet&
{
// 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..
}
}
-
for (size_t i=0; i<localView.size(); i++)
{
auto row = localView.index(i);
@@ -421,58 +411,44 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet&
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<3; i++ )
for (size_t j=0; j<3; j++ )
{
nb.add(phiOffset+i,phiOffset+j); // m*m part of StiffnessMatrix
}
-
}
- //////////////////////////////////////////////////////////////////"localIndex is larger than #shapeFunctions",
+ //////////////////////////////////////////////////////////////////
// setOneBaseFunctionToZero
//////////////////////////////////////////////////////////////////
if(parameterSet.template get<bool>("set_IntegralZero", true)){
FieldVector<int,3> row;
- unsigned int elementCounter = 1;
unsigned int arbitraryLocalIndex = parameterSet.template get<unsigned int>("arbitraryLocalIndex", 0);
unsigned int arbitraryElementNumber = parameterSet.template get<unsigned int>("arbitraryElementNumber", 0);
-
- assert(arbitraryLocalIndex < localView.size() );
- assert(arbitraryElementNumber < basis.gridView().size(0));
+
+
+
+// assert(arbitraryLocalIndex < localView.size() ); // "localIndex is larger than #shapeFunctions" TODO ERROR
+
+
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // macht keinen Unterschied ob hier k oder 0..
+ const auto nSf = localFiniteElement.localBasis().size();
+ assert(arbitraryLocalIndex < nSf );
+
+ std::cout << "arbitraryElementNumber:" << arbitraryElementNumber << std::endl;
+ std::cout << "basis.gridView().size(0:" << basis.gridView().size(0) << std::endl;
+ assert(arbitraryElementNumber < basis.gridView().size(0)); // "arbitraryElementNumber larger than total Number of Elements"
//Determine 3 global indices (one for each component of an arbitrary local FE-function)
row = arbitraryComponentwiseIndices(basis,arbitraryElementNumber,arbitraryLocalIndex);
+ printvector(std::cout, row, "row" , "--");
-// for (const auto& element : elements(basis.gridView()))
-// {
-// localView.bind(element);
-//
-// if(elementCounter == arbitraryElementNumber) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
-// {
-// for (int k = 0; k < 3; k++)
-// {
-// auto rowLocal = localView.tree().child(k).localIndex(arbitraryLocalIndex);
-// row[k] = localView.index(rowLocal);
-// }
-// }
-// // "global assembly"
-// for (int k = 0; k<3; k++)
-// nb.add(row[k],row[k]);
-//
-// elementCounter++;
-// }
- // "global assembly"
for (int k = 0; k<3; k++)
- nb.add(row[k],row[k]);
+ nb.add(row[k],row[k]);
}
-
-
}
@@ -499,7 +475,7 @@ void computeElementLoadVector( const LocalView& localView,
constexpr int dim = Element::dimension;
constexpr int nCompo = dim;
- using VectorRT = FieldVector< double, nCompo>;
+// using VectorRT = FieldVector< double, nCompo>;
using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
// Set of shape functions for a single element
@@ -536,10 +512,11 @@ void computeElementLoadVector( const LocalView& localView,
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
const double integrationElement = geometry.integrationElement(quadPos);
- std::vector<FieldMatrix<double,1,nCompo> > referenceGradients;
+ std::vector<FieldMatrix<double,1,dim> > referenceGradients;
+
localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
- std::vector< VectorRT > gradients(referenceGradients.size());
+ std::vector< FieldVector< double, nCompo> > gradients(referenceGradients.size()); //nComo right? TODO
for (size_t i=0; i< gradients.size(); i++)
jacobian.mv(referenceGradients[i][0], gradients[i]);
@@ -613,171 +590,6 @@ void computeElementLoadVector( const LocalView& localView,
-///////////////////////////////////////////////// TEST ///////////////////////////////////////
-
-/*
-
-// Compute the source term for a single element
-// < L (sym[D_gamma*nabla phi_i] + M_i ), x_3G_alpha >
-template<class LocalView, class LocalFunction1, class LocalFunction2, class Vector, class Force>
-void computeElementLoadTEST( const LocalView& localView,
- LocalFunction1& mu,
- LocalFunction2& lambda,
- const double gamma,
- Vector& elementRhs,
- const Force& forceTerm
- )
-{
- // | 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 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();
-
-
- elementRhs.resize(localView.size() +3);
- elementRhs = 0;
-
-
- // 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.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
- std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
-
-
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1); // für simplex
- const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-
-
- for (const auto& quadPoint : quad)
- {
-
-
- const FieldVector<double,dim>& quadPos = quadPoint.position();
- const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
- const double integrationElement = geometry.integrationElement(quadPos);
-
- 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 strain1 = forceTerm(quadPos);
-
- // "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]);
-
-
-
-
-
-
- // St.Venant-Kirchhoff stress
- MatrixRT stressV(0);
- stressV.axpy(2*mu(quadPos), strain1); //this += 2mu *strainU
-
- double trace = 0.0;
- for (int ii=0; ii<nCompo; ii++)
- trace += strain1[ii][ii];
-
- for (int ii=0; ii<nCompo; ii++)
- stressV[ii][ii] += lambda(quadPos) * trace;
-
-
- // Local energy density: stress times strain
- double energyDensity = 0.0;
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
- energyDensity += stressV[ii][jj] * strainV[ii][jj];
-
-
- size_t row = localView.tree().child(k).localIndex(i);
- elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
- }
-
-
- // "f*m"-part
- for (size_t m=0; m<3; m++)
- {
-
- // St.Venant-Kirchhoff stress
- MatrixRT stressG(0);
- stressG.axpy(2*mu(quadPos), strain1); //this += 2mu *strainU
-
- double traceG = 0;
- for (int ii=0; ii<nCompo; ii++)
- traceG += strain1[ii][ii];
-
- 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] * basisContainer[m][ii][jj];
-
- elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
- }
-
- }
-}
-
-
-*/
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-////////////////////////////////////////////////////////////////////////////////////////////////
-
-
-
-
template<class Basis, class Matrix, class LocalFunction1, class LocalFunction2, class ParameterSet>
void assembleCellStiffness(const Basis& basis,
LocalFunction1& muLocal,
@@ -794,23 +606,15 @@ void assembleCellStiffness(const Basis& basis,
occupationPattern.exportIdx(matrix);
matrix = 0;
-
-
- std::fstream localdebugLog;
- localdebugLog.open("../../outputs/debugLog.txt",std::ios::out);
-
-
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);
-
const int localPhiOffset = localView.size();
// std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
@@ -825,24 +629,20 @@ void assembleCellStiffness(const Basis& basis,
// std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
-
//////////////////////////////////////////////////////////////////////////////
// GLOBAL STIFFNES ASSEMBLY
//////////////////////////////////////////////////////////////////////////////
for (size_t i=0; i<localPhiOffset; i++)
{
-
auto row = localView.index(i);
for (size_t j=0; j<localPhiOffset; j++ )
{
-
auto col = localView.index(j);
matrix[row][col] += elementMatrix[i][j];
}
}
-
for (size_t i=0; i<localPhiOffset; i++)
for(size_t m=0; m<3; m++)
{
@@ -852,7 +652,6 @@ void assembleCellStiffness(const Basis& basis,
matrix[phiOffset+m][row] += elementMatrix[localPhiOffset+m][i];
}
-
for (size_t m=0; m<3; m++ )
for (size_t n=0; n<3; n++ )
{
@@ -863,7 +662,6 @@ void assembleCellStiffness(const Basis& basis,
}
-
}
@@ -876,7 +674,6 @@ void assembleCellLoad(const Basis& basis,
const Force& forceTerm
)
{
-
// std::cout << "assemble load vector." << std::endl;
b.resize(basis.size()+3);
@@ -934,11 +731,6 @@ void assembleCellLoad(const Basis& basis,
-
-
-
-
-
template<class Basis, class LocalFunction1, class LocalFunction2, class MatrixFunction>
auto energy(const Basis& basis,
LocalFunction1& muLocal,
@@ -1037,35 +829,22 @@ void setOneBaseFunctionToZero(const Basis& basis,
Vector& load_alpha2,
Vector& load_alpha3,
ParameterSet& parameterSet
-// FieldVector<int,3>& indexVector // TODO oder input : arbitraryLocalIndex und dann hier function computeArbitraryIndex aufrufen
)
{
std::cout << "Setting one Basis-function to zero" << std::endl;
constexpr int dim = 3;
-// auto localView = basis.localView();
FieldVector<int,3> row;
-
unsigned int arbitraryLocalIndex = parameterSet.template get<unsigned int>("arbitraryLocalIndex", 0);
unsigned int arbitraryElementNumber = parameterSet.template get<unsigned int>("arbitraryElementNumber", 0);
//Determine 3 global indices (one for each component of an arbitrary local FE-function)
row = arbitraryComponentwiseIndices(basis,arbitraryElementNumber,arbitraryLocalIndex);
- // use first element:
-// auto it = basis.gridView().template begin<0>();
-// localView.bind(*it);
-//
-// int arbitraryIndex = 0;
-// FieldVector<int,3> row; //fill with Indices..
-
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_alpha1[row[k]] = 0.0;
load_alpha2[row[k]] = 0.0;
load_alpha3[row[k]] = 0.0;
auto cIt = stiffnessMatrix[row[k]].begin();
@@ -1073,7 +852,6 @@ void setOneBaseFunctionToZero(const Basis& basis,
for (; cIt!=cEndIt; ++cIt)
*cIt = (cIt.index()==row[k]) ? 1.0 : 0.0;
}
-
}
@@ -1089,20 +867,9 @@ auto childToIndexMap(const Basis& basis,
{
// Input : child/component
// Output : determine all Indices that correspond to that component
- constexpr int dim = 3;
-
-
-
- using GridView = typename Basis::GridView;
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
-
-
- using FuncScalar = std::function< FieldVector< double, dim>(const Domain&) >;
auto localView = basis.localView();
-
-
std::vector<int> r = { };
// for (int n: r)
// std::cout << n << ","<< std::endl;
@@ -1111,22 +878,16 @@ auto childToIndexMap(const Basis& basis,
// (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();
-
+ const auto nSf = localFiniteElement.localBasis().size(); //
for(size_t j=0; j<nSf; j++)
{
auto Localidx = localView.tree().child(k).localIndex(j); // local indices
r.push_back(localView.index(Localidx)); // global indices
}
-
-
-
}
-
// Delete duplicate elements
// first remove consecutive (adjacent) duplicates
auto last = std::unique(r.begin(), r.end());
@@ -1136,7 +897,6 @@ auto childToIndexMap(const Basis& basis,
last = std::unique(r.begin(), r.end());
r.erase(last, r.end());
-
return r;
}
@@ -1152,85 +912,36 @@ auto integralMean(const Basis& basis,
)
{
// computes integral mean of given LocalFunction
-
-
-
-
-
- using GridView = typename Basis::GridView;
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- constexpr int dim = 3;
- using FuncScalar = std::function< FieldVector< double, dim>(const Domain&) >;
-
+ constexpr int dim = 3;
auto GVFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis,coeffVector);
auto localfun = localFunction(GVFunction);
-
-
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
+ // Compute Area integral & integral of FE-function
for(const auto& element : elements(basis.gridView()))
{
-
localView.bind(element);
+ localfun.bind(element);
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
- for (const auto& quadPoint : quad)
+ for(const auto& quadPoint : quad)
{
const double integrationElement = element.geometry().integrationElement(quadPoint.position());
area += quadPoint.weight() * integrationElement;
- }
- }
- // std::cout << "Domain-Area: " << area << std::endl;
-
-
- for(const auto& element : elements(basis.gridView()))
- {
-
- localView.bind(element);
- localfun.bind(element);
- const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- int order = 2*(dim*localFiniteElement.localBasis().order()-1);
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), order);
-
- for (const auto& quadPoint : quad)
- {
-
+
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]);
-
- auto value = localfun(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;
-
-
- // r += tmp * quadPoint.weight() * integrationElement;
-
+ r += localfun(quadPos) * quadPoint.weight() * integrationElement;
}
}
-
- // return r/area;
+ // std::cout << "Domain-Area: " << area << std::endl;
return (1.0/area) * r ;
}
@@ -1246,7 +957,7 @@ auto subtractIntegralMean(const Basis& basis,
Vector& coeffVector
)
{
- // subtract correct Integral mean from each associated component function
+ // Substract correct Integral mean from each associated component function
// auto IM = integralMean(basis, fun);
auto IM = integralMean(basis, coeffVector);
@@ -1682,14 +1393,6 @@ for (const auto& e : elements(basis_.gridView()))
-
-
-
-
-
-
-
-
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -1698,24 +1401,11 @@ for (const auto& e : elements(basis_.gridView()))
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-
-
-
-
-
-
-
-
-
-
int main(int argc, char *argv[])
{
MPIHelper::instance(argc, argv);
-
-
ParameterTree parameterSet;
if (argc < 2)
ParameterTreeParser::readINITree("../../inputs/cellsolver.parset", parameterSet);
@@ -1773,7 +1463,6 @@ int main(int argc, char *argv[])
// Generate the grid
///////////////////////////////////
-
//Corrector Problem Domain:
unsigned int cellDomain = parameterSet.get<unsigned int>("cellDomain", 1);
// (shifted) Domain (-1/2,1/2)^3
@@ -1788,7 +1477,6 @@ int main(int argc, char *argv[])
}
std::array<int, dim> nElements = parameterSet.get<std::array<int,dim>>("nElements_Cell", {10, 10, 10});
-
log << "Number of Elements in each direction: " << nElements << std::endl;
using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
@@ -1796,20 +1484,15 @@ int main(int argc, char *argv[])
using GridView = CellGridType::LeafGridView;
const GridView gridView_CE = grid_CE.leafGridView();
-
-
- using ScalarRT = FieldVector< double, 1>;
- using VectorRT = FieldVector< double, nCompo>;
+// 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 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> >;
///////////////////////////////////
@@ -1820,6 +1503,7 @@ int main(int argc, char *argv[])
double mu2 = beta*mu1;
double lambda1 = parameterSet.get<double>("lambda1",0.0);;
double lambda2 = beta*lambda1;
+
///////////////////////////////////
// Create Lambda-Functions for material Parameters depending on microstructure
///////////////////////////////////
@@ -1840,8 +1524,7 @@ int main(int argc, char *argv[])
return lambda2;
};
-
- auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE); // TODO needed here?!
+ auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
auto muLocal = localFunction(muGridF);
auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
auto lambdaLocal = localFunction(lambdaGridF);
@@ -1852,9 +1535,6 @@ int main(int argc, char *argv[])
log << "lambda1: " << lambda1 <<"\nlambda2:" << lambda2 << std::endl;
-
-
-
///////////////////////////////////
// --- Choose Solver ---
// 1 : CG-Solver
@@ -1896,7 +1576,6 @@ int main(int argc, char *argv[])
auto perTest = periodicIndices.indexPairSet(); // TODO remove
-
auto Basis_CE = makeBasis(
gridView_CE,
power<dim>( //lagrange<1>(),
@@ -1912,21 +1591,6 @@ int main(int argc, char *argv[])
const int phiOffset = Basis_CE.size();
debugLog << "phiOffset: "<< phiOffset << std::endl;
-
-// 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;
-
-
-
/////////////////////////////////////////////////////////
// Data structures: Stiffness matrix and right hand side vector
/////////////////////////////////////////////////////////
@@ -1943,10 +1607,10 @@ int main(int argc, char *argv[])
substract_integralMean = true;
}
- debugLog << " set_oneBasisFunction_Zero: " << set_oneBasisFunction_Zero << std::endl;
- debugLog << " substract_integralMean: " << substract_integralMean << std::endl;
-
+ /////////////////////////////////////////////////////////
+ // Define "rhs"-functions
+ /////////////////////////////////////////////////////////
Func2Tensor x3G_1 = [] (const Domain& x) {
return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
@@ -1959,8 +1623,6 @@ int main(int argc, char *argv[])
return MatrixRT{{0.0, 1.0/sqrt(2)*x[2], 0.0}, {1.0/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
-
-
Func2Tensor x3G_1neg = [x3G_1] (const Domain& x) {
return -1.0*x3G_1(x);
};
@@ -1972,7 +1634,8 @@ int main(int argc, char *argv[])
Func2Tensor x3G_3neg = [x3G_3] (const Domain& x) {
return -1.0*x3G_3(x);
};
-
+
+
///////////////////////////////////////////////
// Basis for R_sym^{2x2}
@@ -1981,11 +1644,21 @@ int main(int argc, char *argv[])
MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
-// printmatrix(std::cout, basisContainer[0] , "G_1", "--");
- // printmatrix(std::cout, basisContainer[1] , "G_2", "--");
- // printmatrix(std::cout, basisContainer[2] , "´G_3", "--");
-// log << basisContainer[0] << std::endl;
+// printmatrix(std::cout, basisContainer[0] , "G_1", "--");
+// printmatrix(std::cout, basisContainer[1] , "G_2", "--");
+// printmatrix(std::cout, basisContainer[2] , "´G_3", "--");
+// log << basisContainer[0] << std::endl;
+
+
+ //////////////////////////////////////////////////////////////////////
+ // Determine global indices of components for arbitrary (local) index
+ //////////////////////////////////////////////////////////////////////
+ int arbitraryLocalIndex = parameterSet.get<unsigned int>("arbitraryLocalIndex", 0); // localIdx..assert Number < #ShapeFcts .. also input Element?
+ int arbitraryElementNumber = parameterSet.get<unsigned int>("arbitraryElementNumber", 0);
+ FieldVector<int,3> row;
+ row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLocalIndex);
+ printvector(std::cout, row, "row" , "--");
/////////////////////////////////////////////////////////
// Assemble the system
@@ -2003,28 +1676,13 @@ int main(int argc, char *argv[])
//////////////////////////////////////////////////////////////////////
// Determine global indices of components for arbitrary (local) index
//////////////////////////////////////////////////////////////////////
- auto localView = Basis_CE.localView();
-
- int arbitraryLocalIndex = parameterSet.get<unsigned int>("arbitraryLocalIndex", 0); // localIdx..assert Number < #ShapeFcts .. also input Element?
- int arbitraryElementNumber = parameterSet.get<unsigned int>("arbitraryElementNumber", 0);
-
- FieldVector<int,3> row; //fill with Indices..
- row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLocalIndex);
-// printvector(std::cout, row, "row" , "--");
-
-// // use first element:
-// auto it = Basis_CE.gridView().template begin<0>();
-// localView.bind(*it);
+// int arbitraryLocalIndex = parameterSet.get<unsigned int>("arbitraryLocalIndex", 0); // localIdx..assert Number < #ShapeFcts .. also input Element?
+// int arbitraryElementNumber = parameterSet.get<unsigned int>("arbitraryElementNumber", 0);
+//
+// FieldVector<int,3> row;
+// row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLocalIndex);
+// printvector(std::cout, row, "row" , "--");
//
-// for (int k = 0; k < dim; k++)
-// {
-// auto rowLocal = localView.tree().child(k).localIndex(arbitraryLocalIndex);
-// row[k] = localView.index(rowLocal);
-// std::cout << "row[k]:" << row[k] << std::endl;
-// }
- ///////////////////////////////////////////////////////////
-
-