20 #include <gsl/gsl_sys.h>
24 #define SQR(x) ((x)*(x))
25 #define PointID(triangle_id, vertex_id) Triangles_m[triangle_id][vertex_id]
26 #define Point(triangle_id, vertex_id) Points_m[Triangles_m[triangle_id][vertex_id]]
45 return x1 (0) < x2 (0);
51 return x1 (1) < x2 (1);
57 return x1 (2) < x2 (2);
64 Vector_t get_max_extent (std::vector<Vector_t>& coords) {
66 coords.begin (), coords.end (), VectorLessX ());
68 coords.begin (), coords.end (), VectorLessY ());
70 coords.begin (), coords.end (), VectorLessZ ());
71 return Vector_t (x (0), y (1), z (2));
78 Vector_t get_min_extent (std::vector<Vector_t>& coords) {
80 coords.begin (), coords.end (), VectorLessX ());
82 coords.begin (), coords.end (), VectorLessY ());
84 coords.begin (), coords.end (), VectorLessZ ());
85 return Vector_t (x (0), y (1), z (2));
91 static void write_voxel_mesh (
92 const std::unordered_map<
int, std::unordered_set<int> >& ids,
98 const size_t numpoints = 8 * ids.size ();
100 of.open (std::string (
"data/testBBox.vtk").c_str ());
101 assert (of.is_open ());
104 of <<
"# vtk DataFile Version 2.0" <<
std::endl;
105 of <<
"generated using BoundaryGeometry::computeMeshVoxelization"
108 of <<
"DATASET UNSTRUCTURED_GRID" <<
std::endl;
109 of <<
"POINTS " << numpoints <<
" float" <<
std::endl;
111 const auto nr0_times_nr1 = nr[0] * nr[1];
112 for (
auto& elem: ids) {
113 auto id = elem.first;
114 int k = (
id - 1) / nr0_times_nr1;
115 int rest = (
id - 1) % nr0_times_nr1;
116 int j = rest / nr[0];
117 int i = rest % nr[0];
120 P[0] = i * hr_m[0] + origin[0];
121 P[1] = j * hr_m[1] + origin[1];
122 P[2] = k * hr_m[2] + origin[2];
124 of << P[0] <<
" " << P[1] <<
" " << P[2] <<
std::endl;
125 of << P[0] + hr_m[0] <<
" " << P[1] <<
" " << P[2] <<
std::endl;
126 of << P[0] <<
" " << P[1] + hr_m[1] <<
" " << P[2] <<
std::endl;
127 of << P[0] + hr_m[0] <<
" " << P[1] + hr_m[1] <<
" " << P[2] <<
std::endl;
128 of << P[0] <<
" " << P[1] <<
" " << P[2] + hr_m[2] <<
std::endl;
129 of << P[0] + hr_m[0] <<
" " << P[1] <<
" " << P[2] + hr_m[2] <<
std::endl;
130 of << P[0] <<
" " << P[1] + hr_m[1] <<
" " << P[2] + hr_m[2] <<
std::endl;
131 of << P[0] + hr_m[0] <<
" " << P[1] + hr_m[1] <<
" " << P[2] + hr_m[2] <<
std::endl;
134 const auto num_cells = ids.size ();
135 of <<
"CELLS " << num_cells <<
" " << 9 * num_cells <<
std::endl;
136 for (
size_t i = 0; i < num_cells; i++)
138 << 8 * i <<
" " << 8 * i + 1 <<
" " << 8 * i + 2 <<
" " << 8 * i + 3 <<
" "
139 << 8 * i + 4 <<
" " << 8 * i + 5 <<
" " << 8 * i + 6 <<
" " << 8 * i + 7 << std::endl;
140 of <<
"CELL_TYPES " << num_cells <<
std::endl;
141 for (
size_t i = 0; i < num_cells; i++)
142 of <<
"11" << std::endl;
143 of <<
"CELL_DATA " << num_cells <<
std::endl;
144 of <<
"SCALARS " <<
"cell_attribute_data" <<
" float " <<
"1" <<
std::endl;
145 of <<
"LOOKUP_TABLE " <<
"default" <<
std::endl;
146 for (
size_t i = 0; i < num_cells; i++)
147 of << (
float)(i) << std::endl;
149 of <<
"COLOR_SCALARS " <<
"BBoxColor " << 4 <<
std::endl;
150 for (
size_t i = 0; i < num_cells; i++) {
151 of <<
"1.0" <<
" 1.0 " <<
"0.0 " <<
"1.0" <<
std::endl;
169 #define LERP( A, B, C) ((B)+(A)*((C)-(B)))
170 #define MIN2(a,b) (((a) < (b)) ? (a) : (b))
171 #define MAX2(a,b) (((a) > (b)) ? (a) : (b))
172 #define MIN3(a,b,c) ((((a)<(b))&&((a)<(c))) ? (a) : (((b)<(c)) ? (b) : (c)))
173 #define MAX3(a,b,c) ((((a)>(b))&&((a)>(c))) ? (a) : (((b)>(c)) ? (b) : (c)))
189 inline double v1(
int i)
const {
195 inline double v2(
int i)
const {
201 inline double v3(
int i)
const {
210 pts[0][0] *= scaleby[0];
211 pts[0][1] *= scaleby[1];
212 pts[0][2] *= scaleby[2];
213 pts[1][0] *= scaleby[0];
214 pts[1][1] *= scaleby[1];
215 pts[1][2] *= scaleby[2];
216 pts[2][0] *= scaleby[0];
217 pts[2][1] *= scaleby[1];
218 pts[2][2] *= scaleby[2];
236 int outcode_fcmp = 0;
238 if (gsl_fcmp (p[0], 0.5,
EPS) > 0) outcode_fcmp |= 0x01;
239 if (gsl_fcmp (p[0],-0.5,
EPS) < 0) outcode_fcmp |= 0x02;
240 if (gsl_fcmp (p[1], 0.5,
EPS) > 0) outcode_fcmp |= 0x04;
241 if (gsl_fcmp (p[1],-0.5,
EPS) < 0) outcode_fcmp |= 0x08;
242 if (gsl_fcmp (p[2], 0.5,
EPS) > 0) outcode_fcmp |= 0x10;
243 if (gsl_fcmp (p[2],-0.5,
EPS) < 0) outcode_fcmp |= 0x20;
245 return(outcode_fcmp);
256 int outcode_fcmp = 0;
258 if (gsl_fcmp( p[0] + p[1], 1.0,
EPS) > 0) outcode_fcmp |= 0x001;
259 if (gsl_fcmp( p[0] - p[1], 1.0,
EPS) > 0) outcode_fcmp |= 0x002;
260 if (gsl_fcmp(-p[0] + p[1], 1.0,
EPS) > 0) outcode_fcmp |= 0x004;
261 if (gsl_fcmp(-p[0] - p[1], 1.0,
EPS) > 0) outcode_fcmp |= 0x008;
262 if (gsl_fcmp( p[0] + p[2], 1.0,
EPS) > 0) outcode_fcmp |= 0x010;
263 if (gsl_fcmp( p[0] - p[2], 1.0,
EPS) > 0) outcode_fcmp |= 0x020;
264 if (gsl_fcmp(-p[0] + p[2], 1.0,
EPS) > 0) outcode_fcmp |= 0x040;
265 if (gsl_fcmp(-p[0] - p[2], 1.0,
EPS) > 0) outcode_fcmp |= 0x080;
266 if (gsl_fcmp( p[1] + p[2], 1.0,
EPS) > 0) outcode_fcmp |= 0x100;
267 if (gsl_fcmp( p[1] - p[2], 1.0,
EPS) > 0) outcode_fcmp |= 0x200;
268 if (gsl_fcmp(-p[1] + p[2], 1.0,
EPS) > 0) outcode_fcmp |= 0x400;
269 if (gsl_fcmp(-p[1] - p[2], 1.0,
EPS) > 0) outcode_fcmp |= 0x800;
271 return(outcode_fcmp);
282 int outcode_fcmp = 0;
284 if (gsl_fcmp( p[0] + p[1] + p[2], 1.5,
EPS) > 0) outcode_fcmp |= 0x01;
285 if (gsl_fcmp( p[0] + p[1] - p[2], 1.5,
EPS) > 0) outcode_fcmp |= 0x02;
286 if (gsl_fcmp( p[0] - p[1] + p[2], 1.5,
EPS) > 0) outcode_fcmp |= 0x04;
287 if (gsl_fcmp( p[0] - p[1] - p[2], 1.5,
EPS) > 0) outcode_fcmp |= 0x08;
288 if (gsl_fcmp(-p[0] + p[1] + p[2], 1.5,
EPS) > 0) outcode_fcmp |= 0x10;
289 if (gsl_fcmp(-p[0] + p[1] - p[2], 1.5,
EPS) > 0) outcode_fcmp |= 0x20;
290 if (gsl_fcmp(-p[0] - p[1] + p[2], 1.5,
EPS) > 0) outcode_fcmp |= 0x40;
291 if (gsl_fcmp(-p[0] - p[1] - p[2], 1.5,
EPS) > 0) outcode_fcmp |= 0x80;
293 return(outcode_fcmp);
312 plane_point[0] =
LERP(alpha, p1[0], p2[0]);
313 plane_point[1] =
LERP(alpha, p1[1], p2[1]);
314 plane_point[2] =
LERP(alpha, p1[2], p2[2]);
315 return(face_plane(plane_point) & mask);
329 const int outcode_diff
331 if ((0x01 & outcode_diff) != 0)
332 if (check_point(p1,p2,( .5-p1[0])/(p2[0]-p1[0]),0x3e) ==
INSIDE)
return(
INSIDE);
333 if ((0x02 & outcode_diff) != 0)
334 if (check_point(p1,p2,(-.5-p1[0])/(p2[0]-p1[0]),0x3d) ==
INSIDE)
return(
INSIDE);
335 if ((0x04 & outcode_diff) != 0)
336 if (check_point(p1,p2,( .5-p1[1])/(p2[1]-p1[1]),0x3b) ==
INSIDE)
return(
INSIDE);
337 if ((0x08 & outcode_diff) != 0)
338 if (check_point(p1,p2,(-.5-p1[1])/(p2[1]-p1[1]),0x37) ==
INSIDE)
return(
INSIDE);
339 if ((0x10 & outcode_diff) != 0)
340 if (check_point(p1,p2,( .5-p1[2])/(p2[2]-p1[2]),0x2f) ==
INSIDE)
return(
INSIDE);
341 if ((0x20 & outcode_diff) != 0)
342 if (check_point(p1,p2,(-.5-p1[2])/(p2[2]-p1[2]),0x1f) ==
INSIDE)
return(
INSIDE);
355 return ((A[0] <
EPS) ? 4 : 0 | (A[0] > -
EPS) ? 32 : 0 |
356 (A[1] <
EPS) ? 2 : 0 | (A[1] > -
EPS) ? 16 : 0 |
357 (A[2] <
EPS) ? 1 : 0 | (A[2] > -
EPS) ? 8 : 0);
361 point_triangle_intersection (
386 const int sign12 = SIGN3(cross12_1p);
391 const int sign23 = SIGN3(cross23_2p);
396 const int sign31 = SIGN3(cross31_3p);
416 triangle_intersects_cube (
435 if ((v1_test & v2_test & v3_test) != 0)
return(
OUTSIDE);
440 v1_test |= bevel_2d(t.
v1()) << 8;
441 v2_test |= bevel_2d(t.
v2()) << 8;
442 v3_test |= bevel_2d(t.
v3()) << 8;
443 if ((v1_test & v2_test & v3_test) != 0)
return(
OUTSIDE);
448 v1_test |= bevel_3d(t.
v1()) << 24;
449 v2_test |= bevel_3d(t.
v2()) << 24;
450 v3_test |= bevel_3d(t.
v3()) << 24;
451 if ((v1_test & v2_test & v3_test) != 0)
return(
OUTSIDE);
461 if ((v1_test & v2_test) == 0)
463 if ((v1_test & v3_test) == 0)
465 if ((v2_test & v3_test) == 0)
495 double d = norm[0] * t.
v1(0) + norm[1] * t.
v1(1) + norm[2] * t.
v1(2);
502 if(
fabs(denom=(norm[0] + norm[1] + norm[2]))>
EPS) {
505 if (
fabs(hitpp[0]) <= 0.5)
506 if (point_triangle_intersection(hitpp,t) ==
INSIDE)
return(
INSIDE);
508 if(
fabs(denom=(norm[0] + norm[1] - norm[2]))>
EPS) {
510 hitpn[2] = -(hitpn[0] = hitpn[1] = d / denom);
511 if (
fabs(hitpn[0]) <= 0.5)
512 if (point_triangle_intersection(hitpn,t) ==
INSIDE)
return(
INSIDE);
514 if(
fabs(denom=(norm[0] - norm[1] + norm[2]))>
EPS) {
516 hitnp[1] = -(hitnp[0] = hitnp[2] = d / denom);
517 if (
fabs(hitnp[0]) <= 0.5)
518 if (point_triangle_intersection(hitnp,t) ==
INSIDE)
return(
INSIDE);
520 if(
fabs(denom=(norm[0] - norm[1] - norm[2]))>
EPS) {
522 hitnn[1] = hitnn[2] = -(hitnn[0] = d / denom);
523 if (
fabs(hitnn[0]) <= 0.5)
524 if (point_triangle_intersection(hitnn,t) ==
INSIDE)
return(
INSIDE);
590 pts[0][0] *= scale[0];
591 pts[0][1] *= scale[1];
592 pts[0][2] *= scale[2];
593 pts[1][0] *= scale[0];
594 pts[1][1] *= scale[1];
595 pts[1][2] *= scale[2];
608 if ( (tmin_ > tymax) || (tymin > tmax_) )
616 if ( (tmin_ > tzmax) || (tzmin > tmax_) )
642 t_.
scale (scaleby , v_.
pts[0] + 0.5);
643 return triangle_intersects_cube (t_);
665 static inline Vector_t normalVector (
671 const double magnitude =
sqrt (
SQR (N (0)) +
SQR (N (1)) +
SQR (N (2)));
672 assert (gsl_fcmp (magnitude, 0.0,
EPS) > 0);
673 return N / magnitude;
677 static inline double computeArea (
684 return(0.5 *
sqrt (
dot (AB, AB) *
dot (AC, AC) -
dot (AB, AC) *
dot (AB, AC)));
699 SIZE,
"GEOMETRY",
"The \"GEOMETRY\" statement defines the beam pipe geometry.") {
703 "Specifies the geometry file [h5fed]",
708 "BOX, BOXCORNER, ELLIPTIC if FGEOM is selected topo is over-written ",
713 "Specifies the length of a tube shaped elliptic beam pipe [m]",
718 "Specifies the start of a tube shaped elliptic beam pipe [m]",
723 "Specifies the major semi-axis of a tube shaped elliptic beam pipe [m]",
728 "Specifies the major semi-axis of a tube shaped elliptic beam pipe [m]",
733 "In case of BOXCORNER Specifies first part with hight == B [m]",
738 "In case of BOXCORNER Specifies first second with hight == B-C [m]",
743 "In case of BOXCORNER Specifies hight of corner C [m]",
748 "Distribution to be generated on the surface",
752 "Distribution array to be generated on the surface");
756 "Multiplicative scaling factor for coordinates ",
761 "Multiplicative scaling factor for X coordinates ",
766 "Multiplicative scaling factor for Y coordinates ",
771 "Multiplicative scaling factor for Z coordinates ",
776 "Shift in z direction",
780 (
"APERTURE",
"The element aperture");
801 randGen_m = gsl_rng_alloc(gsl_rng_default);
809 const std::string&
name,
813 randGen_m = gsl_rng_alloc(gsl_rng_default);
858 throw OpalException (
"BoundaryGeometry::find()",
"Geometry \""
859 + name +
"\" not found.");
868 const int triangle_id,
934 const int triangle_id,
950 const double a = -
dot(n,w0);
951 const double b =
dot(n,dir);
952 if (gsl_fcmp (b, 0.0,
EPS) == 0) {
961 const double r = a / b;
969 if (r < 0 || 1.0 < r) {
979 const double uu =
dot(u,u);
980 const double uv =
dot(u,v);
981 const double vv =
dot(v,v);
983 const double wu =
dot(w,u);
984 const double wv =
dot(w,v);
985 const double D = uv * uv - uu * vv;
988 const double s = (uv * wv - vv * wu) / D;
989 if (s < 0.0 || s > 1.0) {
992 const double t = (uv * wu - uu * wv) / D;
993 if (t < 0.0 || (s + t) > 1.0) {
999 }
else if ((0.0 <= r) && (r <= 1.0)) {
1006 static inline double magnitude (
1032 *gmsg <<
"* " << __func__ <<
": "
1033 <<
"reference_pt=" << reference_pt
1034 <<
", P=" << P <<
endl;
1038 const Vector_t v = reference_pt - P;
1046 int triangle_id = -1;
1048 for (
int i = 0; i < N; i++) {
1055 *gmsg <<
"* " << __func__ <<
": "
1056 <<
"result: " << result <<
endl;
1082 *gmsg <<
"* " << __func__ <<
": "
1100 int triangle_id = -1;
1103 I, triangle_id) > 0) ? 1 : 0;
1106 *gmsg <<
"* " << __func__ <<
": "
1107 <<
" result=" << result
1138 #define mapPoint2VoxelIndices(pt, i, j, k) { \
1139 i = floor ((pt[0] - voxelMesh_m.minExtent [0]) / voxelMesh_m.sizeOfVoxel[0]); \
1140 j = floor ((pt[1] - voxelMesh_m.minExtent [1]) / voxelMesh_m.sizeOfVoxel[1]); \
1141 k = floor ((pt[2] - voxelMesh_m.minExtent [2]) / voxelMesh_m.sizeOfVoxel[2]); \
1142 if (!(0 <= i && i < voxelMesh_m.nr_m[0] && \
1143 0 <= j && j < voxelMesh_m.nr_m[1] && \
1144 0 <= k && k < voxelMesh_m.nr_m[2])) { \
1145 *gmsg << "* " << __func__ << ":" \
1146 << " WARNING: pt=" << pt \
1147 << " is outside the bbox" \
1182 for (
unsigned int triangle_id = 0; triangle_id <
Triangles_m.size(); triangle_id++) {
1187 MIN3 (v1[0], v2[0], v3[0]),
1188 MIN3 (v1[1], v2[1], v3[1]),
1189 MIN3 (v1[2], v2[2], v3[2]) };
1191 MAX3 (v1[0], v2[0], v3[0]),
1192 MAX3 (v1[1], v2[1], v3[1]),
1193 MAX3 (v1[2], v2[2], v3[2]) };
1194 int i_min, j_min, k_min;
1195 int i_max, j_max, k_max;
1199 for (
int i = i_min; i <= i_max; i++) {
1200 for (
int j = j_min; j <= j_max; j++) {
1201 for (
int k = k_min; k <= k_max; k++) {
1211 *gmsg <<
"* Mesh voxelization done." <<
endl;
1235 double longest_edge_max_m = 0.0;
1236 for (
unsigned int i = 0; i < bg->
Triangles_m.size(); i++) {
1241 const double length_edge1 =
sqrt (
1242 SQR (x1[0] - x2[0]) +
SQR (x1[1] - x2[1]) +
SQR (x1[2] - x2[2]));
1243 const double length_edge2 =
sqrt (
1244 SQR (x3[0] - x2[0]) +
SQR (x3[1] - x2[1]) +
SQR (x3[2] - x2[2]));
1245 const double length_edge3 =
sqrt (
1246 SQR (x3[0] - x1[0]) +
SQR (x3[1] - x1[1]) +
SQR (x3[2] - x1[2]));
1248 double max = length_edge1;
1249 if (length_edge2 > max) max = length_edge2;
1250 if (length_edge3 > max) max = length_edge3;
1253 if (longest_edge_max_m < max) longest_edge_max_m =
max;
1410 static void computeTriangleNeighbors (
1412 std::vector<std::set<unsigned int>>& neighbors
1414 std::vector<std::set<unsigned int>> adjacencies_to_pt (bg->
Points_m.size());
1417 for (
unsigned int triangle_id = 0; triangle_id < bg->
Triangles_m.size(); triangle_id++) {
1418 for (
unsigned int j = 1; j <= 3; j++) {
1419 auto pt_id = bg->PointID (triangle_id, j);
1420 assert (pt_id < bg->
Points_m.size ());
1421 adjacencies_to_pt [pt_id].insert (triangle_id);
1425 for (
unsigned int triangle_id = 0; triangle_id < bg->
Triangles_m.size(); triangle_id++) {
1426 std::set<unsigned int> to_A = adjacencies_to_pt [bg->PointID (triangle_id, 1)];
1427 std::set<unsigned int> to_B = adjacencies_to_pt [bg->PointID (triangle_id, 2)];
1428 std::set<unsigned int> to_C = adjacencies_to_pt [bg->PointID (triangle_id, 3)];
1430 std::set<unsigned int> intersect;
1431 std::set_intersection (
1432 to_A.begin(), to_A.end(),
1433 to_B.begin(), to_B.end(),
1434 std::inserter(intersect,intersect.begin()));
1435 std::set_intersection(
1436 to_B.begin(), to_B.end(),
1437 to_C.begin(), to_C.end(),
1438 std::inserter(intersect,intersect.begin()));
1439 std::set_intersection(
1440 to_C.begin(), to_C.end(),
1441 to_A.begin(), to_A.end(),
1442 std::inserter(intersect, intersect.begin()));
1443 intersect.erase (triangle_id);
1445 neighbors [triangle_id] = intersect;
1447 *gmsg <<
"* " << __func__ <<
": Computing neighbors done" <<
endl;
1477 std::vector<Vector_t> intersection_points;
1480 for (
unsigned int triangle_id = 0; triangle_id < bg->
Triangles_m.size(); triangle_id++) {
1483 intersection_points.push_back (result);
1488 return ((intersection_points.size () % 2) == 1);
1492 static bool hasInwardPointingNormal (
1494 const int triangle_id
1499 const Vector_t triNormal = normalVector (A, B, C);
1508 const Vector_t P = (A+B+
C)/3 + triNormal * minvoxelmesh;
1518 const bool is_inside = isInside (bg, P);
1519 const double dotPA_N =
dot (P - A, triNormal);
1520 return (is_inside && dotPA_N >= 0) || (!is_inside && dotPA_N < 0);
1524 static void orientTriangle (
BoundaryGeometry* bg,
int ref_id,
int triangle_id) {
1529 for (
int i = 1; i <= 3; i++) {
1530 for (
int j = 1; j <= 3; j++) {
1531 if (bg->PointID (triangle_id, j) == bg->PointID (ref_id, i)) {
1535 if (n == 2)
goto edge_found;
1541 int diff =
id[1] -
id[0];
1542 if ((((ic[1] - ic[0]) == 1) && ((diff == 1) || (diff == -2))) ||
1543 (((ic[1] - ic[0]) == 2) && ((diff == -1) || (diff == 2)))) {
1544 std::swap (bg->PointID (triangle_id,
id[0]), bg->PointID (triangle_id,
id[1]));
1549 std::vector<std::set<unsigned int>> neighbors (bg->
Triangles_m.size());
1551 computeTriangleNeighbors (bg, neighbors);
1554 int triangle_id = 0;
1556 std::vector<unsigned int> triangles (bg->
Triangles_m.size());
1557 std::vector<unsigned int>::size_type queue_cursor = 0;
1558 std::vector<unsigned int>::size_type queue_end = 0;
1559 std::vector <bool> isOriented (bg->
Triangles_m.size(),
false);
1566 while (isOriented [triangle_id])
1570 if (!hasInwardPointingNormal (bg, triangle_id)) {
1571 std::swap (bg->PointID (triangle_id, 2), bg->PointID (triangle_id, 3));
1573 isOriented [triangle_id] =
true;
1576 triangles [queue_end++] = triangle_id;
1578 for (
auto neighbor_id: neighbors [triangle_id]) {
1579 if (isOriented [neighbor_id])
continue;
1580 orientTriangle (bg, triangle_id, neighbor_id);
1581 isOriented [neighbor_id] =
true;
1582 triangles[queue_end++] = neighbor_id;
1584 triangle_id = triangles [++queue_cursor];
1585 }
while (queue_cursor < queue_end);
1589 *gmsg <<
"* " << __func__ <<
": mesh is contiguous." <<
endl;
1591 *gmsg <<
"* " << __func__ <<
": mesh is discontiguous (" << parts <<
") parts." <<
endl;
1593 *gmsg <<
"* Triangle Normal built done." <<
endl;
1599 *gmsg <<
"* Initializing Boundary Geometry..." <<
endl;
1605 *gmsg <<
"* Found additional aperture." <<
endl;
1606 for (
unsigned int i=0; i<
apert_m.size(); i=i+3)
1607 *gmsg <<
"* zmin = " <<
apert_m[i]
1609 <<
" r= " <<
apert_m[i+2] << endl;
1620 *gmsg <<
"* X-scale all points of geometry by " << xscale <<
endl;
1621 *gmsg <<
"* Y-scale all points of geometry by " << yscale <<
endl;
1622 *gmsg <<
"* Z-scale all points of geometry by " << zscale <<
endl;
1623 *gmsg <<
"* Scale all points of geometry by " << xyzscale <<
endl;
1626 #if defined (NDEBUG)
1629 rc = H5SetErrorHandler (H5AbortErrorhandler);
1630 assert (rc != H5_ERR);
1631 H5SetVerbosityLevel (1);
1633 h5_prop_t props = H5CreateFileProp ();
1635 H5SetPropFileMPIOCollective (props, &comm);
1636 h5_file_t f = H5OpenFile (
h5FileName_m.c_str(), H5_O_RDONLY, props);
1637 H5CloseProp (props);
1639 h5t_mesh_t* m = NULL;
1640 H5FedOpenTriangleMesh (f,
"0", &m);
1641 H5FedSetLevel (m, 0);
1643 auto numTriangles = H5FedGetNumElementsTotal (m);
1647 h5_loc_id_t local_id;
1649 h5t_iterator_t* iter = H5FedBeginTraverseEntities (m, 0);
1650 while ((local_id = H5FedTraverseEntities (iter)) >= 0) {
1651 h5_loc_id_t local_vids[4];
1652 H5FedGetVertexIndicesOfEntity (m, local_id, local_vids);
1654 PointID (i, 1) = local_vids[0];
1655 PointID (i, 2) = local_vids[1];
1656 PointID (i, 3) = local_vids[2];
1659 H5FedEndTraverseEntities (iter);
1662 int num_points = H5FedGetNumVerticesTotal (m);
1664 for (i = 0; i < num_points; i++) {
1666 H5FedGetVertexCoordsByIndex (m, i, P);
1668 P[0] * xyzscale * xscale,
1669 P[1] * xyzscale * yscale,
1670 P[2] * xyzscale * zscale + zshift));
1674 *gmsg <<
"* Reading mesh done." <<
endl;
1676 Local::computeGeometryInterval (
this);
1686 Local::makeTriangleNormalInwardPointing (
this);
1702 *gmsg <<
"* Triangle barycent built done." <<
endl;
1704 *gmsg << *
this <<
endl;
1731 *gmsg <<
"* " << __func__ <<
": "
1738 const Ray r =
Ray (P, v_);
1742 MIN2 (P[2], Q[2]) };
1746 MAX2 (P[2], Q[2]) };
1769 std::unordered_set<int> triangle_ids;
1770 for (
int i = i_min; i <= i_max; i++) {
1771 for (
int j = j_min; j <= j_max; j++) {
1772 for (
int k = k_min; k <= k_max; k++) {
1776 if (
debugFlags_m & debug_intersectTinyLineSegmentBoundary) {
1777 *gmsg <<
"* " << __func__ <<
": "
1778 <<
" Test voxel: (" << i <<
", " << j <<
", " << k <<
"), "
1795 const auto triangles_intersecting_voxel =
1797 if (triangles_intersecting_voxel !=
voxelMesh_m.ids.end()) {
1798 triangle_ids.insert (
1799 triangles_intersecting_voxel->second.begin(),
1800 triangles_intersecting_voxel->second.end());
1809 int num_intersections = 0;
1810 int tmp_intersect_result = 0;
1812 for (
auto it = triangle_ids.begin ();
1813 it != triangle_ids.end ();
1822 if (
debugFlags_m & debug_intersectTinyLineSegmentBoundary) {
1823 *gmsg <<
"* " << __func__ <<
": "
1824 <<
" Test triangle: " << *it
1825 <<
" intersect: " << tmp_intersect_result
1832 switch (tmp_intersect_result) {
1840 if (gsl_fcmp (Q[0] - P[0], 0.0,
EPS) != 0) {
1841 t = (tmp_intersect_pt[0] - P[0]) / (Q[0] - P[0]);
1842 }
else if (gsl_fcmp (Q[1] - P[1], 0.0,
EPS) != 0) {
1843 t = (tmp_intersect_pt[1] - P[1]) / (Q[1] - P[1]);
1845 t = (tmp_intersect_pt[2] - P[2]) / (Q[2] - P[2]);
1847 num_intersections++;
1850 if (
debugFlags_m & debug_intersectTinyLineSegmentBoundary) {
1851 *gmsg <<
"* " << __func__ <<
": "
1857 intersect_pt = tmp_intersect_pt;
1858 triangle_id = (*it);
1862 assert (tmp_intersect_result != -1);
1866 return num_intersections;
1884 *gmsg <<
"* " << __func__ <<
": "
1894 int intersect_result = 0;
1896 int i_min, j_min, k_min;
1897 int i_max, j_max, k_max;
1904 MIN2 (P0[2], Q[2]) };
1908 MAX2 (P0[2], Q[2]) };
1911 }
while (( (i_max-i_min+1) * (j_max-j_min+1) * (k_max-k_min+1)) > 27);
1916 for (
int l = 1; l <=
n; l++, P = Q) {
1919 if (triangle_id != -1) {
1924 if (
debugFlags_m & debug_intersectLineSegmentBoundary) {
1925 *gmsg <<
"* " << __func__ <<
": "
1926 <<
" result=" << intersect_result
1927 <<
" intersection pt: " << intersect_pt
1932 return intersect_result;
1956 *gmsg <<
"* " << __func__ <<
": "
1977 int tmp_triangle_id = -1;
1979 if (tmp_triangle_id >= 0) {
1980 intersect_pt = tmp_intersect_pt;
1981 triangle_id = tmp_triangle_id;
1984 else if (Parttype == 1)
1992 *gmsg <<
"* " << __func__ <<
":"
1993 <<
" result=" << ret;
1995 *gmsg <<
" intersetion=" << intersect_pt;
2008 of.open (fn.c_str ());
2009 assert (of.is_open ());
2011 of <<
"# vtk DataFile Version 2.0" <<
std::endl;
2012 of <<
"generated using DataSink::writeGeoToVtk" <<
std::endl;
2013 of <<
"ASCII" << std::endl <<
std::endl;
2014 of <<
"DATASET UNSTRUCTURED_GRID" <<
std::endl;
2016 for (
unsigned int i = 0; i <
Points_m.size (); i++)
2032 of <<
"5" << std::endl;
2034 of <<
"SCALARS " <<
"cell_attribute_data" <<
" float " <<
"1" <<
std::endl;
2035 of <<
"LOOKUP_TABLE " <<
"default" <<
std::endl;
2037 of << (
float)(i) << std::endl;
2044 os <<
"* ************* B O U N D A R Y G E O M E T R Y *********************************** " <<
endl;
2052 if (
getTopology () == std::string (
"BOXCORNER")) {
2057 os <<
"* Total triangle num " <<
Triangles_m.size() <<
'\n'
2058 <<
"* Total points num " <<
Points_m.size () <<
'\n'
2059 <<
"* Geometry bounds(m) Max= " <<
maxExtent_m <<
'\n'
2062 <<
"* Resolution of voxel mesh " <<
voxelMesh_m.nr_m <<
'\n'
2063 <<
"* Size of voxel " <<
voxelMesh_m.sizeOfVoxel <<
'\n'
2064 <<
"* Number of voxels in mesh " <<
voxelMesh_m.ids.size () <<
'\n'
2066 os <<
"* ********************************************************************************** " <<
endl;
2110 const int& triId = itsBunch->
TriID[i];
2111 const double& incQ = itsBunch->
Q[i];
2112 const Vector_t& incMomentum = itsBunch->
P[i];
2113 const double p_sq =
dot (incMomentum, incMomentum);
2123 }
else if (BGtag & BGphysics::SecondaryEmission) {
2128 ERRORMSG (
" cosTheta = " << cosTheta <<
" < 0 (!)" <<
endl <<
2129 " particle position = " << itsBunch->
R[i] <<
endl <<
2130 " incident momentum=" << incMomentum <<
endl <<
2133 " intecoords = " << intecoords <<
endl <<
2134 " triangle ID = " << triId <<
endl <<
2135 " triangle = (" <<
getPoint(triId, 1)
2138 assert(cosTheta>=0);
2141 if (intecoords !=
Point (triId, 1)) {
2174 const int& triId = itsBunch->
TriID[i];
2175 const double& incQ = itsBunch->
Q[i];
2176 const Vector_t& incMomentum = itsBunch->
P[i];
2177 const double p_sq =
dot (incMomentum, incMomentum);
2187 }
else if (BGtag & BGphysics::SecondaryEmission) {
2193 ERRORMSG (
" cosTheta = " << cosTheta <<
" < 0 (!)" <<
endl <<
2194 " particle position = " << itsBunch->
R[i] <<
endl <<
2195 " incident momentum=" << incMomentum <<
endl <<
2198 " intecoords = " << intecoords <<
endl <<
2199 " triangle ID = " << triId <<
endl <<
2201 " Particle No. = (" << i <<
")"
2203 assert(cosTheta>=0);
2206 if (intecoords !=
Point (triId, 1)) {
2249 for (
size_t i = 0; i <
n; i++) {
2267 centroid, itsBunch->
getdT (), E,
B);
2283 for (
size_t i = 0; i < nData; i++) {
2309 double x_low =
minExtent_m (0) + 0.5 * len [0] - 0.49 * len [0];
2315 double x_up =
minExtent_m (0) + 0.5 * len [0] + 0.49 * len [0];
2321 double y_low =
minExtent_m (1) + 0.5 * len [1] - 0.49 * len [1];
2327 double y_up =
minExtent_m (1) + 0.5 * len [1] + 0.49 * len [1];
2329 for (
size_t i = 0; i < n / 2; i++) {
2333 while (zCoord > 0.000001 ||
2334 zCoord < - 0.000001 ||
2355 for (
size_t i = 0; i < n / 2; i++) {
2396 for (
size_t i = 0; i < nData; i++) {
2407 for (
size_t i = 0; i <
n; i++) {
2439 for (
size_t i = 0; i < nData; i++) {
ParticleAttrib< Vector_t > P
std::vector< double > TriFEPartloss_m
Interface for basic beam line object.
void define(Object *newObject)
Define a new object.
The base class for all OPAL definitions.
std::vector< Vector_t > partsr_m
#define PointID(triangle_id, vertex_id)
Triangle(const Vector_t &v1, const Vector_t &v2, const Vector_t &v3)
int intersectLineSegmentBoundary(const Vector_t &P0, const Vector_t &P1, Vector_t &intersection_pt, int &triangle_id)
const Vector_t & v2() const
IpplTimings::TimerRef TfastIsInside_m
void writeGeomToVtk(std::string fn)
PETE_TUTree< FnFabs, typename T::PETE_Expr_t > fabs(const PETE_Expr< T > &l)
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
The base class for all OPAL exceptions.
int intersectLineTriangle(const enum INTERSECTION_TESTS kind, const Vector_t &P0, const Vector_t &P1, const int triangle_id, Vector_t &I)
ParticleAttrib< double > Q
IpplTimings::TimerRef TisInside_m
PETE_TUTree< FnCeil, typename T::PETE_Expr_t > ceil(const PETE_Expr< T > &l)
void scale(const Vector_t &scale)
std::vector< double > apert_m
SecondaryEmissionPhysics sec_phys_m
void createParticlesOnSurface(size_t n, double darkinward, OpalBeamline &itsOpalBeamline, PartBunchBase< double, 3 > *itsBunch)
std::vector< Attribute > itsAttr
The object attributes (see Attribute.hh).
std::vector< Vector_t > Points_m
bool isInside(const Vector_t &P) const
RandomNumberGen IpplRandom
unsigned long getFieldAt(const unsigned int &, const Vector_t &, const long &, const double &, Vector_t &, Vector_t &)
virtual void execute()
Execute the command.
std::vector< Vector_t > TriNormals_m
double dot(const Vector3D &lhs, const Vector3D &rhs)
Vector dot product.
std::vector< double > TriAreas_m
std::vector< std::array< unsigned int, 4 > > Triangles_m
constexpr double m_e
The electron rest mass in GeV.
constexpr double alpha
The fine structure constant, no dimension.
const Vector_t & v1() const
std::vector< short > TriBGphysicstag_m
Attribute makeStringArray(const std::string &name, const std::string &help)
Create a string array attribute.
static OpalData * getInstance()
int emitSecondaryFurmanPivi(const Vector_t &intecoords, const int i, PartBunchBase< double, 3 > *itsBunch, double &seyNum)
const std::string & getOpalName() const
Return object name.
struct BoundaryGeometry::@91 voxelMesh_m
int intersectTinyLineSegmentBoundary(const Vector_t &, const Vector_t &, Vector_t &, int &)
Abstract base class for accelerator geometry classes.
void updateElement(ElementBase *element)
void createPriPart(size_t n, double darkinward, OpalBeamline &itsOpalBeamline, PartBunchBase< double, 3 > *itsBunch)
Voxel(const Vector_t &min, const Vector_t &max)
int emitSecondaryNone(const Vector_t &intecoords, const int &triId)
IpplTimings::TimerRef Tinitialize_m
int partInside(const Vector_t &r, const Vector_t &v, const double dt, int Parttype, const double Qloss, Vector_t &intecoords, int &triId)
ParticleAttrib< int > TriID
Vector_t mapIndices2Voxel(const int, const int, const int)
constexpr double c
The velocity of light in m/s.
Vector3D cross(const Vector3D &lhs, const Vector3D &rhs)
Vector cross product.
void registerOwnership(const AttributeHandler::OwnerType &itsClass) const
static BoundaryGeometry * find(const std::string &name)
static void startTimer(TimerRef t)
std::vector< double > TriSePartloss_m
IpplTimings::TimerRef TRayTrace_m
std::vector< double > getRealArray(const Attribute &attr)
Get array value.
std::vector< double > TriPrPartloss_m
IpplTimings::TimerRef TPartInside_m
Vektor< double, 3 > Vector_t
int fastIsInside(const Vector_t &reference_pt, const Vector_t &P)
static MPI_Comm getComm()
int intersect(const Triangle &t) const
std::vector< Vector_t > partsp_m
Tps< T > sqrt(const Tps< T > &x)
Square root.
virtual ~BoundaryGeometry()
int emitSecondaryVaughan(const Vector_t &intecoords, const int i, PartBunchBase< double, 3 > *itsBunch, double &seyNum)
bool intersect(const Ray &r, double &tmin, double &tmax) const
bool intersect(const Ray &r) const
Vector_t mapPoint2Voxel(const Vector_t &)
Object * find(const std::string &name)
Find entry.
The base class for all OPAL objects.
const Vector_t & v3() const
void scale(const Vector_t &scaleby, const Vector_t &shiftby)
void setOpalName(const std::string &name)
Set object name.
void getMessage(Message &m, T &t)
int intersectTriangleVoxel(const int triangle_id, const int i, const int j, const int k)
#define mapPoint2VoxelIndices(pt, i, j, k)
Inform & printInfo(Inform &os) const
Attribute makeRealArray(const std::string &name, const std::string &help)
Create real array attribute.
Ray(Vector_t o, Vector_t d)
bool builtin
Built-in flag.
void computeMeshVoxelization(void)
std::vector< Vector_t > TriBarycenters_m
std::string::iterator iterator
void putMessage(Message &m, const T &t)
static TimerRef getTimer(const char *nm)
double getReal(const Attribute &attr)
Return real value.
static void stopTimer(TimerRef t)
bool ppdebug
ppdebug flag.
Message * receive_block(int &node, int &tag)
Attribute makeString(const std::string &name, const std::string &help)
Make string attribute.
virtual void update()
Update this object.
static Communicate * Comm
int mapVoxelIndices2ID(const int i, const int j, const int k)
Attribute makeReal(const std::string &name, const std::string &help)
Make real attribute.
virtual bool canReplaceBy(Object *object)
Test if replacement is allowed.
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
virtual int broadcast_all(Message *, int)
void nSec(const double &incEnergy, const double &cosTheta, const int &matNumber, int &seNum, int &seType, const double &incQ, const Vector_t &TriNorm, const Vector_t &inteCoords, const Vector_t &localX, PartBunchBase< double, 3 > *itsBunch, double &seyNum, const double &ppVw, const double &vVThermal, const bool nEmissionMode)
virtual BoundaryGeometry * clone(const std::string &name)
Return a clone.
const Ray & operator=(const Ray &a)=delete
const Vector_t & getPoint(const int triangle_id, const int vertex_id)
std::string getTopology() const
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
int intersectRayBoundary(const Vector_t &P, const Vector_t &v, Vector_t &I)
Inform & endl(Inform &inf)
std::string getString(const Attribute &attr)
Get string value.