Calculate volume using 3d Polyine

[cadgeek33]

Hi All,

I want to calculate cut volume for a pond. I have 2 Surfaces as 3d polylines with elevation in each vertices. Is there any lisp or routine I can follow rather than getting the volume using Civil 3D, because we have more than 500 similar ponds in this project. Sample dwg. and Image Attached for reference.

TIA.

[Sea-Haven]

Hi John the request was for volume of cut so look at my cross section image using water spill point would give incorrect answer this could be a substantial difference in volume. I know of a project where they ran out of dirt needed, the vols suggested an excess of material.

I designed a retard basin exactly that shape where batter chased the up hill slope.

The request is do multiple ponds, a couple of if and buts if all cut then make a super top surface all holes, make a second surface all top and bottom, take the top pline and get surface1 triangles, same again get triangles surface 2 then a VOL calc. The formula I posted adding each triangle down to say z=0 volume, then a VOL1-VOL2. This all falls apart when you have cut and fill.

Can use the ssget WP to find all the triangles. Need a answer for sample dwg correct volume. I got 5892. 

The code is rough as I am having problems with my CIV3d and need to do a full install again. I have checked against a simple pyramid solid and got the same answer for volume.

This is based on 2 surfaces 3dfaces from CIV3D it is the sample dwg C-TINN-VIEW is the 2 plines C-TINN-VIEW2 is the outside only.

; surface volumes 1st go
; the two surfaces must have a common outside boundary 

