From 9f57bf7c1a8609d57d4f0bd5bb34db22b7d9fa4f Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Wed, 21 Mar 2018 18:08:55 +0100
Subject: [PATCH] .
---
src/assemblers_tmpl.cc | 2 +-
src/multi-body-problem-data/mygrids_tmpl.cc | 12 ++--
src/multi-body-problem.cc | 69 +++++++++++--------
src/spatial-solving/fixedpointiterator.cc | 6 +-
src/spatial-solving/fixedpointiterator.hh | 4 +-
src/time-stepping/adaptivetimestepper.cc | 6 +-
src/time-stepping/adaptivetimestepper.hh | 4 +-
src/time-stepping/coupledtimestepper.cc | 12 ++--
src/time-stepping/coupledtimestepper.hh | 4 +-
src/time-stepping/rate.cc | 2 +-
src/time-stepping/rate.hh | 2 +-
src/time-stepping/rate/backward_euler.cc | 6 +-
src/time-stepping/rate/backward_euler.hh | 2 +-
src/time-stepping/rate/newmark.cc | 6 +-
src/time-stepping/rate/newmark.hh | 2 +-
src/time-stepping/rate/rateupdater.cc | 2 +-
src/time-stepping/rate/rateupdater.hh | 8 +--
src/time-stepping/rate_tmpl.cc | 2 +-
src/time-stepping/state.cc | 2 +-
.../state/sliplawstateupdater.cc | 4 +-
src/time-stepping/state_tmpl.cc | 2 +-
src/vtk_tmpl.cc | 4 +-
22 files changed, 90 insertions(+), 73 deletions(-)
diff --git a/src/assemblers_tmpl.cc b/src/assemblers_tmpl.cc
index e068ad5e..0d5ea0c6 100644
--- a/src/assemblers_tmpl.cc
+++ b/src/assemblers_tmpl.cc
@@ -4,4 +4,4 @@
#include "explicitgrid.hh"
-template class MyAssembler<LeafGridView, MY_DIM>;
+template class MyAssembler<DeformedGrid::LeafGridView, MY_DIM>;
diff --git a/src/multi-body-problem-data/mygrids_tmpl.cc b/src/multi-body-problem-data/mygrids_tmpl.cc
index 2fee2d0e..fc8ccba7 100644
--- a/src/multi-body-problem-data/mygrids_tmpl.cc
+++ b/src/multi-body-problem-data/mygrids_tmpl.cc
@@ -8,14 +8,14 @@
template class GridsConstructor<Grid>;
-template struct MyFaces<LeafGridView>;
-template struct MyFaces<LevelGridView>;
+template struct MyFaces<DeformedGrid::LeafGridView>;
+template struct MyFaces<DeformedGrid::LevelGridView>;
-template MyFaces<LeafGridView> GridsConstructor<Grid>::constructFaces(
- LeafGridView const &gridView, CuboidGeometry const &CuboidGeometry_);
+template MyFaces<DeformedGrid::LeafGridView> GridsConstructor<Grid>::constructFaces(
+ DeformedGrid::LeafGridView const &gridView, CuboidGeometry const &CuboidGeometry_);
-template MyFaces<LevelGridView> GridsConstructor<Grid>::constructFaces(
- LevelGridView const &gridView, CuboidGeometry const &CuboidGeometry_);
+template MyFaces<DeformedGrid::LevelGridView> GridsConstructor<Grid>::constructFaces(
+ DeformedGrid::LevelGridView const &gridView, CuboidGeometry const &CuboidGeometry_);
template void refine<Grid, LocalVector>(
Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index d66d8cc1..54e76531 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -127,11 +127,11 @@ int main(int argc, char *argv[]) {
using DefLeafGridView = DeformedGrid::LeafGridView;
using DefLevelGridView = DeformedGrid::LevelGridView;
- using MyAssembler = MyAssembler<DefLeafGridView, dims>;
- using Matrix = MyAssembler::Matrix;
+ using Assembler = MyAssembler<DefLeafGridView, dims>;
+ using Matrix = Assembler::Matrix;
using LocalVector = Vector::block_type;
- using ScalarMatrix = MyAssembler::ScalarMatrix;
- using ScalarVector = MyAssembler::ScalarVector;
+ using ScalarMatrix = Assembler::ScalarMatrix;
+ using ScalarVector = Assembler::ScalarVector;
using field_type = Matrix::field_type;
// set up material properties of all bodies
@@ -263,16 +263,20 @@ int main(int argc, char *argv[]) {
// set up boundary conditions
std::vector<BoundaryPatch<DefLeafGridView>> neumannBoundaries(bodyCount);
- std::vector<BoundaryPatch<DefLeafGridView>> frictionalBoundaries(bodyCount);
std::vector<BoundaryPatch<DefLeafGridView>> surfaces(bodyCount);
+ // friction boundary
+ std::vector<BoundaryPatch<DefLeafGridView>> frictionalBoundaries(bodyCount);
+ std::vector<Dune::BitSetVector<1>> frictionalNodes(bodyCount);
+
// Dirichlet boundary
std::vector<Dune::BitSetVector<dims>> noNodes(bodyCount);
std::vector<Dune::BitSetVector<dims>> dirichletNodes(bodyCount);
// set up functions for time-dependent boundary conditions
using Function = Dune::VirtualFunction<double, double>;
- const Function& velocityDirichletFunction = VelocityDirichletCondition();
+
+ std::vector<const Function*> velocityDirichletFunctions(grids.size());
const Function& neumannFunction = NeumannCondition();
for (size_t i=0; i<grids.size(); i++) {
@@ -294,9 +298,12 @@ int main(int argc, char *argv[]) {
frictionalBoundary = BoundaryPatch<DefLeafGridView>(leafFaces[i]->lower);
frictionalBoundary.addPatch(BoundaryPatch<DefLeafGridView>(leafFaces[i]->upper));
}
+ frictionalNodes[i] = *frictionalBoundary.getVertices();
//surfaces[i] = BoundaryPatch<GridView>(myFaces.upper);
// TODO: Dirichlet Boundary
+ velocityDirichletFunctions[i] = new VelocityDirichletCondition();
+
noNodes[i] = Dune::BitSetVector<dims>(leafVertexCount);
dirichletNodes[i] = Dune::BitSetVector<dims>(leafVertexCount);
auto& gridDirichletNodes = dirichletNodes[i];
@@ -315,16 +322,16 @@ int main(int argc, char *argv[]) {
}
// set up individual assemblers for each body, assemble problem (matrices, forces, rhs)
- std::vector<std::shared_ptr<MyAssembler>> assemblers(bodyCount);
+ std::vector<std::shared_ptr<Assembler>> assemblers(bodyCount);
Matrices<Matrix> matrices(bodyCount);
std::vector<Vector> gravityFunctionals(bodyCount);
std::vector<std::function<void(double, Vector&)>> externalForces(bodyCount);
- std::vector<EnergyNorm<ScalarMatrix, ScalarVector>* > stateEnergyNorms(bodyCount);
+ std::vector<const EnergyNorm<ScalarMatrix, ScalarVector>* > stateEnergyNorms(bodyCount);
for (size_t i=0; i<assemblers.size(); i++) {
auto& assembler = assemblers[i];
- assembler = std::make_shared<MyAssembler>(*leafViews[i]);
+ assembler = std::make_shared<Assembler>(*leafViews[i]);
assembler->assembleElasticity(body.getYoungModulus(), body.getPoissonRatio(), *matrices.elasticity[i]);
assembler->assembleViscosity(body.getShearViscosityField(), body.getBulkViscosityField(), *matrices.damping[i]);
@@ -401,11 +408,11 @@ int main(int argc, char *argv[]) {
: nullptr;
- /* auto restartIO = handleRestarts
+ auto restartIO = handleRestarts
? std::make_unique<RestartIO<MyProgramState>>(
*restartFile, leafVertexCounts)
- : nullptr;*/
-/*
+ : nullptr;
+
if (firstRestart > 0) // automatically adjusts the time and timestep
restartIO->read(firstRestart, programState);
else
@@ -413,7 +420,7 @@ int main(int argc, char *argv[]) {
matrices, assemblers, dirichletNodes,
noNodes, frictionalBoundaries, body);
- */
+
// assemble friction
std::vector<std::shared_ptr<MyGlobalFrictionData<LocalVector>>> frictionInfo(weakPatches.size());
@@ -434,18 +441,19 @@ int main(int argc, char *argv[]) {
vertexCoordinates[vertexMapper.index(v)] = geoToPoint(v.geometry());
}
- using MyVertexBasis = typename MyAssembler::VertexBasis;
+ using MyVertexBasis = typename Assembler::VertexBasis;
auto dataWriter =
writeData ? std::make_unique<
HDF5Writer<MyProgramState, MyVertexBasis, GridView>>(
*dataFile, vertexCoordinates, myAssembler.vertexBasis,
surface, frictionalBoundary, weakPatch)
: nullptr;
- MyVTKWriter<MyVertexBasis, typename MyAssembler::CellBasis> const vtkWriter(
+ MyVTKWriter<MyVertexBasis, typename Assembler::CellBasis> const vtkWriter(
myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
-
+ */
IterationRegister iterationCount;
+ /*
auto const report = [&](bool initial = false) {
if (writeData) {
dataWriter->reportSolution(programState, *myGlobalFriction);
@@ -485,17 +493,18 @@ int main(int argc, char *argv[]) {
using MyUpdater = Updaters<RateUpdater<Vector, Matrix, Function, dims>,
StateUpdater<ScalarVector, Vector>>;
- /*
+
+ std::vector<double> L(bodyCount, parset.get<double>("boundary.friction.L"));
+ std::vector<double> V0(bodyCount, parset.get<double>("boundary.friction.V0"));
MyUpdater current(
initRateUpdater(parset.get<Config::scheme>("timeSteps.scheme"),
- velocityDirichletFunction, dirichletNodes, matrices,
+ velocityDirichletFunctions, dirichletNodes, matrices,
programState.u, programState.v, programState.a),
initStateUpdater<ScalarVector, Vector>(
parset.get<Config::stateModel>("boundary.friction.stateModel"),
- programState.alpha, *frictionalBoundary.getVertices(),
- parset.get<double>("boundary.friction.L"),
- parset.get<double>("boundary.friction.V0")));
- */
+ programState.alpha, frictionalNodes,
+ L, V0));
+
auto const refinementTolerance =
parset.get<double>("timeSteps.refinementTolerance");
@@ -514,16 +523,16 @@ int main(int argc, char *argv[]) {
return energyNorm > refinementTolerance;
};
- /*
+
std::signal(SIGXCPU, handleSignal);
std::signal(SIGINT, handleSignal);
std::signal(SIGTERM, handleSignal);
AdaptiveTimeStepper<NonlinearFactory, MyUpdater,
EnergyNorm<ScalarMatrix, ScalarVector>>
- adaptiveTimeStepper(factory, parset, myGlobalFriction, current,
+ adaptiveTimeStepper(factory, parset, globalFriction, current,
programState.relativeTime, programState.relativeTau,
- computeExternalForces, stateEnergyNorm, mustRefine);
+ externalForces, stateEnergyNorms, mustRefine);
while (!adaptiveTimeStepper.reachedEnd()) {
programState.timeStep++;
@@ -536,14 +545,20 @@ int main(int argc, char *argv[]) {
current.rate_->extractAcceleration(programState.a);
current.state_->extractAlpha(programState.alpha);
- report();
+ for (size_t i=0; i<deformedGridComplex.size() ; i++) {
+ deformedGridComplex.setDeformation(programState.u[i], i);
+ }
+ nBodyAssembler.assembleTransferOperator();
+ nBodyAssembler.assembleObstacle();
+
+ // report();
if (terminationRequested) {
std::cerr << "Terminating prematurely" << std::endl;
break;
}
}
-*/
+
} catch (Dune::Exception &e) {
Dune::derr << "Dune reported error: " << e << std::endl;
} catch (std::exception &e) {
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index 521c513b..18931b36 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -36,7 +36,7 @@ void FixedPointIterationCounter::operator+=(
template <class Factory, class Updaters, class ErrorNorm>
FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
Factory &factory, Dune::ParameterTree const &parset,
- std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, const ErrorNorm& errorNorm)
+ std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, const std::vector<const ErrorNorm* >& errorNorms)
: factory_(factory),
step_(factory_.getStep()),
parset_(parset),
@@ -47,7 +47,7 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
velocityMaxIterations_(parset.get<size_t>("v.solver.maximumIterations")),
velocityTolerance_(parset.get<double>("v.solver.tolerance")),
verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
- errorNorm_(errorNorm) {}
+ errorNorms_(errorNorms) {}
template <class Factory, class Updaters, class ErrorNorm>
FixedPointIterationCounter
@@ -134,7 +134,7 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
bool breakCriterion = true;
for (size_t i=0; i<alpha.size(); i++) {
- if (errorNorm_.diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
+ if (errorNorms_[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
breakCriterion = false;
break;
}
diff --git a/src/spatial-solving/fixedpointiterator.hh b/src/spatial-solving/fixedpointiterator.hh
index e579645b..bae4624f 100644
--- a/src/spatial-solving/fixedpointiterator.hh
+++ b/src/spatial-solving/fixedpointiterator.hh
@@ -37,7 +37,7 @@ class FixedPointIterator {
public:
FixedPointIterator(Factory& factory, const Dune::ParameterTree& parset,
std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
- const ErrorNorm& errorNorm);
+ const std::vector<const ErrorNorm* >& errorNorms);
FixedPointIterationCounter run(Updaters updaters,
const std::vector<Matrix>& velocityMatrices,
@@ -56,6 +56,6 @@ class FixedPointIterator {
size_t velocityMaxIterations_;
double velocityTolerance_;
Solver::VerbosityMode verbosity_;
- ErrorNorm const &errorNorm_;
+ const std::vector<const ErrorNorm* >& errorNorms_;
};
#endif
diff --git a/src/time-stepping/adaptivetimestepper.cc b/src/time-stepping/adaptivetimestepper.cc
index 0808750f..1c36768e 100644
--- a/src/time-stepping/adaptivetimestepper.cc
+++ b/src/time-stepping/adaptivetimestepper.cc
@@ -23,7 +23,7 @@ AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, Updaters ¤t,
double relativeTime, double relativeTau,
std::vector<std::function<void(double, Vector &)>>& externalForces,
- ErrorNorm const &errorNorm,
+ const std::vector<const ErrorNorm* >& errorNorms,
std::function<bool(Updaters &, Updaters &)> mustRefine)
: relativeTime_(relativeTime),
relativeTau_(relativeTau),
@@ -35,7 +35,7 @@ AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
R1_(),
externalForces_(externalForces),
mustRefine_(mustRefine),
- errorNorm_(errorNorm) {}
+ errorNorms_(errorNorms) {}
template <class Factory, class Updaters, class ErrorNorm>
bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::reachedEnd() {
@@ -100,7 +100,7 @@ AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::step(
UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
newUpdatersAndCount.count =
MyCoupledTimeStepper(finalTime_, factory_, parset_, globalFriction_,
- newUpdatersAndCount.updaters, errorNorm_,
+ newUpdatersAndCount.updaters, errorNorms_,
externalForces_)
.step(rTime, rTau);
iterationRegister_.registerCount(newUpdatersAndCount.count);
diff --git a/src/time-stepping/adaptivetimestepper.hh b/src/time-stepping/adaptivetimestepper.hh
index a0bf6113..f7897a29 100644
--- a/src/time-stepping/adaptivetimestepper.hh
+++ b/src/time-stepping/adaptivetimestepper.hh
@@ -34,7 +34,7 @@ class AdaptiveTimeStepper {
Updaters ¤t, double relativeTime,
double relativeTau,
std::vector<std::function<void(double, Vector &)>>& externalForces,
- ErrorNorm const &errorNorm,
+ const std::vector<const ErrorNorm* >& errorNorms,
std::function<bool(Updaters &, Updaters &)> mustRefine);
bool reachedEnd();
@@ -55,7 +55,7 @@ class AdaptiveTimeStepper {
UpdatersWithCount R1_;
std::vector<std::function<void(double, Vector &)>>& externalForces_;
std::function<bool(Updaters &, Updaters &)> mustRefine_;
- ErrorNorm const &errorNorm_;
+ const std::vector<const ErrorNorm* >& errorNorms_;
IterationRegister iterationRegister_;
};
diff --git a/src/time-stepping/coupledtimestepper.cc b/src/time-stepping/coupledtimestepper.cc
index 7501f68a..74e66fc8 100644
--- a/src/time-stepping/coupledtimestepper.cc
+++ b/src/time-stepping/coupledtimestepper.cc
@@ -8,7 +8,7 @@ template <class Factory, class Updaters, class ErrorNorm>
CoupledTimeStepper<Factory, Updaters, ErrorNorm>::CoupledTimeStepper(
double finalTime, Factory &factory, Dune::ParameterTree const &parset,
std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, Updaters updaters,
- ErrorNorm const &errorNorm,
+ const std::vector<const ErrorNorm* >& errorNorms,
std::vector<std::function<void(double, Vector &)>>& externalForces)
: finalTime_(finalTime),
factory_(factory),
@@ -16,7 +16,7 @@ CoupledTimeStepper<Factory, Updaters, ErrorNorm>::CoupledTimeStepper(
globalFriction_(globalFriction),
updaters_(updaters),
externalForces_(externalForces),
- errorNorm_(errorNorm) {}
+ errorNorms_(errorNorms) {}
template <class Factory, class Updaters, class ErrorNorm>
FixedPointIterationCounter
@@ -26,8 +26,10 @@ CoupledTimeStepper<Factory, Updaters, ErrorNorm>::step(double relativeTime,
updaters_.rate_->nextTimeStep();
auto const newRelativeTime = relativeTime + relativeTau;
- std::vector<Vector> ell;
- externalForces_(newRelativeTime, ell);
+ std::vector<Vector> ell(externalForces_.size());
+ for (size_t i=0; i<externalForces_.size(); i++) {
+ externalForces_[i](newRelativeTime, ell[i]);
+ }
std::vector<Matrix> velocityMatrix;
std::vector<Vector> velocityRHS;
@@ -38,7 +40,7 @@ CoupledTimeStepper<Factory, Updaters, ErrorNorm>::step(double relativeTime,
updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS, velocityIterate, velocityMatrix);
FixedPointIterator<Factory, Updaters, ErrorNorm> fixedPointIterator(
- factory_, parset_, globalFriction_, errorNorm_);
+ factory_, parset_, globalFriction_, errorNorms_);
auto const iterations = fixedPointIterator.run(updaters_, velocityMatrix,
velocityRHS, velocityIterate);
return iterations;
diff --git a/src/time-stepping/coupledtimestepper.hh b/src/time-stepping/coupledtimestepper.hh
index 92e3cd4a..97c01282 100644
--- a/src/time-stepping/coupledtimestepper.hh
+++ b/src/time-stepping/coupledtimestepper.hh
@@ -18,7 +18,7 @@ class CoupledTimeStepper {
CoupledTimeStepper(double finalTime, Factory &factory,
Dune::ParameterTree const &parset,
std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
- Updaters updaters, ErrorNorm const &errorNorm,
+ Updaters updaters, const std::vector<const ErrorNorm* >& errorNorms,
std::vector<std::function<void(double, Vector &)>>& externalForces);
FixedPointIterationCounter step(double relativeTime, double relativeTau);
@@ -30,6 +30,6 @@ class CoupledTimeStepper {
std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
Updaters updaters_;
std::vector<std::function<void(double, Vector &)>>& externalForces_;
- ErrorNorm const &errorNorm_;
+ const std::vector<const ErrorNorm* >& errorNorms_;
};
#endif
diff --git a/src/time-stepping/rate.cc b/src/time-stepping/rate.cc
index fe366fe7..ce4050c7 100644
--- a/src/time-stepping/rate.cc
+++ b/src/time-stepping/rate.cc
@@ -9,7 +9,7 @@
template <class Vector, class Matrix, class Function, int dimension>
std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
initRateUpdater(Config::scheme scheme,
- const std::vector<Function>& velocityDirichletFunctions,
+ const std::vector<const Function*>& velocityDirichletFunctions,
const std::vector<Dune::BitSetVector<dimension>>& velocityDirichletNodes,
const Matrices<Matrix>& matrices,
const std::vector<Vector>& u_initial,
diff --git a/src/time-stepping/rate.hh b/src/time-stepping/rate.hh
index 60268240..5d2ae428 100644
--- a/src/time-stepping/rate.hh
+++ b/src/time-stepping/rate.hh
@@ -9,7 +9,7 @@
template <class Vector, class Matrix, class Function, int dimension>
std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
initRateUpdater(Config::scheme scheme,
- const std::vector<Function>& velocityDirichletFunctions,
+ const std::vector<const Function* >& velocityDirichletFunctions,
const std::vector<Dune::BitSetVector<dimension>>& velocityDirichletNodes,
const Matrices<Matrix>& matrices,
const std::vector<Vector>& u_initial,
diff --git a/src/time-stepping/rate/backward_euler.cc b/src/time-stepping/rate/backward_euler.cc
index 2694a970..11b7907f 100644
--- a/src/time-stepping/rate/backward_euler.cc
+++ b/src/time-stepping/rate/backward_euler.cc
@@ -8,7 +8,7 @@ BackwardEuler<Vector, Matrix, Function, dim>::BackwardEuler(
const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
- const std::vector<Function>& _dirichletFunctions)
+ const std::vector<const Function*>& _dirichletFunctions)
: RateUpdater<Vector, Matrix, Function, dim>(
_matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
_dirichletFunctions) {}
@@ -20,7 +20,7 @@ void BackwardEuler<Vector, Matrix, Function, dim>::setup(const std::vector<Vecto
std::vector<Vector>& rhs, std::vector<Vector>& iterate,
std::vector<Matrix>& AM) {
for (size_t i=0; i<this->u.size(); i++) {
- this->dirichletFunctions[i].evaluate(relativeTime, this->dirichletValues[i]);
+ this->dirichletFunctions[i]->evaluate(relativeTime, this->dirichletValue);
this->tau = _tau;
/* We start out with the formulation
@@ -73,7 +73,7 @@ void BackwardEuler<Vector, Matrix, Function, dim>::setup(const std::vector<Vecto
const Dune::BitSetVector<dim>& dirichletNodess = this->dirichletNodes[i];
for (size_t k = 0; k < dirichletNodess.size(); ++k)
for (size_t j = 0; j < dim; ++j)
- if (this->dirichletNodes[k][j])
+ if (dirichletNodess[k][j])
iterate[k][j] = (j == 0) ? this->dirichletValue : 0;
}
}
diff --git a/src/time-stepping/rate/backward_euler.hh b/src/time-stepping/rate/backward_euler.hh
index 245bb74b..169323d6 100644
--- a/src/time-stepping/rate/backward_euler.hh
+++ b/src/time-stepping/rate/backward_euler.hh
@@ -7,7 +7,7 @@ class BackwardEuler : public RateUpdater<Vector, Matrix, Function, dim> {
BackwardEuler(const Matrices<Matrix> &_matrices, const std::vector<Vector> &_u_initial,
const std::vector<Vector> &_v_initial, const std::vector<Vector> &_a_initial,
const std::vector<Dune::BitSetVector<dim> > &_dirichletNodes,
- const std::vector<Function> &_dirichletFunctions);
+ const std::vector<const Function*> &_dirichletFunctions);
void setup(const std::vector<Vector>&, double, double, std::vector<Vector>&, std::vector<Vector>&,
std::vector<Matrix>&) override;
diff --git a/src/time-stepping/rate/newmark.cc b/src/time-stepping/rate/newmark.cc
index 53bf32c4..c88136d7 100644
--- a/src/time-stepping/rate/newmark.cc
+++ b/src/time-stepping/rate/newmark.cc
@@ -8,7 +8,7 @@ Newmark<Vector, Matrix, Function, dim>::Newmark(
const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
- const std::vector<Function>& _dirichletFunctions)
+ const std::vector<const Function*>& _dirichletFunctions)
: RateUpdater<Vector, Matrix, Function, dim>(
_matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
_dirichletFunctions) {}
@@ -20,7 +20,7 @@ void Newmark<Vector, Matrix, Function, dim>::setup(const std::vector<Vector>& el
std::vector<Vector>& rhs, std::vector<Vector>& iterate,
std::vector<Matrix>& AM) {
for (size_t i=0; i<this->u.size(); i++) {
- this->dirichletFunctions[i].evaluate(relativeTime, this->dirichletValues[i]);
+ this->dirichletFunctions[i]->evaluate(relativeTime, this->dirichletValue);
this->tau = _tau;
/* We start out with the formulation
@@ -77,7 +77,7 @@ void Newmark<Vector, Matrix, Function, dim>::setup(const std::vector<Vector>& el
const Dune::BitSetVector<dim>& dirichletNodess = this->dirichletNodes[i];
for (size_t k = 0; k < dirichletNodess.size(); ++k)
for (size_t j = 0; j < dim; ++j)
- if (this->dirichletNodess[k][j])
+ if (dirichletNodess[k][j])
iterate[k][j] = (j == 0) ? this->dirichletValue : 0;
}
}
diff --git a/src/time-stepping/rate/newmark.hh b/src/time-stepping/rate/newmark.hh
index be5cb844..25e014ea 100644
--- a/src/time-stepping/rate/newmark.hh
+++ b/src/time-stepping/rate/newmark.hh
@@ -7,7 +7,7 @@ class Newmark : public RateUpdater<Vector, Matrix, Function, dim> {
Newmark(const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
- const std::vector<Function>& _dirichletFunctions);
+ const std::vector<const Function*>& _dirichletFunctions);
void setup(const std::vector<Vector>&, double, double, std::vector<Vector>&, std::vector<Vector>&,
std::vector<Matrix>&) override;
diff --git a/src/time-stepping/rate/rateupdater.cc b/src/time-stepping/rate/rateupdater.cc
index 23d09536..6c095235 100644
--- a/src/time-stepping/rate/rateupdater.cc
+++ b/src/time-stepping/rate/rateupdater.cc
@@ -9,7 +9,7 @@ RateUpdater<Vector, Matrix, Function, dim>::RateUpdater(
const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
- const std::vector<Function>& _dirichletFunctions)
+ const std::vector<const Function*>& _dirichletFunctions)
: matrices(_matrices),
u(_u_initial),
v(_v_initial),
diff --git a/src/time-stepping/rate/rateupdater.hh b/src/time-stepping/rate/rateupdater.hh
index 5266259d..8b78e1d3 100644
--- a/src/time-stepping/rate/rateupdater.hh
+++ b/src/time-stepping/rate/rateupdater.hh
@@ -13,12 +13,12 @@ class RateUpdater {
RateUpdater(const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
- const std::vector<Function>& _dirichletFunctions);
+ const std::vector<const Function*>& _dirichletFunctions);
public:
void nextTimeStep();
- void virtual setup(Vector const &ell, double _tau, double relativeTime,
- Vector &rhs, Vector &iterate, Matrix &AB) = 0;
+ void virtual setup(const std::vector<Vector>& ell, double _tau, double relativeTime,
+ std::vector<Vector>& rhs, std::vector<Vector>& iterate, std::vector<Matrix>& AB) = 0;
void virtual postProcess(const std::vector<Vector>& iterate) = 0;
void extractDisplacement(std::vector<Vector>& displacements) const;
@@ -33,7 +33,7 @@ class RateUpdater {
const Matrices<Matrix>& matrices;
std::vector<Vector> u, v, a;
const std::vector<Dune::BitSetVector<dim>>& dirichletNodes;
- const std::vector<Function>& dirichletFunctions;
+ const std::vector<const Function*>& dirichletFunctions;
double dirichletValue;
std::vector<Vector> u_o, v_o, a_o;
diff --git a/src/time-stepping/rate_tmpl.cc b/src/time-stepping/rate_tmpl.cc
index 92c15bdd..18593e15 100644
--- a/src/time-stepping/rate_tmpl.cc
+++ b/src/time-stepping/rate_tmpl.cc
@@ -8,7 +8,7 @@ template std::shared_ptr<RateUpdater<Vector, Matrix, Function, MY_DIM>>
initRateUpdater<Vector, Matrix, Function, MY_DIM>(
Config::scheme scheme,
const std::vector<Function>& velocityDirichletFunctions,
- const std::vector<Dune::BitSetVector<dimension>>& velocityDirichletNodes,
+ const std::vector<Dune::BitSetVector<MY_DIM>>& velocityDirichletNodes,
const Matrices<Matrix>& matrices,
const std::vector<Vector>& u_initial,
const std::vector<Vector>& v_initial,
diff --git a/src/time-stepping/state.cc b/src/time-stepping/state.cc
index f07a8033..b8e4349e 100644
--- a/src/time-stepping/state.cc
+++ b/src/time-stepping/state.cc
@@ -8,7 +8,7 @@
template <class ScalarVector, class Vector>
std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
- const std::vector<Dune::BitSetVector<_Tp1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0) {
+ const std::vector<Dune::BitSetVector<1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0) {
switch (model) {
case Config::AgeingLaw:
return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(
diff --git a/src/time-stepping/state/sliplawstateupdater.cc b/src/time-stepping/state/sliplawstateupdater.cc
index 6b322d89..9a36d0a3 100644
--- a/src/time-stepping/state/sliplawstateupdater.cc
+++ b/src/time-stepping/state/sliplawstateupdater.cc
@@ -36,8 +36,8 @@ void SlipLawStateUpdater<ScalarVector, Vector>::solve(
double const V = velocity_fieldi[j].two_norm();
double const mtoL = -tau * V / L[i];
- alphai[j] = (V <= 0) ? alpha_oi[j] : std::expm1(mtVoL) * std::log(V / V0[i]) +
- alpha_oi[j] * std::exp(mtVoL);
+ alphai[j] = (V <= 0) ? alpha_oi[j] : std::expm1(mtoL) * std::log(V / V0[i]) +
+ alpha_oi[j] * std::exp(mtoL);
}
}
}
diff --git a/src/time-stepping/state_tmpl.cc b/src/time-stepping/state_tmpl.cc
index 5703f81a..718150c7 100644
--- a/src/time-stepping/state_tmpl.cc
+++ b/src/time-stepping/state_tmpl.cc
@@ -3,4 +3,4 @@
template std::shared_ptr<StateUpdater<ScalarVector, Vector>>
initStateUpdater<ScalarVector, Vector>(
Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
- const std::vector<Dune::BitSetVector<_Tp1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0);
+ const std::vector<Dune::BitSetVector<1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0);
diff --git a/src/vtk_tmpl.cc b/src/vtk_tmpl.cc
index f6790e5c..ecfecf9a 100644
--- a/src/vtk_tmpl.cc
+++ b/src/vtk_tmpl.cc
@@ -8,8 +8,8 @@
#include <dune/fufem/functionspacebases/p0basis.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-using MyP0Basis = P0Basis<GridView, double>;
-using P1Basis = P1NodalBasis<GridView, double>;
+using MyP0Basis = P0Basis<LeafGridView, double>;
+using P1Basis = P1NodalBasis<LeafGridView, double>;
template class MyVTKWriter<P1Basis, MyP0Basis>;
--
GitLab