Line data Source code
1 : //* This file is part of SALAMANDER: Software for Advanced Large-scale Analysis of MAgnetic confinement for Numerical Design, Engineering & Research,
2 : //* A multiphysics application for modeling plasma facing components
3 : //* https://github.com/idaholab/salamander
4 : //* https://mooseframework.inl.gov/salamander
5 : //*
6 : //* SALAMANDER is powered by the MOOSE Framework
7 : //* https://www.mooseframework.inl.gov
8 : //*
9 : //* Licensed under LGPL 2.1, please see LICENSE for details
10 : //* https://www.gnu.org/licenses/lgpl-2.1.html
11 : //*
12 : //* Copyright 2025, Battelle Energy Alliance, LLC
13 : //* ALL RIGHTS RESERVED
14 : //*
15 :
16 : #include "ElementSampler.h"
17 : #include "MooseRandom.h"
18 : #include "libmesh/enum_to_string.h"
19 : #include "FEProblemBase.h"
20 :
21 : namespace SALAMANDER
22 : {
23 :
24 176 : ElementSampler::ElementSampler(FEProblemBase & problem,
25 : const unsigned int seed,
26 176 : MooseRandom & generator)
27 176 : : _seed(seed),
28 176 : _generator(generator),
29 176 : _arbitrary_qrule(problem.mesh().dimension(), FIRST),
30 176 : _fe(FEBase::build(problem.mesh().dimension(), FEType(CONSTANT, MONOMIAL)))
31 : {
32 176 : _fe->attach_quadrature_rule(&_arbitrary_qrule);
33 : _fe->get_xyz();
34 176 : }
35 :
36 : std::vector<Point>
37 32046 : ElementSampler::sampleElement(const Elem * elem, const unsigned int samples)
38 : {
39 : mooseAssert(samples != 0, "You must request a non-zero number of samples");
40 : // so long as p is a point sampled from a cube with dimensions
41 : // [-1, 1] x [-1, 1] x [-1, 1]
42 : // this lambda will put the point inside of the libmesh reference pyramid
43 464400 : auto put_particle_in_pyramid = [](Point & p)
44 : {
45 : // if the point is not in the pyramids along the x axis
46 : // we need to find which one it is in
47 464400 : if (!((p(0) - std::abs(p(2)) < 0.0) && (-p(0) - std::abs(p(2)) < 0.0) &&
48 231108 : (p(1) - std::abs(p(2)) < 0.0) && (-p(1) - std::abs(p(2)) < 0.0)))
49 : {
50 : // checking out the pyramids with the
51 : // the x axis going through the center
52 310302 : if (((p(0) - p(1) > 0.0) && (-p(0) - p(1) > 0.0)) ||
53 233910 : ((p(0) + p(1) > 0.0) && (-p(0) + p(1) > 0.0)))
54 : {
55 : auto temporary = p(1);
56 154584 : p(1) = p(2);
57 154584 : p(2) = temporary;
58 : }
59 : // all the rest of the points are pyramids with the
60 : // the y axis going through the center
61 : else
62 : {
63 : auto temporary = p(0);
64 155718 : p(0) = p(1);
65 155718 : p(1) = p(2);
66 155718 : p(2) = temporary;
67 : }
68 : }
69 : // now we have all of the pyramids oriented
70 : // along the x axis but we need to put them in
71 : // the correct orientation to match the reference element
72 464400 : if (p(2) > 0)
73 232416 : p(2) = 1 - p(2);
74 : else
75 231984 : p(2) = 1 + p(2);
76 464400 : };
77 :
78 32046 : std::vector<Point> reference_points = std::vector<Point>(samples);
79 :
80 : // reseeding with the element id with an additional
81 : // user selected seed for consistency
82 : // note that this is only consistent across process counts when element ids are
83 : // also consistent across processor counts which in general is not the case
84 32046 : _generator.seed(elem->id() + _seed);
85 :
86 32046 : switch (elem->type())
87 : {
88 : // 1D reference elements x = [-1, 1]
89 : case EDGE2:
90 : {
91 1552 : for (const auto i : make_range(samples))
92 1440 : reference_points[i](0) = 2.0 * _generator.rand() - 1.0;
93 : break;
94 : }
95 : // 2D trianglular element where the vertices are at
96 : // (0,0) (0,1), (0, 1)
97 : case TRI3:
98 : {
99 51200 : for (const auto i : make_range(samples))
100 : {
101 : // sample on a square x = [0, 1] and y = [0, 1]
102 48000 : reference_points[i](0) = _generator.rand();
103 48000 : reference_points[i](1) = _generator.rand();
104 : // if our points are not in the triangle we mirror them into the triangle
105 48000 : if (reference_points[i](1) > 1 - reference_points[i](0))
106 : {
107 23920 : reference_points[i](0) = 1 - reference_points[i](0);
108 23920 : reference_points[i](1) = 1 - reference_points[i](1);
109 : }
110 : }
111 : break;
112 : }
113 : // 2D square reference element where x = [-1, 1] and y = [-1, 1]
114 : case QUAD4:
115 : {
116 25600 : for (const auto i : make_range(samples))
117 : {
118 24000 : reference_points[i](0) = 2.0 * _generator.rand() - 1.0;
119 24000 : reference_points[i](1) = 2.0 * _generator.rand() - 1.0;
120 : }
121 : break;
122 : }
123 : // 3D cubic basis element where x = [-1, 1] and y = [-1, 1] and z = [-1, 1]
124 : case HEX8:
125 : {
126 21704 : for (const auto i : make_range(samples))
127 : {
128 20640 : reference_points[i](0) = 2.0 * _generator.rand() - 1.0;
129 20640 : reference_points[i](1) = 2.0 * _generator.rand() - 1.0;
130 20640 : reference_points[i](2) = 2.0 * _generator.rand() - 1.0;
131 : }
132 : break;
133 : }
134 : case PYRAMID5:
135 : {
136 97668 : for (const auto i : make_range(samples))
137 : {
138 92880 : reference_points[i](0) = 2.0 * _generator.rand() - 1.0;
139 92880 : reference_points[i](1) = 2.0 * _generator.rand() - 1.0;
140 92880 : reference_points[i](2) = 2.0 * _generator.rand() - 1.0;
141 :
142 92880 : put_particle_in_pyramid(reference_points[i]);
143 : }
144 : break;
145 : }
146 : // 3D elements with all triangular faces with nodes at
147 : // (0,0,0), (1,0,0), (0,1,0), (0,0,1)
148 : case TET4:
149 : {
150 : Real temporary;
151 390672 : for (const auto i : make_range(samples))
152 : {
153 371520 : reference_points[i](0) = 2.0 * _generator.rand() - 1.0;
154 371520 : reference_points[i](1) = 2.0 * _generator.rand() - 1.0;
155 371520 : reference_points[i](2) = 2.0 * _generator.rand() - 1.0;
156 :
157 : // start by folding the cube into a pyramid
158 371520 : put_particle_in_pyramid(reference_points[i]);
159 :
160 : // now we are going to fold the pyramid into a single tet
161 371520 : if (reference_points[i](0) < 0.0)
162 185778 : reference_points[i](0) *= -1;
163 :
164 371520 : if (reference_points[i](1) < 0.0)
165 186024 : reference_points[i](1) *= -1;
166 :
167 : // at this point we have two tetrahedrons and we need to Translate them over the proper
168 : // tet
169 371520 : if (reference_points[i](1) > reference_points[i](0))
170 : {
171 : temporary = reference_points[i](0);
172 186054 : reference_points[i](0) = reference_points[i](1);
173 186054 : reference_points[i](1) = temporary;
174 : }
175 :
176 : // now all of our points are in a tet that is bounded by
177 : // (0,0,0), (1,0,0), (1,1,0), (0,0,1)
178 : // to make this our reference tet we perform an affine transformation
179 371520 : reference_points[i](0) -= reference_points[i](1);
180 : }
181 : break;
182 : }
183 : case PRISM6:
184 : {
185 43408 : for (const auto i : make_range(samples))
186 : {
187 : // sample on a rectangular prism x = [0, 1] and y = [0, 1] z = [-1,1]
188 41280 : reference_points[i](0) = _generator.rand();
189 41280 : reference_points[i](1) = _generator.rand();
190 41280 : reference_points[i](2) = 2 * _generator.rand() - 1;
191 : // if our points are not in the triangle we mirror them into the prism
192 41280 : if (reference_points[i](1) > 1 - reference_points[i](0))
193 : {
194 20224 : reference_points[i](0) = 1 - reference_points[i](0);
195 20224 : reference_points[i](1) = 1 - reference_points[i](1);
196 : }
197 : }
198 : break;
199 : }
200 2 : default:
201 4 : mooseError("Particle Initialization has not been implemented for elements of type " +
202 2 : Utility::enum_to_string(elem->type()) + ".\n" +
203 0 : "If your problem requires this element type please reach out to us at\n" +
204 : "https://github.com/idaholab/salamander/discussions");
205 : }
206 :
207 32044 : _arbitrary_qrule.setPoints(reference_points);
208 32044 : _fe->reinit(elem);
209 : // now that all of the particle locations have been placed we need to
210 : // set up the data they will need to be made into actual rays
211 32044 : const auto physical_points = _fe->get_xyz();
212 :
213 32044 : return physical_points;
214 : }
215 : }
|
