Method Interpolate

comp/gridfunction.cpp:285–374

  void GridFunction :: Interpolate (const CoefficientFunction & cf,
                                    const Region * reg, int mdcomp, LocalHeap & clh)
  {
    static Timer t("GridFunction::Interpolate"); RegionTimer r(t);
    shared_ptr<FESpace> fes = GetFESpace();
    shared_ptr<MeshAccess> ma = fes->GetMeshAccess(); 
    int dim   = fes->GetDimension();

    Array<int> cnti(fes->GetNDof());
    cnti = 0;

    auto tempvec = GetVector(mdcomp).CreateVector();
    
    if (reg)
      {
        auto regdofs = make_shared<BitArray> (fes->GetDofs(*reg));
        Projector proj(regdofs, false);
        tempvec = proj * GetVector(mdcomp);
      }
    else
      tempvec = 0.0;
   
    auto vb = reg ? reg->VB() : VOL;
    IterateElements 
      (*fes, vb, clh, 
       [&] (FESpace::Element el, LocalHeap & lh)
       {
         if (reg)
           if (!reg->Mask().Test(el.GetIndex())) return;
         
         if (!fes->DefinedOn(el)) return;

         const FiniteElement & fel = fes->GetFE (el, lh);
         int ndof = fel.GetNDof();
         int fesdim = fespace->GetDimension();
         // int dimcf = cf.Dimension();
         const ElementTransformation & eltrans = ma->GetTrafo (el, lh); 
         
         FlatVector<> elvec(ndof*fesdim, lh), elvec1(ndof*fesdim, lh);

         // GetElementVector (mdcomp, el.GetDofs(), elvec1);
         // cout << "elvec, before: " << endl << elvec1 << endl;
         
         // fel.Interpolate (eltrans, cf, elvec.AsMatrix(ndof/dimcf, dimcf), lh);
         fel.Interpolate (eltrans, cf, elvec.AsMatrix(ndof, fesdim), lh);
         // cout << "elvec: " << endl << elvec << endl;
         
         fes->TransformVec (el, elvec, TRANSFORM_SOL_INVERSE);

         // GetElementVector (mdcomp, el.GetDofs(), elvec1);
         tempvec.GetIndirect (el.GetDofs(), elvec1);
         elvec1 += elvec;
         // SetElementVector (mdcomp, el.GetDofs(), elvec1);
         tempvec.SetIndirect (el.GetDofs(), elvec1);
         
         for (auto d : el.GetDofs())
           if (IsRegularDof(d)) cnti[d]++;
       });