diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C
index 35387607..24d0a23b 100644
--- a/src/ASM/LR/ASMu2D.C
+++ b/src/ASM/LR/ASMu2D.C
@@ -58,62 +58,62 @@ ASMu2D::ASMu2D (const ASMu2D& patch, unsigned char n_f)
 
 bool ASMu2D::read (std::istream& is)
 {
-	if (shareFE) return false;
+  if (shareFE) return false;
 
-	// read inputfile as either an LRSpline file directly or a tensor product B-spline and convert
-	char firstline[256];
-	is.getline(firstline, 256);
-	if(strncmp(firstline, "# LRSPLINE", 10) == 0) {
-		lrspline.reset(new LR::LRSplineSurface());
-		is >> *lrspline;
-	} else { // probably a SplineSurface, so we'll read that and convert
-		tensorspline = new Go::SplineSurface();
-		is >> *tensorspline;
-		lrspline.reset(new LR::LRSplineSurface(tensorspline));
-	}
+  // read inputfile as either an LRSpline file directly or a tensor product B-spline and convert
+  char firstline[256];
+  is.getline(firstline, 256);
+  if(strncmp(firstline, "# LRSPLINE", 10) == 0) {
+    lrspline.reset(new LR::LRSplineSurface());
+    is >> *lrspline;
+  } else { // probably a SplineSurface, so we'll read that and convert
+    tensorspline = new Go::SplineSurface();
+    is >> *tensorspline;
+    lrspline.reset(new LR::LRSplineSurface(tensorspline));
+  }
 
-	// Eat white-space characters to see if there is more data to read
-	char c;
-	while (is.get(c))
-		if (!isspace(c)) {
-			is.putback(c);
-			break;
-		}
+  // Eat white-space characters to see if there is more data to read
+  char c;
+  while (is.get(c))
+    if (!isspace(c)) {
+      is.putback(c);
+      break;
+    }
 
-	if (!is.good() && !is.eof())
-	{
-		std::cerr <<" *** ASMu2D::read: Failure reading spline data"<< std::endl;
-		lrspline.reset();
-		return false;
-	}
-	else if (lrspline->dimension() < 2)
-	{
-		std::cerr <<" *** ASMu2D::read: Invalid spline lrsplineace patch, dim="
-			  << lrspline->dimension() << std::endl;
-		lrspline.reset();
-		return false;
-	}
-	else if (lrspline->dimension() < nsd)
-	{
-		std::cout <<"  ** ASMu2D::read: The dimension of this lrsplineace patch "
-			  << lrspline->dimension() <<" is less than nsd="<< nsd
-			  <<".\n                   Resetting nsd to "<< lrspline->dimension()
-			  <<" for this patch."<< std::endl;
-		nsd = lrspline->dimension();
-	}
+  if (!is.good() && !is.eof())
+  {
+    std::cerr <<" *** ASMu2D::read: Failure reading spline data"<< std::endl;
+    lrspline.reset();
+    return false;
+  }
+  else if (lrspline->dimension() < 2)
+  {
+    std::cerr <<" *** ASMu2D::read: Invalid spline lrsplineace patch, dim="
+        << lrspline->dimension() << std::endl;
+    lrspline.reset();
+    return false;
+  }
+  else if (lrspline->dimension() < nsd)
+  {
+    std::cout <<"  ** ASMu2D::read: The dimension of this lrsplineace patch "
+        << lrspline->dimension() <<" is less than nsd="<< nsd
+        <<".\n                   Resetting nsd to "<< lrspline->dimension()
+        <<" for this patch."<< std::endl;
+    nsd = lrspline->dimension();
+  }
 
-	geo = lrspline.get();
-	return true;
+  geo = lrspline.get();
+  return true;
 }
 
 
 bool ASMu2D::write (std::ostream& os, int) const
 {
-	if (!lrspline) return false;
+  if (!lrspline) return false;
 
-	os << *lrspline;
+  os << *lrspline;
 
-	return os.good();
+  return os.good();
 }
 
 
@@ -143,152 +143,152 @@ bool ASMu2D::checkRightHandSystem ()
 
 bool ASMu2D::cornerRefine (int minBasisfunctions)
 {
-	if (!lrspline) return false;
-	if (shareFE) return true;
+  if (!lrspline) return false;
+  if (shareFE) return true;
 
-	double h = 1.0;
-	int nBasis = lrspline->nBasisFunctions();
-	double unif_step_h = 1.0 / ((minBasisfunctions - nBasis) / 3.0 + 1.0);
-	while(lrspline->nBasisFunctions() < minBasisfunctions) {
-		lrspline->insert_const_u_edge(h-unif_step_h, 0, h);
-		lrspline->insert_const_v_edge(h-unif_step_h, 0, h);
-		h -= unif_step_h;
-	}
+  double h = 1.0;
+  int nBasis = lrspline->nBasisFunctions();
+  double unif_step_h = 1.0 / ((minBasisfunctions - nBasis) / 3.0 + 1.0);
+  while(lrspline->nBasisFunctions() < minBasisfunctions) {
+    lrspline->insert_const_u_edge(h-unif_step_h, 0, h);
+    lrspline->insert_const_v_edge(h-unif_step_h, 0, h);
+    h -= unif_step_h;
+  }
 
-	std::ofstream paramMeshFile("mesh_param.eps");
-	std::ofstream physicalMeshFile("mesh_physical.eps");
-	lrspline->writePostscriptMesh(paramMeshFile);
-	lrspline->writePostscriptMesh(physicalMeshFile);
-	return true;
+  std::ofstream paramMeshFile("mesh_param.eps");
+  std::ofstream physicalMeshFile("mesh_physical.eps");
+  lrspline->writePostscriptMesh(paramMeshFile);
+  lrspline->writePostscriptMesh(physicalMeshFile);
+  return true;
 }
 
 bool ASMu2D::diagonalRefine (int minBasisfunctions)
 {
-	if (!lrspline) return false;
-	if (shareFE) return true;
+  if (!lrspline) return false;
+  if (shareFE) return true;
 
-	double end1 = lrspline->endparam(0);
-	double end2 = lrspline->endparam(1);
-	double h = 1.0;
-	int iter = 0;
-	double u = h/2.0;
-	double v = h/2.0;
-	while(lrspline->nBasisFunctions() < minBasisfunctions) {
-		lrspline->insert_const_u_edge(u, (iter-1<0) ? 0 : (iter-1)*h, ((iter+2)*h>end2) ? end2 : (iter+2)*h);
-		lrspline->insert_const_v_edge(v, (iter-1<0) ? 0 : (iter-1)*h, ((iter+2)*h>end1) ? end1 : (iter+2)*h);
-		u += h;
-		v += h;
-		iter++;
-		if( u>end1 ) {
-			h /= 2.0;
-			iter = 0;
-			u = h/2.0;
-			v = h/2.0;
-		}
-	}
+  double end1 = lrspline->endparam(0);
+  double end2 = lrspline->endparam(1);
+  double h = 1.0;
+  int iter = 0;
+  double u = h/2.0;
+  double v = h/2.0;
+  while(lrspline->nBasisFunctions() < minBasisfunctions) {
+    lrspline->insert_const_u_edge(u, (iter-1<0) ? 0 : (iter-1)*h, ((iter+2)*h>end2) ? end2 : (iter+2)*h);
+    lrspline->insert_const_v_edge(v, (iter-1<0) ? 0 : (iter-1)*h, ((iter+2)*h>end1) ? end1 : (iter+2)*h);
+    u += h;
+    v += h;
+    iter++;
+    if( u>end1 ) {
+      h /= 2.0;
+      iter = 0;
+      u = h/2.0;
+      v = h/2.0;
+    }
+  }
 
-	std::ofstream meshFile("mesh.eps");
-	lrspline->writePostscriptMesh(meshFile);
-	return true;
+  std::ofstream meshFile("mesh.eps");
+  lrspline->writePostscriptMesh(meshFile);
+  return true;
 }
 
 bool ASMu2D::uniformRefine (int minBasisfunctions)
 {
-	if (!lrspline) return false;
-	if (shareFE) return true;
+  if (!lrspline) return false;
+  if (shareFE) return true;
 
-	double end1 = lrspline->endparam(0);
-	double end2 = lrspline->endparam(1);
-	double h = 1.0;
-	bool step_u = true;
-	double u = h/2.0;
-	double v = h/2.0;
-	while(lrspline->nBasisFunctions() < minBasisfunctions) {
-		if(step_u) {
-			lrspline->insert_const_u_edge(u, 0, end2);
-			u += h;
-			if(u > end1) {
-				step_u = !step_u;
-				u = h/4.0;
-			}
-		} else {
-			lrspline->insert_const_v_edge(v, 0, end1);
-			v += h;
-			if(v > end2) {
-				step_u = !step_u;
-				v = h/4.0;
-				h /= 2.0;
-			}
-		}
-	}
+  double end1 = lrspline->endparam(0);
+  double end2 = lrspline->endparam(1);
+  double h = 1.0;
+  bool step_u = true;
+  double u = h/2.0;
+  double v = h/2.0;
+  while(lrspline->nBasisFunctions() < minBasisfunctions) {
+    if(step_u) {
+      lrspline->insert_const_u_edge(u, 0, end2);
+      u += h;
+      if(u > end1) {
+        step_u = !step_u;
+        u = h/4.0;
+      }
+    } else {
+      lrspline->insert_const_v_edge(v, 0, end1);
+      v += h;
+      if(v > end2) {
+        step_u = !step_u;
+        v = h/4.0;
+        h /= 2.0;
+      }
+    }
+  }
 
-	std::ofstream meshFile("mesh.eps");
-	lrspline->writePostscriptMesh(meshFile);
-	return true;
+  std::ofstream meshFile("mesh.eps");
+  lrspline->writePostscriptMesh(meshFile);
+  return true;
 }
 
 bool ASMu2D::uniformRefine (int dir, int nInsert)
 {
-	if (!tensorspline || dir < 0 || dir > 1 || nInsert < 1) return false;
-	if (shareFE) return true;
+  if (!tensorspline || dir < 0 || dir > 1 || nInsert < 1) return false;
+  if (shareFE) return true;
 
-	RealArray extraKnots;
-	RealArray::const_iterator uit = tensorspline->basis(dir).begin();
-	double ucurr, uprev = *(uit++);
-	while (uit != tensorspline->basis(dir).end())
-	{
-		ucurr = *(uit++);
-		if (ucurr > uprev)
-			for (int i = 0; i < nInsert; i++)
-			{
-				double xi = (double)(i+1)/(double)(nInsert+1);
-				extraKnots.push_back(ucurr*xi + uprev*(1.0-xi));
-			}
-		uprev = ucurr;
-	}
+  RealArray extraKnots;
+  RealArray::const_iterator uit = tensorspline->basis(dir).begin();
+  double ucurr, uprev = *(uit++);
+  while (uit != tensorspline->basis(dir).end())
+  {
+    ucurr = *(uit++);
+    if (ucurr > uprev)
+      for (int i = 0; i < nInsert; i++)
+      {
+        double xi = (double)(i+1)/(double)(nInsert+1);
+        extraKnots.push_back(ucurr*xi + uprev*(1.0-xi));
+      }
+    uprev = ucurr;
+  }
 
-	if (dir == 0)
-		tensorspline->insertKnot_u(extraKnots);
-	else
-		tensorspline->insertKnot_v(extraKnots);
+  if (dir == 0)
+    tensorspline->insertKnot_u(extraKnots);
+  else
+    tensorspline->insertKnot_v(extraKnots);
 
-	lrspline.reset(new LR::LRSplineSurface(tensorspline));
-	geo = lrspline.get();
+  lrspline.reset(new LR::LRSplineSurface(tensorspline));
+  geo = lrspline.get();
 
-	return true;
+  return true;
 }
 
 bool ASMu2D::refine (int dir, const RealArray& xi)
 {
-	if (!tensorspline || dir < 0 || dir > 1 || xi.empty()) return false;
-	if (xi.front() < 0.0 || xi.back() > 1.0) return false;
-	if (shareFE) return true;
+  if (!tensorspline || dir < 0 || dir > 1 || xi.empty()) return false;
+  if (xi.front() < 0.0 || xi.back() > 1.0) return false;
+  if (shareFE) return true;
 
-	RealArray extraKnots;
-	RealArray::const_iterator uit = tensorspline->basis(dir).begin();
-	double ucurr, uprev = *(uit++);
-	while (uit != tensorspline->basis(dir).end())
-	{
-		ucurr = *(uit++);
-		if (ucurr > uprev)
-			for (size_t i = 0; i < xi.size(); i++)
-				if (i > 0 && xi[i] < xi[i-1])
-					return false;
-				else
-					extraKnots.push_back(ucurr*xi[i] + uprev*(1.0-xi[i]));
+  RealArray extraKnots;
+  RealArray::const_iterator uit = tensorspline->basis(dir).begin();
+  double ucurr, uprev = *(uit++);
+  while (uit != tensorspline->basis(dir).end())
+  {
+    ucurr = *(uit++);
+    if (ucurr > uprev)
+      for (size_t i = 0; i < xi.size(); i++)
+        if (i > 0 && xi[i] < xi[i-1])
+          return false;
+        else
+          extraKnots.push_back(ucurr*xi[i] + uprev*(1.0-xi[i]));
 
-		uprev = ucurr;
-	}
+    uprev = ucurr;
+  }
 
-	if (dir == 0)
-		tensorspline->insertKnot_u(extraKnots);
-	else
-		tensorspline->insertKnot_v(extraKnots);
+  if (dir == 0)
+    tensorspline->insertKnot_u(extraKnots);
+  else
+    tensorspline->insertKnot_v(extraKnots);
 
-	lrspline.reset(new LR::LRSplineSurface(tensorspline));
-	geo = lrspline.get();
+  lrspline.reset(new LR::LRSplineSurface(tensorspline));
+  geo = lrspline.get();
 
-	return true;
+  return true;
 }
 
 
@@ -480,117 +480,117 @@ bool ASMu2D::generateFEMTopology ()
 
 bool ASMu2D::connectPatch (int edge, ASMu2D& neighbor, int nedge, bool revers)
 {
-	return this->connectBasis(edge,neighbor,nedge,revers);
+  return this->connectBasis(edge,neighbor,nedge,revers);
 }
 
 
 bool ASMu2D::connectBasis (int edge, ASMu2D& neighbor, int nedge, bool revers,
                            int basis, int slave, int master)
 {
-	// Set up the slave node numbers for this surface patch
+  // Set up the slave node numbers for this surface patch
 
-	int n1, n2;
-	if (!this->getSize(n1,n2,basis)) return false;
-	int node = slave+1, i1 = 0;
+  int n1, n2;
+  if (!this->getSize(n1,n2,basis)) return false;
+  int node = slave+1, i1 = 0;
 
-	switch (edge)
-		{
-		case 2: // Positive I-direction
-			node += n1-1;
-		case 1: // Negative I-direction
-			i1 = n1;
-			n1 = n2;
-			break;
+  switch (edge)
+    {
+    case 2: // Positive I-direction
+      node += n1-1;
+    case 1: // Negative I-direction
+      i1 = n1;
+      n1 = n2;
+      break;
 
-		case 4: // Positive J-direction
-			node += n1*(n2-1);
-		case 3: // Negative J-direction
-			i1 = 1;
-			break;
+    case 4: // Positive J-direction
+      node += n1*(n2-1);
+    case 3: // Negative J-direction
+      i1 = 1;
+      break;
 
-		default:
-			std::cerr <<" *** ASMu2D::connectPatch: Invalid slave edge "
-			          << edge << std::endl;
-			return false;
-		}
+    default:
+      std::cerr <<" *** ASMu2D::connectPatch: Invalid slave edge "
+                << edge << std::endl;
+      return false;
+    }
 
-	int i;
-	IntVec slaveNodes(n1,0);
-	for (i = 0; i < n1; i++, node += i1)
-		slaveNodes[i] = node;
+  int i;
+  IntVec slaveNodes(n1,0);
+  for (i = 0; i < n1; i++, node += i1)
+    slaveNodes[i] = node;
 
-	// Set up the master node numbers for the neighboring surface patch
+  // Set up the master node numbers for the neighboring surface patch
 
-	if (!neighbor.getSize(n1,n2,basis)) return false;
-	node = master+1; i1 = 0;
+  if (!neighbor.getSize(n1,n2,basis)) return false;
+  node = master+1; i1 = 0;
 
-	switch (nedge)
-		{
-		case 2: // Positive I-direction
-			node += n1-1;
-		case 1: // Negative I-direction
-			i1 = n1;
-			n1 = n2;
-			break;
+  switch (nedge)
+    {
+    case 2: // Positive I-direction
+      node += n1-1;
+    case 1: // Negative I-direction
+      i1 = n1;
+      n1 = n2;
+      break;
 
-		case 4: // Positive J-direction
-			node += n1*(n2-1);
-		case 3: // Negative J-direction
-			i1 = 1;
-			break;
+    case 4: // Positive J-direction
+      node += n1*(n2-1);
+    case 3: // Negative J-direction
+      i1 = 1;
+      break;
 
-		default:
-			std::cerr <<" *** ASMu2D::connectPatch: Invalid master edge "
-			          << nedge << std::endl;
-			return false;
-		}
+    default:
+      std::cerr <<" *** ASMu2D::connectPatch: Invalid master edge "
+                << nedge << std::endl;
+      return false;
+    }
 
-	if (n1 != (int)slaveNodes.size())
-	{
-		std::cerr <<" *** ASMu2D::connectPatch: Non-matching edges, sizes "
-		          << n1 <<" and "<< slaveNodes.size() << std::endl;
-		return false;
-	}
+  if (n1 != (int)slaveNodes.size())
+  {
+    std::cerr <<" *** ASMu2D::connectPatch: Non-matching edges, sizes "
+              << n1 <<" and "<< slaveNodes.size() << std::endl;
+    return false;
+  }
 
-	const double xtol = 1.0e-4;
-	for (i = 0; i < n1; i++, node += i1)
-	{
-		int k = revers ? n1-i-1 : i;
-		if (!neighbor.getCoord(node).equal(this->getCoord(slaveNodes[k]),xtol))
-		{
-			std::cerr <<" *** ASMu2D::connectPatch: Non-matching nodes "
-			          << node <<": "<< neighbor.getCoord(node)
-			          <<"\n                                          and "
-			          << slaveNodes[k] <<": "<< this->getCoord(slaveNodes[k])
-			          << std::endl;
-			return false;
-		}
-		else
-			ASMbase::collapseNodes(neighbor.MLGN[node-1],MLGN[slaveNodes[k]-1]);
-	}
+  const double xtol = 1.0e-4;
+  for (i = 0; i < n1; i++, node += i1)
+  {
+    int k = revers ? n1-i-1 : i;
+    if (!neighbor.getCoord(node).equal(this->getCoord(slaveNodes[k]),xtol))
+    {
+      std::cerr <<" *** ASMu2D::connectPatch: Non-matching nodes "
+                << node <<": "<< neighbor.getCoord(node)
+                <<"\n                                          and "
+                << slaveNodes[k] <<": "<< this->getCoord(slaveNodes[k])
+                << std::endl;
+      return false;
+    }
+    else
+      ASMbase::collapseNodes(neighbor.MLGN[node-1],MLGN[slaveNodes[k]-1]);
+  }
 
-	return true;
+  return true;
 }
 
 
 void ASMu2D::closeEdges (int dir, int basis, int master)
 {
-	int n1, n2;
-	if (basis < 1) basis = 1;
-	if (!this->getSize(n1,n2,basis)) return;
+  int n1, n2;
+  if (basis < 1) basis = 1;
+  if (!this->getSize(n1,n2,basis)) return;
 
-	switch (dir)
-		{
-		case 1: // Edges are closed in I-direction
-			for (int i2 = 1; i2 <= n2; i2++, master += n1)
-				this->makePeriodic(master,master+n1-1);
-			break;
+  switch (dir)
+    {
+    case 1: // Edges are closed in I-direction
+      for (int i2 = 1; i2 <= n2; i2++, master += n1)
+        this->makePeriodic(master,master+n1-1);
+      break;
 
-		case 2: // Edges are closed in J-direction
-			for (int i1 = 1; i1 <= n1; i1++, master++)
-				this->makePeriodic(master,master+n1*(n2-1));
-			break;
-		}
+    case 2: // Edges are closed in J-direction
+      for (int i1 = 1; i1 <= n1; i1++, master++)
+        this->makePeriodic(master,master+n1*(n2-1));
+      break;
+    }
 }
 */
 
@@ -682,7 +682,7 @@ void ASMu2D::constrainEdge (int dir, bool open, int dof, int code, char basis)
 
 
 size_t ASMu2D::constrainEdgeLocal (int dir, bool open, int dof, int code,
-				   bool project)
+                                   bool project)
 {
   return 0; // TODO...
 }
@@ -716,7 +716,7 @@ int ASMu2D::getCorner(int I, int J, int basis) const
 
   if (edgeFunctions.size() > 1)
     std::cerr <<"  ** ASMu2D::constrainCorner: "<< edgeFunctions.size()
-	      <<" corners returned from LRSplineSurface::getEdgeFunctions()"
+              <<" corners returned from LRSplineSurface::getEdgeFunctions()"
               << std::endl;
 
   size_t ofs = 1;
@@ -744,17 +744,17 @@ void ASMu2D::constrainNode (double xi, double eta, int dof, int code)
 {
   std::cerr <<" *** ASMu2D::constrainNode: Not implemented yet!"<< std::endl;
   /*
-	if (xi  < 0.0 || xi  > 1.0) return;
-	if (eta < 0.0 || eta > 1.0) return;
+  if (xi  < 0.0 || xi  > 1.0) return;
+  if (eta < 0.0 || eta > 1.0) return;
 
-	int n1, n2;
-	if (!this->getSize(n1,n2,1)) return;
+  int n1, n2;
+  if (!this->getSize(n1,n2,1)) return;
 
-	int node = 1;
-	if (xi  > 0.0) node += int(0.5+(n1-1)*xi);
-	if (eta > 0.0) node += n1*int(0.5+(n2-1)*eta);
+  int node = 1;
+  if (xi  > 0.0) node += int(0.5+(n1-1)*xi);
+  if (eta > 0.0) node += n1*int(0.5+(n2-1)*eta);
 
-	this->prescribe(node,dof,code);
+  this->prescribe(node,dof,code);
   */
 }
 
@@ -874,7 +874,7 @@ size_t ASMu2D::getNoBoundaryElms (char lIndex, char ldim) const
   std::vector<LR::Element*> edgeElms;
   lrspline->getEdgeElements(edgeElms, edge);
   
-	return edgeElms.size();
+  return edgeElms.size();
 }
 
 
@@ -964,7 +964,7 @@ bool ASMu2D::integrate (Integrand& integrand,
     {
       this->getGaussPointParameters(gpar[d],d,nGauss,iel,xg);
       if (xr)
-	this->getGaussPointParameters(redpar[d],d,nRed,iel,xr);
+        this->getGaussPointParameters(redpar[d],d,nRed,iel,xr);
     }
 
     if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
@@ -990,33 +990,33 @@ bool ASMu2D::integrate (Integrand& integrand,
       // --- Selective reduced integration loop --------------------------------
 
       for (int j = 0; j < nRed; j++)
-	for (int i = 0; i < nRed; i++)
-	{
-	  // Local element coordinates of current integration point
-	  fe.xi  = xr[i];
-	  fe.eta = xr[j];
+        for (int i = 0; i < nRed; i++)
+        {
+          // Local element coordinates of current integration point
+          fe.xi  = xr[i];
+          fe.eta = xr[j];
 
-	  // Parameter values of current integration point
-	  fe.u = redpar[0][i];
-	  fe.v = redpar[1][j];
+          // Parameter values of current integration point
+          fe.u = redpar[0][i];
+          fe.v = redpar[1][j];
 
-	  // Compute basis function derivatives at current point
-	  Go::BasisDerivsSf spline;
-	  lrspline->computeBasis(fe.u,fe.v,spline,iel-1);
-	  SplineUtils::extractBasis(spline,fe.N,dNdu);
+          // Compute basis function derivatives at current point
+          Go::BasisDerivsSf spline;
+          lrspline->computeBasis(fe.u,fe.v,spline,iel-1);
+          SplineUtils::extractBasis(spline,fe.N,dNdu);
 
-	  // Compute Jacobian inverse and derivatives
-	  fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
+          // Compute Jacobian inverse and derivatives
+          fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
 
-	  // Cartesian coordinates of current integration point
-	  X = Xnod * fe.N;
-	  X.t = time.t;
+          // Cartesian coordinates of current integration point
+          X = Xnod * fe.N;
+          X.t = time.t;
 
-	  // Compute the reduced integration terms of the integrand
-	  fe.detJxW *= 0.25*dA*wr[i]*wr[j];
-	  if (!integrand.reducedInt(*A,fe,X))
-	    return false;
-	}
+          // Compute the reduced integration terms of the integrand
+          fe.detJxW *= 0.25*dA*wr[i]*wr[j];
+          if (!integrand.reducedInt(*A,fe,X))
+            return false;
+        }
     }
 
     // --- Integration loop over all Gauss points in each direction ------------
@@ -1027,57 +1027,57 @@ bool ASMu2D::integrate (Integrand& integrand,
     for (int j = 0; j < nGauss; j++)
       for (int i = 0; i < nGauss; i++, fe.iGP++)
       {
-	// Local element coordinates of current integration point
-	fe.xi  = xg[i];
-	fe.eta = xg[j];
+        // Local element coordinates of current integration point
+        fe.xi  = xg[i];
+        fe.eta = xg[j];
 
-	// Parameter values of current integration point
-	fe.u = gpar[0][i];
-	fe.v = gpar[1][j];
+        // Parameter values of current integration point
+        fe.u = gpar[0][i];
+        fe.v = gpar[1][j];
 
-	// Compute basis function derivatives at current integration point
-	if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES) {
-	  Go::BasisDerivsSf2 spline;
-	  lrspline->computeBasis(fe.u,fe.v,spline,iel-1);
-	  SplineUtils::extractBasis(spline,fe.N,dNdu,d2Ndu2);
-	}
-	else {
-	  Go::BasisDerivsSf spline;
-	  lrspline->computeBasis(fe.u,fe.v,spline, iel-1);
-	  SplineUtils::extractBasis(spline,fe.N,dNdu);
+        // Compute basis function derivatives at current integration point
+        if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES) {
+          Go::BasisDerivsSf2 spline;
+          lrspline->computeBasis(fe.u,fe.v,spline,iel-1);
+          SplineUtils::extractBasis(spline,fe.N,dNdu,d2Ndu2);
+        }
+        else {
+          Go::BasisDerivsSf spline;
+          lrspline->computeBasis(fe.u,fe.v,spline, iel-1);
+          SplineUtils::extractBasis(spline,fe.N,dNdu);
 #if SP_DEBUG > 4
-	  std::cout <<"\nBasis functions at a integration point "
-		    <<" : (u,v) = "<< spline.param[0] <<" "<< spline.param[1]
-		    <<"  left_idx = "<< spline.left_idx[0] <<" "<< spline.left_idx[1];
-	  for (size_t ii = 0; ii < spline.basisValues.size(); ii++)
-	    std::cout <<'\n'<< 1+ii <<'\t' << spline.basisValues[ii] <<'\t'
-		      << spline.basisDerivs_u[ii] <<'\t'<< spline.basisDerivs_v[ii];
-	  std::cout << std::endl;
+          std::cout <<"\nBasis functions at a integration point "
+                    <<" : (u,v) = "<< spline.param[0] <<" "<< spline.param[1]
+                    <<"  left_idx = "<< spline.left_idx[0] <<" "<< spline.left_idx[1];
+          for (size_t ii = 0; ii < spline.basisValues.size(); ii++)
+            std::cout <<'\n'<< 1+ii <<'\t' << spline.basisValues[ii] <<'\t'
+                      << spline.basisDerivs_u[ii] <<'\t'<< spline.basisDerivs_v[ii];
+          std::cout << std::endl;
 #endif
-	}
+        }
 
-	// Compute Jacobian inverse of coordinate mapping and derivatives
-	fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
-	if (fe.detJxW == 0.0) continue; // skip singular points
+        // Compute Jacobian inverse of coordinate mapping and derivatives
+        fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
+        if (fe.detJxW == 0.0) continue; // skip singular points
 
-	// Compute Hessian of coordinate mapping and 2nd order derivatives
-	if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
-	  if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
-	    return false;
+        // Compute Hessian of coordinate mapping and 2nd order derivatives
+        if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
+          if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
+            return false;
 
 #if SP_DEBUG > 4
-	std::cout <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
+        std::cout <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
 #endif
 
-	// Cartesian coordinates of current integration point
-	X = Xnod * fe.N;
-	X.t = time.t;
+        // Cartesian coordinates of current integration point
+        X = Xnod * fe.N;
+        X.t = time.t;
 
-	// Evaluate the integrand and accumulate element contributions
-	fe.detJxW *= 0.25*dA*wg[i]*wg[j];
-	PROFILE3("Integrand::evalInt");
-	if (!integrand.evalInt(*A,fe,time,X))
-	  return false;
+        // Evaluate the integrand and accumulate element contributions
+        fe.detJxW *= 0.25*dA*wg[i]*wg[j];
+        PROFILE3("Integrand::evalInt");
+        if (!integrand.evalInt(*A,fe,time,X))
+          return false;
       }
 
     // Finalize the element quantities
@@ -1163,16 +1163,16 @@ bool ASMu2D::integrate (Integrand& integrand,
       fe.u = elmPts[ip][0];
       fe.v = elmPts[ip][1];
 
-	// Compute basis function derivatives at current integration point
+        // Compute basis function derivatives at current integration point
       if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES) {
-	Go::BasisDerivsSf2 spline;
-	lrspline->computeBasis(fe.u,fe.v,spline,iel-1);
-	SplineUtils::extractBasis(spline,fe.N,dNdu,d2Ndu2);
+        Go::BasisDerivsSf2 spline;
+        lrspline->computeBasis(fe.u,fe.v,spline,iel-1);
+        SplineUtils::extractBasis(spline,fe.N,dNdu,d2Ndu2);
       }
       else {
-	Go::BasisDerivsSf spline;
-	lrspline->computeBasis(fe.u,fe.v,spline,iel-1);
-	SplineUtils::extractBasis(spline,fe.N,dNdu);
+        Go::BasisDerivsSf spline;
+        lrspline->computeBasis(fe.u,fe.v,spline,iel-1);
+        SplineUtils::extractBasis(spline,fe.N,dNdu);
       }
 
       // Compute Jacobian inverse of coordinate mapping and derivatives
@@ -1181,12 +1181,12 @@ bool ASMu2D::integrate (Integrand& integrand,
 
       // Compute Hessian of coordinate mapping and 2nd order derivatives
       if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
-	if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
-	  return false;
+        if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
+          return false;
 
 #if SP_DEBUG > 4
       std::cout <<"\niel, ip = "<< iel <<" "<< ip
-		<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
+                <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
 #endif
 
       // Cartesian coordinates of current integration point
@@ -1197,7 +1197,7 @@ bool ASMu2D::integrate (Integrand& integrand,
       fe.detJxW *= 0.25*dA*elmPts[ip][2];
       PROFILE3("Integrand::evalInt");
       if (!integrand.evalInt(*A,fe,time,X))
-	return false;
+        return false;
     }
 
     // Finalize the element quantities
@@ -1326,7 +1326,7 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
       // Evaluate the integrand and accumulate element contributions
       fe.detJxW *= 0.5*dS*wg[i];
       if (!integrand.evalBou(*A,fe,time,X,normal))
-	return false;
+        return false;
     }
 
     // Finalize the element quantities
@@ -1346,102 +1346,102 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
 
 int ASMu2D::evalPoint (const double* xi, double* param, Vec3& X) const
 {
-	if (!lrspline) return -2;
+  if (!lrspline) return -2;
 
-	param[0] = (1.0-xi[0])*lrspline->startparam(0) + xi[0]*lrspline->endparam(0);
-	param[1] = (1.0-xi[1])*lrspline->startparam(1) + xi[1]*lrspline->endparam(1);
+  param[0] = (1.0-xi[0])*lrspline->startparam(0) + xi[0]*lrspline->endparam(0);
+  param[1] = (1.0-xi[1])*lrspline->startparam(1) + xi[1]*lrspline->endparam(1);
 
-	FiniteElement fe;
-	fe.iel = 1 + lrspline->getElementContaining(param[0],param[1]);
-	fe.u   = param[0];
-	fe.v   = param[1];
+  FiniteElement fe;
+  fe.iel = 1 + lrspline->getElementContaining(param[0],param[1]);
+  fe.u   = param[0];
+  fe.v   = param[1];
 
-	Matrix Xnod;
-	if (!this->getElementCoordinates(Xnod,fe.iel))
-	  return -1;
+  Matrix Xnod;
+  if (!this->getElementCoordinates(Xnod,fe.iel))
+    return -1;
 
-	if (!this->evaluateBasis(fe))
-	  return -1;
+  if (!this->evaluateBasis(fe))
+    return -1;
 
-	X = Xnod * fe.N;
+  X = Xnod * fe.N;
 
-	return 0;
+  return 0;
 }
 
 
 bool ASMu2D::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) const
 {
 #ifdef SP_DEBUG
-	std::cout << "ASMu2D::getGridParameters(  )\n";
+  std::cout << "ASMu2D::getGridParameters(  )\n";
 #endif
 
-	// output is written once for each element resulting in a lot of unnecessary storage
-	// this is preferable to figuring out all element topology information
+  // output is written once for each element resulting in a lot of unnecessary storage
+  // this is preferable to figuring out all element topology information
 
-	std::vector<LR::Element*>::iterator el;
-	for(el=lrspline->elementBegin(); el<lrspline->elementEnd(); el++) {
-		// evaluate element at element corner points
-		double umin = (**el).umin();
-		double umax = (**el).umax();
-		double vmin = (**el).vmin();
-		double vmax = (**el).vmax();
-		for(int iv=0; iv<=nSegPerSpan; iv++) {
-			for(int iu=0; iu<=nSegPerSpan; iu++) {
-				double u = umin + (umax-umin)/nSegPerSpan*iu;
-				double v = vmin + (vmax-vmin)/nSegPerSpan*iv;
-				if(dir==0)
-					prm.push_back(u);
-				else
-					prm.push_back(v);
-			}
-		}
-	}
-	return true;
+  std::vector<LR::Element*>::iterator el;
+  for(el=lrspline->elementBegin(); el<lrspline->elementEnd(); el++) {
+    // evaluate element at element corner points
+    double umin = (**el).umin();
+    double umax = (**el).umax();
+    double vmin = (**el).vmin();
+    double vmax = (**el).vmax();
+    for(int iv=0; iv<=nSegPerSpan; iv++) {
+      for(int iu=0; iu<=nSegPerSpan; iu++) {
+        double u = umin + (umax-umin)/nSegPerSpan*iu;
+        double v = vmin + (vmax-vmin)/nSegPerSpan*iv;
+        if(dir==0)
+          prm.push_back(u);
+        else
+          prm.push_back(v);
+      }
+    }
+  }
+  return true;
 }
 
 #if 0
-	if (!lrspline) return false;
+  if (!lrspline) return false;
 
-	if (nSegPerSpan < 1)
-	{
-		std::cerr <<" *** ASMu2D::getGridParameters: Too few knot-span points "
-		          << nSegPerSpan+1 <<" in direction "<< dir << std::endl;
-		return false;
-	}
+  if (nSegPerSpan < 1)
+  {
+    std::cerr <<" *** ASMu2D::getGridParameters: Too few knot-span points "
+              << nSegPerSpan+1 <<" in direction "<< dir << std::endl;
+    return false;
+  }
 
-	RealArray::const_iterator uit = lrspline->basis(dir).begin();
-	double ucurr = 0.0, uprev = *(uit++);
-	while (uit != lrspline->basis(dir).end())
-	{
-		ucurr = *(uit++);
-		if (ucurr > uprev)
-			if (nSegPerSpan == 1)
-				prm.push_back(uprev);
-			else for (int i = 0; i < nSegPerSpan; i++)
-			{
-				double xg = (double)(2*i-nSegPerSpan)/(double)nSegPerSpan;
-				prm.push_back(0.5*(ucurr*(1.0+xg) + uprev*(1.0-xg)));
-			}
-		uprev = ucurr;
-	}
+  RealArray::const_iterator uit = lrspline->basis(dir).begin();
+  double ucurr = 0.0, uprev = *(uit++);
+  while (uit != lrspline->basis(dir).end())
+  {
+    ucurr = *(uit++);
+    if (ucurr > uprev)
+      if (nSegPerSpan == 1)
+        prm.push_back(uprev);
+      else for (int i = 0; i < nSegPerSpan; i++)
+      {
+        double xg = (double)(2*i-nSegPerSpan)/(double)nSegPerSpan;
+        prm.push_back(0.5*(ucurr*(1.0+xg) + uprev*(1.0-xg)));
+      }
+    uprev = ucurr;
+  }
 
-	if (ucurr > prm.back())
-		prm.push_back(ucurr);
-	return true;
+  if (ucurr > prm.back())
+    prm.push_back(ucurr);
+  return true;
 }
 
 
 bool ASMu2D::getGrevilleParameters (RealArray& prm, int dir) const
 {
-	if (!lrspline) return false;
+  if (!lrspline) return false;
 
-	const Go::BsplineBasis& basis = lrspline->basis(dir);
+  const Go::BsplineBasis& basis = lrspline->basis(dir);
 
-	prm.resize(basis.numCoefs());
-	for (size_t i = 0; i < prm.size(); i++)
-		prm[i] = basis.grevilleParameter(i);
+  prm.resize(basis.numCoefs());
+  for (size_t i = 0; i < prm.size(); i++)
+    prm[i] = basis.grevilleParameter(i);
 
-	return true;
+  return true;
 }
 #endif
 
@@ -1449,64 +1449,64 @@ bool ASMu2D::getGrevilleParameters (RealArray& prm, int dir) const
 bool ASMu2D::tesselate (ElementBlock& grid, const int* npe) const
 {
 #ifdef SP_DEBUG
-	std::cout << "ASMu2D::tesselate(  )\n";
+  std::cout << "ASMu2D::tesselate(  )\n";
 #endif
-	if(!lrspline) return false;
+  if(!lrspline) return false;
 
-	if(npe[0] != npe[1]) {
-		std::cerr << "ASMu2D::tesselate does not support different tesselation resolution in "
-		          << "u- and v-direction. nviz u = " << npe[0] << ", nviz v = " << npe[1] << "\n";
-		return false;
-	}
+  if(npe[0] != npe[1]) {
+    std::cerr << "ASMu2D::tesselate does not support different tesselation resolution in "
+              << "u- and v-direction. nviz u = " << npe[0] << ", nviz v = " << npe[1] << "\n";
+    return false;
+  }
 
-	int nNodesPerElement =  npe[0]   * npe[1];
-	int nSubElPerElement = (npe[0]-1)*(npe[1]-1);
-	int nElements        = lrspline->nElements();
+  int nNodesPerElement =  npe[0]   * npe[1];
+  int nSubElPerElement = (npe[0]-1)*(npe[1]-1);
+  int nElements        = lrspline->nElements();
 
-	// output is written once for each element resulting in a lot of unnecessary storage
-	// this is preferable to figuring out all element topology information
-	grid.unStructResize(nElements * nSubElPerElement,
-	                    nElements * nNodesPerElement);
+  // output is written once for each element resulting in a lot of unnecessary storage
+  // this is preferable to figuring out all element topology information
+  grid.unStructResize(nElements * nSubElPerElement,
+                      nElements * nNodesPerElement);
 
-	std::vector<LR::Element*>::iterator el;
-	int inod = 0;
-	int iel = 0;
-	for(el=lrspline->elementBegin(); el<lrspline->elementEnd(); el++, iel++) {
-		// evaluate element at element corner points
-		double umin = (**el).umin();
-		double umax = (**el).umax();
-		double vmin = (**el).vmin();
-		double vmax = (**el).vmax();
-		for(int iv=0; iv<npe[1]; iv++) {
-			for(int iu=0; iu<npe[0]; iu++) {
-				double u = umin + (umax-umin)/(npe[0]-1)*iu;
-				double v = vmin + (vmax-vmin)/(npe[1]-1)*iv;
-				Go::Point pt;
-				lrspline->point(pt, u,v, iel, iu!=npe[0]-1, iv!=npe[1]-1);
-				for(int dim=0; dim<nsd; dim++)
-					grid.setCoor(inod, dim, pt[dim]);
-				inod++;
-			}
-		}
-	}
+  std::vector<LR::Element*>::iterator el;
+  int inod = 0;
+  int iel = 0;
+  for(el=lrspline->elementBegin(); el<lrspline->elementEnd(); el++, iel++) {
+    // evaluate element at element corner points
+    double umin = (**el).umin();
+    double umax = (**el).umax();
+    double vmin = (**el).vmin();
+    double vmax = (**el).vmax();
+    for(int iv=0; iv<npe[1]; iv++) {
+      for(int iu=0; iu<npe[0]; iu++) {
+        double u = umin + (umax-umin)/(npe[0]-1)*iu;
+        double v = vmin + (vmax-vmin)/(npe[1]-1)*iv;
+        Go::Point pt;
+        lrspline->point(pt, u,v, iel, iu!=npe[0]-1, iv!=npe[1]-1);
+        for(int dim=0; dim<nsd; dim++)
+          grid.setCoor(inod, dim, pt[dim]);
+        inod++;
+      }
+    }
+  }
 
-	int ip = 0;
-	iel = 0;
-	for(int i=0; i<lrspline->nElements(); i++) {
-		int iStart = i*nNodesPerElement;
-		for(int iv=0; iv<npe[1]-1; iv++) {
-			for(int iu=0; iu<npe[0]-1; iu++, iel++) {
-				// enumerate nodes counterclockwise around the quad
-				grid.setNode(ip++, iStart + (iv  )*npe[0] + (iu  ) );
-				grid.setNode(ip++, iStart + (iv  )*npe[0] + (iu+1) );
-				grid.setNode(ip++, iStart + (iv+1)*npe[0] + (iu+1) );
-				grid.setNode(ip++, iStart + (iv+1)*npe[0] + (iu  ) );
-				grid.setElmId(iel+1, i+1);
-			}
-		}
-	}
+  int ip = 0;
+  iel = 0;
+  for(int i=0; i<lrspline->nElements(); i++) {
+    int iStart = i*nNodesPerElement;
+    for(int iv=0; iv<npe[1]-1; iv++) {
+      for(int iu=0; iu<npe[0]-1; iu++, iel++) {
+        // enumerate nodes counterclockwise around the quad
+        grid.setNode(ip++, iStart + (iv  )*npe[0] + (iu  ) );
+        grid.setNode(ip++, iStart + (iv  )*npe[0] + (iu+1) );
+        grid.setNode(ip++, iStart + (iv+1)*npe[0] + (iu+1) );
+        grid.setNode(ip++, iStart + (iv+1)*npe[0] + (iu  ) );
+        grid.setElmId(iel+1, i+1);
+      }
+    }
+  }
 
-	return true;
+  return true;
 }
 
 
