| 93 | |
| 94 | template <class M> |
| 95 | void InitOverlappingComponents( |
| 96 | std::istream& primlist, const Multi<FieldCell<typename M::Scal>*>& fcu, |
| 97 | const Multi<FieldCell<typename M::Scal>*>& fccl, |
| 98 | const GRange<size_t>& layers, const M& m) { |
| 99 | using Scal = typename M::Scal; |
| 100 | using Vect = typename M::Vect; |
| 101 | using Primitive = typename UPrimList<Vect>::Primitive; |
| 102 | constexpr Scal kClNone = -1; // no color |
| 103 | const std::vector<Primitive> primitives = |
| 104 | UPrimList<Vect>::GetPrimitives(primlist, m.flags.edim); |
| 105 | fassert_equal(fcu.size(), layers.size()); |
| 106 | fassert_equal(fccl.size(), layers.size()); |
| 107 | for (auto l : layers) { |
| 108 | fcu[l]->Reinit(m, 0); |
| 109 | fccl[l]->Reinit(m, kClNone); |
| 110 | } |
| 111 | const auto h = m.GetCellSize(); |
| 112 | auto levelset = [&](size_t i, Vect x) { |
| 113 | auto& p = primitives[i]; |
| 114 | return p.ls(x) * (p.mod_minus ? -1 : 1); |
| 115 | }; |
| 116 | for (auto c : m.AllCellsM()) { |
| 117 | const size_t none = -1; |
| 118 | size_t imax0 = none; // primitive with largest levelset in c |
| 119 | size_t imax1 = none; // primitive with second largest levelset in c |
| 120 | Scal lsmax0 = -std::numeric_limits<Scal>::max(); |
| 121 | Scal lsmax1 = -std::numeric_limits<Scal>::max(); |
| 122 | Vect image0(0); // image vector for periodic conditions |
| 123 | Vect image1(0); // image vector for periodic conditions |
| 124 | |
| 125 | const auto wimages = [&]() { // |
| 126 | using MIdx = generic::MIdx<M::dim>; |
| 127 | MIdx begin(0); |
| 128 | MIdx size(1); |
| 129 | for (auto d : m.dirs) { |
| 130 | if (m.flags.is_periodic[d]) { |
| 131 | begin[d] = -1; |
| 132 | size[d] = 3; |
| 133 | } |
| 134 | } |
| 135 | return GBlock<int, M::dim>{begin, size}; |
| 136 | }(); |
| 137 | |
| 138 | for (size_t i = 0; i < primitives.size(); ++i) { |
| 139 | for (auto wimage : wimages) { |
| 140 | const Vect image = Vect(wimage) * m.GetGlobalLength(); |
| 141 | const Scal ls = levelset(i, c.center() + image); |
| 142 | if (ls > lsmax0) { |
| 143 | imax0 = i; |
| 144 | lsmax0 = ls; |
| 145 | image0 = image; |
| 146 | } |
| 147 | } |
| 148 | } |
| 149 | for (size_t i = 0; i < primitives.size(); ++i) { |
| 150 | for (auto wimage : wimages) { |
| 151 | const Vect image = Vect(wimage) * m.GetGlobalLength(); |
| 152 | const Scal ls = levelset(i, c.center() + image); |
no test coverage detected