from moedebug import *
+# When cse() is called for something containing BaseVector, it
+# produces infinite recursion.
+#def cse(x, *a, **kw): return ((), (x,))
+
def dprint(*args):
if not dbg_enabled(): return
print(*args)
print('\n =\n')
pprint(cse(eval(vn)))
+def sqnorm(v): return v & v
+
N = CoordSysCartesian('N')
calculated = False
F0, G0 = vector_symbols('F0 G0')
En, Hn = vector_symbols('En Hn')
-EFl, GHl = symbols('EFl GHl')
+EFlq, HGlq = symbols('EFlq HGlq')
def vector_component(v, ix):
return v.components[N.base_vectors()[ix]]
# x array in numerical algorithm has:
# N x 3 coordinates of points 0..N-3
-# 1 EFl = length parameter |EF| for point 1
-# 1 GHl = length parameter |GH| for point N-2
+# 1 EFlq = sqrt of length parameter |EF| for point 1
+# 1 HGlq = sqrt of length parameter |HG| for point N-2
# fixed array in numerical algorithm has:
# 4 x 3 E, En, H, Hn
#def subst_vect():
+iterations = []
+
+class SomeIteration():
+ def __init__(ar, names, size, expr):
+ ar.names_string = names
+ ar.names = names.split(' ')
+ ar.name = ar.names[0]
+ ar.size = size
+ ar.expr = expr
+ if dbg_enabled():
+ print('\n ' + ar.name + '\n')
+ print(expr)
+ iterations.append(ar)
+
+ def gen_calculate_cost(ar):
+ ar._gen_array()
+ cprint('for (P=0; P<(%s); P++) {' % ar.size)
+ ar._cassign()
+ cprint('}')
+
+class ScalarArray(SomeIteration):
+ def _gen_array(ar):
+ cprint('double A_%s[%s];' % (ar.name, ar.size))
+ def gen_references(ar):
+ for ai in range(0, len(ar.names)):
+ ar._gen_reference(ai, ar.names[ai])
+ def _gen_reference(ar, ai, an):
+ cprintraw('#define %s A_%s[P%+d]' % (an, ar.name, ai))
+ def _cassign(ar):
+ cassign(ar.expr, ar.name, 'tmp_'+ar.name)
+ def s(ar):
+ return symbols(ar.names_string)
+
+class CoordArray(ScalarArray):
+ def _gen_array(ar):
+ cprint('double A_%s[%s][3];' % (ar.name, ar.size))
+ def _gen_reference(ar, ai, an):
+ ScalarArray._gen_reference(ar, ai, an)
+ gen_point_coords_macro(an)
+ def _cassign(ar):
+ cassign_vector(ar.expr, ar.name, 'tmp_'+ar.name)
+ def s(ar):
+ return vector_symbols(ar.names_string)
+
+class CostComponent(SomeIteration):
+ def __init__(cc, size, expr):
+ cc.size = size
+ cc.expr = expr
+ iterations.append(cc)
+ def gen_references(cc): pass
+ def _gen_array(cc): pass
+ def _cassign(cc):
+ cassign(cc.expr, 'P_cost', 'tmp_cost')
+ cprint('cost += P_cost;')
+
def calculate():
global calculated
if calculated: return
- Q = B + 0.5 * (B - A)
- R = C + 0.5 * (C - D)
- QR = R - Q
- BC = C - B
- cost_ABCD = (QR & QR) / (BC & BC)
- global cost_ABCD
- dbg('cost_ABCD')
- dprint(A)
+ # ---------- actual cost computation formulae ----------
global F, G
- F = E + En * EFl
- G = H + Hn * GHl
+ F = E + En * pow(EFlq, 2)
+ G = H + Hn * pow(HGlq, 2)
+
+ global a,b, al,bl, au,bu
+ a, b = CoordArray ('a_ b_', 'NP-1', B-A ).s() # [mm]
+ al, bl = ScalarArray('al bl', 'NP-1', a.magnitude() ).s() # [mm]
+ au, bu = CoordArray ('au bu', 'NP-1', a / al ).s() # [1]
+
+ tan_theta = (au ^ bu) / (au & bu) # [1] bending
- # diff_A_i = cost_ABCD.diff(vector_component(A, 0))
- # global diff_A_i
- # dbg('diff_A_i')
+ global mu, nu
+ mu, nu = CoordArray ('mu nu', 'NP-2', tan_theta ).s() # [1]
+
+ CostComponent('NP-3', sqnorm(mu - nu)) # [1]
+
+ dl2 = pow(al - bl, 2) # [mm^2]
+ CostComponent('NP-2', dl2 / (al*bl)) # [1]
+
+ # ---------- end of cost computation formulae ----------
calculated = True
gen_point_index_macro('E', 'INPUT', '3*0')
gen_point_index_macro('F0', 'INPUT', '3*1')
- gen_point_index_macro('G0', 'INPUT', '3*(N-2)')
- gen_point_index_macro('H', 'INPUT', '3*(N-1)')
- cprintraw( '#define NINPUT ( 3*(N-0) )')
+ gen_point_index_macro('G0', 'INPUT', '3*(NP-2)')
+ gen_point_index_macro('H', 'INPUT', '3*(NP-1)')
+ cprintraw( '#define NINPUT ( 3*(NP-0) )')
gen_point_index_macro('En', 'PREP', '3*0')
gen_point_index_macro('Hn', 'PREP', '3*1')
- cprintraw( '#define PREP_N (3*2)')
+ cprintraw( '#define NPREP (3*2)')
- cprintraw('#define X_N_DIRECT 3*(N-4)')
+ cprintraw('#define NX_DIRECT 3*(NP-4)')
cprint('#define POINT(PP) (')
- cprint(' (PP) == 0 ? E :')
- cprint(' (PP) == 1 ? F :')
- cprint(' (PP) == N-2 ? G :')
- cprint(' (PP) == N-1 ? H :')
+ cprint(' (PP) == 0 ? E :')
+ cprint(' (PP) == 1 ? F :')
+ cprint(' (PP) == NP-2 ? G :')
+ cprint(' (PP) == NP-1 ? H :')
cprint(' &X[3*((PP)-2)]')
cprintraw(')')
- cprintraw('#define EFl X[ X_N_DIRECT + 0 ]')
- cprintraw('#define GHl X[ X_N_DIRECT + 1 ]')
- cprintraw('#define X_N ( X_N_DIRECT + 2 )')
+ cprintraw('#define EFlq X[ NX_DIRECT + 0 ]')
+ cprintraw('#define HGlq X[ NX_DIRECT + 1 ]')
+ cprintraw('#define NX ( NX_DIRECT + 2 )')
for ai in range(0, len(abcd)):
cprintraw('#define %s POINT(P%+d)' % (abcd[ai], ai))
gen_point_coords_macro(abcd[ai])
+ for si in iterations:
+ si.gen_references()
+
cprintraw('')
def gen_prepare():
cprint('#define PREPARE')
- cprint('memcpy(X, &INPUT[3*2], sizeof(double) * X_N_DIRECT)')
- cassign_vector((F0 - E).normalize(), 'En', 'tmp_En')
- cassign_vector((G0 - H).normalize(), 'Hn', 'tmp_Hn')
+ cprint('memcpy(X, &INPUT[3*2], sizeof(double) * NX_DIRECT);')
+ for EH,EHs,FG0,FGs in ((E,'E', F0,'F'),
+ (H,'H', G0,'G')):
+ EFHGv = FG0 - EH
+ EFHGl = EFHGv.magnitude()
+ EFHGlq = sqrt(EFHGl)
+ cassign_vector(EFHGv/EFHGl, EHs+'n', 'tmp_'+EHs)
+ cassign(EFHGlq, EHs+FGs+'lq', 'tmp_l'+EHs)
cprintraw('')
def gen_calculate_FG():
+ cprintraw('#define DECLARE_F_G double F[3], G[3];')
cprint('#define CALCULATE_F_G')
cassign_vector(F,'F','tmp_F')
cassign_vector(G,'G','tmp_G')
def gen_calculate_cost():
cprint('#define CALCULATE_COST')
- cassign(cost_ABCD,'P_cost','tmp_P_cost')
+ cprint('double cost=0, P_cost;')
+ for si in iterations:
+ si.gen_calculate_cost()
cprintraw('')
def gen_C():
~ian / moebius3.git / blobdiff