Method AddSourceElementContribution

bem/ngbem.cpp:1558–1617

  template <typename KERNEL> template <typename T>
  void PotentialCF<KERNEL> ::
  AddSourceElementContribution(const BaseMappedIntegrationPoint & mip,
                               ElementId ei,
                               const IntegrationRule & ir,
                               FlatVector<T> result,
                               T scale,
                               LocalHeap & lh) const
  {
    auto space = this->gf->GetFESpace();
    auto mesh = space->GetMeshAccess();

    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);

    SIMD_IntegrationRule simd_ir(ir);
    Vector<SIMD<T>> simd_result(Dimension());
    simd_result = SIMD<T>(0.0);

    static constexpr int bs = 64;
    for (int k = 0; k < simd_ir.Size(); k += bs)
      {
        HeapReset hr(lh);
        auto simd_ir_range = simd_ir.Range(k, min(simd_ir.Size(), size_t(k+bs)));
        auto & miry = trafo(simd_ir_range, lh);
        FlatMatrix<SIMD<T>> vals(evaluator->Dim(), miry.Size(), lh);

        evaluator->Apply(fel, miry, elvec, vals);
        for (int iy = 0; iy < miry.Size(); iy++)
          {
            Vec<3,SIMD<double>> x = mip.GetPoint();
            Vec<3,SIMD<double>> nx{0.0};
            if constexpr (KERNEL::target_type::needs_normal)
              nx = dynamic_cast<const MappedIntegrationPoint<2,3>&>(mip).GetNV();

            Vec<3,SIMD<double>> y = miry[iy].GetPoint();
            Vec<3,SIMD<double>> ny{0.0};
            if constexpr (KERNEL::source_type::needs_normal)
              {
                if (source_vb != BND)
                  throw Exception("kernel requires boundary source normals");
                ny = static_cast<const SIMD<MappedIntegrationPoint<2,3>>&>(miry[iy]).GetNV();
              }

            auto eval = kernel.Evaluate(x, y, nx, ny);
            for (auto term : kernel.terms)
              {
                auto kernel_ = term.fac * eval(term.kernel_comp);
                simd_result(term.test_comp) += miry[iy].GetWeight() * kernel_ * vals(term.trial_comp, iy);
              }
          }
      }

    for (int i = 0; i < Dimension(); i++)