diff --git a/src/foam/foam.cc b/src/foam/foam.cc
index 324bcd88ad6e55762bad4cbc90b3e8caecfa6a08..02e51281ecb240edbe844c50cce6f6a877f46ad1 100644
--- a/src/foam/foam.cc
+++ b/src/foam/foam.cc
@@ -32,6 +32,7 @@
#include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/fufem/formatstring.hh>
#include <dune/fufem/hdf5/file.hh>
+#include <dune/fufem/assemblers/transferoperatorassembler.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
@@ -41,6 +42,8 @@
#include <dune/contact/common/couplingpair.hh>
#include <dune/contact/projections/normalprojection.hh>
+#include <dune/tnnmg/functionals/quadraticfunctional.hh>
+
#include <dune/tectonic/assemblers.hh>
#include <dune/tectonic/gridselector.hh>
@@ -196,173 +199,173 @@ int main(int argc, char *argv[]) {
programState.setupInitialConditions(parset, contactNetwork);
}
- /*
- // setup initial conditions in program state
- programState.alpha[0] = 0;
- programState.alpha[1] = parset.get<double>("boundary.friction.initialAlpha"); */
-
- /*auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
- frictionalBoundary, frictionData);
- myGlobalNonlinearity->updateLogState(alpha_initial);
-*/
- /*
- // Solving a linear problem with a multigrid solver
- auto linearSolver = makeLinearSolver<ContactNetwork, Vector>(parset, contactNetwork);
const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+ /* auto linearSolver = makeLinearSolver<ContactNetwork, Vector>(parset, contactNetwork);
+ auto nonlinearity = ZeroNonlinearity();
+ // Solving a linear problem with a multigrid solver
auto const solveLinearProblem = [&](
const BitVector& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
const std::vector<Vector>& _rhs, std::vector<Vector>& _x,
Dune::ParameterTree const &_localParset) {
- std::vector<const Matrix*> matrices_ptr(_matrices.size());
- for (size_t i=0; i<matrices_ptr.size(); i++) {
- matrices_ptr[i] = _matrices[i].get();
- }
+ Vector totalX;
+ nBodyAssembler.nodalToTransformed(_x, totalX);
- // assemble full global contact problem
- Matrix bilinearForm;
- nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
-
- Vector totalRhs, oldTotalRhs;
- nBodyAssembler.assembleRightHandSide(_rhs, totalRhs);
- oldTotalRhs = totalRhs;
-
- Vector totalX, oldTotalX;
- nBodyAssembler.nodalToTransformed(_x, totalX);
- oldTotalX = totalX;
-
- // get lower and upper obstacles
- const auto& totalObstacles = nBodyAssembler.totalObstacles_;
- Vector lower(totalObstacles.size());
- Vector upper(totalObstacles.size());
-
- for (size_t j=0; j<totalObstacles.size(); ++j) {
- const auto& totalObstaclesj = totalObstacles[j];
- auto& lowerj = lower[j];
- auto& upperj = upper[j];
- for (size_t d=0; d<dims; ++d) {
- lowerj[d] = totalObstaclesj[d][0];
- upperj[d] = totalObstaclesj[d][1];
- }
- }
+ FunctionalFactory<std::decay_t<decltype(nBodyAssembler)>, decltype(nonlinearity), Matrix, Vector> functionalFactory(nBodyAssembler, nonlinearity, _matrices, _rhs);
+ functionalFactory.build();
+ auto functional = functionalFactory.functional();
- // set up functional
- using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
- Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper);
+ NonlinearSolver<std::remove_reference_t<decltype(*functional)>, BitVector> solver(parset.sub("solver.tnnmg"), linearSolver, functional, _dirichletNodes);
+ solver.solve(_localParset, totalX);
- // set up TNMMG solver
- using Factory = SolverFactory<Functional, BitVector>;
- //Factory factory(parset.sub("solver.tnnmg"), J, linearSolver, _dirichletNodes);
- Factory factory(parset.sub("solver.tnnmg"), J, _dirichletNodes);
- factory.build(linearSolver); */
+ nBodyAssembler.postprocess(totalX, _x);
+ };
- /* std::vector<BitVector> bodyDirichletNodes;
- nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
- for (size_t i=0; i<bodyDirichletNodes.size(); i++) {
- print(bodyDirichletNodes[i], "bodyDirichletNodes_" + std::to_string(i) + ": ");
- }*/
+ programState.timeStep = parset.get<size_t>("initialTime.timeStep");
+ programState.relativeTime = parset.get<double>("initialTime.relativeTime");
+ programState.relativeTau = parset.get<double>("initialTime.relativeTau");
- /*
- auto tnnmgStep = factory.step();
- factory.setProblem(totalX);
+ std::vector<Vector> ell0(bodyCount);
+ for (size_t i=0; i<bodyCount; i++) {
+ programState.u[i] = 0.0;
+ programState.v[i] = 0.0;
+ programState.a[i] = 0.0;
- const EnergyNorm<Matrix, Vector> norm(bilinearForm);
+ ell0[i].resize(programState.u[i].size());
+ ell0[i] = 0.0;
- LoopSolver<Vector> solver(
- *tnnmgStep.get(), _localParset.get<size_t>("maximumIterations"),
- _localParset.get<double>("tolerance"), norm,
- _localParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
+ contactNetwork.body(i)->externalForce()(programState.relativeTime, ell0[i]);
+ }
+ */
+ // Initial displacement: Start from a situation of minimal stress,
+ // which is automatically attained in the case [v = 0 = a].
+ // Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
+ BitVector totalDirichletNodes;
+ contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
+ size_t dof = 0;
+ for (size_t i=0; i<bodyCount; i++) {
+ const auto& body = *contactNetwork.body(i);
+
+ if (body.data()->getYoungModulus() == 0.0) {
+ for (; dof<body.nVertices(); dof++) {
+ totalDirichletNodes[dof] = true;
+ }
+ } else {
+ dof += body.nVertices();
+ }
+ }
- solver.preprocess();
- solver.solve();
+ using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
+ BoundaryNodes dirichletNodes;
+ contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
- nBodyAssembler.postprocess(tnnmgStep->getSol(), _x);
- }; */
+ std::vector<const Dune::BitSetVector<1>*> frictionNodes;
+ contactNetwork.frictionNodes(frictionNodes);
- /*using LinearFactory = SolverFactory<
- dims, BlockNonlinearTNNMGProblem<ConvexProblem<
- ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
- Grid>;
- ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;*/
+ /*
+ std::cout << "solving linear problem for u..." << std::endl;
+
+ {
+ int nVertices = contactNetwork.body(1)->nVertices();
+ BitVector uDirichletNodes(nVertices, false);
+
+ int idx=0;
+ const auto& body = contactNetwork.body(1);
+ using LeafBody = typename ContactNetwork::LeafBody;
+ std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> boundaryConditions;
+ body->boundaryConditions("dirichlet", boundaryConditions);
+
+ if (boundaryConditions.size()>0) {
+ const int idxBackup = idx;
+ for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
+ const auto& nodes = boundaryConditions[bc]->boundaryNodes();
+ for (size_t j=0; j<nodes->size(); j++, idx++)
+ for (int k=0; k<dims; k++)
+ uDirichletNodes[idx][k] = uDirichletNodes[idx][k] or (*nodes)[j][k];
+ idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
+ }
+ }
- /*
- // TODO: clean up once generic lambdas arrive
- auto const solveLinearProblem = [&](
- Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
- Vector const &_rhs, Vector &_x, EnergyNorm<Matrix, Vector> _norm,
- Dune::ParameterTree const &_localParset) {
+ for (auto i=0; i<nVertices; i++) {
+ if ((*frictionNodes[1])[i][0]) {
+ uDirichletNodes[i][1] = true;
+ }
+ }
+ using Norm = EnergyNorm<Matrix, Vector>;
+ using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+ // transfer operators need to be recomputed on change due to setDeformation()
+ using TransferOperator = CompressedMultigridTransfer<Vector>;
+ using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+ const auto& grid = contactNetwork.body(1)->grid();
+ TransferOperators transfer(grid->maxLevel());
+ for (size_t i = 0; i < transfer.size(); ++i)
+ {
+ // create transfer operators from level i to i+1 (note that this will only work for either uniformly refined grids or adaptive grids with RefinementType=COPY)
+ auto t = std::make_shared<TransferOperator>();
+ t->setup(*grid, i, i+1);
+ transfer[i] = t;
+ }
+ auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >(
+ *contactNetwork.matrices().elasticity[1],
+ programState.u[1],
+ ell0[1]);
+ linearMultigridStep->setIgnore(uDirichletNodes);
+ linearMultigridStep->setMGType(1, 3, 3);
+ linearMultigridStep->setSmoother(TruncatedBlockGSStep<Matrix, Vector>());
+ linearMultigridStep->setTransferOperators(transfer);
- typename LinearFactory::ConvexProblem convexProblem(
- 1.0, _matrix, zeroNonlinearity, _rhs, _x);
- typename LinearFactory::BlockProblem problem(parset, convexProblem);
+ const auto& solverParset = parset.sub("u0.solver");
+ Dune::Solvers::LoopSolver<Vector> solver(linearMultigridStep, solverParset.get<size_t>("maximumIterations"),
+ solverParset.get<double>("tolerance"), Norm(*linearMultigridStep),
+ solverParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
- auto multigridStep = factory.getSolver();
- multigridStep->setProblem(_x, problem);
- LoopSolver<Vector> solver(
- multigridStep, _localParset.get<size_t>("maximumIterations"),
- _localParset.get<double>("tolerance"), &_norm,
- _localParset.get<Solver::VerbosityMode>("verbosity"),
- false); // absolute error
solver.preprocess();
solver.solve();
- };
-
- // Solve the stationary problem
- programState.u[0] = 0.0;
- programState.u[1] = 0.0;
-
- solveLinearProblem(dirichletNodes, contactNetwork.matrices().elasticity, ell0, programState.u,
- parset.sub("u0.solver"));
-
- programState.v[0] = 0.0;
- programState.v[1] = 0.0;
-
-
- programState.a[0] = 0.0;
- programState.a[1] = 0.0;
- {
- // Initial acceleration: Computed in agreement with Ma = ell0 - Au
- // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
- std::vector<Vector> accelerationRHS = ell0;
-
- for (size_t i=0; i<bodyCount; i++) {
- const auto& body = contactNetwork.body(i);
- Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, u[i]);
- }
-
- BitVector noNodes(dirichletNodes.size(), false);
- solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, programState.a,
- parset.sub("a0.solver"));
+ }
- }
+ //print(u, "initial u:");
- for (size_t i=0; i<contactNetwork.nCouplings(); i++) {
- const auto& coupling = contactNetwork.coupling(i);
- const auto& contactCoupling = contactNetwork.nBodyAssembler().getContactCouplings()[i];
+ // Initial acceleration: Computed in agreement with Ma = ell0 - Au
+ // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
+ std::vector<Vector> accelerationRHS = ell0;
+ for (size_t i=0; i<bodyCount; i++) {
+ const auto& body = contactNetwork.body(i);
+ Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, programState.u[i]);
+ }
- const auto nonmortarIdx = coupling->gridIdx_[0];
- const auto& body = contactNetwork.body(nonmortarIdx);
+ std::cout << "solving linear problem for a..." << std::endl;
- ScalarVector frictionBoundaryStress(weightedNormalStress[nonmortarIdx].size());
+ BitVector noNodes(totalDirichletNodes.size(), false);
+ solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, programState.a,
+ parset.sub("a0.solver"));
- body->assembler()->assembleWeightedNormalStress(
- contactCoupling->nonmortarBoundary(), frictionBoundaryStress, body->data()->getYoungModulus(),
- body->data()->getPoissonRatio(), u[nonmortarIdx]);
+ // setup initial conditions in program state
+ programState.alpha[0] = 0;
+ programState.alpha[1] = parset.get<double>("boundary.friction.initialAlpha");
- weightedNormalStress[nonmortarIdx] += frictionBoundaryStress;
- }
+ for (size_t i=0; i<contactNetwork.nCouplings(); i++) {
+ const auto& coupling = contactNetwork.coupling(i);
+ const auto& contactCoupling = contactNetwork.nBodyAssembler().getContactCouplings()[i];
+ const auto nonmortarIdx = coupling->gridIdx_[0];
+ const auto& body = contactNetwork.body(nonmortarIdx);
-*/
+ ScalarVector frictionBoundaryStress(programState.weightedNormalStress[nonmortarIdx].size());
+ body->assembler()->assembleWeightedNormalStress(
+ contactCoupling->nonmortarBoundary(), frictionBoundaryStress, body->data()->getYoungModulus(),
+ body->data()->getPoissonRatio(), programState.u[nonmortarIdx]);
+ programState.weightedNormalStress[nonmortarIdx] += frictionBoundaryStress;
+ } */
- auto& nBodyAssembler = contactNetwork.nBodyAssembler();
for (size_t i=0; i<bodyCount; i++) {
contactNetwork.body(i)->setDeformation(programState.u[i]);
}
@@ -387,46 +390,23 @@ int main(int argc, char *argv[]) {
// Set up TNNMG solver
// -------------------
- BitVector totalDirichletNodes;
- contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
-
- /*for (size_t i=0; i<totalDirichletNodes.size(); i++) {
- bool val = false;
- for (size_t d=0; d<dims; d++) {
- val = val || totalDirichletNodes[i][d];
- }
-
- totalDirichletNodes[i] = val;
- for (size_t d=0; d<dims; d++) {
- totalDirichletNodes[i][d] = val;
- }
- }*/
-
- //print(totalDirichletNodes, "totalDirichletNodes:");
-
//using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, field_type>;
using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
using NonlinearFactory = SolverFactory<Functional, BitVector>;
using BoundaryFunctions = typename ContactNetwork::BoundaryFunctions;
- using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
StateUpdater<ScalarVector, Vector>>;
BoundaryFunctions velocityDirichletFunctions;
contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
- BoundaryNodes dirichletNodes;
- contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
-
/*for (size_t i=0; i<dirichletNodes.size(); i++) {
for (size_t j=0; j<dirichletNodes[i].size(); j++) {
print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
}
}*/
- std::vector<const Dune::BitSetVector<1>*> frictionNodes;
- contactNetwork.frictionNodes(frictionNodes);
/*for (size_t i=0; i<frictionNodes.size(); i++) {
print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
@@ -491,7 +471,7 @@ int main(int argc, char *argv[]) {
timeStepper(stepBase, contactNetwork, current,
programState.relativeTime, programState.relativeTau);
- size_t timeSteps = parset.get<size_t>("timeSteps.timeSteps");
+ size_t timeSteps = std::round(parset.get<double>("timeSteps.timeSteps"));
while (!timeStepper.reachedEnd()) {
programState.timeStep++;