A non-iterative method for robustly computing the intersections between a line and a curve or surface

TL;DR

This paper presents a robust, non-iterative linear algebra-based method for computing intersections between lines and high-order curves or surfaces, improving reliability over traditional iterative approaches.

Contribution

The paper introduces a novel, non-iterative intersection computation method using SVD and eigenvalue problems, simplifying implementation and increasing robustness for high-order finite element applications.

Findings

The method determines all intersection points in a single step.

It is more robust than Newton-Raphson iteration.

The approach is applicable to high-order finite elements.

Abstract

The need to compute the intersections between a line and a high-order curve or surface arises in a large number of finite element applications. Such intersection problems are easy to formulate but hard to solve robustly. We introduce a non-iterative method for computing intersections by solving a matrix singular value decomposition (SVD) and an eigenvalue problem. That is, all intersection points and their parametric coordinates are determined in one-shot using only standard linear algebra techniques available in most software libraries. As a result, the introduced technique is far more robust than the widely used Newton-Raphson iteration or its variants. The maximum size of the considered matrices depends on the polynomial degree $q$ of the shape functions and is $2q \times 3q$ for curves and $6 q^2 \times 8 q^2$ for surfaces. The method has its origin in algebraic geometry and has here been considerably simplified with a view to widely used high-order finite elements. In addition, the method is derived from a purely linear algebra perspective without resorting to algebraic geometry terminology. A complete implementation is available from http://bitbucket.org/nitro-project/.

Figures

Equations

$$\mathbf{x}(\theta )=\sum _{i=1}^{q_x+1}{L}_i(\theta ){\mathbf{x}}_i=\sum _{j=1}^{q_x+1}{P}_j(\theta ){\mathbf{\alpha }}_j,$$

$$P_j(\theta )=\sum _{i=1}^{q_x+1}{L}_i(\theta )P_j({\theta }_i),$$

$$l(\theta ,\mathbf{x})=\left(\begin{array}{c}\mathbf{x}\\ 1\end{array}\right)\cdot \mathbf{g}(\theta )=x^1{g}^1(\theta )+x^2{g}^2(\theta )+g^3(\theta )=0,$$

$$\mathbf{g}(\theta )=\sum _{l=1}^{q_g+1}{\tilde{P}}_l(\theta ){\mathbf{g}}_l.$$

$$l(\theta ,\mathbf{x}(\theta ))=\left(\begin{array}{c}\mathbf{x}(\theta )\\ 1\end{array}\right)\cdot \mathbf{g}(\theta )=\left(\begin{array}{c}\mathbf{x}(\theta )\\ 1\end{array}\right)\cdot \left(\sum _{l=1}^{q_g+1}{\tilde{P}}_l(\theta ){\mathbf{g}}_l\right)=0.$$

$$\sum _{l=1}^{q_g+1}\left(\sum _{j=1}^{q_x+1}\left(\begin{array}{c}{\mathbf{\alpha }}_j\\ 1\end{array}\right)P_j(\theta ){\tilde{P}}_l(\theta )\right){\mathbf{g}}_l=0.$$

$$\sum _{k=1}^{q_x+q_q+1}{\hat{P}}_k(\theta )C_{kl}=\sum _{j=1}^{q_x+1}\left(\begin{array}{c}{\mathbf{\alpha }}_j\\ 1\end{array}\right)P_j(\theta ){\tilde{P}}_l(\theta ),$$

$$\sum _{l=1}^{3(q_g+1)}\left(\sum _{k=1}^{q_x+q_g+1}{\hat{P}}_k(\theta )C_{kl}\right)h_l=0,$$

$$\mathbf{h}=\text{setcounter}MaxMatrixCols12{\left(\begin{array}{cccccccccccc}{g}_1^1& g_2^1& \dots & g_{q_g+1}^1& g_1^2& g_2^2& \dots & g_{q_g+1}^2& g_1^3& g_2^3& \dots & g_{q_g+1}^3\end{array}\right)}^{\mathsf{T}}.$$

$$2q_g-q_x+2\ge q_x\Rightarrow q_g\ge q_x-1.$$

$$l^{(i)}(\theta ,\mathbf{x})=\left(\begin{array}{c}\mathbf{x}\\ 1\end{array}\right)\cdot \left(\sum _{l=1}^{q_g+1}{\tilde{P}}_l(\theta ){\mathbf{g}}_l^{(i)}\right)=0$$

