Make compile by using GeometryGrid to implement the deformation of the grid

[[Imported from SVN: r12631]]

nladaptmeasure.cc

@@ -5,6 +5,7 @@
 #include <dune/common/parametertreeparser.hh>
 
 #include <dune/grid/uggrid.hh>
+#include <dune/grid/geometrygrid.hh>
 #include <dune/grid/io/file/amirameshwriter.hh>
 #include <dune/grid/io/file/amirameshreader.hh>
 
@@ -20,6 +21,7 @@
 #include <dune/fufem/assemblers/operatorassembler.hh>
 #include <dune/fufem/functiontools/gridfunctionadaptor.hh>
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/functions/deformationfunction.hh>
 
 #ifdef HAVE_IPOPT
 #include <dune/solvers/solvers/quadraticipopt.hh>
@@ -36,9 +38,9 @@
 #include <dune/fufem/estimators/fractionalmarking.hh>
 
 #include <dune/contact/estimators/hierarchiccontactestimator.hh>
+#include <dune/contact/assemblers/twobodyassembler.hh>
 #include <dune/contact/solvers/nsnewtonmgstep.hh>
 #include <dune/contact/solvers/nsnewtoncontacttransfer.hh>
-#include <dune/contact/assemblers/nonlintwobodyassembler.hh>
 
 using namespace Dune;
 
@@ -51,8 +53,11 @@ typedef double field_type;
 typedef BCRSMatrix<FieldMatrix<double, dim, dim> > MatrixType;
 typedef BlockVector<FieldVector<double, dim> >     VectorType;
 typedef UGGrid<dim> GridType;
+typedef DeformationFunction<GridType::LeafGridView,VectorType> Deformation;
+typedef GeometryGrid<GridType,Deformation> DeformedGridType;
 typedef BoundaryPatch<GridType::LevelGridView> LevelBoundary;
 typedef BoundaryPatch<GridType::LeafGridView> LeafBoundary;
+typedef BoundaryPatch<DeformedGridType::LevelGridView> DefLevelBoundary;
 
 const double E[2]  = {17e6,1e6};
 const double nu[2] = {0.3, 0.3};
