3D-TSV: The 3D Trajectory-based Stress Visualizer

We present the 3D Trajectory-based Stress Visualizer (3D-TSV), a visual analysis tool for the exploration of the principal stress directions in 3D solids under load. 3D-TSV provides a modular and generic implementation of key algorithms required for a trajectory-based visual analysis of principal stress directions, including the automatic seeding of space-filling stress lines, their extraction using numerical schemes, their mapping to an effective renderable representation, and rendering options to convey structures with special mechanical properties. In the design of 3D-TSV, several perceptual challenges have been addressed when simultaneously visualizing three mutually orthogonal stress directions via lines. We present a novel algorithm for generating a space-filling and evenly spaced set of mutually orthogonal lines. The algorithm further considers the locations of lines to obtain a more regular appearance, and enables the extraction of a level-of-detail representation with adjustable sparseness of the trajectories along a certain stress direction. To convey ambiguities in the orientation of the principal stress directions, the user can select a combined visualization of two principal directions via oriented ribbons. Additional depth cues improve the perception of the spatial relationships between trajectories. 3D-TSV is accessible to end users via a C++- and OpenGL-based rendering frontend that is seamlessly connected to a MatLab-based extraction backend. The code (BSD license) of 3D-TSV as well as scripts to make ANSYS and ABAQUS simulation results accessible to the 3D-TSV backend are publicly available.


INTRODUCTION
A S indicated by an informal questionnaire survey involving 9 experienced professors in computational mechanics, techniques for visualizing the three mutually orthogonal principal stress directions in 3D solids under loads are important in a number of use cases. In civil engineering such visualizations are used to develop and assess strategies for steel reinforcement of concrete support structures. In mechanical engineering, where often massive components like engines and pumps are considered, one is interested in how forces "find" their way through these components. The development of lightweight structures which are resistant to external loads is investigated in aerospace engineering, here stress directions provide the first indicators where structures can be hollowed. In topology optimization and bio-mechanics, such techniques are used to show tension and compression pathways simultaneously, and compare different structural designs regarding their mechanical properties [1].
An informative visualization of the stress directions in a 3D solid can be achieved via principal stress lines (PSLs), i.e., trajectories in 3D space along the principal stress directions. These trajectories are effective in communicating the pathways along which external loads are transmitted, and they show the mutual relationships between the different principal stress directions [2], [3]. In computational engineering, PSLs are used in particular to show where and how loads are internally redirected and deflected. There is agreement that such visualizations are necessary for a first qualitative analysis, before a quantitative analysis of certain regions using derived scalar stress measures is commonly performed.
However, our questionnaire has revealed a rather inconsistent use of stress trajectory visualizations in computational mechanics, and provides evidence that no standard tool for such an analysis exists. In 6 of the 9 groups we interviewed, own software packages for showing one particular principal stress direction starting at randomly selected locations are used. In 5 groups, CFD tools for flow visualization are used to show streamlines in a single principal stress direction field. Visualization tools that are able to show all principal stress directions simultaneously are rare, and also available post-processing tools do not offer this functionality.
One reason preventing a wider adoption of such tools is visual clutter and occlusions that are produced when showing the different types of PSLs simultaneously. Due to their mutual orthogonality, the visualizations appear irregular and unstructured, and perceptual coherence breaks up even for sparse sets of trajectories. While this effect can be reduced by starting trajectories from narrow regions and following only a single type of PSLs, this leaves large sub-domains uncovered and does not show the mutual variations of the stress directions. In general, clutter can be reduced by visualizing the single stress directions side-byside, yet juxtaposition makes it difficult to effectively relate the three mutual orthogonal stress directions to each other.

