Function Run

examples/104_embed_diffusion/main.cpp:103–185

}

void Run(M& m, Vars& var) {
  auto sem = m.GetSem();
  struct {
    FieldNode<Scal> fnl;
    FieldCell<Scal> fcu;
    MapEmbed<BCond<Scal>> mebc;
    Scal t = 0;
    Scal diff;
    Scal dt;
    std::unique_ptr<Embed<M>> eb_;
  } * ctx(sem);
  auto& fnl = ctx->fnl;
  auto& fcu = ctx->fcu;
  auto& mebc = ctx->mebc;
  auto& t = ctx->t;
  auto& dt = ctx->dt;
  auto& diff = ctx->diff;

  if (sem("init")) {
    diff = var.Double["diff"];
    dt = var.Double["cfl"] * sqr(m.GetCellSize()[0]) * diff;
    // initial field
    fcu.Reinit(m, 0);
    for (auto c : m.AllCells()) {
      auto x = m.GetCenter(c);
      fcu[c] = ((Vect(0.5, 0.5, 0) - x).norminf() < 0.1);
    }
    // level-set for embedded boundaries
    fnl.Reinit(m, 0);
    for (auto n : m.AllNodes()) {
      auto x = m.GetNode(n);
      fnl[n] = (0.4 - (Vect(0.5, 0.5, 0) - Vect(x[0], x[1], 0.)).norm());
    }
    ctx->eb_.reset(new Embed<M>(m, var.Double["gradlim"]));
    m.Dump(&fcu, "u");
  }
  if (sem.Nested("eb-init")) {
    ctx->eb_->Init(fnl);
  }
  if (sem.Nested("eb-dump")) {
    ctx->eb_->DumpPlaneSection(Vect(0., 0., 1e-3), Vect(0., 0., 1.));
  }
  if (sem("init2")) {
    auto& eb = *ctx->eb_;
    // boundary conditions
    for (auto c : eb.SuCFaces()) {
      auto& bc = mebc[c];
      bc.type = BCondType::dirichlet;
      bc.val = 1;
    }
  }
  sem.LoopBegin();
  if (sem("step")) {
    if (m.IsRoot()) {
      std::cout << "t=" << t << std::endl;
    }
    switch (var.Int["case"]) {
      case 0: