Applied Mathematics concerning CAD Faces

This comprehensive page is a practical collection about CAD faces with layman explanations of mathematics and additional code fragments. It deals with vectors, angles, normals, 2D projections, rotations and more. The best news is that even BricsCAD Classic supports both Lisp and 3D commands like RuleSurf. Even with limited resources you can build great solutions.

I made FlatMesh in an attempt to create blanks from ruled surfaces. While the basics are relative easy – translating 3D coordinates to 2D – it turned out that the exceptions are the real math challenge. I couldn’t find a this page like this, in order to clarify and understand the subject. So this page was born. I hope it is of help. If you think something is wrong, missing or can be done better, just drop a comment.

Definitions

  • Mesh: A mesh is a collection of faces. A mesh forms a surface object. Often, representation in CAD is like a collection of rectangles. The use of rectangles has the advantage of clarity. On the other hand, in 3D, the four corners of each rectangle are almost never in one plane. To fix that, each rectangle consists of two adjacent triangles called faces. The common side of these two faces forms an invisible diagonal of a visible rectangle – the remaining sides of the two triangles.
  • Face: See Mesh. A face is by definition a flat polygon and almost always a triangle in a CAD environment. A face has three edges and vertices. Complex surfaces can be simulated by using triangular faces. Faces with four vertices are not suitable for complex surfaces.
  • Vector: A vector is a direction with a magnitude.
    • A vector V is represented as \vec V
    • with a direction x,y,z
    • and a length or magnitude \|V\|
  • Dot product: An algebraic operation on two vectors \vec a and \vec b with angle \theta in between, resulting in a scalar (number, not a vector) of the area of the parallelogram of \vec a and \vec b, written as
    a \cdot b = \|a\| \|b\| \sin (\theta). More…
  • Cross product: An algebraic operation on two vectors \vec a and \vec b with angle \theta in between, resulting in a vector perpendicular to both \vec a and \vec b, written as
    a \times b = \|a\| \|b\| \sin (\theta) n. The magnitude of the resulting vector equals the area of the parallelogram of \vec a and \vec b. Note that the area – part of both dot and cross product – is valuable for solving equations. n is positive when using the right hand rule. More…
  • Normal: A line, ray or vector perpendicular to something. In this case we mean a vector perpendicular to a plane. The side of the plane on which the normal starts is determined by the right hand rule. Therefore the order of a and b is important, a \times b \ne b \times a and a \times b = - ( b \times a ) . More…

Math solutions

Calculate the magnitude of a vector

For vector \vec{A} with value x,y,z, the magnitude or length is \sqrt{x^2+y^2+z^2} according to the Pythagorean theorem.

For vector \vec{AB} with A=Ax,Ay,Az and B=Bx,By,Bz, the length is

\sqrt{(Bx-Ax)^2+(By-Ay)^2+(Bz-Az)^2}.

Calculate the unit vector of a vector

A unit vector of a vector has the same direction but has a length or magnitude of 1. It can be retrieved by dividing each x, y and z of a vector by the magnitude of that vector, resulting in a new vector.

Calculate the angle between two vectors

In 2D, it is easy to comprehend. 3D is harder. For example grep a pencil in each hand and position and point them randomly before you. What on earth is the angle? Now bring the origins together. At that point, you can imagine a plane through the three points and it becomes clear there is an angle. Then the question becomes: What is the angle? And: Do you want the angle from pencil A to pencil B or the other way around? Finally, you can look at the plane from to sides, which changes the answers too. So there is enough complexity.

Using arccosine

The angle \alpha between \vec {AB} and \vec {AC}

\alpha = \arccos \Big( \frac {\vec {AB} \cdot \vec {AC}} {\|AB\| \|AC\|} \Big)

Example

AB = 1,2,3
AC = 4,5,6
|AB| = SQRT(1^2+2^2+3^2) = 3.74
|AC| = SQRT(4^2+5^2+6^2) = 8.77
AB.AC = 1*4+2*5+3*6 = 32
acos((AB.AC)/(|AB|*|AC|)) = acos(32/(3.74*8.77)) = 0.226 rad = 12.9 degrees

Here, the angle range is 0 to \pi. Although technically correct, that may not always be what you want.

Using arctangent

The same angle again…

\alpha = atan2 \Big( \big( (\vec {AB} \times \vec {AC}) \cdot \vec n ) ,  (\vec {AB} \cdot \vec {AC} \big) \Big)

Is a vector angle obtuse, perpendicular or acute?

Continuing with the previous… Instead of, or in addition to, calculating an angle, the dot product is enough to say whether the angle is obtuse, perpendicular or acute. It is depending on the sign of m = \vec {AB} \cdot \vec {AC}. More specific:

  • If m is positive, the angle is acute.
  • If m is zero, the angle is right, perpendicular.
  • If m is negative, the angle is obtuse.

