2010-12-10 30 views
2

Je veux calculer la masse d'une sphère basée sur une distribution de densité inhomogène discrète tridimensionnelle. Disons qu'un ensemble de 3x3x3 cubes de différentes densités est inscrit par une sphère. Quel est le moyen le plus rapide pour résumer les masses partitionnées en utilisant Python?Comment calculer la masse d'une sphère inhomogène?

I a essayé de calculer le volume sous l'équation mathématique pour une sphère: x^2 + y^2 + z 2 = R^2 pour la plage de l'un des cubes à l'aide de scipy.integrate.dblquad. Cependant, le résultat n'est valide que si les limites sont plus petites que le rayon de la sphère et le calcul répétitif pour disons que 50'000 sphères avec 27 cubes chacune seraient assez lentes. D'autre part, l'équation habituelle pour les évaluations CoM ne pourrait pas être utilisée à mon avis, en raison de la distribution de masse plutôt grossière et discrète.

+1

« plus petit que le rayon » devrait être « inférieur au rayon ». Utilisez "que" pour les comparaisons, "puis" pour la séquence temporelle. –

+0

Si les cubes sont homogènes ... il suffit d'additionner la masse des cubes.Si les cubes ne sont pas homogènes, voyons la formule de densité. Vous n'avez pas besoin de la CoM pour calculer la masse. –

+0

Je pense que ce n'est pas une vraie question à moins que la densité des cubes ne soit spécifiée. –

Répondre

0

Je ne peux pas obtenir votre signification exacte d'inscrit par une sphère. Aussi, je n'ai pas essayé le scipy.integrate. Toutefois, voici quelques-uns:

Définissez un cube 3x3x3 à la densité de l'unité. Ensuite, prenez l'intégration pour chaque cube, donc vous devriez avoir le cube de volumeV_ijk ici. Maintenant pour chaque sphère, vous pouvez obtenir la masse de chaque sphère en additionnant V_ijk*D_ijk, où D_ijk est la densité de la sphère.

Cela devrait être beaucoup plus rapide car vous n'avez pas besoin de faire l'intégration maintenant.

0

Vous pouvez obtenir une formule analytique pour le volume d'intersection entre un cube (ou un prisme rectangulaire) et une sphère. Ce ne sera pas facile, mais cela devrait être possible. Je l'ai fait pour un triangle arbitraire et un cercle en 2D. L'idée de base est de décomposer l'intersection en parties plus simples, comme les secteurs de tétraèdres et de triangles sphériques volumétriques, pour lesquels des formules de volume relativement simples sont connues. Le principal difficile est de considérer tous les cas possibles d'intersections. Heureusement, les deux objets sont convexes, ce qui vous garantit un seul volume d'intersection convexe.

Une méthode approximative pourrait être de simplement subdiviser les cubes jusqu'à ce que votre algorithme d'intégration numérique approximatif fonctionne; cela devrait toujours être relativement rapide. Savez-vous à propos de Pick's Theorem? Cela ne fonctionne qu'en 2D, mais il y a, je crois, 3D generalizations.

1

Timing Expérience

Vous n'avez pas spécifié vos contraintes de temps, donc je l'ai fait une petite expérience avec un joli paquet d'intégration.

Sans optimisation, chaque intégrale en coordonnées sphériques peut être évaluée en 0.005 secondes dans un ordinateur portable standard si les densités de cubes sont des fonctions simples.

Tout comme une référence, c'est le programme Mathematica:

[email protected]; 
(* Define a cuboid as density function *) 
iP = IntegerPart; 
f[{x_, y_, z_}, {lx_, ly_, lz_}] := iP[x - lx] + iP[y - ly] + iP[z - lz] /; 
    lx <= x <= lx + 3 && ly <= y <= ly + 3 && lz <= z <= lz + 3; 

f[{x_, y_, z_}, {lx_, ly_, lz_}] := Break[] /; True; 

Timing[Table[s = RandomReal[{0, 3}, 3]; (*sphere center random*) 
    sphereRadius = Min[Union[s, 3 - s]]; (*max radius inside cuboid *) 
    NIntegrate[(f[{x, y, z} - s, -s] /. (*integrate in spherical coords *) 
     {x -> r [email protected] [email protected], 
     y -> r [email protected] [email protected], 
     z -> r [email protected]}) r^2 [email protected], 
     {r, 0, sphereRadius}, {th, 0, 2 Pi}, {phi, 0, Pi}], 
     {10000}]][[1]] 

Le résultat est 52 secondes pour 10^4 itérations.

besoin donc peut-être vous ne pas optimiser beaucoup ...