Diffusion solver with redistribution to neighbor cells and LoopFaces() to avoid code duplicatoin.
| 75 | // Diffusion solver with redistribution to neighbor cells |
| 76 | // and LoopFaces() to avoid code duplicatoin. |
| 77 | void Diffusion2( |
| 78 | FieldCell<Scal>& fcu, const MapEmbed<BCond<Scal>>& mebc, Scal diff, Scal dt, |
| 79 | const Embed<M>& eb) { |
| 80 | const auto feg = UEmbed<M>::Gradient(fcu, mebc, eb); |
| 81 | // Compute flux. |
| 82 | FieldEmbed<Scal> fed(eb, 0); |
| 83 | eb.LoopFaces([&](auto cf) { // lambda-function applied to faces and |
| 84 | // embedded faces |
| 85 | fed[cf] = feg[cf] * diff * eb.GetArea(cf); |
| 86 | }); |
| 87 | // Compute the change at one time step. |
| 88 | FieldCell<Scal> fct(eb, 0); |
| 89 | for (auto c : eb.Cells()) { |
| 90 | Scal sum = fed[c]; |
| 91 | for (auto q : eb.Nci(c)) { |
| 92 | sum += fed[eb.GetFace(c, q)] * eb.GetOutwardFactor(c, q); |
| 93 | } |
| 94 | fct[c] = sum * dt; |
| 95 | } |
| 96 | fct = UEmbed<M>::RedistributeCutCells(fct, eb); |
| 97 | // Advance in time. |
| 98 | for (auto c : eb.Cells()) { |
| 99 | fcu[c] += fct[c] / eb.GetVolume(c); |
| 100 | } |
| 101 | } |
| 102 | |
| 103 | void Run(M& m, Vars& var) { |
| 104 | auto sem = m.GetSem(); |