add invert matrix by analytical expression with low dimension

NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp

@@ -260,14 +260,63 @@ void fullMatrixOperation::transpose(const fullMatrix<double>& A, fullMatrix<doub
     }
   }
 }
-bool fullMatrixOperation::invertMatrix(const fullMatrix<double>& A, fullMatrix<double>& invA, bool stiff, hyperFullMatrix& DinvADA)
+bool fullMatrixOperation::invertMatrix(const fullMatrix<double>& a, fullMatrix<double>& inva, bool stiff, hyperFullMatrix& DinvaDa)
 {
-  int N = A.size1();
-  allocateAndMakeZero(invA,N,N);
-  bool ok = A.invert(invA);
+  int N = a.size1();
+  allocateAndMakeZero(inva,N,N);
+  bool ok = true;
+  if (N==1)
+  {
+    inva(0,0) = 1./a(0,0);
+  }
+  else if (N == 2)
+  {
+    double det = a(0,0)*a(1,1) -a(0,1)*a(1,0); 
+    if(fabs(det) > 0.) 
+    {
+      double ud = 1. / det;
+      inva(0,0) = a(1,1) * ud;
+      inva(1,0) = -a(1,0) * ud;
+      inva(0,1) = -a(0,1) * ud;
+      inva(1,1) = a(0,0) * ud;
+    }
+    else 
+    {
+      ok = false;
+    }
+  }
+  else if (N == 3)
+  {
+    double detA =  (a(0,0) * (a(1,1) * a(2,2) - a(1,2) * a(2,1)) -
+            a(0,1) * (a(1,0) * a(2,2) - a(1,2) * a(2,0)) +
+            a(0,2) * (a(1,0) * a(2,1) - a(1,1) * a(2,0)));
+            
+    if (fabs(detA) > 0.)
+    {
+      double udet = 1./detA;
+      inva(0,0) =  (a(1,1) * a(2,2) - a(1,2) * a(2,1))* udet;
+      inva(1,0) = -(a(1,0) * a(2,2) - a(1,2) * a(2,0))* udet;
+      inva(2,0) =  (a(1,0) * a(2,1) - a(1,1) * a(2,0))* udet;
+      inva(0,1) = -(a(0,1) * a(2,2) - a(0,2) * a(2,1))* udet;
+      inva(1,1) =  (a(0,0) * a(2,2) - a(0,2) * a(2,0))* udet;
+      inva(2,1) = -(a(0,0) * a(2,1) - a(0,1) * a(2,0))* udet;
+      inva(0,2) =  (a(0,1) * a(1,2) - a(0,2) * a(1,1))* udet;
+      inva(1,2) = -(a(0,0) * a(1,2) - a(0,2) * a(1,0))* udet;
+      inva(2,2) =  (a(0,0) * a(1,1) - a(0,1) * a(1,0))* udet;
+    }
+    else
+    {
+      ok = false;
+    }
+  }
+  else
+  {
+    ok = a.invert(inva);
+  }
+  
   if (stiff)
   {
-    DinvADA.allocate(N,N,0.);
+    DinvaDa.allocate(N,N,0.);
     for (int k=0; k<N; k++)
     {
       for (int l=0; l<N; l++)
@@ -276,7 +325,7 @@ bool fullMatrixOperation::invertMatrix(const fullMatrix<double>& A, fullMatrix<d
         {
           for (int j=0; j<N; j++)
           {
-            DinvADA(i,j,k,l) = -0.5*(invA(i,k)*invA(j,l)+invA(i,l)*invA(j,k));
+            DinvaDa(i,j,k,l) = -0.5*(inva(i,k)*inva(j,l)+inva(i,l)*inva(j,k));
           }
         }
       }