Skip to content
Snippets Groups Projects
Commit 4badc347 authored by Klaus Böhnlein's avatar Klaus Böhnlein
Browse files

deleted unnecessary Files

parent 91c99a53
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/microstructure/cellassembler.hh deleted 100644 → 0
+ 0
− 746
View file @ 91c99a53
#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; };
*/
dune/microstructure/cellsolver.hh deleted 100644 → 0
+ 0
− 272
View file @ 91c99a53
#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
dune/microstructure/effectivemodelcomputer.hh deleted 100644 → 0
+ 0
− 228
View file @ 91c99a53
#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
dune/microstructure/microstructuredrodgrid.hh deleted 100644 → 0
+ 0
− 115
View file @ 91c99a53
#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
dune/microstructure/testclass.hh deleted 100644 → 0
+ 0
− 68
View file @ 91c99a53
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)
{}
};
*/
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment