| 10 | |
| 11 | template <class M> |
| 12 | void InitVfList( |
| 13 | FieldCell<typename M::Scal>& fc, std::istream& primlist, int approx, |
| 14 | size_t edim, const M& m, bool verbose) { |
| 15 | using Scal = typename M::Scal; |
| 16 | using Vect = typename M::Vect; |
| 17 | using Primitive = typename UPrimList<Vect>::Primitive; |
| 18 | const std::vector<Primitive> ppa = |
| 19 | UPrimList<Vect>::GetPrimitives(primlist, edim); |
| 20 | if (verbose && m.IsRoot()) { |
| 21 | std::cerr << "Read " << ppa.size() << " primitives" << std::endl; |
| 22 | } |
| 23 | |
| 24 | const Vect h = m.GetCellSize(); |
| 25 | // filter to bounding box |
| 26 | auto& bc = m.GetSuBlockCells(); |
| 27 | Rect<Vect> rect(Vect(bc.GetBegin()) * h, Vect(bc.GetEnd()) * h); |
| 28 | |
| 29 | std::vector<Primitive> pp; |
| 30 | for (auto& p : ppa) { |
| 31 | if (p.inter(rect)) { |
| 32 | pp.push_back(p); |
| 33 | } |
| 34 | } |
| 35 | |
| 36 | if (pp.empty()) { |
| 37 | fc.Reinit(m, 0.); |
| 38 | } else { |
| 39 | auto lsmax0 = [&pp](Vect x) -> std::pair<Scal, size_t> { |
| 40 | Scal lmax = -std::numeric_limits<Scal>::max(); // maximum level-set |
| 41 | size_t imax0 = 0; // index of maximum |
| 42 | for (size_t i = 0; i < pp.size(); ++i) { |
| 43 | auto& p = pp[i]; |
| 44 | Scal li = p.ls(x); |
| 45 | if (p.mod_minus) { |
| 46 | li = -li; |
| 47 | } |
| 48 | if (p.mod_and) { |
| 49 | lmax = std::min(lmax, li); |
| 50 | } else { |
| 51 | if (li > lmax) { |
| 52 | lmax = li; |
| 53 | imax0 = i; |
| 54 | } |
| 55 | } |
| 56 | } |
| 57 | return {lmax, imax0}; |
| 58 | }; |
| 59 | if (approx == 0) { // stepwise |
| 60 | for (auto c : m.Cells()) { |
| 61 | fc[c] = (lsmax0(m.GetCenter(c)).first >= 0 ? 1 : 0); |
| 62 | } |
| 63 | } else if (approx == 1) { // level set |
| 64 | for (auto c : m.Cells()) { |
| 65 | const auto x = m.GetCenter(c); |
| 66 | const auto& p = pp[lsmax0(x).second]; |
| 67 | fc[c] = GetLevelSetVolume<Scal>(p.ls, x, h); |
| 68 | } |
| 69 | } else if (approx == 2) { // overlap |
no test coverage detected