$$\mathbf{r}(\xi )=\left(\begin{array}{c}{r}^1(\xi )\\ r^2(\xi )\end{array}\right)=\left(\begin{array}{c}{c}_1^1\\ c_1^2\end{array}\right)\xi +\left(\begin{array}{c}{c}_0^1\\ c_0^2\end{array}\right)$$

$$l^{(i)}(\theta ,\mathbf{r}(\xi ))=\left(\begin{array}{c}\mathbf{r}(\xi )\\ 1\end{array}\right)\cdot \left(\sum _{l=1}^{q_g+1}{\tilde{P}}_l(\theta ){\mathbf{g}}_l^{(i)}\right)=0.$$

$$\sum _{l=1}^{q_g+1}{\tilde{P}}_l(\theta )\left[\left(\begin{array}{c}\mathbf{r}(\xi )\\ 1\end{array}\right)\cdot \left({\mathbf{g}}_l^{(i)}\right)\right]=0.$$

$$\sum _{l=1}^{q_g+1}{\tilde{P}}_l(\theta )\left[A_{li}-\xi B_{li}\right]=0.$$

$$\mathbf{\varphi }\left(\mathbf{A}-\xi \mathbf{B}\right)=0.$$

$${\theta }^{(j)}=\frac{{\tilde{P}}_{i+1}}{{\tilde{P}}_i}=\frac{{\varphi }_{i+1}^{(j)}}{{\varphi }_i^{(j)}}.$$

$$({\mathbf{\Delta }}_{i+1}-\theta {\mathbf{\Delta }}_i)\mathbf{\psi }=0.$$

$$\mathbf{x}(\theta )=\sum _{j=1}^4{P}_j(\theta ){\mathbf{\alpha }}_j=\left(\begin{array}{cccc}0& 4& -4.5& 4.5\\ 0& 11.25& -31.5& 20.25\end{array}\right)\left(\begin{array}{c}{P}_1(\theta )\\ P_2(\theta )\\ P_3(\theta )\\ P_4(\theta )\end{array}\right).$$

$$\mathbf{r}(\xi )=\left(\begin{array}{c}0\\ 1\end{array}\right)+\xi \left(\begin{array}{c}4\\ -2\end{array}\right)=\left(\begin{array}{c}4\xi \\ 1-2\xi \end{array}\right)$$

$$\mathbf{g}(\theta )=\sum _{l=1}^4{\tilde{P}}_l(\theta ){\mathbf{g}}_l.$$

$$\mathbf{C}=\text{setcounter}MaxMatrixCols12\left(\begin{array}{cccccccccccc}[r]0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 4& 0& 0& 0& 11.25& 0& 0& 0& 0& 1& 0& 0\\ -4.5& 4& 0& 0& -31.5& 11.25& 0& 0& 0& 0& 1& 0\\ 4.5& -4.5& 4& 0& 20.25& -31.5& 11.25& 0& 0& 0& 0& 1\\ 0& 4.5& -4.5& 4& 0& 20.25& -31.5& 11.25& 0& 0& 0& 0\\ 0& 0& 4.5& -4.5& 0& 0& 20.25& -31.5& 0& 0& 0& 0\\ 0& 0& 0& 4.5& 0& 0& 0& 20.25& 0& 0& 0& 0\end{array}\right).$$

$$\mathbf{A}=\left(\begin{array}{ccccc}[r]0.08710& 0.08740& 0.04298& 0.05767& 0.01708\\ -0.2222& -0.02814& 0.8999& -0.01814& 0.03355\\ 0.04655& -0.07068& -0.03786& 0.9989& 0.01874\\ -0.03786& 0.1289& -0.01752& -0.001102& 1.010\end{array}\right)$$

$$\mathbf{B}=\left(\begin{array}{ccccc}[r]1.046& 1.101& 1.499& 0.7322& 0.2201\\ -3.351& -1.347& -0.1816& 0.04373& 0.1881\\ 2.390& -2.402& -0.1053& 0.2525& 0.3777\\ -0.1302& 2.987& -0.3424& 0.1006& 0.3987\end{array}\right).$$

$${\mathbf{\varphi }}^{\square }\left({\mathbf{A}}^{\square }-\xi {\mathbf{B}}^{\square }\right)=0.$$

