No commit message

No commit message

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

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()
 {
   

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);