diff --git a/cmake/Modules/FindIFEMDeps.cmake b/cmake/Modules/FindIFEMDeps.cmake
index bfb5c909..cdb81d1e 100644
--- a/cmake/Modules/FindIFEMDeps.cmake
+++ b/cmake/Modules/FindIFEMDeps.cmake
@@ -113,18 +113,18 @@ ENDIF(IFEM_USE_SUPERLU OR IFEM_USE_SUPERLU_MT)
 IF(IFEM_USE_LRSPLINES)
   FIND_PACKAGE(LRSpline)
   IF(LRSpline_LIBRARIES AND LRSpline_INCLUDE_DIRS)
-    IF(LRSPLINE_VERSION_MINOR LESS 4 AND LRSPLINE_VERSION_MAJOR LESS 1)
+    IF(LRSPLINE_VERSION_MINOR LESS 5 AND LRSPLINE_VERSION_MAJOR LESS 1)
       MESSAGE(STATUS "LR-spline library is out of date.
-                      Found ${LRSPLINE_VERSION_MAJOR}.${LRSPLINE_VERSION_MINOR}.${LRSPLINE_VERSION_PATCH}, need at least 0.4.0.
+                      Found ${LRSPLINE_VERSION_MAJOR}.${LRSPLINE_VERSION_MINOR}.${LRSPLINE_VERSION_PATCH}, need at least 0.5.0.
                       Support not enabled")
       SET(LRSpline_LIBRARIES "")
       SET(LRSpline_INCLUDE_DIRS "")
-    ELSE(LRSPLINE_VERSION_MINOR LESS 4 AND LRSPLINE_VERSION_MAJOR LESS 1)
+    ELSE(LRSPLINE_VERSION_MINOR LESS 5 AND LRSPLINE_VERSION_MAJOR LESS 1)
       SET(IFEM_DEPLIBS ${IFEM_DEPLIBS} ${LRSpline_LIBRARIES})
       SET(IFEM_DEPINCLUDES ${IFEM_DEPINCLUDES} ${LRSpline_INCLUDE_DIRS})
       SET(IFEM_BUILD_CXX_FLAGS "${IFEM_BUILD_CXX_FLAGS} -DHAS_LRSPLINE=1 ${LRSpline_DEFINITIONS}")
       SET(IFEM_CXX_FLAGS "${IFEM_CXX_FLAGS} -DHAS_LRSPLINE=1 ${LRSpline_DEFINITIONS}")
-    ENDIF(LRSPLINE_VERSION_MINOR LESS 4 AND LRSPLINE_VERSION_MAJOR LESS 1)
+    ENDIF(LRSPLINE_VERSION_MINOR LESS 5 AND LRSPLINE_VERSION_MAJOR LESS 1)
   ENDIF(LRSpline_LIBRARIES AND LRSpline_INCLUDE_DIRS)
 ENDIF(IFEM_USE_LRSPLINES)
 
diff --git a/src/ASM/ASMunstruct.C b/src/ASM/ASMunstruct.C
index ce6a4d43..3001f06e 100644
--- a/src/ASM/ASMunstruct.C
+++ b/src/ASM/ASMunstruct.C
@@ -26,7 +26,7 @@ int ASMunstruct::gNod = 0;
 
 
 ASMunstruct::ASMunstruct (unsigned char n_p, unsigned char n_s,
-			  unsigned char n_f)
+                          unsigned char n_f)
   : ASMbase(n_p,n_s,n_f)
 {
   geo = NULL;
@@ -53,15 +53,24 @@ ASMunstruct::~ASMunstruct ()
 bool ASMunstruct::refine (const RealArray&, const IntVec&, Vectors*, const char*)
 {
   std::cerr <<" *** ASMunstruct::refine: Not available without LR-Spline"
-	    << std::endl;
+            << std::endl;
   return false;
 }
 
 bool ASMunstruct::refine (const IntVec&, const IntVec&, Vectors*, const char*)
 {
   std::cerr <<" *** ASMunstruct::refine: Not available without LR-Spline"
-	    << std::endl;
+            << std::endl;
   return false;
 }
 
+
+Go::BsplineBasis ASMunstruct::getBezierBasis (int p)
+{
+  std::vector<double> knot(2*p);
+  std::fill(knot.begin(),   knot.begin()+p, -1.0);
+  std::fill(knot.begin()+p, knot.end()    ,  1.0);
+  return Go::BsplineBasis(p,p,knot.begin());
+}
+
 #endif
diff --git a/src/ASM/ASMunstruct.h b/src/ASM/ASMunstruct.h
index c29a501d..2c57940d 100644
--- a/src/ASM/ASMunstruct.h
+++ b/src/ASM/ASMunstruct.h
@@ -15,6 +15,7 @@
 #define _ASM_UNSTRUCT_H
 
 #include "ASMbase.h"
+#include "GoTools/geometry/BsplineBasis.h"
 
 namespace LR {
   class LRSpline;
@@ -79,6 +80,11 @@ public:
   //! \param[in] integrand Object with problem-specific data and methods
   virtual LR::LRSpline* evalSolution(const IntegrandBase& integrand) const = 0;
 
+  /*!
+    \brief Returns a Bezier basis of order \a p.
+  */
+  static Go::BsplineBasis getBezierBasis (int p);
+
 protected:
   LR::LRSpline* geo; //!< Pointer to the actual spline geometry object
 
diff --git a/src/ASM/LR/ASMLRSpline.C b/src/ASM/LR/ASMLRSpline.C
index 3b34a01b..8e6d00df 100644
--- a/src/ASM/LR/ASMLRSpline.C
+++ b/src/ASM/LR/ASMLRSpline.C
@@ -126,8 +126,10 @@ bool ASMunstruct::refine (const RealArray& elementError,
 		bool isLinIndep = geo->isLinearIndepByOverloading(false);
 		if(!isLinIndep) {
 			std::cout << "Inconclusive..." << std::endl;
+#ifdef HAS_BOOST
 			std::cout << "Testing for linear independence by full tensor expansion " << std::endl;
 			isLinIndep = geo->isLinearIndepByMappingMatrix(false);
+#endif
 		}
 		if(!isLinIndep)
 		{
@@ -259,8 +261,10 @@ bool ASMunstruct::refine (const IntVec& elements,
 		bool isLinIndep = geo->isLinearIndepByOverloading(false);
 		if(!isLinIndep) {
 			std::cout << "Inconclusive..." << std::endl;
+#ifdef HAS_BOOST
 			std::cout << "Testing for linear independence by full tensor expansion " << std::endl;
 			isLinIndep = geo->isLinearIndepByMappingMatrix(false);
+#endif
 		}
 		if(!isLinIndep)
 		{
@@ -286,3 +290,11 @@ bool ASMunstruct::refine (const IntVec& elements,
 
 	return true;
 }
+
+Go::BsplineBasis ASMunstruct::getBezierBasis (int p)
+{
+	double knot[2*p];
+	std::fill(knot,   knot+p,  -1.0);
+	std::fill(knot+p, knot+2*p, 1.0);
+	return Go::BsplineBasis(p,p,knot);
+}
diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C
index 5a555d63..586e3cda 100644
--- a/src/ASM/LR/ASMu2D.C
+++ b/src/ASM/LR/ASMu2D.C
@@ -297,6 +297,82 @@ bool ASMu2D::raiseOrder (int ru, int rv)
 }
 
 
+void ASMu2D::evaluateBasis(FiniteElement &fe, int derivs) const
+{
+	// skipping all error testing. Its probably really really smart
+
+	LR::Element *el = lrspline->getElement(fe.iel-1);
+	fe.xi  = 2.0*(fe.u - el->umin()) / (el->umax() - el->umin()) - 1.0; // sets value between -1 and 1
+	fe.eta = 2.0*(fe.v - el->vmin()) / (el->vmax() - el->vmin()) - 1.0;
+	std::vector<double> Nu = bezier_u.computeBasisValues(fe.xi,  derivs);
+	std::vector<double> Nv = bezier_v.computeBasisValues(fe.eta, derivs);
+
+	const Matrix &C = bezierExtract[fe.iel-1];
+
+	int p1 = lrspline->order(0);
+	int p2 = lrspline->order(1);
+	int p  = p1*p2;
+	Vector B(p);  // Bezier basis functions (vector of p1 x p2 components)
+	Vector Bu(p); // Bezier basis functions differentiated wrt u
+	Vector Bv(p); // Bezier basis functions differentiated wrt v
+	size_t k = 0;
+	for(size_t j=0; j<Nv.size(); j+=(derivs+1)) {
+		for(size_t i=0; i<Nu.size(); i+=(derivs+1), k++) {
+			B[k] = Nu[i]*Nv[j];
+			if(derivs > 0) {
+				Bu[k] = Nu[i+1]*Nv[ j ];
+				Bv[k] = Nu[ i ]*Nv[j+1];
+			}
+		}
+	}
+	fe.N = C*B;
+
+	int N    = el->nBasisFunctions();
+	int allP = p1*p2;
+	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 #" << 1 << std::endl;
+		exit(123);
+	}
+	sum = 0;
+	for(int qq=1; qq<=allP; qq++) sum+= B(qq);
+	if(fabs(sum-1) > 1e-10) {
+		std::cerr << "Bezier basis not sums to one at integration point #" << 1 << std::endl;
+		exit(123);
+	}
+
+	if(derivs > 0) {
+		Matrix dNdu(el->nBasisFunctions(), 2);
+		dNdu.fillColumn(1, C*Bu);
+		dNdu.fillColumn(2, C*Bv);
+		Matrix Xnod, Jac;
+		getElementCoordinates(Xnod, fe.iel);
+		fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
+
+		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 #" << 1 << 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 #" << 1 << std::endl;
+			exit(123);
+		}
+
+		sum = 0;
+		for(int qq=1; qq<=allP; qq++) sum+= Bu(qq);
+		if(fabs(sum) > 1e-10) {
+			std::cerr << "Bezier derivatives not sums to zero at integration point #" << 1 << std::endl;
+			exit(123);
+		}
+	}
+
+}
+
 bool ASMu2D::generateFEMTopology ()
 {
 	// At this point we are through with the tensor spline object.
@@ -322,24 +398,52 @@ bool ASMu2D::generateFEMTopology ()
 	else if (shareFE)
 		return true;
 
+	const int p1 = lrspline->order(0);
+	const int p2 = lrspline->order(1);
+
 	myMLGN.resize(nBasis);
 	myMLGE.resize(nElements);
 	myMNPC.resize(nElements);
+	bezierExtract.resize(nElements);
 	lrspline->generateIDs();
 
+	std::vector<double> extrMat ;
 	std::vector<LR::Element*>::iterator el_it = lrspline->elementBegin();
 	for (int iel=0; iel<nElements; iel++, el_it++)
 	{
+		LR::Element *el = *el_it;
+		int nSupportFunctions = el->nBasisFunctions();
+
 		myMLGE[iel] = ++gEl; // global element number over all patches
-		myMNPC[iel].resize((*el_it)->nBasisFunctions());
+		myMNPC[iel].resize(el->nBasisFunctions());
 
 		int lnod = 0;
-		for (LR::Basisfunction *b : (*el_it)->support())
+		for (LR::Basisfunction *b : el->support())
 			myMNPC[iel][lnod++] = b->getId();
+
+		{
+		PROFILE("Bezier extraction");
+
+		// get bezier extraction matrix 
+		lrspline->getBezierExtraction(iel, extrMat);
+		int width  = p1*p2;
+		int height = nSupportFunctions;
+
+		// wrap the extraction data into a Matrix class
+		Matrix newM(height, width);
+		for(int c=1; c<=width; c++) 
+			newM.fillColumn(c, &extrMat[(c-1)*height]);
+
+		// keep it for use later
+		bezierExtract[iel] = newM;
+		}
 	}
 
 	for (int inod = 0; inod < nBasis; inod++)
 		myMLGN[inod] = ++gNod;
+	
+	bezier_u = getBezierBasis(p1);
+	bezier_v = getBezierBasis(p2);
 
 	nnod = gNod;
 	return true;
@@ -1189,10 +1293,18 @@ int ASMu2D::evalPoint (const double* xi, double* param, Vec3& X) const
 	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);
 
-	Go::Point X0;
-	lrspline->point(X0,param[0],param[1]);
-	for (unsigned char d = 0; d < nsd; d++)
-		X[d] = X0[d];
+	int iel = lrspline->getElementContaining(param[0],param[1]);
+	FiniteElement fe(lrspline->getElement(iel)->nBasisFunctions());
+	fe.iel = iel + 1;
+	fe.u   = param[0];
+	fe.v   = param[1];
+
+	Matrix Xnod;
+	getElementCoordinates(Xnod, fe.iel);
+
+	evaluateBasis(fe);
+	X = Xnod * fe.N;
+
 	return 0;
 }
 
@@ -1385,6 +1497,10 @@ bool ASMu2D::evalSolution (Matrix& sField, const Vector& locSol,
     // Fetch element containing evaluation point.
     // Sadly, points are not always ordered in the same way as the elements.
     int iel = lrspline->getElementContaining(gpar[0][i],gpar[1][i]);
+    FiniteElement fe(lrspline->getElement(iel)->nBasisFunctions());
+    fe.iel = iel + 1;
+    fe.u   = gpar[0][i];
+    fe.v   = gpar[1][i];
     utl::gather(MNPC[iel],nComp,locSol,eSol);
 
     // Set up control point (nodal) coordinates for current element
@@ -1396,16 +1512,14 @@ bool ASMu2D::evalSolution (Matrix& sField, const Vector& locSol,
     switch (deriv) {
 
     case 0: // Evaluate the solution
-      lrspline->computeBasis(gpar[0][i],gpar[1][i],spline0,iel);
-      eSol.multiply(spline0.basisValues,ptSol);
+      evaluateBasis(fe, deriv);
+      ptSol = eSol * fe.N;
       sField.fillColumn(1+i,ptSol);
       break;
 
     case 1: // Evaluate first derivatives of the solution
-      lrspline->computeBasis(gpar[0][i],gpar[1][i],spline1,iel);
-      SplineUtils::extractBasis(spline1,ptSol,dNdu);
-      utl::Jacobian(Jac,dNdX,Xnod,dNdu);
-      ptDer.multiply(eSol,dNdX);
+      evaluateBasis(fe, deriv);
+      ptDer.multiply(eSol,fe.dNdX);
       sField.fillColumn(1+i,ptDer);
       break;
 
diff --git a/src/ASM/LR/ASMu2D.h b/src/ASM/LR/ASMu2D.h
index 33b4ff94..ae074058 100644
--- a/src/ASM/LR/ASMu2D.h
+++ b/src/ASM/LR/ASMu2D.h
@@ -15,6 +15,7 @@
 #define _ASM_U2D_H
 
 #include "ASMunstruct.h"
+#include "FiniteElement.h"
 #include "ASM2D.h"
 
 #include <memory>
@@ -383,6 +384,9 @@ protected:
   //! \param[out] XC Coordinates of the element corners
   void getElementCorners(int iel, std::vector<Vec3>& XC) const;
 
+  //! \brief Evaluate all basis functions and \a derivs number of derivatives on one element
+  virtual void evaluateBasis(FiniteElement &el, int derivs=0) const;
+
 public:
   //! \brief Returns the number of elements on a boundary.
   virtual size_t getNoBoundaryElms(char lIndex, char ldim) const;
@@ -414,6 +418,8 @@ protected:
 
   std::shared_ptr<LR::LRSplineSurface> lrspline; //!< Pointer to the LR-spline surface object
 
+  std::vector<Matrix> bezierExtract; //!< Bezier extraction matrices
+
   Go::SplineSurface* tensorspline; //!< Pointer to original tensor spline object
   // The tensor spline object is kept for backward compatability with the REFINE
   // and RAISEORDER key-words, although we take note that there is a possibility
@@ -422,6 +428,9 @@ protected:
 
   //! Inhomogeneous Dirichlet boundary condition data
   std::vector<DirichletEdge> dirich;
+  
+  Go::BsplineBasis bezier_u; // we're keeping one of these objects to do evaluation on
+  Go::BsplineBasis bezier_v; 
 };
 
 #endif
diff --git a/src/ASM/LR/ASMu3D.C b/src/ASM/LR/ASMu3D.C
index 761d413f..55b3438c 100644
--- a/src/ASM/LR/ASMu3D.C
+++ b/src/ASM/LR/ASMu3D.C
@@ -816,18 +816,6 @@ void ASMu3D::getElementCorners (int iEl, Vec3Vec& XC) const
       }
 }
 
-/*!
-  \brief Returns a Bezier basis of order \a p.
-*/
-
-static Go::BsplineBasis getBezierBasis (int p)
-{
-	double knot[2*p];
-	std::fill(knot,   knot+p,  -1.0);
-	std::fill(knot+p, knot+2*p, 1.0);
-	return Go::BsplineBasis(p,p,knot);
-}
-
 
 void ASMu3D::evaluateBasis (FiniteElement &el, Matrix &dNdu) const
 {
diff --git a/src/SIM/SIMgeneric.C b/src/SIM/SIMgeneric.C
index 1348bbe0..6ebbc704 100644
--- a/src/SIM/SIMgeneric.C
+++ b/src/SIM/SIMgeneric.C
@@ -46,3 +46,14 @@ Vector SIMgeneric::getSolution (const Vector& psol, const double* par,
 
   return tmpVal.getColumn(1);
 }
+
+int SIMgeneric::evalPoint(const double* xi, Vec3& X, double* param, int patch) const
+{
+  ASMbase *pch = getPatch(patch-1);
+
+  double dummy[2];
+  if(param == NULL)
+    return pch->evalPoint(xi, dummy, X);
+  else
+    return pch->evalPoint(xi, param, X);
+}
diff --git a/src/SIM/SIMgeneric.h b/src/SIM/SIMgeneric.h
index 7733bbbb..bb4a2374 100644
--- a/src/SIM/SIMgeneric.h
+++ b/src/SIM/SIMgeneric.h
@@ -45,6 +45,15 @@ public:
   //! \return Evaluated solution values
   Vector getSolution(const Vector& psol, const double* par,
                      int deriv = 0, int patch = 1) const;
+
+  //! \brief Evaluates the mapping of the geometry at the given point.
+  //! \param[in] xi Dimensionless parameters in range [0,1] of the point
+  //! \param[out] X The Cartesian coordinates of the point
+  //! \param[out] param The parameters of the point in the knot-span domain
+  //! \param[in] patch The patch to evaluate
+  //! \return 0 if the evaluation went good
+  int evalPoint(const double* xi, Vec3& X, 
+                double* param=NULL, int patch = 1) const;
 };
 
 #endif