$${\xi }^{(1)}=0.8112,{\xi }^{(2)}=0.3594,{\xi }^{(3)}=0.08875,{\xi }^{(4)}=0.05326,$$

$${\xi }^{(1)}=28.05,{\xi }^{(2)}=0.8112,{\xi }^{(3)}=0.3594,{\xi }^{(4)}=0.08875.$$

$${\mathbf{x}}^{(1)}=(3.245,-0.6225{)}^{\mathsf{T}},{\mathbf{x}}^{(2)}=(1.438,0.2813{)}^{\mathsf{T}},{\mathbf{x}}^{(3)}=(0.3550,0.8225{)}^{\mathsf{T}}.$$

$${\mathbf{x}}^{(1)}$$

$${\mathbf{x}}^{(4)}$$

$${\theta }^{(1)}=0.9014,{\theta }^{(2)}=0.5,{\theta }^{(3)}=0.09861,{\theta }^{(4)}=0.03932.$$

Full text

1 Introduction

There has recently been an increased academic and industrial interest in high-order finite elements due to their efficiency advantages over classical low-order elements, see e.g.[1, 2]. To achieve their full potential, high-order methods require the curved domain boundaries to be approximated with non-planar elements. The intersection between the curved elements and a given line is required in a wide range of applications, including contact [3], mesh generation [4, 5] and immersed finite elements [6]. As known, intersection computations lead to an easy to formulate, but hard to solve nonlinear root-finding problem. The prevalent technique in computational mechanics for solving such problems is Newton-Raphson iteration, which is in general not very robust, especially when no good starting points are available or multiple intersections are present.

In computer-aided geometric design (CAD) and manufacturing (CAM) intersection computation is a recurring task and, to this end, a number of ingenious methods have been developed [7]. Especially promising are the non-iterative methods with an origin in algebraic geometry, which are for the most part unknown in computational mechanics. Algebraic geometry deals with systems of polynomial equations and geometric objects defined by them and provides the most rigorous framework for intersection computations [8]. The specific technique considered in this paper is the implicitisation technique proposed by Busé [9], which shares some commonalities with the method of moving lines/planes introduced by Sederberg and Chen [10]. Different from the original work, in this paper we derive the method from a purely linear algebra viewpoint. To follow the presented derivations, it is sufficient to only know the notion of the algebraic degree of a curve or surface. The algebraic degree of a curve or surface is defined as its number of intersections with a line. Counting all intersections (real, complex, multiple ones and ones at infinity), the algebraic degree of a polynomial curve of degree $q$ and the corresponding tensor-product surface are $q$ and $2q^{2}$, respectively. As will become clear, adopting a linear algebra viewpoint has the added benefit that many linear algebra techniques, like preconditioning and various matrix decompositions, become readily available for intersection computations.

In the following, we first discuss the intersection between a planar Lagrange curve and a line and provide an easy to follow illustrative example. This simple case is sufficient to introduce and discuss the key aspects of the proposed non-iterative technique. Its extension to the surface case is straightforward and is discussed last.

2 Intersection of lines with curves

2.1 Moving lines and intersections

Let $\bm{x}(\theta)=\left(x^{1}(\theta),\,x^{2}(\theta)\right)^{\mathsf{T}}$ be a planar parametric curve, with $\bm{x}(\theta)\in\mathbb{R}^{2}$, of degree $q_{x}$ given either in Lagrange basis $L_{i}(\theta)$ or power (monomial) basis $P_{j}(\theta)$ with

[TABLE]

where $\bm{x}_{i}\in\mathbb{R}^{2}$ are the nodal coordinates and $\bm{\alpha}_{j}\in\mathbb{R}^{2}$ are the coefficients in the power basis. As usual, the power basis $P_{j}(\theta)$ contains the consecutive powers of $\theta$ from [math] up to $q_{x}$. The two basis are related by

[TABLE]

where $P_{j}(\theta_{i})$ is the Vandermonde matrix and $\theta_{i}$ is the parametric coordinate of the $i$-th Lagrange node. Following a similar approach a curve given in any other polynomial basis can be re-expressed in the power basis.

To define a point $\bm{x}(\theta)$ on the curve as the intersection of several moving lines, or pencils of lines, consider the line

[TABLE]

