diff --git a/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc
index 7fb4bd3023aab974d967f8c373ec023f84ae0967..303e1ea0e8ae9fcc5245e7e4bd462c9feb90c027 100644
--- a/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc
+++ b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc
@@ -149,7 +149,10 @@ void DualMortarCoupling<field_type, GridType0, GridType1>::assembleNonmortarLagr
// 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);
+ if (v == -1)
+ continue;
+
+ Dune::MatrixVector::addToDiagonal(nonmortarLagrangeMatrix_[v][v],sfI);
}
}
@@ -178,9 +181,13 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
const auto& vertices = *nonmortarBoundary_.getVertices();
globalToLocal_.resize(vertices.size());
+ hasObstacle_.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;
+ for (size_t i=0; i<globalToLocal_.size(); i++) {
+ bool nonmortarVertex = vertices[i][0] and crosspoints_.count(i)==0;
+ globalToLocal_[i] = nonmortarVertex ? idx++ : -1;
+ hasObstacle_[i][0] = nonmortarVertex ? true : false;
+ }
// Assemble the diagonal matrix coupling the nonmortar side with the lagrange multiplyers there
assembleNonmortarLagrangeMatrix();
@@ -249,8 +256,6 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
// loop over all intersections and compute the matrix entries
for (const auto& rIs : intersections(*glue_)) {
-
-
const auto& inside = rIs.inside();
if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
@@ -291,9 +296,8 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
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++) {
+ for (int i=0; i<domainRefElement.size(rIs.indexInInside(),1,dim); i++) {
int faceIdxi = domainRefElement.subEntity(rIs.indexInInside(), 1, i, dim);
int globalIdx = indexSet0.subIndex(inside,faceIdxi,dim);
@@ -326,35 +330,59 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
if (isCrosspointFace) {
constantDualFE.localBasis().evaluateFunction(nonmortarQuadPos,dualQuadValues);
+
+ // loop over all Lagrange multiplier shape functions
+ for (size_t j=0; j<nonmortarFaceNodes.size(); j++) {
+
+ int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
+ int rowIdx = globalToLocal_[globalDomainIdx];
+
+ weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
+ * dualQuadValues[0] * (gapVector*avNormals[globalDomainIdx]);
+
+ // loop over all mortar shape functions
+ for (int k=0; k<noOfMortarVec; k++) {
+
+ int colIdx = indexSet1.subIndex(outside, k, dim);
+ if (!mortarBoundary_.containsVertex(colIdx))
+ continue;
+
+ // Integrate over the product of two shape functions
+ field_type mortarEntry = integrationElement* quadPt.weight()* dualQuadValues[0]* mortarQuadValues[k];
+
+ Dune::MatrixVector::addToDiagonal(mortarLagrangeMatrix_[rowIdx][colIdx], mortarEntry);
+
+ }
+
+ }
} else {
dualCache->evaluateFunction(nonmortarQuadPos,dualQuadValues);
- }
- // loop over all Lagrange multiplier shape functions
- for (int j=0; j<nNonmortarFaceNodes; j++) {
+ // loop over all Lagrange multiplier shape functions
+ for (size_t j=0; j<nonmortarFaceNodes.size(); j++) {
- int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
- int rowIdx = globalToLocal_[globalDomainIdx];
+ int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
+ int rowIdx = globalToLocal_[globalDomainIdx];
- weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
- * dualQuadValues[nonmortarFaceNodes[j]] * (gapVector*avNormals[globalDomainIdx]);
+ weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
+ * dualQuadValues[nonmortarFaceNodes[j]] * (gapVector*avNormals[globalDomainIdx]);
- // loop over all mortar shape functions
- for (int k=0; k<noOfMortarVec; k++) {
+ // loop over all mortar shape functions
+ for (int k=0; k<noOfMortarVec; k++) {
- int colIdx = indexSet1.subIndex(outside, k, dim);
- if (!mortarBoundary_.containsVertex(colIdx))
- continue;
+ int colIdx = indexSet1.subIndex(outside, k, dim);
+ if (!mortarBoundary_.containsVertex(colIdx))
+ continue;
- // Integrate over the product of two shape functions
- field_type mortarEntry = integrationElement* quadPt.weight()* dualQuadValues[nonmortarFaceNodes[j]]* mortarQuadValues[k];
+ // Integrate over the product of two shape functions
+ field_type mortarEntry = integrationElement* quadPt.weight()* dualQuadValues[nonmortarFaceNodes[j]]* mortarQuadValues[k];
- Dune::MatrixVector::addToDiagonal(mortarLagrangeMatrix_[rowIdx][colIdx], mortarEntry);
+ Dune::MatrixVector::addToDiagonal(mortarLagrangeMatrix_[rowIdx][colIdx], mortarEntry);
- }
+ }
+ }
}
-
}
}
diff --git a/dune/tectonic/spatial-solving/contact/dualmortarcoupling.hh b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.hh
index 04782aab79b773e05c2b8e0cbe7af6931becc0bc..daa93708d0e36ad5b01085b1819a369e4a8b191b 100644
--- a/dune/tectonic/spatial-solving/contact/dualmortarcoupling.hh
+++ b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.hh
@@ -144,6 +144,10 @@ class DualMortarCoupling {
return weakObstacle_;
}
+ const Dune::BitSetVector<1>* hasObstacle() const {
+ return &hasObstacle_;
+ }
+
const std::vector<int>& globalToLocal() const {
return globalToLocal_;
}
@@ -181,6 +185,8 @@ class DualMortarCoupling {
/** \brief The weak obstacles. */
ObstacleVector weakObstacle_;
+ Dune::BitSetVector<1> hasObstacle_;
+
/** \brief Allow some small penetration of the bodies. */
field_type overlap_;
diff --git a/dune/tectonic/spatial-solving/contact/nbodyassembler.cc b/dune/tectonic/spatial-solving/contact/nbodyassembler.cc
index 13cee85783346c57bc2f0da1ba5c53fc632538e6..0c0d7ca9fdad7009f979fcfb2c4a1f5dafd30226 100644
--- a/dune/tectonic/spatial-solving/contact/nbodyassembler.cc
+++ b/dune/tectonic/spatial-solving/contact/nbodyassembler.cc
@@ -15,7 +15,7 @@
template <class GridType, class VectorType>
void NBodyAssembler<GridType, VectorType>::setCrosspoints() {
- std::vector<Dune::BitSetVector<1>> bodywiseBoundaries(nGrids());
+ /* std::vector<Dune::BitSetVector<1>> bodywiseBoundaries(nGrids());
for (size_t i=0; i<nGrids(); i++) {
bodywiseBoundaries[i] = *dirichletVertices_[i];
}
@@ -36,6 +36,29 @@ void NBodyAssembler<GridType, VectorType>::setCrosspoints() {
bodywiseBoundary[j][0] = bodywiseBoundary[j][0] or nmBoundaryVertices[j][0];
}
+ }*/
+
+ for (int 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];
+
+ auto& crosspoints = crosspoints_[nonmortarIdx];
+
+ Dune::BitSetVector<1> nonmortarPatchBoundary;
+ contactCoupling->nonmortarBoundary().getPatchBoundaryVertices(nonmortarPatchBoundary);
+
+ for (size_t j=0; j<nonmortarPatchBoundary.size(); j++) {
+ if (nonmortarPatchBoundary[j][0]) {
+ crosspoints.emplace(j);
+ }
+ }
+
+ }
+
+
+ for (size_t i=0; i<crosspoints_.size(); i++) {
+ print(crosspoints_[i], "crosspoints " + std::to_string(i));
}
}
@@ -73,13 +96,16 @@ void NBodyAssembler<GridType, VectorType>::assembleTransferOperator()
setCrosspoints();
- for (size_t i=0; i<nGrids(); i++) {
+ for (int 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();
+
+ print(contactCoupling_[i]->nonmortarLagrangeMatrix() , "nonmortarLagrangeMatrix " + std::to_string(i) + ": ");
+ print(contactCoupling_[i]->mortarLagrangeMatrix() , "mortarLagrangeMatrix " + std::to_string(i) + ": ");
}
// setup Householder reflections
@@ -135,14 +161,15 @@ void NBodyAssembler<GridType, VectorType>::assembleHouseholderReflections()
for (int k=0;k<grid0Idx;k++)
idx += grids_[k]->size(dim);
- const auto& nonmortarBoundary = contactCoupling_[j]->nonmortarBoundary();
const auto& globalToLocal = contactCoupling_[j]->globalToLocal();
// 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(globalToLocal[k] > -1)
+ if(globalToLocal[k] > -1 ) //or crosspoints.count(k) > 0)
this->localCoordSystems_[idx+k] = coordinateSystems[j][k];
}
+
+ //print(this->localCoordSystems_, "localCoordSystems:");
}
template <class GridType, class VectorType>
@@ -162,7 +189,6 @@ void NBodyAssembler<GridType, VectorType>::assembleObstacle()
for (int j=0; j<nCouplings(); j++) {
int grid0Idx = coupling_[j].gridIdx_[0];
- const auto& nonmortarBoundary = contactCoupling_[j]->nonmortarBoundary();
const auto& globalToLocal = globalToLocals[j];
// compute offset
@@ -243,6 +269,11 @@ void NBodyAssembler<GridType, VectorType>::assembleObstacle()
}
}
+ //totalHasObstacle_[35] = true;
+ //totalObstacles_[35].lower(0) = 0;
+ // totalObstacles_[35].lower(1) = 0;
+ //totalObstacles_[35].upper(0) = 0;
+ // totalObstacles_[35].upper(1) = 0;
}
template <class GridType, class VectorType>
diff --git a/dune/tectonic/spatial-solving/contact/nbodycontacttransfer.cc b/dune/tectonic/spatial-solving/contact/nbodycontacttransfer.cc
index 347e2a22bd70770701aaf325d7a56f198129e505..0bdea11fa1c4be6ce142282e9064d769db5fc0ab 100644
--- a/dune/tectonic/spatial-solving/contact/nbodycontacttransfer.cc
+++ b/dune/tectonic/spatial-solving/contact/nbodycontacttransfer.cc
@@ -81,7 +81,7 @@ void NBodyContactTransfer<ContactNetwork, VectorType>::setup(const ContactNetwor
for (size_t i=0; i<nCouplings; i++) {
mortarTransferOperators[i] = &contactCouplings[i]->mortarLagrangeMatrix();
- fineHasObstacle[i] = contactCouplings[i]->nonmortarBoundary().getVertices();
+ fineHasObstacle[i] = contactCouplings[i]->hasObstacle();
gridIdx[i] = nBodyAssembler.getCoupling(i).gridIdx_;
}
diff --git a/src/foam/foam-2D.cfg b/src/foam/foam-2D.cfg
index c22b8a7224393ab58613ba0f4a075e142f412bfc..1d8c5160eb785e415ba26301883786d2b4465758 100644
--- a/src/foam/foam-2D.cfg
+++ b/src/foam/foam-2D.cfg
@@ -1,9 +1,9 @@
# -*- mode:conf -*-
[body0]
-smallestDiameter = 0.01 # 2e-3 [m]
+smallestDiameter = 0.02 # 2e-3 [m]
[body1]
-smallestDiameter = 0.01 # 2e-3 [m]
+smallestDiameter = 0.02 # 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-5 # 1e-5
diff --git a/src/foam/foam.cc b/src/foam/foam.cc
index ce0f8d3a9ab4f1c188ea71afbb1b4dc9853d3d48..9e4c7e56c6b7bb5df22dcf9fb1a3f304315f31c0 100644
--- a/src/foam/foam.cc
+++ b/src/foam/foam.cc
@@ -198,9 +198,6 @@ int main(int argc, char *argv[]) {
else
programState.setupInitialConditions(parset, contactNetwork);
-
- //DUNE_THROW(Dune::Exception, "Just need to stop here!");
-
auto& nBodyAssembler = contactNetwork.nBodyAssembler();
for (size_t i=0; i<bodyCount; i++) {
contactNetwork.body(i)->setDeformation(programState.u[i]);
@@ -343,6 +340,8 @@ int main(int argc, char *argv[]) {
print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
}
+ DUNE_THROW(Dune::Exception, "Just need to stop here!");
+
Updaters current(
initRateUpdater(
parset.get<Config::scheme>("timeSteps.scheme"),
diff --git a/src/foam/foam.cfg b/src/foam/foam.cfg
index ee78e2fc57ca31b4e0f2eca7f668fb5098664e29..c8604d621faee5052e095bf051016b3a884d133a 100644
--- a/src/foam/foam.cfg
+++ b/src/foam/foam.cfg
@@ -2,19 +2,19 @@
gravity = 9.81 # [m/s^2]
[body0]
-length = 0.04 # [m]
+length = 0.08 # [m]
height = 0.04 # [m]
depth = 0.04 # [m]
bulkModulus = 0.5e5 # [Pa] #2190
poissonRatio = 0.30 # [1] #0.11
[body0.elastic]
density = 900 # [kg/m^3] #750
-shearViscosity = 1e3 # [Pas]
-bulkViscosity = 1e3 # [Pas]
+shearViscosity = 1e20 # [Pas]
+bulkViscosity = 1e20 # [Pas]
[body0.viscoelastic]
density = 1000 # [kg/m^3]
-shearViscosity = 1e4 # [Pas]
-bulkViscosity = 1e4 # [Pas]
+shearViscosity = 1e20 # [Pas]
+bulkViscosity = 1e20 # [Pas]
[body1]
length = 0.04 # [m]
@@ -24,12 +24,12 @@ bulkModulus = 0.5e5 # [Pa]
poissonRatio = 0.30 # [1]
[body1.elastic]
density = 900 # [kg/m^3]
-shearViscosity = 1e3 # [Pas]
-bulkViscosity = 1e3 # [Pas]
+shearViscosity = 1e20 # [Pas]
+bulkViscosity = 1e20 # [Pas]
[body1.viscoelastic]
density = 1000 # [kg/m^3]
-shearViscosity = 1e4 # [Pas]
-bulkViscosity = 1e4 # [Pas]
+shearViscosity = 1e20 # [Pas]
+bulkViscosity = 1e20 # [Pas]
[boundary.friction]
C = 10 # [Pa]
diff --git a/src/multi-body-problem/multi-body-problem-2D.cfg b/src/multi-body-problem/multi-body-problem-2D.cfg
index 641078a70876d7286c6691333b63f81740c6898d..cfee414057863bdf3b609cd5bd751f2371cb48bc 100644
--- a/src/multi-body-problem/multi-body-problem-2D.cfg
+++ b/src/multi-body-problem/multi-body-problem-2D.cfg
@@ -1,6 +1,6 @@
# -*- mode:conf -*-
[boundary.friction]
-smallestDiameter = 0.05 # 2e-3 [m]
+smallestDiameter = 0.1 # 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-5 # 1e-5
diff --git a/src/multi-body-problem/multi-body-problem.cc b/src/multi-body-problem/multi-body-problem.cc
index 30fc02b349bf8299a6087c2910820503269d52e8..eb2fd28fff6e668ee0e982e432dd30f622ab426c 100644
--- a/src/multi-body-problem/multi-body-problem.cc
+++ b/src/multi-body-problem/multi-body-problem.cc
@@ -137,20 +137,21 @@ int main(int argc, char *argv[]) {
const size_t bodyCount = contactNetwork.nBodies();
for (size_t i=0; i<contactNetwork.nLevels(); i++) {
- // printDofLocation(contactNetwork.body(i)->gridView());
+ //printDofLocation(contactNetwork.body(i)->gridView());
//Vector def(contactNetwork.deformedGrids()[i]->size(dims));
//def = 1;
//deformedGridComplex.setDeformation(def, i);
- const auto& level = *contactNetwork.level(i);
+ /*const auto& level = *contactNetwork.level(i);
for (size_t j=0; j<level.nBodies(); j++) {
writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
- }
+ }*/
}
for (size_t i=0; i<bodyCount; i++) {
+ printDofLocation(contactNetwork.body(i)->gridView());
writeToVTK(contactNetwork.body(i)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
}
@@ -199,8 +200,7 @@ int main(int argc, char *argv[]) {
else
programState.setupInitialConditions(parset, contactNetwork);
-
- //DUNE_THROW(Dune::Exception, "Just need to stop here!");
+ //print(programState.u, "u:");
auto& nBodyAssembler = contactNetwork.nBodyAssembler();
for (size_t i=0; i<bodyCount; i++) {
@@ -296,6 +296,8 @@ int main(int argc, char *argv[]) {
};
report(true);
+ //DUNE_THROW(Dune::Exception, "Just need to stop here!");
+
// -------------------
// Set up TNNMG solver
// -------------------
@@ -303,7 +305,7 @@ int main(int argc, char *argv[]) {
BitVector totalDirichletNodes;
contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
- /*for (size_t i=0; i<totalDirichletNodes.size(); i++) {
+ 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];
@@ -313,9 +315,9 @@ int main(int argc, char *argv[]) {
for (size_t d=0; d<dims; d++) {
totalDirichletNodes[i][d] = val;
}
- }*/
+ }
- //print(totalDirichletNodes, "totalDirichletNodes:");
+ 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>;
@@ -332,18 +334,18 @@ int main(int argc, char *argv[]) {
BoundaryNodes dirichletNodes;
contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
- /*for (size_t i=0; i<dirichletNodes.size(); i++) {
+ 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++) {
+ for (size_t i=0; i<frictionNodes.size(); i++) {
print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
- }*/
+ }
Updaters current(
initRateUpdater(