OpenNL: make the API thread safe by always passing context.

Previously two laplacian smooth or deform modifiers executing
simultaneously could crash.

7 files changed, 282 insertions, 308 deletions

diff --git a/source/blender/bmesh/operators/bmo_smooth_laplacian.c b/source/blender/bmesh/operators/bmo_smooth_laplacian.c
index 661b1c21892..624fbd93818 100644
--- a/source/blender/bmesh/operators/bmo_smooth_laplacian.c
+++ b/source/blender/bmesh/operators/bmo_smooth_laplacian.c
@@ -333,9 +333,9 @@ static void fill_laplacian_matrix(LaplacianSystem *sys)
 					w4 = w4 / 4.0f;
 
 					if (!vert_is_boundary(vf[j]) && sys->zerola[idv1] == 0) {
-						nlMatrixAdd(idv1, idv2, w2 * sys->vweights[idv1]);
-						nlMatrixAdd(idv1, idv3, w3 * sys->vweights[idv1]);
-						nlMatrixAdd(idv1, idv4, w4 * sys->vweights[idv1]);
+						nlMatrixAdd(sys->context, idv1, idv2, w2 * sys->vweights[idv1]);
+						nlMatrixAdd(sys->context, idv1, idv3, w3 * sys->vweights[idv1]);
+						nlMatrixAdd(sys->context, idv1, idv4, w4 * sys->vweights[idv1]);
 					}
 				}
 			}
@@ -346,16 +346,16 @@ static void fill_laplacian_matrix(LaplacianSystem *sys)
 				/* Is ring if number of faces == number of edges around vertice*/
 				i = BM_elem_index_get(f);
 				if (!vert_is_boundary(vf[0]) && sys->zerola[idv1] == 0) {
-					nlMatrixAdd(idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]);
-					nlMatrixAdd(idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]);
+					nlMatrixAdd(sys->context, idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]);
+					nlMatrixAdd(sys->context, idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]);
 				}
 				if (!vert_is_boundary(vf[1]) && sys->zerola[idv2] == 0) {
-					nlMatrixAdd(idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]);
-					nlMatrixAdd(idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]);
+					nlMatrixAdd(sys->context, idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]);
+					nlMatrixAdd(sys->context, idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]);
 				}
 				if (!vert_is_boundary(vf[2]) && sys->zerola[idv3] == 0) {
-					nlMatrixAdd(idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]);
-					nlMatrixAdd(idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]);
+					nlMatrixAdd(sys->context, idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]);
+					nlMatrixAdd(sys->context, idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]);
 				}
 			}
 		}
@@ -368,8 +368,8 @@ static void fill_laplacian_matrix(LaplacianSystem *sys)
 			idv2 = BM_elem_index_get(e->v2);
 			if (sys->zerola[idv1] == 0 && sys->zerola[idv2] == 0) {
 				i = BM_elem_index_get(e);
-				nlMatrixAdd(idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
-				nlMatrixAdd(idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
+				nlMatrixAdd(sys->context, idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
+				nlMatrixAdd(sys->context, idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
 			}
 		}
 	}
@@ -434,12 +434,12 @@ static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez
 		idv2 = BM_elem_index_get(e->v2);
 		vi1 = e->v1->co;
 		vi2 =  e->v2->co;
-		ve1[0] = nlGetVariable(0, idv1);
-		ve1[1] = nlGetVariable(1, idv1);
-		ve1[2] = nlGetVariable(2, idv1);
-		ve2[0] = nlGetVariable(0, idv2);
-		ve2[1] = nlGetVariable(1, idv2);
-		ve2[2] = nlGetVariable(2, idv2);
+		ve1[0] = nlGetVariable(sys->context, 0, idv1);
+		ve1[1] = nlGetVariable(sys->context, 1, idv1);
+		ve1[2] = nlGetVariable(sys->context, 2, idv1);
+		ve2[0] = nlGetVariable(sys->context, 0, idv2);
+		ve2[1] = nlGetVariable(sys->context, 1, idv2);
+		ve2[2] = nlGetVariable(sys->context, 2, idv2);
 		leni = len_v3v3(vi1, vi2);
 		lene = len_v3v3(ve1, ve2);
 		if (lene > leni * SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE || lene < leni * SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE) {
@@ -455,13 +455,13 @@ static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez
 		m_vertex_id = BM_elem_index_get(v);
 		if (sys->zerola[m_vertex_id] == 0) {
 			if (usex) {
-				v->co[0] =  nlGetVariable(0, m_vertex_id);
+				v->co[0] =  nlGetVariable(sys->context, 0, m_vertex_id);
 			}
 			if (usey) {
-				v->co[1] =  nlGetVariable(1, m_vertex_id);
+				v->co[1] =  nlGetVariable(sys->context, 1, m_vertex_id);
 			}
 			if (usez) {
-				v->co[2] =  nlGetVariable(2, m_vertex_id);
+				v->co[2] =  nlGetVariable(sys->context, 2, m_vertex_id);
 			}
 		}
 	}
@@ -501,33 +501,32 @@ void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op)
 	preserve_volume = BMO_slot_bool_get(op->slots_in, "preserve_volume");
 
 
-	nlNewContext();
-	sys->context = nlGetCurrent();
+	sys->context = nlNewContext();
 
-	nlSolverParameteri(NL_NB_VARIABLES, bm->totvert);
-	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
-	nlSolverParameteri(NL_NB_ROWS, bm->totvert);
-	nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 3);
+	nlSolverParameteri(sys->context, NL_NB_VARIABLES, bm->totvert);
+	nlSolverParameteri(sys->context, NL_LEAST_SQUARES, NL_TRUE);
+	nlSolverParameteri(sys->context, NL_NB_ROWS, bm->totvert);
+	nlSolverParameteri(sys->context, NL_NB_RIGHT_HAND_SIDES, 3);
 
-	nlBegin(NL_SYSTEM);
+	nlBegin(sys->context, NL_SYSTEM);
 	for (i = 0; i < bm->totvert; i++) {
-		nlLockVariable(i);
+		nlLockVariable(sys->context, i);
 	}
 	BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) {
 		m_vertex_id = BM_elem_index_get(v);
-		nlUnlockVariable(m_vertex_id);
-		nlSetVariable(0, m_vertex_id, v->co[0]);
-		nlSetVariable(1, m_vertex_id, v->co[1]);
-		nlSetVariable(2, m_vertex_id, v->co[2]);
+		nlUnlockVariable(sys->context, m_vertex_id);
+		nlSetVariable(sys->context, 0, m_vertex_id, v->co[0]);
+		nlSetVariable(sys->context, 1, m_vertex_id, v->co[1]);
+		nlSetVariable(sys->context, 2, m_vertex_id, v->co[2]);
 	}
 
-	nlBegin(NL_MATRIX);
+	nlBegin(sys->context, NL_MATRIX);
 	init_laplacian_matrix(sys);
 	BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) {
 		m_vertex_id = BM_elem_index_get(v);
-		nlRightHandSideAdd(0, m_vertex_id, v->co[0]);
-		nlRightHandSideAdd(1, m_vertex_id, v->co[1]);
-		nlRightHandSideAdd(2, m_vertex_id, v->co[2]);
+		nlRightHandSideAdd(sys->context, 0, m_vertex_id, v->co[0]);
+		nlRightHandSideAdd(sys->context, 1, m_vertex_id, v->co[1]);
+		nlRightHandSideAdd(sys->context, 2, m_vertex_id, v->co[2]);
 		i = m_vertex_id;
 		if (sys->zerola[i] == 0) {
 			w = sys->vweights[i] * sys->ring_areas[i];
@@ -536,22 +535,22 @@ void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op)
 			sys->vlengths[i] = (w == 0.0f) ? 0.0f : -lambda_border  * 2.0f / w;
 
 			if (!vert_is_boundary(v)) {
-				nlMatrixAdd(i, i,  1.0f + lambda_factor / (4.0f * sys->ring_areas[i]));
+				nlMatrixAdd(sys->context, i, i,  1.0f + lambda_factor / (4.0f * sys->ring_areas[i]));
 			}
 			else {
-				nlMatrixAdd(i, i,  1.0f + lambda_border * 2.0f);
+				nlMatrixAdd(sys->context, i, i,  1.0f + lambda_border * 2.0f);
 			}
 		}
 		else {
-			nlMatrixAdd(i, i, 1.0f);
+			nlMatrixAdd(sys->context, i, i, 1.0f);
 		}
 	}
 	fill_laplacian_matrix(sys);
 
