Mode Surfaces of Symmetric Tensor Fields: Topological Analysis and Seamless Extraction
©© 2020 IEEE. This is the author’s version of the article that has been published in IEEE Transactions on Visualization andComputer Graphics. The final version of this record is available at: xx.xxxx/TVCG.201x.xxxxxxx/
Mode Surfaces of Symmetric Tensor Fields:Topological Analysis and Seamless Extraction
Botong Qu, Lawrence Roy, Yue Zhang,
Member, IEEE , and Eugene Zhang,
Senior Member, IEEE (a) µ = − . and µ = ± . (b) µ = − . (c) µ = . and µ = − . Fig. 1: Mode surfaces of a stress tensor field for a block under compression. When the mode value µ is close to ± µ is close 0(c), the mode surface converges onto the neutral surfaces (chartreuse). By observing the change in mode surfaces’ geometry andtopology ((a)-(c)), we gain insight into the interaction between degenerate curves and neutral surfaces, the two constituents of tensorfield topology. In addition, some mode surfaces can contain interesting features not present in degenerate curves and neutral surfaces,such as the bottom layer with a hole (b). Finally, the change in mode surfaces’ topology when mode values change, such as thecontraction of the mode surface from (b) to a vascular structure in (a) can also provide important insight into the underlying physics. Abstract —Mode surfaces are the generalization of degenerate curves and neutral surfaces, which constitute 3D symmetric tensor fieldtopology. Efficient analysis and visualization of mode surfaces can provide additional insight into not only degenerate curves and neutralsurfaces, but also how these features transition into each other. Moreover, the geometry and topology of mode surfaces can helpdomain scientists better understand the tensor fields in their applications. Existing mode surface extraction methods can miss featuresin the surfaces. Moreover, the mode surfaces extracted from neighboring cells have gaps, which make their subsequent analysisdifficult. In this paper, we provide novel analysis on the topological structures of mode surfaces, including a common parameterizationof all mode surfaces of a tensor field using 2D asymmetric tensors. This allows us to not only better understand the structures in modesurfaces and their interactions with degenerate curves and neutral surfaces, but also develop an efficient algorithm to seamlesslyextract mode surfaces, including neutral surfaces. The seamless mode surfaces enable efficient analysis of their geometric structures,such as the principal curvature directions. We apply our analysis and visualization to a number of solid mechanics data sets.
Index Terms —Tensor field visualization, tensor field topology, traceless tensors, degenerate curve extraction, neutral surface extraction,mode surface extraction.
NTRODUCTION • B. Qu is with the School of Electrical Engineering and Computer Science,Oregon State University. E-mail: [email protected].• L. Roy is with the School of Electrical Engineering and Computer Science,Oregon State University. E-mail: [email protected].• Y. Zhang is an Associate Professor with the School of Electrical Engineeringand Computer Science, Oregon State University. E-mail:[email protected].• E. Zhang is a Professor with the School of Electrical Engineering andComputer Science, Oregon State University. E-mail:[email protected] received xx xxx. 201x; accepted xx xxx. 201x. Date of Publicationxx xxx. 201x; date of current version xx xxx. 201x. For information onobtaining reprints of this article, please send e-mail to: [email protected] Object Identifier: xx.xxxx/TVCG.201x.xxxxxxx
Symmetric tensor fields have a wide range of applications in science, en-gineering, and medical domains. The diffusion tensor field analysis inmedical imaging plays a key role in diagnosing and treatment planningfor brain cancers. The stress and strain tensors in continuum mechanicsenable the predictions of structural failures. Topology-driven analy-sis and visualization of 3D symmetric tensor fields have made muchprogress in recent years, of which the focus is on high-quality extractionand visualization of the tensor field topology such as degenerate curves (where the tensors have repeating eigenvalues) and neutral surfaces (where the major, medium, and minor eigenvalues of the tensors forman arithmetic sequence). On the other hand, degenerate curves andneutral surfaces are often treated as unrelated objects and interpretedseparately.In fact, both are a level set of the mode function of the tensor field [6],which ranges from − ± a r X i v : . [ c s . G R ] S e p et. In solid mechanics [6], a mode − uniaxialcompression and a mode 1 tensor corresponds to uniaxial extension . Incontrast, a mode 0 tensor corresponds to pure shear . At other modevalues, we observe biaxial extension and compression.Degenerate curves and neutral surfaces transition into each other, andunderstanding their interaction can provide more insight than interpret-ing them separately. Such an interaction can be better understood andvisualized through level set surfaces whose iso-values are between 0and ±
1. These surfaces are referred to as mode surfaces such as thegold-colored surfaces shown in Figure 1. Note how mode surfaces withmode values close to ± (a) A-patches Method [27](b) Our method Fig. 2: The holes and gaps in mode surfaces extracted using the A-patches method [27] lead to a mesh with seams, impeding the computa-tion of important surface properties such as normal and curvature (top:LIC texture lacks clarity in showing the major principal direction in themode surface). Our method generates seamless meshes, which lead torobust curvature computation (bottom).In this paper, we address the aforementioned challenges by providing aparameterization for all mode surfaces given a 3D linear tensor field.This common parameterization allows us to extract mode surfaces atany accuracy without the need for mesh subdivision. This leads tofewer missing pieces in the extracted surfaces than from the A-patchesmethod. Furthermore, our algorithm extracts mode surfaces inside eachface in the mesh, which are then used to find mode surfaces insideeach tetrahedron. This makes it straightforward to stitch mode surfacesfrom adjacent tetrahedra without gaps, resulting in a seamless modesurface on which differential properties can be computed and visualized(Figure 2 (b)). Moreover, a mode surface can be extracted under fiveseconds for our simulation data sets.Our pipeline also applies to the extraction of neutral surfaces, whichis a non-orientable surface [2] characterized by a degree-three poly-nomial [32]. To our knowledge, the algorithm and the pipeline wedevelop for mode and neutral surface extraction, in conjunction withthe degenerate curve extraction method of Roy et al. [32], is the first unified framework for all isosurfaces of the mode function, includingneutral surfaces (mode 0) and degenerate curves (mode ± (a) The A-patches method [27] (b) Our method(c) The method of [32] (d) Our method Fig. 3: Existing mode surface extraction methods such as [27] (top-left)can have holes due to non-convergence in the extraction process. Inaddition, these methods [27, 32] can have gaps between mode surfacesextracted from different tets (top-left and bottom-left). These holesand gaps can not only mislead interpretation but also make it difficultto compute differential properties of the mode surfaces such as theircurvature tensors. Our method addresses these challenges with a unifiedframework in which mode surfaces (including neutral surfaces) canbe extracted more accurately, faster, and in a seamless fashion (rightcolumn).To demonstrate the utility of our approach, we apply our tensor fieldanalysis and visualization to solid mechanics applications.
ELATED W ORK
Tensor field visualization has advanced much in the last decades [4, 19].Delmarcelle and Hesselink [9] introduce the notion of topology for2D symmetric tensor fields, which consists of degenerate points. Thetopological features of 3D symmetric tensor fields are first studied byHesselink et al. [15]. Zheng and Pang [37] point out that degeneratepoints form curves under structurally stable conditions, i.e. the structurepersists under arbitrarily small perturbations [7]. Several methods havebeen proposed to extract degenerate curves [27, 32, 35, 39]. Palacios etal. [26] introduce editing operations for degenerate curves, such as de-generate curve removal, degenerate curve deformation, and degeneratecurve reconnection.Besides feature curves, another type of topological features is surfaces.Zobel and Scheuermann [40] introduce the notion of extremal sur-faces for 3D symmetric tensor fields. Raith et al. [30] extract fibersurfaces of tensor fields by linearly interpolating tensor invariants ineach tetrahedron. Palacios et al. [27] introduce the notion of neutralsurfaces.Extracting implicit surfaces is a well-researched area [8]. The mostpopular technique, Marching Cubes [20] and its variants, focus ontrilinear functions that are degree-three polynomials, while mode sur-faces are of degree-six. Using surface extraction methods designedfor a lower-degree polynomial, even with a guarantee for topologi-cal correctness [13, 24, 31], can still miss important topological andgeometric features of isosurfaces corresponding to a higher-degreepolynomial. Implicit surface extraction methods that can handle more general functions with topological guarantees [3, 25, 34] usually re-quire C functions. Since the mode function in our case is a piecewisedegree-six polynomial and C at the cell boundaries, it is not clear howto adapt these techniques to the mode function.The A-patches method [21] is a technique to extract algebraic surfaces,which is adapted by Palacios et al. [27] to extract mode surfaces. How-ever, to our knowledge, the extracted surfaces can contain seams at thefaces of the mesh.Roy et al. [32] provide algorithms to extract degenerate curves andneutral surfaces from linear tensor fields based on parameterizationsof features in the field by their eigenvectors (medium eigenvectorsfor neutral surfaces and dominant eigenvectors for degenerate curves).These parameterizations lead to a more accurate extraction of degener-ate curves and neutral surfaces than previous techniques based on theA-patches algorithm [27]. However, while their extraction of degener-ate curves is seamless, their extraction of neutral surfaces is not. In thispaper, we provide a unified parameterization of all mode surfaces ofa linear tensor field, which enables a unified pipeline for the seamlessextraction of mode surfaces including neutral surfaces and degeneratecurves.Our analysis makes use results from 2D asymmetric tensor fields.Zheng and Pang [38] introduce the notion of real domains and complexdomains for 2D asymmetric tensor fields. Zhang et al. [36] providetopological analysis of asymmetric tensor fields on surfaces, whichKhan et al. [18] extend to a multi-scale framework. Chen et al. [29]visualize asymmetric tensor fields on surfaces with a hybrid approach:glyphs for the complex domain and hyperstreamlines for the real do-main. ENSOR B ACKGROUND
In this section, we review the relevant math background on tensors andproperties of 3D linear symmetric tensor fields. An n -dimensional tensor T can be expressed as an n × n matrix under agiven orthonormal basis.The trace of a tensor T = ( T i j ) is the sum of its diagonal elements.When the trace is zero, the tensor is referred to as being traceless .A tensor T can be uniquely decomposed as the sum of the tensor D = trace Tn I (a multiple of the identity matrix) and a traceless tensor A = T − D (referred to as the deviator of T ). Note that T and A havethe same set of eigenvectors. The set of all n × n tensors form a linearspace, on which the following inner product of two tensors R and S canbe introduced [33]: (cid:104) R , S (cid:105) = n ∑ i = n ∑ j = R i j S i j = trace ( S T R ) . (1).With this product, one can define the magnitude of a tensor T as || T || = (cid:112) (cid:104) T , T (cid:105) . Another important quantity of a given tensor is its determinant | T | , which is the product of its eigenvalues.A tensor T is symmetric if it is equal to its transpose. Otherwise, itis asymmetric . The eigenvalues of a symmetric tensor are guaranteedto be real-valued, while the eigenvalues of an asymmetric tensor canbe either real-valued or complex-valued. Furthermore, the eigenvec-tors belonging to different eigenvalues of a symmetric tensor forman orthonormal basis. For asymmetric tensors, even when the eigen-3alues are real-valued, their respective eigenvectors are not mutuallyperpendicular.Given our focus on 3D symmetric tensors and occasional mention of 2Dasymmetric tensors, in the remainder of the paper we will drop the word“symmetric” for symmetric tensors and keep the word “asymmetric”for asymmetric tensors. Moreover, we only consider 3D traceless(symmetric) tensors and therefore will also omit the word “traceless”for 3D tensors. In contrast, when discussing 2D asymmetric tensors,we do not assume that they are traceless. A 3 × T has three eigenvalues λ ≥ λ ≥ λ , which are re-ferred to respectively as its major eigenvalue , medium eigenvalue , and minor eigenvalue . Eigenvectors corresponding to T ’s major eigenvalueare referred to its major eigenvectors . We can define T ’s mediumeigenvectors and minor eigenvectors in a similar fashion. T is degenerate if it has repeating eigenvalues. Under structurally sta-ble conditions, a degenerate tensor T has two eigenvalues being thesame (referred to as the repeating eigenvalue ). The third eigenvalueis the dominant eigenvalue . Furthermore, if the dominant eigenvalueis larger than the repeating eigenvalue, T is referred to as being lineardegenerate . If the dominant eigenvalue is smaller than the repeatingeigenvalue, T is referred to as being planar degenerate . The eigenvec-tors corresponding to the dominant eigennvalue are referred to as the dominant eigenvectors .A 3 × T is neutral if its medium eigenvalue is the averageof its major and minor eigenvalues. The dominant eigenvalue andeigenvectors are not well-defined for neutral tensors.The mode of a 3D (traceless, symmetric) tensor T is µ ( T ) = √ det ( T ) (cid:107) T (cid:107) ,with a range of [ − , ] . As special instances, neutral tensors are mode0 tensors, while linear degenerate tensors and planar degenerate tensorscorrespond to mode 1 and mode − µ as follows [6, 16, 23]: λ = (cid:114)
23 sin (
13 arcsin ( − µ ) + π ) , λ = (cid:114)
23 sin (
13 arcsin ( − µ )) , λ = (cid:114)
23 sin (
13 arcsin ( − µ ) − π ) . (2)A 3D tensor field is a continuous tensor-valued function. A degeneratepoint and a neutral point are where the tensor values are degenerate and neutral , respectively. Under structurally stable conditions, degeneratepoints form curves ( degenerate curves ) [37], and neutral points formsurfaces ( neutral surfaces ) [27]. In general, the µ level set of the modefunction is a surface when − < µ <
1. Such a level set is referred toas a mode - µ surface [27]. Note that both degenerate curves and neutralsurfaces are special level sets of the mode µ . We focus on 3D linear tensor fields, which can be written in the formof T ( x , y , z ) = T + xT x + yT y + zT z where T , T x , T y , and T z are linearlyindependent 3D tensors. Let U be the set of 3D (traceless, symmetric)tensors, which is a five-dimensional space. Under structurally stableconditions, there exists a 3D tensor T that satisfies the following: (cid:104) T , T (cid:105) = (cid:104) T , T x (cid:105) = (cid:104) T , T y (cid:105) = (cid:104) T , T z (cid:105) = , (cid:104) T , T (cid:105) = , det ( T ) ≤ . (3)Note that T plays an important role in the behavior of the tensorfield [32]. We refer to T as the characteristic tensor of the lineartensor field.The set of degenerate points of a 3D linear tensor field can be parameter-ized by a topological circle [32]. That is, the union of the set of mode 1points and mode − RP (referred to as the medium eigenvector manifold ) except for two lines,each of which corresponds to a single point in the medium eigenvectormanifold (referred to as a singularity in the parameterization). Dueto the existence of the two singularities, the set of neutral points of a3D linear tensor field is homeomorphic to RP attached with a handle(thus non-orientable) [2].In this paper, we provide analysis on the topology of mode surfaces aswell as efficient algorithms to extract them. ODE S URFACES A NALYSIS
Given a number µ ∈ ( − , ) (cid:83) ( , ) , we wish to seamlessly extractthe mode µ surface from a piecewise linear tensor field defined on atetrahedral mesh. To do so, we provide a unified framework in the samespirit of the degenerate curve extraction method of Roy et al. [32]. Thatis, we first extract mode curves on each triangular face in the mesh.Next, we extend the mode curves from the four faces of each tet toextract the mode surface inside the tet. Finally, we stitch the modesurface from adjacent tets across their common faces to generate aseamless mode surface in the whole mesh.This requires the ability to extract mode surfaces inside a tet at high-quality. Recall that inside each tet of the tetrahedral mesh, the tensorfield is linear. In the remainder of this section, we will describe ournovel analysis of mode surfaces for 3D linear tensor fields, whichleads to a parameterization of such surfaces that enables high-qualityextraction. We will state the results of our analysis in the paper andprovide their proofs in Appendix A.Similar to the case of degenerate curves [32], we will consider the set ofmode ± µ points together for our mode surface analysis and extraction.These points satisfy the following degree-six equation:54 ( det ( T )) − µ ( T ) (cid:107) T (cid:107) = , (4)and we refer to the collection of such points as the generalized mode µ surface . As in the case of neutral surfaces, we show that a generalizedmode µ surface can also be parameterized by its medium eigenvectors(Theorem 1 in Appendix A). This parameterization is based on thefollowing 2D asymmetric tensor field defined on the unit sphere: A ( v ) = R θ + π T (cid:48) ( v ) R θ − π (5)in which v is a unit vector, T (cid:48) ( v ) is the projection of T onto theplane whose normal is v , θ = arcsin ( √ ( arcsin ( µ ))) , and R φ = (cid:18) cos φ − sin φ sin φ cos φ (cid:19) .4 Figure 4 illustrates the asymmetric tensor field A with an example 3Dlinear tensor field created manually. There are four types of regionson the aforementioned sphere (the medium eigenvector manifold): (1)grey, (2) cyan, (3) magenta, and (4) blue. (a) Medium eigenvector manifold (b) Generalized mode µ surfaces Fig. 4: A generalized mode µ surface (right) can be parameterized bythe medium eigenvector manifold (left). Each pair of antipodal points inthe blue region of the medium eigenvector manifold (left: red dot) givesrise to two points in the generalized mode µ surfaces (right: red dots),one with a positive mode value (right: the red dot on the teal surface)and one with a negative mode value (right: the red dot on the goldsurface). The dominant eigenvector directions at these points (right:shown with LIC textures in the planes) are given by the eigenvectors ofan asymmetric tensor field in the medium eigenvector manifold (left:LIC texture directions). Points in the cyan and magenta regions in themedium eigenvector manifold correspond to two positive-mode pointsand two negative-mode points, respectively. Points in the grey regiondo not correspond to any point in the generalized mode µ surface.The grey region is the complex domain of the asymmetric tensor field A ( v ) , i.e. with complex eigenvalues. Unit vectors in this region cannotappear as the medium eigenvector of any tensor in the 3D linear tensorfield. That is, there are no points in the generalized mode µ surfacethat correspond to these unit vectors. Note that if a unit vector v is inthe complex domain, so is − v . Therefore, the complex domain of theasymmetric field respects the antipodal symmetry.The cyan, magenta, and blue regions together form the real domain of the asymmetric tensor field. Each pair of antipodal points in themedium eigenvector manifold inside these regions correspond to twopoints in the generalized mode µ surface. If both points have thepositive mode µ , we color the original pair in the sphere with cyan. Ifboth points have the mode − µ , we color the pair in the sphere magenta.If one point has the mode µ and other − µ , we color the correspondingpair in the sphere blue.Note that the major and minor eigenvectors of the asymmetric tensorfield A give rise to the dominant eigenvectors of the 3D tensor fieldat the corresponding points in the generalized mode µ surface. Anexample is shown in Figure 4. In (a), a unit vector v in the sphere(highlighted by a red dot) corresponds to two points in the generalizedmode µ surface (b), each of which also highlighted by a red dot. Thepoint in the teal surface (b) has a positive mode while the point in thegold surface has a negative mode. The planes normal to the mediumeigenvectors are shown for both points (b). Notice that the two planesare parallel, i.e. the medium eigenvectors are the same for both points,which correspond to the vector in the medium eigenvector manifold(the red dot in (a)). The dominant eigenvector directions at the pointsin the generalized mode µ surface ((b): the LIC textures in the twoplanes) together match the major and minor eigenvector directions of A ( v ) ((a): the LIC directions at the red dot).Due to the symmetry in the tensor field, − v corresponds to the sametwo points in the generalized mode µ surface (Theorem 1 in Ap- pendix A). However, the minor eigenvector of A ( v ) becomes the majoreigenvector of A ( − v ) and the major eigenvector of A ( v ) becomes theminor eigenvector of A ( − v ) . To make our parameterization easier forsubsequent processing, we choose the unit vector from v and − v sothat its major eigenvector gives rise to the dominant eigenvector of thecorresponding point in the generalized mode µ surface. This schemeremoves the ambiguity in our parameterization by converting the two-to-two correspondence ( ± v to the two points in the generalized mode µ surface) to a one-to-one correspondence.Finally, points on the boundary between the real domain and com-plex domain ( complex domain boundary ) have one real eigenvalue ofmultiplicity of two. The complex domain boundary also satisfies theantipodal symmetry. That is, if a vector v is in the complex domainboundary, so is − v . Moreover, v and − v correspond to exactly onepoint in the generalized mode µ surface.Given a 3D linear tensor field, all of its generalized mode µ surface canbe parameterized by the same sphere. The following equation charac-terizes the complex domain boundary (Theorem 2 in Appendix A):12 − v T T v +
14 cos θ ( v T T v ) = θ = arcsin ( √ ( arcsin ( µ ))) is the same as that in Equa-tion 5. We illustrate the changes in the geometry and topology ofgeneralized mode µ surfaces with an example tensor field in Figure 5.Three generalized mode µ surfaces (teal and gold surfaces) are shownin (a), (b), and (c), respectively. Their corresponding eigenvector man-ifolds are also shown (the three spheres). In addition, we show thedegenerate curves (the yellow and green curves in (a)) and the neutralsurfaces (the chartreuse surface in (c)).When µ =
1, the mode surface is essentially the degenerate curves. Thecomplex domain in the corresponding medium eigenvector manifold(not shown) consists of two connected components, satisfying theantipodal symmetry. When µ decreases, the complex domain shrinksin size (a). At this stage, the complex domain boundary still consistsof two curves satisfying the antipodal symmetry (one is visible in (a)).Since each pair of antipodal points on the complex domain boundarycorresponds to exactly one point in the generalized mode µ surface,the latter can be topologically constructed by removing the complexdomains from the sphere (punching two holes) and gluing the surfacealong the complex domain boundary based on the antipodal symmetry.This results in a torus.Note that as µ decreases, the complex domain grows smaller (Corol-lary 3 in Appendix A). When µ reaches (cid:113) − µ where µ is the modeof T , the complex domain boundary touches itself (Figure 5 (b)). Thegeneralized mode µ surface at this point is a non-manifold.After this, the complex domain continues to shrink and now has fourconnected components. Performing the same topological surgery (re-moving the complex domains and gluing along the complex domainboundary) results in a double-torus, i.e. a surface with two handles(Figure 5 (c)).Finally, when µ =
0, the complex domain disappears, and the complexdomain boundary degenerates into two pairs of antipodal points in themedium eigenvector manifold. They are precisely the singularities inthe parameterization [32], each of which corresponds to a straight linein the neutral surface. Using the eigenvectors of T as the coordinatesystem for the medium eigenvector manifold, we have the singularitiesbeing (cid:18) ± (cid:114) λ − λ λ − λ ± (cid:114) λ − λ λ − λ (cid:19) .We refer the reader to Theorem 5 in Appendix A for the proof of the5ig. 5: Given a 3D linear tensor field, the positive mode surfaces (teal) converge to linear degenerate curves (green) and the negative modesurfaces (gold) converge to planar degenerate curves (yellow) when µ approaches ± µ surface has a topology of a torus. In contrast, when µ approaches 0 (c), the positive mode and negative mode surfaces together converge onto the neutral surface (chartreuse). The complex domain nowconsists of two pairs of antipodal loops, and the generalized mode µ surface has a topology of a double-torus. The bifurcation between the twocases occurs (b) when the complex domain resembles a figure-eight and the generalized mode µ surface has a non-manifold point.above analysis. Note that neutral surfaces are the only non-orientablemode surfaces, which highlights their topological significance.The complex domain boundary is characterized by a degree-six polyno-mial, which makes its extraction challenging. Fortunately, similar todegenerate curves in a 3D linear tensor field which can be character-ized by an elliptical loop, the complex domain boundary can also beparameterized as illustrated in Figure 6. Due to the four-fold symmetryin the complex domain boundary, it is sufficient to focus on one quarterof the curve, which, in the coordinate system α = v T T v , (7) β = v T T v (8)is characterized by 0 = − β +
12 cos θ α . (9)This is the equation of a parabola, which can be parameterized by α (Lemma 4 in Appendix A). Each α gives one corresponding β fromEquation 9 and thus one point in each of the four quarters of the complexdomain boundary (Figure 6). XTRACTION OF M ODE S URFACES
In this section, we describe our mode surface extraction algorithm, theinput of which is a tetrahedral mesh, at whose vertices tensors are given.These tensor values are linearly interpolated into the faces and interiorsof each tet. This results in a piecewise linear 3D tensor field, i.e. insideeach tet, the tensor field is linear. Similarly, inside each face, we have a2D linear tensor field.The pipeline of our unified mode surface (including neutral surfaces)extraction is illustrated in Figure 7 with an example tensor field. First,we extract mode curves inside each face in the mesh, including internalloops. Next, for each tet with known mode curves in any of its faces, weextract the mode surface inside the tet. Finally, mode surfaces extractedfrom the tets will be stitched together along shared faces. We describethe detail of each step next. Fig. 6: The complex domain boundary has a four-fold symmetry, witheach quarter parameterizable. A quarter segment can intersect the hori-zontal axis and the vertical axis at one point each, making the complexdomain a single region (the outermost loop). This corresponds to thecase shown in Figure 5 (a). In contrast, a quarter segment can also in-tersect the horizontal axis only (inner-most case), splitting the complexdomain into two regions (the innermost loop). This corresponds to thecase shown in Figure 5 (c). The bifurcation occurs when the quartersegments pass through the origin (the loop between the outermost andinnermost loops), which corresponds to Figure 5 (b) .
Mode curves inside a plane can consist of a number of open curvesand loops. If an open curve does not intersect any edge of a trianglein our mesh, the curve must be entirely outside the triangle and is notpart of mode curves for the triangle. In contrast, there exist loops thatare entirely inside the triangle (thus valid) but do not intersect any ofthe edges of the triangle. Thus, detecting mode curves for a triangle byonly detecting their intersections with the boundary edges can miss theinner loops.To overcome this problem, we use a property from Morse theory [22],which states that for any level set loop of a function there must be atleast one local maximum or minimum of that function enclosed by theloop. Based on this insight, we first extract the critical points of themode function inside the plane containing the triangle. Next, for eachextremum (a critical point that is not a saddle) inside the triangle, weinsert a new vertex at the point and subdivide the triangle into threetriangles. In Theorem 6 (Appendix A) we show that there are at mostfour extrema inside a plane P containing the triangle. Only the ones that6 Fig. 7: The pipeline of our mode extraction algorithm for a given tet:starting from mode curves extracted from the faces of a tet (a), weidentify the region in the medium eigenvector manifold bounded bythe aforementioned curves (b). The region is then triangulated (c) andmapped back to the
XY Z space to give the mode surface in the tet (d).are inside the triangle are valid. Thus, the aforementioned subdivisionprocess is executed at most four times for each face in the mesh. Tofind the critical points of the mode function inside a plane P , we makeuse of the fact that they must be solutions to the following system ofpolynomial equations: vh ( u , v , w ) = wg ( u , v , w ) , (10) w f ( u , v , w ) = uh ( u , v , w ) , (11) u + v + w = ∇ det ( uT + vT + wT ) = ( f ( u , v , w ) , g ( u , v , w ) , h ( u , v , w )) . Here, T , T , and T is an orthonormal basis for the set of tensor values in-side P given by the tensor field. In addition, ( u , v , w ) is a unit vectorsuch that T ( x , y ) || T ( x , y ) || = u ( x , y ) T + v ( x , y ) T + w ( x , y ) T . Note that Equa-tions 10 and 11 are both homogeneous cubic polynomials. Moreover,if ( u , v , w ) is a solution, so is ( − u , − v , − w ) . Furthermore, both ( u , v , w ) and ( − u , − v , − w ) correspond to the same tensor, i.e. the same criticalpoint. Let u (cid:48) = uw and v (cid:48) = vw . The original system of equations istransformed into the following system of two cubic equations: v (cid:48) h ( u (cid:48) , v (cid:48) ) = g ( u (cid:48) , v (cid:48) ) , f ( u (cid:48) , v (cid:48) ) = u (cid:48) h ( u (cid:48) , v (cid:48) ) (13)where f , g , and h are non-homogeneous cubic polynomials derivedrespectively from f , g , and h with the change from u , v , w to u (cid:48) , v (cid:48) .According to B´ezout’s theorem [10], this system of cubic equationsis equivalent to a degree-nine polynomial which has nine solutions.However, we show that two of the solutions are spurious and four othersolutions correspond to degenerate points in the plane P . Moreover, thespurious solutions can be found by solving h ( u , v , ) = u + v =
1, while the degenerate points can be found using the method of Royet al. [32]. Once these solutions are factored out, we obtain a cubicpolynomial whose roots, when real-valued, are additional critical points(local extrema and saddles). We compute these roots using the Eigen library [14]. For each solution in the form of u , v , and w values, werecover T = uT + vT + wT which we use to find ( x , y ) such that T ( x , y ) = T .We can extract the intersection points of the mode curves with theedges in the triangles (including the subdivided triangles). Startingfrom these intersection points, we perform numerical tracing in thedirection perpendicular to the gradient of the mode function. Thisguarantees that all internal loops in the mode curves are found. To traceout the face intersections we use a numerical ODE integrator. Thisleads to mode curves whose ends are on the edges of the face. Tracingis finished when we exit the (possibly subdivided) triangle. We thenrun a numerical root-finding algorithm [1] on the interpolation functionproduced by the ODE integrator to accurately find the point where itexits, which is connected to the closest edge intersection point. Once we have extracted mode curves from all the faces in the mesh, weproceed to extract mode surfaces from inside each tet.First, we gather the mode curves extracted from the four faces of thetet and stitch open curves (those intersecting the edges of the tet) intoloops (Figure 7 (a): blue curves).Next, for each sample point on the aforementioned mode curves, wefind the corresponding point on the medium eigenvector manifoldfor the 3D linear tensor field in the tet. Therefore, the mode curves(represented as polylines) on the boundary of the tet lead to a setof curves (also represented as polylines) in the medium eigenvectormanifold (Figure 7 (b): blue curves).Recall that a pair of antipodal points on the complex domain bound-ary corresponds to the same point in the mode surface. Therefore, aconnected component in the polyline in the tet’s boundary can becomedisconnected in the medium eigenvector manifold at the complex do-main boundary (green points in Figure 7 (a) and (b)). To overcome thisproblem, we need to not only find the exact intersection points of thepolylines with the complex domain boundary but also generate curvesconnecting such points in the medium eigenvector manifold, whichcorrespond to segments of the complex domain boundary (Figure 7 (b):green curves).Detecting the existence of such points is relatively straightforwardas each such point corresponds to an antipodal pair in the mediumeigenvector manifold. Therefore, along a mode curve, such a pointdivides the curve into two disconnected components, whose mediumeigenvectors are situated on opposite hemispheres of the medium eigen-vector manifold. Consequently, for every segment in the polyline, wecheck the dot product between the medium eigenvectors of the two endpoints. If the sign of the dot product is negative, we mark the segmentas intersecting the complex domain boundary and proceed to find theexact location of the intersection point through numerical root-findingusing Equation 6. This segment is then split into two with the insertionof the newly found intersection point.Once we have found all the complex domain boundary points in thefaces of the tet, we compute their corresponding antipodal point pairson the medium eigenvector manifold. To decide which segments on thecomplex domain boundary are part of the mode surface, we make useof the parameterization for the complex domain boundary (Figure 6).This parameterization allows us to select a new point on the complexdomain boundary and check whether its corresponding point on themode surface is inside or outside the tet. This information is then usedto decide which segments along the complex domain boundary are partof the mode surfaces inside the tet. Notice that our method is similar tothe method of Roy et al. [32] for computing degenerate curves, whichare parameterizable by an ellipse.7he segments along the complex domain boundary (Figure 7 (b): greencurves) and the polylines representing the mode curves from the faces(Figure 7 (b): blue curves) bound the region in the medium eigenvectormanifold that corresponds to the mode surface inside the tet (Figure 7(b): yellow region). To find the interior of this region, we perform aconstrained Delaunay triangulation on the set of sample points (frommode curves and complex domain boundary segments). Since Delaunaytriangulation is normally defined on a plane, we take the stereographicprojection of the vertices as the stereographic projection maps circlesto circles and thus preserves the Delaunay condition. This leads toa triangulation of the region with rather poor aspect ratios for thetriangles.To improve the quality of the tessellation, we add more points insidethe region on the medium eigenvector manifold. Given our unifiedapproach of using the medium eigenvector manifold for mode surfaceextraction (including neutral surfaces), we make use of the followingquad parameterization of the medium eigenvector manifold based onthe singularities in the neutral surfaces. x = cos θ (cid:113) f cos η + , y = cos θ sin η , z = cos η , (14)where f = λ ( T ) − λ ( T ) − λ ( T ) λ ( T ) − λ ( T ) , 0 ≤ η < π , and 0 ≤ θ < π . Thisparameterization leads to a perfect quadrangulation of the mediumeigenvector manifold (Figure 7 (b): the quad grid) with the only irregu-lar vertices in the quadrangulation being the singularities in the mediumeigenvector manifold (Figure 7 (b): yellow dots). After the grid pointsinside the region are added, we perform a second constrained Delaunaytriangulation with both boundary and interior sample points. This leadsto improved aspect ratios in the triangulation (Figure 7 (c)).Finally, we need to map the sample points from the medium eigenvec-tor manifold back to the XY Z space while preserving the connectivityamong them to construct the mode surfaces inside the tetrahedron. Tomap a point v in the medium eigenvector manifold to its correspondingpoint in the XY Z space, we need to identify the tensor T whose mediumeigenvector is v . The eigenvalues of T can be computed using Equa-tion 2 given the mode value µ . The dominant eigenvector of T can becomputed from the asymmetric tensor A ( v ) (Equation 5). This givesus the tensor T . We can find its corresponding point in the XY Z spaceby solving a system of linear equations [32].
To get a seamless mode surface, we need to stitch the sheets fromdifferent tets across their common faces.This step is relatively straightforward as we have computed modecurves in the faces first. Therefore, mode surfaces from neighboringtets already have matching polylines with the same vertices.
We use the same pipeline for both mode surfaces and neutral surfaces,with the following two differences.First, neutral surfaces and curves are characterized by degree-threepolynomials. Therefore, when finding the intersection with an edge,we can directly compute them using the cubic formula.Second, since there is no complex domain in the medium eigenvector manifold for neutral surfaces, we do not need to detect the intersectionof mode curves with the complex domain boundary. Instead, we needto find the singularities in the medium eigenvector manifold, whichcorrespond to straight lines in the
XY Z space. Their intersection withthe boundary faces can be easily computed. Then the intersectionpoints corresponding to the same singularity are connected using a linesegment.Other than these differences, our framework handles mode surfaces andneutral surfaces in a unified fashion.
ERFORMANCE
Our methods can extract mode surfaces and neutral surfaces with higherquality, i.e. the extracted surfaces are seamless and more accurate thanexisting mode surface and neutral surface extraction methods [27, 32](Figure 3). The seamlessness of the extracted surfaces enables addi-tional information about mode surfaces to be computed more robustly,such as principal curvature directions (Figure 2).In addition, our technique is faster than existing techniques. On av-erage, our neutral surface extraction method is about 1 . . .
40 GHz, 64GB of RAM, and an NVIDIA QuadroP620 GPU. We used four test data sets: a block with a single compres-sion force (Figure 1: 480000 tets), a block with a compression forceand an extension force (Figure 9: 450000 tets), a block with three com-pression forces (Figure 10: 384000 tets), and a block with a twistingcompression force (Figure 13 in Appendix B: 286416 tets).
PPLICATIONS
The application of compressive loads is ubiquitous in engineering wheresolids are extruded to designed geometries and liquids are compressedin combustion engines, and in medical research where human or an-imal organs are compressed to enable successful imaging processes.On the other hand, compression is a challenging state to model andelucidate. In fluid mechanics, Navier-Stokes modeling [5] applies toincompressible fluids and needs a completely different solution schemewhen the fluids become compressible. In solid mechanics, a scalarquantity referred to as the Poissons ratio [11] is used to dial frompolymers, which hardly compress, to materials that change volumeduring compression. Here we examine four scenarios of compressionin a solid block by visualizing the mode surfaces of the stress tensorfields. We vary boundary conditions on this block to reveal the materialbehavior embedded in the stress tensors. Three scenarios (Figure 8) arediscussed in this section while the fourth is covered in Appendix B.In the first scenario (Figure 8 (a)), there is a compression force onthe top of a cubic block. The resulting mode surfaces are displayedin Figure 1. Notice that the volume is dominated by negative modesurfaces (gold: indicating compression-dominant). In addition, thenetwork of degenerate curves leads to an interesting vascular structurethat is not commonly observed. This indicates that despite a simpleboundary condition, the shape of the block can also play an important (a) (b) (c)
Fig. 8: Three scenarios of compression in a solid block.8
Fig. 9: Mode surfaces of stress tensor fields for a rectangle block beingpushed and pulled simultaneously. The left side is compressed whilethe right side is extended through boundary conditions applied at thetop.role in deciding where uniaxial compression can occur.We contrast compression against extension in our second scenario(Figure 8 (b)), in which one side of the block is pushed down andthe other side is pulled up forming a dome-shaped dent and a dome,respectively. Our visualization (Figure 9) confirms this, as all planardegenerate curves (yellow: uniaxial compression) appear in the lefthalf and all linear degenerate curves (green: uniaxial extension) appearin the right half of the block. In addition, we show the generalizedmode 0 .
97 surface. Notice the symmetry between the mode surfaces forcompression (left side) and extension (right side). While compressionand extension are usually conceptualized as two different types ofdeformations (respectively volume loss and gain), to our knowledgethe symmetry between them in the generalized mode surfaces is newto the research and application communities. In addition, Figure 2(b) shows the major principal curvature directions for the generalizedmode 0 .
99 surface. Notice that in the negative mode part of this surface(gold), the major principal curvature direction is mostly aligned withthe width of the block except in the middle where it is aligned with thelength of the block. The sudden change in the directions is a reflectionof the fact that due to the boundary condition, the middle part of thissurface is being pushed out more along the block than across the block.The collision of compression directions leads to some trisectors on thesurface near the middle. Due to symmetry, similar deformations of thesurface (this time extension) can be observed on the positive mode partof the surface (teal). The number and location of the umbilical points(e.g. trisectors) can indicate uniform compression or extension andwill be investigated further. Note that such observations are impossiblewithout the ability to generate seamless meshes for mode surface (e.g.compare this to Figure 2 (a)).Our third example studies spatially consecutive compressive loads(Figure 8 (c)) which can be seen on highway bridges where fleets ofheavy trucks are parked due to traffic lights or in structures wherepairs of nuts and bolts are installed at a number of evenly spacedlocations. Intuitively, we can see that these compressive “holding-down” points create a wave form in the underlying geometry. In ourexample where we have three consecutive compressive loads, we notethat through analyzing the mode surfaces of the stress tensor field(Figure 10 (a)), this wave shape can be observed at mode value − . ONCLUSION AND F UTURE W ORK
In this paper, we provide to our knowledge the first approach to theseamless extraction of mode surfaces (including neutral surfaces). The (a) µ = ± . µ = ± . Fig. 10: Mode surfaces of stress tensor fields of a rectangular blockbeing pushed down at three locations.method is faster and of higher quality than existing methods of mode(neutral) surface extraction. Together with the degenerate curve extrac-tion algorithm of Roy et. [32], our approach represents the first unifiedapproach for the seamless extraction of mode surfaces of any mode.At the core of the approach is our novel topological analysis of modesurfaces (including degenerate curves and neutral surfaces) using thesame medium eigenvector manifold.In addition, we apply our technique to a number of data sets to fo-cus on analyzing stress tensor fields for compressive behavior that isfundamental in solid mechanics. We present a new demonstration ofcompression in volumes to indicate different behavior from differentregions of the material.Our current sampling strategy for the medium eigenvector manifolddoes not guarantee best quality of the triangle mesh such as the aspectratios of the triangles. We plan to explore other sampling patterns toimprove on this.There are other types of feature surfaces in 3D symmetric tensor fields,such as extremal surfaces [40] as well as magnitude surfaces andanisotropy index surfaces [27]. Adapting our parameterization ap-proach to those surfaces is also a future research avenue.Extending tensor field analysis and visualization to more complex phys-ical behaviors that are hybrids of different kinds of deformation can alsobe fruitful future research directions. Furthermore, correlating patternsobserved in the geometry of mode surfaces to the load configurationshas the potential of increasing our understanding and thus control inapplications such as tire performance and bridge maintenance. We planto explore these avenues in our future research. A CKNOWLEDGMENTS
We wish to thank our anonymous reviewers for their valuable sugges-tions. We thank Kyle Hiebel for making the voice recording of ourvideo. This research is partially supported by NSF awards ( R EFERENCES [1] G. E. Alefeld, F. A. Potra, and Y. Shi. Algorithm 748: Enclosing zerosof continuous functions.
ACM Trans. Math. Softw. , 21(3):327–344, Sept.1995. doi: 10.1145/210089.210111[2] M. Armstrong.
Basic topology . McGraw-Hill Book Co., 1979.
3] J.-D. Boissonnat and S. Oudot. Provably good sampling and meshing ofsurfaces.
Graphical Models , 67(5):405 – 451, 2005. Solid Modeling andApplications. doi: 10.1016/j.gmod.2005.01.004[4] L. Cammoun, C. A. Castano-Moraga, E. Munoz-Moreno, D. Sosa-Cabrera,B. Acar, M. Rodriguez-Florido, A. Brun, H. Knutsson, J. Thiran, S. Aja-Fernandez, R. de Luis Garcia, D. Tao, and X. Li.
Tensors in Image Pro-cessing and Computer Vision . Advances in Pattern Recognition. SpringerLondon, London, 2009.[5] A. J. Chorin, J. E. Marsden, and J. E. Marsden.
A mathematical introduc-tion to fluid mechanics , vol. 3. Springer, 1990.[6] J. C. Criscione, J. D. Humphrey, A. S. Douglas, and W. C. Hunter. Aninvariant basis for natural strain which yields orthogonal stress responseterms in isotropic hyperelasticity.
Journal of the Mechanics and Physics ofSolids , 48(12):2445 – 2465, 2000. doi: 10.1016/S0022-5096(00)00023-5[7] J. Damon. Generic structure of two-dimensional images under gaussianblurring.
SIAM Journal on Applied Mathematics , 59(1):97–138, 1998.[8] B. De-Arajo, D. Lopes, P. Jepp, J. Jorge, and B. Wyvill. A survey onimplicit surface polygonization.
ACM Computing Surveys , 47:1–39, 052015. doi: 10.1145/2732197[9] T. Delmarcelle and L. Hesselink. Visualizing second-order tensor fieldswith hyperstream lines.
IEEE Computer Graphics and Applications ,13(4):25–33, July 1993.[10] W. Fulton.
Algebraic Curves . Mathematics Lecture Note Series. W.A.Benjamin, 1974.[11] G. N. Greaves, A. Greer, R. S. Lakes, and T. Rouxel. Poisson’s ratio andmodern materials.
Nature materials , 10(11):823–837, 2011.[12] J. Greene. ”traces of matrix products.
Electronic Journal of LinearAlgebra. , 27, 2018. doi: 10.13001/1081-3810.1999[13] R. Grosso. An asymptotic decider for robust and topologically correcttriangulation of isosurfaces: Topologically correct isosurfaces. In
Pro-ceedings of the Computer Graphics International Conference , CGI 17.Association for Computing Machinery, New York, NY, USA, 2017. doi:10.1145/3095140.3095179[14] G. Guennebaud, B. Jacob, et al. Eigen v3. http://eigen.tuxfamily.org,2010.[15] L. Hesselink, Y. Levy, and Y. Lavin. The topology of symmetric, second-order 3D tensor fields.
IEEE Transactions on Visualization and ComputerGraphics , 3(1):1–11, Mar. 1997.[16] R. S. Irving.
Integers, polynomials, and rings: a course in algebra .Springer Science & Business Media, 2003.[17] C. Johnson and C. Hansen.
Visualization Handbook . Academic Press,Inc., USA, 2004.[18] F. Khan, L. Roy, E. Zhang, B. Qu, S. Hung, H. Yeh, R. S. Laramee, andY. Zhang. Multi-scale topological analysis of asymmetric tensor fields onsurfaces.
IEEE Transactions on Visualization and Computer Graphics ,26(1):270–279, 2020.[19] A. Kratz, C. Auer, M. Stommel, and I. Hotz. Visualization and analysis ofsecond-order tensors: Moving beyond the symmetric positive-definite case.
Computer Graphics Forum , 32(1):49–74, 2013. doi: 10.1111/j.1467-8659.2012.03231.x[20] W. E. Lorensen and H. E. Cline. Marching cubes: A high resolution 3dsurface construction algorithm. In
Proceedings of the 14th Annual Confer-ence on Computer Graphics and Interactive Techniques , SIGGRAPH 87,p. 163169. Association for Computing Machinery, New York, NY, USA,1987. doi: 10.1145/37401.37422[21] C. Luk and S. Mann. Tessellating algebraic curves and surfaces usinga-patches. In
GRAPP , pp. 82–89, 2009.[22] Y. Matsumoto.
An introduction to Morse theory , vol. 208. AmericanMathematical Soc., 2002.[23] R. Nickalls. Viete, descartes and the cubic equation.
The MathematicalGazette , 90(518):203–208, 2006.[24] G. M. Nielson. Dual marching cubes. In
Proceedings of the Conferenceon Visualization 04 , VIS 04, p. 489496. IEEE Computer Society, USA,2004. doi: 10.1109/VISUAL.2004.28[25] A. Paiva, T. Lewiner, Luiz, and H. D. Figueiredo. Robust adaptive meshesfor implicit surfaces. In
Computer Graphics and Image Processing, Brazil-ian Symposium on 0 , pp. 205–212, 2006.[26] J. Palacios, L. Roy, P. Kumar, C. Hsu, W. Chen, C. Ma, L. Wei,and E. Zhang. Tensor field design in volumes.
ACM Trans. Graph. ,36(6):188:1–188:15, 2017. doi: 10.1145/3130800.3130844[27] J. Palacios, H. Yeh, W. Wang, Y. Zhang, R. S. Laramee, R. Sharma,T. Schultz, and E. Zhang. Feature surfaces in symmetric tensor fieldsbased on eigenvalue manifold.
IEEE Transactions on Visualization and Computer Graphics , 22(3):1248–1260, Mar. 2016. doi: 10.1109/TVCG.2015.2484343[28] J. Palacios and E. Zhang. Interactive visualization of rotational symmetryfields on surfaces.
IEEE Trans. Vis. Comput. Graph. , 17(7):947–955, 2011.doi: 10.1109/TVCG.2010.121[29] D. Palke, Z. Lin, G. Chen, H. Yeh, P. Vincent, R. Laramee, and E. Zhang.Asymmetric tensor field visualization for surfaces.
IEEE Transactions onVisualization and Computer Graphics , 17(12):1979–1988, Dec. 2011. doi:10.1109/TVCG.2011.170[30] F. Raith, C. Blecha, T. Nagel, F. Parisio, O. Kolditz, F. G¨unther, M. Stom-mel, and G. Scheuermann. Tensor field visualization using fiber surfacesof invariant space.
IEEE transactions on visualization and computergraphics , 25(1):1122–1131, 2019.[31] X. Renbo, L. Weijun, and Y. Wang. A robust and topological correctmarching cube algorithm without look-up table. pp. 565 – 569, 10 2005.doi: 10.1109/CIT.2005.44[32] L. Roy, P. Kumar, Y. Zhang, and E. Zhang. Robust and fast extraction of3d symmetric tensor field topology.
IEEE Trans. Vis. Comput. Graph. ,25(1):1102–1111, 2019. doi: 10.1109/TVCG.2018.2864768[33] L. E. Spence, A. J. Insel, and S. H. Friedberg.
Elementary linear algebra .Prentice Hall, 2000.[34] B. T. Stander and J. C. Hart. Guaranteeing the topology of an implicitsurface polygonization for interactive modeling. In
Proceedings of the 24thAnnual Conference on Computer Graphics and Interactive Techniques ,SIGGRAPH 97, p. 279286. ACM Press/Addison-Wesley Publishing Co.,USA, 1997. doi: 10.1145/258734.258868[35] X. Tricoche, G. Kindlmann, and C.-F. Westin. Invariant crease lines fortopological and structural analysis of tensor fields.
IEEE Transactions onVisualization and Computer Graphics , 14(6):1627–1634, 2008. doi: 10.1109/TVCG.2008.148[36] E. Zhang, H. Yeh, Z. Lin, and R. S. Laramee. Asymmetric tensor analysisfor flow visualization.
IEEE Transactions on Visualization and ComputerGraphics , 15(1):106–122, 2009.[37] X. Zheng and A. Pang. Topological lines in 3d tensor fields. In
ProceedingsIEEE Visualization 2004 , VIS ’04, pp. 313–320. IEEE Computer Society,Washington, DC, USA, 2004. doi: 10.1109/VISUAL.2004.105[38] X. Zheng and A. Pang. 2D asymmetric tensor analysis.
IEEE Proceedingson Visualization , pp. 3–10, Oct 2005.[39] X. Zheng, B. N. Parlett, and A. Pang. Topological lines in 3d tensor fieldsand discriminant hessian factorization.
IEEE Transactions on Visualizationand Computer Graphics , 11(4):395–407, July 2005.[40] V. Zobel and G. Scheuermann. Extremal curves and surfaces in symmetrictensor fields.
The Visual Computer , Oct 2017. doi: 10.1007/s00371-017-1450-1 A T
HEORETICAL R ESULTS AND P ROOFS
Theorem 1.
Given a 3D linear tensor field T ( x , y , z ) = T + xT x + yT y + zT z , a mode value < µ < , and a unit vector v , the num-ber of points on the generalized mode µ surface with ± v as itsmedium eigenvector is the same as the number of real eigenvaluesof the 2D asymmetric tensor A = R θ + π T (cid:48) R θ − π where T (cid:48) is the pro-jection of characteristic tensor T onto the plane P with normal v , θ = arcsin ( √ ( arcsin ( µ ))) , and R φ = (cid:18) cos φ − sin φ sin φ cos φ (cid:19) . Thereal-valued eigenvectors of A give rise to the dominant eigenvectors ofthe corresponding points in the generalized mode µ surface.Proof. As pointed out in [32], given a tensor t = λ v v T + λ v v T + λ v v T in the linear tensor field T ( x , y , z ) = T + xT x + yT y + zT z , where v , v , and v are respectively the major, medium, and minor eigenvec-tors, we have v T T v + v T T v + v T T v = . (15)In addition, since T is the characteristic tensor of T ( x , y , z ) ,0 = (cid:68) T , λ v v T + λ v v T + λ v v T (cid:69) = λ trace ( T v v T ) + λ trace ( T v v T ) + λ trace ( T v v T )= λ v T T v + λ v T T v + λ v T T v . (16)according to the cyclic property of trace [12].Combining this equation with Equation 15, we have0 = ( λ − λ ) v T T v − ( λ − λ ) v T T v . (17)Notice that the major eigenvector v and minor eigenvector v must beinside P , the plane that contains the point where t occurs in the fieldand whose normal is v . Let v (cid:48) and v (cid:48) be v and v expressed in thecoordinate system of P . Consequently, Equation 17 can be rewritten asthe following:0 = ( λ − λ ) v (cid:48) T T (cid:48) v (cid:48) − ( λ − λ ) v (cid:48) T T (cid:48) v (cid:48) . (18)We first consider the case when µ ( t ) > v is the horizontal axis and v is the verticalaxis.For simplification purposes, we define u = kv (cid:48) + lv (cid:48) and w = kv (cid:48) − lv (cid:48) where k = (cid:113) λ − λ λ − λ and l = (cid:113) λ − λ λ − λ . Therefore, v (cid:48) = u + w k and v (cid:48) = u − w l .It is straightforward to verify that u and w both have unit length. More-over, since both v (cid:48) and v (cid:48) are unit vectors and have the same length, itcan be verified that u · w = k − l = λ + λ − λ λ − λ . (19)Since λ + λ + λ =
0, we have λ + λ − λ = − λ = √ (
13 arcsin ( µ )) . (20)Similarly, it can be shown that λ − λ = (cid:113) − λ = √ (
13 arcsin ( µ )) . (21)Consequently, u · w = λ + λ − λ λ − λ = √ (
13 arcsin ( µ )) = sin θ . (22)Notice that since k ≥ l ≥ w is always counterclockwise of u .Therefore, w is obtained by rotating u counterclockwise by π − θ . Notethat Equation 18 can be rewritten in terms of u and w as w T T (cid:48) u =
0, orequivalently w · ( T (cid:48) u ) =
0. That is, w ⊥ T (cid:48) u .As v (cid:48) is a bisector of u and w , v (cid:48) can be obtained by either rotating u counterclockwise by π − θ or rotating w counterclockwise by θ − π .That is, u = R θ − π v (cid:48) and w = R π − θ v (cid:48) . Therefore, R π − θ v (cid:48) ⊥ T (cid:48) R θ − π v (cid:48) . (23)Since v (cid:48) ⊥ v (cid:48) , we have v (cid:48) = R − π v (cid:48) . Consequently, R π − θ R − π v (cid:48) ⊥ T (cid:48) R θ − π v (cid:48) , which is equivalent to R − π − θ v (cid:48) ⊥ T (cid:48) R θ − π v (cid:48) . Rotatingboth sides conterclockwise by π + θ , we have v (cid:48) ⊥ R π + θ T (cid:48) R θ − π v (cid:48) . (24)Recall that A = R θ + π T (cid:48) R θ − π . Then v (cid:48) ⊥ Av (cid:48) .This means that Av (cid:48) must be orthogonal to a vector that is orthogonal to v (cid:48) . As v (cid:48) is in two dimensions, this implies that Av (cid:48) is a scalar multipleof v (cid:48) . Consequently, v (cid:48) is an eigenvector of A .In the case where µ ( t ) <
0, we use the right-handed coordinate systemof P such that v is now the horizontal axis and v is the vertical axis.We again define u = kv (cid:48) + lv (cid:48) but negate w , i.e. w = − kv (cid:48) + lv (cid:48) . We canstill show that u and w are unit vectors and u · w = √ ( arcsin ( µ )) .From here, we follow the same argument for the case where µ ( t ) > v (cid:48) and v (cid:48) . This leads to that v (cid:48) is an eigenvector of A = R θ + π T (cid:48) R θ − π .We now consider a number of cases. First, when A has zero real eigen-values, its eigenvectors must be complex-valued. Since the eigenvectorsof a 3D symmetric tensor cannot be complex-valued, there are no pointson the generalized mode µ surface with v as its medium eigenvectors.When A has only one real eigenvalue, there is only one point in the gen-eralized mode µ surface whose medium eigenvector is v . Dependingon the sign of the mode value of this point, either its major eigenvector v ( µ ( t ) >
0) or minor eigenvector v ( µ ( t ) <
0) is given by the cor-responding eigenvector of A . Recall that the dominant eigenvector isthe major eigenvector v when µ ( t ) = µ > v when µ ( t ) = − µ <
0. Therefore, the eigenvector of A gives thedominant eigenvector of the point in the generalized mode µ surface.11o see the uniqueness, we note that if there is another point in thegeneralized mode µ surface with the same combination of v , v , v and µ ( t ) , then the tensor must be a multiple of t . However, as pointedout by Roy et al. [32], if a tensor t appears in a 3D linear tensor field,then none of its multiples can appear in the same field. Consequently,there is only one point in the generalized mode µ surface with v as itsmedium eigenvector.On the other hand, when A has two real eigenvalues, the major eigen-vector and the minor eigenvector of A each corresponds to a point inthe generalized mode µ surface whose medium eigenvectors are givenby v . If the point has a positive mode value, its major eigenvector v is given by the corresponding eigenvector of A . In contrast, if thepoint has a negative mode value, its minor eigenvector v is given bythe corresponding eigenvector of A .What remains to be shown is that using − v as the medium eigenvector,we arrive at the same set of points as using v . Again, we have the twocases µ ( t ) > µ ( t ) < A arearound v , and they are reversed when − v is chosen as the mediumeigenvector. Assume that A ’s eigenvalues are ˆ λ and ˆ λ , and that v (cid:48) isthe eigenvector corresponding to ˆ λ . Then we find the eigenvalue onthe opposite side of the sphere with R − ( θ + π ) T (cid:48) R − ( θ − π ) v (cid:48) = A T v (cid:48) = R π ( R − π A T R π ) R − π v (cid:48) = R π ad j ( A ) R − π v (cid:48) = R π ad j ( A ) v (cid:48) = R π ˆ λ v (cid:48) = ˆ λ R π v (cid:48) = ˆ λ v (cid:48) (25)where we have used a property of asymmetric 2x2 tensors [33] that if A = (cid:20) a bc d (cid:21) is a 2D asymmetric tensor and R θ is the two-dimensionalcounterclockwise rotation matrix of angle θ , R − π A T R π = (cid:20) d − b − c a (cid:21) = adj ( A ) (26)is the adjugate of A , whose eigenvectors are the same as A ’s and the rolesof whose major and minor eigenvalues are swapped. Consequently, thepoint on the generalized mode µ surface with v given by the majoreigenvector of A ( v ) is the same point whose v is given by the minoreigenvector of A ( − v ) .A similar argument applies when µ ( t ) < v and − v correspond to the same set of points in thegeneralized mode µ surface. Moreover, when there are two pointswith v and − v as medium eigenvectors, each of A ( v ) and A ( − v ) corresponds to exactly one point in the generalized mode µ surfacewhose dominant eigenvector (as a 3D symmetric tensor) is given by themajor eigenvector of A (as an asymmetric tensor). Theorem 2.
Given a 3D linear tensor field T ( x , y , z ) = T + xT x + yT y + zT z and a mode value µ , the real domain in the medium eigen-vector manifold corresponding to µ is characterized by − v T T (cid:48) v + cos θ ( v T T (cid:48) v ) ≥ where θ = arcsin ( √ ( arcsin ( µ ))) . Theboundary between the real and complex domain occurs if and only ifthe equal sign holds. Proof. Let T (cid:48) be the projection of the characteristic tensor T of thetensor field onto the plane perpendicular to unit medium eigenvectors v .We consider the asymmetric tensor field A ( v ) = R θ + π T (cid:48) ( v ) R θ − π .As mentioned in [36], a 2 × A can be uniquelydecomposed as follows: A = γ d (cid:20) (cid:21) + γ r (cid:20) −
11 0 (cid:21) + γ s (cid:20) cos τ sin τ sin τ − cos τ (cid:21) (27)where γ d , γ r , and γ r are the isotropic, rotational, and anisotropic compo-nents, respectively. Note that τ gives rise to the eigenvector informationof A . Depending on the discriminant of ∆ = γ s − γ r , A has either tworeal eigenvalues ( ∆ >
0) or two complex-valued eigenvalues ( ∆ < ∆ = A has a pair of repeating eigenvalues. For our asymmetrictensor field, this implies that v is on the complex domain boundary inthe medium eigenvector manifold. Notice that this condition is bothnecessary and sufficient. To compute ∆ , we need to compute both γ d and γ r . γ d =
12 trace ( A )=
12 trace ( R θ + π T (cid:48) R θ − π )=
12 trace ( R θ − π R θ + π T (cid:48) )=
12 trace ( R θ T (cid:48) )=
12 trace ( cos θ I T (cid:48) ) +
12 trace (cid:18) sin θ (cid:20) −
11 0 (cid:21) T (cid:48) (cid:19) = cos θ trace ( T (cid:48) ) T (cid:48) is symmetric.Since trace is the same in any basis, we havetrace ( T (cid:48) ) = v (cid:48) T T (cid:48) v (cid:48) + v (cid:48) T T (cid:48) v (cid:48) . Recall that T (cid:48) , v (cid:48) , and v (cid:48) are respectively the projection of T , v , and v onto the plane perpendicular to v . Thus, v (cid:48) T T (cid:48) v (cid:48) + v (cid:48) T T (cid:48) v (cid:48) = v T T v + v T T v = − v T T v (29)in which the second equality above comes from rearranging Equa-tion 15.Consequently, γ d = − cos θ v T T v . (30)A similar calculation shows that γ r = sin θ trace ( T (cid:48) ) = − sin θ v T T v . (31)Notice that γ d + γ r = ( v T Tv ) . To compute γ s it is convenient to findthe squared magnitude of A , which is12 trace ( A T A ) = trace (( R θ + π T (cid:48) R θ − π ) T R θ + π T (cid:48) R θ − π ) = trace ( T (cid:48) ) . (32)Since v has unit length, the matrix I − v v T is a projection matrix andis thus equal to its square. Using this to project the tensor T , we havetrace ( T (cid:48) ) = trace (( I − v v T ) T ( I − v v T ) T ( I − v v T ))= trace ( T ( I − v v T ) T ( I − v v T ) )= trace ( T ( I − v v T ) T ( I − v v T ))= trace ( T ) − trace ( T v v T T ) − trace ( T v v T ) + trace ( T v v T T v v T )= − v T T v + ( v T T v ) . (33)Because A ’s magnitude is 2 ( γ r + γ s + γ d ) , we can use this to find γ s . γ s =
12 trace ( A T A ) − γ d − γ r = − v T T v + ( v T T v ) − (cid:18) v T T v (cid:19) = − v T T v + ( v T T v ) . (34)The real domain is characterized by γ s ≥ γ r , i.e.12 − v T T v + ( v T T v ) ≥ ( − sin θ v T T v ) (35)or equivalently, − v T T v + cos θ ( v T T v ) ≥
0. The complexdomain boundary is thus12 − v T T v +
14 cos θ ( v T T v ) = Corollary 3.
Given a 3D linear tensor field T ( x , y , z ) = T + xT x + yT y + zT z and two mode values µ and µ such that ≤ µ < µ ≤ ,the complex domain on the medium eigenvector manifold correspondingto µ is a proper subset of that corresponding to µ .Proof. The complex domain is characterized by the following inequal-ity opposite that for the real domain, i.e.,12 − v T T v +
14 cos θ ( v T T v ) < . (37)Recall that θ = arcsin ( √ ( arcsin ( µ ))) (Theorem 1), which im-plies θ is monotonically increasing with respect to mode values µ .Consider the range of θ which is [ , π ] , cos θ is a monotonicallydecreasing function with respect to µ . That is, for µ < µ , theircorresponding θ ’s satisfy cos θ > cos θ .Consequently, a unit vector v that satisfies12 − v T T v +
14 cos θ ( v T T v ) < − v T T v +
14 cos θ ( v T T v ) < . (39)That is, if v is in the complex domain for µ then it must also residein the complex domain µ . The reverse is not true. Consequently, thecomplex domain for µ is a proper subset of that of µ . Lemma 4.
Given a 3D linear tensor field T ( x , y , z ) = T + xT x + yT y + zT z and a mode value µ , the complex domain boundary correspondingto µ is parameterizable by α = v T v .Proof. Let λ ≥ λ ≥ λ be the eigenvalues of T and v , v , and v be the corresponding eigenvectors. For convenience, we choose thecoordinate system { v , v , v } in which T is a diagonal matrix. Let v = abc and m = m x m y m z = a b c . Equation 36 implies0 = − ( m x λ + m y λ + m z λ )+
12 cos θ ( m x λ + m y λ + m z λ ) . (40)This equation’s dependence on only m x = a , m y = b , and m z = c implies an eight-fold symmetry for the medium eigenvector manifold,and thus the complex domain boundary. That is, if a unit vector ( a , b , c ) is on the complex domain boundary, so is ( ± a , ± b , ± c ) .Therefore, we only need to parameterize the segment of the complexdomain boundary where a , b , c ≥
0. The other seven segments can beparameterized in a similar fashion.Notice that m is always on the plane m x + m y + m z = (cid:107) v (cid:107) =
1. In thecoordinate system α = λ m x + λ m y + λ m z = v T T v (41) β = λ m x + λ m y + λ m z = v T T v (42)Equation 40 becomes 0 = − β +
12 cos θ α (43)which is the equation of a parabola that can be parameterized by α .Each α gives one corresponding β from Equation 43 and thus one pointon the segment of the complex domain boundary where a , b , c ≥ Theorem 5.
Given a 3D linear tensor field T ( x , y , z ) = T + xT x + yT y + zT z , let µ is the mode of the characteristic tensor T and µ = (cid:113) − µ . The topology of the generalized mode µ surface isa topological torus when µ > µ and a topological double-torus when µ < µ .Proof. We reuse the expressions m x , m y , m z , α , and β from Lemma 4.Note that m x ≥ m y ≥ m z ≥
0, and m x + m y + m z =
1. Therefore,points satisfying these conditions form an equilateral triangle in theplane m x + m y + m z =
1, which is illustrated in Figure 11.Recall that T is a unit, traceless tensor with a non-positive determinant(Equation 3). We have13 a) Before bifurcation (b) Bifurcation (c) After bifurcation Fig. 11: In the triangle bounded by m x = m y =
0, and m z = m x + m y + m z =
1, the corners of the triangle are as follows: ( , , ) (red), ( , , ) (green), and ( , , ) (cyan). When µ decreases,the complex domain (grey regions) reduce in sizes. The complexdomain boundary initially intersects each of the m x = m y = m y = µ surface depends on whether there is a solutionwhen m x = λ + λ + λ = λ + λ + λ = λ ≥ λ ≥ ≥ λ (46)Plugging in λ = − λ − λ into Equation 45, we have λ + λ λ + λ = . This implies that λ + λ < . Consequently, λ < , λ < ,and λ > .For degenerate curves, i.e. µ =
1, cos θ = = − β , (47)which further reduces to 1 − λ i at the corners of the triangle. Con-sequently, this function is negative at ( , , ) (complex domain) andpositive at ( , , ) and ( , , ) (real domain). Because the function islinear, its zeroth levelset (complex domain boundary) must intersecteach of the edges m x = m y = ( , , ) is the complexdomain corresponding to µ = µ decreases, the complex domain (grey regions in Figure 11) re-duce in size (Lemma 4), and the parabola characterized by Equation 43moves lower while still intersecting m x = m y = c >
0. Since the original parabolic segment touches both the m x = m y = c <
0) also form a single loop. Moreover, this loop is to beidentified with the loop where c >
0. The real domain in the mediumeigenvector manifold is thus the part of the sphere between these twoloops. Gluing the two loops based on the antipodal symmetry results ina space without a boundary, the torus, which is homeomorphic to thegeneralized mode µ surface. The above situation changes when the parabola intersects the cornerof m x = m y = m y = µ = (cid:113) − µ , which isthe bifurcation point. When this happens, the four segments where c > µ surface in this case is a non-manifold surface.Decreasing µ further, the parabolic segment only intersects the m y = c >
0) in the medium eigenvector manifold form two loops(Figure 6: the innermost loop). Similarly, the four segments ( c <
0) alsoform two loops. The real domain is the part of the sphere between thefour loops. Gluing the four loops pairwise according to the antipodalsymmetry results in a sphere with two handles attached, i.e. double-torus. This is the topology of generalized mode µ surface when µ < µ .When µ =
0, it is straightforward to verify that the parabolic segmentdegenerates to a single point, and the aforementioned four loops shrinkto two pairs of antipodal points: (cid:18) ± (cid:114) λ − λ λ − λ ± (cid:114) λ − λ λ − λ (cid:19) (48)Each of the pairs corresponds to one of the two singularities in themedium eigenvector manifold for neutral surfaces [32]. Theorem 6.
Given a 3D linear tensor field T ( x , y , z ) = T + xT x + yT y + zT z and a plane P, the critical points of the mode function on theplane P consist of at most four extrema.Proof. On the plane P , the tensor field T ( x , y , z ) is still a linear tensorfield. Therefore, without the loss of generality, we assume the plane P to be the XY plane. If this is not the case, we can simply perform aspace transformation so that P becomes the XY plane under the newcoordinate system.In the XY plane, the tensor field has the form T ( x , y ) = T + xT x + yT y .We define a map χ ( x , y ) = T ( x , y ) (cid:107) T ( x , y ) (cid:107) from the plane P to the set ofunit tensors. Since the tensors on the plane is a combination of T , T x , and T y , the image of this map together with its negation forms atwo-dimensional sphere in the linear subspace spanned by T , T x , and T y . The map χ is injective because if a tensor has already appearedin a plane, its multiples cannot. Furthermore, χ is locally surjective.Therefore, the set of critical points of the mode function in the plane P has a one-to-one correspondence to the critical points of the modefunction on the aforementioned two-dimensional sphere.The mode of a unit tensor T is 3 √ ( T ) , which means that the criticalpoints of the mode function on the sphere can be found by computingthe critical points of the determinant function det ( T ) on the sphere.For convenience, we choose T , T , T to be a orthonormal basis for thespace spanned by T , T x , and T y .The critical points of a function defined on the sphere are the pointswhere its gradient is colinear to the sphere’s normal. Thus, the criticalpoints of the determinant function satisfy ∇ det ( uT + vT + wT ) × uvw = u + v + w = . (50)14 Let ∇ det ( uT + vT + wT ) = f ( u , v , w ) g ( u , v , w ) h ( u , v , w ) (51)where f ( u , v , w ) , g ( u , v , w ) , and h ( u , v , w ) are quadratic polynomials.Thus, Equation 49 consists of the following three cubic equations: vh ( u , v , w ) = wg ( u , v , w ) w f ( u , v , w ) = uh ( u , v , w ) ug ( u , v , w ) = v f ( u , v , w ) . (52)Together with the condition u + v + w =
1, we have an over-determined system of three cubic equations and one quadratic equation.Fortunately, the three cubic equations (Equation 52) are almost redun-dant. To see this, notice that multiplying the first two of these equationsgives rises to vw f ( u , v , w ) h ( u , v , w ) = uwg ( u , v , w ) h ( u , v , w ) (53)When w (cid:54) = h ( u , v , w ) (cid:54) =
0, the above equation reduces to the thirdcubic equation ug ( u , v , w ) = v f ( u , v , w ) . This shows the redundancy ofover-determined system.After removing the redundant third equation, we get the system vh ( u , v , w ) = wg ( u , v , w ) (54) w f ( u , v , w ) = uh ( u , v , w ) (55) u + v + w = , (56)which has up to 18 complex solutions based on B´ezout’s theo-rem [10]. However, there are some spurious solutions. Note thatfor vw f ( u , v , w ) h ( u , v , w ) = uwg ( u , v , w ) h ( u , v , w ) to be equivalent to ug ( u , v , w ) = v f ( u , v , w ) , we require that w (cid:54) = h ( u , v , w ) (cid:54) =
0. As-suming that w =
0, then we must have vh ( u , v , w ) = = uh ( u , v , w ) .Since u and v cannot be both zeros (or the vector ( u , v , w ) is the zerovector), we must have h ( u , v , w ) = h ( u , v , ) = u + v = − =
14 complex solutions, thus 14 real solutions atmost. Note that if ( u , v , w ) is a solution, so is ( − u , − v , − w ) . Moreover,they represent the same tensor. Consequently, there are only up to sevenreal solutions, i.e. up to seven critical points in the mode function inthe plane.The Euler Characteristic of a sphere is two [2], which means that thenumber of extrema is always two more than the number of saddlesaccording to Morse theory [22]. Together with the fact that there areat most 14 critical points, we can conclude that there are at most eightextrema on the sphere. Since the antipodal points on the sphere givethe same point on the plane, we have at most four extrema of the modefunction on the plane P . B O NE M ORE A PPLICATION S CENARIO
We provide one more scenario of compression on a block (Figure 12)in which the compressive load is misaligned with the natural geometricorientation of the block due to the addition of a slanted slab on thetop of the block. This slanted slab presses down on the block. At thebottom of the block, we also add a full-length layer with a more rigid (a)
Fig. 12: Another scenario of compression in a solid block. (a) Neutral Surface(b) µ = ± . Fig. 13: Mode surfaces of stress tensor fields of a rectangular blockbeing pushed down by a slanted slab.material. This set-up aims to provide a simplified representation of atire tread block anchored on some steel belts that hold the tire togetherwhile driving on hard pavement. Nonetheless, this type of boundarycondition is not limited to tire design. It can be due to misuse of thestructure. In particular, both cases exist on tires, i.e. some designs angletread blocks from the circumferential direction of the tire to facilitaterainwater drainage while irregular wear such as tire cupping inducesuneven tread block wear on the shoulders of the tire.The neutral surface of the stress is shown in Figure 13 (a), which,interestingly, reflects the existence of the slab. In addition, there aretwo sheets of neutral surfaces attached to the top face, two sheetsattached to the sidewalls, and two tubes that connect the front and backfaces. The location of these surfaces which indicate pure shear regionsmatch our expectation due to our design of the boundary conditions.In addition, other mode surfaces can also provide meaningful insight.For example, we show the generalized mode 0 ..