/
FindValueOnLine.C
166 lines (145 loc) · 4.66 KB
/
FindValueOnLine.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
#include "FindValueOnLine.h"
// MOOSE includes
#include "MooseMesh.h"
#include "MooseUtils.h"
#include "MooseVariable.h"
registerMooseObject("MooseApp", FindValueOnLine);
template <>
InputParameters
validParams<FindValueOnLine>()
{
InputParameters params = validParams<GeneralPostprocessor>();
params.addClassDescription("Find a specific target value along a sampling line. The variable "
"values along the line should change monotonically. The target value "
"is searched using a bisection algorithm.");
params.addParam<Point>("start_point", "Start point of the sampling line.");
params.addParam<Point>("end_point", "End point of the sampling line.");
params.addParam<Real>("target", "Target value to locate.");
params.addParam<unsigned int>("depth", 36, "Maximum number of bisections to perform.");
params.addParam<Real>(
"tol",
1e-10,
"Stop search if a value is found that is equal to the target with this tolerance applied.");
params.addCoupledVar("v", "Variable to inspect");
return params;
}
FindValueOnLine::FindValueOnLine(const InputParameters & parameters)
: GeneralPostprocessor(parameters),
Coupleable(this, false),
_start_point(getParam<Point>("start_point")),
_end_point(getParam<Point>("end_point")),
_length((_end_point - _start_point).norm()),
_target(getParam<Real>("target")),
_depth(getParam<unsigned int>("depth")),
_tol(getParam<Real>("tol")),
_coupled_var(*getVar("v", 0)),
_position(0.0),
_mesh(_subproblem.mesh()),
_point_vec(1)
{
}
void
FindValueOnLine::initialize()
{
// We do this here just in case it's been destroyed and recreated becaue of mesh adaptivity.
_pl = _mesh.getPointLocator();
_pl->enable_out_of_mesh_mode();
}
void
FindValueOnLine::execute()
{
Real s;
Real s_left = 0.0;
Real left = getValueAtPoint(_start_point);
Real s_right = 1.0;
Real right = getValueAtPoint(_end_point);
/**
* Here we determine the direction of the solution. i.e. the left might be the high value
* while the right might be the low value.
*/
bool left_to_right = left < right;
// Initial bounds check
if ((left_to_right && _target < left) || (!left_to_right && _target < right))
mooseError("Target value \"",
_target,
"\" is less than the minimum sampled value \"",
std::min(left, right),
"\"");
if ((left_to_right && _target > right) || (!left_to_right && _target > left))
mooseError("Target value \"",
_target,
"\" is greater than the maximum sampled value \"",
std::max(left, right),
"\"");
bool found_it = false;
Real value = 0;
for (unsigned int i = 0; i < _depth; ++i)
{
// find midpoint
s = (s_left + s_right) / 2.0;
Point p = s * (_end_point - _start_point) + _start_point;
// sample value
value = getValueAtPoint(p);
// have we hit the target value yet?
if (MooseUtils::absoluteFuzzyEqual(value, _target, _tol))
{
found_it = true;
break;
}
// bisect
if ((left_to_right && _target < value) || (!left_to_right && _target > value))
// to the left
s_right = s;
else
// to the right
s_left = s;
}
if (!found_it)
mooseError("Target value \"",
std::setprecision(10),
_target,
"\" not found on line within tolerance, last sample: ",
value,
".");
_position = s * _length;
}
Real
FindValueOnLine::getValueAtPoint(const Point & p)
{
const Elem * elem = (*_pl)(p);
processor_id_type elem_proc_id = elem ? elem->processor_id() : DofObject::invalid_processor_id;
_communicator.min(elem_proc_id);
if (elem_proc_id == DofObject::invalid_processor_id)
{
// there is no element
mooseError("No element found at the current search point. Please make sure the sampling line "
"stays inside the mesh completely.");
}
Real value = 0;
if (elem)
{
if (elem->processor_id() == processor_id())
{
// element is local
_point_vec[0] = p;
_subproblem.reinitElemPhys(elem, _point_vec, 0);
value = _coupled_var.sln()[0];
}
}
// broadcast value
_communicator.broadcast(value, elem_proc_id);
return value;
}
PostprocessorValue
FindValueOnLine::getValue()
{
return _position;
}