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

fixed local Element assembly

parent da401b8c
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/dune-microstructure.cc
+ 268
− 279
View file @ d7e19c9b
......@@ -56,225 +56,219 @@ using namespace Dune;
template<class LocalView, class Matrix, class localFunction1, class localFunction2>
void computeElementStiffnessMatrixOLD(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();
constexpr int dim = Element::dimension;
constexpr int nCompo = dim;
using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
// using Domain = typename LocalView::GridView::Codim<0>::Geometry::GlobalCoordinate;
// using FuncScalar = std::function< ScalarRT(const Domain&) >;
// using Func2Tensor = std::function< 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;
// template<class LocalView, class Matrix, class localFunction1, class localFunction2>
// void computeElementStiffnessMatrixOLD(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();
//
// constexpr int dim = Element::dimension;
// constexpr int nCompo = dim;
//
// using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
// // using Domain = typename LocalView::GridView::Codim<0>::Geometry::GlobalCoordinate;
// // using FuncScalar = std::function< ScalarRT(const Domain&) >;
// // using Func2Tensor = std::function< 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++)
// energyDensityPhiG += stressU[ii][jj] * basisContainer[m][ii][jj]; // "phi*m"-part
// {
// 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] * strainV[ii][jj]; // "phi*phi"-part
//
// /* 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;
// }
//
//
//
// }
//
//
//
// }
//
//
// }
//
//
//
//
// }
//
//
// }
//
// 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>
......@@ -408,13 +402,13 @@ void computeElementStiffnessMatrix(const LocalView& localView,
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];
energyDensity += stressU[ii][jj] * strainV[ii][jj]; // "phi*phi"-part
// energyDensity += stressU[ii][jj] * strainU[ii][jj];
/* size_t row = localView.tree().child(k).localIndex(i); // kann auf Unterscheidung zwischen k & l verzichten?!
size_t col = localView.tree().child(l).localIndex(j); */ // siehe DUNE-book p.394
// size_t col = localView.tree().child(k).localIndex(j); // Indizes mit k=l genügen .. Kroenecker-Delta_kl NEIN???
size_t col = localView.tree().child(k).localIndex(i); // kann auf Unterscheidung zwischen k & l verzichten?!
// size_t row = localView.tree().child(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;
......@@ -424,7 +418,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
}
// "m*phi"-part
// "m*phi" & "phi*m" -part
for (size_t l=0; l< nCompo; l++)
for (size_t j=0; j < nSf; j++ )
{
......@@ -472,7 +466,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
elementMatrix[row][localPhiOffset+m] += value;
// elementMatrix[localPhiOffset+m][row] += value; // ---- reicht das ? --- TODO
elementMatrix[localPhiOffset+m][row] += value; // ---- reicht das ? --- TODO
}
......@@ -481,60 +475,59 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// "phi*m"-part
for (size_t k=0; k < nCompo; k++)
for (size_t i=0; i < nSf; i++)
{
// (scaled) Deformation gradient of the ansatz basis function
MatrixRT defGradientU(0);
defGradientU[k][0] = gradients[i][0]; // Y
defGradientU[k][1] = gradients[i][1]; //X2
defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
// symmetric Gradient (Elastic Strains)
MatrixRT strainU;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
{
strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
}
// St.Venant-Kirchhoff stress
MatrixRT stressU(0);
stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
double trace = 0;
for (int ii=0; ii<nCompo; ii++)
trace += strainU[ii][ii];
for (int ii=0; ii<nCompo; ii++)
stressU[ii][ii] += lambda(quadPos) * trace;
for( size_t n=0; n<3; n++)
{
// < L sym[D_gamma*nabla phi_i], G_j >
double energyDensityPhiG = 0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
energyDensityPhiG += stressU[ii][jj] * basisContainer[n][ii][jj]; // "phi*m"-part
size_t col = localView.tree().child(k).localIndex(i);
elementMatrix[localPhiOffset+n][col] += energyDensityPhiG * quadPoint.weight() * integrationElement;
}
}
// for (size_t k=0; k < nCompo; k++)
// for (size_t i=0; i < nSf; i++)
// {
//
// // (scaled) Deformation gradient of the ansatz basis function
// MatrixRT defGradientU(0);
// defGradientU[k][0] = gradients[i][0]; // Y
// defGradientU[k][1] = gradients[i][1]; //X2
// defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
//
//
// // symmetric Gradient (Elastic Strains)
// MatrixRT strainU;
// for (int ii=0; ii<nCompo; ii++)
// for (int jj=0; jj<nCompo; jj++)
// {
// strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
// }
//
// // St.Venant-Kirchhoff stress
// MatrixRT stressU(0);
// stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
//
// double trace = 0;
// for (int ii=0; ii<nCompo; ii++)
// trace += strainU[ii][ii];
//
// for (int ii=0; ii<nCompo; ii++)
// stressU[ii][ii] += lambda(quadPos) * trace;
//
// for( size_t n=0; n<3; n++)
// {
//
// // < L sym[D_gamma*nabla phi_i], G_j >
// double energyDensityPhiG = 0;
// for (int ii=0; ii<nCompo; ii++)
// for (int jj=0; jj<nCompo; jj++)
// energyDensityPhiG += stressU[ii][jj] * basisContainer[n][ii][jj]; // "phi*m"-part
//
// size_t col = localView.tree().child(k).localIndex(i);
//
// elementMatrix[localPhiOffset+n][col] += energyDensityPhiG * quadPoint.weight() * integrationElement;
//
// }
//
//
//
// }
// "m*m"-part
for(size_t m=0; m<3; m++)
for(size_t n=0; n<3; n++)
{
......@@ -560,10 +553,6 @@ void computeElementStiffnessMatrix(const LocalView& localView,
}
}
......@@ -868,9 +857,9 @@ void assembleCellProblem(const Basis& basis,
computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
Dune::Matrix<double> TestelementMatrix;
computeElementStiffnessMatrixOLD(localView, TestelementMatrix, muLocal, lambdaLocal, gamma);
printmatrix(std::cout, TestelementMatrix, "TESTElementMatrix", "--");
// Dune::Matrix<double> TestelementMatrix;
// computeElementStiffnessMatrixOLD(localView, TestelementMatrix, muLocal, lambdaLocal, gamma);
// printmatrix(std::cout, TestelementMatrix, "TESTElementMatrix", "--");
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment