www.mooseframework.org
LayeredBase.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
21 template <>
24 {
26  MooseEnum directions("x y z");
27 
28  params.addRequiredParam<MooseEnum>("direction", directions, "The direction of the layers.");
29  params.addParam<unsigned int>("num_layers", "The number of layers.");
30  params.addParam<std::vector<Real>>("bounds",
31  "The 'bounding' positions of the layers i.e.: '0, "
32  "1.2, 3.7, 4.2' will mean 3 layers between those "
33  "positions.");
34 
35  MooseEnum sample_options("direct interpolate average", "direct");
36  params.addParam<MooseEnum>("sample_type",
37  sample_options,
38  "How to sample the layers. 'direct' means get the value of the layer "
39  "the point falls in directly (or average if that layer has no value). "
40  " 'interpolate' does a linear interpolation between the two closest "
41  "layers. 'average' averages the two closest layers.");
42 
43  params.addParam<unsigned int>("average_radius",
44  1,
45  "When using 'average' sampling this is how "
46  "the number of values both above and below "
47  "the layer that will be averaged.");
48 
49  params.addParam<bool>(
50  "cumulative",
51  false,
52  "When true the value in each layer is the sum of the values up to and including that layer");
53 
54  params.addParam<std::vector<SubdomainName>>(
55  "block", "The list of block ids (SubdomainID) that this object will be applied");
56 
57  params.addParam<std::vector<SubdomainName>>("layer_bounding_block",
58  "List of block ids (SubdomainID) that are used to "
59  "determine the upper and lower geometric bounds for "
60  "all layers. If this is not specified, the ids "
61  "specified in 'block' are used for this purpose.");
62 
63  return params;
64 }
65 
67  : Restartable(parameters.getCheckedPointerParam<SubProblem *>("_subproblem")->getMooseApp(),
68  parameters.get<std::string>("_object_name") + "_layered_base",
69  "LayeredBase",
70  parameters.get<THREAD_ID>("_tid")),
71  _layered_base_name(parameters.get<std::string>("_object_name")),
72  _layered_base_params(parameters),
73  _direction_enum(parameters.get<MooseEnum>("direction")),
74  _direction(_direction_enum),
75  _sample_type(parameters.get<MooseEnum>("sample_type")),
76  _average_radius(parameters.get<unsigned int>("average_radius")),
77  _using_displaced_mesh(_layered_base_params.get<bool>("use_displaced_mesh")),
78  _layer_values(declareRestartableData<std::vector<Real>>("layer_values")),
79  _layer_has_value(declareRestartableData<std::vector<int>>("layer_has_value")),
80  _layered_base_subproblem(*parameters.getCheckedPointerParam<SubProblem *>("_subproblem")),
81  _cumulative(parameters.get<bool>("cumulative")),
82  _layer_bounding_blocks()
83 {
84  if (_layered_base_params.isParamValid("num_layers") &&
86  mooseError("'bounds' and 'num_layers' cannot both be set in ", _layered_base_name);
87 
88  if (_layered_base_params.isParamValid("num_layers"))
89  {
90  _num_layers = _layered_base_params.get<unsigned int>("num_layers");
91  _interval_based = true;
92  }
93  else if (_layered_base_params.isParamValid("bounds"))
94  {
95  _interval_based = false;
96 
97  _layer_bounds = _layered_base_params.get<std::vector<Real>>("bounds");
98 
99  // Make sure the bounds are sorted - we're going to depend on this
100  std::sort(_layer_bounds.begin(), _layer_bounds.end());
101 
102  _num_layers = _layer_bounds.size() - 1; // Layers are only in-between the bounds
103  }
104  else
105  mooseError("One of 'bounds' or 'num_layers' must be specified for ", _layered_base_name);
106 
107  if (!_interval_based && _sample_type == 1)
108  mooseError("'sample_type = interpolate' not supported with 'bounds' in ", _layered_base_name);
109 
110  if (_layered_base_params.isParamValid("layer_bounding_block"))
112  _layered_base_params.get<std::vector<SubdomainName>>("layer_bounding_block"));
113  else if (_layered_base_params.isParamValid("block"))
115  _layered_base_params.get<std::vector<SubdomainName>>("block"));
116 
117  _layer_values.resize(_num_layers);
119 
120  getBounds();
121 }
122 
123 Real
125 {
126  unsigned int layer = getLayer(p);
127 
128  int higher_layer = -1;
129  int lower_layer = -1;
130 
131  for (unsigned int i = layer; i < _layer_values.size(); i++)
132  {
133  if (_layer_has_value[i])
134  {
135  higher_layer = i;
136  break;
137  }
138  }
139 
140  for (int i = layer - 1; i >= 0; i--)
141  {
142  if (_layer_has_value[i])
143  {
144  lower_layer = i;
145  break;
146  }
147  }
148 
149  if (higher_layer == -1 && lower_layer == -1)
150  return 0; // TODO: We could error here but there are startup dependency problems
151 
152  switch (_sample_type)
153  {
154  case 0: // direct
155  {
156  if (higher_layer == -1) // Didn't find a higher layer
157  return _layer_values[lower_layer];
158 
159  if (unsigned(higher_layer) == layer) // constant in a layer
160  return _layer_values[higher_layer];
161 
162  if (lower_layer == -1) // Didn't find a lower layer
163  return _layer_values[higher_layer];
164 
165  return (_layer_values[higher_layer] + _layer_values[lower_layer]) / 2;
166  }
167  case 1: // interpolate
168  {
169  if (higher_layer == -1) // Didn't find a higher layer
170  return _layer_values[lower_layer];
171 
172  Real layer_length = (_direction_max - _direction_min) / _num_layers;
173  Real lower_coor = _direction_min;
174  Real lower_value = 0;
175  if (lower_layer != -1)
176  {
177  lower_coor += (lower_layer + 1) * layer_length;
178  lower_value = _layer_values[lower_layer];
179  }
180 
181  // Interpolate between the two points
182  Real higher_value = _layer_values[higher_layer];
183 
184  // Linear interpolation
185  return lower_value +
186  (higher_value - lower_value) * (p(_direction) - lower_coor) / layer_length;
187  }
188  case 2: // average
189  {
190  Real total = 0;
191  unsigned int num_values = 0;
192 
193  if (higher_layer != -1)
194  {
195  for (unsigned int i = 0; i < _average_radius; i++)
196  {
197  int current_layer = higher_layer + i;
198 
199  if ((size_t)current_layer >= _layer_values.size())
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  if (lower_layer != -1)
211  {
212  for (unsigned int i = 0; i < _average_radius; i++)
213  {
214  int current_layer = lower_layer - i;
215 
216  if (current_layer < 0)
217  break;
218 
219  if (_layer_has_value[current_layer])
220  {
221  total += _layer_values[current_layer];
222  num_values += 1;
223  }
224  }
225  }
226 
227  return total / num_values;
228  }
229  default:
230  mooseError("Unknown sample type!");
231  }
232 }
233 
234 Real
235 LayeredBase::getLayerValue(unsigned int layer) const
236 {
237  if (layer >= _layer_values.size())
238  mooseError("Layer '", layer, "' not found in '", _layered_base_name, "'.");
239  return _layer_values[layer];
240 }
241 
242 void
244 {
246  getBounds();
247 
248  for (unsigned int i = 0; i < _layer_values.size(); i++)
249  {
250  _layer_values[i] = 0.0;
251  _layer_has_value[i] = false;
252  }
253 }
254 
255 void
257 {
260 
261  if (_cumulative)
262  {
263  Real value = 0;
264 
265  for (unsigned i = 0; i < _num_layers; ++i)
266  {
267  value += getLayerValue(i);
268  setLayerValue(i, value);
269  }
270  }
271 }
272 
273 void
275 {
276  const LayeredBase & lb = dynamic_cast<const LayeredBase &>(y);
277  for (unsigned int i = 0; i < _layer_values.size(); i++)
278  if (lb.layerHasValue(i))
280 }
281 
282 unsigned int
283 LayeredBase::getLayer(Point p) const
284 {
285  Real direction_x = p(_direction);
286 
287  if (direction_x < _direction_min)
288  return 0;
289 
290  if (_interval_based)
291  {
292  unsigned int layer =
293  std::floor(((direction_x - _direction_min) / (_direction_max - _direction_min)) *
294  static_cast<Real>(_num_layers));
295 
296  if (layer >= _num_layers)
297  layer = _num_layers - 1;
298 
299  return layer;
300  }
301  else // Figure out what layer we are in from the bounds
302  {
303  // This finds the first entry in the vector that is larger than what we're looking for
304  std::vector<Real>::const_iterator one_higher =
305  std::upper_bound(_layer_bounds.begin(), _layer_bounds.end(), direction_x);
306 
307  if (one_higher == _layer_bounds.end())
308  {
309  return static_cast<unsigned int>(
310  _layer_bounds.size() -
311  2); // Just return the last layer. -2 because layers are "in-between" bounds
312  }
313  else if (one_higher == _layer_bounds.begin())
314  return 0; // Return the first layer
315  else
316  // The -1 is because the interval that we fall in is just _before_ the number that is bigger
317  // (which is what we found
318  return static_cast<unsigned int>(std::distance(_layer_bounds.begin(), one_higher - 1));
319  }
320 }
321 
322 void
323 LayeredBase::setLayerValue(unsigned int layer, Real value)
324 {
325  _layer_values[layer] = value;
326  _layer_has_value[layer] = true;
327 }
328 
329 void
331 {
332  if (_layer_bounding_blocks.size() == 0)
333  {
334  BoundingBox bounding_box = MeshTools::create_bounding_box(_layered_base_subproblem.mesh());
335  _direction_min = bounding_box.min()(_direction);
336  _direction_max = bounding_box.max()(_direction);
337  }
338  else
339  {
340  _direction_min = std::numeric_limits<Real>::infinity();
341  _direction_max = -std::numeric_limits<Real>::infinity();
342 
344 
345  for (auto & elem_ptr : *mesh.getActiveLocalElementRange())
346  {
347  auto subdomain_id = elem_ptr->subdomain_id();
348 
349  if (std::find(_layer_bounding_blocks.begin(), _layer_bounding_blocks.end(), subdomain_id) ==
351  continue;
352 
353  for (auto & node : elem_ptr->node_ref_range())
354  {
355  _direction_min = std::min(_direction_min, node(_direction));
356  _direction_max = std::max(_direction_max, node(_direction));
357  }
358  }
359 
360  mesh.comm().min(_direction_min);
361  mesh.comm().max(_direction_max);
362  }
363 }
virtual MooseMesh & mesh()=0
std::vector< Real > & _layer_values
Value of the integral for each layer.
Definition: LayeredBase.h:123
bool _using_displaced_mesh
true if this object operates on the displaced mesh, otherwise false
Definition: LayeredBase.h:116
ConstElemRange * getActiveLocalElementRange()
Return pointers to range objects for various types of ranges (local nodes, boundary elems...
Definition: MooseMesh.C:737
A class for creating restricted objects.
Definition: Restartable.h:29
std::string _layered_base_name
Name of this object.
Definition: LayeredBase.h:89
SubProblem & _layered_base_subproblem
Subproblem for the child object.
Definition: LayeredBase.h:129
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:207
bool _interval_based
Whether or not this object is based on equally spaced intervals or "bounds".
Definition: LayeredBase.h:101
bool _cumulative
Whether the values are cumulative over the layers.
Definition: LayeredBase.h:132
void getBounds()
Compute bounds, restricted to blocks if given.
Definition: LayeredBase.C:330
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
std::vector< int > & _layer_has_value
Whether or not each layer has had any value summed into it.
Definition: LayeredBase.h:126
const InputParameters & _layered_base_params
Params for this object.
Definition: LayeredBase.h:92
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:274
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:98
std::vector< SubdomainID > _layer_bounding_blocks
List of SubdomainIDs, if given.
Definition: LayeredBase.h:135
unsigned int _average_radius
How many layers both above and below the found layer will be used in the average. ...
Definition: LayeredBase.h:113
bool layerHasValue(unsigned int layer) const
Whether or not a layer has a value.
Definition: LayeredBase.h:81
virtual Real integralValue(Point p) const
Given a Point return the integral value associated with the layer that point falls in...
Definition: LayeredBase.C:124
virtual Real getLayerValue(unsigned int layer) const
Get the value for a given layer.
Definition: LayeredBase.C:235
void setLayerValue(unsigned int layer, Real value)
Set the value for a particular layer.
Definition: LayeredBase.C:323
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
Definition: MooseMesh.h:74
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:31
virtual void initialize()
Definition: LayeredBase.C:243
unsigned int _num_layers
Number of layers to split the mesh into.
Definition: LayeredBase.h:104
InputParameters validParams< LayeredBase >()
Definition: LayeredBase.C:23
virtual void finalize()
Definition: LayeredBase.C:256
std::vector< Real > _layer_bounds
The boundaries of the layers.
Definition: LayeredBase.h:107
Generic class for solving transient nonlinear problems.
Definition: SubProblem.h:59
virtual unsigned int getLayer(Point p) const
Helper function to return the layer the point lies in.
Definition: LayeredBase.C:283
This UserObject (????? this isn&#39;t a userobject) computes volume integrals of a variable storing parti...
Definition: LayeredBase.h:40
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:118
std::vector< SubdomainID > getSubdomainIDs(const std::vector< SubdomainName > &subdomain_name) const
Get the associated subdomainIDs for the subdomain names that are passed in.
Definition: MooseMesh.C:1090
unsigned int _sample_type
How to sample the values.
Definition: LayeredBase.h:110
LayeredBase(const InputParameters &parameters)
Definition: LayeredBase.C:66
Real _direction_max
Definition: LayeredBase.h:119
Base class for user-specific data.
Definition: UserObject.h:37
unsigned int THREAD_ID
Definition: MooseTypes.h:161
bool isParamValid(const std::string &name) const
This method returns parameters that have been initialized in one fashion or another, i.e.