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 106 : AverageSectionValueSampler::validParams()
20 : {
21 106 : InputParameters params = GeneralVectorPostprocessor::validParams();
22 106 : params.addClassDescription("Compute the section's variable average in three-dimensions "
23 : "given a user-defined definition of the cross section.");
24 212 : 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 212 : params.addRequiredParam<Point>("axis_direction", "Direction of the structural component's axis");
29 212 : params.addRequiredParam<Point>("reference_point",
30 : "Structural component reference starting point from which the "
31 : "input parameter 'lengths' applies.");
32 212 : params.addParam<Real>("cross_section_maximum_radius",
33 212 : 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 212 : params.addRequiredParam<std::vector<VariableName>>(
39 : "variables", "Variables for the cross section output. These variables must be nodal.");
40 212 : 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 212 : params.addParam<std::vector<Real>>(
46 : "positions",
47 : {},
48 : "Positions along axis from reference_point at which to compute average values.");
49 106 : params.addParam<RealVectorValue>(
50 : "symmetry_plane",
51 106 : 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 212 : params.addParam<Real>("tolerance",
57 212 : 1.0e-6,
58 : "Maximum axial distance of nodes from the specified axial lengths to "
59 : "consider them in the cross-section average");
60 212 : params.addParam<bool>(
61 : "require_equal_node_counts",
62 212 : true,
63 : "Whether to require the number of nodes at each axial location to be equal");
64 106 : return params;
65 0 : }
66 :
67 55 : AverageSectionValueSampler::AverageSectionValueSampler(const InputParameters & parameters)
68 : : GeneralVectorPostprocessor(parameters),
69 55 : _mesh(_app.actionWarehouse().mesh()),
70 165 : _variables(getParam<std::vector<VariableName>>("variables")),
71 110 : _direction(getParam<Point>("axis_direction")),
72 110 : _reference_point(getParam<Point>("reference_point")),
73 220 : _positions(isParamSetByUser("lengths") ? getParam<std::vector<Real>>("lengths")
74 165 : : getParam<std::vector<Real>>("positions")),
75 110 : _have_symmetry_plane(isParamSetByUser("symmetry_plane")),
76 110 : _symmetry_plane(getParam<RealVectorValue>("symmetry_plane")),
77 55 : _automatically_locate_positions(_positions.size() == 0),
78 110 : _tolerance(getParam<Real>("tolerance")),
79 55 : _number_of_nodes(_positions.size()),
80 110 : _cross_section_maximum_radius(getParam<Real>("cross_section_maximum_radius")),
81 110 : _require_equal_node_counts(getParam<bool>("require_equal_node_counts")),
82 55 : _need_mesh_initializations(true)
83 : {
84 55 : if (parameters.isParamSetByUser("lengths") && parameters.isParamSetByUser("positions"))
85 0 : paramError("lengths", "The 'lengths' and 'positions' parameters cannot both be set.");
86 :
87 110 : if (parameters.isParamSetByUser("lengths") && _positions.size() == 0)
88 0 : paramError("lengths", "If 'lengths' is specified, at least one value must be provided");
89 :
90 110 : if (parameters.isParamSetByUser("positions") && _positions.size() == 0)
91 2 : paramError("positions", "If 'positions' is specified, at least one value must be provided");
92 :
93 53 : if (_mesh->dimension() != 3)
94 0 : mooseError("The AverageSectionValueSampler postprocessor can only be used with three "
95 : "dimensional meshes.");
96 :
97 53 : if (!MooseUtils::absoluteFuzzyEqual(_direction.norm_sq(), 1.0))
98 0 : paramError("axis_direction", "Axis direction must be a unit vector.");
99 :
100 53 : _output_vector.resize(_variables.size() + 2);
101 208 : for (const auto j : make_range(_variables.size()))
102 155 : _output_vector[j] = &declareVector(_variables[j]);
103 : const auto pos_idx = _variables.size();
104 53 : _output_vector[pos_idx] = &declareVector("axial_position");
105 53 : _output_vector[pos_idx + 1] = &declareVector("node_count");
106 :
107 206 : for (const auto j : make_range(_variables.size()))
108 : {
109 155 : if (!_sys.hasVariable(_variables[j]))
110 2 : paramError("variables",
111 : "One or more of the specified variables do not exist in this system.");
112 153 : const MooseVariable & variable = _sys.getFieldVariable<Real>(_tid, _variables[j]);
113 153 : if (!variable.isNodal())
114 0 : paramError("variables",
115 : "One or more of the specified variables are not defined at the nodes.");
116 153 : _var_numbers.push_back(variable.number());
117 : }
118 :
119 51 : if (_have_symmetry_plane)
120 : {
121 7 : const auto symm_plane_norm = _symmetry_plane.norm();
122 7 : if (MooseUtils::absoluteFuzzyEqual(symm_plane_norm, 0.0))
123 0 : mooseError("Vector defined in 'symmetry_plane' cannot have a norm of zero");
124 7 : _symmetry_plane /= _symmetry_plane.norm();
125 7 : 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 7 : mooseInfo("In AverageSectionValueSampler " + name() +
130 14 : ": \nTreating the variables as vector components (" + _variables[0] + ", " +
131 21 : _variables[1] + ", " + _variables[2] +
132 : ") and zeroing out the component normal to the prescribed 'symmetry_plane'\n");
133 : }
134 51 : }
135 :
136 : void
137 51 : AverageSectionValueSampler::initialSetup()
138 : {
139 51 : }
140 :
141 : void
142 9 : AverageSectionValueSampler::meshChanged()
143 : {
144 9 : _need_mesh_initializations = true;
145 9 : }
146 :
147 : void
148 67 : AverageSectionValueSampler::initialize()
149 : {
150 67 : if (_need_mesh_initializations)
151 : {
152 : // Option 1: locate positions first, then check counts
153 60 : if (_automatically_locate_positions)
154 39 : automaticallyLocatePositions();
155 60 : _need_mesh_initializations = false;
156 : }
157 :
158 402 : for (const auto j : make_range(_variables.size() + 2))
159 : {
160 335 : _output_vector[j]->clear();
161 335 : _output_vector[j]->resize(_positions.size());
162 : }
163 67 : _number_of_nodes.clear();
164 67 : _number_of_nodes.resize(_positions.size());
165 67 : }
166 :
167 : void
168 67 : AverageSectionValueSampler::finalize()
169 : {
170 268 : for (const auto j : make_range(_variables.size()))
171 201 : _communicator.sum(*(_output_vector[j]));
172 :
173 67 : _communicator.sum(_number_of_nodes);
174 :
175 597 : for (const auto li : index_range(_positions))
176 530 : if (_number_of_nodes[li] < 1)
177 0 : mooseError("No nodes were found in AverageSectionValueSampler postprocessor.");
178 :
179 597 : for (const auto li : index_range(_positions))
180 2120 : for (const auto j : make_range(_variables.size()))
181 1590 : (*_output_vector[j])[li] /= _number_of_nodes[li];
182 :
183 : const auto pos_idx = _variables.size();
184 67 : (*_output_vector[pos_idx]) = _positions;
185 67 : std::vector<Real> num_nodes_real(_number_of_nodes.begin(), _number_of_nodes.end());
186 67 : (*_output_vector[pos_idx + 1]) = num_nodes_real;
187 :
188 67 : if (_require_equal_node_counts)
189 : {
190 404 : for (const auto li : index_range(_number_of_nodes))
191 : {
192 354 : 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 65 : }
207 :
208 : void
209 67 : AverageSectionValueSampler::execute()
210 : {
211 134 : _block_ids = _mesh->getSubdomainIDs(getParam<std::vector<SubdomainName>>("block"));
212 :
213 67 : auto * active_nodes = _mesh->getActiveSemiLocalNodeRange();
214 :
215 : std::vector<std::vector<Real>> output_vector_partial(_variables.size(),
216 67 : std::vector<Real>(_positions.size()));
217 :
218 67 : const NumericVector<Number> & _sol = *_sys.currentSolution();
219 :
220 34751 : for (const auto & node : *active_nodes)
221 : {
222 34684 : const std::set<SubdomainID> & node_blk_ids = _mesh->getNodeBlockIds(*node);
223 :
224 279028 : for (const auto li : index_range(_positions))
225 : {
226 468636 : for (const auto i : _block_ids)
227 : {
228 244344 : if (node_blk_ids.find(i) != node_blk_ids.end())
229 : {
230 : // Check if node is within tolerance of user-prescribed plane
231 244344 : if (std::abs(axialPosition(*node) - _positions[li]) > _tolerance)
232 224292 : continue;
233 20500 : if ((*node).processor_id() != processor_id())
234 448 : continue;
235 :
236 : // Retrieve nodal variables
237 20052 : std::vector<Real> this_node_vars(_var_numbers.size());
238 80208 : for (const auto j : make_range(_var_numbers.size()))
239 60156 : this_node_vars[j] = _sol(node->dof_number(_sys.number(), _var_numbers[j], 0));
240 :
241 20052 : 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 1716 : this_node_vec_var -= _symmetry_plane * this_node_vec_var * _symmetry_plane;
246 1716 : this_node_vars[0] = this_node_vec_var(0);
247 1716 : this_node_vars[1] = this_node_vec_var(1);
248 1716 : this_node_vars[2] = this_node_vec_var(2);
249 : }
250 :
251 80208 : for (const auto j : make_range(_var_numbers.size()))
252 60156 : output_vector_partial[j][li] += this_node_vars[j];
253 :
254 20052 : _number_of_nodes[li]++;
255 : break;
256 20052 : }
257 : }
258 : }
259 : }
260 :
261 268 : for (const auto j : make_range(_variables.size()))
262 201 : *_output_vector[j] = output_vector_partial[j];
263 67 : }
264 :
265 : Real
266 260798 : AverageSectionValueSampler::axialPosition(const Node & node) const
267 : {
268 : // Compute node location w.r.t. structural component length
269 : const Point relative_distance{
270 260798 : 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 260798 : (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 260798 : if (in_plane_distance > _cross_section_maximum_radius)
279 : return std::numeric_limits<Real>::max();
280 :
281 248126 : return axial_position;
282 : }
283 :
284 : void
285 39 : AverageSectionValueSampler::automaticallyLocatePositions()
286 : {
287 39 : _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 78 : _block_ids = _mesh->getSubdomainIDs(getParam<std::vector<SubdomainName>>("block"));
297 39 : auto * nodes = _mesh->getLocalNodeRange();
298 :
299 16493 : for (const auto node : *nodes)
300 : {
301 16454 : const std::set<SubdomainID> & node_blk_ids = _mesh->getNodeBlockIds(*node);
302 :
303 32908 : for (const auto i : _block_ids)
304 : {
305 16454 : if (node_blk_ids.find(i) != node_blk_ids.end())
306 : {
307 : bool found_match = false;
308 16454 : Real axial_position = axialPosition(*node);
309 16454 : if (axial_position == std::numeric_limits<Real>::max()) // Node not in this component
310 0 : continue;
311 98724 : for (auto & pos : pos_vec)
312 : {
313 98286 : const Real dist_to_plane = std::abs(axial_position - pos);
314 98286 : if (dist_to_plane <= _tolerance)
315 : {
316 : found_match = true;
317 : break;
318 : }
319 : }
320 16454 : if (!found_match)
321 438 : pos_vec.emplace_back(axial_position);
322 : }
323 : }
324 : }
325 :
326 39 : _communicator.allgather(pos_vec);
327 :
328 701 : for (const auto & posv : pos_vec)
329 : {
330 : bool found_match = false;
331 4050 : for (auto & poss : pos_set)
332 : {
333 3612 : if (std::abs(posv - poss) < _tolerance)
334 : {
335 : found_match = true;
336 : break;
337 : }
338 : }
339 662 : if (!found_match)
340 438 : pos_set.insert(posv);
341 : }
342 477 : for (const auto & poss : pos_set)
343 438 : _positions.emplace_back(poss);
344 39 : }
|
