Method T_Evaluate_impl

comp/interpolate.cpp:183–295

    template <typename MIR, typename T, ORDERING ORD>
    void T_Evaluate_impl (const MIR & ir, BareSliceMatrix<T,ORD> values) const
    {
      // #ifdef FIRSTDRAFT
      LocalHeapMem<2000000> lh("interpolate");

      // static Timer t("interpolate");
      // RegionTracer reg(TaskManager::GetThreadId(), t);    

      const ElementTransformation & trafo = ir.GetTransformation();
      // const MeshAccess & ma = *static_cast<const MeshAccess*> (trafo.GetMesh());
      ElementId ei = trafo.GetElementId();
      auto & fel = fes->GetFE(ei, lh);
      // int dim   = fes->GetDimension();
      int dim = func->Dimension();

      
      // cout << " eval for ei " << ei << endl;
      // cout << " ndof = " << fel.GetNDof() << endl;

      // if (dim != 1)
      // { throw Exception("Dim != 1 porbably does not work (yet)"); }

      /** func * dual_shape **/
      FlatVector<T> elflux(fel.GetNDof(), lh);

      FlatVector<T> elfluxadd(fel.GetNDof(), lh);  elflux = 0; // non-SIMD version
      for (auto el_vb : fes->GetDualShapeNodes(trafo.VB()))
	{
          if (el_vb == VOL)
            {
              IntegrationRule ir(fel.ElementType(), 2*fel.Order()+bonus_intorder);
              auto & mir = trafo(ir, lh);
              FlatMatrix<T> mflux(ir.Size(), dim, lh);
              func->Evaluate (mir, mflux);
              for (size_t j : Range(mir))
                mflux.Row(j) *= mir[j].GetWeight();
              dual_diffop -> ApplyTrans (fel, mir, mflux, elfluxadd, lh);
              elflux += elfluxadd;
            }
          else
            {
              Facet2ElementTrafo f2el (fel.ElementType(), el_vb);
              for (int locfnr : Range(f2el.GetNFacets()))
                {
                  // SIMD does not work yet
                  // SIMD_IntegrationRule irfacet(f2el.FacetType(locfnr), 2 * fel.Order());
                  IntegrationRule irfacet(f2el.FacetType(locfnr), 2*fel.Order()+bonus_intorder);
                  auto & irvol = f2el(locfnr, irfacet, lh);
                  auto & mir = trafo(irvol, lh);
                  mir.ComputeNormalsAndMeasure(fel.ElementType(), locfnr);
                  
                  // FlatMatrix<T,ORD> mflux(dim, irfacet.Size(), lh);
                  // func->Evaluate (mir, mflux);
                  // for (size_t j : Range(mir))
                  // 	mflux.Col(j) *= mir[j].GetWeight();
                  // SIMD only
                  // dual_diffop->AddTrans (fel, mir, mflux, elflux);
                  // NON-simd version