diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset index b98b4f8781b32d2d57dbf2dd9e9f4c2326319bdd..f096251875be01605c9ca6628a47baa316c7eef5 100644 --- a/inputs/cellsolver.parset +++ b/inputs/cellsolver.parset @@ -68,8 +68,8 @@ alpha = 2.0 #--- Lame-Parameters mu1=1.0 -lambda1=0.0 -#lambda1=1.0 +#lambda1=0.0 +lambda1=1.0 # ---volume fraction (default value = 1.0/4.0) diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc index dbb84e936e410afad1a61000662d03606bbbd905..53b9688b190cbd8e3da43f8216d33ee7fb04c08b 100644 --- a/src/Cell-Problem.cc +++ b/src/Cell-Problem.cc @@ -194,12 +194,18 @@ void computeElementStiffnessMatrix(const LocalView& localView, for (size_t i=0; i<gradients.size(); i++) jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]); - for (size_t k=0; k < dimWorld; k++) - for (size_t i=0; i < nSf; i++) + for (size_t l=0; l< dimWorld; l++) + for (size_t j=0; j < nSf; j++ ) { + size_t row = localView.tree().child(l).localIndex(j); + // (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 // "phi*phi"-part - for (size_t l=0; l< dimWorld; l++) - for (size_t j=0; j < nSf; j++ ) + for (size_t k=0; k < dimWorld; k++) + for (size_t i=0; i < nSf; i++) { // (scaled) Deformation gradient of the ansatz basis function MatrixRT defGradientU(0); @@ -207,29 +213,24 @@ void computeElementStiffnessMatrix(const LocalView& localView, defGradientU[k][1] = gradients[i][1]; //X2 defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3 // printmatrix(std::cout, defGradientU , "defGradientU", "--"); - // (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 - + double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); // double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); // also works.. size_t col = localView.tree().child(k).localIndex(i); - size_t row = localView.tree().child(l).localIndex(j); + elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement; + } - // "m*phi" & "phi*m" - part - for( size_t m=0; m<3; m++) - { - double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV); - - auto value = energyDensityGphi * quadPoint.weight() * integrationElement; - elementMatrix[row][localPhiOffset+m] += value; - elementMatrix[localPhiOffset+m][row] += value; + // "m*phi" & "phi*m" - part + for( size_t m=0; m<3; m++) + { + double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV); + auto value = energyDensityGphi * quadPoint.weight() * integrationElement; + elementMatrix[row][localPhiOffset+m] += value; + elementMatrix[localPhiOffset+m][row] += value; } - } + } // "m*m"-part for(size_t m=0; m<3; m++) @@ -418,10 +419,13 @@ void assembleCellStiffness(const Basis& basis, muLocal.bind(element); lambdaLocal.bind(element); const int localPhiOffset = localView.size(); - +// Dune::Timer Time; //std::cout << "localPhiOffset : " << localPhiOffset << std::endl; Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix; computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma); + +// std::cout << "local assembly time:" << Time.elapsed() << std::endl; + //printmatrix(std::cout, elementMatrix, "ElementMatrix", "--"); //std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl; //std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl; @@ -848,6 +852,8 @@ double computeL2Norm(const Basis& basis, int main(int argc, char *argv[]) { MPIHelper::instance(argc, argv); + + Dune::Timer globalTimer; ParameterTree parameterSet; if (argc < 2) @@ -1875,4 +1881,6 @@ int main(int argc, char *argv[]) log << std::string(tableWidth*7 + 3*7, '-') << "\n"; log.close(); + +// std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl; }