diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt index bd74ac25eabdedb8d7d01294e6cc6221ec90a2a0..9f64e02cb5a07a16ecdef2b6b810ccef08046efe 100644 --- a/Solver/CMakeLists.txt +++ b/Solver/CMakeLists.txt @@ -10,15 +10,15 @@ set(SRC dofManager.cpp groupOfElements.cpp elasticityTerm.cpp - elasticitySolver.cpp +elasticitySolver.cpp SElement.cpp eigenSolver.cpp multiscaleLaplace.cpp - functionSpace.cpp +functionSpace.cpp filters.cpp sparsityPattern.cpp STensor43.cpp - terms.cpp +terms.cpp ) file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h) diff --git a/Solver/STensor43.h b/Solver/STensor43.h index 7ea38c1f24b5f46b2265c4cf9200ee5ec0d3e8ed..351c39657fae03aa4b94327e976f469d3069ef43 100644 --- a/Solver/STensor43.h +++ b/Solver/STensor43.h @@ -69,6 +69,67 @@ class STensor43 { return *this; }*/ void print(const char *) const; + + STensor43 transpose (int n, int m) const + { + STensor43 ithis; + if ((n==0 && m==1) || (n==1 && m==0)) + { + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) + ithis(i,j,k,l) = (*this)(j,i,k,l); + return ithis; + } + if ((n==0 && m==2) || (n==2 && m==0)) + { + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) + ithis(i,j,k,l) = (*this)(k,j,i,l); + return ithis; + } + if ((n==0 && m==3) || (n==3 && m==0)) + { + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) + ithis(i,j,k,l) = (*this)(l,j,k,i); + return ithis; + } + if ((n==1 && m==2) || (n==2 && m==1)) + { + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) + ithis(i,j,k,l) = (*this)(i,k,j,l); + return ithis; + } + if ((n==1 && m==3) || (n==3 && m==1)) + { + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) + ithis(i,j,k,l) = (*this)(i,l,k,j); + return ithis; + } + if ((n==2 && m==3) || (n==3 && m==2)) + { + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) + ithis(i,j,k,l) = (*this)(i,j,l,k); + return ithis; + } + return ithis+=(*this); + } + }; // tensor product diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp index 6da24682378d648457a044348f0581df96db9764..97dcc481359271f6bfd068720b00b0c76c81c27e 100644 --- a/Solver/elasticitySolver.cpp +++ b/Solver/elasticitySolver.cpp @@ -65,6 +65,36 @@ void elasticitySolver::setMesh(const std::string &meshFileName) LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpace(_tag+1); } + +void elasticitySolver::exportKb() +{ + + FILE *f = fopen ( "K.txt", "w" ); + double valeur; + std::string sysname = "A"; + for ( int i = 0 ; i < pAssembler->sizeOfR() ; i ++ ) + { + for ( int j = 0 ; j < pAssembler->sizeOfR() ; j ++ ) + { + pAssembler->getLinearSystem ( sysname )->getFromMatrix ( i,j, valeur ); + fprintf ( f,"%+e ",valeur ) ; + } + fprintf ( f,"\n" ); + } + + fclose ( f ); + + f = fopen ( "b.txt", "w" ); + for ( int i = 0 ; i < pAssembler->sizeOfR() ; i ++ ) + { + pAssembler->getLinearSystem ( sysname )->getFromRightHandSide ( i,valeur ); + fprintf ( f,"%+e\n",valeur ) ; + } + + fclose ( f ); + +} + void elasticitySolver::solve() { diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h index 0731c4d2f9a547df3bffc8b91f0452909b90448b..3ae30a7bdc8263decf963699cbdc5153b24611cf 100644 --- a/Solver/elasticitySolver.h +++ b/Solver/elasticitySolver.h @@ -93,6 +93,7 @@ class elasticitySolver virtual void setMesh(const std::string &meshFileName); void solve(); void postSolve(); + void exportKb(); void getSolutionOnElement(MElement *el, fullMatrix<double> &sol); virtual PView *buildDisplacementView(const std::string postFileName); virtual PView *buildStressesView(const std::string postFileName);