From 53dcb10d447d593ddf28369ce190b045a4cdafd9 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Mon, 13 Sep 2021 09:59:36 +0200
Subject: [PATCH] deleted unnecessary Files
---
dune/microstructure/cellassembler.hh | 746 ------------------
dune/microstructure/cellsolver.hh | 272 -------
dune/microstructure/effectivemodelcomputer.hh | 228 ------
dune/microstructure/microstructuredrodgrid.hh | 115 ---
dune/microstructure/testclass.hh | 68 --
5 files changed, 1429 deletions(-)
delete mode 100644 dune/microstructure/cellassembler.hh
delete mode 100644 dune/microstructure/cellsolver.hh
delete mode 100644 dune/microstructure/effectivemodelcomputer.hh
delete mode 100644 dune/microstructure/microstructuredrodgrid.hh
delete mode 100644 dune/microstructure/testclass.hh
diff --git a/dune/microstructure/cellassembler.hh b/dune/microstructure/cellassembler.hh
deleted file mode 100644
index 15d76355..00000000
--- a/dune/microstructure/cellassembler.hh
+++ /dev/null
@@ -1,746 +0,0 @@
-#ifndef DUNE_REDUCEDROD_CELLASSEMBLER_HH
-#define DUNE_REDUCEDROD_CELLASSEMBLER_HH
-
-#include <dune/common/parametertree.hh>
-
-#include <dune/istl/matrixindexset.hh>
-
-#include <dune/functions/functionspacebases/interpolate.hh>
-#include <dune/functions/gridfunctions/gridviewfunction.hh>
-#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
-
-using namespace Dune;
-using namespace Functions;
-
-
-template <class Basis> //, class LocalScalar, class Local2Tensor> // LocalFunction derived from basis?
-class CellAssembler {
-
-public:
- static const int nCompo = 3; //GridView::dimensionworld;
- static const int dim = 3; // typename Basis::GridView::dimension
- // enum {dim=Basis::LocalView::GridView::dimension};
-
- using GridView = typename Basis::GridView;
-
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
-
- using ScalarRT = FieldVector< double, 1>;
- using VectorRT = FieldVector< double, nCompo>;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
-
- 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 ElementMatrixCT = Matrix<FieldMatrix<double,1,1> >;
-
-
-protected:
-//private:
-
- const Basis& basis_; // Basis kopieren?
-
- std::fstream& log_;
- const ParameterTree& parameterSet_;
-
- const FuncScalar mu_; //ref? const?
- const FuncScalar lambda_;
- double gamma_;
-
- std::map<int,int> perIndexMap_;
- int sizeBasisTorus_;
-
- // kann fest im objekt enthalten sein
- /*
- Grid2Tensor E_a_;
- Grid2Tensor E_K1_;
- Grid2Tensor E_K2_;
- Grid2Tensor E_K3_; */
-
-
-public:
-
- // Constructor
- CellAssembler( const Basis& basis,
- const FuncScalar& mu,
- const FuncScalar& lambda,
- double gamma,
- std::fstream& log,
- const ParameterTree& parameterSet)
- : basis_(basis),
- mu_(mu),
- lambda_(lambda),
- gamma_(gamma),
- log_(log),
- parameterSet_(parameterSet)
- {
- //periodic FE index mapper
- sizeBasisTorus_ = assemblePeriodicIndexMap();
- }
-
- // getter
-
- //Basis getBasis() const {return basis_;}
- const Basis* getBasis() {return &basis_;}
-
- std::map<int,int> getIndexMapPeriodicBoundary(){return perIndexMap_;}
-
- ParameterTree getParameterSet() const {return parameterSet_;}
-
- std::fstream* getLog(){return &log_;}
-
- int getSizeTorusElements(){return sizeBasisTorus_;}
-
- // setter
-
- template <class LocalScalar1, class LocalScalar2>
- void computeLocalStiffnessMatrix(const typename Basis::LocalView& localView, const LocalScalar1& mu, const LocalScalar2& lambda, ElementMatrixCT& elementMatrix)
- {// Get the grid element from the local FE basis view
- auto element = localView.element();
- auto geometry = element.geometry();
-
- const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- const auto nSf = localFiniteElement.localBasis().size();
-
- elementMatrix.setSize(localView.size(), localView.size());
- elementMatrix = 0;
- int order = 2*(dim*localFiniteElement.localBasis().order()-1); // für simplex
- const QuadratureRule< double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), order);
-
- for (const auto& quadPoint : quad){
- const Domain& 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]);
-
- for (size_t i=0; i < nSf; i++)
- for (size_t j=0; j < nSf; j++ )
- for (size_t k=0; k < nCompo; k++)
- for (size_t l=0; l< nCompo; l++)
- {
- // Deformation gradient of the ansatz basis function, Domain and Range components are swaped
- MatrixRT defGradientU(0);
- defGradientU[k][0] = gradients[i][0]; // Y
- defGradientU[k][1] = 1.0/gamma_*gradients[i][1]; //X2
- defGradientU[k][2] = 1.0/gamma_*gradients[i][2]; //X3
-
- // Deformation gradient of the test basis function, Domain and Range components are swaped
- MatrixRT defGradientV(0);
- defGradientV[l][0] = gradients[j][0]; // Y
- defGradientV[l][1] = 1.0/gamma_*gradients[j][1]; //X2
- defGradientV[l][2] = 1.0/gamma_*gradients[j][2]; //X3
-
- // 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]); // same ?
- }
-
- // 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;
-
- // 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];
-
- size_t row = localView.tree().child(k).localIndex(i);
- size_t col = localView.tree().child(l).localIndex(j);
- elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
- }
- }
- }
-
-template<class LocalScalar1, class LocalScalar2, class Local2Tensor>
-void computeLocalLoad(const typename Basis::LocalView& localView, const LocalScalar1& mu, const LocalScalar2& lambda, const Local2Tensor& strainE, VectorCT& elementRhs)
-{// Get the grid element from the local FE basis view
-
- elementRhs.resize(localView.size());
- elementRhs = 0;
- const auto element = localView.element();
- const auto geometry = element.geometry();
-
- const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- const auto nSf = localFiniteElement.localBasis().size();
-
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1); // für simplex
- //log << "orderQR: " << orderQR << std::endl;
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
-
- for (const auto& quadPoint : quad){
- const Domain& 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]);
-
- for (size_t i=0; i < nSf; i++)
- for (size_t k=0; k < nCompo; k++)
- {
- // Deformation gradient of the ansatz basis function, no: Domain and Range components are swaped
- MatrixRT defGradientU(0);
- defGradientU[k][0] = gradients[i][0]; // Y
- defGradientU[k][1] = 1.0/gamma_*gradients[i][1]; //X2
- defGradientU[k][2] = 1.0/gamma_*gradients[i][2]; //X3
-
- //log_ << defGradientU << std::endl;
- // Elastic Strain
- MatrixRT strainU, strainV; //strainV never used?
- 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]);
-
- // 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;
- //log_ << "stressU:\n" << stressU << std::endl;
- // 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] * strainE(quadPos)[ii][jj];
- //log_ << "strainE:\n" << strainE(quadPoint.position()) << std::endl;
-
- size_t row = localView.tree().child(k).localIndex(i);
- elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
- //log_ << "\nk: " << k << "\ni: " << i << "\nrow: " << row << std::endl;
-
- //log << "energyDensity quadPoint.weight() integrationElement: " << energyDensity << " " << quadPoint.weight() << " " << integrationElement << std::endl;
- } //end for k
- } // end for quadPoint
-
-}// end method
-
-
-
- void assembleStiffnessMatrix(MatrixCT& matrix)
- {
- std::cout << "assemble stiffnessmatrix begins.\n";
- matrix = 0;
- auto localView = basis_.localView();
-
- auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
- auto muLocal = localFunction(muGridF);
- auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
- auto lambdaLocal = localFunction(lambdaGridF);
-
- unsigned int printEachElement = parameterSet_.get<unsigned int>("print_each_element");
- unsigned int counter =1;
- for (const auto& e : elements(basis_.gridView()))
- {
- if (counter % printEachElement == 1 )
- std::cout << "assemble stiffness matrix - element: " << counter << std::endl;
- localView.bind(e);
- muLocal.bind(e);
- auto B_Term = prestrainImp
- lambdaLocal.bind(e);
-
- Matrix<FieldMatrix<double,1,1> > elementMatrix; // type ?
- computeLocalStiffnessMatrix(localView, muLocal, lambdaLocal, elementMatrix);
-
- for (size_t i=0; i<elementMatrix.N(); i++){
- auto row = localView.index(i);
- for (size_t j=0; j<elementMatrix.M(); j++ ){
- auto col = localView.index(j);
- matrix[perIndexMap_[row]][perIndexMap_[col]] += elementMatrix[i][j];
- }
- }
- counter += 1;
- }
- }
-
-
- void assembleLoadVector(const Func2Tensor& matrixFieldFunc, VectorCT& rhs)
- {
- std::cout << "assemble load vector.\n";
-
- //rhs.resize(basis_.size());
- rhs.resize(sizeBasisTorus_);
-
- auto localView = basis_.localView();
- auto loadGVF = makeGridViewFunction(matrixFieldFunc, basis_.gridView());
- auto loadFunctional = localFunction(loadGVF);
-
- auto muGVF = makeGridViewFunction(mu_, basis_.gridView());
- auto muLocal = localFunction(muGVF);
- auto lambdaGVF = makeGridViewFunction(lambda_, basis_.gridView());
- auto lambdaLocal = localFunction(lambdaGVF);
-
- unsigned int printEachElement = parameterSet_.get<unsigned int>("print_each_element");
- unsigned int counter =1;
- for (const auto& e : elements(basis_.gridView()))
- {
- if (counter % printEachElement == 1)
- std::cout << "assemble load vector: next element. " << counter << std::endl;
- // Bind the local FE basis view to the current element
- localView.bind(e);
- loadFunctional.bind(e);
- muLocal.bind(e);
- lambdaLocal.bind(e);
-
- VectorCT elementRhs;
- computeLocalLoad(localView, muLocal, lambdaLocal, loadFunctional, elementRhs);
-
- for (size_t i=0; i<elementRhs.size(); i++) {
- auto row = localView.index(i);
- rhs[perIndexMap_[row]] += elementRhs[i];
- /*log_ << "\ni: " << i
- << "\nrow: " << row
- << "\nperIndexMap_[row]: " << perIndexMap_[row]
- << "\nelementRhs[i]: " << elementRhs[i] << std::endl;*/
- }
- counter += 1;
- }
- }
-
-
- ScalarRT energy(const Func2Tensor& matrixFieldFunc1,
- const Func2Tensor& matrixFieldFunc2)
- {
- ScalarRT energy = 0.0;
-
- auto localView = basis_.localView();
-
- auto matrixField1GVF = makeGridViewFunction(matrixFieldFunc1, basis_.gridView());
- auto matrixField1 = localFunction(matrixField1GVF);
- auto matrixField2GVF = makeGridViewFunction(matrixFieldFunc2, basis_.gridView());
- auto matrixField2 = localFunction(matrixField2GVF);
-
- auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
- auto muLocal = localFunction(muGridF);
- auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
- auto lambdaLocal = localFunction(lambdaGridF);
-
- for (const auto& e : elements(basis_.gridView()))
- {
- localView.bind(e);
- matrixField1.bind(e);
- matrixField2.bind(e);
- muLocal.bind(e);
- lambdaLocal.bind(e);
-
- FieldVector<double,1> elementEnergy(0);
-
- 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);
-
- for (size_t pt=0; pt < quad.size(); pt++) {
- const Domain& quadPos = quad[pt].position();
- const double integrationElement = geometry.integrationElement(quadPos);
-
- auto strain1 = matrixField1(quadPos);
-
- // 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];
-
- for (int ii=0; ii<nCompo; ii++)
- stressU[ii][ii] += lambdaLocal(quadPos) * trace;
-
- auto strain2 = matrixField2(quadPos);
-
- //log_ << "strain2= \n";
- //for (int iii = 0; iii < 3; iii++)
- //log_ << strain2[iii] << std::endl;
-
- // 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];
-
- //elementConstantTerm += linElasTensor4(localMatrixField1(quadPos), localMatrixField2(quadPos), mu(quadPos), lambda(quadPos)) * quad[pt].weight() * integrationElement;
- elementEnergy += energyDensity * quad[pt].weight() * integrationElement;
- }
- energy += elementEnergy;
- }
- return energy;
- }
-
-
- int assemblePeriodicIndexMap() // see Sander-DuneBook p.313..determine positions of Lagrange-nodes
- {
- std::vector<double> verticesL;
- auto fL = [](const Domain& x){return FieldVector<double,1>{x[0]};};
- Functions::interpolate(basis_, verticesL, fL);
-
- std::vector<double> verticesX2;
- auto fAcrossFromX2 = [](const Domain& x){return FieldVector<double,1>{x[1]};};
- Functions::interpolate(basis_, verticesX2, fAcrossFromX2);
-
- std::vector<double> verticesX3;
- auto fAcrossFromX3 = [](const Domain& x){return FieldVector<double,1>{x[2]};};
- Functions::interpolate(basis_, verticesX3, fAcrossFromX3);
-
- int sizeBasisTorus = 0;
- for (int d=0; d < nCompo; d++)
- for (int i = d*basis_.size()/nCompo; i < (d+1)*verticesL.size()/nCompo; i++)
- if (verticesL[i]!=1.0){
- perIndexMap_[i]=sizeBasisTorus;
- sizeBasisTorus++;
- }else{
- auto iL = oppositeVertice(d, verticesL, verticesX2, verticesX3, i);
- perIndexMap_[i]=perIndexMap_[iL];
- }
- log_ << "sizeBasisTorus: " << sizeBasisTorus << std::endl;
- std::cout << "sizeBasisTorus: " << sizeBasisTorus << std::endl;
- if (parameterSet_.get<bool>("plot_periodic_index_mapper")) {
- log_ << "\nplot_periodic_index_mapper:\n ";
- for (int i = 0; i < perIndexMap_.size(); i++){
- log_ << i << " " << perIndexMap_[i] << std::endl;
- }
- }
-
- return sizeBasisTorus;
- }
-
-
- int oppositeVertice(const int d,
- const std::vector<double>& verticesL,
- const std::vector<double>& verticesAcrossFromX2 ,
- const std::vector<double>& verticesAcrossFromX3,
- int i)
- {
- int j=0;
- for (int iL = d*verticesL.size()/nCompo; iL < (d+1)*verticesL.size()/nCompo; iL++)
- if (verticesL[iL]==0.0 && verticesAcrossFromX2[i]==verticesAcrossFromX2[iL] && verticesAcrossFromX3[i]==verticesAcrossFromX3[iL])
- j = iL;
- return j;
- // test nodeL[i]==1.0
- // test d < nCompo
- }
-
-
- void getOccupationPattern(MatrixIndexSet& occupationPattern)
- {
- // in flat interleaved ist occupationPattern ein großer block statt 9, . erste 3x3 block in occupationPattern wuerde lauten:
- // 1x1x 1x1y 1x1z
- // 1y1x 1y1y 1y1z
- // 1z1x 1z1y 1z1z
- // Machen aber lexicographic
- std::cout << "setting non zero entries for stiffness matrix begins.\n";
-
- occupationPattern.resize(sizeBasisTorus_,sizeBasisTorus_);
-
- auto localView = basis_.localView();
- for(const auto& e : elements(basis_.gridView()))
- {
- localView.bind(e);
- for (size_t i =0; i<localView.size(); i++) {
- auto iIdx = localView.index(i);
- for (size_t j=0; j<localView.size(); j++) {
- auto jIdx = localView.index(j);
- occupationPattern.add(perIndexMap_[iIdx],perIndexMap_[jIdx]);
- }
- }
- }
-
- auto nE = basis_.size();
- auto nNodeBasis = nE/nCompo;
-
- if(parameterSet_.get<bool>("set_one_base_function_to_0" )){ // warum wird dieser Teil benötigt?
- for (int d=0; d < nCompo; d++)
- for (int j = 0; j < nNodeBasis; j++)
- occupationPattern.add(perIndexMap_[0 + d*nNodeBasis], perIndexMap_[0 + d*nNodeBasis]); // j ????!!
- }
-
- if(parameterSet_.get<bool>("set_integral_0")){
- for (int d=0; d < nCompo; d++)
- for (int j = 0; j < nNodeBasis; j++)
- occupationPattern.add(perIndexMap_[0 + d*nNodeBasis], perIndexMap_[j + d*nNodeBasis]);
- }
-
- if(parameterSet_.get<bool>("set_first_moment_0"))
- for (int j = 0; j < nE; j++)
- occupationPattern.add(perIndexMap_[nNodeBasis] + 1, perIndexMap_[j]);
-
- std::cout << "finished occupation pattern\n";
- }
-
-
- void setOneBaseFunctionToZero (MatrixCT& stiffnessMatrix,
- VectorCT& load_a,
- VectorCT& load_K1, VectorCT& load_K2, VectorCT& load_K3)
- {
- if (parameterSet_.get<bool>("set_one_base_function_to_0")){
- log_ << "\nassembling: one base function to zero\n";
- std::cout << "\nassembling: one base function to zero\n";
- auto nE = basis_.size();
- auto nENode = nE/nCompo;
-
- //arbitrary index < nENode
- int arbitraryIndex = 0;
- for (int d = 0; d < dim; d++){
- auto row = perIndexMap_[arbitraryIndex + d * nENode]; // in jeder Komponente auf 0 gesetzt
- load_a[row] = 0.0;
- load_K1[row] = 0.0;
- load_K2[row] = 0.0;
- load_K3[row] = 0.0;
- auto cIt = stiffnessMatrix[row].begin();
- auto cEndIt = stiffnessMatrix[row].end();
- for (; cIt!=cEndIt; ++cIt)
- *cIt = (cIt.index()==row) ? 1.0 : 0.0;
- }
- }
-
-
- }
-
-
-
- void setIntegralZero (MatrixCT& stiffnessMatrix,
- VectorCT& load_a,
- VectorCT& load_K1, VectorCT& load_K2, VectorCT& load_K3)
- {
-
- //
- // 0. moment invariance, 3 dofs
- //
- if (parameterSet_.get<bool>("set_integral_0")){
- log_ << "\nassembling: integral to zero\n";
- std::cout << "\nassembling: integral to zero\n";
- auto n = basis_.size();
- auto nNodeBasis = n/nCompo;
- auto localView = basis_.localView();
-
- int arbitraryIndex = 0;
- for (int d = 0; d < dim; d++) {
- auto row = perIndexMap_[arbitraryIndex + d*nNodeBasis];
- stiffnessMatrix[row] = 0;
- load_a[row] = 0.0;
- load_K1[row] = 0.0;
- load_K2[row] = 0.0;
- load_K3[row] = 0.0;
- }
-
- unsigned int elementCounter = 1;
- for (const auto& e : elements(basis_.gridView()))
- {
- if (elementCounter % 50 == 1) // ?????
- std::cout << "integral = zero - treatment - element: " << elementCounter << std::endl;
- localView.bind(e);
-
- BlockVector<FieldVector<double,1> > elementColumns;
- elementColumns.resize(localView.size());
- 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(e.type(), order);
-
- for (size_t pt=0; pt < quad.size(); pt++) {
-
- const FieldVector<double,dim>& quadPos = quad[pt].position();
- const double integrationElement = e.geometry().integrationElement(quadPos);
-
- std::vector<FieldVector<double,1> > shapeFunctionValues;
- localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
-
- for (size_t k=0; k<nSf; k++)
- for (size_t d=0; d<nCompo; d++)
- elementColumns[ k + d*nSf] += shapeFunctionValues[k] * quad[pt].weight() * integrationElement;
- }
-
- for (int d = 0; d < dim; d++) {
- auto row = perIndexMap_[arbitraryIndex + d*nNodeBasis];
- for (size_t j=0; j<nSf; j++){
- auto col = perIndexMap_[ (localView.index(j+ d*nSf)) ];
- stiffnessMatrix[row][col] += elementColumns[j+d*nSf];
- }
- }
- elementCounter ++;
- }//end for each element
- }
-
- }
-
-
- void setFirstMomentZero (MatrixCT& stiffnessMatrix,
- VectorCT& load_a,
- VectorCT& load_K1, VectorCT& load_K2, VectorCT& load_K3)
- {
- if (parameterSet_.get<bool>("set_first_moment_0")){
- log_<< "\nassembling: first moment = zero\n";
- std::cout<< "\nassembling: first moment = zero\n";
-
- auto zerox2x3Term = [] (const Domain& z) { return Domain({0, z[1], z[2]}); }; //Range components = Domain compo Y X2 X3 // ???
- auto zerox2x3 = Functions::makeGridViewFunction(zerox2x3Term, basis_.gridView());
-
- auto n = basis_.size();
- auto nNodeBasis = n/nCompo;
- auto localView = basis_.localView();
- auto zerox2x3Local = localFunction(zerox2x3);
-
- //arbitrary index < nENode
- int arbitraryIndex = nNodeBasis; //setOneBaseFunctionTo0 takes index= 0, we hav to modify rows which affect 2. or 3. component
- auto row = perIndexMap_[arbitraryIndex] + 1; // ???
- stiffnessMatrix[row] = 0;
- load_a[row] = 0;
- load_K1[row] = 0;
- load_K2[row] = 0;
- load_K3[row] = 0;
-
- unsigned int elementCounter = 1;
- for (const auto& e : elements(basis_.gridView()))
- {
- if (elementCounter % 50 == 1) // ????
- std::cout << "first moment = zero - element: " << elementCounter << std::endl;
-
- localView.bind(e);
- zerox2x3Local.bind(e);
-
- BlockVector<FieldVector<double,1> > elementColumns;
- elementColumns.resize(localView.size());
- 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(e.type(), order);
-
- for (size_t pt=0; pt < quad.size(); pt++) {
-
- const FieldVector<double,dim>& quadPos = quad[pt].position();
- const double integrationElement = e.geometry().integrationElement(quadPos);
- std::vector<FieldVector<double,1> > shapeFunctionValues;
- localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
-
- for (size_t k=0; k<nSf; k++)
- for (size_t d=0; d<dim; d++){
- FieldVector<double,3> sFVVector(0);
- sFVVector[d] = shapeFunctionValues[k];
- size_t ll = localView.tree().child(d).localIndex(k);
- elementColumns[ll] += (sFVVector * zerox2x3Local(quadPos)) * quad[pt].weight() * integrationElement;
- }
- }
-
- for (int j = 0; j < localView.size(); j++){
- auto col = perIndexMap_[(localView.index(j))];
- stiffnessMatrix[row][col] += elementColumns[j]; //arg1 = 1...secound vertice, arg2 = 0...first dimension
- }
- elementCounter ++;
- }//end for each element
-
- } // end if
- }
-
-
-
-
-
- // main method which uses the previous ones
-
- void assemble(MatrixCT& stiffnessMatrix,
- VectorCT& load_a, VectorCT& load_K1,
- VectorCT& load_K2, VectorCT& load_K3)
- {
-
- MatrixIndexSet occupationPattern;
- getOccupationPattern(occupationPattern);
- occupationPattern.exportIdx(stiffnessMatrix);
-
- assembleStiffnessMatrix(stiffnessMatrix);
-
- // Forces for cell problems = strain basis functions
- // Basis for the limiting Biot strain:
- // E_h = 1/h (sqrt(F^T F) - Id) -> E(U,K) + G
- // E(U,K) = sym( K * (0,x2,x3)^T + U ot e1)
-
- // lateral stretch strain E_U1 = e1 ot e1
- Func2Tensor E_a = [] (const Domain& z) {
- return MatrixRT{{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; };
- //{{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}};
-
- // x2 0 0 | x3 0 0 | 0 1/2 x3 -1/2 x2
- // 0 0 0 | 0 0 0 | 1/2 x3 0 0
- // 0 0 0 | 0 0 0 | -1/2 x2 0 0
- //
- // Ki.. basis of Skew(3) ... rod curvature
- //
- // twist in x2-x3 plane E_K1 = sym (e1 x (0,x2,x3)^T ot e1)
- Func2Tensor E_K1 = [] (const Domain& z) {
- return MatrixRT{{0.0 , -0.5*z[2], 0.5 * z[1]},
- {-0.5*z[2], 0.0 , 0.0 },
- {0.5*z[1], 0.0 , 0.0}}; };
- //{{0.0, 0.5*z[1],-0.5*z[0]}, {0.5*z[1],0.0, 0.0}, {-0.5*z[0], 0.0, 0.0}};
-
- // bend strain in x1-x2 plane E_K2 = sym (e2 x (0,x2,x3)^T ot e1)
- Func2Tensor E_K2 = [] (const Domain& z) {
- return MatrixRT{{z[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0}}; };
- //{{z[0], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
-
- // bend strain in x1-x3 plane E_K3 = sym (e3 x (0,x2,x3)^T ot e1)
- Func2Tensor E_K3 = [] (const Domain& z) {
- return MatrixRT{{-z[1], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; };
- //{{z[1], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
-
- assembleLoadVector(E_a, load_a);
- //log_ << "\n\n\n\n\n\n";
- assembleLoadVector(E_K1, load_K1);
- assembleLoadVector(E_K2, load_K2);
- assembleLoadVector(E_K3, load_K3);
-
- // invariants treatments
- setOneBaseFunctionToZero(stiffnessMatrix, load_a, load_K1, load_K2, load_K3);
- setFirstMomentZero(stiffnessMatrix, load_a, load_K1, load_K2, load_K3);
- setIntegralZero(stiffnessMatrix, load_a, load_K1, load_K2, load_K3);
- };
-
-}; // end class
-
-
-#endif
-
-
-
-// shear strain in x2 E_U2 = sym (e2 ot e1)
- /*
- F E_U2Term = [] (const Domain& z) {
- MatrixRT ret = {{0.0, 0.5, 0.0}, {0.5, 0.0, 0.0}, {0.0, 0.0, 0.0}};
- return ret; };
-
- // shear strain in x2 E_U2 = sym (e3 ot e1)
- F E_U3Term = [] (const Domain& z) {
- MatrixRT ret = {{0.0, 0.0, 0.5}, {0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}};
- return ret; };
- */
diff --git a/dune/microstructure/cellsolver.hh b/dune/microstructure/cellsolver.hh
deleted file mode 100644
index 340029eb..00000000
--- a/dune/microstructure/cellsolver.hh
+++ /dev/null
@@ -1,272 +0,0 @@
-#ifndef DUNE_REDUCEDROD_CELLSOLVER_HH
-#define DUNE_REDUCEDROD_CELLSOLVER_HH
-
-#include <dune/reducedrod/cellassembler.hh>
-#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
-
-#include <dune/istl/operators.hh>
-#include <dune/istl/preconditioners.hh>
-#include <dune/istl/solvers.hh>
-#include <dune/istl/spqr.hh>
-
-
-
-template <class Basis>
-class CellSolver {
-
-public:
-
- static const int nCompo = 3;
-
- using VectorRT = typename CellAssembler<Basis>::VectorRT;
-
- using VectorCT = typename CellAssembler<Basis>::VectorCT;
- using MatrixCT = typename CellAssembler<Basis>::MatrixCT;
-
- using HierarchicVectorView = Dune::Functions::HierarchicVectorWrapper<VectorCT, double>;
-
-protected:
-
- CellAssembler<Basis>& cellAssembler_;
-
- MatrixCT stiffnessMatrix_;
- VectorCT load_a_, load_K1_, load_K2_, load_K3_;
-
-public:
-
- // Constructor
- CellSolver(CellAssembler<Basis>& cellAssembler)
- : cellAssembler_(cellAssembler)
- {
- cellAssembler.assemble(stiffnessMatrix_,load_a_, load_K1_, load_K2_, load_K3_); //assemble hier rein?
- }
-
- //return ref ???
- CellAssembler<Basis> getCellAssembler(){return cellAssembler_;}
-
- MatrixCT * getStiffnessMatrix(){return &stiffnessMatrix_;}
- VectorCT * getLoad_a(){return &load_a_;}
- VectorCT * getLoad_K1(){return &load_K1_;}
- VectorCT * getLoad_K2(){return &load_K2_;}
- VectorCT * getLoad_K3(){return &load_K3_;}
-
- void solveCorrectorSystems(VectorCT& x_a, VectorCT& x_K1, VectorCT& x_K2, VectorCT& x_K3)
- {
-
- std::cout << "CellSolver: Solve linear systems begin.\n";
- ParameterTree parameterSet = cellAssembler_.getParameterSet();
- //auto logPointer = cellAssembler_.getLog();
- //auto & log = *logPointer;
- auto & log = *(cellAssembler_.getLog());
-
- int sizeBasisTorus_ = cellAssembler_.getSizeTorusElements();
- x_a.resize(sizeBasisTorus_);
- x_K1.resize(sizeBasisTorus_);
- x_K2.resize(sizeBasisTorus_);
- x_K3.resize(sizeBasisTorus_);
- x_a = 0;
- x_K1 = 0;
- x_K2 = 0;
- x_K3 = 0;
-
- if (parameterSet.get<int>("iterate_with_load") == 1){
- x_a = load_a_;
- x_K1 = load_K1_;
- x_K2 = load_K2_;
- x_K3 = load_K3_;
- }
-
- if (parameterSet.get<int>("useCG") == 1){
- std::cout << "We use CG solver.\n";
-
- MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_);
-
- std::cout << "MatrixAdapter op.\n";
-
- // Sequential incomplete LU decomposition as the preconditioner
- SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_,1.0);
- int iter = parameterSet.get<double>("iterations_CG", 999);
- // Preconditioned conjugate-gradient solver
- CGSolver<VectorCT> solver(op,
- ilu0, //NULL,
- 1e-8, // desired residual reduction factor
- iter, // maximum number of iterations
- 2); // verbosity of the solver
- InverseOperatorResult statistics;
- std::cout << "solve linear system for xa.\n";
- solver.apply(x_a, load_a_, statistics);
- std::cout << "solve linear system for xK1.\n";
- solver.apply(x_K1, load_K1_, statistics);
- std::cout << "solve linear system for xK2.\n";
- solver.apply(x_K2, load_K2_, statistics);
- std::cout << "solve linear system for xK3.\n";
- solver.apply(x_K3, load_K3_, statistics);
-
- }
-
- if (parameterSet.get<int>("useQR") == 1){
-
- std::cout << "We use QR solver.\n";
- log << "solveLinearSystems: We use QR solver.\n";
- //TODO install suitesparse
- SPQR<MatrixCT> sPQR(stiffnessMatrix_);
- sPQR.setVerbosity(1);
- InverseOperatorResult statisticsQR;
-
- sPQR.apply(x_a, load_a_, statisticsQR);
- sPQR.apply(x_K1, load_K1_, statisticsQR);
- sPQR.apply(x_K2, load_K2_, statisticsQR);
- sPQR.apply(x_K3, load_K3_, statisticsQR);
-
- }
-
- // Write solutions in logs
- if (parameterSet.get<int>("write_solutions_corrector_problems") == 1){
- log << "\nSolution of Corrector problems:\n";
-
- auto sizeComp = x_a.size()/nCompo;
- std::vector<VectorRT> x_a_vec(sizeComp);
- std::vector<VectorRT> x_K1_vec(sizeComp);
- std::vector<VectorRT> x_K2_vec(sizeComp);
- std::vector<VectorRT> x_K3_vec(sizeComp);
-
- for (int i=0; i < x_a_vec.size(); i++){
- x_a_vec[i] = VectorRT{x_a[i], x_a[i + sizeComp], x_a[i + 2*sizeComp]};
- x_K1_vec[i] = VectorRT{x_K1[i], x_K1[i + sizeComp], x_K1[i + 2*sizeComp]};
- x_K2_vec[i] = VectorRT{x_K2[i], x_K2[i + sizeComp], x_K2[i + 2*sizeComp]};
- x_K3_vec[i] = VectorRT{x_K3[i], x_K3[i + sizeComp], x_K3[i + 2*sizeComp]};
- }
-
- log << "\nxa:\n";
- for (int i=0; i < x_a_vec.size(); i++)
- log << i << ": " << x_a_vec[i] << std::endl;
-
- log << "\nxK1:\n";
- for (int i=0; i < x_K1_vec.size(); i++)
- log << i << ": " << x_K1_vec[i] << std::endl;
-
- log << "\nxK2:\n";
- for (int i=0; i < x_K2_vec.size(); i++)
- log << i << ": " << x_K2_vec[i] << std::endl;
-
- log << "\nxK3:\n";
- for (int i=0; i < x_K3_vec.size(); i++)
- log << i << ": " << x_K3_vec[i] << std::endl;
-
- /*
- log << "\nxa:\n";
- for (int i=0; i < x_a.size(); i++)
- log << i << " " << x_a[i] << std::endl;
-
- log << "\nxK1:\n";
- for (int i=0; i < x_K1.size(); i++)
- log << i << " " << x_K1[i] << std::endl;
-
- log << "\nxK2:\n";
- for (int i=0; i < x_K2.size(); i++)
- log << i << " " << x_K2[i] << std::endl;
-
- log << "\nxK3:\n";
- for (int i=0; i < x_K3.size(); i++)
- log << i << " " << x_K3[i] << std::endl;
- */
- }
-
- }
-
-
- void writeSystemMatrixAndRHSs()
- {
- ParameterTree parameterSet = cellAssembler_.getParameterSet();
- auto logP = cellAssembler_.getLog();
- auto & log = *logP;
-
- if (parameterSet.get<int>("write_stiffness_matrix_Cell") == 1){
- std::cout << "write stiffness matrix of cell problem in logfile\n";
-
- log << "\nwriteSystemMatrixAndRHSs:\n";
- printSparseMatrix(log, stiffnessMatrix_, "stiffnessMatrix with periodic bc", "");
- //log << "\n";
- //printmatrix(log, stiffnessMatrix_, "stiffnessMatrix with periodic bc", "");
- //writeMatrixToMatlab(stiffnessMatrix, "../outputs/reducedrod_output/postprocessing/stiffness_matrix_for_octave.txt");
- }
-
- if (parameterSet.get<int>("write_loads_Cell") == 1){
-
- std::cout << "write rhs of cell problem in logfile\n";
-
- auto sizeComp = load_a_.size()/nCompo;
- std::vector<VectorRT> load_a_vec(sizeComp);
- std::vector<VectorRT> load_K1_vec(sizeComp);
- std::vector<VectorRT> load_K2_vec(sizeComp);
- std::vector<VectorRT> load_K3_vec(sizeComp);
-
- for (int i=0; i < load_a_vec.size(); i++){
- load_a_vec[i] = VectorRT{load_a_[i], load_a_[i + sizeComp], load_a_[i + 2*sizeComp]};
- load_K1_vec[i] = VectorRT{load_K1_[i], load_K1_[i + sizeComp], load_K1_[i + 2*sizeComp]};
- load_K2_vec[i] = VectorRT{load_K2_[i], load_K2_[i + sizeComp], load_K2_[i + 2*sizeComp]};
- load_K3_vec[i] = VectorRT{load_K3_[i], load_K3_[i + sizeComp], load_K3_[i + 2*sizeComp]};
- }
-
- log << "\n RHS stretch strain - load_a:\n";
- for (int i=0; i < load_a_vec.size(); i++)
- log << i << " " << load_a_vec[i] << std::endl;
-
- log << "\n RHS twist strain - load_K1:\n";
- for (int i=0; i < load_K1_vec.size(); i++)
- log << i << " " << load_K1_vec[i] << std::endl;
-
- log << "\n RHS bend1 strain - load_K2:\n";
- for (int i=0; i < load_K2_vec.size(); i++)
- log << i << " " << load_K2_vec[i] << std::endl;
-
- log << "\n RHS bend2 strain - load_K3:\n";
- for (int i=0; i < load_K3_vec.size(); i++)
- log << i << " " << load_K3_vec[i] << std::endl;
-
- }
- }
-
-
- void plotLoads()
- {
- auto basis = cellAssembler_.getBasis();
- auto perIndexMap = cellAssembler_.getIndexMapPeriodicBoundary();
- auto size = basis->size();
-
- VectorCT load_a_Ce(size);
- VectorCT load_K1_Ce(size);
- VectorCT load_K2_Ce(size);
- VectorCT load_K3_Ce(size);
- //index_PeriodicFE = indexmap(index_FE)
- for (int i = 0; i < size; i++){
- auto idx = perIndexMap[i];
- load_a_Ce[i] = load_a_[idx];
- load_K1_Ce[i] = load_K1_[idx];
- load_K2_Ce[i] = load_K2_[idx];
- load_K3_Ce[i] = load_K3_[idx];
- }
-
- std::cout << "plot loads.\n";
-
- auto stretchF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(Functions::subspaceBasis(*basis), HierarchicVectorView(load_a_Ce));
- auto twistF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(Functions::subspaceBasis(*basis), HierarchicVectorView(load_K1_Ce));
- auto bend1F = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(Functions::subspaceBasis(*basis), HierarchicVectorView(load_K2_Ce));
- auto bend2F = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(Functions::subspaceBasis(*basis), HierarchicVectorView(load_K3_Ce));
-
- // Write result to VTK file, Maybe we write for every gamma?
- VTKWriter<typename Basis::GridView> vtkWriter(basis->gridView());
- vtkWriter.addVertexData(stretchF, VTK::FieldInfo("load_a_stretch", VTK::FieldInfo::Type::vector, nCompo));
- vtkWriter.addVertexData(twistF, VTK::FieldInfo("load_K1_twist", VTK::FieldInfo::Type::vector, nCompo));
- vtkWriter.addVertexData(bend1F, VTK::FieldInfo("load_K2_bend1", VTK::FieldInfo::Type::vector, nCompo));
- vtkWriter.addVertexData(bend2F, VTK::FieldInfo("load_K3_bend2", VTK::FieldInfo::Type::vector, nCompo));
-
- std::string outputPath = cellAssembler_.getParameterSet().get("outputPath", "../outputs/clampedprestressedrod_output");
- vtkWriter.write(outputPath + "/loads_cell_problems");
- }
-
-}; // end class
-
-
-
-#endif
diff --git a/dune/microstructure/effectivemodelcomputer.hh b/dune/microstructure/effectivemodelcomputer.hh
deleted file mode 100644
index 0020a959..00000000
--- a/dune/microstructure/effectivemodelcomputer.hh
+++ /dev/null
@@ -1,228 +0,0 @@
-#ifndef DUNE_REDUCEDROD_EFFECTIVEMODELCOMPUTER_HH
-#define DUNE_REDUCEDROD_EFFECTIVEMODELCOMPUTER_HH
-
-
-#include <dune/reducedrod/cellsolver.hh>
-
-
-template <class Basis>
-class EffectiveModelComputer {
-
-public:
-
- static const int nCompo = 3;
- static const int nCells = 4;
-
- using Domain = typename CellAssembler<Basis>::Domain;
-
- using VectorRT = typename CellSolver<Basis>::VectorRT;
- using MatrixRT = typename CellAssembler<Basis>::MatrixRT;
-
- using Func2Tensor = typename CellAssembler<Basis>::Func2Tensor;
-
- using VectorCT = typename CellSolver<Basis>::VectorCT;
-
- using HierarchicVectorView = typename CellSolver<Basis>::HierarchicVectorView;
-
-protected:
-
- CellSolver<Basis>& cellSolver_;
- Func2Tensor prestrain_;
-
- VectorCT x_a_, x_K1_, x_K2_, x_K3_; // in torus basis
- FieldMatrix<double, nCells, nCells> M_Ce_;
- FieldVector<double, nCells> aeffKeff_;
-
-public:
-
- // Constructor
- EffectiveModelComputer(CellSolver<Basis>& cellSolver, Func2Tensor prestrain)
- : cellSolver_(cellSolver), prestrain_(prestrain)
- {
- cellSolver.solveCorrectorSystems(x_a_, x_K1_, x_K2_, x_K3_ );
- computeEffectiveModel();
- //computeCorrector();
- }
-
- VectorCT * getCorr_a(){return &x_a_;}
- VectorCT * getCorr_K1(){return &x_K1_;}
- VectorCT * getCorr_K2(){return &x_K2_;}
- VectorCT * getCorr_K3(){return &x_K3_;}
-
- FieldMatrix<double, nCells, nCells> getEffMaterial(){return M_Ce_;}
- FieldVector<double, nCells> getEffPreStrain(){return aeffKeff_;}
-
- CellSolver<Basis> getCellSolver(){return cellSolver_;}
- const Basis* getCellBasis() {
- auto cellAssembler = cellSolver_.getCellAssembler();
- return cellAssembler.getBasis();}
-
- void plotCorrectors()
- {
- std::cout << "Write corrector solutions.\n";
-
- auto cellAssembler = cellSolver_.getCellAssembler();
- auto basis = getCellBasis();
- auto perIndexMap = cellAssembler.getIndexMapPeriodicBoundary();
- auto size = basis->size();
-
- VectorCT x_a_Ce(size);
- VectorCT x_K1_Ce(size);
- VectorCT x_K2_Ce(size);
- VectorCT x_K3_Ce(size);
- //index_PeriodicFE = indexmap(index_FE)
- for (int i = 0; i < size; i++){
- auto idx = perIndexMap[i];
- x_a_Ce[i] = x_a_[idx];
- x_K1_Ce[i] = x_K1_[idx];
- x_K2_Ce[i] = x_K2_[idx];
- x_K3_Ce[i] = x_K3_[idx];
- }
-
- VectorCT x_corrector(size);
- computeCorrector(x_corrector);
-
- auto stretchF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(Functions::subspaceBasis(*basis), HierarchicVectorView(x_a_Ce));
- auto twistF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(Functions::subspaceBasis(*basis), HierarchicVectorView(x_K1_Ce));
- auto bend1F = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(Functions::subspaceBasis(*basis), HierarchicVectorView(x_K2_Ce));
- auto bend2F = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(Functions::subspaceBasis(*basis), HierarchicVectorView(x_K3_Ce));
- auto correctorF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(Functions::subspaceBasis(*basis), HierarchicVectorView(x_corrector));
-
- // Write result to VTK file, Maybe we write for every gamma?
- VTKWriter<typename Basis::GridView> vtkWriter(basis->gridView());
- vtkWriter.addVertexData(stretchF, VTK::FieldInfo("phi_a_stretch", VTK::FieldInfo::Type::vector, nCompo));
- vtkWriter.addVertexData(twistF, VTK::FieldInfo("phi_K1_twist", VTK::FieldInfo::Type::vector, nCompo));
- vtkWriter.addVertexData(bend1F, VTK::FieldInfo("phi_K2_bend1", VTK::FieldInfo::Type::vector, nCompo));
- vtkWriter.addVertexData(bend2F, VTK::FieldInfo("phi_K3_bend2", VTK::FieldInfo::Type::vector, nCompo));
- vtkWriter.addVertexData(correctorF, VTK::FieldInfo("phi_full", VTK::FieldInfo::Type::vector, nCompo));
-
- std::string outputPath = cellAssembler.getParameterSet().get("outputPath", "../outputs/clampedprestressedrod_output");
- vtkWriter.write(outputPath + "/corrector_solutions");
-
- //add corrector
- }
-
-
- void computeCorrector(VectorCT& x_corrector)
- {
- VectorCT x_Corrector_Torus(x_a_);
- x_Corrector_Torus *= aeffKeff_[0];
- auto tmp_xK1 = x_K1_;
- auto tmp_xK2 = x_K2_;
- auto tmp_xK3 = x_K3_;
- tmp_xK1 *= aeffKeff_[1];
- tmp_xK2 *= aeffKeff_[2];
- tmp_xK3 *= aeffKeff_[3];
- x_Corrector_Torus += tmp_xK1;
- x_Corrector_Torus += tmp_xK2;
- x_Corrector_Torus += tmp_xK3;
-
- auto cellAssembler = cellSolver_.getCellAssembler();
- auto basis = cellAssembler.getBasis();
- auto perIndexMap = cellAssembler.getIndexMapPeriodicBoundary();
-
- x_corrector.resize(basis->size());
- for (int i = 0; i < basis->size(); i++)
- x_corrector[i] = x_Corrector_Torus[perIndexMap[i]];
- }
-
-
-private:
- void computeEffectiveModel()
- {
- auto cellAssembler = cellSolver_.getCellAssembler();
- ParameterTree parameterSet = cellAssembler.getParameterSet();
- auto logP = cellAssembler.getLog();
- auto & log = *logP;
-
- // lateral stretch strain E_U1 = e1 ot e1
- Func2Tensor E_a = [] (const Domain& z) //extra Klasse
- {
- return MatrixRT{{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
- };
-
- // twist in x2-x3 plane E_K1 = sym (e1 x (0,x2,x3)^T ot e1)
- Func2Tensor E_K1 = [] (const Domain& z)
- {
- return MatrixRT{{0.0,-0.5*z[2], 0.5*z[1]}, {-0.5*z[2], 0.0 , 0.0 }, {0.5*z[1],0.0,0.0}};
- };
-
- // bend strain in x1-x2 plane E_K2 = sym (e2 x (0,x2,x3)^T ot e1)
- Func2Tensor E_K2 = [] (const Domain& z)
- {
- return MatrixRT{{z[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0}};
- };
-
- // bend strain in x1-x3 plane E_K3 = sym (e3 x (0,x2,x3)^T ot e1)
- Func2Tensor E_K3 = [] (const Domain& z)
- {
- return MatrixRT{{-z[1], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
- };
-
- VectorCT *xCorrectors[nCells] = {&x_a_, &x_K1_, &x_K2_, &x_K3_};
- VectorCT *loadsCell[nCells] = {cellSolver_.getLoad_a(), cellSolver_.getLoad_K1(), cellSolver_.getLoad_K2(), cellSolver_.getLoad_K3()};
-
- //using GridViewFunction_Ce = decltype(E_a);
- const std::array<Func2Tensor*, nCells> strainsE = {&E_a, &E_K1, &E_K2, &E_K3};
-
- //wir brauchen sym grad phi als gridFunction!
-
- FieldMatrix<double, nCells, nCells> QEE(0);
-
- // Assemble 4x4 matrix //symmetrisch? reicht hälfte zu füllen?
- for (unsigned int k = 0; k < nCells; k++)
- for (unsigned int l = 0; l < nCells; l++){
- auto energyEkEl = cellAssembler.energy(*strainsE[k], *strainsE[l]); //<LE,E>
- QEE[k][l] = energyEkEl;
- M_Ce_[k][l]= 3.0*((*xCorrectors[k])*(*loadsCell[l])) + energyEkEl; //<LChi, Chi>= <E,Chi>, 2<LChi, E> // 3* ?
- }// oder braucht man nicht periodische objekte?
-
- /*
- FieldMatrix<double, nCells, nCells> QChiChi(0);
- FieldMatrix<double, nCells, nCells> QEChi(0);
- if (parameterSet.get<bool>("material_computation_2"))
- for (unsigned int k = 0; k < nCells; k++)
- for (unsigned int l = 0; l < nCells; l++){
- auto energyEE = cellAssembler.computeScalarTerm(*strainsE[k], *strainsE[l]); //<LE,E>
- auto energyChiChi = cellAssembler.computeScalarTerm(*xCorrectors[k], *xCorrectors[l]);
- auto energyEChi = cellAssembler.computeScalarTerm(*strainsE[k], *xCorrectors[l]);
- QEE[k][l] = energyEE;
- QChiChi[k][l] = energyChiChi;
- QEChi[k][l] = energyEChi;
- M_Ce_[k][l]= energyEE + energyChiChi + 2*energyEChi; //<LChi, Chi>= <E,Chi>, 2<LChi, E> // eher so ??
- }
- */
-
-
-
- VectorCT load_B;
- cellAssembler.assembleLoadVector(prestrain_, load_B);
-
- // Assemble rhs by projecting prestrain to strain basis
- FieldVector<double, nCells> b;
- for (unsigned int k = 0; k < nCells; k++){
- auto strainsEkB = cellAssembler.energy(*strainsE[k], prestrain_); // <LB, E>
- b[k]= (*xCorrectors[k]) * load_B + strainsEkB; // < B, chi>
- }
-
- M_Ce_.solve( aeffKeff_ ,b);
-
- log << "\nQEE= \n";
- for (int i = 0; i < nCells; i++)
- log << QEE[i] << std::endl;
- log << "\nM_Ce= \n";
- for (int i = 0; i < nCells; i++)
- log << M_Ce_[i] << std::endl;
- log << "\nb= " << b << std::endl;
- log << "\nstretch twist bend1 bend2\naeffKeff = "<< aeffKeff_ << std::endl;
- std::cout << "\naeffKeff: stretch twist bend1 bend2 = "<< aeffKeff_ << std::endl;
-
- }
-
-}; // end class
-
-
-
-
-
-#endif
diff --git a/dune/microstructure/microstructuredrodgrid.hh b/dune/microstructure/microstructuredrodgrid.hh
deleted file mode 100644
index 252af1f0..00000000
--- a/dune/microstructure/microstructuredrodgrid.hh
+++ /dev/null
@@ -1,115 +0,0 @@
-
-#ifndef DUNE_REDUCEDROD_MICROSTRUCTUREDRODGRID_HH
-#define DUNE_REDUCEDROD_MICROSTRUCTUREDRODGRID_HH
-
-#include <dune/grid/yaspgrid.hh>
-#include <dune/grid/onedgrid.hh>
-#include <dune/functions/functionspacebases/powerbasis.hh>
-#include <dune/functions/functionspacebases/lagrangebasis.hh>
-
-using namespace Dune;
-using namespace Functions::BasisBuilder;
-
-
-class MicrostructuredRodGrid {
-
-public:
-
- //static const int orderFE_Ce = 1;
- static const int dim = 3;
- static const int nCompo = 3;
-
- using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
- using Rod1DGridType = OneDGrid;
- using Rod3DGridType = CellGridType;
-
- using IsometricRodBasisType = Functions::LagrangeBasis<Rod1DGridType::LeafGridView, 1>;
-
-protected:
-
- const double width_;
- const double len_;
- const double eps_;
- const double h_; //thickness h könnte auch in plotmethoden vorkommen und Querschnitt bleibt skaliert
-
- const std::array<int,dim> nE_Ce_;
- const FieldVector<double,dim> bboxL_Ce_;
- const FieldVector<double,dim> bboxR_Ce_;
-
- const int nE_R1D_;
-
- const std::array<int,dim> nE_R3D_;
- const FieldVector<double,dim> bboxL_R3D_;
- const FieldVector<double,dim> bboxR_R3D_;
-
- const CellGridType grid_Ce_;
- const Rod1DGridType grid_R1D_;
- const Rod3DGridType grid_R3D_;
-
- const IsometricRodBasisType basis_IsoRod_;
-
-
-
-public:
-
- // Constructor
- MicrostructuredRodGrid( const double width,
- const double len,
- const double eps,
- const double h,
- const std::array<int,dim> nE_Ce):
- width_(width),
- len_(len),
- eps_(eps),
- h_(h),
-
- nE_Ce_(nE_Ce),
- bboxL_Ce_({0.0, -width/2.0, -width/2.0}),
- bboxR_Ce_({1.0, width/2.0, width/2.0}),
-
- nE_R1D_(int( (len *nE_Ce[0])/ eps) ),
-
-
- //nE_R3D_({nE_R1D_, int(double(nE_Ce[1])/h), int(double(nE_Ce[2])/h)}),
- nE_R3D_({nE_R1D_, nE_Ce[1], nE_Ce[2]}),
- bboxL_R3D_({0, h*bboxL_Ce_[1], h*bboxL_Ce_[2] }),
- bboxR_R3D_({len, h*bboxR_Ce_[1], h*bboxR_Ce_[2] }),
-
- grid_Ce_(bboxL_Ce_, bboxR_Ce_, nE_Ce),
- grid_R1D_(nE_R1D_, 0, len),
- grid_R3D_(bboxL_R3D_, bboxR_R3D_, nE_R3D_),
-
- basis_IsoRod_(grid_R1D_.leafGridView())
-
- {}
-
-
- // Getter
- double getWidth() {return width_;}
- double getLen() {return len_;}
- double getEps() {return eps_;}
- double getH() {return h_;}
-
- std::array<int,dim> getNECell() {return nE_Ce_;}
- FieldVector<double,dim> getBboxLCell() {return bboxL_Ce_;}
- FieldVector<double,dim> getBboxRCell() {return bboxR_Ce_;}
-
- const int getNERod1D(){return nE_R1D_;}
-
- std::array<int,dim> getNERod3D() {return nE_R3D_;}
- FieldVector<double,dim> getBboxLRod3D() {return bboxL_R3D_;}
- FieldVector<double,dim> getBboxRRod3D() {return bboxR_R3D_;}
-
- const CellGridType & getCellGrid() {return grid_Ce_;}
- const Rod1DGridType & getRod1DGrid() {return grid_R1D_;}
- const Rod3DGridType & getRod3DGrid() {return grid_R3D_;}
-
- const IsometricRodBasisType* getIsometricRodBasis() {return &basis_IsoRod_;}
-
- //return gridViews
-
- //Rod1DGridType copyRod1DGrid() {Rod1DGridType grid = grid_R1D_; return grid;}
-
-}; // end class
-
-#endif
\ No newline at end of file
diff --git a/dune/microstructure/testclass.hh b/dune/microstructure/testclass.hh
deleted file mode 100644
index 7f843176..00000000
--- a/dune/microstructure/testclass.hh
+++ /dev/null
@@ -1,68 +0,0 @@
-
-template<class Basis>
-class TestClass {
-
-public:
- static const int nCompo = 3;
-
- using GridView = typename Basis::GridView;
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRange = FieldMatrix< double, nCompo, nCompo>;
- using Func2Tensor = std::function< MatrixRange(const Domain&) >;
-
- const Basis& basis_;
- Func2Tensor prestress_;
-
- MatrixRange evaluateZero()
- {
- auto preStressGrid = Functions::makeGridViewFunction(prestress_, basis_.gridView());
- return prestress_({0,0,0});
- }
-
-
-
- TestClass(const Basis& basis, Func2Tensor prestress) : basis_(basis), prestress_(prestress) {}
-};
-
-
-
-
-
-
-
-
-
-
-
-
-/*
-template <class Basis, class GridF1, class GridF2> //, class GridScalar> //, class LocalScalar, class Local2Tensor> // LocalFunction derived from basis?
-class TestClass {
-
- using GridView = typename Basis::GridView;
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using ScalarRange = FieldVector< double, 1>;
- using FuncScalar = std::function< ScalarRange(const Domain&) >;
-
-private:
-
- //const int dim = 3;
- const Basis basis_;
- const ParameterTree& parameterSet_;
- const GridF1& mu1_;
- const GridF2& mu2_;
-
- std::map<int,int> perIndexMap_;
-
-public:
-
- TestClass(const Basis& basis,
- const GridF1& mu
- const ParameterTree& parameterSet)
- : basis_(basis),
- mu1_(mu),
- mu2_(mu)
- {}
-
-};
-*/
--
GitLab