compute relative velocities

src/spatial-solving/fixedpointiterator.cc

@@ -119,8 +119,14 @@ std::ostream &operator<<(std::ostream &stream,
 
 template <class Factory, class Updaters, class ErrorNorm>
 void FixedPointIterator<Factory, Updaters, ErrorNorm>::relativeVelocities(std::vector<Vector>& v_m) const {
+  // needs assemblers to obtain basis
+  std::vector<std::shared_ptr<MyAssembler>> assemblers(bodyCount);
+
   using field_type = typename Factory::Matrix::field_type;
 
+
+
+
   // adaptation of DualMortarCoupling::setup()
 
   const size_t dim = DeformedGrid::dimension;
@@ -134,11 +140,11 @@ void FixedPointIterator<Factory, Updaters, ErrorNorm>::relativeVelocities(std::v
   using DualCache = Dune::Contact::DualBasisAdapter<GridView, field_type>;
   std::unique_ptr<DualCache> dualCache;
   dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView, field_type> >();
-/*
+
   // define FE grid functions
-  std::vector<BasisGridFunction<> > gridFunctions(nBodyAssembler_.nGrids());
+  std::vector<BasisGridFunction<VertexBasis, Vector>* > gridFunctions(v_m.size());
   for (size_t i=0; i<gridFunctions.size(); i++) {
-
+    gridFunctions[i] = new BasisGridFunction<MyAssembler::VertexBasis, Vector>(assemblers[i]->vertexBasis, v_m[i]);
   }
 
   const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
@@ -245,63 +251,6 @@ void FixedPointIterator<Factory, Updaters, ErrorNorm>::relativeVelocities(std::v
 
   }
 
-      ///////////////////////////////////
-      //  reducing nonmortar boundary
-      /////////////////////////////////
-
-      // Get all fine grid boundary segments that are totally covered by the grid-glue segments
-      typedef std::pair<int,int> Pair;
-      std::map<Pair,ctype> coveredArea, fullArea;
-
-      // initialize with area of boundary faces
-      for (const auto& bIt : nonmortarBoundary_) {
-          const Pair p(indexSet0.index(bIt.inside()),bIt.indexInInside());
-          fullArea[p] = bIt.geometry().volume();
-          coveredArea[p] = 0;
-      }
-
-      // sum up the remote intersection areas to find out which are totally covered
-      for (const auto& rIs : intersections(glue))
-          coveredArea[Pair(indexSet0.index(rIs.inside()),rIs.indexInInside())] += rIs.geometry().volume();
-
-      // add all fine grid faces that are totally covered by the contact mapping
-      for (const auto& bIt : nonmortarBoundary_) {
-          const auto& inside = bIt.inside();
-          if(coveredArea[Pair(indexSet0.index(inside),bIt.indexInInside())]/
-              fullArea[Pair(indexSet0.index(inside),bIt.indexInInside())] >= coveredArea_)
-              boundaryWithMapping.addFace(inside, bIt.indexInInside());
-      }
-
-      //writeBoundary(boundaryWithMapping,debugPath_ + "relevantNonmortar");
-
-
-      // \todo replace by all fine grid segments which are totally covered by the RemoteIntersections.
-      //for (const auto& rIs : intersections(glue))
-      //    boundaryWithMapping.addFace(rIs.inside(),rIs.indexInInside());
-
-      printf("contact mapping could be built for %d of %d boundary segments.\n",
-             boundaryWithMapping.numFaces(), nonmortarBoundary_.numFaces());
-
-      nonmortarBoundary_ = boundaryWithMapping;
-      mortarBoundary_.setup(gridView1);
-      for (const auto& rIs : intersections(glue))
-          if (nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
-              mortarBoundary_.addFace(rIs.outside(),rIs.indexInOutside());
-
-
-      // Assemble the diagonal matrix coupling the nonmortar side with the lagrange multiplyers there
-      assembleNonmortarLagrangeMatrix();
-
-      // The weak obstacle vector
-      weakObstacle_.resize(nonmortarBoundary_.numVertices());
-      weakObstacle_ = 0;
-
-      // ///////////////////////////////////////////////////////////
-      //   Get the occupation structure for the mortar matrix
-      // ///////////////////////////////////////////////////////////
-
-      // 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;
@@ -340,46 +289,7 @@ void FixedPointIterator<Factory, Updaters, ErrorNorm>::relativeVelocities(std::v
           }
       }
 
-      mortarIndices.exportIdx(mortarLagrangeMatrix_);
-
-      // 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> >();
-
-      std::vector<Dune::FieldVector<ctype,dim> > avNormals;
-      avNormals = nonmortarBoundary_.getNormals();
-
-
-      }
-
-      // ///////////////////////////////////////
-      //    Compute M = D^{-1} \hat{M}
-      // ///////////////////////////////////////
-
-      Dune::BCRSMatrix<MatrixBlock>& M  = mortarLagrangeMatrix_;
-      Dune::BDMatrix<MatrixBlock>& D    = nonmortarLagrangeMatrix_;
-
-      // First compute D^{-1}
-      D.invert();
-
-      // Then the matrix product D^{-1} \hat{M}
-      for (auto rowIt = M.begin(); rowIt != M.end(); ++rowIt) {
-          const auto rowIndex = rowIt.index();
-          for (auto& entry : *rowIt)
-              entry.leftmultiply(D[rowIndex][rowIndex]);
-      }
-
-      // weakObstacles in transformed basis = D^{-1}*weakObstacle_
-      for(size_t rowIdx=0; rowIdx<weakObstacle_.size(); rowIdx++)
-          weakObstacle_[rowIdx] *= D[rowIdx][rowIdx][0][0];
 
-      gridGlueBackend_->clear(); */
 }
 
 #include "fixedpointiterator_tmpl.cc"