Results of my master's thesis
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

baker.py 3.7KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147
  1. from grid import exact_func
  2. import numpy as np
  3. import sys
  4. def get_phis(X, r):
  5. """
  6. The get_phis function is used to get barycentric coordonites for a point on a triangle.
  7. X -- the destination point (2D)
  8. X = [0,0]
  9. r -- the three points that make up the triangle (2D)
  10. r = [[-1, -1], [0, 2], [1, -1]]
  11. this will return [0.333, 0.333, 0.333]
  12. """
  13. # baker: eq 7
  14. A = np.array([
  15. [ 1, 1, 1],
  16. [r[0][0], r[1][0], r[2][0]],
  17. [r[0][1], r[1][1], r[2][1]],
  18. ])
  19. b = np.array([ 1,
  20. X[0],
  21. X[1]
  22. ])
  23. try:
  24. phi = np.linalg.solve(A,b)
  25. except:
  26. print >> sys.stderr, "warning: calculation of phis yielded a linearly dependant system"
  27. phi = np.dot(np.linalg.pinv(A), b)
  28. return phi
  29. def get_phis_3D(X, r):
  30. """
  31. The get_phis function is used to get barycentric coordonites for a point on a triangle.
  32. X -- the destination point (3D)
  33. X = [0,0,0]
  34. r -- the four points that make up the tetrahedron (3D)
  35. r = [[-1, -1], [0, 2], [1, -1]]
  36. this will return [0.333, 0.333, 0.333]
  37. """
  38. # baker: eq 7
  39. A = np.array([
  40. [ 1, 1, 1, 1 ],
  41. [r[0][0], r[1][0], r[2][0], r[3][0]],
  42. [r[0][1], r[1][1], r[2][1], r[3][1]],
  43. [r[0][2], r[1][2], r[2][2], r[3][2]],
  44. ])
  45. b = np.array([ 1,
  46. X[0],
  47. X[1],
  48. X[2]
  49. ])
  50. try:
  51. phi = np.linalg.solve(A,b)
  52. except:
  53. print >> sys.stderr, "warning: calculation of phis yielded a linearly dependant system"
  54. phi = np.dot(np.linalg.pinv(A), b)
  55. return phi
  56. def qlinear(X, r, q):
  57. """
  58. this calculates the linear portion of q from X to r
  59. X = destination point
  60. r = simplex points
  61. q = CFD quantities of interest at the simplex points
  62. """
  63. phis = get_phis(X, r)
  64. qlin = sum([q_i * phi_i for q_i, phi_i in zip(q[:len(phis)], phis)])
  65. return qlin
  66. def qlinear_3D(X, r, q):
  67. """
  68. this calculates the linear portion of q from X to r
  69. X = destination point
  70. r = simplex points
  71. q = CFD quantities of interest at the simplex points(r)
  72. """
  73. phis = get_phis_3D(X, r)
  74. qlin = sum([q_i * phi_i for q_i, phi_i in zip(q[:len(phis)], phis)])
  75. return qlin
  76. def run_baker(X, g, tree, extra_points = 3, verbose = False):
  77. """
  78. This is the main function to call to get an interpolation to X from the tree
  79. X -- the destination point (2D)
  80. X = [0,0]
  81. g -- the grid object
  82. tree -- the kdtree search object (built from the g mesh)
  83. """
  84. (dist, indicies) = tree.query(X, 3 + extra_points)
  85. nn = [g.points[i] for i in indicies]
  86. nq = [g.q[i] for i in indicies]
  87. # calculate values only for the triangle
  88. phi = get_phis(X, nn[:3])
  89. qlin = qlinear(X, nn[:3], nq[:3])# nq[0] * phi[0] + nq[1] * phi[1] + nq[2] * phi[2]
  90. error_term = 0.0
  91. if extra_points != 0:
  92. B = [] # baker eq 9
  93. w = [] # baker eq 11
  94. for index in indicies[3:]:
  95. (phi1,phi2,phi3) = get_phis(g.points[index], nn)
  96. B.append([phi1 * phi2, phi2*phi3, phi3*phi1])
  97. w.append(g.q[index] - qlinear(g.points[index], nn, nq))
  98. B = np.array(B)
  99. w = np.array(w)
  100. A = np.dot(B.T, B)
  101. b = np.dot(B.T, w)
  102. # baker solve eq 10
  103. try:
  104. (a, b, c) = np.linalg.solve(A,b)
  105. except:
  106. print >> sys.stderr, "warning: linear calculation went bad, resorting to np.linalg.pinv"
  107. (a, b, c) = np.dot(np.linalg.pinv(A), b)
  108. error_term = a * phi[0] * phi[1]\
  109. + b * phi[1] * phi[2]\
  110. + c * phi[2] * phi[0]
  111. q_final = qlin + error_term
  112. return qlin, error_term, q_final