diff --git a/dune/tectonic/assemblers.cc b/dune/tectonic/assemblers.cc index 37c8cb39aa662cfd2d78d5c0b095d1b7ffa9aabe..0f5ac9906d032c31ad7b71bb8a5212e712614ec5 100644 --- a/dune/tectonic/assemblers.cc +++ b/dune/tectonic/assemblers.cc @@ -144,6 +144,7 @@ template <class GridView, int dimension> void MyAssembler<GridView, dimension>::assembleWeightedNormalStress( BoundaryPatch<GridView> const &frictionalBoundary, ScalarVector &weightedNormalStress, + ScalarVector &weights, double youngModulus, double poissonRatio, Vector const &displacement) const { @@ -159,7 +160,6 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress( auto const nodalTractionAverage = interpolateP0ToP1(frictionalBoundary, traction); - ScalarVector weights; { NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type> frictionalBoundaryAssembler( @@ -168,6 +168,7 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress( vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler, weights, frictionalBoundary); } + auto const normals = frictionalBoundary.getNormals(); for (size_t i = 0; i < vertexBasis.size(); ++i) weightedNormalStress[i] = diff --git a/dune/tectonic/assemblers.hh b/dune/tectonic/assemblers.hh index e5354141074f9646a6a0407b41d66c978018e707..d7ca0d8826450924558cde30afe8f3a478b24bb4 100644 --- a/dune/tectonic/assemblers.hh +++ b/dune/tectonic/assemblers.hh @@ -89,6 +89,7 @@ class MyAssembler { void assembleWeightedNormalStress( BoundaryPatch<GridView> const &frictionalBoundary, ScalarVector &weightedNormalStress, + ScalarVector &weights, double youngModulus, double poissonRatio, Vector const &displacement) const; diff --git a/dune/tectonic/data-structures/program_state.hh b/dune/tectonic/data-structures/program_state.hh index 4bc169ece1892c173eca33965861f7fdd3acf246..d4573f9bcdbc7de50c8e7ad5ee69267142322c86 100644 --- a/dune/tectonic/data-structures/program_state.hh +++ b/dune/tectonic/data-structures/program_state.hh @@ -74,6 +74,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState { v(bodyCount_), a(bodyCount_), alpha(bodyCount_), + weights(bodyCount_), weightedNormalStress(bodyCount_) { for (size_t i=0; i<bodyCount_; i++) { size_t leafVertexCount = leafVertexCounts[i]; @@ -82,6 +83,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState { v[i].resize(leafVertexCount), a[i].resize(leafVertexCount), alpha[i].resize(leafVertexCount), + weights[i].resize(leafVertexCount), weightedNormalStress[i].resize(leafVertexCount), bodyStates[i] = new BodyState(&u[i], &v[i], &a[i], &alpha[i], &weightedNormalStress[i]); @@ -205,12 +207,14 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState { const auto& body = contactNetwork.body(nonmortarIdx); ScalarVector frictionBoundaryStress(weightedNormalStress[nonmortarIdx].size()); + ScalarVector nmWeights(frictionBoundaryStress); body->assembler()->assembleWeightedNormalStress( - contactCoupling->nonmortarBoundary(), frictionBoundaryStress, body->data()->getYoungModulus(), + contactCoupling->nonmortarBoundary(), frictionBoundaryStress, nmWeights, body->data()->getYoungModulus(), body->data()->getPoissonRatio(), u[nonmortarIdx]); weightedNormalStress[nonmortarIdx] += frictionBoundaryStress; + weights[nonmortarIdx] += nmWeights; } std::cout << "solving linear problem for a..." << std::endl; @@ -231,6 +235,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState { std::vector<Vector> v; std::vector<Vector> a; std::vector<ScalarVector> alpha; + std::vector<ScalarVector> weights; std::vector<ScalarVector> weightedNormalStress; double relativeTime; double relativeTau;