Water drop path over TIN

[autoid374ceb4990]

I am not sure if this is the correct forum for this post, perhaps I should in the ObjectARX or Civil3D forum, but here goes. I am attempting to implement a water drop function on a TIN surface model. I know that the newer AutoCAD versions have such a function, but from what I have read on the Civil3D forum it has several shortfalls plus I am still using R14 and I have written my own software for contouring, earthwork, and other civil engineering and land surveying applications.

My question is does anyone know (in general) the approach that the AutoCad software ( or any other similar software) uses to implement the water drop function? Not looking for source code, just an idea of the methodology used to trace the path of the water across the surface of a TIN.

Thanks,

Charles

[leeminardi · Accepted solution]

A TIN can be considered as a collection of 3DFACES.  You can determine the maximum downward slope (the direction a drop of water will travel) as follows:

1. Give  a 3DFace with vertices A, B, and C and a point P on the face determine the normal vector to the face:
   nrm = (B-A) cross (C-A)
2. Check that the normal has a positive Z component:
   if ( nrm dot [0,0,1] < 0 then multiply the nrm by -1.
3. Take the cross product of the normal and the z direction
   v = nrm cross [0,0,1]
4. Take the cross product of the normal and v
   w = nrm cross v
   w is the vector pointing in the direction of the water drop travel.
5. Find the point Q on the edge of the face where the drop leaves by looking at the three intersection points of a line going from p in the direction w. Note, two of the intersection points wil be on the finite line segment AB, BC or CA.  Disregard the intersection point that is not on the line segment and choose the intersection point with the lowest z coordinate of the other two. 
   Q1 = intersection of line p (p + w) and A B
   Q2 = intersection of line p (p + w) and B C
   Q3 = intersection of line p (p + w) and C A
6. A line from P to Q is the path of the water drop
7. Use Q one the adjacent face as the P for the next face and repeat!

[autoid374ceb4990]

Leemardi:

Thanks for your explanation, exactly what I was looking for.

Have you implemented this solution yourself?

[leeminardi]

Years ago I wrote a LISP program to create a geodesic curve given an array of 3dfaces. It had similar requirements to the water drop problem.  I'll see what I can put together.

[leeminardi]

@autoid374ceb4990 Here's a program that will project a give point vertically to a selected 3dface and determine the downward slope point of exit for a drop of water. The program determines if the projected given point is within the bounds of the face. It igores faces that are vertical although I am not sure this is a needed limiation.
To test it create a 3DFACE that is not parallel to the WCS XY plane, go to the top view and specify a point.

The next step is to enable the selection of multiple 3dfaces and compute the slope lines from one face to the next. Is this something you would use?

; Calculates the point on a 3Dface vertically aboce a given point and the
; the point on its edge downslope from the point on the 3dface.
; L. Minardi  5/20/2024

(defun c:FaceSlope (/ ss en edata pxy p A B C N Apxy verticalFace tk x1 x2 x3 qlist w q1 q2 q3 q)
  (princ "\nPlease select 3DFACE and press ENTER.")
  (setq	ss    (ssget '((0 . "3DFACE") ))			; get data for 3Dface
	en    (ssname ss 0)
	edata (entget en)
  )
  (setq	pxy   (getpoint "\nEnter point under face:")	; a point
	p   (mapcar '+ pxy '(0 0 1))  ; point above pxy
  )
  (setq	A (cdr (assoc 10 edata))	;three corners of the 3DFACE
	B (cdr (assoc 11 edata))
	C (cdr (assoc 12 edata))
  )
;;;  (setq	AB (mapcar '- B A)	;define vectors along edges of face
;;;	BC (mapcar '- C B)
;;;  )
  (setq N (cross (mapcar '- B A) (mapcar '- C B)))		; normal to face
  (setq	Apxy (mapcar '- pxy A)	; vector from face to pxy
  )
;;;; determine if face is vertical
(if (< (abs (caddr n)) 0.00001)
  (setq verticalFace T)
  (setq verticalFace nil)
)
  (if verticalFace
    (princ
      "\nVertical face."
    )
    (progn ; face not vertical, continue
      (setq tk (- (/ (dot N Apxy) (dot N '(0 0 1))))
	; value of parameter t at intersection
      )		; intersection point of line and face using parametric definition of a line  
      (setq P (mapcar '+ ; redefine p as point of face 
			 pxy
			 (mapcar '* (mapcar '- p pxy) (list tk tk tk))
		 )
      )
      (command "point" "_non" P)
      ; determine if point is inside face boundary
(setq x1 (dot N (cross (mapcar '- P A) AB))
      x2 (dot N (cross (mapcar '- P B) BC))
      x3 (dot N (cross (mapcar '- P C) CA))
)
(if  ;test inside
  (or
    (and (>= x1 0) (>= x2 0) (>= x3 0))
    (and (<= x1 0) (<= x2 0) (<= x3 0))
  )
(progn
  ;(princ "\nThe intersection point is INSIDE the Face")
  (setq qlist nil) ; list of valid intersections of slope line from p
  (setq w (cross N (cross N '(0 0 1)))) ; downslope vector
  (setq q1 (inters A B P (mapcar '+ P W) nil))
  (if (onSegment A B q1) 
    (setq qlist (append qlist (list q1)))
  )
  (setq q2 (inters B C P (mapcar '+ P W) nil))
  (if (onSegment B C q2)
    (setq qlist (append qlist (list q2)))
  )
  (setq q3 (inters C A P (mapcar '+ P W) nil))
  (if (onSegment C A q3)
    (setq qlist (append qlist (list q3)))
  )
; which intersection point has the lowest z coordinate
  (if (< (caddr (car qlist)) (caddr (cadr qlist)))
	 (setq q (car qlist))
    (setq q (cadr qlist))) ; point on face edge downslope from p
  (command "_point" "_non" q)
  (command "_line" "_non" p "_non" q "")
)					; end progn	
   (princ "\nThe intersection point is OUTSIDE the Face")
)  ;end if test inside outside
)
)	; end if vertical face

  (princ)
)					;end Face-Intr



; Determines if the projection of point pt is on the
; line segment from point p1 to p2. True if it is,
; nil if it isn't.
(defun onSegment (p1 p2 pt / x v ti)
  (setq x (distance p1 p2))
  (setq v (mapcar '/ (mapcar '- p2 p1) (list x x x)))
  (setq	ti (/ (dot (mapcar '- pt p1) v)
	      (dot (mapcar '- p2 p1) v)
	   )
  )
  (if (and (>= ti 0) (<= ti 1.))
    (setq x T)
    (setq x nil)
  )
)
  
;;; Compute the cross product of vectors a and b

(defun cross (a b / crs)
  (setq	crs (list
	      (- (* (nth 1 a) (nth 2 b))
		 (* (nth 1 b) (nth 2 a))
	      )
	      (- (* (nth 0 b) (nth 2 a))
		 (* (nth 0 a) (nth 2 b))
	      )
	      (- (* (nth 0 a) (nth 1 b))
		 (* (nth 0 b) (nth 1 a))
	      )
	    )				;end list
  )					;end setq c
)					;end cross
;;; Compute the dot product of 2 vectors a and b
(defun dot (a b / dd)
  (setq dd (mapcar '* a b))
  (setq dd (+ (nth 0 dd) (nth 1 dd) (nth 2 dd)))
)					;end of dot
; Unit vector of v
(defun unitvec (v / x)
  (setq	x (distance '(0 0 0) v)
	x (mapcar '/ v (list x x x))
  )
)

[autoid374ceb4990]

Lee:

Thanks for the LISP code.  When I load it in to AutoCAD I get an error: "ACAD error: invalid character".  I found some LISP code in the AutoLISP forum  that purports to do something similar to yours, but it gave me the same error.

Perhaps this error is because I am using an ancient (R14) version.  Most of my code is written in "C" and I haven't written any LISP code for many years, so I will have to try to convert your code to "C" .  I already have functions for manipulating the 3DFACE's in a TIN, like calculating the X,Y,Z coordinates by picking a random point, plotting the profile of a POLYLINE path across the TIN,  and calculating the volumes between TIN's.  I have "C" code to calculate the face normal vector of a a given 3DFACE, so hopefully all I need is the cross product code to calculate the water drop path.

Thanks again,

Charles

[leeminardi]

Charles,

I do not know what is causing the invaid character problem.

Attached is a copy of the code that I created by copyling it to Notepad which should only allow standard ASCII characters.  See if it will load.

Please let me know if you have anyquesions about my code.

Lee

[autoid374ceb4990]

Well loading the code from NOTEPAD cured part of the problem.  I had been using WORDPAD and when I load the code in WORDPAD, save it, and run the code, I get the same "invalid character"error???.   Anyway with the WORDPAD version I am able to  select a 3DFACE and when I pick a point I get a POINT drawn at the pick point and then: ACAD error: invalid argument

[leeminardi]

Well, at the last minute I decided to do some editing to clean up the code.  I thought I had checked after I did.  Apparently I was too agressive in my edits and careless in not checking the results.

Please try this version.

[autoid374ceb4990]

That seems to have fixed it.

Thanks,

Charles

[autoid374ceb4990]

"The next step is to enable the selection of multiple 3dfaces and compute the slope lines from one face to the next. Is this something you would use?"

Yes, having the program automatically move from face to face is my ultimate goal.  I have some TIN's created from LIDAR points that are 500,000 to 700,000 faces with many triangles less than 1 sq.ft. in area.

[leeminardi]

The main challenge in implementing a routine to trace the path of a drop of water flowing dow a TIN is making the search for the adjacent face effieicient.   If I were to write the code it would just do a linear search trhough all the faces until one is found that shares the edge with the know edge point Q. With hundreds of thousands of faces this could take a while to compute.  Also, it may be more computationally efficient to work with a TIN than a file with individual faces.  

Can you post a TIN with a few hundred triangles I could use as a test?

Have you looked at Civil 3D.  It has watershed analysis functions.

https://help.autodesk.com/view/CIV3D/2022/ENU/?guid=GUID-E021850D-0749-42CC-A49E-6C8BB6CC1DC8

[autoid374ceb4990]

I thought about the “locate the adjacent triangle” problem also. My TIN’s are saved in a binary file and I am loading all the triangles into memory to make searches faster. I also have a function that locates a triangle given any X,Y point, so I can pick a point and instantaneously load the correct triangle and calculate the elevation on the triangle surface. So I thought I would be able to locate the adjacent triangle by first locating the downhill point of the trace line on the first triangle, then extend the slope line by say 0.1 unit, calculate the coordinates of the extended point, and use those X,Y coordinates to identify the adjacent triangle.

Using LISP I thought you might be able to use the downhill point and select the two 3DFACES that cross that point, filter out the uphill face, and use the downhill face to continue the trace.

Using either method my require some filtering in case the downhill point lands on a triangle corner or a flat face.

I have attached a TIN drawing file for you to play with. It has about 1000 3DFACE's and also some contour lines.  I cut the data from a larger drawing that was created from LIDAR points.  It is an R14 drawing.  Let me know if you can read it.

I have not tried Civil3D since I can do everything I need to do in R14 plus I do not care for AutoDesk's "rent a program" policy and I did not like the idea of having to buy a new C compiler every time AutoCAD was updated.  I am also old and retired so I just program to keep my mind working.

[leeminardi]

If the point P is on an edge of a face but not at a corner then a simple search for another face that shares that edge can be done by comparing the three vertices of the face under examnation with the vertces of the edge containing P.  It gets more complcated if P is at a vertex of a face.  Then an examination of the slope for each of the faces that share that vertex has to be made.  Messy! 

For example.

Note that starting at a random point 1 yields the following flow.  The two faces at point 2 have an upward slope so the water flows down the edge to point 3. The slope of the four faces that share point 3 (two can be ignored) have to be compared before moving on.

[autoid374ceb4990]

Where is point P?

I think the major question is: How do you do you select the adjacent (next) triangle?

I would think that if you reach a point where the trace is running uphill or level you could just stop searching.

[leeminardi]

Point P is initially at #1 on the surface of the 3DFACE. Where its path leaves the face is point Q.  Point Q becomes point P for the adjacent face and a new Q is calculated.  The process continues until there's nowhere lower to go from the last point P.  As noted, the adjacet faces for when P is at #2 slope upwards, the path from there follow the edge (crevice) to the corner at #3.