Changed "most aniso" mesh metric intersection

Geo/STensor3.cpp

@@ -64,22 +64,76 @@ SMetric3 intersection_conserveM1 (const SMetric3 &m1, const SMetric3 &m2)
   return iv;
 }
 
+namespace {
+
+double anisoRatio2D (double v0x, double v0y, double v0z,
+                     double v1x, double v1y, double v1z,
+                     double v2x, double v2y, double v2z,
+                     double s0, double s1, double s2) {
+
+  static const double eps = 1.e-8;
+
+  double ratio;
+  if (v0x*v0x+v0y*v0y+(v0z-1.)*(v0z-1.) < eps)        // If 1st eigenvect corresponds to z dir. ...
+    ratio = std::min(fabs(s1/s2),fabs(s2/s1));        // ... consider ratio of 2nd and 3rd eigenval.
+  else if (v1x*v1x+v1y*v1y+(v1z-1.)*(v1z-1.) < eps)   // If 2nd eigenvect corresponds to z dir. ...
+    ratio = std::min(fabs(s0/s2),fabs(s2/s0));        // ... consider ratio of 1st and 3rd eigenval.
+  else if (v2x*v2x+v2y*v2y+(v2z-1.)*(v2z-1.) < eps)   // If 3rd eigenvect corresponds to z dir. ...
+    ratio = std::min(fabs(s0/s1),fabs(s1/s0));        // ... consider ratio of 1st and 2nd eigenval.
+  else {
+    printf("Error in anisoRatio2D: z direction not found!\n");
+    ratio = 0.;
+  }
+
+  return ratio;
+
+}
+
+}
+
 // preserve orientation of the most anisotropic metric !!!
 SMetric3 intersection_conserve_mostaniso (const SMetric3 &m1, const SMetric3 &m2)
 {
   fullMatrix<double> V1(3,3);
   fullVector<double> S1(3);
   m1.eig(V1,S1,true);
-  double lambda1_min = std::min(std::min(fabs(S1(0)),fabs(S1(1))),fabs(S1(2)));
+  double ratio1 = fabs(S1(0)/S1(2));  // Minimum ratio because we take sorted eigenvalues
   fullMatrix<double> V2(3,3);
   fullVector<double> S2(3);
   m2.eig(V2,S2,true);
-  double lambda2_min = std::min(std::min(fabs(S2(0)),fabs(S2(1))),fabs(S2(2)));
+  double ratio2 = fabs(S2(0)/S2(2));  // Minimum ratio because we take sorted eigenvalues
+
+  if (ratio1 < ratio2)
+    return intersection_conserveM1(m1, m2);
+  else
+    return intersection_conserveM1(m2, m1);
+}
+
+// preserve orientation of the most anisotropic metric in 2D!!!
+SMetric3 intersection_conserve_mostaniso_2d (const SMetric3 &m1, const SMetric3 &m2)
+{
+
+  fullMatrix<double> V1(3,3);
+  fullVector<double> S1(3);
+  m1.eig(V1,S1,false);
+  double ratio1 = anisoRatio2D(V1(0,0),V1(1,0),V1(2,0),
+                               V1(0,1),V1(1,1),V1(2,1),
+                               V1(0,2),V1(1,2),V1(2,2),
+                               S1(0),S1(1),S1(2));
+
+  fullMatrix<double> V2(3,3);
+  fullVector<double> S2(3);
+  m2.eig(V2,S2,false);
+  double ratio2 = anisoRatio2D(V2(0,0),V2(1,0),V2(2,0),
+                               V2(0,1),V2(1,1),V2(2,1),
+                               V2(0,2),V2(1,2),V2(2,2),
+                               S2(0),S2(1),S2(2));
 
-  if (lambda1_min < lambda2_min)
+  if (ratio1 < ratio2)
     return intersection_conserveM1(m1, m2);
   else
     return intersection_conserveM1(m2, m1);
+
 }
 
 // (1-t) * m1 + t * m2

Geo/STensor3.h

@@ -169,6 +169,7 @@ SMetric3 intersection_conserveM1 (const SMetric3 &m1,
 
 // preserve orientation of the most anisotropic metric
 SMetric3 intersection_conserve_mostaniso (const SMetric3 &m1, const SMetric3 &m2);
+SMetric3 intersection_conserve_mostaniso_2d (const SMetric3 &m1, const SMetric3 &m2);
 // compute the largest inscribed ellipsoid...
 SMetric3 intersection (const SMetric3 &m1,
                        const SMetric3 &m2);

Mesh/meshMetric.cpp

@@ -95,7 +95,9 @@ void meshMetric::updateMetrics()
 
     if (setOfMetrics.size() > 1)
       for (unsigned int i=1;i<setOfMetrics.size();i++){
-        _nodalMetrics[ver] = intersection_conserve_mostaniso(_nodalMetrics[ver],setOfMetrics[i][ver]);
+        _nodalMetrics[ver] = (_dim == 3) ?
+          intersection_conserve_mostaniso(_nodalMetrics[ver],setOfMetrics[i][ver]) :
+          intersection_conserve_mostaniso_2d(_nodalMetrics[ver],setOfMetrics[i][ver]);
         _nodalSizes[ver] = std::min(_nodalSizes[ver],setOfSizes[i][ver]);
       }
   }