-	nlEnd(NL_MATRIX);
-	nlEnd(NL_SYSTEM);
+	nlEnd(sys->context, NL_MATRIX);
+	nlEnd(sys->context, NL_SYSTEM);
 
-	if (nlSolve(NL_TRUE) ) {
+	if (nlSolve(sys->context, NL_TRUE) ) {
 		validate_solution(sys, usex, usey, usez, preserve_volume);
 	}
 
diff --git a/source/blender/editors/armature/meshlaplacian.c b/source/blender/editors/armature/meshlaplacian.c
index 956cbe25c3b..3f27a58930d 100644
--- a/source/blender/editors/armature/meshlaplacian.c
+++ b/source/blender/editors/armature/meshlaplacian.c
@@ -182,18 +182,18 @@ static void laplacian_triangle_weights(LaplacianSystem *sys, int f, int i1, int
 	t2 = cotangent_tri_weight_v3(v2, v3, v1) / laplacian_edge_count(sys->edgehash, i3, i1);
 	t3 = cotangent_tri_weight_v3(v3, v1, v2) / laplacian_edge_count(sys->edgehash, i1, i2);
 
-	nlMatrixAdd(i1, i1, (t2 + t3) * varea[i1]);
-	nlMatrixAdd(i2, i2, (t1 + t3) * varea[i2]);
-	nlMatrixAdd(i3, i3, (t1 + t2) * varea[i3]);
+	nlMatrixAdd(sys->context, i1, i1, (t2 + t3) * varea[i1]);
+	nlMatrixAdd(sys->context, i2, i2, (t1 + t3) * varea[i2]);
+	nlMatrixAdd(sys->context, i3, i3, (t1 + t2) * varea[i3]);
 
-	nlMatrixAdd(i1, i2, -t3 * varea[i1]);
-	nlMatrixAdd(i2, i1, -t3 * varea[i2]);
+	nlMatrixAdd(sys->context, i1, i2, -t3 * varea[i1]);
+	nlMatrixAdd(sys->context, i2, i1, -t3 * varea[i2]);
 
-	nlMatrixAdd(i2, i3, -t1 * varea[i2]);
-	nlMatrixAdd(i3, i2, -t1 * varea[i3]);
+	nlMatrixAdd(sys->context, i2, i3, -t1 * varea[i2]);
+	nlMatrixAdd(sys->context, i3, i2, -t1 * varea[i3]);
 
-	nlMatrixAdd(i3, i1, -t2 * varea[i3]);
-	nlMatrixAdd(i1, i3, -t2 * varea[i1]);
+	nlMatrixAdd(sys->context, i3, i1, -t2 * varea[i3]);
+	nlMatrixAdd(sys->context, i1, i3, -t2 * varea[i1]);
 
 	if (sys->storeweights) {
 		sys->fweights[f][0] = t1 * varea[i1];
@@ -219,12 +219,10 @@ static LaplacianSystem *laplacian_system_construct_begin(int totvert, int totfac
 	sys->storeweights = 0;
 
 	/* create opennl context */
-	nlNewContext();
-	nlSolverParameteri(NL_NB_VARIABLES, totvert);
+	sys->context = nlNewContext();
+	nlSolverParameteri(sys->context, NL_NB_VARIABLES, totvert);
 	if (lsq)
-		nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
-
-	sys->context = nlGetCurrent();
+		nlSolverParameteri(sys->context, NL_LEAST_SQUARES, NL_TRUE);
 
 	return sys;
 }
@@ -274,7 +272,7 @@ static void laplacian_system_construct_end(LaplacianSystem *sys)
 
 		/* for heat weighting */
 		if (sys->heat.H)
-			nlMatrixAdd(a, a, sys->heat.H[a]);
+			nlMatrixAdd(sys->context, a, a, sys->heat.H[a]);
 	}
 
 	if (sys->storeweights)
@@ -312,41 +310,41 @@ void laplacian_begin_solve(LaplacianSystem *sys, int index)
 	int a;
 
 	if (!sys->nlbegun) {
-		nlBegin(NL_SYSTEM);
+		nlBegin(sys->context, NL_SYSTEM);
 
 		if (index >= 0) {
 			for (a = 0; a < sys->totvert; a++) {
 				if (sys->vpinned[a]) {
-					nlSetVariable(0, a, sys->verts[a][index]);
-					nlLockVariable(a);
+					nlSetVariable(sys->context, 0, a, sys->verts[a][index]);
+					nlLockVariable(sys->context, a);
 				}
 			}
 		}
 
-		nlBegin(NL_MATRIX);
+		nlBegin(sys->context, NL_MATRIX);
 		sys->nlbegun = 1;
 	}
 }
 
-void laplacian_add_right_hand_side(LaplacianSystem *UNUSED(sys), int v, float value)
+void laplacian_add_right_hand_side(LaplacianSystem *sys, int v, float value)
 {
-	nlRightHandSideAdd(0, v, value);
+	nlRightHandSideAdd(sys->context, 0, v, value);
 }
 
 int laplacian_system_solve(LaplacianSystem *sys)
 {
-	nlEnd(NL_MATRIX);
-	nlEnd(NL_SYSTEM);
+	nlEnd(sys->context, NL_MATRIX);
+	nlEnd(sys->context, NL_SYSTEM);
 	sys->nlbegun = 0;
 
-	//nlPrintMatrix();
+	//nlPrintMatrix(sys->context, );
 
-	return nlSolve(NL_TRUE);
+	return nlSolve(sys->context, NL_TRUE);
 }
 
-float laplacian_system_get_solution(int v)
+float laplacian_system_get_solution(LaplacianSystem *sys, int v)
 {
-	return nlGetVariable(0, v);
+	return nlGetVariable(sys->context, 0, v);
 }
 
 /************************* Heat Bone Weighting ******************************/
@@ -723,7 +721,7 @@ void heat_bone_weighting(Object *ob, Mesh *me, float (*verts)[3], int numsource,
 				if (mask && !mask[a])
 					continue;
 
-				solution = laplacian_system_get_solution(a);
+				solution = laplacian_system_get_solution(sys, a);
 				
 				if (bbone) {
 					if (solution > 0.0f)
@@ -1286,7 +1284,7 @@ static float meshdeform_boundary_total_weight(MeshDeformBind *mdb, int x, int y,
 	return totweight;
 }
 
-static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, int x, int y, int z)
+static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, NLContext *context, int x, int y, int z)
 {
 	MDefBoundIsect *isect;
 	float weight, totweight;
@@ -1296,7 +1294,7 @@ static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, int x, int y, int z)
 	if (mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR)
 		return;
 
-	nlMatrixAdd(mdb->varidx[acenter], mdb->varidx[acenter], 1.0f);
+	nlMatrixAdd(context, mdb->varidx[acenter], mdb->varidx[acenter], 1.0f);
 	
 	totweight = meshdeform_boundary_total_weight(mdb, x, y, z);
 	for (i = 1; i <= 6; i++) {
@@ -1307,12 +1305,12 @@ static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, int x, int y, int z)
 		isect = mdb->boundisect[acenter][i - 1];
 		if (!isect) {
 			weight = (1.0f / mdb->width[0]) / totweight;
-			nlMatrixAdd(mdb->varidx[acenter], mdb->varidx[a], -weight);
+			nlMatrixAdd(context, mdb->varidx[acenter], mdb->varidx[a], -weight);
 		}
 	}
 }
 