@@ -1532,7 +1532,7 @@ bool ASMu2D::evalSolution (Matrix& sField, const Vector& locSol,
 {
 #ifdef SP_DEBUG
   std::cout <<"ASMu2D::evalSolution(Matrix&,const Vector&,const RealArray*,"
-	    <<"bool,int)"<< std::endl;
+            <<"bool,int)"<< std::endl;
 #endif
   size_t nComp = locSol.size() / this->getNoNodes();
   if (nComp*this->getNoNodes() != locSol.size())
@@ -1625,26 +1625,26 @@ bool ASMu2D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
     // Compute parameter values of the result sampling points
     std::array<RealArray,2> gpar;
     if (this->getGridParameters(gpar[0],0,npe[0]-1) &&
-	this->getGridParameters(gpar[1],1,npe[1]-1))
+        this->getGridParameters(gpar[1],1,npe[1]-1))
     {
       if (!project)
-	// Evaluate the secondary solution directly at all sampling points
-	return this->evalSolution(sField,integrand,gpar.data());
+        // Evaluate the secondary solution directly at all sampling points
+        return this->evalSolution(sField,integrand,gpar.data());
       else if (s)
       {
-	// Evaluate the projected field at the result sampling points
-	Go::Point p;
-	sField.resize(s->dimension(),gpar[0].size()*gpar[1].size());
+        // Evaluate the projected field at the result sampling points
+        Go::Point p;
+        sField.resize(s->dimension(),gpar[0].size()*gpar[1].size());
 
-	int iel = 0; // evaluation points are always structured in element order
-	for (size_t i = 0; i < gpar[0].size(); i++)
-	{
-	  if ((i+1)%npe[0] == 0) iel++;
-	  s->point(p,gpar[0][i],gpar[1][i],iel);
-	  sField.fillColumn(i+1,p.begin());
-	}
-	delete s;
-	return true;
+        int iel = 0; // evaluation points are always structured in element order
+        for (size_t i = 0; i < gpar[0].size(); i++)
+        {
+          if ((i+1)%npe[0] == 0) iel++;
+          s->point(p,gpar[0][i],gpar[1][i],iel);
+          sField.fillColumn(i+1,p.begin());
+        }
+        delete s;
+        return true;
       }
     }
     else if (s)
@@ -1789,13 +1789,13 @@ bool ASMu2D::updateDirichlet (const std::map<int,RealFunc*>& func,
     else
     {
       std::cerr <<" *** ASMu2D::updateDirichlet: Code "<< dirich[i].code
-		<<" is not associated with any function."<< std::endl;
+                <<" is not associated with any function."<< std::endl;
       return false;
     }
     if (!dcrv)
     {
       std::cerr <<" *** ASMu2D::updateDirichlet: Projection failure."
-		<< std::endl;
+                << std::endl;
       return false;
     }
 
@@ -1803,23 +1803,23 @@ bool ASMu2D::updateDirichlet (const std::map<int,RealFunc*>& func,
     for (nit = dirich[i].nodes.begin(); nit != dirich[i].nodes.end(); ++nit)
       for (int dofs = dirich[i].dof; dofs > 0; dofs /= 10)
       {
-	int dof = dofs%10;
-	// Find the constraint equation for current (node,dof)
-	MPC pDOF(MLGN[nit->second-1],dof);
-	MPCIter mit = mpcs.find(&pDOF);
-	if (mit == mpcs.end()) continue; // probably a deleted constraint
+        int dof = dofs%10;
+        // Find the constraint equation for current (node,dof)
+        MPC pDOF(MLGN[nit->second-1],dof);
+        MPCIter mit = mpcs.find(&pDOF);
+        if (mit == mpcs.end()) continue; // probably a deleted constraint
 
-	// Find index to the control point value for this (node,dof) in dcrv
-	RealArray::const_iterator cit = dcrv->coefs_begin();
-	if (dcrv->dimension() > 1) // A vector field is specified
-	  cit += (nit->first-1)*dcrv->dimension() + (dof-1);
-	else // A scalar field is specified at this dof
-	  cit += (nit->first-1);
+        // Find index to the control point value for this (node,dof) in dcrv
+        RealArray::const_iterator cit = dcrv->coefs_begin();
+        if (dcrv->dimension() > 1) // A vector field is specified
+          cit += (nit->first-1)*dcrv->dimension() + (dof-1);
+        else // A scalar field is specified at this dof
+          cit += (nit->first-1);
 
-	// Now update the prescribed value in the constraint equation
-	(*mit)->setSlaveCoeff(*cit);
+        // Now update the prescribed value in the constraint equation
+        (*mit)->setSlaveCoeff(*cit);
 #if SP_DEBUG > 1
-	std::cout <<"Updated constraint: "<< **mit;
+        std::cout <<"Updated constraint: "<< **mit;
 #endif
       }
     delete dcrv;
diff --git a/src/ASM/LR/ASMu2D.h b/src/ASM/LR/ASMu2D.h
index 97dc608d..b996fa5a 100644
--- a/src/ASM/LR/ASMu2D.h
+++ b/src/ASM/LR/ASMu2D.h
@@ -164,7 +164,7 @@ public:
   //! \param[in] project If \e true, the local axis directions are projected
   //! \return Number of additional nodes added due to local axis constraints
   virtual size_t constrainEdgeLocal(int dir, bool open, int dof, int code,
-				    bool project = false);
+                                    bool project = false);
 
   //! \brief Constrains a corner node identified by the two parameter indices.
   //! \param[in] I Parameter index in u-direction
@@ -358,8 +358,8 @@ public:
   //! \param[in] integrand Object with problem-specific data and methods
   //! \param[in] continuous If \e true, a continuous L2-projection is used
   virtual bool globalL2projection(Matrix& sField,
-				  const IntegrandBase& integrand,
-				  bool continuous = false) const;
+                                  const IntegrandBase& integrand,
+                                  bool continuous = false) const;
 
   //! \brief Transfers Gauss point variables from old basis to this patch.
   //! \param[in] oldBasis The LR-spline basis to transfer from
diff --git a/src/ASM/LR/ASMu2Dmx.C b/src/ASM/LR/ASMu2Dmx.C
index 9e374496..05206ae1 100644
--- a/src/ASM/LR/ASMu2Dmx.C
+++ b/src/ASM/LR/ASMu2Dmx.C
@@ -39,7 +39,7 @@
 
 
 ASMu2Dmx::ASMu2Dmx (unsigned char n_s,
-		    const std::vector<unsigned char>& n_f)
+                    const std::vector<unsigned char>& n_f)
   : ASMu2D(n_s), ASMmxBase(n_f)
 {
 }
@@ -149,14 +149,14 @@ void ASMu2Dmx::initMADOF (const int* sysMadof)
 
 
 void ASMu2Dmx::extractNodeVec (const Vector& globRes, Vector& nodeVec,
-			       unsigned char, int basis) const
+                               unsigned char, int basis) const
 {
   this->extractNodeVecMx(globRes,nodeVec,basis);
 }
 
 
 bool ASMu2Dmx::injectNodeVec (const Vector& nodeRes, Vector& globRes,
-			      unsigned char, int basis) const
+                              unsigned char, int basis) const
 {
   this->injectNodeVecMx(globRes,nodeRes,basis);
   return true;
@@ -164,7 +164,7 @@ bool ASMu2Dmx::injectNodeVec (const Vector& nodeRes, Vector& globRes,
 
 
 bool ASMu2Dmx::getSolution (Matrix& sField, const Vector& locSol,
-			    const IntVec& nodes) const
+                            const IntVec& nodes) const
 {
   return this->getSolutionMx(sField,locSol,nodes);
 }
@@ -244,8 +244,8 @@ bool ASMu2Dmx::generateFEMTopology ()
 
 
 bool ASMu2Dmx::integrate (Integrand& integrand,
-			  GlobalIntegral& glInt,
-			  const TimeDomain& time)
+                          GlobalIntegral& glInt,
+                          const TimeDomain& time)
 {
   if (m_basis.empty())
     return true; // silently ignore empty patches
@@ -309,38 +309,38 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
     for (int j = 0; j < nGauss; j++)
       for (int i = 0; i < nGauss; i++, fe.iGP++)
       {
-	// Local element coordinates of current integration point
-	fe.xi  = xg[i];
-	fe.eta = xg[j];
+        // Local element coordinates of current integration point
+        fe.xi  = xg[i];
+        fe.eta = xg[j];
 
-	// Parameter values of current integration point
-	fe.u = gpar[0][i];
-	fe.v = gpar[1][j];
+        // Parameter values of current integration point
+        fe.u = gpar[0][i];
+        fe.v = gpar[1][j];
 
-	// Compute basis function derivatives at current integration point
+        // Compute basis function derivatives at current integration point
         std::vector<Go::BasisDerivsSf> splinex(m_basis.size());
         for (size_t i=0; i < m_basis.size(); ++i) {
           m_basis[i]->computeBasis(fe.u, fe.v, splinex[i], els[i]-1);
           SplineUtils::extractBasis(splinex[i],fe.basis(i+1),dNxdu[i]);
         }
 
-	// Compute Jacobian inverse of coordinate mapping and derivatives
+        // Compute Jacobian inverse of coordinate mapping and derivatives
         // basis function derivatives w.r.t. Cartesian coordinates
         fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1]);
-	if (fe.detJxW == 0.0) continue; // skip singular points
+        if (fe.detJxW == 0.0) continue; // skip singular points
         for (size_t b = 0; b < m_basis.size(); ++b)
           if (b != (size_t)geoBasis-1)
             fe.grad(b+1).multiply(dNxdu[b],Jac);
 
-	// Cartesian coordinates of current integration point
+        // Cartesian coordinates of current integration point
         X = Xnod * fe.basis(geoBasis);
-	X.t = time.t;
+        X.t = time.t;
 
-	// Evaluate the integrand and accumulate element contributions
-	fe.detJxW *= 0.25*dA*wg[i]*wg[j];
-	PROFILE3("Integrand::evalInt");
-	if (!integrand.evalIntMx(*A,fe,time,X))
-	  return false;
+        // Evaluate the integrand and accumulate element contributions
+        fe.detJxW *= 0.25*dA*wg[i]*wg[j];
+        PROFILE3("Integrand::evalInt");
+        if (!integrand.evalIntMx(*A,fe,time,X))
+          return false;
       }
 
     // Finalize the element quantities
@@ -359,8 +359,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
 
 
 bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
-			  GlobalIntegral& glInt,
-			  const TimeDomain& time)
+                          GlobalIntegral& glInt,
+                          const TimeDomain& time)
 {
   if (!m_basis[0] || !m_basis[1])
     return true; // silently ignore empty patches
@@ -490,7 +490,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
       // Evaluate the integrand and accumulate element contributions
       fe.detJxW *= 0.5*dS*wg[i];
       if (!integrand.evalBouMx(*A,fe,time,X,normal))
-	return false;
+        return false;
     }
 
     // Finalize the element quantities
@@ -569,7 +569,7 @@ bool ASMu2Dmx::evalSolution (Matrix& sField, const Vector& locSol,
 
 
 bool ASMu2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
-			     const RealArray* gpar, bool regular) const
+                             const RealArray* gpar, bool regular) const
 {
   return evalSolution(sField,
                      const_cast<IntegrandBase&>(integrand).getSolution(0),
diff --git a/src/ASM/LR/ASMu2Dmx.h b/src/ASM/LR/ASMu2Dmx.h
index df3bebc1..19b1f495 100644
--- a/src/ASM/LR/ASMu2Dmx.h
+++ b/src/ASM/LR/ASMu2Dmx.h
@@ -88,7 +88,7 @@ public:
   //! \param glbInt The integrated quantity
   //! \param[in] time Parameters for nonlinear/time-dependent simulations
   virtual bool integrate(Integrand& integrand,
-			 GlobalIntegral& glbInt, const TimeDomain& time);
+                         GlobalIntegral& glbInt, const TimeDomain& time);
 
   //! \brief Evaluates a boundary integral over a patch edge.
   //! \param integrand Object with problem-specific data and methods
@@ -96,7 +96,7 @@ public:
   //! \param glbInt The integrated quantity
   //! \param[in] time Parameters for nonlinear/time-dependent simulations
   virtual bool integrate(Integrand& integrand, int lIndex,
-			 GlobalIntegral& glbInt, const TimeDomain& time);
+                         GlobalIntegral& glbInt, const TimeDomain& time);
 
 
   // Post-processing methods
@@ -107,7 +107,7 @@ public:
   //! \param[in] locSol Solution vector local to current patch
   //! \param[in] nodes 1-based local node numbers to extract solution for
   virtual bool getSolution(Matrix& sField, const Vector& locSol,
-			   const IntVec& nodes) const;
+                           const IntVec& nodes) const;
 
   //! \brief Evaluates the primary solution field at the given points.
   //! \param[out] sField Solution field
@@ -139,29 +139,29 @@ public:
   //! Otherwise, we assume that it contains the \a u and \a v parameters
   //! directly for each sampling point.
   virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
-			    const RealArray* gpar, bool regular = true) const;
+                            const RealArray* gpar, bool regular = true) const;
 
   //! \brief Extracts nodal results for this patch from the global vector.
   //! \param[in] globVec Global solution vector in DOF-order
   //! \param[out] nodeVec Nodal result vector for this patch
   //! \param[in] basis Which basis (or 0 for both) to extract nodal values for
   virtual void extractNodeVec(const Vector& globVec, Vector& nodeVec,
-			      unsigned char = 0, int basis = 0) const;
+                              unsigned char = 0, int basis = 0) const;
 
   //! \brief Inject nodal results for this patch into a global vector.
   //! \param[in] nodeVec Nodal result vector for this patch
   //! \param[out] globVec Global solution vector in DOF-order
   //! \param[in] basis Which basis (or 0 for both) to extract nodal values for
   virtual bool injectNodeVec(const Vector& nodeVec, Vector& globVec,
-			     unsigned char = 0, int basis = 0) const;
+                             unsigned char = 0, int basis = 0) const;
 
   //! \brief Projects the secondary solution using a discrete global L2-norm.
   //! \param[out] sField Secondary solution field control point values
   //! \param[in] integrand Object with problem-specific data and methods
   //! \param[in] continuous If \e true, a continuous L2-projection is used
   virtual bool globalL2projection(Matrix& sField,
-				  const IntegrandBase& integrand,
-				  bool continuous = false) const;
+                                  const IntegrandBase& integrand,
+                                  bool continuous = false) const;
 
   using ASMu2D::refine;
   //! \brief Refines the mesh adaptively.
diff --git a/src/ASM/LR/ASMu2Dmxrecovery.C b/src/ASM/LR/ASMu2Dmxrecovery.C
index 0c2a021c..839797f7 100644
--- a/src/ASM/LR/ASMu2Dmxrecovery.C
+++ b/src/ASM/LR/ASMu2Dmxrecovery.C
@@ -42,8 +42,8 @@ static void expandTensorGrid (const RealArray* in, RealArray* out)
 
 
 bool ASMu2Dmx::globalL2projection (Matrix& sField,
-				   const IntegrandBase& integrand,
-				   bool continuous) const
+                                   const IntegrandBase& integrand,
+                                   bool continuous) const
 {
   if (!m_basis[0] || !m_basis[1]) return true; // silently ignore empty patches
 
@@ -158,14 +158,14 @@ bool ASMu2Dmx::globalL2projection (Matrix& sField,
               int jnod = MNPC[iel-1][jj+el_ofs]+1;
               for (size_t k=1;k<=nfx[b];++k) {
                 A((inod-nod_ofs)*nfx[b]+k+eq_ofs,(jnod-nod_ofs)*nfx[b]+k+eq_ofs) += phi[b][ii]*phi[b][jj]*dJw;
-	   }
+           }
             }
             for (size_t k=1;k<=nfx[b];++k)
               B((inod-nod_ofs)*nfx[b]+k+eq_ofs) += phi[b][ii]*sField(k,ip+1)*dJw;
           }
           el_ofs += elem_sizes[b];
           eq_ofs += nb[b]*nfx[b];
-	  nod_ofs += nb[b];
+          nod_ofs += nb[b];
         }
       }
   }
diff --git a/src/ASM/LR/ASMu2Drecovery.C b/src/ASM/LR/ASMu2Drecovery.C
index dffd359a..e146911b 100644
--- a/src/ASM/LR/ASMu2Drecovery.C
+++ b/src/ASM/LR/ASMu2Drecovery.C
@@ -96,8 +96,8 @@ LR::LRSpline* ASMu2D::evalSolution (const IntegrandBase& integrand) const
 
 
 bool ASMu2D::globalL2projection (Matrix& sField,
-				 const IntegrandBase& integrand,
-				 bool continuous) const
+                                 const IntegrandBase& integrand,
+                                 bool continuous) const
 {
   if (!lrspline) return true; // silently ignore empty patches
 
@@ -339,29 +339,29 @@ LR::LRSplineSurface* ASMu2D::scRecovery (const IntegrandBase& integrand) const
       // Loop over the Gauss points in current knot-span
       int i, j, ig = 1;
       for (j = 0; j < ng2; j++)
-	for (i = 0; i < ng1; i++, ig++)
-	{
-	  // Evaluate the polynomial expansion at current Gauss point
-	  lrspline->point(X,gaussPt[0][i],gaussPt[1][j]);
-	  evalMonomials(n1,n2,X[0]-G[0],X[1]-G[1],P);
+        for (i = 0; i < ng1; i++, ig++)
+        {
+          // Evaluate the polynomial expansion at current Gauss point
+          lrspline->point(X,gaussPt[0][i],gaussPt[1][j]);
+          evalMonomials(n1,n2,X[0]-G[0],X[1]-G[1],P);
 #if SP_DEBUG > 2
-	  std::cout <<"Greville     point: "<< G
-		    <<"\nGauss        point: "<< gaussPt[0][i] <<", "<< gaussPt[0][j]
-		    <<"\nMapped gauss point: "<< X
-		    <<"\nP-matrix:"<< P <<"--------------------\n"<< std::endl;
+          std::cout <<"Greville     point: "<< G
+                    <<"\nGauss        point: "<< gaussPt[0][i] <<", "<< gaussPt[0][j]
+                    <<"\nMapped gauss point: "<< X
+                    <<"\nP-matrix:"<< P <<"--------------------\n"<< std::endl;
 #endif
 
-	  for (k = 1; k <= nPol; k++)
-	  {
-	    // Accumulate the projection matrix, A += P^t * P
-	    for (l = 1; l <= nPol; l++)
-	      A(k,l) += P(k)*P(l);
+          for (k = 1; k <= nPol; k++)
+          {
+            // Accumulate the projection matrix, A += P^t * P
+            for (l = 1; l <= nPol; l++)
+              A(k,l) += P(k)*P(l);
 
-	    // Accumulate the right-hand-side matrix, B += P^t * sigma
-	    for (l = 1; l <= nCmp; l++)
-	      B(k,l) += P(k)*sField(l,ig);
-	  }
-	}
+            // Accumulate the right-hand-side matrix, B += P^t * sigma
+            for (l = 1; l <= nCmp; l++)
+              B(k,l) += P(k)*sField(l,ig);
+          }
+        }
     }
 
 #if SP_DEBUG > 2
