From 41d79203f561d9072901fd54e88e90e62ec163e0 Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de>
Date: Wed, 17 Jun 2020 22:57:57 +0200
Subject: [PATCH] Read in deformation function on each process - parallel
assembly works again!
---
dune/gfe/parallel/globalmapper.hh | 9 +++---
dune/gfe/riemanniantrsolver.cc | 6 ++--
src/film-on-substrate.cc | 50 ++++++++++++++++++++-----------
3 files changed, 40 insertions(+), 25 deletions(-)
diff --git a/dune/gfe/parallel/globalmapper.hh b/dune/gfe/parallel/globalmapper.hh
index 554a70fd..4cf617dd 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 e70e6649..1b4431a7 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/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;
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 471ea478..19f7b967 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -217,11 +217,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;
@@ -364,33 +364,47 @@ 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);
+
// 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::unordered_map<std::string, FieldVector<double,3>> deformationMap;
+ std::string line, displacement, entry;
+ //TODO: Read this using MPI_File_open and MPI_File_read
+ std::ifstream file(pathToGridDeformationFile + parameterSet.get<std::string>("gridDeformationFile") + std::to_string(mpiHelper.rank()));
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 (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);
--
GitLab