-static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, int x, int y, int z, int cagevert)
+static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, NLContext *context, int x, int y, int z, int cagevert)
 {
 	MDefBoundIsect *isect;
 	float rhs, weight, totweight;
@@ -1333,7 +1331,7 @@ static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, int x, int y, int z,
 		if (isect) {
 			weight = (1.0f / isect->len) / totweight;
 			rhs = weight * meshdeform_boundary_phi(mdb, isect, cagevert);
-			nlRightHandSideAdd(0, mdb->varidx[acenter], rhs);
+			nlRightHandSideAdd(context, 0, mdb->varidx[acenter], rhs);
 		}
 	}
 }
@@ -1406,39 +1404,38 @@ static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind
 	progress_bar(0, "Starting mesh deform solve");
 
 	/* setup opennl solver */
-	nlNewContext();
-	context = nlGetCurrent();
+	context = nlNewContext();
 
-	nlSolverParameteri(NL_NB_VARIABLES, totvar);
-	nlSolverParameteri(NL_NB_ROWS, totvar);
-	nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 1);
+	nlSolverParameteri(context, NL_NB_VARIABLES, totvar);
+	nlSolverParameteri(context, NL_NB_ROWS, totvar);
+	nlSolverParameteri(context, NL_NB_RIGHT_HAND_SIDES, 1);
 
-	nlBegin(NL_SYSTEM);
-	nlBegin(NL_MATRIX);
+	nlBegin(context, NL_SYSTEM);
+	nlBegin(context, NL_MATRIX);
 
 	/* build matrix */
 	for (z = 0; z < mdb->size; z++)
 		for (y = 0; y < mdb->size; y++)
 			for (x = 0; x < mdb->size; x++)
