diff --git a/intern/opennl/extern/ONL_opennl.h b/intern/opennl/extern/ONL_opennl.h
index 7f501629290..be76aa95eac 100644
--- a/intern/opennl/extern/ONL_opennl.h
+++ b/intern/opennl/extern/ONL_opennl.h
@@ -79,20 +79,16 @@ typedef void* NLContext;
 
 #define NL_SYSTEM  0x0
 #define NL_MATRIX  0x1
-#define NL_ROW     0x2
 
 /* Solver Parameters */
 
-#define NL_SOLVER           0x100
-#define NL_NB_VARIABLES     0x101
-#define NL_LEAST_SQUARES    0x102
-#define NL_SYMMETRIC        0x106
-#define NL_ERROR            0x108
-
-/* Row parameters */
-
-#define NL_RIGHT_HAND_SIDE 0x500
-#define NL_ROW_SCALING     0x501
+#define NL_SOLVER              0x100
+#define NL_NB_VARIABLES        0x101
+#define NL_LEAST_SQUARES       0x102
+#define NL_SYMMETRIC           0x106
+#define NL_ERROR               0x108
+#define NL_NB_ROWS             0x110
+#define NL_NB_RIGHT_HAND_SIDES 0x112 /* 4 max */
 
 /* Contexts */
 
@@ -106,9 +102,6 @@ NLContext nlGetCurrent(void);
 void nlSolverParameterf(NLenum pname, NLfloat param);
 void nlSolverParameteri(NLenum pname, NLint param);
 
-void nlRowParameterf(NLenum pname, NLfloat param);
-void nlRowParameteri(NLenum pname, NLint param);
-
 void nlGetBooleanv(NLenum pname, NLboolean* params);
 void nlGetFloatv(NLenum pname, NLfloat* params);
 void nlGetIntergerv(NLenum pname, NLint* params);
@@ -119,8 +112,8 @@ NLboolean nlIsEnabled(NLenum pname);
 
 /* Variables */
 
-void nlSetVariable(NLuint index, NLfloat value);
-NLfloat nlGetVariable(NLuint index);
+void nlSetVariable(NLuint rhsindex, NLuint index, NLfloat value);
+NLfloat nlGetVariable(NLuint rhsindex, NLuint index);
 void nlLockVariable(NLuint index);
 void nlUnlockVariable(NLuint index);
 NLboolean nlVariableIsLocked(NLuint index);
@@ -129,14 +122,16 @@ NLboolean nlVariableIsLocked(NLuint index);
 
 void nlBegin(NLenum primitive);
 void nlEnd(NLenum primitive);
-void nlCoefficient(NLuint index, NLfloat value);
 
-/* Setting random elements matrix/vector - not supported for
-   least squares! */
+/* Setting elements in matrix/vector */
 
 void nlMatrixAdd(NLuint row, NLuint col, NLfloat value);
-void nlRightHandSideAdd(NLuint index, NLfloat value);
-void nlRightHandSideSet(NLuint index, NLfloat value);
+void nlRightHandSideAdd(NLuint rhsindex, NLuint index, NLfloat value);
+void nlRightHandSideSet(NLuint rhsindex, NLuint index, NLfloat value);
+
+/* Multiply */
+
+void nlMatrixMultiply(NLfloat *x, NLfloat *y);
 
 /* Solve */
 
diff --git a/intern/opennl/intern/opennl.c b/intern/opennl/intern/opennl.c
index 836c9d01e60..2d30da075d3 100644
--- a/intern/opennl/intern/opennl.c
+++ b/intern/opennl/intern/opennl.c
@@ -207,10 +207,6 @@ static void __nlRowColumnAppend(__NLRowColumn* c, NLint index, NLfloat value) {
 	c->size++;
 }
 
-static void __nlRowColumnZero(__NLRowColumn* c) {
-	c->size = 0;
-}
-
 static void __nlRowColumnClear(__NLRowColumn* c) {
 	c->size	 = 0;
 	c->capacity = 0;
@@ -432,13 +428,62 @@ static void __nlSparseMatrixMult(__NLSparseMatrix* A, NLfloat* x, NLfloat* y) {
 	}
 }
 