Contribution
This paper presents the 3D Trajectory-based Stress Visualizer (3D-TSV), a system and methodology for the visual analysis of the PSLs in 3D stress fields. Fig. 1 gives an overview of the visualization options provided by 3D-TSV. With 3D-TSV, we release, to the best of our knowledge, the first system supporting a comprehensive trajectory-based analysis of 3D stress fields. To achieve this, 3D-TSV builds upon existing techniques for line seeding in vector fields [4], [5], and it extends them towards the specific use case by considering simultaneously the three principal stress directions in the seeding process. To improve the regularity of the extracted arXiv:2112.09202v1 [cs.GR] 16 Dec 2021 PSLs, in the sequential seeding process every new seed point is located on an existing PSL belonging to a different principal stress direction. As proposed for streamlines in [4], [5], the seeding process is parameterized using different distance thresholds for each type of PSL, which allows controlling separately the sparseness of the PSLs of each type. We use this possibility to enable a Focus+Context (F+C) visualization that provides a fine-granular representation of a selected PSL type and simultaneously conveys the respective other PSLs as context information.
To ease integration into existing systems and accessibility to end users typically from engineering science, we have split the system into an OpenGL rendering frontend and a backend implemented in MATLAB. The frontend has been designed to reduce clutter and occlusions by keeping the number of rendered primitives low. This is achieved by converting one set of context PSLs into ribbons that follow either the major or minor principal stress direction and align with the medium principal stress direction. This approach enables to show-at least locally-also the third PSL type, and it conveys twists in the PSL-frame which cannot be seen when only two PSL-types are shown by lines.
The backend extracts trajectories from a given stress field using parameters that are either specified via the GUI that is built into the renderer or a configuration file. We have chosen a MATLAB backend especially due to the popularity of MATLAB in mechanical engineering. This enables engineers to easily integrate new model representations and algorithms. For instance, we provide MATLAB code for trilinear and inverse distance-based interpolation of stress tensors. If other types of basis functions are used, the corresponding MATLAB functions can easily be exchanged by specialized routines. Currently, the backend supports the extraction of stress trajectories in Cartesian and deformed hexahedral grids using tri-linear interpolation and inverse distance weighting, respectively. Due to the cell adjacency structure that is built internally to efficiently find the next cell during trajectory integration, other cell types can be supported with only minor additional effort.
To summarize, the contributions of this work are • a publicly available tool for trajectory-based stress tensor visualization supporting stress fields on arbitrary hexahedral grids, • the adaptation of evenly-spaced line seeding to create a space-filling set of PSLs with improved regularity, • a focus and context visualization using varying PSL density and visual mappings to lines and ribbons, • a GPU renderer that renders a ribbon as a flat elliptic tube to avoid flickering at oblique angles.
In a number of experiments using datasets with different shapes and stress states we demonstrate the application of 3D-TSV. To demonstrate the practical relevance of the proposed system, we focus on hexahedral simulation meshes, as they are in particular used in stress analysis and topology optimization. The code is made publicly available under a BSD license, and published on https://github.com/chrismile/LineVis (frontend) and https://github. com/Junpeng-Wang-TUM/3D-TSV (backend).

RELATED WORK
Stress Tensor Field Visualization. Stress tensor field visualization can be classified into trajectory-, glyph-and topology-based methods [1], [6]. Trajectory-based methods choose the PSLs as visual abstractions of the stress field, focusing on the directional structure of the principal stresses. Delmarcelle and Hesselink [7] introduced the concept of hyperstreamlines, a visual mapping of the medium and minor principal stresses onto a tube surface with a single selected major PSL as centerline. We propose a mapping of the medium principal stresses onto a ribbon with the major or minor PSL as centerline. Dick et al. [2] trace the major and minor PSLs from randomly distributed seed points in the loading area of the solid object, and different types of stress state like tension and compression are distinguished by color. In order to identify and visualize regions where stress trajectories are of rotational or hyperbolic behavior, Oster et al. [8] proposed the concept of tensor core lines in 3D second-order tensor fields. Hotz et al. [9] smear out dye along the PSLs using line integral convolution. In this way, a density field is generated that resembles a grid-like structure. This approach provides a global overview of a 2D stress distribution, yet an extension to 3D is problematic due to the generation of a dense volumetric field. Glyph-based methods, on the other hand, depict the stress field by a set of welldesigned geometric primitives -so-called tensor glyphs. Tensor glyphs were originally designed for glyph-based diffusion tensor visualization [10], and later adapted to visualize positive definite tensors [11], general symmetric tensors [12], as well as asymmetric tensors [13], [14]. Glyph-based techniques are problematic when used to visualize 3D stress fields, due to their inherent occlusion effects. Specific placement strategies can be used to reduce the number of glyphs and occlusions thereof [15], [16]. Tensor glyphs are effective in showing the local stress states, but they cannot effectively communicate the global structure of stress lines. Patel and Laidlaw [17] proposed to guide the placement of glyphs by principal trajectories in the underlying field, and thus to provide a better understanding of the global relationships in this field. Topology-based approaches for stress tensor visualization abstract is from the depiction of stress directions and focus on revealing specific topological characteristics of the tensor field. Delmarcelle and Hesselink [18], [19] studied the topology of symmetric 2D and 3D tensor fields, and introduced the fundamental concepts of degenerate points and topological skeletons. Zheng and Pang [20], and later Roy et al. [21], discussed the robust extraction of these topological features. Zobel and Scheuermann proposed the notion of extremal points to analyze the complete invariant part of the tensor [22]. Raith et al. presented a general approach for the generation of separating surfaces in the invariant space [23]. Palacios et al. introduced the eigenvalue manifold and visualized the 3D eigenvectors as curve surfaces [24]. Qu et al. [25] further generalized the concepts of degenerate curves and neutral surfaces to a unified framework called mode surfaces.
Streamline Seeding. Seeding strategies to control the density and placement of trajectories in vector fields are widely used in flow visualization. Turk and Banks [26] and Jobard and Lefer [4] were the first to introduce seeding strategies for generating evenlyspaced sets of streamlines in 2D vector field. Numerous extensions and improvements of these concepts have been proposed since then. In particular, Vilanova et al. [27] proposed an extension of the approach by Jobard and Lefer to diffusion tensor fields, which detects the distance between the new streamline and the existing ones during the tracing process. They demonstrate the generation of evenly distributed streamlines, however, the approach suffers from 'unfinished' streamlines that are caused by an artificial stopping criterion and only considers a single eigenvector field at a time. For 3D flow visualization, dedicated approaches have been developed to reduce the visual clutter and occlusion of densely distributed streamlines in 3D fields [28], [29], [30], [31]. However, these techniques do not fit our goal of visualizing PSLs and their mutual relationships, which requires considering three sets of orthogonal PSLs simultaneously.
Streamline Visualization. Illuminated streamlines are often used as a means of visualizing streamlines in a 3D environment. The streamlines are mapped to tubes (sometimes also called generalized cylinders [32]) with circular profile and then shaded, e.g., using the Blinn-Phong shading model [33]. Early work on illuminated streamlines was done by Zöckler et al. [34]. Mattausch et al. [5] further introduced the use of depth cues and halos at the screen space edges of streamlines. In their work, halos are created by rendering the streamlines a second time with black color, increased radius and adapted depth testing. Stoll et al. [35] extended this work by introducing stylized line primitives, rendered by a hybrid CPU-GPU renderer which performs line splatting on the GPU and line primitive tessellation on the CPU. We generalize the rendering of illuminated streamlines to tubes with an elliptic profile and perform all operations entirely on the GPU and in one single rendering pass. In order to identify the position of the halos, we compute the silhouette points of the elliptic tube cross-section, which generalizes the computation of silhouette points for circular tube cross-section by Blinn [36].
Hexahedral Meshing. An alternative approach to PSL-based stress field visualization is to generate a frame field from the principal stress field first and employ field-aligned hexahedral meshing to produce orthogonal edges that follow PSLs. The edges of such hex-meshes can follow the directions of PSLs excellently in situations where degenerate points are not present and the stress lines show low degrees of convergence and divergence. However, when guided with frame fields corresponding to realistic load situations, yet still much more benign than those demonstrated in this work, it is an unsolved problem to reliably produce an all-hex mesh. Hexahedral-dominant meshing has been resorted as an intermediate solution. For instance, Wu et al. [37] propose a conforming stress-guided lattice structure by combining topology optimization with the field-guided polyhedral meshing algorithm from [38]. Arora et al. [39] generate similar structural designs via the guidance of the principal stress field, where they modify the stress field to get a smooth frame field. However, hexahedraldominant meshes often contain either T-junctions or non-hexahedral elements with non-orthgonal edges, significantly deviating from the PSLs and are, thus, not applicable for stress field analysis either.

STRESS TENSOR DIRECTIONS
At each point in a 3D solid under load, the stress state is fully described by the stress vectors for three mutually orthogonal orientations. The second-order stress tensor contains these vectors for the axes of a Cartesian coordinate system. T is symmetric since the shear stresses given by the offdiagonal elements in T are equal on mutually orthogonal planes. The principal stress directions of the stress tensor indicate the three mutually orthogonal directions along which the shear stresses vanish. These directions are given by the eigenvectors of T , with magnitudes given by the corresponding eigenvalues. The signs of the principal stress magnitudes classify the stresses into tension (positive sign) or compression (negative sign). However, since there are three principal stresses acting at each point, the classification is with respect to a specific direction.
In descending order, the three eigenvalues of T represent the major σ 1 , medium σ 2 and minor σ 3 principal stresses, with the corresponding eigenvectors indicating the principal stress directions at each point in the 3D solid. The trajectories along these directions are called the principal stress lines (PSLs). They are computed by numerically integrating massless particles in each single (normalized) eigenvector field.
In general, σ 1 , σ 2 and σ 3 are mutually unequal, and the eigenvectors are linearly independent and even mutually orthogonal due to the symmetry of T . However, so-called degenerate points can exist where two or more eigenvalues are equal. In the vicinity of these points, which are classified by σ 1 = σ 2 > σ 3 or σ 1 > σ 2 = σ 3 1 , the PSL direction cannot be decided. Therefore, when tracing along a principal stress direction, we test whether the eigenvalue σ i corresponding to this direction is too close to another eigenvalue If this is the case and the angle between the current and next PSL tangent is too large and indicates an ambiguous situation, the integration is stopped. Furthermore, we provide the option to map deg to color via a color table (see Sec. 5), and color the PSLs so that the proximity to a degenerate point is indicated. The integration is also stopped when the next integration point is located on a boundary face, the point is closer to a previous point on the same trajectory than a predefined distance threshold, or the number of integration steps reaches the pre-defined threshold.
The integration of PSLs requires to select seed points from which they start until they arrive at a degenerate point or the boundary. In this work we employ a Runge-Kutta RK2 scheme with fixed integration step size δ for numerical integration. In each integration step, we interpolate the stress tensor T and compute the eigenvalues and eigenvectors from the interpolated tensor. If none of the mentioned stopping criteria holds, the next step is performed in the direction with the least deviation from the previous direction.

PSL SEEDING AND LEVEL OF DETAIL
Finding a set of PSLs that effectively convey the principal stress directions in 3D stress fields requires to consider perceptual issues related to the visualization of large sets of trajectories. While in principle the PSLs of a single type, i.e., major, medium, or minor, can be visualized separately using techniques from flow visualization, in a stress field the different types of PSLs need to be shown simultaneously to understand their mutual interplay. However, an effective and efficient visual analysis is hindered by the mutual orthogonality of the different types, which is perceived as a disordered state even when a low number of PSLs is shown. Our proposed seeding strategy cannot completely avoid this problem, but it has some built-in regularity due to enforced PSL intersections.

