1 Introduction
Techniques for visualizing the three mutually orthogonal principal stress directions in 3D solids under load are important in a number of use cases in computational mechanics. In civil engineering such visualizations are used to develop and assess strategies for steel reinforcement of concrete support structures [tam2015stress]. 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 load bearing structures is investigated in e.g., aerospace engineering, here stress directions provide the first indicators where structures can be hollowed [kratz2014tensor, kwok2016structural, daynes2017optimisation]. In biomechanics, such techniques are used to show tension and compression pathways simultaneously, and compare different structural designs regarding their mechanical properties [Dick2009Stress]
. For an overview of stress tensor visualization, we refer to the recent review article by Hergl et al.
[hergl2021visualization].An informative visualization of the stress directions in a 3D solid can be achieved via principal stress lines (PSLs), i.e., integral curves in 3D space along the principal stress directions. PSLs are effective in communicating the pathways along which external loads are transmitted, and they show the mutual relationships between the different principal stress directions [Dick2009Stress, Wang2020Globally]. In computational engineering, PSLs are used in particular to show where and how loads are internally redirected and deflected. 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, in computational mechanics stress trajectory visualizations are used in a rather inconsistent way, and, to the best of our knowledge, no standard tool for such an analysis exists. In many research groups in computational mechanics, own software packages for showing one particular principal stress direction starting at randomly selected locations are used. Often, 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 postprocessing 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 subdomains uncovered and does not show the mutual variations of the stress directions. In general, clutter can be reduced by visualizing the single stress directions sidebyside, yet juxtaposition makes it difficult to effectively relate the three mutual orthogonal stress directions to each other.
Contribution
This paper presents the 3D Trajectorybased Stress Visualizer (3DTSV), a system and methodology for the visual analysis of the PSLs in 3D stress fields. Figure 1
gives an overview of the visualization options provided by 3DTSV. With 3DTSV, we release a system that supports a comprehensive integral linebased analysis of 3D stress fields. To achieve this, 3DTSV builds upon existing techniques for line seeding in vector fields
[Jobard1997creating, IlluminatedStreamlines], and it extends them towards the specific use case by considering simultaneously the three principal stress directions in the seeding process. 3DTSV is designed to achieve improved regularity of the extracted PSLs, i.e., it aims for a gridlike structure where PSLs roughly intersect, uniformly cover the domain, and reveal symmetries in the underlying fields. To achieve this, 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 [Jobard1997creating, IlluminatedStreamlines], 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 levelofdetail (LoD) visualization that combines a dense seeding of a selected PSL type with a seeding at a userselected sparseness level of the respective other PSLs.To ease integration into existing systems and accessibility to end users, 3DTSV is implemented as a clientserver tool connecting a MatLab PSL extraction backend with an OpenGL rendering frontend. 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 due to the popularity of MatLab in mechanical engineering, and, thus, to enable engineers to easily integrate new model representations and algorithms. Currently, 3DTSV works with hexahedral simulation grids, including MatLab code for trilinear and inverse distancebased interpolation of stress tensors in such grids. If other types of basis functions are used, the corresponding MatLab functions simply need to be exchanged. Due to the cell adjacency structure that is built internally to efficiently find the next cell during trajectory integration in deformed hexahedral grids, other cell types can be supported with only minor additional effort.
The frontend renders whatever set of lines that is sent from the backend using advanced rendering options such as depth cues, outlines, as well as ambient occlusion effects to improve the perception of the spatial relationships between trajectories. Furthermore, the user can select to visualize one pair of stress directions via ribbons. Ribbons follow one of the selected directions and twist according to the other one, and they can effectively convey regions where the assignment of the eigenvector directions to the type of PSL (i.e., major, medium, or minor) changes.
To summarize, the contributions of this work are

an advanced and publicly available tool for trajectorybased stress tensor visualization supporting stress fields on arbitrary hexahedral grids,

the adaptation of evenly spaced line seeding to create a spacefilling set of PSLs with improved regularity,

