add new function

NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp

@@ -1215,7 +1215,170 @@ void LaminateHomogenization::compute(const std::vector<fullMatrix<double> >& Cin
       DCeffDcin[0]=AA;
     }
   }
+};
+
+
+void LaminateHomogenization::compute(const std::vector<fullMatrix<double> >& Cin, const std::vector<double>& f, 
+                        const std::vector<SVector3>& normal,
+                        fullMatrix<double>& Ceff, 
+                        bool stiff, 
+                        std::vector<hyperFullMatrix>& DCeffDcin, 
+                        std::vector<fullMatrix<double> >& DCeffDf, 
+                        std::vector<std::vector<fullMatrix<double> > >&DCeffDnormal) const
+{
+  int numPhase = Cin.size();
+  if (stiff)
+  {
+    if (DCeffDcin.size() != numPhase)
+    {
+      DCeffDcin.resize(numPhase);
+    }
+    if (DCeffDf.size() != numPhase)
+    {
+      DCeffDf.resize(numPhase);
+    }
+    if (DCeffDnormal.size() != normal.size())
+    {
+      DCeffDnormal.resize(normal.size(), std::vector<fullMatrix<double> >(3));
+    }
+  };
+  
+  if (numPhase ==1) 
+  {
+    Ceff = Cin[0];
+    if (stiff)
+    {
+      DCeffDcin[0].allocate(Ceff.size1(),Cin[0].size1(),1.);
+      fullMatrixOperation::allocateAndMakeZero(DCeffDf[0],Cin[0].size1(),Cin[0].size2());
+      for (int j=0; j< normal.size(); j++)
+      {
+        fullMatrixOperation::allocateAndMakeZero(DCeffDnormal[j][0],Cin[0].size1(),Cin[0].size2());
+        fullMatrixOperation::allocateAndMakeZero(DCeffDnormal[j][1],Cin[0].size1(),Cin[0].size2());
+        fullMatrixOperation::allocateAndMakeZero(DCeffDnormal[j][2],Cin[0].size1(),Cin[0].size2());        
+      }
+    }
+  }
+  else if (numPhase == 2)
+  {
+    computeBimaterial(Cin[0],Cin[1],f[0],f[1],normal[0],Ceff,
+                    stiff,&DCeffDcin[0],&DCeffDcin[1],
+                    &DCeffDf[0],&DCeffDf[1],&DCeffDnormal[0][0],
+                    &DCeffDnormal[0][1],&DCeffDnormal[0][2]);
   }
+  else
+  {
+    //
+    auto getSumFaction = [](const std::vector<double>& ff, int i, fullVector<double>& DoutDf)
+    {
+      DoutDf.resize(ff.size());
+      DoutDf.setAll(0.);
+      double v= 0;
+      for (int k=0; k< i+1; k++)
+      {
+        v += ff[k];
+        DoutDf(k) = 1.;
+      }
+      return v;
+    };
+    
+    static std::vector<hyperFullMatrix> A, B;
+    static std::vector<fullMatrix<double> > CC, F0, F1, N0, N1, N2;
+    static fullMatrix<double> Df0iDf, Df1iDf;
+    
+    if (CC.size() != numPhase)
+    {
+      CC.resize(numPhase);
+    };
+    //
+    if (stiff)
+    {
+      // resallocate for static variable
+      if (A.size() != numPhase-1) A.resize(numPhase-1);
+      if (B.size() != numPhase-1) B.resize(numPhase-1);
+      if (F0.size() != numPhase-1) F0.resize(numPhase-1);
+      if (F1.size() != numPhase-1) F1.resize(numPhase-1);
+      if (N0.size() != numPhase-1) N0.resize(numPhase-1);
+      if (N1.size() != numPhase-1) N1.resize(numPhase-1);
+      if (N2.size() != numPhase-1) N2.resize(numPhase-1);
+      
+      fullMatrixOperation::allocateAndMakeZero(Df0iDf,numPhase-1,numPhase);
+      fullMatrixOperation::allocateAndMakeZero(Df1iDf,numPhase-1,numPhase);
+      
+      for (int j=0; j< numPhase; j++)
+      {
+        fullMatrixOperation::allocateAndMakeZero(DCeffDf[j],Cin[0].size1(),Cin[0].size2());
+      }
+    }
+    //
+    CC[0] = Cin[0];
+    for (int i=0; i< numPhase-1; i++)
+    {
+      static fullVector<double> DftolDf, DftolPrevDf;
+      double ftol = getSumFaction(f,i+1,DftolDf);
+      double ftolPrev = getSumFaction(f,i,DftolPrevDf);
+      
+      double f0i = ftolPrev/ftol;
+      double f1i = 1.-f0i;
+      if (stiff)
+      {     
+        for (int r=0; r< numPhase; r++)
+        {
+          Df0iDf(i,r) = (-ftolPrev/ftol/ftol)*DftolDf(r) + (1./ftol)*DftolPrevDf(r);
+          Df1iDf(i,r) = -Df0iDf(i,r);
+        };
+      }
+      // two phase
+      computeBimaterial(CC[i],Cin[i+1],f0i,f1i,normal[i],CC[i+1],stiff,&A[i],&B[i],&F0[i],&F1[i],&N0[i],&N1[i],&N2[i]);
+    };
+    //
+    Ceff = CC[numPhase-1];
+    //
+    if (stiff)
+    {
+      
+      DCeffDcin[numPhase-1] = B[numPhase-2];
+      //
+      for (int j=0; j< numPhase; j++)
+      {
+        DCeffDf[j].axpy(F0[numPhase-2],Df0iDf(numPhase-2,j));
+        DCeffDf[j].axpy(F1[numPhase-2],Df1iDf(numPhase-2,j));
+      }
+      //
+      DCeffDnormal[numPhase-2][0] = N0[numPhase-2];
+      DCeffDnormal[numPhase-2][1] = N1[numPhase-2];
+      DCeffDnormal[numPhase-2][2] = N2[numPhase-2];
+      
+      static hyperFullMatrix AA;
+      AA = A[numPhase-2];
+      for (int i=numPhase-2; i>0; i--)
+      {        
+        fullMatrixOperation::mult(AA,B[i-1],DCeffDcin[i],1.);
+        
+        static fullMatrix<double> AAF0, AAF1;
+        fullMatrixOperation::mult(AA,F0[i-1],AAF0);
+        fullMatrixOperation::mult(AA,F1[i-1],AAF1);
+        //
+        for (int j=0; j< numPhase; j++)
+        {
+          DCeffDf[j].axpy(AAF0,Df0iDf(i-1,j));
+          DCeffDf[j].axpy(AAF1,Df1iDf(i-1,j));
+        }
+        //
+        
+        fullMatrixOperation::mult(AA,N0[i-1],DCeffDnormal[i-1][0]);
+        fullMatrixOperation::mult(AA,N1[i-1],DCeffDnormal[i-1][1]);
+        fullMatrixOperation::mult(AA,N2[i-1],DCeffDnormal[i-1][2]);
+        //
+        hyperFullMatrix temp;
+        fullMatrixOperation::mult(AA,A[i-1],temp,1.);
+        AA = temp;
+      }
+      DCeffDcin[0]=AA;
+    }
+  }
+};
+                        
+                        
 void LaminateHomogenization::printInfos() const
 {
   printf("using Voigt, Reuss with laminate \n");

NonLinearSolver/modelReduction/DeepMaterialNetworks.h

@@ -167,6 +167,14 @@ class LaminateHomogenization : public BimaterialHomogenization
                         std::vector<fullMatrix<double>>& DCeffDf, 
                         std::vector<fullMatrix<double>>&DCeffDnormal) const;
                         
+    virtual void compute(const std::vector<fullMatrix<double> >& Cin, const std::vector<double>& f, 
+                        const std::vector<SVector3>& normal,
+                        fullMatrix<double>& Ceff, 
+                        bool stiff, 
+                        std::vector<hyperFullMatrix>& DCeffDcin, 
+                        std::vector<fullMatrix<double> >& DCeffDf, 
+                        std::vector<std::vector<fullMatrix<double> > >&DCeffDnormal) const; 
+    
     virtual void printInfos() const;
     #endif //SWIG
 };