diff --git a/src/ASM/LR/ASMu3D.C b/src/ASM/LR/ASMu3D.C
index 0aa39a72..b32be6ce 100644
--- a/src/ASM/LR/ASMu3D.C
+++ b/src/ASM/LR/ASMu3D.C
@@ -55,55 +55,55 @@ ASMu3D::ASMu3D (const ASMu3D& patch, unsigned char n_f)
 
 bool ASMu3D::read (std::istream& is)
 {
-	if (shareFE) return true;
+  if (shareFE) return true;
         lrspline.reset();
 
-	// read inputfile as either an LRSpline file directly or a tensor product B-spline and convert
-	char firstline[256];
-	is.getline(firstline, 256);
-	if(strncmp(firstline, "# LRSPLINE", 10) == 0) {
+  // read inputfile as either an LRSpline file directly or a tensor product B-spline and convert
+  char firstline[256];
+  is.getline(firstline, 256);
+  if(strncmp(firstline, "# LRSPLINE", 10) == 0) {
                 lrspline.reset(new LR::LRSplineVolume());
-		is >> *lrspline;
-	} else { // probably a SplineVolume, so we'll read that and convert
-		tensorspline = new Go::SplineVolume();
-		is >> *tensorspline;
+    is >> *lrspline;
+  } else { // probably a SplineVolume, so we'll read that and convert
+    tensorspline = new Go::SplineVolume();
+    is >> *tensorspline;
                 lrspline.reset(new LR::LRSplineVolume(tensorspline));
-	}
+  }
 
-	// Eat white-space characters to see if there is more data to read
-	char c;
-	while (is.get(c))
-		if (!isspace(c)) {
-			is.putback(c);
-			break;
-		}
+  // Eat white-space characters to see if there is more data to read
+  char c;
+  while (is.get(c))
+    if (!isspace(c)) {
+      is.putback(c);
+      break;
+    }
 
-	if (!is.good() && !is.eof())
-	{
-		std::cerr <<" *** ASMu3D::read: Failure reading spline data"<< std::endl;
+  if (!is.good() && !is.eof())
+  {
+    std::cerr <<" *** ASMu3D::read: Failure reading spline data"<< std::endl;
                 lrspline.reset();
-		return false;
-	}
-	else if (lrspline->dimension() < 3)
-	{
-		std::cerr <<" *** ASMu3D::read: Invalid spline volume patch, dim="
-			  << lrspline->dimension() << std::endl;
+    return false;
+  }
+  else if (lrspline->dimension() < 3)
+  {
+    std::cerr <<" *** ASMu3D::read: Invalid spline volume patch, dim="
+              << lrspline->dimension() << std::endl;
                 lrspline.reset();
-		return false;
-	}
+    return false;
+  }
 
         geo = lrspline.get();
-	return true;
+  return true;
 }
 
 
 bool ASMu3D::write (std::ostream& os, int) const
 {
-	if (!lrspline) return false;
+  if (!lrspline) return false;
 
-	os << *lrspline;
+  os << *lrspline;
 
-	return os.good();
+  return os.good();
 }
 
 
@@ -133,67 +133,67 @@ bool ASMu3D::checkRightHandSystem ()
 
 bool ASMu3D::refine (int dir, const RealArray& xi)
 {
-	if (!tensorspline || dir < 0 || dir > 2 || xi.empty()) return false;
-	if (xi.front() < 0.0 || xi.back() > 1.0) return false;
-	if (shareFE) return true;
+  if (!tensorspline || dir < 0 || dir > 2 || xi.empty()) return false;
+  if (xi.front() < 0.0 || xi.back() > 1.0) return false;
+  if (shareFE) return true;
 
-	RealArray extraKnots;
-	RealArray::const_iterator uit = tensorspline->basis(dir).begin();
-	double ucurr, uprev = *(uit++);
-	while (uit != tensorspline->basis(dir).end())
-	{
-		ucurr = *(uit++);
-		if (ucurr > uprev)
-			for (size_t i = 0; i < xi.size(); i++)
-	if (i > 0 && xi[i] < xi[i-1])
-	  return false;
-	else
-	  extraKnots.push_back(ucurr*xi[i] + uprev*(1.0-xi[i]));
+  RealArray extraKnots;
+  RealArray::const_iterator uit = tensorspline->basis(dir).begin();
+  double ucurr, uprev = *(uit++);
+  while (uit != tensorspline->basis(dir).end())
+  {
+    ucurr = *(uit++);
+    if (ucurr > uprev)
+      for (size_t i = 0; i < xi.size(); i++)
+        if (i > 0 && xi[i] < xi[i-1])
+          return false;
+        else
+          extraKnots.push_back(ucurr*xi[i] + uprev*(1.0-xi[i]));
 
-		uprev = ucurr;
-	}
+    uprev = ucurr;
+  }
 
-	tensorspline->insertKnot(dir,extraKnots);
+  tensorspline->insertKnot(dir,extraKnots);
         lrspline.reset(new LR::LRSplineVolume(tensorspline));
         geo = lrspline.get();
-	return true;
+  return true;
 }
 
 bool ASMu3D::uniformRefine (int dir, int nInsert)
 {
-	if (!tensorspline || dir < 0 || dir > 2 || nInsert < 1) return false;
-	if (shareFE) return true;
+  if (!tensorspline || dir < 0 || dir > 2 || nInsert < 1) return false;
+  if (shareFE) return true;
 
-	RealArray extraKnots;
-	RealArray::const_iterator uit = tensorspline->basis(dir).begin();
-	double ucurr, uprev = *(uit++);
-	while (uit != tensorspline->basis(dir).end())
-	{
-		ucurr = *(uit++);
-		if (ucurr > uprev)
-			for (int i = 0; i < nInsert; i++)
-			{
-				double xi = (double)(i+1)/(double)(nInsert+1);
-				extraKnots.push_back(ucurr*xi + uprev*(1.0-xi));
-			}
-		uprev = ucurr;
-	}
+  RealArray extraKnots;
+  RealArray::const_iterator uit = tensorspline->basis(dir).begin();
+  double ucurr, uprev = *(uit++);
+  while (uit != tensorspline->basis(dir).end())
+  {
+    ucurr = *(uit++);
+    if (ucurr > uprev)
+      for (int i = 0; i < nInsert; i++)
+      {
+        double xi = (double)(i+1)/(double)(nInsert+1);
+        extraKnots.push_back(ucurr*xi + uprev*(1.0-xi));
+      }
+    uprev = ucurr;
+  }
 
-	tensorspline->insertKnot(dir,extraKnots);
-        lrspline.reset(new LR::LRSplineVolume(tensorspline));
-        geo = lrspline.get();
-	return true;
+  tensorspline->insertKnot(dir,extraKnots);
+  lrspline.reset(new LR::LRSplineVolume(tensorspline));
+  geo = lrspline.get();
+  return true;
 }
 
 bool ASMu3D::raiseOrder (int ru, int rv, int rw)
 {
-	if (!tensorspline) return false;
-	if (shareFE) return true;
+  if (!tensorspline) return false;
+  if (shareFE) return true;
 
-	tensorspline->raiseOrder(ru,rv,rw);
-        lrspline.reset(new LR::LRSplineVolume(tensorspline));
-        geo = lrspline.get();
-	return true;
+  tensorspline->raiseOrder(ru,rv,rw);
+  lrspline.reset(new LR::LRSplineVolume(tensorspline));
+  geo = lrspline.get();
+  return true;
 }
 
 bool ASMu3D::generateFEMTopology ()
@@ -269,222 +269,222 @@ bool ASMu3D::connectPatch (int face, ASMu3D& neighbor, int nface, int norient)
 
 
 bool ASMu3D::connectBasis (int face, ASMu3D& neighbor, int nface, int norient,
-				                   int basis, int slave, int master)
+                           int basis, int slave, int master)
 {
-	if (shareFE && neighbor.shareFE)
-		return true;
-	else if (shareFE || neighbor.shareFE)
-	{
-		std::cerr <<" *** ASMu3D::connectPatch: Logic error, cannot"
-	      <<" connect a sharedFE patch with an unshared one"<< std::endl;
-		return false;
-	}
+  if (shareFE && neighbor.shareFE)
+    return true;
+  else if (shareFE || neighbor.shareFE)
+  {
+    std::cerr <<" *** ASMu3D::connectPatch: Logic error, cannot"
+        <<" connect a sharedFE patch with an unshared one"<< std::endl;
+    return false;
+  }
 
-	// Set up the slave node numbers for this volume patch
+  // Set up the slave node numbers for this volume patch
 
-	int n1, n2, n3;
-	if (!this->getSize(n1,n2,n3,basis)) return false;
-	int node = slave+1, i1 = 0, i2 = 0;
+  int n1, n2, n3;
+  if (!this->getSize(n1,n2,n3,basis)) return false;
+  int node = slave+1, i1 = 0, i2 = 0;
 
-	switch (face)
-		{
-		case 2: // Positive I-direction
-			node += n1-1;
-		case 1: // Negative I-direction
-			i1 = n1;
-			n1 = n2;
-			n2 = n3;
-			break;
+  switch (face)
+    {
+    case 2: // Positive I-direction
+      node += n1-1;
+    case 1: // Negative I-direction
+      i1 = n1;
+      n1 = n2;
+      n2 = n3;
+      break;
 
-		case 4: // Positive J-direction
-			node += n1*(n2-1);
-		case 3: // Negative J-direction
-			i2 = n1*(n2-1);
-			i1 = 1;
-			n2 = n3;
-			break;
+    case 4: // Positive J-direction
+      node += n1*(n2-1);
+    case 3: // Negative J-direction
+      i2 = n1*(n2-1);
+      i1 = 1;
+      n2 = n3;
+      break;
 
-		case 6: // Positive K-direction
-			node += n1*n2*(n3-1);
-		case 5: // Negative K-direction
-			i1 = 1;
-			break;
+    case 6: // Positive K-direction
+      node += n1*n2*(n3-1);
+    case 5: // Negative K-direction
+      i1 = 1;
+      break;
 
-		default:
-			std::cerr <<" *** ASMu3D::connectPatch: Invalid slave face "
-		<< face << std::endl;
-			return false;
-		}
+    default:
+      std::cerr <<" *** ASMu3D::connectPatch: Invalid slave face "
+    << face << std::endl;
+      return false;
+    }
 
-	int i, j;
-	IntMat slaveNodes(n1,IntVec(n2,0));
-	for (j = 0; j < n2; j++, node += i2)
-		for (i = 0; i < n1; i++, node += i1)
-			slaveNodes[i][j] = node;
+  int i, j;
+  IntMat slaveNodes(n1,IntVec(n2,0));
+  for (j = 0; j < n2; j++, node += i2)
+    for (i = 0; i < n1; i++, node += i1)
+      slaveNodes[i][j] = node;
 
-	// Set up the master node numbers for the neighboring volume patch
+  // Set up the master node numbers for the neighboring volume patch
 
-	if (!neighbor.getSize(n1,n2,n3,basis)) return false;
-	node = master+1; i1 = i2 = 0;
+  if (!neighbor.getSize(n1,n2,n3,basis)) return false;
+  node = master+1; i1 = i2 = 0;
 
-	switch (nface)
-		{
-		case 2: // Positive I-direction
-			node += n1-1;
-		case 1: // Negative I-direction
-			i1 = n1;
-			n1 = n2;
-			n2 = n3;
-			break;
+  switch (nface)
+    {
+    case 2: // Positive I-direction
+      node += n1-1;
+    case 1: // Negative I-direction
+      i1 = n1;
+      n1 = n2;
+      n2 = n3;
+      break;
 
-		case 4: // Positive J-direction
-			node += n1*(n2-1);
-		case 3: // Negative J-direction
-			i2 = n1*(n2-1);
-			i1 = 1;
-			n2 = n3;
-			break;
+    case 4: // Positive J-direction
+      node += n1*(n2-1);
+    case 3: // Negative J-direction
+      i2 = n1*(n2-1);
+      i1 = 1;
+      n2 = n3;
+      break;
 
-		case 6: // Positive K-direction
-			node += n1*n2*(n3-1);
-		case 5: // Negative K-direction
-			i1 = 1;
-			break;
+    case 6: // Positive K-direction
+      node += n1*n2*(n3-1);
+    case 5: // Negative K-direction
+      i1 = 1;
+      break;
 
-		default:
-			std::cerr <<" *** ASMu3D::connectPatch: Invalid master face "
-		<< nface << std::endl;
-			return false;
-		}
+    default:
+      std::cerr <<" *** ASMu3D::connectPatch: Invalid master face "
+    << nface << std::endl;
+      return false;
+    }
 
-	if (norient < 0 || norient > 7)
-	{
-		std::cerr <<" *** ASMu3D::connectPatch: Orientation flag "
-	      << norient <<" is out of range [0,7]"<< std::endl;
-		return false;
-	}
+  if (norient < 0 || norient > 7)
+  {
+    std::cerr <<" *** ASMu3D::connectPatch: Orientation flag "
+        << norient <<" is out of range [0,7]"<< std::endl;
+    return false;
+  }
 
-	int m1 = slaveNodes.size();
-	int m2 = slaveNodes.front().size();
-	if (norient < 4 ? (n1 != m1 || n2 != m2) : (n2 != m1 || n1 != m2))
-	{
-		std::cerr <<" *** ASMu3D::connectPatch: Non-matching faces, sizes "
-	      << n1 <<","<< n2 <<" and "<< m1 <<","<< m2 << std::endl;
-		return false;
-	}
+  int m1 = slaveNodes.size();
+  int m2 = slaveNodes.front().size();
+  if (norient < 4 ? (n1 != m1 || n2 != m2) : (n2 != m1 || n1 != m2))
+  {
+    std::cerr <<" *** ASMu3D::connectPatch: Non-matching faces, sizes "
+        << n1 <<","<< n2 <<" and "<< m1 <<","<< m2 << std::endl;
+    return false;
+  }
 
-	const double xtol = 1.0e-4;
-	for (j = 0; j < n2; j++, node += i2)
-		for (i = 0; i < n1; i++, node += i1)
-		{
-			int k = i, l = j;
-			switch (norient)
-	{
-	case  1: k =    i  ; l = n2-j-1; break;
-	case  2: k = n1-i-1; l =    j  ; break;
-	case  3: k = n1-i-1; l = n2-j-1; break;
-	case  4: k =    j  ; l =    i  ; break;
-	case  5: k =    j  ; l = n1-i-1; break;
-	case  6: k = n2-j-1; l =    i  ; break;
-	case  7: k = n2-j-1; l = n1-i-1; break;
-	default: k =    i  ; l = j     ;
-	}
+  const double xtol = 1.0e-4;
+  for (j = 0; j < n2; j++, node += i2)
+    for (i = 0; i < n1; i++, node += i1)
+    {
+      int k = i, l = j;
+      switch (norient)
+  {
+  case  1: k =    i  ; l = n2-j-1; break;
+  case  2: k = n1-i-1; l =    j  ; break;
+  case  3: k = n1-i-1; l = n2-j-1; break;
+  case  4: k =    j  ; l =    i  ; break;
+  case  5: k =    j  ; l = n1-i-1; break;
+  case  6: k = n2-j-1; l =    i  ; break;
+  case  7: k = n2-j-1; l = n1-i-1; break;
+  default: k =    i  ; l = j     ;
+  }
 
-			int slave = slaveNodes[k][l];
-			if (!neighbor.getCoord(node).equal(this->getCoord(slave),xtol))
-			{
-	std::cerr <<" *** ASMu3D::connectPatch: Non-matching nodes "
-		  << node <<": "<< neighbor.getCoord(node)
-		  <<"\n                                          and "
-		  << slave <<": "<< this->getCoord(slave) << std::endl;
-	return false;
-			}
-			else
-	ASMbase::collapseNodes(neighbor,node,*this,slave);
-		}
+      int slave = slaveNodes[k][l];
+      if (!neighbor.getCoord(node).equal(this->getCoord(slave),xtol))
+      {
+  std::cerr <<" *** ASMu3D::connectPatch: Non-matching nodes "
+      << node <<": "<< neighbor.getCoord(node)
+      <<"\n                                          and "
+      << slave <<": "<< this->getCoord(slave) << std::endl;
+  return false;
+      }
+      else
+  ASMbase::collapseNodes(neighbor,node,*this,slave);
+    }
 
-	return true;
+  return true;
 }
 
 
 void ASMu3D::closeFaces (int dir, int basis, int master)
 {
-	int n1, n2, n3;
-	if (basis < 1) basis = 1;
-	if (!this->getSize(n1,n2,n3,basis)) return;
+  int n1, n2, n3;
+  if (basis < 1) basis = 1;
+  if (!this->getSize(n1,n2,n3,basis)) return;
 
-	switch (dir)
-		{
-		case 1: // Faces are closed in I-direction
-			for (int i3 = 1; i3 <= n3; i3++)
-	for (int i2 = 1; i2 <= n2; i2++, master += n1)
-	  this->makePeriodic(master,master+n1-1);
-			break;
+  switch (dir)
+    {
+    case 1: // Faces are closed in I-direction
+      for (int i3 = 1; i3 <= n3; i3++)
+  for (int i2 = 1; i2 <= n2; i2++, master += n1)
+    this->makePeriodic(master,master+n1-1);
+      break;
 
-		case 2: // Faces are closed in J-direction
-			for (int i3 = 1; i3 <= n3; i3++, master += n1*(n2-1))
-	for (int i1 = 1; i1 <= n1; i1++, master++)
-	  this->makePeriodic(master,master+n1*(n2-1));
-			break;
+    case 2: // Faces are closed in J-direction
+      for (int i3 = 1; i3 <= n3; i3++, master += n1*(n2-1))
+  for (int i1 = 1; i1 <= n1; i1++, master++)
+    this->makePeriodic(master,master+n1*(n2-1));
+      break;
 
-		case 3: // Faces are closed in K-direction
-			for (int i2 = 1; i2 <= n2; i2++)
-	for (int i1 = 1; i1 <= n1; i1++, master++)
-	  this->makePeriodic(master,master+n1*n2*(n3-1));
-			break;
-		}
+    case 3: // Faces are closed in K-direction
+      for (int i2 = 1; i2 <= n2; i2++)
+  for (int i1 = 1; i1 <= n1; i1++, master++)
+    this->makePeriodic(master,master+n1*n2*(n3-1));
+      break;
+    }
 }
 */
 
 
 /*!
-	A negative \a code value implies direct evaluation of the Dirichlet condition
-	function at the control point. Positive \a code implies projection onto the
-	spline basis representing the boundary surface (needed for curved faces and/or
-	non-constant functions).
+  A negative \a code value implies direct evaluation of the Dirichlet condition
+  function at the control point. Positive \a code implies projection onto the
+  spline basis representing the boundary surface (needed for curved faces and/or
+  non-constant functions).
 */
 
 void ASMu3D::constrainFace (int dir, bool open, int dof, int code, char)
 {
-	if(open)
-		std::cerr << "\nWARNING: ASMu3D::constrainFace, open boundary conditions not supported yet. Treating it as closed" << std::endl;
+  if(open)
+    std::cerr << "\nWARNING: ASMu3D::constrainFace, open boundary conditions not supported yet. Treating it as closed" << std::endl;
 
-	int bcode = code;
-	if (code > 0) {// Dirichlet projection will be performed
-		std::cerr << "\nWARNING: Projective (nonhomogenuous) dirichlet boundary conditions not implemented. ";
-		std::cerr << "\n         Performing variational diminishing approximation instead" << std::endl;
-		// dirich.push_back(DirichletFace(this->getBoundary(dir),dof,code));
-	}
-	else if (code < 0)
-		bcode = -code;
+  int bcode = code;
+  if (code > 0) {// Dirichlet projection will be performed
+    std::cerr << "\nWARNING: Projective (nonhomogenuous) dirichlet boundary conditions not implemented. ";
+    std::cerr << "\n         Performing variational diminishing approximation instead" << std::endl;
+    // dirich.push_back(DirichletFace(this->getBoundary(dir),dof,code));
+  }
+  else if (code < 0)
+    bcode = -code;
 
-	// get all the boundary functions from the LRspline object
-	std::vector<LR::Basisfunction*> thisEdge;
-	if(dir == -1) 
-		lrspline->getEdgeFunctions(thisEdge, LR::WEST, 1);
-	else if(dir == 1) 
-		lrspline->getEdgeFunctions(thisEdge, LR::EAST, 1);
-	else if(dir == -2) 
-		lrspline->getEdgeFunctions(thisEdge, LR::SOUTH, 1);
-	else if(dir == 2) 
-		lrspline->getEdgeFunctions(thisEdge, LR::NORTH, 1);
-	else if(dir == -3) 
-		lrspline->getEdgeFunctions(thisEdge, LR::BOTTOM, 1);
-	else if(dir == 3) 
-		lrspline->getEdgeFunctions(thisEdge, LR::TOP, 1);
+  // get all the boundary functions from the LRspline object
+  std::vector<LR::Basisfunction*> thisEdge;
+  if(dir == -1) 
+    lrspline->getEdgeFunctions(thisEdge, LR::WEST, 1);
+  else if(dir == 1) 
+    lrspline->getEdgeFunctions(thisEdge, LR::EAST, 1);
+  else if(dir == -2) 
+    lrspline->getEdgeFunctions(thisEdge, LR::SOUTH, 1);
+  else if(dir == 2) 
+    lrspline->getEdgeFunctions(thisEdge, LR::NORTH, 1);
+  else if(dir == -3) 
+    lrspline->getEdgeFunctions(thisEdge, LR::BOTTOM, 1);
+  else if(dir == 3) 
+    lrspline->getEdgeFunctions(thisEdge, LR::TOP, 1);
 
-	std::cout << "\nNumber of constraints: " << thisEdge.size() << std::endl;
+  std::cout << "\nNumber of constraints: " << thisEdge.size() << std::endl;
 
-	// enforce the boundary conditions
-	for(LR::Basisfunction* b  : thisEdge)
-		this->prescribe(b->getId()+1,dof,bcode);
+  // enforce the boundary conditions
+  for(LR::Basisfunction* b  : thisEdge)
+    this->prescribe(b->getId()+1,dof,bcode);
 }
 
 size_t ASMu3D::constrainFaceLocal(int dir, bool open, int dof, int code, bool project, char T1)
 {
-	std::cerr << "ASMu3D::constrainFaceLocal not implemented properly yet" << std::endl;
-	exit(776654);
-	return 0;
+  std::cerr << "ASMu3D::constrainFaceLocal not implemented properly yet" << std::endl;
+  exit(776654);
+  return 0;
 }
 
 
