Function CalcPotential

src/test/embed_potential/main.cpp:27–70

template <class M>
void CalcPotential(
    const MapEmbed<BCond<typename M::Scal>>& mebc, M& m, const Embed<M>& eb,
    FieldCell<typename M::Scal>& fcp, FieldEmbed<typename M::Scal>& feg,
    FieldCell<typename M::Scal>& fc_residual,
    std::shared_ptr<linear::Solver<M>> linsolver,
    FieldCell<typename M::Expr>& fc_system) {
  using ExprFace = typename M::ExprFace;
  using Expr = typename M::Expr;

  auto sem = m.GetSem();
  struct {
    FieldEmbed<ExprFace> feg;
  } * ctx(sem);
  auto& t = *ctx;

  if (sem("init")) {
    t.feg = UEmbed<M>::GradientImplicit(mebc, eb);
    fc_system.Reinit(m, Expr::GetUnit(0));
    for (auto c : eb.Cells()) {
      Expr sum(0);
      eb.LoopNci(c, [&](auto q) {
        const auto cf = eb.GetFace(c, q);
        eb.AppendExpr(
            sum, t.feg[cf] * eb.GetArea(cf) * eb.GetOutwardFactor(c, q), q);
      });
      fc_system[c] = sum;
    }
  }
  if (sem.Nested("solve")) {
    linsolver->Solve(fc_system, &fcp, fcp, m);
  }
  if (sem("flux")) {
    feg.Reinit(m, 0);
    eb.LoopFaces([&](auto cf) { //
      feg[cf] = UEmbed<M>::Eval(t.feg[cf], cf, fcp, eb);
    });
    fc_residual.Reinit(m, 0);
    for (auto c : m.Cells()) {
      fc_residual[c] = UEmbed<M>::Eval(fc_system[c], c, fcp, eb);
    }
  }
  if (sem()) {
  }
}

template <class M>
void Run(M& m, Vars& var) {