www.mooseframework.org
LayeredBase.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 #include "LayeredBase.h"
16 
17 // MOOSE includes
18 #include "MooseEnum.h"
19 #include "MooseMesh.h"
20 #include "SubProblem.h"
21 #include "UserObject.h"
22 
23 #include "libmesh/mesh_tools.h"
24 #include "libmesh/point.h"
25 
26 template <>
29 {
31  MooseEnum directions("x y z");
32 
33  params.addRequiredParam<MooseEnum>("direction", directions, "The direction of the layers.");
34  params.addParam<unsigned int>("num_layers", "The number of layers.");
35  params.addParam<std::vector<Real>>("bounds",
36  "The 'bounding' positions of the layers i.e.: '0, "
37  "1.2, 3.7, 4.2' will mean 3 layers between those "
38  "positions.");
39 
40  MooseEnum sample_options("direct interpolate average", "direct");
41  params.addParam<MooseEnum>("sample_type",
42  sample_options,
43  "How to sample the layers. 'direct' means get the value of the layer "
44  "the point falls in directly (or average if that layer has no value). "
45  " 'interpolate' does a linear interpolation between the two closest "
46  "layers. 'average' averages the two closest layers.");
47 
48  params.addParam<unsigned int>("average_radius",
49  1,
50  "When using 'average' sampling this is how "
51  "the number of values both above and below "
52  "the layer that will be averaged.");
53 
54  params.addParam<bool>(
55  "cumulative",
56  false,
57  "When true the value in each layer is the sum of the values up to and including that layer");
58 
59  return params;
60 }
61 
63  : _layered_base_name(parameters.get<std::string>("_object_name")),
64  _layered_base_params(parameters),
65  _direction_enum(parameters.get<MooseEnum>("direction")),
66  _direction(_direction_enum),
67  _sample_type(parameters.get<MooseEnum>("sample_type")),
68  _average_radius(parameters.get<unsigned int>("average_radius")),
69  _layered_base_subproblem(*parameters.get<SubProblem *>("_subproblem")),
70  _cumulative(parameters.get<bool>("cumulative"))
71 {
72  if (_layered_base_params.isParamValid("num_layers") &&
74  mooseError("'bounds' and 'num_layers' cannot both be set in ", _layered_base_name);
75 
76  if (_layered_base_params.isParamValid("num_layers"))
77  {
78  _num_layers = _layered_base_params.get<unsigned int>("num_layers");
79  _interval_based = true;
80  }
81  else if (_layered_base_params.isParamValid("bounds"))
82  {
83  _interval_based = false;
84 
85  _layer_bounds = _layered_base_params.get<std::vector<Real>>("bounds");
86 
87  // Make sure the bounds are sorted - we're going to depend on this
88  std::sort(_layer_bounds.begin(), _layer_bounds.end());
89 
90  _num_layers = _layer_bounds.size() - 1; // Layers are only in-between the bounds
91  }
92  else
93  mooseError("One of 'bounds' or 'num_layers' must be specified for ", _layered_base_name);
94 
95  if (!_interval_based && _sample_type == 1)
96  mooseError("'sample_type = interpolate' not supported with 'bounds' in ", _layered_base_name);
97 
98  BoundingBox bounding_box = MeshTools::create_bounding_box(_layered_base_subproblem.mesh());
99  _layer_values.resize(_num_layers);
101 
102  _direction_min = bounding_box.min()(_direction);
103  _direction_max = bounding_box.max()(_direction);
104 }
105 
106 Real
108 {
109  unsigned int layer = getLayer(p);
110 
111  int higher_layer = -1;
112  int lower_layer = -1;
113 
114  for (unsigned int i = layer; i < _layer_values.size(); i++)
115  {
116  if (_layer_has_value[i])
117  {
118  higher_layer = i;
119  break;
120  }
121  }
122 
123  for (int i = layer - 1; i >= 0; i--)
124  {
125  if (_layer_has_value[i])
126  {
127  lower_layer = i;
128  break;
129  }
130  }
131 
132  if (higher_layer == -1 && lower_layer == -1)
133  return 0; // TODO: We could error here but there are startup dependency problems
134 
135  switch (_sample_type)
136  {
137  case 0: // direct
138  {
139  if (higher_layer == -1) // Didn't find a higher layer
140  return _layer_values[lower_layer];
141 
142  if (unsigned(higher_layer) == layer) // constant in a layer
143  return _layer_values[higher_layer];
144 
145  if (lower_layer == -1) // Didn't find a lower layer
146  return _layer_values[higher_layer];
147 
148  return (_layer_values[higher_layer] + _layer_values[lower_layer]) / 2;
149  }
150  case 1: // interpolate
151  {
152  if (higher_layer == -1) // Didn't find a higher layer
153  return _layer_values[lower_layer];
154 
155  Real layer_length = (_direction_max - _direction_min) / _num_layers;
156  Real lower_coor = _direction_min;
157  Real lower_value = 0;
158  if (lower_layer != -1)
159  {
160  lower_coor += (lower_layer + 1) * layer_length;
161  lower_value = _layer_values[lower_layer];
162  }
163 
164  // Interpolate between the two points
165  Real higher_value = _layer_values[higher_layer];
166 
167  // Linear interpolation
168  return lower_value +
169  (higher_value - lower_value) * (p(_direction) - lower_coor) / layer_length;
170  }
171  case 2: // average
172  {
173  Real total = 0;
174  unsigned int num_values = 0;
175 
176  if (higher_layer != -1)
177  {
178  for (unsigned int i = 0; i < _average_radius; i++)
179  {
180  int current_layer = higher_layer + i;
181 
182  if ((size_t)current_layer >= _layer_values.size())
183  break;
184 
185  if (_layer_has_value[current_layer])
186  {
187  total += _layer_values[current_layer];
188  num_values += 1;
189  }
190  }
191  }
192 
193  if (lower_layer != -1)
194  {
195  for (unsigned int i = 0; i < _average_radius; i++)
196  {
197  int current_layer = lower_layer - i;
198 
199  if (current_layer < 0)
200  break;
201 
202  if (_layer_has_value[current_layer])
203  {
204  total += _layer_values[current_layer];
205  num_values += 1;
206  }
207  }
208  }
209 
210  return total / num_values;
211  }
212  default:
213  mooseError("Unknown sample type!");
214  }
215 }
216 
217 Real
218 LayeredBase::getLayerValue(unsigned int layer) const
219 {
220  if (layer >= _layer_values.size())
221  mooseError("Layer '", layer, "' not found in '", _layered_base_name, "'.");
222  return _layer_values[layer];
223 }
224 
225 void
227 {
228  for (unsigned int i = 0; i < _layer_values.size(); i++)
229  {
230  _layer_values[i] = 0.0;
231  _layer_has_value[i] = false;
232  }
233 }
234 
235 void
237 {
240 
241  if (_cumulative)
242  {
243  Real value = 0;
244 
245  for (unsigned i = 0; i < _num_layers; ++i)
246  {
247  value += getLayerValue(i);
248  setLayerValue(i, value);
249  }
250  }
251 }
252 
253 void
255 {
256  const LayeredBase & lb = dynamic_cast<const LayeredBase &>(y);
257  for (unsigned int i = 0; i < _layer_values.size(); i++)
258  if (lb.layerHasValue(i))
260 }
261 
262 unsigned int
263 LayeredBase::getLayer(Point p) const
264 {
265  Real direction_x = p(_direction);
266 
267  if (direction_x < _direction_min)
268  return 0;
269 
270  if (_interval_based)
271  {
272  unsigned int layer =
273  std::floor(((direction_x - _direction_min) / (_direction_max - _direction_min)) *
274  static_cast<Real>(_num_layers));
275 
276  if (layer >= _num_layers)
277  layer = _num_layers - 1;
278 
279  return layer;
280  }
281  else // Figure out what layer we are in from the bounds
282  {
283  // This finds the first entry in the vector that is larger than what we're looking for
284  std::vector<Real>::const_iterator one_higher =
285  std::upper_bound(_layer_bounds.begin(), _layer_bounds.end(), direction_x);
286 
287  if (one_higher == _layer_bounds.end())
288  {
289  return static_cast<unsigned int>(
290  _layer_bounds.size() -
291  2); // Just return the last layer. -2 because layers are "in-between" bounds
292  }
293  else if (one_higher == _layer_bounds.begin())
294  return 0; // Return the first layer
295  else
296  // The -1 is because the interval that we fall in is just _before_ the number that is bigger
297  // (which is what we found
298  return static_cast<unsigned int>(std::distance(_layer_bounds.begin(), one_higher - 1));
299  }
300 }
301 
302 void
303 LayeredBase::setLayerValue(unsigned int layer, Real value)
304 {
305  _layer_values[layer] = value;
306  _layer_has_value[layer] = true;
307 }
virtual MooseMesh & mesh()=0
std::string _layered_base_name
Name of this object.
Definition: LayeredBase.h:88
SubProblem & _layered_base_subproblem
Subproblem for the child object.
Definition: LayeredBase.h:125
bool layerHasValue(unsigned int layer) const
Whether or not a layer has a value.
Definition: LayeredBase.h:85
void mooseError(Args &&...args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:182
bool _interval_based
Whether or not this object is based on equally spaced intervals or "bounds".
Definition: LayeredBase.h:100
virtual unsigned int getLayer(Point p) const
Helper function to return the layer the point lies in.
Definition: LayeredBase.C:263
bool _cumulative
Whether the values are cumulative over the layers.
Definition: LayeredBase.h:128
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
const InputParameters & _layered_base_params
Params for this object.
Definition: LayeredBase.h:91
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()
virtual void threadJoin(const UserObject &y)
Definition: LayeredBase.C:254
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:97
virtual Real getLayerValue(unsigned int layer) const
Get the value for a given layer.
Definition: LayeredBase.C:218
std::vector< Real > _layer_values
Value of the integral for each layer.
Definition: LayeredBase.h:119
unsigned int _average_radius
How many layers both above and below the found layer will be used in the average. ...
Definition: LayeredBase.h:112
bool isParamValid(const std::string &name) const
This method returns parameters that have been initialized in one fashion or another, i.e.
void setLayerValue(unsigned int layer, Real value)
Set the value for a particular layer.
Definition: LayeredBase.C:303
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:37
virtual void initialize()
Definition: LayeredBase.C:226
unsigned int _num_layers
Number of layers to split the mesh into.
Definition: LayeredBase.h:103
InputParameters validParams< LayeredBase >()
Definition: LayeredBase.C:28
virtual void finalize()
Definition: LayeredBase.C:236
std::vector< Real > _layer_bounds
The boundaries of the layers.
Definition: LayeredBase.h:106
Generic class for solving transient nonlinear problems.
Definition: SubProblem.h:53
virtual Real integralValue(Point p) const
Given a Point return the integral value associated with the layer that point falls in...
Definition: LayeredBase.C:107
This UserObject computes volume integrals of a variable storing partial sums for the specified number...
Definition: LayeredBase.h:44
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an option parameter and a documentation string to the InputParameters object...
Real _direction_min
Definition: LayeredBase.h:114
unsigned int _sample_type
How to sample the values.
Definition: LayeredBase.h:109
LayeredBase(const InputParameters &parameters)
Definition: LayeredBase.C:62
Real _direction_max
Definition: LayeredBase.h:115
Base class for user-specific data.
Definition: UserObject.h:42
std::vector< bool > _layer_has_value
Whether or not each layer has had any value summed into it.
Definition: LayeredBase.h:122