Evenly-spaced PSL Seeding
The proposed seeding strategy builds upon the evenly-spaced streamline seeding approach by Jobard and Lefer [4], and extends 1. We do not consider triple degenerate points with σ 1 = σ 2 = σ 3 , since they do not exist under structurally stable conditions [20]. this approach in several ways to account for the application to PSLs. For the sake of clarity, we describe the strategy in the context of 2D stress fields, yet it will become clear that the extension to 3D is straightforward. However, when applied in 3D, the resulting PSL structures show a fundamental difference. Unlike in 2D, where due to the intersections between major and minor PSLs a fairly regular grid-like structure is generated, such intersections are rare or do not exist at all when seeding PSLs in 3D. This counteracts the impression of a consistent grid-like structure and results in a rather disordered appearance. We propose a seeding strategy that weakens this effect, but it needs to be considered that due to the nature of PSLs in 3D stress fields a globally consistent 3D grid-like structure is impossible to achieve in general.
Our method builds upon the selection of new seed points in the spirit of Jobard and Lefer, where the potential candidates are those points which are at least a prescribed distance away from any already extracted PSL. Of these candidates, the one with minimum distance is selected and a new trajectory is started at that point. In contrast, in our approach the distance is always wrt. the initial seed point, so that the PSLs grow around that point instead of being seeded at vastly different locations.
To adapt the seeding strategy to the situation of different types of PSLs, we first introduce the concept of seed valence. In 2D, the seed valence ϑ is a 2 × 1 binary array, which is associated to each seed point to indicate whether and of which type PSLs have been traced from this point. ϑ can take on four different bit combinations, i.e., empty seed [0 0] (passed by no PSL), solid seed [1 1] (passed by both major and minor PSLs) and semi-empty seed [1 0] (only passed by major PSL) or [0 1] (only passed by minor PSL). The sampling process is repeated until all valences of all possible seed points become solid [1 1]. With this definition of seed valence, the sampling process is performed iteratively, by using the seed valence to characterize the state of each seed point at a specific iteration. To ensure that the generated PSLs are space-filling, the initial candidate seed points (with ϑ = [0 0]) are located at the vertices of a space-filling Cartesian grid (step 0 in Fig. 2).
Seeding starts by selecting one of the candidate seed points and tracing the major and minor PSLs from it (Step 1 in Fig. 2), setting ϑ = [1 1] at this point. Per default, the system starts with the seed point closest to the center of the bounding box of the domain, to preserve an existing plane symmetry of the stress field in the PSLs (see Fig. 11, Fig. 12 and Fig. 13). Then, all candidate  Fig. 10 for the simulated load conditions). Major (ocher), medium (green) and minor (blue) PSLs generated by (top) evenly-spaced streamline seeding [4] in each PSL field separately, and (bottom) by our method. Note that since the stress field is not strictly symmetric, the PSL set shows some asymmetry. seed points with ϑ not equal to [1 1] are re-classified with respect to the currently existing PSLs. To exclude candidates too close to an existing major or minor PSL, ϑ of these candidates is set to [1 0] or [0 1], respectively. If a point is classified as [1 0] or [0 1] and closer to a minor or major PSL, respectively, its valence is set to [1 1]. The distance between a point and a PSL is computed as the minimum distance between the point and any of the integration points on the PSL. Proximity is decided via a distance threshold ε, which also controls the density of the extracted PSLs.
To obtain a more regular PSL structure, each re-classified candidate point is re-located (i.e., snapped) to the position of the closest integration point on the PSL causing its classification. This creates an "empty" band around the PSLs where no candidate seed point exists. The snapping operation enforces that newly selected seed points lie on an existing PSL, so that the final PSL structure appears more regular and less cluttered (see Fig. 3 for a comparison to the seeding approach by Jobard and Lefer). By placing the initial seed point in a region deemed important, the user can specifically enforce regularity in this region.
If the last computed PSL was a major or a minor PSL, then the next seed point is selected from the set of candidates with ϑ = [1 0] or [0 1], respectively. Thus, we alternate the order of major and minor PSL extraction to obtain a uniform distribution of both types. Of all these, the one closest to the initial seed point is selected as the new seed point, and the respectively transverse PSL is computed. The entire procedure is then restarted until no more candidate is available (see steps 2-5 in Fig. 2).
We further consider the situation where some empty seed points may get too close (measured by ε) to the other type of existing PSLs after they are merged to the current PSL, e.g., the seed valence ϑ of some empty seed points become [1 0] after merging them to the newly traced major PSL. However, it can also happen that some of these merged seed points might be close to some of the existing minor PSLs, which would unavoidably cause inappropriate placement of minor PSLs in the final visualization. Given this, we identify those semi-empty seed points after snapping, and compute the distances of them to the corresponding type of PSLs. If there are distances less than ε, the valences of these seed points are set to [1 1]. By simply making ϑ a binary array with three elements referring to the major, medium and minor PSL, the proposed seeding strategy can be lifted to 3D.

PSL LoD Structure
To change the density of the generated PSLs, the seeding process can simply be re-run with an appropriately set distance threshold ε. The larger this threshold is, the less PSLs are extracted. However, the different sets of PSLs that are generated for different thresholds are not nested, i.e., the PSLs at a coarser representation with lower PSL density are not a subset of the PSLs at a representation with higher density. Therefore, in an exploration session where the user interactively selects different PSL LoDs, there are abrupt changes when transitioning from one level to another. To avoid this, we propose to generate a nested PSL hierarchy.
The basic idea underlying the construction of a nested hierarchy is to let the PSLs at a level with higher PSL density 'grow out' sequentially from the PSLs at a lower density level. As a side effect, this enables saving computations by progressively computing a new level from the previous coarser level. For a given set of PSLs that have been generated with distance threshold ε 0 , the refined set of PSLs according to a distance threshold ε 1 < ε 0 is computed as follows: Firstly, the candidate seed points are reset to their initial positions. Secondly, the candidate seed points are snapped to the existing PSLs according to ε 1 , to create "empty" bands around the existing PSLs. The valences are updated accordingly to [1 0], [0 1] or [1 1] depending on the types of PSL they are snapped to. After this, some (semi-)empty or empty seeds are left, because ε 0 is larger than ε 1 . With these seeds the seeding is subsequently performed, including the iteration of seed point selection, PSL computation, and re-classification as described in subsection 4.1.
To generate a full LoD PSL hierarchy, the user defines the minimum distance threshold ε and the number M of levels to construct. Then, the distance thresholds of each level are computed as 2 (M−k) ε, k = 1 : M from coarse to fine, and the hierarchy is constructed progressively from the coarsest resolution level (see 1st and 2nd row in Fig. 4). To compute a PSL structure with different types of PSLs at different LoDs, the distance thresholds for each PSL type are first selected by the user, and then the multi-type LoD is computed by alternatively considering the different PSL types with their respective distances.

RIBBON-BASED FOCUS+CONTEXT
Regardless the efforts one puts into the creation of more regular PSL structure, showing all three types of PSLs together can lead to a rather cluttered view. Clutter can be reduced by showing PSLs at coarser LoDs, yet at a sparseness level that allows for an uncluttered view, many important principal structures can be missed. To address this problem, the system enables users to focus their visualization efforts on one selected PSL type, while preserving a surrounding context that conveys important directional cues provided by the respective other types. While in principle this can be achieved by showing the focus and context PSLs at a fine and coarse LoD, respectively, when visualizing PSLs via lines the few context lines quickly get lost in the many focus lines. Our solution is to show a selected principal direction as focus via lines, and utilize a contextual ribbon-shaped geometry [40] that is aligned along another direction and twists according to the remaining one (see Fig. 5 (a)(b)). The mapping of two principal stress directions to a ribbon geometry is conceptually similar to the well-known hyperstreamlines [7], i.e., a mapping of two principal stress directions to a tube centered at the remaining directions. The constructed ribbon doesn't necessarily coincide with a streamsurface that is traced from a PSL along one other stress direction, and which might not even exist [23], yet it enables an effective qualitative analysis of the mutual variations of two stress directions. Compared to a hyperstreamline, a ribbon causes less occlusions and allows to access the remaining third stress direction in relation.
A beneficial property of a ribbon-shaped geometry is that it can effectively convey changes in the assignment of the eigenvector directions to the type of PSL in the vicinity of degenerate points. When a ribbon is traced along a principal stress direction and twists according to another direction, flips can occur in the vicinity of a degenerate point (see Fig. 5 (c)). This is because the two directions can exchange their classification as major, medium, and minor, since this depends only on their position in the sorted sequence of eigenvalues. Thus, ribbons provide a strong visual cue to indicate topological changes of the force pathways at degenerate points. Fig. 6 compares the options to visualize principal stress directions via ribbons and lines, and combine them into a F+C visualization. As can be seen, twists in the ribbon geometry effectively hint to regions where degenerate points might exist. For lines, 3D-TSV can map the degeneracy measure introduced in Sec. 3 to color. An interesting observation is that high degeneracy and flips thereof frequently occur close to the object boundaries when Cartesian simulation meshes are used. These flips occur due to the well-known inaccuracies at curved boundaries that are represented by hexahedral simulation elements in a Cartesian grid.

RIBBON-BASED PSL RENDERING
The rendering frontend enables the user to interactively explore the stress field. In order to make real-time exploration possible, all computations performed after data loading, which itself is performed in a separate background thread, are parallelized on the GPU by utilizing OpenGL.
GPU rendering of tubes and ribbons can be performed efficiently by using a geometry shader to map lines to circular tubes and flat quads. Circular tubes are constructed on-the-fly in a geometry shader by generating for each line vertex a set of vertices and forming triangles by connecting these sets of vertices across line segments. For the latter, each vertex is equipped with an extrusion direction so that a quad can be spanned on either side of the initial line. However, in this way ribbons become extremely thin and invisible when viewed at oblique angles. This turns out to be a major weakness due to flickering that is introduced by ribbons which suddenly disappear and popup again in a next view (cf. Fig. 7).
To avoid the aforementioned effect, we propose a rendering technique for ribbons that builds upon the mapping of lines onto elliptic tubes (cf. Fig. 8). While circular tube geometry has been used in previous work for stylized line rendering [35], [36], we are not aware of any work that addresses the generalization to stylized rendering of elliptic tubes and its use for ribbon rendering including halos (i.e., screen space outlines).
An elliptic tube of a curve x(τ) with radius r and dilation factor w in normal direction n(τ) is defined as the set of points p(τ, ϕ) = x(τ) + rF(τ) w · cos(ϕ) sin(ϕ) 0 T . (3) is a frame matrix based on the Frenet-Serret frame [41], [42] with t(τ) being the tangent, n(τ) being the normal, and b(τ) being the binormal of the curve x at τ. ϕ is an angle in the range [0, 2π). For discretized curves in 3D space, the tangent is approximated by the direction of the line segment between two adjacent points, and the normal and binormal are chosen as two vectors orthonormal to the tangent (cf. Stoll et al. [35]). The cross-section of a tube, which is shown in Fig. 9 left, is defined as the intersection of the tube with a normal plane of its underlying curve. The normal plane is a plane containing the normal and binormal of the curve at a certain point. Lines can be mapped to circular or elliptic tubes by adapting the dilation factor w. For w = 1, the cross section of the tube is a circle and lines are rendered as circular tubes. Ribbons are rendered as elliptic tubes (i.e., the cross section of the tube is an ellipse). By using a dilation factor w < 1 of the ellipse in the normal direction of the ribbon, the torsion of the band can be encoded due to the anisotropy of the ellipse. For lim w→0 , a perfectly flat ribbon geometry is obtained. This solves the problem that ribbons can become extremely thin, as w > 0 ensures that they have a minimum thickness when viewed from the side. However, a problem that arises when rendering elliptic tubes is that the computation of stylized screen space outlines is not as straight-forward as for circular tubes.
In order to compute the outline, we need to project the position of a fragment processed in a shader to an underlying screenoriented silhouette and compute the ratio of the distance of the projected point to the center of the silhouette and the total length of the silhouette. Then, black outlines can be rendered when the absolute silhouette position ratio is close to one. In the following paragraphs, we assume that all ellipses have unit extent in one direction, and extent w in the other direction. The extents will normally be modulated with the ribbon width. The simpler special case of circular cross sections has first been described by Blinn [36].
Without loss of generality, we further assume that the center of the elliptic cross-section lies at (0, 0) T and the coordinate axes are aligned to the tangent frame of the curve. The world-space camera position c world (assuming that the coordinate system is normalized to have its origin at the center of the ellipse) can be projected onto the cross-section normal plane by computing The quadratic form to express that a point in homogeneous coordinates p = (x, y, 1) T lies on the ellipse in Fig. 9 right is This quadratic form can also be expressed as p T Ap = 0 with An ellipse is just a special case of a conic section in two dimensions. Given the position c of the camera in the normal plane, Fig. 8. From left to right: Band rendering, circular tube rendering with outlines and depth cues, elliptic tube rendering with outlines and depth cues. Bands can become very thin, and circular tubes cannot convey twists. Elliptic tubes solve both issues.
we want to compute the two points p max and p max on the conic section which are touched by the tangents of the conic section wrt. the camera position. For this, we can compute the polar l = Ac of the camera point, which intersects the conic section in these two points [43]. In order to compute the intersection of l with the conic section, the algorithm described by Richter-Gebert [44] can be employed. Then, with τ being a non-zero entry in l, we can compute Let C = B + αM l and (i, j) be the index of a non-zero element C i, j of C. Then the two intersection points p max and p max are the ith row and jth column of C, respectively. Finally, we can compute p = meet(l, join(c, p)) = l × (c × p). The meet operation applied on two lines yields their intersection point, and the join operation applied on two points yields their connecting line (cf. [45]). Then, the absolute silhouette position is equal to pos = p −p max 2 p max −p max 2 · 2 − 1 .