-				meshdeform_matrix_add_cell(mdb, x, y, z);
+				meshdeform_matrix_add_cell(mdb, context, x, y, z);
 
 	/* solve for each cage vert */
 	for (a = 0; a < mdb->totcagevert; a++) {
 		if (a != 0) {
-			nlBegin(NL_SYSTEM);
-			nlBegin(NL_MATRIX);
+			nlBegin(context, NL_SYSTEM);
+			nlBegin(context, NL_MATRIX);
 		}
 
 		/* fill in right hand side and solve */
 		for (z = 0; z < mdb->size; z++)
 			for (y = 0; y < mdb->size; y++)
 				for (x = 0; x < mdb->size; x++)
-					meshdeform_matrix_add_rhs(mdb, x, y, z, a);
+					meshdeform_matrix_add_rhs(mdb, context, x, y, z, a);
 
-		nlEnd(NL_MATRIX);
-		nlEnd(NL_SYSTEM);
+		nlEnd(context, NL_MATRIX);
+		nlEnd(context, NL_SYSTEM);
 
-		if (nlSolve(NL_TRUE)) {
+		if (nlSolve(context, NL_TRUE)) {
 			for (z = 0; z < mdb->size; z++)
 				for (y = 0; y < mdb->size; y++)
 					for (x = 0; x < mdb->size; x++)
@@ -1451,7 +1448,7 @@ static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind
 
 			for (b = 0; b < mdb->size3; b++) {
 				if (mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR)
-					mdb->phi[b] = nlGetVariable(0, mdb->varidx[b]);
+					mdb->phi[b] = nlGetVariable(context, 0, mdb->varidx[b]);
 				mdb->totalphi[b] += mdb->phi[b];
 			}
 
diff --git a/source/blender/editors/armature/meshlaplacian.h b/source/blender/editors/armature/meshlaplacian.h
index 1412136c1a8..bba8c739abf 100644
--- a/source/blender/editors/armature/meshlaplacian.h
+++ b/source/blender/editors/armature/meshlaplacian.h
@@ -48,7 +48,7 @@ void laplacian_add_triangle(LaplacianSystem *sys, int v1, int v2, int v3);
 void laplacian_begin_solve(LaplacianSystem *sys, int index);
 void laplacian_add_right_hand_side(LaplacianSystem *sys, int v, float value);
 int laplacian_system_solve(LaplacianSystem *sys);
-float laplacian_system_get_solution(int v);
+float laplacian_system_get_solution(LaplacianSystem *sys, int v);
 
 /* Heat Weighting */
 
diff --git a/source/blender/editors/armature/reeb.c b/source/blender/editors/armature/reeb.c
index e10ae4bcf3f..d53e5b857a0 100644
--- a/source/blender/editors/armature/reeb.c
+++ b/source/blender/editors/armature/reeb.c
@@ -2497,7 +2497,7 @@ int weightFromLoc(EditMesh *em, int axis)
 	return 1;
 }
 
-static void addTriangle(EditVert *v1, EditVert *v2, EditVert *v3, int e1, int e2, int e3)
+static void addTriangle(NLContext *context, EditVert *v1, EditVert *v2, EditVert *v3, int e1, int e2, int e3)
 {
 	/* Angle opposite e1 */
 	float t1 = cotangent_tri_weight_v3(v1->co, v2->co, v3->co) / e2;
@@ -2512,22 +2512,23 @@ static void addTriangle(EditVert *v1, EditVert *v2, EditVert *v3, int e1, int e2
 	int i2 = indexData(v2);
 	int i3 = indexData(v3);
 	
-	nlMatrixAdd(i1, i1, t2 + t3);
-	nlMatrixAdd(i2, i2, t1 + t3);
-	nlMatrixAdd(i3, i3, t1 + t2);
+	nlMatrixAdd(context, i1, i1, t2 + t3);
+	nlMatrixAdd(context, i2, i2, t1 + t3);
+	nlMatrixAdd(context, i3, i3, t1 + t2);
 
-	nlMatrixAdd(i1, i2, -t3);
-	nlMatrixAdd(i2, i1, -t3);
+	nlMatrixAdd(context, i1, i2, -t3);
+	nlMatrixAdd(context, i2, i1, -t3);
 
-	nlMatrixAdd(i2, i3, -t1);
-	nlMatrixAdd(i3, i2, -t1);
+	nlMatrixAdd(context, i2, i3, -t1);
+	nlMatrixAdd(context, i3, i2, -t1);
 
-	nlMatrixAdd(i3, i1, -t2);
-	nlMatrixAdd(i1, i3, -t2);
+	nlMatrixAdd(context, i3, i1, -t2);
+	nlMatrixAdd(context, i1, i3, -t2);
 }
 
 int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
 {
+	NLContext *context;
 	NLboolean success;
 	EditVert *eve;
 	EditEdge *eed;
@@ -2543,11 +2544,11 @@ int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
 
 	/* Solve with openNL */
 	
-	nlNewContext();
+	context = nlNewContext();
 
-	nlSolverParameteri(NL_NB_VARIABLES, totvert);
+	nlSolverParameteri(context, NL_NB_VARIABLES, totvert);
 
-	nlBegin(NL_SYSTEM);
+	nlBegin(context, NL_SYSTEM);
 	
 	/* Find local extrema */
 	for (index = 0, eve = em->verts.first; eve; index++, eve = eve->next) {
@@ -2582,8 +2583,8 @@ int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
 			if (maximum || minimum) {
 				float w = weightData(eve);
 				eve->f1 = 0;
-				nlSetVariable(0, index, w);
-				nlLockVariable(index);
+				nlSetVariable(context, 0, index, w);
+				nlLockVariable(context, index);
 			}
 			else {
 				eve->f1 = 1;
@@ -2591,7 +2592,7 @@ int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
 		}
 	}
 	
-	nlBegin(NL_MATRIX);
+	nlBegin(context, NL_MATRIX);
 
 	/* Zero edge weight */
 	for (eed = em->edges.first; eed; eed = eed->next) {
@@ -2615,32 +2616,32 @@ int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
 	for (efa = em->faces.first; efa; efa = efa->next) {
 		if (efa->h == 0) {
 			if (efa->v4 == NULL) {
-				addTriangle(efa->v1, efa->v2, efa->v3, efa->e1->tmp.l, efa->e2->tmp.l, efa->e3->tmp.l);
+				addTriangle(context, efa->v1, efa->v2, efa->v3, efa->e1->tmp.l, efa->e2->tmp.l, efa->e3->tmp.l);
 			}
 			else {
-				addTriangle(efa->v1, efa->v2, efa->v3, efa->e1->tmp.l, efa->e2->tmp.l, 2);
-				addTriangle(efa->v3, efa->v4, efa->v1, efa->e3->tmp.l, efa->e4->tmp.l, 2);
+				addTriangle(context, efa->v1, efa->v2, efa->v3, efa->e1->tmp.l, efa->e2->tmp.l, 2);
+				addTriangle(context, efa->v3, efa->v4, efa->v1, efa->e3->tmp.l, efa->e4->tmp.l, 2);
 			}
 		}
 	}
 	
-	nlEnd(NL_MATRIX);
+	nlEnd(context, NL_MATRIX);
 
-	nlEnd(NL_SYSTEM);
+	nlEnd(context, NL_SYSTEM);
 
-	success = nlSolve(NL_TRUE);
+	success = nlSolve(context, NL_TRUE);
 
 	if (success) {
 		rval = 1;
 		for (index = 0, eve = em->verts.first; eve; index++, eve = eve->next) {
-			weightSetData(eve, nlGetVariable(0, index));
+			weightSetData(eve, nlGetVariable(context, 0, index));
 		}
 	}
 	else {
 		rval = 0;
 	}
 
-	nlDeleteContext(nlGetCurrent());
+	nlDeleteContext(context);
 
 	return rval;
 }
diff --git a/source/blender/editors/uvedit/uvedit_parametrizer.c b/source/blender/editors/uvedit/uvedit_parametrizer.c
index cb387b4af22..bb233919089 100644
--- a/source/blender/editors/uvedit/uvedit_parametrizer.c
+++ b/source/blender/editors/uvedit/uvedit_parametrizer.c
@@ -2471,16 +2471,17 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 	PEdge *e;
 	int i, j, ninterior = sys->ninterior, nvar = 2 * sys->ninterior;
 	PBool success;
+	NLContext *context;
 
-	nlNewContext();
-	nlSolverParameteri(NL_NB_VARIABLES, nvar);
+	context = nlNewContext();
+	nlSolverParameteri(context, NL_NB_VARIABLES, nvar);
 
-	nlBegin(NL_SYSTEM);
+	nlBegin(context, NL_SYSTEM);
 
-	nlBegin(NL_MATRIX);
+	nlBegin(context, NL_MATRIX);
 
 	for (i = 0; i < nvar; i++)
-		nlRightHandSideAdd(0, i, sys->bInterior[i]);
+		nlRightHandSideAdd(context, 0, i, sys->bInterior[i]);
 
 	for (f = chart->faces; f; f = f->nextlink) {
 		float wi1, wi2, wi3, b, si, beta[3], j2[3][3], W[3][3];
@@ -2526,8 +2527,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 			sys->J2dt[e2->u.id][0] = j2[1][0] = p_abf_compute_sin_product(sys, v1, e2->u.id) * wi2;
 			sys->J2dt[e3->u.id][0] = j2[2][0] = p_abf_compute_sin_product(sys, v1, e3->u.id) * wi3;
 
-			nlRightHandSideAdd(0, v1->u.id, j2[0][0] * beta[0]);
-			nlRightHandSideAdd(0, ninterior + v1->u.id, j2[1][0] * beta[1] + j2[2][0] * beta[2]);
+			nlRightHandSideAdd(context, 0, v1->u.id, j2[0][0] * beta[0]);
+			nlRightHandSideAdd(context, 0, ninterior + v1->u.id, j2[1][0] * beta[1] + j2[2][0] * beta[2]);
 
 			row1[0] = j2[0][0] * W[0][0];
 			row2[0] = j2[0][0] * W[1][0];
@@ -2546,8 +2547,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 			sys->J2dt[e2->u.id][1] = j2[1][1] = 1.0f * wi2;
 			sys->J2dt[e3->u.id][1] = j2[2][1] = p_abf_compute_sin_product(sys, v2, e3->u.id) * wi3;
 
-			nlRightHandSideAdd(0, v2->u.id, j2[1][1] * beta[1]);
-			nlRightHandSideAdd(0, ninterior + v2->u.id, j2[0][1] * beta[0] + j2[2][1] * beta[2]);
+			nlRightHandSideAdd(context, 0, v2->u.id, j2[1][1] * beta[1]);
+			nlRightHandSideAdd(context, 0, ninterior + v2->u.id, j2[0][1] * beta[0] + j2[2][1] * beta[2]);
 
 			row1[1] = j2[1][1] * W[0][1];
 			row2[1] = j2[1][1] * W[1][1];
@@ -2566,8 +2567,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 			sys->J2dt[e2->u.id][2] = j2[1][2] = p_abf_compute_sin_product(sys, v3, e2->u.id) * wi2;
 			sys->J2dt[e3->u.id][2] = j2[2][2] = 1.0f * wi3;
 
-			nlRightHandSideAdd(0, v3->u.id, j2[2][2] * beta[2]);
-			nlRightHandSideAdd(0, ninterior + v3->u.id, j2[0][2] * beta[0] + j2[1][2] * beta[1]);
+			nlRightHandSideAdd(context, 0, v3->u.id, j2[2][2] * beta[2]);
+			nlRightHandSideAdd(context, 0, ninterior + v3->u.id, j2[0][2] * beta[0] + j2[1][2] * beta[1]);
 
 			row1[2] = j2[2][2] * W[0][2];
 			row2[2] = j2[2][2] * W[1][2];
@@ -2591,29 +2592,29 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 					continue;
 
 				if (i == 0)
-					nlMatrixAdd(r, c, j2[0][i] * row1[j]);
+					nlMatrixAdd(context, r, c, j2[0][i] * row1[j]);
 				else
-					nlMatrixAdd(r + ninterior, c, j2[0][i] * row1[j]);
+					nlMatrixAdd(context, r + ninterior, c, j2[0][i] * row1[j]);
 
 				if (i == 1)
-					nlMatrixAdd(r, c, j2[1][i] * row2[j]);
+					nlMatrixAdd(context, r, c, j2[1][i] * row2[j]);
 				else
-					nlMatrixAdd(r + ninterior, c, j2[1][i] * row2[j]);
+					nlMatrixAdd(context, r + ninterior, c, j2[1][i] * row2[j]);
 
 
 				if (i == 2)
-					nlMatrixAdd(r, c, j2[2][i] * row3[j]);
+					nlMatrixAdd(context, r, c, j2[2][i] * row3[j]);
 				else
-					nlMatrixAdd(r + ninterior, c, j2[2][i] * row3[j]);
+					nlMatrixAdd(context, r + ninterior, c, j2[2][i] * row3[j]);
 			}
 		}
 	}
 
-	nlEnd(NL_MATRIX);
+	nlEnd(context, NL_MATRIX);
 
-	nlEnd(NL_SYSTEM);
+	nlEnd(context, NL_SYSTEM);
 
-	success = nlSolve(NL_FALSE);
+	success = nlSolve(context, NL_FALSE);
 
 	if (success) {
 		for (f = chart->faces; f; f = f->nextlink) {
@@ -2624,24 +2625,24 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 			pre[0] = pre[1] = pre[2] = 0.0;
 
 			if (v1->flag & PVERT_INTERIOR) {
-				float x = nlGetVariable(0, v1->u.id);
-				float x2 = nlGetVariable(0, ninterior + v1->u.id);
+				float x = nlGetVariable(context, 0, v1->u.id);
+				float x2 = nlGetVariable(context, 0, ninterior + v1->u.id);
 				pre[0] += sys->J2dt[e1->u.id][0] * x;
 				pre[1] += sys->J2dt[e2->u.id][0] * x2;
 				pre[2] += sys->J2dt[e3->u.id][0] * x2;
 			}
 
 			if (v2->flag & PVERT_INTERIOR) {
-				float x = nlGetVariable(0, v2->u.id);
-				float x2 = nlGetVariable(0, ninterior + v2->u.id);
+				float x = nlGetVariable(context, 0, v2->u.id);
+				float x2 = nlGetVariable(context, 0, ninterior + v2->u.id);
 				pre[0] += sys->J2dt[e1->u.id][1] * x2;
 				pre[1] += sys->J2dt[e2->u.id][1] * x;
 				pre[2] += sys->J2dt[e3->u.id][1] * x2;
 			}
 
 			if (v3->flag & PVERT_INTERIOR) {
-				float x = nlGetVariable(0, v3->u.id);
-				float x2 = nlGetVariable(0, ninterior + v3->u.id);
+				float x = nlGetVariable(context, 0, v3->u.id);
+				float x2 = nlGetVariable(context, 0, ninterior + v3->u.id);
 				pre[0] += sys->J2dt[e1->u.id][2] * x2;
 				pre[1] += sys->J2dt[e2->u.id][2] * x2;
 				pre[2] += sys->J2dt[e3->u.id][2] * x;
@@ -2672,12 +2673,12 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 		}
 
 		for (i = 0; i < ninterior; i++) {
-			sys->lambdaPlanar[i] += (float)nlGetVariable(0, i);
-			sys->lambdaLength[i] += (float)nlGetVariable(0, ninterior + i);
+			sys->lambdaPlanar[i] += (float)nlGetVariable(context, 0, i);
+			sys->lambdaLength[i] += (float)nlGetVariable(context, 0, ninterior + i);
 		}
 	}
 
-	nlDeleteContext(nlGetCurrent());
+	nlDeleteContext(context);
 
 	return success;
 }
@@ -3003,11 +3004,12 @@ static void p_chart_extrema_verts(PChart *chart, PVert **pin1, PVert **pin2)
 
 static void p_chart_lscm_load_solution(PChart *chart)
 {
+	NLContext *context = chart->u.lscm.context;
 	PVert *v;
 
 	for (v = chart->verts; v; v = v->nextlink) {
-		v->uv[0] = nlGetVariable(0, 2 * v->u.id);
-		v->uv[1] = nlGetVariable(0, 2 * v->u.id + 1);
+		v->uv[0] = nlGetVariable(context, 0, 2 * v->u.id);
+		v->uv[1] = nlGetVariable(context, 0, 2 * v->u.id + 1);
 	}
 }
 
@@ -3062,17 +3064,16 @@ static void p_chart_lscm_begin(PChart *chart, PBool live, PBool abf)
 		for (v = chart->verts; v; v = v->nextlink)
 			v->u.id = id++;
 
-		nlNewContext();
-		nlSolverParameteri(NL_NB_VARIABLES, 2 * chart->nverts);
-		nlSolverParameteri(NL_NB_ROWS, 2 * chart->nfaces);
-		nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
-
-		chart->u.lscm.context = nlGetCurrent();
+		chart->u.lscm.context = nlNewContext();
+		nlSolverParameteri(chart->u.lscm.context, NL_NB_VARIABLES, 2 * chart->nverts);
+		nlSolverParameteri(chart->u.lscm.context, NL_NB_ROWS, 2 * chart->nfaces);
+		nlSolverParameteri(chart->u.lscm.context, NL_LEAST_SQUARES, NL_TRUE);
 	}
 }
 
 static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 {
+	NLContext *context = chart->u.lscm.context;
 	PVert *v, *pin1 = chart->u.lscm.pin1, *pin2 = chart->u.lscm.pin2;
 	PFace *f;
 	float *alpha = chart->u.lscm.abf_alpha;
@@ -3080,9 +3081,7 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 	bool flip_faces;
 	int row;
 
-	nlMakeCurrent(chart->u.lscm.context);
-
-	nlBegin(NL_SYSTEM);
+	nlBegin(context, NL_SYSTEM);
 
 #if 0
 	/* TODO: make loading pins work for simplify/complexify. */
@@ -3093,25 +3092,25 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 			p_vert_load_pin_select_uvs(handle, v);  /* reload for live */
 
 	if (chart->u.lscm.pin1) {
-		nlLockVariable(2 * pin1->u.id);
-		nlLockVariable(2 * pin1->u.id + 1);
-		nlLockVariable(2 * pin2->u.id);
-		nlLockVariable(2 * pin2->u.id + 1);
+		nlLockVariable(context, 2 * pin1->u.id);
+		nlLockVariable(context, 2 * pin1->u.id + 1);
+		nlLockVariable(context, 2 * pin2->u.id);
+		nlLockVariable(context, 2 * pin2->u.id + 1);
 
-		nlSetVariable(0, 2 * pin1->u.id, pin1->uv[0]);
-		nlSetVariable(0, 2 * pin1->u.id + 1, pin1->uv[1]);
-		nlSetVariable(0, 2 * pin2->u.id, pin2->uv[0]);
-		nlSetVariable(0, 2 * pin2->u.id + 1, pin2->uv[1]);
+		nlSetVariable(context, 0, 2 * pin1->u.id, pin1->uv[0]);
+		nlSetVariable(context, 0, 2 * pin1->u.id + 1, pin1->uv[1]);
+		nlSetVariable(context, 0, 2 * pin2->u.id, pin2->uv[0]);
+		nlSetVariable(context, 0, 2 * pin2->u.id + 1, pin2->uv[1]);
 	}
 	else {
 		/* set and lock the pins */
 		for (v = chart->verts; v; v = v->nextlink) {
 			if (v->flag & PVERT_PIN) {
-				nlLockVariable(2 * v->u.id);
-				nlLockVariable(2 * v->u.id + 1);
+				nlLockVariable(context, 2 * v->u.id);
+				nlLockVariable(context, 2 * v->u.id + 1);
 
-				nlSetVariable(0, 2 * v->u.id, v->uv[0]);
-				nlSetVariable(0, 2 * v->u.id + 1, v->uv[1]);
+				nlSetVariable(context, 0, 2 * v->u.id, v->uv[0]);
+				nlSetVariable(context, 0, 2 * v->u.id + 1, v->uv[1]);
 			}
 		}
 	}
@@ -3138,7 +3137,7 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 
 	/* construct matrix */
 
-	nlBegin(NL_MATRIX);
+	nlBegin(context, NL_MATRIX);
 
 	row = 0;
 	for (f = chart->faces; f; f = f->nextlink) {
@@ -3186,44 +3185,26 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 		cosine = cosf(a1) * ratio;
 		sine = sina1 * ratio;
 
-#if 0
-		nlBegin(NL_ROW);
-		nlCoefficient(2 * v1->u.id,   cosine - 1.0);
-		nlCoefficient(2 * v1->u.id + 1, -sine);
-		nlCoefficient(2 * v2->u.id,   -cosine);
-		nlCoefficient(2 * v2->u.id + 1, sine);
-		nlCoefficient(2 * v3->u.id,   1.0);
-		nlEnd(NL_ROW);
-
-		nlBegin(NL_ROW);
-		nlCoefficient(2 * v1->u.id,   sine);
-		nlCoefficient(2 * v1->u.id + 1, cosine - 1.0);
-		nlCoefficient(2 * v2->u.id,   -sine);
-		nlCoefficient(2 * v2->u.id + 1, -cosine);
-		nlCoefficient(2 * v3->u.id + 1, 1.0);
-		nlEnd(NL_ROW);
-#else
-		nlMatrixAdd(row, 2 * v1->u.id,   cosine - 1.0f);
-		nlMatrixAdd(row, 2 * v1->u.id + 1, -sine);
-		nlMatrixAdd(row, 2 * v2->u.id,   -cosine);
-		nlMatrixAdd(row, 2 * v2->u.id + 1, sine);
-		nlMatrixAdd(row, 2 * v3->u.id,   1.0);
+		nlMatrixAdd(context, row, 2 * v1->u.id,   cosine - 1.0f);
+		nlMatrixAdd(context, row, 2 * v1->u.id + 1, -sine);
+		nlMatrixAdd(context, row, 2 * v2->u.id,   -cosine);
+		nlMatrixAdd(context, row, 2 * v2->u.id + 1, sine);
+		nlMatrixAdd(context, row, 2 * v3->u.id,   1.0);
 		row++;
 
-		nlMatrixAdd(row, 2 * v1->u.id,   sine);
-		nlMatrixAdd(row, 2 * v1->u.id + 1, cosine - 1.0f);
-		nlMatrixAdd(row, 2 * v2->u.id,   -sine);
-		nlMatrixAdd(row, 2 * v2->u.id + 1, -cosine);
-		nlMatrixAdd(row, 2 * v3->u.id + 1, 1.0);
+		nlMatrixAdd(context, row, 2 * v1->u.id,   sine);
+		nlMatrixAdd(context, row, 2 * v1->u.id + 1, cosine - 1.0f);
+		nlMatrixAdd(context, row, 2 * v2->u.id,   -sine);
+		nlMatrixAdd(context, row, 2 * v2->u.id + 1, -cosine);
+		nlMatrixAdd(context, row, 2 * v3->u.id + 1, 1.0);
 		row++;
-#endif
 	}
 
-	nlEnd(NL_MATRIX);
+	nlEnd(context, NL_MATRIX);
 
-	nlEnd(NL_SYSTEM);
+	nlEnd(context, NL_SYSTEM);
 
-	if (nlSolve(NL_TRUE)) {
+	if (nlSolve(context, NL_TRUE)) {
 		p_chart_lscm_load_solution(chart);
 		return P_TRUE;
 	}
diff --git a/source/blender/modifiers/intern/MOD_laplaciandeform.c b/source/blender/modifiers/intern/MOD_laplaciandeform.c
index 916f8b8a36d..fdaacc7cd9e 100644
--- a/source/blender/modifiers/intern/MOD_laplaciandeform.c
+++ b/source/blender/modifiers/intern/MOD_laplaciandeform.c
@@ -283,9 +283,9 @@ static void initLaplacianMatrix(LaplacianSystem *sys)
 			sys->delta[idv[0]][1] -= v3[1] * w3;
 			sys->delta[idv[0]][2] -= v3[2] * w3;
 
-			nlMatrixAdd(idv[0], idv[1], -w2);
-			nlMatrixAdd(idv[0], idv[2], -w3);
-			nlMatrixAdd(idv[0], idv[0], w2 + w3);
+			nlMatrixAdd(sys->context, idv[0], idv[1], -w2);
+			nlMatrixAdd(sys->context, idv[0], idv[2], -w3);
+			nlMatrixAdd(sys->context, idv[0], idv[0], w2 + w3);
 		}
 	}
 }