where $\bm{x}\in\mathbb{R}^{2}$ and  $\bm{g}(\theta)=\left(g^{1}(\theta),\,g^{2}(\theta),\,g^{3}(\theta)\right)^{\mathsf{T}}$ is an auxiliary vector collecting the parameters of the line. For a fixed $\theta$ equation (3) describes a line and the line moves with the parameter $\theta$, see Figure 1a . The three parameters of the line are assumed to be polynomial functions given by

[TABLE]

The degree $q_{g}$ of the power basis $\widetilde{P}_{l}(\theta)$ has to be chosen sufficiently high in order to be able to compute all the intersection points (real, complex, multiple ones and ones at infinity). The number of intersection points is equivalent to the algebraic degree of the curve $\bm{x}(\theta)$. A curve $\bm{x}(\theta)$ of degree $q_{x}$ has $q_{x}$ intersection points with a line. As will become clear, the number of intersections implies a constraint on the minimum possible degree $q_{g}$ for $\widetilde{P}_{l}(\theta)$.

Next, the aim is to find several moving lines of the form (3), or more specifically their coefficients $\bm{g}_{l}$ in (4), with a common intersection point which is a point on the curve $\bm{x}(\theta)$. It is required that each line satisfies at the common intersection point

[TABLE]

After introducing the definition of the curve (1) this yields

[TABLE]

The bracketed term can be expressed in a new power basis $\widehat{P}_{k}$ of dimension $q_{x}+q_{g}+1$ with

[TABLE]

where the matrix components $C_{kl}$ contain the known coefficients $\bm{\alpha}_{j}$. Equation (6) can now be rewritten as

[TABLE]

where the array $\bm{h}$ contains the components of the yet unknown vectors $\bm{g}_{l}$ sorted (by choice) in the following way

[TABLE]

It is required that (8) is always satisfied irrespective of $\theta$, which is the case for the right null vectors of the matrix $C_{kl}$. The right null vectors are determined with a SVD, see e.g. [11], yielding the set of null vectors $\bm{g}_{l}^{(i)}$, where the index $(i)$ denotes the number of the null vector.

The number of null vectors of $C_{kl}$ depends on the degrees $q_{x}$ and $q_{g}$ of the basis $P_{j}(\theta)$ and $\widetilde{P}_{l}(\theta)$, and the coefficients $\bm{\alpha}_{j}$ of the specific curve considered. For subsequent computations the number of null vectors must be more than the number of intersections of the curve with a line (or its algebraic degree). The non-square matrix $C_{kl}$ has $q_{x}+q_{g}+1$ rows and $3(q_{g}+1)$ columns. Hence, its number of right null vectors must be equal or greater than $3(q_{g}+1)-(q_{x}+q_{g}+1)=2q_{g}-q_{x}+2$.111 The number of right null vectors is larger when, for instance, a quadratic curve is described with a cubic polynomial (i.e. $q_{x}=3$ but $\alpha_{4}=0$). More precisely, if $n$ is the largest integer such that $\bm{\alpha}_{n}\neq 0$ then the number of right null vectors of $C_{kl}$ is exactly $2q_{g}-n+3$.

In order to obtain the $q_{x}$ intersections it is necessary to have

[TABLE]

The set of null vectors denoted with $\bm{g}_{l}^{(i)}$ introduced in (3) yields a set of moving lines

[TABLE]

with a common intersection point on the curve $\bm{x}(\theta)$. As an example, in Figure 1a the description of a quadratic curve with $q_{x}=2$ by two moving lines $l^{(1)}(\theta,\,\bm{x})$ and $l^{(2)}(\theta,\,\bm{x})$ with $q_{g}=1$ is shown.

Next, the intersection of a given parametric line

[TABLE]

with the curve $\bm{x}(\theta)$ is considered, where $\bm{c}_{1}=(c_{1}^{1},\,c_{1}^{2})^{\mathsf{T}}$ and $\bm{c}_{0}=(c_{0}^{1},\,c_{0}^{2})^{\mathsf{T}}$ are two given vectors. We require that the determined set of moving lines $\{l^{(i)}(\theta,\,\bm{x})=0\}$ and $\bm{r}(\xi)$ intersect at the same point to be yet determined. According to (5), at the common intersection point it is required that

[TABLE]

These equations describe the intersection of each moving line $l^{(i)}(\theta,\,\bm{x})=0$ with the given line $\bm{r}(\xi)$ and can be rewritten as

