🕊️ ✍️ 😀 Algoritma untuk mencari volume dan pusat massa polihedron 🙏🏿 👩🏻‍✈️ 🧜

Mungkin semua orang tahu algoritma ini, tetapi "pihak berwenang bersembunyi dari saya." Saya menemukan deskripsi verbal di halaman ketiga dari mesin pencari di arsip terjemahan otomatis dari forum berbahasa Inggris. Tampaknya bagi saya bahwa uraian terperinci (dan dengan kode) layak untuk habrosta.

Jadi, misalnya, Anda perlu menghasilkan massa untuk mainan dan di suatu tempat dalam proses untuk menyingkirkan mereka yang tidak berdiri di atas kaki mereka. Untuk melakukan ini, Anda perlu menemukan pusat massa gerombolan (dan ini hampir sama dengan menemukan volumenya) dan pastikan bahwa itu ada di suatu tempat di atas kaki gerombolan.

Massa adalah polyhedron, untuk kesederhanaan kami menganggap bahwa polyhedron hanya terdiri dari segitiga (dalam algoritma di dalamnya terdapat rumus area Gaussian , sehingga Anda dapat mengembangkannya untuk polyhedron apa pun, tetapi mengapa ...). Selain itu, polyhedron tidak boleh memiliki persimpangan sendiri dan membatasi volume tertutup, sebagaimana layaknya polyhedron yang layak.

(yah, seperti itu)

Sebuah UPD kecil yang menjelaskan mengapa pada KDPV massa kanan tidak baik-baik saja, tetapi kiri OK:
Gambar yang tepat tidak OK karena massa akan jatuh ke depan, karena pusat massanya diperluas melampaui area pendukung. Area pendukung dari poligon yang berdiri di permukaan didefinisikan sebagai poligon minimum di dalam semua titik pada permukaan. Dalam kasus kiri, area penopang bergeser ke pusat massa dan lebih banyak (karena cakar dinosaurus lebih besar), dan pada gambar kanan area itu sendiri lebih kecil dan lebih dekat ke ekor.
Rasio area referensi dengan pusat massa akan menjadi seperti ini:

Saya akan segera mulai dengan kode pencarian volume (Python, input - daftar poin dan matriks transisi):

beberapa kode

def RecSetDirsTriangles(para, Connects, TR): """ ,           """ #1.   , ,     for i in range(0,len(Connects)): if i != para[0] and i != para[1] and Connects[i][para[0]] and Connects[i][para[1]]: #  ! fl = 1 for T in TR: if i in T and para[1] in T and para[0] in T: fl = 0 #    break if fl: # ! TR += [(para[1],para[0],i)] Recc((para[0], i) , Connects, TR) Recc((i, para[1]) , Connects, TR) def FindV(dots, Connects): """ .   - dots     [x, y, z], Connects -  , Connects[i][j]=1      i, j,  =0 """ #1.      TR = [] for i in range(1,len(Connects)):#        - if Connects[i][0]: for j in range(i+1, len(Connects)): if Connects[0][j] and Connects[i][j]: TR += [(0,i,j)] break RecSetDirsTriangles((0,i),Connects, TR) break print(" : ", len(TR)) #2.        V = 0 for T in TR: ''' : x1y2 x2y3 x3y1 x2y1 x3y2 x1y3''' S = 0.5 * (dots[T[0]][0]*dots[T[1]][1] + dots[T[1]][0]*dots[T[2]][1] + dots[T[2]][0]*dots[T[0]][1] - dots[T[1]][0]*dots[T[0]][1] - dots[T[2]][0]*dots[T[1]][1] - dots[T[0]][0]*dots[T[2]][1]) #S   +  -    ,    V += S*(dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2])/3 #    ... return math.fabs(V)

Inti dari algoritma ini adalah untuk mempertimbangkan volume gambar yang membentuk wajah polyhedron yang “jatuh” ke bidang xy. Untuk melakukan ini, Anda perlu mengetahui area proyeksi segitiga dan tanda yang digunakan untuk menambahkan volume gambar (prisma terpotong). Bahkan, jika segitiga dipesan terlebih dahulu, volume dan tanda dikurangi menjadi satu perhitungan.

Oleh karena itu, hal pertama yang fungsi rekursif mengumpulkan segitiga dari adalah input. Merakit sedemikian rupa sehingga ketika melihat "luar" pada sebuah polyhedron, arah berkeliling segitiga adalah sama (idealnya berlawanan arah jarum jam; jika Anda mengambil arah searah jarum jam, hasilnya akan benar, tetapi negatif - oleh karena itu, volume modulus diberikan untuk pengembalian).

Ini sangat sederhana untuk dicapai - ambil segitiga (poin a1, a2, a3), cari tetangganya dan buat daftar dua simpul yang cocok dalam urutan terbalik (misalnya, seperti ini: a2, a1, b1).
Ternyata sesuatu seperti ini:

Sekarang, jika kita memproyeksikan segitiga seperti itu ke bidang xy, maka urutan traversal untuk proyeksi segitiga "atas" akan bertepatan dengan yang awalnya dipilih, dan urutan traversal untuk proyeksi segitiga "bawah" akan mengubah arahnya. Akibatnya, itu akan mengubah tanda dan luas segitiga ini, dihitung dengan rumus Gauss. Di sini, segitiga "bawah" - konsep kondisional - berarti volume tepat di bawahnya tidak termasuk dalam volume polihedron. Segitiga "bawah" dari polyhedron non-cembung bisa lebih tinggi dari yang "atas".