@@ -338,9 +338,9 @@ static void rotateDifferentialCoordinates(LaplacianSystem *sys)
 		beta = dot_v3v3(uij, di);
 		gamma = dot_v3v3(e2, di);
 
-		pi[0] = nlGetVariable(0, i);
-		pi[1] = nlGetVariable(1, i);
-		pi[2] = nlGetVariable(2, i);
+		pi[0] = nlGetVariable(sys->context, 0, i);
+		pi[1] = nlGetVariable(sys->context, 1, i);
+		pi[2] = nlGetVariable(sys->context, 2, i);
 		zero_v3(ni);
 		num_fni = 0;
 		num_fni = sys->ringf_map[i].count;
@@ -349,9 +349,9 @@ static void rotateDifferentialCoordinates(LaplacianSystem *sys)
 			fidn = sys->ringf_map[i].indices;
 			vin = sys->tris[fidn[fi]];
 			for (j = 0; j < 3; j++) {
-				vn[j][0] = nlGetVariable(0, vin[j]);
-				vn[j][1] = nlGetVariable(1, vin[j]);
-				vn[j][2] = nlGetVariable(2, vin[j]);
+				vn[j][0] = nlGetVariable(sys->context, 0, vin[j]);
+				vn[j][1] = nlGetVariable(sys->context, 1, vin[j]);
+				vn[j][2] = nlGetVariable(sys->context, 2, vin[j]);
 				if (vin[j] == sys->unit_verts[i]) {
 					copy_v3_v3(pj, vn[j]);
 				}
@@ -372,14 +372,14 @@ static void rotateDifferentialCoordinates(LaplacianSystem *sys)
 		fni[2] = alpha * ni[2] + beta * uij[2] + gamma * e2[2];
 
 		if (len_squared_v3(fni) > FLT_EPSILON) {
-			nlRightHandSideSet(0, i, fni[0]);
-			nlRightHandSideSet(1, i, fni[1]);
-			nlRightHandSideSet(2, i, fni[2]);
+			nlRightHandSideSet(sys->context, 0, i, fni[0]);
+			nlRightHandSideSet(sys->context, 1, i, fni[1]);
+			nlRightHandSideSet(sys->context, 2, i, fni[2]);
 		}
 		else {
-			nlRightHandSideSet(0, i, sys->delta[i][0]);
-			nlRightHandSideSet(1, i, sys->delta[i][1]);
-			nlRightHandSideSet(2, i, sys->delta[i][2]);
+			nlRightHandSideSet(sys->context, 0, i, sys->delta[i][0]);
+			nlRightHandSideSet(sys->context, 1, i, sys->delta[i][1]);
+			nlRightHandSideSet(sys->context, 2, i, sys->delta[i][2]);
 		}
 	}
 }