@@ -545,73 +545,73 @@ void ASMu3D::constrainEdge (int lEdge, bool open, int dof, int code, char)
 void ASMu3D::constrainLine (int fdir, int ldir, double xi, int dof,
                             int code, char)
 {
-	std::cerr << "ASMu3D::constrainLine not implemented properly yet" << std::endl;
-	exit(776654);
+  std::cerr << "ASMu3D::constrainLine not implemented properly yet" << std::endl;
+  exit(776654);
 #if 0
-	if (xi < 0.0 || xi > 1.0) return;
+  if (xi < 0.0 || xi > 1.0) return;
 
-	int n1, n2, n3, node = 1;
-	if (!this->getSize(n1,n2,n3,1)) return;
+  int n1, n2, n3, node = 1;
+  if (!this->getSize(n1,n2,n3,1)) return;
 
-	switch (fdir)
-		{
-		case  1: // Right face (positive I-direction)
-			node += n1-1;
-		case -1: // Left face (negative I-direction)
-			if (ldir == 2)
-			{
-	// Line goes in J-direction
-	node += n1*n2*int(0.5+(n3-1)*xi);
-	for (int i2 = 1; i2 <= n2; i2++, node += n1)
-	  this->prescribe(node,dof,code);
-			}
-			else if (ldir == 3)
-			{
-	// Line goes in K-direction
-	node += n1*int(0.5+(n2-1)*xi);
-	for (int i3 = 1; i3 <= n3; i3++, node += n1*n2)
-	  this->prescribe(node,dof,code);
-			}
-			break;
+  switch (fdir)
+    {
+    case  1: // Right face (positive I-direction)
+      node += n1-1;
+    case -1: // Left face (negative I-direction)
+      if (ldir == 2)
+      {
+  // Line goes in J-direction
+  node += n1*n2*int(0.5+(n3-1)*xi);
+  for (int i2 = 1; i2 <= n2; i2++, node += n1)
+    this->prescribe(node,dof,code);
+      }
+      else if (ldir == 3)
+      {
+  // Line goes in K-direction
+  node += n1*int(0.5+(n2-1)*xi);
+  for (int i3 = 1; i3 <= n3; i3++, node += n1*n2)
+    this->prescribe(node,dof,code);
+      }
+      break;
 
-		case  2: // Back face (positive J-direction)
-			node += n1*(n2-1);
-		case -2: // Front face (negative J-direction)
-			if (ldir == 1)
-			{
-	// Line goes in I-direction
-	node += n1*n2*int(0.5+(n3-1)*xi);
-	for (int i1 = 1; i1 <= n1; i1++, node++)
-	  this->prescribe(node,dof,code);
-			}
-			else if (ldir == 3)
-			{
-	// Line goes in K-direction
-	node += int(0.5+(n1-1)*xi);
-	for (int i3 = 1; i3 <= n3; i3++, node += n1*n2)
-	  this->prescribe(node,dof,code);
-			}
-			break;
+    case  2: // Back face (positive J-direction)
+      node += n1*(n2-1);
+    case -2: // Front face (negative J-direction)
+      if (ldir == 1)
+      {
+  // Line goes in I-direction
+  node += n1*n2*int(0.5+(n3-1)*xi);
+  for (int i1 = 1; i1 <= n1; i1++, node++)
+    this->prescribe(node,dof,code);
+      }
+      else if (ldir == 3)
+      {
+  // Line goes in K-direction
+  node += int(0.5+(n1-1)*xi);
+  for (int i3 = 1; i3 <= n3; i3++, node += n1*n2)
+    this->prescribe(node,dof,code);
+      }
+      break;
 
-		case  3: // Top face (positive K-direction)
-			node += n1*n2*(n3-1);
-		case -3: // Bottom face (negative K-direction)
-			if (ldir == 1)
-			{
-	// Line goes in I-direction
-	node += n1*int(0.5+(n2-1)*xi);
-	for (int i1 = 1; i1 <= n1; i1++, node++)
-	  this->prescribe(node,dof,code);
-			}
-			else if (ldir == 2)
-			{
-	// Line goes in J-direction
-	node += int(0.5+(n1-1)*xi);
-	for (int i2 = 1; i2 <= n2; i2++, node += n1)
-	  this->prescribe(node,dof,code);
-			}
-			break;
-		}
+    case  3: // Top face (positive K-direction)
+      node += n1*n2*(n3-1);
+    case -3: // Bottom face (negative K-direction)
+      if (ldir == 1)
+      {
+  // Line goes in I-direction
+  node += n1*int(0.5+(n2-1)*xi);
+  for (int i1 = 1; i1 <= n1; i1++, node++)
+    this->prescribe(node,dof,code);
+      }
+      else if (ldir == 2)
+      {
+  // Line goes in J-direction
+  node += int(0.5+(n1-1)*xi);
+  for (int i2 = 1; i2 <= n2; i2++, node += n1)
+    this->prescribe(node,dof,code);
+      }
+      break;
+    }
 #endif
 }
 
@@ -635,24 +635,24 @@ void ASMu3D::constrainCorner (int I, int J, int K, int dof, int code, char)
 
 
 void ASMu3D::constrainNode (double xi, double eta, double zeta,
-			    int dof, int code, char)
+                            int dof, int code, char)
 {
-	std::cerr << "ASMu3D::constrainNode not implemented properly yet" << std::endl;
-	exit(776654);
-	if (xi   < 0.0 || xi   > 1.0) return;
-	if (eta  < 0.0 || eta  > 1.0) return;
-	if (zeta < 0.0 || zeta > 1.0) return;
+  std::cerr << "ASMu3D::constrainNode not implemented properly yet" << std::endl;
+  exit(776654);
+  if (xi   < 0.0 || xi   > 1.0) return;
+  if (eta  < 0.0 || eta  > 1.0) return;
+  if (zeta < 0.0 || zeta > 1.0) return;
 
 #if 0
-	int n1, n2, n3;
-	if (!this->getSize(n1,n2,n3,1)) return;
+  int n1, n2, n3;
+  if (!this->getSize(n1,n2,n3,1)) return;
 
-	int node = 1;
-	if (xi   > 0.0) node += int(0.5+(n1-1)*xi);
-	if (eta  > 0.0) node += n1*int(0.5+(n2-1)*eta);
-	if (zeta > 0.0) node += n1*n2*int(0.5+(n3-1)*zeta);
+  int node = 1;
+  if (xi   > 0.0) node += int(0.5+(n1-1)*xi);
+  if (eta  > 0.0) node += n1*int(0.5+(n2-1)*eta);
+  if (zeta > 0.0) node += n1*n2*int(0.5+(n3-1)*zeta);
 
-	this->prescribe(node,dof,code);
+  this->prescribe(node,dof,code);
 #endif
 }
 
