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 "FindValueOnLine.h"
11 :
12 : // MOOSE includes
13 : #include "MooseMesh.h"
14 : #include "MooseUtils.h"
15 : #include "MooseVariable.h"
16 :
17 : registerMooseObject("MooseApp", FindValueOnLine);
18 :
19 : InputParameters
20 14417 : FindValueOnLine::validParams()
21 : {
22 14417 : InputParameters params = GeneralPostprocessor::validParams();
23 14417 : params.addClassDescription("Find a specific target value along a sampling line. The variable "
24 : "values along the line should change monotonically. The target value "
25 : "is searched using a bisection algorithm.");
26 14417 : params.addParam<Point>("start_point", "Start point of the sampling line.");
27 14417 : params.addParam<Point>("end_point", "End point of the sampling line.");
28 14417 : params.addParam<Real>("target", "Target value to locate.");
29 43251 : params.addParam<bool>(
30 : "error_if_not_found",
31 28834 : true,
32 : "If true, stop with error if target value is not found on the line. If false, "
33 : "return default_value.");
34 43251 : params.addParam<Real>("default_value",
35 28834 : -1,
36 : "Value to return if target value is not found on line and "
37 : "error_if_not_found is false.");
38 14417 : params.addParam<unsigned int>("depth", 36, "Maximum number of bisections to perform.");
39 43251 : params.addParam<Real>(
40 : "tol",
41 28834 : 1e-10,
42 : "Stop search if a value is found that is equal to the target with this tolerance applied.");
43 14417 : params.addCoupledVar("v", "Variable to inspect");
44 14417 : return params;
45 0 : }
46 :
47 76 : FindValueOnLine::FindValueOnLine(const InputParameters & parameters)
48 : : GeneralPostprocessor(parameters),
49 : Coupleable(this, false),
50 76 : _start_point(getParam<Point>("start_point")),
51 76 : _end_point(getParam<Point>("end_point")),
52 76 : _length((_end_point - _start_point).norm()),
53 76 : _target(getParam<Real>("target")),
54 76 : _error_if_not_found(getParam<bool>("error_if_not_found")),
55 76 : _default_value(getParam<Real>("default_value")),
56 76 : _depth(getParam<unsigned int>("depth")),
57 76 : _tol(getParam<Real>("tol")),
58 76 : _coupled_var(*getVar("v", 0)),
59 76 : _position(0.0),
60 76 : _mesh(_subproblem.mesh()),
61 152 : _point_vec(1)
62 : {
63 76 : }
64 :
65 : void
66 346 : FindValueOnLine::initialize()
67 : {
68 : // We do this here just in case it's been destroyed and recreated becaue of mesh adaptivity.
69 346 : _pl = _mesh.getPointLocator();
70 346 : _pl->enable_out_of_mesh_mode();
71 346 : }
72 :
73 : void
74 346 : FindValueOnLine::execute()
75 : {
76 : Real s;
77 346 : Real s_left = 0.0;
78 346 : Real left = getValueAtPoint(_start_point);
79 346 : Real s_right = 1.0;
80 346 : Real right = getValueAtPoint(_end_point);
81 :
82 : // Helper for erroring; won't error if error_if_not_found=false
83 332 : const auto error = [this](auto &&... message)
84 : {
85 166 : if (_error_if_not_found)
86 12 : mooseError(std::scientific, std::setprecision(4), message...);
87 154 : _position = _default_value;
88 496 : };
89 :
90 : /**
91 : * Here we determine the direction of the solution. i.e. the left might be the high value
92 : * while the right might be the low value.
93 : */
94 342 : bool left_to_right = left < right;
95 : // Initial bounds check
96 342 : if ((left_to_right && _target < left) || (!left_to_right && _target < right))
97 : {
98 136 : error("Target value ",
99 136 : _target,
100 : " is less than the minimum sampled value ",
101 : std::min(left, right));
102 154 : return;
103 : }
104 206 : if ((left_to_right && _target > right) || (!left_to_right && _target > left))
105 : {
106 26 : error("Target value ",
107 26 : _target,
108 : " is greater than the maximum sampled value ",
109 : std::max(left, right));
110 22 : return;
111 : }
112 :
113 180 : bool found_it = false;
114 180 : Real value = 0;
115 5982 : for (unsigned int i = 0; i < _depth; ++i)
116 : {
117 : // find midpoint
118 5978 : s = (s_left + s_right) / 2.0;
119 5978 : Point p = s * (_end_point - _start_point) + _start_point;
120 :
121 : // sample value
122 5978 : value = getValueAtPoint(p);
123 :
124 : // have we hit the target value yet?
125 5978 : if (MooseUtils::absoluteFuzzyEqual(value, _target, _tol))
126 : {
127 176 : found_it = true;
128 176 : break;
129 : }
130 :
131 : // bisect
132 5802 : if ((left_to_right && _target < value) || (!left_to_right && _target > value))
133 : // to the left
134 2707 : s_right = s;
135 : else
136 : // to the right
137 3095 : s_left = s;
138 : }
139 :
140 : // Return error if target value (within tol) was not found within depth bisections
141 180 : if (!found_it)
142 : {
143 4 : error("Target value ",
144 4 : _target,
145 : " not found on line within tolerance.\nLast sample: ",
146 : value,
147 : ", difference from target: ",
148 4 : std::abs(_target - value));
149 0 : return;
150 : }
151 :
152 176 : _position = s * _length;
153 : }
154 :
155 : Real
156 6670 : FindValueOnLine::getValueAtPoint(const Point & p)
157 : {
158 6670 : const Elem * elem = (*_pl)(p);
159 :
160 : processor_id_type elem_proc_id =
161 6670 : elem ? elem->processor_id() : libMesh::DofObject::invalid_processor_id;
162 6670 : _communicator.min(elem_proc_id);
163 :
164 6670 : if (elem_proc_id == libMesh::DofObject::invalid_processor_id)
165 : {
166 : // there is no element
167 4 : mooseError("No element found at the current search point. Please make sure the sampling line "
168 : "stays inside the mesh completely.");
169 : }
170 :
171 6666 : Real value = 0;
172 :
173 6666 : if (elem)
174 : {
175 6115 : if (elem->processor_id() == processor_id())
176 : {
177 : // element is local
178 4872 : _point_vec[0] = p;
179 4872 : _subproblem.reinitElemPhys(elem, _point_vec, 0);
180 4872 : value = _coupled_var.sln()[0];
181 : }
182 : }
183 :
184 : // broadcast value
185 6666 : _communicator.broadcast(value, elem_proc_id);
186 6666 : return value;
187 : }
188 :
189 : PostprocessorValue
190 330 : FindValueOnLine::getValue() const
191 : {
192 330 : return _position;
193 : }
|