an adaptive levelofdetail visualization using varying PSL density and visual mappings to lines and ribbons.
The application of 3DTSV is demonstrated in a number of experiments using datasets with different shapes and stress states. The code of 3DTSV is made publicly available under a BSD license, and published on https://github.com/JunpengWangTUM/3DTSV. In video1^{1}^{1}1https://youtu.be/lN9CxgvfgNY, the seeding of trajectories by 3DTSV is compared to the seeding of trajectories separately in each principal stress direction field via evenly spaced seeding [Jobard1997creating]. 3DTSV can be used as clientserver system as described (see video2^{2}^{2}2https://youtu.be/h7BzP7Jg_o), or as standalone tool solely in MatLab providing rudimentary visualization options (see video3^{3}^{3}3https://youtu.be/99Jn938ZoVk). Also the frontend can be used standalone, reading the PSL specific information from "psl.dat" files (see video4^{4}^{4}4https://youtu.be/zafBOAt9Xvs). Thus, any other backend can be used to generate PSLs and let the frontend visualize them. We also provide a script written in the ANSYS builtin language APDL, which automatically converts the result of an ANSYS finite element stress analysis into the format required by the 3DTSV backend (see video5^{5}^{5}5https://youtu.be/Yri_B7m3AWU). To support the output from ABAQUS, the mesh information needs to be read from the ABAQUS input file (".inp"), and the stress data can be acquired from the result file (".rpt"). We provide datasets, description and configuration files, as well as scripts for all use cases of 3DTSV on the publicly available GitHub repository.
2 Related work
Stress Tensor Field Visualization.
Stress tensor field visualization can be classified into trajectory, glyph and topologybased methods
[Kratz2013Visualization, hergl2021visualization]. Trajectorybased methods choose the PSLs as visual abstractions of the stress field, focusing on the directional structure of the principal stresses. Delmarcelle and Hesselink [Delmarcelle1993Visualizing] 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. Dick et al. [Dick2009Stress] 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. [Oster2018Core] proposed the concept of tensor core lines in 3D second‐order tensor fields. Hotz et al. [Hotz2006Tensor] smear out dye along the PSLs using line integral convolution. In this way, a density field is generated that resembles a gridlike 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.It’s worth noting that even though stresses are frequently simulated and analysed in engineering applications, the use of trajectorybased visualizations that consider the whole stress field as a tensor field instead of several scalar fields are not commonplace. In particular, such functionality seems neither provided by any of the wellestablished software packages for stress simulation, like ABAQUS and ANSYS, nor by dedicated environments for visualizing finiteelement simulation results [lee2008femvrml, weng2011web].
Glyphbased methods, on the other hand, depict the stress field by a set of welldesigned geometric primitives – socalled tensor glyphs. Tensor glyphs were originally designed for glyphbased diffusion tensor visualization [Kindlmann2004Superquadric], and later adapted to visualize positive definite tensors [Kindlmann2008Quantification], general symmetric tensors [Schultz2010Superquadric], as well as asymmetric tensors[Seltzer2016Glyphs, Gerrits2017Glyphs]. Glyphbased 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 [Kindlmann2006Diffusion, Hlawitschka2007Interactive]. Tensor glyphs are effective in showing the local stress states, but they cannot effectively communicate the global structure of stress lines. Patel and Laidlaw [patel2020visualization] 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.
Topologybased approaches for stress tensor visualization abstract from the depiction of stress directions and focus on revealing specific topological characteristics of the tensor field. Delmarcelle and Hesselink [Delmarcelle1994Topology, Hesselink1997Topology] studied the topology of symmetric 2D and 3D tensor fields, and introduced the fundamental concepts of degenerate points and topological skeletons. Zheng and Pang [Zheng2004Topological], and later Roy et al. [Roy2018Robust], 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 [Zobel2018Extremal]. Raith et al. presented a general approach for the generation of separating surfaces in the invariant space [Raith2018Tensor]
. Palacios et al. introduced the eigenvalue manifold and visualized the 3D eigenvectors as curve surfaces
[Palacios2015Feature]. Qu et al. [Qu2020Mode] 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 [TurkStream] and Jobard and Lefer [Jobard1997creating] were the first to introduce seeding strategies for generating evenly spaced sets of streamlines in 2D vector field. Numerous extensions and improvements of these concepts have been proposed since then. In particular, Vilanova et al. [Vilanova:2004:DTI] 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 and frameworks have been developed to reduce the visual clutter and occlusion of densely distributed streamlines in 3D fields [Ye2005Strategy, chen2007similarity, Yu2011Hierarchical, Kanzler2016Line]. 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 and then shaded, e.g., using the BlinnPhong shading model [BlinnPhong]. Early work on illuminated streamlines was done by Zöckler et al. [IllumStreamlines] and Mattausch et al. [IlluminatedStreamlines]. Stoll et al. [StylizedLines] extended this work by introducing stylized line primitives, rendered by a hybrid CPUGPU renderer. Liu [liu2019prototype] presented the DOXIV, a prototype framework for highperformance visual analysis of large flow data. Volpe [StreamribbonsInitial] first introduced the concept of streamribbons for flow field visualization.
Hexahedral Meshing. An alternative approach to PSLbased stress field visualization is to generate a frame field from the principal stress field first and employ fieldaligned hexahedral meshing to produce orthogonal edges that follow PSLs. The edges of such hexmeshes 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 allhex mesh. Hexahedraldominant meshing has been resorted as an intermediate solution. For instance, Wu et al. [Wu2021TVCG] propose a conforming stressguided lattice structure by combining topology optimization with the fieldguided polyhedral meshing algorithm from [Gao2017Robust]. Arora et al. [Arora2018DEsigning] 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 Tjunctions or nonhexahedral elements with nonorthgonal edges, significantly deviating from the PSLs and are, thus, not applicable for stress field analysis either.
3 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 secondorder stress tensor
(1) 
contains these vectors for the axes of a Cartesian coordinate system.
is symmetric since the shear stresses given by the offdiagonal elements in 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 , 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 represent the major , medium and minor 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, , and are mutually unequal, and the eigenvectors are linearly independent and even mutually orthogonal due to the symmetry of . However, socalled degenerate points can exist where two or more eigenvalues are equal. In the vicinity of these points, which are classified by or ^{6}^{6}6We do not consider triple degenerate points with , since they do not exist under structurally stable conditions [Zheng2004Topological]., the PSL direction cannot be decided. Therefore, when tracing along a principal stress direction, we test whether the eigenvalue corresponding to this direction is too close to another eigenvalue , i.e., . If this is the case and the angle between the PSL tangents at the current and next integration point is too large, the integration is stopped. Furthermore, we provide the option to map to the color of a PSL via a color table (see Sec. 4.3), so that the proximity to a degenerate point is indicated. PSL 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 (i.e., to avoid running into closed orbits), or the number of integration steps reaches the predefined threshold.
The integration of PSLs requires to select seed points from which they start until they arrive at a degenerate point or the boundary. 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 vertices where loads are applied. Furthermore, different integration schemes can be used for PSL tracing, including the 1storder Euler method, and the 2nd and 4thorder RungeKutta methods, where the fixed integration step size is used for Cartesian meshes, and an adaptive for unstructured hexahedral meshes. In each integration step, the stress tensor is interpolated, and the eigenvalues and eigenvectors are computed 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.
4 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 builtin regularity due to enforced PSL intersections.
4.1 Evenly Spaced PSL Seeding
The proposed seeding strategy builds upon the evenly spaced streamline seeding approach by Jobard and Lefer [Jobard1997creating], and extends 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 gridlike 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 gridlike 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 gridlike 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 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 (passed by no PSL), solid seed (passed by both major and minor PSLs) and semiempty seed (only passed by major PSL) or (only passed by minor PSL). The sampling process is repeated until all valences of all possible seed points become solid . 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 spacefilling, the initial candidate seed points (with ) are located at the vertices of a spacefilling Cartesian grid (step 0 in Figure 2).
Seeding starts by selecting one of the candidate seed points and tracing the major and minor PSLs from it (Step 1 in Figure 2), setting 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 Figure 10 and Figure 11). Then, all candidate seed points with not equal to are reclassified 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 or , respectively. If a point is classified as or and closer to a minor or major PSL, respectively, its valence is set to . 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 reclassified candidate point is relocated (i.e., merged) 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 merging 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 Figure 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 or
, 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 25 in
Figure 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 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 semiempty seed points after merging, 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 . 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.
4.2 PSL LoD Structure
To change the density of the generated PSLs, the seeding process can simply be rerun 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 , the refined set of PSLs according to a distance threshold is computed as follows: Firstly, the candidate seed points are reset to their initial positions. Secondly, the candidate seed points are merged to the existing PSLs according to , to create “empty” bands around the existing PSLs. The valences are updated accordingly to , or depending on the types of PSL they are merged to. After this, some nonsolid seeds are left, because is larger than . With these seeds the seeding is subsequently performed, including the iteration of seed point selection, PSL computation, and reclassification as described in subsection 4.1.
To generate a full LoD PSL hierarchy, the user defines the minimum distance threshold and the number of levels to construct. Then, the distance thresholds of each level are computed as from coarse to fine, and the hierarchy is constructed progressively from the coarsest resolution level (see 1st and 2nd rows in Figure 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 multitype LoD is computed by alternatively considering the different PSL types with their respective distances.
4.3 Ribbonbased Stress Visualization
Instead of rendering lines, the user can select a PSL type (i.e., major, medium, minor) and visualize ribbonshaped geometry [Streamribbons] that is centered at the PSLs of the selected type and twists according to the direction of another stress type (see Figure 5 a,b). At each integration point along a PSL of the selected type, two lines with adjustable length are traced forward and backward along the other direction. The lines’ endpoints at subsequent integration points are connected to form a ribbon. It is worth noting that the constructed ribbons don’t coincide with streamsurfaces that are integrated from a PSL along one other stress direction. As shown by Raith et al. [Raith2018Tensor], such surface might not even exist, i.e., when integrating from two points on the same PSL over a certain length along another stress direction, the two endpoints are not lying on a PSL in general. The mapping of two principal stress directions to a ribbon geometry is conceptually similar to the wellknown hyperstreamlines [Delmarcelle1993Visualizing], i.e., a mapping of two principal stress directions to a tube centered at the PSL along the third direction.
We let the user select a visualization using ribbons to convey changes in the assignment of the eigenvector directions to the type of PSL in the vicinity of degenerate points. When a ribbon is formed as described, flips often occur in the vicinity of a degenerate point (see Figure 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 an additional visual cue to indicate topological changes of the PSLs in the vicinity of degenerate points.
Figure 6 compares the options to visualize principal stress directions via ribbons and lines, and combine them into a single visualization. As can be seen, twists in the ribbon geometry effectively hint to regions where degenerate points might exist. For lines, 3DTSV 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 wellknown inaccuracies at curved boundaries that are represented by hexahedral simulation elements in a Cartesian grid.
5 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. 3DTSV relies on the requestreply 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.
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 multithreaded 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 hardwareaccelerated features from the C++ frontend, like depth cues or ambient occlusion effects. 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.
5.1 Numerical PSL Integration
3DTSV is designed to support the visualization of PSLs in solids discretized by hexahedral grids, where the stress tensors are given at the grid vertices. When computing PSLs in Cartesian grids, componentwise trilinear interpolation of the tensors is used during numerical line integration. In deformed hexahedral cells, tensor interpolation is performed via inverse distance weighting [shepard1968two].
To integrate PSLs in Cartesian grids, the system provides fixedstep integration schemes with user adjustable stepsize of at least half the cell diameter. In deformed hexahedral grids, a different approach is taken since the size of the simulation elements can vary, and with a constant stepsize the risk increases that multiple cells smaller than this size are missed in one single integration step. To reduce this risk, the integration stepsize is automatically adapted to the size (i.e., the length of the shortest edge) of the cell at the current integration point . These values are precomputed and stored per cell. In each integration step, the size of the current cell is read and multiplied by a user selected scaling factor . can be made smaller than 1 to obtain more accurate PSLs. With the stepsize , the PSL is integrated from the current point in cell to the new point . Then, the integration process is restarted with and the cell containing .
To find , it is first tested whether is still contained in . The following inout criterion is used to test whether a point is located in a hexahedral cell: Given a hexahedral element with the centers and outfacing normal of its 6 faces and . Any point in the interior or on the boundary of the element satisfies , see Figure 7a. In practice, the criterion is slightly relaxed to , to account for nonplanar cell faces, i.e., a slight variation of the normal vectors across the faces.
If does not contain , the cell needs to be determined. To this end, we further test whether lies in any of the adjacent cells of . For each cell, the set of adjacent cells as well as the adjacency type, i.e., face, edge, and vertexadjacency, is precomputed and stored. In case is not within or , we scale down the stepsize via a dichotomy strategy, i.e., , until is located in or it’s adjacent cells .
In the case where and are connected by a single edge or vertex, it may still happen that cells are skipped when going from to . In this situation, stepsize refinement is performed multiple times until the cell shares a face with or is below a userselected threshold. The latter situation is encountered when the PSL goes through a cell vertex or edge, so that faceadjacency cannot be determined. In Figure 8, for the given mesh two PSLs that have been extracted without and with additional stepsize refinement are compared. As can be seen, cells that would be skipped when using only facetoface adjacency are now determined and considered in the integration.
5.2 Rendering
The line and ribbon primitives are rendered in a stylized fashion similar to the techniques by Zöckler et al. [IllumStreamlines], Stoll et al. [StylizedLines] and Mattausch et al. [IlluminatedStreamlines], using default colors, halos and depth cues as shown in the first three images in Figure 1. Focus PSLs and contextual ribbons are rendered in ocher and blue, respectively. The base color is modulated using BlinnPhong shading [BlinnPhong, IllumStreamlines], 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 Figure 1 using the scalar von Mises stress measure. The scalar values are issued via the backend as pervertex attributes. The standard color scheme we use for the different principal stress directions (blue, green, ocher) is the ‘3class Set2’ transfer function from ColorBrewer^{7}^{7}7https://colorbrewer2.org/#type=qualitative&scheme=Set2&n=3. 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. A translucent simulation mesh outline hull can be rendered together with the stress field data in order to hint at the extents of the simulation domain.
5.3 3DTSV Settings
Data Set  #Cells  #Seeds  #PSLs  Time (s)  
Cantilever  250K  2K  1/5  1  85  0.4 
Rod  536K  18K  1/5  1  174  2.1 
Femur  696K  10K  1/18  3  823  9.0 
Bracket  650K  9K  1/12  3  293  5.4 
Bearing  189K  55K  1/18  3  1,364  33.4 
Parts1  253K  46K  1/20  3  1,557  27.9 
3DTSV 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 introduced in subsection 4.1 and subsection 4.2, respectively. Another set of parameters enables a userguided 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 two PSL types that are used to generate ribbons. Via a dropdown menu, the user can select a scalar stress measures that are 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.
6 Results
In all of our experiments, 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 onchip memory. The rendering times are always below 10 milliseconds. The data sets we use in our experiments are shown in Figure 9. The stress fields are simulated by a finite element method (FEM), using the solid objects under the shown load conditions. Table 1 lists the numbers of simulation elements of each of the data sets, the seed points that are used to generate the PSLs, the number of generated PSLs, and the time required for PSL generation.
For the three models ’Bridge’, ’Cantilever’ and ’Rod’, we demonstrate the improvements of the proposed seeding strategy over evenly spaced streamline seeding. 3DTSV is used to visually analyze the stress fields in ’Femur’ and ’Bracket’. These two data sets that are frequently seen in structural design and optimization [Wu2018TVCG]. Finally, we consider the two mechanical parts ’Bearing’ and ’Parts1’ to demonstrate the application of 3DTSV to unstructured hexahedral simulation meshes.
Figs. 10 and 11 emphasize the improvements by the proposed seeding strategy regarding the regularity of the extracted set of PSLs. 3DTSV generates a fairly uniform spacefilling PSL structure, which, in particular, maintains the symmetry of the stress field in ’Cantilever’. Evenly spaced streamline seeding, on the other hand, generates a far less regular design which introduces severe visual clutter.
The visualization also highlights the importance of showing different PSL types simultaneously. In the analyzed tensor field, the signs of the eigenvalues along the major and minor PSLs are mostly positive and negative, respectively. This means that the major PSLs are mainly under tension and the minor PSLs mainly under compression. Thus, either of both effects could be shown by visualizing one PSL type, but not both.
Figure 12 (top) shows the spacefilling PSLs in the stress field in the interior of ’Bracket’. From the boundary condition in Figure 9, we see that the structure is mainly under tension. Thus, we choose to show the major PSLs at the higher level of detail () and the minor PSLs at lower level (see Figure 12 (bottom)). The minor PSLs are shown via ribbons, with the medium principal stress direction indicating the twist. This enables a fine granular analysis of the major principal stress directions, and simultaneously provide a coarse representation of the other principal directions. A similar setting has been selected to visualize the stress directions in ’Femur’ (see Figure 1).
3DTSV works with Cartesian meshes and deformed hexahedral meshes, which are both frequently used in mechanical engineering applications. Here we use the stress fields due to external loads in the interior of ‘Bearing’ and ‘Parts1’, to demonstrate the capability of 3DTSV. As shown in Figure 9, especially in ‘Bearing’ the element sizes change considerably over the 3D domain. The distribution of PSLs of ‘Bearing’ is shown in Figure 13 (top), and the bottom image shows the combination of major at the third level of detail () and minor at , where the minor PSLs are shown via ribbons. The full distribution of PSLs of ‘Parts1’ can be seen in the Figure 14 (left), on the right the minor PSLs at and major PSLs at are shown simultaneously, where the major PSLs are rendered via ribbons.
7 Conclusion and Future Work
In this paper, we have introduced 3DTSV, a tool for visualizing the principal stress directions in 3D solids under load. 3DTSV makes use of a novel seeding strategy, to generate a spacefilling 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 multiresolution hierarchy is formed, which can be utilized to show different PSL types with different resolutions simultaneously. 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 3DTSV 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 multithreading architectures. Furthermore, we are interested in using spacefilling 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 gridlike structure, we aim at combining our seeding strategy with the automatic growth process underlying topology optimization.
Acknowledgment
This work was supported in part by a grant from German Research Foundation (DFG) under grant number WE 2754/101. We acknowledge the help of Chunxiao Meng at Northwestern Polytechnical University and Yingjian Liu at The University of Texas at Dallas in adapting 3DTSV to ANSYS and ABAQUS, respectively.