Method GetTraceOperator

comp/l2hofespace.cpp:1960–2042

  shared_ptr<BaseMatrix> L2HighOrderFESpace ::
  GetTraceOperator (shared_ptr<FESpace> tracespace, bool avg) const
  {
    LocalHeap lh(1000000);
    Array<short> classnr(ma->GetNE());
    ma->IterateElements
      (VOL, lh, [&] (auto el, LocalHeap & llh)
       {
         classnr[el.Nr()] =
           SwitchET<ET_TRIG,ET_TET>
           (el.GetType(),
            [el] (auto et) { return ET_trait<et.ElementType()>::GetClassNr(el.Vertices()); });
       });

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

    shared_ptr<BaseMatrix> sum;

    // size_t ne = ma->GetNE();

    shared_ptr<VVector<double>> cnt;
    if (avg)
      {
        cnt = make_shared<VVector<double>>(tracespace->GetNDof());
        *cnt = 0;
      }


    for (auto elclass_inds : table)
      {
        if (elclass_inds.Size() == 0) continue;

        ElementId ei(VOL,elclass_inds[0]);
        auto & felx = GetFE (ei, lh);
        //auto & trafo = GetMeshAccess()->GetTrafo(ei, lh);

        Matrix<> trace_op_x = GetTraceMatrix (felx);


        Table<DofId> xdofs(elclass_inds.Size(), felx.GetNDof()),
          ydofs(elclass_inds.Size(), trace_op_x.Height());

        Array<DofId> dnumsx, dnumsy;
        for (auto i : Range(elclass_inds))
          {
            ElementId ei(VOL, elclass_inds[i]);
            GetDofNrs(ei, dnumsx);
            tracespace->GetDofNrs(ei, dnumsy);
            xdofs[i] = dnumsx;
            ydofs[i] = dnumsy;

            if (avg)
              for (auto d : dnumsy)
                (*cnt)(d) += 1;