Line data Source code
1 : //
2 : // Fill a domain based on a list of points to form a closed 2D shape.
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 for a single closed polgyon
5 :
6 : #ifndef IC_POINTLIST_H_
7 : #define IC_POINTLIST_H_
8 :
9 : #include "Set/Set.H"
10 : #include "IC/IC.H"
11 : #include "Util/Util.H"
12 : using namespace std;
13 : #include <iostream>
14 : #include <vector>
15 : #include <algorithm>
16 : #include <iterator>
17 : #include <mpi.h>
18 : #include "IO/ParmParse.H"
19 :
20 : namespace IC
21 : {
22 : class PointList : public IC<Set::Scalar>
23 : {
24 : public:
25 : static constexpr const char* name = "pointlist";
26 :
27 : enum Type
28 : {
29 : Partition,
30 : Values
31 : };
32 :
33 : PointList(amrex::Vector<amrex::Geometry>& _geom) : IC(_geom) {}
34 :
35 1 : PointList(amrex::Vector<amrex::Geometry>& _geom, IO::ParmParse& pp, std::string name) : IC(_geom)
36 : {
37 1 : pp.queryclass(name, *this);
38 1 : }
39 : void Define() {
40 :
41 : };
42 :
43 14 : void Add(const int& lev, Set::Field<Set::Scalar>& a_phi, Set::Scalar)
44 : {
45 262 : for (amrex::MFIter mfi(*a_phi[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
46 : {
47 : if (AMREX_SPACEDIM != 2)
48 : {
49 : amrex::Abort("This code only supports 2D (AMREX_SPACEDIM must be 2)");
50 : }
51 :
52 248 : amrex::Box bx;
53 248 : amrex::IndexType type = a_phi[lev]->ixType();
54 496 : if (type == amrex::IndexType::TheCellType()) bx = mfi.growntilebox();
55 0 : else if (type == amrex::IndexType::TheNodeType()) bx = mfi.grownnodaltilebox();
56 0 : else Util::Abort(INFO, "Unkonwn index type");
57 248 : amrex::Array4<Set::Scalar> const& phi = a_phi[lev]->array(mfi);
58 :
59 248 : if (X.back() != X.front())
60 : // The first and last point must match, if they do not add another point to the end
61 : {
62 1 : X.push_back(X.front());
63 : }
64 :
65 248 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
66 : {
67 61080 : Set::Vector x = Set::Position(i, j, k, geom[lev], type);
68 : Set::Scalar x_cross;
69 61080 : Set::Scalar min_dist_to_edge = std::numeric_limits<Set::Scalar>::max(); // Mininmum distance from point to nearest edge
70 : Set::Scalar dist_to_edge; // Distance from point to nearest edge
71 : Set::Scalar x1;
72 : Set::Scalar y1;
73 : Set::Scalar x2;
74 : Set::Scalar y2;
75 61080 : int num_cross = 0;
76 :
77 : Set::Scalar px;
78 : Set::Scalar py;
79 : Set::Scalar dx;
80 : Set::Scalar dy;
81 : Set::Scalar len2;
82 : Set::Scalar t;
83 61080 : Set::Vector x_int; // x and y position of the intercept of the mesh point and the edge, assuming edge is infinite
84 :
85 1954560 : for (unsigned int n = 0; n < X.size()-1; n++)
86 : {
87 :
88 1893480 : x1 = X[n](0);
89 1893480 : y1 = X[n](1);
90 1893480 : x2 = X[n+1](0);
91 1893480 : y2 = X[n+1](1);
92 :
93 1893480 : if ((y1 > x(1)) != (y2 > x(1)))
94 : {
95 94680 : x_cross = x1 + (x(1) - y1) * (x2 - x1) / (y2 - y1);
96 :
97 94680 : if (x_cross > x(0))
98 45592 : num_cross++;
99 : }
100 :
101 : // Computes the closest point on a line segment (X[n] → X[n+1]) to a query point x,
102 : // using vector projection onto the segment.
103 : //
104 : // Let A = X[n], B = X[n+1], and P = x.
105 : //
106 : // The parameter t is the normalized projection of P onto the infinite line AB:
107 : //
108 : // t = ((P - A) · (B - A)) / |B - A|^2
109 : //
110 : // If t ∈ [0, 1], the perpendicular projection lies within the segment.
111 : // If t < 0, the closest point is A.
112 : // If t > 1, the closest point is B.
113 : //
114 : // The resulting closest point x_int is then used to compute the Euclidean
115 : // distance from P to the segment, and the minimum distance over all edges is tracked.
116 1893480 : px = x(0);
117 1893480 : py = x(1);
118 1893480 : dx = x2 - x1;
119 1893480 : dy = y2 - y1;
120 1893480 : len2 = dx*dx + dy*dy;
121 1893480 : t = ((px - x1)*dx + (py - y1)*dy) / len2;
122 :
123 1893480 : if (t < 0.0)
124 : {
125 904728 : x_int(0) = x1;
126 904728 : x_int(1) = y1;
127 : }
128 988752 : else if (t > 1.0)
129 : {
130 881024 : x_int(0) = x2;
131 881024 : x_int(1) = y2;
132 : }
133 : else
134 : {
135 107728 : x_int(0) = x1 + t*dx;
136 107728 : x_int(1) = y1 + t*dy;
137 : }
138 :
139 1893480 : dist_to_edge = dist(x(0), x(1), x_int(0), x_int(1));
140 1893480 : min_dist_to_edge = std::min(min_dist_to_edge, dist_to_edge);
141 :
142 : }
143 61080 : if (num_cross % 2 == 0)
144 : {
145 52280 : min_dist_to_edge = 1*min_dist_to_edge; // If the number of crossings is even, the point is inside the shape
146 : } else
147 : {
148 8800 : min_dist_to_edge = -1*min_dist_to_edge;
149 : }
150 :
151 61080 : phi(i,j,k) = 0.5 * (1.0 - std::tanh(min_dist_to_edge /(std::sqrt(2.0) * eps)));
152 :
153 61080 : if (invert) {
154 0 : phi(i,j,k) = 1 - phi(i,j,k);
155 : }
156 :
157 61080 : });
158 14 : }
159 14 : }
160 :
161 1 : static void Parse(PointList& value, IO::ParmParse& pp)
162 : {
163 1 : std::string filename;
164 1 : int verbose = 0;
165 : // Diffuseness of the solid boundary
166 5 : pp.query_default("eps", value.eps, "0.0", Unit::Length());
167 :
168 4 : pp.forbid("filename", "use file.name instead");
169 : // Location of .xy file
170 2 : pp.query_file("file.name", filename);
171 :
172 : // Verbosity (used in parser only)
173 1 : pp.query_default("verbose", verbose, 0);
174 :
175 : // Unitless coordinate multiplier, applied after the rotation
176 4 : pp.query_default("mult",value.mult, "1.0", Unit::Less());
177 :
178 : // Whether to invert (make IC 1 outside of solid instead of inside solid)
179 2 : pp.query_default("invert",value.invert,false);
180 :
181 : // X-offsed, applied after the rotation
182 4 : pp.queryarr_default("x0",value.x0,"0.0 0.0 0.0", Unit::Length());
183 :
184 : // if (pp.contains("x0"))
185 : // pp_queryarr("x0",value.x0); // Coordinate offset
186 :
187 : // Center of rigid-body rotation
188 4 : pp.queryarr_default("rotation.center",value.rotation.center,"0.0 0.0 0.0", Unit::Length());
189 :
190 : // Amount of rotation about rotation center (clockwise)
191 4 : pp.query_default("rotation.angle",value.rotation.angle,"0.0", Unit::Angle());
192 :
193 :
194 1 : Unit unit = Unit::Length();
195 : // Units of length in the file
196 1 : pp.queryunit("file.unit",unit);
197 7 : Util::AssertException( INFO,TEST(unit.isType(Unit::Length())),
198 2 : "Unit must be of type length but got unit ",unit.normalized_unitstring());
199 :
200 1 : std::ifstream datafile(filename);
201 1 : std::string line;
202 1 : if (datafile.is_open())
203 : {
204 1 : value.X.clear();
205 :
206 32 : while (getline(datafile, line))
207 : {
208 31 : std::istringstream in(line);
209 :
210 31 : std::string strx, stry, strz, strR;
211 31 : in >> strx >> stry >> strz >> strR;
212 :
213 31 : Set::Scalar x = (std::stod(strx) * unit).normalized_value();
214 31 : Set::Scalar y = (std::stod(stry) * unit).normalized_value();
215 : #if AMREX_SPACEDIM > 2
216 0 : Set::Scalar z = (std::stod(strz) * unit).normalized_value();
217 : #endif
218 :
219 31 : Set::Vector X(AMREX_D_DECL(x,y,z));
220 :
221 31 : Set::Scalar cx = value.rotation.center[0];
222 31 : Set::Scalar cy = value.rotation.center[1];
223 :
224 31 : Set::Scalar theta = value.rotation.angle;
225 :
226 31 : Set::Scalar cos_t = std::cos(theta);
227 31 : Set::Scalar sin_t = std::sin(theta);
228 :
229 31 : Set::Scalar x_shifted = x - cx;
230 31 : Set::Scalar y_shifted = y - cy;
231 :
232 31 : X[0] = cx + x_shifted * cos_t - y_shifted * sin_t;
233 31 : X[1] = cy + x_shifted * sin_t + y_shifted * cos_t;
234 :
235 31 : X = value.x0 + value.mult*X;
236 :
237 31 : value.X.push_back(X);
238 31 : if (verbose > 0)
239 0 : Util::Message(INFO, "x=", value.X.back().transpose());
240 31 : }
241 1 : datafile.close();
242 : }
243 : else
244 : {
245 0 : Util::Abort(INFO, "Unable to open file ", filename);
246 : }
247 1 : }
248 :
249 : private:
250 : std::vector<Set::Vector> X; // Big X is a vector of vectors of each coordinate read into the file
251 : Set::Scalar eps;
252 : Set::Scalar mult = 1.0;
253 : Set::Vector x0 = Set::Vector::Zero();
254 : bool invert = false;
255 :
256 : struct{
257 : Set::Vector center;
258 : Set::Scalar angle;
259 : } rotation;
260 :
261 1893480 : double dist(double x1, double y1, double x2, double y2)
262 : {
263 : // Find distance between 2 points
264 1893480 : double dx = x2 - x1;
265 1893480 : double dy = y2 - y1;
266 :
267 1893480 : return std::sqrt(dx*dx + dy*dy);
268 : }
269 : };
270 : }
271 : #endif
|