@@ -662,30 +662,30 @@ void ASMu3D::constrainNode (double xi, double eta, double zeta,
 double ASMu3D::getParametricArea (int iel, int dir) const
 {
 #ifdef INDEX_CHECK
-	if (iel < 1 || (size_t)iel > MNPC.size())
-	{
-		std::cerr <<" *** ASMu3D::getParametricArea: Element index "<< iel
-	      <<" out of range [1,"<< MNPC.size() <<"]."<< std::endl;
-		return DERR;
-	}
+  if (iel < 1 || (size_t)iel > MNPC.size())
+  {
+    std::cerr <<" *** ASMu3D::getParametricArea: Element index "<< iel
+        <<" out of range [1,"<< MNPC.size() <<"]."<< std::endl;
+    return DERR;
+  }
 #endif
-	if (MNPC[iel-1].empty())
-		return 0.0;
+  if (MNPC[iel-1].empty())
+    return 0.0;
 
-	LR::Element *el = lrspline->getElement(iel-1);
-	double du = el->getParmax(0) - el->getParmin(0);
-	double dv = el->getParmax(1) - el->getParmin(1);
-	double dw = el->getParmax(2) - el->getParmin(2);
-	switch (dir)
-		{
-		case 1: return dv * dw;
-		case 2: return du * dw;
-		case 3: return du * dv;
-		}
+  LR::Element *el = lrspline->getElement(iel-1);
+  double du = el->getParmax(0) - el->getParmin(0);
+  double dv = el->getParmax(1) - el->getParmin(1);
+  double dw = el->getParmax(2) - el->getParmin(2);
+  switch (dir)
+    {
+    case 1: return dv * dw;
+    case 2: return du * dw;
+    case 3: return du * dv;
+    }
 
-	std::cerr <<" *** ASMu3D::getParametricArea: Invalid face direction "
-	    << dir << std::endl;
-	return DERR;
+  std::cerr <<" *** ASMu3D::getParametricArea: Invalid face direction "
+            << dir << std::endl;
+  return DERR;
 }
 
 double ASMu3D::getParametricVolume (int iel) const
@@ -761,29 +761,29 @@ bool ASMu3D::updateCoords (const Vector& displ)
 
 size_t ASMu3D::getNoBoundaryElms (char lIndex, char ldim) const
 {
-	if (!lrspline) return 0;
+  if (!lrspline) return 0;
 
-	LR::parameterEdge edge;
-	switch(lIndex)
-	{
-	case 1: edge = LR::WEST;   break;
-	case 2: edge = LR::EAST;   break;
-	case 3: edge = LR::SOUTH;  break;
-	case 4: edge = LR::NORTH;  break;
-	case 5: edge = LR::BOTTOM; break;
-	case 6: edge = LR::TOP;    break;
-	default:edge = LR::NONE;
-	}
+  LR::parameterEdge edge;
+  switch(lIndex)
+  {
+  case 1: edge = LR::WEST;   break;
+  case 2: edge = LR::EAST;   break;
+  case 3: edge = LR::SOUTH;  break;
+  case 4: edge = LR::NORTH;  break;
+  case 5: edge = LR::BOTTOM; break;
+  case 6: edge = LR::TOP;    break;
+  default:edge = LR::NONE;
+  }
 
-	std::vector<LR::Element*> edgeElms;
-	lrspline->getEdgeElements(edgeElms, edge);
+  std::vector<LR::Element*> edgeElms;
+  lrspline->getEdgeElements(edgeElms, edge);
 
-	return edgeElms.size();
+  return edgeElms.size();
 }
 
 
 void ASMu3D::getGaussPointParameters (RealArray& uGP, int dir, int nGauss,
-				      int iEl, const double* xi) const
+                                      int iEl, const double* xi) const
 {
   LR::Element* el = lrspline->getElement(iEl);
   double start = el->getParmin(dir);
@@ -810,415 +810,415 @@ void ASMu3D::getElementCorners (int iEl, Vec3Vec& XC) const
     for (int j = 0; j < 2; j++)
       for (int i = 0; i < 2; i++)
       {
-	lrspline->point(pt,u[i],v[j],w[k],iEl);
-	XC.push_back(SplineUtils::toVec3(pt));
+        lrspline->point(pt,u[i],v[j],w[k],iEl);
+        XC.push_back(SplineUtils::toVec3(pt));
       }
 }
 
 
 void ASMu3D::evaluateBasis (FiniteElement &el, Matrix &dNdu) const
 {
-	PROFILE2("Spline evaluation");
-	size_t nBasis = lrspline->getElement(el.iel-1)->nBasisFunctions();
+  PROFILE2("Spline evaluation");
+  size_t nBasis = lrspline->getElement(el.iel-1)->nBasisFunctions();
 
-	std::vector<RealArray> result;
-	lrspline->computeBasis(el.u, el.v, el.w, result, 1, el.iel-1);
+  std::vector<RealArray> result;
+  lrspline->computeBasis(el.u, el.v, el.w, result, 1, el.iel-1);
 
-	el.N.resize(nBasis);
-	dNdu.resize(nBasis,3);
-	size_t jp, n = 1;
-	for (jp = 0; jp < nBasis; jp++, n++) {
-		el.N(n)   = result[jp][0];
-		dNdu(n,1) = result[jp][1];
-		dNdu(n,2) = result[jp][2];
-		dNdu(n,3) = result[jp][3];
-	}
+  el.N.resize(nBasis);
+  dNdu.resize(nBasis,3);
+  size_t jp, n = 1;
+  for (jp = 0; jp < nBasis; jp++, n++) {
+    el.N(n)   = result[jp][0];
+    dNdu(n,1) = result[jp][1];
+    dNdu(n,2) = result[jp][2];
+    dNdu(n,3) = result[jp][3];
+  }
 }
 
 void ASMu3D::evaluateBasis (FiniteElement &el, Matrix &dNdu, Matrix &C, Matrix &B) const
 {
-	PROFILE2("BeSpline evaluation");
-	Matrix N = C*B;
-	size_t nBasis = lrspline->getElement(el.iel-1)->nBasisFunctions();
-	el.N = N.getColumn(1);
-	dNdu.resize(nBasis,3);
-	dNdu.fillColumn(1,N.getColumn(2));
-	dNdu.fillColumn(2,N.getColumn(3));
-	dNdu.fillColumn(3,N.getColumn(4));
+  PROFILE2("BeSpline evaluation");
+  Matrix N = C*B;
+  size_t nBasis = lrspline->getElement(el.iel-1)->nBasisFunctions();
+  el.N = N.getColumn(1);
+  dNdu.resize(nBasis,3);
+  dNdu.fillColumn(1,N.getColumn(2));
+  dNdu.fillColumn(2,N.getColumn(3));
+  dNdu.fillColumn(3,N.getColumn(4));
 }
 
 void ASMu3D::evaluateBasis (FiniteElement &el, Matrix &dNdu, Matrix3D &d2Ndu2) const
 {
-	PROFILE2("Spline evaluation");
-	size_t nBasis = lrspline->getElement(el.iel-1)->nBasisFunctions();
+  PROFILE2("Spline evaluation");
+  size_t nBasis = lrspline->getElement(el.iel-1)->nBasisFunctions();
 
-	std::vector<RealArray> result;
-	lrspline->computeBasis(el.u, el.v, el.w, result, 2, el.iel-1 );
+  std::vector<RealArray> result;
+  lrspline->computeBasis(el.u, el.v, el.w, result, 2, el.iel-1 );
 
-	el.N.resize(nBasis);
-	dNdu.resize(nBasis,3);
-	d2Ndu2.resize(nBasis,3,3);
-	size_t jp, n = 1;
-	for (jp = 0; jp < nBasis; jp++, n++) {
-		el.N   (n)    = result[jp][0];
-		dNdu (n,1)    = result[jp][1];
-		dNdu (n,2)    = result[jp][2];
-		dNdu (n,3)    = result[jp][3];
-		d2Ndu2(n,1,1) = result[jp][4];
-		d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = result[jp][5];
-		d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = result[jp][6];
-		d2Ndu2(n,2,2) = result[jp][7];
-		d2Ndu2(n,2,3) = d2Ndu2(n,3,2) = result[jp][8];
-		d2Ndu2(n,3,3) = result[jp][9];
-	}
+  el.N.resize(nBasis);
+  dNdu.resize(nBasis,3);
+  d2Ndu2.resize(nBasis,3,3);
+  size_t jp, n = 1;
+  for (jp = 0; jp < nBasis; jp++, n++) {
+    el.N   (n)    = result[jp][0];
+    dNdu (n,1)    = result[jp][1];
+    dNdu (n,2)    = result[jp][2];
+    dNdu (n,3)    = result[jp][3];
+    d2Ndu2(n,1,1) = result[jp][4];
+    d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = result[jp][5];
+    d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = result[jp][6];
+    d2Ndu2(n,2,2) = result[jp][7];
+    d2Ndu2(n,2,3) = d2Ndu2(n,3,2) = result[jp][8];
+    d2Ndu2(n,3,3) = result[jp][9];
+  }
 }
 
 void ASMu3D::evaluateBasis (FiniteElement &el, int derivs) const
 {
-	PROFILE2("Spline evaluation");
-	size_t nBasis = lrspline->getElement(el.iel-1)->nBasisFunctions();
+  PROFILE2("Spline evaluation");
+  size_t nBasis = lrspline->getElement(el.iel-1)->nBasisFunctions();
 
-	std::vector<RealArray> result;
-	lrspline->computeBasis(el.u, el.v, el.w, result, derivs, el.iel-1 );
+  std::vector<RealArray> result;
+  lrspline->computeBasis(el.u, el.v, el.w, result, derivs, el.iel-1 );
 
-	el.N.resize(nBasis);
-	el.dNdX.resize(nBasis,3);
-	el.d2NdX2.resize(nBasis,3,3);
-	size_t jp, n = 1;
-	for (jp = 0; jp < nBasis; jp++, n++) {
-		el.N   (n)      = result[jp][0];
-		if(derivs > 0) {  
-			el.dNdX (n,1)    = result[jp][1];
-			el.dNdX (n,2)    = result[jp][2];
-			el.dNdX (n,3)    = result[jp][3];
-		}
-		if(derivs > 1) {
-			el.d2NdX2(n,1,1) = result[jp][4];
-			el.d2NdX2(n,1,2) = el.d2NdX2(n,2,1) = result[jp][5];
-			el.d2NdX2(n,1,3) = el.d2NdX2(n,3,1) = result[jp][6];
-			el.d2NdX2(n,2,2) = result[jp][7];
-			el.d2NdX2(n,2,3) = el.d2NdX2(n,3,2) = result[jp][8];
-			el.d2NdX2(n,3,3) = result[jp][9];
-		}
-	}
+  el.N.resize(nBasis);
+  el.dNdX.resize(nBasis,3);
+  el.d2NdX2.resize(nBasis,3,3);
+  size_t jp, n = 1;
+  for (jp = 0; jp < nBasis; jp++, n++) {
+    el.N   (n)      = result[jp][0];
+    if(derivs > 0) {  
+      el.dNdX (n,1)    = result[jp][1];
+      el.dNdX (n,2)    = result[jp][2];
+      el.dNdX (n,3)    = result[jp][3];
+    }
+    if(derivs > 1) {
+      el.d2NdX2(n,1,1) = result[jp][4];
+      el.d2NdX2(n,1,2) = el.d2NdX2(n,2,1) = result[jp][5];
+      el.d2NdX2(n,1,3) = el.d2NdX2(n,3,1) = result[jp][6];
+      el.d2NdX2(n,2,2) = result[jp][7];
+      el.d2NdX2(n,2,3) = el.d2NdX2(n,3,2) = result[jp][8];
+      el.d2NdX2(n,3,3) = result[jp][9];
+    }
+  }
 }
 
 bool ASMu3D::integrate (Integrand& integrand,
                         GlobalIntegral& glInt,
                         const TimeDomain& time)
 {
-	if (!lrspline) return true; // silently ignore empty patches
+  if (!lrspline) return true; // silently ignore empty patches
 
-	PROFILE2("ASMu3D::integrate(I)");
+  PROFILE2("ASMu3D::integrate(I)");
 
-	// Get Gaussian quadrature points and weights
-	const double* xg = GaussQuadrature::getCoord(nGauss);
-	const double* wg = GaussQuadrature::getWeight(nGauss);
-	if (!xg || !wg) return false;
+  // Get Gaussian quadrature points and weights
+  const double* xg = GaussQuadrature::getCoord(nGauss);
+  const double* wg = GaussQuadrature::getWeight(nGauss);
+  if (!xg || !wg) return false;
 
-	// evaluate all gauss points on the bezier patch (-1, 1)
-	int p1 = lrspline->order(0);
-	int p2 = lrspline->order(1);
-	int p3 = lrspline->order(2);
-	Go::BsplineBasis basis1 = getBezierBasis(p1);
-	Go::BsplineBasis basis2 = getBezierBasis(p2);
-	Go::BsplineBasis basis3 = getBezierBasis(p3);
+  // evaluate all gauss points on the bezier patch (-1, 1)
+  int p1 = lrspline->order(0);
+  int p2 = lrspline->order(1);
+  int p3 = lrspline->order(2);
+  Go::BsplineBasis basis1 = getBezierBasis(p1);
+  Go::BsplineBasis basis2 = getBezierBasis(p2);
+  Go::BsplineBasis basis3 = getBezierBasis(p3);
 
-	Matrix BN(   p1*p2*p3, nGauss*nGauss*nGauss);
-	Matrix BdNdu(p1*p2*p3, nGauss*nGauss*nGauss);
-	Matrix BdNdv(p1*p2*p3, nGauss*nGauss*nGauss);
-	Matrix BdNdw(p1*p2*p3, nGauss*nGauss*nGauss);
-	int ig=1; // gauss point iterator
-	for(int zeta=0; zeta<nGauss; zeta++) {
-		for(int eta=0; eta<nGauss; eta++) {
-			for(int xi=0; xi<nGauss; xi++, ig++) {
-				double u[2*p1];
-				double v[2*p2];
-				double w[2*p3];
-				basis1.computeBasisValues(xg[xi],   u, 1); 
-				basis2.computeBasisValues(xg[eta],  v, 1);
-				basis3.computeBasisValues(xg[zeta], w, 1);
-				int ib=1; // basis function iterator
-				double sum = 0;
-				for(int k=0; k<p3; k++) {
-					for(int j=0; j<p2; j++) {
-						for(int i=0; i<p1; i++, ib++) {
-							BN(ib,ig)    = u[2*i  ]*v[2*j  ]*w[2*k  ];
-							BdNdu(ib,ig) = u[2*i+1]*v[2*j  ]*w[2*k  ];
-							BdNdv(ib,ig) = u[2*i  ]*v[2*j+1]*w[2*k  ];
-							BdNdw(ib,ig) = u[2*i  ]*v[2*j  ]*w[2*k+1];
-							sum += BN(ib,ig);
-						}
-					}
-				}
-			}
-		}
-	}
+  Matrix BN(   p1*p2*p3, nGauss*nGauss*nGauss);
+  Matrix BdNdu(p1*p2*p3, nGauss*nGauss*nGauss);
+  Matrix BdNdv(p1*p2*p3, nGauss*nGauss*nGauss);
+  Matrix BdNdw(p1*p2*p3, nGauss*nGauss*nGauss);
+  int ig=1; // gauss point iterator
+  for(int zeta=0; zeta<nGauss; zeta++) {
+    for(int eta=0; eta<nGauss; eta++) {
+      for(int xi=0; xi<nGauss; xi++, ig++) {
+        double u[2*p1];
+        double v[2*p2];
+        double w[2*p3];
+        basis1.computeBasisValues(xg[xi],   u, 1); 
+        basis2.computeBasisValues(xg[eta],  v, 1);
+        basis3.computeBasisValues(xg[zeta], w, 1);
+        int ib=1; // basis function iterator
+        double sum = 0;
+        for(int k=0; k<p3; k++) {
+          for(int j=0; j<p2; j++) {
+            for(int i=0; i<p1; i++, ib++) {
+              BN(ib,ig)    = u[2*i  ]*v[2*j  ]*w[2*k  ];
+              BdNdu(ib,ig) = u[2*i+1]*v[2*j  ]*w[2*k  ];
+              BdNdv(ib,ig) = u[2*i  ]*v[2*j+1]*w[2*k  ];
+              BdNdw(ib,ig) = u[2*i  ]*v[2*j  ]*w[2*k+1];
+              sum += BN(ib,ig);
+            }
+          }
+        }
+      }
+    }
+  }
 
-	// Get the reduced integration quadrature points, if needed
-	const double* xr = nullptr;
-	const double* wr = nullptr;
-	int nRed = integrand.getReducedIntegration(nGauss);
-	if (nRed > 0)
-	{
-		xr = GaussQuadrature::getCoord(nRed);
-		wr = GaussQuadrature::getWeight(nRed);
-		if (!xr || !wr) return false;
-	}
-	else if (nRed < 0)
-		nRed = nGauss; // The integrand needs to know nGauss
+  // Get the reduced integration quadrature points, if needed
+  const double* xr = nullptr;
+  const double* wr = nullptr;
+  int nRed = integrand.getReducedIntegration(nGauss);
+  if (nRed > 0)
+  {
+    xr = GaussQuadrature::getCoord(nRed);
+    wr = GaussQuadrature::getWeight(nRed);
+    if (!xr || !wr) return false;
+  }
+  else if (nRed < 0)
+    nRed = nGauss; // The integrand needs to know nGauss
 
 
-	// === Assembly loop over all elements in the patch ==========================
+  // === Assembly loop over all elements in the patch ==========================
 
-	bool ok=true;
-	for(LR::Element *el : lrspline->getAllElements()) {
-		int iEl = el->getId();
-		int nBasis = el->nBasisFunctions();
-		FiniteElement fe(nBasis);
-		fe.iel = iEl+1;
+  bool ok=true;
+  for(LR::Element *el : lrspline->getAllElements()) {
+    int iEl = el->getId();
+    int nBasis = el->nBasisFunctions();
+    FiniteElement fe(nBasis);
+    fe.iel = iEl+1;
 
-		Matrix   C = bezierExtract[iEl];
-		Matrix   dNdu, Xnod, Jac;
-		Matrix3D d2Ndu2, Hess;
-		double   dXidu[3];
-		Vec4     X;
-		// Get element volume in the parameter space
-		double du = el->umax() - el->umin();
-		double dv = el->vmax() - el->vmin();
-		double dw = el->wmax() - el->wmin();
-		double vol = el->volume();
-		if (vol < 0.0)
-		{
-			ok = false; // topology error (probably logic error)
-			break;
-		}
+    Matrix   C = bezierExtract[iEl];
+    Matrix   dNdu, Xnod, Jac;
+    Matrix3D d2Ndu2, Hess;
+    double   dXidu[3];
+    Vec4     X;
+    // Get element volume in the parameter space
+    double du = el->umax() - el->umin();
+    double dv = el->vmax() - el->vmin();
+    double dw = el->wmax() - el->wmin();
+    double vol = el->volume();
+    if (vol < 0.0)
+    {
+      ok = false; // topology error (probably logic error)
+      break;
+    }
 
-		// Set up control point (nodal) coordinates for current element
-		if (!this->getElementCoordinates(Xnod,iEl+1))
-		{
-			ok = false;
-			break;
-		}
+    // Set up control point (nodal) coordinates for current element
+    if (!this->getElementCoordinates(Xnod,iEl+1))
+    {
+      ok = false;
+      break;
+    }
 
-		// Compute parameter values of the Gauss points over the whole element
-		std::array<RealArray,3> gpar, redpar;
-		for (int d = 0; d < 3; d++)
-		{
-			this->getGaussPointParameters(gpar[d],d,nGauss,iEl,xg);
-			if (xr)
-				this->getGaussPointParameters(redpar[d],d,nRed,iEl,xr);
-		}
+    // Compute parameter values of the Gauss points over the whole element
+    std::array<RealArray,3> gpar, redpar;
+    for (int d = 0; d < 3; d++)
+    {
+      this->getGaussPointParameters(gpar[d],d,nGauss,iEl,xg);
+      if (xr)
+        this->getGaussPointParameters(redpar[d],d,nRed,iEl,xr);
+    }
 
 
-		if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
-			this->getElementCorners(iEl, fe.XC);
+    if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
+      this->getElementCorners(iEl, fe.XC);
 
-		if (integrand.getIntegrandType() & Integrand::G_MATRIX)
-		{
-			// Element size in parametric space
-			dXidu[0] = el->getParmin(0);
-			dXidu[1] = el->getParmin(1);
-			dXidu[2] = el->getParmin(2);
-		}
-		else if (integrand.getIntegrandType() & Integrand::AVERAGE)
-		{
-			// --- Compute average value of basis functions over the element -----
+    if (integrand.getIntegrandType() & Integrand::G_MATRIX)
+    {
+      // Element size in parametric space
+      dXidu[0] = el->getParmin(0);
+      dXidu[1] = el->getParmin(1);
+      dXidu[2] = el->getParmin(2);
+    }
+    else if (integrand.getIntegrandType() & Integrand::AVERAGE)
+    {
+      // --- Compute average value of basis functions over the element -----
 
-			fe.Navg.resize(nBasis,true);
-			double vol = 0.0;
-			for (int k = 0; k < nGauss; k++)
-				for (int j = 0; j < nGauss; j++)
-					for (int i = 0; i < nGauss; i++)
-					{
-						// Fetch basis function derivatives at current integration point
-						evaluateBasis(fe, dNdu);
+      fe.Navg.resize(nBasis,true);
+      double vol = 0.0;
+      for (int k = 0; k < nGauss; k++)
+        for (int j = 0; j < nGauss; j++)
+          for (int i = 0; i < nGauss; i++)
+          {
+            // Fetch basis function derivatives at current integration point
+            evaluateBasis(fe, dNdu);
 
-						// Compute Jacobian determinant of coordinate mapping
-						// and multiply by weight of current integration point
-						double detJac = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu,false);
-						double weight = 0.125*vol*wg[i]*wg[j]*wg[k];
+            // Compute Jacobian determinant of coordinate mapping
+            // and multiply by weight of current integration point
+            double detJac = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu,false);
+            double weight = 0.125*vol*wg[i]*wg[j]*wg[k];
 
-						// Numerical quadrature
-						fe.Navg.add(fe.N,detJac*weight);
-						vol += detJac*weight;
-			}
+            // Numerical quadrature
+            fe.Navg.add(fe.N,detJac*weight);
+            vol += detJac*weight;
+      }
 
-			// Divide by element volume
-			fe.Navg /= vol;
-		}
+      // Divide by element volume
+      fe.Navg /= vol;
+    }
 
-		else if (integrand.getIntegrandType() & Integrand::ELEMENT_CENTER)
-		{
-			// Compute the element center
-			Go::Point X0;
-			double u0 = 0.5*(el->getParmin(0) + el->getParmax(0));
-			double v0 = 0.5*(el->getParmin(1) + el->getParmax(1));
-			double w0 = 0.5*(el->getParmin(2) + el->getParmax(2));
-			lrspline->point(X0,u0,v0,w0);
-			X = SplineUtils::toVec3(X0);
-		}
+    else if (integrand.getIntegrandType() & Integrand::ELEMENT_CENTER)
+    {
+      // Compute the element center
+      Go::Point X0;
+      double u0 = 0.5*(el->getParmin(0) + el->getParmax(0));
+      double v0 = 0.5*(el->getParmin(1) + el->getParmax(1));
+      double w0 = 0.5*(el->getParmin(2) + el->getParmax(2));
+      lrspline->point(X0,u0,v0,w0);
+      X = SplineUtils::toVec3(X0);
+    }
 
-		// Initialize element quantities
-		LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel);
-		if (!integrand.initElement(MNPC[iEl],fe,X,nRed*nRed*nRed,*A))
-		{
-			A->destruct();
-			ok = false;
-			break;
-		}
+    // Initialize element quantities
+    LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel);
+    if (!integrand.initElement(MNPC[iEl],fe,X,nRed*nRed*nRed,*A))
+    {
+      A->destruct();
+      ok = false;
+      break;
+    }
 
-		if (xr)
-		{
-			std::cerr << "Haven't really figured out what this part does yet\n";
-			exit(42142);
+    if (xr)
+    {
+      std::cerr << "Haven't really figured out what this part does yet\n";
+      exit(42142);
 #if 0
-			// --- Selective reduced integration loop ----------------------------
+      // --- Selective reduced integration loop ----------------------------
 
-			int ip = (((i3-p3)*nRed*nel2 + i2-p2)*nRed*nel1 + i1-p1)*nRed;
-			for (int k = 0; k < nRed; k++, ip += nRed*(nel2-1)*nRed*nel1)
-				for (int j = 0; j < nRed; j++, ip += nRed*(nel1-1))
-					for (int i = 0; i < nRed; i++, ip++)
-					{
-						// Local element coordinates of current integration point
-						fe.xi   = xr[i];
-						fe.eta  = xr[j];
-						fe.zeta = xr[k];
+      int ip = (((i3-p3)*nRed*nel2 + i2-p2)*nRed*nel1 + i1-p1)*nRed;
+      for (int k = 0; k < nRed; k++, ip += nRed*(nel2-1)*nRed*nel1)
+        for (int j = 0; j < nRed; j++, ip += nRed*(nel1-1))
+          for (int i = 0; i < nRed; i++, ip++)
+          {
+            // Local element coordinates of current integration point
+            fe.xi   = xr[i];
+            fe.eta  = xr[j];
+            fe.zeta = xr[k];
 
-						// Parameter values of current integration point
-						fe.u = redpar[0](i+1,i1-p1+1);
-						fe.v = redpar[1](j+1,i2-p2+1);
-						fe.w = redpar[2](k+1,i3-p3+1);
+            // Parameter values of current integration point
+            fe.u = redpar[0](i+1,i1-p1+1);
+            fe.v = redpar[1](j+1,i2-p2+1);
+            fe.w = redpar[2](k+1,i3-p3+1);
 
-						// Fetch basis function derivatives at current point
-						evaluateBasis(fe, 1); 
+            // Fetch basis function derivatives at current point
+            evaluateBasis(fe, 1); 
 
-						// Compute Jacobian inverse and derivatives
-						fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
+            // Compute Jacobian inverse and derivatives
+            fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
 
-						// Cartesian coordinates of current integration point
-						X = Xnod * fe.N;
-						X.t = time.t;
+            // Cartesian coordinates of current integration point
+            X = Xnod * fe.N;
+            X.t = time.t;
 
-						// Compute the reduced integration terms of the integrand
-						fe.detJxW *= 0.125*dV*wr[i]*wr[j]*wr[k];
-						if (!integrand.reducedInt(*A,fe,X))
-							ok = false;
-			}
+            // Compute the reduced integration terms of the integrand
+            fe.detJxW *= 0.125*dV*wr[i]*wr[j]*wr[k];
+            if (!integrand.reducedInt(*A,fe,X))
+              ok = false;
+      }
 #endif
-		}
+    }
 
 
-		// --- Integration loop over all Gauss points in each direction --------
+    // --- Integration loop over all Gauss points in each direction --------
 
-		fe.iGP = iEl*nGauss*nGauss; // Global integration point counter
+    fe.iGP = iEl*nGauss*nGauss; // Global integration point counter
 
-		Matrix B(p1*p2*p3, 4); // Bezier evaluation points and derivatives
-		ig = 1;
-		for (int k = 0; k < nGauss; k++) 
-			for (int j = 0; j < nGauss; j++) 
-				for (int i = 0; i < nGauss; i++, fe.iGP++, ig++)
-				{
-					// Local element coordinates of current integration point
-					fe.xi   = xg[i];
-					fe.eta  = xg[j];
-					fe.zeta = xg[k];
+    Matrix B(p1*p2*p3, 4); // Bezier evaluation points and derivatives
+    ig = 1;
+    for (int k = 0; k < nGauss; k++) 
+      for (int j = 0; j < nGauss; j++) 
+        for (int i = 0; i < nGauss; i++, fe.iGP++, ig++)
+        {
+          // Local element coordinates of current integration point
+          fe.xi   = xg[i];
+          fe.eta  = xg[j];
+          fe.zeta = xg[k];
 
-					// Parameter values of current integration point
-					fe.u = gpar[0][i];
-					fe.v = gpar[1][j];
-					fe.w = gpar[2][k];
+          // Parameter values of current integration point
+          fe.u = gpar[0][i];
+          fe.v = gpar[1][j];
+          fe.w = gpar[2][k];
 
-					// Extract bezier basis functions
-					B.fillColumn(1, BN.getColumn(ig));
-					B.fillColumn(2, BdNdu.getColumn(ig)*2.0/du);
-					B.fillColumn(3, BdNdv.getColumn(ig)*2.0/dv);
-					B.fillColumn(4, BdNdw.getColumn(ig)*2.0/dw);
+          // Extract bezier basis functions
+          B.fillColumn(1, BN.getColumn(ig));
+          B.fillColumn(2, BdNdu.getColumn(ig)*2.0/du);
+          B.fillColumn(3, BdNdv.getColumn(ig)*2.0/dv);
+          B.fillColumn(4, BdNdw.getColumn(ig)*2.0/dw);
 
-					// Fetch basis function derivatives at current integration point
-					if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
-						evaluateBasis(fe, dNdu, d2Ndu2);
-					else
-						evaluateBasis(fe, dNdu, C, B) ;
+          // Fetch basis function derivatives at current integration point
+          if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
+            evaluateBasis(fe, dNdu, d2Ndu2);
+          else
+            evaluateBasis(fe, dNdu, C, B) ;
 
-					// look for errors in bezier extraction
-					/*
-					int N    = nBasis;
-					int allP = p1*p2*p3;
-					double sum = 0;
-					for(int qq=1; qq<=N; qq++) sum+= fe.N(qq);
-					if(fabs(sum-1) > 1e-10) {
-						std::cerr << "fe.N not sums to one at integration point #" << ig << std::endl;
-						exit(123);
-					}
-					sum = 0;
-					for(int qq=1; qq<=N; qq++) sum+= dNdu(qq,1);
-					if(fabs(sum) > 1e-10) {
-						std::cerr << "dNdu not sums to zero at integration point #" << ig << std::endl;
-						exit(123);
-					}
-					sum = 0;
-					for(int qq=1; qq<=N; qq++) sum+= dNdu(qq,2);
-					if(fabs(sum) > 1e-10) {
-						std::cerr << "dNdv not sums to zero at integration point #" << ig << std::endl;
-						exit(123);
-					}
-					sum = 0;
-					for(int qq=1; qq<=N; qq++) sum+= dNdu(qq,3);
-					if(fabs(sum) > 1e-10) {
-						std::cerr << "dNdw not sums to zero at integration point #" << ig << std::endl;
-						exit(123);
-					}
-					sum = 0;
-					for(int qq=1; qq<=allP; qq++) sum+= B(qq,1);
-					if(fabs(sum-1) > 1e-10) {
-						std::cerr << "Bezier basis not sums to one at integration point #" << ig << std::endl;
-						exit(123);
-					}
-					sum = 0;
-					for(int qq=1; qq<=allP; qq++) sum+= B(qq,2);
-					if(fabs(sum) > 1e-10) {
-						std::cerr << "Bezier derivatives not sums to zero at integration point #" << ig << std::endl;
-						exit(123);
-					}
-					*/
+          // look for errors in bezier extraction
+          /*
+          int N    = nBasis;
+          int allP = p1*p2*p3;
+          double sum = 0;
+          for(int qq=1; qq<=N; qq++) sum+= fe.N(qq);
+          if(fabs(sum-1) > 1e-10) {
+            std::cerr << "fe.N not sums to one at integration point #" << ig << std::endl;
+            exit(123);
+          }
+          sum = 0;
+          for(int qq=1; qq<=N; qq++) sum+= dNdu(qq,1);
+          if(fabs(sum) > 1e-10) {
+            std::cerr << "dNdu not sums to zero at integration point #" << ig << std::endl;
+            exit(123);
+          }
+          sum = 0;
+          for(int qq=1; qq<=N; qq++) sum+= dNdu(qq,2);
+          if(fabs(sum) > 1e-10) {
+            std::cerr << "dNdv not sums to zero at integration point #" << ig << std::endl;
+            exit(123);
+          }
+          sum = 0;
+          for(int qq=1; qq<=N; qq++) sum+= dNdu(qq,3);
+          if(fabs(sum) > 1e-10) {
+            std::cerr << "dNdw not sums to zero at integration point #" << ig << std::endl;
+            exit(123);
+          }
+          sum = 0;
+          for(int qq=1; qq<=allP; qq++) sum+= B(qq,1);
+          if(fabs(sum-1) > 1e-10) {
+            std::cerr << "Bezier basis not sums to one at integration point #" << ig << std::endl;
+            exit(123);
+          }
+          sum = 0;
+          for(int qq=1; qq<=allP; qq++) sum+= B(qq,2);
+          if(fabs(sum) > 1e-10) {
+            std::cerr << "Bezier derivatives not sums to zero at integration point #" << ig << std::endl;
+            exit(123);
+          }
+          */
 
-					// Compute Jacobian inverse of coordinate mapping and derivatives
-					fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
-					if (fe.detJxW == 0.0) continue; // skip singular points
+          // Compute Jacobian inverse of coordinate mapping and derivatives
+          fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
+          if (fe.detJxW == 0.0) continue; // skip singular points
 
-					// Compute Hessian of coordinate mapping and 2nd order derivatives
-					if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
-						if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
-							ok = false;
+          // Compute Hessian of coordinate mapping and 2nd order derivatives
+          if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
+            if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
+              ok = false;
 
-					// Compute G-matrix
-					if (integrand.getIntegrandType() & Integrand::G_MATRIX)
-						utl::getGmat(Jac,dXidu,fe.G);
+          // Compute G-matrix
+          if (integrand.getIntegrandType() & Integrand::G_MATRIX)
+            utl::getGmat(Jac,dXidu,fe.G);
 
-					// Cartesian coordinates of current integration point
-					X   = Xnod * fe.N;
-					X.t = time.t;
+          // Cartesian coordinates of current integration point
+          X   = Xnod * fe.N;
+          X.t = time.t;
 
-					// Evaluate the integrand and accumulate element contributions
-					fe.detJxW *= 0.125*vol*wg[i]*wg[j]*wg[k];
-					if (!integrand.evalInt(*A,fe,time,X))
-						ok = false;
+          // Evaluate the integrand and accumulate element contributions
+          fe.detJxW *= 0.125*vol*wg[i]*wg[j]*wg[k];
+          if (!integrand.evalInt(*A,fe,time,X))
+            ok = false;
 
-		} // end gauss integrand
+    } // end gauss integrand
 
-		// Finalize the element quantities
-		if (ok && !integrand.finalizeElement(*A,time,0))
-			ok = false;
+    // Finalize the element quantities
+    if (ok && !integrand.finalizeElement(*A,time,0))
+      ok = false;
 
-		// Assembly of global system integral
-		if (ok && !glInt.assemble(A->ref(),fe.iel))
-			ok = false;
+    // Assembly of global system integral
+    if (ok && !glInt.assemble(A->ref(),fe.iel))
+      ok = false;
 
-		A->destruct();
-	}
+    A->destruct();
+  }
 
-	return ok;
+  return ok;
 }
 
 
@@ -1226,178 +1226,178 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex,
                         GlobalIntegral& glInt,
                         const TimeDomain& time)
 {
-	if (!lrspline) return true; // silently ignore empty patches
+  if (!lrspline) return true; // silently ignore empty patches
 
-	PROFILE2("ASMu3D::integrate(B)");
+  PROFILE2("ASMu3D::integrate(B)");
 
-	// Get Gaussian quadrature points and weights
-	int nGP = integrand.getBouIntegrationPoints(nGauss);
-	const double* xg = GaussQuadrature::getCoord(nGP);
-	const double* wg = GaussQuadrature::getWeight(nGP);
-	if (!xg || !wg) return false;
+  // Get Gaussian quadrature points and weights
+  int nGP = integrand.getBouIntegrationPoints(nGauss);
+  const double* xg = GaussQuadrature::getCoord(nGP);
+  const double* wg = GaussQuadrature::getWeight(nGP);
+  if (!xg || !wg) return false;
 
-	// Find the parametric direction of the face normal {-3,-2,-1, 1, 2, 3}
-	const int faceDir = (lIndex+1)/(lIndex%2 ? -2 : 2);
+  // Find the parametric direction of the face normal {-3,-2,-1, 1, 2, 3}
+  const int faceDir = (lIndex+1)/(lIndex%2 ? -2 : 2);
 
-	const int t1 = 1 + abs(faceDir)%3; // first tangent direction
-	const int t2 = 1 + t1%3;           // second tangent direction
+  const int t1 = 1 + abs(faceDir)%3; // first tangent direction
+  const int t2 = 1 + t1%3;           // second tangent direction
 
   std::map<char,size_t>::const_iterator iit = firstBp.find(lIndex);
   size_t firstp = iit == firstBp.end() ? 0 : iit->second;
-	
+  
   LR::parameterEdge edge;
-	switch(lIndex)
-	{
-	case 1: edge = LR::WEST;   break;
-	case 2: edge = LR::EAST;   break;
-	case 3: edge = LR::SOUTH;  break;
-	case 4: edge = LR::NORTH;  break;
-	case 5: edge = LR::BOTTOM; break;
-	case 6: edge = LR::TOP;    break;
-	default:edge = LR::NONE;
-	}
-	
-	// fetch all elements along the chosen edge
-	std::vector<LR::Element*> edgeElms;
-	lrspline->getEdgeElements(edgeElms, (LR::parameterEdge) edge);
+  switch(lIndex)
+  {
+  case 1: edge = LR::WEST;   break;
+  case 2: edge = LR::EAST;   break;
+  case 3: edge = LR::SOUTH;  break;
+  case 4: edge = LR::NORTH;  break;
+  case 5: edge = LR::BOTTOM; break;
+  case 6: edge = LR::TOP;    break;
+  default:edge = LR::NONE;
+  }
+  
+  // fetch all elements along the chosen edge
+  std::vector<LR::Element*> edgeElms;
+  lrspline->getEdgeElements(edgeElms, (LR::parameterEdge) edge);
 
-	// iterate over all edge elements
-	bool ok = true;
-	for(LR::Element *el : edgeElms) {
-		int iEl = el->getId();
-		int nBasis = el->nBasisFunctions();
-		FiniteElement fe(nBasis);
-		fe.iel = iEl+1;
+  // iterate over all edge elements
+  bool ok = true;
+  for(LR::Element *el : edgeElms) {
+    int iEl = el->getId();
+    int nBasis = el->nBasisFunctions();
+    FiniteElement fe(nBasis);
+    fe.iel = iEl+1;
 
-		// Compute parameter values of the Gauss points over the whole element
-		std::array<Vector,3> gpar;
-		for (int d = 0; d < 3; d++)
-			if (-1-d == faceDir)
-			{
-				gpar[d].resize(1);
-				gpar[d].fill(lrspline->startparam(d));
-			}
-			else if (1+d == faceDir)
-			{
-				gpar[d].resize(1);
-				gpar[d].fill(lrspline->endparam(d));
-			}
-			else
-				this->getGaussPointParameters(gpar[d],d,nGP,iEl,xg);
+    // Compute parameter values of the Gauss points over the whole element
+    std::array<Vector,3> gpar;
+    for (int d = 0; d < 3; d++)
+      if (-1-d == faceDir)
+      {
+        gpar[d].resize(1);
+        gpar[d].fill(lrspline->startparam(d));
+      }
+      else if (1+d == faceDir)
+      {
+        gpar[d].resize(1);
+        gpar[d].fill(lrspline->endparam(d));
+      }
+      else
+        this->getGaussPointParameters(gpar[d],d,nGP,iEl,xg);
 
-		fe.xi = fe.eta = fe.zeta = faceDir < 0 ? -1.0 : 1.0;
-		fe.u = gpar[0](1);
-		fe.v = gpar[1](1);
-		fe.w = gpar[2](1);
+    fe.xi = fe.eta = fe.zeta = faceDir < 0 ? -1.0 : 1.0;
+    fe.u = gpar[0](1);
+    fe.v = gpar[1](1);
+    fe.w = gpar[2](1);
 
-		Matrix dNdu, Xnod, Jac;
-		Vec4   X;
-		Vec3   normal;
-		double dXidu[3];
+    Matrix dNdu, Xnod, Jac;
+    Vec4   X;
+    Vec3   normal;
+    double dXidu[3];
 
-		// Get element face area in the parameter space
-		double dA = this->getParametricArea(iEl+1,abs(faceDir));
-		if (dA < 0.0) // topology error (probably logic error)
-		{
-		  ok = false;
-		  break;
-		}
+    // Get element face area in the parameter space
+    double dA = this->getParametricArea(iEl+1,abs(faceDir));
+    if (dA < 0.0) // topology error (probably logic error)
+    {
+      ok = false;
+      break;
+    }
 
-		// Set up control point coordinates for current element
-		if (!this->getElementCoordinates(Xnod,iEl+1))
-		{
-		  ok = false;
-		  break;
-		}
+    // Set up control point coordinates for current element
+    if (!this->getElementCoordinates(Xnod,iEl+1))
+    {
+      ok = false;
+      break;
+    }
 
-		if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
-		  this->getElementCorners(iEl,fe.XC);
+    if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
+      this->getElementCorners(iEl,fe.XC);
 
-		if (integrand.getIntegrandType() & Integrand::G_MATRIX)
-		{
-		  // Element size in parametric space
-		  dXidu[0] = el->getParmax(0) - el->getParmin(0);   
-		  dXidu[1] = el->getParmax(1) - el->getParmin(1); 
-		  dXidu[2] = el->getParmax(2) - el->getParmin(2);
-		}
+    if (integrand.getIntegrandType() & Integrand::G_MATRIX)
+    {
+      // Element size in parametric space
+      dXidu[0] = el->getParmax(0) - el->getParmin(0);   
+      dXidu[1] = el->getParmax(1) - el->getParmin(1); 
+      dXidu[2] = el->getParmax(2) - el->getParmin(2);
+    }
 
-        // Initialize element quantities
-        LocalIntegral* A = integrand.getLocalIntegral(nBasis,fe.iel,true);
-        if (!integrand.initElementBou(MNPC[iEl],*A))
+    // Initialize element quantities
+    LocalIntegral* A = integrand.getLocalIntegral(nBasis,fe.iel,true);
+    if (!integrand.initElementBou(MNPC[iEl],*A))
+    {
+      A->destruct();
+      ok = false;
+      break;
+    }
+
+    // --- Integration loop over all Gauss points in each direction --------
+
+    fe.iGP = firstp; // Global integration point counter
+    int k1,k2,k3;
+    for (int j = 0; j < nGP; j++)
+      for (int i = 0; i < nGP; i++, fe.iGP++)
+      {
+        // Local element coordinates and parameter values
+        // of current integration point
+        switch (abs(faceDir))
         {
-          A->destruct();
-          ok = false;
-          break;
+          case 1: k2 = i; k3 = j; k1 = 0; break;
+          case 2: k1 = i; k3 = j; k2 = 0; break;
+          case 3: k1 = i; k2 = j; k3 = 0; break;
+          default: k1 = k2 = k3 = 0;
+        }
+        if (gpar[0].size() > 1)
+        {
+          fe.xi = xg[k1];
+          fe.u = gpar[0](k1+1);
+        }
+        if (gpar[1].size() > 1)
+        {
+          fe.eta = xg[k2];
+          fe.v = gpar[1](k2+1);
+        }
+        if (gpar[2].size() > 1)
+        {
+          fe.zeta = xg[k3];
+          fe.w = gpar[2](k3+1);
         }
 
-		// --- Integration loop over all Gauss points in each direction --------
+        // Fetch basis function derivatives at current integration point
+        evaluateBasis(fe, dNdu);
 
-		fe.iGP = firstp; // Global integration point counter
-		int k1,k2,k3;
-		for (int j = 0; j < nGP; j++)
-			for (int i = 0; i < nGP; i++, fe.iGP++)
-			{
-				// Local element coordinates and parameter values
-				// of current integration point
-				switch (abs(faceDir))
-				{
-					case 1: k2 = i; k3 = j; k1 = 0; break;
-					case 2: k1 = i; k3 = j; k2 = 0; break;
-					case 3: k1 = i; k2 = j; k3 = 0; break;
-					default: k1 = k2 = k3 = 0;
-				}
-				if (gpar[0].size() > 1)
-				{
-					fe.xi = xg[k1];
-					fe.u = gpar[0](k1+1);
-				}
-				if (gpar[1].size() > 1)
-				{
-					fe.eta = xg[k2];
-					fe.v = gpar[1](k2+1);
-				}
-				if (gpar[2].size() > 1)
-				{
-					fe.zeta = xg[k3];
-					fe.w = gpar[2](k3+1);
-				}
+        // Compute basis function derivatives and the face normal
+        fe.detJxW = utl::Jacobian(Jac,normal,fe.dNdX,Xnod,dNdu,t1,t2);
+        if (fe.detJxW == 0.0) continue; // skip singular points
 
-				// Fetch basis function derivatives at current integration point
-				evaluateBasis(fe, dNdu);
+        if (faceDir < 0) normal *= -1.0;
 
-				// Compute basis function derivatives and the face normal
-				fe.detJxW = utl::Jacobian(Jac,normal,fe.dNdX,Xnod,dNdu,t1,t2);
-				if (fe.detJxW == 0.0) continue; // skip singular points
+        // Compute G-matrix
+        if (integrand.getIntegrandType() & Integrand::G_MATRIX)
+          utl::getGmat(Jac,dXidu,fe.G);
 
-				if (faceDir < 0) normal *= -1.0;
+        // Cartesian coordinates of current integration point
+        X = Xnod * fe.N;
+        X.t = time.t;
 
-				// Compute G-matrix
-				if (integrand.getIntegrandType() & Integrand::G_MATRIX)
-					utl::getGmat(Jac,dXidu,fe.G);
-
-				// Cartesian coordinates of current integration point
-				X = Xnod * fe.N;
-				X.t = time.t;
-
-				// Evaluate the integrand and accumulate element contributions
-				fe.detJxW *= 0.25*dA*wg[i]*wg[j];
-				if (!integrand.evalBou(*A,fe,time,X,normal))
-					ok = false;
-		}
+        // Evaluate the integrand and accumulate element contributions
+        fe.detJxW *= 0.25*dA*wg[i]*wg[j];
+        if (!integrand.evalBou(*A,fe,time,X,normal))
+          ok = false;
+    }
 
     // Finalize the element quantities
     if (ok && !integrand.finalizeElementBou(*A,fe,time))
       ok = false;
 
-		// Assembly of global system integral
-		if (ok && !glInt.assemble(A->ref(),fe.iel))
-			ok = false;
-		A->destruct();
+    // Assembly of global system integral
+    if (ok && !glInt.assemble(A->ref(),fe.iel))
+      ok = false;
+    A->destruct();
 
-		firstp += nGP*nGP;
-	}
+    firstp += nGP*nGP;
+  }
 
-	return ok;
+  return ok;
 }
 
 
@@ -1405,326 +1405,326 @@ bool ASMu3D::integrateEdge (Integrand& integrand, int lEdge,
                             GlobalIntegral& glInt,
                             const TimeDomain& time)
 {
-	std::cerr << "ASMu3D::integrateEdge(...) is not properly implemented yet :(" << std::endl;
-	exit(776654);
+  std::cerr << "ASMu3D::integrateEdge(...) is not properly implemented yet :(" << std::endl;
+  exit(776654);
 #if 0
-	if (!lrspline) return true; // silently ignore empty patches
+  if (!lrspline) return true; // silently ignore empty patches
 
-	PROFILE2("ASMu3D::integrate(E)");
+  PROFILE2("ASMu3D::integrate(E)");
 
-	// Get Gaussian quadrature points and weights
-	const double* xg = GaussQuadrature::getCoord(nGauss);
-	const double* wg = GaussQuadrature::getWeight(nGauss);
-	if (!xg || !wg) return false;
+  // Get Gaussian quadrature points and weights
+  const double* xg = GaussQuadrature::getCoord(nGauss);
+  const double* wg = GaussQuadrature::getWeight(nGauss);
+  if (!xg || !wg) return false;
 
-	// Compute parameter values of the Gauss points along the whole patch edge
-	std::array<Matrix,3> gpar;
-	for (int d = 0; d < 3; d++)
-		if (lEdge < d*4+1 || lEdge >= d*4+5)
-		{
-			gpar[d].resize(1,1);
-			if (lEdge%4 == 1)
-	gpar[d].fill(lrspline->startparam(d));
-			else if (lEdge%4 == 0)
-	gpar[d].fill(lrspline->endparam(d));
-			else if (lEdge == 6 || lEdge == 10)
-	gpar[d].fill(d == 0 ? lrspline->endparam(d) : lrspline->startparam(d));
-			else if (lEdge == 2 || lEdge == 11)
-	gpar[d].fill(d == 1 ? lrspline->endparam(d) : lrspline->startparam(d));
-			else if (lEdge == 3 || lEdge == 7)
-	gpar[d].fill(d == 2 ? lrspline->endparam(d) : lrspline->startparam(d));
-		}
-		else
-		{
-			int pm1 = lrspline->order(d) - 1;
-			RealArray::const_iterator uit = lrspline->basis(d).begin() + pm1;
-			double ucurr, uprev = *(uit++);
-			int nCol = lrspline->numCoefs(d) - pm1;
-			gpar[d].resize(nGauss,nCol);
-			for (int j = 1; j <= nCol; uit++, j++)
-			{
-	ucurr = *uit;
-	for (int i = 1; i <= nGauss; i++)
-	  gpar[d](i,j) = 0.5*((ucurr-uprev)*xg[i-1] + ucurr+uprev);
-	uprev = ucurr;
-			}
-		}
+  // Compute parameter values of the Gauss points along the whole patch edge
+  std::array<Matrix,3> gpar;
+  for (int d = 0; d < 3; d++)
+    if (lEdge < d*4+1 || lEdge >= d*4+5)
+    {
+      gpar[d].resize(1,1);
+      if (lEdge%4 == 1)
+  gpar[d].fill(lrspline->startparam(d));
+      else if (lEdge%4 == 0)
+  gpar[d].fill(lrspline->endparam(d));
+      else if (lEdge == 6 || lEdge == 10)
+  gpar[d].fill(d == 0 ? lrspline->endparam(d) : lrspline->startparam(d));
+      else if (lEdge == 2 || lEdge == 11)
+  gpar[d].fill(d == 1 ? lrspline->endparam(d) : lrspline->startparam(d));
+      else if (lEdge == 3 || lEdge == 7)
+  gpar[d].fill(d == 2 ? lrspline->endparam(d) : lrspline->startparam(d));
+    }
+    else
+    {
+      int pm1 = lrspline->order(d) - 1;
+      RealArray::const_iterator uit = lrspline->basis(d).begin() + pm1;
+      double ucurr, uprev = *(uit++);
+      int nCol = lrspline->numCoefs(d) - pm1;
+      gpar[d].resize(nGauss,nCol);
+      for (int j = 1; j <= nCol; uit++, j++)
+      {
+  ucurr = *uit;
+  for (int i = 1; i <= nGauss; i++)
+    gpar[d](i,j) = 0.5*((ucurr-uprev)*xg[i-1] + ucurr+uprev);
+  uprev = ucurr;
+      }
+    }
 
-	// Evaluate basis function derivatives at all integration points
-	std::vector<Go::BasisDerivs> spline;
-	{
-		PROFILE2("Spline evaluation");
-		lrspline->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline);
-	}
+  // Evaluate basis function derivatives at all integration points
+  std::vector<Go::BasisDerivs> spline;
+  {
+    PROFILE2("Spline evaluation");
+    lrspline->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline);
+  }
 
-	const int n1 = lrspline->numCoefs(0);
-	const int n2 = lrspline->numCoefs(1);
-	const int n3 = lrspline->numCoefs(2);
+  const int n1 = lrspline->numCoefs(0);
+  const int n2 = lrspline->numCoefs(1);
+  const int n3 = lrspline->numCoefs(2);
 
-	const int p1 = lrspline->order(0);
-	const int p2 = lrspline->order(1);
-	const int p3 = lrspline->order(2);
+  const int p1 = lrspline->order(0);
+  const int p2 = lrspline->order(1);
+  const int p3 = lrspline->order(2);
 
-	std::map<char,size_t>::const_iterator iit = firstBp.find(lEdge);
-	size_t firstp = iit == firstBp.end() ? 0 : iit->second;
+  std::map<char,size_t>::const_iterator iit = firstBp.find(lEdge);
+  size_t firstp = iit == firstBp.end() ? 0 : iit->second;
 