[TABLE]

As discussed the number of moving lines satisfies $\max(i)\geq q_{x}$ and all of them can be combined in one homogenous equation system

[TABLE]

with $A_{li}$ and $B_{li}$ representing the components of two matrices $\bm{A}$ and $\bm{B}$. Both matrices have $q_{g}+1$ rows and $2q_{g}-q_{x}+2$ or more columns. To obtain $\widetilde{P}_{l}(\theta)$ and $\xi$ that satisfy (15) the following generalised eigenvalue problem is considered

[TABLE]

That is, the eigenvalues $\xi^{(j)}$ are the parametric coordinates of the intersection points on the line $\bm{r}(\xi)$ and the eigenvectors are (up to a multiplicative constant) the basis functions $\widetilde{P}_{l}(\theta^{(j)})$ evaluated at the intersection points $\theta^{(j)}$. Unfortunately, the matrices $\bm{A}$ and $\bm{B}$ are not always square and computing the values $\xi^{(j)}$ that satisfy this equation requires non-standard linear algebra techniques. However, as will be discussed in Section 2.2, for the purposes of intersection computation it is sufficient to consider a square eigenvalue problem obtained from (16) by taking only some of its columns. The non-complex eigenvalues of this square eigenvalue problem contain all the intersection points between the given line and the curve. Some of these non-complex eigenvalues may not be actual intersection points, but they can easily be identified.

A non-complex eigenvalue $\xi^{(j)}$ of (16), or its corresponding square eigenvalue problem, gives the potential intersection point $\bm{r}(\xi^{(j)})$. The respective unknown parameter value $\theta^{(j)}$ on the curve satisfies the equation $\bm{x}(\theta^{(j)})=\bm{r}(\xi^{(j)})$, which is a hard to solve nonlinear problem. According to (15) and (16), however, the left null vector $\bm{\phi}^{(j)}$ corresponding $\xi^{(j)}$ is proportional to the vector $\widetilde{P}_{l}(\theta^{(j)})$, which is exploited to determine $\theta^{(j)}$. More precisely, if there is a single parameter value $\theta^{(j)}$, the ratio of any two consecutive components yields

[TABLE]

It is assumed here that the monomials in $\widetilde{P}_{i}$ are sorted in increasing order so that the ratio of two consecutive entries is simply $\theta^{(j)}$. If there are $p$ parameter values $(\theta^{(j,1)},\,\dotsc,\theta^{(j,p)})$ corresponding to the eigenvalue $\xi^{(j)}$, then the corresponding left null vector space is given by a matrix $\bm{K}$ with $p$ rows. Keeping the previous assumption on the ordering of the monomials in $\widetilde{P}_{i}$, we define the $p\times p$ matrices $\bm{\Delta}_{i}$ by taking the columns $i$ to $p+i$ of $\bm{K}$. Then, the parameter values $(\theta^{(j,1)},\,\dotsc,\theta^{(j,p)})$ are obtained by solving a generalised eigenvalue problem

[TABLE]

2.2 Illustrative example

The intersection of a cubic Lagrange curve with a line depicted in Figure 2 is considered next. The cubic curve $\bm{x}(\theta)$ interpolating the points $\bm{x}_{1}=(0,\,0)^{\mathsf{T}}$, $\bm{x}_{2}=(1,\,1)^{\mathsf{T}}$, $\bm{x}_{3}=(2,\,-0.5)^{\mathsf{T}}$ and $\bm{x}_{4}=(4,\,0)^{\mathsf{T}}$ is expressed in power basis with

[TABLE]

The intersections between the given line

[TABLE]

and the curve $\bm{x}(\theta)$ are sought.

As discussed, the degree $q_{g}$ of the auxiliary polynomial needs to satisfy $q_{g}\geq q_{x}-1=2$. In this example choosing $q_{g}=2$ yields matrices $\bm{A}$ and $\bm{B}$ that have the dimensions $3\times 3$ and it is straightforward to compute the eigenvalues $\xi$ of the generalised eigenvalue problem (16). To illustrate the more challenging case with non-square matrices $\bm{A}$ and $\bm{B}$, which turns out to be inevitable in the case of surfaces, we choose here $q_{g}=3$ such that

[TABLE]

The right null vectors $\bm{g}_{l}^{(i)}$ that define the moving lines are obtained from (8) with $\bm{C}$ and the components

