Returns the integral of the function over a specified domain.
135{
136
137
138 const T *x1_ptr, *x2_ptr;
139 bool switch_bounds;
142 if (MooseUtils::absoluteFuzzyEqual(rawA, rawB))
143 return 0.0;
144 else if (rawB > rawA)
145 {
146 x1_ptr = &xA;
147 x2_ptr = &xB;
148 raw1 = rawA;
149 raw2 = rawB;
150 switch_bounds = false;
151 }
152 else
153 {
154 x1_ptr = &xB;
155 x2_ptr = &xA;
156 raw1 = rawB;
157 raw2 = rawA;
158 switch_bounds = true;
159 }
160 const auto & x1 = *x1_ptr;
161 const auto & x2 = *x2_ptr;
162
163
164 T integral = 0.0;
165
166 const auto n =
_x.size();
167 const unsigned int i1 =
168 raw1 <=
_x[n - 1] ? std::distance(
_x.begin(), std::upper_bound(
_x.begin(),
_x.end(), raw1))
169 : n;
170 const unsigned int i2 =
171 raw2 <=
_x[n - 1] ? std::distance(
_x.begin(), std::upper_bound(
_x.begin(),
_x.end(), raw2))
172 : n;
173 unsigned int i = i1;
174 while (i <= i2)
175 {
176 if (i == 0)
177 {
178
179 T integral1, integral2;
181 {
182 const auto dydx = (
_y[1] -
_y[0]) / (
_x[1] -
_x[0]);
183 const auto y1 =
_y[0] + dydx * (x1 -
_x[0]);
184 integral1 = 0.5 * (y1 +
_y[0]) * (
_x[0] - x1);
185 if (i2 == i)
186 {
187 const auto y2 =
_y[0] + dydx * (x2 -
_x[0]);
188 integral2 = 0.5 * (y2 +
_y[0]) * (
_x[0] - x2);
189 }
190 else
191 integral2 = 0.0;
192 }
193 else
194 {
195 integral1 =
_y[0] * (
_x[0] - x1);
196 if (i2 == i)
197 integral2 =
_y[0] * (
_x[0] - x2);
198 else
199 integral2 = 0.0;
200 }
201
202 integral += integral1 - integral2;
203 }
204 else if (i == n)
205 {
206
207 T integral1, integral2;
209 {
210 const auto dydx = (
_y[n - 1] -
_y[n - 2]) / (
_x[n - 1] -
_x[n - 2]);
211 const auto y2 =
_y[n - 1] + dydx * (x2 -
_x[n - 1]);
212 integral2 = 0.5 * (y2 +
_y[n - 1]) * (x2 -
_x[n - 1]);
213 if (i1 == n)
214 {
215 const auto y1 =
_y[n - 1] + dydx * (x1 -
_x[n - 1]);
216 integral1 = 0.5 * (y1 +
_y[n - 1]) * (x1 -
_x[n - 1]);
217 }
218 else
219 integral1 = 0.0;
220 }
221 else
222 {
223 integral2 =
_y[n - 1] * (x2 -
_x[n - 1]);
224 if (i1 == n)
225 integral1 =
_y[n - 1] * (x1 -
_x[n - 1]);
226 else
227 integral1 = 0.0;
228 }
229
230 integral += integral2 - integral1;
231 }
232 else
233 {
234 T integral1;
235 if (i == i1)
236 {
237 const auto dydx = (
_y[i] -
_y[i - 1]) / (
_x[i] -
_x[i - 1]);
238 const auto y1 =
_y[i - 1] + dydx * (x1 -
_x[i - 1]);
239 integral1 = 0.5 * (y1 +
_y[i - 1]) * (x1 -
_x[i - 1]);
240 }
241 else
242 integral1 = 0.0;
243
244 T integral2;
245 if (i == i2)
246 {
247 const auto dydx = (
_y[i] -
_y[i - 1]) / (
_x[i] -
_x[i - 1]);
248 const auto y2 =
_y[i - 1] + dydx * (x2 -
_x[i - 1]);
249 integral2 = 0.5 * (y2 +
_y[i - 1]) * (x2 -
_x[i - 1]);
250 }
251 else
252 integral2 = 0.5 * (
_y[i] +
_y[i - 1]) * (
_x[i] -
_x[i - 1]);
253
254 integral += integral2 - integral1;
255 }
256
257 i++;
258 }
259
260
261 if (switch_bounds)
262 return -1.0 * integral;
263 else
264 return integral;
265}