require 'matrix'
require 'bigdecimal'
require 'mathn'
print("regular mesh\n")
Delaunay = [[[0.0, 0.0, 10.0], [0.0, 800.0, 10.0], [160.0, 280.0, 200.0]], [[800.0, 800.0, 10.0], [800.0, 0.0, 10.0], [600.0, 280.0, 100.0]], [[800.0, 0.0, 10.0], [0.0, 0.0, 10.0], [600.0, 280.0, 100.0]], [[0.0, 800.0, 10.0], [800.0, 800.0, 10.0], [440.0, 600.0, 300.0]], [[440.0, 600.0, 300.0], [600.0, 280.0, 100.0], [160.0, 280.0, 200.0]], [[600.0, 280.0, 100.0], [0.0, 0.0, 10.0], [160.0, 280.0, 200.0]], [[800.0, 800.0, 10.0], [600.0, 280.0, 100.0], [440.0, 600.0, 300.0]], [[160.0, 280.0, 200.0], [0.0, 800.0, 10.0], [440.0, 600.0, 300.0]]]
num_tri = Delaunay.length
def point_seg(s1, s2, p)
u = [s2[0] - s1[0], s2[1] - s1[1]] #s1s2
v = [p[0] - s1[0], p[1] - s1[0]] #s1p
#print "u = ", u, " v = ", v, "\n"
cross = u[0] * v[1] - v[0] * u[1]
if cross == 0
if (((s1[0] - p[0]) * (s2[0] - p[0])) <= 0) && (((s1[1] - p[1]) * (s2[1] - p[1])) <= 0)
#print "on edge\n"
return true
else
return false
end
else
return false
end
end
def inside_tri(tri, p)
a = tri[0]
b = tri[1]
c = tri[2]
m = (a[0] - p[0]) * (b[1] - a[1]) - (b[0] - a[0]) * (a[1] - p[1])
n = (b[0] - p[0]) * (c[1] - b[1]) - (c[0] - b[0]) * (b[1] - p[1])
k = (c[0] - p[0]) * (a[1] - c[1]) - (a[0] - c[0]) * (c[1] - p[1])
#print "m: ", m, "n: ", n, "k: ", k, "\n"
if ((m > 0) && (n > 0) && (k > 0)) || ((m < 0) && (n < 0) && (k < 0)) || (m == 0) || (n == 0) || (k == 0)
#print "inside\n"
return true
end
return false
end
def in_tri(tri, p)
#print tri, "\n"
if point_seg(tri[0], tri[1], p)
#point on seg 0-1
return [true, 0]
elsif point_seg(tri[1], tri[2], p)
#point on seg 1-2
return [true, 1]
elsif point_seg(tri[0], tri[2], p)
#point on seg 0-2
return [true, 2]
elsif inside_tri(tri, p)
return [true, 3]
end
return [false, -1]
end
def dist(x0, y0, x1, y1)
return Math.sqrt((x1 - x0) ** 2 + (y1 - y0) ** 2)
end
def polsum(x0, y0, x1, y1, x2, y2)
a = dist(x0, y0, x1, y1)
b = dist(x1, y1, x2, y2)
c = dist(x0, y0, x2, y2)
return (a + b + c)/2.0
end
def sq(x0, y0, x1, y1, p1, p2)
pols = polsum(x0, y0, x1, y1, p1, p2)
a = dist(x0, y0, x1, y1)
b = dist(x1, y1, p1, p2)
c = dist(x0, y0, p1, p2)
return Math.sqrt(pols * (pols - a) * (pols - b) * (pols - c))
end
def ro(x0, y0, x1, y1, p1, p2)
sqr = sq(x0, y0, x1, y1, p1, p2)
a = dist(x0, y0, x1, y1)
return (2 * sqr) / a
end
def count_rec(tr_ind, p, triang)
#tr = triang[tr_ind]
tr = Marshal.load(Marshal.dump(triang[tr_ind]))
v0 = tr[0]
v1 = tr[1]
v2 = tr[2]
#count distance between point and each edge
r0 = ro(v0[0], v0[1], v1[0], v1[1], p[0], p[1])
r1 = ro(v1[0], v1[1], v2[0], v2[1], p[0], p[1])
r2 = ro(v0[0], v0[1], v2[0], v2[1], p[0], p[1])
rarr = [r0, r1, r2]
#min less then delta found height, else find max and divide
minind = 0
minval = rarr[minind]
maxind = 0
maxval = rarr[maxind]
for i in (1 .. 2)
if rarr[i] < minval
minind = i
minval = rarr[i]
end
if rarr[i] > maxval
maxind = i
maxval = rarr[i]
end
end
if (minval < 10)
#count height
else
#delete old triangle, build 2 new, find in which triangle point lie and run search from tis triangle
triang.delete_at(tr_ind)
if (maxind == 0)
tr1 = [v0, v2, []]
end
def count_tri(tr_ind, p, triang)
interp = count_rec(tr_ind, p, triang)
#return value in p
return (first * secondi * third).element(0, 0)
end
#create grid
LOW = 0
HIGH = 800
kol = 8
step = (HIGH - LOW)/kol
print "step: ", step, "\n"
grid = Hash.new
npx = []
npy = []
npz = []
#i = 0
#j = 0
for i in (0 .. kol)
for j in (0 .. kol)
#print "i: ", i, "j: ", j, " istep ", i * step, " jstep ", j * step, "\n"
#print "will write: ", [i * step, j * step, 0], "\n"
grid[[i, j]] = [i * step, j * step]
#print "current grid: ", grid, "\n"
#j += 1
end
#i += 1
end
print "end grid: \n", grid, "\n"
for i in (0 .. kol)
for j in (0 .. kol)
#проходим в цикле по треугольникам, проверяем лежит ли точка внутри, если нет,
#проверяем лежит ли на каком то ребре,
#когда нашли где находится точка, то выполняем расчет
cur_p = grid[[i, j]]
for k in (0 .. (Delaunay.length - 1))
#print "ex cur_p: ", cur_p, "\n"
#print "ex in k ", k, "\n"
cur_tri = Delaunay[k]
#print "cur_p: ", cur_p, "\n"
res = in_tri(cur_tri, cur_p)
if res[0]
#посчитать высоту
#print "in k ", k, " triangle: ", cur_tri, "\n"
grid[[i, j]] = [grid[[i, j]][0], grid[[i, j]][1], count_tri(k, cur_p, Marshal.load(Marshal.dump(Delaunay)))]
npx += [grid[[i, j]][0]]
npy += [grid[[i, j]][1]]
npz += [grid[[i, j]][2]]
break
end
end
end
end
#print "grid with z: ", grid, "\n"
#grid.each { |key,value| print "vertex(", value[0], ", ", value[1], ", ", value[2], ")\n" }
#for i in (0 .. 66)
# print "vertex(", npx[i], ", ", npy[i], ", ", npz[i], ");\n"
#end
#print npx.length, "\n", npy.length, "\n", npz.length, "\n"
print npx, "\n"
print npy, "\n"
print npz, "\n"
#q = Marshal.load(Marshal.dump(Delaunay))
#print "oldD ", Delaunay, "\n"
#q[0][0][0] = 10000
#print "newq ", q, "\n"
#print "newD ", Delaunay, "\n"
=begin
print "test plot\n"
npxx = []
npyy = []
npzz = []
for i in Delaunay
npxx += i[0]
npyy += i[1]
npzz += i[2]
end
print npxx.length, npyy.length, npzz.length, "\n"
print npxx, "\n"
print npyy, "\n"
print npzz, "\n"
=end