Calcul du moment d'inertie d'un polygone 2D concave par rapport à son origine

Question

Je souhaite calculer le moment d'inertie d'un polygone concave (2D). J'ai trouvé ça sur internet. Mais je ne suis pas très sûr de savoir comment interpréter la formule ...

Formula http://img101.imageshack.us/img101/8141/92175941c14cadeeb956d8f.gif

1) Cette formule est correcte?
2) Si oui, ma conversion en C++ est-elle correcte?

float sum (0); 
for (int i = 0; i < N; i++) // N = number of vertices 
{ 
    int j = (i + 1) % N; 
    sum += (p[j].y - p[i].y) * (p[j].x + p[i].x) * (pow(p[j].x, 2) + pow(p[i].x, 2)) - (p[j].x - p[i].x) * (p[j].y + p[i].y) * (pow(p[j].y, 2) + pow(p[i].y, 2)); 
} 
float inertia = (1.f/12.f * sum) * density; 

Martijn

Answer 1

#include <math.h> //for abs 
float dot (vec a, vec b) { 
    return (a.x*b.x + a.y*b.y); 
} 
float lengthcross (vec a, vec b) { 
    return (abs(a.x*b.y - a.y*b.x)); 
} 
... 
do stuff 
... 
float sum1=0; 
float sum2=0; 
for (int n=0;n<N;++n) { //equivalent of the Σ 
    sum1 += lengthcross(P[n+1],P[n])* 
      (dot(P[n+1],P[n+1]) + dot(P[n+1],P[n]) + dot(P[n],P[n])); 
    sum2 += lengthcross(P[n+1],P[n]); 
} 
return (m/6*sum1/sum2); 

Edit: Beaucoup de petites maths change

Answer 2

Je pense que vous avez plus de travail à faire que le simple fait traduire des formules en code. Vous devez comprendre exactement ce que cette formule signifie. Lorsque vous avez un polygone 2D, vous avez trois moments d'inertie que vous pouvez calculer par rapport à un système de coordonnées donné: le moment autour de x, le moment autour de y et le moment d'inertie polaire. Il y a un théorème d'axe parallèle qui vous permet de traduire d'un système de coordonnées à un autre. Savez-vous exactement quel moment et quel système de coordonnées s'applique à cette formule?

est ici un code qui pourrait vous aider, ainsi que d'un test JUnit pour prouver que cela fonctionne:

import java.awt.geom.Point2D; 

/** 
* PolygonInertiaCalculator 
* User: Michael 
* Date: Jul 25, 2010 
* Time: 9:51:47 AM 
*/ 
public class PolygonInertiaCalculator 
{ 
    private static final int MIN_POINTS = 2; 

    public static double dot(Point2D u, Point2D v) 
    { 
     return u.getX()*v.getX() + u.getY()*v.getY(); 
    } 

    public static double cross(Point2D u, Point2D v) 
    { 
     return u.getX()*v.getY() - u.getY()*v.getX(); 
    } 

    /** 
    * Calculate moment of inertia about x-axis 
    * @param poly of 2D points defining a closed polygon 
    * @return moment of inertia about x-axis 
    */ 
    public static double ix(Point2D [] poly) 
    { 
     double ix = 0.0; 

     if ((poly != null) && (poly.length > MIN_POINTS)) 
     { 
      double sum = 0.0; 
      for (int n = 0; n < (poly.length-1); ++n) 
      { 
       double twiceArea = poly[n].getX()*poly[n+1].getY() - poly[n+1].getX()*poly[n].getY(); 
       sum += (poly[n].getY()*poly[n].getY() + poly[n].getY()*poly[n+1].getY() + poly[n+1].getY()*poly[n+1].getY())*twiceArea; 
      } 

      ix = sum/12.0; 
     } 

     return ix; 
    } 

    /** 
    * Calculate moment of inertia about y-axis 
    * @param poly of 2D points defining a closed polygon 
    * @return moment of inertia about y-axis 
    * @link http://en.wikipedia.org/wiki/Second_moment_of_area 
    */ 
    public static double iy(Point2D [] poly) 
    { 
     double iy = 0.0; 

     if ((poly != null) && (poly.length > MIN_POINTS)) 
     { 
      double sum = 0.0; 
      for (int n = 0; n < (poly.length-1); ++n) 
      { 
       double twiceArea = poly[n].getX()*poly[n+1].getY() - poly[n+1].getX()*poly[n].getY(); 
       sum += (poly[n].getX()*poly[n].getX() + poly[n].getX()*poly[n+1].getX() + poly[n+1].getX()*poly[n+1].getX())*twiceArea; 
      } 

      iy = sum/12.0; 
     } 

     return iy; 
    } 

