code cleanup: make bmesh operator names more consistant since python has access to...
[blender.git] / source / blender / bmesh / operators / bmo_smooth_laplacian.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * Contributor(s): Alexander Pinzon
19  *
20  * ***** END GPL LICENSE BLOCK *****
21  */
22
23 /** \file blender/bmesh/operators/bmo_smooth_laplacian.c
24  *  \ingroup bmesh
25  */
26
27 #include "MEM_guardedalloc.h"
28
29 #include "DNA_meshdata_types.h"
30
31 #include "BLI_array.h"
32 #include "BLI_heap.h"
33 #include "BLI_math.h"
34 #include "BLI_math_geom.h"
35 #include "BLI_smallhash.h"
36
37 #include "BKE_customdata.h"
38 #include "BKE_mesh.h"
39
40 #include "bmesh.h"
41
42 #include "ONL_opennl.h"
43
44 #include "intern/bmesh_operators_private.h" /* own include */
45
46 // #define SMOOTH_LAPLACIAN_AREA_FACTOR 4.0f  /* UNUSED */
47 // #define SMOOTH_LAPLACIAN_EDGE_FACTOR 2.0f  /* UNUSED */
48 #define SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE 1.8f
49 #define SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE 0.15f
50
51 struct BLaplacianSystem {
52         float *eweights;        /* Length weights per Edge */
53         float (*fweights)[3];   /* Cotangent weights per face */
54         float *ring_areas;      /* Total area per ring*/
55         float *vlengths;        /* Total sum of lengths(edges) per vertice*/
56         float *vweights;        /* Total sum of weights per vertice*/
57         int numEdges;           /* Number of edges*/
58         int numFaces;           /* Number of faces*/
59         int numVerts;           /* Number of verts*/
60         short *zerola;          /* Is zero area or length*/
61
62         /* Pointers to data*/
63         BMesh *bm;
64         BMOperator *op;
65         NLContext *context;
66
67         /*Data*/
68         float min_area;
69 };
70 typedef struct BLaplacianSystem LaplacianSystem;
71
72 static float compute_volume(BMesh *bm);
73 static float cotan_weight(float *v1, float *v2, float *v3);
74 static int vert_is_boundary(BMVert *v);
75 static LaplacianSystem *init_laplacian_system(int a_numEdges, int a_numFaces, int a_numVerts);
76 static void init_laplacian_matrix(LaplacianSystem *sys);
77 static void delete_laplacian_system(LaplacianSystem *sys);
78 static void delete_void_pointer(void *data);
79 static void fill_laplacian_matrix(LaplacianSystem *sys);
80 static void memset_laplacian_system(LaplacianSystem *sys, int val);
81 static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez, int preserve_volume);
82 static void volume_preservation(BMOperator *op, float vini, float vend, int usex, int usey, int usez);
83
84 static void delete_void_pointer(void *data)
85 {
86         if (data) {
87                 MEM_freeN(data);
88                 data = NULL;
89         }
90 }
91
92 static void delete_laplacian_system(LaplacianSystem *sys)
93 {
94         delete_void_pointer(sys->eweights);
95         delete_void_pointer(sys->fweights);
96         delete_void_pointer(sys->ring_areas);
97         delete_void_pointer(sys->vlengths);
98         delete_void_pointer(sys->vweights);
99         delete_void_pointer(sys->zerola);
100         if (sys->context) {
101                 nlDeleteContext(sys->context);
102         }
103         sys->bm = NULL;
104         sys->op = NULL;
105         MEM_freeN(sys);
106 }
107
108 static void memset_laplacian_system(LaplacianSystem *sys, int val)
109 {
110         memset(sys->eweights,     val, sizeof(float) * sys->numEdges);
111         memset(sys->fweights,     val, sizeof(float) * sys->numFaces * 3);
112         memset(sys->ring_areas,   val, sizeof(float) * sys->numVerts);
113         memset(sys->vlengths,     val, sizeof(float) * sys->numVerts);
114         memset(sys->vweights,     val, sizeof(float) * sys->numVerts);
115         memset(sys->zerola,       val, sizeof(short) * sys->numVerts);
116 }
117
118 static LaplacianSystem *init_laplacian_system(int a_numEdges, int a_numFaces, int a_numVerts)
119 {
120         LaplacianSystem *sys;
121         sys = MEM_callocN(sizeof(LaplacianSystem), "ModLaplSmoothSystem");
122         sys->numEdges = a_numEdges;
123         sys->numFaces = a_numFaces;
124         sys->numVerts = a_numVerts;
125
126         sys->eweights =  MEM_callocN(sizeof(float) * sys->numEdges, "ModLaplSmoothEWeight");
127         if (!sys->eweights) {
128                 delete_laplacian_system(sys);
129                 return NULL;
130         }
131
132         sys->fweights =  MEM_callocN(sizeof(float) * 3 * sys->numFaces, "ModLaplSmoothFWeight");
133         if (!sys->fweights) {
134                 delete_laplacian_system(sys);
135                 return NULL;
136         }
137
138         sys->ring_areas =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothRingAreas");
139         if (!sys->ring_areas) {
140                 delete_laplacian_system(sys);
141                 return NULL;
142         }
143
144         sys->vlengths =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVlengths");
145         if (!sys->vlengths) {
146                 delete_laplacian_system(sys);
147                 return NULL;
148         }
149
150         sys->vweights =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVweights");
151         if (!sys->vweights) {
152                 delete_laplacian_system(sys);
153                 return NULL;
154         }
155
156         sys->zerola =  MEM_callocN(sizeof(short) * sys->numVerts, "ModLaplSmoothZeloa");
157         if (!sys->zerola) {
158                 delete_laplacian_system(sys);
159                 return NULL;
160         }
161
162         return sys;
163 }
164
165 /* Compute weigth between vertice v_i and all your neighbors
166  * weight between v_i and v_neighbor
167  * Wij = cot(alpha) + cot(beta) / (4.0 * total area of all faces  * sum all weight)
168  *        v_i *
169  *          / | \
170  *         /  |  \
171  *  v_beta*   |   * v_alpha
172  *         \  |  /
173  *          \ | /
174  *            * v_neighbor
175  */
176
177 static void init_laplacian_matrix(LaplacianSystem *sys)
178 {
179         float areaf;
180         float *v1, *v2, *v3, *v4;
181         float w1, w2, w3, w4;
182         int i, j;
183         int has_4_vert;
184         unsigned int idv1, idv2, idv3, idv4, idv[4];
185         BMEdge *e;
186         BMFace *f;
187         BMIter eiter;
188         BMIter fiter;
189         BMIter vi;
190         BMVert *vn;
191         BMVert *vf[4];
192
193         BM_ITER_MESH_INDEX (e, &eiter, sys->bm, BM_EDGES_OF_MESH, j) {
194                 if (!BM_elem_flag_test(e, BM_ELEM_SELECT) && BM_edge_is_boundary(e)) {
195                         v1 = e->v1->co;
196                         v2 =  e->v2->co;
197                         idv1 = BM_elem_index_get(e->v1);
198                         idv2 = BM_elem_index_get(e->v2);
199
200                         w1 = len_v3v3(v1, v2);
201                         if (w1 > sys->min_area) {
202                                 w1 = 1.0f / w1;
203                                 i = BM_elem_index_get(e);
204                                 sys->eweights[i] = w1;
205                                 sys->vlengths[idv1] += w1;
206                                 sys->vlengths[idv2] += w1;
207                         }
208                         else {
209                                 sys->zerola[idv1] = 1;
210                                 sys->zerola[idv2] = 1;
211                         }
212                 }
213         }
214
215         BM_ITER_MESH (f, &fiter, sys->bm, BM_FACES_OF_MESH) {
216                 if (BM_elem_flag_test(f, BM_ELEM_SELECT)) {
217
218                         BM_ITER_ELEM_INDEX (vn, &vi, f, BM_VERTS_OF_FACE, i) {
219                                 vf[i] = vn;
220                         }
221                         has_4_vert = (i == 4) ? 1 : 0;
222                         idv1 = BM_elem_index_get(vf[0]);
223                         idv2 = BM_elem_index_get(vf[1]);
224                         idv3 = BM_elem_index_get(vf[2]);
225                         idv4 = has_4_vert ? BM_elem_index_get(vf[3]) : 0;
226
227                         v1 = vf[0]->co;
228                         v2 = vf[1]->co;
229                         v3 = vf[2]->co;
230                         v4 = has_4_vert ? vf[3]->co : 0;
231
232                         if (has_4_vert) {
233                                 areaf = area_quad_v3(v1, v2, v3, v4);
234                         }
235                         else {
236                                 areaf = area_tri_v3(v1, v2, v3);
237                         }
238
239                         if (fabsf(areaf) < sys->min_area) {
240                                 sys->zerola[idv1] = 1;
241                                 sys->zerola[idv2] = 1;
242                                 sys->zerola[idv3] = 1;
243                                 if (has_4_vert) sys->zerola[idv4] = 1;
244                         }
245
246                         sys->ring_areas[idv1] += areaf;
247                         sys->ring_areas[idv2] += areaf;
248                         sys->ring_areas[idv3] += areaf;
249                         if (has_4_vert) sys->ring_areas[idv4] += areaf;
250
251                         if (has_4_vert) {
252
253                                 idv[0] = idv1;
254                                 idv[1] = idv2;
255                                 idv[2] = idv3;
256                                 idv[3] = idv4;
257
258                                 for (j = 0; j < 4; j++) {
259                                         idv1 = idv[j];
260                                         idv2 = idv[(j + 1) % 4];
261                                         idv3 = idv[(j + 2) % 4];
262                                         idv4 = idv[(j + 3) % 4];
263
264                                         v1 = vf[j]->co;
265                                         v2 = vf[(j + 1) % 4]->co;
266                                         v3 = vf[(j + 2) % 4]->co;
267                                         v4 = vf[(j + 3) % 4]->co;
268
269                                         w2 = cotan_weight(v4, v1, v2) + cotan_weight(v3, v1, v2);
270                                         w3 = cotan_weight(v2, v3, v1) + cotan_weight(v4, v1, v3);
271                                         w4 = cotan_weight(v2, v4, v1) + cotan_weight(v3, v4, v1);
272
273                                         sys->vweights[idv1] += (w2 + w3 + w4) / 4.0f;
274                                 }
275                         }
276                         else {
277                                 i = BM_elem_index_get(f);
278
279                                 w1 = cotan_weight(v1, v2, v3);
280                                 w2 = cotan_weight(v2, v3, v1);
281                                 w3 = cotan_weight(v3, v1, v2);
282
283                                 sys->fweights[i][0] += w1;
284                                 sys->fweights[i][1] += w2;
285                                 sys->fweights[i][2] += w3;
286
287                                 sys->vweights[idv1] += w2 + w3;
288                                 sys->vweights[idv2] += w1 + w3;
289                                 sys->vweights[idv3] += w1 + w2;
290                         }
291                 }
292         }
293 }
294
295 static void fill_laplacian_matrix(LaplacianSystem *sys)
296 {
297         float *v1, *v2, *v3, *v4;
298         float w2, w3, w4;
299         int i, j;
300         int has_4_vert;
301         unsigned int idv1, idv2, idv3, idv4, idv[4];
302
303         BMEdge *e;
304         BMFace *f;
305         BMIter eiter;
306         BMIter fiter;
307         BMIter vi;
308         BMVert *vn;
309         BMVert *vf[4];
310
311         BM_ITER_MESH (f, &fiter, sys->bm, BM_FACES_OF_MESH) {
312                 if (BM_elem_flag_test(f, BM_ELEM_SELECT)) {
313                         BM_ITER_ELEM_INDEX (vn, &vi, f, BM_VERTS_OF_FACE, i) {
314                                 vf[i] = vn;
315                         }
316                         has_4_vert = (i == 4) ? 1 : 0;
317                         if (has_4_vert) {
318                                 idv[0] = BM_elem_index_get(vf[0]);
319                                 idv[1] = BM_elem_index_get(vf[1]);
320                                 idv[2] = BM_elem_index_get(vf[2]);
321                                 idv[3] = BM_elem_index_get(vf[3]);
322                                 for (j = 0; j < 4; j++) {
323                                         idv1 = idv[j];
324                                         idv2 = idv[(j + 1) % 4];
325                                         idv3 = idv[(j + 2) % 4];
326                                         idv4 = idv[(j + 3) % 4];
327
328                                         v1 = vf[j]->co;
329                                         v2 = vf[(j + 1) % 4]->co;
330                                         v3 = vf[(j + 2) % 4]->co;
331                                         v4 = vf[(j + 3) % 4]->co;
332
333                                         w2 = cotan_weight(v4, v1, v2) + cotan_weight(v3, v1, v2);
334                                         w3 = cotan_weight(v2, v3, v1) + cotan_weight(v4, v1, v3);
335                                         w4 = cotan_weight(v2, v4, v1) + cotan_weight(v3, v4, v1);
336
337                                         w2 = w2 / 4.0f;
338                                         w3 = w3 / 4.0f;
339                                         w4 = w4 / 4.0f;
340
341                                         if (!vert_is_boundary(vf[j]) && sys->zerola[idv1] == 0) {
342                                                 nlMatrixAdd(idv1, idv2, w2 * sys->vweights[idv1]);
343                                                 nlMatrixAdd(idv1, idv3, w3 * sys->vweights[idv1]);
344                                                 nlMatrixAdd(idv1, idv4, w4 * sys->vweights[idv1]);
345                                         }
346                                 }
347                         }
348                         else {
349                                 idv1 = BM_elem_index_get(vf[0]);
350                                 idv2 = BM_elem_index_get(vf[1]);
351                                 idv3 = BM_elem_index_get(vf[2]);
352                                 /* Is ring if number of faces == number of edges around vertice*/
353                                 i = BM_elem_index_get(f);
354                                 if (!vert_is_boundary(vf[0]) && sys->zerola[idv1] == 0) {
355                                         nlMatrixAdd(idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]);
356                                         nlMatrixAdd(idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]);
357                                 }
358                                 if (!vert_is_boundary(vf[1]) && sys->zerola[idv2] == 0) {
359                                         nlMatrixAdd(idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]);
360                                         nlMatrixAdd(idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]);
361                                 }
362                                 if (!vert_is_boundary(vf[2]) && sys->zerola[idv3] == 0) {
363                                         nlMatrixAdd(idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]);
364                                         nlMatrixAdd(idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]);
365                                 }
366                         }
367                 }
368         }
369         BM_ITER_MESH (e, &eiter, sys->bm, BM_EDGES_OF_MESH) {
370                 if (!BM_elem_flag_test(e, BM_ELEM_SELECT) && BM_edge_is_boundary(e) ) {
371                         v1 = e->v1->co;
372                         v2 =  e->v2->co;
373                         idv1 = BM_elem_index_get(e->v1);
374                         idv2 = BM_elem_index_get(e->v2);
375                         if (sys->zerola[idv1] == 0 && sys->zerola[idv2] == 0) {
376                                 i = BM_elem_index_get(e);
377                                 nlMatrixAdd(idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
378                                 nlMatrixAdd(idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
379                         }
380                 }
381         }
382 }
383
384 static float cotan_weight(float *v1, float *v2, float *v3)
385 {
386         float a[3], b[3], c[3], clen;
387
388         sub_v3_v3v3(a, v2, v1);
389         sub_v3_v3v3(b, v3, v1);
390         cross_v3_v3v3(c, a, b);
391
392         clen = len_v3(c);
393
394         if (clen == 0.0f)
395                 return 0.0f;
396
397         return dot_v3v3(a, b) / clen;
398 }
399
400 static int vert_is_boundary(BMVert *v)
401 {
402         BMEdge *ed;
403         BMFace *f;
404         BMIter ei;
405         BMIter fi;
406         BM_ITER_ELEM (ed, &ei, v, BM_EDGES_OF_VERT) {
407                 if (BM_edge_is_boundary(ed)) {
408                         return 1;
409                 }
410         }
411         BM_ITER_ELEM (f, &fi, v, BM_FACES_OF_VERT) {
412                 if (!BM_elem_flag_test(f, BM_ELEM_SELECT)) {
413                         return 1;
414                 }
415         }
416         return 0;
417 }
418
419 static float compute_volume(BMesh *bm)
420 {
421         float vol = 0.0f;
422         float x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4;
423         int i;
424         BMFace *f;
425         BMIter fiter;
426         BMIter vi;
427         BMVert *vn;
428         BMVert *vf[4];
429
430         BM_ITER_MESH (f, &fiter, bm, BM_FACES_OF_MESH) {
431                 BM_ITER_ELEM_INDEX (vn, &vi, f, BM_VERTS_OF_FACE, i) {
432                         vf[i] = vn;
433                 }
434                 x1 = vf[0]->co[0];
435                 y1 = vf[0]->co[1];
436                 z1 = vf[0]->co[2];
437
438                 x2 = vf[1]->co[0];
439                 y2 = vf[1]->co[1];
440                 z2 = vf[1]->co[2];
441
442                 x3 = vf[2]->co[0];
443                 y3 = vf[2]->co[1];
444                 z3 = vf[2]->co[2];
445
446                 vol += (1.0f / 6.0f) * (0.0f - x3 * y2 * z1 + x2 * y3 * z1 + x3 * y1 * z2 - x1 * y3 * z2 - x2 * y1 * z3 + x1 * y2 * z3);
447
448                 if (i == 4) {
449                         x4 = vf[3]->co[0];
450                         y4 = vf[3]->co[1];
451                         z4 = vf[3]->co[2];
452                         vol += (1.0f / 6.0f) * (x1 * y3 * z4 - x1 * y4 * z3 - x3 * y1 * z4 + x3 * z1 * y4 + y1 * x4 * z3 - x4 * y3 * z1);
453                 }
454         }
455         return fabs(vol);
456 }
457
458 static void volume_preservation(BMOperator *op, float vini, float vend, int usex, int usey, int usez)
459 {
460         float beta;
461         BMOIter siter;
462         BMVert *v;
463
464         if (vend != 0.0f) {
465                 beta  = pow(vini / vend, 1.0f / 3.0f);
466                 BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) {
467                         if (usex) {
468                                 v->co[0] *= beta;
469                         }
470                         if (usey) {
471                                 v->co[1] *= beta;
472                         }
473                         if (usez) {
474                                 v->co[2] *= beta;
475                         }
476
477                 }
478         }
479 }
480
481 static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez, int preserve_volume)
482 {
483         int m_vertex_id;
484         float leni, lene;
485         float vini, vend;
486         float *vi1, *vi2, ve1[3], ve2[3];
487         unsigned int idv1, idv2;
488         BMOIter siter;
489         BMVert *v;
490         BMEdge *e;
491         BMIter eiter;
492
493         BM_ITER_MESH  (e, &eiter, sys->bm, BM_EDGES_OF_MESH) {
494                 idv1 = BM_elem_index_get(e->v1);
495                 idv2 = BM_elem_index_get(e->v2);
496                 vi1 = e->v1->co;
497                 vi2 =  e->v2->co;
498                 ve1[0] = nlGetVariable(0, idv1);
499                 ve1[1] = nlGetVariable(1, idv1);
500                 ve1[2] = nlGetVariable(2, idv1);
501                 ve2[0] = nlGetVariable(0, idv2);
502                 ve2[1] = nlGetVariable(1, idv2);
503                 ve2[2] = nlGetVariable(2, idv2);
504                 leni = len_v3v3(vi1, vi2);
505                 lene = len_v3v3(ve1, ve2);
506                 if (lene > leni * SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE || lene < leni * SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE) {
507                         sys->zerola[idv1] = 1;
508                         sys->zerola[idv2] = 1;
509                 }
510         }
511
512         if (preserve_volume) {
513                 vini = compute_volume(sys->bm);
514         }
515         BMO_ITER (v, &siter, sys->op->slots_in, "verts", BM_VERT) {
516                 m_vertex_id = BM_elem_index_get(v);
517                 if (sys->zerola[m_vertex_id] == 0) {
518                         if (usex) {
519                                 v->co[0] =  nlGetVariable(0, m_vertex_id);
520                         }
521                         if (usey) {
522                                 v->co[1] =  nlGetVariable(1, m_vertex_id);
523                         }
524                         if (usez) {
525                                 v->co[2] =  nlGetVariable(2, m_vertex_id);
526                         }
527                 }
528         }
529         if (preserve_volume) {
530                 vend = compute_volume(sys->bm);
531                 volume_preservation(sys->op, vini, vend, usex, usey, usez);
532         }
533
534 }
535
536 void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op)
537 {
538         int i;
539         int m_vertex_id;
540         int usex, usey, usez, preserve_volume;
541         float lambda, lambda_border;
542         float w;
543         BMOIter siter;
544         BMVert *v;
545         LaplacianSystem *sys;
546
547         sys = init_laplacian_system(bm->totedge, bm->totface, bm->totvert);
548         if (!sys) return;
549         sys->bm = bm;
550         sys->op = op;
551
552         memset_laplacian_system(sys, 0);
553
554         BM_mesh_elem_index_ensure(bm, BM_VERT);
555         lambda = BMO_slot_float_get(op->slots_in, "lambda");
556         lambda_border = BMO_slot_float_get(op->slots_in, "lambda_border");
557         sys->min_area = 0.00001f;
558         usex = BMO_slot_bool_get(op->slots_in, "use_x");
559         usey = BMO_slot_bool_get(op->slots_in, "use_y");
560         usez = BMO_slot_bool_get(op->slots_in, "use_z");
561         preserve_volume = BMO_slot_bool_get(op->slots_in, "preserve_volume");
562
563
564         nlNewContext();
565         sys->context = nlGetCurrent();
566
567         nlSolverParameteri(NL_NB_VARIABLES, bm->totvert);
568         nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
569         nlSolverParameteri(NL_NB_ROWS, bm->totvert);
570         nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 3);
571
572         nlBegin(NL_SYSTEM);
573         for (i = 0; i < bm->totvert; i++) {
574                 nlLockVariable(i);
575         }
576         BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) {
577                 m_vertex_id = BM_elem_index_get(v);
578                 nlUnlockVariable(m_vertex_id);
579                 nlSetVariable(0, m_vertex_id, v->co[0]);
580                 nlSetVariable(1, m_vertex_id, v->co[1]);
581                 nlSetVariable(2, m_vertex_id, v->co[2]);
582         }
583
584         nlBegin(NL_MATRIX);
585         init_laplacian_matrix(sys);
586         BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) {
587                 m_vertex_id = BM_elem_index_get(v);
588                 nlRightHandSideAdd(0, m_vertex_id, v->co[0]);
589                 nlRightHandSideAdd(1, m_vertex_id, v->co[1]);
590                 nlRightHandSideAdd(2, m_vertex_id, v->co[2]);
591                 i = m_vertex_id;
592                 if (sys->zerola[i] == 0) {
593                         w = sys->vweights[i] * sys->ring_areas[i];
594                         sys->vweights[i] = (w == 0.0f) ? 0.0f : -lambda  / (4.0f * w);
595                         w = sys->vlengths[i];
596                         sys->vlengths[i] = (w == 0.0f) ? 0.0f : -lambda_border  * 2.0f / w;
597
598                         if (!vert_is_boundary(v)) {
599                                 nlMatrixAdd(i, i,  1.0f + lambda / (4.0f * sys->ring_areas[i]));
600                         }
601                         else {
602                                 nlMatrixAdd(i, i,  1.0f + lambda_border * 2.0f);
603                         }
604                 }
605                 else {
606                         nlMatrixAdd(i, i, 1.0f);
607                 }
608         }
609         fill_laplacian_matrix(sys);
610
611         nlEnd(NL_MATRIX);
612         nlEnd(NL_SYSTEM);
613
614         if (nlSolveAdvanced(NULL, NL_TRUE) ) {
615                 validate_solution(sys, usex, usey, usez, preserve_volume);
616         }
617
618         delete_laplacian_system(sys);
619 }