diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 5e3e54231ce4383fb6700c1fa16800cfd2128c8b..a352da55283f59f131df1d88a7c04bf16e820273 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -81,20 +81,6 @@
size_t const dims = MY_DIM;
-template <class GridView>
-void writeToVTK(const GridView& gridView, const std::string path, const std::string name) {
- Dune::VTKWriter<GridView> vtkwriter(gridView);
-
- //std::ofstream lStream( "garbage.txt" );
- std::streambuf* lBufferOld = std::cout.rdbuf();
- //std::cout.rdbuf(lStream.rdbuf());
-
- vtkwriter.pwrite(name, path, path);
-
- std::cout.rdbuf( lBufferOld );
-}
-
-
Dune::ParameterTree getParameters(int argc, char *argv[]) {
Dune::ParameterTree parset;
Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem.cfg", parset);
@@ -160,7 +146,9 @@ int main(int argc, char *argv[]) {
}
#endif
+ // ------------
// set up grids
+ // ------------
GridsConstructor<Grid> gridConstructor(cuboidGeometries);
auto& grids = gridConstructor.getGrids();
@@ -174,8 +162,6 @@ int main(int argc, char *argv[]) {
weakPatch = getWeakPatch<LocalVector>(parset.sub("boundary.friction.weakening"), cuboidGeometry);
refine(*grids[i], weakPatch, parset.get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
- writeToVTK(grids[i]->leafGridView(), "", "grid" + std::to_string(i));
-
// determine minDiameter and maxDiameter
double minDiameter = std::numeric_limits<double>::infinity();
double maxDiameter = 0.0;
@@ -197,18 +183,9 @@ int main(int argc, char *argv[]) {
const auto deformedGrids = deformedGridComplex.gridVector();
- /*
- // test deformed grids
- for (size_t i=0; i<deformedGrids.size(); i++) {
- Vector def(deformedGrids[i]->size(dims));
- def = 1;
- deformedGridComplex.setDeformation(def, i);
-
- writeToVTK(deformedGrids[i]->leafGridView(), "", "deformedGrid" + std::to_string(i));
- }
- // test deformed grids
- */
-
+ // ------------------------
+ // set up contact couplings
+ // ------------------------
std::vector<std::shared_ptr<DefLeafGridView>> leafViews(deformedGrids.size());
std::vector<std::shared_ptr<DefLevelGridView>> levelViews(deformedGrids.size());
std::vector<size_t> leafVertexCounts(deformedGrids.size());
@@ -247,12 +224,9 @@ int main(int argc, char *argv[]) {
couplings[i]->set(i, i+1, nonmortarPatch, mortarPatch, 0.1, Dune::Contact::CouplingPairBase::STICK_SLIP, contactProjection, backend);
}
+ // ----------------------------------
// set up dune-contact nBodyAssembler
- /*std::vector<const DeformedGrid*> const_grids(bodyCount);
- for (size_t i=0; i<const_grids.size(); i++) {
- const_grids[i] = &deformedGrids.grid("body" + std::to_string(i));
- }*/
-
+ // ----------------------------------
Dune::Contact::NBodyAssembler<DeformedGrid, Vector> nBodyAssembler(bodyCount, bodyCount-1);
nBodyAssembler.setGrids(deformedGrids);
@@ -263,7 +237,9 @@ int main(int argc, char *argv[]) {
nBodyAssembler.assembleTransferOperator();
nBodyAssembler.assembleObstacle();
+ // --------------------------
// set up boundary conditions
+ // --------------------------
std::vector<BoundaryPatch<DefLeafGridView>> neumannBoundaries(bodyCount);
std::vector<BoundaryPatch<DefLeafGridView>> surfaces(bodyCount);
@@ -273,7 +249,7 @@ int main(int argc, char *argv[]) {
// Dirichlet boundary
std::vector<Dune::BitSetVector<dims>> noNodes(bodyCount);
- std::vector<Dune::BitSetVector<dims>> dirichletNodes(bodyCount);
+ std::vector<std::vector<Dune::BitSetVector<dims>>> dirichletNodes(bodyCount);
// set up functions for time-dependent boundary conditions
using Function = Dune::VirtualFunction<double, double>;
@@ -306,8 +282,10 @@ int main(int argc, char *argv[]) {
velocityDirichletFunctions[i] = new VelocityDirichletCondition();
noNodes[i] = Dune::BitSetVector<dims>(leafVertexCount);
- dirichletNodes[i] = Dune::BitSetVector<dims>(leafVertexCount);
+
auto& gridDirichletNodes = dirichletNodes[i];
+ gridDirichletNodes = Dune::BitSetVector<dims>(leafVertexCount);
+
for (int j=0; j<leafVertexCount; j++) {
if (leafFaces[i]->right.containsVertex(j))
gridDirichletNodes[j][0] = true;
@@ -322,7 +300,42 @@ int main(int argc, char *argv[]) {
}
}
+ // extend Dirichlet information to all grid levels
+ /* BoundaryPatchProlongator<GridType>::prolong(dirichletBoundary[i]);
+ dirichletNodes[i].resize(toplevel+1);
+
+ for (int j=0; j<=toplevel; j++) {
+
+ int fSSize = grids[i]->size(j,dim);
+ dirichletNodes[i][j].resize(fSSize);
+ for (int k=0; k<fSSize; k++)
+ for (int l=0; l<dim; l++)
+ dirichletNodes[i][j][k][l] = dirichletBoundary[i][j].containsVertex(k);
+
+ }
+
+ sampleOnBitField(*grids[i], dirichletValues[i], dirichletNodes[i]);*/
+
+ const int toplevel = deformedGrids[0]->maxLevel();
+ std::vector<Dune::BitSetVector<dims> > totalDirichletNodes(toplevel+1);
+
+ for (int i=0; i<=toplevel; i++) {
+ int totalSize = 0;
+ for (int j=0; j<deformedGrids.size(); j++)
+ totalSize += deformedGrids[j]->size(i, dims);
+
+ totalDirichletNodes[i].resize(totalSize);
+
+ int idx=0;
+ for (int j=0; j<deformedGrids.size(); j++)
+ for (size_t k=0; k<dirichletNodes[j][i].size(); k++)
+ for (int l=0; l<dims; l++)
+ totalDirichletNodes[i][idx++][l] = dirichletNodes[j][i][k][l];
+ }
+
+ // ------------------------------------------------------------------------------------
// set up individual assemblers for each body, assemble problem (matrices, forces, rhs)
+ // ------------------------------------------------------------------------------------
std::vector<std::shared_ptr<Assembler>> assemblers(bodyCount);
Matrices<Matrix> matrices(bodyCount);
@@ -355,40 +368,9 @@ int main(int argc, char *argv[]) {
};
}
- /* Jonny 2bcontact
- // make dirichlet bitfields containing dirichlet information for both grids
- int size = rhs[0].size() + rhs[1].size();
-
- Dune::BitSetVector<dims> totalDirichletNodes(size);
-
- for (size_t i=0; i<dirichletNodes[0].size(); i++)
- for (int j=0; j<dims; j++)
- totalDirichletNodes[i][j] = dirichletNodes[0][i][j];
-
- int offset = rhs[0].size();
- for (size_t i=0; i<dirichletNodes[1].size(); i++)
- for (int j=0; j<dims; j++)
- totalDirichletNodes[offset + i][j] = dirichletNodes[1][i][j];
-
- // assemble separate linear elasticity problems
- std::array<MatrixType,2> stiffnessMatrix;
- std::array<const MatrixType*, 2> submat;
-
- for (size_t i=0; i<2; i++) {
- OperatorAssembler<P1Basis,P1Basis> globalAssembler(*p1Basis[i],*p1Basis[i]);
-
- double s = (i==0) ? E0 : E1;
- StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement> localAssembler(s, nu);
-
- globalAssembler.assemble(localAssembler, stiffnessMatrix[i]);
-
- submat[i] = &stiffnessMatrix[i];
- }
-
- MatrixType bilinearForm;
- contactAssembler.assembleJacobian(submat, bilinearForm);
- */
-
+ // -----------------
+ // init input/output
+ // -----------------
using MyProgramState = ProgramState<Vector, ScalarVector>;
MyProgramState programState(leafVertexCounts);
auto const firstRestart = parset.get<size_t>("io.restarts.first");
@@ -418,11 +400,10 @@ int main(int argc, char *argv[]) {
restartIO->read(firstRestart, programState);
else
programState.setupInitialConditions(parset, nBodyAssembler, externalForces,
- matrices, assemblers, dirichletNodes,
+ matrices, assemblers, totalDirichletNodes,
noNodes, frictionalBoundaries, body);
-
// assemble friction
std::vector<std::shared_ptr<MyGlobalFrictionData<LocalVector>>> frictionInfo(weakPatches.size());
std::vector<std::shared_ptr<GlobalFriction<Matrix, Vector>>> globalFriction(weakPatches.size());
@@ -432,7 +413,6 @@ int main(int argc, char *argv[]) {
globalFriction[i]->updateAlpha(programState.alpha[i]);
}
-
using MyVertexBasis = typename Assembler::VertexBasis;
using MyCellBasis = typename Assembler::CellBasis;
std::vector<Vector> vertexCoordinates(leafVertexCounts.size());
@@ -499,9 +479,11 @@ int main(int argc, char *argv[]) {
};
report(true);
+ // -------------------
// Set up TNNMG solver
+ // -------------------
using NonlinearFactory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
- NonlinearFactory factory(parset.sub("solver.tnnmg"), nBodyAssembler, dirichletNodes);
+ NonlinearFactory factory(parset.sub("solver.tnnmg"), nBodyAssembler, totalDirichletNodes);
using MyUpdater = Updaters<RateUpdater<Vector, Matrix, Function, dims>,
StateUpdater<ScalarVector, Vector>>;
@@ -518,8 +500,7 @@ int main(int argc, char *argv[]) {
L, V0));
- auto const refinementTolerance =
- parset.get<double>("timeSteps.refinementTolerance");
+ auto const refinementTolerance = parset.get<double>("timeSteps.refinementTolerance");
auto const mustRefine = [&](MyUpdater &coarseUpdater,
MyUpdater &fineUpdater) {
std::vector<ScalarVector> coarseAlpha;
@@ -540,11 +521,11 @@ int main(int argc, char *argv[]) {
std::signal(SIGINT, handleSignal);
std::signal(SIGTERM, handleSignal);
- AdaptiveTimeStepper<NonlinearFactory, MyUpdater,
- EnergyNorm<ScalarMatrix, ScalarVector>>
+ AdaptiveTimeStepper<NonlinearFactory, MyUpdater, EnergyNorm<ScalarMatrix, ScalarVector>>
adaptiveTimeStepper(factory, parset, globalFriction, current,
programState.relativeTime, programState.relativeTau,
externalForces, stateEnergyNorms, mustRefine);
+
while (!adaptiveTimeStepper.reachedEnd()) {
programState.timeStep++;
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index f9220a826a7e5e9a5aae422ea8def92933421458..726592df93cd52a2356843051ecee65f5fb33eb8 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -16,7 +16,7 @@
template <class DeformedGrid, class Nonlinearity, class Matrix, class Vector>
SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::SolverFactory(
const Dune::ParameterTree& parset, const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
- const std::vector<Dune::BitSetVector<DeformedGrid::dimension>>& ignoreNodes)
+ const std::vector<std::vector<Dune::BitSetVector<DeformedGrid::dimension>>>& ignoreNodes)
: nBodyAssembler_(nBodyAssembler),
baseEnergyNorm_(baseSolverStep_),
baseSolver_(&baseSolverStep_,
@@ -28,11 +28,11 @@ SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::SolverFactory(
multigridStep_ = std::make_shared<Step>();
// tnnmg iteration step
multigridStep_->setMGType(parset.get<int>("main.multi"), parset.get<int>("main.pre"), parset.get<int>("main.post"));
- //multigridStep_->setIgnore(ignoreNodes);
+ multigridStep_->setIgnore(ignoreNodes);
multigridStep_->setBaseSolver(baseSolver_);
multigridStep_->setSmoother(&presmoother_, &postsmoother_);
- // multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
- // multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
+ multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
+ multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
// create the transfer operators
const int nCouplings = nBodyAssembler_.nCouplings();