A rectangle contains two adjacent faces. We can say something about the angle between the faces by first calculating the normal vectors of each face.

Calculate the normal of a face

The normal vector is the cross product of vector AB times vector AC. According to the right hand rule, it is written as \vec {AB} \times \vec {AC}

Example

AB = (1,2,3)
AC = (4,5,6)
AB🞨AC = (2*6-3*5,1*6-3*4,1*5-2*4) = (-3,6,-3)

The opposite normal vector is AC🞨AB, being AB🞨AC with opposite signs, i.e. (3,-6,3).

The magnitude of normal AB🞨AC is area \|AB\| \|AC\| \sin (\angle CAB). So you can tell something about angle CAB if you have the magnitudes of all three vectors.

Calculate 2D projection of a face – Pythagorean

With three coordinates A, B and C known of a 3D face, we can calculate a projection in 2D. Alternative, law of cosines can be used here too.

Calculating C based on A is (0,0) and B is (p,0)

Summary

Known:
Distances p, q, r
Constraints:
Point A = (0,0)
Point B = (p,0)
Point C is above axis AB
First Cy is positive
z coordinates always 0
Requested:
Coordinates of C.
Answer:
C = ((p^2+r^2-q^2)/(2*p),sqrt(r^2-t^2))
Subsequent Cy values can be negative too, i.e. -sqrt(r^2-t^2)

Rationale

