Method CreateMatrixFMM

bem/ngbem.cpp:106–472

  template <typename KERNEL>
  shared_ptr<BaseMatrix> GenericIntegralOperator<KERNEL> ::
  CreateMatrixFMM(LocalHeap & lh) const
  {
    static Timer tall("ngbem fmm setup"); RegionTimer r(tall);
    Array<Vec<3>> xpts, ypts, xnv, ynv;
    // IntegrationRule ir(ET_TRIG, intorder);
    IntegrationRule ir_trig(ET_TRIG, intorder);
    IntegrationRule ir_quad(ET_QUAD, intorder);
    auto trial_mesh = trial_space->GetMeshAccess();
    auto test_mesh = test_space->GetMeshAccess();

    Array<int> compress_trial_els(trial_mesh->GetNE(BND));
    Array<int> compress_test_els(test_mesh->GetNE(BND));
    compress_trial_els = -1;
    compress_test_els = -1;

    int cnt = 0;
    for (auto el : trial_mesh->Elements(BND))
      if (trial_space->DefinedOn(el))
        if (!trial_definedon || (*trial_definedon).Mask().Test(trial_mesh->GetElIndex(el)))
          {
            HeapReset hr(lh);
            auto & trafo = trial_mesh->GetTrafo(el, lh);
            // auto & mir = static_cast<MappedIntegrationRule<2,3>&>(trafo(ir, lh));
            auto & ir = (trafo.GetElementType() == ET_QUAD) ? ir_quad : ir_trig;
            auto & mir = static_cast<MappedIntegrationRule<2,3>&>(trafo(ir, lh));
            for (auto & mip : mir)
              {
                xpts.Append(mip.GetPoint());
                xnv.Append(mip.GetNV());
              }
            compress_trial_els[el.Nr()] = cnt++;
          }

    cnt = 0;
    for (auto el : test_mesh->Elements(BND))
      if (test_space->DefinedOn(el))
        if (!test_definedon || (*test_definedon).Mask().Test(test_mesh->GetElIndex(el)))
          {
            HeapReset hr(lh);
            auto & trafo = test_mesh->GetTrafo(el, lh);
            // auto & mir = static_cast<MappedIntegrationRule<2,3>&>(trafo(ir, lh));
            auto & ir = (trafo.GetElementType() == ET_QUAD) ? ir_quad : ir_trig;
            auto & mir = static_cast<MappedIntegrationRule<2,3>&>(trafo(ir, lh));
            for (auto & mip : mir)
              {
                ypts.Append(mip.GetPoint());
                ynv.Append(mip.GetNV());
              }
            compress_test_els[el.Nr()] = cnt++;
          }

    auto create_eval = [&](const FESpace & fes,
                           const Array<int> & compress_els,
                           const DifferentialOperator & evaluator)
    {
      auto mesh = fes.GetMeshAccess();
      Array<short> classnr(mesh->GetNE(BND));