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