A: (x,y,z)
C: (x',y',z')

Variant 1 \vee Variant 2 (\vee is OR)
\angle CAB <= pi()/2 \vee  \angle CAB > pi()/2 =>
x <= x' \vee x > x'

r^2 = t^2+s^2 =>
s^2 = r^2-t^2

q^2 = (p-t)^2+s^2 \vee (p+t)^2+s^2
= (p-t)^2+r^2-t^2 \vee (p+t)^2+r^2-t^2
= p^2-2pt+t^2+r^2-t^2 \vee p^2+2pt+t^2+r^2-t^2
= p^2-2pt+r^2 \vee p^2+2pt+r^2 =>
2pt = p^2+r^2-q^2 \vee -p^2-r^2+q^2 =>
t = (p^2+r^2-q^2)/(2*p) \vee (p^2+r^2-q^2)/(-2*p)
s^2=r^2-t^2
s=sqrt(r^2-t^2) \vee s=-sqrt(r^2-t^2)

With two answers for t it is tempting to just code both values. Annoying and itching is the fact I cannot say anything beforehand about CAB.

Hmm, is that so? Climbing from the bottom – the resulting formula’s – up, there is fragment p^2+r^2-q^2. Wait, Pythagoras! That should be 0 when \angle CAB is 90 ^{\circ}!

More specific:

  • p^2+r^2-q^2 < 0 \Rightarrow \angle CAB > 90 ^{\circ}
  • p^2+r^2-q^2 = 0 \Rightarrow \angle CAB = 90 ^{\circ}
  • p^2+r^2-q^2 > 0 \Rightarrow \angle CAB < 90 ^{\circ}

One problem solved, one problem created! Let’s see, If p^2+r^2-q^2 < 0then t = (p^2+r^2-q^2)/(-2*p) and else t = (p^2+r^2-q^2)/(2*p). So t is always a positive value and whether it is on the left or right side of x=0 (point A) depends on the value of p^2+r^2-q^2.

More specific, as an x-coordinate, t = (p^2+r^2-q^2)/(2*p). And that is always true! Great, and even my old brains survived this.

A final note about s. It is a square root so the argument can be both positive and negative. It won’t bother us for the assumed first C, but it will bother us for subsequent C’s and all D’s.

Calculate 2D projection of an adjacent face

Getting the coordinate D is similar to finding C.

One segment of the RuleSurf result is build out of two triangles or faces. Using FlatMesh means in most cases: Start with A and B and calculate C and D, then treat those resulting C and D as a new set of A and B and start allover, until the complete mesh is processed. In most cases, the first picture is in play. However, sometimes the second picture can be in play. In addition, even C can be under axis AB but focus is on D for proper explanation.

A,B,C and D are 3D points. There is a plane ABC. D is somewhere relative to this plane. We need to know where D lies. More specific, in a planar view to ABC, is D under or above axis BC? To get a better impression, look at it from a different angle:

This question seems not easy to answer, but, as long as we talk about angles, we can say that:

  • If the angle between the face normals is acute, D lies above BC.
  • If the angle between the face normals is obtuse, D lies under BC.

After calculating the normals as explained before, the dot product is all we need to answer the question. Just remember to use the right hand rule.

That was actually the hard part. D is more or less calculated the same way as C.

However, we need to use the axis CB (not AB) as the base, as the mirror line for when D gets under CB. Finally, C and B were rotated around A, D should also be rotated properly. That is where \alpha starts to play a role. Rotating…

Calculate rotation of a coordinate in 2D

Seriously written for mathematicians, so don’t read https://en.wikipedia.org/wiki/Rotation_(mathematics). Pff… Okay, it boils down to this:

Given is coordinate pair x,y and angle α. Get x',y' by using:

x'=x*cos(α)-y*sin(α)
y'=y*cos(α)+x*sin(α)

In order not to get unexpected results, angle α is positive when counter-clockwise – in this example it is a positive value.

Translations of coordinates

This one is easy, translations like command Move is just adding static values to x, y and z values.

Programming solutions

The purpose of this code is educational. It works and aims to show the structure. Unless you experience slow processing, it should be sufficient.

Lisp lacking Arcsine and Arccosine in AutoCAD

This code can be used for (acos ...):

1
(defun acos (r / ) (atan (sqrt (- 1 (expt r 2))) r))

This code can be used for (asin ...):

1
(defun asin (r / ) (atan r (sqrt (- 1 (expt r 2)))))

BricsCAD uses the native Lisp function. Tip: Ditch AutoCAD for Lisp programming, BricsCAD CAD-Lisp is way faster and is constantly under development with for example not a two decades old IDE, search “Blade BricsCAD”.

Distance between point a and b

In Lisp there is a function (distance a b) but it is tempting to use your own function because (distance a b) treats all points 2D when it encounters one 2D point. An alternative function (dist3d a b) works always 3D and calculates “the square root of the sum of the squares of the delta x, y and z values”. It could be as simple as this one-liner :

1
2
3
4
5
6
7
8
9
(defun dist3d (a b / )
   (sqrt
      (+
         (expt (- (car b) (car a)) 2)
         (expt (- (cadr b) (cadr a)) 2)
         (expt (- (caddr b) (caddr a)) 2)
      )
   )
)

With the distance between a and b, a is x,y,z and b is x',y',z', the formula in spreadsheet notation is:
dist3d=sqrt(((x'-x)^2)+((y'-y)^2)+((z'-z)^2))

Vector between point a and b

For a vector coordinate set from a to b, subtract coordinates of a from b. In Lisp:

1
2
3
4
5
6
7
(defun vector (a b / )
   (list
      (- (car b) (car a))
      (- (cadr b) (cadr a))
      (- (caddr b) (caddr a))
   )
)

Cross product of vector a and b

Take care of the order of a and b, it affects the sign of the coordinates. In Lisp:

1
2
3
4
5
6
7
(defun cross-product (a b / )
   (list
      (- (* (cadr a) (caddr b)) (* (cadr b) (caddr a)))
      (- (* (car a) (caddr b)) (* (car b) (caddr a)))
      (- (* (car a) (cadr b)) (* (car b) (cadr a)))
   )
)

Dot product of vector a and b

In Lisp:

1
2
3
4
5
6
7
(defun dot-product (a b / )
   (+
      (* (car a) (car b))
      (* (cadr a) (cadr b))
      (* (caddr a) (caddr b))
   )
)

Magnitude of vector a

In Lisp:

1
2
3
4
5
6
7
8
9
(defun magnitude (a / )
   (sqrt
      (+
         (expt (car a))
         (expt (cadr a))
         (expt (caddr a))
      )
   )
)

Unit vector of vector a

In Lisp:

1
2
3
4
5
6
7
8
(defun unit-vector (a / l)
   (setq l (magnitude a)
   (list
      (/ (car a) l)
      (/ (cadr a) l)
      (/ (caddr a) l)
   )
)

Return value: A list containing the x, y and z values of the unit vector.

Acute angle between vector a and b

In Lisp, is an angle acute or obtuse? This code treats right angles as acute which makes the function very fast but not always suitable, check first.

1
2
3
(defun acute (a b / )
   (if (&gt;= (dot-product a b) 0.0) T)
)

Return values: If the angle is right or acute, T (True) is returned. the dot product is zero or positive, T (True) is returned. If the angle is obtuse, nil is returned.

Vector angle between a and b (acos)

In Lisp, see angle between vectors based on acos.

Return value: Angle in radians.

1
2
3
(defun vector-ang (a b / )
   (acos (/ (dot-product a b) (* (magnitude a) (magnitude b))))
)

Rotation of point a over angle g

In Lisp. Angle g in radians, sign is positive when counter clock wise. Point a: only x and y are evaluated.

Return value: A list containing the calculated x and y values after rotation.

1
2
3
4
5
6
7
(defun rot-pt (a g / sg cg x y)
   (setq sg (sin g) cg (cos g) x (car a) y (cadr a))
   (list
      (- (* x cg) (* y sg))
      (+ (* y cg) (* x sg))
   )
)
Posted in CAD

Leave a Reply

Your email address will not be published. Required fields are marked *