[TABLE]

This matrix has the dimensions $7\times 12$, its rank is seven and the dimension of its right null space is five. Each of the corresponding five independent null vectors $\bm{g}_{l}^{(i)}$ defines one moving line. Figure 1b depicts the first four moving lines $l^{(1)}(\theta,\,\bm{x})=0$, $l^{(2)}(\theta,\,\bm{x})=0$, $l^{(3)}(\theta,\,\bm{x})=0$ and $l^{(4)}(\theta,\,\bm{x})=0$ at the parameter values $\theta=0.4$ and $\theta=0.8$.

Substituting the line (20) as in (15) gives the two non-square matrices (with four significant digits)

[TABLE]

and

[TABLE]

One approach to obtaining the generalised eigenvalues of matrices $\bm{A}$ and $\bm{B}$ is to use pencil reduction, see [12, 13, 14], which is not a widely used linear algebra operation and may introduce additional numerical issues because of several numerical rank estimations. Alternatively, the eigenvalue problems defined by square submatrices $\bm{A}^{\Box}$ and $\bm{B}^{\Box}$ of the largest size, e.g. the first four columns of $\bm{A}$ and $\bm{B}$, can be considered

[TABLE]

Although there can only be three intersection points for a cubic curve, this problem has four eigenvalues and eigenvectors. Notice that each column $i$ of $\bm{A}-\xi\bm{B}$, i.e. $\sum_{l}\phi_{l}(A_{li}-\xi B_{li})=0$, represents the intersection between the line $\bm{r}(\xi)$ and the moving line $l^{(i)}(\theta,\,\bm{x})=0$. Taking a different couple of submatrices $\bm{A}^{\Box}$ and $\bm{B}^{\Box}$ of largest size, like the last four columns of $\bm{A}$ and $\bm{B}$, yields a different set of intersection points. Both sets of intersection points contain the three true intersection points in addition to one fictitious intersection point.

To demonstrate this observation, the intersection points with two different pairs of largest square submatrices $\bm{A}^{\Box}$ and $\bm{B}^{\Box}$ are computed. The eigenvalues for the problem defined by the first four columns of $\bm{A}$ and $\bm{B}$ are

[TABLE]

and the eigenvalues for the problem defined by the last four columns $\bm{A}$ and $\bm{B}$ are

[TABLE]

It is evident that $\xi^{(4)}$ for the first problem and $\xi^{(1)}$ for the second problem correspond to fictitious intersection points while the other three eigenvalues correspond to the true intersection points, see Figure 2. The coordinates of the three true intersection points are computed by introducing $\xi^{(j)}$ in the line equation (20) yielding

[TABLE]

The fictitious points can also be detected without computing several eigenvalue problems and comparing their eigenvalues; that is, different from above, by computing only one single eigenvalue problem. This is accomplished by determining the parametric coordinates of the intersection points $\theta^{(j)}$ on the curve $\bm{x}(\theta)$. As indicated in (17), these parametric coordinates are computed using the eigenvectors. As an example, consider the eigenvalue problem (25) defined by first four columns of $\bm{A}$ and $\bm{B}$. Its eigenvalues are given in (26) and the coordinates of the intersection points are

[TABLE]

The respective parameters $\theta^{(j)}$ of the intersection points are

[TABLE]

A point is a true intersection point if and only if $\bm{x}(\theta^{(j)})=\bm{x}^{(j)}$; otherwise, it is a fictitious point. It can easily be found $\bm{x}(\theta^{(4)})=(0.1506,\,0.3949)\neq\bm{x}^{(4)}$ such that the fourth point is not an intersection point. Thus, in practice the intersection points are computed from a single couple of square matrices $\bm{A}^{\Box}$ and $\bm{B}^{\Box}$.

3 Intersection of lines with surfaces

The extension of the introduced method to surfaces is straightforward. Let $\bm{x}(\bm{\theta})=(x^{1}(\bm{\theta}),\,x^{2}(\bm{\theta}),\,x^{3}(\bm{\theta}))^{\mathsf{T}}$ be a parametric surface, with $\bm{x}(\bm{\theta})\in\mathbb{R}^{3}$, of bi-degree $(q_{x}^{1},\,q_{x}^{2})$ given either in Lagrange basis $L_{\bm{i}}(\theta^{1},\,\theta^{2})$ or power basis $P_{\bm{j}}(\theta^{1},\,\theta^{2})$ with

