Fixes parallel assembly:

- read in the deformation file on each process
- modify trustregionsolver so it can run in parallel

dune/gfe/localgeodesicfeadolcstiffness.hh

@@ -106,13 +106,13 @@ energy(const typename Basis::LocalView& localView,
     try {
         energy = localEnergy_->energy(localView,localASolution);
     } catch (Dune::Exception &e) {
-        trace_off(rank);
+        trace_off();
         throw e;
     }
 
     energy >>= pureEnergy;
 
-    trace_off(rank);
+    trace_off();
 #if 0
     size_t tape_stats[STAT_SIZE];
     tapestats(rank,tape_stats);             // reading of tape statistics

dune/gfe/parallel/globalmapper.hh

@@ -28,12 +28,11 @@ namespace Dune {
       Master m = Dune::ParMG::entityMasterRank(basis.gridView(), [&](int, int codim) -> bool {
               return dofmap.codimSet().test(codim);
             });
-      DofMaster dofmaster = Dune::ParMG::dofMaster(basis, m);
 
-      coarseGlobalDof_ = globalDof(basis,m, dofmap);
+      globalDof_ = globalDof(basis,m, dofmap);
 
       // total number of degrees of freedom
-      size_ = coarseGlobalDof_.size;
+      size_ = globalDof_.size;
     }
 #else
     GlobalMapper(const Basis& basis)
@@ -46,7 +45,7 @@ namespace Dune {
     /** \brief Given a local index, retrieve its index globally unique over all processes. */
     Index index(const int& localIndex) const {
 #if HAVE_DUNE_PARMG
-      return coarseGlobalDof_.globalDof[localIndex];
+      return globalDof_.globalDof[localIndex];
 #else
       return localIndex;
 #endif
@@ -58,7 +57,7 @@ namespace Dune {
     }
 
 #if HAVE_DUNE_PARMG
-    ParMG::GlobalDof coarseGlobalDof_;
+    ParMG::GlobalDof globalDof_;
 #endif
 
     std::size_t size_;

dune/gfe/riemanniantrsolver.cc

@@ -418,8 +418,10 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
                 if ((*mgStep->ignoreNodes_)[j][k])  // global Dirichlet nodes set
                   gradient[j][k] = 0;
 
-            if (this->verbosity_ == Solver::FULL and rank==0)
-              std::cout << "Gradient norm: " << l2Norm_->operator()(gradient) << std::endl;
+            if (this->verbosity_ == Solver::FULL and rank==0) {
+              std::cout << "Gradient operator norm: " << l2Norm_->operator()(gradient) << std::endl;
+              std::cout << "Gradient norm: " << gradient.two_norm() << std::endl;
+            }
 
             if (this->verbosity_ == Solver::FULL)
               std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
@@ -460,15 +462,14 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
             } catch (Dune::Exception &e) {
                 std::cerr << "Error while solving: " << e << std::endl;
                 solved = false;
-                corr_global = 0;
             }
             std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
             totalSolverTime += solutionTimer.elapsed();
 
-            if (mgStep && solved)
+            if (mgStep && solved) {
                 corr_global = mgStep->getSol();
-
-            //std::cout << "Correction: " << std::endl << corr_global << std::endl;
+                std::cout << "Two norm of the correction: " << corr_global.two_norm() << std::endl;
+            }
         }
 
         // Distribute solution
@@ -476,7 +477,13 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
             std::cout << "Transfer solution back to root process ..." << std::endl;
 
 #if HAVE_MPI
-        corr = vectorComm.scatter(corr_global);
+        solved = grid_->comm().min(solved);
+        if (solved) {
+            corr = vectorComm.scatter(corr_global);
+        } else  {
+            corr_global = 0;
+            corr = 0;
+        }
 #else
         corr = corr_global;
 #endif
@@ -582,11 +589,14 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
             } catch (Dune::Exception &e) {
                 std::cerr << "Error while computing the energy of the new Iterate: " << e << std::endl;
                 std::cerr << "Redoing trust region step with smaller radius..." << std::endl;
-                newIterate = x_;
                 solved = false;
-                energy = oldEnergy;
             }