-	FiniteElement fe(p1*p2*p3);
-	fe.u = gpar[0](1,1);
-	fe.v = gpar[1](1,1);
-	fe.w = gpar[2](1,1);
-	if (gpar[0].size() == 1) fe.xi = fe.u == lrspline->startparam(0) ? -1.0 : 1.0;
-	if (gpar[1].size() == 1) fe.eta = fe.v == lrspline->startparam(1) ? -1.0 : 1.0;
-	if (gpar[2].size() == 1) fe.zeta = fe.w == lrspline->startparam(2) ? -1.0 : 1.0;
+  FiniteElement fe(p1*p2*p3);
+  fe.u = gpar[0](1,1);
+  fe.v = gpar[1](1,1);
+  fe.w = gpar[2](1,1);
+  if (gpar[0].size() == 1) fe.xi = fe.u == lrspline->startparam(0) ? -1.0 : 1.0;
+  if (gpar[1].size() == 1) fe.eta = fe.v == lrspline->startparam(1) ? -1.0 : 1.0;
+  if (gpar[2].size() == 1) fe.zeta = fe.w == lrspline->startparam(2) ? -1.0 : 1.0;
 
-	Matrix dNdu, Xnod, Jac;
-	Vec4   X;
-	Vec3   tang;
+  Matrix dNdu, Xnod, Jac;
+  Vec4   X;
+  Vec3   tang;
 
 
-	// === Assembly loop over all elements on the patch edge =====================
+  // === Assembly loop over all elements on the patch edge =====================
 
-	int iel = 1;
-	for (int i3 = p3; i3 <= n3; i3++)
-		for (int i2 = p2; i2 <= n2; i2++)
-			for (int i1 = p1; i1 <= n1; i1++, iel++)
-			{
-	fe.iel = MLGE[iel-1];
-	if (fe.iel < 1) continue; // zero-volume element
+  int iel = 1;
+  for (int i3 = p3; i3 <= n3; i3++)
+    for (int i2 = p2; i2 <= n2; i2++)
+      for (int i1 = p1; i1 <= n1; i1++, iel++)
+      {
+  fe.iel = MLGE[iel-1];
+  if (fe.iel < 1) continue; // zero-volume element
 
-	// Skip elements that are not on current boundary edge
-	bool skipMe = false;
-	switch (lEdge)
-	  {
-	  case  1: if (i2 > p2 || i3 > p3) skipMe = true; break;
-	  case  2: if (i2 < n2 || i3 > p3) skipMe = true; break;
-	  case  3: if (i2 > p2 || i3 < n3) skipMe = true; break;
-	  case  4: if (i2 < n2 || i3 < n3) skipMe = true; break;
-	  case  5: if (i1 > p1 || i3 > p3) skipMe = true; break;
-	  case  6: if (i1 < n1 || i3 > p3) skipMe = true; break;
-	  case  7: if (i1 > p1 || i3 < n3) skipMe = true; break;
-	  case  8: if (i1 < n1 || i3 < n3) skipMe = true; break;
-	  case  9: if (i1 > p1 || i2 > p2) skipMe = true; break;
-	  case 10: if (i1 < n1 || i2 > p2) skipMe = true; break;
-	  case 11: if (i1 > p1 || i2 < n2) skipMe = true; break;
-	  case 12: if (i1 < n1 || i2 < n2) skipMe = true; break;
-	  }
-	if (skipMe) continue;
+  // Skip elements that are not on current boundary edge
+  bool skipMe = false;
+  switch (lEdge)
+    {
+    case  1: if (i2 > p2 || i3 > p3) skipMe = true; break;
+    case  2: if (i2 < n2 || i3 > p3) skipMe = true; break;
+    case  3: if (i2 > p2 || i3 < n3) skipMe = true; break;
+    case  4: if (i2 < n2 || i3 < n3) skipMe = true; break;
+    case  5: if (i1 > p1 || i3 > p3) skipMe = true; break;
+    case  6: if (i1 < n1 || i3 > p3) skipMe = true; break;
+    case  7: if (i1 > p1 || i3 < n3) skipMe = true; break;
+    case  8: if (i1 < n1 || i3 < n3) skipMe = true; break;
+    case  9: if (i1 > p1 || i2 > p2) skipMe = true; break;
+    case 10: if (i1 < n1 || i2 > p2) skipMe = true; break;
+    case 11: if (i1 > p1 || i2 < n2) skipMe = true; break;
+    case 12: if (i1 < n1 || i2 < n2) skipMe = true; break;
+    }
+  if (skipMe) continue;
 
-	// Get element edge length in the parameter space
-	double dS = 0.0;
-	int ip = MNPC[iel-1].back();
+  // Get element edge length in the parameter space
+  double dS = 0.0;
+  int ip = MNPC[iel-1].back();
 #ifdef INDEX_CHECK
-	if (ip < 0 || (size_t)ip >= nodeInd.size()) return false;
+  if (ip < 0 || (size_t)ip >= nodeInd.size()) return false;
 #endif
-	if (lEdge < 5)
-	{
-	  dS = lrspline->knotSpan(0,nodeInd[ip].I);
-	  ip = (i1-p1)*nGauss;
-	}
-	else if (lEdge < 9)
-	{
-	  dS = lrspline->knotSpan(1,nodeInd[ip].J);
-	  ip = (i2-p2)*nGauss;
-	}
-	else if (lEdge < 13)
-	{
-	  dS = lrspline->knotSpan(2,nodeInd[ip].K);
-	  ip = (i3-p3)*nGauss;
-	}
+  if (lEdge < 5)
+  {
+    dS = lrspline->knotSpan(0,nodeInd[ip].I);
+    ip = (i1-p1)*nGauss;
+  }
+  else if (lEdge < 9)
+  {
+    dS = lrspline->knotSpan(1,nodeInd[ip].J);
+    ip = (i2-p2)*nGauss;
+  }
+  else if (lEdge < 13)
+  {
+    dS = lrspline->knotSpan(2,nodeInd[ip].K);
+    ip = (i3-p3)*nGauss;
+  }
 
-	// Set up control point coordinates for current element
-	if (!this->getElementCoordinates(Xnod,iel)) return false;
+  // Set up control point coordinates for current element
+  if (!this->getElementCoordinates(Xnod,iel)) return false;
 
-	// Initialize element quantities
-	LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel,true);
-	bool ok = integrand.initElementBou(MNPC[iel-1],*A);
+  // Initialize element quantities
+  LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel,true);
+  bool ok = integrand.initElementBou(MNPC[iel-1],*A);
 
 
-	// --- Integration loop over all Gauss points along the edge -----------
+  // --- Integration loop over all Gauss points along the edge -----------
 
-	fe.iGP = firstp + ip; // Global integration point counter
+  fe.iGP = firstp + ip; // Global integration point counter
 
-	for (int i = 0; i < nGauss && ok; i++, ip++, fe.iGP++)
-	{
-	  // Parameter values of current integration point
-	  if (gpar[0].size() > 1) fe.u = gpar[0](i+1,i1-p1+1);
-	  if (gpar[1].size() > 1) fe.v = gpar[1](i+1,i2-p2+1);
-	  if (gpar[2].size() > 1) fe.w = gpar[2](i+1,i3-p3+1);
+  for (int i = 0; i < nGauss && ok; i++, ip++, fe.iGP++)
+  {
+    // Parameter values of current integration point
+    if (gpar[0].size() > 1) fe.u = gpar[0](i+1,i1-p1+1);
+    if (gpar[1].size() > 1) fe.v = gpar[1](i+1,i2-p2+1);
+    if (gpar[2].size() > 1) fe.w = gpar[2](i+1,i3-p3+1);
 
-	  // Fetch basis function derivatives at current integration point
-	  evaluateBasis(fe, dNdu);
+    // Fetch basis function derivatives at current integration point
+    evaluateBasis(fe, dNdu);
 
-	  // Compute basis function derivatives and the edge tang
-	  fe.detJxW = utl::Jacobian(Jac,tang,fe.dNdX,Xnod,dNdu,1+(lEdge-1)/4);
-	  if (fe.detJxW == 0.0) continue; // skip singular points
+    // Compute basis function derivatives and the edge tang
+    fe.detJxW = utl::Jacobian(Jac,tang,fe.dNdX,Xnod,dNdu,1+(lEdge-1)/4);
+    if (fe.detJxW == 0.0) continue; // skip singular points
 
-	  // Cartesian coordinates of current integration point
-	  X = Xnod * fe.N;
-	  X.t = time.t;
+    // Cartesian coordinates of current integration point
+    X = Xnod * fe.N;
+    X.t = time.t;
 
-	  // Evaluate the integrand and accumulate element contributions
-	  fe.detJxW *= 0.5*dS*wg[i];
-				  ok = integrand.evalBou(*A,fe,time,X,tang);
-	}
+    // Evaluate the integrand and accumulate element contributions
+    fe.detJxW *= 0.5*dS*wg[i];
+          ok = integrand.evalBou(*A,fe,time,X,tang);
+  }
 
-	// Assembly of global system integral
-	if (ok && !glInt.assemble(A->ref(),fe.iel))
-	  ok = false;
+  // Assembly of global system integral
+  if (ok && !glInt.assemble(A->ref(),fe.iel))
+    ok = false;
 
-	A->destruct();
+  A->destruct();
 
-	if (!ok) return false;
-	}
+  if (!ok) return false;
+  }
 
-	return true;
+  return true;
 #endif
-	return false;
+  return false;
 }
 
 
 int ASMu3D::evalPoint (const double* xi, double* param, Vec3& X) const
 {
-	std::cerr << "ASMu3D::evalPoint(...) is not properly implemented yet :(" << std::endl;
-	exit(776654);
+  std::cerr << "ASMu3D::evalPoint(...) is not properly implemented yet :(" << std::endl;
+  exit(776654);
 #if 0
-	if (!lrspline) return -3;
+  if (!lrspline) return -3;
 
-	int i;
-	for (i = 0; i < 3; i++)
-		param[i] = (1.0-xi[i])*lrspline->startparam(i) + xi[i]*lrspline->endparam(i);
+  int i;
+  for (i = 0; i < 3; i++)
+    param[i] = (1.0-xi[i])*lrspline->startparam(i) + xi[i]*lrspline->endparam(i);
 
-	Go::Point X0;
-	lrspline->point(X0,param[0],param[1],param[2]);
-	for (i = 0; i < 3 && i < lrspline->dimension(); i++)
-		X[i] = X0[i];
+  Go::Point X0;
+  lrspline->point(X0,param[0],param[1],param[2]);
+  for (i = 0; i < 3 && i < lrspline->dimension(); i++)
+    X[i] = X0[i];
 
-	// Check if this point matches any of the control points (nodes)
-	Vec3 Xnod;
-	size_t inod = 1;
-	RealArray::const_iterator cit = lrspline->coefs_begin();
-	for (i = 0; cit != lrspline->coefs_end(); cit++, i++)
-	{
-		if (i < 3) Xnod[i] = *cit;
-		if (i+1 == lrspline->dimension())
-			if (X.equal(Xnod,0.001))
-	return inod;
-			else
-			{
-	inod++;
-	i = -1;
-			}
-	}
+  // Check if this point matches any of the control points (nodes)
+  Vec3 Xnod;
+  size_t inod = 1;
+  RealArray::const_iterator cit = lrspline->coefs_begin();
+  for (i = 0; cit != lrspline->coefs_end(); cit++, i++)
+  {
+    if (i < 3) Xnod[i] = *cit;
+    if (i+1 == lrspline->dimension())
+      if (X.equal(Xnod,0.001))
+  return inod;
+      else
+      {
+  inod++;
+  i = -1;
+      }
+  }
 
-	return 0;
+  return 0;
 #endif
-	return -1;
+  return -1;
 }
 
 
 bool ASMu3D::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) const
 {
-	if (!lrspline) return false;
+  if (!lrspline) return false;
 
-	// output is written once for each element resulting in a lot of unnecessary storage
-	// this is preferable to figuring out all element topology information
+  // output is written once for each element resulting in a lot of unnecessary storage
+  // this is preferable to figuring out all element topology information
 
-	for(LR::Element *el : lrspline->getAllElements() ) {
-		// evaluate element at element corner points
-		double umin = el->umin();
-		double umax = el->umax();
-		double vmin = el->vmin();
-		double vmax = el->vmax();
-		double wmin = el->wmin();
-		double wmax = el->wmax();
-		for(int iw=0; iw<=nSegPerSpan; iw++) {
-			for(int iv=0; iv<=nSegPerSpan; iv++) {
-				for(int iu=0; iu<=nSegPerSpan; iu++) {
-					double u = umin + (umax-umin)/nSegPerSpan*iu;
-					double v = vmin + (vmax-vmin)/nSegPerSpan*iv;
-					double w = wmin + (wmax-wmin)/nSegPerSpan*iw;
-					if(dir==0)
-						prm.push_back(u);
-					else if(dir==1)
-						prm.push_back(v);
-					else
-						prm.push_back(w);
-				}
-			}
-		}
-	}
-	return true;
+  for(LR::Element *el : lrspline->getAllElements() ) {
+    // evaluate element at element corner points
+    double umin = el->umin();
+    double umax = el->umax();
+    double vmin = el->vmin();
+    double vmax = el->vmax();
+    double wmin = el->wmin();
+    double wmax = el->wmax();
+    for(int iw=0; iw<=nSegPerSpan; iw++) {
+      for(int iv=0; iv<=nSegPerSpan; iv++) {
+        for(int iu=0; iu<=nSegPerSpan; iu++) {
+          double u = umin + (umax-umin)/nSegPerSpan*iu;
+          double v = vmin + (vmax-vmin)/nSegPerSpan*iv;
+          double w = wmin + (wmax-wmin)/nSegPerSpan*iw;
+          if(dir==0)
+            prm.push_back(u);
+          else if(dir==1)
+            prm.push_back(v);
+          else
+            prm.push_back(w);
+        }
+      }
+    }
+  }
+  return true;
 }
 
 bool ASMu3D::tesselate (ElementBlock& grid, const int* npe) const
 {
-	if(!lrspline) return false;
+  if(!lrspline) return false;
 
-	if(npe[0] != npe[1] || npe[0] != npe[2]) {
-		std::cerr << "ASMu2D::tesselate does not support different tesselation resolution in "
-		          << "u- and v-direction. nviz u = " << npe[0] << ", nviz v = " << npe[1] 
-		          << ", nviz w = " << npe[2] << std::endl;
-		return false;
-	}
+  if(npe[0] != npe[1] || npe[0] != npe[2]) {
+    std::cerr << "ASMu2D::tesselate does not support different tesselation resolution in "
+              << "u- and v-direction. nviz u = " << npe[0] << ", nviz v = " << npe[1] 
+              << ", nviz w = " << npe[2] << std::endl;
+    return false;
+  }
 
-	int nNodesPerElement =  npe[0]   * npe[1]   * npe[2];
-	int nSubElPerElement = (npe[0]-1)*(npe[1]-1)*(npe[2]-1);
-	int nElements        = lrspline->nElements();
+  int nNodesPerElement =  npe[0]   * npe[1]   * npe[2];
+  int nSubElPerElement = (npe[0]-1)*(npe[1]-1)*(npe[2]-1);
+  int nElements        = lrspline->nElements();
 
-	// output is written once for each element resulting in a lot of unnecessary storage
-	// this is preferable to figuring out all element topology information
-	grid.unStructResize(nElements * nSubElPerElement,
-	                    nElements * nNodesPerElement);
+  // output is written once for each element resulting in a lot of unnecessary storage
+  // this is preferable to figuring out all element topology information
+  grid.unStructResize(nElements * nSubElPerElement,
+                      nElements * nNodesPerElement);
 
-	std::vector<LR::Element*>::iterator el;
-	int inod = 0;
-	int iel = 0;
-	for(el=lrspline->elementBegin(); el<lrspline->elementEnd(); el++, iel++) {
-		// evaluate element at element corner points
-		double umin = (**el).umin();
-		double umax = (**el).umax();
-		double vmin = (**el).vmin();
-		double vmax = (**el).vmax();
-		double wmin = (**el).wmin();
-		double wmax = (**el).wmax();
-		for(int iw=0; iw<npe[2]; iw++) {
-			for(int iv=0; iv<npe[1]; iv++) {
-				for(int iu=0; iu<npe[0]; iu++) {
-					double u = umin + (umax-umin)/(npe[0]-1)*iu;
-					double v = vmin + (vmax-vmin)/(npe[1]-1)*iv;
-					double w = wmin + (wmax-wmin)/(npe[2]-1)*iw;
-					Go::Point pt;
-					lrspline->point(pt, u,v,w, iel, iu!=npe[0]-1, iv!=npe[1]-1, iw!=npe[2]-1);
-					for(int dim=0; dim<nsd; dim++)
-						grid.setCoor(inod, dim, pt[dim]);
-					inod++;
-				}
-			}
-		}
-	}
+  std::vector<LR::Element*>::iterator el;
+  int inod = 0;
+  int iel = 0;
+  for(el=lrspline->elementBegin(); el<lrspline->elementEnd(); el++, iel++) {
+    // evaluate element at element corner points
+    double umin = (**el).umin();
+    double umax = (**el).umax();
+    double vmin = (**el).vmin();
+    double vmax = (**el).vmax();
+    double wmin = (**el).wmin();
+    double wmax = (**el).wmax();
+    for(int iw=0; iw<npe[2]; iw++) {
+      for(int iv=0; iv<npe[1]; iv++) {
+        for(int iu=0; iu<npe[0]; iu++) {
+          double u = umin + (umax-umin)/(npe[0]-1)*iu;
+          double v = vmin + (vmax-vmin)/(npe[1]-1)*iv;
+          double w = wmin + (wmax-wmin)/(npe[2]-1)*iw;
+          Go::Point pt;
+          lrspline->point(pt, u,v,w, iel, iu!=npe[0]-1, iv!=npe[1]-1, iw!=npe[2]-1);
+          for(int dim=0; dim<nsd; dim++)
+            grid.setCoor(inod, dim, pt[dim]);
+          inod++;
+        }
+      }
+    }
+  }
 
-	int ip = 0;
-	iel = 0;
-	for(int i=0; i<lrspline->nElements(); i++) {
-		int iStart = i*nNodesPerElement;
-		for(int iw=0; iw<npe[2]-1; iw++) {
-			for(int iv=0; iv<npe[1]-1; iv++) {
-				for(int iu=0; iu<npe[0]-1; iu++, iel++) {
-					// enumerate nodes counterclockwise around the hex
-					grid.setNode(ip++, iStart + (iw  )*npe[0]*npe[1] + (iv  )*npe[0] + (iu  ) );
-					grid.setNode(ip++, iStart + (iw  )*npe[0]*npe[1] + (iv  )*npe[0] + (iu+1) );
-					grid.setNode(ip++, iStart + (iw  )*npe[0]*npe[1] + (iv+1)*npe[0] + (iu+1) );
-					grid.setNode(ip++, iStart + (iw  )*npe[0]*npe[1] + (iv+1)*npe[0] + (iu  ) );
-					grid.setNode(ip++, iStart + (iw+1)*npe[0]*npe[1] + (iv  )*npe[0] + (iu  ) );
-					grid.setNode(ip++, iStart + (iw+1)*npe[0]*npe[1] + (iv  )*npe[0] + (iu+1) );
-					grid.setNode(ip++, iStart + (iw+1)*npe[0]*npe[1] + (iv+1)*npe[0] + (iu+1) );
-					grid.setNode(ip++, iStart + (iw+1)*npe[0]*npe[1] + (iv+1)*npe[0] + (iu  ) );
-					grid.setElmId(iel+1, i+1);
-				}
-			}
-		}
-	}
+  int ip = 0;
+  iel = 0;
+  for(int i=0; i<lrspline->nElements(); i++) {
+    int iStart = i*nNodesPerElement;
+    for(int iw=0; iw<npe[2]-1; iw++) {
+      for(int iv=0; iv<npe[1]-1; iv++) {
+        for(int iu=0; iu<npe[0]-1; iu++, iel++) {
+          // enumerate nodes counterclockwise around the hex
+          grid.setNode(ip++, iStart + (iw  )*npe[0]*npe[1] + (iv  )*npe[0] + (iu  ) );
+          grid.setNode(ip++, iStart + (iw  )*npe[0]*npe[1] + (iv  )*npe[0] + (iu+1) );
+          grid.setNode(ip++, iStart + (iw  )*npe[0]*npe[1] + (iv+1)*npe[0] + (iu+1) );
+          grid.setNode(ip++, iStart + (iw  )*npe[0]*npe[1] + (iv+1)*npe[0] + (iu  ) );
+          grid.setNode(ip++, iStart + (iw+1)*npe[0]*npe[1] + (iv  )*npe[0] + (iu  ) );
+          grid.setNode(ip++, iStart + (iw+1)*npe[0]*npe[1] + (iv  )*npe[0] + (iu+1) );
+          grid.setNode(ip++, iStart + (iw+1)*npe[0]*npe[1] + (iv+1)*npe[0] + (iu+1) );
+          grid.setNode(ip++, iStart + (iw+1)*npe[0]*npe[1] + (iv+1)*npe[0] + (iu  ) );
+          grid.setElmId(iel+1, i+1);
+        }
+      }
+    }
+  }
 
-	return true;
+  return true;
 }