Line data Source code
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 "AverageSectionValueSampler.h"
11 : #include "MooseMesh.h"
12 : #include "SystemBase.h"
13 : #include "Conversion.h"
14 : #include <limits>
15 :
16 : registerMooseObject("SolidMechanicsApp", AverageSectionValueSampler);
17 :
18 : InputParameters
19 92 : AverageSectionValueSampler::validParams()
20 : {
21 92 : InputParameters params = GeneralVectorPostprocessor::validParams();
22 92 : params.addClassDescription("Compute the section's variable average in three-dimensions "
23 : "given a user-defined definition of the cross section.");
24 184 : params.addParam<std::vector<SubdomainName>>(
25 : "block",
26 : "The list of blocks in which to search for cross sectional nodes to compute the variable "
27 : "average.");
28 184 : params.addRequiredParam<Point>("axis_direction", "Direction of the structural component's axis");
29 184 : params.addRequiredParam<Point>("reference_point",
30 : "Structural component reference starting point from which the "
31 : "input parameter 'lengths' applies.");
32 184 : params.addParam<Real>("cross_section_maximum_radius",
33 184 : std::numeric_limits<Real>::max(),
34 : "Radial distance with respect to the body axis within which nodes are "
35 : "considered to belong to this "
36 : "structural component. Used to disambiguate multiple components that share "
37 : "the same mesh block.");
38 184 : params.addRequiredParam<std::vector<VariableName>>(
39 : "variables", "Variables for the cross section output. These variables must be nodal.");
40 184 : params.addDeprecatedParam<std::vector<Real>>(
41 : "lengths",
42 : {},
43 : "Positions along axis from reference_point at which to compute average values.",
44 : "Use the 'positions' parameter instead");
45 184 : params.addParam<std::vector<Real>>(
46 : "positions",
47 : {},
48 : "Positions along axis from reference_point at which to compute average values.");
49 92 : params.addParam<RealVectorValue>(
50 : "symmetry_plane",
51 92 : RealVectorValue(0, 0, 0),
52 : "Vector normal to a symmetry plane, used if the section has a symmetry plane through it. "
53 : "Causes the variables to be treated as components of a vector and zeros out the compomnent "
54 : "in the direction normal to symmetry_plane. Three variables must be defined in 'variables' "
55 : "to use this option.");
56 184 : params.addParam<Real>("tolerance",
57 184 : 1.0e-6,
58 : "Maximum axial distance of nodes from the specified axial lengths to "
59 : "consider them in the cross-section average");
60 184 : params.addParam<bool>(
61 : "require_equal_node_counts",
62 184 : true,
63 : "Whether to require the number of nodes at each axial location to be equal");
64 92 : return params;
65 0 : }
66 :
67 48 : AverageSectionValueSampler::AverageSectionValueSampler(const InputParameters & parameters)
68 : : GeneralVectorPostprocessor(parameters),
69 48 : _mesh(_app.actionWarehouse().mesh()),
70 96 : _variables(getParam<std::vector<VariableName>>("variables")),
71 96 : _direction(getParam<Point>("axis_direction")),
72 96 : _reference_point(getParam<Point>("reference_point")),
73 192 : _positions(isParamSetByUser("lengths") ? getParam<std::vector<Real>>("lengths")
74 144 : : getParam<std::vector<Real>>("positions")),
75 96 : _have_symmetry_plane(isParamSetByUser("symmetry_plane")),
76 96 : _symmetry_plane(getParam<RealVectorValue>("symmetry_plane")),
77 48 : _automatically_locate_positions(_positions.size() == 0),
78 96 : _tolerance(getParam<Real>("tolerance")),
79 48 : _number_of_nodes(_positions.size()),
80 96 : _cross_section_maximum_radius(getParam<Real>("cross_section_maximum_radius")),
81 96 : _require_equal_node_counts(getParam<bool>("require_equal_node_counts")),
82 96 : _need_mesh_initializations(true)
83 : {
84 48 : if (parameters.isParamSetByUser("lengths") && parameters.isParamSetByUser("positions"))
85 0 : paramError("lengths", "The 'lengths' and 'positions' parameters cannot both be set.");
86 :
87 96 : if (parameters.isParamSetByUser("lengths") && _positions.size() == 0)
88 0 : paramError("lengths", "If 'lengths' is specified, at least one value must be provided");
89 :
90 96 : if (parameters.isParamSetByUser("positions") && _positions.size() == 0)
91 2 : paramError("positions", "If 'positions' is specified, at least one value must be provided");
92 :
93 46 : if (_mesh->dimension() != 3)
94 0 : mooseError("The AverageSectionValueSampler postprocessor can only be used with three "
95 : "dimensional meshes.");
96 :
97 46 : if (!MooseUtils::absoluteFuzzyEqual(_direction.norm_sq(), 1.0))
98 0 : paramError("axis_direction", "Axis direction must be a unit vector.");
99 :
100 46 : _output_vector.resize(_variables.size() + 2);
101 180 : for (const auto j : make_range(_variables.size()))
102 134 : _output_vector[j] = &declareVector(_variables[j]);
103 : const auto pos_idx = _variables.size();
104 46 : _output_vector[pos_idx] = &declareVector("axial_position");
105 46 : _output_vector[pos_idx + 1] = &declareVector("node_count");
106 :
107 178 : for (const auto j : make_range(_variables.size()))
108 : {
109 134 : if (!_sys.hasVariable(_variables[j]))
110 2 : paramError("variables",
111 : "One or more of the specified variables do not exist in this system.");
112 132 : const MooseVariable & variable = _sys.getFieldVariable<Real>(_tid, _variables[j]);
113 132 : if (!variable.isNodal())
114 0 : paramError("variables",
115 : "One or more of the specified variables are not defined at the nodes.");
116 132 : _var_numbers.push_back(variable.number());
117 : }
118 :
119 44 : if (_have_symmetry_plane)
120 : {
121 6 : const auto symm_plane_norm = _symmetry_plane.norm();
122 6 : if (MooseUtils::absoluteFuzzyEqual(symm_plane_norm, 0.0))
123 0 : mooseError("Vector defined in 'symmetry_plane' cannot have a norm of zero");
124 6 : _symmetry_plane /= _symmetry_plane.norm();
125 6 : if (_variables.size() != 3)
126 0 : paramError("variables",
127 : "If 'symmetry_plane' is prescribed, three variables must be provided (for the x, "
128 : "y, z vector components, in that order)");
129 6 : mooseInfo("In AverageSectionValueSampler " + name() +
130 12 : ": \nTreating the variables as vector components (" + _variables[0] + ", " +
131 18 : _variables[1] + ", " + _variables[2] +
132 : ") and zeroing out the component normal to the prescribed 'symmetry_plane'\n");
133 : }
134 44 : }
135 :
136 : void
137 44 : AverageSectionValueSampler::initialSetup()
138 : {
139 44 : }
140 :
141 : void
142 8 : AverageSectionValueSampler::meshChanged()
143 : {
144 8 : _need_mesh_initializations = true;
145 8 : }
146 :
147 : void
148 52 : AverageSectionValueSampler::initialize()
149 : {
150 52 : if (_need_mesh_initializations)
151 : {
152 : // Option 1: locate positions first, then check counts
153 52 : if (_automatically_locate_positions)
154 34 : automaticallyLocatePositions();
155 52 : _need_mesh_initializations = false;
156 : }
157 :
158 312 : for (const auto j : make_range(_variables.size() + 2))
159 : {
160 260 : _output_vector[j]->clear();
161 260 : _output_vector[j]->resize(_positions.size());
162 : }
163 52 : _number_of_nodes.clear();
164 52 : _number_of_nodes.resize(_positions.size());
165 52 : }
166 :
167 : void
168 52 : AverageSectionValueSampler::finalize()
169 : {
170 208 : for (const auto j : make_range(_variables.size()))
171 156 : _communicator.sum(*(_output_vector[j]));
172 :
173 52 : _communicator.sum(_number_of_nodes);
174 :
175 470 : for (const auto li : index_range(_positions))
176 418 : if (_number_of_nodes[li] < 1)
177 0 : mooseError("No nodes were found in AverageSectionValueSampler postprocessor.");
178 :
179 470 : for (const auto li : index_range(_positions))
180 1672 : for (const auto j : make_range(_variables.size()))
181 1254 : (*_output_vector[j])[li] /= _number_of_nodes[li];
182 :
183 : const auto pos_idx = _variables.size();
184 52 : (*_output_vector[pos_idx]) = _positions;
185 52 : std::vector<Real> num_nodes_real(_number_of_nodes.begin(), _number_of_nodes.end());
186 52 : (*_output_vector[pos_idx + 1]) = num_nodes_real;
187 :
188 52 : if (_require_equal_node_counts)
189 : {
190 314 : for (const auto li : index_range(_number_of_nodes))
191 : {
192 276 : if (_number_of_nodes[li] != _number_of_nodes[0])
193 : {
194 2 : std::ostringstream pos_out;
195 26 : for (const auto li2 : index_range(_number_of_nodes))
196 24 : pos_out << std::setw(10) << _positions[li2] << std::setw(10) << _number_of_nodes[li2]
197 24 : << "\n";
198 4 : mooseError(
199 : "Node counts do not match for all axial positions. If this behavior "
200 : "is desired to accommodate nonuniform meshes, set "
201 2 : "'require_equal_node_counts=false'\n Axial Node\n Position Count\n" +
202 0 : pos_out.str());
203 0 : }
204 : }
205 : }
206 50 : }
207 :
208 : void
209 52 : AverageSectionValueSampler::execute()
210 : {
211 104 : _block_ids = _mesh->getSubdomainIDs(getParam<std::vector<SubdomainName>>("block"));
212 :
213 52 : auto * active_nodes = _mesh->getActiveSemiLocalNodeRange();
214 :
215 : std::vector<std::vector<Real>> output_vector_partial(_variables.size(),
216 52 : std::vector<Real>(_positions.size()));
217 :
218 52 : const NumericVector<Number> & _sol = *_sys.currentSolution();
219 :
220 24852 : for (const auto & node : *active_nodes)
221 : {
222 24800 : const std::set<SubdomainID> & node_blk_ids = _mesh->getNodeBlockIds(*node);
223 :
224 207076 : for (const auto li : index_range(_positions))
225 : {
226 349680 : for (const auto i : _block_ids)
227 : {
228 182276 : if (node_blk_ids.find(i) != node_blk_ids.end())
229 : {
230 : // Check if node is within tolerance of user-prescribed plane
231 182276 : if (std::abs(axialPosition(*node) - _positions[li]) > _tolerance)
232 167404 : continue;
233 15320 : if ((*node).processor_id() != processor_id())
234 448 : continue;
235 :
236 : // Retrieve nodal variables
237 14872 : std::vector<Real> this_node_vars(_var_numbers.size());
238 59488 : for (const auto j : make_range(_var_numbers.size()))
239 44616 : this_node_vars[j] = _sol(node->dof_number(_sys.number(), _var_numbers[j], 0));
240 :
241 14872 : if (_have_symmetry_plane)
242 : {
243 : RealVectorValue this_node_vec_var(
244 : this_node_vars[0], this_node_vars[1], this_node_vars[2]);
245 1144 : this_node_vec_var -= _symmetry_plane * this_node_vec_var * _symmetry_plane;
246 1144 : this_node_vars[0] = this_node_vec_var(0);
247 1144 : this_node_vars[1] = this_node_vec_var(1);
248 1144 : this_node_vars[2] = this_node_vec_var(2);
249 : }
250 :
251 59488 : for (const auto j : make_range(_var_numbers.size()))
252 44616 : output_vector_partial[j][li] += this_node_vars[j];
253 :
254 14872 : _number_of_nodes[li]++;
255 : break;
256 : }
257 : }
258 : }
259 : }
260 :
261 208 : for (const auto j : make_range(_variables.size()))
262 156 : *_output_vector[j] = output_vector_partial[j];
263 52 : }
264 :
265 : Real
266 195996 : AverageSectionValueSampler::axialPosition(const Node & node) const
267 : {
268 : // Compute node location w.r.t. structural component length
269 : const Point relative_distance{
270 195996 : node(0) - _reference_point(0), node(1) - _reference_point(1), node(2) - _reference_point(2)};
271 :
272 : const Real axial_position = _direction * relative_distance;
273 : const Real in_plane_distance =
274 195996 : (relative_distance - relative_distance * _direction * _direction).norm();
275 :
276 : // If the in-plane distance is greater than the specified cross-section radius, the point is not
277 : // in this component
278 195996 : if (in_plane_distance > _cross_section_maximum_radius)
279 : return std::numeric_limits<Real>::max();
280 :
281 187548 : return axial_position;
282 : }
283 :
284 : void
285 34 : AverageSectionValueSampler::automaticallyLocatePositions()
286 : {
287 34 : _positions.clear();
288 :
289 : // Data structure used to collect the locations with nodes on each processor,
290 : // and gather in parallel.
291 : std::vector<Real> pos_vec;
292 :
293 : // Data structure used to store parallel-gathered locations
294 : std::set<Real> pos_set;
295 :
296 68 : _block_ids = _mesh->getSubdomainIDs(getParam<std::vector<SubdomainName>>("block"));
297 34 : auto * nodes = _mesh->getLocalNodeRange();
298 :
299 13754 : for (const auto node : *nodes)
300 : {
301 13720 : const std::set<SubdomainID> & node_blk_ids = _mesh->getNodeBlockIds(*node);
302 :
303 27440 : for (const auto i : _block_ids)
304 : {
305 13720 : if (node_blk_ids.find(i) != node_blk_ids.end())
306 : {
307 : bool found_match = false;
308 13720 : Real axial_position = axialPosition(*node);
309 13720 : if (axial_position == std::numeric_limits<Real>::max()) // Node not in this component
310 0 : continue;
311 82320 : for (auto & pos : pos_vec)
312 : {
313 81938 : const Real dist_to_plane = std::abs(axial_position - pos);
314 81938 : if (dist_to_plane <= _tolerance)
315 : {
316 : found_match = true;
317 : break;
318 : }
319 : }
320 13720 : if (!found_match)
321 382 : pos_vec.emplace_back(axial_position);
322 : }
323 : }
324 : }
325 :
326 34 : _communicator.allgather(pos_vec);
327 :
328 640 : for (const auto & posv : pos_vec)
329 : {
330 : bool found_match = false;
331 3708 : for (auto & poss : pos_set)
332 : {
333 3326 : if (std::abs(posv - poss) < _tolerance)
334 : {
335 : found_match = true;
336 : break;
337 : }
338 : }
339 606 : if (!found_match)
340 382 : pos_set.insert(posv);
341 : }
342 416 : for (const auto & poss : pos_set)
343 382 : _positions.emplace_back(poss);
344 34 : }
|
