This are the equality constraints in the format required by NLOpt. It is basically multiplying out the matrices and taking the squared sum norm (Frobenius norm). I cannot figure out how to encode the sum of multiplications in the `hpijkr` format you provided, any help would be appreciated (I hope this code makes sense)..

constraints.cpp

//Identity constraint: `p . q = id`, where `p` and `q` are opposite projections on the same edge.
double identity_constraint(unsigned n, const double *x, double *, void *data)
{
	surface_data* d = (surface_data*) data;
	auto surface = d->surface;
	auto neighbors = d->neighbors;
	double error = 0;
	visit_edges(*neighbors, [&](hpuint p, hpuint i) {
			//`q` is CW!! in this case
			auto q = *(neighbors->begin() + 3*p + i);

			hpuint j;
			visit_triplet(*neighbors, q, [&](hpuint n0, hpuint n1, hpuint n2) { j = (p == n0) ? 0 : (p == n1) ? 1 : 2; });

			//3x3 matrix entries
			//Encoding: every patch identifies 3 3x3 matrices. Use neighbor index as "stride" of the 9 entries
			double ap = x[27*p + (9*i) + 0];
			double bp = x[27*p + (9*i) + 1];
			double cp = x[27*p + (9*i) + 2];
			double dp = x[27*p + (9*i) + 3];
			double ep = x[27*p + (9*i) + 4];
			double fp = x[27*p + (9*i) + 5];
			double gp = x[27*p + (9*i) + 6];
			double hp = x[27*p + (9*i) + 7];
			double ip = x[27*p + (9*i) + 8];

			double aq = x[27*q + (9*j) + 0];
			double bq = x[27*q + (9*j) + 1];
			double cq = x[27*q + (9*j) + 2];
			double dq = x[27*q + (9*j) + 3];
			double eq = x[27*q + (9*j) + 4];
			double fq = x[27*q + (9*j) + 5];
			double gq = x[27*q + (9*j) + 6];
			double hq = x[27*q + (9*j) + 7];
			double iq = x[27*q + (9*j) + 8];

			//Matrix multiplication and square diff to the identity
			error += norm(ap * aq + bp * dq + cp * gq - 1);
			error += norm(ap * bq + bp * eq + cp * hq);
			error += norm(ap * cq + bp * fq + cp * iq);
			error += norm(dp * aq + ep * dq + fp * gq);
			error += norm(dp * bq + ep * eq + fp * hq - 1);
			error += norm(dp * cq + ep * fq + fp * iq);
			error += norm(gp * aq + hp * dq + ip * gq);
			error += norm(gp * bq + hp * eq + ip * hq);
			error += norm(gp * cq + hp * fq + ip * iq - 1);
		});

	return error;
}


//Compatibility constraint: `p_1 . p_2 . p_3 = id`, where `p_i` are the three projections inside a given patch. Note that this can be rewritten as `p_1 . p_2 = p_3^(-1)` where the inverse denotes the projection running opposite the edge.
double compatibility_constraint(unsigned n, const double *x, double *, void *data)
{
	surface_data* d = (surface_data*) data;
	auto surface = d->surface;
	auto neighbors = d->neighbors;
	double error = 0;
	for(auto p = 0lu; p < surface->getNumberOfPatches(); p++) {
		//`q` is CW!! in this case
		auto q = *(neighbors->begin() + 3*p + 2); //Choose an arbitrary neighbor. TODO: check for null

		hpuint j;
		visit_triplet(*neighbors, q, [&](hpuint n0, hpuint n1, hpuint n2) { j = (p == n0) ? 0 : (p == n1) ? 1 : 2; });

		//3x3 matrix entries
		//Encoding: every patch identifies 3 3x3 matrices. Use neighbor index as "stride" of the 9 entries
		double a1 = x[27*p + (9*0) + 0];
		double b1 = x[27*p + (9*0) + 1];
		double c1 = x[27*p + (9*0) + 2];
		double d1 = x[27*p + (9*0) + 3];
		double e1 = x[27*p + (9*0) + 4];
		double f1 = x[27*p + (9*0) + 5];
		double g1 = x[27*p + (9*0) + 6];
		double h1 = x[27*p + (9*0) + 7];
		double i1 = x[27*p + (9*0) + 8];

		double a2 = x[27*p + (9*1) + 0];
		double b2 = x[27*p + (9*1) + 1];
		double c2 = x[27*p + (9*1) + 2];
		double d2 = x[27*p + (9*1) + 3];
		double e2 = x[27*p + (9*1) + 4];
		double f2 = x[27*p + (9*1) + 5];
		double g2 = x[27*p + (9*1) + 6];
		double h2 = x[27*p + (9*1) + 7];
		double i2 = x[27*p + (9*1) + 8];

		double aq = x[27*q + (9*j) + 0];
		double bq = x[27*q + (9*j) + 1];
		double cq = x[27*q + (9*j) + 2];
		double dq = x[27*q + (9*j) + 3];
		double eq = x[27*q + (9*j) + 4];
		double fq = x[27*q + (9*j) + 5];
		double gq = x[27*q + (9*j) + 6];
		double hq = x[27*q + (9*j) + 7];
		double iq = x[27*q + (9*j) + 8];

		//Matrix multiplication and square diff
		error += norm(a1 * a2 + b1 * d2 + c1 * g2 - aq);
		error += norm(a1 * b2 + b1 * e2 + c1 * h2 - bq);
		error += norm(a1 * c2 + b1 * f2 + c1 * i2 - cq);
		error += norm(d1 * a2 + e1 * d2 + f1 * g2 - dq);
		error += norm(d1 * b2 + e1 * e2 + f1 * h2 - eq);
		error += norm(d1 * c2 + e1 * f2 + f1 * i2 - fq);
		error += norm(g1 * a2 + h1 * d2 + i1 * g2 - gq);
		error += norm(g1 * b2 + h1 * e2 + i1 * h2 - hq);
		error += norm(g1 * c2 + h1 * f2 + i1 * i2 - iq);
	};

	return error;
}