# coding: utf-8
"""
Defines the DEQATN class and sub-functions.
The capitalization of the sub-functions is important.
"""
from __future__ import annotations
import re
from typing import Union, Any, TYPE_CHECKING
import numpy as np
from numpy import (
cos, sin, tan, log, log10, mean, exp, sqrt, square, mod, abs, sum,
arcsin as asin, arccos as acos, arctan as atan, arctan2 as atan2,
arcsinh as asinh, arccosh as acosh, arctanh as atanh)
# atan2h
from numpy.linalg import norm # type: ignore
from pyNastran.bdf.cards.base_card import BaseCard
if TYPE_CHECKING: # pragma: no cover
from pyNastran.bdf.bdf import BDF, DTABLE
[docs]
def pi(num):
"""weird way to multiply π by a number"""
return np.pi * num
[docs]
def avg(*args):
"""average"""
return np.mean(args)
[docs]
def ssq(*args):
"""sum of squares"""
return np.square(args).sum()
[docs]
def logx(x, y):
"""log base_x(y)"""
return log(y**x) / log(x)
[docs]
def dim(x, y):
"""positive difference"""
return x - min(x, y)
[docs]
def db(p, pref):
"""sound pressure in decibels"""
return 20. * log(p / pref)
#def _Log(z):
[docs]
def atan2h(x, y):
"""
Hyperbolic arctangent
>>> arctanh(z) = 1/2 Log((1+z)/(1-z))
real:
the following must be true: |x1| > |x2| and x2 ≠ 0.
complex:
x1 = a + bi
x2 = b + di
a = b = 0 and (sign of c) = (sign of d):
the result is 0.
a = b = 0 and (sign of c) ≠ (sign of d):
the result is π.
c = d = 0 and (sign of a) = (sign of b):
the result is π/2.
c = d = 0 and (sign of a) ≠ (sign of b):
the result is −π/2
"""
#integer_float_types = (int, np.int32, float, np.float32)
#if isinstance(x, integer_float_types):
#assert x >= 0, 'x=%s y=%s' % (x, y)
#assert y > 0, 'x=%s y=%s' % (x, y)
#return np.arctanh(x, y)
raise NotImplementedError()
[docs]
def invdb(dbi: float, pref: float) -> float:
"""inverse Db"""
return 10. ** (dbi / 20. + log(pref))
[docs]
def dba(p: float, pref: float, f: float) -> float:
"""
sound pressure in decibels (perceived)
Parameters
----------
p : float
structural responses or acoustic pressure
f : float
forcing frequency
pref : float
reference pressure
Returns
-------
dbi : float
acoustic pressure in Decibels
"""
ta1, ta2 = _get_ta(f)
return 20. * log(p / pref) + 10 * log(ta1) + 10. * log(ta2)
[docs]
def invdba(dbai: float, pref: float, f: float) -> float:
"""
Inverse Dba
Parameters
----------
dbai : float
acoustic pressure in Decibels (perceived)
f : float
forcing frequency
pref : float
reference pressure
Returns
-------
p : float
structural responses or acoustic pressure
"""
ta1, ta2 = _get_ta(f)
#dbai = dba(p, pref, f)
return 10. ** ((dbai - 10. * log(ta1) - 10. * log(ta2))/20)
[docs]
def _get_ta(f: float) -> float:
"""gets the factors for dba, invdba"""
k1 = 2.242882e16
k3 = 1.562339
p1 = 20.598997
p2 = 107.65265
p3 = 737.86223
p4 = 12194.22
ta1 = k3 * f**4 / ((f**2 + p2**2) * (f**2 + p3**2))
ta2 = k1 * f**4 / ((f**2 + p1**2)**2 * (f**2 + p4**2)**2)
return ta1, ta2
# we'll add _ to the beginning of these variables
BUILTINS = ['del', 'eval', 'yield', 'async', 'await', 'property',
'slice', 'filter', 'map']
[docs]
class DEQATN(BaseCard): # needs work...
"""
Design Equation Definition
Defines one or more equations for use in design sensitivity analysis.
+--------+------+-----+-----+-----+-----+-------+-----+
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
+========+======+=====+=====+=====+=====+=======+=====+
| DEQATN | EQID | EQUATION |
+--------+------+-------------------------------------+
| | EQUATION (cont.) |
+--------+--------------------------------------------+
"""
type = 'DEQATN'
_properties = ['dtable']
def __init__(self, equation_id: int, eqs: list[str], comment=''):
"""
Creates a DEQATN card
Parameters
----------
equation_id : int
the id of the equation
eqs : list[str]
the equations, which may overbound the field
split them by a semicolon (;)
comment : str; default=''
a comment for the card
DEQATN 41 F1(A,B,C,D,R) = A+B *C–(D**3 + 10.0) + sin(PI(1) * R)
+ A**2 / (B - C); F = A + B - F1 * D
def F1(A, B, C, D, R):
F1 = A+B *C-(D**3 + 10.0) + sin(PI(1) * R) + A**2 / (B – C)
F = A + B - F1 * D
return F
eqs = [
'F1(A,B,C,D,R) = A+B *C–(D**3 + 10.0) + sin(PI(1) * R) + A**2 / (B – C)',
'F = A + B – F1 * D',
]
>>> deqatn = DEQATN(41, eq, comment='')
"""
if comment:
self.comment = comment
self.dtable = None
self.dtable_ref = None
self.equation_id = equation_id
self.eqs = eqs
self.func = None
self.func_str = ''
self.func_name = ''
self.nargs = 0
[docs]
@classmethod
def _init_from_empty(cls):
equation_id = 1
eqs = []
return DEQATN(equation_id, eqs, comment='')
[docs]
@classmethod
def add_card(cls, card: list[str], comment: str=''):
"""
Adds a DEQATN card from ``BDF.add_card(...)``
Parameters
----------
card : list[str]
this card is special and is not a ``BDFCard`` like other cards
comment : str; default=''
a comment for the card
"""
#print(card)
line0 = card[0]
if '\t' in line0:
line0 = line0.expandtabs()
name_eqid = line0[:16]
#print('name_eqid = %r' % name_eqid)
assert ',' not in name_eqid, name_eqid
try:
name, eq_id = name_eqid.split()
assert name.strip().upper() == 'DEQATN', card
except ValueError:
msg = 'cannot split %r\n' % name_eqid
msg += "Expected data of the form 'DEQATN 100'\n"
msg += 'card=%s' % card
raise ValueError(msg)
equation_id = int(eq_id)
# combine the equations into a single organized block
line0_eq = line0[16:]
eqs_temp = [line0_eq] + card[1:]
#eqs_temp2 = [line.replace(';;', ';') for line in eqs_temp]
#for line in eqs_temp2:
#print(line)
eqs = lines_to_eqs(eqs_temp)
return DEQATN(equation_id, eqs, comment=comment)
[docs]
def _setup_equation(self) -> None:
"""
creates an executable equation object from self.eqs
x = 10.
>>> deqatn.func(x)
42.0
>>> deqatn.func_str
def stress(x):
x = float(x)
return x + 32.
"""
#default_values = {}
#if self.dtable is not None:
#default_values = self.dtable_ref.default_values
#func_name, nargs, func_str = fortran_to_python(
#self.equation_id, self.eqs, default_values, str(self))
func_name, nargs, func_str = _setup_deqatn(
self.equation_id, self.eqs,
self.dtable_ref,
str(self),
)
self.func_str = func_str
self.func_name = func_name
try:
exec(func_str)
except SyntaxError:
print(func_str)
raise
#print(locals().keys())
func = locals()[func_name]
setattr(self, func_name, func)
#print(func)
self.func = func
self.nargs = nargs
[docs]
def cross_reference(self, model: BDF) -> None:
"""
Cross links the card so referenced cards can be extracted directly
Parameters
----------
model : BDF()
the BDF object
"""
# TODO: get defaults from DTABLE
# TODO: get limits from DCONSTR
self.dtable = model.dtable
self.dtable_ref = self.dtable
self._setup_equation()
[docs]
def safe_cross_reference(self, model: BDF) -> None:
self.cross_reference(model)
[docs]
def uncross_reference(self) -> None:
"""Removes cross-reference links"""
del self.func
#del self.f
#del getattr(self, self.func_name)
setattr(self, self.func_name, None)
del self.func_name
del self.nargs
del self.dtable, self.dtable_ref
def _verify(self, xref: bool) -> None:
pass
[docs]
def evaluate(self, *args) -> float:
"""Makes a call to self.func"""
#args2 = args[:self.nargs]
#print('args =', args2)
if len(args) > self.nargs:
msg = 'len(args) > nargs\n'
msg += 'nargs=%s len(args)=%s; func_name=%s' % (
self.nargs, len(args), self.func_name)
raise RuntimeError(msg)
return self.func(*args)
#self.func(*args)
[docs]
def raw_fields(self) -> list[str]:
return [self.write_card()]
[docs]
def repr_fields(self) -> list[str]:
return self.raw_fields()
[docs]
def write_card(self, size: int=8, is_double: bool=False) -> str:
msg = write_deqatn(self.equation_id, self.eqs,
size=size, is_double=is_double)
return msg
[docs]
def write_deqatn(equation_id: int, eqs: list[str],
size: int=8, is_double: bool=False) -> str:
"""directly write a DEQATN"""
#self.evaluate(1, 2)
split_eqs = split_equations(eqs)
equation_line0 = split_eqs[0]
#assert len(equation_line0) <= 56, equation_line0
msg = 'DEQATN %-8i%-56s\n' % (equation_id, equation_line0)
assert len(equation_line0) <= 56, equation_line0
for eq in split_eqs[1:]:
msg += ' %-64s\n' % eq
assert len(eq) <= 64, eq
#print(msg)
return msg
[docs]
def lines_to_eqs(eqs_in: list[str]) -> list[str]:
"""splits the equations"""
eqs_wrapped = _split_equations_by_semicolon(eqs_in)
eqs = _join_wrapped_equation_lines(eqs_in, eqs_wrapped)
assert len(eqs) > 0, eqs
return eqs
[docs]
def _split_equations_by_semicolon(eqs_in: list[str]) -> list[str]:
"""helper for ``lines_to_eqs``"""
eqs_temp_out = []
nchars = 72 - 16
for iline, eq in enumerate(eqs_in):
if iline == 0:
eq2 = eq[:nchars].strip(' \t\n')
else:
eq2 = eq.expandtabs()[8:nchars].strip(' \t\n')
semicolon = ';' if eq2.rstrip().endswith(';') else ''
eq2 = eq2.rstrip(' \t;')
#nline = len(eq.rstrip('; \n')) + 16
#print('eq2=%r' % eq2)
if ';' in eq2:
eq2s = eq2.split(';')
eq_tempi = [eqi.strip() + ';' for eqi in eq2s if eqi.strip()]
#for check_line in eq2s:
#print(check_line)
#_check_for_valid_line(check_line, eq)
#print('eq_tempi = %r' % eq_tempi)
#eq_tempi[-1] += semicolon
eqs_temp_out += eq_tempi
else:
check_line = eq2 + semicolon
#print('check_line = %r' % (check_line))
#_check_for_valid_line(check_line, eq)
eqs_temp_out.append(check_line)
nchars = 72
return eqs_temp_out
#def _check_for_valid_line(check_line, full_line):
#if '=' not in check_line:
#msg = (
#'expected an equal sign (the first 8 characters are removed)\n'
#'line =%r\n'
#'full_line=%r' % (check_line, full_line))
#raise SyntaxError(msg)
[docs]
def _join_wrapped_equation_lines(unused_eqs_temp_in, eqs_temp: list[str]) -> list[str]:
"""helper for ``lines_to_eqs``"""
eqs = []
neqs = len(eqs_temp)
is_join = False
eqi = ''
for i, eq in enumerate(eqs_temp):
#print(f'i={i} join={is_join} eq={eq!r}')
if is_join:
eq = eqi.rstrip() + eq.lstrip()
eqi = eq.strip().replace(' ', '')
if i == 0 and eqi == '':
#self.eqs.append(eqi)
continue
if i == 0:
# first line
if eqi.endswith(';'):
eqi = eqi[:-1]
assert not eqi.endswith(';'), eq
else:
is_join = True
assert len(eqi) <= 56, eqi
elif i != neqs-1:
# mid line
#assert len(eqi) <= 64, 'len(eqi)=%s eq=%r' % (len(eqi), eqi)
if eqi.endswith(';'):
eqi = eqi[:-1]
is_join = False
assert not eqi.endswith(';'), eq
else:
is_join = True
else:
# last line
is_join = False
if not is_join:
if '=' not in eqi:
raise SyntaxError('line=%r expected an equal sign' % eqi)
eqs.append(eqi)
#print(i, eqi)
#assert not is_join
if is_join:
eqs.append(eqi)
return eqs
[docs]
def split_equations(lines: list[str]) -> list[str]:
"""takes an overbounded DEQATN card and shortens it"""
# first line must be < 56
# second line may be < 64
lines2 = []
for i, line in enumerate(lines):
#print('-------------------------')
# we'll add ; to the end of each line
if i == 0:
lines2 += _split_equation([], line.strip() + ';', 56)
else:
lines2 += _split_equation([], line.strip() + ';', 64)
# remove the trailing semicolon
lines2[-1] = lines2[-1][:-1]
return lines2
[docs]
def _split_equation(lines_out: list[str], line: str, n: int,
isplit: int=0) -> list[str]:
"""
Takes an overbounded DEQATN line and shortens it using recursion
Parameters
----------
lines_out : list[str]
len(lines) = 0 : first iteration
len(lines) = 1 : second iteration
line : str
the line to split
n : int
the maximum number of characters allowed
the first line of the DEQATN has a different number of fields
allowed vs. subsequent lines
isplit : int; default=0
the number of levels deep in the recursive function we are
Returns
-------
lines_out : list[str]
the long line broken into shorter lines
"""
#print('n=%s -> line=%r len=%s' % (n, line, len(line)))
if len(line) <= n:
lines_out.append(line.strip())
return lines_out
# equation must be split
line0 = line[:n][::-1].replace('**', '^')
# fore, aft = line0.split('+-()*', 1)
#print('line0 = %r; len=%s' % (str(line0[::-1]), len(line0)))
out = {}
for operator in ('+', '*', '^', '-', ')', ',', '='):
if operator in line0:
i = line0.index(operator)
out[i] = operator
try:
imin = min(out)
except ValueError:
msg = "Couldn't find an operator ()+-/*= in %r\n" % line[n:]
msg += 'line = %r' % line
raise ValueError(msg)
operator = out[imin]
#print('operator = %r' % operator)
unused_fore, aft = line0.split(operator, 1)
i = len(aft) + 1
line_out = line[:i]
#print('appending %r; len=%s' % (line_out, len(line_out)))
#print('fore = %r' % fore[::-1])
#print('aft = %r' % aft[::-1])
lines_out.append(line_out.replace('^', '**').strip())
isplit += 1
if isplit > 360:
raise RuntimeError('Recursion issue splitting line; isplit=%i' % isplit)
lines_out = _split_equation(lines_out, line[i:], n, isplit+1)
return lines_out
[docs]
def fortran_to_python_short(line: str, unused_default_values: Any) -> Any:
"""the function used by the DRESP2"""
func_str = 'def func(args):\n'
func_str += ' return %s(args)\n' % line.strip()
local_dict = {}
exec(func_str, globals(), local_dict)
return local_dict['func']
[docs]
def split_to_equations(lines: list[str]) -> list[str]:
"""
Splits a line like::
b = a + z; c = 42
into::
b = a + z
c = 42
"""
equation_lines = []
for line in lines:
line = line.rstrip(' ;')
if ';' in line:
lines2 = line.split(';')
equation_lines.extend(lines2)
else:
equation_lines.append(line)
return equation_lines
[docs]
def _setup_deqatn(equation_id: int,
eqs: list[str],
dtable_ref: Optional[DTABLE],
comment: str):
default_values = {}
if dtable_ref is not None:
default_values = dtable_ref.default_values
func_name, nargs, func_str = fortran_to_python(
equation_id, eqs, default_values, comment)
return func_name, nargs, func_str
[docs]
def fortran_to_python(deqatn_id: int,
lines: list[str],
default_values: dict[str, Union[float, np.ndarray]],
comment: str='') -> tuple[str, int, str]:
"""
Creates the python function
Parameters
----------
lines : list[str]
the equations to write broken up by statement
default_values : dict[name] = value
the default values from the DTABLE card
Returns
-------
func_name : str
the name of the function
nargs : int
the number of variables to the function
func_msg : str
the python function
def f(x, y=10.):
'''
$ deqatn
DEQATN 1000 f(x,y) = x+y
'''
try:
if isinstance(x, (int, float, str)):
x = float(x)
if isinstance(y, (int, float, str)):
y = float(y)
except Exception:
print(locals())
raise
f = x + y
return f
"""
func_msg = ''
variables = []
assert len(lines) > 0, lines
equation_lines = split_to_equations(lines)
for i, line in enumerate(equation_lines):
#print('--------------------')
line = line.lower()
#func_msg += '#i=%s\n' % i
assert ';' not in line, line
try:
# f(x, y) = 10.
# f(x, y) = abs(x) + y
# f = 42.
f, eq = line.split('=')
except ValueError:
if '=' not in line:
raise SyntaxError('= not found in %r' % (line))
msg = 'only 1 = sign may be found a line\n'
msg += 'line = %r\n' % line
if len(lines) > 1:
msg += 'lines:\n%s' % '\n'.join(lines)
raise SyntaxError(msg)
f = f.strip()
eq = eq.strip().rstrip(';')
#print('f=%r eq=%r' % (f, eq))
for builtin in BUILTINS:
if builtin == f:
f = f.replace(builtin, builtin + '_')
if builtin == eq:
eq = eq.replace(builtin, builtin + '_')
if i == 0:
func_name, func_msg, variables = write_function_header(
f, eq, default_values, comment)
f = func_name # return the value...
func_msg += ' # i=0 write_function_header\n'
#print(func_msg)
else:
out = f
func_msg += ' %s = %s\n' % (out, eq)
#print('out = %r' % out)
#print('eq = %r' % eq)
func_msg += ' return %s' % f
#print(func_msg)
if func_name in variables:
raise RuntimeError(f'The function name {func_name!r} for DEQATN,{deqatn_id:d} '
f'must not also be an argument; arguments={variables}')
nargs = len(variables)
return func_name, nargs, func_msg
[docs]
def _write_function_line(func_name: str, variables: list[str],
default_values: dict[str, float]) -> str:
"""writes the ``def f(x, y, z=1.):`` part of the function"""
vals = []
is_default = False
for var in variables:
if var in BUILTINS:
var += '_'
if var in default_values:
vals.append('%s=%s' % (var, default_values[var]))
is_default = True
else:
vals.append('%s' % (var))
if is_default:
msg = 'default variables must be set at the end of the function\n'
msg += 'variables = %s\n' % variables
msg += 'default_values = %s' % default_values
raise RuntimeError(msg)
vals2 = ', '.join(vals)
msg = 'def %s(%s):\n' % (func_name, vals2)
return msg
[docs]
def _write_variables(variables: list[str]) -> str:
"""type checks the inputs"""
msg = ' try:\n'
for var in variables:
if var in BUILTINS:
var += '_'
#msg += " assert isinstance(%s, float), '%s is not a float; type(%s)=%s' % (%s)")
#msg += ' %s = float(%s)\n' % (var, var)
msg += ' if isinstance(%s, (int, str)):\n' % var
msg += ' %s = float(%s)\n' % (var, var)
msg += ' except Exception:\n'
msg += ' print(locals())\n'
msg += ' raise\n'
return msg