https://mooseframework.inl.gov
LayeredBase.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "LayeredBase.h"
11 
12 // MOOSE includes
13 #include "MooseEnum.h"
14 #include "MooseMesh.h"
15 #include "SubProblem.h"
16 #include "UserObject.h"
17 
18 #include "libmesh/mesh_tools.h"
19 #include "libmesh/point.h"
20 
23 {
25  MooseEnum directions("x y z");
26 
27  params.addRequiredParam<MooseEnum>("direction", directions, "The direction of the layers.");
28  params.addParam<unsigned int>("num_layers", "The number of layers.");
29  params.addParam<std::vector<Real>>("bounds",
30  "The 'bounding' positions of the layers i.e.: '0, "
31  "1.2, 3.7, 4.2' will mean 3 layers between those "
32  "positions.");
33  params.addRangeCheckedParam<unsigned int>("bound_uniform_splits",
34  "bound_uniform_splits > 0",
35  "The number of times the bins specified in 'bounds' "
36  "should be split uniformly.");
37 
38  MooseEnum sample_options("direct interpolate average", "direct");
39  params.addParam<MooseEnum>("sample_type",
40  sample_options,
41  "How to sample the layers. 'direct' means get the value of the layer "
42  "the point falls in directly (or average if that layer has no value). "
43  " 'interpolate' does a linear interpolation between the two closest "
44  "layers. 'average' averages the two closest layers.");
45 
46  params.addParam<unsigned int>("average_radius",
47  1,
48  "When using 'average' sampling this is how "
49  "the number of values both above and below "
50  "the layer that will be averaged.");
51 
52  params.addParam<bool>(
53  "cumulative",
54  false,
55  "When true the value in each layer is the sum of the values up to and including that layer");
56  params.addParam<bool>(
57  "positive_cumulative_direction",
58  true,
59  "When 'cumulative' is true, whether the direction for summing the cumulative value "
60  "is the positive direction or negative direction");
61 
62  params.addParam<std::vector<SubdomainName>>(
63  "block", "The list of block ids (SubdomainID) that this object will be applied");
64 
65  params.addParam<std::vector<SubdomainName>>("layer_bounding_block",
66  "List of block ids (SubdomainID) that are used to "
67  "determine the upper and lower geometric bounds for "
68  "all layers. If this is not specified, the ids "
69  "specified in 'block' are used for this purpose.");
70 
71  params.addParam<Real>("direction_min",
72  "Minimum coordinate along 'direction' that bounds the layers");
73  params.addParam<Real>("direction_max",
74  "Maximum coordinate along 'direction' that bounds the layers");
75  params.addParamNamesToGroup(
76  "direction num_layers bounds direction_min direction_max bound_uniform_splits",
77  "Layers extent and definition");
78  params.addParamNamesToGroup("sample_type average_radius cumulative positive_cumulative_direction",
79  "Value sampling / aggregating");
80  return params;
81 }
82 
84  : Restartable(parameters.getCheckedPointerParam<SubProblem *>("_subproblem")->getMooseApp(),
85  parameters.get<std::string>("_object_name") + "_layered_base",
86  "LayeredBase",
87  parameters.get<THREAD_ID>("_tid")),
88  _layered_base_name(parameters.get<std::string>("_object_name")),
89  _layered_base_params(parameters),
90  _direction_enum(parameters.get<MooseEnum>("direction")),
91  _direction(_direction_enum),
92  _sample_type(parameters.get<MooseEnum>("sample_type")),
93  _average_radius(parameters.get<unsigned int>("average_radius")),
94  _using_displaced_mesh(_layered_base_params.get<bool>("use_displaced_mesh")),
95  _layer_values(declareRestartableData<std::vector<Real>>("layer_values")),
96  _layer_has_value(declareRestartableData<std::vector<int>>("layer_has_value")),
97  _cumulative(parameters.get<bool>("cumulative")),
98  _positive_cumulative_direction(parameters.get<bool>("positive_cumulative_direction")),
99  _layered_base_subproblem(*parameters.getCheckedPointerParam<SubProblem *>("_subproblem")),
100  _layer_bounding_blocks(),
101  _has_direction_max_min(false)
102 {
103  if (_layered_base_params.isParamValid("num_layers") &&
105  mooseError("'bounds' and 'num_layers' cannot both be set");
106 
107  if (!_cumulative && parameters.isParamSetByUser("positive_cumulative_direction"))
108  mooseWarning(
109  "The 'positive_cumulative_direction' parameter is unused when 'cumulative' is false");
110 
111  if (_layered_base_params.isParamValid("num_layers"))
112  {
113  _num_layers = _layered_base_params.get<unsigned int>("num_layers");
114  _interval_based = true;
115  }
116  else if (_layered_base_params.isParamValid("bounds"))
117  {
118  _interval_based = false;
119 
120  _layer_bounds = _layered_base_params.get<std::vector<Real>>("bounds");
121 
122  if (_layer_bounds.size() < 2)
123  mooseError("At least two boundaries must be provided in 'bounds' to form layers!");
124 
125  // Make sure the bounds are sorted - we're going to depend on this
126  std::sort(_layer_bounds.begin(), _layer_bounds.end());
127 
128  // If requested, we uniformly split the layers.
129  if (_layered_base_params.isParamValid("bound_uniform_splits"))
130  {
131  const auto splits = _layered_base_params.get<unsigned int>("bound_uniform_splits");
132  for (unsigned int s = 0; s < splits; ++s)
133  {
134  std::vector<Real> new_bnds;
135  new_bnds.reserve(2 * (_layer_bounds.size() - 1) + 1);
136  for (unsigned int i = 0; i < _layer_bounds.size() - 1; ++i)
137  {
138  new_bnds.emplace_back(_layer_bounds[i]);
139  new_bnds.emplace_back(0.5 * (_layer_bounds[i] + _layer_bounds[i + 1]));
140  }
141  _layer_bounds = new_bnds;
142  }
143  }
144 
145  _num_layers = _layer_bounds.size() - 1; // Layers are only in-between the bounds
146  _direction_min = _layer_bounds.front();
147  _direction_max = _layer_bounds.back();
148  _has_direction_max_min = true;
149  }
150  else
151  mooseError("One of 'bounds' or 'num_layers' must be specified");
152 
153  if (!_interval_based && _sample_type == 1)
154  mooseError("'sample_type = interpolate' not supported with 'bounds'");
155 
156  bool has_layer_bounding_block = _layered_base_params.isParamValid("layer_bounding_block");
157  bool has_block = _layered_base_params.isParamValid("block");
158  bool has_direction_min = _layered_base_params.isParamValid("direction_min");
159  bool has_direction_max = _layered_base_params.isParamValid("direction_max");
160 
161  if (_has_direction_max_min && has_direction_min)
162  mooseWarning("'direction_min' is unused when providing 'bounds'");
163 
164  if (_has_direction_max_min && has_direction_max)
165  mooseWarning("'direction_max' is unused when providing 'bounds'");
166 
167  // can only specify one of layer_bounding_block or the pair direction_max/min
168  if (has_layer_bounding_block && (has_direction_min || has_direction_max))
169  mooseError("Only one of 'layer_bounding_block' and the pair 'direction_max' and "
170  "'direction_min' can be provided");
171 
172  // if either one of direction_min or direction_max is specified, must provide the other one
173  if (has_direction_min != has_direction_max)
174  mooseError("If providing the layer max/min directions, both 'direction_max' and "
175  "'direction_min' must be specified.");
176 
177  if (has_layer_bounding_block)
179  _layered_base_params.get<std::vector<SubdomainName>>("layer_bounding_block"));
180  else if (has_block)
182  _layered_base_params.get<std::vector<SubdomainName>>("block"));
183 
184  // specifying the direction max/min overrides anything set with the 'block'
185  if (has_direction_min && has_direction_max)
186  {
187  _direction_min = parameters.get<Real>("direction_min");
188  _direction_max = parameters.get<Real>("direction_max");
189  _has_direction_max_min = true;
190 
192  mooseError("'direction_max' must be larger than 'direction_min'");
193  }
194 
195  _layer_values.resize(_num_layers);
197 
198  // if we haven't already figured out the max/min in specified direction
199  // (either with the 'bounds' or explicit specification from the user), do so
201  getBounds();
202 
204 }
205 
206 Real
207 LayeredBase::integralValue(const Point & p) const
208 {
209  unsigned int layer = getLayer(p);
210 
211  int higher_layer = -1;
212  int lower_layer = -1;
213 
214  for (unsigned int i = layer; i < _layer_values.size(); i++)
215  {
216  if (_layer_has_value[i])
217  {
218  higher_layer = i;
219  break;
220  }
221  }
222 
223  for (int i = layer - 1; i >= 0; i--)
224  {
225  if (_layer_has_value[i])
226  {
227  lower_layer = i;
228  break;
229  }
230  }
231 
232  if (higher_layer == -1 && lower_layer == -1)
233  return 0; // TODO: We could error here but there are startup dependency problems
234 
235  switch (_sample_type)
236  {
237  case 0: // direct
238  {
239  if (higher_layer == -1) // Didn't find a higher layer
240  return _layer_values[lower_layer];
241 
242  if (unsigned(higher_layer) == layer) // constant in a layer
243  return _layer_values[higher_layer];
244 
245  if (lower_layer == -1) // Didn't find a lower layer
246  return _layer_values[higher_layer];
247 
248  return (_layer_values[higher_layer] + _layer_values[lower_layer]) / 2;
249  }
250  case 1: // interpolate
251  {
252  if (higher_layer == -1) // Didn't find a higher layer
253  return _layer_values[lower_layer];
254 
255  Real layer_length = (_direction_max - _direction_min) / _num_layers;
256  Real lower_coor = _direction_min;
257  Real lower_value = 0;
258  if (lower_layer != -1)
259  {
260  lower_coor += (lower_layer + 1) * layer_length;
261  lower_value = _layer_values[lower_layer];
262  }
263 
264  // Interpolate between the two points
265  Real higher_value = _layer_values[higher_layer];
266 
267  // Linear interpolation
268  return lower_value +
269  (higher_value - lower_value) * (p(_direction) - lower_coor) / layer_length;
270  }
271  case 2: // average
272  {
273  Real total = 0;
274  unsigned int num_values = 0;
275 
276  if (higher_layer != -1)
277  {
278  for (const auto i : make_range(_average_radius))
279  {
280  int current_layer = higher_layer + i;
281 
282  if ((size_t)current_layer >= _layer_values.size())
283  break;
284 
285  if (_layer_has_value[current_layer])
286  {
287  total += _layer_values[current_layer];
288  num_values += 1;
289  }
290  }
291  }
292 
293  if (lower_layer != -1)
294  {
295  for (const auto i : make_range(_average_radius))
296  {
297  int current_layer = lower_layer - i;
298 
299  if (current_layer < 0)
300  break;
301 
302  if (_layer_has_value[current_layer])
303  {
304  total += _layer_values[current_layer];
305  num_values += 1;
306  }
307  }
308  }
309 
310  return total / num_values;
311  }
312  default:
313  mooseError("Unknown sample type!");
314  }
315 }
316 
317 Real
318 LayeredBase::getLayerValue(unsigned int layer) const
319 {
320  if (layer >= _layer_values.size())
321  mooseError("Layer '", layer, "' not found in '", _layered_base_name, "'.");
322  return _layer_values[layer];
323 }
324 
325 void
327 {
329  getBounds();
330 
331  for (const auto i : index_range(_layer_values))
332  {
333  _layer_values[i] = 0.0;
334  _layer_has_value[i] = false;
335  }
336 }
337 
338 void
340 {
343 
344  if (_cumulative)
345  {
346  Real value = 0;
347 
349  for (const auto i : make_range(_num_layers))
350  {
351  value += getLayerValue(i);
352  setLayerValue(i, value);
353  }
354  else
355  for (int i = _num_layers - 1; i >= 0; i--)
356  {
357  value += getLayerValue(i);
358  setLayerValue(i, value);
359  }
360  }
361 }
362 
363 void
365 {
366  const auto & lb = dynamic_cast<const LayeredBase &>(y);
367  for (const auto i : index_range(_layer_values))
368  if (lb.layerHasValue(i))
369  setLayerValue(i, getLayerValue(i) + lb._layer_values[i]);
370 }
371 
372 unsigned int
373 LayeredBase::getLayer(const Point & p) const
374 {
375  Real direction_x = p(_direction);
376 
377  if (direction_x < _direction_min)
378  return 0;
379 
380  if (_interval_based)
381  {
382  unsigned int layer =
383  std::floor(((direction_x - _direction_min) / (_direction_max - _direction_min)) *
384  static_cast<Real>(_num_layers));
385 
386  if (layer >= _num_layers)
387  layer = _num_layers - 1;
388 
389  return layer;
390  }
391  else // Figure out what layer we are in from the bounds
392  {
393  // This finds the first entry in the vector that is larger than what we're looking for
394  std::vector<Real>::const_iterator one_higher =
395  std::upper_bound(_layer_bounds.begin(), _layer_bounds.end(), direction_x);
396 
397  if (one_higher == _layer_bounds.end())
398  {
399  return static_cast<unsigned int>(
400  _layer_bounds.size() -
401  2); // Just return the last layer. -2 because layers are "in-between" bounds
402  }
403  else if (one_higher == _layer_bounds.begin())
404  return 0; // Return the first layer
405  else
406  // The -1 is because the interval that we fall in is just _before_ the number that is bigger
407  // (which is what we found
408  return static_cast<unsigned int>(std::distance(_layer_bounds.begin(), one_higher - 1));
409  }
410 }
411 
412 void
414 {
415  _layer_centers.resize(_num_layers);
416 
417  if (_interval_based)
418  {
420 
421  for (const auto i : make_range(_num_layers))
422  _layer_centers[i] = (i + 0.5) * dx;
423  }
424  else
425  {
426  for (const auto i : make_range(_num_layers))
427  _layer_centers[i] = 0.5 * (_layer_bounds[i + 1] + _layer_bounds[i]);
428  }
429 }
430 
431 void
432 LayeredBase::setLayerValue(unsigned int layer, Real value)
433 {
434  _layer_values[layer] = value;
435  _layer_has_value[layer] = true;
436 }
437 
438 void
440 {
441  if (_layer_bounding_blocks.size() == 0)
442  {
443  BoundingBox bounding_box = MeshTools::create_bounding_box(_layered_base_subproblem.mesh());
444  _direction_min = bounding_box.min()(_direction);
445  _direction_max = bounding_box.max()(_direction);
446  }
447  else
448  {
449  _direction_min = std::numeric_limits<Real>::infinity();
450  _direction_max = -std::numeric_limits<Real>::infinity();
451 
453 
454  for (auto & elem_ptr : *mesh.getActiveLocalElementRange())
455  {
456  auto subdomain_id = elem_ptr->subdomain_id();
457 
458  if (std::find(_layer_bounding_blocks.begin(), _layer_bounding_blocks.end(), subdomain_id) ==
460  continue;
461 
462  for (auto & node : elem_ptr->node_ref_range())
463  {
466  }
467  }
468 
469  mesh.comm().min(_direction_min);
470  mesh.comm().max(_direction_max);
471  }
472 }
virtual MooseMesh & mesh()=0
std::vector< Real > & _layer_values
Value of the integral for each layer.
Definition: LayeredBase.h:139
bool _using_displaced_mesh
true if this object operates on the displaced mesh, otherwise false
Definition: LayeredBase.h:130
A class for creating restricted objects.
Definition: Restartable.h:28
std::string _layered_base_name
Name of this object.
Definition: LayeredBase.h:103
SubProblem & _layered_base_subproblem
Subproblem for the child object.
Definition: LayeredBase.h:152
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
std::vector< std::pair< R1, R2 > > get(const std::string &param1, const std::string &param2) const
Combine two vector parameters into a single vector of pairs.
void mooseWarning(Args &&... args)
Emit a warning message with the given stringified, concatenated args.
Definition: MooseError.h:336
T * get(const std::unique_ptr< T > &u)
The MooseUtils::get() specializations are used to support making forwards-compatible code changes fro...
Definition: MooseUtils.h:1155
bool _interval_based
Whether or not this object is based on equally spaced intervals or "bounds".
Definition: LayeredBase.h:115
bool _cumulative
Whether the values are cumulative over the layers.
Definition: LayeredBase.h:145
MeshBase & mesh
void getBounds()
Compute bounds, restricted to blocks if given.
Definition: LayeredBase.C:439
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
const Parallel::Communicator & comm() const
std::vector< int > & _layer_has_value
Whether or not each layer has had any value summed into it.
Definition: LayeredBase.h:142
static InputParameters validParams()
Definition: LayeredBase.C:22
const InputParameters & _layered_base_params
Params for this object.
Definition: LayeredBase.h:106
std::vector< Real > _layer_centers
center coordinates of each layer
Definition: LayeredBase.h:133
std::vector< SubdomainID > getSubdomainIDs(const std::vector< SubdomainName > &subdomain_names) const
Get the associated subdomainIDs for the subdomain names that are passed in.
Definition: MooseMesh.C:1734
void addRequiredParam(const std::string &name, const std::string &doc_string)
This method adds a parameter and documentation string to the InputParameters object that will be extr...
InputParameters emptyInputParameters()
auto max(const L &left, const R &right)
virtual void threadJoin(const UserObject &y)
Definition: LayeredBase.C:364
unsigned int _direction
The component direction the layers are going in. We cache this for speed (so we&#39;re not always going t...
Definition: LayeredBase.h:112
std::vector< SubdomainID > _layer_bounding_blocks
List of SubdomainIDs, if given.
Definition: LayeredBase.h:155
unsigned int _average_radius
How many layers both above and below the found layer will be used in the average. ...
Definition: LayeredBase.h:127
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
virtual Real getLayerValue(unsigned int layer) const
Get the value for a given layer.
Definition: LayeredBase.C:318
void computeLayerCenters()
Compute the center points for each layer.
Definition: LayeredBase.C:413
void setLayerValue(unsigned int layer, Real value)
Set the value for a particular layer.
Definition: LayeredBase.C:432
virtual Real integralValue(const Point &p) const
Given a Point return the integral value associated with the layer that point falls in...
Definition: LayeredBase.C:207
const bool _positive_cumulative_direction
Whether the cumulative values should be summed in the positive or negative direction.
Definition: LayeredBase.h:148
virtual unsigned int getLayer(const Point &p) const
Helper function to return the layer the point lies in.
Definition: LayeredBase.C:373
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
Definition: MooseMesh.h:88
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:33
virtual void initialize()
Definition: LayeredBase.C:326
unsigned int _num_layers
Number of layers to split the mesh into.
Definition: LayeredBase.h:118
virtual void finalize()
Definition: LayeredBase.C:339
std::vector< Real > _layer_bounds
The boundaries of the layers.
Definition: LayeredBase.h:121
bool isParamSetByUser(const std::string &name) const
Method returns true if the parameter was set by the user.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Generic class for solving transient nonlinear problems.
Definition: SubProblem.h:78
void max(const T &r, T &o, Request &req) const
This base class computes volume integrals of a variable storing partial sums for the specified number...
Definition: LayeredBase.h:36
IntRange< T > make_range(T beg, T end)
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an optional parameter and a documentation string to the InputParameters object...
Real _direction_min
Definition: LayeredBase.h:135
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
unsigned int _sample_type
How to sample the values.
Definition: LayeredBase.h:124
LayeredBase(const InputParameters &parameters)
Definition: LayeredBase.C:83
Real _direction_max
Definition: LayeredBase.h:136
auto min(const L &left, const R &right)
void ErrorVector unsigned int
auto index_range(const T &sizable)
Base class for user-specific data.
Definition: UserObject.h:40
unsigned int THREAD_ID
Definition: MooseTypes.h:209
bool _has_direction_max_min
whether the max/min coordinate in the direction is known from user input
Definition: LayeredBase.h:158
void addParamNamesToGroup(const std::string &space_delim_names, const std::string group_name)
This method takes a space delimited list of parameter names and adds them to the specified group name...
bool isParamValid(const std::string &name) const
This method returns parameters that have been initialized in one fashion or another, i.e.