Function CalcFluxProject

comp/postproc.cpp:64–211

  
  template <class SCAL>
  void CalcFluxProject (const S_GridFunction<SCAL> & u,
			S_GridFunction<SCAL> & flux,
			shared_ptr<BilinearFormIntegrator> bli,
			bool applyd, const BitArray & domains, LocalHeap & clh)
  {
    static Timer timer("CalcFluxProject");
    RegionTimer reg (timer);
    
    auto fes = u.GetFESpace();
    auto fesflux = flux.GetFESpace();

    shared_ptr<MeshAccess> ma = fesflux->GetMeshAccess();

    BaseStatusHandler::PushStatus ("Post-processing");
    

    auto vb = bli->VB();

    int ne      = ma->GetNE(vb);
    int dim     = fes->GetDimension();
    int dimflux = fesflux->GetDimension();
    int dimfluxvec = bli->DimFlux(); 

    shared_ptr<BilinearFormIntegrator> fluxbli = fesflux->GetIntegrator(vb);
    shared_ptr<BilinearFormIntegrator> single_fluxbli = fluxbli;
    if (dynamic_pointer_cast<BlockBilinearFormIntegrator> (single_fluxbli))
      single_fluxbli = dynamic_pointer_cast<BlockBilinearFormIntegrator> (single_fluxbli)->BlockPtr();
    
    auto flux_evaluator = fesflux->GetEvaluator(vb);
    if (!fluxbli)
      {
        cout << IM(5) << "make a symbolic integrator for CalcFluxProject" << endl;
        auto single_evaluator =  flux_evaluator;
        if (dynamic_pointer_cast<BlockDifferentialOperator>(single_evaluator))
          single_evaluator = dynamic_pointer_cast<BlockDifferentialOperator>(single_evaluator)->BaseDiffOp();
        
        auto trial = make_shared<ProxyFunction>(fesflux, false, false, single_evaluator,
                                                nullptr, nullptr, nullptr, nullptr, nullptr);
        auto test  = make_shared<ProxyFunction>(fesflux, true, false, single_evaluator,
                                                nullptr, nullptr, nullptr, nullptr, nullptr);
        fluxbli = make_shared<SymbolicBilinearFormIntegrator> (InnerProduct(trial,test), vb, VOL);
        single_fluxbli = fluxbli;
        // throw Exception ("no integrator available");
      }

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

    flux.GetVector() = 0.0;

    ProgressOutput progress (ma, "postprocessing element", ne);

    IterateElements  
      (*fesflux, vb, clh, 
       [&] (Ngs_Element ei, LocalHeap & lh)
       {
         HeapReset hr(lh);
         progress.Update ();