@@ -395,71 +395,70 @@ static void laplacianDeformPreview(LaplacianSystem *sys, float (*vertexCos)[3])
 #endif
 
 	if (!sys->is_matrix_computed) {
-		nlNewContext();
-		sys->context = nlGetCurrent();
-
-		nlSolverParameteri(NL_NB_VARIABLES, n);
-		nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
-		nlSolverParameteri(NL_NB_ROWS, n + na);
-		nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 3);
-		nlBegin(NL_SYSTEM);
+		sys->context = nlNewContext();
+
+		nlSolverParameteri(sys->context, NL_NB_VARIABLES, n);
+		nlSolverParameteri(sys->context, NL_LEAST_SQUARES, NL_TRUE);
+		nlSolverParameteri(sys->context, NL_NB_ROWS, n + na);
+		nlSolverParameteri(sys->context, NL_NB_RIGHT_HAND_SIDES, 3);
+		nlBegin(sys->context, NL_SYSTEM);
 		for (i = 0; i < n; i++) {
-			nlSetVariable(0, i, sys->co[i][0]);
-			nlSetVariable(1, i, sys->co[i][1]);
-			nlSetVariable(2, i, sys->co[i][2]);
+			nlSetVariable(sys->context, 0, i, sys->co[i][0]);
+			nlSetVariable(sys->context, 1, i, sys->co[i][1]);
+			nlSetVariable(sys->context, 2, i, sys->co[i][2]);
 		}
 		for (i = 0; i < na; i++) {
 			vid = sys->index_anchors[i];
-			nlSetVariable(0, vid, vertexCos[vid][0]);
-			nlSetVariable(1, vid, vertexCos[vid][1]);
-			nlSetVariable(2, vid, vertexCos[vid][2]);
+			nlSetVariable(sys->context, 0, vid, vertexCos[vid][0]);
+			nlSetVariable(sys->context, 1, vid, vertexCos[vid][1]);
+			nlSetVariable(sys->context, 2, vid, vertexCos[vid][2]);
 		}
-		nlBegin(NL_MATRIX);
+		nlBegin(sys->context, NL_MATRIX);
 
 		initLaplacianMatrix(sys);
 		computeImplictRotations(sys);
 
 		for (i = 0; i < n; i++) {
-			nlRightHandSideSet(0, i, sys->delta[i][0]);
-			nlRightHandSideSet(1, i, sys->delta[i][1]);
-			nlRightHandSideSet(2, i, sys->delta[i][2]);
+			nlRightHandSideSet(sys->context, 0, i, sys->delta[i][0]);
+			nlRightHandSideSet(sys->context, 1, i, sys->delta[i][1]);
+			nlRightHandSideSet(sys->context, 2, i, sys->delta[i][2]);
 		}
 		for (i = 0; i < na; i++) {
 			vid = sys->index_anchors[i];
-			nlRightHandSideSet(0, n + i, vertexCos[vid][0]);
-			nlRightHandSideSet(1, n + i, vertexCos[vid][1]);
-			nlRightHandSideSet(2, n + i, vertexCos[vid][2]);
-			nlMatrixAdd(n + i, vid, 1.0f);
+			nlRightHandSideSet(sys->context, 0, n + i, vertexCos[vid][0]);
+			nlRightHandSideSet(sys->context, 1, n + i, vertexCos[vid][1]);
+			nlRightHandSideSet(sys->context, 2, n + i, vertexCos[vid][2]);
+			nlMatrixAdd(sys->context, n + i, vid, 1.0f);
 		}
