From d7e19c9b94f9995cb9229d49f607b54626e5af46 Mon Sep 17 00:00:00 2001
From: Klaus Boehnlein <klaus.boehnlein@tu-dresden.de>
Date: Tue, 15 Jun 2021 07:45:36 +0200
Subject: [PATCH] fixed local Element assembly
---
src/dune-microstructure.cc | 547 ++++++++++++++++++-------------------
1 file changed, 268 insertions(+), 279 deletions(-)
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 5fb0bec0..1edc7730 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -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", "--");
--
GitLab