-            if (solved) {
+            solved = grid_->comm().min(solved);
+
+            if (!solved) {
+                newIterate = x_;
+                energy = oldEnergy;
+            } else {
                 energy = grid_->comm().sum(energy);
 
                 // compute the model decrease

src/film-on-substrate.cc

@@ -68,9 +68,6 @@
 #include <dune/solvers/solvers/iterativesolver.hh>
 #include <dune/solvers/norms/energynorm.hh>
 
-#include <iostream>
-#include <fstream>
-
 // grid dimension
 #ifndef WORLD_DIM
 #  define WORLD_DIM 3
@@ -217,11 +214,11 @@ int main (int argc, char *argv[]) try
 
     grid->adapt();
 
-    grid->loadBalance();
-
     numLevels--;
   }
 
+  grid->loadBalance();
+
   if (mpiHelper.rank()==0)
     std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl;
 
@@ -262,10 +259,8 @@ int main (int argc, char *argv[]) try
   auto neumannBoundary = std::make_shared<BoundaryPatch<GridView>>(gridView, neumannVertices);
   BoundaryPatch<GridView> surfaceShellBoundary(gridView, surfaceShellVertices);
 
-  if (mpiHelper.rank()==0) {
-    std::cout << "Neumann boundary has " << neumannBoundary->numFaces() << " faces\n";
-    std::cout << "Shell boundary has " << surfaceShellBoundary.numFaces() << " faces\n";
-  }
+  std::cout << "On rank " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary->numFaces() << " faces\n";
+  std::cout << "On rank " << mpiHelper.rank() << ": Shell boundary has " << surfaceShellBoundary.numFaces() << " faces\n";
 
 
   BitSetVector<1> dirichletNodes(feBasis.size(), false);
@@ -364,33 +359,50 @@ int main (int argc, char *argv[]) try
   // Read in the grid deformation
   if (startFromFile) {
     const std::string pathToGridDeformationFile = parameterSet.get("pathToGridDeformationFile", "");
-
     // for this, we create a basis of order 1 in order to deform the geometries on the surface shell boundary
-    typedef Dune::Functions::LagrangeBasis<GridView, 1> FEBasisOrder1;
-    FEBasisOrder1 feBasisOrder1(gridView);
-  
+    auto feBasisOrder1 = makeBasis(
+      gridView,
+      power<dim>(
+        lagrange<1>(),
+        blockedInterleaved()
+    ));
+    GlobalIndexSet<GridView> globalVertexIndexSet(gridView,dim);
+
+    std::unordered_map<std::string, FieldVector<double,3>> deformationMap;
+    std::string line, displacement, entry;
+    if (mpiHelper.rank() == 0) 
+      std::cout << "Reading in deformation file: " << pathToGridDeformationFile + parameterSet.get<std::string>("gridDeformationFile") << std::endl;
     // Read grid deformation information from the file specified in the parameter set via gridDeformationFile
-    BlockVector<FieldVector<double,3> > gridDeformationFromFile(feBasisOrder1.size());
-    std::string line;
-    std::ifstream file(pathToGridDeformationFile + parameterSet.get<std::string>("gridDeformationFile"));
+    std::ifstream file(pathToGridDeformationFile + parameterSet.get<std::string>("gridDeformationFile"), std::ios::in);
     if (file.is_open()) {
-      size_t i = 0;
       while (std::getline(file, line)) {
         size_t j = 0;
-        std::stringstream entries(line);
-        std::string entry;
-        FieldVector<double,3> coord(0);
-        while(entries >> entry) {
-          coord[j++] = std::stod(entry);
-        }
-        gridDeformationFromFile[i++] = coord;
+        size_t pos = line.find(":");
+        displacement = line.substr(pos + 1);
+        line.erase(pos);
+        std::stringstream entries(displacement);
+        FieldVector<double,3> displacementVector(0);
+        while(entries >> entry)
+          displacementVector[j++] = std::stod(entry);
+        deformationMap.insert({line,displacementVector});
       }
-      if (i != feBasisOrder1.size())
+      if (mpiHelper.rank() == 0)
+        std::cout << "... done: The grid has " << globalVertexIndexSet.size(dim) << " vertices and the defomation file has " << deformationMap.size() << " entries." << std::endl;
+      if (deformationMap.size() != globalVertexIndexSet.size(dim))
         DUNE_THROW(Exception, "Error: Grid and deformation vector do not match!");
       file.close();
     } else {
       DUNE_THROW(Exception, "Error: Could not open the file containing the deformation vector!");
     }
+  
+    BlockVector<FieldVector<double,dim>> gridDeformationFromFile;
+    Dune::Functions::interpolate(feBasisOrder1, gridDeformationFromFile, [](FieldVector<double,dim> x){ return x; });
+
+    for (auto& entry : gridDeformationFromFile) {
+      std::stringstream stream;
+      stream << entry;
+      entry = deformationMap.at(stream.str()); //Look up the deformation for this vertex in the deformationMap
+    }
 
     auto gridDeformationFromFileFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasisOrder1, gridDeformationFromFile);
     auto localGridDeformationFromFileFunction = localFunction(gridDeformationFromFileFunction);