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 "Units.h"
11 : #include "libmesh/int_range.h"
12 :
13 : #include <stack>
14 : #include <algorithm>
15 : #include <sstream>
16 :
17 : const std::map<std::string, Real> MooseUnits::_si_prefix = {
18 : {"Y", 1e24}, {"Z", 1e21}, {"E", 1e18}, {"P", 1e15}, {"T", 1e12}, {"G", 1e9}, {"M", 1e6},
19 : {"k", 1e3}, {"h", 1e2}, {"da", 1e1}, {"d", 1e-1}, {"c", 1e-2}, {"m", 1e-3}, {"mu", 1e-6},
20 : {"n", 1e-9}, {"p", 1e-12}, {"f", 1e-15}, {"a", 1e-18}, {"z", 1e-21}, {"y", 1e-24}};
21 :
22 : const std::vector<std::pair<std::string, MooseUnits>> MooseUnits::_unit_table = {
23 : // avoid matching meters
24 : {"Ohm",
25 : {
26 : 1, // conversion factor
27 : 0, // additive shift (currently only used for Celsius and Fahrenheit)
28 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
29 : {MooseUnits::BaseUnit::METER, 2},
30 : {MooseUnits::BaseUnit::SECOND, -3},
31 : {MooseUnits::BaseUnit::AMPERE, -2}},
32 : }},
33 : {"atm",
34 : {101.325e3,
35 : 0,
36 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
37 : {MooseUnits::BaseUnit::METER, -1},
38 : {MooseUnits::BaseUnit::SECOND, -2}}}}, // Standard atmosphere
39 : // avoid matching volts
40 : {"eV",
41 : {
42 : 1.602176634e-19,
43 : 0,
44 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
45 : {MooseUnits::BaseUnit::METER, 2},
46 : {MooseUnits::BaseUnit::SECOND, -2}},
47 : }}, // electronvolt
48 : // avoid matching grams
49 : {"erg",
50 : {
51 : 1e-7,
52 : 0,
53 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
54 : {MooseUnits::BaseUnit::METER, 2},
55 : {MooseUnits::BaseUnit::SECOND, -2}},
56 : }}, // erg
57 :
58 : // avoid matching Farad or Coulomb
59 : {"degC", {1, 273.15, {{MooseUnits::BaseUnit::KELVIN, 1}}}}, // Celsius
60 : {"degF", {5.0 / 9.0, 459.67, {{MooseUnits::BaseUnit::KELVIN, 1}}}}, // Fahrenheit
61 : {"degR", {5.0 / 9.0, 0, {{MooseUnits::BaseUnit::KELVIN, 1}}}}, // Rankine
62 :
63 : {"Ang", {1e-10, 0, {{MooseUnits::BaseUnit::METER, 1}}}}, // Angstrom
64 :
65 : // Base units
66 : {"m", {1, 0, {{MooseUnits::BaseUnit::METER, 1}}}},
67 : {"g", {0.001, 0, {{MooseUnits::BaseUnit::KILOGRAM, 1}}}},
68 : {"s", {1, 0, {{MooseUnits::BaseUnit::SECOND, 1}}}},
69 : {"A", {1, 0, {{MooseUnits::BaseUnit::AMPERE, 1}}}},
70 : {"K", {1, 0, {{MooseUnits::BaseUnit::KELVIN, 1}}}},
71 : {"mol", {6.02214076e23, 0, {{MooseUnits::BaseUnit::COUNT, 1}}}},
72 : {"cd", {1, 0, {{MooseUnits::BaseUnit::CANDELA, 1}}}},
73 :
74 : // Derived units
75 : {"N",
76 : {1,
77 : 0,
78 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
79 : {MooseUnits::BaseUnit::METER, 1},
80 : {MooseUnits::BaseUnit::SECOND, -2}}}}, // Newton
81 : {"Pa",
82 : {1,
83 : 0,
84 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
85 : {MooseUnits::BaseUnit::METER, -1},
86 : {MooseUnits::BaseUnit::SECOND, -2}}}}, // Pascal
87 : {"J",
88 : {1,
89 : 0,
90 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
91 : {MooseUnits::BaseUnit::METER, 2},
92 : {MooseUnits::BaseUnit::SECOND, -2}}}}, // Joule
93 : {"cal",
94 : {4.184,
95 : 0,
96 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
97 : {MooseUnits::BaseUnit::METER, 2},
98 : {MooseUnits::BaseUnit::SECOND, -2}}}}, // Calorie
99 : {"W",
100 : {1,
101 : 0,
102 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
103 : {MooseUnits::BaseUnit::METER, 2},
104 : {MooseUnits::BaseUnit::SECOND, -3}}}}, // Watt
105 : {"C",
106 : {1, 0, {{MooseUnits::BaseUnit::AMPERE, 1}, {MooseUnits::BaseUnit::SECOND, 1}}}}, // Coulomb
107 : {"V",
108 : {1,
109 : 0,
110 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
111 : {MooseUnits::BaseUnit::METER, 2},
112 : {MooseUnits::BaseUnit::SECOND, -3},
113 : {MooseUnits::BaseUnit::AMPERE, -1}}}},
114 : {"F",
115 : {1,
116 : 0,
117 : {{MooseUnits::BaseUnit::KILOGRAM, -1},
118 : {MooseUnits::BaseUnit::METER, -2},
119 : {MooseUnits::BaseUnit::SECOND, 4},
120 : {MooseUnits::BaseUnit::AMPERE, 2}}}},
121 : {"S",
122 : {1,
123 : 0,
124 : {{MooseUnits::BaseUnit::KILOGRAM, -1},
125 : {MooseUnits::BaseUnit::METER, -2},
126 : {MooseUnits::BaseUnit::SECOND, 3},
127 : {MooseUnits::BaseUnit::AMPERE, 2}}}}, // Siemens = 1/Ohm (electrical conductance)
128 : {"Wb",
129 : {1,
130 : 0,
131 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
132 : {MooseUnits::BaseUnit::METER, 2},
133 : {MooseUnits::BaseUnit::SECOND, -2},
134 : {MooseUnits::BaseUnit::AMPERE, -1}}}}, // Weber (magnetic flux)
135 : {"T",
136 : {1,
137 : 0,
138 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
139 : {MooseUnits::BaseUnit::SECOND, -2},
140 : {MooseUnits::BaseUnit::AMPERE, -1}}}}, // Tesla (magnetic flux density)
141 : {"H",
142 : {1,
143 : 0,
144 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
145 : {MooseUnits::BaseUnit::METER, 2},
146 : {MooseUnits::BaseUnit::SECOND, -2},
147 : {MooseUnits::BaseUnit::AMPERE, -2}}}}, // Henry (inductance)
148 :
149 : // cgs units
150 : {"Ba",
151 : {0.1,
152 : 0,
153 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
154 : {MooseUnits::BaseUnit::METER, -1},
155 : {MooseUnits::BaseUnit::SECOND,
156 : -2}}}}, // barye (unit of pressure - not to be confused with bar)
157 : {"dyn",
158 : {1e-5,
159 : 0,
160 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
161 : {MooseUnits::BaseUnit::METER, 1},
162 : {MooseUnits::BaseUnit::SECOND, -2}}}}, // dyne
163 :
164 : // Freedom units
165 : {"ft", {0.3048, 0, {{MooseUnits::BaseUnit::METER, 1}}}}, // foot
166 : {"in", {25.4e-3, 0, {{MooseUnits::BaseUnit::METER, 1}}}}, // inch
167 : {"lb", {0.45359237, 0, {{MooseUnits::BaseUnit::KILOGRAM, 1}}}}, // pound of mass
168 : {"lbf",
169 : {4.4482216152605,
170 : 0,
171 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
172 : {MooseUnits::BaseUnit::METER, 1},
173 : {MooseUnits::BaseUnit::SECOND, -2}}}}, // pound of force
174 : {"psi",
175 : {6.894757e3,
176 : 0,
177 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
178 : {MooseUnits::BaseUnit::METER, -1},
179 : {MooseUnits::BaseUnit::SECOND, -2}}}}, // pound-force per square inch
180 :
181 : // misc.
182 : {"BTU",
183 : {1055.06,
184 : 0,
185 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
186 : {MooseUnits::BaseUnit::METER, 2},
187 : {MooseUnits::BaseUnit::SECOND, -2}}}}, // ISO 31-4 British thermal unit
188 : {"bar",
189 : {1e5,
190 : 0,
191 : {{MooseUnits::BaseUnit::KILOGRAM, 1},
192 : {MooseUnits::BaseUnit::METER, -1},
193 : {MooseUnits::BaseUnit::SECOND, -2}}}}, // bar (unit of pressure)
194 : {"h", {60 * 60, 0, {{MooseUnits::BaseUnit::SECOND, 1}}}}, // hour
195 : {"day", {60 * 60 * 24, 0, {{MooseUnits::BaseUnit::SECOND, 1}}}}, // day
196 : {"year",
197 : {60 * 60 * 24 * 365.25,
198 : 0,
199 : {{MooseUnits::BaseUnit::SECOND, 1}}}}, // year (equal to an annum, we did not use a for the
200 : // symbol as 1 Pa would be both 1 Pascal and 1e15 annums)
201 : {"l", {1e-3, 0, {{MooseUnits::BaseUnit::METER, 3}}}}, // liter
202 : {"u",
203 : {1.66053906660e-27, 0, {{MooseUnits::BaseUnit::KILOGRAM, 3}}}}, // unified atomic mass unit
204 : {"at", {1, 0, {{MooseUnits::BaseUnit::COUNT, 1}}}} // 1 single count (atom)
205 : };
206 :
207 204728 : MooseUnits::MooseUnits(const std::string & unit_string) : _factor(1.0), _shift(0.0), _base()
208 : {
209 : // parse the passed in unit string
210 204728 : parse(unit_string);
211 204728 : }
212 :
213 : bool
214 78 : MooseUnits::conformsTo(const MooseUnits & u) const
215 : {
216 78 : if (_base.size() != u._base.size())
217 0 : return false;
218 247 : for (std::size_t i = 0; i < _base.size(); ++i)
219 169 : if (_base[i] != u._base[i])
220 0 : return false;
221 :
222 78 : return true;
223 : }
224 :
225 : Real
226 48 : MooseUnits::convert(Real from_value, const MooseUnits & from_unit) const
227 : {
228 48 : if (!conformsTo(from_unit))
229 0 : mooseError("Cannot convert between non-conforming units '", *this, "' and '", from_unit, "'.");
230 48 : return ((from_value + from_unit._shift) * from_unit._factor) / _factor - _shift;
231 : }
232 :
233 : void
234 204728 : MooseUnits::parse(const std::string & unit_string)
235 : {
236 : // we need to understand * / ^int ( ) string
237 : // m^2*kg/(N*s)^2
238 : // no numerical expressions are permitted (e.g. no 1/2 0.5 4*2)
239 :
240 : // parse stack
241 204728 : std::stack<std::vector<std::pair<MooseUnits::BaseUnit, int>>> stack;
242 204728 : std::stack<char> op_stack;
243 204728 : std::stack<MooseUnits> out;
244 :
245 204728 : auto it = unit_string.begin();
246 204728 : const auto end = unit_string.end();
247 :
248 : do
249 : {
250 : // skip whitespace
251 611785 : if (*it == ' ')
252 : {
253 2 : it++;
254 2 : continue;
255 : }
256 : // opening parenthesis
257 611783 : else if (*it == '(')
258 7 : op_stack.push(*(it++));
259 : // multiply or divide
260 611776 : else if (*it == '*' || *it == '/')
261 : {
262 : // pop */ off
263 203525 : while (!op_stack.empty() && (op_stack.top() == '*' || op_stack.top() == '/'))
264 : {
265 12 : auto top = op_stack.top();
266 12 : op_stack.pop();
267 :
268 12 : if (out.size() < 2)
269 0 : parseError(
270 : unit_string, it, "Applying ", top, " but don't have enough items on the stack.");
271 12 : auto rhs = out.top();
272 12 : out.pop();
273 12 : if (top == '*')
274 6 : out.top() = out.top() * rhs;
275 : else
276 6 : out.top() = out.top() / rhs;
277 12 : }
278 :
279 203513 : op_stack.push(*(it++));
280 : }
281 : // closing parenthesis
282 408263 : else if (*it == ')')
283 : {
284 : do
285 : {
286 13 : if (op_stack.empty())
287 0 : parseError(unit_string, it, "Mismatching right parenthesis.");
288 :
289 13 : auto top = op_stack.top();
290 13 : op_stack.pop();
291 :
292 13 : if (top == '(')
293 7 : break;
294 :
295 6 : if (top == '*' || top == '/')
296 : {
297 6 : if (out.size() < 2)
298 0 : parseError(
299 : unit_string, it, "Applying ", top, " but don't have enough items on the stack.");
300 6 : auto rhs = out.top();
301 6 : out.pop();
302 6 : if (top == '*')
303 4 : out.top() = out.top() * rhs;
304 : else
305 2 : out.top() = out.top() / rhs;
306 6 : }
307 6 : } while (true);
308 7 : it++;
309 : }
310 : // exponent
311 408256 : else if (*it == '^')
312 : {
313 : // check for sign
314 15 : int sign = 1;
315 15 : ++it;
316 15 : if (it != end && *it == '-')
317 : {
318 0 : sign = -1;
319 0 : it++;
320 : }
321 :
322 15 : int num = 0;
323 15 : bool valid = false;
324 30 : while (*it >= '0' && *it <= '9')
325 : {
326 15 : num *= 10;
327 15 : num += *(it++) - '0';
328 15 : valid = true;
329 : }
330 :
331 : // error
332 15 : if (!valid)
333 : {
334 0 : if (it == end)
335 0 : parseError(unit_string, it, "Expected a number but found end of string.");
336 : else
337 0 : parseError(unit_string, it, "Expected a number but got '", (*it), "'.");
338 : }
339 :
340 : // exponents do not get pushed on the stack, but are applied immediately
341 15 : if (out.empty())
342 0 : parseError(unit_string, it, "Applying an exponent, but found no units in front of it.");
343 15 : out.top() = std::pow(out.top(), sign * num);
344 : }
345 : else
346 : {
347 408241 : std::string unit;
348 1227424 : while (it != end && ((*it >= 'a' && *it <= 'z') || (*it >= 'A' && *it <= 'Z') ||
349 407817 : (*it >= '0' && *it <= '9') || *it == '.' || *it == '-'))
350 411366 : unit += *(it++);
351 :
352 : // check if a number prefix is present
353 : Real si_factor;
354 408241 : std::stringstream ss(unit);
355 408241 : if ((ss >> si_factor).fail() || !std::isfinite(si_factor))
356 204789 : si_factor = 1.0;
357 : else
358 : {
359 203452 : if (ss.eof())
360 : {
361 203443 : out.push(MooseUnits(si_factor));
362 203443 : continue;
363 : }
364 9 : ss >> unit;
365 : }
366 :
367 : // parse the unit and a potential SI prefix
368 204798 : auto len = unit.length();
369 204798 : if (len == 0)
370 : {
371 0 : if (it == end)
372 0 : parseError(unit_string, it, "Expected unit but found end of string.");
373 : else
374 0 : parseError(unit_string, it, "Expected unit but found '", *it, "'.");
375 : }
376 :
377 204798 : unsigned int i = 0;
378 1857354 : for (; i < _unit_table.size(); ++i)
379 : {
380 1857354 : const auto & s = _unit_table[i].first;
381 1857354 : const auto slen = s.length();
382 1857354 : if (slen <= len && unit.substr(len - slen) == s)
383 204798 : break;
384 : }
385 :
386 : // invalid unit
387 204798 : if (i == _unit_table.size())
388 0 : parseError(unit_string, it, "Unknown unit '", unit, "'.");
389 :
390 : // get prefix
391 204798 : const auto & s = _unit_table[i].first;
392 204798 : const auto slen = s.length();
393 204798 : if (slen < len)
394 : {
395 918 : const auto prefix = unit.substr(0, len - slen);
396 918 : auto jt = _si_prefix.find(prefix);
397 918 : if (jt != _si_prefix.end())
398 918 : si_factor = jt->second;
399 : else
400 0 : parseError(unit_string, it, "Unknown SI prefix '", prefix, "' on unit '", unit, "'.");
401 918 : }
402 :
403 : // construct a unit object and push to the output stack
404 204798 : MooseUnits d = _unit_table[i].second * si_factor;
405 204798 : out.push(d);
406 611684 : }
407 611785 : } while (it != end);
408 :
409 : // unwind operator stack
410 408223 : while (!op_stack.empty())
411 : {
412 203495 : auto top = op_stack.top();
413 203495 : op_stack.pop();
414 :
415 203495 : if (top == '(' || top == ')')
416 0 : parseError(unit_string, it, "Unit string contains unmatched parenthesis.");
417 :
418 203495 : if (top == '*' || top == '/')
419 : {
420 203495 : if (out.size() < 2)
421 0 : parseError(unit_string,
422 : it,
423 : "Applying ",
424 : top,
425 : " but don't have enough items on the stack.",
426 : out.size());
427 203495 : auto rhs = out.top();
428 203495 : out.pop();
429 203495 : if (top == '*')
430 203400 : out.top() = out.top() * rhs;
431 : else
432 95 : out.top() = out.top() / rhs;
433 203495 : }
434 : }
435 :
436 : // consistency check
437 204728 : if (out.size() != 1)
438 0 : mooseError("Internal parse error. Remaining output stack size should be 1, but is ",
439 0 : out.size());
440 :
441 204728 : *this = out.top();
442 204728 : simplify();
443 204728 : }
444 :
445 : MooseUnits
446 204798 : MooseUnits::operator*(const Real f) const
447 : {
448 204798 : MooseUnits u = *this;
449 204798 : u._factor *= f;
450 204798 : u._shift *= f;
451 204798 : return u;
452 : }
453 :
454 : MooseUnits
455 203452 : MooseUnits::operator*(const MooseUnits & rhs) const
456 : {
457 203452 : MooseUnits u = rhs;
458 203452 : u._factor *= _factor;
459 203452 : u._base.insert(u._base.end(), _base.begin(), _base.end());
460 203452 : u._shift = 0.0;
461 203452 : u.simplify();
462 203452 : return u;
463 0 : }
464 :
465 : MooseUnits
466 1174 : MooseUnits::operator/(const MooseUnits & rhs) const
467 : {
468 1174 : MooseUnits u = rhs;
469 1174 : u._factor = _factor / u._factor;
470 3427 : for (auto & b : u._base)
471 2253 : b.second *= -1;
472 1174 : u._base.insert(u._base.end(), _base.begin(), _base.end());
473 1174 : u._shift = 0.0;
474 1174 : u.simplify();
475 1174 : return u;
476 0 : }
477 :
478 : bool
479 1 : MooseUnits::operator==(const MooseUnits & rhs) const
480 : {
481 1 : if (_factor != rhs._factor || _shift != rhs._shift)
482 0 : return false;
483 1 : return conformsTo(rhs);
484 : }
485 :
486 : bool
487 0 : MooseUnits::operator==(const Real rhs) const
488 : {
489 0 : return (_factor == rhs && _shift == 0.0 && _base.empty());
490 : }
491 :
492 1115 : MooseUnits::operator Real() const
493 : {
494 1115 : if (!_base.empty())
495 0 : mooseError("Unit '", *this, "' is not a pure number.");
496 1115 : return _factor;
497 : }
498 :
499 : void
500 409354 : MooseUnits::simplify()
501 : {
502 409354 : std::map<BaseUnit, int> base_map;
503 823586 : for (auto & b : _base)
504 : {
505 414232 : auto j = base_map.find(b.first);
506 414232 : if (j == base_map.end())
507 412010 : base_map.insert(b);
508 : else
509 2222 : j->second += b.second;
510 : }
511 :
512 : // replace base vector with map contents
513 409354 : _base.clear();
514 821364 : for (auto & u : base_map)
515 412010 : if (u.second != 0)
516 409799 : _base.push_back(u);
517 409354 : }
518 :
519 : bool
520 140 : MooseUnits::isBase(const MooseUnits::BaseUnit base) const
521 : {
522 : // must have only one base unit
523 140 : if (_base.size() != 1)
524 0 : return false;
525 :
526 : // its exponent must be one
527 140 : if (_base[0].second != 1)
528 0 : return false;
529 :
530 : // and it has to be the right base unit
531 140 : return (_base[0].first == base);
532 : }
533 :
534 : namespace std
535 : {
536 : MooseUnits
537 15 : pow(const MooseUnits & u, int e)
538 : {
539 15 : MooseUnits r = u;
540 15 : r._factor = std::pow(u._factor, e);
541 15 : r._shift = 0.0;
542 29 : for (auto & b : r._base)
543 14 : b.second *= e;
544 15 : return r;
545 : }
546 : }
547 :
548 : std::ostream &
549 0 : MooseUnits::latex(std::ostream & os)
550 : {
551 0 : os.iword(geti()) = 1;
552 0 : return os;
553 : }
554 : std::ostream &
555 0 : MooseUnits::text(std::ostream & os)
556 : {
557 0 : os.iword(geti()) = 0;
558 0 : return os;
559 : }
560 :
561 : int
562 18 : MooseUnits::geti()
563 : {
564 18 : static int i = std::ios_base::xalloc();
565 18 : return i;
566 : }
567 :
568 : std::ostream &
569 18 : operator<<(std::ostream & os, const MooseUnits & u)
570 : {
571 18 : bool latex = (os.iword(MooseUnits::geti()) == 1);
572 26 : static const std::vector<std::string> base_unit = {"m", "kg", "s", "A", "K", "at", "cd"};
573 :
574 18 : if (u._factor != 1.0 || u._base.empty())
575 : {
576 14 : os << u._factor;
577 14 : if (!u._base.empty())
578 11 : os << ' ';
579 : }
580 :
581 57 : for (const auto i : index_range(u._base))
582 : {
583 39 : os << (i ? (latex ? "\\cdot " : "*") : "");
584 :
585 39 : const auto & unit = base_unit[int(u._base[i].first)];
586 39 : if (latex)
587 0 : os << "\\text{" << unit << "}";
588 : else
589 39 : os << unit;
590 :
591 39 : if (u._base[i].second != 1)
592 : {
593 25 : if (latex)
594 0 : os << "^{" << u._base[i].second << '}';
595 : else
596 25 : os << '^' << u._base[i].second;
597 : }
598 : }
599 :
600 18 : return os;
601 2 : }
|