[TABLE]

where $\bm{i}=(i^{1},\,i^{2})$ and $\bm{j}=(j^{1},\,j^{2})$ are multi-indices, $\bm{\theta}=(\theta^{1},\,\theta^{2})$ are the parametric surface coordinates, and $\bm{x}_{\bm{i}}\in\mathbb{R}^{3}$ and $\bm{\alpha}_{\bm{j}}\in\mathbb{R}^{3}$ are the coefficients in the two basis. Usually, in finite element applications the degrees $q_{x}^{1}$ and $q_{x}^{2}$ of the surface $\bm{x}(\bm{\theta})$ are the same.

In line with the curve case, a point $\bm{x}(\bm{\theta})$ on the surface is defined as the intersection of several moving planes of the form

[TABLE]

where $\bm{g}(\bm{\theta})=\left(g^{1}(\bm{\theta}),\,g^{2}(\bm{\theta}),\,g^{3}(\bm{\theta}),\,g^{4}(\bm{\theta})\right)^{\mathsf{T}}$ is an auxiliary vector collecting the parameters of the plane. Although $l(\bm{\theta},\,\bm{x})=0$ is now a plane instead of a line, it is still denoted with the same symbol to keep the notation simple. The parameters $\bm{g}(\bm{\theta})$ are assumed to be of the following form

[TABLE]

The bi-degree $(q_{g}^{1},\,q_{g}^{2})$ of the power basis $\widetilde{P}(\bm{\theta})$ has to be sufficiently high to describe all the intersection points of the surface $\bm{x}(\bm{\theta})$ with a line.

The planes describing the surface $\bm{x}(\bm{\theta})$ have to satisfy

[TABLE]

which can be rearranged to

[TABLE]

with a new power basis $\widehat{P}_{\bm{k}}$ of bi-degree $(q_{x}^{1}+q_{g}^{1}+1,\,q_{x}^{2}+q_{g}^{2}+1)$ and the array $\bm{h}$ containing the sorted components of the vectors $\bm{g}_{\bm{l}}$. The matrix $C_{\bm{k}\bm{l}}$ has $(q_{x}^{1}+q_{g}^{1}+1)(q_{x}^{2}+q_{g}^{2}+1)$ rows and $4(q_{g}^{1}+1)(q_{g}^{2}+1)$ columns. It has at least $4(q_{g}^{1}+1)(q_{g}^{2}+1)-(q_{x}^{1}+q_{g}^{1}+1)(q_{x}^{2}+q_{g}^{2}+1)$ right null vectors, i.e. the difference between the number of columns and rows. The surface $\bm{x}(\bm{\theta})$ has the algebraic degree $2q_{x}^{1}q_{x}^{2}$, which is equal to its number of intersections with a line. In order to obtain all the intersections the following condition has to be satisfied

[TABLE]

As in the curve case, the degree of $\widetilde{P}(\bm{\theta})$ along the $\theta^{2}$ direction can be chosen as $q_{g}^{2}\geq q_{x}^{2}-1$ yielding

[TABLE]

By symmetry, it is also valid to choose $q_{g}^{1}\geq q_{x}^{1}-1$ and $q_{g}^{2}\geq 2q_{x}^{2}-1$. After the right null vectors of $C_{\bm{k}\bm{l}}$ are computed the subsequent steps in computing the intersections are identical to the curve case.

Finally, for the sake of completeness, we mention that the introduced intersection algorithm also applies to triangular finite elements. If $\bm{x}(\bm{\theta})$ is a triangular parametric surface of degree $q_{x}$ then the degree $q_{g}$ of the auxiliary vector $\bm{g}(\bm{\theta})$ has to be chosen to satisfy $q_{g}\geq 2(q_{x}-1)$; see also [9, §3].

4 Conclusions

The introduced method is able to determine in one-shot, simultaneously without iterating, all the intersection points between a line and a curve or surface. It can be applied to curves or surfaces given in any polynomial basis, like the Lagrange or Bernstein, after a straightforward conversion to the power basis. Aspects requiring further research include cases for which the eigenvalue problem (18) is degenerate and the preconditioning of the eigenvalue problem (25). In the accompanying implementation, several choices have been considered for both.