Method T_Integrate

comp/postproc.cpp:1531–1828

  
  template <typename TSCAL>
  TSCAL Integral :: T_Integrate (const ngcomp::MeshAccess & ma,
                               FlatVector<TSCAL> element_wise)
  {
    LocalHeap glh(10000000, "integrate-lh");
    bool use_simd = true;
    TSCAL sum = 0.0;

    BitArray defon;

    if (dx.definedon)
      {
        if (auto definedon_bitarray = get_if<BitArray> (&*dx.definedon))
          defon = *definedon_bitarray;
        if (auto definedon_string = get_if<string> (&*dx.definedon))
          {
            shared_ptr<MeshAccess> spma(const_cast<MeshAccess*>(&ma), NOOP_Deleter);
            Region reg(spma, dx.vb, *definedon_string);
            defon = reg.Mask();
          }
      }

    constexpr int max_element_type = std::max(element_types);
    std::array<IntegrationRule*, max_element_type+1> user_ir;
    for (auto & ir_ptr : user_ir)
      ir_ptr = nullptr;
    std::array<unique_ptr<SIMD_IntegrationRule>, max_element_type+1> user_simd_ir;

    for (auto [type, val] : dx.userdefined_intrules)
      {
        user_ir[type] = val.get();
        user_simd_ir[type] = make_unique<SIMD_IntegrationRule>(*user_ir[type]);
      }
    
    if (dx.element_vb == VOL)
      {
        ma.IterateElements
          (this->dx.vb, glh, [&] (Ngs_Element el, LocalHeap & lh)
           {
             if (this->dx.definedonelements && !this->dx.definedonelements->Test(el.Nr())) return;
             // if(!mask.Test(el.GetIndex())) return;
             auto & trafo1 = ma.GetTrafo (el, lh);
             auto & trafo = trafo1.AddDeformation(this->dx.deformation.get(), lh);

             if (defon.Size() && !defon.Test(el.GetIndex()))
               return;
             
             TSCAL hsum = 0.0;
             
             bool this_simd = use_simd;
             int order = 5+this->dx.bonus_intorder;
             
             if (this_simd)
               {
                 try
                   {
                     SIMD_IntegrationRule ir(trafo.GetElementType(), order);
                     auto myir = user_simd_ir[trafo.GetElementType()] ? user_simd_ir[trafo.GetElementType()].get() : &ir;
                     auto & mir = trafo(*myir, lh);