    /** 
    * Calculate polar moment of inertia xy 
    * @param poly of 2D points defining a closed polygon 
    * @return polar moment of inertia xy 
    * @link http://en.wikipedia.org/wiki/Second_moment_of_area 
    */ 
    public static double ixy(Point2D [] poly) 
    { 
     double ixy = 0.0; 

     if ((poly != null) && (poly.length > MIN_POINTS)) 
     { 
      double sum = 0.0; 
      for (int n = 0; n < (poly.length-1); ++n) 
      { 
       double twiceArea = poly[n].getX()*poly[n+1].getY() - poly[n+1].getX()*poly[n].getY(); 
       sum += (poly[n].getX()*poly[n+1].getY() + 2.0*poly[n].getX()*poly[n].getY() + 2.0*poly[n+1].getX()*poly[n+1].getY() + poly[n+1].getX()*poly[n].getY())*twiceArea; 
      } 

      ixy = sum/24.0; 
     } 

     return ixy; 
    } 

    /** 
    * Calculate the moment of inertia of a 2D concave polygon 
    * @param poly array of 2D points defining the perimeter of the polygon 
    * @return moment of inertia 
    * @link http://www.physicsforums.com/showthread.php?t=43071 
    * @link http://www.physicsforums.com/showthread.php?t=25293 
    * @link http://stackoverflow.com/questions/3329383/help-me-with-converting-latex-formula-to-code 
    */ 
    public static double inertia(Point2D[] poly) 
    { 
     double inertia = 0.0; 

     if ((poly != null) && (poly.length > MIN_POINTS)) 
     { 
      double numer = 0.0; 
      double denom = 0.0; 
      double scale; 
      double mag; 
      for (int n = 0; n < (poly.length-1); ++n) 
      { 
       scale = dot(poly[n + 1], poly[n + 1]) + dot(poly[n + 1], poly[n]) + dot(poly[n], poly[n]); 
       mag = Math.sqrt(cross(poly[n], poly[n+1])); 
       numer += mag * scale; 
       denom += mag; 
      } 

      inertia = numer/denom/6.0; 
     } 

     return inertia; 
    } 
} 

Voici le test JUnit pour l'accompagner:

import org.junit.Test; 

import java.awt.geom.Point2D; 

import static org.junit.Assert.assertEquals; 

/** 
* PolygonInertiaCalculatorTest 
* User: Michael 
* Date: Jul 25, 2010 
* Time: 10:16:04 AM 
*/ 
public class PolygonInertiaCalculatorTest 
{ 
    @Test 
    public void testTriangle() 
    { 
     Point2D[] poly = 
     { 
      new Point2D.Double(0.0, 0.0), 
      new Point2D.Double(1.0, 0.0), 
      new Point2D.Double(0.0, 1.0) 
     }; 


     // Moment of inertia about the y1 axis 
     // http://www.efunda.com/math/areas/triangle.cfm 
     double expected = 1.0/3.0; 
     double actual = PolygonInertiaCalculator.inertia(poly); 

     assertEquals(expected, actual, 1.0e-6); 
    } 

    @Test 
    public void testSquare() 
    { 
     Point2D[] poly = 
     { 
      new Point2D.Double(0.0, 0.0), 
      new Point2D.Double(1.0, 0.0), 
      new Point2D.Double(1.0, 1.0), 
      new Point2D.Double(0.0, 1.0) 
     }; 

     // Polar moment of inertia about z axis 
     // http://www.efunda.com/math/areas/Rectangle.cfm 
     double expected = 2.0/3.0; 
     double actual = PolygonInertiaCalculator.inertia(poly); 

     assertEquals(expected, actual, 1.0e-6); 
    } 

    @Test 
    public void testRectangle() 
    { 
     // This gives the moment of inertia about the y axis for a coordinate system 
     // through the centroid of the rectangle 
     Point2D[] poly = 
     { 
      new Point2D.Double(0.0, 0.0), 
      new Point2D.Double(4.0, 0.0), 
      new Point2D.Double(4.0, 1.0), 
      new Point2D.Double(0.0, 1.0) 
     }; 

     double expected = 5.0 + 2.0/3.0; 
     double actual = PolygonInertiaCalculator.inertia(poly); 

     assertEquals(expected, actual, 1.0e-6); 

     double ix = PolygonInertiaCalculator.ix(poly); 
     double iy = PolygonInertiaCalculator.iy(poly); 
     double ixy = PolygonInertiaCalculator.ixy(poly); 

     assertEquals(ix, (1.0 + 1.0/3.0), 1.0e-6); 
     assertEquals(iy, (21.0 + 1.0/3.0), 1.0e-6); 
     assertEquals(ixy, 4.0, 1.0e-6); 
    } 
} 

Answer 3

Pour référence, voici un mise en œuvre mutable 2D org.gcs.kinetic.Vector et une implémentation org.jscience.mathematics.vector plus polyvalente et immuable. Cet article sur Calculating a 2D Vector’s Cross Product est également utile.

Answer 4

Je l'ai fait avec Tesselation. Et prendre les MOI tous ensemble.