From 47411592a01ae03027f666eba11c73c0bb1633bc Mon Sep 17 00:00:00 2001
From: Eric Bechet <eric.bechet@ulg.ac.be>
Date: Fri, 20 Jan 2012 12:43:32 +0000
Subject: [PATCH]

---
 Solver/CMakeLists.txt       |  6 ++--
 Solver/STensor43.h          | 61 +++++++++++++++++++++++++++++++++++++
 Solver/elasticitySolver.cpp | 30 ++++++++++++++++++
 Solver/elasticitySolver.h   |  1 +
 4 files changed, 95 insertions(+), 3 deletions(-)

diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index bd74ac25ea..9f64e02cb5 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 7ea38c1f24..351c39657f 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 6da2468237..97dcc48135 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 0731c4d2f9..3ae30a7bdc 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);
-- 
GitLab