soft constraints compiles but buggy.

include/directional/polyvector_field.h

@@ -192,12 +192,12 @@ namespace directional
       }
       singleReducRhs(realN-currFaceNumDof-1) = -constVectorComplex;
       
-      cout<<"localFaceReducMats[pvData.constFaces(i)]: "<< localFaceReducMats[pvData.constFaces(i)]<<endl;
+      /*cout<<"localFaceReducMats[pvData.constFaces(i)]: "<< localFaceReducMats[pvData.constFaces(i)]<<endl;
       cout<<"localFaceReducRhs[pvData.constFaces(i)]: "<< localFaceReducRhs[pvData.constFaces(i)]<<endl;
       
       cout<<"singleReducMat: "<<singleReducMat<<endl;
       cout<<"singleReducRhs: "<<singleReducRhs<<endl;
-      cout<<"currFaceNumDof: "<<currFaceNumDof<<endl;
+      cout<<"currFaceNumDof: "<<currFaceNumDof<<endl;*/
       
       localFaceReducRhs[pvData.constFaces(i)] = localFaceReducMats[pvData.constFaces(i)]*singleReducRhs + localFaceReducRhs[pvData.constFaces(i)];
       localFaceReducMats[pvData.constFaces(i)] = localFaceReducMats[pvData.constFaces(i)]*singleReducMat;
@@ -269,8 +269,9 @@ namespace directional
       MatrixXcd invSingleReducMat = singleReducMat.completeOrthogonalDecomposition().pseudoInverse();
       MatrixXcd IAiA = MatrixXcd::Identity(realN,realN) - singleReducMat*invSingleReducMat;
       singleReducRhs = IAiA*singleReducRhs;
-      for (int j=0;j<pvData.N;j++)
-        for (int k=0;k<pvData.N;k+=jump)
+      //This assumes the sign symmetry is dealt by a reduction matrix, so acts on realN
+      for (int j=0;j<IAiA.rows();j++)
+        for (int k=0;k<IAiA.cols();k++)
           alignTriplets.push_back(Triplet<complex<double>>(rowCounter+j, k*F.rows()+i, IAiA(j,k)));
       
       alignRhsList.push_back(singleReducRhs);
@@ -470,6 +471,8 @@ namespace directional
                                    const Eigen::MatrixXi& F,
                                    const Eigen::VectorXi& constFaces,
                                    const Eigen::MatrixXd& constVectors,
+                                   const double smoothWeight,
+                                   const double roSyWeight,
                                    const Eigen::VectorXd& alignWeights,
                                    const int N,
                                    Eigen::MatrixXcd& polyVectorField)
@@ -490,6 +493,8 @@ namespace directional
     pvData.constFaces=constFaces;
     pvData.constVectors=constVectors;
     pvData.wAlignment = alignWeights;
+    pvData.wSmooth = smoothWeight;
+    pvData.wRoSy = roSyWeight;
     igl::edge_topology(V, F, EV, xi, EF);
     polyvector_precompute(V,F,EV,EF, B1,B2, N, pvData);
     polyvector_field(pvData, polyVectorField);

tutorial/303_FlexibleConstraints/main.cpp

@@ -17,12 +17,13 @@ Eigen::MatrixXi F, EV, EF, FE;
 Eigen::MatrixXd V;
 Eigen::MatrixXd normals,constVectors;
 Eigen::MatrixXd rawFieldConstraints, rawFieldHard, rawFieldSoft;
-Eigen::VectorXd w;
+Eigen::VectorXd alignWeights;
 Eigen::MatrixXcd pvFieldHard, pvFieldSoft;
+double smoothWeight, roSyWeight;
 
 directional::DirectionalViewer viewer;
 
-int N = 4;
+int N = 3;
 
 typedef enum {CONSTRAINTS, HARD_PRESCRIPTION, SOFT_PRESCRIPTION} ViewingModes;
 ViewingModes viewingMode=CONSTRAINTS;
@@ -31,9 +32,9 @@ bool normalized=false;
 
 void recompute_field()
 {
-  directional::polyvector_field(V, F, constFaces, constVectors, Eigen::VectorXd::Constant(constFaces.size(),-1), N, pvFieldHard);
+  directional::polyvector_field(V, F, constFaces, constVectors, smoothWeight, roSyWeight, Eigen::VectorXd::Constant(constFaces.size(),-1), N, pvFieldHard);
   //std::cout<<"Done computing field!"<<std::endl;
-  //directional::polyvector_field(V, F, b, bc, w, N, pvFieldSoft);
+  directional::polyvector_field(V, F, constFaces, constVectors, smoothWeight, roSyWeight, alignWeights, N, pvFieldSoft);
 }
 
 void update_visualization()
@@ -65,7 +66,7 @@ void update_visualization()
   
   if (viewingMode==SOFT_PRESCRIPTION){
     viewer.set_selected_faces(Eigen::VectorXi());
-    directional::polyvector_to_raw(V, F, pvFieldSoft, N, rawFieldSoft, true);
+    directional::polyvector_to_raw(V, F, pvFieldSoft, N, rawFieldSoft, N%2==0);
     if (normalized)
       for(int n = 0; n < N; n++)
         rawFieldSoft.middleCols(n*3, 3).rowwise().normalize();
@@ -85,23 +86,39 @@ bool key_down(igl::opengl::glfw::Viewer& viewer, int key, int modifiers)
   {
       // Select vector
     case '1':
-      w.array()*=2.0;
+      roSyWeight+=0.1;
+      std::cout<<"roSyWeight: "<<roSyWeight<<std::endl;
       recompute_field();
       break;
     case '2':
-      w.array()/=2.0;
+      roSyWeight-=0.1;
+      roSyWeight = (roSyWeight < 0.0 ? 0.0 : roSyWeight);
+      std::cout<<"roSyWeight: "<<roSyWeight<<std::endl;
       recompute_field();
       break;
       
     case '3':
-      viewingMode=CONSTRAINTS;
+      alignWeights.array()+=0.1;
+      std::cout<<"alignWeights: "<<alignWeights(1)<<std::endl;
+      recompute_field();
       break;
       
     case '4':
-      viewingMode=SOFT_PRESCRIPTION;
+      alignWeights.array()-=0.1;
+      alignWeights = (alignWeights(1) < 0.0 ? Eigen::VectorXd::Zero(alignWeights.size()) : alignWeights);
+      std::cout<<"alignWeights: "<<alignWeights(1)<<std::endl;
+      recompute_field();
       break;
       
     case '5':
+      viewingMode=CONSTRAINTS;
+      break;
+      
+    case '6':
+      viewingMode=SOFT_PRESCRIPTION;
+      break;
+      
+    case '7':
       viewingMode=HARD_PRESCRIPTION;
       break;
     
@@ -178,6 +195,10 @@ int main()
     posInFace(constFaces(i))++;
   }
   
+  smoothWeight = 1.0;
+  roSyWeight = 0.0;
+  alignWeights = Eigen::VectorXd::Constant(constFaces.size(),10000.0);
+  
   //triangle mesh setup
   viewer.set_mesh(V, F);
   recompute_field();