Mesh Deform Modifier

====================

The MeshDeform modifier can deform a mesh with another 'cage' mesh.
It is similar to a lattice modifier, but instead of being restricted
to the regular grid layout of a lattice, the cage mesh can be modeled
to fit the mesh better.

http://www.blender.org/development/current-projects/changes-since-244/modifiers/

Implementation Notes:
- OpenNL has been refactored a bit to allow least squares matrices to
  be built without passing the matrix row by row, but instead with
  random access. MDef doesn't need this actually, but it's using this
  version of OpenNL so I'm just committing it now.
- Mean value weights for polygons have been added to arithb.c, a type
  of barycentric coordinates for polygons with >= 3 vertices. This
  might be useful for other parts of blender too.

17 files changed, 1582 insertions, 277 deletions

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,32 +640,21 @@ 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;		
-	}
-	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;
+	case NL_NB_RIGHT_HAND_SIDES: {
+		__nlCurrentContext->nb_rhs = (NLuint)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;
+	default: {
+		__nl_assert_not_reached;
 	} break;
 	}
 }
@@ -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++;
 		}
 
-		context->n = n;
-
-		/* a least squares problem results in a symmetric matrix */
-		if(context->least_squares)
-			context->symmetric = NL_TRUE;
+		m = (context->nb_rows == 0)? n: context->nb_rows;
 
-		if(context->symmetric)
-			storage = (storage | __NL_SYMMETRIC);
-
-		/* SuperLU storage does not support symmetric storage */
-		storage = (storage & ~__NL_SYMMETRIC);
+		context->m = m;
+		context->n = n;
 
-		__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
-}
-
-static void __nlBeginRow() {
-	__nlTransition(__NL_STATE_MATRIX, __NL_STATE_ROW);
-	__nlRowColumnZero(&__nlCurrentContext->af);
-	__nlRowColumnZero(&__nlCurrentContext->al);
+	if(context->least_squares)
+		__nlSparseMatrix_transpose_mult_rows(&context->M, b, Mtb);
 }
 
-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
-					);
-				}
-			}
-		}
-
-		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;
+	__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;
+
+			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);