#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/contact/projections/normalprojection.hh>
#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
#include <dune/fufem/functions/constantfunction.hh>
#include "../friction/globalratestatefriction.hh"
#include "../friction/frictionpotential.hh"
#include "contactnetwork.hh"
#include "../../assemblers.hh"
#include "../../utils/tobool.hh"
template <class HostGridType, class VectorType>
ContactNetwork<HostGridType, VectorType>::ContactNetwork(
size_t nBodies,
size_t nCouplings,
size_t nLevels) :
levelContactNetworks_(nLevels),
bodies_(nBodies),
couplings_(nCouplings),
frictionBoundaries_(nBodies),
stateEnergyNorms_(nBodies),
nBodyAssembler_(nBodies, nCouplings, 0.0)
{}
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::addLevel(std::shared_ptr<LevelContactNetwork> levelContactNetwork) {
size_t level = levelContactNetwork->level();
if (level >= levelContactNetworks_.size()) {
levelContactNetworks_.resize(level);
}
levelContactNetworks_[level] = levelContactNetwork;
}
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::addLevel(const std::vector<size_t>& bodyLevels, const size_t level) {
assert(bodyLevels.size() == nBodies());
if (level >= levelContactNetworks_.size()) {
levelContactNetworks_.resize(level+1);
}
auto& levelContactNetwork = levelContactNetworks_[level];
levelContactNetwork = std::make_shared<LevelContactNetwork>(nBodies(), nCouplings(), level);
std::vector<std::shared_ptr<Body>> bodies(nBodies());
for (size_t id=0; id<nBodies(); id++) {
const auto& body = this->body(id);
bodies[id] = std::make_shared<Body>(body->data(), body->grid(), std::min((size_t) body->grid()->maxLevel(), bodyLevels[id]));
}
levelContactNetwork->setBodies(bodies);
levelContactNetwork->setCouplings(this->couplings_);
levelContactNetwork->build();
}
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::assemble() {
std::vector<const GridType*> grids(nBodies());
for (size_t i=0; i<nBodies(); i++) {
grids[i] = bodies_[i]->grid().get();
bodies_[i]->assemble();
frictionBoundaries_[i] = std::make_unique<BoundaryPatch>(bodies_[i]->gridView(), false);
}
// set up dune-contact nBodyAssembler
nBodyAssembler_.setGrids(grids);
for (size_t i=0; i<nCouplings(); i++) {
nBodyAssembler_.setCoupling(*couplings_[i], i);
}
std::vector<std::shared_ptr<Dune::BitSetVector<1>>> dirichletVertices(nBodies());
for (size_t i=0; i<nBodies(); i++) {
const auto& body = this->body(i);
std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> boundaryConditions;
body->boundaryConditions("dirichlet", boundaryConditions);
dirichletVertices[i] = std::make_shared<Dune::BitSetVector<1>>(body->nVertices());
auto& vertices = *dirichletVertices[i];
vertices.unsetAll();
if (boundaryConditions.size()<1)
continue;
for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
const auto& boundaryVertices = *boundaryConditions[bc]->boundaryPatch()->getVertices();
for (size_t j=0; j<boundaryVertices.size(); j++) {
vertices[j][0] = vertices[j][0] or boundaryVertices[j][0];
}
}
}
nBodyAssembler_.setDirichletVertices(dirichletVertices);
nBodyAssembler_.assembleTransferOperator();
nBodyAssembler_.assembleObstacle();
for (size_t i=0; i<nCouplings(); i++) {
auto& coupling = couplings_[i];
const auto& contactCoupling = nBodyAssembler_.getContactCouplings()[i]; // dual mortar object holding boundary patches
const auto nonmortarIdx = coupling->gridIdx_[0];
//const auto mortarIdx = coupling->gridIdx_[1];
frictionBoundaries_[nonmortarIdx]->addPatch(contactCoupling->nonmortarBoundary());
//frictionBoundaries_[mortarIdx]->addPatch(contactCoupling->mortarBoundary());
}
// assemble state energy norm
frictionBoundaryMass_.resize(nBodies());
for (size_t i=0; i<nBodies(); i++) {
frictionBoundaryMass_[i] = std::make_unique<ScalarMatrix>();
bodies_[i]->assembler()->assembleFrictionalBoundaryMass(*frictionBoundaries_[i], *frictionBoundaryMass_[i]);
*frictionBoundaryMass_[i] /= frictionBoundaries_[i]->area(); // TODO: weight by individual friction patches?
stateEnergyNorms_[i] = std::make_unique<typename LeafBody::StateEnergyNorm>(*frictionBoundaryMass_[i]);
}
}
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::assembleFriction(
const Config::FrictionModel& frictionModel,
const std::vector<ScalarVector>& weightedNormalStress) {
assert(weightedNormalStress.size() == bodies_.size());
const size_t nBodies = bodies_.size();
// Lumping of the nonlinearity
std::vector<ScalarVector> weights(nBodies);
{
NeumannBoundaryAssembler<GridType, typename ScalarVector::block_type>
frictionalBoundaryAssembler(std::make_shared<
ConstantFunction<LocalVector, typename ScalarVector::block_type>>(
1));
for (size_t i=0; i<nBodies; i++) {
bodies_[i]->assembler()->assembleBoundaryFunctional(frictionalBoundaryAssembler, weights[i], *frictionBoundaries_[i]);
}
}
/* globalFriction_ = std::make_shared<GlobalRateStateFriction<
Matrix, Vector, ZeroFunction, DeformedGridType>>(
nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress); */
switch (frictionModel) {
case Config::Truncated:
globalFriction_ = std::make_shared<GlobalRateStateFriction<
Matrix, VectorType, TruncatedRateState, GridType>>(
nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
break;
case Config::Regularised:
globalFriction_ = std::make_shared<GlobalRateStateFriction<
Matrix, VectorType, RegularisedRateState, GridType>>(
nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
break;
case Config::Tresca:
globalFriction_ = std::make_shared<GlobalRateStateFriction<
Matrix, VectorType, Tresca, GridType>>(
nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
break;
case Config::None:
globalFriction_ = std::make_shared<GlobalRateStateFriction<
Matrix, VectorType, ZeroFunction, GridType>>(
nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
break;
default:
assert(false);
break;
}
}
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::setBodies(const std::vector<std::shared_ptr<LeafBody>> bodies) {
assert(nBodies()==bodies.size());
bodies_ = bodies;
matrices_.elasticity.resize(nBodies());
matrices_.damping.resize(nBodies());
matrices_.mass.resize(nBodies());
for (size_t i=0; i<nBodies(); i++) {
matrices_.elasticity[i] = bodies_[i]->matrices().elasticity;
matrices_.damping[i] = bodies_[i]->matrices().damping;
matrices_.mass[i] = bodies_[i]->matrices().mass;
}
}
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings) {
assert(this->nCouplings()==couplings.size());
couplings_ = couplings;
}
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::constructBody(
const std::shared_ptr<BodyData<dim>>& bodyData,
const std::shared_ptr<HostGridType>& grid,
std::shared_ptr<LeafBody>& body) const {
body = std::make_shared<LeafBody>(bodyData, grid);
}
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::setDeformation(const std::vector<VectorType>& totalDeformation) {
assert(totalDeformation.size() == nBodies());
for (size_t i=0; i<nBodies(); i++) {
body(i)->setDeformation(totalDeformation[i]);
}
nBodyAssembler_.assembleTransferOperator();
nBodyAssembler_.assembleObstacle();
for (size_t i=0; i<levelContactNetworks_.size(); i++) {
levelContactNetworks_[i]->build();
}
}
template <class HostGridType, class VectorType>
template <class LevelBoundaryPatch>
void ContactNetwork<HostGridType, VectorType>::constructCoupling(
int nonMortarBodyIdx, int mortarBodyIdx,
const std::shared_ptr<LevelBoundaryPatch>& nonmortarPatch,
const std::shared_ptr<LevelBoundaryPatch>& mortarPatch,
const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch,
const std::shared_ptr<GlobalFrictionData<dim>>& frictionData,
std::shared_ptr<FrictionCouplingPair>& coupling) const {
coupling = std::make_shared<FrictionCouplingPair>();
auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<BoundaryPatch>>();
std::shared_ptr<typename FrictionCouplingPair::BackEndType> backend = nullptr;
coupling->set(nonMortarBodyIdx, mortarBodyIdx, nonmortarPatch, mortarPatch, 0.1, Dune::Contact::CouplingPairBase::STICK_SLIP, contactProjection, backend);
coupling->setWeakeningPatch(weakeningPatch);
coupling->setFrictionData(frictionData);
}
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::level(size_t level) const
-> const std::shared_ptr<LevelContactNetwork>& {
return levelContactNetworks_[level];
}
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::nLevels() const
-> size_t {
return levelContactNetworks_.size();
}
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::nBodies() const
-> size_t {
return bodies_.size();
}
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::nCouplings() const
-> size_t {
return couplings_.size();
}
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::body(int i) const
-> const std::shared_ptr<LeafBody>& {
return bodies_[i];
}
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::coupling(int i) const
-> const std::shared_ptr<FrictionCouplingPair>& {
return couplings_[i];
}
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::couplings() const
-> const std::vector<std::shared_ptr<FrictionCouplingPair>>& {
return couplings_;
}
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::nBodyAssembler()
-> NBodyAssembler& {
return nBodyAssembler_;
}
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::nBodyAssembler() const
-> const NBodyAssembler& {
return nBodyAssembler_;
}
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::matrices() const
-> const Matrices<Matrix,2>& {
return matrices_;
}
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::stateEnergyNorms() const
-> const StateEnergyNorms& {
return stateEnergyNorms_;
}
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::globalFriction() const
-> GlobalFriction& {
return *globalFriction_;
}
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::externalForces(ExternalForces& externalForces) const {
externalForces.resize(nBodies());
for (size_t i=0; i<nBodies(); i++) {
externalForces[i] = std::make_unique<typename LeafBody::ExternalForce>(bodies_[i]->externalForce());
}
}
// collects all leaf boundary nodes from the different bodies identified by "tag" in a single BitSetVector
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::totalNodes(
const std::string& tag,
Dune::BitSetVector<dim>& totalNodes) const {
int totalSize = 0;
for (size_t i=0; i<nBodies(); i++) {
totalSize += this->body(i)->nVertices();
}
totalNodes.resize(totalSize);
int idx=0;
for (size_t i=0; i<nBodies(); i++) {
const auto& body = this->body(i);
std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> boundaryConditions;
body->boundaryConditions(tag, 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<dim; k++)
totalNodes[idx][k] = totalNodes[idx][k] or (*nodes)[j][k];
idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
}
} else {
idx += body->nVertices();
}
}
}
// collects all leaf boundary nodes from the different bodies identified by "tag" in a vector of BitSetVector
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::boundaryPatchNodes(
const std::string& tag,
BoundaryPatchNodes& nodes) const {
nodes.resize(nBodies());
for (size_t i=0; i<nBodies(); i++) {
this->body(i)->boundaryPatchNodes(tag, nodes[i]);
}
}
// collects all leaf boundary nodes from the different bodies identified by "tag" in a vector of BitSetVector
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::boundaryNodes(
const std::string& tag,
BoundaryNodes& nodes) const {
nodes.resize(nBodies());
for (size_t i=0; i<nBodies(); i++) {
this->body(i)->boundaryNodes(tag, nodes[i]);
}
}
// collects all leaf boundary functions from the different bodies identified by "tag"
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::boundaryFunctions(
const std::string& tag,
BoundaryFunctions& functions) const {
functions.resize(nBodies());
for (size_t i=0; i<nBodies(); i++) {
this->body(i)->boundaryFunctions(tag, functions[i]);
}
}
// collects all leaf boundary patches from the different bodies identified by "tag"
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::boundaryPatches(
const std::string& tag,
BoundaryPatches& patches) const {
patches.resize(nBodies());
for (size_t i=0; i<nBodies(); i++) {
this->body(i)->boundaryPatches(tag, patches[i]);
}
}
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::frictionNodes(std::vector<const Dune::BitSetVector<1>*>& nodes) const {
nodes.resize(nBodies());
for (size_t i=0; i<nBodies(); i++) {
nodes[i] = frictionBoundaries_[i]->getVertices();
}
}
template <class HostGridType, class VectorType>
void ContactNetwork<HostGridType, VectorType>::frictionPatches(std::vector<const BoundaryPatch*>& patches) const {
patches.resize(nBodies());
for (size_t i=0; i<nBodies(); i++) {
patches[i] = frictionBoundaries_[i].get();
}
}
#include "contactnetwork_tmpl.cc"