+/* ****************** Routines for least squares ******************* */
+
+static void __nlSparseMatrix_square(
+	__NLSparseMatrix* AtA, __NLSparseMatrix *A
+) {
+	NLuint m = A->m;
+	NLuint n = A->n;
+	NLuint i, j0, j1;
+	__NLRowColumn *Ri = NULL;
+	__NLCoeff *c0 = NULL, *c1 = NULL;
+	float value;
+
+	__nlSparseMatrixConstruct(AtA, n, n, A->storage);
+
+	for(i=0; i<m; i++) {
+		Ri = &(A->row[i]);
+
+		for(j0=0; j0<Ri->size; j0++) {
+			c0 = &(Ri->coeff[j0]);
+			for(j1=0; j1<Ri->size; j1++) {
+				c1 = &(Ri->coeff[j1]);
+
+				value = c0->value*c1->value;
+				__nlSparseMatrixAdd(AtA, c0->index, c1->index, value);
+			}
+		}
+	}
+}
+
+static void __nlSparseMatrix_transpose_mult_rows(
+	__NLSparseMatrix* A, NLfloat* x, NLfloat* y
+) {
+	NLuint m = A->m;
+	NLuint n = A->n;
+	NLuint i,ij;
+	__NLRowColumn* Ri = NULL;
+	__NLCoeff* c = NULL;
+
+	__NL_CLEAR_ARRAY(NLfloat, y, n);
+
+	for(i=0; i<m; i++) {
+		Ri = &(A->row[i]);
+		for(ij=0; ij<Ri->size; ij++) {
+			c = &(Ri->coeff[ij]);
+			y[c->index] += c->value * x[i];
+		}
+	}
+}
+
 /************************************************************************************/
 /* NLContext data structure */
 
 typedef void(*__NLMatrixFunc)(float* x, float* y);
 
 typedef struct {
-	NLfloat  value;
+	NLfloat  value[4];
 	NLboolean locked;
 	NLuint	index;
 	__NLRowColumn *a;
@@ -447,32 +492,32 @@ typedef struct {
 #define __NL_STATE_INITIAL				0
 #define __NL_STATE_SYSTEM				1
 #define __NL_STATE_MATRIX				2
-#define __NL_STATE_ROW					3
-#define __NL_STATE_MATRIX_CONSTRUCTED	4
-#define __NL_STATE_SYSTEM_CONSTRUCTED	5
-#define __NL_STATE_SYSTEM_SOLVED		7
+#define __NL_STATE_MATRIX_CONSTRUCTED	3
+#define __NL_STATE_SYSTEM_CONSTRUCTED	4
+#define __NL_STATE_SYSTEM_SOLVED		5
 
 typedef struct {
-	NLenum		   state;
-	NLuint		   n;
+	NLenum			state;
+	NLuint			n;
+	NLuint			m;
 	__NLVariable*	variable;
 	NLfloat*		b;
+	NLfloat*		Mtb;
 	__NLSparseMatrix M;
-	__NLRowColumn	af;
-	__NLRowColumn	al;
+	__NLSparseMatrix MtM;
 	NLfloat*		x;
-	NLfloat		 right_hand_side;
-	NLuint		   nb_variables;
-	NLuint		   current_row;
+	NLuint			nb_variables;
+	NLuint			nb_rows;
 	NLboolean		least_squares;
 	NLboolean		symmetric;
+	NLuint			nb_rhs;
 	NLboolean		solve_again;
 	NLboolean		alloc_M;
-	NLboolean		alloc_af;
-	NLboolean		alloc_al;
+	NLboolean		alloc_MtM;
 	NLboolean		alloc_variable;
 	NLboolean		alloc_x;
 	NLboolean		alloc_b;
+	NLboolean		alloc_Mtb;
 	NLfloat			error;
 	__NLMatrixFunc	matrix_vector_prod;
 
@@ -494,8 +539,8 @@ static void __nlMatrixVectorProd_default(NLfloat* x, NLfloat* y) {
 NLContext nlNewContext(void) {
 	__NLContext* result	  = __NL_NEW(__NLContext);
 	result->state			= __NL_STATE_INITIAL;
-	result->right_hand_side  = 0.0;
 	result->matrix_vector_prod = __nlMatrixVectorProd_default;
+	result->nb_rhs = 1;
 	nlMakeCurrent(result);
 	return result;
 }
@@ -512,11 +557,8 @@ void nlDeleteContext(NLContext context_in) {
 	if(context->alloc_M) {
 		__nlSparseMatrixDestroy(&context->M);
 	}
-	if(context->alloc_af) {
-		__nlRowColumnDestroy(&context->af);
-	}
-	if(context->alloc_al) {
-		__nlRowColumnDestroy(&context->al);
+	if(context->alloc_MtM) {
+		__nlSparseMatrixDestroy(&context->MtM);
 	}
 	if(context->alloc_variable) {
 		for(i=0; i<context->nb_variables; i++) {
@@ -529,6 +571,9 @@ void nlDeleteContext(NLContext context_in) {
 	if(context->alloc_b) {
 		__NL_DELETE_ARRAY(context->b);
 	}
+	if(context->alloc_Mtb) {
+		__NL_DELETE_ARRAY(context->Mtb);
+	}
 	if(context->alloc_x) {
 		__NL_DELETE_ARRAY(context->x);
 	}
@@ -569,12 +614,19 @@ void nlSolverParameterf(NLenum pname, NLfloat param) {
 		__nl_assert(param > 0);
 		__nlCurrentContext->nb_variables = (NLuint)param;
 	} break;
+	case NL_NB_ROWS: {
+		__nl_assert(param > 0);
+		__nlCurrentContext->nb_rows = (NLuint)param;
+	} break;
 	case NL_LEAST_SQUARES: {
 		__nlCurrentContext->least_squares = (NLboolean)param;
 	} break;
 	case NL_SYMMETRIC: {
 		__nlCurrentContext->symmetric = (NLboolean)param;		
-	}
+	} break;
+	case NL_NB_RIGHT_HAND_SIDES: {
+		__nlCurrentContext->nb_rhs = (NLuint)param;
+	} break;
 	default: {
 		__nl_assert_not_reached;
 	} break;
@@ -588,36 +640,25 @@ void nlSolverParameteri(NLenum pname, NLint param) {
 		__nl_assert(param > 0);
 		__nlCurrentContext->nb_variables = (NLuint)param;
 	} break;
+	case NL_NB_ROWS: {
+		__nl_assert(param > 0);
+		__nlCurrentContext->nb_rows = (NLuint)param;
+	} break;
 	case NL_LEAST_SQUARES: {
 		__nlCurrentContext->least_squares = (NLboolean)param;
 	} break;
 	case NL_SYMMETRIC: {
 		__nlCurrentContext->symmetric = (NLboolean)param;		
-	}
+	} break;
+	case NL_NB_RIGHT_HAND_SIDES: {
+		__nlCurrentContext->nb_rhs = (NLuint)param;
+	} break;
 	default: {
 		__nl_assert_not_reached;
 	} break;
 	}
 }
 
-void nlRowParameterf(NLenum pname, NLfloat param) {
-	__nlCheckState(__NL_STATE_MATRIX);
-	switch(pname) {
-	case NL_RIGHT_HAND_SIDE: {
-		__nlCurrentContext->right_hand_side = param;
-	} break;
-	}
-}
-
-void nlRowParameteri(NLenum pname, NLint param) {
-	__nlCheckState(__NL_STATE_MATRIX);
-	switch(pname) {
-	case NL_RIGHT_HAND_SIDE: {
-		__nlCurrentContext->right_hand_side = (NLfloat)param;
-	} break;
-	}
-}
-
 void nlGetBooleanv(NLenum pname, NLboolean* params) {
 	switch(pname) {
 	case NL_LEAST_SQUARES: {
@@ -637,6 +678,9 @@ void nlGetFloatv(NLenum pname, NLfloat* params) {
 	case NL_NB_VARIABLES: {
 		*params = (NLfloat)(__nlCurrentContext->nb_variables);
 	} break;
+	case NL_NB_ROWS: {
+		*params = (NLfloat)(__nlCurrentContext->nb_rows);
+	} break;
 	case NL_LEAST_SQUARES: {
 		*params = (NLfloat)(__nlCurrentContext->least_squares);
 	} break;
@@ -657,6 +701,9 @@ void nlGetIntergerv(NLenum pname, NLint* params) {
 	case NL_NB_VARIABLES: {
 		*params = (NLint)(__nlCurrentContext->nb_variables);
 	} break;
+	case NL_NB_ROWS: {
+		*params = (NLint)(__nlCurrentContext->nb_rows);
+	} break;
 	case NL_LEAST_SQUARES: {
 		*params = (NLint)(__nlCurrentContext->least_squares);
 	} break;
@@ -700,16 +747,16 @@ NLboolean nlIsEnabled(NLenum pname) {
 /************************************************************************************/
 /* Get/Set Lock/Unlock variables */
 
-void nlSetVariable(NLuint index, NLfloat value) {
+void nlSetVariable(NLuint rhsindex, NLuint index, NLfloat value) {
 	__nlCheckState(__NL_STATE_SYSTEM);
 	__nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1);
-	__nlCurrentContext->variable[index].value = value;	
+	__nlCurrentContext->variable[index].value[rhsindex] = value;	
 }
 
-NLfloat nlGetVariable(NLuint index) {
+NLfloat nlGetVariable(NLuint rhsindex, NLuint index) {
 	__nl_assert(__nlCurrentContext->state != __NL_STATE_INITIAL);
 	__nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1);
-	return __nlCurrentContext->variable[index].value;
+	return __nlCurrentContext->variable[index].value[rhsindex];
 }
 
 void nlLockVariable(NLuint index) {
@@ -734,31 +781,41 @@ NLboolean nlVariableIsLocked(NLuint index) {
 /* System construction */
 
 static void __nlVariablesToVector() {
-	NLuint i;
+	__NLContext *context = __nlCurrentContext;
+	NLuint i, j, nb_rhs;
 
-	__nl_assert(__nlCurrentContext->alloc_x);
-	__nl_assert(__nlCurrentContext->alloc_variable);
+	__nl_assert(context->alloc_x);
+	__nl_assert(context->alloc_variable);
 
-	for(i=0; i<__nlCurrentContext->nb_variables; i++) {
-		__NLVariable* v = &(__nlCurrentContext->variable[i]);
+	nb_rhs= context->nb_rhs;
+
+	for(i=0; i<context->nb_variables; i++) {
+		__NLVariable* v = &(context->variable[i]);
 		if(!v->locked) {
-			__nl_assert(v->index < __nlCurrentContext->n);
-			__nlCurrentContext->x[v->index] = v->value;
+			__nl_assert(v->index < context->n);
+
+			for(j=0; j<nb_rhs; j++)
+				context->x[context->n*j + v->index] = v->value[j];
 		}
 	}
 }
 
 static void __nlVectorToVariables() {
-	NLuint i;
+	__NLContext *context = __nlCurrentContext;
+	NLuint i, j, nb_rhs;
 
-	__nl_assert(__nlCurrentContext->alloc_x);
-	__nl_assert(__nlCurrentContext->alloc_variable);
+	__nl_assert(context->alloc_x);
+	__nl_assert(context->alloc_variable);
 
-	for(i=0; i<__nlCurrentContext->nb_variables; i++) {
-		__NLVariable* v = &(__nlCurrentContext->variable[i]);
+	nb_rhs= context->nb_rhs;
+
+	for(i=0; i<context->nb_variables; i++) {
+		__NLVariable* v = &(context->variable[i]);
 		if(!v->locked) {
-			__nl_assert(v->index < __nlCurrentContext->n);
-			v->value = __nlCurrentContext->x[v->index];
+			__nl_assert(v->index < context->n);
+
+			for(j=0; j<nb_rhs; j++)
+				v->value[j] = context->x[context->n*j + v->index];
 		}
 	}
 }
@@ -783,8 +840,8 @@ static void __nlEndSystem() {
 }
 
 static void __nlBeginMatrix() {
-	NLuint i, j;
-	NLuint n = 0;
+	NLuint i;
+	NLuint m = 0, n = 0;
 	NLenum storage = __NL_ROWS;
 	__NLContext *context = __nlCurrentContext;
 
@@ -801,57 +858,37 @@ static void __nlBeginMatrix() {
 				context->variable[i].index = n++;
 		}
 
+		m = (context->nb_rows == 0)? n: context->nb_rows;
+
+		context->m = m;
 		context->n = n;
 
-		/* a least squares problem results in a symmetric matrix */
-		if(context->least_squares)
-			context->symmetric = NL_TRUE;
-
-		if(context->symmetric)
-			storage = (storage | __NL_SYMMETRIC);
-
-		/* SuperLU storage does not support symmetric storage */
-		storage = (storage & ~__NL_SYMMETRIC);
-
-		__nlSparseMatrixConstruct(&context->M, n, n, storage);
+		__nlSparseMatrixConstruct(&context->M, m, n, storage);
 		context->alloc_M = NL_TRUE;
 
-		context->b = __NL_NEW_ARRAY(NLfloat, n);
+		context->b = __NL_NEW_ARRAY(NLfloat, m*context->nb_rhs);
 		context->alloc_b = NL_TRUE;
 
-		context->x = __NL_NEW_ARRAY(NLfloat, n);
+		context->x = __NL_NEW_ARRAY(NLfloat, n*context->nb_rhs);
 		context->alloc_x = NL_TRUE;
 	}
 	else {
 		/* need to recompute b only, A is not constructed anymore */
-		__NL_CLEAR_ARRAY(NLfloat, context->b, context->n);
+		__NL_CLEAR_ARRAY(NLfloat, context->b, context->m*context->nb_rhs);
 	}
 
 	__nlVariablesToVector();
-
-	__nlRowColumnConstruct(&context->af);
-	context->alloc_af = NL_TRUE;
-	__nlRowColumnConstruct(&context->al);
-	context->alloc_al = NL_TRUE;
-
-	context->current_row = 0;
 }
 
-static void __nlEndMatrix() {
+static void __nlEndMatrixRHS(NLuint rhs) {
 	__NLContext *context = __nlCurrentContext;
 	__NLVariable *variable;
 	__NLRowColumn *a;
-	NLfloat *b;
+	NLfloat *b, *Mtb;
 	NLuint i, j;
 
-	__nlTransition(__NL_STATE_MATRIX, __NL_STATE_MATRIX_CONSTRUCTED);	
-	
-	__nlRowColumnDestroy(&context->af);
-	context->alloc_af = NL_FALSE;
-	__nlRowColumnDestroy(&context->al);
-	context->alloc_al = NL_FALSE;
-	
-	b = context->b;
+	b = context->b + context->m*rhs;
+	Mtb = context->Mtb + context->n*rhs;
 
 	for(i=0; i<__nlCurrentContext->nb_variables; i++) {
 		variable = &(context->variable[i]);
@@ -860,73 +897,34 @@ static void __nlEndMatrix() {
 			a = variable->a;
 
 			for(j=0; j<a->size; j++) {
-				b[a->coeff[j].index] -= a->coeff[j].value*variable->value;
+				b[a->coeff[j].index] -= a->coeff[j].value*variable->value[rhs];
 			}
 		}
 	}
 
-#if 0
-	if(!context->least_squares) {
-		__nl_assert(
-			context->current_row == 
-			context->n
-		);
-	}
-#endif
+	if(context->least_squares)
+		__nlSparseMatrix_transpose_mult_rows(&context->M, b, Mtb);
 }
 
-static void __nlBeginRow() {
-	__nlTransition(__NL_STATE_MATRIX, __NL_STATE_ROW);
-	__nlRowColumnZero(&__nlCurrentContext->af);
-	__nlRowColumnZero(&__nlCurrentContext->al);
-}
-
-static void __nlEndRow() {
-	__NLRowColumn*	af = &__nlCurrentContext->af;
-	__NLRowColumn*	al = &__nlCurrentContext->al;
-	__NLSparseMatrix* M  = &__nlCurrentContext->M;
-	NLfloat* b		= __nlCurrentContext->b;
-	NLuint nf		  = af->size;
-	NLuint nl		  = al->size;
-	NLuint current_row = __nlCurrentContext->current_row;
+static void __nlEndMatrix() {
+	__NLContext *context = __nlCurrentContext;
 	NLuint i;
-	NLuint j;
-	NLfloat S;
-	__nlTransition(__NL_STATE_ROW, __NL_STATE_MATRIX);
 
-	if(__nlCurrentContext->least_squares) {
-		if (!__nlCurrentContext->solve_again) {
-			for(i=0; i<nf; i++) {
-				for(j=0; j<nf; j++) {
-					__nlSparseMatrixAdd(
-						M, af->coeff[i].index, af->coeff[j].index,
-						af->coeff[i].value * af->coeff[j].value
-					);
-				}
-			}
-		}
+	__nlTransition(__NL_STATE_MATRIX, __NL_STATE_MATRIX_CONSTRUCTED);	
+	
+	if(context->least_squares) {
+		if(!__nlCurrentContext->solve_again) {
+			__nlSparseMatrix_square(&context->MtM, &context->M);
+			context->alloc_MtM = NL_TRUE;
 
-		S = -__nlCurrentContext->right_hand_side;
-		for(j=0; j<nl; j++)
-			S += al->coeff[j].value;
-
-		for(i=0; i<nf; i++)
-			b[ af->coeff[i].index ] -= af->coeff[i].value * S;
-	} else {
-		if (!__nlCurrentContext->solve_again) {
-			for(i=0; i<nf; i++) {
-				__nlSparseMatrixAdd(
-					M, current_row, af->coeff[i].index, af->coeff[i].value
-				);
-			}
-		}
-		b[current_row] = -__nlCurrentContext->right_hand_side;
-		for(i=0; i<nl; i++) {
-			b[current_row] -= al->coeff[i].value;
+			context->Mtb =
+				__NL_NEW_ARRAY(NLfloat, context->n*context->nb_rhs);
+			context->alloc_Mtb = NL_TRUE;
 		}
 	}
-	__nlCurrentContext->current_row++;
-	__nlCurrentContext->right_hand_side = 0.0;	
+
+	for(i=0; i<context->nb_rhs; i++)
+		__nlEndMatrixRHS(i);
 }
 
 void nlMatrixAdd(NLuint row, NLuint col, NLfloat value)
@@ -934,68 +932,70 @@ void nlMatrixAdd(NLuint row, NLuint col, NLfloat value)
 	__NLContext *context = __nlCurrentContext;
 
 	__nlCheckState(__NL_STATE_MATRIX);
-	__nl_assert(!context->least_squares);
 
-	if (context->variable[row].locked);
+	if(context->solve_again)
+		return;
+
+	if (!context->least_squares && context->variable[row].locked);
 	else if (context->variable[col].locked) {
-		row = context->variable[row].index;
+		if(!context->least_squares)
+			row = context->variable[row].index;
 		__nlRowColumnAppend(context->variable[col].a, row, value);
 	}
 	else {
 		__NLSparseMatrix* M  = &context->M;
-
-		row = context->variable[row].index;
+		
+		if(!context->least_squares)
+			row = context->variable[row].index;
 		col = context->variable[col].index;
 		
-		__nl_range_assert(row, 0, context->n - 1);
+		__nl_range_assert(row, 0, context->m - 1);
 		__nl_range_assert(col, 0, context->n - 1);
 
 		__nlSparseMatrixAdd(M, row, col, value);
 	}
 }
 
-void nlRightHandSideAdd(NLuint index, NLfloat value)
+void nlRightHandSideAdd(NLuint rhsindex, NLuint index, NLfloat value)
 {
 	__NLContext *context = __nlCurrentContext;
 	NLfloat* b = context->b;
 
 	__nlCheckState(__NL_STATE_MATRIX);
-	__nl_assert(!context->least_squares);
 
-	if(!context->variable[index].locked) {
-		index = context->variable[index].index;
-		__nl_range_assert(index, 0, context->n - 1);
+	if(context->least_squares) {
+		__nl_range_assert(index, 0, context->m - 1);
+		b[rhsindex*context->m + index] += value;
+	}
+	else {
+		if(!context->variable[index].locked) {
+			index = context->variable[index].index;
+			__nl_range_assert(index, 0, context->m - 1);
 
-		b[index] += value;
+			b[rhsindex*context->m + index] += value;
+		}
 	}
 }
 
-void nlRightHandSideSet(NLuint index, NLfloat value)
+void nlRightHandSideSet(NLuint rhsindex, NLuint index, NLfloat value)
 {
 	__NLContext *context = __nlCurrentContext;
 	NLfloat* b = context->b;
 
 	__nlCheckState(__NL_STATE_MATRIX);
-	__nl_assert(!context->least_squares);
 
-	if(!context->variable[index].locked) {
-		index = context->variable[index].index;
-		__nl_range_assert(index, 0, context->n - 1);
-
-		b[index] = value;
+	if(context->least_squares) {
+		__nl_range_assert(index, 0, context->m - 1);
+		b[rhsindex*context->m + index] = value;
 	}
-}
+	else {
+		if(!context->variable[index].locked) {
+			index = context->variable[index].index;
+			__nl_range_assert(index, 0, context->m - 1);
 
-void nlCoefficient(NLuint index, NLfloat value) {
-	__NLVariable* v;
-	unsigned int zero= 0;
-	__nlCheckState(__NL_STATE_ROW);
-	__nl_range_assert(index, zero, __nlCurrentContext->nb_variables - 1);
-	v = &(__nlCurrentContext->variable[index]);
-	if(v->locked)
-		__nlRowColumnAppend(&(__nlCurrentContext->al), 0, value*v->value);
-	else
-		__nlRowColumnAppend(&(__nlCurrentContext->af), v->index, value);
+			b[rhsindex*context->m + index] = value;
+		}
+	}
 }
 
 void nlBegin(NLenum prim) {
@@ -1006,9 +1006,6 @@ void nlBegin(NLenum prim) {
 	case NL_MATRIX: {
 		__nlBeginMatrix();
 	} break;
-	case NL_ROW: {
-		__nlBeginRow();
-	} break;
 	default: {
 		__nl_assert_not_reached;
 	}
@@ -1023,9 +1020,6 @@ void nlEnd(NLenum prim) {
 	case NL_MATRIX: {
 		__nlEndMatrix();
 	} break;
-	case NL_ROW: {
-		__nlEndRow();
-	} break;
 	default: {
 		__nl_assert_not_reached;
 	}
@@ -1040,7 +1034,7 @@ void nlEnd(NLenum prim) {
 static NLboolean __nlFactorize_SUPERLU(__NLContext *context, NLint *permutation) {
 
 	/* OpenNL Context */
-	__NLSparseMatrix* M = &(context->M);
+	__NLSparseMatrix* M = (context->least_squares)? &context->MtM: &context->M;
 	NLuint n = context->n;
 	NLuint nnz = __nlSparseMatrixNNZ(M); /* number of non-zero coeffs */
 
@@ -1057,7 +1051,6 @@ static NLboolean __nlFactorize_SUPERLU(__NLContext *context, NLint *permutation)
 	superlu_options_t options;
 
 	/* Temporary variables */
-	__NLRowColumn* Ri = NULL;
 	NLuint i, jj, count;
 	
 	__nl_assert(!(M->storage & __NL_SYMMETRIC));
@@ -1136,31 +1129,33 @@ static NLboolean __nlFactorize_SUPERLU(__NLContext *context, NLint *permutation)
 static NLboolean __nlInvert_SUPERLU(__NLContext *context) {
 
 	/* OpenNL Context */
-	NLfloat* b = context->b;
+	NLfloat* b = (context->least_squares)? context->Mtb: context->b;
 	NLfloat* x = context->x;
-	NLuint n = context->n;
+	NLuint n = context->n, j;
 
 	/* SuperLU variables */
 	SuperMatrix B;
 	NLint info;
 
-	/* Create superlu array for B */
-	sCreate_Dense_Matrix(
-		&B, n, 1, b, n, 
-		SLU_DN, /* Fortran-type column-wise storage */
-		SLU_S,  /* floats						  */
-		SLU_GE  /* general						  */
-	);
+	for(j=0; j<context->nb_rhs; j++, b+=n, x+=n) {
+		/* Create superlu array for B */
+		sCreate_Dense_Matrix(
+			&B, n, 1, b, n, 
+			SLU_DN, /* Fortran-type column-wise storage */
+			SLU_S,  /* floats						  */
+			SLU_GE  /* general						  */
+		);
 
-	/* Forward/Back substitution to compute x */
-	sgstrs(TRANS, &(context->slu.L), &(context->slu.U),
-		context->slu.perm_c, context->slu.perm_r, &B,
-		&(context->slu.stat), &info);
+		/* Forward/Back substitution to compute x */
+		sgstrs(TRANS, &(context->slu.L), &(context->slu.U),
+			context->slu.perm_c, context->slu.perm_r, &B,
+			&(context->slu.stat), &info);
 
-	if(info == 0)
-		memcpy(x, ((DNformat*)B.Store)->nzval, sizeof(*x)*n);
+		if(info == 0)
+			memcpy(x, ((DNformat*)B.Store)->nzval, sizeof(*x)*n);
 
-	Destroy_SuperMatrix_Store(&B);
+		Destroy_SuperMatrix_Store(&B);
+	}
 
 	return (info == 0);
 }
@@ -1179,15 +1174,18 @@ static void __nlFree_SUPERLU(__NLContext *context) {
 }
 
 void nlPrintMatrix(void) {
-	__NLSparseMatrix* M  = &(__nlCurrentContext->M);
-	float *b = __nlCurrentContext->b;
+	__NLContext *context = __nlCurrentContext;
+	__NLSparseMatrix* M  = &(context->M);
+	__NLSparseMatrix* MtM  = &(context->MtM);
+	float *b = context->b;
 	NLuint i, jj, k;
-	NLuint n = __nlCurrentContext->n;
+	NLuint m = context->m;
+	NLuint n = context->n;
 	__NLRowColumn* Ri = NULL;
-	float *value = malloc(sizeof(*value)*n);
+	float *value = malloc(sizeof(*value)*(n+m));
 
 	printf("A:\n");
-	for(i=0; i<n; i++) {
+	for(i=0; i<m; i++) {
 		Ri = &(M->row[i]);
 
 		memset(value, 0.0, sizeof(*value)*n);
@@ -1199,10 +1197,35 @@ void nlPrintMatrix(void) {
 		printf("\n");
 	}
 
-	printf("b:\n");
-	for(i=0; i<n; i++)
-		printf("%f ", b[i]);
-	printf("\n");
+	for(k=0; k<context->nb_rhs; k++) {
+		printf("b (%d):\n", k);
+		for(i=0; i<n; i++)
+			printf("%f ", b[context->n*k + i]);
+		printf("\n");
+	}
+
+	if(context->alloc_MtM) {
+		printf("AtA:\n");
+		for(i=0; i<n; i++) {
+			Ri = &(MtM->row[i]);
+
+			memset(value, 0.0, sizeof(*value)*m);
+			for(jj=0; jj<Ri->size; jj++)
+				value[Ri->coeff[jj].index] = Ri->coeff[jj].value;
+
+			for (k = 0; k<n; k++)
+				printf("%.3f ", value[k]);
+			printf("\n");
+		}
+
+		for(k=0; k<context->nb_rhs; k++) {
+			printf("Mtb (%d):\n", k);
+			for(i=0; i<n; i++)
+				printf("%f ", context->Mtb[context->n*k + i]);
+			printf("\n");
+		}
+		printf("\n");
+	}
 
 	free(value);
 }
diff --git a/source/blender/blenkernel/BKE_bad_level_calls.h b/source/blender/blenkernel/BKE_bad_level_calls.h
index 32128c68af5..b942add80e2 100644
--- a/source/blender/blenkernel/BKE_bad_level_calls.h
+++ b/source/blender/blenkernel/BKE_bad_level_calls.h
@@ -225,5 +225,11 @@ void antialias_tagbuf(int xsize, int ysize, char *rectmove);
 /* imagetexture.c */
 void ibuf_sample(struct ImBuf *ibuf, float fx, float fy, float dx, float dy, float *result);
 
+/* modifier.c */
+struct MeshDeformModifierData;
+
+void harmonic_coordinates_bind(struct MeshDeformModifierData *mmd,
+	float (*vertexcos)[3], int totvert, float cagemat[][4]);
+
 #endif
 
diff --git a/source/blender/blenkernel/BKE_mesh.h b/source/blender/blenkernel/BKE_mesh.h
index e9fd059f833..a8e4969ad43 100644
--- a/source/blender/blenkernel/BKE_mesh.h
+++ b/source/blender/blenkernel/BKE_mesh.h
@@ -90,6 +90,7 @@ void mesh_calc_normals(struct MVert *mverts, int numVerts, struct MFace *mfaces,
 	/* Return a newly MEM_malloc'd array of all the mesh vertex locations
 	 * (_numVerts_r_ may be NULL) */
 float (*mesh_getVertexCos(struct Mesh *me, int *numVerts_r))[3];
+float (*mesh_getRefKeyCos(struct Mesh *me, int *numVerts_r))[3];
 
 /* map from uv vertex to face (for select linked, stitch, uv suburf) */
 
diff --git a/source/blender/blenkernel/bad_level_call_stubs/stubs.c b/source/blender/blenkernel/bad_level_call_stubs/stubs.c
index aca8b916323..63bc23a71bf 100644
--- a/source/blender/blenkernel/bad_level_call_stubs/stubs.c
+++ b/source/blender/blenkernel/bad_level_call_stubs/stubs.c
@@ -332,3 +332,7 @@ void BIF_filelist_freelib(struct FileList* filelist) {};
 /* edittime.c stub */
 TimeMarker *get_frame_marker(int frame){return 0;};
 
+/* modifier.c stub */
+void harmonic_coordinates_bind(struct MeshDeformModifierData *mmd,
+	float (*vertexcos)[3], int totvert, float cagemat[][4]) {}
+
diff --git a/source/blender/blenkernel/intern/DerivedMesh.c b/source/blender/blenkernel/intern/DerivedMesh.c
index 04b111ffe2f..4e6d4a31173 100644
--- a/source/blender/blenkernel/intern/DerivedMesh.c
+++ b/source/blender/blenkernel/intern/DerivedMesh.c
@@ -1917,7 +1917,10 @@ static void mesh_calc_modifiers(Object *ob, float (*inputVertexCos)[3],
 	} else {
 		if(!fluidsimMeshUsed) {
 			/* default behaviour for meshes */
-			deformedVerts = inputVertexCos;
+			if(inputVertexCos)
+				deformedVerts = inputVertexCos;
+			else
+				deformedVerts = mesh_getRefKeyCos(me, &numVerts);
 		} else {
 			/* the fluid sim mesh might have more vertices than the original 
 			 * one, so inputVertexCos shouldnt be used
diff --git a/source/blender/blenkernel/intern/mesh.c b/source/blender/blenkernel/intern/mesh.c
index e6e6ba49066..4e551e28885 100644
--- a/source/blender/blenkernel/intern/mesh.c
+++ b/source/blender/blenkernel/intern/mesh.c
@@ -486,16 +486,7 @@ static float *make_orco_mesh_internal(Object *ob, int render)
 		/* Get appropriate vertex coordinates */
 
 	if(me->key && me->texcomesh==0 && me->key->refkey) {
-		KeyBlock *kb= me->key->refkey;
-		float *fp= kb->data;
-		totvert= MIN2(kb->totelem, me->totvert);
-		vcos = MEM_callocN(sizeof(*vcos)*me->totvert, "orco mesh");
-
-		for(a=0; a<totvert; a++, fp+=3) {
-			vcos[a][0]= fp[0];
-			vcos[a][1]= fp[1];
-			vcos[a][2]= fp[2];
-		}
+		vcos= mesh_getRefKeyCos(me, &totvert);
 	}
 	else {
 		MultiresLevel *lvl = NULL;
@@ -1120,9 +1111,8 @@ float (*mesh_getVertexCos(Mesh *me, int *numVerts_r))[3]
 		float (*cos)[3] = MEM_mallocN(sizeof(*cos)*numVerts, "vertexcos1");
         
 		if (numVerts_r) *numVerts_r = numVerts;
-		for (i=0; i<numVerts; i++) {
+		for (i=0; i<numVerts; i++)
 			VECCOPY(cos[i], me->mvert[i].co);
-		}
         
 		return cos;
 #ifdef WITH_VERSE
@@ -1130,6 +1120,25 @@ float (*mesh_getVertexCos(Mesh *me, int *numVerts_r))[3]
 #endif
 }
 
+float (*mesh_getRefKeyCos(Mesh *me, int *numVerts_r))[3]
+{
+	KeyBlock *kb;
+	float (*cos)[3] = NULL;
+	int totvert;
+	
+	if(me->key && me->key->refkey) {
+		if(numVerts_r) *numVerts_r= me->totvert;
+		cos= MEM_mallocN(sizeof(*cos)*me->totvert, "vertexcos1");
+
+		kb= me->key->refkey;
+		totvert= MIN2(kb->totelem, me->totvert);
+
+		memcpy(cos, kb->data, sizeof(*cos)*totvert);
+	}
+
+	return cos;
+}
+
 UvVertMap *make_uv_vert_map(struct MFace *mface, struct MTFace *tface, unsigned int totface, unsigned int totvert, int selected, float *limit)
 {
 	UvVertMap *vmap;
diff --git a/source/blender/blenkernel/intern/modifier.c b/source/blender/blenkernel/intern/modifier.c
index 98ba46d49d5..243a24bc3b8 100644
--- a/source/blender/blenkernel/intern/modifier.c
+++ b/source/blender/blenkernel/intern/modifier.c
@@ -47,6 +47,7 @@
 #include "BLI_linklist.h"
 #include "BLI_edgehash.h"
 #include "BLI_ghash.h"
+#include "BLI_memarena.h"
 
 #include "MEM_guardedalloc.h"
 
@@ -4921,6 +4922,223 @@ static DerivedMesh *booleanModifier_applyModifier(
 	return derivedData;
 }
 
+/* MeshDeform */
+
+static void meshdeformModifier_initData(ModifierData *md)
+{
+	MeshDeformModifierData *mmd = (MeshDeformModifierData*) md;
+
+	mmd->gridsize= 5;
+}
+
+static void meshdeformModifier_freeData(ModifierData *md)
+{
+	MeshDeformModifierData *mmd = (MeshDeformModifierData*) md;
+
+	if (mmd->bindweights) MEM_freeN(mmd->bindweights);
+	if (mmd->bindcos) MEM_freeN(mmd->bindcos);
+}
+
+static void meshdeformModifier_copyData(ModifierData *md, ModifierData *target)
+{
+	MeshDeformModifierData *mmd = (MeshDeformModifierData*) md;
+	MeshDeformModifierData *tmmd = (MeshDeformModifierData*) target;
+
+	tmmd->gridsize = mmd->gridsize;
+	tmmd->object = mmd->object;
+}
+
+CustomDataMask meshdeformModifier_requiredDataMask(ModifierData *md)
+{	
+	MeshDeformModifierData *mmd = (MeshDeformModifierData *)md;
+	CustomDataMask dataMask = 0;
+
+	/* ask for vertexgroups if we need them */
+	if(mmd->defgrp_name[0]) dataMask |= (1 << CD_MDEFORMVERT);
+
+	return dataMask;
+}
+
+static int meshdeformModifier_isDisabled(ModifierData *md)
+{
+	MeshDeformModifierData *mmd = (MeshDeformModifierData*) md;
+
+	return !mmd->object;
+}
+
+static void meshdeformModifier_foreachObjectLink(
+                ModifierData *md, Object *ob,
+                void (*walk)(void *userData, Object *ob, Object **obpoin),
+                void *userData)
+{
+	MeshDeformModifierData *mmd = (MeshDeformModifierData*) md;
+
+	walk(userData, ob, &mmd->object);
+}
+
+static void meshdeformModifier_updateDepgraph(
+                ModifierData *md, DagForest *forest, Object *ob,
+                DagNode *obNode)
+{
+	MeshDeformModifierData *mmd = (MeshDeformModifierData*) md;
+
+	if (mmd->object) {
+		DagNode *curNode = dag_get_node(forest, mmd->object);
+
+		dag_add_relation(forest, curNode, obNode,
+		                 DAG_RL_DATA_DATA|DAG_RL_OB_DATA|DAG_RL_DATA_OB|DAG_RL_OB_OB);
+	}
+}
+
+static void meshdeformModifier_do(
+                ModifierData *md, Object *ob, DerivedMesh *dm,
+                float (*vertexCos)[3], int numVerts)
+{
+	MeshDeformModifierData *mmd = (MeshDeformModifierData*) md;
+	float imat[4][4], cagemat[4][4], icagemat[4][4], icmat[3][3];
+	float weight, totweight, fac, co[3], *weights, (*dco)[3], (*bindcos)[3];
+	int a, b, totvert, totcagevert, defgrp_index;
+	DerivedMesh *tmpdm, *cagedm;
+	MDeformVert *dvert = NULL;
+	MDeformWeight *dw;
+	MVert *cagemvert;
+
+	if(!mmd->object || (!mmd->bindweights && !mmd->needbind))
+		return;
+	
+	/* get cage derivedmesh */
+	if(mmd->object == G.obedit) {
+		tmpdm= editmesh_get_derived_cage_and_final(&cagedm, 0);
+		if(tmpdm)
+			tmpdm->release(tmpdm);
+	}
+	else
+		cagedm= mesh_get_derived_final(mmd->object, CD_MASK_BAREMESH);
+	
+	/* TODO: this could give inifinite loop for circular dependency */
+	if(!cagedm)
+		return;
+
+ 	/* compute matrices to go in and out of cage object space */
+	Mat4Invert(imat, mmd->object->obmat);
+	Mat4MulMat4(cagemat, ob->obmat, imat);
+	Mat4Invert(icagemat, cagemat);
+	Mat3CpyMat4(icmat, icagemat);
+
+	/* bind weights if needed */
+	if(!mmd->bindweights)
+		harmonic_coordinates_bind(mmd, vertexCos, numVerts, cagemat);
+
+	/* verify we have compatible weights */
+	totvert= numVerts;
+	totcagevert= cagedm->getNumVerts(cagedm);
+
+	if(mmd->totvert!=totvert || mmd->totcagevert!=totcagevert || !mmd->bindweights) {
+		cagedm->release(cagedm);
+		return;
+	}
+	
+	/* setup deformation data */
+	cagemvert= cagedm->getVertArray(cagedm);
+	weights= mmd->bindweights;
+	bindcos= (float(*)[3])mmd->bindcos;
+
+	dco= MEM_callocN(sizeof(*dco)*totcagevert, "MDefDco");
+	for(a=0; a<totcagevert; a++) {
+		VECCOPY(co, cagemvert[a].co);
+		Mat4MulVecfl(mmd->object->obmat, co);
+		VECSUB(dco[a], co, bindcos[a]);
+	}
+
+	defgrp_index = -1;
+
+	if(mmd->defgrp_name[0]) {
+		bDeformGroup *def;
+
+		for(a=0, def=ob->defbase.first; def; def=def->next, a++) {
+			if(!strcmp(def->name, mmd->defgrp_name)) {
+				defgrp_index= a;
+				break;
+			}
+		}
+
+		if (defgrp_index >= 0)
+			dvert= dm->getVertDataArray(dm, CD_MDEFORMVERT);
+	}
+
+	/* do deformation */
+	for(b=0; b<totvert; b++) {
+		totweight= 0.0f;
+		co[0]= co[1]= co[2]= 0.0f;
+
+		for(a=0; a<totcagevert; a++) {
+			weight= weights[a + b*totcagevert];
+			co[0]+= weight*dco[a][0];
+			co[1]+= weight*dco[a][1];
+			co[2]+= weight*dco[a][2];
+			totweight += weight;
+		}
+
+		if(totweight > 0.0f) {
+			if(dvert) {
+				for(dw=NULL, a=0; a<dvert[b].totweight; a++) {
+					if(dvert[b].dw[a].def_nr == defgrp_index) {
+						dw = &dvert[b].dw[a];
+						break;
+					}
+				}
+				if(!dw) continue;
+
+				fac= dw->weight;
+			}
+			else
+				fac= 1.0f;
+
+			VecMulf(co, fac/totweight);
+			Mat3MulVecfl(icmat, co);
+			VECADD(vertexCos[b], vertexCos[b], co);
+		}
+	}
+
+	/* release cage derivedmesh */
+	MEM_freeN(dco);
+	cagedm->release(cagedm);
+}
+
+static void meshdeformModifier_deformVerts(
+                ModifierData *md, Object *ob, DerivedMesh *derivedData,
+                float (*vertexCos)[3], int numVerts)
+{
+	DerivedMesh *dm;
+
+	if(derivedData) dm = CDDM_copy(derivedData);
+	else dm = CDDM_from_mesh(ob->data, ob);
+
+	CDDM_apply_vert_coords(dm, vertexCos);
+	CDDM_calc_normals(dm);
+
+	meshdeformModifier_do(md, ob, dm, vertexCos, numVerts);
+
+	dm->release(dm);
+}
+
+static void meshdeformModifier_deformVertsEM(
+                ModifierData *md, Object *ob, EditMesh *editData,
+                DerivedMesh *derivedData, float (*vertexCos)[3], int numVerts)
+{
+	DerivedMesh *dm;
+
+	if(derivedData) dm = CDDM_copy(derivedData);
+	else dm = CDDM_from_editmesh(editData, ob->data);
+
+	CDDM_apply_vert_coords(dm, vertexCos);
+	CDDM_calc_normals(dm);
+
+	meshdeformModifier_do(md, ob, dm, vertexCos, numVerts);
+
+	dm->release(dm);
+}
+
 /***/
 
 static ModifierTypeInfo typeArr[NUM_MODIFIER_TYPES];
@@ -5149,6 +5367,20 @@ ModifierTypeInfo *modifierType_getInfo(ModifierType type)
 		mti->foreachObjectLink = booleanModifier_foreachObjectLink;
 		mti->updateDepgraph = booleanModifier_updateDepgraph;
 
+		mti = INIT_TYPE(MeshDeform);
+		mti->type = eModifierTypeType_OnlyDeform;
+		mti->flags = eModifierTypeFlag_AcceptsCVs
+		             | eModifierTypeFlag_SupportsEditmode;
+		mti->initData = meshdeformModifier_initData;
+		mti->freeData = meshdeformModifier_freeData;
+		mti->copyData = meshdeformModifier_copyData;
+		mti->requiredDataMask = meshdeformModifier_requiredDataMask;
+		mti->isDisabled = meshdeformModifier_isDisabled;
+		mti->foreachObjectLink = meshdeformModifier_foreachObjectLink;
+		mti->updateDepgraph = meshdeformModifier_updateDepgraph;
+		mti->deformVerts = meshdeformModifier_deformVerts;
+		mti->deformVertsEM = meshdeformModifier_deformVertsEM;
+
 		typeArrInit = 0;
 #undef INIT_TYPE
 	}
@@ -5523,3 +5755,4 @@ int modifiers_isDeformed(Object *ob)
 	}
 	return 0;
 }
+
diff --git a/source/blender/blenlib/BLI_arithb.h b/source/blender/blenlib/BLI_arithb.h
index fc132250fc6..2942439504c 100644
--- a/source/blender/blenlib/BLI_arithb.h
+++ b/source/blender/blenlib/BLI_arithb.h
@@ -282,7 +282,8 @@ extern short IsectLL2Ds(short *v1, short *v2, short *v3, short *v4);
 
 /* interpolation weights of point in a triangle or quad, v4 may be NULL */
 void InterpWeightsQ3Dfl(float *v1, float *v2, float *v3, float *v4, float *co, float *w);
-
+/* interpolation weights of point in a polygon with >= 3 vertices */
+void MeanValueWeights(float v[][3], int n, float *co, float *w);
 
 void i_lookat(
 	float vx, float vy, 
diff --git a/source/blender/blenlib/intern/arithb.c b/source/blender/blenlib/intern/arithb.c
index f95d102763a..721df3a1a0c 100644
--- a/source/blender/blenlib/intern/arithb.c
+++ b/source/blender/blenlib/intern/arithb.c
@@ -2506,7 +2506,53 @@ void InterpWeightsQ3Dfl(float *v1, float *v2, float *v3, float *v4, float *co, f
 		else
 			BarycentricWeights(v1, v2, v3, co, n, w);
 	}
-} 
+}
+
+/* Mean value weights - smooth interpolation weights for polygons with
+ * more than 3 vertices */
+static float MeanValueHalfTan(float *v1, float *v2, float *v3)
+{
+	float d2[3], d3[3], cross[3], area, dot, len;
+
+	VecSubf(d2, v2, v1);
+	VecSubf(d3, v3, v1);
+	Crossf(cross, d2, d3);
+
+	area= VecLength(cross);
+	dot= Inpf(d2, d3);
+	len= VecLength(d2)*VecLength(d3);
+
+	if(area == 0.0f)
+		return 0.0f;
+	else
+		return (len - dot)/area;
+}
+
+void MeanValueWeights(float v[][3], int n, float *co, float *w)
+{
+	float totweight, t1, t2, len, *vmid, *vprev, *vnext;
+	int i;
+
+	totweight= 0.0f;
+
+	for(i=0; i<n; i++) {
+		vmid= v[i];
+		vprev= (i == 0)? v[n-1]: v[i-1];
+		vnext= (i == n-1)? v[0]: v[i+1];
+
+		t1= MeanValueHalfTan(co, vprev, vmid);
+		t2= MeanValueHalfTan(co, vmid, vnext);
+
+		len= VecLenf(co, vmid);
+		w[i]= (t1+t2)/len;
+		totweight += w[i];
+	}
+
+	if(totweight != 0.0f)
+		for(i=0; i<n; i++)
+			w[i] /= totweight;
+}
+
 
 /* ************ EULER *************** */
 
diff --git a/source/blender/blenloader/intern/readfile.c b/source/blender/blenloader/intern/readfile.c
index 82923ca1aa5..384a6d93480 100644
--- a/source/blender/blenloader/intern/readfile.c
+++ b/source/blender/blenloader/intern/readfile.c
@@ -2889,6 +2889,21 @@ static void direct_link_modifiers(FileData *fd, ListBase *lb)
 				}
 			}
 		}
+		else if (md->type==eModifierType_MeshDeform) {
+			MeshDeformModifierData *mmd = (MeshDeformModifierData*) md;
+
+			mmd->bindweights= newdataadr(fd, mmd->bindweights);
+			mmd->bindcos= newdataadr(fd, mmd->bindcos);
+
+			if(fd->flags & FD_FLAGS_SWITCH_ENDIAN) {
+				int a;
+
+				for(a=0; a<mmd->totcagevert*mmd->totvert; a++)
+					SWITCH_INT(mmd->bindweights[a])
+				for(a=0; a<mmd->totcagevert*3; a++)
+					SWITCH_INT(mmd->bindcos[a])
+			}
+		}
 	}
 }
 
diff --git a/source/blender/blenloader/intern/writefile.c b/source/blender/blenloader/intern/writefile.c
index 45621c5fcd0..43907a30ac2 100644
--- a/source/blender/blenloader/intern/writefile.c
+++ b/source/blender/blenloader/intern/writefile.c
@@ -785,6 +785,14 @@ static void write_modifiers(WriteData *wd, ListBase *modbase)
 			
 			writedata(wd, DATA, sizeof(int)*hmd->totindex, hmd->indexar);
 		}
+		else if (md->type==eModifierType_MeshDeform) {
+			MeshDeformModifierData *mmd = (MeshDeformModifierData*) md;
+
+			writedata(wd, DATA, sizeof(float)*mmd->totvert*mmd->totcagevert,
+				mmd->bindweights);
+			writedata(wd, DATA, sizeof(float)*3*mmd->totcagevert,
+				mmd->bindcos);
+		}
 	}
 }
 
diff --git a/source/blender/include/BIF_meshlaplacian.h b/source/blender/include/BIF_meshlaplacian.h
index 1e10336723f..74e4fef0937 100644
--- a/source/blender/include/BIF_meshlaplacian.h
+++ b/source/blender/include/BIF_meshlaplacian.h
@@ -39,6 +39,7 @@
 struct Object;
 struct Mesh;
 struct bDeformGroup;
+struct MeshDeformModifierData;
 
 #ifdef RIGID_DEFORM
 struct EditMesh;
@@ -77,5 +78,10 @@ void rigid_deform_iteration(void);
 void rigid_deform_end(int cancel);
 #endif
 
+/* Harmonic Coordinates */
+
+void harmonic_coordinates_bind(struct MeshDeformModifierData *mmd,
+	float (*vertexcos)[3], int totvert, float cagemat[][4]);
+
 #endif
 
diff --git a/source/blender/makesdna/DNA_modifier_types.h b/source/blender/makesdna/DNA_modifier_types.h
index 71e850e4368..c8dcba4fae7 100644
--- a/source/blender/makesdna/DNA_modifier_types.h
+++ b/source/blender/makesdna/DNA_modifier_types.h
@@ -28,6 +28,7 @@ typedef enum ModifierType {
 	eModifierType_UVProject,
 	eModifierType_Smooth,
 	eModifierType_Cast,
+	eModifierType_MeshDeform,
 	NUM_MODIFIER_TYPES
 } ModifierType;
 
@@ -347,4 +348,17 @@ typedef struct BooleanModifierData {
 	int operation, pad;
 } BooleanModifierData;
 
+typedef struct MeshDeformModifierData {
+	ModifierData modifier;
+
+	struct Object *object;			/* mesh object */
+	char defgrp_name[32];			/* optional vertexgroup name */
+
+	float *bindweights, *bindcos;	/* computed binding weights */
+	short gridsize, needbind;
+	int pad;
+
+	int totvert, totcagevert;
+} MeshDeformModifierData;
+
 #endif
diff --git a/source/blender/src/buttons_editing.c b/source/blender/src/buttons_editing.c
index a9a75126f51..40e2643a390 100644
--- a/source/blender/src/buttons_editing.c
+++ b/source/blender/src/buttons_editing.c
@@ -1465,6 +1465,35 @@ void set_uvproject_uvlayer(void *arg1, void *arg2)
 	strcpy(umd->uvlayer_name, layer->name);
 }
 
+static void modifiers_bindMeshDeform(void *ob_v, void *md_v)
+{
+	MeshDeformModifierData *mmd = (MeshDeformModifierData*) md_v;
+	Object *ob = (Object*)ob_v;
+
+	if(mmd->bindweights) {
+		MEM_freeN(mmd->bindweights);
+		MEM_freeN(mmd->bindcos);
+		mmd->bindweights= NULL;
+		mmd->bindcos= NULL;
+		mmd->totvert= 0;
+		mmd->totcagevert= 0;
+	}
+	else {
+		DerivedMesh *dm;
+		int mode= mmd->modifier.mode;
+
+		/* force modifier to run, it will call binding routine */
+		mmd->needbind= 1;
+		mmd->modifier.mode |= eModifierMode_Realtime;
+
+		dm= mesh_create_derived_view(ob, 0);
+		dm->release(dm);
+
+		mmd->needbind= 0;
+		mmd->modifier.mode= mode;
+	}
+}
+
 static void draw_modifier(uiBlock *block, Object *ob, ModifierData *md, int *xco, int *yco, int index, int cageIndex, int lastCageIndex)
 {
 	ModifierTypeInfo *mti = modifierType_getInfo(md->type);
@@ -1605,6 +1634,8 @@ static void draw_modifier(uiBlock *block, Object *ob, ModifierData *md, int *xco
 			height = 48;
 		} else if (md->type==eModifierType_Array) {
 			height = 211;
+		} else if (md->type==eModifierType_MeshDeform) {
+			height = 73;
 		} 
 		
 							/* roundbox 4 free variables: corner-rounding, nop, roundbox type, shade */
@@ -2084,7 +2115,27 @@ static void draw_modifier(uiBlock *block, Object *ob, ModifierData *md, int *xco
 			                     &amd->end_cap,
 			                     "Mesh object to use as end cap");
 			uiButSetCompleteFunc(but, autocomplete_meshob, (void *)ob);
+		} else if (md->type==eModifierType_MeshDeform) {
+			MeshDeformModifierData *mmd = (MeshDeformModifierData*) md;
+
+			uiBlockBeginAlign(block);
+			uiDefIDPoinBut(block, test_meshobpoin_but, ID_OB, B_CHANGEDEP, "Ob: ", lx, (cy-=19), buttonWidth,19, &mmd->object, "Mesh object to be use as cage"); 
+			but=uiDefBut(block, TEX, B_MODIFIER_RECALC, "VGroup: ",				  lx, (cy-=19), buttonWidth,19, &mmd->defgrp_name, 0.0, 31.0, 0, 0, "Vertex Group name to control overall meshdeform influence");
+			uiButSetCompleteFunc(but, autocomplete_vgroup, (void *)ob);
+
+			uiBlockBeginAlign(block);
+			if(mmd->bindweights) {
+				but= uiDefBut(block, BUT, B_MODIFIER_RECALC, "Unbind", lx,(cy-24), buttonWidth,19, 0, 0, 0, 0, 0, "Unbind mesh from cage");
+				uiButSetFunc(but,modifiers_bindMeshDeform,ob,md);
+			}
+			else {
+				but= uiDefBut(block, BUT, B_MODIFIER_RECALC, "Bind", lx,(cy-24), buttonWidth/2,19, 0, 0, 0, 0, 0, "Bind mesh to cage");
+				uiButSetFunc(but,modifiers_bindMeshDeform,ob,md);
+				uiDefButS(block, NUM, B_NOP, "Precision:", lx+(buttonWidth+1)/2,(cy-=24), buttonWidth/2,19, &mmd->gridsize, 2, 10, 0.5, 0, "The grid size for binding");
+			}
+			uiBlockEndAlign(block);
 		}
+
 		uiBlockEndAlign(block);
 
 		y-=height;
diff --git a/source/blender/src/meshlaplacian.c b/source/blender/src/meshlaplacian.c
index a660b42996e..c27694ae51e 100644
--- a/source/blender/src/meshlaplacian.c
+++ b/source/blender/src/meshlaplacian.c
@@ -31,6 +31,7 @@
  * meshlaplacian.c: Algorithms using the mesh laplacian.
  */
 
+#include <math.h>
 #include <string.h>
 
 #include "MEM_guardedalloc.h"
@@ -39,17 +40,23 @@
 #include "DNA_object_types.h"
 #include "DNA_mesh_types.h"
 #include "DNA_meshdata_types.h"
+#include "DNA_modifier_types.h"
 
 #include "BLI_arithb.h"
 #include "BLI_edgehash.h"
+#include "BLI_memarena.h"
 
+#include "BKE_DerivedMesh.h"
 #include "BKE_utildefines.h"
 
 #include "BIF_editdeform.h"
 #include "BIF_meshlaplacian.h"
 #include "BIF_meshtools.h"
+#include "BIF_screen.h"
 #include "BIF_toolbox.h"
 
+#include "BSE_headerbuttons.h"
+
 #ifdef RIGID_DEFORM
 #include "BLI_editVert.h"
 #include "BLI_polardecomp.h"
@@ -336,7 +343,7 @@ void laplacian_begin_solve(LaplacianSystem *sys, int index)
 		if(index >= 0) {
 			for(a=0; a<sys->totvert; a++) {
 				if(sys->vpinned[a]) {
-					nlSetVariable(a, sys->verts[a][index]);
+					nlSetVariable(0, a, sys->verts[a][index]);
 					nlLockVariable(a);
 				}
 			}
@@ -349,7 +356,7 @@ void laplacian_begin_solve(LaplacianSystem *sys, int index)
 
 void laplacian_add_right_hand_side(LaplacianSystem *sys, int v, float value)
 {
-	nlRightHandSideAdd(v, value);
+	nlRightHandSideAdd(0, v, value);
 }
 
 int laplacian_system_solve(LaplacianSystem *sys)
@@ -365,7 +372,7 @@ int laplacian_system_solve(LaplacianSystem *sys)
 
 float laplacian_system_get_solution(int v)
 {
-	return nlGetVariable(v);
+	return nlGetVariable(0, v);
 }
 
 /************************* Heat Bone Weighting ******************************/
@@ -453,6 +460,7 @@ static int heat_ray_bone_visible(LaplacianSystem *sys, int vertex, int bone)
 		return 1;
 
 	/* setup isec */
+	memset(&isec, 0, sizeof(isec));
 	isec.mode= RE_RAY_SHADOW;
 	isec.lay= -1;
 	isec.face_last= NULL;
@@ -904,3 +912,866 @@ void rigid_deform_end(int cancel)
 }
 #endif
 
+/************************** Harmonic Coordinates ****************************/
+/* From "Harmonic Coordinates for Character Articulation",
+	Pushkar Joshi, Mark Meyer, Tony DeRose, Brian Green and Tom Sanocki,
+	SIGGRAPH 2007. */
+
+#define EPSILON 0.0001f
+
+#define MESHDEFORM_TAG_UNTYPED  0
+#define MESHDEFORM_TAG_BOUNDARY 1
+#define MESHDEFORM_TAG_INTERIOR 2
+#define MESHDEFORM_TAG_EXTERIOR 3
+
+#define MESHDEFORM_LEN_THRESHOLD 1e-6
+
+static int MESHDEFORM_OFFSET[7][3] =
+		{{0,0,0}, {1,0,0}, {-1,0,0}, {0,1,0}, {0,-1,0}, {0,0,1}, {0,0,-1}};
+
+typedef struct MDefBoundIsect {
+	float co[3], uvw[4];
+	int nvert, v[4], facing;
+	float len;
+} MDefBoundIsect;
+
+typedef struct MeshDeformBind {
+	/* grid dimensions */
+	float min[3], max[3];
+	float width[3], halfwidth[3];
+	int size, size3;
+
+	/* meshes */
+	DerivedMesh *cagedm;
+	float (*cagecos)[3];
+	float (*vertexcos)[3];
+	int totvert, totcagevert;
+
+	/* grids */
+	MemArena *memarena;
+	MDefBoundIsect *(*boundisect)[6];
+	int *semibound;
+	int *tag;
+	float *phi, *totalphi;
+
+	/* mesh stuff */
+	int *inside;
+	float *weights;
+	float cagemat[4][4];
+
+	/* direct solver */
+	int *varidx;
+
+	/* raytrace */
+	RayTree *raytree;
+} MeshDeformBind;
+
+/* ray intersection */
+
+/* our own triangle intersection, so we can fully control the epsilons and
+ * prevent corner case from going wrong*/
+static int meshdeform_tri_intersect(float orig[3], float end[3], float vert0[3],
+    float vert1[3], float vert2[3], float *isectco, float *uvw)
+{
+	float edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
+	float det,inv_det, u, v, dir[3], isectdir[3];
+
+	VECSUB(dir, end, orig);
+
+	/* find vectors for two edges sharing vert0 */
+	VECSUB(edge1, vert1, vert0);
+	VECSUB(edge2, vert2, vert0);
+
+	/* begin calculating determinant - also used to calculate U parameter */
+	Crossf(pvec, dir, edge2);
+
+	/* if determinant is near zero, ray lies in plane of triangle */
+	det = INPR(edge1, pvec);
+
+	if (det == 0.0f)
+	  return 0;
+	inv_det = 1.0f / det;
+
+	/* calculate distance from vert0 to ray origin */
+	VECSUB(tvec, orig, vert0);
+
+	/* calculate U parameter and test bounds */
+	u = INPR(tvec, pvec) * inv_det;
+	if (u < -EPSILON || u > 1.0f+EPSILON)
+	  return 0;
+
+	/* prepare to test V parameter */
+	Crossf(qvec, tvec, edge1);
+
+	/* calculate V parameter and test bounds */
+	v = INPR(dir, qvec) * inv_det;
+	if (v < -EPSILON || u + v > 1.0f+EPSILON)
+	  return 0;
+
+	isectco[0]= (1.0f - u - v)*vert0[0] + u*vert1[0] + v*vert2[0];
+	isectco[1]= (1.0f - u - v)*vert0[1] + u*vert1[1] + v*vert2[1];
+	isectco[2]= (1.0f - u - v)*vert0[2] + u*vert1[2] + v*vert2[2];
+
+	uvw[0]= 1.0 - u - v;
+	uvw[1]= u;
+	uvw[2]= v;
+
+	/* check if it is within the length of the line segment */
+	VECSUB(isectdir, isectco, orig);
+
+	if(INPR(dir, isectdir) < -EPSILON)
+		return 0;
+	
+	if(INPR(dir, dir) + EPSILON < INPR(isectdir, isectdir))
+		return 0;
+
+	return 1;
+}
+
+/* blender's raytracer is not use now, even though it is much faster. it can
+ * give problems with rays falling through, so we use our own intersection 
+ * function above with tweaked epsilons */
+
+#if 0
+static MeshDeformBind *MESHDEFORM_BIND = NULL;
+
+static void meshdeform_ray_coords_func(RayFace *face, float **v1, float **v2, float **v3, float **v4)
+{
+	MFace *mface= (MFace*)face;
+	float (*cagecos)[3]= MESHDEFORM_BIND->cagecos;
+
+	*v1= cagecos[mface->v1];
+	*v2= cagecos[mface->v2];
+	*v3= cagecos[mface->v3];
+	*v4= (mface->v4)? cagecos[mface->v4]: NULL;
+}
+
+static int meshdeform_ray_check_func(Isect *is, RayFace *face)
+{
+	return 1;
+}
+
+static void meshdeform_ray_tree_create(MeshDeformBind *mdb)
+{
+	MFace *mface;
+	float min[3], max[3];
+	int a, totface;
+
+	/* create a raytrace tree from the mesh */
+	INIT_MINMAX(min, max);
+
+	for(a=0; a<mdb->totcagevert; a++)
+		DO_MINMAX(mdb->cagecos[a], min, max)
+
+	MESHDEFORM_BIND= mdb;
+
+	mface= mdb->cagedm->getFaceArray(mdb->cagedm);
+	totface= mdb->cagedm->getNumFaces(mdb->cagedm);
+
+	mdb->raytree= RE_ray_tree_create(64, totface, min, max,
+		meshdeform_ray_coords_func, meshdeform_ray_check_func);
+
+	for(a=0; a<totface; a++, mface++)
+		RE_ray_tree_add_face(mdb->raytree, mface);
+
+	RE_ray_tree_done(mdb->raytree);
+}
+
+static void meshdeform_ray_tree_free(MeshDeformBind *mdb)
+{
+	MESHDEFORM_BIND= NULL;
+	RE_ray_tree_free(mdb->raytree);
+}
+#endif
+
+static int meshdeform_intersect(MeshDeformBind *mdb, Isect *isec)
+{
+	MFace *mface;
+	float face[4][3], co[3], uvw[3], len, nor[3];
+	int f, hit, is= 0, totface;
+
+	isec->labda= 1e10;
+
+	mface= mdb->cagedm->getFaceArray(mdb->cagedm);
+	totface= mdb->cagedm->getNumFaces(mdb->cagedm);
+
+	for(f=0; f<totface; f++, mface++) {
+		VECCOPY(face[0], mdb->cagecos[mface->v1]);
+		VECCOPY(face[1], mdb->cagecos[mface->v2]);
+		VECCOPY(face[2], mdb->cagecos[mface->v3]);
+
+		if(mface->v4) {
+			VECCOPY(face[3], mdb->cagecos[mface->v4]);
+			hit= meshdeform_tri_intersect(isec->start, isec->end, face[0], face[1], face[2], co, uvw);
+
+			if(hit) {
+				CalcNormFloat(face[0], face[1], face[2], nor);
+			}
+			else {
+				hit= meshdeform_tri_intersect(isec->start, isec->end, face[0], face[2], face[3], co, uvw);
+				CalcNormFloat(face[0], face[2], face[3], nor);
+			}
+		}
+		else {
+			hit= meshdeform_tri_intersect(isec->start, isec->end, face[0], face[1], face[2], co, uvw);
+			CalcNormFloat(face[0], face[1], face[2], nor);
+		}
+
+		if(hit) {
+			len= VecLenf(isec->start, co)/VecLenf(isec->start, isec->end);
+			if(len < isec->labda) {
+				isec->labda= len;
+				isec->face= mface;
+				isec->isect= (INPR(isec->vec, nor) <= 0.0f);
+				is= 1;
+			}
+		}
+	}
+
+	return is;
+}
+
+static MDefBoundIsect *meshdeform_ray_tree_intersect(MeshDeformBind *mdb, float *co1, float *co2)
+{
+	MDefBoundIsect *isect;
+	Isect isec;
+	float (*cagecos)[3];
+	MFace *mface;
+	float vert[4][3], len;
+	static float epsilon[3]= {0, 0, 0}; //1e-4, 1e-4, 1e-4};
+
+	/* setup isec */
+	memset(&isec, 0, sizeof(isec));
+	isec.mode= RE_RAY_MIRROR; /* we want the closest intersection */
+	isec.lay= -1;
+	isec.face_last= NULL;
+	isec.faceorig= NULL;
+	isec.labda= 1e10f;
+
+	VECADD(isec.start, co1, epsilon);
+	VECADD(isec.end, co2, epsilon);
+	VECSUB(isec.vec, isec.end, isec.start);
+
+#if 0
+	/*if(RE_ray_tree_intersect(mdb->raytree, &isec)) {*/
+#endif
+
+	if(meshdeform_intersect(mdb, &isec)) {
+		len= isec.labda;
+		mface= isec.face;
+
+		/* create MDefBoundIsect */
+		isect= BLI_memarena_alloc(mdb->memarena, sizeof(*isect));
+
+		/* compute intersection coordinate */
+		isect->co[0]= co1[0] + isec.vec[0]*len;
+		isect->co[1]= co1[1] + isec.vec[1]*len;
+		isect->co[2]= co1[2] + isec.vec[2]*len;
+
+		isect->len= VecLenf(co1, isect->co);
+		if(isect->len < MESHDEFORM_LEN_THRESHOLD)
+			isect->len= MESHDEFORM_LEN_THRESHOLD;
+
+		isect->v[0]= mface->v1;
+		isect->v[1]= mface->v2;
+		isect->v[2]= mface->v3;
+		isect->v[3]= mface->v4;
+		isect->nvert= (mface->v4)? 4: 3;
+
+		isect->facing= isec.isect;
+
+		/* compute mean value coordinates for interpolation */
+		cagecos= mdb->cagecos;
+		VECCOPY(vert[0], cagecos[mface->v1]);
+		VECCOPY(vert[1], cagecos[mface->v2]);
+		VECCOPY(vert[2], cagecos[mface->v3]);
+		if(mface->v4) VECCOPY(vert[3], cagecos[mface->v4]);
+		MeanValueWeights(vert, isect->nvert, isect->co, isect->uvw);
+
+		return isect;
+	}
+
+	return NULL;
+}
+
+static int meshdeform_inside_cage(MeshDeformBind *mdb, float *co)
+{
+	MDefBoundIsect *isect;
+	float outside[3], start[3], dir[3];
+	int i, counter;
+
+	for(i=1; i<=6; i++) {
+		counter = 0;
+
+		outside[0] = co[0] + (mdb->max[0] - mdb->min[0] + 1.0f)*MESHDEFORM_OFFSET[i][0];
+		outside[1] = co[1] + (mdb->max[1] - mdb->min[1] + 1.0f)*MESHDEFORM_OFFSET[i][1];
+		outside[2] = co[2] + (mdb->max[2] - mdb->min[2] + 1.0f)*MESHDEFORM_OFFSET[i][2];
+
+		VECSUB(dir, outside, start);
+		Normalize(dir);
+		VECCOPY(start, co);
+		
+		isect = meshdeform_ray_tree_intersect(mdb, start, outside);
+		if(isect && !isect->facing)
+			return 1;
+	}
+
+	return 0;
+}
+
+/* solving */
+
+static int meshdeform_index(MeshDeformBind *mdb, int x, int y, int z, int n)
+{
+	int size= mdb->size;
+	
+	x += MESHDEFORM_OFFSET[n][0];
+	y += MESHDEFORM_OFFSET[n][1];
+	z += MESHDEFORM_OFFSET[n][2];
+
+	if(x < 0 || x >= mdb->size)
+		return -1;
+	if(y < 0 || y >= mdb->size)
+		return -1;
+	if(z < 0 || z >= mdb->size)
+		return -1;
+
+	return x + y*size + z*size*size;
+}
+
+static void meshdeform_cell_center(MeshDeformBind *mdb, int x, int y, int z, int n, float *center)
+{
+	x += MESHDEFORM_OFFSET[n][0];
+	y += MESHDEFORM_OFFSET[n][1];
+	z += MESHDEFORM_OFFSET[n][2];
+
+	center[0]= mdb->min[0] + x*mdb->width[0] + mdb->halfwidth[0];
+	center[1]= mdb->min[1] + y*mdb->width[1] + mdb->halfwidth[1];
+	center[2]= mdb->min[2] + z*mdb->width[2] + mdb->halfwidth[2];
+}
+
+static void meshdeform_add_intersections(MeshDeformBind *mdb, int x, int y, int z)
+{
+	MDefBoundIsect *isect;
+	float center[3], ncenter[3];
+	int i, a;
+
+	a= meshdeform_index(mdb, x, y, z, 0);
+	meshdeform_cell_center(mdb, x, y, z, 0, center);
+
+	/* check each outgoing edge for intersection */
+	for(i=1; i<=6; i++) {
+		if(meshdeform_index(mdb, x, y, z, i) == -1)
+			continue;
+
+		meshdeform_cell_center(mdb, x, y, z, i, ncenter);
+
+		isect= meshdeform_ray_tree_intersect(mdb, center, ncenter);
+		if(isect) {
+			mdb->boundisect[a][i-1]= isect;
+			mdb->tag[a]= MESHDEFORM_TAG_BOUNDARY;
+		}
+	}
+}
+
+static void meshdeform_bind_floodfill(MeshDeformBind *mdb)
+{
+	int *stack, *tag= mdb->tag;
+	int a, b, i, xyz[3], stacksize, size= mdb->size;
+
+	stack= MEM_callocN(sizeof(int)*mdb->size3, "MeshDeformBindStack");
+
+	/* we know lower left corner is EXTERIOR because of padding */
+	tag[0]= MESHDEFORM_TAG_EXTERIOR;
+	stack[0]= 0;
+	stacksize= 1;
+
+	/* floodfill exterior tag */
+	while(stacksize > 0) {
+		a= stack[--stacksize];
+
+		xyz[2]= a/(size*size);
+		xyz[1]= (a - xyz[2]*size*size)/size;
+		xyz[0]= a - xyz[1]*size - xyz[2]*size*size;
+
+		for(i=1; i<=6; i++) {
+			b= meshdeform_index(mdb, xyz[0], xyz[1], xyz[2], i);
+
+			if(b != -1) {
+				if(tag[b] == MESHDEFORM_TAG_UNTYPED ||
+				   (tag[b] == MESHDEFORM_TAG_BOUNDARY && !mdb->boundisect[a][i-1])) {
+					tag[b]= MESHDEFORM_TAG_EXTERIOR;
+					stack[stacksize++]= b;
+				}
+			}
+		}
+	}
+
+	/* other cells are interior */
+	for(a=0; a<size*size*size; a++)
+		if(tag[a]==MESHDEFORM_TAG_UNTYPED)
+			tag[a]= MESHDEFORM_TAG_INTERIOR;
+
+#if 0
+	{
+		int tb, ti, te, ts;
+		tb= ti= te= ts= 0;
+		for(a=0; a<size*size*size; a++)
+			if(tag[a]==MESHDEFORM_TAG_BOUNDARY)
+				tb++;
+			else if(tag[a]==MESHDEFORM_TAG_INTERIOR)
+				ti++;
+			else if(tag[a]==MESHDEFORM_TAG_EXTERIOR) {
+				te++;
+
+				if(mdb->semibound[a])
+					ts++;
+			}
+		
+		printf("interior %d exterior %d boundary %d semi-boundary %d\n", ti, te, tb, ts);
+	}
+#endif
+
+	MEM_freeN(stack);
+}
+
+static float meshdeform_boundary_phi(MeshDeformBind *mdb, MDefBoundIsect *isect, int cagevert)
+{
+	int a;
+
+	for(a=0; a<isect->nvert; a++)
+		if(isect->v[a] == cagevert)
+			return isect->uvw[a];
+	
+	return 0.0f;
+}
+
+static float meshdeform_interp_w(MeshDeformBind *mdb, float *gridvec, float *vec, int cagevert)
+{
+	float dvec[3], ivec[3], wx, wy, wz, result=0.0f;
+	float weight, totweight= 0.0f;
+	int i, a, x, y, z;
+
+	for(i=0; i<3; i++) {
+		ivec[i]= (int)gridvec[i];
+		dvec[i]= gridvec[i] - ivec[i];
+	}
+
+	for(i=0; i<8; i++) {
+		if(i & 1) { x= ivec[0]+1; wx= dvec[0]; }
+		else { x= ivec[0]; wx= 1.0f-dvec[0]; } 
+
+		if(i & 2) { y= ivec[1]+1; wy= dvec[1]; }
+		else { y= ivec[1]; wy= 1.0f-dvec[1]; } 
+
+		if(i & 4) { z= ivec[2]+1; wz= dvec[2]; }
+		else { z= ivec[2]; wz= 1.0f-dvec[2]; } 
+
+		CLAMP(x, 0, mdb->size-1);
+		CLAMP(y, 0, mdb->size-1);
+		CLAMP(z, 0, mdb->size-1);
+
+		a= meshdeform_index(mdb, x, y, z, 0);
+		weight= wx*wy*wz;
+		result += weight*mdb->phi[a];
+		totweight += weight;
+	}
+
+	if(totweight > 0.0f)
+		result /= totweight;
+
+	return result;
+}
+
+static void meshdeform_check_semibound(MeshDeformBind *mdb, int x, int y, int z)
+{
+	int i, a;
+
+	a= meshdeform_index(mdb, x, y, z, 0);
+	if(mdb->tag[a] != MESHDEFORM_TAG_EXTERIOR)
+		return;
+
+	for(i=1; i<=6; i++)
+		if(mdb->boundisect[a][i-1]) 
+			mdb->semibound[a]= 1;
+}
+
+static float meshdeform_boundary_total_weight(MeshDeformBind *mdb, int x, int y, int z)
+{
+	float weight, totweight= 0.0f;
+	int i, a;
+
+	a= meshdeform_index(mdb, x, y, z, 0);
+
+	/* count weight for neighbour cells */
+	for(i=1; i<=6; i++) {
+		if(meshdeform_index(mdb, x, y, z, i) == -1)
+			continue;
+
+		if(mdb->boundisect[a][i-1])
+			weight= 1.0f/mdb->boundisect[a][i-1]->len;
+		else if(!mdb->semibound[a])
+			weight= 1.0f/mdb->width[0];
+		else
+			weight= 0.0f;
+
+		totweight += weight;
+	}
+
+	return totweight;
+}
+
+static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, int x, int y, int z)
+{
+	MDefBoundIsect *isect;
+	float weight, totweight;
+	int i, a, acenter;
+
+	acenter= meshdeform_index(mdb, x, y, z, 0);
+	if(mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR)
+		return;
+
+	nlMatrixAdd(mdb->varidx[acenter], mdb->varidx[acenter], 1.0f);
+	
+	totweight= meshdeform_boundary_total_weight(mdb, x, y, z);
+	for(i=1; i<=6; i++) {
+		a= meshdeform_index(mdb, x, y, z, i);
+		if(a == -1 || mdb->tag[a] == MESHDEFORM_TAG_EXTERIOR)
+			continue;
+
+		isect= mdb->boundisect[acenter][i-1];
+		if (!isect) {
+			weight= (1.0f/mdb->width[0])/totweight;
+			nlMatrixAdd(mdb->varidx[acenter], mdb->varidx[a], -weight);
+		}
+	}
+}
+
+static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, int x, int y, int z, int cagevert)
+{
+	MDefBoundIsect *isect;
+	float rhs, weight, totweight;
+	int i, a, acenter;
+
+	acenter= meshdeform_index(mdb, x, y, z, 0);
+	if(mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR)
+		return;
+
+	totweight= meshdeform_boundary_total_weight(mdb, x, y, z);
+	for(i=1; i<=6; i++) {
+		a= meshdeform_index(mdb, x, y, z, i);
+		if(a == -1)
+			continue;
+
+		isect= mdb->boundisect[acenter][i-1];
+
+		if (isect) {
+			weight= (1.0f/isect->len)/totweight;
+			rhs= weight*meshdeform_boundary_phi(mdb, isect, cagevert);
+			nlRightHandSideAdd(0, mdb->varidx[acenter], rhs);
+		}
+	}
+}
+
+static void meshdeform_matrix_add_semibound_phi(MeshDeformBind *mdb, int x, int y, int z, int cagevert)
+{
+	MDefBoundIsect *isect;
+	float rhs, weight, totweight;
+	int i, a;
+
+	a= meshdeform_index(mdb, x, y, z, 0);
+	if(!mdb->semibound[a])
+		return;
+	
+	mdb->phi[a]= 0.0f;
+
+	totweight= meshdeform_boundary_total_weight(mdb, x, y, z);
+	for(i=1; i<=6; i++) {
+		isect= mdb->boundisect[a][i-1];
+
+		if (isect) {
+			weight= (1.0f/isect->len)/totweight;
+			rhs= weight*meshdeform_boundary_phi(mdb, isect, cagevert);
+			mdb->phi[a] += rhs;
+		}
+	}
+}
+
+static void meshdeform_matrix_add_exterior_phi(MeshDeformBind *mdb, int x, int y, int z, int cagevert)
+{
+	float phi, totweight;
+	int i, a, acenter;
+
+	acenter= meshdeform_index(mdb, x, y, z, 0);
+	if(mdb->tag[acenter] != MESHDEFORM_TAG_EXTERIOR || mdb->semibound[acenter])
+		return;
+
+	phi= 0.0f;
+	totweight= 0.0f;
+	for(i=1; i<=6; i++) {
+		a= meshdeform_index(mdb, x, y, z, i);
+
+		if(a != -1 && mdb->semibound[a]) {
+			phi += mdb->phi[a];
+			totweight += 1.0f;
+		}
+	}
+
+	if(totweight != 0.0f)
+		mdb->phi[acenter]= phi/totweight;
+}
+
+static void meshdeform_matrix_solve(MeshDeformBind *mdb)
+{
+	NLContext *context;
+	float vec[3], gridvec[3];
+	int a, b, x, y, z, totvar;
+	char message[1024];
+
+	/* setup variable indices */
+	mdb->varidx= MEM_callocN(sizeof(int)*mdb->size3, "MeshDeformDSvaridx");
+	for(a=0, totvar=0; a<mdb->size3; a++)
+		mdb->varidx[a]= (mdb->tag[a] == MESHDEFORM_TAG_EXTERIOR)? -1: totvar++;
+
+	if(totvar == 0) {
+		MEM_freeN(mdb->varidx);
+		return;
+	}
+
+	progress_bar(0, "Starting mesh deform solve");
+
+	/* setup opennl solver */
+	nlNewContext();
+	context= nlGetCurrent();
+
+	nlSolverParameteri(NL_NB_VARIABLES, totvar);
+	nlSolverParameteri(NL_NB_ROWS, totvar);
+	nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 1);
+
+	nlBegin(NL_SYSTEM);
+	nlBegin(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);
+
+	/* solve for each cage vert */
+	for(a=0; a<mdb->totcagevert; a++) {
+		if(a != 0) {
+			nlBegin(NL_SYSTEM);
+			nlBegin(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);
+
+		nlEnd(NL_MATRIX);
+		nlEnd(NL_SYSTEM);
+
+#if 0
+		nlPrintMatrix();
+#endif
+
+		if(nlSolveAdvanced(NULL, NL_TRUE)) {
+			for(z=0; z<mdb->size; z++)
+				for(y=0; y<mdb->size; y++)
+					for(x=0; x<mdb->size; x++)
+						meshdeform_matrix_add_semibound_phi(mdb, x, y, z, a);
+
+			for(z=0; z<mdb->size; z++)
+				for(y=0; y<mdb->size; y++)
+					for(x=0; x<mdb->size; x++)
+						meshdeform_matrix_add_exterior_phi(mdb, x, y, z, a);
+
+			for(b=0; b<mdb->size3; b++) {
+				if(mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR)
+					mdb->phi[b]= nlGetVariable(0, mdb->varidx[b]);
+				mdb->totalphi[b] += mdb->phi[b];
+			}
+
+			/* compute weights for each vertex */
+			for(b=0; b<mdb->totvert; b++) {
+				if(mdb->inside[b]) {
+					VECCOPY(vec, mdb->vertexcos[b]);
+					Mat4MulVecfl(mdb->cagemat, vec);
+					gridvec[0]= (vec[0] - mdb->min[0] - mdb->halfwidth[0])/mdb->width[0];
+					gridvec[1]= (vec[1] - mdb->min[1] - mdb->halfwidth[1])/mdb->width[1];
+					gridvec[2]= (vec[2] - mdb->min[2] - mdb->halfwidth[2])/mdb->width[2];
+
+					mdb->weights[b*mdb->totcagevert + a]= meshdeform_interp_w(mdb, gridvec, vec, a);
+				}
+			}
+		}
+		else {
+			error("Mesh Deform: failed to find solution.");
+			break;
+		}
+
+		sprintf(message, "Mesh deform solve %d / %d       |||", a+1, mdb->totcagevert);
+		progress_bar((float)(a+1)/(float)(mdb->totcagevert), message);
+	}
+
+#if 0
+	/* sanity check */
+	for(b=0; b<mdb->size3; b++)
+		if(mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR)
+			if(fabs(mdb->totalphi[b] - 1.0f) > 1e-4)
+				printf("totalphi deficiency [%s|%d] %d: %.10f\n",
+					(mdb->tag[b] == MESHDEFORM_TAG_INTERIOR)? "interior": "boundary", mdb->semibound[b], mdb->varidx[b], mdb->totalphi[b]);
+#endif
+	
+	/* free */
+	MEM_freeN(mdb->varidx);
+
+	nlDeleteContext(context);
+}
+
+void harmonic_coordinates_bind(MeshDeformModifierData *mmd, float (*vertexcos)[3], int totvert, float cagemat[][4])
+{
+	MeshDeformBind mdb;
+	MVert *mvert;
+	float center[3], vec[3], maxwidth;
+	int a, x, y, z, totinside;
+
+	waitcursor(1);
+	start_progress_bar();
+
+	/* free exisiting weights */
+	if(mmd->bindweights) {
+		MEM_freeN(mmd->bindweights);
+		MEM_freeN(mmd->bindcos);
+		mmd->bindweights= NULL;
+		mmd->bindcos= NULL;
+	}
+
+	memset(&mdb, 0, sizeof(MeshDeformBind));
+
+	/* get mesh and cage mesh */
+	mdb.vertexcos= vertexcos;
+	mdb.totvert= totvert;
+	
+	mdb.cagedm= mesh_create_derived_no_deform(mmd->object, NULL, CD_MASK_BAREMESH);
+	mdb.totcagevert= mdb.cagedm->getNumVerts(mdb.cagedm);
+	mdb.cagecos= MEM_callocN(sizeof(*mdb.cagecos)*mdb.totcagevert, "MeshDeformBindCos");
+	Mat4CpyMat4(mdb.cagemat, cagemat);
+
+	mvert= mdb.cagedm->getVertArray(mdb.cagedm);
+	for(a=0; a<mdb.totcagevert; a++)
+		VECCOPY(mdb.cagecos[a], mvert[a].co)
+
+	/* compute bounding box of the cage mesh */
+	INIT_MINMAX(mdb.min, mdb.max);
+
+	for(a=0; a<mdb.totcagevert; a++)
+		DO_MINMAX(mdb.cagecos[a], mdb.min, mdb.max);
+
+	/* allocate memory */
+	mdb.size= (2<<(mmd->gridsize-1)) + 2;
+	mdb.size3= mdb.size*mdb.size*mdb.size;
+	mdb.tag= MEM_callocN(sizeof(int)*mdb.size3, "MeshDeformBindTag");
+	mdb.phi= MEM_callocN(sizeof(float)*mdb.size3, "MeshDeformBindPhi");
+	mdb.totalphi= MEM_callocN(sizeof(float)*mdb.size3, "MeshDeformBindTotalPhi");
+	mdb.boundisect= MEM_callocN(sizeof(*mdb.boundisect)*mdb.size3, "MDefBoundIsect");
+	mdb.semibound= MEM_callocN(sizeof(int)*mdb.size3, "MDefSemiBound");
+
+	mdb.weights= MEM_callocN(sizeof(float)*mdb.totvert*mdb.totcagevert, "MDefWeights");
+	mdb.inside= MEM_callocN(sizeof(int)*mdb.totvert, "MDefInside");
+
+	mdb.memarena= BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE);
+	BLI_memarena_use_calloc(mdb.memarena);
+
+	/* make bounding box equal size in all directions, add padding, and compute
+	 * width of the cells */
+	maxwidth = -1.0f;
+	for(a=0; a<3; a++)
+		if(mdb.max[a]-mdb.min[a] > maxwidth)
+			maxwidth= mdb.max[a]-mdb.min[a];
+
+	for(a=0; a<3; a++) {
+		center[a]= (mdb.min[a]+mdb.max[a])*0.5f;
+		mdb.min[a]= center[a] - maxwidth*0.5f;
+		mdb.max[a]= center[a] + maxwidth*0.5f;
+
+		mdb.width[a]= (mdb.max[a]-mdb.min[a])/(mdb.size-4);
+		mdb.min[a] -= 2.1f*mdb.width[a];
+		mdb.max[a] += 2.1f*mdb.width[a];
+
+		mdb.width[a]= (mdb.max[a]-mdb.min[a])/mdb.size;
+		mdb.halfwidth[a]= mdb.width[a]*0.5f;
+	}
+
+	progress_bar(0, "Setting up mesh deform system");
+
+#if 0
+	/* create ray tree */
+	meshdeform_ray_tree_create(&mdb);
+#endif
+
+	totinside= 0;
+	for(a=0; a<mdb.totvert; a++) {
+		VECCOPY(vec, mdb.vertexcos[a]);
+		Mat4MulVecfl(mdb.cagemat, vec);
+		mdb.inside[a]= meshdeform_inside_cage(&mdb, vec);
+		if(mdb.inside[a])
+			totinside++;
+	}
+
+	/* free temporary MDefBoundIsects */
+	BLI_memarena_free(mdb.memarena);
+	mdb.memarena= BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE);
+
+	/* start with all cells untyped */
+	for(a=0; a<mdb.size3; a++)
+		mdb.tag[a]= MESHDEFORM_TAG_UNTYPED;
+	
+	/* detect intersections and tag boundary cells */
+	for(z=0; z<mdb.size; z++)
+		for(y=0; y<mdb.size; y++)
+			for(x=0; x<mdb.size; x++)
+				meshdeform_add_intersections(&mdb, x, y, z);
+
+#if 0
+	/* free ray tree */
+	meshdeform_ray_tree_free(&mdb);
+#endif
+
+	/* compute exterior and interior tags */
+	meshdeform_bind_floodfill(&mdb);
+
+	for(z=0; z<mdb.size; z++)
+		for(y=0; y<mdb.size; y++)
+			for(x=0; x<mdb.size; x++)
+				meshdeform_check_semibound(&mdb, x, y, z);
+
+	/* solve */
+	meshdeform_matrix_solve(&mdb);
+
+	/* assign results */
+	mmd->bindweights= mdb.weights;
+	mmd->bindcos= (float*)mdb.cagecos;
+	mmd->totvert= mdb.totvert;
+	mmd->totcagevert= mdb.totcagevert;
+
+	/* transform bindcos to world space */
+	for(a=0; a<mdb.totcagevert; a++)
+		Mat4MulVecfl(mmd->object->obmat, mmd->bindcos+a*3);
+
+	/* free */
+	mdb.cagedm->release(mdb.cagedm);
+	MEM_freeN(mdb.tag);
+	MEM_freeN(mdb.phi);
+	MEM_freeN(mdb.totalphi);
+	MEM_freeN(mdb.boundisect);
+	MEM_freeN(mdb.semibound);
+	MEM_freeN(mdb.inside);
+	BLI_memarena_free(mdb.memarena);
+
+	end_progress_bar();
+	waitcursor(0);
+}
+
diff --git a/source/blender/src/parametrizer.c b/source/blender/src/parametrizer.c
index 62a07a437b0..d0a51027ad3 100644
--- a/source/blender/src/parametrizer.c
+++ b/source/blender/src/parametrizer.c
@@ -2213,7 +2213,7 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 	nlBegin(NL_MATRIX);
 
 	for (i = 0; i < nvar; i++)
-		nlRightHandSideAdd(i, sys->bInterior[i]);
+		nlRightHandSideAdd(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];
@@ -2259,8 +2259,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(v1->u.id, j2[0][0]*beta[0]);
-			nlRightHandSideAdd(ninterior + v1->u.id, j2[1][0]*beta[1] + j2[2][0]*beta[2]);
+			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]);
 
 			row1[0] = j2[0][0]*W[0][0];
 			row2[0] = j2[0][0]*W[1][0];
@@ -2279,8 +2279,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 			sys->J2dt[e2->u.id][1] = j2[1][1] = 1.0*wi2;
 			sys->J2dt[e3->u.id][1] = j2[2][1] = p_abf_compute_sin_product(sys, v2, e3->u.id)*wi3;
 
-			nlRightHandSideAdd(v2->u.id, j2[1][1]*beta[1]);
-			nlRightHandSideAdd(ninterior + v2->u.id, j2[0][1]*beta[0] + j2[2][1]*beta[2]);
+			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]);
 
 			row1[1] = j2[1][1]*W[0][1];
 			row2[1] = j2[1][1]*W[1][1];
@@ -2299,8 +2299,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.0*wi3;
 
-			nlRightHandSideAdd(v3->u.id, j2[2][2]*beta[2]);
-			nlRightHandSideAdd(ninterior + v3->u.id, j2[0][2]*beta[0] + j2[1][2]*beta[1]);
+			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]);
 
 			row1[2] = j2[2][2]*W[0][2];
 			row2[2] = j2[2][2]*W[1][2];
@@ -2357,24 +2357,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(v1->u.id);
-				float x2 = nlGetVariable(ninterior + v1->u.id);
+				float x = nlGetVariable(0, v1->u.id);
+				float x2 = nlGetVariable(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(v2->u.id);
-				float x2 = nlGetVariable(ninterior + v2->u.id);
+				float x = nlGetVariable(0, v2->u.id);
+				float x2 = nlGetVariable(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(v3->u.id);
-				float x2 = nlGetVariable(ninterior + v3->u.id);
+				float x = nlGetVariable(0, v3->u.id);
+				float x2 = nlGetVariable(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;
@@ -2405,8 +2405,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 		}
 
 		for (i = 0; i < ninterior; i++) {
-			sys->lambdaPlanar[i] += nlGetVariable(i);
-			sys->lambdaLength[i] += nlGetVariable(ninterior + i);
+			sys->lambdaPlanar[i] += nlGetVariable(0, i);
+			sys->lambdaLength[i] += nlGetVariable(0, ninterior + i);
 		}
 	}
 
@@ -2738,8 +2738,8 @@ static void p_chart_lscm_load_solution(PChart *chart)
 	PVert *v;
 
 	for (v=chart->verts; v; v=v->nextlink) {
-		v->uv[0] = nlGetVariable(2*v->u.id);
-		v->uv[1] = nlGetVariable(2*v->u.id + 1);
+		v->uv[0] = nlGetVariable(0, 2*v->u.id);
+		v->uv[1] = nlGetVariable(0, 2*v->u.id + 1);
 	}
 }
 
@@ -2796,6 +2796,7 @@ static void p_chart_lscm_begin(PChart *chart, PBool live, PBool abf)
 
 		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();
@@ -2807,6 +2808,7 @@ static PBool p_chart_lscm_solve(PChart *chart)
 	PVert *v, *pin1 = chart->u.lscm.pin1, *pin2 = chart->u.lscm.pin2;
 	PFace *f;
 	float *alpha = chart->u.lscm.abf_alpha;
+	int row;
 
 	nlMakeCurrent(chart->u.lscm.context);
 
@@ -2826,10 +2828,10 @@ static PBool p_chart_lscm_solve(PChart *chart)
 		nlLockVariable(2*pin2->u.id);
 		nlLockVariable(2*pin2->u.id + 1);
 	
-		nlSetVariable(2*pin1->u.id, pin1->uv[0]);
-		nlSetVariable(2*pin1->u.id + 1, pin1->uv[1]);
-		nlSetVariable(2*pin2->u.id, pin2->uv[0]);
-		nlSetVariable(2*pin2->u.id + 1, pin2->uv[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]);
 	}
 	else {
 		/* set and lock the pins */
@@ -2838,8 +2840,8 @@ static PBool p_chart_lscm_solve(PChart *chart)
 				nlLockVariable(2*v->u.id);
 				nlLockVariable(2*v->u.id + 1);
 
-				nlSetVariable(2*v->u.id, v->uv[0]);
-				nlSetVariable(2*v->u.id + 1, v->uv[1]);
+				nlSetVariable(0, 2*v->u.id, v->uv[0]);
+				nlSetVariable(0, 2*v->u.id + 1, v->uv[1]);
 			}
 		}
 	}
@@ -2848,6 +2850,7 @@ static PBool p_chart_lscm_solve(PChart *chart)
 
 	nlBegin(NL_MATRIX);
 
+	row = 0;
 	for (f=chart->faces; f; f=f->nextlink) {
 		PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
 		PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
@@ -2890,6 +2893,7 @@ static PBool p_chart_lscm_solve(PChart *chart)
 		cosine = cos(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);
@@ -2905,6 +2909,21 @@ static PBool p_chart_lscm_solve(PChart *chart)
 		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.0);
+		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);
+		row++;
+
+		nlMatrixAdd(row, 2*v1->u.id,   sine);
+		nlMatrixAdd(row, 2*v1->u.id+1, cosine - 1.0);
+		nlMatrixAdd(row, 2*v2->u.id,   -sine);
+		nlMatrixAdd(row, 2*v2->u.id+1, -cosine);
+		nlMatrixAdd(row, 2*v3->u.id+1, 1.0);
+		row++;
+#endif
 	}
 
 	nlEnd(NL_MATRIX);