From 09929cc35af78fe163de54381ccd391d9b2e8c52 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Fri, 6 Sep 2019 19:24:02 +0200
Subject: [PATCH] towards crosspoints
---
.../friction/globalratestatefriction.hh | 4 +-
.../data-structures/network/contactnetwork.cc | 28 ++-
.../data-structures/network/contactnetwork.hh | 5 +-
.../tectonic/data-structures/program_state.hh | 2 +-
.../factories/stackedblocksfactory.hh | 2 +-
dune/tectonic/factories/threeblocksfactory.cc | 8 +-
dune/tectonic/factories/threeblocksfactory.hh | 2 +-
.../contact/dualmortarcoupling.cc | 178 +++++++++++-------
.../contact/dualmortarcoupling.hh | 31 ++-
.../spatial-solving/contact/nbodyassembler.cc | 70 +++++--
.../spatial-solving/contact/nbodyassembler.hh | 41 ++--
.../spatial-solving/fixedpointiterator.cc | 2 +-
.../fixedpointiterator_tmpl.cc | 6 +-
.../multilevelpatchpreconditioner.hh | 2 +-
.../time-stepping/adaptivetimestepper_tmpl.cc | 4 +-
.../time-stepping/coupledtimestepper_tmpl.cc | 6 +-
dune/tectonic/time-stepping/state_tmpl.cc | 3 +-
src/multi-body-problem/multi-body-problem.cc | 5 +-
18 files changed, 256 insertions(+), 143 deletions(-)
diff --git a/dune/tectonic/data-structures/friction/globalratestatefriction.hh b/dune/tectonic/data-structures/friction/globalratestatefriction.hh
index 0d007d8e..ec42936b 100644
--- a/dune/tectonic/data-structures/friction/globalratestatefriction.hh
+++ b/dune/tectonic/data-structures/friction/globalratestatefriction.hh
@@ -10,7 +10,7 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
-#include <dune/contact/assemblers/dualmortarcoupling.hh>
+#include "../../spatial-solving/contact/dualmortarcoupling.hh"
#include "globalfrictiondata.hh"
#include "globalfriction.hh"
@@ -33,7 +33,7 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
using typename Base::LocalVectorType;
using FrictionCoupling = FrictionCouplingPair<GridType, LocalVectorType, field_type>;
- using ContactCoupling = Dune::Contact::DualMortarCoupling<field_type, GridType>;
+ using ContactCoupling = DualMortarCoupling<field_type, GridType>;
size_t bodyIndex(const size_t globalIdx) {
size_t i=0;
diff --git a/dune/tectonic/data-structures/network/contactnetwork.cc b/dune/tectonic/data-structures/network/contactnetwork.cc
index e369614c..d9779a74 100644
--- a/dune/tectonic/data-structures/network/contactnetwork.cc
+++ b/dune/tectonic/data-structures/network/contactnetwork.cc
@@ -26,7 +26,7 @@ ContactNetwork<HostGridType, VectorType>::ContactNetwork(
couplings_(nCouplings),
frictionBoundaries_(nBodies),
stateEnergyNorms_(nBodies),
- nBodyAssembler_(nBodies, nCouplings, false, 0.0)
+ nBodyAssembler_(nBodies, nCouplings, 0.0)
{}
template <class HostGridType, class VectorType>
@@ -79,6 +79,28 @@ void ContactNetwork<HostGridType, VectorType>::assemble() {
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();
@@ -263,14 +285,14 @@ auto ContactNetwork<HostGridType, VectorType>::couplings() const
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::nBodyAssembler()
--> Dune::Contact::NBodyAssembler<GridType, VectorType>& {
+-> NBodyAssembler& {
return nBodyAssembler_;
}
template <class HostGridType, class VectorType>
auto ContactNetwork<HostGridType, VectorType>::nBodyAssembler() const
--> const Dune::Contact::NBodyAssembler<GridType, VectorType>& {
+-> const NBodyAssembler& {
return nBodyAssembler_;
}
diff --git a/dune/tectonic/data-structures/network/contactnetwork.hh b/dune/tectonic/data-structures/network/contactnetwork.hh
index fc9f554c..87b01365 100644
--- a/dune/tectonic/data-structures/network/contactnetwork.hh
+++ b/dune/tectonic/data-structures/network/contactnetwork.hh
@@ -5,8 +5,7 @@
#include <dune/common/promotiontraits.hh>
-#include <dune/contact/assemblers/nbodyassembler.hh>
-
+#include "../../spatial-solving/contact/nbodyassembler.hh"
#include "../../assemblers.hh"
#include "../enums.hh"
#include "../matrices.hh"
@@ -54,7 +53,7 @@ class ContactNetwork {
using StateEnergyNorms = std::vector<std::unique_ptr<typename LeafBody::StateEnergyNorm> >;
using ExternalForces = std::vector<std::unique_ptr<const typename LeafBody::ExternalForce> >;
- using NBodyAssembler = Dune::Contact::NBodyAssembler<GridType, VectorType>;
+ using NBodyAssembler = NBodyAssembler<GridType, VectorType>;
using LevelContactNetwork = LevelContactNetwork<GridType, FrictionCouplingPair, field_type>;
using Body = typename LevelContactNetwork::Body;
diff --git a/dune/tectonic/data-structures/program_state.hh b/dune/tectonic/data-structures/program_state.hh
index 05198b41..9e1f63f5 100644
--- a/dune/tectonic/data-structures/program_state.hh
+++ b/dune/tectonic/data-structures/program_state.hh
@@ -29,7 +29,7 @@
#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
-#include "../spatial-solving/preconditioners/nbodycontacttransfer.hh"
+#include "../spatial-solving/contact/nbodycontacttransfer.hh"
template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class BodyState {
public:
diff --git a/dune/tectonic/factories/stackedblocksfactory.hh b/dune/tectonic/factories/stackedblocksfactory.hh
index e8f9454f..fb21774a 100644
--- a/dune/tectonic/factories/stackedblocksfactory.hh
+++ b/dune/tectonic/factories/stackedblocksfactory.hh
@@ -54,7 +54,7 @@ template <class HostGridType, class VectorType> class StackedBlocksFactory : pub
public:
StackedBlocksFactory(const Dune::ParameterTree& parset) :
Base(parset, parset.get<size_t>("problem.bodyCount"), parset.get<size_t>("problem.bodyCount")-1),
- bodyData_(std::make_shared<MyBodyData<dim>>(this->parset_, zenith_())),
+ bodyData_(std::make_shared<MyBodyData<dim>>(this->parset_.sub("body"), this->parset_.template get<double>("gravity"), zenith_())),
cuboidGeometries_(this->bodyCount_),
leafFaces_(this->bodyCount_),
levelFaces_(this->bodyCount_),
diff --git a/dune/tectonic/factories/threeblocksfactory.cc b/dune/tectonic/factories/threeblocksfactory.cc
index 61275c63..6c395565 100644
--- a/dune/tectonic/factories/threeblocksfactory.cc
+++ b/dune/tectonic/factories/threeblocksfactory.cc
@@ -116,10 +116,10 @@ void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setCouplings() {
std::array<std::array<size_t, 2>, 3> couplingIndices;
- nonmortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->upper);
- mortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->lower);
- weakPatches_[0] = std::make_shared<WeakeningRegion>(cuboidGeometries_[0]->weakeningRegions()[0]);
- couplingIndices[0] = {0, 1};
+ nonmortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->lower);
+ mortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->upper);
+ weakPatches_[0] = std::make_shared<WeakeningRegion>(cuboidGeometries_[1]->weakeningRegions()[0]);
+ couplingIndices[0] = {1, 0};
nonmortarPatch_[1] = std::make_shared<LevelBoundaryPatch>(levelFaces_[2]->lower);
mortarPatch_[1] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->upper);
diff --git a/dune/tectonic/factories/threeblocksfactory.hh b/dune/tectonic/factories/threeblocksfactory.hh
index e72c7f1c..b846f9f3 100644
--- a/dune/tectonic/factories/threeblocksfactory.hh
+++ b/dune/tectonic/factories/threeblocksfactory.hh
@@ -55,7 +55,7 @@ class ThreeBlocksFactory : public ContactNetworkFactory<HostGridType, VectorType
public:
ThreeBlocksFactory(const Dune::ParameterTree& parset) :
Base(parset, 3, 3),
- bodyData_(std::make_shared<MyBodyData<dim>>(this->parset_, zenith_())),
+ bodyData_(std::make_shared<MyBodyData<dim>>(this->parset_.sub("body"), this->parset_.template get<double>("gravity"), zenith_())),
cuboidGeometries_(3),
leafFaces_(3),
levelFaces_(3),
diff --git a/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc
index 9f4daa82..7fb4bd30 100644
--- a/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc
+++ b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc
@@ -19,62 +19,11 @@
#include <dune/contact/common/dualbasisadapter.hh>
#include <dune/matrix-vector/addtodiagonal.hh>
-namespace Dune {
-namespace Contact {
-
-template <class field_type, class GridType0, class GridType1>
-void DualMortarCoupling<field_type, GridType0, GridType1>::assembleNonmortarLagrangeMatrix()
-{
- // Create mapping from the global set of block dofs to the ones on the contact boundary
- std::vector<int> globalToLocal;
- nonmortarBoundary_.makeGlobalToLocal(globalToLocal);
-
- // clear matrix
- nonmortarLagrangeMatrix_ = Dune::BDMatrix<MatrixBlock>(nonmortarBoundary_.numVertices());
- nonmortarLagrangeMatrix_ = 0;
-
- const auto& indexSet = grid0_->leafIndexSet();
-
- // loop over all faces of the boundary patch
- for (const auto& nIt : nonmortarBoundary_) {
-
- const auto& geometry = nIt.geometry();
- const field_type numCorners = geometry.corners();
-
- ctype intElem = geometry.integrationElement(Dune::FieldVector<ctype,dim-1>(0));
-
- field_type sfI = (numCorners==3) ? intElem/6.0 : intElem/numCorners;
-
- // turn scalar element mass matrix into vector-valued one
- // and add element mass matrix to the global matrix
-
- // Get global node ids
- const auto& inside = nIt.inside();
- const auto& refElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
-
- // Add localMatrix to nonmortarLagrangeMat
- for (int i=0; i<refElement.size(nIt.indexInInside(),1,dim); i++) {
- // we can use subEntity here because we add all indices anyway
- int v = globalToLocal[indexSet.subIndex(inside,refElement.subEntity(nIt.indexInInside(),1,
- i,dim),dim)];
- Dune::MatrixVector::addToDiagonal(nonmortarLagrangeMatrix_[v][v],sfI);
- }
-
- }
-
-}
-
template <class field_type, class GridType0, class GridType1>
-void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
-{
+void DualMortarCoupling<field_type, GridType0, GridType1>::reducePatches() {
typedef typename GridType0::LeafGridView GridView0;
typedef typename GridType1::LeafGridView GridView1;
- typedef Dune::PQkLocalFiniteElementCache<typename GridType1::ctype, field_type, GridType1::dimension,1> FiniteElementCache1;
-
- // cache for the dual functions on the boundary
- using DualCache = Dune::Contact::DualBasisAdapter<GridView0, field_type>;
-
using Element0 = typename GridView0::template Codim<0>::Entity;
using Element1 = typename GridView1::template Codim<0>::Entity;
@@ -150,13 +99,95 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
for (const auto& rIs : intersections(glue))
if (nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
mortarBoundary_.addFace(rIs.outside(),rIs.indexInOutside());
+}
+
+template <class field_type, class GridType0, class GridType1>
+void DualMortarCoupling<field_type, GridType0, GridType1>::setCrosspoints(const std::set<size_t>& allCrosspoints) {
+ crosspoints_.clear();
+
+ if (!allCrosspoints.empty()) {
+ Dune::BitSetVector<1> nonmortarPatchBoundary;
+ nonmortarBoundary_.getPatchBoundaryVertices(nonmortarPatchBoundary);
+
+ for (auto p : allCrosspoints) {
+ if (nonmortarPatchBoundary[p][0]) {
+ crosspoints_.emplace(p);
+ }
+ }
+ }
+}
+
+template <class field_type, class GridType0, class GridType1>
+void DualMortarCoupling<field_type, GridType0, GridType1>::assembleNonmortarLagrangeMatrix()
+{
+ // clear matrix
+ const auto nNonmortarVertices = nonmortarBoundary_.numVertices() - crosspoints_.size();
+ nonmortarLagrangeMatrix_ = Dune::BDMatrix<MatrixBlock>(nNonmortarVertices);
+ nonmortarLagrangeMatrix_ = 0;
+
+ const auto& indexSet = grid0_->leafIndexSet();
+
+ // loop over all faces of the boundary patch
+ for (const auto& nIt : nonmortarBoundary_) {
+
+ const auto& geometry = nIt.geometry();
+ const field_type numCorners = geometry.corners();
+
+ ctype intElem = geometry.integrationElement(Dune::FieldVector<ctype,dim-1>(0));
+
+ field_type sfI = (numCorners==3) ? intElem/6.0 : intElem/numCorners;
+
+ // turn scalar element mass matrix into vector-valued one
+ // and add element mass matrix to the global matrix
+
+ // Get global node ids
+ const auto& inside = nIt.inside();
+ const auto& refElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+
+ // Add localMatrix to nonmortarLagrangeMat
+ for (int i=0; i<refElement.size(nIt.indexInInside(),1,dim); i++) {
+ // we can use subEntity here because we add all indices anyway
+ int v = globalToLocal_[indexSet.subIndex(inside,refElement.subEntity(nIt.indexInInside(),1,
+ i,dim),dim)];
+ Dune::MatrixVector::addToDiagonal(nonmortarLagrangeMatrix_[v][v],sfI);
+ }
+
+ }
+}
+
+template <class field_type, class GridType0, class GridType1>
+void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
+{
+ typedef typename GridType0::LeafGridView GridView0;
+ typedef typename GridType1::LeafGridView GridView1;
+
+ using FiniteElementCache1 = Dune::PQkLocalFiniteElementCache<typename GridType1::ctype, field_type, GridType1::dimension,1>;
+
+ // cache for the dual functions on the boundary
+ using DualCache = Dune::Contact::DualBasisAdapter<GridView0, field_type>;
+ using ConstantDualCache = Dune::PQkLocalFiniteElementCache<typename GridType0::ctype, field_type, GridType0::dimension, 0>;
+
+ GridView0 gridView0 = grid0_->leafGridView();
+ GridView1 gridView1 = grid1_->leafGridView();
+
+ const auto& indexSet0 = gridView0.indexSet();
+ const auto& indexSet1 = gridView1.indexSet();
+
+ // Create mapping from the global set of block dofs to the ones on the contact boundary
+ const auto& vertices = *nonmortarBoundary_.getVertices();
+
+ globalToLocal_.resize(vertices.size());
+ int idx = 0;
+ for (size_t i=0; i<globalToLocal_.size(); i++)
+ globalToLocal_[i] = (vertices[i][0] and crosspoints_.count(i)==0) ? idx++ : -1;
// Assemble the diagonal matrix coupling the nonmortar side with the lagrange multiplyers there
assembleNonmortarLagrangeMatrix();
// The weak obstacle vector
- weakObstacle_.resize(nonmortarBoundary_.numVertices());
+ const auto nNonmortarVertices = nonmortarBoundary_.numVertices()-crosspoints_.size();
+ weakObstacle_.resize(nNonmortarVertices);
weakObstacle_ = 0;
// ///////////////////////////////////////////////////////////
@@ -164,14 +195,10 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
// ///////////////////////////////////////////////////////////
/** \todo Also restrict mortar indices and don't use the whole grid level. */
- Dune::MatrixIndexSet mortarIndices(nonmortarBoundary_.numVertices(), grid1_->size(dim));
-
- // Create mapping from the global set of block dofs to the ones on the contact boundary
- std::vector<int> globalToLocal;
- nonmortarBoundary_.makeGlobalToLocal(globalToLocal);
+ Dune::MatrixIndexSet mortarIndices(nNonmortarVertices, grid1_->size(dim));
// loop over all intersections
- for (const auto& rIs : intersections(glue)) {
+ for (const auto& rIs : intersections(*glue_)) {
if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
continue;
@@ -187,9 +214,9 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
for (int j=0; j<nDomainVertices; j++) {
- int localDomainIdx = globalToLocal[indexSet0.subIndex(inside,j,dim)];
+ int localDomainIdx = globalToLocal_[indexSet0.subIndex(inside,j,dim)];
- // if the vertex is not contained in the restricted contact boundary then dismiss it
+ // if the vertex is not contained in the restricted contact boundary or a crosspoint, then dismiss it
if (localDomainIdx == -1)
continue;
@@ -208,18 +235,20 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
// Clear it
mortarLagrangeMatrix_ = 0;
-
//cache of local bases
FiniteElementCache1 cache1;
std::unique_ptr<DualCache> dualCache;
dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView0, field_type> >();
+ // adapt dual basis close to crosspoints using ConstantDualCache
+ ConstantDualCache constantDualCache;
+
std::vector<Dune::FieldVector<ctype,dim> > avNormals;
avNormals = nonmortarBoundary_.getNormals();
// loop over all intersections and compute the matrix entries
- for (const auto& rIs : intersections(glue)) {
+ for (const auto& rIs : intersections(*glue_)) {
const auto& inside = rIs.inside();
@@ -250,6 +279,8 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
//const typename FiniteElementCache0::FiniteElementType& nonmortarFiniteElement = cache0.get(nonmortarEType);
const auto& mortarFiniteElement = cache1.get(mortarEType);
+
+ const auto& constantDualFE = constantDualCache.get(nonmortarEType);
dualCache->bind(inside, rIs.indexInInside());
std::vector<Dune::FieldVector<field_type,1> > mortarQuadValues, dualQuadValues;
@@ -259,11 +290,18 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
const auto& rGeomInInside = rIs.geometryInInside();
const auto& rGeomInOutside = rIs.geometryInOutside();
+ bool isCrosspointFace = false;
int nNonmortarFaceNodes = domainRefElement.size(rIs.indexInInside(),1,dim);
std::vector<int> nonmortarFaceNodes;
for (int i=0; i<nNonmortarFaceNodes; i++) {
int faceIdxi = domainRefElement.subEntity(rIs.indexInInside(), 1, i, dim);
- nonmortarFaceNodes.push_back(faceIdxi);
+ int globalIdx = indexSet0.subIndex(inside,faceIdxi,dim);
+
+ if (crosspoints_.count(globalIdx)>0) {
+ isCrosspointFace = true;
+ } else {
+ nonmortarFaceNodes.push_back(faceIdxi);
+ }
}
for (const auto& quadPt : quadRule) {
@@ -285,13 +323,18 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
//evaluate all shapefunctions at the quadrature point
//nonmortarFiniteElement.localBasis().evaluateFunction(nonmortarQuadPos,nonmortarQuadValues);
mortarFiniteElement.localBasis().evaluateFunction(mortarQuadPos,mortarQuadValues);
- dualCache->evaluateFunction(nonmortarQuadPos,dualQuadValues);
+
+ if (isCrosspointFace) {
+ constantDualFE.localBasis().evaluateFunction(nonmortarQuadPos,dualQuadValues);
+ } else {
+ dualCache->evaluateFunction(nonmortarQuadPos,dualQuadValues);
+ }
// loop over all Lagrange multiplier shape functions
for (int j=0; j<nNonmortarFaceNodes; j++) {
int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
- int rowIdx = globalToLocal[globalDomainIdx];
+ int rowIdx = globalToLocal_[globalDomainIdx];
weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
* dualQuadValues[nonmortarFaceNodes[j]] * (gapVector*avNormals[globalDomainIdx]);
@@ -339,6 +382,3 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
gridGlueBackend_->clear();
}
-
-} /* namespace Contact */
-} /* namespace Dune */
diff --git a/dune/tectonic/spatial-solving/contact/dualmortarcoupling.hh b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.hh
index cfd8ed5f..04782aab 100644
--- a/dune/tectonic/spatial-solving/contact/dualmortarcoupling.hh
+++ b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.hh
@@ -1,7 +1,7 @@
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
-#ifndef DUNE_CONTACT_ASSEMBLERS_DUAL_MORTAR_COUPLING_HH
-#define DUNE_CONTACT_ASSEMBLERS_DUAL_MORTAR_COUPLING_HH
+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_CONTACT_DUAL_MORTAR_COUPLING_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_CONTACT_DUAL_MORTAR_COUPLING_HH
#include <memory>
@@ -16,10 +16,9 @@
#include <dune/grid-glue/extractors/codim1extractor.hh>
#include <dune/grid-glue/merging/merger.hh>
-namespace Dune {
-namespace Contact {
-
-/** \brief Assembles the transfer operator for two-body contact
+/** \brief Assembles the transfer operator for two-body contact.
+ * Extension of Dune::Contact::DualMortarCoupling by crosspoints.
+ * Note: so war only works in 2D!
*/
template<class field_type, class GridType0, class GridType1=GridType0>
class DualMortarCoupling {
@@ -66,6 +65,10 @@ class DualMortarCoupling {
/* Nothing. */
}
+ void reducePatches();
+
+ void setCrosspoints(const std::set<size_t>& allCrosspoints);
+
/** \brief Sets up the contact coupling operator on the grid leaf level */
virtual void setup();
@@ -83,7 +86,7 @@ class DualMortarCoupling {
}
/** \brief Get the non-mortar matrix. */
- const BDMatrix<MatrixBlock>& nonmortarLagrangeMatrix() const {
+ const Dune::BDMatrix<MatrixBlock>& nonmortarLagrangeMatrix() const {
return nonmortarLagrangeMatrix_;
}
@@ -141,6 +144,14 @@ class DualMortarCoupling {
return weakObstacle_;
}
+ const std::vector<int>& globalToLocal() const {
+ return globalToLocal_;
+ }
+
+ /* const auto& crosspoints() const {
+ return crosspoints_;
+ }*/
+
// /////////////////////////////////////////////
// Data members
// /////////////////////////////////////////////
@@ -158,6 +169,9 @@ class DualMortarCoupling {
/** \brief The crosspoints of the nonmortar body. */
std::set<size_t> crosspoints_;
+ /** \brief Nonmortar boundary patch global to local map. */
+ std::vector<int> globalToLocal_;
+
/** \brief The matrix coupling mortar side and Lagrange multipliers. */
MatrixType mortarLagrangeMatrix_;
@@ -184,9 +198,6 @@ class DualMortarCoupling {
field_type coveredArea_;
};
-} /* namespace Contact */
-} /* namespace Dune */
-
#include "dualmortarcoupling.cc"
#endif
diff --git a/dune/tectonic/spatial-solving/contact/nbodyassembler.cc b/dune/tectonic/spatial-solving/contact/nbodyassembler.cc
index 2d75582f..13cee857 100644
--- a/dune/tectonic/spatial-solving/contact/nbodyassembler.cc
+++ b/dune/tectonic/spatial-solving/contact/nbodyassembler.cc
@@ -10,8 +10,34 @@
#include <dune/fufem/boundarypatchprolongator.hh>
#include <dune/contact/projections/normalprojection.hh>
-namespace Dune {
-namespace Contact {
+
+#include "../../utils/debugutils.hh"
+
+template <class GridType, class VectorType>
+void NBodyAssembler<GridType, VectorType>::setCrosspoints() {
+ std::vector<Dune::BitSetVector<1>> bodywiseBoundaries(nGrids());
+ for (size_t i=0; i<nGrids(); i++) {
+ bodywiseBoundaries[i] = *dirichletVertices_[i];
+ }
+
+ for (size_t i=0; i<nCouplings(); i++) {
+ auto& coupling = coupling_[i];
+ const auto& contactCoupling = contactCoupling_[i]; // dual mortar object holding boundary patches
+ const auto nonmortarIdx = coupling.gridIdx_[0];
+ //const auto mortarIdx = coupling->gridIdx_[1];
+
+ auto& bodywiseBoundary = bodywiseBoundaries[nonmortarIdx];
+ auto& crosspoints = crosspoints_[nonmortarIdx];
+ const auto& nmBoundaryVertices = *contactCoupling->nonmortarBoundary().getVertices();
+ for (size_t j=0; j<nmBoundaryVertices.size(); j++) {
+ if (bodywiseBoundary[j][0] and nmBoundaryVertices[j][0]) {
+ crosspoints.emplace(j);
+ }
+ bodywiseBoundary[j][0] = bodywiseBoundary[j][0] or nmBoundaryVertices[j][0];
+ }
+
+ }
+}
template <class GridType, class VectorType>
void NBodyAssembler<GridType, VectorType>::assembleTransferOperator()
@@ -42,6 +68,17 @@ void NBodyAssembler<GridType, VectorType>::assembleTransferOperator()
contactCoupling_[i]->setupContactPatch(*coupling_[i].patch0(),*coupling_[i].patch1());
contactCoupling_[i]->gridGlueBackend_ = coupling_[i].backend();
contactCoupling_[i]->setCoveredArea(coveredArea_);
+ contactCoupling_[i]->reducePatches();
+ }
+
+ setCrosspoints();
+
+ for (size_t i=0; i<nGrids(); i++) {
+ print(crosspoints_[i], "crosspoints " + std::to_string(i));
+ }
+
+ for (size_t i=0; i<contactCoupling_.size(); i++) {
+ contactCoupling_[i]->setCrosspoints(crosspoints_[coupling_[i].gridIdx_[0]]);
contactCoupling_[i]->setup();
}
@@ -91,7 +128,6 @@ void NBodyAssembler<GridType, VectorType>::assembleHouseholderReflections()
this->localCoordSystems_[i] = id;
for (int j=0; j<nCouplings(); j++) {
-
int grid0Idx = coupling_[j].gridIdx_[0];
// compute offset
@@ -100,10 +136,11 @@ void NBodyAssembler<GridType, VectorType>::assembleHouseholderReflections()
idx += grids_[k]->size(dim);
const auto& nonmortarBoundary = contactCoupling_[j]->nonmortarBoundary();
+ const auto& globalToLocal = contactCoupling_[j]->globalToLocal();
- // if a body is non-mortar more then once, one has to be careful with not overwriting entries
+ // if a body is non-mortar more than once, one has to be careful with not overwriting entries
for (std::size_t k=0; k<coordinateSystems[j].size(); k++)
- if(nonmortarBoundary.containsVertex(k))
+ if(globalToLocal[k] > -1)
this->localCoordSystems_[idx+k] = coordinateSystems[j][k];
}
}
@@ -114,7 +151,7 @@ void NBodyAssembler<GridType, VectorType>::assembleObstacle()
std::vector<std::vector<int> > globalToLocals(nCouplings());
for (std::size_t i = 0; i < globalToLocals.size(); ++i)
- contactCoupling_[i]->nonmortarBoundary().makeGlobalToLocal(globalToLocals[i]);
+ globalToLocals[i] = contactCoupling_[i]->globalToLocal();
///////////////////////////////////////
// Set the obstacle values
@@ -126,6 +163,7 @@ void NBodyAssembler<GridType, VectorType>::assembleObstacle()
int grid0Idx = coupling_[j].gridIdx_[0];
const auto& nonmortarBoundary = contactCoupling_[j]->nonmortarBoundary();
+ const auto& globalToLocal = globalToLocals[j];
// compute offset
int idx = 0;
@@ -133,7 +171,7 @@ void NBodyAssembler<GridType, VectorType>::assembleObstacle()
idx += grids_[k]->size(dim);
for (int k=0; k<grids_[grid0Idx]->size(dim); k++)
- if (nonmortarBoundary.containsVertex(k))
+ if (globalToLocal[k] > -1)
totalHasObstacle_[idx+k] = true;
}
@@ -181,23 +219,23 @@ void NBodyAssembler<GridType, VectorType>::assembleObstacle()
// Set the obstacle
switch (coupling_[i].type_) {
- case CouplingPairBase::CONTACT:
+ case Dune::Contact::CouplingPairBase::CONTACT:
totalObstacles_[idx+vIdx].upper(0) = obs[globalToLocal[vIdx]];
break;
- case CouplingPairBase::STICK_SLIP:
+ case Dune::Contact::CouplingPairBase::STICK_SLIP:
totalObstacles_[idx+vIdx].lower(0) = 0;
totalObstacles_[idx+vIdx].upper(0) = 0;
break;
- case CouplingPairBase::GLUE:
+ case Dune::Contact::CouplingPairBase::GLUE:
for (int j=0; j<dim; j++) {
totalObstacles_[idx+vIdx].lower(j) = 0;
totalObstacles_[idx+vIdx].upper(j) = 0;
}
break;
- case CouplingPairBase::NONE:
+ case Dune::Contact::CouplingPairBase::NONE:
break;
default:
DUNE_THROW(Dune::NotImplemented, "Coupling type not handled yet!");
@@ -438,6 +476,7 @@ computeTransformationMatrix()
for (int i=0; i<nCouplings(); i++) {
const auto& nonmortarBoundary = contactCoupling_[i]->nonmortarBoundary();
+ const auto& globalToLocal = contactCoupling_[i]->globalToLocal();
// If the contact mapping could not be built at all then skip this coupling
if (nonmortarBoundary.numVertices() == 0)
continue;
@@ -450,8 +489,9 @@ computeTransformationMatrix()
nonmortarToGlobal[i].resize(nonmortarBoundary.numVertices());
int idx = 0;
for (int j=0; j<grids_[grid0Idx]->size(dim); j++)
- if (nonmortarBoundary.containsVertex(j))
+ if (globalToLocal[j] > -1)
nonmortarToGlobal[i][idx++] = j;
+ nonmortarToGlobal[i].resize(idx);
// the off-diagonal part
const MatrixType& M = contactCoupling_[i]->mortarLagrangeMatrix();
@@ -485,6 +525,7 @@ computeTransformationMatrix()
for (int i=0; i<nCouplings(); i++) {
const auto& nonmortarBoundary = contactCoupling_[i]->nonmortarBoundary();
+ const auto& globalToLocal = contactCoupling_[i]->globalToLocal();
if (nonmortarBoundary.numVertices() == 0)
continue;
@@ -510,10 +551,7 @@ computeTransformationMatrix()
// modify non-mortar dofs and rotate them in normal and tangential coordinates
for (int k=0; k<grids_[grid0Idx]->size(dim); k++)
- if(nonmortarBoundary.containsVertex(k))
+ if(globalToLocal[k] > -1)
BT_[offset + k][offset+k] =this->localCoordSystems_[offset + k];
}
}
-
-} /* namespace Contact */
-} /* namespace Dune */
diff --git a/dune/tectonic/spatial-solving/contact/nbodyassembler.hh b/dune/tectonic/spatial-solving/contact/nbodyassembler.hh
index 55a5d633..e197497e 100644
--- a/dune/tectonic/spatial-solving/contact/nbodyassembler.hh
+++ b/dune/tectonic/spatial-solving/contact/nbodyassembler.hh
@@ -1,7 +1,7 @@
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
-#ifndef DUNE_CONTACT_ASSEMBLERS_N_BODY_MORTAR_ASSEMBLER_HH
-#define DUNE_CONTACT_ASSEMBLERS_N_BODY_MORTAR_ASSEMBLER_HH
+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_CONTACT_N_BODY_MORTAR_ASSEMBLER_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_CONTACT_N_BODY_MORTAR_ASSEMBLER_HH
#include <dune/common/bitsetvector.hh>
#include <dune/common/promotiontraits.hh>
@@ -9,18 +9,16 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/solvers/common/boxconstraint.hh>
#include <dune/fufem/boundarypatch.hh>
-#include <dune/contact/assemblers/dualmortarcouplinghierarchy.hh>
-#include <dune/contact/assemblers/dualmortarcoupling.hh>
#include <dune/contact/assemblers/contactassembler.hh>
#include <dune/contact/common/couplingpair.hh>
#include <dune/contact/projections/contactprojection.hh>
+#include "dualmortarcoupling.hh"
+
#ifdef HAVE_DUNE_GRID_GLUE
#include <dune/grid-glue/merging/merger.hh>
#endif
-namespace Dune {
-namespace Contact {
/** \brief Assembler for mortar discretized contact problems with arbitrary number of bodies.
*
@@ -28,16 +26,15 @@ namespace Contact {
* \tparam VectorType The vector type for the displacements.
*/
template <class GridType, class VectorType>
-class NBodyAssembler : public ContactAssembler<GridType::dimension, typename VectorType::field_type>
+class NBodyAssembler : public Dune::Contact::ContactAssembler<GridType::dimension, typename VectorType::field_type>
{
protected:
enum {dim = GridType::dimension};
- using field_type = typename PromotionTraits<typename VectorType::field_type,
+ using field_type = typename Dune::PromotionTraits<typename VectorType::field_type,
typename GridType::ctype>::PromotedType;
using CouplingType = DualMortarCoupling<field_type, GridType>;
- using HierarchyCouplingType = DualMortarCouplingHierarchy<field_type, GridType>;
using MatrixBlock = Dune::FieldMatrix<field_type, dim, dim>;
using MatrixType = Dune::BCRSMatrix<MatrixBlock>;
@@ -52,21 +49,19 @@ class NBodyAssembler : public ContactAssembler<GridType::dimension, typename Vec
*
* \param nGrids The number of involved bodies.
* \param nCouplings The number of couplings.
- * \param hierarchyCoupling This boolean determines if the whole hierarchy of mortar matrices should be setup
*/
- NBodyAssembler(int nGrids, int nCouplings, bool hierarchyCoupling=false, field_type coveredArea = 0.99) :
+ NBodyAssembler(int nGrids, int nCouplings, field_type coveredArea = 0.99) :
coveredArea_(coveredArea)
{
grids_.resize(nGrids);
+ crosspoints_.resize(nGrids);
+
coupling_.resize(nCouplings);
contactCoupling_.resize(nCouplings);
// initialise the couplings with DualMortarCoupling
for (int i=0; i<nCouplings; i++)
- if (hierarchyCoupling)
- contactCoupling_[i] = std::make_shared<HierarchyCouplingType>();
- else
- contactCoupling_[i] = std::make_shared<CouplingType>();
+ contactCoupling_[i] = std::make_shared<CouplingType>();
}
/** \brief Get the number of couplings.*/
@@ -83,6 +78,8 @@ class NBodyAssembler : public ContactAssembler<GridType::dimension, typename Vec
return n;
}
+ void setCrosspoints();
+
/** \brief Setup the patches for the contact boundary and compute the obstacles. */
void assembleObstacle();
@@ -159,13 +156,16 @@ class NBodyAssembler : public ContactAssembler<GridType::dimension, typename Vec
void setGrids(const std::vector<const GridType*>& grids) {grids_ = grids;}
/** \brief Set grids. */
- void setCoupling(const CouplingPair<GridType, GridType, field_type>& coupling, size_t i) {
+ void setCoupling(const Dune::Contact::CouplingPair<GridType, GridType, field_type>& coupling, size_t i) {
coupling_[i] = coupling;
}
/** \brief Get reference to i'th coupling. */
const auto& getCoupling(size_t i) const {return coupling_[i];}
+ void setDirichletVertices(std::vector<std::shared_ptr<Dune::BitSetVector<1>>> dirichletVertices) {
+ dirichletVertices_ = dirichletVertices;
+ }
protected:
/** \brief Compute the transposed mortar transformation matrix. */
void computeTransformationMatrix();
@@ -174,17 +174,21 @@ class NBodyAssembler : public ContactAssembler<GridType::dimension, typename Vec
void assembleHouseholderReflections();
public:
+ std::vector<std::shared_ptr<Dune::BitSetVector<1>>> dirichletVertices_;
+
/** \brief Bitvector for all bodies whose flags are set if a dof has an obstacle.*/
Dune::BitSetVector<dim> totalHasObstacle_;
/** \brief Obstacles for all bodies on the leafview.*/
std::vector<BoxConstraint<field_type,dim> > totalObstacles_;
+ std::vector<std::set<size_t>> crosspoints_;
+
/** \brief The mortar couplings between grid pairs */
std::vector<std::shared_ptr<CouplingType> > contactCoupling_;
/** \brief The coupling pairs. */
- std::vector<CouplingPair<GridType,GridType,field_type> > coupling_;
+ std::vector<Dune::Contact::CouplingPair<GridType,GridType,field_type> > coupling_;
/** \brief Vector containing the involved grids. */
std::vector<const GridType*> grids_;
@@ -198,9 +202,6 @@ class NBodyAssembler : public ContactAssembler<GridType::dimension, typename Vec
field_type coveredArea_;
};
-} /* namespace Contact */
-} /* namespace Dune */
-
#include "nbodyassembler.cc"
#endif
diff --git a/dune/tectonic/spatial-solving/fixedpointiterator.cc b/dune/tectonic/spatial-solving/fixedpointiterator.cc
index 298e2f94..39610f48 100644
--- a/dune/tectonic/spatial-solving/fixedpointiterator.cc
+++ b/dune/tectonic/spatial-solving/fixedpointiterator.cc
@@ -34,7 +34,7 @@
#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
-#include "../spatial-solving/preconditioners/nbodycontacttransfer.hh"
+#include "../spatial-solving/contact/nbodycontacttransfer.hh"
#include "tnnmg/functional.hh"
#include "tnnmg/zerononlinearity.hh"
diff --git a/dune/tectonic/spatial-solving/fixedpointiterator_tmpl.cc b/dune/tectonic/spatial-solving/fixedpointiterator_tmpl.cc
index bae2cc32..c434a9b2 100644
--- a/dune/tectonic/spatial-solving/fixedpointiterator_tmpl.cc
+++ b/dune/tectonic/spatial-solving/fixedpointiterator_tmpl.cc
@@ -29,13 +29,13 @@ using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using LinearSolver = Dune::Solvers::LoopSolver<Vector>;
using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
-using NBodyAssembler = typename MyContactNetwork::NBodyAssembler;
+using MyNBodyAssembler = typename MyContactNetwork::NBodyAssembler;
using MyGlobalFriction = GlobalFriction<Matrix, Vector>;
using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Vector&, double>;
using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
-template class FixedPointIterator<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>;
+template class FixedPointIterator<MySolverFactory, MyNBodyAssembler, MyUpdaters, ErrorNorms>;
-template FixedPointIterationCounter FixedPointIterator<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>::run<LinearSolver>(
+template FixedPointIterationCounter FixedPointIterator<MySolverFactory, MyNBodyAssembler, MyUpdaters, ErrorNorms>::run<LinearSolver>(
MyUpdaters, std::shared_ptr<LinearSolver>&, const std::vector<Matrix>&, const std::vector<Vector>&, std::vector<Vector>&);
diff --git a/dune/tectonic/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh b/dune/tectonic/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
index 88de3c32..e67dfe55 100644
--- a/dune/tectonic/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
+++ b/dune/tectonic/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
@@ -15,7 +15,7 @@
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
-#include "nbodycontacttransfer.hh"
+#include "../contact/nbodycontacttransfer.hh"
#include "levelpatchpreconditioner.hh"
template <class ContactNetwork, class MatrixType, class VectorType>
diff --git a/dune/tectonic/time-stepping/adaptivetimestepper_tmpl.cc b/dune/tectonic/time-stepping/adaptivetimestepper_tmpl.cc
index 0e7cd4f9..64155c96 100644
--- a/dune/tectonic/time-stepping/adaptivetimestepper_tmpl.cc
+++ b/dune/tectonic/time-stepping/adaptivetimestepper_tmpl.cc
@@ -29,7 +29,7 @@ using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using LinearSolver = Dune::Solvers::LoopSolver<Vector>;
using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
-using NBodyAssembler = typename MyContactNetwork::NBodyAssembler;
+using MyNBodyAssembler = typename MyContactNetwork::NBodyAssembler;
using MyGlobalFriction = GlobalFriction<Matrix, Vector>;
using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Vector&, double>;
@@ -39,4 +39,4 @@ template class AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters
template typename AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>::UpdatersWithCount
AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>::step<LinearSolver>(
- const MyUpdaters&, const NBodyAssembler&, std::shared_ptr<LinearSolver>&, double, double);
+ const MyUpdaters&, const MyNBodyAssembler&, std::shared_ptr<LinearSolver>&, double, double);
diff --git a/dune/tectonic/time-stepping/coupledtimestepper_tmpl.cc b/dune/tectonic/time-stepping/coupledtimestepper_tmpl.cc
index e09f3f89..a2b7fc84 100644
--- a/dune/tectonic/time-stepping/coupledtimestepper_tmpl.cc
+++ b/dune/tectonic/time-stepping/coupledtimestepper_tmpl.cc
@@ -29,12 +29,12 @@ using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using LinearSolver = Dune::Solvers::LoopSolver<Vector>;
using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
-using NBodyAssembler = typename MyContactNetwork::NBodyAssembler;
+using MyNBodyAssembler = typename MyContactNetwork::NBodyAssembler;
using MyGlobalFriction = GlobalFriction<Matrix, Vector>;
using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Vector&, double>;
using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
-template class CoupledTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>;
+template class CoupledTimeStepper<MySolverFactory, MyNBodyAssembler, MyUpdaters, ErrorNorms>;
-template FixedPointIterationCounter CoupledTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>::step<LinearSolver>(std::shared_ptr<LinearSolver>&, double, double);
+template FixedPointIterationCounter CoupledTimeStepper<MySolverFactory, MyNBodyAssembler, MyUpdaters, ErrorNorms>::step<LinearSolver>(std::shared_ptr<LinearSolver>&, double, double);
diff --git a/dune/tectonic/time-stepping/state_tmpl.cc b/dune/tectonic/time-stepping/state_tmpl.cc
index d4d99038..89834f69 100644
--- a/dune/tectonic/time-stepping/state_tmpl.cc
+++ b/dune/tectonic/time-stepping/state_tmpl.cc
@@ -5,11 +5,12 @@
#include <dune/contact/assemblers/dualmortarcoupling.hh>
#include "../data-structures/friction/frictioncouplingpair.hh"
+#include "../spatial-solving/contact/dualmortarcoupling.hh"
using field_type = typename Dune::PromotionTraits<typename Vector::field_type,
typename DeformedGrid::ctype>::PromotedType;
-using MyContactCoupling = Dune::Contact::DualMortarCoupling<field_type, DeformedGrid>;
+using MyContactCoupling = DualMortarCoupling<field_type, DeformedGrid>;
using MyFrictionCouplingPair = FrictionCouplingPair<DeformedGrid, LocalVector, field_type>;
template std::shared_ptr<StateUpdater<ScalarVector, Vector>>
diff --git a/src/multi-body-problem/multi-body-problem.cc b/src/multi-body-problem/multi-body-problem.cc
index 715acca3..30fc02b3 100644
--- a/src/multi-body-problem/multi-body-problem.cc
+++ b/src/multi-body-problem/multi-body-problem.cc
@@ -110,7 +110,7 @@ int main(int argc, char *argv[]) {
std::cout << argv[0] << std::endl;
}
- std::ofstream out("../log.txt");
+ std::ofstream out("multi-body-problem.log");
std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
std::cout.rdbuf(out.rdbuf()); //redirect std::cout to log.txt
@@ -132,7 +132,7 @@ int main(int argc, char *argv[]) {
using ContactNetwork = typename ThreeBlocksFactory<Grid, Vector>::ContactNetwork;
threeBlocksFactory.build();
- ContactNetwork& contactNetwork = threeBlocksFactory.contactNetwork(); */
+ ContactNetwork& contactNetwork = threeBlocksFactory.contactNetwork();*/
const size_t bodyCount = contactNetwork.nBodies();
@@ -224,6 +224,7 @@ int main(int argc, char *argv[]) {
std::vector<const MyCellBasis* > cellBases(bodyCount);
auto& wPatches = stackedBlocksFactory.weakPatches();
+ //auto& wPatches = threeBlocksFactory.weakPatches();
std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches(bodyCount);
--
GitLab