1 /* 2 * Copyright (c) 2012-2016, Bruno Levy 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms, with or without 6 * modification, are permitted provided that the following conditions are met: 7 * 8 * * Redistributions of source code must retain the above copyright notice, 9 * this list of conditions and the following disclaimer. 10 * * Redistributions in binary form must reproduce the above copyright notice, 11 * this list of conditions and the following disclaimer in the documentation 12 * and/or other materials provided with the distribution. 13 * * Neither the name of the ALICE Project-Team nor the names of its 14 * contributors may be used to endorse or promote products derived from this 15 * software without specific prior written permission. 16 * 17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 27 * POSSIBILITY OF SUCH DAMAGE. 28 * 29 * If you modify this software, you should include a notice giving the 30 * name of the person performing the modification, the date of modification, 31 * and the reason for such modification. 32 * 33 * Contact: Bruno Levy 34 * 35 * Bruno.Levy@inria.fr 36 * http://www.loria.fr/~levy 37 * 38 * ALICE Project 39 * LORIA, INRIA Lorraine, 40 * Campus Scientifique, BP 239 41 * 54506 VANDOEUVRE LES NANCY CEDEX 42 * FRANCE 43 * 44 */ 45 46 #include <geogram/parameterization/mesh_PGP_2d.h> 47 #include <geogram/mesh/mesh.h> 48 #include <geogram/mesh/mesh_geometry.h> 49 #include <geogram/NL/nl.h> 50 #include <geogram/bibliography/bibliography.h> 51 52 #include <stack> 53 54 namespace { 55 using namespace GEO; 56 57 58 /** 59 * \brief computes gradients in a triangle in 2D. 60 */ 61 class ParamTrglGradient { 62 public: ParamTrglGradient(const vec2 & p1,const vec2 & p2,const vec2 & p3)63 ParamTrglGradient( 64 const vec2& p1, const vec2& p2, const vec2& p3 65 ) { 66 vertex_[0] = p1 ; 67 vertex_[1] = p2 ; 68 vertex_[2] = p3 ; 69 70 double x1 = p1.x ; 71 double y1 = p1.y ; 72 double x2 = p2.x ; 73 double y2 = p2.y ; 74 double x3 = p3.x ; 75 double y3 = p3.y ; 76 77 double d = x2*y3 - y2*x3 + x3*y1 - y3*x1 + x1*y2 - y1*x2 ; 78 79 if(fabs(d) < 1e-10) { 80 d = 1.0 ; 81 is_flat_ = true ; 82 } else { 83 is_flat_ = false ; 84 } 85 86 TX_[0] = (y2 - y3)/d ; 87 TX_[1] = (y3 - y1)/d ; 88 TX_[2] = (y1 - y2)/d ; 89 90 TY_[0] = -(x2 - x3)/d ; 91 TY_[1] = -(x3 - x1)/d ; 92 TY_[2] = -(x1 - x2)/d ; 93 } 94 TX(int i) const95 double TX(int i) const { 96 geo_debug_assert(i<3); 97 return TX_[i]; 98 } 99 TY(int i) const100 double TY(int i) const { 101 geo_debug_assert(i<3); 102 return TY_[i]; 103 } 104 is_flat() const105 bool is_flat() const { 106 return is_flat_ ; 107 } 108 109 private: 110 double TX_[3] ; 111 double TY_[3] ; 112 vec2 vertex_[3] ; 113 bool is_flat_ ; 114 } ; 115 116 117 /** 118 * \brief Retrieves a coordinate from an angle computed by PGP. 119 * \param[in] alpha the input variable. 120 * \param[in] ref the reference variable. 121 * \return a number congruent to \p alpha modulo 2 pi in the inverval 122 * [ \p ref - M_PI, \p ref + M_PI] 123 */ normalize_periodic_variable(double alpha,double ref)124 double normalize_periodic_variable( 125 double alpha, double ref 126 ) { 127 int count = 0; 128 if(Numeric::is_nan(alpha)) { 129 return 0.0 ; 130 } 131 if(Numeric::is_nan(ref)) { 132 return 0.0 ; 133 } 134 double result = alpha ; 135 count = 0 ; 136 while(ref - result > M_PI) { 137 result += 2.0 * M_PI ; 138 count ++ ; 139 if(count > 100) { 140 return 0.0 ; 141 } 142 } 143 count = 0 ; 144 while(result - ref > M_PI) { 145 result -= 2.0 * M_PI ; 146 count ++ ; 147 if(count > 100) { 148 return 0.0 ; 149 } 150 } 151 return result ; 152 } 153 } 154 155 namespace GEO { 156 157 namespace GlobalParam2d { 158 PGP(Mesh * mesh,Attribute<vec3> & B,Attribute<vec2> & U,double scaling,bool constrain_hard_edges,bool use_direct_solver,double maximum_scaling_correction)159 void PGP( 160 Mesh* mesh, 161 Attribute<vec3>& B, Attribute<vec2>& U, 162 double scaling, bool constrain_hard_edges, 163 bool use_direct_solver, 164 double maximum_scaling_correction 165 ) { 166 geo_cite("DBLP:journals/tog/RayLLSA06"); 167 168 bool do_brush = true; 169 170 // We use f = c/3 171 // (could be fixed in the future, using a c2f array...). 172 geo_assert(mesh->facets.are_simplices()); 173 174 // Step 0: Preparation 175 176 scaling *= 2.0; 177 scaling *= surface_average_edge_length(*mesh); 178 179 // Step 0.1: Preparation / Brushing 180 if(do_brush) { 181 Internal::brush(mesh,B); 182 } 183 184 // Step 0.2: Preparation / Compute relative rotation of B between 185 // pairs of adjacent facets 186 Attribute<index_t> R_ff(mesh->facet_corners.attributes(),"R"); 187 Internal::compute_R_ff(mesh,B,R_ff); 188 Attribute<index_t> R_fv(mesh->facet_corners.attributes(),"R_fv"); 189 Internal::compute_R_fv(mesh,R_ff,R_fv); 190 191 Attribute<double> CC; 192 if(maximum_scaling_correction != 1.0) { 193 Logger::out("PGP") << "Computing scaling correction" 194 << std::endl; 195 CC.bind(mesh->vertices.attributes(), "CC"); 196 Attribute<vec3> Bv(mesh->vertices.attributes(), "B"); 197 Internal::transfer_B_to_vertices(mesh, B, Bv, R_fv); 198 curl_correction( 199 mesh, Bv, R_fv, CC, 200 use_direct_solver, maximum_scaling_correction 201 ); 202 } 203 204 205 Logger::out("PGP") << "Solving for PGP" << std::endl; 206 207 208 // Step 1: Determine the structure of the problem: 209 // - There are four variables per vertex (cu, su, cv, sv) 210 // - Each vertex has a reference corner (v2c[v]) 211 // - Each corner c knows the number of rotations Rc[c] required 212 // to make the B of its triangle match the B of the triangle 213 // of the reference corner attached to its vertex (clear 214 // enough ?) 215 216 // Rot[ R_ff[c] ][i][j] corresponds to the transform to 217 // be applied to the variables associated with a 218 // vertex to express their coordinates in the frame 219 // of corner c (it corresponds to the same Rot[][][] 220 // matrices as QuadCover, with the exception that each 221 // coefficient 1 is replaced with the 2x2 identity matrix 222 // and each coefficient -1 with Mat2x2(1,0,0,-1) (there is 223 // only one "-1" coefficient because cos(-x) = cos(x) !!). 224 // Note: in the PGP paper, the wrong coefficient is negated. 225 226 static const double Rot[4][4][4] = { 227 {{ 1, 0, 0, 0}, 228 { 0, 1, 0, 0}, 229 { 0, 0, 1, 0}, 230 { 0, 0, 0, 1}}, 231 232 {{ 0, 0, 1, 0}, 233 { 0, 0, 0, 1}, 234 { 1, 0, 0, 0}, 235 { 0,-1, 0, 0}}, 236 237 {{ 1, 0, 0, 0}, 238 { 0,-1, 0, 0}, 239 { 0, 0, 1, 0}, 240 { 0, 0, 0,-1}}, 241 242 {{ 0, 0, 1, 0}, 243 { 0, 0, 0,-1}, 244 { 1, 0, 0, 0}, 245 { 0, 1, 0, 0}} 246 }; 247 248 static double Rotc2[4][4]; 249 250 // Step 2: setup and solve linear system 251 { 252 253 nlNewContext(); 254 255 nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE); 256 nlSolverParameteri( 257 NL_NB_VARIABLES, NLint(mesh->vertices.nb()*4) 258 ); 259 260 if(use_direct_solver) { 261 if(nlInitExtension("CHOLMOD")) { 262 nlSolverParameteri(NL_SOLVER, NL_CHOLMOD_EXT); 263 } else if(nlInitExtension("SUPERLU")) { 264 nlSolverParameteri(NL_SOLVER, NL_PERM_SUPERLU_EXT); 265 } else { 266 Logger::warn("PGP") 267 << "Could not initialize direct sovlver" 268 << std::endl; 269 Logger::warn("PGP") 270 << "Falling back to Jacobi pre-CG" 271 << std::endl; 272 use_direct_solver = false; 273 } 274 } 275 276 if(!use_direct_solver) { 277 // With the iterative solver, a very small threshold is 278 // needed, because of the very high scaling on the 279 // solution due to the varying modulus of the complex 280 // numbers far away from the constrained point. 281 nlSolverParameterd(NL_THRESHOLD, 1e-20); 282 } 283 284 nlBegin(NL_SYSTEM); 285 286 if(constrain_hard_edges) { 287 for(index_t c: mesh->facet_corners) { 288 index_t cnstr = 289 Internal::get_edge_constraints(mesh, c, B); 290 index_t v = mesh->facet_corners.vertex(c); 291 // Inverse R, because when the axis turns 292 // clockwise, coordinates turn anticlockwise. 293 index_t Rcc = Internal::inverse_R(R_fv[c]); 294 if(cnstr & Internal::CNSTR_U) { 295 if(Rcc == 0 || Rcc == 2) { 296 nlLockVariable(4*v); 297 nlSetVariable(4*v,1e4); 298 nlLockVariable(4*v+1); 299 nlSetVariable(4*v+1,0.0); 300 } else { 301 nlLockVariable(4*v+2); 302 nlSetVariable(4*v+2,1e4); 303 nlLockVariable(4*v+3); 304 nlSetVariable(4*v+3,0.0); 305 } 306 } 307 if(cnstr & Internal::CNSTR_V) { 308 if(Rcc == 0 || Rcc == 2) { 309 nlLockVariable(4*v+2); 310 nlSetVariable(4*v+2,1e4); 311 nlLockVariable(4*v+3); 312 nlSetVariable(4*v+3,0.0); 313 } else { 314 nlLockVariable(4*v); 315 nlSetVariable(4*v,1e4); 316 nlLockVariable(4*v+1); 317 nlSetVariable(4*v+1,0.0); 318 } 319 } 320 } 321 } 322 323 // Lock at least one variable per connected component. 324 { 325 std::vector<bool> f_visited(mesh->facets.nb(),false); 326 for(index_t f: mesh->facets) { 327 if(!f_visited[f]) { 328 index_t nb_locked=0; 329 index_t first_v = mesh->facets.vertex(f,0); 330 std::stack<index_t> S; 331 f_visited[f] = true; 332 S.push(f); 333 while(!S.empty()) { 334 index_t f_top = S.top(); 335 S.pop(); 336 FOR(le,mesh->facets.nb_vertices(f_top)) { 337 index_t v = mesh->facets.vertex(f_top,le); 338 if( 339 nlVariableIsLocked(4*v) || 340 nlVariableIsLocked(4*v+2)) { 341 ++nb_locked; 342 } 343 index_t f_neigh = 344 mesh->facets.adjacent(f_top,le); 345 346 if(f_neigh != NO_FACET && 347 !f_visited[f_neigh] 348 ) { 349 f_visited[f_neigh] = true; 350 S.push(f_neigh); 351 } 352 } 353 } 354 355 // Lock one of the points in each 356 // connected component to make sure that 357 // the minimum is well defined. 358 if(nb_locked == 0) { 359 nlLockVariable(4*first_v); 360 nlSetVariable(4*first_v,1e4); 361 362 nlLockVariable(4*first_v+1); 363 nlSetVariable(4*first_v+1,0.0); 364 365 nlLockVariable(4*first_v+2); 366 nlSetVariable(4*first_v+2,1e4); 367 368 nlLockVariable(4*first_v+3); 369 nlSetVariable(4*first_v+3,0.0); 370 } 371 } 372 } 373 } 374 375 nlBegin(NL_MATRIX); 376 377 // This one will be replaced in-place 378 // with the rotation that encodes the 379 // delta u and delta v along each edge. 380 381 double RotDelta[4][4]; 382 FOR(i,4) { 383 FOR(j,4) { 384 RotDelta[i][j] = ((i==j) ? 1.0 : 0.0); 385 } 386 } 387 388 for(index_t f: mesh->facets) { 389 390 vec3 Bf, BTf; 391 392 for(index_t c1=mesh->facets.corners_begin(f); 393 c1 < mesh->facets.corners_end(f); ++c1) { 394 395 Internal::get_B_on_edge(mesh, B, R_ff, f, c1, Bf, BTf); 396 397 index_t c2 = 398 mesh->facets.next_corner_around_facet(f,c1); 399 index_t v1 = mesh->facet_corners.vertex(c1); 400 index_t v2 = mesh->facet_corners.vertex(c2); 401 402 vec3 E = vec3(mesh->vertices.point_ptr(v2)) - 403 vec3(mesh->vertices.point_ptr(v1)); 404 double delta_u = 2.0 * M_PI * dot(E,Bf)/scaling; 405 double delta_v = 2.0 * M_PI * dot(E,BTf)/scaling; 406 407 if(CC.is_bound()) { 408 double s = (0.5*(CC[v1] + CC[v2])); 409 delta_u /= s; 410 delta_v /= s; 411 } 412 413 double sdu = sin(delta_u); 414 double cdu = cos(delta_u); 415 double sdv = sin(delta_v); 416 double cdv = cos(delta_v); 417 418 RotDelta[0][0] = cdu; 419 RotDelta[0][1] = -sdu; 420 RotDelta[1][0] = sdu; 421 RotDelta[1][1] = cdu; 422 423 RotDelta[2][2] = cdv; 424 RotDelta[2][3] = -sdv; 425 RotDelta[3][2] = sdv; 426 RotDelta[3][3] = cdv; 427 428 // Inverse R, because when the axis turns 429 // clockwise, coordinates turn anticlockwise. 430 index_t Rc1 = Internal::inverse_R(R_fv[c1]); 431 index_t Rc2 = Internal::inverse_R(R_fv[c2]); 432 433 // Compute the product of the "delta u, delta v" 434 // rotation with the Rc2 "90 degrees rotation" matrix, 435 // exactly like in the PGP article. 436 437 for(index_t i=0; i<4; ++i) { 438 for(index_t j=0; j<4; ++j) { 439 Rotc2[i][j] = 0.0; 440 for(index_t k=0; k<4; ++k) { 441 Rotc2[i][j] += 442 RotDelta[i][k] * Rot[Rc2][k][j]; 443 } 444 } 445 } 446 447 for(index_t i=0; i<4; ++i) { 448 nlBegin(NL_ROW); 449 for(index_t j=0; j<4; ++j) { 450 double a1 = Rot[Rc1][i][j]; 451 if(a1 != 0.0) { 452 nlCoefficient(v1*4+j, a1); 453 } 454 } 455 for(index_t j=0; j<4; ++j) { 456 double a2 = Rotc2[i][j]; 457 if(a2 != 0.0) { 458 nlCoefficient(v2*4+j, -a2); 459 } 460 } 461 nlEnd(NL_ROW); 462 } 463 } 464 } 465 466 nlEnd(NL_MATRIX); 467 nlEnd(NL_SYSTEM); 468 469 nlSolve(); 470 471 Attribute<double> PGP; 472 PGP.bind_if_is_defined(mesh->facet_corners.attributes(),"PGP"); 473 if(!PGP.is_bound()) { 474 PGP.create_vector_attribute( 475 mesh->facet_corners.attributes(), "PGP", 4 476 ); 477 } 478 479 for(index_t c: mesh->facet_corners) { 480 index_t v = mesh->facet_corners.vertex(c); 481 // Inverse R, because when the axis turns 482 // clockwise, coordinates turn anticlockwise. 483 index_t Rcc = Internal::inverse_R(R_fv[c]); 484 double vars[4]; 485 for(index_t i=0; i<4; ++i) { 486 vars[i] = 0.0; 487 for(index_t j=0; j<4; ++j) { 488 vars[i] += Rot[Rcc][i][j] * nlGetVariable(4*v+j); 489 } 490 } 491 double s = sqrt(vars[0]*vars[0]+vars[1]*vars[1]); 492 vars[0] /= s; 493 vars[1] /= s; 494 495 s = sqrt(vars[2]*vars[2]+vars[3]*vars[3]); 496 vars[2] /= s; 497 vars[3] /= s; 498 499 U[c].x = atan2(vars[1], vars[0]); 500 U[c].y = atan2(vars[3], vars[2]); 501 502 FOR(i,4) { 503 PGP[4*c+i] = vars[i]; 504 } 505 } 506 507 Attribute<index_t> singular( 508 mesh->facets.attributes(),"is_singular" 509 ); 510 511 for(index_t f: mesh->facets) { 512 singular[f] = false; 513 for(index_t c1=mesh->facets.corners_begin(f); 514 c1<mesh->facets.corners_end(f); ++c1) { 515 index_t c2 = 516 mesh->facets.next_corner_around_facet(f,c1); 517 518 index_t v1 = mesh->facet_corners.vertex(c1); 519 index_t v2 = mesh->facet_corners.vertex(c2); 520 vec3 E = vec3(mesh->vertices.point_ptr(v2)) - 521 vec3(mesh->vertices.point_ptr(v1)); 522 523 vec3 Bf, BTf; 524 Internal::get_B_on_edge(mesh, B, R_ff, f, c1, Bf, BTf); 525 526 double delta_u = 2.0 * M_PI * dot(E,Bf)/scaling; 527 double delta_v = 2.0 * M_PI * dot(E,BTf)/scaling; 528 529 if(CC.is_bound()) { 530 double s = (0.5*(CC[v1] + CC[v2])); 531 delta_u /= s; 532 delta_v /= s; 533 } 534 535 double expected_u = U[c1].x-delta_u; 536 double expected_v = U[c1].y-delta_v; 537 vec2 uv( 538 normalize_periodic_variable(U[c2].x,expected_u), 539 normalize_periodic_variable(U[c2].y,expected_v) 540 ); 541 if(c2 == mesh->facets.corners_begin(f)) { 542 singular[f] = 543 (length(uv - U[c2]) > 1e-20) ? 1 : 0; 544 } else { 545 U[c2] = uv; 546 } 547 548 } 549 } 550 551 for(index_t c: mesh->facet_corners) { 552 U[c].x /= M_PI; 553 U[c].y /= M_PI; 554 Internal::snap_tex_coord(U[c].x); 555 Internal::snap_tex_coord(U[c].y); 556 } 557 558 nlDeleteContext(nlGetCurrent()); 559 } 560 561 562 } 563 564 curl_correction(Mesh * mesh,Attribute<vec3> & Bv,Attribute<index_t> & R_fv,Attribute<double> & CC,bool use_direct_solver,double max_scaling_correction)565 void curl_correction( 566 Mesh* mesh, Attribute<vec3>& Bv, 567 Attribute<index_t>& R_fv, Attribute<double>& CC, 568 bool use_direct_solver, double max_scaling_correction 569 ) { 570 geo_assert(Bv.manager() == &mesh->vertices.attributes()); 571 nlNewContext(); 572 573 nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE); 574 nlSolverParameteri( 575 NL_NB_VARIABLES, NLint(mesh->vertices.nb()*4) 576 ); 577 578 if(use_direct_solver) { 579 if(nlInitExtension("CHOLMOD")) { 580 nlSolverParameteri(NL_SOLVER, NL_CHOLMOD_EXT); 581 } else if(nlInitExtension("SUPERLU")) { 582 nlSolverParameteri(NL_SOLVER, NL_PERM_SUPERLU_EXT); 583 } else { 584 Logger::warn("PGP") 585 << "Could not initialize direct solver" 586 << std::endl; 587 Logger::warn("PGP") 588 << "Falling back to Jacobi pre-CG" 589 << std::endl; 590 use_direct_solver = false; 591 } 592 } 593 594 if(!use_direct_solver) { 595 // With the iterative solver, a very small threshold is 596 // needed, because of the very high scaling on the 597 // solution due to the varying modulus of the complex 598 // numbers far away from the constrained point. 599 nlSolverParameterd(NL_THRESHOLD, 1e-20); 600 } 601 602 const double locked_value = 1.0; 603 const double solver_scale = 1e3; 604 605 nlBegin(NL_SYSTEM); 606 nlLockVariable(0u); 607 nlSetVariable(0u,locked_value); 608 nlBegin(NL_MATRIX); 609 for(index_t f: mesh->facets) { 610 611 index_t va = mesh->facets.vertex(f,0); 612 index_t vb = mesh->facets.vertex(f,1); 613 index_t vc = mesh->facets.vertex(f,2); 614 615 vec3 N = normalize(Geom::mesh_facet_normal(*mesh,f)); 616 617 vec3 fieldA3d = normalize(Bv[va]); 618 vec3 fieldB3d = normalize(Bv[vb]); 619 vec3 fieldC3d = normalize(Bv[vc]); 620 621 FOR(i, Internal::inverse_R(R_fv[f*3])) { 622 fieldA3d = cross(N,fieldA3d); 623 } 624 625 FOR(i, Internal::inverse_R(R_fv[f*3+1])) { 626 fieldB3d = cross(N,fieldB3d); 627 } 628 629 FOR(i, Internal::inverse_R(R_fv[f*3+2])) { 630 fieldC3d = cross(N,fieldC3d); 631 } 632 633 vec3 U = normalize(cross(cross(N, fieldA3d),N)); 634 vec3 V = normalize(cross(N, U)); 635 636 637 // fieldx is the 2d field at point x. The fields are 638 // rotated in a coherent way to take the modulus into 639 // account notice that it will be better if the field 640 // is rotated to the facet instead of being projected 641 vec2 fieldA(1.0,0.0); 642 vec2 fieldB(dot(fieldB3d, U), dot(fieldB3d, V)) ; 643 fieldB = normalize(fieldB); 644 645 vec2 fieldC (dot(fieldC3d, U), dot(fieldC3d, V)) ; 646 fieldC = normalize(fieldC); 647 648 vec3 A3d(mesh->vertices.point_ptr(va)); 649 vec3 B3d(mesh->vertices.point_ptr(vb)); 650 vec3 C3d(mesh->vertices.point_ptr(vc)); 651 652 vec3 AB3d = B3d-A3d; 653 vec3 AC3d = C3d-A3d; 654 655 vec2 A(0.0,0.0); 656 vec2 B(dot(AB3d,U), dot(AB3d,V)); 657 vec2 C(dot(AC3d,U), dot(AC3d,V)); 658 659 // hummm... flat triangles ? 660 double a=0; 661 double b=0; 662 ParamTrglGradient trg(A,B,C); 663 if (trg.is_flat() || 664 length(cross(AB3d,AC3d))<1e-10 || 665 ::fabs(fieldB.x) < 1e-20 || 666 ::fabs(fieldC.x) < 1e-20 667 ) { 668 std::cerr<<"bad triangle ..." << std::endl ; 669 670 nlBegin(NL_ROW); 671 nlCoefficient(va,1e-2) ; 672 nlCoefficient(vb,-1e-2) ; 673 nlEnd(NL_ROW); 674 675 nlBegin(NL_ROW); 676 nlCoefficient(va,1e-2) ; 677 nlCoefficient(vc,-1e-2) ; 678 nlEnd(NL_ROW); 679 680 nlBegin(NL_ROW); 681 nlCoefficient(vc,1e-2) ; 682 nlCoefficient(vb,-1e-2) ; 683 nlEnd(NL_ROW); 684 continue ; 685 } else { 686 687 // the direction field is represented by angles 688 // double alpha_A = 0; 689 // double alpha_B = atan(fieldB.y/fieldB.x); 690 // double alpha_C = atan(fieldC.y/fieldC.x); 691 692 // Order 1 Taylor Expansion .... 693 double alpha_A = 0; 694 double alpha_B = fieldB.y/fieldB.x; 695 double alpha_C = fieldC.y/fieldC.x; 696 697 vec2 grad_alpha ( 698 trg.TX(0)*alpha_A + 699 trg.TX(1)*alpha_B + 700 trg.TX(2)*alpha_C , 701 trg.TY(0)*alpha_A + 702 trg.TY(1)*alpha_B + 703 trg.TY(2)*alpha_C 704 ); 705 a = grad_alpha.x; 706 b = grad_alpha.y; 707 } 708 709 double sqrt_area = ::sqrt(Geom::mesh_facet_area(*mesh,f)) ; 710 711 if( 712 Numeric::is_nan(sqrt_area) || 713 Numeric::is_nan(trg.TX(0)) || 714 Numeric::is_nan(trg.TX(1)) || 715 Numeric::is_nan(trg.TX(2)) || 716 Numeric::is_nan(trg.TY(0)) || 717 Numeric::is_nan(trg.TY(1)) || 718 Numeric::is_nan(trg.TY(2)) || 719 Numeric::is_nan(a) || 720 Numeric::is_nan(b) 721 ) { 722 std::cerr << "Found NAN !!" << std::endl ; 723 } 724 725 726 nlRowScaling(sqrt_area) ; 727 nlBegin(NL_ROW); 728 nlCoefficient(va, trg.TX(0)) ; 729 nlCoefficient(vb, trg.TX(1)) ; 730 nlCoefficient(vc, trg.TX(2)) ; 731 nlRightHandSide(solver_scale*b); 732 nlEnd(NL_ROW); 733 734 nlRowScaling(sqrt_area) ; 735 nlBegin(NL_ROW); 736 nlCoefficient(va, trg.TY(0)) ; 737 nlCoefficient(vb, trg.TY(1)) ; 738 nlCoefficient(vc, trg.TY(2)) ; 739 nlRightHandSide(solver_scale*-a); 740 nlEnd(NL_ROW); 741 } 742 nlEnd(NL_MATRIX); 743 nlEnd(NL_SYSTEM); 744 745 nlSolve(); 746 double min_val = Numeric::max_float64(); 747 double max_val = Numeric::min_float64(); 748 for(index_t v: mesh->vertices) { 749 CC[v] = ::exp( 750 (nlGetVariable(v)-locked_value)/solver_scale) 751 ; 752 min_val = std::min(min_val,CC[v]); 753 max_val = std::max(max_val,CC[v]); 754 } 755 double avg_val = 0.5 * (min_val + max_val); 756 double min_limit = 1.0 / ::sqrt(max_scaling_correction); 757 double max_limit = ::sqrt(max_scaling_correction); 758 for(index_t v: mesh->vertices) { 759 CC[v] = CC[v] / avg_val; 760 geo_clamp(CC[v], min_limit, max_limit); 761 } 762 nlDeleteContext(nlGetCurrent()); 763 } 764 765 } // namespace GlobalParam2d 766 } // namespace OGF 767