@@ -105,19 +107,22 @@ void meshMetric::updateMetrics()
 void meshMetric::exportInfo(const char * fileendname)
 {
   if (needMetricUpdate) updateMetrics();
-  std::stringstream sg,sm,sl,sh;
+  std::stringstream sg,sm,sl,sh,shm;
   sg << "meshmetric_gradients_" << fileendname;
   sm << "meshmetric_metric_" << fileendname;
   sl << "meshmetric_levelset_" << fileendname;
   sh << "meshmetric_hessian_" << fileendname;
+  shm << "meshmetric_hessianmat_" << fileendname;
   std::ofstream out_grad(sg.str().c_str());
   std::ofstream out_metric(sm.str().c_str());
   std::ofstream out_ls(sl.str().c_str());
   std::ofstream out_hess(sh.str().c_str());
+  std::ofstream out_hessmat(shm.str().c_str());
   out_grad << "View \"ls_gradient\"{" << std::endl;
   out_metric << "View \"metric\"{" << std::endl;
   out_ls << "View \"ls\"{" << std::endl;
   out_hess << "View \"hessian\"{" << std::endl;
+  out_hessmat << "View \"hessian_mat\"{" << std::endl;
   std::vector<MElement*>::iterator itelem = _elements.begin();
   std::vector<MElement*>::iterator itelemen = _elements.end();
   for (;itelem!=itelemen;itelem++){
@@ -127,12 +132,14 @@ void meshMetric::exportInfo(const char * fileendname)
       out_grad << "VT(";
       out_ls << "ST(";
       out_hess << "ST(";
+      out_hessmat << "TT(";
     }
     else {
       out_metric << "TS(";
       out_grad << "VS(";
       out_ls << "SS(";
       out_hess << "SS(";
+      out_hessmat << "TS(";
     }
     for ( int i = 0; i < e->getNumVertices(); i++) {
       MVertex *ver = e->getVertex(i);
@@ -140,17 +147,20 @@ void meshMetric::exportInfo(const char * fileendname)
       out_grad << ver->x() << "," << ver->y() << "," << ver->z();
       out_ls << ver->x() << "," << ver->y() << "," << ver->z();
       out_hess << ver->x() << "," << ver->y() << "," << ver->z();
+      out_hessmat << ver->x() << "," << ver->y() << "," << ver->z();
       if (i!=e->getNumVertices()-1){
         out_metric << ",";
         out_grad << ",";
         out_ls << ",";
         out_hess << ",";
+        out_hessmat << ",";
       }
       else{
         out_metric << "){";
         out_grad << "){";
         out_ls << "){";
         out_hess << "){";
+        out_hessmat << "){";
       }
     }
     for ( int i = 0; i < e->getNumVertices(); i++) {
@@ -174,8 +184,15 @@ void meshMetric::exportInfo(const char * fileendname)
         else out_grad << ",";
         for (int l=0;l<3;l++){
           out_metric << _nodalMetrics[ver](k,l);
-          if ((k == 2) && (l == 2) && (i == (e->getNumVertices() - 1))) out_metric << "};" << std::endl;
-          else out_metric << ",";
+          out_hessmat << dgrads[k][ver](l);
+          if ((k == 2) && (l == 2) && (i == (e->getNumVertices() - 1))) {
+            out_metric << "};" << std::endl;
+            out_hessmat << "};" << std::endl;
+          }
+          else {
+            out_metric << ",";
+            out_hessmat << ",";
+          }
         }
       }
     }
@@ -184,10 +201,12 @@ void meshMetric::exportInfo(const char * fileendname)
   out_metric << "};" << std::endl;
   out_ls << "};" << std::endl;
   out_hess << "};" << std::endl;
+  out_hessmat << "};" << std::endl;
   out_grad.close();
   out_metric.close();
   out_ls.close();
   out_hess.close();
+  out_hessmat.close();
 }
 
 meshMetric::~meshMetric()
@@ -360,10 +379,9 @@ void meshMetric::computeMetricLevelSet(MVertex *ver, SMetric3 &hessian,  SMetric
 
 void meshMetric::computeMetricHessian(MVertex *ver, SMetric3 &hessian,  SMetric3 &metric, double &size, double x, double y, double z)
 {
-  double signed_dist;
+
   SVector3 gradudx, gradudy,gradudz,gr;
   if(ver != NULL){
-    signed_dist = vals[ver];
     gr = grads[ver];
     gradudx = dgrads[0][ver];
     gradudy = dgrads[1][ver];
@@ -373,7 +391,6 @@ void meshMetric::computeMetricHessian(MVertex *ver, SMetric3 &hessian,  SMetric3
     hessian(2,0) = gradudz(0); hessian(2,1) = gradudz(1); hessian(2,2) = gradudz(2);
   }
   else if (ver == NULL){
-    signed_dist = (*_fct)(x,y,z); 
     _fct->gradient(x,y,z,gr(0),gr(1),gr(2));
     _fct->hessian(x,y,z, hessian(0,0),hessian(0,1),hessian(0,2),
 		  hessian(1,0),hessian(1,1),hessian(1,2),
@@ -394,7 +411,7 @@ void meshMetric::computeMetricHessian(MVertex *ver, SMetric3 &hessian,  SMetric3
   
   SVector3 t1 (V(0,0),V(1,0),V(2,0));
   SVector3 t2 (V(0,1),V(1,1),V(2,1));
-  SVector3 t3 (V(0,2),V(1,2),V(2,2));
+  SVector3 t3  = (_dim == 3) ? SVector3(V(0,2),V(1,2),V(2,2)) : SVector3(0.,0.,1.);
   
   size =  std::min(std::min(1/sqrt(lambda1),1/sqrt(lambda2)),1/sqrt(lambda3));
   metric = SMetric3(lambda1,lambda2,lambda3,t1,t2,t3);