Method BuildLocalExpansion

bem/ngbem.cpp:1266–1413

  template <typename KERNEL>
  void PotentialCF<KERNEL> ::
  BuildLocalExpansion(const Region & reg)
  {
    LocalHeapMem<100000> lh("PotentialCF::BuildLocalExpansion");

    auto space = this->gf->GetFESpace();
    auto mesh = space->GetMeshAccess();


    Vec<3> smax(-1e99, -1e99, -1e99);
    Vec<3> smin(1e99, 1e99, 1e99);

    for (size_t i = 0; i < mesh->GetNE(source_vb); i++)
      {
        HeapReset hr(lh);
        ElementId ei(source_vb, i);
        if (!space->DefinedOn(ei)) continue;
        if (definedon && !(*definedon).Mask().Test(mesh->GetElIndex(ei))) continue;

        // const FiniteElement &fel = space->GetFE(ei, lh);
        const ElementTransformation &trafo = mesh->GetTrafo(ei, lh);
        IntegrationRule ir(trafo.GetElementType(), intorder);
        auto & miry = trafo(ir, lh);

        for (int k = 0; k < miry.Size(); k++)
          for (int j = 0; j < 3; j++)
            {
              smin(j) = min(smin(j), miry[k].GetPoint()(j));
              smax(j) = max(smax(j), miry[k].GetPoint()(j));
            }
      }

    Vec<3> cs = 0.5*(smin+smax);
    double rs = MaxNorm(smax-smin);

    // cout << "cs = " << cs << ", rs = " << rs << endl;
    auto singmp = kernel.source.CreateMultipoleExpansion(cs, rs, io_params);

    typedef typename KERNEL::value_type T;
    for (size_t i = 0; i < mesh->GetNE(source_vb); i++)
      {
        HeapReset hr(lh);
        ElementId ei(source_vb, i);

        if (!space->DefinedOn(ei)) continue;
        if (definedon && !(*definedon).Mask().Test(mesh->GetElIndex(ei))) continue;

        const FiniteElement &fel = space->GetFE(ei, lh);
        const ElementTransformation &trafo = mesh->GetTrafo(ei, lh);

        Array<DofId> dnums(fel.GetNDof(), lh);
        space->GetDofNrs(ei, dnums);
        FlatVector<T> elvec(fel.GetNDof(), lh);
        gf->GetElementVector(dnums, elvec);

        IntegrationRule ir(fel.ElementType(), intorder);
        auto & miry = trafo(ir, lh);
        FlatMatrix<T> vals(miry.Size(), evaluator->Dim(), lh);