-		nlEnd(NL_MATRIX);
-		nlEnd(NL_SYSTEM);
-		if (nlSolve(NL_TRUE)) {
+		nlEnd(sys->context, NL_MATRIX);
+		nlEnd(sys->context, NL_SYSTEM);
+		if (nlSolve(sys->context, NL_TRUE)) {
 			sys->has_solution = true;
 
 			for (j = 1; j <= sys->repeat; j++) {
-				nlBegin(NL_SYSTEM);
-				nlBegin(NL_MATRIX);
+				nlBegin(sys->context, NL_SYSTEM);
+				nlBegin(sys->context, NL_MATRIX);
 				rotateDifferentialCoordinates(sys);
 
 				for (i = 0; i < na; i++) {
 					vid = sys->index_anchors[i];
-					nlRightHandSideSet(0, n + i, vertexCos[vid][0]);
-					nlRightHandSideSet(1, n + i, vertexCos[vid][1]);
-					nlRightHandSideSet(2, n + i, vertexCos[vid][2]);
+					nlRightHandSideSet(sys->context, 0, n + i, vertexCos[vid][0]);
+					nlRightHandSideSet(sys->context, 1, n + i, vertexCos[vid][1]);
+					nlRightHandSideSet(sys->context, 2, n + i, vertexCos[vid][2]);
 				}
 
-				nlEnd(NL_MATRIX);
-				nlEnd(NL_SYSTEM);
-				if (!nlSolve(NL_FALSE)) {
+				nlEnd(sys->context, NL_MATRIX);
+				nlEnd(sys->context, NL_SYSTEM);
+				if (!nlSolve(sys->context, NL_FALSE)) {
 					sys->has_solution = false;
 					break;
 				}
 			}
 			if (sys->has_solution) {
 				for (vid = 0; vid < sys->total_verts; vid++) {
-					vertexCos[vid][0] = nlGetVariable(0, vid);
-					vertexCos[vid][1] = nlGetVariable(1, vid);
-					vertexCos[vid][2] = nlGetVariable(2, vid);
+					vertexCos[vid][0] = nlGetVariable(sys->context, 0, vid);
+					vertexCos[vid][1] = nlGetVariable(sys->context, 1, vid);
+					vertexCos[vid][2] = nlGetVariable(sys->context, 2, vid);
 				}
 			}
 			else {
@@ -474,51 +473,49 @@ static void laplacianDeformPreview(LaplacianSystem *sys, float (*vertexCos)[3])
 
 	}
 	else if (sys->has_solution) {
-		nlMakeCurrent(sys->context);
-
-		nlBegin(NL_SYSTEM);
-		nlBegin(NL_MATRIX);
+		nlBegin(sys->context, NL_SYSTEM);
+		nlBegin(sys->context, NL_MATRIX);
 
 		for (i = 0; i < n; i++) {
-			nlRightHandSideSet(0, i, sys->delta[i][0]);
-			nlRightHandSideSet(1, i, sys->delta[i][1]);
-			nlRightHandSideSet(2, i, sys->delta[i][2]);
+			nlRightHandSideSet(sys->context, 0, i, sys->delta[i][0]);
+			nlRightHandSideSet(sys->context, 1, i, sys->delta[i][1]);
+			nlRightHandSideSet(sys->context, 2, i, sys->delta[i][2]);
 		}
 		for (i = 0; i < na; i++) {
 			vid = sys->index_anchors[i];
-			nlRightHandSideSet(0, n + i, vertexCos[vid][0]);
-			nlRightHandSideSet(1, n + i, vertexCos[vid][1]);
-			nlRightHandSideSet(2, n + i, vertexCos[vid][2]);
-			nlMatrixAdd(n + i, vid, 1.0f);
+			nlRightHandSideSet(sys->context, 0, n + i, vertexCos[vid][0]);
+			nlRightHandSideSet(sys->context, 1, n + i, vertexCos[vid][1]);
+			nlRightHandSideSet(sys->context, 2, n + i, vertexCos[vid][2]);
+			nlMatrixAdd(sys->context, n + i, vid, 1.0f);
 		}
 
-		nlEnd(NL_MATRIX);
-		nlEnd(NL_SYSTEM);
-		if (nlSolve(NL_FALSE)) {
+		nlEnd(sys->context, NL_MATRIX);
+		nlEnd(sys->context, NL_SYSTEM);
+		if (nlSolve(sys->context, NL_FALSE)) {
 			sys->has_solution = true;
 			for (j = 1; j <= sys->repeat; j++) {
-				nlBegin(NL_SYSTEM);
-				nlBegin(NL_MATRIX);
+				nlBegin(sys->context, NL_SYSTEM);
+				nlBegin(sys->context, NL_MATRIX);
 				rotateDifferentialCoordinates(sys);
 
 				for (i = 0; i < na; i++) {
 					vid = sys->index_anchors[i];
-					nlRightHandSideSet(0, n + i, vertexCos[vid][0]);
-					nlRightHandSideSet(1, n + i, vertexCos[vid][1]);
-					nlRightHandSideSet(2, n + i, vertexCos[vid][2]);
+					nlRightHandSideSet(sys->context, 0, n + i, vertexCos[vid][0]);
+					nlRightHandSideSet(sys->context, 1, n + i, vertexCos[vid][1]);
+					nlRightHandSideSet(sys->context, 2, n + i, vertexCos[vid][2]);
 				}
-				nlEnd(NL_MATRIX);
-				nlEnd(NL_SYSTEM);
-				if (!nlSolve(NL_FALSE)) {
+				nlEnd(sys->context, NL_MATRIX);
+				nlEnd(sys->context, NL_SYSTEM);
+				if (!nlSolve(sys->context, NL_FALSE)) {
 					sys->has_solution = false;
 					break;
 				}
 			}
 			if (sys->has_solution) {
 				for (vid = 0; vid < sys->total_verts; vid++) {
-					vertexCos[vid][0] = nlGetVariable(0, vid);
-					vertexCos[vid][1] = nlGetVariable(1, vid);
-					vertexCos[vid][2] = nlGetVariable(2, vid);
+					vertexCos[vid][0] = nlGetVariable(sys->context, 0, vid);
+					vertexCos[vid][1] = nlGetVariable(sys->context, 1, vid);
+					vertexCos[vid][2] = nlGetVariable(sys->context, 2, vid);
 				}
 			}
 			else {
diff --git a/source/blender/modifiers/intern/MOD_laplaciansmooth.c b/source/blender/modifiers/intern/MOD_laplaciansmooth.c
index aef28d24a51..189ceb11d08 100644
--- a/source/blender/modifiers/intern/MOD_laplaciansmooth.c
+++ b/source/blender/modifiers/intern/MOD_laplaciansmooth.c
@@ -300,16 +300,16 @@ static void fill_laplacian_matrix(LaplacianSystem *sys)
 
 			/* Is ring if number of faces == number of edges around vertice*/
 			if (sys->numNeEd[l_curr->v] == sys->numNeFa[l_curr->v] && sys->zerola[l_curr->v] == 0) {
-				nlMatrixAdd(l_curr->v, l_next->v, sys->fweights[l_curr_index][2] * sys->vweights[l_curr->v]);
-				nlMatrixAdd(l_curr->v, l_prev->v, sys->fweights[l_curr_index][1] * sys->vweights[l_curr->v]);
+				nlMatrixAdd(sys->context, l_curr->v, l_next->v, sys->fweights[l_curr_index][2] * sys->vweights[l_curr->v]);
+				nlMatrixAdd(sys->context, l_curr->v, l_prev->v, sys->fweights[l_curr_index][1] * sys->vweights[l_curr->v]);
 			}
 			if (sys->numNeEd[l_next->v] == sys->numNeFa[l_next->v] && sys->zerola[l_next->v] == 0) {
-				nlMatrixAdd(l_next->v, l_curr->v, sys->fweights[l_curr_index][2] * sys->vweights[l_next->v]);
-				nlMatrixAdd(l_next->v, l_prev->v, sys->fweights[l_curr_index][0] * sys->vweights[l_next->v]);
+				nlMatrixAdd(sys->context, l_next->v, l_curr->v, sys->fweights[l_curr_index][2] * sys->vweights[l_next->v]);
+				nlMatrixAdd(sys->context, l_next->v, l_prev->v, sys->fweights[l_curr_index][0] * sys->vweights[l_next->v]);
 			}
 			if (sys->numNeEd[l_prev->v] == sys->numNeFa[l_prev->v] && sys->zerola[l_prev->v] == 0) {
-				nlMatrixAdd(l_prev->v, l_curr->v, sys->fweights[l_curr_index][1] * sys->vweights[l_prev->v]);
-				nlMatrixAdd(l_prev->v, l_next->v, sys->fweights[l_curr_index][0] * sys->vweights[l_prev->v]);
+				nlMatrixAdd(sys->context, l_prev->v, l_curr->v, sys->fweights[l_curr_index][1] * sys->vweights[l_prev->v]);
+				nlMatrixAdd(sys->context, l_prev->v, l_next->v, sys->fweights[l_curr_index][0] * sys->vweights[l_prev->v]);
 			}
 		}
 	}
@@ -323,8 +323,8 @@ static void fill_laplacian_matrix(LaplacianSystem *sys)
 		    sys->zerola[idv1] == 0 &&
 		    sys->zerola[idv2] == 0)
 		{
-			nlMatrixAdd(idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
-			nlMatrixAdd(idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
+			nlMatrixAdd(sys->context, idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
+			nlMatrixAdd(sys->context, idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
 		}
 	}
 }
@@ -342,13 +342,13 @@ static void validate_solution(LaplacianSystem *sys, short flag, float lambda, fl
 		if (sys->zerola[i] == 0) {
 			lam = sys->numNeEd[i] == sys->numNeFa[i] ? (lambda >= 0.0f ? 1.0f : -1.0f) : (lambda_border >= 0.0f ? 1.0f : -1.0f);
 			if (flag & MOD_LAPLACIANSMOOTH_X) {
-				sys->vertexCos[i][0] += lam * ((float)nlGetVariable(0, i) - sys->vertexCos[i][0]);
+				sys->vertexCos[i][0] += lam * ((float)nlGetVariable(sys->context, 0, i) - sys->vertexCos[i][0]);
 			}
 			if (flag & MOD_LAPLACIANSMOOTH_Y) {
-				sys->vertexCos[i][1] += lam * ((float)nlGetVariable(1, i) - sys->vertexCos[i][1]);
+				sys->vertexCos[i][1] += lam * ((float)nlGetVariable(sys->context, 1, i) - sys->vertexCos[i][1]);
 			}
 			if (flag & MOD_LAPLACIANSMOOTH_Z) {
-				sys->vertexCos[i][2] += lam * ((float)nlGetVariable(2, i) - sys->vertexCos[i][2]);
+				sys->vertexCos[i][2] += lam * ((float)nlGetVariable(sys->context, 2, i) - sys->vertexCos[i][2]);
 			}
 		}
 	}
@@ -390,21 +390,20 @@ static void laplaciansmoothModifier_do(
 	modifier_opennl_lock();
 #endif
 
-	nlNewContext();
-	sys->context = nlGetCurrent();
-	nlSolverParameteri(NL_NB_VARIABLES, numVerts);
-	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
-	nlSolverParameteri(NL_NB_ROWS, numVerts);
-	nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 3);
+	sys->context = nlNewContext();
+	nlSolverParameteri(sys->context, NL_NB_VARIABLES, numVerts);
+	nlSolverParameteri(sys->context, NL_LEAST_SQUARES, NL_TRUE);
+	nlSolverParameteri(sys->context, NL_NB_ROWS, numVerts);
+	nlSolverParameteri(sys->context, NL_NB_RIGHT_HAND_SIDES, 3);
 
 	init_laplacian_matrix(sys);
 
 	for (iter = 0; iter < smd->repeat; iter++) {
-		nlBegin(NL_SYSTEM);
+		nlBegin(sys->context, NL_SYSTEM);
 		for (i = 0; i < numVerts; i++) {
-			nlSetVariable(0, i, vertexCos[i][0]);
-			nlSetVariable(1, i, vertexCos[i][1]);
-			nlSetVariable(2, i, vertexCos[i][2]);
+			nlSetVariable(sys->context, 0, i, vertexCos[i][0]);
+			nlSetVariable(sys->context, 1, i, vertexCos[i][1]);
+			nlSetVariable(sys->context, 2, i, vertexCos[i][2]);
 			if (iter == 0) {
 				add_v3_v3(sys->vert_centroid, vertexCos[i]);
 			}
@@ -413,12 +412,12 @@ static void laplaciansmoothModifier_do(
 			mul_v3_fl(sys->vert_centroid, 1.0f / (float)numVerts);
 		}
 
-		nlBegin(NL_MATRIX);
+		nlBegin(sys->context, NL_MATRIX);
 		dv = dvert;
 		for (i = 0; i < numVerts; i++) {
-			nlRightHandSideSet(0, i, vertexCos[i][0]);
-			nlRightHandSideSet(1, i, vertexCos[i][1]);
-			nlRightHandSideSet(2, i, vertexCos[i][2]);
+			nlRightHandSideSet(sys->context, 0, i, vertexCos[i][0]);
+			nlRightHandSideSet(sys->context, 1, i, vertexCos[i][1]);
+			nlRightHandSideSet(sys->context, 2, i, vertexCos[i][2]);
 			if (iter == 0) {
 				if (dv) {
 					wpaint = defvert_find_weight(dv, defgrp_index);
@@ -435,10 +434,10 @@ static void laplaciansmoothModifier_do(
 						w = sys->vlengths[i];
 						sys->vlengths[i] = (w == 0.0f) ? 0.0f : -fabsf(smd->lambda_border) * wpaint * 2.0f / w;
 						if (sys->numNeEd[i] == sys->numNeFa[i]) {
-							nlMatrixAdd(i, i,  1.0f + fabsf(smd->lambda) * wpaint);
+							nlMatrixAdd(sys->context, i, i,  1.0f + fabsf(smd->lambda) * wpaint);
 						}
 						else {
-							nlMatrixAdd(i, i,  1.0f + fabsf(smd->lambda_border) * wpaint * 2.0f);
+							nlMatrixAdd(sys->context, i, i,  1.0f + fabsf(smd->lambda_border) * wpaint * 2.0f);
 						}
 					}
 					else {
@@ -448,15 +447,15 @@ static void laplaciansmoothModifier_do(
 						sys->vlengths[i] = (w == 0.0f) ? 0.0f : -fabsf(smd->lambda_border) * wpaint * 2.0f / w;
 
 						if (sys->numNeEd[i] == sys->numNeFa[i]) {
-							nlMatrixAdd(i, i,  1.0f + fabsf(smd->lambda) * wpaint / (4.0f * sys->ring_areas[i]));
+							nlMatrixAdd(sys->context, i, i,  1.0f + fabsf(smd->lambda) * wpaint / (4.0f * sys->ring_areas[i]));
 						}
 						else {
-							nlMatrixAdd(i, i,  1.0f + fabsf(smd->lambda_border) * wpaint * 2.0f);
+							nlMatrixAdd(sys->context, i, i,  1.0f + fabsf(smd->lambda_border) * wpaint * 2.0f);
 						}
 					}
 				}
 				else {
-					nlMatrixAdd(i, i, 1.0f);
+					nlMatrixAdd(sys->context, i, i, 1.0f);
 				}
 			}
 		}
@@ -465,10 +464,10 @@ static void laplaciansmoothModifier_do(
 			fill_laplacian_matrix(sys);
 		}
 
-		nlEnd(NL_MATRIX);
-		nlEnd(NL_SYSTEM);
+		nlEnd(sys->context, NL_MATRIX);
+		nlEnd(sys->context, NL_SYSTEM);
 
-		if (nlSolve(NL_TRUE)) {
+		if (nlSolve(sys->context, NL_TRUE)) {
 			validate_solution(sys, smd->flag, smd->lambda, smd->lambda_border);
 		}
 	}