Line data Source code
1 : //
2 : // Fill a domain based on a list of points to form one or more closed 2D shapes.
3 : // The distance from each point to the nearest edge is found and an error function is then used to compute
4 : // the value of phi between 0 (outside) and 1 (inside) the domain. Only works in 2 dimensions.
5 : // Points belonging to different polygons are distinguished using the ObjNum column of the input file:
6 : // whenever ObjNum changes by more than obj_num_threshold between consecutive lines, a new polygon begins.
7 : // Each polygon's field is computed into phi_temp. Polygons are then combined into either the "solid" field
8 : // (pointwise max over all non-void polygons) or the "void" field (pointwise max over all void polygons,
9 : // as selected by the per-polygon invert list). The void field is inverted (zero inside, one outside the
10 : // void polygons) and multiplied into the solid field, punching holes wherever a void polygon overlaps it.
11 : // invert_all is then applied once at the end to flip the whole finished field if desired.
12 :
13 : #ifndef IC_POINTLIST_H_
14 : #define IC_POINTLIST_H_
15 :
16 : #include "Set/Set.H"
17 : #include "IC/IC.H"
18 : #include "Util/Util.H"
19 : using namespace std;
20 : #include <iostream>
21 : #include <vector>
22 : #include <algorithm>
23 : #include <iterator>
24 : #include <mpi.h>
25 : #include "IO/ParmParse.H"
26 :
27 : namespace IC
28 : {
29 : class PointList : public IC<Set::Scalar>
30 : {
31 : public:
32 : static constexpr const char* name = "pointlist";
33 :
34 : enum Type
35 : {
36 : Partition,
37 : Values
38 : };
39 :
40 : PointList(amrex::Vector<amrex::Geometry>& _geom) : IC(_geom) {}
41 :
42 2 : PointList(amrex::Vector<amrex::Geometry>& _geom, IO::ParmParse& pp, std::string name) : IC(_geom)
43 : {
44 2 : pp.queryclass(name, *this);
45 2 : }
46 : void Define() {
47 :
48 : };
49 :
50 22 : void Add(const int& lev, Set::Field<Set::Scalar>& a_phi, Set::Scalar)
51 : {
52 : if (AMREX_SPACEDIM != 2)
53 : {
54 : amrex::Abort("This code only supports 2D (AMREX_SPACEDIM must be 2)");
55 : }
56 :
57 22 : if ((int)phi_temp.size() != a_phi.size()) phi_temp.resize(a_phi.size());
58 22 : if ((int)phi_void.size() != a_phi.size()) phi_void.resize(a_phi.size());
59 36 : if (!phi_temp[lev] ||
60 36 : phi_temp[lev]->boxArray() != a_phi[lev]->boxArray() ||
61 11 : phi_temp[lev]->DistributionMap() != a_phi[lev]->DistributionMap())
62 : {
63 22 : phi_temp[lev] = std::make_unique<amrex::MultiFab>(a_phi[lev]->boxArray(), a_phi[lev]->DistributionMap(),
64 33 : a_phi[lev]->nComp(), a_phi[lev]->nGrow());
65 : }
66 36 : if (!phi_void[lev] ||
67 36 : phi_void[lev]->boxArray() != a_phi[lev]->boxArray() ||
68 11 : phi_void[lev]->DistributionMap() != a_phi[lev]->DistributionMap())
69 : {
70 22 : phi_void[lev] = std::make_unique<amrex::MultiFab>(a_phi[lev]->boxArray(), a_phi[lev]->DistributionMap(),
71 33 : a_phi[lev]->nComp(), a_phi[lev]->nGrow());
72 : }
73 :
74 22 : bool first_solid = true;
75 22 : bool first_void = true;
76 22 : bool any_void = false;
77 :
78 68 : for (unsigned int p = 0; p < polygons.size(); p++)
79 : {
80 46 : std::vector<Set::Vector> polygon = polygons[p];
81 46 : if (polygon.back() != polygon.front())
82 : // The first and last point must match, if they do not add another point to the end
83 : {
84 46 : polygon.push_back(polygon.front());
85 : }
86 46 : bool is_void = (p < invert.size() && invert[p] != 0);
87 46 : bool is_first = is_void ? first_void : first_solid;
88 :
89 342 : for (amrex::MFIter mfi(*a_phi[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
90 : {
91 296 : amrex::Box bx;
92 296 : amrex::IndexType type = a_phi[lev]->ixType();
93 592 : if (type == amrex::IndexType::TheCellType()) bx = mfi.growntilebox();
94 0 : else if (type == amrex::IndexType::TheNodeType()) bx = mfi.grownnodaltilebox();
95 0 : else Util::Abort(INFO, "Unkonwn index type");
96 296 : amrex::Array4<Set::Scalar> const& phi = a_phi[lev]->array(mfi);
97 296 : amrex::Array4<Set::Scalar> const& phi_t = phi_temp[lev]->array(mfi);
98 296 : amrex::Array4<Set::Scalar> const& phi_v = phi_void[lev]->array(mfi);
99 :
100 296 : std::vector<Set::Vector> X = polygon;
101 :
102 296 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
103 : {
104 203288 : Set::Vector x = Set::Position(i, j, k, geom[lev], type);
105 : Set::Scalar x_cross;
106 203288 : Set::Scalar min_dist_to_edge = std::numeric_limits<Set::Scalar>::max(); // Mininmum distance from point to nearest edge
107 : Set::Scalar dist_to_edge; // Distance from point to nearest edge
108 : Set::Scalar x1;
109 : Set::Scalar y1;
110 : Set::Scalar x2;
111 : Set::Scalar y2;
112 203288 : int num_cross = 0;
113 :
114 : Set::Scalar px;
115 : Set::Scalar py;
116 : Set::Scalar dx;
117 : Set::Scalar dy;
118 : Set::Scalar len2;
119 : Set::Scalar t;
120 203288 : Set::Vector x_int; // x and y position of the intercept of the mesh point and the edge, assuming edge is infinite
121 :
122 2594496 : for (unsigned int n = 0; n < X.size()-1; n++)
123 : {
124 :
125 2391208 : x1 = X[n](0);
126 2391208 : y1 = X[n](1);
127 2391208 : x2 = X[n+1](0);
128 2391208 : y2 = X[n+1](1);
129 :
130 2391208 : if ((y1 > x(1)) != (y2 > x(1)))
131 : {
132 154664 : x_cross = x1 + (x(1) - y1) * (x2 - x1) / (y2 - y1);
133 :
134 154664 : if (x_cross > x(0))
135 76272 : num_cross++;
136 : }
137 :
138 : // Computes the closest point on a line segment (X[n] → X[n+1]) to a query point x,
139 : // using vector projection onto the segment.
140 : //
141 : // Let A = X[n], B = X[n+1], and P = x.
142 : //
143 : // The parameter t is the normalized projection of P onto the infinite line AB:
144 : //
145 : // t = ((P - A) · (B - A)) / |B - A|^2
146 : //
147 : // If t ∈ [0, 1], the perpendicular projection lies within the segment.
148 : // If t < 0, the closest point is A.
149 : // If t > 1, the closest point is B.
150 : //
151 : // The resulting closest point x_int is then used to compute the Euclidean
152 : // distance from P to the segment, and the minimum distance over all edges is tracked.
153 2391208 : px = x(0);
154 2391208 : py = x(1);
155 2391208 : dx = x2 - x1;
156 2391208 : dy = y2 - y1;
157 2391208 : len2 = dx*dx + dy*dy;
158 2391208 : t = ((px - x1)*dx + (py - y1)*dy) / len2;
159 :
160 2391208 : if (t < 0.0)
161 : {
162 1098168 : x_int(0) = x1;
163 1098168 : x_int(1) = y1;
164 : }
165 1293040 : else if (t > 1.0)
166 : {
167 1059140 : x_int(0) = x2;
168 1059140 : x_int(1) = y2;
169 : }
170 : else
171 : {
172 233900 : x_int(0) = x1 + t*dx;
173 233900 : x_int(1) = y1 + t*dy;
174 : }
175 :
176 2391208 : dist_to_edge = dist(x(0), x(1), x_int(0), x_int(1));
177 2391208 : min_dist_to_edge = std::min(min_dist_to_edge, dist_to_edge);
178 :
179 : }
180 203288 : if (num_cross % 2 == 0)
181 : {
182 184512 : min_dist_to_edge = 1*min_dist_to_edge; // If the number of crossings is even, the point is inside the shape
183 : } else
184 : {
185 18776 : min_dist_to_edge = -1*min_dist_to_edge;
186 : }
187 :
188 203288 : phi_t(i,j,k) = 0.5 * (1.0 - std::tanh(min_dist_to_edge /(std::sqrt(2.0) * eps)));
189 :
190 : // Combine this polygon's field into either the running solid union or the running
191 : // void union (pointwise max), depending on whether it is marked as a void polygon.
192 274392 : if (is_void) phi_v(i,j,k) = is_first ? phi_t(i,j,k) : std::max(phi_v(i,j,k), phi_t(i,j,k));
193 574312 : else phi(i,j,k) = is_first ? phi_t(i,j,k) : std::max(phi(i,j,k), phi_t(i,j,k));
194 203288 : });
195 342 : }
196 :
197 46 : if (is_void) { first_void = false; any_void = true; }
198 38 : else first_solid = false;
199 46 : }
200 :
201 22 : if (any_void)
202 : // Punch holes into the solid field wherever a void polygon overlaps it: phi_void is zero inside
203 : // the void polygons and one outside them, so multiplying it in zeros out the void regions.
204 : {
205 20 : for (amrex::MFIter mfi(*a_phi[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
206 : {
207 12 : amrex::Box bx;
208 12 : amrex::IndexType type = a_phi[lev]->ixType();
209 24 : if (type == amrex::IndexType::TheCellType()) bx = mfi.growntilebox();
210 0 : else if (type == amrex::IndexType::TheNodeType()) bx = mfi.grownnodaltilebox();
211 0 : else Util::Abort(INFO, "Unkonwn index type");
212 12 : amrex::Array4<Set::Scalar> const& phi = a_phi[lev]->array(mfi);
213 12 : amrex::Array4<Set::Scalar> const& phi_v = phi_void[lev]->array(mfi);
214 :
215 12 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
216 : {
217 106656 : phi(i,j,k) = phi(i,j,k) * (1.0 - phi_v(i,j,k));
218 35552 : });
219 8 : }
220 : }
221 :
222 22 : if (invert_all)
223 : {
224 20 : for (amrex::MFIter mfi(*a_phi[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
225 : {
226 12 : amrex::Box bx;
227 12 : amrex::IndexType type = a_phi[lev]->ixType();
228 24 : if (type == amrex::IndexType::TheCellType()) bx = mfi.growntilebox();
229 0 : else if (type == amrex::IndexType::TheNodeType()) bx = mfi.grownnodaltilebox();
230 0 : else Util::Abort(INFO, "Unkonwn index type");
231 12 : amrex::Array4<Set::Scalar> const& phi = a_phi[lev]->array(mfi);
232 :
233 12 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
234 : {
235 71104 : phi(i,j,k) = 1 - phi(i,j,k);
236 35552 : });
237 8 : }
238 : }
239 22 : }
240 :
241 2 : static void Parse(PointList& value, IO::ParmParse& pp)
242 : {
243 2 : std::string filename;
244 2 : int verbose = 0;
245 : // Diffuseness of the solid boundary
246 10 : pp.query_default("eps", value.eps, "0.0", Unit::Length());
247 :
248 8 : pp.forbid("filename", "use file.name instead");
249 : // Location of .xy file
250 4 : pp.query_file("file.name", filename);
251 :
252 : // Verbosity (used in parser only)
253 2 : pp.query_default("verbose", verbose, 0);
254 :
255 : // Unitless coordinate multiplier, applied after the rotation
256 8 : pp.query_default("mult",value.mult, "1.0", Unit::Less());
257 :
258 : // Whether to invert the finished field as a whole (make IC 1 outside of solid instead of inside solid)
259 4 : pp.query_default("invert_all",value.invert_all,false);
260 :
261 : // X-offsed, applied after the rotation
262 8 : pp.queryarr_default("x0",value.x0,"0.0 0.0 0.0", Unit::Length());
263 :
264 : // Center of rigid-body rotation
265 8 : pp.queryarr_default("rotation.center",value.rotation.center,"0.0 0.0 0.0", Unit::Length());
266 :
267 : // Amount of rotation about rotation center (clockwise)
268 8 : pp.query_default("rotation.angle",value.rotation.angle,"0.0", Unit::Angle());
269 :
270 : // ObjNum jump (between consecutive points in the file) larger than this starts a new polygon
271 8 : pp.query_default("obj_num_threshold",value.objnum_threshold,"0.5", Unit::Less());
272 :
273 :
274 2 : Unit unit = Unit::Length();
275 : // Units of length in the file
276 2 : pp.queryunit("file.unit",unit);
277 14 : Util::AssertException( INFO,TEST(unit.isType(Unit::Length())),
278 4 : "Unit must be of type length but got unit ",unit.normalized_unitstring());
279 :
280 2 : std::ifstream datafile(filename);
281 2 : std::string line;
282 2 : if (datafile.is_open())
283 : {
284 2 : value.polygons.clear();
285 2 : bool have_prev_obj_num = false;
286 2 : Set::Scalar prev_obj_num = 0.0;
287 :
288 47 : while (getline(datafile, line))
289 : {
290 45 : std::istringstream in(line);
291 :
292 45 : std::string strx, stry, strz, strObjNum;
293 45 : in >> strx >> stry >> strz;
294 :
295 45 : Set::Scalar x = (std::stod(strx) * unit).normalized_value();
296 45 : Set::Scalar y = (std::stod(stry) * unit).normalized_value();
297 : #if AMREX_SPACEDIM > 2
298 0 : Set::Scalar z = (std::stod(strz) * unit).normalized_value();
299 : #endif
300 :
301 : // ObjNum is optional; if the file has no fourth column, treat every point as one polygon
302 45 : Set::Scalar ObjNum = have_prev_obj_num ? prev_obj_num : 0.0;
303 45 : if (in >> strObjNum) ObjNum = std::stod(strObjNum);
304 :
305 45 : Set::Vector X(AMREX_D_DECL(x,y,z));
306 :
307 45 : Set::Scalar cx = value.rotation.center[0];
308 45 : Set::Scalar cy = value.rotation.center[1];
309 :
310 45 : Set::Scalar theta = value.rotation.angle;
311 :
312 45 : Set::Scalar cos_t = std::cos(theta);
313 45 : Set::Scalar sin_t = std::sin(theta);
314 :
315 45 : Set::Scalar x_shifted = x - cx;
316 45 : Set::Scalar y_shifted = y - cy;
317 :
318 45 : X[0] = cx + x_shifted * cos_t - y_shifted * sin_t;
319 45 : X[1] = cy + x_shifted * sin_t + y_shifted * cos_t;
320 :
321 45 : X = value.x0 + value.mult*X;
322 :
323 : // A jump in ObjNum larger than the threshold marks the start of a new polygon
324 45 : if (!have_prev_obj_num || std::abs(ObjNum - prev_obj_num) > value.objnum_threshold)
325 : {
326 5 : value.polygons.push_back(std::vector<Set::Vector>());
327 : }
328 45 : have_prev_obj_num = true;
329 45 : prev_obj_num = ObjNum;
330 :
331 45 : value.polygons.back().push_back(X);
332 45 : if (verbose > 0)
333 0 : Util::Message(INFO, "x=", value.polygons.back().back().transpose());
334 45 : }
335 2 : datafile.close();
336 : }
337 : else
338 : {
339 0 : Util::Abort(INFO, "Unable to open file ", filename);
340 : }
341 :
342 : // Per-polygon void flag: nonzero marks a polygon as a void/negative region that is cut out of
343 : // the solid field instead of being unioned into it. Must have one entry per polygon if given.
344 4 : pp.queryarr("invert", value.invert);
345 :
346 14 : Util::AssertException( INFO, TEST(value.invert.size() == value.polygons.size()),
347 2 : "'invert' must have one entry per polygon: expected ", value.polygons.size(),
348 2 : " but got ", value.invert.size());
349 2 : }
350 :
351 : private:
352 : std::vector<std::vector<Set::Vector>> polygons; // Each entry is the ordered list of vertices for one closed polygon
353 : Set::Scalar eps;
354 : Set::Scalar mult = 1.0;
355 : Set::Vector x0 = Set::Vector::Zero();
356 : Set::Field<Set::Scalar> phi_temp;
357 : Set::Field<Set::Scalar> phi_void;
358 : Set::Scalar objnum_threshold = 0.5;
359 : std::vector<int> invert; // per-polygon void flag: nonzero = void/negative region
360 : bool invert_all = false;
361 :
362 : struct{
363 : Set::Vector center;
364 : Set::Scalar angle;
365 : } rotation;
366 :
367 2391208 : double dist(double x1, double y1, double x2, double y2)
368 : {
369 : // Find distance between 2 points
370 2391208 : double dx = x2 - x1;
371 2391208 : double dy = y2 - y1;
372 :
373 2391208 : return std::sqrt(dx*dx + dy*dy);
374 : }
375 : };
376 : }
377 : #endif
|
