|
|
@@ -0,0 +1,1179 @@
|
|
|
+// This file is part of libigl, a simple c++ geometry processing library.
|
|
|
+//
|
|
|
+// Copyright (C) 2015 Olga Diamanti <olga.diam@gmail.com>
|
|
|
+//
|
|
|
+// This Source Code Form is subject to the terms of the Mozilla Public License
|
|
|
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
|
|
|
+// obtain one at http://mozilla.org/MPL/2.0/.
|
|
|
+
|
|
|
+#include <igl/integrable_polyvector_fields.h>
|
|
|
+#include <igl/field_local_global_conversions.h>
|
|
|
+#include <igl/parallel_transport_angles.h>
|
|
|
+#include <igl/local_basis.h>
|
|
|
+#include <igl/edge_topology.h>
|
|
|
+#include <igl/sparse.h>
|
|
|
+#include <igl/sort.h>
|
|
|
+#include <igl/slice.h>
|
|
|
+#include <igl/slice_into.h>
|
|
|
+#include <igl/sort_vectors_ccw.h>
|
|
|
+#include <iostream>
|
|
|
+
|
|
|
+IGL_INLINE igl::integrable_polyvector_fields_parameters::integrable_polyvector_fields_parameters():
|
|
|
+numIter(5),
|
|
|
+wBarrier(0.1),
|
|
|
+sBarrier(0.9),
|
|
|
+wCurl(10),
|
|
|
+wQuotCurl(10),
|
|
|
+wSmooth(1.),
|
|
|
+wCloseUnconstrained(1e-3),
|
|
|
+wCloseConstrained(100),
|
|
|
+redFactor_wsmooth(.8),
|
|
|
+gamma(0.1),
|
|
|
+tikh_gamma(1e-8),
|
|
|
+do_partial(false)
|
|
|
+{}
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+namespace igl {
|
|
|
+ template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+ class IntegrableFieldSolver
|
|
|
+ {
|
|
|
+ private:
|
|
|
+
|
|
|
+ IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC> &data;
|
|
|
+ //Symbolic calculations
|
|
|
+ IGL_INLINE void rj_barrier_face(const Eigen::RowVectorXd &vec2D_a,
|
|
|
+ const double &s,
|
|
|
+ Eigen::VectorXd &residuals,
|
|
|
+ bool do_jac = false,
|
|
|
+ Eigen::MatrixXd &J = *(Eigen::MatrixXd*)NULL);
|
|
|
+ IGL_INLINE void rj_polycurl_edge(const Eigen::RowVectorXd &vec2D_a,
|
|
|
+ const Eigen::RowVector2d &ea,
|
|
|
+ const Eigen::RowVectorXd &vec2D_b,
|
|
|
+ const Eigen::RowVector2d &eb,
|
|
|
+ Eigen::VectorXd &residuals,
|
|
|
+ bool do_jac = false,
|
|
|
+ Eigen::MatrixXd &Jac = *(Eigen::MatrixXd*)NULL);
|
|
|
+ IGL_INLINE void rj_quotcurl_edge_polyversion(const Eigen::RowVectorXd &vec2D_a,
|
|
|
+ const Eigen::RowVector2d &ea,
|
|
|
+ const Eigen::RowVectorXd &vec2D_b,
|
|
|
+ const Eigen::RowVector2d &eb,
|
|
|
+ Eigen::VectorXd &residuals,
|
|
|
+ bool do_jac = false,
|
|
|
+ Eigen::MatrixXd &Jac = *(Eigen::MatrixXd*)NULL);
|
|
|
+ IGL_INLINE void rj_smoothness_edge(const Eigen::RowVectorXd &vec2D_a,
|
|
|
+ const Eigen::RowVectorXd &vec2D_b,
|
|
|
+ const double &k,
|
|
|
+ const int nA,
|
|
|
+ const int nB,
|
|
|
+ Eigen::VectorXd &residuals,
|
|
|
+ bool do_jac = false,
|
|
|
+ Eigen::MatrixXd &Jac = *(Eigen::MatrixXd*)NULL);
|
|
|
+
|
|
|
+ public:
|
|
|
+ IGL_INLINE IntegrableFieldSolver(IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC> &cffsoldata);
|
|
|
+
|
|
|
+ IGL_INLINE bool solve(integrable_polyvector_fields_parameters ¶ms,
|
|
|
+ Eigen::PlainObjectBase<DerivedFF>& current_field,
|
|
|
+ bool current_field_is_not_ccw);
|
|
|
+
|
|
|
+ IGL_INLINE void solveGaussNewton(integrable_polyvector_fields_parameters ¶ms,
|
|
|
+ const Eigen::VectorXd &x_initial,
|
|
|
+ Eigen::VectorXd &x);
|
|
|
+
|
|
|
+ //Compute residuals and Jacobian for Gauss Newton
|
|
|
+ IGL_INLINE double RJ(const Eigen::VectorXd &x,
|
|
|
+ const Eigen::VectorXd &x0,
|
|
|
+ const integrable_polyvector_fields_parameters ¶ms,
|
|
|
+ bool doJacs = false);
|
|
|
+
|
|
|
+ IGL_INLINE void RJ_Smoothness(const Eigen::MatrixXd &sol2D,
|
|
|
+ const double &wSmoothSqrt,
|
|
|
+ const int startRowInJacobian,
|
|
|
+ bool doJacs = false,
|
|
|
+ const int startIndexInVectors = 0);
|
|
|
+ IGL_INLINE void RJ_Barrier(const Eigen::MatrixXd &sol2D,
|
|
|
+ const double &s,
|
|
|
+ const double &wBarrierSqrt,
|
|
|
+ const int startRowInJacobian,
|
|
|
+ bool doJacs = false,
|
|
|
+ const int startIndexInVectors = 0);
|
|
|
+ IGL_INLINE void RJ_Closeness(const Eigen::MatrixXd &sol2D,
|
|
|
+ const Eigen::MatrixXd &sol02D,
|
|
|
+ const double &wCloseUnconstrainedSqrt,
|
|
|
+ const double &wCloseConstrainedSqrt,
|
|
|
+ const bool do_partial,
|
|
|
+ const int startRowInJacobian,
|
|
|
+ bool doJacs = false,
|
|
|
+ const int startIndexInVectors = 0);
|
|
|
+ IGL_INLINE void RJ_Curl(const Eigen::MatrixXd &sol2D,
|
|
|
+ const double &wCASqrt,
|
|
|
+ const double &wCBSqrt,
|
|
|
+ const int startRowInJacobian,
|
|
|
+ bool doJacs = false,
|
|
|
+ const int startIndexInVectors = 0);
|
|
|
+ IGL_INLINE void RJ_QuotCurl(const Eigen::MatrixXd &sol2D,
|
|
|
+ const double &wQuotCurlSqrt,
|
|
|
+ const int startRowInJacobian,
|
|
|
+ bool doJacs = false,
|
|
|
+ const int startIndexInVectors = 0);
|
|
|
+
|
|
|
+
|
|
|
+ };
|
|
|
+};
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::IntegrableFieldSolverData()
|
|
|
+{}
|
|
|
+
|
|
|
+
|
|
|
+template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+precomputeMesh(const Eigen::PlainObjectBase<DerivedV> &_V,
|
|
|
+ const Eigen::PlainObjectBase<DerivedF> &_F)
|
|
|
+{
|
|
|
+ numV = _V.rows();
|
|
|
+ numF = _F.rows();
|
|
|
+ numVariables = 2*2*numF;
|
|
|
+ //Mesh stuff
|
|
|
+ igl::edge_topology(_V,_F,E,F2E,E2F);
|
|
|
+ numE = E.rows();
|
|
|
+ igl::local_basis(_V,_F,B1,B2,FN);
|
|
|
+ computeInteriorEdges();
|
|
|
+ igl::parallel_transport_angles(_V, _F, FN, E2F, F2E, K);
|
|
|
+ EVecNorm.setZero(numE,3);
|
|
|
+ for (int k = 0; k<numE; ++k)
|
|
|
+ EVecNorm.row(k) = (_V.row(E(k,1))-_V.row(E(k,0))).normalized();
|
|
|
+}
|
|
|
+
|
|
|
+template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+initializeConstraints(const Eigen::VectorXi& b,
|
|
|
+ const Eigen::PlainObjectBase<DerivedC>& bc,
|
|
|
+ const Eigen::VectorXi& constraint_level)
|
|
|
+{
|
|
|
+ //save constrained
|
|
|
+ Eigen::VectorXi iSorted, constrained_unsorted;
|
|
|
+ constrained_unsorted.resize(2*2*b.size());
|
|
|
+ is_constrained_face.setZero(numF, 1);
|
|
|
+ int ind = 0;
|
|
|
+ indInConstrained.setConstant(numF,1,-1);
|
|
|
+ for (int k =0; k<b.size(); ++k)
|
|
|
+ {
|
|
|
+ is_constrained_face[b[k]] = constraint_level[k];
|
|
|
+ for (int i=0; i<2;++i)
|
|
|
+ {
|
|
|
+ int xi = 2*2*b[k] + 2*i +0;
|
|
|
+ int yi = 2*2*b[k] + 2*i +1;
|
|
|
+ constrained_unsorted[ind++] = xi;
|
|
|
+ constrained_unsorted[ind++] = yi;
|
|
|
+ }
|
|
|
+ indInConstrained[b[k]] = k;
|
|
|
+ }
|
|
|
+ //sort in descending order (so removing rows will work)
|
|
|
+ igl::sort(constrained_unsorted, 1, false, constrained, iSorted);
|
|
|
+ constrained_vec3 = bc.template cast<double>();
|
|
|
+
|
|
|
+}
|
|
|
+
|
|
|
+template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+makeFieldCCW(Eigen::MatrixXd &sol3D)
|
|
|
+{
|
|
|
+ //sort ccw
|
|
|
+ Eigen::RowVectorXd t;
|
|
|
+ Eigen::RowVectorXd all(1,2*2*3);
|
|
|
+ Eigen::VectorXi order;
|
|
|
+ for (int fi=0; fi<numF; ++fi)
|
|
|
+ {
|
|
|
+ //take all 4 vectors (including opposites) and pick two that are in ccw order
|
|
|
+ all << sol3D.row(fi), -sol3D.row(fi);
|
|
|
+ igl::sort_vectors_ccw(all, FN.row(fi).eval(), order, true, t);
|
|
|
+ //if we are in a constrained face, we need to make sure that the first vector is always the same vector as in the constraints
|
|
|
+ if(is_constrained_face[fi])
|
|
|
+ {
|
|
|
+ const Eigen::RowVector3d &constraint = constrained_vec3.block(indInConstrained[fi], 0, 1, 3);;
|
|
|
+ int best_i = -1; double best_score = std::numeric_limits<double>::max();
|
|
|
+ for (int i = 0; i<2*2; ++i)
|
|
|
+ {
|
|
|
+ double score = (t.segment(i*3,3) - constraint).norm();
|
|
|
+ if (score<best_score)
|
|
|
+ {
|
|
|
+ best_score = score;
|
|
|
+ best_i = i;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ //do a circshift
|
|
|
+ Eigen::RowVectorXd temp = t.segment(0, 3*best_i);
|
|
|
+ int s1 = 3*best_i;
|
|
|
+ int n2 = 3*best_i;
|
|
|
+ int n1 = 3*2*2-n2;
|
|
|
+ t.segment(0,n1) = t.segment(s1,n1);
|
|
|
+ t.segment(n1, n2) = temp;
|
|
|
+
|
|
|
+ }
|
|
|
+ sol3D.row(fi) = t.segment(0, 2*3);
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+initializeOriginalVariable(const Eigen::PlainObjectBase<DerivedFF>& original_field)
|
|
|
+{
|
|
|
+ Eigen::MatrixXd sol2D;
|
|
|
+ Eigen::MatrixXd sol3D = original_field.template cast<double>();
|
|
|
+ makeFieldCCW(sol3D);
|
|
|
+ igl::global2local(B1, B2, sol3D, sol2D);
|
|
|
+ xOriginal.setZero(numVariables);
|
|
|
+ for (int i =0; i<numF; i++)
|
|
|
+ xOriginal.segment(i*2*2, 2*2) = sol2D.row(i);
|
|
|
+
|
|
|
+}
|
|
|
+
|
|
|
+template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+computeInteriorEdges()
|
|
|
+{
|
|
|
+ Eigen::VectorXi isBorderEdge;
|
|
|
+ // Flag border edges
|
|
|
+ numInteriorEdges = 0;
|
|
|
+ isBorderEdge.setZero(numE,1);
|
|
|
+ indFullToInterior = -1.*Eigen::VectorXi::Ones(numE,1);
|
|
|
+
|
|
|
+ for(unsigned i=0; i<numE; ++i)
|
|
|
+ {
|
|
|
+ if ((E2F(i,0) == -1) || ((E2F(i,1) == -1)))
|
|
|
+ isBorderEdge[i] = 1;
|
|
|
+ else
|
|
|
+ {
|
|
|
+ indFullToInterior[i] = numInteriorEdges;
|
|
|
+ numInteriorEdges++;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ E2F_int.resize(numInteriorEdges, 2);
|
|
|
+ indInteriorToFull.setZero(numInteriorEdges,1);
|
|
|
+ int ii = 0;
|
|
|
+ for (int k=0; k<numE; ++k)
|
|
|
+ {
|
|
|
+ if (isBorderEdge[k])
|
|
|
+ continue;
|
|
|
+ E2F_int.row(ii) = E2F.row(k);
|
|
|
+ indInteriorToFull[ii] = k;
|
|
|
+ ii++;
|
|
|
+ }
|
|
|
+
|
|
|
+}
|
|
|
+
|
|
|
+template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+add_Jacobian_to_svector(const int &toplace,
|
|
|
+ const Eigen::MatrixXd &tJac,
|
|
|
+ Eigen::VectorXd &SS_Jac)
|
|
|
+{
|
|
|
+ int numInnerRows = tJac.rows();
|
|
|
+ int numInnerCols = tJac.cols();
|
|
|
+ int ind = toplace;
|
|
|
+ for (int j=0; j<numInnerRows; ++j)
|
|
|
+ for (int k=0; k<numInnerCols; ++k, ++ind)
|
|
|
+ SS_Jac(ind) = tJac(j,k);
|
|
|
+}
|
|
|
+
|
|
|
+template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+add_jac_indices_face(const int numInnerRows,
|
|
|
+ const int numInnerCols,
|
|
|
+ const int startRowInJacobian,
|
|
|
+ const int startIndexInVectors,
|
|
|
+ Eigen::VectorXi &Rows,
|
|
|
+ Eigen::VectorXi &Columns)
|
|
|
+{
|
|
|
+ for (int fi=0; fi<numF; ++fi)
|
|
|
+ {
|
|
|
+
|
|
|
+ int startRow = startRowInJacobian+numInnerRows*fi;
|
|
|
+ int startIndex = startIndexInVectors+numInnerRows*numInnerCols*fi;
|
|
|
+ face_Jacobian_indices(startRow, startIndex, fi, 2, numInnerRows, numInnerCols, Rows, Columns);
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+face_Jacobian_indices(const int &startRow,
|
|
|
+ const int &toplace,
|
|
|
+ const int& fi,
|
|
|
+ const int& half_degree,
|
|
|
+ const int &numInnerRows,
|
|
|
+ const int &numInnerCols,
|
|
|
+ Eigen::VectorXi &rows,
|
|
|
+ Eigen::VectorXi &columns)
|
|
|
+{
|
|
|
+ int ind = toplace;
|
|
|
+ for (int j=0; j<numInnerRows; ++j)
|
|
|
+ {
|
|
|
+ for (int k=0; k<numInnerCols; ++k, ++ind)
|
|
|
+ {
|
|
|
+ int iv = k/2;//which vector (0..half_degree-1) am i at
|
|
|
+ int ixy = k%2;//which part (real/imag) am i at
|
|
|
+ rows(ind) = startRow+j;
|
|
|
+ columns(ind) = 2*half_degree*fi + 2*iv +ixy;
|
|
|
+ }
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+add_jac_indices_edge(const int numInnerRows,
|
|
|
+ const int numInnerCols,
|
|
|
+ const int startRowInJacobian,
|
|
|
+ const int startIndexInVectors,
|
|
|
+ Eigen::VectorXi &Rows,
|
|
|
+ Eigen::VectorXi &Columns)
|
|
|
+{
|
|
|
+ for (int ii=0; ii<numInteriorEdges; ++ii)
|
|
|
+ {
|
|
|
+ // the two faces of the flap
|
|
|
+ int a = E2F_int(ii,0);
|
|
|
+ int b = E2F_int(ii,1);
|
|
|
+
|
|
|
+
|
|
|
+ int startRow = startRowInJacobian+numInnerRows*ii;
|
|
|
+ int startIndex = startIndexInVectors+numInnerRows*numInnerCols*ii;
|
|
|
+
|
|
|
+ edge_Jacobian_indices(startRow, startIndex, a, b, 2, numInnerRows, numInnerCols, Rows, Columns);
|
|
|
+
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+edge_Jacobian_indices(const int &startRow,
|
|
|
+ const int &toplace,
|
|
|
+ const int& a,
|
|
|
+ const int& b,
|
|
|
+ const int& half_degree,
|
|
|
+ const int &numInnerRows,
|
|
|
+ const int &numInnerCols,
|
|
|
+ Eigen::VectorXi &rows,
|
|
|
+ Eigen::VectorXi &columns)
|
|
|
+{
|
|
|
+ int ind = toplace;
|
|
|
+ for (int j=0; j<numInnerRows; ++j)
|
|
|
+ {
|
|
|
+ for (int k=0; k<numInnerCols; ++k, ++ind)
|
|
|
+ {
|
|
|
+ int f = (k<2*half_degree)?a:b;//which face i am at
|
|
|
+ int iv = (k%(2*half_degree))/2;//which vector (0..half_degree-1) am i at
|
|
|
+ int ixy = k%2;//which part (real/imag) am i at
|
|
|
+ rows(ind) = startRow+j;
|
|
|
+ columns(ind) = 2*half_degree*f + 2*iv +ixy;
|
|
|
+ }
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+computeJacobianPattern()
|
|
|
+{
|
|
|
+ num_residuals_smooth = 4*numInteriorEdges;
|
|
|
+ num_residuals_close = 4*numF;
|
|
|
+ num_residuals_polycurl = 2*numInteriorEdges;
|
|
|
+ num_residuals_quotcurl = numInteriorEdges;
|
|
|
+ num_residuals_barrier = numF;
|
|
|
+
|
|
|
+ num_residuals = num_residuals_smooth + num_residuals_close + num_residuals_polycurl + num_residuals_quotcurl + num_residuals_barrier;
|
|
|
+
|
|
|
+ residuals.setZero(num_residuals,1);
|
|
|
+
|
|
|
+ //per edge
|
|
|
+ numJacElements_smooth = numInteriorEdges*numInnerJacCols_edge*numInnerJacRows_smooth;
|
|
|
+ numJacElements_polycurl = numInteriorEdges*numInnerJacCols_edge*numInnerJacRows_polycurl;
|
|
|
+ numJacElements_quotcurl = numInteriorEdges*numInnerJacCols_edge*numInnerJacRows_quotcurl;
|
|
|
+
|
|
|
+ //per face
|
|
|
+ numJacElements_barrier = numF*numInnerJacCols_face*numInnerJacRows_barrier;
|
|
|
+ numJacElements_close = numF*numInnerJacCols_face*numInnerJacRows_close;
|
|
|
+
|
|
|
+ numJacElements = (numJacElements_smooth +numJacElements_polycurl + numJacElements_quotcurl) + (numJacElements_barrier +numJacElements_close);
|
|
|
+ //allocate
|
|
|
+ II_Jac.setZero(numJacElements);
|
|
|
+ JJ_Jac.setZero(numJacElements);
|
|
|
+ SS_Jac.setOnes(numJacElements);
|
|
|
+
|
|
|
+ //set stuff (attention: order !)
|
|
|
+ int startRowInJacobian = 0;
|
|
|
+ int startIndexInVectors = 0;
|
|
|
+
|
|
|
+ //smoothness
|
|
|
+ add_jac_indices_edge(numInnerJacRows_smooth,
|
|
|
+ numInnerJacCols_edge,
|
|
|
+ startRowInJacobian,
|
|
|
+ startIndexInVectors,
|
|
|
+ II_Jac,
|
|
|
+ JJ_Jac);
|
|
|
+ startRowInJacobian += num_residuals_smooth;
|
|
|
+ startIndexInVectors += numJacElements_smooth;
|
|
|
+
|
|
|
+ //closeness
|
|
|
+ add_jac_indices_face(numInnerJacRows_close,
|
|
|
+ numInnerJacCols_face,
|
|
|
+ startRowInJacobian,
|
|
|
+ startIndexInVectors,
|
|
|
+ II_Jac,
|
|
|
+ JJ_Jac);
|
|
|
+ startRowInJacobian += num_residuals_close;
|
|
|
+ startIndexInVectors += numJacElements_close;
|
|
|
+
|
|
|
+ //barrier
|
|
|
+ add_jac_indices_face(numInnerJacRows_barrier,
|
|
|
+ numInnerJacCols_face,
|
|
|
+ startRowInJacobian,
|
|
|
+ startIndexInVectors,
|
|
|
+ II_Jac,
|
|
|
+ JJ_Jac);
|
|
|
+ startRowInJacobian += num_residuals_barrier;
|
|
|
+ startIndexInVectors += numJacElements_barrier;
|
|
|
+
|
|
|
+ //polycurl
|
|
|
+ add_jac_indices_edge(numInnerJacRows_polycurl,
|
|
|
+ numInnerJacCols_edge,
|
|
|
+ startRowInJacobian,
|
|
|
+ startIndexInVectors,
|
|
|
+ II_Jac,
|
|
|
+ JJ_Jac);
|
|
|
+ startRowInJacobian += num_residuals_polycurl;
|
|
|
+ startIndexInVectors += numJacElements_polycurl;
|
|
|
+
|
|
|
+ //quotcurl
|
|
|
+ add_jac_indices_edge(numInnerJacRows_quotcurl,
|
|
|
+ numInnerJacCols_edge,
|
|
|
+ startRowInJacobian,
|
|
|
+ startIndexInVectors,
|
|
|
+ II_Jac,
|
|
|
+ JJ_Jac);
|
|
|
+ igl::sparse(II_Jac, JJ_Jac, SS_Jac, Jac);
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+computeHessianPattern()
|
|
|
+{
|
|
|
+ //II_Jac is sorted in ascending order already
|
|
|
+ int starti = 0;
|
|
|
+ int currI = II_Jac(0);
|
|
|
+ for (int ii = 0; ii<II_Jac.rows(); ++ii)
|
|
|
+ {
|
|
|
+ if(currI != II_Jac(ii))
|
|
|
+ {
|
|
|
+ starti = ii;
|
|
|
+ currI = II_Jac(ii);
|
|
|
+ }
|
|
|
+ int k1 = II_Jac(ii);
|
|
|
+ for (int jj = starti; jj<II_Jac.rows(); ++jj)
|
|
|
+ {
|
|
|
+ int k2 = II_Jac(jj);
|
|
|
+ if (k1 !=k2)
|
|
|
+ break;
|
|
|
+ indInSS_Hess_1_vec.push_back(ii);
|
|
|
+ indInSS_Hess_2_vec.push_back(jj);
|
|
|
+ Hess_triplets.push_back(Eigen::Triplet<double> (JJ_Jac(ii),
|
|
|
+ JJ_Jac(jj),
|
|
|
+ SS_Jac(ii)*SS_Jac(jj)
|
|
|
+ )
|
|
|
+ );
|
|
|
+ }
|
|
|
+ }
|
|
|
+ Hess.resize(Jac.cols(),Jac.cols());
|
|
|
+ Hess.setFromTriplets(Hess_triplets.begin(), Hess_triplets.end());
|
|
|
+ Hess.makeCompressed();
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+computeNewHessValues()
|
|
|
+{
|
|
|
+ for (int i =0; i<Hess_triplets.size(); ++i)
|
|
|
+ Hess_triplets[i] = Eigen::Triplet<double>(Hess_triplets[i].row(),
|
|
|
+ Hess_triplets[i].col(),
|
|
|
+ SS_Jac(indInSS_Hess_1_vec[i])*SS_Jac(indInSS_Hess_2_vec[i])
|
|
|
+ );
|
|
|
+
|
|
|
+ Hess.setFromTriplets(Hess_triplets.begin(), Hess_triplets.end());
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC>::IntegrableFieldSolver( IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC> &cffsoldata):
|
|
|
+data(cffsoldata)
|
|
|
+{ };
|
|
|
+
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE bool igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+solve(igl::integrable_polyvector_fields_parameters ¶ms,
|
|
|
+ Eigen::PlainObjectBase<DerivedFF>& current_field,
|
|
|
+ bool current_field_is_not_ccw)
|
|
|
+{
|
|
|
+ Eigen::MatrixXd sol2D;
|
|
|
+ Eigen::MatrixXd sol3D = current_field.template cast<double>();
|
|
|
+ if (current_field_is_not_ccw)
|
|
|
+ data.makeFieldCCW(sol3D);
|
|
|
+
|
|
|
+ igl::global2local(data.B1, data.B2, sol3D, sol2D);
|
|
|
+ Eigen::VectorXd x;
|
|
|
+ x.setZero(data.numVariables);
|
|
|
+ for (int i =0; i<data.numF; i++)
|
|
|
+ x.segment(i*2*2, 2*2) = sol2D.row(i);
|
|
|
+
|
|
|
+ //get x
|
|
|
+ solveGaussNewton(params, data.xOriginal, x);
|
|
|
+ //get output from x
|
|
|
+ for (int i =0; i<data.numF; i++)
|
|
|
+ sol2D.row(i) = x.segment(i*2*2, 2*2);
|
|
|
+ igl::local2global(data.B1, data.B2, sol2D, sol3D);
|
|
|
+ current_field = sol3D.cast<typename DerivedFF::Scalar>();
|
|
|
+ return true;
|
|
|
+}
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+solveGaussNewton(integrable_polyvector_fields_parameters ¶ms,
|
|
|
+ const Eigen::VectorXd &x_initial,
|
|
|
+ Eigen::VectorXd &x)
|
|
|
+{
|
|
|
+ bool converged = false;
|
|
|
+
|
|
|
+ double F;
|
|
|
+ Eigen::VectorXd xprev = x;
|
|
|
+ Eigen::VectorXd xc = igl::slice(x_initial, data.constrained, 1);
|
|
|
+ // double ESmooth, EClose, ECurl, EQuotCurl, EBarrier;
|
|
|
+ for (int innerIter = 0; innerIter<params.numIter; ++innerIter)
|
|
|
+ {
|
|
|
+
|
|
|
+ //set constrained entries to those of the initial
|
|
|
+ igl::slice_into(xc, data.constrained, 1, xprev);
|
|
|
+
|
|
|
+
|
|
|
+ //get function, gradients and Hessians
|
|
|
+ F = RJ(x, xprev, params, true);
|
|
|
+
|
|
|
+ printf("IntegrableFieldSolver -- Iteration %d\n", innerIter);
|
|
|
+
|
|
|
+ if((data.residuals.array() == std::numeric_limits<double>::infinity()).any())
|
|
|
+ {
|
|
|
+ std::cerr<<"IntegrableFieldSolver -- residuals: got infinity somewhere"<<std::endl;
|
|
|
+ exit(1);
|
|
|
+ };
|
|
|
+ if((data.residuals.array() != data.residuals.array()).any())
|
|
|
+ {
|
|
|
+ std::cerr<<"IntegrableFieldSolver -- residuals: got infinity somewhere"<<std::endl;
|
|
|
+ exit(1);
|
|
|
+ };
|
|
|
+
|
|
|
+ converged = false;
|
|
|
+
|
|
|
+ Eigen::VectorXd rhs = data.Jac.transpose()*data.residuals;
|
|
|
+
|
|
|
+ bool success;
|
|
|
+ data.solver.factorize(data.Hess);
|
|
|
+ success = data.solver.info() == Eigen::Success;
|
|
|
+
|
|
|
+ if(!success)
|
|
|
+ std::cerr<<"IntegrableFieldSolver -- Could not do LU"<<std::endl;
|
|
|
+
|
|
|
+ Eigen::VectorXd direction;
|
|
|
+
|
|
|
+ double error;
|
|
|
+ direction = data.solver.solve(rhs);
|
|
|
+ error = (data.Hess*direction - rhs).cwiseAbs().maxCoeff();
|
|
|
+ if(error> 1e-4)
|
|
|
+ {
|
|
|
+ std::cerr<<"IntegrableFieldSolver -- Could not solve"<<std::endl;
|
|
|
+ }
|
|
|
+
|
|
|
+ // adaptive backtracking
|
|
|
+ bool repeat = true;
|
|
|
+ int run = 0;
|
|
|
+ Eigen::VectorXd cx;
|
|
|
+ Eigen::VectorXd tRes;
|
|
|
+ double newF;
|
|
|
+ while(repeat)
|
|
|
+ {
|
|
|
+ cx = x - params.gamma*direction;
|
|
|
+ newF = RJ(cx, xprev, params);
|
|
|
+ if(newF < F)
|
|
|
+ {
|
|
|
+ repeat = false;
|
|
|
+ if(run == 0)
|
|
|
+ params.gamma *= 2;
|
|
|
+ }
|
|
|
+ else
|
|
|
+ {
|
|
|
+ params.gamma *= 0.5f;
|
|
|
+ if(params.gamma<1e-30)
|
|
|
+ {
|
|
|
+ repeat = false;
|
|
|
+ converged = true;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ run++;
|
|
|
+ }
|
|
|
+
|
|
|
+
|
|
|
+ if (!converged)
|
|
|
+ {
|
|
|
+ xprev = x;
|
|
|
+ x = cx;
|
|
|
+ }
|
|
|
+ else
|
|
|
+ {
|
|
|
+ std::cerr<<"IntegrableFieldSolver -- Converged"<<std::endl;
|
|
|
+ break;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+
|
|
|
+}
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE double igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+RJ(const Eigen::VectorXd &x,
|
|
|
+ const Eigen::VectorXd &x0,
|
|
|
+ const integrable_polyvector_fields_parameters ¶ms,
|
|
|
+ bool doJacs)
|
|
|
+{
|
|
|
+ Eigen::MatrixXd sol2D(data.numF,4), sol02D(data.numF,4);
|
|
|
+ for (int i =0; i<data.numF; i++)
|
|
|
+ sol2D.row(i) = x.segment(i*2*2, 2*2);
|
|
|
+ for (int i =0; i<data.numF; i++)
|
|
|
+ sol02D.row(i) = x0.segment(i*2*2, 2*2);
|
|
|
+
|
|
|
+ data.SS_Jac.setZero(data.numJacElements);
|
|
|
+
|
|
|
+ //set stuff (attention: order !)
|
|
|
+ int startRowInJacobian = 0;
|
|
|
+ int startIndexInVectors = 0;
|
|
|
+
|
|
|
+ //smoothness
|
|
|
+ RJ_Smoothness(sol2D, sqrt(params.wSmooth), startRowInJacobian, doJacs, startIndexInVectors);
|
|
|
+ startRowInJacobian += data.num_residuals_smooth;
|
|
|
+ startIndexInVectors += data.numJacElements_smooth;
|
|
|
+
|
|
|
+ //closeness
|
|
|
+ RJ_Closeness(sol2D, sol02D, sqrt(params.wCloseUnconstrained), sqrt(params.wCloseConstrained), params.do_partial, startRowInJacobian, doJacs, startIndexInVectors);
|
|
|
+ startRowInJacobian += data.num_residuals_close;
|
|
|
+ startIndexInVectors += data.numJacElements_close;
|
|
|
+
|
|
|
+ //barrier
|
|
|
+ RJ_Barrier(sol2D, params.sBarrier, sqrt(params.wBarrier), startRowInJacobian, doJacs, startIndexInVectors);
|
|
|
+ startRowInJacobian += data.num_residuals_barrier;
|
|
|
+ startIndexInVectors += data.numJacElements_barrier;
|
|
|
+
|
|
|
+ //polycurl
|
|
|
+ RJ_Curl(sol2D, params.wCurl, powl(params.wCurl, 2), startRowInJacobian, doJacs, startIndexInVectors);
|
|
|
+ startRowInJacobian += data.num_residuals_polycurl;
|
|
|
+ startIndexInVectors += data.numJacElements_polycurl;
|
|
|
+
|
|
|
+ //quotcurl
|
|
|
+ RJ_QuotCurl(sol2D, sqrt(params.wQuotCurl), startRowInJacobian, doJacs, startIndexInVectors);
|
|
|
+
|
|
|
+ if(doJacs)
|
|
|
+ {
|
|
|
+ for (int i =0; i<data.numJacElements; ++i)
|
|
|
+ data.Jac.coeffRef(data.II_Jac(i), data.JJ_Jac(i)) = data.SS_Jac(i);
|
|
|
+ data.computeNewHessValues();
|
|
|
+ }
|
|
|
+
|
|
|
+ return data.residuals.transpose()*data.residuals;
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+rj_smoothness_edge(const Eigen::RowVectorXd &vec2D_a,
|
|
|
+ const Eigen::RowVectorXd &vec2D_b,
|
|
|
+ const double &k,
|
|
|
+ const int nA,
|
|
|
+ const int nB,
|
|
|
+ Eigen::VectorXd &residuals,
|
|
|
+ bool do_jac,
|
|
|
+ Eigen::MatrixXd &Jac)
|
|
|
+{
|
|
|
+ const Eigen::RowVector2d &ua = vec2D_a.segment(0, 2);
|
|
|
+ const Eigen::RowVector2d &va = vec2D_a.segment(2, 2);
|
|
|
+ const Eigen::RowVector2d &ub = vec2D_b.segment(0, 2);
|
|
|
+ const Eigen::RowVector2d &vb = vec2D_b.segment(2, 2);
|
|
|
+ const double &xua=ua[0], &yua=ua[1], &xva=va[0], &yva=va[1];
|
|
|
+ const double &xub=ub[0], &yub=ub[1], &xvb=vb[0], &yvb=vb[1];
|
|
|
+
|
|
|
+ double xua_2 = xua*xua;
|
|
|
+ double xva_2 = xva*xva;
|
|
|
+ double yua_2 = yua*yua;
|
|
|
+ double yva_2 = yva*yva;
|
|
|
+ double xub_2 = xub*xub;
|
|
|
+ double xvb_2 = xvb*xvb;
|
|
|
+ double yub_2 = yub*yub;
|
|
|
+ double yvb_2 = yvb*yvb;
|
|
|
+
|
|
|
+ double sA = sin(nA*k);
|
|
|
+ double cA = cos(nA*k);
|
|
|
+ double sB = sin(nB*k);
|
|
|
+ double cB = cos(nB*k);
|
|
|
+
|
|
|
+ double t1 = xua*yua;
|
|
|
+ double t2 = xva*yva;
|
|
|
+ double t3 = xub*xvb;
|
|
|
+ double t4 = yub*yvb;
|
|
|
+ double t5 = xua*xva;
|
|
|
+ double t6 = xub*yub;
|
|
|
+ double t7 = yua*yva;
|
|
|
+ double t8 = xvb*yvb;
|
|
|
+
|
|
|
+ double t9 = xva_2 - yva_2;
|
|
|
+ double t10 = xua_2 - yua_2;
|
|
|
+ double t11 = xvb_2 - yvb_2;
|
|
|
+ double t12 = xub_2 - yub_2;
|
|
|
+
|
|
|
+ double t13 = 2*t1 + 2*t2;
|
|
|
+
|
|
|
+ double t17 = (2*t1*t9 + 2*t2*t10);
|
|
|
+ double t19 = (t10*t9 - 4*t5*t7);
|
|
|
+
|
|
|
+
|
|
|
+ residuals.resize(4, 1);
|
|
|
+ residuals <<
|
|
|
+ cA*(t10 + t9) - sA*(t13) - t12 - t11,
|
|
|
+ sA*(t10 + t9) - 2*t8 - 2*t6 + cA*(t13),
|
|
|
+ cB*t19 - (t12)*(t11) - sB*t17 + 4*t3*t4,
|
|
|
+ cB*t17 + sB*t19 - 2*t6*(t11) - 2*t8*(t12);
|
|
|
+
|
|
|
+
|
|
|
+ if (do_jac)
|
|
|
+ {
|
|
|
+ double t20 = 2*yua*t9 + 4*xua*t2;
|
|
|
+ double t21 = 2*xua*t9 - 4*xva*t7;
|
|
|
+ double t22 = 2*yva*t10 + 4*t5*yua;
|
|
|
+ double t23 = 2*xva*t10 - 4*t1*yva;
|
|
|
+
|
|
|
+ Jac.resize(4,8);
|
|
|
+ Jac << 2*xua*cA - 2*yua*sA, - 2*yua*cA - 2*xua*sA, 2*xva*cA - 2*yva*sA, - 2*yva*cA - 2*xva*sA, -2*xub, 2*yub, -2*xvb, 2*yvb,
|
|
|
+ 2*yua*cA + 2*xua*sA, 2*xua*cA - 2*yua*sA, 2*yva*cA + 2*xva*sA, 2*xva*cA - 2*yva*sA, -2*yub, -2*xub, -2*yvb, -2*xvb,
|
|
|
+ cB*(t21) - sB*(t20), - cB*(t20) - sB*(t21), cB*(t23) - sB*(t22), - cB*(t22) - sB*(t23), 4*xvb*t4 - 2*xub*t11, 2*yub*t11 + 4*t3*yvb, 4*xub*t4 - 2*xvb*t12, 2*yvb*t12 + 4*t3*yub,
|
|
|
+ cB*(t20) + sB*(t21), cB*(t21) - sB*(t20), cB*(t22) + sB*(t23), cB*(t23) - sB*(t22), - 2*yub*t11 - 4*t3*yvb, 4*xvb*t4 - 2*xub*t11, - 2*yvb*t12 - 4*t3*yub, 4*xub*t4 - 2*xvb*t12;
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+RJ_Smoothness(const Eigen::MatrixXd &sol2D,
|
|
|
+ const double &wSmoothSqrt,
|
|
|
+ const int startRowInJacobian,
|
|
|
+ bool doJacs,
|
|
|
+ const int startIndexInVectors)
|
|
|
+{
|
|
|
+ if (wSmoothSqrt ==0)
|
|
|
+ return;
|
|
|
+ for (int ii=0; ii<data.numInteriorEdges; ++ii)
|
|
|
+ {
|
|
|
+ // the two faces of the flap
|
|
|
+ int a = data.E2F_int(ii,0);
|
|
|
+ int b = data.E2F_int(ii,1);
|
|
|
+
|
|
|
+ int k = data.indInteriorToFull[ii];
|
|
|
+
|
|
|
+ Eigen::MatrixXd tJac;
|
|
|
+ Eigen::VectorXd tRes;
|
|
|
+ rj_smoothness_edge(sol2D.row(a),
|
|
|
+ sol2D.row(b),
|
|
|
+ data.K[k],
|
|
|
+ 2*(0+1), //degree of first coefficient
|
|
|
+ 2*(1+1), //degree of second coefficient
|
|
|
+ tRes,
|
|
|
+ doJacs,
|
|
|
+ tJac);
|
|
|
+
|
|
|
+ int startRow = startRowInJacobian+data.numInnerJacRows_smooth*ii;
|
|
|
+ data.residuals.segment(startRow,data.numInnerJacRows_smooth) = wSmoothSqrt*tRes;
|
|
|
+
|
|
|
+ if(doJacs)
|
|
|
+ {
|
|
|
+ int startIndex = startIndexInVectors+data.numInnerJacRows_smooth*data.numInnerJacCols_edge*ii;
|
|
|
+ data.add_Jacobian_to_svector(startIndex, wSmoothSqrt*tJac,data.SS_Jac);
|
|
|
+ }
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+rj_barrier_face(const Eigen::RowVectorXd &vec2D_a,
|
|
|
+ const double &s,
|
|
|
+ Eigen::VectorXd &residuals,
|
|
|
+ bool do_jac,
|
|
|
+ Eigen::MatrixXd &Jac)
|
|
|
+{
|
|
|
+
|
|
|
+ const Eigen::RowVector2d &ua = vec2D_a.segment(0, 2);
|
|
|
+ const Eigen::RowVector2d &va = vec2D_a.segment(2, 2);
|
|
|
+ const double &xua=ua[0], &yua=ua[1], &xva=va[0], &yva=va[1];
|
|
|
+
|
|
|
+
|
|
|
+ double xva_2 = xva*xva;
|
|
|
+ double yua_2 = yua*yua;
|
|
|
+ double xua_2 = xua*xua;
|
|
|
+ double yva_2 = yva*yva;
|
|
|
+
|
|
|
+ double s_2 = s*s;
|
|
|
+ double s_3 = s*s_2;
|
|
|
+
|
|
|
+ double t00 = xua*yva;
|
|
|
+ double t01 = xva*yua;
|
|
|
+ double t05 = t00 - t01;
|
|
|
+ double t05_2 = t05*t05;
|
|
|
+ double t05_3 = t05*t05_2;
|
|
|
+
|
|
|
+ if (do_jac)
|
|
|
+ Jac.setZero(1,4);
|
|
|
+ residuals.resize(1, 1);
|
|
|
+ if (t05>=s)
|
|
|
+ residuals << 0;
|
|
|
+ else if (t05<0)
|
|
|
+ residuals << std::numeric_limits<double>::infinity();
|
|
|
+ else
|
|
|
+ {
|
|
|
+ residuals << 1/((3*t00 - 3*t01)/s - (3*t05_2)/s_2 + t05_3/s_3) - 1;
|
|
|
+ double t03 = (s - t05)*(s - t05);
|
|
|
+ double t06 = (3*s_2 - 3*s*t00 + 3*s*t01 + xua_2*yva_2 - 2*xua*t01*yva + xva_2*yua_2);
|
|
|
+ double t04 = t06*t06;
|
|
|
+ if (do_jac)
|
|
|
+ Jac<<
|
|
|
+ -(3*s_3*yva*t03)/(t05_2*t04),
|
|
|
+ (3*s_3*xva*t03)/(t05_2*t04),
|
|
|
+ (3*s_3*yua*t03)/(t05_2*t04),
|
|
|
+ -(3*s_3*xua*t03)/(t05_2*t04);
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+RJ_Barrier(const Eigen::MatrixXd &sol2D,
|
|
|
+ const double &s,
|
|
|
+ const double &wBarrierSqrt,
|
|
|
+ const int startRowInJacobian,
|
|
|
+ bool doJacs,
|
|
|
+ const int startIndexInVectors)
|
|
|
+{
|
|
|
+ if (wBarrierSqrt ==0)
|
|
|
+ return;
|
|
|
+
|
|
|
+ for (int fi=0; fi<data.numF; ++fi)
|
|
|
+ {
|
|
|
+ Eigen::MatrixXd tJac;
|
|
|
+ Eigen::VectorXd tRes;
|
|
|
+ rj_barrier_face(sol2D.row(fi),
|
|
|
+ s,
|
|
|
+ tRes,
|
|
|
+ doJacs,
|
|
|
+ tJac);
|
|
|
+
|
|
|
+ int startRow = startRowInJacobian+ data.numInnerJacRows_barrier * fi;
|
|
|
+ data.residuals.segment(startRow,data.numInnerJacRows_barrier) = wBarrierSqrt*tRes;
|
|
|
+
|
|
|
+ if(doJacs)
|
|
|
+ {
|
|
|
+ int startIndex = startIndexInVectors+data.numInnerJacRows_barrier*data.numInnerJacCols_face*fi;
|
|
|
+ data.add_Jacobian_to_svector(startIndex, wBarrierSqrt*tJac,data.SS_Jac);
|
|
|
+ }
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+RJ_Closeness(const Eigen::MatrixXd &sol2D,
|
|
|
+ const Eigen::MatrixXd &sol02D,
|
|
|
+ const double &wCloseUnconstrainedSqrt,
|
|
|
+ const double &wCloseConstrainedSqrt,
|
|
|
+ const bool do_partial,
|
|
|
+ const int startRowInJacobian,
|
|
|
+ bool doJacs,
|
|
|
+ const int startIndexInVectors)
|
|
|
+{
|
|
|
+ if (wCloseUnconstrainedSqrt ==0 && wCloseConstrainedSqrt ==0)
|
|
|
+ return;
|
|
|
+ for (int fi=0; fi<data.numF; ++fi)
|
|
|
+ {
|
|
|
+ Eigen::Vector4d weights;
|
|
|
+ if (!data.is_constrained_face[fi])
|
|
|
+ weights.setConstant(wCloseUnconstrainedSqrt);
|
|
|
+ else
|
|
|
+ {
|
|
|
+ if (do_partial&&data.is_constrained_face[fi]==1)
|
|
|
+ //only constrain the first vector
|
|
|
+ weights<<wCloseConstrainedSqrt,wCloseConstrainedSqrt,wCloseUnconstrainedSqrt,wCloseUnconstrainedSqrt;
|
|
|
+ else
|
|
|
+ //either not partial, or partial with 2 constraints
|
|
|
+ weights.setConstant(wCloseConstrainedSqrt);
|
|
|
+ }
|
|
|
+
|
|
|
+ Eigen::MatrixXd tJac;
|
|
|
+ Eigen::VectorXd tRes;
|
|
|
+ tJac = Eigen::MatrixXd::Identity(4,4);
|
|
|
+ tRes.resize(4, 1); tRes<<(sol2D.row(fi)-sol02D.row(fi)).transpose();
|
|
|
+ int startRow = startRowInJacobian+data.numInnerJacRows_close*fi;
|
|
|
+ data.residuals.segment(startRow,data.numInnerJacRows_close) = weights.array()*tRes.array();
|
|
|
+
|
|
|
+ if(doJacs)
|
|
|
+ {
|
|
|
+ int startIndex = startIndexInVectors+data.numInnerJacRows_close*data.numInnerJacCols_face*fi;
|
|
|
+ data.add_Jacobian_to_svector(startIndex, weights.asDiagonal()*tJac,data.SS_Jac);
|
|
|
+ }
|
|
|
+
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+rj_polycurl_edge(const Eigen::RowVectorXd &vec2D_a,
|
|
|
+ const Eigen::RowVector2d &ea,
|
|
|
+ const Eigen::RowVectorXd &vec2D_b,
|
|
|
+ const Eigen::RowVector2d &eb,
|
|
|
+ Eigen::VectorXd &residuals,
|
|
|
+ bool do_jac,
|
|
|
+ Eigen::MatrixXd &Jac)
|
|
|
+{
|
|
|
+ const Eigen::RowVector2d &ua = vec2D_a.segment(0, 2);
|
|
|
+ const Eigen::RowVector2d &va = vec2D_a.segment(2, 2);
|
|
|
+ const Eigen::RowVector2d &ub = vec2D_b.segment(0, 2);
|
|
|
+ const Eigen::RowVector2d &vb = vec2D_b.segment(2, 2);
|
|
|
+ const double &xua=ua[0], &yua=ua[1], &xva=va[0], &yva=va[1];
|
|
|
+ const double &xub=ub[0], &yub=ub[1], &xvb=vb[0], &yvb=vb[1];
|
|
|
+ const double &xea=ea[0], &yea=ea[1];
|
|
|
+ const double &xeb=eb[0], &yeb=eb[1];
|
|
|
+
|
|
|
+ const double dua = (xea*xua + yea*yua);
|
|
|
+ const double dub = (xeb*xub + yeb*yub);
|
|
|
+ const double dva = (xea*xva + yea*yva);
|
|
|
+ const double dvb = (xeb*xvb + yeb*yvb);
|
|
|
+
|
|
|
+ const double dua_2 = dua*dua;
|
|
|
+ const double dva_2 = dva*dva;
|
|
|
+ const double dub_2 = dub*dub;
|
|
|
+ const double dvb_2 = dvb*dvb;
|
|
|
+
|
|
|
+ residuals.resize(2, 1);
|
|
|
+ residuals << dua_2 - dub_2 + dva_2 - dvb_2,
|
|
|
+ dua_2*dva_2 - dub_2*dvb_2 ;
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+ if (do_jac)
|
|
|
+ {
|
|
|
+
|
|
|
+ Jac.resize(2,8);
|
|
|
+ Jac << 2*xea*dua, 2*yea*dua, 2*xea*dva, 2*yea*dva, -2*xeb*dub, -2*yeb*dub, -2*xeb*dvb, -2*yeb*dvb,
|
|
|
+ 2*xea*dua*dva_2, 2*yea*dua*dva_2, 2*xea*dua_2*dva, 2*yea*dua_2*dva, -2*xeb*dub*dvb_2, -2*yeb*dub*dvb_2, -2*xeb*dub_2*dvb, -2*yeb*dub_2*dvb;
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+RJ_Curl(const Eigen::MatrixXd &sol2D,
|
|
|
+ const double &wCASqrt,
|
|
|
+ const double &wCBSqrt,
|
|
|
+ const int startRowInJacobian,
|
|
|
+ bool doJacs,
|
|
|
+ const int startIndexInVectors)
|
|
|
+{
|
|
|
+ if((wCASqrt==0) &&(wCBSqrt==0))
|
|
|
+ return;
|
|
|
+ for (int ii=0; ii<data.numInteriorEdges; ++ii)
|
|
|
+ {
|
|
|
+ // the two faces of the flap
|
|
|
+ int a = data.E2F_int(ii,0);
|
|
|
+ int b = data.E2F_int(ii,1);
|
|
|
+
|
|
|
+ int k = data.indInteriorToFull[ii];
|
|
|
+
|
|
|
+ // the common edge, a 3 vector
|
|
|
+ const Eigen::RowVector3d &ce = data.EVecNorm.row(k);
|
|
|
+
|
|
|
+ // the common edge expressed in local coordinates in the two faces
|
|
|
+ // x/y denotes real/imaginary
|
|
|
+ double xea = data.B1.row(a).dot(ce);
|
|
|
+ double yea = data.B2.row(a).dot(ce);
|
|
|
+ Eigen::RowVector2d ea; ea<<xea, yea;
|
|
|
+ double xeb = data.B1.row(b).dot(ce);
|
|
|
+ double yeb = data.B2.row(b).dot(ce);
|
|
|
+ Eigen::RowVector2d eb; eb<<xeb, yeb;
|
|
|
+
|
|
|
+
|
|
|
+ Eigen::MatrixXd tJac;
|
|
|
+ Eigen::VectorXd tRes;
|
|
|
+ rj_polycurl_edge(sol2D.row(a),
|
|
|
+ ea,
|
|
|
+ sol2D.row(b),
|
|
|
+ eb,
|
|
|
+ tRes,
|
|
|
+ doJacs,
|
|
|
+ tJac);
|
|
|
+
|
|
|
+ tRes[0] = tRes[0]*wCASqrt;
|
|
|
+ tRes[1] = tRes[1]*wCBSqrt;
|
|
|
+
|
|
|
+
|
|
|
+ int startRow = startRowInJacobian+data.numInnerJacRows_polycurl*ii;
|
|
|
+ data.residuals.segment(startRow,data.numInnerJacRows_polycurl) = tRes;
|
|
|
+
|
|
|
+ if(doJacs)
|
|
|
+ {
|
|
|
+ tJac.row(0) = tJac.row(0)*wCASqrt;
|
|
|
+ tJac.row(1) = tJac.row(1)*wCBSqrt;
|
|
|
+ int startIndex = startIndexInVectors+data.numInnerJacRows_polycurl*data.numInnerJacCols_edge*ii;
|
|
|
+ data.add_Jacobian_to_svector(startIndex, tJac,data.SS_Jac);
|
|
|
+ }
|
|
|
+
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+rj_quotcurl_edge_polyversion(const Eigen::RowVectorXd &vec2D_a,
|
|
|
+ const Eigen::RowVector2d &ea,
|
|
|
+ const Eigen::RowVectorXd &vec2D_b,
|
|
|
+ const Eigen::RowVector2d &eb,
|
|
|
+ Eigen::VectorXd &residuals,
|
|
|
+ bool do_jac,
|
|
|
+ Eigen::MatrixXd &Jac)
|
|
|
+{
|
|
|
+ const Eigen::RowVector2d &ua = vec2D_a.segment(0, 2);
|
|
|
+ const Eigen::RowVector2d &va = vec2D_a.segment(2, 2);
|
|
|
+ const Eigen::RowVector2d &ub = vec2D_b.segment(0, 2);
|
|
|
+ const Eigen::RowVector2d &vb = vec2D_b.segment(2, 2);
|
|
|
+ const double &xua=ua[0], &yua=ua[1], &xva=va[0], &yva=va[1];
|
|
|
+ const double &xub=ub[0], &yub=ub[1], &xvb=vb[0], &yvb=vb[1];
|
|
|
+ const double &xea=ea[0], &yea=ea[1];
|
|
|
+ const double &xeb=eb[0], &yeb=eb[1];
|
|
|
+
|
|
|
+ double dua = (xea*xua + yea*yua);
|
|
|
+ double dva = (xea*xva + yea*yva);
|
|
|
+ double dub = (xeb*xub + yeb*yub);
|
|
|
+ double dvb = (xeb*xvb + yeb*yvb);
|
|
|
+
|
|
|
+ double dua_2 = dua * dua;
|
|
|
+ double dva_2 = dva * dva;
|
|
|
+ double dub_2 = dub * dub;
|
|
|
+ double dvb_2 = dvb * dvb;
|
|
|
+ double t00 = (dub_2 - dvb_2);
|
|
|
+ double t01 = (dua_2 - dva_2);
|
|
|
+
|
|
|
+
|
|
|
+ residuals.resize(1, 1);
|
|
|
+ residuals << dua*dva*t00 - dub*dvb*t01;
|
|
|
+
|
|
|
+ if (do_jac)
|
|
|
+ {
|
|
|
+ Jac.resize(1,8);
|
|
|
+ Jac << xea*dva*t00 - 2*xea*dua*dub*dvb, yea*dva*t00 - 2*yea*dua*dub*dvb, xea*dua*t00 + 2*xea*dub*dva*dvb, yea*dua*t00 + 2*yea*dub*dva*dvb, 2*xeb*dua*dub*dva - xeb*dvb*t01, 2*yeb*dua*dub*dva - yeb*dvb*t01, - xeb*dub*t01 - 2*xeb*dua*dva*dvb, - yeb*dub*t01 - 2*yeb*dua*dva*dvb;
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC>::
|
|
|
+RJ_QuotCurl(const Eigen::MatrixXd &sol2D,
|
|
|
+ const double &wQuotCurlSqrt,
|
|
|
+ const int startRowInJacobian,
|
|
|
+ bool doJacs,
|
|
|
+ const int startIndexInVectors)
|
|
|
+{
|
|
|
+ for (int ii=0; ii<data.numInteriorEdges; ++ii)
|
|
|
+ {
|
|
|
+ // the two faces of the flap
|
|
|
+ int a = data.E2F_int(ii,0);
|
|
|
+ int b = data.E2F_int(ii,1);
|
|
|
+
|
|
|
+ int k = data.indInteriorToFull[ii];
|
|
|
+
|
|
|
+ // the common edge, a 3 vector
|
|
|
+ const Eigen::RowVector3d &ce = data.EVecNorm.row(k);
|
|
|
+
|
|
|
+ // the common edge expressed in local coordinates in the two faces
|
|
|
+ // x/y denotes real/imaginary
|
|
|
+ double xea = data.B1.row(a).dot(ce);
|
|
|
+ double yea = data.B2.row(a).dot(ce);
|
|
|
+ Eigen::RowVector2d ea; ea<<xea, yea;
|
|
|
+ double xeb = data.B1.row(b).dot(ce);
|
|
|
+ double yeb = data.B2.row(b).dot(ce);
|
|
|
+ Eigen::RowVector2d eb; eb<<xeb, yeb;
|
|
|
+
|
|
|
+
|
|
|
+ Eigen::MatrixXd tJac;
|
|
|
+ Eigen::VectorXd tRes;
|
|
|
+ rj_quotcurl_edge_polyversion(sol2D.row(a),
|
|
|
+ ea,
|
|
|
+ sol2D.row(b),
|
|
|
+ eb,
|
|
|
+ tRes,
|
|
|
+ doJacs,
|
|
|
+ tJac);
|
|
|
+
|
|
|
+ int startRow = startRowInJacobian+ data.numInnerJacRows_quotcurl*ii;
|
|
|
+ data.residuals.segment(startRow,data.numInnerJacRows_quotcurl) = wQuotCurlSqrt*tRes;
|
|
|
+
|
|
|
+ if(doJacs)
|
|
|
+ {
|
|
|
+ int startIndex = startIndexInVectors+data.numInnerJacRows_quotcurl*data.numInnerJacCols_edge*ii;
|
|
|
+ data.add_Jacobian_to_svector(startIndex, wQuotCurlSqrt*tJac,data.SS_Jac);
|
|
|
+ }
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::integrable_polyvector_fields_precompute(
|
|
|
+ const Eigen::PlainObjectBase<DerivedV>& V,
|
|
|
+ const Eigen::PlainObjectBase<DerivedF>& F,
|
|
|
+ const Eigen::VectorXi& b,
|
|
|
+ const Eigen::PlainObjectBase<DerivedC>& bc,
|
|
|
+ const Eigen::VectorXi& constraint_level,
|
|
|
+ const Eigen::PlainObjectBase<DerivedFF>& original_field,
|
|
|
+ igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC> &data)
|
|
|
+{
|
|
|
+ data.precomputeMesh(V,F);
|
|
|
+
|
|
|
+ data.computeJacobianPattern();
|
|
|
+ data.computeHessianPattern();
|
|
|
+ data.solver.analyzePattern(data.Hess);
|
|
|
+
|
|
|
+ data.initializeConstraints(b,bc,constraint_level);
|
|
|
+ data.initializeOriginalVariable(original_field);
|
|
|
+};
|
|
|
+
|
|
|
+
|
|
|
+template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
|
|
|
+IGL_INLINE void igl::integrable_polyvector_fields_solve(igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC> &cffsoldata,
|
|
|
+ integrable_polyvector_fields_parameters ¶ms,
|
|
|
+ Eigen::PlainObjectBase<DerivedFF>& current_field,
|
|
|
+ bool current_field_is_not_ccw)
|
|
|
+{
|
|
|
+ igl::IntegrableFieldSolver<DerivedV, DerivedF, DerivedFF, DerivedC> cffs(cffsoldata);
|
|
|
+ cffs.solve(params, current_field, current_field_is_not_ccw);
|
|
|
+};
|
|
|
+
|
|
|
+#ifdef IGL_STATIC_LIBRARY
|
|
|
+// Explicit template specialization
|
|
|
+template igl::IntegrableFieldSolverData<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >::IntegrableFieldSolverData();
|
|
|
+template void igl::integrable_polyvector_fields_solve<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(igl::IntegrableFieldSolverData<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, igl::integrable_polyvector_fields_parameters&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, bool);
|
|
|
+template void igl::integrable_polyvector_fields_precompute<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, igl::IntegrableFieldSolverData<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
|
|
|
+#endif
|
Olga Diamanti