(defun trivols ( )
(setq vol 0.0 totarea 0.0)
(repeat (setq x (sslength ss))
(setq tri (vlax-ename->vla-object (ssname ss (setq x (- x 1)))))
(setq lst (vlax-get tri 'coordinates))
(setq j 1 x 1)
(repeat 3
(set (read (strcat "p" (rtos x 2 0) "x")) (nth (- j 1) lst))
(set (read (strcat "p" (rtos x 2 0) "y")) (nth j lst))
(set (read (strcat "p" (rtos x 2 0) "z")) (nth (+ j 1) lst))
(setq j (+ j 3) x (+ 1 x))
)
(setq lst2 '())
(setq lst2 (cons (list p1x p1y) lst2))
(setq lst2 (cons (list p2x p2y) lst2))
(setq lst2 (cons (list p3x p3y) lst2))
(setq lst2 (reverse lst2))
(setq lst2 (cons (last lst2) lst2))
(setq k 0 tot 0.0)
(repeat (- (length lst2) 1)
  (setq xy (- (* (car (nth k lst2))(cadr (nth (+ k 1) lst2))) (* (car (nth (+ k 1) lst2))(cadr (nth k lst2)))))
  (setq tot (+ tot xy))
  (setq k (+ k 1))
)
(setq tot (abs (/ tot 2.0)))
(setq newvol (* tot (/ (+ p1z p2z p3z) 3.0)))
(setq vol (+ vol newvol))
)
(princ)
)

(defun c:trivol ()
(setq ss (ssget (list (cons 0 "3DFACE")(cons 8 "C-TINN-VIEW"))))
(trivols)
(setq vol1 vol)
(setq ss (ssget (list (cons 0 "3DFACE")(cons 8 "C-TINN-VIEW2"))))
(trivols)
(alert (strcat "total volume is " (rtos  (- vol vol1) 2 2)))
(princ)
)

(c:trivol)

This is just a start so much more to add. 

[john.uhden]

@cadgeek33 ,

@Sea-Haven brought up an important question...

Are you looking for the total cut volume to excavate the pond or just the storage capacity of the pond?

If it's for the cut volume, I presume the upper 3Dpoly is drawn at existing grade, right?

[Sea-Haven]

Like John's comment we need that clarified, if all cut then provide a sample dwg with say 5 ponds need 2layers, layer1 has the 3d faces for the existing surface, layer2 has the bottoms added to the existing surface. The plines need to be on 2 layers as well to automate would look at outer layer and find the 3dfaces within it for the 2 surfaces.

Look in surfaces CIV3D there is a export option for the 3dfaces. Change the ouput of the triangle layers.

Test shape for idea created a solid then used top bottom lines with Z. Compared massprop to volume.

[ronjonp]

---

@john.uhden wrote:
Thanks, @hak_vz,
but I already have the interpolation function, and I don't like using
commands unless I have to. It's all about manipulating data.

---

Agreed 🍺

[cadgeek33]

yes, its for the cut volume, as you said upper line indicates the existing grade. Thanks

[hak_vz]

My current progress with the code to calculate volume from two base 3dpolylines using method of parallel vertical cuts (slices, sections). I'll try to make it universal so it can use as many objects in 3d.

[Sea-Haven]

Good answer fairly simple to make work, need a bounding box sort of answer  to interpret the best angle for the slices.

The more accurate answer would still intersect triangles look at ends. As the OP has CIV3D can make triangle faces. Then use the outer pline to select triangles. A iterative approach using ssget "F" so section line chooses correct 3dfaces. I did something drew a line over size and used intersectwith which returns 2 points ie edge of pline.

There is something weird in sample dwg as I can make 2 surfaces can see triangles and contours but they have a surface area of 0.0 no idea need more time to work out why not working. 

[john.uhden · Accepted solution]

I think I've got it.

Someone just needs to turn the 'items into a table.

(defun c:3DVOL ( / *error* Doc vars vals @Anonymous @group @closest @getpointatZ @enhance e obj plist1 plist2
                   dmax okay lines pts p1 p2 z zmin zmax done IncVol Vol A1 Z1 A2 Z2 data items title uline)
  ;; Program attempts to compute volumes from two (2) 3DPolylines
  ;;   via a specified elevation interval and the truncated prism method,
  ;;   or as @Sea-Haven (Alan Houston) insists the truncated pyramid method.
  ;; The Storage volume is computed from the lowest elevation of the bottom poly
  ;;   to the lowest elevation of the top poly.
  ;; The Cut volume is computed from the lowest elevation of the bottom poly
  ;;   to the highest elevation of the top poly.
  ;; Note that it draws the contours whether you want them or not
  ;;   in order to get the area of each.
  (gc)
  (vl-load-com)
  (prompt "\n3DVOL.lsp v1.0 (c)2021, John F. Uhden")
  (defun *error* (err)
    (mapcar 'setvar vars vals)
    (vla-endundomark Doc)
    (cond
      ((not err))
      ((wcmatch (strcase err) "*CANCEL*,*QUIT*")
         (vl-exit-with-error "\r                                              ")
      )
      (1 (vl-exit-with-error (strcat "\r*ERROR*: " err)))
    )
    (princ)
  )
  ;;----------------------------------------
  ;; Initialize
  ;;
  (setq Doc (vlax-get (vlax-get-acad-object) 'ActiveDocument))
  (vla-endundomark Doc)
  (vla-startundomark Doc)
  (setq vars '("cmdecho" "dimzin" "highlight"))
  (setq vals (mapcar 'getvar vars))
  (mapcar 'setvar vars '(0 1 0))
  (command "_.expert" (getvar "expert")) ;; dummy command
  (defun @Anonymous (p)(list (car p)(cadr p)))
  ;; Function to group a list of items into a list of
  ;; multiple lists, each of length N, e.g.
  ;; '(A B C D E F G H I) -> '((A B C)(D E F)(G H I))
  (defun @group (lst n / item new)
    (foreach element (reverse lst)
      (setq item (cons element item))
      (if (= (length item) n)
        (setq new (cons item new) item nil)
      )
    )
    new
  )
  (defun @closest (p plist / d dmin closest)
    (setq dmin 1000000.0)
    (foreach pt plist
      (if (< (setq d (distance p pt)) dmin)
        (setq dmin d closest pt)
      )
    )
    closest
  )
  (defun @getpoly (what / e obj plist)
    (and
      (setq e (car (entsel (strcat "\nSelect " what " 3DPoly: "))))
      (setq obj (vlax-ename->vla-object e))
      (or
        (= (vlax-get obj 'objectName) "AcDb3dPolyline")
        (prompt "  Object selected is NOT a 3DPoly.")
      )
      (or
        (not (zerop (vlax-get obj 'Closed)))
        (prompt "  3DPoly is NOT closed.")
      )
      (setq plist (@group (vlax-get obj 'Coordinates) 3))
      (setq plist (append plist (list (car plist))))
    )
    plist
  )
  (defun @getpointatZ (line z / p1 p2 z1 z2 m p)
    (setq p1 (car line)
          z1 (last p1)
          p2 (cadr line)
          z2 (last p2)
    )
    (cond
      ((= z z1)(setq p p1))
      ((= z z2)(setq p p2))
      ((or (< z1 z z2)(< z2 z z1))
        (setq m (/ (- z z1)(- z2 z1))
              p (mapcar '+ p1 (mapcar '* (list m m m)(mapcar '- p2 p1)))
        )
      )
    )
    p
  )
  (defun @enhance (old / i p1 p2 p d z1 z2 z n new)
    (setq i 0)
    (while (< i (1- (length old)))
      (setq p1 (nth i old)
	     i (1+ i)
	    p2 (nth i old)
            z1 (last p1)
	    z2 (last p2)
	     z z1
	     d (distance p1 p2)
	   new (append new (list p1))
      )
      (and
	(> d dmax)
	(setq n (1+ (read (rtos (/ d dmax) 2 0)))
	     dz (/ (- z2 z1) n)
        )
	(repeat n
	  (setq z (+ z dz)
	        p (@getpointatZ (list p1 p2) z)
	      new (append new (list p))
          )
        )
      )
    )
    (setq new (vl-remove nil new))
    (if (equal (last new)(last old))
      new
      (append new (list (last old)))
    )
  )
  (defun @makePline (plist elev layer closed / plist2)  ; closed must be 1 or 0
    (setq plist2 (mapcar '@2d plist)
          plist2 (mapcar '(lambda (x)(cons 10 x)) plist2)
    )
    (entmakex
      (append
        (list '(0 . "LWPOLYLINE")
              '(100 . "AcDbEntity")
              '(100 . "AcDbPolyline")
               (cons 8 layer)
               (cons 90 (length plist))
               (cons 70 closed)                    ; 1 for closed 0 otherwise
               (cons 38 elev)
             ; (cons 210 zdir)
        )
        plist2
      )
    )
  )
  (initget 1 "Cut Storage")
  (setq voltype (getkword "\nVolume type, Cut/Storage: "))
  (initget 7)
  (setq dmax (getdist "\nMaximum distance between tessellations: "))
  (and
    (setq plist1 (@enhance (@getpoly "bottom")))
    (setq plist2 (@enhance (@getpoly "top")))
   ; (foreach p plist1 (entmake (list '(0 . "POINT")(cons 10 p)'(8 . "3DPOINT"))))
   ; (foreach p plist2 (entmake (list '(0 . "POINT")(cons 10 p)'(8 . "3DPOINT"))))
    (setq zmin (apply 'min (mapcar 'last plist1))
          zmax (if (= voltype "Cut")(apply 'max (mapcar 'last plist2))(apply 'min (mapcar 'last plist2)))
          z zmin
          Vol 0.0
    )
    (while (not okay)
      (not (initget 7))
      (setq inc (getdist "\nVertical increment: "))
      (if (<= inc (- Zmax Zmin))
        (setq okay 1)
        (prompt "  Increment is too large.")
      )
    )
    (not (initget "Yes No"))
    (or (setq draw (getkword "\nDraw tessellations?, Yes/<No>: ")) 1)
    (or (setq draw (if (= draw "Yes") 1 nil)) 1)
    (foreach p2 plist2
      (setq p1 (@closest p2 plist1))
      (if draw (entmakex (list '(0 . "LINE")'(8 . "3DLINE")(cons 10 p1)(cons 11 p2))))
      (setq lines (cons (list p1 p2) lines))
    )
    (while (not done)
      (setq pts (vl-remove nil (mapcar '(lambda (x)(@getpointatZ x z)) lines)))
    ;  (print z)(print pts)
      (setq e (@makepline pts z "CONTOURS" 1))
      (setq area (vlax-get (vlax-ename->vla-object e) 'Area))
      (setq data (cons (cons z area) data))
      (cond
        ((= z zmax)(setq done 1))
        ((> (+ z inc) zmax)(setq z zmax))
        ((setq z (+ z inc)))
      )
    )
    (setq data (reverse data))
    (setq items '(("Elev" "Area" "Inc. Vol." "****. Vol.")))
    (setq items (cons (list (caar data)(cdar data) 0.0 0.0) items))
    (repeat (1- (length data))
      (setq z1 (caar data) z2 (caadr data) a1 (cdar data) a2 (cdadr data) data (cdr data))
      (setq IncVol (* (+ A1 A2 (sqrt (* A1 A2)))(/ (- z2 z1) 3.0)))
      (setq Vol (+ Vol IncVol))
      (setq items (cons (list z2 a2 IncVol Vol) Items)) 
    )
    (princ (setq title (strcat "\n" voltype " Volume Tabulation")))
    (setq uline "\n")(repeat (strlen title)(setq uline (strcat uline "-")))
    (princ uline)
    (mapcar 'print (setq items (reverse items)))
  )
  (*error* nil)
)

[hak_vz · Accepted solution]

@cadgeek33 

Here is my early version to calculate volume of the cut defined by two 3dpolyline using vertical sections. Select upper and lower base and two points on different side of lower base, let say two points on diagonal.

Because of some hiccups regarding area command on 3dpoly I'll probably have to make changes.

@john.uhden  Nice work!! I like it.

(defun c:cutvol( / 
adoc create_value_vector plane_line_intersection take pointlist3d normal_vector vector unit_vect vect_dot vect_cross
mod_vect create_value_vector acos asin angles_vect angles_two_vect poly vol upper lower eu el a b bb to_delete ang ang1 ang2 ang3 nsect ceu cel zlmin zumax i j p1 p2 p3 k pts_l pts_u int_pts tmp m n mp tp e an
to_delete ent f g d1 d2 step vtotal *error*
)
(defun *error* ( msg )
	(if (not (member msg '("Function cancelled" "quit / exit abort")))
		(princ)
	)
	(setvar 'cmdecho 1)
	(if (and adoc) (vla-endundomark adoc))
	(princ)
)
(defun create_value_vector (n val / r ) (repeat n (setq r (cons val r))) r)
(defun plane_line_intersection (pa pb pc la lb / m n l p d m1 m2 ip);
(setq
	n (unit_vect(normal_vector pa pb pc))
	m (vector pa la)
	l (vector lb la)
	p (* 1.0(vect_dot l n))
)
(if (/= p 0.0)(setq d (/ (* 1.0 (vect_dot m n)) p)))
(cond 
	((and d )
		(setq ip (mapcar '- la (mapcar '* l (list d d d))))
		(setq l (mapcar 'cadr (list ip la lb)))
		(setq m1 (car(vl-sort l '< )) m2 (car(vl-sort l '> )))
		(if (or (= m1 (cadr ip))(= m2 (cadr ip))) (setq ip nil)) 
	)
)
IP
)
(defun take (amount lst / ret)(repeat amount (setq ret (cons (car lst) (take (1- amount) (cdr lst))))))
(defun pointlist3d (lst / ret) (while lst (setq	ret (cons (take 3 lst) ret) lst (cdddr lst))) (reverse ret))
(defun normal_vector (p1 p2 p3) (vect_cross(vector p1 p2) (vector p1 p3)))
(defun vector (p1 p2) (mapcar '- p2 p1))
(defun unit_vect (v)(mapcar '* v (create_value_vector (length v) (/ 1 (mod_vect v)))))
(defun vect_dot (v1 v2)(apply '+ (mapcar '* v1 v2)))
(defun vect_cross (v1 v2 / a b c d e f)
(mapcar 'set '(a b c d e f) (append v1 v2))
(list 
	(-(* b f)(* c e))
	(* -1(-(* a f )(* c d)))
	(-(* a e)(* b d))
)
)
(defun mod_vect (v)(sqrt(vect_dot v v)))
(defun create_value_vector (n val / r ) (repeat n (setq r (cons val r))) r)

(defun acos (x) (cond ((and(>= x -1.0)(<= x 1.0)) (-(* pi 0.5) (asin x)))))
(defun asin (x) (cond ((and(> x -1.0)(< x 1.0)) (atan (/ x (sqrt (- 1.0 (* x x)))))) ((= x -1.0) (* -1.0 (/ pi 2))) ((= x  1) (/ pi 2)) ) )
(defun angles_vect (v)(mapcar 'acos (unit_vect v)))
(defun angles_two_vect (v1 v2 / ang xa ya za xb yb zb) 
	(mapcar 'set '(xa ya za) v1)
	(mapcar 'set '(xb yb zb) v2)
	(setq ang 
		(acos 
			(/
				(+ (* xa xb)(* ya yb)(* za zb)) 
				(* (sqrt (+(* xa xa)(* ya ya)(* za za)))(sqrt (+(* xb xb)(* yb yb)(* zb zb))))
			)
		)
	)
	(if (< zb za) (setq ang (- (* pi 2.0) ang))) ; to test angle relative to horizontal vector in plane givs angle in full circle
	ang
)

(defun poly (lst)
(entmake '((0 . "POLYLINE")(70 . 9)))
(foreach l lst
(setq  k (append '((10)) l))
(entmake (list(cons 0 "VERTEX")(cons 70 32) (apply 'append k))) 
)
(entmake '((0 . "SEQEND")))
)
(defun vol (a b d ) (setq a (sqrt a) b (sqrt b))  (* (/ d 3.0)(+ (* a a)(* b b)(* a b))))
    (setq adoc (vla-get-activedocument (vlax-get-acad-object))) 
	(setq upper (car(entsel "\nSelect upper base >")))
	(setq lower (car(entsel "\nSelect lower base >")))
	(setq eu (vlax-ename->vla-object upper))
	(setq el (vlax-ename->vla-object lower))
	(setq a (take 2(getpoint "\nSelect first section point on smaller base >")))
	(setq b (take 2(getpoint "\nSelect second section point on smaler base >")))
	(setq bb b)
	(setq to_delete (ssadd))
	(setq area_list (list))
	(setq ang1 (angle a b) ang2(- ang1 (/ pi 2)) ang3(+ ang1 (/ pi 2)))
	(setq a (polar a ang1 0.01))
	(setq b (polar b (angle b a) 0.01))
	(setq nsect 100)
	(setq step (/ (distance a b) (+ nsect 1)))
	(setq ceu (pointlist3d (vlax-get eu 'Coordinates)))
	(setq cel (pointlist3d (vlax-get el 'Coordinates)))
	(setq zumax (+ (apply 'max(mapcar 'last ceu)) 1.0))
	(setq zlmin (- (apply 'min(mapcar 'last cel)) 1.0))
	(setq i -1)
	(setvar 'cmdecho 0)
	(vla-endundomark adoc)
	(vla-startundomark adoc)
	(while (< (setq i (1+ i)) (1+ nsect))
		(setq b (polar a ang1 (* i step)))
		(setq p1 (append (polar b ang3 10.0) (list zlmin)))
		(setq p2 (append (polar b ang2 10.0) (list zlmin)))
		(setq p3 (append (polar b ang2 10.0) (list zumax)))
		;(command "_.3dpoly" "_non" p1 "_non" p2 "_non" p3 "")
		(setq j -1 k (length ceu) pts_u nil pts_l nil)
		(setq ceu (append ceu (list (car ceu))))
		(while (< (setq j (1+ j)) k)
			(if (and (setq m (plane_line_intersection p1 p2 p3 (nth j ceu)(nth (1+ j) ceu))))
				(setq pts_u (cons m pts_u))
			)
		)
		(setq j -1 k (length cel))
		(setq cel (append cel (list (car cel))))
		(while (< (setq j (1+ j)) k)
			(if (and (setq m (plane_line_intersection p1 p2 p3 (nth j cel)(nth (1+ j) cel))))
				(setq pts_l (cons m pts_l))
			)
		)
		(setq int_pts (append pts_l pts_u))
		(cond 
			((and int_pts (>= (length int_pts) 3))
				(setq tmp int_pts m (car tmp) tmp (cdr tmp))
				(while tmp
					(setq m (mapcar '+ m (car tmp)) tmp (cdr tmp))
				)
				(setq mp (mapcar '/ m  (create_value_vector (length int_pts)(length int_pts))))
				(setq tp (polar mp ang2 5.0))
				(setq e -1)
				(while (< (setq e (1+ e)) (length int_pts))
					(setq an (angles_two_vect (vector mp tp)(vector mp (nth e int_pts))))
					(setq tmp (cons (list an (nth e int_pts)) tmp))
				)
				(poly(mapcar 'cdr(vl-sort tmp '(lambda (a b) (< (car a)(car b))))))
				
				(setq to_delete (ssadd (entlast) to_delete))
				(command "_.copy" (entlast) "" '(0 0 0)'(0 0 0))
				(command "_.rotate3d" (entlast) "" "2p" "_non" p1 "_non" p2 90.0)
				(command "_.region" (entlast) "")
				(setq area_list (append area_list (list (vlax-get (vlax-ename->vla-object (entlast)) 'Area))))
				(entdel (entlast))
			)
		)		
	)
	(command "_.line" "_non" a "_non" bb "")
	(command "_.extend" upper "" (entlast) "_non" a "_non" bb "")
	(setq ent (entget (entlast)))
	(setq f (take 2(cdr (assoc 10 ent))) g (take 2 (cdr (assoc 11 ent))))
	(setq d1 (distance a f) d2 (distance b g))
	(entdel (entlast))
	(setq vtotal (vol 0.0 (car area_list) d1))
	(setq vtotal (+ vtotal (vol 0.0 (last area_list) d2)))
	(setq i -1)
	(while (< (setq i (1+ i)) (- (length area_list)1))
		(setq vtotal (+ vtotal (vol (nth i area_list)(nth (1+ i) area_list) step)))
	)
	(setvar 'cmdecho 1)	
	(vla-endundomark adoc)
	(command "_.erase" to_delete "")
	(princ (strcat "\nAproximate volume = " (rtos vtotal 2 2)))
	(princ)
)

[john.uhden]

@hak_vz 

I am honored to receive compliments from you.

I am equally impressed with your vector manipulation and the slicing show it puts on.

However, I drew a couple of 3DPolys and got very different answers...

CUTVOL = 2991

3DVOL = 3378

Yours appeared to start taking slices a ways in from the upper poly, so I'm wondering if it lost some volume there.

Wait a second... my first run was with an elevation interval of 0.2.  I ran it again using 0.5 and got a cut volume of 3360.  I should add a filter so that it can be run with a tight interval but report at say every 1 unit of depth.

[cadgeek33]

@john.uhden @hak_vz 

Thanks alot guys, finally got what I wanted. Much appreciated. Keep up the good work.