Topic: Math problem - calculating a volume given points

[VovKa]

Alan, it means that my drawing skills are so bad, that no one is going to believe from the first sight

[dgorsman]

Taking the easy part first - volume of the middle chunk.

V_middle = A * L, where

L = min (P_Fz P_Gz P_Hz) - max ( P_Az P_Bz P_Cz)


Area is bit trickier, we need to create a set of intermediary points.  Using the three points at the lower end we get

I_z = max ( P_Az P_Bz P_Cz)

P_AM = (P_Ax, P_Ay, I_z)
P_BM = (P_Bx, P_By, I_z)
P_CM = (P_Cx, P_Cy, I_z)

Area of a triangle is 0.5 * base * perpendicular height.  We can use any of the three sides as the base, but the height needs to be trig'd.  So using line P_AM P_BM as the base, the distance d is perpendicular from line P_AM P_BM to P_CM:

A = 0.5 * | P_AM P_BM | * ( | P_BM P_CM | * sin P_AM-P_BM-P_CM )


( LISP to come next week, unless somebody is feeling industrious )

[VovKa]

my 3d experience equals to zero, so i had to spend a couple of hours restudying stereometry
here's what i got so far

Code: [Select]

(defun vk_GetArea (CoordsList)
  (abs
    (/ (apply '+
     (mapcar (function
(lambda (p1 p2)
 (* (+ (car p1) (car p2)) (- (cadr p1) (cadr p2)))
)
     )
     CoordsList
     (cons (last CoordsList) CoordsList)
     )
       )
       2.0
    )
  )
)
(defun vk_IsPointOnLine (p1 p2 p / r)
  (equal (+ (distance p p1) (distance p p2))
(distance p1 p2)
1e-8
  )
)
(defun vk_GetOrt (v / k)
  (setq k (/ 1.0 (distance '(0 0 0) v)))
  (mapcar (function (lambda (c) (* c k))) v)
)
(defun vk_GetCrossProduct (v1 v2)
  (list (- (* (cadr v1) (caddr v2)) (* (caddr v1) (cadr v2)))
(- (* (caddr v1) (car v2)) (* (car v1) (caddr v2)))
(- (* (car v1) (cadr v2)) (* (cadr v1) (car v2)))
  )
)
(defun vk_Get3pNormal (p1 p2 p3 /)
;;;  (c:cal "nor ('p1','p2','p3')")
  (vk_GetOrt
    (vk_GetCrossProduct (mapcar '- p1 p2) (mapcar '- p3 p2))
  )
)
(defun vk_GetPlaneNormal (PlaneLst / p1 p2 p3)
  (setq p1 (car PlaneLst)
p2 (cadr PlaneLst)
PlaneLst (cddr PlaneLst)
  )
  (while (and (setq p3 (car PlaneLst)) (vk_IsPointOnLine p1 p3 p2))
    (setq PlaneLst (cdr PlaneLst))
  )
  (if p3
    (vk_Get3pNormal p1 p2 p3)
  )
)
(defun vk_TransPlane (PlaneLst N / pn)
  (setq pn (vk_GetPlaneNormal PlaneLst))
  (mapcar (function (lambda (c) (trans c N pn))) PlaneLst)
)
(defun vk_GetPointToPlaneDist (p PlaneLst / n)
;;;(c:cal "dpp ('p','p1','p2','p3')")
  (setq n (vk_GetPlaneNormal PlaneLst))
  (abs (- (last (trans (car PlaneLst) '(0 0 1) n))
 (last (trans p '(0 0 1) n))
       )
  )
)
(defun vk_GetFrustumVolume (BottomLst TopLst)
  ((lambda (s1 s2)
     (/ (* (vk_GetPointToPlaneDist (car TopLst) BottomLst)
  (+ s1 s2 (sqrt (* s1 s2)))
)
3.0
     )
   )
    (vk_GetArea (vk_TransPlane BottomLst '(0 0 1)))
    (if (cdr TopLst)
      (vk_GetArea (vk_TransPlane TopLst '(0 0 1)))
      0.0
    )
  )
)
(defun vk_GetPyramidVolume (BaseLst apex)
  (vk_GetFrustumVolume BaseLst (list apex))
)
(defun vk_GetPrismVolume (BottomLst TopLst)
  (vk_GetFrustumVolume BottomLst TopLst)
)
all functions are meant to be generic
still to be done: split the solid into three frustums. one of them will be a prism and two other - pyramids

[Kerry]

Taking the easy part first - volume of the middle chunk.

V_middle = A * L, where

L = min (P_Fz P_Gz P_Hz) - max ( P_Az P_Bz P_Cz)


< ..... >

This assumes that the body is vertical.

Is this a given for your problem.

Regards
Kerry

[dgorsman]

Taking the easy part first - volume of the middle chunk.

V_middle = A * L, where

L = min (P_Fz P_Gz P_Hz) - max ( P_Az P_Bz P_Cz)


< ..... >

This assumes that the body is vertical.

Is this a given for your problem.

Regards
Kerry

Yup.  From the first post, P_Axy = P_Fxy, P_Bxy = P_Gxy, and P_Cxy = P_Dxy.  If this wasn't a case, then an arbitrary coordinate system could be set up to make it so, with all initial points converted to that system.

[dgorsman]

Never fails - I get started on something and three other things pop up which need doing *now*.

Code: [Select]

;;
;; math_lib:volumeOfSolid4
;;
;; arguments: plane1 - list of three points (as list) describing the lower
;;   plane
;;   plane2 - list of three points (as list) describing the upper
;;   plane
;;
;; return value: real with shape volume
;;
;; comments: calculates the volume of shape defined by two triangular planes.
;; This function is intented to work in consistent metric units i.e.
;; if the drawing is in millimeters the volume will be in mm^3.
;;

(defun math_lib:volumeOfSolid4 ( plane1
                                       plane2 / return_volume
                                      vol_middle vol_lower vol_upper
                                      PA PB PC PF PG PH
                                      Iz1 Iz2
                                      PAM PBM PCM
                                      area_section )

;; Local defun: x value of a coordinate

   (defun |x| ( input ) (car input))

;; Local defun: y value of a coordinate

   (defun |y| ( input ) (cadr input))

;; Local defun: z value of a coordinate

   (defun |z| ( input ) (caddr input))

;; Local defun: angle from three points, as interior and exterior values

   (defun get< ( pt1 vertex pt2 / return_angle angle1 angle2 )
(if
         (<
            (setq return_angle (- (angle vertex pt1) (angle vertex pt2)))
            0
         )
      (setq return_angle (+ (* pi 2) return_angle))
   )

      (list
         (max (- (* pi 2) return_angle) return_angle)
         (min (- (* pi 2) return_angle) return_angle)
       )
)


   
;; Split up the arguments into the lower plane coordinates (PA PB PC) and the
;; upper plane coordinates (PF PG PH)

   (setq
      PA (car plane1)
      PB (cadr plane1)
      PC (caddr plane1)
   )

   (setq
      PF (car plane2)
      PG (cadr plane2)
      PH (caddr plane2)
   )

;; We break the solid into three parts - a triangular solid with parallel end
;; faces and two pyramids.

;; Step 1: getting the volume of the triangular solid

;; Construct a plane based on the max z of the lower plane

   (setq Iz1 (max (|z| PA) (|z| PB) (|z| PC)))
   (setq PAM (list (|x| PA) (|y| PA) Iz1))
   (setq PBM (list (|x| PB) (|y| PB) Iz1))
   (setq PCM (list (|x| PC) (|y| PC) Iz1))

;; Get the area of the end planes, using 0.5 * base * height, or
;; 0.5 * |PAM PBM| * (|PBM PCM| * sin < PAM-PBM-PCM)

   (setq
      area_section
        (*
           0.5
           (distance PAM PBM)
           (distance PBM PCM)
           (sin (cadr (get< PAM PBM PCM)))
        )
   )

;; Calculate the volume using the distance from the max z of the lower plane
;; and the min z of the upper plane

   (setq Iz2 (min (|z| PF) (|z| PG) (|z| PH)))

   (setq
      vol_middle
        (*
           area_section
           (distance
              PAM
              (list (|x| PF) (|y| PF) Iz2)
           )
        )
   )

;; Step 2: get the volume of the lower pyramid

;; This is just a placeholder for the moment
   
   (setq vol_lower 0)

;; Step 2: get the volume of the upper pyramid

;; This is just a placeholder for the moment
   
   (setq vol_upper 0)

   (+ vol_middle vol_lower vol_upper)
)

Apologies for the gibbled formatting - I'm a little tight on time.