Setelah langkah-langkah awal ini, untuk menghitung volume total polyhedron, Anda hanya perlu menambahkan (dengan memperhitungkan tanda, yang diperoleh "dengan sendirinya") semua volume prisma terpotong yang dikumpulkan dari wajah dan proyeksi wajah-wajah ini pada bidang xy. Dan volume prisma dianggap sebagai produk daerah (Gaussian, dengan tanda) dan rata-rata aritmatika z-koordinat dari simpul segitiga.

Jika polyhedron memotong bidang xy, maka ketika menghitung volume, semua tanda akan saling membatalkan dan hasilnya tetap benar (Anda hanya perlu mengambil ketinggian prisma tanpa modul).

(entah bagaimana prisma terpotong tampak seperti)

Dengan pencarian pusat massa, semuanya hampir sama. Demikian pula, kita perlu menemukan pusat massa untuk setiap prisma terpotong dan meringkasnya secara bersamaan, mengalikannya dengan volume prisma (diasumsikan bahwa massa didistribusikan secara merata di seluruh volume dan satu dapat diganti dengan yang lain). Untuk menemukan pusat massa prisma terpotong, perlu untuk menghitung pusat massa dua tetrahedra (fungsi +1) dan satu prisma biasa. Algoritma juga "tidak memburuk" jika polyhedron melintasi bidang xy (dan di sini bisa menjadi reproduksi Magritte).

(kedua tetrahedra ini, bertanda merah dan merah, bersama dengan prisma segitiga (di bawah tetrahedron merah) membentuk prisma terpotong yang diinginkan. Kita perlu menemukan pusat massa dan volume ketiga figur. Penandaan kira-kira sesuai dengan penunjukan dalam kode)

Kode yang menghitung ini dan itu:

sedikit lebih banyak kode

 def RecSetDirsTriangles(para, Connects, TR): #1.   , ,     for i in range(0,len(Connects)): if i != para[0] and i != para[1] and Connects[i][para[0]] and Connects[i][para[1]]: #  ! fl = 1 for T in TR: if i in T and para[1] in T and para[0] in T: fl = 0 break if fl: # ! TR += [(para[1],para[0],i)] Recc((para[0], i) , Connects, TR) Recc((i, para[1]) , Connects, TR) def TetrV(mas):#dot1, dot2, dot3, dot4): """   """ M = np.zeros((3,3),float) for i in range(1,4): for j in range(0,3): M[i-1][j] = mas[i][j] - mas[0][j] #print(M) return math.fabs(np.linalg.det(M)/6) def FindVandCM(dots, Connects): """     """ #1.      TR = [] for i in range(1,len(Connects)): #        - if Connects[i][0]: for j in range(i+1, len(Connects)): if Connects[0][j] and Connects[i][j]: TR += [(0,i,j)] break RecSetDirsTriangles((0,i),Connects, TR) break print(" : ", len(TR)) #2.   ,          V = 0 CM = [0, 0, 0] for T in TR: ''' : x1y2 x2y3 x3y1 x2y1 x3y2 x1y3''' S = 0.5 * (dots[T[0]][0]*dots[T[1]][1] + dots[T[1]][0]*dots[T[2]][1] + dots[T[2]][0]*dots[T[0]][1] - dots[T[1]][0]*dots[T[0]][1] - dots[T[2]][0]*dots[T[1]][1] - dots[T[0]][0]*dots[T[2]][1]) #S   +  -    ,    V += S*(dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2])/3 #    ... #c        c1 = ((dots[T[0]][0] + dots[T[1]][0] + dots[T[2]][0])/3, (dots[T[0]][1]+ dots[T[1]][1]+ dots[T[2]][1])/3) #    hm = min([dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]]) hM = max([dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]]) indM = [dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]].index(hM) indm = [dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]].index(hm) V3 = S * hm if indM == indm: # ! CM[0] += V3*c1[0] CM[1] += V3*c1[1] CM[2] += V3*hm/2 continue L = [0,1,2] L.remove(indM) L.remove(indm) indmidle = L[0] dots1 = [dots[T[0]], dots[T[1]], dots[T[2]], (dots[T[indM]][0], dots[T[indM]][1] , hm)] #  V1 = TetrV(dots1) if S < 0: V1 = -V1 V2 = S * ( dots[T[indmidle]][2] - hm)/3 #V3 = S * hm CM[0] += V1*((dots[T[0]][0] + dots[T[1]][0] + dots[T[2]][0] + dots[T[indM]][0])/4) + V2*((dots[T[0]][0] + dots[T[1]][0] + dots[T[2]][0] + dots[T[indmidle]][0])/4) + V3*c1[0] CM[1] += V1*((dots[T[0]][1] + dots[T[1]][1] + dots[T[2]][1] + dots[T[indM]][1])/4) + V2*((dots[T[0]][1] + dots[T[1]][1] + dots[T[2]][1] + dots[T[indmidle]][1])/4) + V3*c1[1] CM[2] += V1*((dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2] + hm)/4) + V2*((dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2] + hm)/4) + V3*hm/2 CM[0] = CM[0]/V CM[1] = CM[1]/V CM[2] = CM[2]/V return (math.fabs(V), CM)

Sepotong algoritma di mana arah segitiga dipertimbangkan dan digunakan untuk memahami volume eksternal dan internal adalah langkah yang sangat kuat, dapat digunakan banyak saat Anda bekerja dengan polyhedra. Misalnya, jika Anda perlu menghitung arah normals "out" - cukup untuk mengetahui arah "berlawanan arah jarum jam" untuk satu wajah - dan voila!

(tebak filmnya!)

Algoritma untuk mencari volume dan pusat massa polihedron

More articles: