Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Functions
Garfield::Polygon Namespace Reference

Functions

void Inside (const std::vector< double > &xpl, const std::vector< double > &ypl, const double x, const double y, bool &inside, bool &edge)
 
double Area (const std::vector< double > &xp, const std::vector< double > &yp)
 Determine the (signed) area of a polygon.
 
bool NonTrivial (const std::vector< double > &xp, const std::vector< double > &yp)
 Check whether a set of points builds a non-trivial polygon.
 
void EliminateButterflies (std::vector< double > &xp, std::vector< double > &yp, std::vector< double > &zp)
 

Function Documentation

◆ Area()

double Garfield::Polygon::Area ( const std::vector< double > &  xp,
const std::vector< double > &  yp 
)

Determine the (signed) area of a polygon.

Definition at line 274 of file Polygon.cc.

274 {
275
276 double f = 0.;
277 const unsigned int n = xp.size();
278 for (unsigned int i = 0; i < n; ++i) {
279 const unsigned int ii = i < n - 1 ? i + 1 : 0;
280 f += xp[i] * yp[ii] - xp[ii] * yp[i];
281 }
282 return 0.5 * f;
283}

Referenced by Garfield::ComponentNeBem2d::AddRegion().

◆ EliminateButterflies()

void Garfield::Polygon::EliminateButterflies ( std::vector< double > &  xp,
std::vector< double > &  yp,
std::vector< double > &  zp 
)

Try to eliminate "butterflies" (crossing of two adjacent segments
of a polygon), by point exchanges.

Definition at line 355 of file Polygon.cc.

356 {
357 //----------------------------------------------------------------------
358 // BUTFLD - Tries to eliminate "butterflies", i.e. the crossing of 2
359 // adjacent segments of a polygon, by point exchanges.
360 //----------------------------------------------------------------------
361 unsigned int np = xp.size();
362 if (np <= 3) return;
363 // Compute range.
364 const double xmin = *std::min_element(std::begin(xp), std::end(xp));
365 const double xmax = *std::max_element(std::begin(xp), std::end(xp));
366 const double ymin = *std::min_element(std::begin(yp), std::end(yp));
367 const double ymax = *std::max_element(std::begin(yp), std::end(yp));
368 const double zmin = *std::min_element(std::begin(zp), std::end(zp));
369 const double zmax = *std::max_element(std::begin(zp), std::end(zp));
370 // Set tolerances.
371 const double epsx = std::max(1.e-10 * std::abs(xmax - xmin), 1.e-6);
372 const double epsy = std::max(1.e-10 * std::abs(ymax - ymin), 1.e-6);
373 const double epsz = std::max(1.e-10 * std::abs(zmax - zmin), 1.e-6);
374 double xsurf = 0.;
375 double ysurf = 0.;
376 double zsurf = 0.;
377 for (unsigned int i = 2; i < np; ++i) {
378 const double x1 = xp[i - 1] - xp[0];
379 const double y1 = yp[i - 1] - yp[0];
380 const double z1 = zp[i - 1] - zp[0];
381 xsurf += std::abs((yp[i] - yp[0]) * z1 - y1 * (zp[i] - zp[0]));
382 ysurf += std::abs((xp[i] - xp[0]) * z1 - x1 * (zp[i] - zp[0]));
383 zsurf += std::abs((xp[i] - xp[0]) * y1 - x1 * (yp[i] - yp[0]));
384 }
385 // Eliminate points appearing twice, initialise marks.
386 std::vector<bool> mark(np, false);
387 // Scan the list.
388 for (unsigned int i = 0; i < np; ++i) {
389 if (mark[i]) continue;
390 for (unsigned int j = i + 1; j < np; ++j) {
391 if (std::abs(xp[i] - xp[j]) <= epsx && std::abs(yp[i] - yp[j]) <= epsy &&
392 std::abs(zp[i] - zp[j]) <= epsz) {
393 mark[j] = true;
394 }
395 }
396 }
397 // And remove the duplicate points.
398 unsigned int nNew = 0;
399 for (unsigned int i = 0; i < np; ++i) {
400 if (mark[i]) continue;
401 xp[nNew] = xp[i];
402 yp[nNew] = yp[i];
403 zp[nNew] = zp[i];
404 ++nNew;
405 }
406 // std::cout << "ElminateButterflies: old/new number of points: " << np
407 // << "/" << nNew << "\n";
408 // Update the number of points.
409 np = nNew;
410 xp.resize(np);
411 yp.resize(np);
412 zp.resize(np);
413 // No risk of having a butterfly with less than 4 points.
414 if (np <= 3) return;
415 // Select the axis with the largest norm.
416 unsigned int iaxis = 0;
417 if (xsurf > ysurf && xsurf > zsurf) {
418 iaxis = 1;
419 } else if (ysurf > zsurf) {
420 iaxis = 2;
421 } else {
422 iaxis = 3;
423 }
424 // Set number of passes to avoid endless loop.
425 unsigned int nPass = 0;
426 bool repass = true;
427 while (repass) {
428 // Make a pass.
429 ++nPass;
430 repass = false;
431 for (unsigned int i = 0; i < np; ++i) {
432 const unsigned int ii = (i + 1) % np;
433 for (unsigned int j = i + 2; j < np; ++j) {
434 const unsigned int jj = (j + 1) % np;
435 if (j + 1 >= np && jj >= i) continue;
436 // Check for a crossing.
437 if (iaxis == 1 && !Crossing(yp[i], zp[i], yp[ii], zp[ii],
438 yp[j], zp[j], yp[jj], zp[jj])) {
439 continue;
440 } else if (iaxis == 2 && !Crossing(xp[i], zp[i], xp[ii], zp[ii],
441 xp[j], zp[j], xp[jj], zp[jj])) {
442 continue;
443 } else if (iaxis == 3 && !Crossing(xp[i], yp[i], xp[ii], yp[ii],
444 xp[j], yp[j], xp[jj], yp[jj])) {
445 continue;
446 }
447 // If there is a crossing, exchange the portion in between.
448 for (unsigned int k = 0; k < (j - i) / 2; ++k) {
449 const unsigned int k1 = (i + k + 1) % np;
450 const unsigned int k2 = (j - k) % np;
451 std::swap(xp[k1], xp[k2]);
452 std::swap(yp[k1], yp[k2]);
453 std::swap(zp[k1], zp[k2]);
454 }
455 // And remember we have to do another pass after this.
456 repass = true;
457 }
458 }
459 // See whether we have to do another pass.
460 if (repass && nPass > np) {
461 std::cerr << "Butterfly: unable to eliminate the internal crossings;\n"
462 << " plot probably incorrect.\n";
463 break;
464 }
465 }
466}

Referenced by Garfield::SolidBox::Cut(), Garfield::SolidExtrusion::Cut(), Garfield::SolidRidge::Cut(), Garfield::SolidSphere::Cut(), and Garfield::SolidTube::Cut().

◆ Inside()

void Garfield::Polygon::Inside ( const std::vector< double > &  xpl,
const std::vector< double > &  ypl,
const double  x,
const double  y,
bool &  inside,
bool &  edge 
)

Determine whether the point (x, y) is located inside of the polygon (xpl, ypl).

Definition at line 187 of file Polygon.cc.

188 {
189 // Initial settings.
190 inside = false;
191 edge = false;
192 const unsigned int npl = xpl.size();
193 if (ypl.size() != npl) return;
194 // Special treatment for few points.
195 if (npl < 2) {
196 return;
197 } else if (npl == 2) {
198 edge = OnLine(xpl[0], ypl[0], xpl[1], ypl[1], x, y);
199 return;
200 }
201 // Determine the range of the data.
202 const double xmin = *std::min_element(std::begin(xpl), std::end(xpl));
203 const double xmax = *std::max_element(std::begin(xpl), std::end(xpl));
204 const double ymin = *std::min_element(std::begin(ypl), std::end(ypl));
205 const double ymax = *std::max_element(std::begin(ypl), std::end(ypl));
206
207 // Set tolerances.
208 double epsx = 1.e-8 * std::max(fabs(xmin), fabs(xmax));
209 double epsy = 1.e-8 * std::max(fabs(ymin), fabs(ymax));
210 epsx = std::max(epsx, 1.e-8);
211 epsy = std::max(epsy, 1.e-8);
212
213 // Ensure that we have a range.
214 if (fabs(xmax - xmin) <= epsx) {
215 if (y >= ymin - epsy && y <= ymax + epsy &&
216 fabs(xmax + xmin - 2 * x) <= epsx) {
217 edge = true;
218 } else {
219 edge = false;
220 }
221 } else if (fabs(ymax - ymin) <= epsy) {
222 if (x >= xmin - epsx && x <= xmax + epsx &&
223 fabs(ymax + ymin - 2 * y) <= epsy) {
224 edge = true;
225 } else {
226 edge = false;
227 }
228 }
229 // Choose a point at "infinity".
230 double xinf = xmin - fabs(xmax - xmin);
231 double yinf = ymin - fabs(ymax - ymin);
232
233 unsigned int nIter = 0;
234 bool ok = false;
235 while (!ok && nIter < 100) {
236 ok = true;
237 // Loop over the edges counting intersections.
238 unsigned int nCross = 0;
239 for (unsigned int j = 0; j < npl; ++j) {
240 const unsigned int jj = j < npl - 1 ? j + 1 : 0;
241 // Flag points located on one of the edges.
242 if (OnLine(xpl[j], ypl[j], xpl[jj], ypl[jj], x, y)) {
243 edge = true;
244 return;
245 }
246 // Count mid-line intersects.
247 double xc = 0., yc = 0.;
248 if (Crossing(x, y, xinf, yinf, xpl[j], ypl[j], xpl[jj], ypl[jj], xc,
249 yc)) {
250 ++nCross;
251 }
252 // Ensure that the testing line doesn't cross a corner.
253 if (OnLine(x, y, xinf, yinf, xpl[j], ypl[j])) {
254 xinf = xmin - ::Garfield::RndmUniform() * fabs(xmax - xinf);
255 yinf = ymin - ::Garfield::RndmUniform() * fabs(ymax - yinf);
256 ok = false;
257 break;
258 }
259 }
260 if (ok) {
261 // Set the INSIDE flag.
262 if (nCross != 2 * (nCross / 2)) inside = true;
263 return;
264 }
265 ++nIter;
266 }
267
268 std::cerr << "Polygon::Inside:\n Warning. Unable to verify "
269 << "whether a point is internal; setting to edge.\n";
270 inside = false;
271 edge = true;
272}
Garfield::RndmUniform
double RndmUniform()
Draw a random number uniformly distributed in the range [0, 1).
Definition: Random.hh:14
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

