Help Needed lisp Routine for Creating 3D Polylines with Specified Slope"

[libarraSCXGZ]

Hi everyone,

I need to develop a LISP routine that allows creating a 3D polyline with a user-specified slope. The goal is for the user to select a starting point, and then the polyline should continue along the surface, following the slope until it reaches a point where the maximum slope at that point is less than the design slope.

I have to carry out feasibility studies for roads and irrigation canals, and it is tedious to do it manually.

Here is a reference of what I need to do:

https://filesnj.carlsonsw.com/mirror/manuals/Carlson_2007/online/index.html?page=source%2FSite_Road_...

Could you help me with some basic LISP code to help me study and move forward with this? My knowledge of LISP is basic, but any help would be appreciated.

[leeminardi]

@autoid374ceb4990 

The attached version 2 creates a 3dpoly of the path and fixes the user input of the design slope.

The program will occasionally prematurely stop due to it not selecting both 3daces that share the point Q. I tried  specifying Q twice and creating Lower Left and Uppe right points slightly offset from Q.  E.g.,

  (setq pLL (mapcar '- Q '(0.002 0.002 0.002))) ; Lower Left of select cross
  (setq pUR (mapcar '+ Q '(0.002 0.002 0.002))) ; Upper Right of select cross
  (setq ss (ssget "_c" pur pLL '((0 . "3dface"))))

  This did not help!  When a selection set does not include the 2 faces an error message is output and a small circle is drawn at Q.  The user can reissue the commans and use the end of the 3dpoly as the start for a new path.

When time allows I will write a routine to search faces and look to see if Q is on one of its edges.  I suspect that this will slow down the calculation significantly.

[autoid374ceb4990]

@leeminardi

When I ran your latest code I initially had a "bad argument type" error.  I traced it to the 3dpoly command.  I changed your line "(while (= (getvar 'cmdactive) 1)" to "(while (= (getvar "cmdactive") 1)" and that fixed that problem.  Then I noticed that the 3DPOLY was following the triangle sides, because I had a snap mode active, so I added "(setvar "OSMODE" 0)" just after the "(command "_3dpoly) line and that fixed it.

Nice work! 

[leeminardi]

@autoid374ceb4990 

Thanks for doing some debugging of my code.  Versionn 3 (attached) has the changes you made although instead of toggling "osmode"  I added "_non" in front of points.  I prefer doing this rather than storing crrent osmode settings and then reinstating them.

It's clear to me that some additonal debugging is required in addition to a better method for finding the next 3dface.

(defun c:SpecSlope (/)  ; needs list of variables here.
; Creates a path along a collection of 3d faces of user defined slope.
; The path is defined by individual line segments.
; User input:
;   Design slope, default = 0.2 (20%)
;   Point on edge of a 3dface
;   Selection of the 3dface with the point
; Fixes in v.02
;   1. Fixed bug that restricted slope to 0.2
;   2.  A 3dpoly is created.  
; Known problems:
;   1. Sometimes fails because selecting cross at q does not find both faces.
;   2. Does not handle case of path leading to the vertex of several 3dfaces.
;   3. User control of which of 2 valid directions to take is not provided.
; L. Minardi 8/31/2024  v3
(setq designSlope (getreal "\nEnter slope, (e.g. .2 = 20%): "))
(if (not designslope)
  (setq designslope 0.2)
)
(setq p (getpoint "\nEnter start point on the edge of a face. "))
(command "point" "_non" p)
(princ "\n Select the 3dface of the stating point: ")
(setq ss (ssget '((0 . "3DFACE"))) ; get data for 3Dface
      faceName (ssname ss 0)
      Qlist (list p)
)

(while faceName

  (setq	faceData    (get_face_data faceName)
	A	    (nth 0 faceData)
	B	    (nth 1 faceData)
	C	    (nth 2 faceData)
	Normal	    (nth 3 faceData)
	slopeVector (nth 4 faceData)
	slopeFace   (nth 5 faceData)
	err	    (nth 6 faceData)
  )
  (if (= err 1)
		(progn
		  (princ  "\nNoSolution, 3dFace slope is less than maximum slope.\n
				   Slope  = "
		  )
		  (princ slopeface)
		  (setq faceName nil)
		)
(progn
  (setq gamma (acos (/ designSlope slopeFace)))
  (setq q (FindQ p A B C slopeVector designSlope gamma))

  (if (< (caddr q) (caddr p))		; is q lower
    (progn				; try minus gamma
      (setq gamma (- gamma))
      (setq q (FindQ p A B C slopeVector designSlope gamma))
    )
  )
;  (command "_point" "_non" Q)
 ; (command "_line" "_non" P "_non" Q "")
  (setq qList (append qList (list q)))
;;;     (setq ss (ssget "_c" q q '((0 . "3dface"))))

;;;  (setq Q2D (list (car Q) (cadr Q)))
;;;  (setq pLL (mapcar '- Q2D '(0.005 0.005)))
;;;  (setq pUR (mapcar '+ Q2D '(0.005 0.005)))
  (setq pLL (mapcar '- Q '(0.002 0.002 0.002))) ; Lower Left of select cross
  (setq pUR (mapcar '+ Q '(0.002 0.002 0.002))) ; Upper Right of select cross
  (setq ss (ssget "_c" pur pLL '((0 . "3dface"))))
;;(setq ss (ssget "_c" q q '((0 . "3dface"))))
  (if (not ss)
    (progn
      (princ "\nNo selection at Q.")
      (command "_circle" "_non" q 0.2)		; draw circle to show error location
      (setq facename nil)		; to terminate while
    )
    (progn
      (setq face1 (ssname ss 0)
	    face2 (ssname ss 1)
      )

      (if (equal face1 facename)
	(setq faceName face2)
	(setq faceName face1)
      )
      (setq p q)
    )					; end progn no errors
  )					; end if
)					; end progn, error not 1
    )  ; end if
  ); end while
 (command "_3dpoly"  )			;start 3dpoly command
  (while (= (getvar "cmdactive") 1)	;while a command is running loop (it is because we haven't entered any points)
    (repeat (setq x (length qList))	
      (command "_non" (nth (setq x (- x 1)) qList))
    )					;end repeat
    (command "")			;end 3dpoly
  )


    
    (princ)
)


(defun FindQ ( p v1 v2 v3 slopeVector designSlope gamma / q s q1 q2 q3 )
; determine up hill point at slope = designSlope
; p = start point
; v1, v2, v3 = 3dface vertices
; slopeVector = slope vector of 3dface that points in the uphill direction
; designSlope = desired slope of line  (e.g., 0.2 = 20%)
; gamma = angle between the projection of slopeVector to xy plane and the
;         direction of the designSlope line projected to the XY plane  
  (setq	Q
	 (mapcar
	   '+
	   p
	   (list
	     (-
	       (* (car slopevector) (cos gamma))
	       (* (cadr slopevector) (sin gamma))
	     )
	     (+
	       (* (car slopevector) (sin gamma))
	       (* (cadr slopevector) (cos gamma))
	     )
	     (*	designSlope
		(sqrt
		  (+ (expt (car slopevector) 2) (expt (cadr slopevector) 2))
		)
	     )
	   )
	 )
  )
; get the three intersections of design slope vector with face edges
(IF (equal p q 0.0001)
  (setq s (mapcar '+ p slopevector))
(setq s (mapcar '+ p (mapcar '* (mapcar '- q p) '(1000 1000 1000))))
)					; point in direction of slopevector from p
(setq q1 (inters v1 v2 p s nil)
      q2 (inters v2 v3 p s nil)
      q3 (inters v3 v1 p s nil)
)
(princ "\n q1 = ")
(princ q1   )
(princ "\n q2 = ")
(princ q2   )
(princ "\n q3 = ")
(princ q3   )
  
; check if q1, q2, q3, are on an edge segment of face
; set q to nil if not on line segment
  (if (not (equal (distance v1 v2)
		  (+ (distance v1 q1) (distance q1 v2))
		  1e-5
	   )
      )	; q1 on line segment?
    (setq q1 nil)
  )
  (if (not (equal (distance v2 v3)
		  (+ (distance v2 q2) (distance q2 v3))
		  1e-5
	   )
      )					; q2 on line segment?
    (setq q2 nil)
  )
  (if (not (equal (distance v3 v1)
		  (+ (distance v3 q3) (distance q3 v1))
		  1e-5
	   )
      )					; q3 on line segment?
    (setq q3 nil)
  )
					;set q to nil if q coincident with p
  (if (equal q1 p 0.001)
    (setq q1 nil)
  )
  (if (equal q2 p 0.001)
    (setq q2 nil)
  )
  (if (equal q3 p 0.001)
    (setq q3 nil)
  )
					;  set q to only non-nil q
  (cond
    (q1 (setq q q1))
    (q2 (setq q q2))
    (q3 (setq q q3))
  )

)


;;;;;;;;;;;;;;;;;;;
(defun get_face_data (faceName / A B C Normal slopeVector data)
; Given a face name determine its vertices, normal vector, slope vector, slope, and error if any.
; Returns a list with the data
; err = nil, no errors
; err = 1, no solution, the slope of the face is less than the design slope  
  (setq faceData (entget faceName))
  (setq	A      (cdr (assoc 10 faceData)) ;three corners of the 3DFACE
	B      (cdr (assoc 11 faceData))
	C      (cdr (assoc 12 faceData))
	Normal (cross (mapcar '- B A) (mapcar '- C B))
	err    nil
  )
  (setq slopeVector (cross Normal (cross Normal '(0 0 1))))
  (if (< (caddr slopeVector) 0.)	; make sure normal has a + z coorindate
    (setq slopeVector (mapcar '* '(-1 -1 -1) slopeVector))
  )
  (setq	slopeFace
	 (/
	   (caddr slopeVector)
	   (distance '(0 0 0)
		     (list (car slopeVector) (cadr slopeVector) 0.0)
	   )
	 )
  )
  (if (< slopeface designSlope)
    (setq err 1)
  )
  (setq data (list A B C Normal slopeVector slopeFace err))
)


; Unit vector of v
(defun unitvec (v / x)
  (setq	x (distance '(0 0 0) v)
	x (mapcar '/ v (list x x x))
  )
)

;; ArcCosine  -  Lee Mac
;; Args: -1 <= x <= 1

(defun acos ( x )
    (if (<= -1.0 x 1.0)
        (atan (sqrt (- 1.0 (* x x))) x)
    )
)
(defun dot (a b / dd)
  (setq dd (mapcar '* a b))
  (setq dd (+ (nth 0 dd) (nth 1 dd) (nth 2 dd)))
);end of dot

;;; 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


;; Tangent  -  Lee Mac
;; Args: x - real

(defun tan ( x )
    (if (not (equal 0.0 (cos x) 1e-10))
        (/ (sin x) (cos x))
    )
)
; test if 3 points are collinear. 
;; Collinear-p  -  Lee Mac
;; Returns T if p1,p2,p3 are collinear
(defun Collinear-p (p1 p2 p3)
  (
   (lambda (a b c)
     (or
       (equal (+ a b) c 1e-8)
       (equal (+ b c) a 1e-8)
       (equal (+ c a) b 1e-8)
     )
   )
    (distance p1 p2)
    (distance p2 p3)
    (distance p1 p3)
  )
)
; some function to help debugginh
(defun c:faceName (/)
  (setq	ss	 (ssget '((0 . "3DFACE"))) ; get data for 3Dface
	faceName (ssname ss 0)
  )
  (entget faceName)
  (princ facename)
  (princ)
)

(defun c:GetSlope (/ s p1 p2 data d slope)
  (setq	s     (ssget '((0 . "line")))
	data  (entget (ssname s 0))
	p1    (cdr (assoc 10 data))
	p2    (cdr (assoc 11 data))
	d     (distance	(list (car p1) (cadr p1))
			(list (car p2) (cadr p2))
	      )
	slope (/ (- (caddr p2) (caddr p1)) d)
  )
  (princ "\nSlope = ")
  (princ slope)
  (princ)
)

[john.uhden]

@leeminardi ,

I'm just not getting this.

Using AutoLisp (without vector math) I can find my way down a slope (the old water drop path), but how the hello do you specify a slope if a face doesn't have that slope anywhere across it, like if it's flat?

[john.uhden]

@leeminardi , @autoid374ceb4990 ,

I'm getting the idea that it's better to use my profile widget because there are too may dead ends trying to force existing to match your desires.

Anyway, it looks to me as though he is trying to find the best path for a grand slalom run.  But I'm afraid it will require some land grading to achieve the desired slopes.  I think I can speed up my elevation detection if that process (Working >>>>>) is too big a speed bump (like missing a gate and flying off the course, which Lee might enjoy  but not Lindsay Von).

[leeminardi]

@john.uhden 

I'm not sure I understand your question.  Are you asking:

1.  How to find the direction on a face witha specified slope?

or 

2. How to find the slope of a 3dface without using vector math?

If the latter then look at the following:

ABC  are the vertices, pAB the intersection of AB with the gorund plane, pAC ...

Point D is the projection of A onto the line from pAB to pAC.

The slope = AE/DE 

Remember, when working in 3D, vectors are your friend!

Lee

[john.uhden]

@leeminardi ,

That graphic you made is very revealing.  Maybe it will help me to learn.

BUT, my main concern/question remains how can you find where a 20% slope is headed when the 3Dface (plane) is at a slope of say 10%?!

[leeminardi]

---

@john.uhden wrote:

@leeminardi ,

That graphic you made is very revealing.  Maybe it will help me to learn.

BUT, my main concern/question remains how can you find where a 20% slope is headed when the 3Dface (plane) is at a slope of say 10%?!

---

@john.uhden 

If the slope of a 3dface is 10% there is no direction you can go on that 3dface where the slope will be 20%.  The task is to find a less steep slope on a more steep sloped hill.  E.g., find a direction on a 20% sloped incline that has a 10% slope.

I tried to explains how I do it in my post #53.

First I find the slope of the 3dface.  I named it slopeVector.  You've implied that vectors freak you out but they are not that complex and a spatial guy like you should like them.  Consider SlopeVector to be a arrow pointing uphill with the same slope as the hill. The desired slope that we want to find is DesignSlope and can be show as another arrow (although as specified, it's a scalar)!  We can define this arrow (vector) as being below the SlopeVector. Note, (and this may be a source for confusion) DesignSlope in the drawing below implies a vector but in my code it is a scalar (e.g., 0.0 for a 20% slope). Its rise/run is equal to the user specified design slope. But, as drawn, the DesignSlope arrow is below the surface of the hill (3dface). That is, its not as steep as the slope of the hill. We want to rotate this vector about the world Z axis such that it becomes flush wth the 3dface.  I named this angle gamma and derived the following expression for gamma:

gamma = acos ( (design_slope) / (slope_of_the_3DFace))

where:

   design_slope = user specifed

   slope_of_the_3DFace = ( z coordinate of the slope vector) /

                                            (the length of the slope vector projected to the XY plane)

Let me know if you'd like an explanation of the derivation of gamma.

There are two possible solutions for the direction of the less steep path up the hill.  One is pointing uphill to the left of going straight up the hill (the direction of slopeVector) and the other pointing to the right of going straight up the hill.

Here are top views looking down at the 3dface showing plu and minus values for gamma.  The "Q" points show potential exit points from the 3dface at the design slope.

Please review the latest revision of the code in post #63 and don't hesitate to ask if you have more questions .

Lee

[autoid374ceb4990]

When I try your version 3 (and previous versions) I find that the search for the next adjoining 3DFACE stops when the “ssget _c " reaches the edge of the screen. I assume that this is because the ssget requires the points to be visible on the screen. If I zoom out and try the same start point and 3DFACE the polyline goes further but still stops at the edge of the screen and a circle is placed on the last point. When I zoom out and get a dense pattern of 3DFACES I get a “bad argument type error”. So, I am wondering would it be possible to zoom or pan as the 3Dfaces are being selected so that the ssget points are always in the visible window?

[leeminardi]

@autoid374ceb4990 

I did a quick test of adding the following after the lines:

(setq pLL...

(setq pUR...

and before:

(setq ss...

(command "_zoom" "c" q 5)

and it seemed to help!  Try it out.

I am still working on improving the reliability of how to find the nect face.  A point size or small crss window does not consistently work. I plan to make a larger select cross window that will select several faces and then sort through them to determine which is adjacent to Q but is not the previous face.   I'm juggling a few other tasks right now.

[EDIT]  I just did a little more testing with the Zoom addition and it seems to solve the "ssget cross" problem!  Let me know if you get similar results.

[john.uhden]

@autoid374ceb4990 ,

Yes, I have had to include zooming as part of the code so that AutoCAD can "see" the target.

Of course I am using 2002.  More recent releases have some SELECTOFFSCREEN (or something like that) variable that is supposed to compensate for targets being offscreen.  But you are using R14, so you'll have to include the zooming.

[john.uhden]

@leeminardi ,

Try using a short fence selection across the previous edge.  With a 3Dface filter you should get only 2 members, one of which will be the previous 3Dface.  I'm not sure how you would handle searching from a corner.

[leeminardi]

@john.uhden  Zooming before ssget seems to have fixed the problem. 

Did my explanation of how to find the design slope direction answer your  question?

[autoid374ceb4990]

I added the ZOOM line and so far it seems to work well. I have tried it and created a 90 vertex POLYLINE all the way to the edge of my TIN with 6000 faces.

Now just add a VIEW save at the start and a VIEW restore at the end.

[john.uhden]

I'm sorry, @leeminardi .

Despite all your efforts I haven't studied your code at all.  Right now I'm trying to convince myself that it's no more than learning trigonometric identities, which reminds me of my first real job, with the Ocean County, New Jersey, Engineering Dept.  Dick Lane, the County Engineer, had no one on his staff that could solve trigonometry.  We were to calculate all the horizontal geometry for curves and right-of-way acquisitions for a 3 mile road widening project with jughandles galore.  When he found that I knew trigonometry, I inherited the assignment.  And all we had was a crank calculator (+ - * /) with no memory and books of sines, cosines, etc., so there were pages and pages of handwritten calculations with every letter in the Greek alphabet.  Of course it started as a summer job, so all they could pay me was $1.66 per hour as a laborer.  But it was fun and challenging and he taught me a lot about engineering.

[john.uhden]

@libarraSCXGZ , @leeminardi , @autoid374ceb4990 ,

I've been thinking (which can be dangerous) that they might be better off with a slope range map.  We used to be able to make them all the way back in DCA (before Softdesk and Land Desktop and C3D).

They would color code all the triangles that were in different slope ranges that you could specify, like 0%-2%, 2%-5%, 5%-10%, 10%-20%, 20%+, plus it would report the total area of each range.

That would give them a snapshot of where they might want to study first.

[autoid374ceb4990]

‎

@leeminardi

 In your post #66 you stated:

"ABC are the vertices, pAB the intersection of AB with the ground plane, pAC ...

Point D is the projection of A onto the line from pAB to pAC.

The slope = AE/DE "

How do you "project" A onto the line from pAB to pAC?

Report

[autoid374ceb4990]

OK, I have a general vector question concerning triangular 3DFACEs. I can find the face normal vector by calculating the cross product of two sides.  Is the slope (not slope vector) of the triangle with respect to a horizontal plane calculated by taking the acos of the dot product of the face normal and the normal to the horizontal plane (0,0,1)?

[leeminardi]

---

@autoid374ceb4990 wrote:

‎

@leeminardi

 In your post #66 you stated:

"ABC are the vertices, pAB the intersection of AB with the ground plane, pAC ...

Point D is the projection of A onto the line from pAB to pAC.

The slope = AE/DE "

How do you "project" A onto the line from pAB to pAC?

Report

---

@autoid374ceb4990 

To project point A onto a line defined by the two points pAB and pCA we first create a unit vector in the direction from pAB to pCA.   Let's call it uAB_CA.

To find D, the projection of A onto the line, we can take the dot product of a vector from anywhere on the line to A with the unit vector   and add it to the beginning of the line pAB.

D = pAB + (dot   uAB_CA  (A - pAB))

In lisp it woud look like:

(setq d (distance pAB pCA))
(setq uAB_CA (mapcar '/
		     (mapcar '- pCA pAB)
		     (list d d d)
	     )
)
(setq D	(mapcar	'+
		pAB
		(dot (mapcar '- A pAB) uAB_CA)
	)
)

[leeminardi]

---

@autoid374ceb4990 wrote:

OK, I have a general vector question concerning triangular 3DFACEs. I can find the face normal vector by calculating the cross product of two sides.  Is the slope (not slope vector) of the triangle with respect to a horizontal plane calculated by taking the acos of the dot product of the face normal and the normal to the horizontal plane (0,0,1)?

---

The slope can be expressd as a ratio of rise/run or by an angle.  If we want the former then we need to divide the vertical component of the slope vector by the length of the vector in the horizontal (XY) plane. Line #209 of my post #63 does that via:

  (setq	slopeFace
	 (/
	   (caddr slopeVector)
	   (distance '(0 0 0)
		     (list (car slopeVector) (cadr slopeVector) 0.0)
	   )
	 )
  )

If you want the angle then use the atan function with the z component divided by the length of the vector in the horizontal (XY) plane (similar to the above code).

Of course, since the normal vector is perpendicular to the slope you could find the slope of the normal and take its recipricol for the slope of the face.

Note, at the end of my code in post #63 are functions for atan, dot product, cross product and unit vector.

BTW, I've found some bugs in my latest post of code (#63) and am in the process of debugging and cleaning up the code and adding a feature to allow left or right bias for finding the path at the specified slope.