Method ConvertL2Operator

comp/fespace.cpp:2352–2484

  

  shared_ptr<BaseMatrix> FESpace :: ConvertL2Operator (shared_ptr<FESpace> l2space) const
  {
    LocalHeap lh(100000000);
    
    Array<int> classnr(ma->GetNE());
    
    FlatArray<int> vertex_map;
    if (const PeriodicFESpace * periodic = dynamic_cast<const PeriodicFESpace*> (this))
      vertex_map.Assign (periodic->GetVertexMap());
    
    ma->IterateElements
      (VOL, lh, [&] (auto el, LocalHeap & llh)
       {
         if (!DefinedOn (el) || !l2space->DefinedOn (el)) { classnr[el.Nr()] = -1; return; }
         classnr[el.Nr()] = 
           SwitchET<ET_TRIG,ET_QUAD,ET_TET,ET_HEX>
           (el.GetType(),
            [el, vertex_map] (auto et) {
             if (vertex_map.Size())
               return ET_trait<et.ElementType()>::GetClassNr(vertex_map[el.Vertices()]);
             else
               return ET_trait<et.ElementType()>::GetClassNr(el.Vertices());
            });
       });

    TableCreator<size_t> creator;
    for ( ; !creator.Done(); creator++)
      for (auto i : Range(classnr))
        if(classnr[i] >= 0)
          creator.Add (classnr[i], i);
    Table<size_t> table = creator.MoveTable();

    shared_ptr<BaseMatrix> sum;
  
    // size_t ne = ma->GetNE();

    // scaling for multiple dofs (if actually not an L2 space)
    Array<int> cnt(l2space->GetNDof());
    cnt = 0;
    
    for (auto elclass_inds : table)
      {
        HeapReset hr(lh);
        if (elclass_inds.Size() == 0) continue;
        
        ElementId ei(VOL,elclass_inds[0]);
        auto & fel = GetFE (ei, lh);
        auto & fel_l2 = l2space->GetFE (ei, lh);

        auto & trafo = GetFEElementTransformation(fel.ElementType());
        MixedFiniteElement fel_mixed(fel, fel_l2);
        auto evaluator = GetEvaluator(VOL);
        auto l2evaluator = l2space->GetEvaluator(VOL);

        auto trial = make_shared<ProxyFunction>(dynamic_pointer_cast<FESpace>(const_cast<FESpace*>(this)->shared_from_this()),
                                                false, false, evaluator,
                                                nullptr, nullptr, nullptr, nullptr, nullptr);
        auto trial_l2 = make_shared<ProxyFunction>(l2space,