diff --git a/examples/biceps_activeMuscle_3D.cpp b/examples/biceps_activeMuscle_3D.cpp
index 64d0592..6ffb6da 100644
--- a/examples/biceps_activeMuscle_3D.cpp
+++ b/examples/biceps_activeMuscle_3D.cpp
@@ -133,6 +133,7 @@ int main(int argc, char* argv[]){
     // setting Newton's method
     gsIterative<real_t> solver(assembler);
     solver.options().setInt("Verbosity",solver_verbosity::all);
+    solver.options().setReal("RelTol",1e-4); // Should be 1e-9
     solver.options().setInt("Solver",linear_solver::LDLT);
 
     gsInfo << "Solving...\n";
diff --git a/examples/muscleBeam_activeMuscle_3Dt.cpp b/examples/muscleBeam_activeMuscle_3Dt.cpp
index 90234ef..5cadcb3 100644
--- a/examples/muscleBeam_activeMuscle_3Dt.cpp
+++ b/examples/muscleBeam_activeMuscle_3Dt.cpp
@@ -40,7 +40,7 @@ int main(int argc, char* argv[]){
     gsVector<> fiberDirection(3);
     fiberDirection << 1.,0.,0.;
     // space discretization
-    index_t numUniRefDirX = 2;
+    index_t numUniRefDirX = 0;
     index_t numUniRef = 0;
     index_t numDegElev = 0;
     bool subgridOrTaylorHood = false;
diff --git a/src/gsVisitorElasticityNeumann.h b/src/gsVisitorElasticityNeumann.h
index 87305cb..992b5a9 100644
--- a/src/gsVisitorElasticityNeumann.h
+++ b/src/gsVisitorElasticityNeumann.h
@@ -36,7 +36,7 @@ class gsVisitorElasticityNeumann
                     const gsOptionList & options,
                     gsQuadRule<T> & rule)
     {
-        GISMO_UNUSED(patchIndex);
+        currentPatch = patchIndex;
         // parametric dimension of the first displacement component
         dim = basisRefs.front().dim();
         // a quadrature rule is defined by the basis for the first displacement component.
@@ -59,8 +59,8 @@ class gsVisitorElasticityNeumann
         md.flags = NEED_VALUE | NEED_MEASURE;
         // Compute image of the quadrature points plus gradient, jacobian and other necessary data
         geo.computeMap(md);
-        // Evaluate the Neumann functon on the images of the quadrature points
-        neumannFunction_ptr->eval_into(md.values[0], neumannValues);
+
+        neumannFunction_ptr->piece(currentPatch).eval_into(md.values[0], neumannValues);
         // find local indices of the displacement basis functions active on the element
         basisRefs.front().active_into(quNodes.col(0),localIndicesDisp);
         N_D = localIndicesDisp.rows();
@@ -107,6 +107,7 @@ class gsVisitorElasticityNeumann
     // problem info
     short_t dim;
     const gsFunctionSet<T> * neumannFunction_ptr;
+    index_t currentPatch; //gsLOOKUPFUNCTION
     T forceScaling;
     boxSide patchSide;
     // geometry mapping
diff --git a/src/gsVisitorThermoBoundary.h b/src/gsVisitorThermoBoundary.h
index 988af25..f3c9819 100644
--- a/src/gsVisitorThermoBoundary.h
+++ b/src/gsVisitorThermoBoundary.h
@@ -69,10 +69,7 @@ class gsVisitorThermoBoundary
         // Compute image of the quadrature points plus gradient, jacobian and other necessary data
         geo.computeMap(md);
         // Evaluate temperature
-        if (paramTemp) // evaluate temperature in the parametric domain
-            temperatureField.piece(patch).eval_into(quNodes,tempValues);
-        else           // evaluate temperature in the physical domain
-            temperatureField.eval_into(md.values[0],tempValues);
+        temperatureField.piece(patch).eval_into(quNodes,tempValues);
         // find local indices of the displacement basis functions active on the element
         basisRefs.front().active_into(quNodes.col(0),localIndicesDisp);
         N_D = localIndicesDisp.rows();