diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc index 1edc7730716059823efc9a8300a8c5d3d50f9ae4..f4ee1a89401ebc8ca9e9a2470b3b33f2b30bc045 100644 --- a/src/dune-microstructure.cc +++ b/src/dune-microstructure.cc @@ -56,220 +56,6 @@ 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] * 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; -// } -// -// -// -// } -// -// -// -// } -// -// -// } -// -// -// -// -// } -// -// -// } -// -// - template<class LocalView, class Matrix, class localFunction1, class localFunction2> void computeElementStiffnessMatrix(const LocalView& localView, @@ -617,13 +403,14 @@ void computeElementStiffnessMatrix(const LocalView& localView, template<class Basis> void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb) { - + // OccupationPattern: + // | phi*phi m*phi | + // | phi *m m*m | auto gridView = basis.gridView(); auto localView = basis.localView(); - int maxCounter = 0; - const int maxPhiIndex = basis.dimension(); + const int phiOffset = basis.dimension(); nb.resize(basis.size()+3, basis.size()+3); // ??? @@ -641,45 +428,28 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb) nb.add(row[0],col[0]); // nun würde auch nb.add(row,col) gehen.. } - - if (row > maxCounter) - maxCounter = row; - } - - - - // Extend the OccupationPattern: - // | phi*phi m*phi | - // | phi *m m*m | - for (size_t i=0; i<localView.size(); i++) { - auto PhiIndex = localView.index(i); + auto row = localView.index(i); for (size_t j=0; j<3; j++ ) { - nb.add(PhiIndex,maxPhiIndex+j); // m*phi part of StiffnessMatrix + nb.add(row,phiOffset+j); // m*phi part of StiffnessMatrix - nb.add(maxPhiIndex+j,PhiIndex); // phi*m part of StiffnessMatrix + nb.add(phiOffset+j,row); // phi*m part of StiffnessMatrix } } - for (size_t i=0; i<3; i++ ) for (size_t j=0; j<3; j++ ) { - nb.add(maxPhiIndex+i,maxPhiIndex+j); // m*m part of StiffnessMatrix + nb.add(phiOffset+i,phiOffset+j); // m*m part of StiffnessMatrix } } - - - std::cout << "maxCounter: " << maxCounter << std::endl; - std::cout << "maxPhiIndex: " << maxPhiIndex << std::endl; - } @@ -820,7 +590,6 @@ void assembleCellProblem(const Basis& basis, { std::cout << "assemble stiffnessmatrix begins." << std::endl; - MatrixIndexSet occupationPattern; getOccupationPattern(basis, occupationPattern); occupationPattern.exportIdx(matrix); @@ -836,7 +605,8 @@ void assembleCellProblem(const Basis& basis, auto localView = basis.localView(); // 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? @@ -851,31 +621,68 @@ void assembleCellProblem(const Basis& basis, muLocal.bind(element); lambdaLocal.bind(element); + const int localPhiOffset = localView.size(); + std::cout << "localPhiOffset : " << localPhiOffset << std::endl; Dune::Matrix<double> elementMatrix; //eher das ?? // Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix; computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma); printmatrix(std::cout, elementMatrix, "ElementMatrix", "--"); + +// std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl; +// std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl; + // Dune::Matrix<double> TestelementMatrix; // computeElementStiffnessMatrixOLD(localView, TestelementMatrix, muLocal, lambdaLocal, gamma); // printmatrix(std::cout, TestelementMatrix, "TESTElementMatrix", "--"); - - - for (size_t i=0; i<elementMatrix.N(); i++) + ////////////////////////////////////////////////////////////////////////////// + // GLOBAL ASSEMBLY ..(hier dreht sich i, j..) + ////////////////////////////////////////////////////////////////////////////// + for (size_t i=0; i<localPhiOffset i++) { // The global index of the i-th local degree of freedom of the current element auto row = localView.index(i); - for (size_t j=0; j<elementMatrix.M(); j++ ) + + for (size_t j=0; j<localPhiOffset; j++ ) { // The global index of the j-th local degree of freedom of the current element auto col = localView.index(j); - matrix[row][col] += elementMatrix[i][j]; // keine periodicMap nötig..wird durch Basis behandelt.. + matrix[row][col] += elementMatrix[i][j]; // keine periodicMap nötig..wird durch Basis behandelt.. } } + for (size_t i=0; i<localPhiOffset; i++) + for(size_t m=0; m<3; m++) + { + auto row = localView.index(i); + + matrix[row][phiOffset+m] += elementMatrix[i][localPhiOffset+m]; + matrix[phiOffset+m][row] += elementMatrix[localPhiOffset+m][i]; + } + + + for (size_t m=0; m<3; m++ ) + for (size_t n=0; n<3; n++ ) + { + matrix[phiOffset+m][phiOffset+n] += elementMatrix[localPhiOffset+m][localPhiOffset+n]; + } + + + + + + + + + + + /* + + + ///////////////////////////////////////////////////////////////////////////////////// // Now get the local contribution to the right-hand side vector // BlockVector<double> localB; BlockVector<FieldVector<double,1>> elementRhs; @@ -888,6 +695,8 @@ void assembleCellProblem(const Basis& basis, b[row] += elementRhs[p]; } + */ + } @@ -928,6 +737,7 @@ void assembleCellProblem(const Basis& basis, return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0])) and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1])) and (FloatCmp::eq(x[2],y[2]))); +// and (FloatCmp::eq(x[2],y[2]) or FloatCmp::eq(x[2]+1,y[2]) or FloatCmp::eq(x[2]-1,y[2]))); }; @@ -996,7 +806,7 @@ int main(int argc, char *argv[]) FieldVector<double,dim> lower({0.0, 0.0, -1.0/2.0}); FieldVector<double,dim> upper({1.0, 1.0, 1.0/2.0}); - int nE = 1; + int nE = 2; std::array<int,dim> nElements={nE,nE,nE}; //#Elements in each direction @@ -1138,7 +948,7 @@ Func2Tensor x3G_1 = [] (const Domain& z) { printmatrix(std::cout, basisContainer[0] , "G_1", "--"); printmatrix(std::cout, basisContainer[1] , "G_2", "--"); - printmatrix(std::cout, basisContainer[2] , "G_3", "--"); + printmatrix(std::cout, basisContainer[2] , "´G_3", "--"); @@ -1148,7 +958,7 @@ assembleCellProblem(Basis_CE, muLocal, lambdaLocal,gamma, stiffnessMatrix_CE ,l -// printmatrix(std::cout, stiffnessMatrix, "StiffnessMatrix before B.C", "--"); +printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix before B.C", "--");