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

found error in local assembly

parent 4fb04581
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
.kdev4/dune-microstructure.kdev4
+ 3
− 0
View file @ bb147190
...@@ -35,3 +35,6 @@ Project Target=dune-microstructure,src,dune-microstructure ...@@ -35,3 +35,6 @@ Project Target=dune-microstructure,src,dune-microstructure
Use External Terminal=false Use External Terminal=false
Working Directory=file:///home/klaus/Desktop/DUNE2.8/dune-microstructure/build-cmake/src Working Directory=file:///home/klaus/Desktop/DUNE2.8/dune-microstructure/build-cmake/src
isExecutable=false isExecutable=false
[Project]
VersionControlSupport=kdevgit
build-cmake/CMakeFiles/CMakeError.log
+ 22
− 0
View file @ bb147190
...@@ -519,3 +519,25 @@ Determing location of SIONLIB failed: ...@@ -519,3 +519,25 @@ Determing location of SIONLIB failed:
Include directory: Include directory:
Library directory: Library directory:
Determing location of ARPACK failed:
Libraries to link against:
Determing location of ARPACK++ failed:
Include directory:
Libraries to link against:
Determing location of SIONLIB failed:
Include directory:
Library directory:
Determing location of ARPACK failed:
Libraries to link against:
Determing location of ARPACK++ failed:
Include directory:
Libraries to link against:
Determing location of SIONLIB failed:
Include directory:
Library directory:
build-cmake/CMakeFiles/CMakeOutput.log
+ 8
− 0
View file @ bb147190
...@@ -862,3 +862,11 @@ Determining location of SuperLU succeeded: ...@@ -862,3 +862,11 @@ Determining location of SuperLU succeeded:
Include directory: /usr/include/superlu Include directory: /usr/include/superlu
Library directory: /usr/lib/x86_64-linux-gnu/libsuperlu.so Library directory: /usr/lib/x86_64-linux-gnu/libsuperlu.so
Determining location of SuperLU succeeded:
Include directory: /usr/include/superlu
Library directory: /usr/lib/x86_64-linux-gnu/libsuperlu.so
Determining location of SuperLU succeeded:
Include directory: /usr/include/superlu
Library directory: /usr/lib/x86_64-linux-gnu/libsuperlu.so
build-cmake/DartConfiguration.tcl
+ 2
− 2
View file @ bb147190
...@@ -57,9 +57,9 @@ P4UpdateOptions: ...@@ -57,9 +57,9 @@ P4UpdateOptions:
P4UpdateCustom: P4UpdateCustom:
# Generic update command # Generic update command
UpdateCommand: UpdateCommand: /usr/bin/git
UpdateOptions: UpdateOptions:
UpdateType: UpdateType: git
# Compiler info # Compiler info
Compiler: /usr/bin/c++ Compiler: /usr/bin/c++
......
src/.dune-microstructure.cc.kate-swp 0 → 100644
+ 0
− 0
View file @ bb147190
File added
src/dune-microstructure.cc
+ 410
− 78
View file @ bb147190
...@@ -54,40 +54,265 @@ ...@@ -54,40 +54,265 @@
using namespace Dune; using namespace Dune;
template<class LocalView, class Matrix, class localFunction1, class localFunction2> template<class LocalView, class Matrix, class localFunction1, class localFunction2>
void computeElementStiffnessMatrix(const LocalView& localView, void computeElementStiffnessMatrixOLD(const LocalView& localView,
// Matrix<FieldMatrix<double,1,1> >& elementMatrix,
Matrix& elementMatrix, Matrix& elementMatrix,
const localFunction1& mu, const localFunction1& mu,
const localFunction2& lambda, const localFunction2& lambda,
const double gamma, const double gamma
const double phiOffset
) )
{ {
// Local StiffnessMatrix of the form: // Local StiffnessMatrix of the form:
// | phi*phi m*phi | // | phi*phi m*phi |
// | phi *m m*m | // | phi *m m*m |
using Element = typename LocalView::Element; using Element = typename LocalView::Element;
const Element element = localView.element(); const Element element = localView.element();
auto geometry = element.geometry();
constexpr int dim = Element::dimension; constexpr int dim = Element::dimension;
constexpr int nCompo = dim; constexpr int nCompo = dim;
using MatrixRT = FieldMatrix< double, nCompo, nCompo>; using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
// using Domain = typename LocalView::GridView::Codim<0>::Geometry::GlobalCoordinate; // using Domain = typename LocalView::GridView::Codim<0>::Geometry::GlobalCoordinate;
// using FuncScalar = std::function< ScalarRT(const Domain&) >; // using FuncScalar = std::function< ScalarRT(const Domain&) >;
// using Func2Tensor = std::function< MatrixRT(const Domain&) >; // using Func2Tensor = std::function< MatrixRT(const Domain&) >;
elementMatrix.setSize(localView.size()+3, localView.size()+3);
elementMatrix = 0;
// LocalBasis-Offset
const int localPhiOffset = localView.size();
const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
const auto nSf = localFiniteElement.localBasis().size();
///////////////////////////////////////////////
// Basis for R_sym^{2x2} // wird nicht als Funktion benötigt da konstant...
//////////////////////////////////////////////
MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
MatrixRT G_3 {{0.0, 0.5, 0.0}, {0.5, 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
//print:
printmatrix(std::cout, basisContainer[0] , "G_1", "--");
printmatrix(std::cout, basisContainer[1] , "G_2", "--");
printmatrix(std::cout, basisContainer[2] , "G_3", "--");
////////////////////////////////////////////////////
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
for (const auto& quadPoint : quad)
{
const auto& quadPos = quadPoint.position();
const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
const auto integrationElement = geometry.integrationElement(quadPoint.position());
// std::vector<FieldMatrix<double,1,dim> > referenceGradients; // old
std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
localFiniteElement.localBasis().evaluateJacobian(
quadPoint.position(),
referenceGradients);
// Compute the shape function gradients on the grid element
std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
// std::vector< VectorRT> gradients(referenceGradients.size());
for (size_t i=0; i<gradients.size(); i++)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
for (size_t k=0; k < nCompo; k++)
for (size_t l=0; l< nCompo; l++)
{
for (size_t i=0; i < nSf; i++)
for (size_t j=0; j < nSf; j++ )
{
// (scaled) Deformation gradient of the ansatz basis function
MatrixRT defGradientU(0);
defGradientU[k][0] = gradients[i][0]; // Y
defGradientU[k][1] = gradients[i][1]; //X2
defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
// (scaled) Deformation gradient of the test basis function
MatrixRT defGradientV(0);
defGradientV[l][0] = gradients[j][0]; // Y
defGradientV[l][1] = gradients[j][1]; //X2
defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
// symmetric Gradient (Elastic Strains)
MatrixRT strainU, strainV;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
// strainV[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // same ? genügt strainU
}
// St.Venant-Kirchhoff stress
// < L sym[D_gamma*nabla phi_i], sym[D_gamma*nabla phi_j] >
// stressU*strainV
MatrixRT stressU(0);
stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
double trace = 0;
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] * strainU[ii][jj]; // "phi*phi"-part
energyDensity += stressU[ii][jj] * strainV[ii][jj];
/* size_t row = localView.tree().child(k).localIndex(i); // kann auf Unterscheidung zwischen k & l verzichten?!
size_t col = localView.tree().child(l).localIndex(j); */ // siehe DUNE-book p.394
// size_t col = localView.tree().child(k).localIndex(j); // Indizes mit k=l genügen .. Kroenecker-Delta_kl NEIN???
size_t col = localView.tree().child(k).localIndex(i); // kann auf Unterscheidung zwischen k & l verzichten?!
size_t row = localView.tree().child(l).localIndex(j);
elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
for( size_t m=0; m<3; m++)
{
// < L G_i, sym[D_gamma*nabla phi_j] >
// L G_i* strainV
// St.Venant-Kirchhoff stress
MatrixRT stressG(0);
stressG.axpy(2*mu(quadPos), basisContainer[m]); //this += 2mu *strainU
double traceG = 0;
for (int ii=0; ii<nCompo; ii++)
traceG += basisContainer[m][ii][ii];
for (int ii=0; ii<nCompo; ii++)
stressG[ii][ii] += lambda(quadPos) * traceG;
double energyDensityGphi = 0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
energyDensityGphi += stressG[ii][jj] * strainV[ii][jj]; // "m*phi"-part
auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
elementMatrix[row][localPhiOffset+m] += value;
elementMatrix[localPhiOffset+m][row] += value; // ---- reicht das ? --- TODO
// // equivalent? :
// // < L sym[D_gamma*nabla phi_i], G_j >
// double energyDensityPhiG = 0;
// for (int ii=0; ii<nCompo; ii++)
// for (int jj=0; jj<nCompo; jj++)
// energyDensityPhiG += stressU[ii][jj] * basisContainer[m][ii][jj]; // "phi*m"-part
//
// elementMatrix[localPhiOffset+m][row] += energyDensityPhiG * quadPoint.weight() * integrationElement;
//
// St.Venant-Kirchhoff stress
// < L G_alpha, G_alpha >
for(size_t n=0; n<3; n++)
{
double energyDensityGG = 0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
energyDensityGG += stressG[ii][jj] * basisContainer[n][ii][jj]; // "m*m"-part
elementMatrix[localPhiOffset+m][localPhiOffset+n]= energyDensityGG * quadPoint.weight() * integrationElement;
}
}
}
}
}
}
template<class LocalView, class Matrix, class localFunction1, class localFunction2>
void computeElementStiffnessMatrix(const LocalView& localView,
Matrix& elementMatrix,
const localFunction1& mu,
const localFunction2& lambda,
const double gamma
)
{
// Local StiffnessMatrix of the form:
// | phi*phi m*phi |
// | phi *m m*m |
using Element = typename LocalView::Element;
const Element element = localView.element();
auto geometry = element.geometry(); auto geometry = element.geometry();
constexpr int dim = Element::dimension;
constexpr int nCompo = dim;
using MatrixRT = FieldMatrix< double, nCompo, nCompo>; using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
// using Domain = typename LocalView::GridView::Codim<0>::Geometry::GlobalCoordinate;
// using FuncScalar = std::function< ScalarRT(const Domain&) >;
// using Func2Tensor = std::function< MatrixRT(const Domain&) >;
elementMatrix.setSize(localView.size()+3, localView.size()+3); elementMatrix.setSize(localView.size()+3, localView.size()+3);
elementMatrix = 0; elementMatrix = 0;
const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
// LocalBasis-Offset
const int localPhiOffset = localView.size();
const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
const auto nSf = localFiniteElement.localBasis().size(); const auto nSf = localFiniteElement.localBasis().size();
...@@ -95,8 +320,7 @@ void computeElementStiffnessMatrix(const LocalView& localView, ...@@ -95,8 +320,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
/////////////////////////////////////////////// ///////////////////////////////////////////////
// Basis for R_sym^{2x2} // wird nicht als Funktion benötigt da konstant... // Basis for R_sym^{2x2} // wird nicht als Funktion benötigt da konstant...
// //////////////////////////////////////////////
MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}}; MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}}; MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
MatrixRT G_3 {{0.0, 0.5, 0.0}, {0.5, 0.0, 0.0}, {0.0, 0.0, 0.0}}; MatrixRT G_3 {{0.0, 0.5, 0.0}, {0.5, 0.0, 0.0}, {0.0, 0.0, 0.0}};
...@@ -107,7 +331,7 @@ void computeElementStiffnessMatrix(const LocalView& localView, ...@@ -107,7 +331,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
printmatrix(std::cout, basisContainer[0] , "G_1", "--"); printmatrix(std::cout, basisContainer[0] , "G_1", "--");
printmatrix(std::cout, basisContainer[1] , "G_2", "--"); printmatrix(std::cout, basisContainer[1] , "G_2", "--");
printmatrix(std::cout, basisContainer[2] , "G_3", "--"); printmatrix(std::cout, basisContainer[2] , "G_3", "--");
////////////////////////////////////////////////////
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1); int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR); const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
...@@ -135,112 +359,215 @@ void computeElementStiffnessMatrix(const LocalView& localView, ...@@ -135,112 +359,215 @@ void computeElementStiffnessMatrix(const LocalView& localView,
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]); jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
// "phi*phi"-part
for (size_t k=0; k < nCompo; k++)
for (size_t l=0; l< nCompo; l++)
{
for (size_t i=0; i < nSf; i++)
for (size_t j=0; j < nSf; j++ )
{
// (scaled) Deformation gradient of the ansatz basis function
MatrixRT defGradientU(0);
defGradientU[k][0] = gradients[i][0]; // Y
defGradientU[k][1] = gradients[i][1]; //X2
defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
// (scaled) Deformation gradient of the test basis function
MatrixRT defGradientV(0);
defGradientV[l][0] = gradients[j][0]; // Y
defGradientV[l][1] = gradients[j][1]; //X2
defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
// symmetric Gradient (Elastic Strains)
MatrixRT strainU, strainV;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
// strainV[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // same ? genügt strainU
}
// St.Venant-Kirchhoff stress
// < L sym[D_gamma*nabla phi_i], sym[D_gamma*nabla phi_j] >
// stressU*strainV
MatrixRT stressU(0);
stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
double trace = 0;
for (int ii=0; ii<nCompo; ii++)
trace += strainU[ii][ii];
for (int ii=0; ii<nCompo; ii++)
stressU[ii][ii] += lambda(quadPos) * trace;
for (size_t i=0; i < nSf; i++) // Local energy density: stress times strain
for (size_t j=0; j < nSf; j++ ) double energyDensity = 0;
for (size_t k=0; k < nCompo; k++) for (int ii=0; ii<nCompo; ii++)
// for (size_t l=0; l< nCompo; l++) for (int jj=0; jj<nCompo; jj++)
{ // energyDensity += stressU[ii][jj] * strainU[ii][jj]; // "phi*phi"-part
energyDensity += stressU[ii][jj] * strainV[ii][jj];
// scale Gradient /* size_t row = localView.tree().child(k).localIndex(i); // kann auf Unterscheidung zwischen k & l verzichten?!
MatrixRT defGradientU(0); size_t col = localView.tree().child(l).localIndex(j); */ // siehe DUNE-book p.394
defGradientU[k][0] = gradients[i][0]; // Y // size_t col = localView.tree().child(k).localIndex(j); // Indizes mit k=l genügen .. Kroenecker-Delta_kl NEIN???
defGradientU[k][1] = gradients[i][1]; //X2 size_t col = localView.tree().child(k).localIndex(i); // kann auf Unterscheidung zwischen k & l verzichten?!
defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3 size_t row = localView.tree().child(l).localIndex(j);
elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
// symmetric Gradient (Elastic Strains) }
MatrixRT strainU, strainV;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
strainV[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // same ? genügt strainU
// strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
}
// St.Venant-Kirchhoff stress }
// < L sym[D_gamma*nabla phi_i], sym[D_gamma*nabla phi_i] >
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++) // "m*phi"-part
stressU[ii][ii] += lambda(quadPos) * trace; for (size_t l=0; l< nCompo; l++)
for (size_t j=0; j < nSf; j++ )
{
// 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] * strainU[ii][jj]; // "phi*phi"-part
// energyDensity += stressU[ii][jj] * strainV[ii][jj];
size_t row = localView.tree().child(k).localIndex(i); // kann auf Unterscheidung zwischen k & l verzichten?! // (scaled) Deformation gradient of the test basis function
// size_t col = localView.tree().child(l).localIndex(j); // siehe DUNE-book p.394 MatrixRT defGradientV(0);
size_t col = localView.tree().child(k).localIndex(j); // Indizes mit k=l genügen .. Kroenecker-Delta_kl defGradientV[l][0] = gradients[j][0]; // Y
defGradientV[l][1] = gradients[j][1]; //X2
defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
// symmetric Gradient (Elastic Strains)
MatrixRT strainV;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
{
strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
}
elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
for( size_t j=0; j<3; j++) for( size_t m=0; m<3; m++)
{ {
// < L G_alpha, sym[D_gamma*nabla phi_i] > // < L G_i, sym[D_gamma*nabla phi_j] >
// St.Venant-Kirchhoff stress // L G_i* strainV
MatrixRT stressG(0);
stressG.axpy(2*mu(quadPos), basisContainer[j]); //this += 2mu *strainU
double traceG = 0; // St.Venant-Kirchhoff stress
for (int ii=0; ii<nCompo; ii++) MatrixRT stressG(0);
traceG += basisContainer[j][ii][ii]; stressG.axpy(2*mu(quadPos), basisContainer[m]); //this += 2mu *strainU
for (int ii=0; ii<nCompo; ii++) double traceG = 0;
stressG[ii][ii] += lambda(quadPos) * traceG; for (int ii=0; ii<nCompo; ii++)
traceG += basisContainer[m][ii][ii];
double energyDensityGphi = 0; for (int ii=0; ii<nCompo; ii++)
for (int ii=0; ii<nCompo; ii++) stressG[ii][ii] += lambda(quadPos) * traceG;
for (int jj=0; jj<nCompo; jj++)
energyDensityGphi += stressG[ii][jj] * strainU[ii][jj]; // "phi*m"-part
double energyDensityGphi = 0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
energyDensityGphi += stressG[ii][jj] * strainV[ii][jj]; // "m*phi"-part
size_t row = localView.tree().child(l).localIndex(j);
auto value = energyDensityGphi * quadPoint.weight() * integrationElement; auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
elementMatrix[row][phiOffset+j] += value; elementMatrix[row][localPhiOffset+m] += value;
elementMatrix[phiOffset+j][row] += value; // elementMatrix[localPhiOffset+m][row] += value; // ---- reicht das ? --- TODO
}
}
// "phi*m"-part
for (size_t k=0; k < nCompo; k++)
for (size_t i=0; i < nSf; i++)
{
// (scaled) Deformation gradient of the ansatz basis function
MatrixRT defGradientU(0);
defGradientU[k][0] = gradients[i][0]; // Y
defGradientU[k][1] = gradients[i][1]; //X2
defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
// symmetric Gradient (Elastic Strains)
MatrixRT strainU;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
}
// St.Venant-Kirchhoff stress // St.Venant-Kirchhoff stress
// < L G_alpha, G_alpha > MatrixRT stressU(0);
for(size_t m=0; m<3; m++) stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
double trace = 0;
for (int ii=0; ii<nCompo; ii++)
trace += strainU[ii][ii];
for (int ii=0; ii<nCompo; ii++)
stressU[ii][ii] += lambda(quadPos) * trace;
for( size_t n=0; n<3; n++)
{ {
double energyDensityGG = 0;
// < L sym[D_gamma*nabla phi_i], G_j >
double energyDensityPhiG = 0;
for (int ii=0; ii<nCompo; ii++) for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++) for (int jj=0; jj<nCompo; jj++)
energyDensityGG += stressG[ii][jj] * basisContainer[m][ii][jj]; // "m*m"-part energyDensityPhiG += stressU[ii][jj] * basisContainer[n][ii][jj]; // "phi*m"-part
size_t col = localView.tree().child(k).localIndex(i);
elementMatrix[localPhiOffset+n][col] += energyDensityPhiG * quadPoint.weight() * integrationElement;
elementMatrix[phiOffset+m][phiOffset+j]= energyDensityGG * quadPoint.weight() * integrationElement;
} }
} }
}
// "m*m"-part
for(size_t m=0; m<3; m++)
for(size_t n=0; n<3; n++)
{
} // St.Venant-Kirchhoff stress
MatrixRT stressG(0);
stressG.axpy(2*mu(quadPos), basisContainer[m]); //this += 2mu *strainU
} double traceG = 0;
for (int ii=0; ii<nCompo; ii++)
traceG += basisContainer[m][ii][ii];
for (int ii=0; ii<nCompo; ii++)
stressG[ii][ii] += lambda(quadPos) * traceG;
double energyDensityGG = 0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
energyDensityGG += stressG[ii][jj] * basisContainer[n][ii][jj]; // "m*m"-part
elementMatrix[localPhiOffset+m][localPhiOffset+n]= energyDensityGG * quadPoint.weight() * integrationElement;
}
}
}
...@@ -520,12 +847,11 @@ void assembleCellProblem(const Basis& basis, ...@@ -520,12 +847,11 @@ void assembleCellProblem(const Basis& basis,
auto localView = basis.localView(); auto localView = basis.localView();
// std::cout << localView.maxSize() << std::endl; // std::cout << localView.maxSize() << std::endl;
const int phiOffset = basis.dimension(); // const int phiOffset = basis.dimension();
// Transform G_alpha's to GridViewFunctions/LocalFunctions -- why exactly? // Transform G_alpha's to GridViewFunctions/LocalFunctions -- why exactly?
auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView()); auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
auto loadFunctional = localFunction(loadGVF); // Dune::Functions:: notwendig? auto loadFunctional = localFunction(loadGVF); // Dune::Functions:: notwendig?
...@@ -539,9 +865,15 @@ void assembleCellProblem(const Basis& basis, ...@@ -539,9 +865,15 @@ void assembleCellProblem(const Basis& basis,
Dune::Matrix<double> elementMatrix; //eher das ?? Dune::Matrix<double> elementMatrix; //eher das ??
// Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix; // Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma, phiOffset); computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
printmatrix(std::cout, elementMatrix, "ElementMatrix", "--"); printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
Dune::Matrix<double> TestelementMatrix;
computeElementStiffnessMatrixOLD(localView, TestelementMatrix, muLocal, lambdaLocal, gamma);
printmatrix(std::cout, TestelementMatrix, "TESTElementMatrix", "--");
for (size_t i=0; i<elementMatrix.N(); i++) for (size_t i=0; i<elementMatrix.N(); i++)
{ {
// The global index of the i-th local degree of freedom of the current element // The global index of the i-th local degree of freedom of the current element
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment