diff --git a/NonLinearSolver/materialLaw/mlawJ2VMSmallStrain.cpp b/NonLinearSolver/materialLaw/mlawJ2VMSmallStrain.cpp
index 09ec824d6d808e4d5982078526cca864d1a4c323..7f6d1fce75bcdfa0837a6fed3dada2a1d1ce1441 100644
--- a/NonLinearSolver/materialLaw/mlawJ2VMSmallStrain.cpp
+++ b/NonLinearSolver/materialLaw/mlawJ2VMSmallStrain.cpp
@@ -322,8 +322,6 @@ void mlawJ2VMSmallStrain::predictorCorector(const STensor3& F0,         // initi
   Ep = q0->_j2lepsp; // plastic strain tensor (equal I if no plastic occurs)
   gamma = q0->_j2lepspbarre; // equivalent plastic strain
   
-  eigenstrain = Ep;
-  
   static const STensor3 I(1.); // identity tensor
 
   // small strain
@@ -417,7 +415,6 @@ void mlawJ2VMSmallStrain::predictorCorector(const STensor3& F0,         // initi
 
     Ep.daxpy(N,deps);
     q1->_Ee.daxpy(N,-deps);
-    eigenstrain = Ep;
     
     double fact = -2.*_mu*deps;
     P.daxpy(N,fact);
@@ -498,6 +495,9 @@ void mlawJ2VMSmallStrain::predictorCorector(const STensor3& F0,         // initi
       }
     }
   }
+  
+  eigenstrain = Ep;
+  eigenstrain -= I;
 
   q1->_elasticEnergy = deformationEnergy(q1->_Ee);
 
diff --git a/NonLinearSolver/modelReduction/Clustering.cpp b/NonLinearSolver/modelReduction/Clustering.cpp
index 5e01a5f430e12dbe8752307696029dc3ba314361..68f61e4be3653c3e82b51e09053ef4ae6d5c9780 100644
--- a/NonLinearSolver/modelReduction/Clustering.cpp
+++ b/NonLinearSolver/modelReduction/Clustering.cpp
@@ -3321,6 +3321,197 @@ bool Clustering::computeClustersPlasticity(int minNumberElementEachCluster, int
   return true;
 };
 
+bool Clustering::computeClustersLocalOrientation(int minNumberElementEachCluster, int maxNbIterations,
+                                                 bool saveToFile, std::string clusterFileName,
+                                                 bool saveSummaryFile, std::string summaryFileName,
+                                                 bool saveClusterMeanFiberOrientationFile, const std::string saveClusterMeanFiberOrientationFileName
+                                                )
+{
+  // init cluster 
+  clearAllClusterData();
+  for (std::map<int,int>::const_iterator it = _materialMap.begin(); it != _materialMap.end(); it++)
+  {
+    _clusterIndex[it->first] = -1;
+  }
+  
+  //
+  _totalNumberOfClusters = 0;
+  for (std::map<int, std::map<int, fullMatrix<double> > >::const_iterator  itMat = _concentrationTensorMap.begin(); 
+                      itMat != _concentrationTensorMap.end(); itMat++)
+  {
+    int matNum = itMat->first;
+    int nbCluster = 0;
+    std::map<int,int>::const_iterator itnb = _numberOfClustersPerLaw.find(matNum);
+    if (itnb !=_numberOfClustersPerLaw.end())
+    {
+      nbCluster = itnb->second;
+    }
+    if (nbCluster > 0)
+    {
+      const std::map<int, fullMatrix<double> >& dataInMat = itMat->second;
+      const std::map<int, fullVector<double> >& dataInMatInelasticStrainConc = _inelasticStrainConcentrationVectorMap[itMat->first];      
+      const std::map<int, fullVector<double> >& dataInMatLocOri = _localOrientationVectorMap[itMat->first];
+      
+      std::map<int, fullVector<double> >& dataInMatInelasticStrainConcNormalized = _inelasticStrainConcentrationNormalizedVectorMap[itMat->first];
+      clutering_impls::vecNormalizeByOverallNorm(dataInMatInelasticStrainConc,_mapGaussPoints,dataInMatInelasticStrainConcNormalized);
+   
+      //Msg::Info("start clustering for material %d \n number of clusters = %d, min number in each cluster = %d",matNum, nbCluster,minNumberElementEachCluster);
+      std::map<int, std::map<int, const fullMatrix<double>*> > clusterInMat;
+      std::map<int, std::map<int, const fullVector<double>*> > clusterInMatInelasticTensor;      
+      std::map<int, std::map<int, const fullVector<double>*> > clusterInMatLocOri;
+
+      bool ok=false;
+                                                             
+      //////// for woven /////////
+      if(matNum==11)
+      {                                                        
+      ok= clutering_impls::kmean_cluster_FullVectorBased(dataInMatInelasticStrainConcNormalized,dataInMat,nbCluster,minNumberElementEachCluster,maxNbIterations,_mapGaussPoints,
+                                                              clusterInMatInelasticTensor,clusterInMat);
+      }
+      if(matNum==12 or matNum==13)
+      {                                                        
+      //ok= clutering_impls::kmean_cluster_orientationBased(dataInMatLocOri,dataInMat,nbCluster,minNumberElementEachCluster,maxNbIterations,_mapGaussPoints,
+        //                                                       clusterInMatLocOri,clusterInMat);
+        
+      ok= clutering_impls::kmean_cluster_hierarchical_orientation_fullVectorBased(dataInMatLocOri,dataInMatInelasticStrainConcNormalized,dataInMat,
+                                                                                  nbCluster,minNumberElementEachCluster,maxNbIterations,_mapGaussPoints,
+                                                                                  clusterInMatLocOri,clusterInMatInelasticTensor,clusterInMat);
+                                                                                  
+      //ok= clutering_impls::kmean_cluster_FullVectorBased(dataInMatInelasticStrainConcNormalized,dataInMat,nbCluster,minNumberElementEachCluster,maxNbIterations,_mapGaussPoints,
+        //                                                      clusterInMatInelasticTensor,clusterInMat);
+      }
+      
+      if (!ok)
+      {
+        return false;
+      }
+      for (std::map<int, std::map<int, const fullMatrix<double>*> >::const_iterator it = clusterInMat.begin(); it!= clusterInMat.end(); it++)
+      {
+        int clusterNum = it->first + _totalNumberOfClusters;
+        const std::map<int, const fullMatrix<double>*>& cluster = it->second;
+        //
+        for (std::map<int, const fullMatrix<double>*>::const_iterator itset = cluster.begin(); itset != cluster.end(); itset++)
+        {
+          // update cluster index at each material point
+          _clusterIndex[itset->first] = clusterNum;
+        }
+        // average concentration tensor in this cluster
+        static fullMatrix<double> Amean(6,6);
+        _clusterVolumeMap[clusterNum] = clutering_impls::meanCluster(cluster,_mapGaussPoints,Amean);
+        STensorOperation::fromFullMatrixToSTensor43(Amean,_clusterAverageStrainConcentrationTensorMap[clusterNum]);
+                                            
+        _clusterMaterialMap[clusterNum] = matNum;
+        
+        Msg::Info("cluster %d, matNum = %d, volume = %e, stdev %e",clusterNum,matNum,_clusterVolumeMap[clusterNum],0.);
+        Amean.print("average strain concentration tensor");
+                
+        fullVector<double> cluster_fiberOrientation(3);
+        double E1_mean(0.), E2_mean(0.), E3_mean(0.);
+        clutering_impls::meanClusterFiberOrientation(cluster,_mapGaussPoints,_mapLocalOrientation,E1_mean,E2_mean,E3_mean);
+        cluster_fiberOrientation(0) = E1_mean;
+        cluster_fiberOrientation(1) = E2_mean;
+        cluster_fiberOrientation(2) = E3_mean;
+        _clusterMeanFiberOrientationMap[clusterNum] = cluster_fiberOrientation;
+        Msg::Info("mean fiber orientation of cluster = %d is %.16g,%.16g,%.16g ",clusterNum,E1_mean,E2_mean,E3_mean);
+               
+      }
+      _totalNumberOfClusters += nbCluster;
+      Msg::Info("done clustering for material %d",matNum);
+    }
+    else
+    {
+      Msg::Info("material %d is not clustered",matNum);
+    }
+  }
+  
+  if (saveClusterMeanFiberOrientationFile)
+  {
+    std::vector<std::string> splitFileName = SplitFileName(saveClusterMeanFiberOrientationFileName);
+    std::string fname = splitFileName[1]+"_nbClusters_"+std::to_string(_totalNumberOfClusters)+splitFileName[2];
+    FILE* f = fopen(fname.c_str(),"w");
+    printf("start writing cluster mean fiber orientation to file: %s\n",fname.c_str());
+    std::string str ="CLUSTER_ID MAT_NUM E1 E2 E3";
+    fprintf(f,"%s\n",str.c_str());
+    for (std::map<int, fullVector<double> >::const_iterator it = _clusterMeanFiberOrientationMap.begin();  it!= _clusterMeanFiberOrientationMap.end(); it++)
+    { 
+      int clusterId = it->first;
+      int matNum = _clusterMaterialMap[clusterId];
+      fullVector<double> angles = it->second;
+      double E1 = angles(0);
+      double E2 = angles(1);
+      double E3 = angles(2);
+      fprintf(f, "%d %d %.16g %.16g %.16g\n",clusterId,matNum,E1,E2,E3);
+    }
+    fclose(f);
+    printf("done writing cluster mean fiber orientation to file: %s \n",fname.c_str());
+  }
+  
+  if (saveToFile)
+  {
+    std::vector<std::string> splitFileName = SplitFileName(clusterFileName);
+    std::string fname = splitFileName[1]+"_nbClusters_"+std::to_string(_totalNumberOfClusters)+splitFileName[2];
+    FILE* f = fopen(fname.c_str(),"w");
+    printf("start writing cluster data to file: %s\n",fname.c_str());
+    std::string str ="ELE_NUM GP_NUM CLUSTER_ID MAT_NUM";
+    fprintf(f,"%s\n",str.c_str());
+    // write data
+    for (std::map<int, int>::const_iterator it = _clusterIndex.begin();  it!= _clusterIndex.end(); it++)
+    { 
+      int type = it->first;
+      int clusterId = it->second;
+      int ele, gp;
+      numericalMaterialBase::getTwoIntsFromType(type,ele,gp);
+      fprintf(f, "%d %d %d %d\n",ele,gp,clusterId,_materialMap[type]);
+    }
+    fclose(f);
+    printf("done writing cluster data to file: %s \n",fname.c_str());
+  }
+  
+  if (saveSummaryFile)
+  {
+    std::vector<std::string> splitFileName = SplitFileName(summaryFileName);
+    std::string fname = splitFileName[1]+"_nbClusters_"+std::to_string(_totalNumberOfClusters)+splitFileName[2];
+    FILE* f = fopen(fname.c_str(),"w");
+    printf("start writing cluster summary to file: %s\n",fname.c_str());
+    std::string str ="CLUSTER_ID MAT_NUM VOLUME STD";
+    // write data
+    for (int i=0; i<6; i++)
+    {
+      for (int j=0; j<6; j++)
+      {
+        str += " A_"+std::to_string(i)+std::to_string(j);
+      }
+    }
+    fprintf(f,"%s\n",str.c_str());
+    // write data
+    for (std::map<int, STensor43>::const_iterator it = _clusterAverageStrainConcentrationTensorMap.begin();  it!= _clusterAverageStrainConcentrationTensorMap.end(); it++)
+    { 
+      int cluster = it->first;
+      int matNum = _clusterMaterialMap[cluster];
+      double volume = _clusterVolumeMap[cluster];
+      double stdev = 0.;
+      fprintf(f,"%d %d %.16g %.16g",cluster, matNum, volume, stdev);
+      const STensor43& Amean = it->second;
+      static fullMatrix<double> mat(6,6);
+      STensorOperation::fromSTensor43ToFullMatrix(Amean,mat);
+      //printf("writing data for cluster %d\n",cluster);
+      //mat.print("average strain concentration tensor");
+      for (int i=0; i<6; i++)
+      {
+        for (int j=0; j<6; j++)
+        {
+          fprintf(f, " %.16g",mat(i,j));
+        }
+      }
+      fprintf(f, "\n");
+    }
+    fclose(f);
+    printf("done writing cluster summary to file: %s \n",fname.c_str());
+  }
+
+  return true;
+};
+
 
 bool Clustering::computeClustersPlasticityHierarchical(int minNumberElementEachCluster, int maxNbIterations,
                                                        bool saveToFile, std::string clusterFileName,
diff --git a/NonLinearSolver/modelReduction/Clustering.h b/NonLinearSolver/modelReduction/Clustering.h
index 114b25853f995f2a996a915c6a102e12207593ca..0288a3466037e3153372782f29010be22bbeab4a 100644
--- a/NonLinearSolver/modelReduction/Clustering.h
+++ b/NonLinearSolver/modelReduction/Clustering.h
@@ -272,6 +272,12 @@ class Clustering
                                  bool saveClusterMeanFiberOrientationFile=true, const std::string saveClusterMeanFiberOrientationFileName = "ClusterMeanFiberOrientation.csv"
                                  );
                                  
+    bool computeClustersLocalOrientation(int minNumberElementEachCluster, int maxNbIterations=10000,
+                                         bool saveToFile=true, const std::string clusterFileName = "ClusterData.csv", 
+                                         bool saveSummaryFile=true, const std::string summaryFileName = "ClusterDataSummary.csv",
+                                         bool saveClusterMeanFiberOrientationFile=true, const std::string saveClusterMeanFiberOrientationFileName = "ClusterMeanFiberOrientation.csv"
+                                        );
+                                 
     bool computeClustersPlasticityHierarchical(int minNumberElementEachCluster, int maxNbIterations=10000,
                                                       bool saveToFile=true, const std::string clusterFileName = "ClusterDataHierarchical.csv", 
                                                       bool saveSummaryFile=true, const std::string summaryFileName = "ClusterDataSummaryHierarchical.csv"