Referenced by Garfield::ComponentNeBem2d::GetMedium(), Garfield::SolidExtrusion::IsInside(), Garfield::SolidHole::IsInside(), and Garfield::SolidTube::IsInside().

◆ NonTrivial()

bool Garfield::Polygon::NonTrivial ( const std::vector< double > &  xp,
const std::vector< double > &  yp 
)

Check whether a set of points builds a non-trivial polygon.

Definition at line 285 of file Polygon.cc.

286 {
287
288 // -----------------------------------------------------------------------
289 // PLACHK - Checks whether a set of points builds a non-trivial
290 // polygon in the (x,y) plane.
291 // -----------------------------------------------------------------------
292
293 // First check number of points.
294 if (xp.size() != yp.size()) return false;
295 const unsigned int np = xp.size();
296 if (xp.size() < 3) return false;
297 // Find a second point at maximum distance of the first.
298 double d1 = 0.;
299 double x1 = 0.;
300 double y1 = 0.;
301 unsigned int i1 = 0;
302 double xmin = xp[0];
303 double ymin = yp[0];
304 double xmax = xp[0];
305 double ymax = yp[0];
306 for (unsigned int i = 1; i < np; ++i) {
307 xmin = std::min(xmin, xp[i]);
308 ymin = std::min(ymin, yp[i]);
309 xmax = std::max(xmax, xp[i]);
310 ymax = std::max(ymax, yp[i]);
311 const double dx = xp[i] - xp[0];
312 const double dy = yp[i] - yp[0];
313 const double d = dx * dx + dy * dy;
314 if (d > d1) {
315 x1 = dx;
316 y1 = dy;
317 i1 = i;
318 d1 = d;
319 }
320 }
321 // Set tolerances.
322 double epsx = 1.e-6 * (std::abs(xmax) + std::abs(xmin));
323 double epsy = 1.e-6 * (std::abs(ymax) + std::abs(ymin));
324 if (epsx <= 0) epsx = 1.e-6;
325 if (epsy <= 0) epsy = 1.e-6;
326 // See whether there is a range at all.
327 if (std::abs(xmax - xmin) <= epsx && std::abs(ymax - ymin) <= epsy) {
328 // Single point.
329 return false;
330 }
331 // See whether there is a second point.
332 if (d1 <= epsx * epsx + epsy * epsy || i1 == 0) return false;
333
334 // Find a third point maximising the external product.
335 double d2 = 0.;
336 unsigned int i2 = 0;
337 for (unsigned int i = 1; i < np; ++i) {
338 if (i == i1) continue;
339 const double dx = xp[i] - xp[0];
340 const double dy = yp[i] - yp[0];
341 const double d = std::abs(x1 * dy - y1 * dx);
342 if (d > d2) {
343 d2 = d;
344 i2 = i;
345 }
346 }
347 if (d2 <= epsx * epsy || i2 <= 0) {
348 // No third point.
349 return false;
350 }
351 // Seems to be OK.
352 return true;
353}

Referenced by Garfield::SolidExtrusion::SetProfile().