@@ -103,12 +108,13 @@ void setTrustRegionObstacles(double trustRegionRadius,
 
 
 template <int numGrids>
-void compute (const array<GridType, numGrids>& grid,
+void compute (const array<DeformedGridType*, numGrids>& grid,
+              array<Deformation*,numGrids>& deformation,
               array<VectorType,numGrids>& x,
               const array<VectorType,numGrids>& coarseDirichletValues, 
               const array<BitSetVector<1>,numGrids>& coarseDirichletNodes,
               const array<LevelBoundary, numGrids>& coarseDirichletBoundary,
-              const LevelBoundary& obsPatch,
+              const DefLevelBoundary& obsPatch,
               double obsDistance,
               double trustRegionRadius,
               int maxTrustRegionSteps,
@@ -120,38 +126,38 @@ void compute (const array<GridType, numGrids>& grid,
     array<std::vector<BitSetVector<dim> >, 2> dirichletNodes;
     array<std::vector<LevelBoundary>, 2> dirichletBoundary;
         
-    array<VectorType, 2> rhs;
+    std::vector<VectorType> rhs(numGrids);
     
     for (int i=0; i<numGrids; i++) {
         
-        dirichletValues[i].resize(grid[i].maxLevel()+1);
+        dirichletValues[i].resize(grid[i]->maxLevel()+1);
         dirichletValues[i][0] = coarseDirichletValues[i];
         
         /** \todo We don't need the whole hierarchy, do we? */
-        dirichletNodes[i].resize(grid[i].maxLevel()+1);
+        dirichletNodes[i].resize(grid[i]->maxLevel()+1);
         
-        dirichletBoundary[i].resize(grid[i].maxLevel()+1);
-        dirichletBoundary[i][0].setup(grid[i].levelView(0), coarseDirichletNodes[i]);
+        dirichletBoundary[i].resize(grid[i]->maxLevel()+1);
+        dirichletBoundary[i][0].setup(grid[i]->hostGrid().levelGridView(0), coarseDirichletNodes[i]);
         
         // Initial solution
-        assert(x[i].size() == grid[i].size(dim));
+        assert(x[i].size() == grid[i]->size(dim));
         
         BoundaryPatchProlongator<GridType>::prolong(dirichletBoundary[i]);
         
-        for (int j=0; j<=grid[i].maxLevel(); j++) {
+        for (int j=0; j<=grid[i]->maxLevel(); j++) {
             
-            int fSSize = grid[i].size(j,dim);
+            int fSSize = grid[i]->size(j,dim);
             dirichletNodes[i][j].resize(fSSize);
             for (int k=0; k<fSSize; k++)
                 dirichletNodes[i][j][k] = dirichletBoundary[i][j].containsVertex(k);
             
         }
         
-        sampleOnBitField(grid[i], dirichletValues[i], dirichletNodes[i]);
+        sampleOnBitField(*grid[i], dirichletValues[i], dirichletNodes[i]);
         
         //    Create rhs vectors
         
-        rhs[i].resize(grid[i].size(grid[i].maxLevel(),dim));
+        rhs[i].resize(grid[i]->size(grid[i]->maxLevel(),dim));
         rhs[i] = 0;
         
         // Copy Dirichlet information into the solution vector
@@ -164,7 +170,7 @@ void compute (const array<GridType, numGrids>& grid,
         
     }
 
-    int toplevel = std::max(grid[0].maxLevel(), grid[1].maxLevel());
+    int toplevel = std::max(grid[0]->maxLevel(), grid[1]->maxLevel());
     
     
     // make dirichlet bitfields containing dirichlet information for both grids
@@ -172,8 +178,8 @@ void compute (const array<GridType, numGrids>& grid,
     
     for (int colevel=0; colevel<=toplevel; colevel++) {
         
-        int levelGrid0 = std::max(0, grid[0].maxLevel()-colevel);
-        int levelGrid1 = std::max(0, grid[1].maxLevel()-colevel);
+        int levelGrid0 = std::max(0, grid[0]->maxLevel()-colevel);
+        int levelGrid1 = std::max(0, grid[1]->maxLevel()-colevel);
         
         int offset = dirichletValues[0][levelGrid0].size();
         
@@ -189,17 +195,13 @@ void compute (const array<GridType, numGrids>& grid,
 
 
     // Assemble contact problem
-    NonlinearTwoBodyAssembler<GridType, VectorType, false> contactAssembler(grid[0],
-                                                                            grid[1],
-                                                                            obsPatch,
-                                                                            obsDistance,
-                                                                            CouplingPairBase::CONTACT
-#ifdef EXPLICIT_CONTACT_DIRECTIONS
-                                                                            , &contactDirection
-#endif
-                                                                            );
-    BitSetVector<1> hasObstacle(grid[0].size(dim) + grid[1].size(dim), false);
-    std::vector<BoxConstraint<field_type,dim> > trustRegionObstacles(grid[0].size(dim) + grid[1].size(dim));
+    DefLevelBoundary mortarPatch(grid[1]->levelGridView(0),true);
+    TwoBodyAssembler<DeformedGridType, VectorType, false> contactAssembler(*grid[0], *grid[1],
+                                                                   obsPatch, mortarPatch, obsDistance,
+                                                                   CouplingPairBase::CONTACT);
+
+    BitSetVector<1> hasObstacle(grid[0]->size(dim) + grid[1]->size(dim), false);
+    std::vector<BoxConstraint<field_type,dim> > trustRegionObstacles(grid[0]->size(dim) + grid[1]->size(dim));
 
 
     // //////////////////////////////////////
@@ -228,7 +230,7 @@ void compute (const array<GridType, numGrids>& grid,
     multigridStep.verbosity_         = Solver::QUIET;
 
     // Create the transfer operators
-    multigridStep.mgTransfer_.resize(grid[0].maxLevel());
+    multigridStep.mgTransfer_.resize(grid[0]->maxLevel());
     for (int i=0; i<multigridStep.mgTransfer_.size(); i++)
         multigridStep.mgTransfer_[i] = NULL;
 
@@ -240,19 +242,6 @@ void compute (const array<GridType, numGrids>& grid,
                                                    &energyNorm,
                                                    Solver::QUIET);
 
-    // ///////////////////////////////////////////////////
-    //   Do a homotopy of the Dirichlet boundary data
-    // ///////////////////////////////////////////////////
-    typedef OgdenAssembler<GridType> Assembler;
-    //typedef NeoHookeAssembler<GridType, 3> Assembler;
-
-    array<Assembler,2> ogdenAssembler;
-    
-    for (int i=0; i<2; i++) {
-        ogdenAssembler[i].grid_ = &grid[i];
-        ogdenAssembler[i].setEandNu(E[i], nu[i]);
-        ogdenAssembler[i].d_ = d[i];
-    }
     
     // //////////////////////////////////////////////////////////
     //   Assemble linear elasticity matrix for the energy norm
@@ -260,8 +249,8 @@ void compute (const array<GridType, numGrids>& grid,
     // //////////////////////////////////////////////////////////
     // Assemble separate linear elasticity problems
     typedef P1NodalBasis<GridType::LeafGridView> P1Basis;
-    P1Basis p1Basis0(grid[0].leafView());
-    P1Basis p1Basis1(grid[1].leafView());
+    P1Basis p1Basis0(grid[0]->hostGrid().leafView());
+    P1Basis p1Basis1(grid[1]->hostGrid().leafView());
     
     OperatorAssembler<P1Basis,P1Basis> globalAssembler0(p1Basis0,p1Basis0);
     OperatorAssembler<P1Basis,P1Basis> globalAssembler1(p1Basis1,p1Basis1);
@@ -273,17 +262,28 @@ void compute (const array<GridType, numGrids>& grid,
     globalAssembler1.assemble(localAssembler, stiffnessMatrix1);
         
     array<const MatrixType*, 2> linearElasticityMatrices;
-    submat[0] = &stiffnessMatrix0;
-    submat[1] = &stiffnessMatrix1;
+    linearElasticityMatrices[0] = &stiffnessMatrix0;
+    linearElasticityMatrices[1] = &stiffnessMatrix1;
 
+    typedef OgdenAssembler<P1Basis, P1Basis> Assembler;
+    //typedef NeoHookeAssembler<P1Basis, P1Basis> Assembler;
+
+    array<Assembler*,2> ogdenAssembler;
+    ogdenAssembler[0] = new Assembler(p1Basis0,p1Basis0);
+    ogdenAssembler[1] = new Assembler(p1Basis1,p1Basis1);
+    
+    for (int i=0; i<2; i++) {
+        ogdenAssembler[i]->setEandNu(E[i], nu[i]);
+        ogdenAssembler[i]->d_ = d[i];
+    }
     // /////////////////////////////////////////////////////
     //   Trust-Region Solver
     // /////////////////////////////////////////////////////
 
     for (int i=0; i<maxTrustRegionSteps; i++) {
         
-        double oldEnergy = ogdenAssembler[0].computeEnergy(x[0]) 
-            + ogdenAssembler[1].computeEnergy(x[1]);
+        double oldEnergy = ogdenAssembler[0]->computeEnergy(x[0])
+            + ogdenAssembler[1]->computeEnergy(x[1]);
         
         if (std::isnan(oldEnergy))
             DUNE_THROW(UnfeasableInitialIterateException, "Unfeasable initial iterate!");
@@ -292,20 +292,24 @@ void compute (const array<GridType, numGrids>& grid,
         std::cout << "      Trust-Region Step Number: " << i << ",   Energy: " << oldEnergy << std::endl;
         std::cout << "------------------------------------------------------" << std::endl;
 
+        // Deform the grid
+        deformation[0]->setDeformation(x[0]);
+        deformation[1]->setDeformation(x[1]);
+
         // /////////////////////////////////////////////////////
         //   Assemble the obstacle and the transfer operator 
         // /////////////////////////////////////////////////////
-        contactAssembler.assembleObstacle(x);
-        contactAssembler.assembleTransferOperator(x);
+        contactAssembler.assembleObstacle();
+        contactAssembler.assembleTransferOperator();
 
         // /////////////////////////////////////////////////////
         //   Check if initial iterate is inadmissible
         // /////////////////////////////////////////////////////
         bool initialIterateIsAdmissible = true;
-        for (int j=0; j<contactAssembler.obstacles_[grid[0].maxLevel()].size() && initialIterateIsAdmissible; j++)
+        for (int j=0; j<contactAssembler.obstacles_[grid[0]->maxLevel()].size() && initialIterateIsAdmissible; j++)
             for (int k=0; k<dim; k++)
-                if (contactAssembler.obstacles_[grid[0].maxLevel()][j].lower(k) > 0
-                    || contactAssembler.obstacles_[grid[0].maxLevel()][j].upper(k) < 0) {
+                if (contactAssembler.obstacles_[grid[0]->maxLevel()][j].lower(k) > 0
+                    || contactAssembler.obstacles_[grid[0]->maxLevel()][j].upper(k) < 0) {
                     initialIterateIsAdmissible = false;
                     break;
                 }
@@ -317,8 +321,8 @@ void compute (const array<GridType, numGrids>& grid,
         // //////////////////////////////////////////////////////
         
         hasObstacle.unsetAll();
-        for (int j=0; j<grid[0].size(dim); j++)
-            hasObstacle[j] = contactAssembler.nonmortarBoundary_[0][grid[0].maxLevel()].containsVertex(j);
+        for (int j=0; j<grid[0]->size(dim); j++)
+            hasObstacle[j] = contactAssembler.nonmortarBoundary_[0][grid[0]->maxLevel()].containsVertex(j);
 
         // //////////////////////////////////////////////////////////////////////
         //   The MG transfer operators involve the mortar coupling operator
@@ -327,9 +331,9 @@ void compute (const array<GridType, numGrids>& grid,
         for (int j=0; j<multigridStep.mgTransfer_.size(); j++) {
         
             NonSmoothNewtonContactTransfer<VectorType>* newTransferOp = new NonSmoothNewtonContactTransfer<VectorType>;
-            newTransferOp->setup(grid[0], grid[1], j,
-                                 contactAssembler.contactCoupling_.mortarLagrangeMatrix_,
-                                 contactAssembler.localCoordSystems_[grid[0].maxLevel()-j], 
+            newTransferOp->setup(grid[0]->hostGrid(), grid[1]->hostGrid(), j,
+                                 contactAssembler.contactCoupling_[0].maxLevelMortarLagrangeMatrix(),
+                                 contactAssembler.localCoordSystems_[grid[0]->maxLevel()-j],
                                  hasObstacle);
             
             if (multigridStep.mgTransfer_[toplevel-j-1])
@@ -348,7 +352,7 @@ void compute (const array<GridType, numGrids>& grid,
         MatrixType transformedLinearMatrix;
         contactAssembler.assemble(linearElasticityMatrices, transformedLinearMatrix);
         
-        dynamic_cast<EnergyNorm<MatrixType,VectorType>*>(solver.errorNorm_)->setMatrix(&transformedLinearMatrix);
+        energyNorm.setMatrix(&transformedLinearMatrix);
     
 
         // ///////////////////////////////////////////////////////////////////////////////
@@ -357,14 +361,18 @@ void compute (const array<GridType, numGrids>& grid,
         rhs[0] = 0;
         rhs[1] = 0;
 
-        ogdenAssembler[0].assembleMatrix(x[0], rhs[0]);
-        ogdenAssembler[1].assembleMatrix(x[1], rhs[1]);
-        
+        std::vector<const MatrixType*> elastMatrix(numGrids);
+        for (int i=0; i<numGrids; i++) {
+            elastMatrix[i] = new MatrixType;
+            ogdenAssembler[i]->assembleProblem(*const_cast<MatrixType*>(elastMatrix[i]), x[i], rhs[i]);
+        }
 
-        contactAssembler.assembleJacobian(ogdenAssembler[0].getMatrix(), ogdenAssembler[1].getMatrix());
+        // Transform the nonlinear elasticity matrix into the mortar basis
+        MatrixType transformedMatrix;
+        contactAssembler.assembleJacobian(elastMatrix,transformedMatrix);
         
         VectorType totalRhs(rhs[0].size() + rhs[1].size());
-        contactAssembler.assembleRightHandSide(totalRhs, rhs);
+        contactAssembler.assembleRightHandSide(rhs, totalRhs);
         
         // The correction for both grids together in local coordinates
         VectorType totalCorr(x[0].size() + x[1].size());
@@ -384,11 +392,11 @@ void compute (const array<GridType, numGrids>& grid,
             // Create trust-region obstacle on grid0.maxLevel()
             setTrustRegionObstacles(trustRegionRadius,
                                     trustRegionObstacles,
-                                    contactAssembler.obstacles_[grid[0].maxLevel()],
-                                    totalDirichletNodes[grid[0].maxLevel()]);
+                                    contactAssembler.obstacles_[grid[0]->maxLevel()],
+                                    totalDirichletNodes[grid[0]->maxLevel()]);
             
             if (dynamic_cast<MultigridStep<MatrixType,VectorType>*>(solver.iterationStep_))
-                dynamic_cast<MultigridStep<MatrixType,VectorType>*>(solver.iterationStep_)->setProblem(*contactAssembler.mat_, totalCorr, totalRhs, grid[0].maxLevel()+1);
+                dynamic_cast<MultigridStep<MatrixType,VectorType>*>(solver.iterationStep_)->setProblem(transformedMatrix, totalCorr, totalRhs, grid[0]->maxLevel()+1);
             else {
                 DUNE_THROW(Exception, "You need a multigrid step!");
                 //solver.iterationStep_->setProblem(*contactAssembler.mat_, totalCorr, totalRhs);
@@ -422,8 +430,8 @@ void compute (const array<GridType, numGrids>& grid,
             VectorType newIterate0 = x[0];  newIterate0 += corr[0];
             VectorType newIterate1 = x[1];  newIterate1 += corr[1];
             
-            double energy = ogdenAssembler[0].computeEnergy(newIterate0) 
-                + ogdenAssembler[1].computeEnergy(newIterate1);
+            double energy = ogdenAssembler[0]->computeEnergy(newIterate0)
+                + ogdenAssembler[1]->computeEnergy(newIterate1);
             
             // //////////////////////////////////////////////////////
             //   Compute the model decrease
@@ -434,7 +442,7 @@ void compute (const array<GridType, numGrids>& grid,
             for (size_t j=0; j<2; j++) {
                 VectorType tmp(corr[j].size());
                 tmp = 0;
-                ogdenAssembler[j].getMatrix()->umv(corr[j], tmp);
+                elastMatrix[j]->umv(corr[j], tmp);
                 modelDecrease += (rhs[j]*corr[j]) - 0.5 * (corr[j]*tmp);
             }
 
@@ -486,6 +494,11 @@ void compute (const array<GridType, numGrids>& grid,
         x[0] += corr[0];
         x[1] += corr[1];
 
+        // cleanup
+        for (int i=0; i<numGrids;i++) {
+            delete(elastMatrix[i]);
+            delete(ogdenAssembler[i]);
+        }
     }
 
 }
@@ -564,13 +577,6 @@ int main(int argc, char *argv[]) try
     readBoundaryPatch<GridType>(coarseDirichletBoundary[0], path + dirichletFile0);
     readBoundaryPatch<GridType>(coarseDirichletBoundary[1], path + dirichletFile1);
 
-    // /////////////////////////////////////////////////////
-    // read field describing the obstacles
-    // /////////////////////////////////////////////////////
-
-    LevelBoundary obsPatch(uniformGrid[0].levelView(0));
-    readBoundaryPatch<GridType>(obsPatch, path + obsFilename);
-
     // ////////////////////////////////////////////////////////////////////
     //   Solve the problem on the coarsest level (using path following,
     //   if necessary), to get a good start iterate.
@@ -582,6 +588,21 @@ int main(int argc, char *argv[]) try
     uniformX[0] = 0;
     uniformX[1] = 0;
 
+    // create deformed grids
+    array<Deformation*,2> deformation;
+    array<DeformedGridType*,2> defUniformGrid;
+    for (int i=0;i<2;i++) {
+        deformation[i]=new Deformation(uniformGrid[i].leafGridView(),uniformX[i]);
+        defUniformGrid[i]= new DeformedGridType(&uniformGrid[i],deformation[i]);
+    }
+
+    // /////////////////////////////////////////////////////
+    // read field describing the obstacles
+    // /////////////////////////////////////////////////////
+
+    DefLevelBoundary obsPatch(defUniformGrid[0]->levelView(0));
+    readBoundaryPatch<DeformedGridType>(obsPatch, path + obsFilename);
+
     double loadFactor    = 0;
     double loadIncrement = 1;
 
@@ -617,7 +638,7 @@ int main(int argc, char *argv[]) try
                 std::cout << "####################################################" << std::endl;
 
                 try {
-                    compute<2>(uniformGrid, uniformX, scaledDirichletValues, dirichletNodes, coarseDirichletBoundary,
+                    compute<2>(defUniformGrid, deformation, uniformX, scaledDirichletValues, dirichletNodes, coarseDirichletBoundary,
                                obsPatch, obsDistance, trustRegionRadius, maxTrustRegionSteps, tolerance, numIt);
                 } catch (UnfeasableInitialIterateException e) {
                     std::cout << "Reducing path following parameter\n";
@@ -660,6 +681,10 @@ int main(int argc, char *argv[]) try
             
             if (paramBoundaries)
                 improveGrid(uniformGrid[j], 10);
+
+            deformation[j]->setDeformation(uniformX[j]);
+            deformation[j]->setGridView(uniformGrid[j].leafGridView());
+            defUniformGrid[j]->update();
         }
     }
 
@@ -668,8 +693,9 @@ int main(int argc, char *argv[]) try
     // /////////////////////////////////////////////////////////
 
     if (maxLevel!=0)
-        compute<2>(uniformGrid, uniformX, dirichletValues, dirichletNodes, coarseDirichletBoundary,
-                   obsPatch, obsDistance, trustRegionRadius, maxTrustRegionSteps, tolerance, numIt);
+        compute<2>(defUniformGrid, deformation, uniformX, dirichletValues, dirichletNodes,
+                    coarseDirichletBoundary, obsPatch, obsDistance, trustRegionRadius,
+                    maxTrustRegionSteps, tolerance, numIt);
 
     // Output result
     LeafAmiraMeshWriter<GridType>::writeGrid(uniformGrid[0], "0uniformGrid");
@@ -677,18 +703,11 @@ int main(int argc, char *argv[]) try
     
     if (dim==2) {
         
-        for (int i=0; i<2; i++) {
-
-            GridType::Codim<dim>::LeafIterator vIt    = uniformGrid[i].leafbegin<dim>();
-            GridType::Codim<dim>::LeafIterator vEndIt = uniformGrid[i].leafend<dim>();
-            
-            for (; vIt!=vEndIt; ++vIt)
-                uniformGrid[i].setPosition(vIt, vIt->geometry()[0] + uniformX[i][uniformGrid[i].leafIndexSet().index(*vIt)]);
-                
-        }
+        for (int i=0; i<2; i++)
+            deformation[i]->setDeformation(uniformX[i]);
 
-        LeafAmiraMeshWriter<GridType>::writeGrid(uniformGrid[0], "0uniformGridDef");
-        LeafAmiraMeshWriter<GridType>::writeGrid(uniformGrid[1], "1uniformGridDef");
+        LeafAmiraMeshWriter<DeformedGridType>::writeGrid(*defUniformGrid[0], "0uniformGridDef");
+        LeafAmiraMeshWriter<DeformedGridType>::writeGrid(*defUniformGrid[1], "1uniformGridDef");
 
     }
     // /////////////////////////////////////////////////////////
@@ -707,13 +726,24 @@ int main(int argc, char *argv[]) try
         AmiraMeshReader<GridType>::read(adaptiveGrid[1], path + object1Name);
     }
 
+    // create the deformed grids
+    array<DeformedGridType*,2> defAdaptiveGrid;
+    for (int i=0;i<2;i++) {
+        VectorType displace(adaptiveGrid[i].size(dim));
+        displace = 0;
+
+        deformation[i]->setDeformation(displace);
+        deformation[i]->setGridView(adaptiveGrid[i].leafGridView());
+        defAdaptiveGrid[i]= new DeformedGridType(&adaptiveGrid[i],deformation[i]);
+    }
+
     // /////////////////////////////////////////////////////
     // read field describing the obstacles again, so this time
     // it is attached to the adaptive grid.
     // /////////////////////////////////////////////////////
 
-    obsPatch.setup(adaptiveGrid[0].levelView(0));
-    readBoundaryPatch<GridType>(obsPatch, path + obsFilename);
+    obsPatch.setup(defAdaptiveGrid[0]->levelGridView(0));
+    readBoundaryPatch<DeformedGridType>(obsPatch, path + obsFilename);
 
 
     // /////////////////////////////////////////////////////////
@@ -730,14 +760,13 @@ int main(int argc, char *argv[]) try
         //   Solve problem on locally refined grid
         // /////////////////////////////////////////////////////////////
         if (i!=0)
-            compute<2>(adaptiveGrid, adaptiveX, dirichletValues, dirichletNodes, coarseDirichletBoundary,
+            compute<2>(defAdaptiveGrid, deformation, adaptiveX, dirichletValues, dirichletNodes, coarseDirichletBoundary,
                        obsPatch, obsDistance, trustRegionRadius, maxTrustRegionSteps, tolerance, numIt);
         
         // /////////////////////////////////////////////////////////////
         //   Compute 'actual' error by comparing the the solution
         //   on the uniformly refined grid
         // /////////////////////////////////////////////////////////////
-        array<LinearElasticityLocalStiffness<GridType::LeafGridView,double>,2> linearLocalStiffness;
         array<StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement>*,2> linearLocalStiffness;
         linearLocalStiffness[0]=new StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement> (E[0],nu[0]);
         linearLocalStiffness[1]=new StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement> (E[1],nu[1]);
@@ -759,7 +788,8 @@ int main(int argc, char *argv[]) try
         // ////////////////////////////////////////////////////
         HierarchicContactEstimator<GridType,field_type> estimator(adaptiveGrid[0], adaptiveGrid[1]);
 
-        estimator.couplings_[0].obsPatch_    = &obsPatch;
+        LevelBoundary undefObsPatch(adaptiveGrid[0].levelGridView(0),*obsPatch.getVertices());
+        estimator.couplings_[0].obsPatch_    = &undefObsPatch;
         estimator.couplings_[0].obsDistance_ = obsDistance;
         estimator.couplings_[0].type_        = CouplingPairBase::CONTACT;
 
@@ -781,8 +811,8 @@ int main(int argc, char *argv[]) try
         }
         BoundaryPatchProlongator<GridType>::prolong(coarseDirichletBoundary[0], leafDirichletBoundary[0]);
         BoundaryPatchProlongator<GridType>::prolong(coarseDirichletBoundary[1], leafDirichletBoundary[1]);
-
-        array<OgdenMaterialLocalStiffness<GridType::LeafGridView,double>,2> ogdenLocalStiffness;
+        typedef P2NodalBasis<GridType::LeafGridView> P2Basis;
+        array<OgdenMaterialLocalStiffness<GridType,typename P2Basis::LocalFiniteElement, typename P2Basis::LocalFiniteElement>,2> ogdenLocalStiffness;
         ogdenLocalStiffness[0].setEandNu(E[0], nu[0], d[0]);
         ogdenLocalStiffness[1].setEandNu(E[1], nu[1], d[1]);
 
@@ -834,6 +864,10 @@ int main(int argc, char *argv[]) try
             if (paramBoundaries)
                 improveGrid(adaptiveGrid[j],10);
 
+            deformation[j]->setDeformation(adaptiveX[j]);
+            deformation[j]->setGridView(adaptiveGrid[j].leafGridView());
+            defAdaptiveGrid[j]->update();
+
         }
 
         std::cout << "########################################################" << std::endl;
@@ -861,18 +895,11 @@ int main(int argc, char *argv[]) try
 
     if (dim==2) {
         
-        for (int i=0; i<2; i++) {
-
-            GridType::Codim<dim>::LeafIterator vIt    = adaptiveGrid[i].leafbegin<dim>();
-            GridType::Codim<dim>::LeafIterator vEndIt = adaptiveGrid[i].leafend<dim>();
-            
-            for (; vIt!=vEndIt; ++vIt)
-                adaptiveGrid[i].setPosition(vIt, vIt->geometry()[0] + adaptiveX[i][adaptiveGrid[i].leafIndexSet().index(*vIt)]);
-                
-        }
+        for (int i=0; i<2; i++)
+            deformation[i]->setDeformation(adaptiveX[i]);
 
-        LeafAmiraMeshWriter<GridType>::writeGrid(adaptiveGrid[0], "0resultGridDef");
-        LeafAmiraMeshWriter<GridType>::writeGrid(adaptiveGrid[1], "1resultGridDef");
+        LeafAmiraMeshWriter<DeformedGridType>::writeGrid(*defAdaptiveGrid[0], "0resultGridDef");
+        LeafAmiraMeshWriter<DeformedGridType>::writeGrid(*defAdaptiveGrid[1], "1resultGridDef");
 
     }