diff --git a/dune/gfe/localgeodesicfeadolcstiffness.hh b/dune/gfe/localgeodesicfeadolcstiffness.hh
index 148f2c4ca79ae435782d3c4a690832a497a184cd..f7c4d27c629fd9b6e61344f0e524bf039a05bd59 100644
--- a/dune/gfe/localgeodesicfeadolcstiffness.hh
+++ b/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
diff --git a/dune/gfe/parallel/globalmapper.hh b/dune/gfe/parallel/globalmapper.hh
index 554a70fd67c42d9ea71fabb4e9267738be5bbab1..4cf617dd8b219ad5de4ac32ab927e603d6340e8d 100644
--- a/dune/gfe/parallel/globalmapper.hh
+++ b/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_;
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index e70e6649f60d03a8e3e69ce99285c9812fbaad00..b8604dfc6c94f033342cbc61addc33454fd577fa 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -403,6 +403,7 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
*hessianMatrix_,
i==0 // assemble occupation pattern only for the first call
);
+ std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
rhs *= -1; // The right hand side is the _negative_ gradient
@@ -418,11 +419,12 @@ 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 operator norm: " << l2Norm_->operator()(gradient) << std::endl;
+ std::cout << "Gradient norm: " << gradient.two_norm() << std::endl;
+ }
if (this->verbosity_ == Solver::FULL and rank==0)
- std::cout << "Gradient norm: " << l2Norm_->operator()(gradient) << std::endl;
-
- if (this->verbosity_ == Solver::FULL)
- std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
+ std::cout << "Oveall assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
totalAssemblyTime += gradientTimer.elapsed();
// Transfer matrix data
@@ -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
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 471ea478564f912ead8f141f3acc0c03caf597d9..3d3fa4d8def6a22dd602b2bcc046a767a8f30463 100644
--- a/src/film-on-substrate.cc
+++ b/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);