SYSTEM IMPLEMENTATION
To implement the communication between the C++ visualization frontend and the MATLAB extraction backend, the messaging library ZeroMQ is utilized, which can be used for communication over a wide variety of protocols, like TCP/IP. 3D-TSV relies on the request-reply pattern implemented in ZeroMQ, where the frontend issues a new request to the backend when the user changes simulation settings in the graphical user interface, and the backend sends back a reply as soon as the simulation is finished in order to notify the frontend of the availability of new data. Fig. 10. The solid objects used in this work and the applied external loads. Red and blue arrows indicate the loading positions and directions, black regions indicate fixed boundaries. A finite-element-based elasticity analysis has been used to compute the stress field for each model under the predicted loads. The unstructured hexahedral meshes 'BunnyHex' and 'KittenHex' are courtesy of [46] and [47], respectively. All other meshes are Cartesian meshes. 'Arched Bridge', 'Parts' and 'Rod' are courtesy of [48], [46] and [47], respectively. All simulated stress fields are made publicly available.
The reason why we turned to MATLAB instead of C++ for the implementation of the backend is, on the one hand, that the sampling method is an inherently sequential algorithm. Thus, it cannot benefit significantly from multi-threaded PSL tracing or GPU parallelization. On the other hand, MATLAB is widely spread in engineering, where most of our collaborators regarding stress visualization come from, and the engineers tend to use mainstream commercial software they are already familiar with to finish the design iteration quickly. In this case, they can run the MATLAB backend independently without any complicated compilation and setup process. To this end, we also provide a slim MATLAB visualization implementation, which can provide users a fast and easy way to explore the stress field, while discarding some more complex hardware-accelerated features from the C++ frontend, like depth cues or elliptic tubes. It is worth noting that also the rendering frontend can be used standalone, by reading trajectories from a file specifying the exchange format regarding PSL type and LoD representation.

Numerical PSL Integration
3D-TSV is designed to support a visual analysis of the PSLs in solids discretized by hexahedral grids, where the stress tensors are given at the grid vertex. When computing PSLs in Cartesian grids, component-wise trilinear interpolation of the tensors is used during numerical integration. For arbitrarily deformed hexhedral grids, we use inverse distance weighting for element interpolation [49]. Note here that we assume each face of a hexahedral cell to be planar, as it is common in finite cell and element simulation methods. For celllocation, i.e., locating the cell which contains the next integration point, we perform a two-step strategy: Firstly, elements potentially containing the point are quickly located by means of a hierarchical spatial search data structure [50]. Secondly, the element containing the point is determined by taking into account the following inout criterion. Given a hexahedral element with the centers and outer normal of its six faces C i and V i , i ∈ {1, · · · , 6}, any point B satisfies max(arccos( BC i ,V i )) ≤ π 2 , i ∈ {1, · · · , 6} iff it is inside of this element or lies on its boundary surface. It is worth noting that sometimes the normal of each element face might not be unique, since the 4 vertices of each element face do not exactly lie on one single plane. To deal with this, on the one hand, we approximate the quadrilateral element face by its two diagonals. On the other hand, the in-out criterion is slightly relaxed to max(arccos( BC i ,V i )) ≤ 91π 180 , i ∈ {1, · · · , 6}. With a constant integration step, it can further happen that tiny elements are skipped when adjacent elements are small compared to the step size. In principle, this can be avoided by traversing from cell boundary to cell boundary along the current stress direction, and adapting the direction at each cell until the step size is reached. However, due to numerical issues regarding the accuracy of line-face intersections in the vicinity of face edges, we refrained from considering this strategy.

Rendering
The line and ribbon primitives are rendered in a stylized fashion visually similar to the techniques by Zöckler et al. [34], Stoll et al. [35] and Mattausch et al. [5], using default colors, halos and depth cues as shown in the first three images in Fig. 1. However, unlike Stoll et al. [35], the discretized geometry of the tubes is created solely in a geometry shader on the GPU and no CPU-GPU communication is necessary. In contrast to the halo rendering introduced by Mattausch et al. [5], the black outlines are computed in one rendering pass (cf. section 6) and duplicated rendering of the geometry is not necessary. Focus PSLs and contextual ribbons are rendered in ocher and blue, respectively. The base color is modulated using Blinn-Phong shading [33], [34], which assumes a point light source at the world space position of the viewer (i.e., a head light).
The user can interactively change the color mapping-also separately for each PSL type-and can in particular switch to a mapping of some scalar quantity to color, as indicated in the last image in Fig. 1 using the scalar von Mises stress measure. The scalar values are issued via the backend as per-vertex attributes. The standard color scheme we use for the different principal stress directions (blue, green, ocher) is the '3-class Set2' transfer function from ColorBrewer. It is colorblind safe and print friendly.
For enhanced depth perception, depth cues are added, i.e., with increasing distance to the camera, fragments are increasingly desaturated. Therefore, the program recomputes the depth range each time the camera moves or the displayed data is changed. This is done solely in compute shaders on the GPU by projecting all vertices into view space and then computing the minimum and maximum depth using parallel reduction operations. The reduction operations are highly efficient due to the implementation utilizing on-chip shared memory (cf. [51]), and they avoid costly synchronization and copy operations between the CPU and GPU.
Additionally, a translucent simulation mesh outline hull is rendered together with the stress field data in order to hint at the extents of the simulation domain. In order to avoid the desaturation of rendered line and band primitives when being overlaid by the mesh hull, the opacity of the hull is modulated by one minus the cosine of the angle between the surface normal and the viewing direction.

