Method Update

comp/contact.cpp:260–335

  template<int DIM>
  void T_GapFunction<DIM> :: Update(shared_ptr<GridFunction> displacement_, int intorder2, double h_, bool both_sides)
  {
    h = h_;
    this->both_sides = both_sides;

    displacement = displacement_;

    LocalHeap lh(1000000, "T_GapFunction::Update");
    netgen::Box<DIM> bbox{netgen::Box<DIM>::EMPTY_BOX};
    double maxh = 0;
    for (Ngs_Element el2 : ma->Elements(BND))
      {
        HeapReset hr(lh);
        auto & mask = other.Mask();
        if (!mask.Test(el2.GetIndex())) continue;

        auto & trafo2 = ma->GetTrafo (el2, lh);
        auto & trafo2_def = trafo2.AddDeformation(displacement.get(), lh);

        IntegrationRule ir2_(trafo2.GetElementType(), intorder2);
        for(auto ir2 : ir2_.Split())
        {
            HeapReset hr(lh);
            MappedIntegrationRule<DIM-1, DIM> mir2(ir2, trafo2, lh);
            MappedIntegrationRule<DIM-1, DIM> mir2_def(ir2, trafo2_def, lh);

            netgen::Box<DIM> el_box{netgen::Box<DIM>::EMPTY_BOX};
            for (auto & mip : mir2_def)
              {
                netgen::Point<DIM> p;
                for (int j = 0; j < DIM; j++)
                  p(j) = mip.GetPoint()(j);
                bbox.Add(p);
                if(h==0.0)
                  el_box.Add(p);
              }
            maxh = max(maxh, el_box.Diam());
        }
      }

    // Default-value for h is 2 * maximum_element_diameter
    if(h==0.0)
      h = 2*maxh;

    bbox.Scale(2); // make sure we don't add boxes outside of tree bounding box
    searchtree = make_unique<netgen::BoxTree<DIM, int>>(bbox);
    for (Ngs_Element el2 : ma->Elements(BND))
      {
        HeapReset hr(lh);
        auto & mask = other.Mask();
        if (!mask.Test(el2.GetIndex())) continue;

        auto & trafo2 = ma->GetTrafo (el2, lh);
        auto & trafo2_def = trafo2.AddDeformation(displacement.get(), lh);

        IntegrationRule ir2_(trafo2.GetElementType(), intorder2);
        netgen::Box<DIM> elbox{netgen::Box<DIM>::EMPTY_BOX};
        for(auto ir2 : ir2_.Split())