3D-TSV Settings
3D-TSV provides a number of parameters that can be changed by the user to control the generation of PSLs. These parameters include the merging threshold ε and the number of levels M introduced in subsection 4.1 and subsection 4.2, respectively. Another set of parameters enables a user-guided interaction with the PSL distribution, including sliders for controlling the LoD resolution of major, medium and minor PSLs. In addition, the user can select the type of PSL in focus and context, and the LoD combination. Via a drop-down menu, the user can select a scalar stress component field that is mapped to PSL color using a transfer function. The backend provides different stress components, such as the principal stress amplitudes, von Mises stress, and the six Cartesian stress components.
After discussing 3D-TSV with domain experts, some additional options have been realized to enable a more expert-like use of 3D-TSV: While uniform seeding in the entire domain is used as the default option, the user can select seeding from the boundary vertices as well as the loaded and fixed vertices. Furthermore, different integration schemes can be used for PSL tracing, including the 1st-order Euler method, and the 2nd-and 4th-order Runge-Kutta methods.

ADDITIONAL RESULTS
In all of our examples, PSL generation is performed on the CPU, i.e., a workstation running Ubuntu 20.04 with an AMD Ryzen 9 3900X @3.80GHz CPU and 32GB RAM. Rendering is done on an NVIDIA RTX 2070 SUPER GPU with 8GB of on-chip memory. The rendering times are always below 10 milliseconds. The stress fields in the solid objects under different load conditions in Fig. 10 are simulated by a finite element method (FEM). Table 1 lists the corresponding number of hex-elements of each of the test data sets, the seed points that are used to generate the PSLs, the number of generated PSLs, and the time required for PSL generation.
In this section, we first consider three representative examples to show the improvements of our proposed PSL seeding strategy over the existing evenly-spaced streamline seeding approach. Then, additional examples show the combination of the F+C and ribbon rendering as provided by 3D-TSV. Finally, we demonstrate the application of 3D-TSV for PSL extraction and visualization in unstructured hexahedral meshes.
First, we consider the symmetric stress field in 'Cantilever' (Fig. 10). As shown in Fig. 11left, 3D-TSV generates a spacefilling and well-structured visualization, and the symmetry of the stress field is well revealed. Fig. 11right shows the result of evenly spaced streamline seeding, demonstrating the improved regularity and reduced visual clutter by our proposed strategy.
The visualization also highlights the importance of showing different PSL types simultaneously, i.e., otherwise one physical  Fig. 12 and Fig. 13 further demonstrate this property of our seeding strategy. Especially, the closeup views reveal more PSL intersections and improved visual quality thereof. Fig. 14top shows the space-filling PSLs in the 'Bracket' stress field. By analyzing the corresponding boundary condition in Fig. 10, we see that the structure is mainly under tension, so we choose the major PSLs at L2 as focus and the minor PSLs at L1 as context in Fig. 14bottom. This enables a fine granular analysis of the major principal stress directions.
We further consider the 'Parts1' model, which has more geometric details, hence, the resolution of the PSL distribution should be so high that no important features are missed in the visualization, see Fig. 15left. The high resolution of PSLs transmits the overall impression of the stress field, but accompanied by serious occlusion problem. The LoD enables an effective way to explore this type of data set by interactively adjusting the density of PSLs, according to specific analysis requirements. In the middle of Fig. 15, the F+C setting major PSLs at L2 and minor ones at L1 is used to highlight the major principal stress field, and the opposite settings are used for Fig. 15right.
3D-TSV works with both Cartesian meshes and deformed hexahedral meshes. Here we use the stress fields from 'BunnyHex' and 'KittenHex' as test cases to demonstrate the capability of our method. Especially in the 'KittenHex' model, the element sizes change considerably as shown in Fig. 10. The full distribution of PSLs is shown in the first column of Fig. 16. The F+C visualizations are shown in the second and third columns, where, for the bunny, the combinations of major PSLs at L1 and minor PSLs at L3 (middle), major PSLs at L2 and minor PSLs at L1 (right) are used, for the kitten, we use major at L3 and minor at L1 (middle), major at L1 and minor at L2 (right), respectively.

CONCLUSION AND FUTURE WORK
In this paper, we have introduced 3D-TSV, a tool for visualizing the principal stress directions in 3D solids under load. 3D-TSV makes use of a novel seeding strategy, to generate a space-filling and evenly-spaced set of PSLs. By considering all three types of PSLs simultaneously in the construction process, the regularity of the resulting PSL structure is improved. By incorporating different merging thresholds for each PSL type into the construction process, a consistent multi-resolution hierarchy is formed. This possibility has been utilized in the development of a focus+context visualization, which simultaneously shows a PSL type in focus via lines and the remaining PSLs as context information using ribbons. Efficient rendering options for lines and ribbons on the GPU enable interactive analysis of large sets of PSLs.
In the future, we intend to couple 3D-TSV with load simulation processes, so that dynamic changes of the stress field can be instantly monitored. Therefore, we will analyze whether the intrinsically iterative parts of the algorithm can be parallelized on modern multi-threading architectures. Furthermore, we are interested in using space-filling evenly-spaced seeding to guide the material growth in topology optimization. Topology optimization seeks to distribute material in a way that makes the object resistant to external loads. To automatically generate support structures that follow the major stress directions and eventually can form a 3D grid-like structure, we aim at combining our seeding strategy with the automatic growth process underlying topology optimization.   Table 1. Bottom: F+C visualization with major PSLs at L3 in focus and minor PSLs at L1 as context.