Interpolation methods for orthotropic fourth-order fiber orientation tensors in context of virtual composites manufacturing

The objective of this work is to develop and investigate interpolation methods for fourth-order fiber orientation tensors. The developed methods aim to minimize information loss during mapping in the virtual manufacturing process for fiber-reinforced composites. The strategy of decomposing features that describe tensor shape and orientation, followed by separate interpolation, has been found to be a suitable approach for interpolating second-order fiber orientation tensors. Thus, in this work, decomposition-based approaches are also applied to fourth-order fiber orientation tensors. The decomposition of orthotropic fourth-order tensors into shape-and orientation-describing features with respect to three-dimensional space is presented. Various proposals for shape interpolation are described and evaluated. The developed interpolation methods are then applied to the numerical solution of a continuous problem from fluid mechanics. The continuous solution is used to generate data for the associated discrete problem and to evaluate the interpolation results of the different methods. The evaluation results demonstrate that the proposed interpolation methods are clearly superior to the conventional interpolation methods. Separating shape-and orientation-describing features with respect to three-dimensional space also proves to be useful for fourth-order tensors. The developed methods can significantly reduce the loss of information during mapping in virtual manufacturing process chains.


Notation
Arbitrary sized numerical arrays, such as components in a specific coordinate system, are denoted by indices (e.g.,   ,   ,   ) where the range of index values corresponds to their dimensionality.Spatial components in the three-dimensional space are denoted by Latin lower case letters , , … ∈ {1; 2; 3}.Iterators in a set of discrete values are denoted by upper case Latin letters, e.g.,  , … ∈ {1; 2; 3; ...}.Unless otherwise indicated, Einstein's convention for summation holds.This means that indices appearing twice in a single expression imply summation.The transition from symbolic to index notation is represented by the operator = under the assumption of orthogonal basis vectors.
Symbolic tensor notation is preferred throughout this work.Scalars are denoted by standard Latin and Greek letters, e.g., ,  .First-order tensors are represented by bold lower case letters, e.g., , whereas upper case Greek or Latin letters are used for second-order tensors such as , .Fourth-order tensors are denoted by C, S. Tensors of higher-order than four are denoted using superscript indices, e.g.C ⟨⟩ , where  represents the order of the tensor.The composition of two equal-order tensors is denoted without special operator symbols, e.g. =     , AB =     .The composition of two tensors of different order is denoted by a circle, e.g.•B =     , A• =     .All contractions of higher-order tensors by a lower-order tensor, The trace of a tensor of arbitrary order is defined as the projection of the tensor on its according unit tensor (tr() =  ⋅  and tr(A) = A ⋅ I).The identity of symmetric second-order tensors is referred to as I S .The rotation of second-or fourth-order tensors is denoted using the Rayleigh Product  ⋆  and  ⋆ A, where the second-order tensor  is an element of the orthogonal group ℎ ⟨3⟩ .Components of second-or fourth-order tensors that are displayed in Mandel notation [1] are denoted using a hat and Greek letters as indices, e.g.Â , Ĉ .For Mandel notation, Einstein's convention for summation is applied for all indices appearing twice with , , … ∈ {1; … ; 6}.All conventions concerning the Mandel notation refer to Bauer and Böhlke [2].
Throughout this work various sets of tensors are utilized.All considered first-order tensors are part of the three-dimensional vector space .The set of all symmetric second-order tensors is called .If a second-order tensors  satisfies the condition   = , it is called orthogonal and is an element of the orthogonal group ℎ ⟨3⟩ .A fourth-order tensor that possesses both the major symmetry (A = A   ) and minor symmetries (A = A   = A   ) is called Hooke's tensor.The set of all Hooke's tensors is referred to as .The set of all fourth-order tensors that provide full index symmetry (  =   =   =   = … ) is referred to as  ⟨4⟩ .Fourth-order tensor R that satisfy the condition RR   = I S are called orthotropic.

Motivation
This work addresses a methodological problem that is relevant, but not limited to, the holistic modeling of discontinuous fiber reinforced polymers (DiCoFRP).Numerous contributions have documented the crucial influence of process-induced fiber orientation on mechanical properties [3][4][5].To this end, component design is usually accompanied by a CAE-chain interlinking process simulation and structural simulation.Deviating numerical challenges and requirements necessitate different spatial discretizations in the individual simulation modules (cf.Fig. 1 left).Especially if the results are transferred from a coarser source mesh onto a higher-resolution target mesh, as indicated by blue respectively green colors in Fig. 1, interpolation schemes have to be applied to sustain smoothness.In nearly all available macroscopic process simulation approaches, the spatio-temporal evolution of fiber orientations is modeled via direction-dependent second-order fiber orientation tensors following the fundamental work of Advani and Tucker [6] and Kanatani [7].In contrast to scalar fields, their interpolation harbors additional challenges.Treating tensor components as independent scalar fields leads to systematic overestimation of isotropy and more sophisticated decomposition-based schemes are favorable [8][9][10].
Within the interface between process simulation and structural simulation, the fiber orientation tensors are used to determine effective macroscopic properties via an orientation average.For elastic properties, the knowledge of the second-order fiber orientation tensor field is not sufficient [5].Instead, a fourth-order fiber orientation tensor field has to be estimated algebraically by means of a closure operation.Taking only the second-order fiber orientation tensor into account, most closure approaches preserve the material symmetry of its argument.Therefore, the weakest material symmetry a closed fourth-order fiber orientation tensor can possess is orthotropy [2].Closure approaches have been and continue to be subject of research, cf.[11] for a comprehensive overview.
Recently, Kraußet al. [12] gave strong reasoning that the closure of averaged second-order fiber orientation tensor may yield invalid results and, instead, the closed fourth-order fiber orientation tensor should be averaged.The study in Kraußet al. [12] was conducted under the assumption of a continuous, piecewise constant fiber orientation tensor field, i.e. for any spatial point the specific value of the fiber orientation tensor is known.Yet, in various numerical schemes the spatial domain is discretized and an approximated solution is only available at specific coordinates with a finite resolution.If source and target locations mismatch, interpolation techniques become necessary.Ultimately, this leads to the questions, for which this work seeks to give answer to: How should fourth-order interpolation techniques be constructed to preserve meaningful physical quantities, and how do they perform in qualitative and quantitative comparison with a naive interpolation of tensor components?All fourth-order fiber orientation tensors acting as basic values in the interpolation are assumed to be result of a valid closure scheme.Therefore, this work is limited to orthotropic fiber orientation tensors.

Related work 2.2.1. Interpolation of second-order tensors
The largest application area for tensor interpolation is medicine.Here, it is mainly used for diffusion-weighted magnetic resonance imaging (DW-MRI).DW-MRI was first mentioned by Le Bihan et al. [13] and is a noninvasive diagnostic imaging technique that can map the diffusive motion of water molecules in the brain.The local diffusivity is thereby mapped by a diffusion tensor  DW [14].By evaluating the diffusion movements, physicians can draw conclusions about diseases of the central nervous system [15][16][17].Using DW-MRI, the diffusion tensor  DW is only evaluated at specific spatial points in the brain.However, the physician requires high resolution information on the orientation and shape of  DW .Thus, the diffusion tensors must be interpolated between these low resolution points to draw accurate diagnoses [8].
In the mechanical context, tensor interpolation has been considered less frequently.In context of CAE chains, it is needed to transfer tensor-valued information from a source mesh to a target mesh [18,19].To minimize information loss during transfer, Krauß and Kärger [9] examined various interpolation techniques for symmetric second-order tensors.In industry, commercial software is commonly used for these tasks.Examples include Converse from PART Engineering GmbH, MpCCI MAPPER from Fraunhofer SCAI, and Digimat-MAP from MSC Software GmbH.
In the current paper the considered second-order tensors are second-order fiber orientation tensors .They were first introduced by Kanatani [7] and are a statistical measure to describe the orientation distribution of fibers at a central moment.Second-order fiber orientation tensors  are the second moment of the orientation distribution function  () Here,  is a unit vector that characterizes the orientation of a specific fiber.Since  () is a probability density function, it is normalized over all possible orientation states.Thus, the trace of the second-order fiber orientation tensor is normalized tr() = 1.Also,  () is always non-negative.This means that all eigenvalues of  must be non-negative; respectively,  is positive-semidefinite, so the condition ⋅ ≥ 0 holds for all  ≠ .Note that the unit vectors  can be permutated arbitrarily in Eq. ( 1).Hence,  must be symmetric.If  fibers with known orientation  are weighted equally by means of the distribution function  (), the second-order orientation tensor  can be calculated empirically Interpolation methods for tensors are divided into global and decomposition-based interpolation techniques.The simplest global interpolation technique is the Euclidean interpolation (EU).The result of the EU is the weighted arithmetic averaging of the individual tensor components in the global coordinate system Here, the tensor components  , are considered as independent scalar fields of the  tensors.The factors   are the weights of the supporting points.The interpolation is called Euclidean because ĀEU  minimizes the Euclidean distance with respect to the Frobenius norm [9].An alternative global interpolation technique is the log-Euclidean interpolation, which was introduced by Arsigny et al. [20].Here, the averaged tensor is computed by Log-Euclidean interpolation is a generalization of the geometric mean for scalars.This interpolation method is only valid for positive definite second-order tensors [21].
In contrast to global interpolation techniques, decomposition-based interpolation techniques separate the tensor orientation from the tensor shape.Afterwards, orientation-describing and shape-describing features are interpolated separately and finally reassembled into the final interpolation result.The shape-describing features must be tensor invariants.This means that they are invariant with respect to all orthogonal transformations  ∈ ℎ ⟨3⟩ .The best known example for tensor invariants are the eigenvalues  ⟨⟩  , which result from the eigenvalue problem for second-order tensors (Section 3.1).Krauß and Kärger [9] showed that the use of orthogonal invariants, introduced by Ennis and Kindlmann [22], are useful to average the shape-describing features.In this context, the invariants   were selected to fulfill the orthogonality condition [23] After the linear interpolation of the chosen invariants, the interpolated tensor invariants can be used to infer the eigenvalues λ⟨⟩  of the associated tensor.The calculation rules needed for this can be found in Gahm et al. [8] and Gahm et al. [24].The averaged shape is described by the diagonal tensor Λ. Shape averaging using the orthogonal invariants leads to a monotonic change of the tensor shape between the grid points, which cannot be achieved with the global interpolation techniques [9].For orientation averaging of second-order tensors, the orthogonal tensor  ∈ ℎ ⟨3⟩ obtained from the spectral decomposition (Section 3.1) can be used.Here,  ∈ ℎ ⟨3⟩ describes the orientation of the eigensystem of  relative to the considered global coordinate system.The orientation information contained in  can be interpolated.To this end, Gahm and Ennis [25] suggest two options.The first one is orientation averaging using quaternions.Quaternions provide the possibility to distinctly describe rotations using a tuple of four real numbers.An iterative algorithm for averaging quaternions is described by Krauß and Kärger [9].An alternative procedure for averaging tensor orientations is based on the eigenprojector representation of second-order tensors [e.g.26].The projectors of the orthogonal tensors are arithmetically averaged at the supporting points.Using the polar decomposition [e.g.27], the rotational part can be filtered.For extended information on orientation averaging using eigenprojectors, see Gahm and Ennis [25].The averaged orientation is described by the rotation R ∈ ℎ ⟨3⟩ .Krauß and Kärger [9] reassemble the averaged orientation R and the averaged shape Λ to the final interpolation result using Note that Eq. ( 6) corresponds to the inverse spectral decomposition for second-order tensors.The main result of the investigations by Krauß and Kärger [9] was that a separation of shape-and orientation-describing features is a suitable approach to interpolate tensor-valued properties.Moreover, the averaging of orthogonal invariants provides a monotonic change of the tensor shape, which cannot be achieved using simple arithmetic averaging of tensor components (e.g.Eq. ( 3)).
The decomposition-based approach was used by Blarr et al. [10,28] to generate high resolution orientation tensor fields from scarce input data.The input data was determined from microscopic X-ray computed tomography scans.The main result was that the implemented methods provide good macroscopic interpolation results without using high amounts of resource-consuming micro-computed tomography (CT) scans.However, it is assumed that the decomposition-based method provides poorer results for extrapolation tasks.

Interpolation of fourth-order tensors
Evaluation of DW-MRI in complex brain structures often leads to inadequate results [14].Complex brain structures include fiber bundles, fiber crossings, and branched fibers [29].At these locations, more information on local diffusivity is needed than is contained in the second-order diffusion tensor  DW .Methods are needed to interpolate this additional information.In the field of DW-MRI, there is little knowledge on how to meaningfully interpolate fourth-order tensors.Currently, research is still focusing on the description of complex diffusion states.The problem of describing complex diffusion states was first addressed by Tuch [30].Diffusion states are described using a fiber orientation distribution function, which describes the local orientation state of brain fibers [31,32].The method is called ''high angular resolution diffusion tensor imaging'' (HARDI) and is the topic of various recent publications [e.g.33,34].
In context of HARDI, an -th-degree tensor describes a homogeneous th degree polynomial whose coefficients approximate the fiber orientation function  ().Therefore, in context of HARDI, the probability density function  () should not be interpreted as a Fourier series (as by Advani and Tucker [6]).It rather corresponds to a truncated Taylor series, which depends on the coordinates of the longitudinal fiber axis .
Most publications on spatial tensor interpolation are limited to second-order tensors (Section 2.2.1).Interpolation methods applicable to higher-order tensors have been considered rarely.Two decomposition-based interpolation techniques, that use the Tucker decomposition and the canonical decomposition, are considered by Cardona et al. [35] and Vargas-Cardona et al. [36] in the context of DW-MRI.The tensors from those decompositions are high-dimensional feature-tensors.Data from the decompositions are used as training data for a Markov chain Monte Carlo procedure.The canonical decomposition corresponds to an eigenprojector decomposition for fourth-order tensors and is superior to the Tucker decomposition in terms of interpolation results.To evaluate the results, the two decomposition-based interpolation techniques were compared to the Euclidean averaging of the tensor components (see Eq. ( 3)) The Frobenius norm was used as the error metric.The main result of Cardona et al. [35] and Vargas-Cardona et al. [36] was that both decomposition-based interpolation techniques are superior to the linear interpolation of tensor coefficients (Eq.( 7)).Bauer et al. [37] recently investigated the interpolation of fourth-order fiber orientation tensors.They proposed an interpolation method that separates the information of the tensors's shape and the orientation by means of a parametrization that is based on tensor components and a unique eigenvalue system.They also discussed a convention for defining a unique eigensystem in the absence of material symmetries.
In this paper the considered fourth-order tensors are the fourth-order fiber orientation tensors A. Similarly to the second-order fiber orientation tensors  (Eq.( 1)), they were first introduced by Kanatani [7].They are a statistical measure to describe the orientation distribution of fibers at a central moment.Fourth-order fiber orientation tensors A are the fourth moment of the orientation distribution function  () The function  () is always non-negative.Thus, the tensor A is positive-semidefinite and the condition  ⋅ A[] ≥ 0 holds for all  ≠ .Note that the unit vectors  can be permutated arbitrarily in Eq. ( 8).Hence, A must be fully symmetric.If  fibers with known orientation  are weighted equally by means of the distribution function  (), the fourth-order orientation tensor A can be calculated empirically Higher-order fiber orientation tensors are linked with lower-order fiber orientation tensors by Consequently, all information of a lower-order orientation tensor can be derived from the higher-order orientation tensors which describe the same orientation state.

Scope of this work
This paper focuses on interpolation techniques for fourth-order fully symmetric tensor fields, such as fiber orientation tensors.The techniques presented were developed in context of virtual manufacturing workflows for composites, but they are not limited to this domain.Note that the interpolation schemes presented are limited to orthogonal material symmetries.This is due to the fact that most fiber orientation tensors related to the virtual manufacturing process are at least orthotropic.Injection molding simulation usually exclusively provides second-order fiber orientation tensors.Hence, the associated fourth-order fiber orientation tensors are constructed using closure approximations which only provide orthotropic material symmetry.
The aim of this work is to minimize systematical errors when transferring direction-dependent field data between deviating discretizations while also maintaining various mathematical and physical conditions.The developed interpolation schemes must satisfy the following conditions: • the full index symmetry is preserved; • the trace conditions tr() = 1 and tr(A) = 1 are preserved; • the interpolation of two supporting points that only differ in orientation-describing features must preserve the tensor shape; • the interpolation of two supporting points that share material symmetry must preserve the material symmetry.This paper is structured as follows: Section 3 focuses on orthogonal decompositions for fourth-order tensors that are utilized to develop decomposition-based interpolation techniques.The underlying mathematics and simplifications for fiber orientation tensors are briefly described and applied.Section 4 focuses on the development of decomposition-based interpolation techniques for orthotropic fourth-order fiber orientation tensors.For this, the extraction of orientation-and shape-describing features of orthotropic fourth-order fiber orientation tensors is described.Afterwards, multiple interpolation schemes to interpolate the shape-describing features are presented, while the interpolation schemes for orientation-describing features can be adopted from Krauß and Kärger [9].In Section 5, a continuous solution for a theoretical example from fluid dynamics is derived.This continuous solution is used to generate data for the associated discrete problem on which the developed interpolation techniques are applied.Finally, the performances of the interpolation techniques are compared using various error metrics.

Spectral decomposition
For symmetric second-order tensors  ∈  the spectral decomposition is a well-known orthogonal decomposition.The spectral decomposition for second-order tensors is based on the solution of the eigenvalue problem.The eigenvalue problem for symmetric second-order tensors  reads (no summation over  = 1, 2, 3) Here,  denotes the zero vector.The vectors   are called eigenvectors of  [26].If  is displayed in the basis system which is constructed by its eigenvectors   , the eigenvalues  ⟨⟩  are on the diagonal of the associated component matrix The representation in Eq. ( 12) is called spectral decomposition and can be denoted in an alternative way by where  is a diagonal second-order tensor Note that the component matrices in Eq. ( 12) and Eq. ( 14) are identical while the basis systems differ.Hence, the tensor  is an orthogonal second-order tensor which describes the rotation from the global basis system { 1 ,  2 ,  3 } to the (normalized) eigenvector system { 1 ,  2 ,  3 }.However, the rotation tensor  is not distinct if no convention for the sequence of the eigenvalues  ⟨⟩  is applied.Since the trace of a second-order orientation tensor  is normalized, the sum of the eigenvalues always equals one For fourth-order tensors A, the spectral decomposition is less well-known.It is defined for Hooke's tensors A ∈  that are positive-semidefinite (no summation over  = 1 … 6) [38,39] Here,  denotes the second-order zero tensor.The characteristic polynomial for the eigenvalue problem in Eq. ( 15) is defined In this work the determinant of a fourth-order tensor A ∈  is equal to the determinant of its six-by-sixmatrix in Mandel notation det(A) = det( Â ).Thus, the values  ⟨A⟩  that solve the characteristic equation (A) = 0 can be calculated in Mandel notation [39] and are the eigenvalues of A. The values  ⟨A⟩  are also referred to as Kelvin moduli.They were first described by Thomson [40].The symmetric second-order tensors   that are associated to the eigenvalues  ⟨A⟩  are the eigentensors of A. In the mechanical context the tensors   are also referred to as eigenstates [e.g.41,42].The eigentensors   are pairwise orthogonal if all associated eigenvalues  ⟨A⟩  are unique.For fourth-order tensors A ∈  the spectral decomposition can be defined analogously to Eq. ( 12) For the eigenvalues of an orientation tensor A the additional condition holds.The equation that is analogous to Eq. ( 13) can be formulated in Mandel notation The diagonal of the matrix Λ is populated with the six eigenvalues of A. In the following the transformation of R to tensor notation is referred to as R ⟨SP⟩ with where   are the basis vectors of the Mandel notation [e.g.1,43,44].
demonstrates that R ⟨SP⟩ is an orthogonal tensor.

Harmonic decomposition
The decomposition of symmetric second-order tensors  ∈  in isotropic and symmetric deviatoric parts is well-known.Since the trace of a second-order fiber orientation tensor  is normalized, its isotropic part is constant  • = 1  3 .The symmetric deviatoric part is traceless (tr( ′ ) = 0) and reads  ′ =  −  • .Both parts are symmetric and orthogonal ( The decomposition of fourth-order tensors in its isotropic and its deviatoric parts is called harmonic decomposition.It is often applied for Hooke tensors A ∈  [e.g.[44][45][46][47][48].The general form of the harmonic decomposition reads [47] where the three parts A Iso , A Dev 2 and A Dev 4 are biorthogonal and distinct.Further decompositions of the general harmonic decomposition are outlined by Rychlewski [47].Since fiber orientation tensors are fully symmetric and the condition  ⋅ A[] = 1 holds, the harmonic decomposition can be simplified.The general form of the harmonic decomposition for fourth-order fiber orientation tensors A is given in irreducible form [2] Eq. (22) shows that the harmonic decomposition of fourth-order fiber orientation tensors can be described with only one secondorder deviator  ′ and one fourth-order deviator A Dev 4 .Additionally, the isotropic part A Iso is constant for all fourth-order fiber orientation tensors.Further decompositions addressing the harmonic decomposition for fiber orientation tensors are presented by Bauer and Böhlke [2].

Application of interpolation techniques for second-order tensors on fourth-order tensors using spectral decomposition 4.1.1. Calculation rule
Krauß and Kärger [9] demonstrated the advantages of a decomposition-based interpolation method for second-order tensors.To achieve this, they decomposed the tensors into orientation-and shape-describing features, interpolated the two features separately, and finally reassembled them.The averaged orientation is described by the orthogonal tensor R ∈ ℎ ⟨3⟩ .The averaged shape is described by the diagonal tensor Λ, which contains the eigenvalues of the interpolated tensor.The averaged tensor is obtained by reassembling using the inverse spectral decomposition (Eq.( 6)).
In the following, a procedure is presented to directly extend the decomposition-based interpolation approach from second-order tensors [9] to fourth-order tensors using Mandel notation.For the interpolation, each of the  tensors at the supporting points is decomposed using the spectral decomposition for fourth-order tensors (Eq.( 18)).The coefficient matrices Λ, are diagonal and are intended to describe the shape of the tensors A  in the six-dimensional Mandel space.They represent the extension of the diagonal tensor  from Eq. ( 13).On the diagonal of the matrix Λ, are the Kelvin moduli of the according fourth-order tensor A  .The transformation of the matrix R, into tensor space is denoted by R ⟨SP⟩  .The matrices R, are used in this interpolation method as an extension of the rotation tensor  from equation Eq. (13).The supporting points are weighted using the scalar weights   .
The interpolation of the tensor shapes can be performed by the linear interpolation of the Kelvin moduli Λ, .The averaged tensor shape Λ is calculated in Mandel notation The matrices Λ and Λ, both are symmetric and positive-semidefinite for fiber orientation tensors A. The sum of the eigenvalues (respectively the trace) of the interpolated shape Λ is still normalized.However, since the matrices Λ and Λ, are diagonal, their representation in tensor notation does not provide full index symmetry.Thus, Λ and Λ, do not represent fourth-order fiber orientation tensors A in Mandel notation.
To perform the orientation interpolation, the dyadic orientation averaging for orthogonal second-order tensors  ∈ ℎ ⟨3⟩ [25] can be extended to the orthogonal matrices R, in Mandel notation.The only difference to the method for  ∈ ℎ ⟨3⟩ is that all the tensor operations are performed in the six-dimensional space.The singular value decomposition that is used for the dyadic orientation average can also be evaluated in the six-dimensional space.After the dyadic orientation averaging is performed, the resulting interpolated orientation R is given.Finally, the interpolated tensor can be obtained (similarly to Eq. ( 6)) by reassembling the orientation average R and the shape average Λ in Mandel notation The matrix Ā can be transformed to tensor notation ĀSP .Since this method mainly uses the spectral method for the shape interpolation it is called "spectral-method" (SP) in the following.All results using the SP-method are denoted as ĀSP .

Results and discussion
In the following the described SP-method is used to average fourth-order fiber orientation tensors A. To visualize the tensors the fiber alignment is used.The fiber alignment  4 () is plotted as a surface in the three-dimensional space (see Fig. 2).This graphical representation is similar to that of the Youngs's modulus from Böhlke and Brüggemann [48].The coordinate axes  1 and  2 are marked using red arrows.Now a special case to apply the interpolation method is considered.Here only two fiber orientation tensors A 1 and A 2 are averaged.Both tensors are orthotropic.Using Mandel notation, the first supporting point A 1 is described by The visualization of A 1 is shown in Fig. 2(a).The representation is projected into the (  1 ,  2 ) -plane.The principal axes of A 1 coincide with the global basis system { 1 ,  2 ,  3 }.The orthotropic material symmetry is evident in the graphical representation, as both the (  1 ,  3 ) -plane and the (  2 ,  3 ) -plane serve as planes of symmetry.Note that the first principal axis (in  1 -direction) is a local maximum of  4 () while the second principal axis (in  2 -direction) is a local minimum of  4 ().The second supporting point, A 2 (not visualized), is obtained by rotating A 1 by an angle of 45 • around the  3 -axis.Consequently, A 1 and A 2 share the same Kelvin moduli and hence the same shape.They only differ in their relative orientation to the global basis system.
To interpolate the tensors, the SP-method is employed with weights  1 = 0.6 and  2 = 0.4.The resulting interpolated tensor, denoted as ĀSP , shares the same Kelvin moduli as the two supporting points A 1 and A 2 .The tensor ĀSP is positive-semidefinite and the condition tr( ĀSP ) = 1 holds.The visualization of ĀSP is shown in Fig. 2(b).The local maximum of  4 near the  1 -direction is marked as  max  1 .The local minimum of  4 near the  2 -direction is marked as  min 2 and almost coincides with  2 .It is important to note that the two directions of the extrema ( max 1 and  min  2 ) are no longer orthogonal, which results in the absence of symmetry planes in Fig. 2(b).Furthermore, the interpolated fourth-order fiber orientation tensor ĀSP is no longer fully symmetric.Thus, ĀSP is not an admissible fourth-order fiber orientation tensor.
In conclusion, the SP-method is not an appropriate interpolation method because the interpolation process loses the full index symmetry of fiber orientation tensors.The interpolation procedure should result in admissible fourth-order fiber orientation tensors.Furthermore, this indicates that the spectral decomposition is not a suitable method for decomposing orientation-and shape-describing features for fourth-order fiber orientation tensors in three-dimensional space.Therefore, to separately describe the orientation-and shape-defining features of fourth-order fiber orientation tensors in three-dimensional space, alternative decompositions based on the harmonic decomposition (cf.Section 3.2) are pursued.

Framework for the new decomposition-based interpolation 4.2.1. Shape-describing features of orthotropic fourth-order fiber orientation tensors
Krauß and Kärger [9] demonstrated that the decomposition of orientation-and shape-describing features, followed by their separate interpolation, is a suitable technique for interpolating symmetric second-order tensors.The direct application of methods for second-order tensors was rejected (see Section 4.1) because the SP-method did not allow for a complete separation of shape and orientation describing features in three-dimensional space.Thus, this section defines the features that describe the shape and orientation of orthotropic fourth-order fiber orientation tensors.
To describe an arbitrary fourth-order fiber orientation tensor A, 15 independent coefficients are required e.g.[2].The orientation of the considered orthotropic tensor is described using a rotation tensor  ∈ ℎ ⟨3⟩ .To describe the tensor's orientation in threedimensional space, three of the 15 coefficients are needed [49].The remaining 12 coefficients, which describe the tensor's shape, can be expressed using 12 invariants [50,51].If the considered fourth-order fiber orientation tensor A exhibits material symmetries, fewer scalar values are required to fully describe its tensor shape.An orthotropic fourth-order fiber orientation tensor, A Orth , can be displayed in Mandel notation [2] To describe the component matrix in Eq. ( 27) only six scalar values are required.However, the special matrix form in Eq. ( 27) is only achieved if a certain basis system {ẽ 1 , ẽ2 , ẽ3 } is used.This basis system coincides with the principal axes of A Orth [52].Therefore, the basis system of the principal axes {ẽ 1 , ẽ2 , ẽ3 } is shape-defining.Note that if the basis system {ẽ 1 , ẽ2 , ẽ3 } does not coincide with the orthotropic axes of A Orth , the component matrix in Eq. ( 27) can be completely populated.In this case, it may not be possible to describe the component matrix of A Orth with only six scalar values.For fourth-order fiber orientation tensors with weaker symmetry than orthotropy (e.g., arbitrary triclinic fiber orientation tensors), it is not possible to identify an orthonormal basis system {ẽ 1 , ẽ2 , ẽ3 } to construct the component matrix of A as outlined in Eq. (27).It is therefore impossible to describe the shape of an arbitrary triclinic fiber orientation tensors A using only six scalar values within the approach presented in this work.In summary, an orthotropic fourth-order fiber orientation tensor A Orth requires nine coefficients for a distinct description.Three coefficients are required to define the rotation,  ∈ ℎ ⟨3⟩ , from the global basis system to the principal axes of A Orth .The remaining six coefficients are used to describe the shape of A Orth , as shown in Eq. ( 27).

Decomposition of orientation-and shape-describing features
To impose the form in Eq. ( 27), the principal axes of A Orth must coincide with the chosen basis system.Following Cintra and Tucker [52], the principal axes of A Orth align with the eigenvectors of the corresponding second-order fiber orientation tensor  = A Orth [].Therefore, the orientation-describing feature of an orthotropic fourth-order fiber orientation tensor can be calculated by solving the eigenvalue problem for .Rearranging Eq. ( 13) the spectral decomposition for  reads The orthogonal tensor  ∈ ℎ ⟨3⟩ describes the rotation of  to its eigensystem.The resulting fourth-order rotation tensor R □ =  □   transforms the principal axes of A Orth to the global basis system.Consequently, the resulting rotated tensor A Orth, * , that describes the shape of A Orth , is calculated Eq. ( 29) decomposes orientation-and shape-describing features for orthotropic fourth-order fiber orientation tensors A Orth .All shape information of A Orth is contained in the tensor A Orth, * , while the rotation tensor  ∈ ℎ ⟨3⟩ from Eq. ( 28) contains all the orientation information of A Orth in the three-dimensional space.Analogous algorithms can be utilized to find the principal axes of an orthotropic stiffness tensor C Orth [see 45,53].
Since the principal axes can be arbitrarily permutated or multiplied by (−1) the found rotation tensor  ∈ ℎ ⟨3⟩ in Eq. ( 28) is not distinct [53].In these cases the coefficient matrix of the resulting A Orth, * still has its special form from Eq. (27).In total, there are 24 different rotation tensors , that rotate the principal axes of A Orth to the global basis system.To limit this arbitrariness, the following two conditions are applied to  and A Orth, * det() = 1 and Eq. ( 30) postulates a right-handed coordinate system and introduces a convention for the sequence of the principal axes {ẽ 1 , ẽ2 , ẽ3 }.
By minimizing the rotation angle as an intrinsic metric of ℎ ⟨3⟩ , the shape of the tensor A Orth, * is uniquely determined.Note that there are additional cases for non-distinct orientation information caused by material symmetries.If an isotropic fourthorder fiber orientation tensor A Iso is considered, every basis system will align with the principal axes.In case of a transversal isotropic material symmetry, an infinite number of principal basis systems can be obtained by rotating a suitable basis system around the transversal isotropic axis.However, for triclinic or monoclinic material symmetries no orthogonal basis system of principal axes can be found.Therefore, for triclinic or monoclinic fiber orientation tensors the representation from Eq. ( 27) does not exist and no suitable  ∈ ℎ ⟨3⟩ can be found.

Reassembly
In the previous sections the decomposition of orthotropic fourth-order fiber orientation tensors A Orth in shape-and orientationdescribing features was presented.The tensor's shapes are fully defined by A Orth, *  (Eq.( 27)), whereas the tensor's orientations are fully described by   ∈ ℎ ⟨3⟩ (Eq.( 28) and Eq. ( 29)) The features describing the shape and orientation are interpolated separately.
The shape interpolation is accomplished by utilizing the shape-describing tensors A Orth, *  , the weights   , and one of the presented shape interpolation methods  S in Section 4.3 The principal axes of ĀOrth, * and A Orth, *  coincide with the global coordinate axes { 1 ,  2 ,  3 }.To execute the orientation interpolation  O of the rotation tensors   ∈ ℎ ⟨3⟩ , either the dyadic orientation interpolation [24] or the unit quaternion interpolation method [9] is used Subsequently, the resulting orientation average R ∈ ℎ ⟨3⟩ is available.Finally, the resulting interpolated fourth-order fiber orientation tensor ĀOrth can be calculated (compare to Eq. ( 29))

Shape interpolation for fourth-order orientation tensors 4.3.1. Decomposition-based component interpolation
In the following the decomposition-based component interpolation (DBC) is presented.All results addressing the DBC-method are denoted by (⋅) DBC .The orientation of the supporting points A  is extracted in a previous step (refer to Section 4.2.2) and interpolated separately.For all supporting points A  , the resulting shape-describing tensors A *  are calculated.Shape interpolation using the DBC-method does not involve tensor invariants.Instead, the interpolation is carried out by computing a weighted mean of the tensor components of the shape-describing tensor A *  in the global basis system.The interpolated tensor shape ĀDBC, * is calculated by The only difference compared to the classical EU-method (Eq.( 7)) is the decomposition into orientation and shape describing features and their separate interpolation.All tensors A *  are positive-semidefinite, fully symmetric and satisfy the condition Thus, all A *  are admissible fourth-order fiber orientation tensors.The sum of positive-semidefinite tensors stays positive-semidefinite.Similarly, the sum of fully symmetric tensors remains fully symmetric.The condition stated in Eq. ( 35) remains satisfied after summation, Thus, the interpolated tensor shape ĀDBC, * and the interpolated tensor ĀDBC (Eq.( 33)) are both admissible fourth-order fiber orientation tensors.Note that this shape interpolation is similar to the one proposed by Bauer et al. [37] for orthotropic material symmetry.

Deviator invariants interpolation
In the following the shape interpolation using deviator invariants (DI) is presented.All results addressing the DI-method are denoted by (⋅) DI .First, the used invariants are described.The DI-method utilizes the harmonic decomposition for fourth-order fiber orientation tensors A (see Section 3.2) In order to develop an invariant-based interpolation scheme, the orthogonality of the parts A Iso , A Dev 2 and A Dev 4 is utilized.The constant term, A Iso , is always equal for all admissible fourth-order fiber orientation tensors.Therefore, A Iso can be excluded from the interpolation.The part A Dev 2 only depends on .Therefore, A Dev 2 can be fully described using the orthogonal invariants  ⟨2⟩  , which were introduced by Ennis and Kindlmann [22] Using A[] = , the three invariants  ⟨2⟩  can be displayed as functions of A. They are referred to as  DI  ( = 1, 2, 3) The parameters  1 ,  2 and  3 are multiples of the Kelvin moduli of A Dev 4 .Due to this fact, they are invariant with respect to all rotations  ∈ ℎ ⟨3⟩ .Thus, they are suitable invariants for the shape interpolation The principal axes of the shape-describing tensors A Dev 4 , * and the global basis system always coincide, which allows calculation of the invariants  DI 4 ,  DI 5 and  DI 6 using the projections (without summation over  = 4, 5, 6) The shape of the fourth-order deviator A Dev 4 , * is fully described with the invariants  DI 4 ,  DI 5 and  DI 6 To perform the shape interpolation, a weighted mean of the tensor invariants is evaluated The first three interpolated invariants KDI 1 , KDI 2 and KDI 3 describe the shape of A Dev 2 , * (Eq.( 37) and Eq. ( 38)).The remaining three interpolated invariants KDI 4 , KDI 5 and KDI 6 describe the shape of A Dev 4 , * (Eq.( 40) and Eq. ( 43)).The interpolation of the invariants KDI  corresponds to the separate interpolation of the individual parts of the harmonic decomposition Due to the fact that ĀDev 2 , * is only dependent on  * , the interpolation of ĀDev 2 , * can be performed by interpolating  *  .Therefore, the invariant-based interpolation method presented by Krauß and Kärger [9] is used.Since Krauß and Kärger [9] use the same invariants (Eq.( 39)) as the DI-method, Eq. ( 44) is satisfied for KDI 1 , KDI 2 and KDI 3 .Consequently, only KDI 4 , KDI 5 and KDI 6 have to be calculated using Eq. ( 44).The interpolated fourth-order deviator ĀDev 4 , * can be calculated by reinserting d1 , d2 and d3 into Eq.(40).

Orthotropic invariants interpolation
In the following the shape interpolation using orthotropic invariants (OI) is presented.All results addressing the OI-method are denoted by (⋅) OI .First, the used invariants are described.A separate eigenvalue decomposition for the second-and fourth-order fiber orientation tensor  and A is performed.For this purpose, it is used that orthotropic fourth-order fiber orientation tensors A have the form in the basis system of their principal axes.This representation shows that the parameters 2 2233 , 2 1133 and 2 1122 are eigenvalues (respectively Kelvin moduli) of A and thus are invariants with respect to  ∈ ℎ ⟨3⟩ Since the principal axes of A and the global basis system align, the according second-order fiber orientation tensor Hence, the parameters  11 ,  22 and  33 are eigenvalues of  and consequently are invariants of A with respect to  ∈ ℎ ⟨3⟩ The components  1111 ,  2222 and  3333 in Eq. ( 46) can be calculated as a function of the Kelvin moduli  ⟨A⟩  (Eq.( 47)) and the eigenvalues  ⟨⟩  of the according second-order fiber orientation tensor  (Eq.( 49)).If A is given in its principal axes system the following conditions hold [52] These relations are often used to derive closure approximations e.g.[11,54].
The invariants  ⟨2⟩  , which were introduced by Ennis and Kindlmann [22], provide suitable shape interpolation for second-order tensors  [9].In Appendix A it is shown that the gradients of the invariants  ⟨4⟩  ( = 1, 2, 3), depending on A, still provide pairwise orthogonal gradients.Thus, those invariants are used for the shape interpolation of the OI-method [9,22] Since  OI  ( = 1, 2, 3) can be used for shape interpolation of second-order tensors , the invariants from Eq. ( 51) are used to describe the eigenvalues  ⟨⟩ 1 , To perform the shape interpolation, the tensors A *  are used.For those tensors the orthotropic system and the global basis system always coincide.Thus, the invariants  OI  4 ,  OI  5 and  OI 6 can be calculated using the projections (no summation over  = 4, 5, 6) To perform the shape interpolation a weighted mean of the tensor invariants is evaluated Since Krauß and Kärger [9]  are calculated using Eq. ( 54).Finally, the resulting interpolated shape ĀOI, * is given with the results for λ⟨A⟩ The interpolated coefficients Ā1111 , Ā2222 and Ā3333 can be calculated by rearranging Eq. ( 50)

Problem statement
The validation is performed using the numerical solution of a theoretical problem from fluid mechanics.The relevant equations are the evolution equations for the fiber orientation tensors A described by Advani and Tucker [6].In the originally proposed tensor-valued differential equation (DE) for the fourth-order fiber orientation tensors A [6], the right-hand side is not necessarily fully symmetrized.As a result, the solutions of the DE are not necessarily positive-semidefinite.Therefore, Jack and Smith [55]  presented a DE with a fully-symmetrized second term on the right-hand side.The evolution of the fourth-order fiber orientation tensor A depends on itself, the sixth-order fiber orientation tensor A ⟨6⟩ , the velocity field , its gradient  and the material parameters  e and  I .The property  e describes the geometry of the fibers.For elliptical particles,  e is calculated by  e = ( 2 e − 1)∕( 2 e + 1) where the quantity  e is the ratio of the fiber length to the fiber diameter (aspect ratio).The quantity  I is a material property that is a measure for the degree of fiber-fiber interaction in the melt [6].
The solution to the DE is derived for a special case that significantly simplifies the relevant equations.The problem under consideration involves the point sprue of a plate with infinite diameter, allowing the melt to spread radially unhindered (refer to Fig. 3(a)).The problem is assumed to be stationary.Therefore, all considered fields are only dependent on the location  but not on the time
The absence of a velocity component in tangential direction   and thickness direction  z is assumed, i.e. () =  r () r .Consequently, mass transport in the thickness and tangential direction is not considered in this approach.The velocity field  is assumed to be rotationally symmetric () =  r (, ) r .To investigate the influence of a non-constant velocity profile over the plate thickness , the velocity field in radial direction is approached as a quadratic function of the thickness coordinate  (see Fig. 3(b)).
The mid-plane of the cavity is located at  = 0, while the wall of the cavity is located at  =  0 .Using the symmetry condition  r  (,  = 0) = 0 and the boundary conditions  r (,  =  0 ) =  min () and  r (,  = 0) =  max (), the -dependency of  is defined (see Fig. 3(b)).Furthermore, if the initial values  max ( =  0 ) =  max,0 and  min ( =  0 ) =  min,0 are specified as well, the velocity field is completely described by Here,  describes the radial distance to the sprue and  0 the radial distance where the initial velocity value  0 () is specified.The values  max,0 ,  min,0 , ,  0 ,  and  0 are always positive and the conditions  ≤  0 and  ≥  0 hold.This particular case closely resembles the one considered by Krauß and Kärger [9] for second-order fiber orientation tensors .Therefore, some parts of the solution can be derived analogously.However, it addresses the DE for a fourth-order fiber orientation tensor A and considers a non-constant velocity profile  across the plate thickness .The solution to the problem is developed in cylindrical coordinates.
To evaluate the DE for the fiber orientation tensor A, the velocity gradient  is needed.It is calculated by Its symmetric part  and its skew-symmetric part  can be computed using Fig. 4. Courses of the tensor components of A ana in dependency on the normalized radius ∕ 0 .The dimensionless parameters  min,0 ∕ max,0 = 1 and  e = 0.3 as well as the initial value A ana 0 = A Iso are chosen.The solution is not dependent on ∕ 0 .
Note that if the velocity function  r (, ) is constant over the plate thickness (e.g.,  min,0 =  max,0 ), the velocity gradient is purely symmetric ( = ) and the skew-symmetric part vanishes ( = ).This case corresponds to a plug-flow.
Since the considered problem is stated to be stationary, the local change of A on the left-hand side of the DE vanishes [55].Assuming no interaction between the fibers ( I = 0), the DE is finally simplified to Note that the operator sym (⋅) applies the index symmetry for all indices.To approximate the sixth-order fiber orientation tensor the quadratic closure is utilized A ⟨6⟩ = A ⊗  [6].After evaluating Eq. (61) in cylindrical coordinates, a coupled system of 15 scalar DEs is derived.Since the velocity vector (, ) is directed in radial direction  r (see Eq. ( 58)), only partial derivatives with respect to the radius  can be found on the left-hand side of Eq. (61).As a result, all considered DEs are in fact ordinary DEs.Therefore, the thickness coordinate  can be treated as a parameter that can be specified in advance, and thus the problem can be treated as planar for any given thickness coordinate .The solutions for the tensor components of A in the cylindrical basis system { r ,   ,  z } were computed numerically using an explicit fourth-order Runge-Kutta method, which requires a suitable initial value A 0 for the fiber orientation tensor A at the initial radius  0 .
Note that the solutions for the tensor components of A can be calculated analytically if the velocity profile (, ) is uniform over the plate thickness (e.g.,  min,0 =  max,0 ).The resulting explicit equations to compute the tensor components were obtained using a computer algebra system.The formulas are provided in Appendix B.

Results of the thickness-independent problem
First, the special case  max,0 =  min,0 (plug flow) is considered.The velocity gradient  is purely symmetric ( = ) and its skew-symmetric part vanishes  = .Therefore, the analytical functions given in Appendix B can be used to calculate the resulting tensor components.The trajectories of the resulting non-zero components of A ana over the normalized radius ∕ 0 are depicted in Fig. 4. The solution is calculated for  e = 0.3 and the initial fiber orientation A ana 0 = A Iso .Note that in this special case, the velocity vector  is not dependent on the thickness  anymore but only on the radius .Therefore, the normalized thickness ∕ 0 can be chosen arbitrarily (see Eq. ( 58)).
Despite the problem being planar, all non-zero components of A ana 0 change with increasing normalized radii ∕ 0 .Independent of the choice of the initial value A ana 0 and the material parameter  e , the magnitude of  ana rrrr ,  ana rr ,  ana rrzz and  ana zzzz decay with increasing normalized radius ∕ 0 .The value of  ana zz initially increases and then decreases at a slower rate than the other tensor components.Only the  ana  component converges towards the value of one.As the normalized radius ∕ 0 increases, the orientation state becomes more anisotropic.At the limit  → ∞, the solution of the continuous point sprue problem is unidirectional (lim →∞ A =   ⊗   ⊗   ⊗   ) for all admissible initial values.The analogous observation was made by Krauß and Kärger [9] for second-order fiber orientation tensors (lim →∞  =   ⊗   ).
All solutions A ana are fully symmetric, positive semi-definite, and the trace is normalized.Therefore, all solutions A ana are admissible fourth-order fiber orientation tensors.However, the material symmetry of solutions differ depending on the initial fiber orientation tensor A ana 0 .If the initial value A ana 0 is chosen isotropic (as in Fig. 4) an orthotropic fourth-order fiber orientation tensor Fig. 5. Courses of the tensor components of A num in dependency on the normalized radius ∕ 0 .The dimensionless parameters ∕ 0 = 0.9,  min,0 ∕ max,0 = 0.3 and  e = 0.3 as well as the initial value A 0 = A Iso are chosen here.
is obtained for ∕ 0 > 1, whose principal axes correspond to the coordinate axes { r ,   ,  z }.If the initial value A ana 0 is chosen orthotropic, the solution for ∕ 0 > 1 remains orthotropic if the principal axes of the initial value A ana 0 align with the coordinate axes { r ,   ,  z }.If an orthotropic initial value A ana 0 is chosen whose principal axes do not align with the coordinate axes { r ,   ,  z }, a triclinic fiber orientation tensor is obtained for ∕ 0 > 1.If the initial value A ana 0 is chosen triclinic, the resulting tensors A ana are triclinic for any normalized radius ∕ 0 ≥ 1, except in the limit  → ∞ where the solution becomes unidirectional.

Results of the thickness-dependent problem
Now, the condition  max,0 =  min,0 is dismissed, and the radial velocity component  r from Eq. ( 58) is indeed described with a quadratic function.Consequently, the resulting velocity gradient is not purely symmetric anymore ( ≠  for all ∕ 0 ≠ 0).The parameters are chosen ∕ 0 = 0.9,  min,0 ∕ max,0 = 0.3 and  e = 0.3.The initial fiber orientation tensor is again considered isotropic A num 0 = A Iso .The trajectories of the non-zero tensor components for the now considered case are displayed in Fig. 5. Unlike the previous example (Section 5.2.1), the courses of the tensor components are oscillating.The amplitude of the oscillations increases if ∕ 0 or  min,0 ∕ max,0 is raised, or if  e is decreased.Increasing ∕ 0 or  min,0 ∕ max,0 causes an increasing skew-symmetric part of the velocity gradient  .Additionally, a decreasing  e decreases the influence of the symmetric part  on the right-hand side of the DE (Eq.( 61)).Since no oscillations were observed in the previous example (Section 5.2.1), it can be concluded that the oscillations are directly caused by the skew-symmetric part of the velocity gradient  .Similarly to the previous example, the magnitude of each tensor component decays with increasing normalized radius ∕ 0 .Again, the only exception to this is  num  .This behavior is observed for arbitrary initial values A num 0 and arbitrary parameter sets.Note that the components  num rrrrz ,  num rz and  num zzrz are non-zero in this particular case, which was not observed in the thickness-independent case from Section 5.2.1 (refer to Fig. 4).
Since the solutions are always fully symmetric, positive-semidefinite, and the trace is normalized, all solutions A num are admissible fourth-order fiber orientation tensors.However, the material symmetry of solutions differ depending on the initial fiber orientation tensor A num 0 .If the initial value A num 0 is chosen isotropic (as in Fig. 4 and Fig. 5), an orthotropic fourth-order fiber orientation tensor is obtained for ∕ 0 > 1.

Procedure
To assess the performance of the different interpolation methods, they are applied to an associated discrete problem.The continuous point sprue problem from Section 5.1 is employed to generate data for the discrete problem.Note that the developed interpolation methods from Section 4 are designed for orthotropic fiber orientation tensors; therefore, the isotropic fiber orientation tensor is used as an initial value, A num 0 = A Iso .The schematic visualization of the associated discrete problem is presented in Fig. 6.The hexahedral cell considered, as shown in Fig. 6(b), has a square base with an edge length of √ 2  and a height of 2  .The aspect ratio of the cell is denoted as  C =   ∕( √ 2  ).The location of the middle of the cell can be described using (62) Here,  M and  M describe the position of  M relative to the origin  in the ( 1 ,  2 )-plane using cylindrical coordinates (Fig. 6(a)).
The coordinate  M describes the according component in thickness ( 3 ) direction (Fig. 6(b)).The locations of the supporting points   with  = 1..8 can be calculated analogously.Since the numerical solutions of the continuous problem only depend on the radii   and the thickness coordinates  I , the fiber orientation tensors in the vertices A  = A num (  ) can be computed.Similarly, the fiber orientation tensor in the middle of the cell can be computed A M = A num ( M ).In the following, the calculated fiber orientation tensors in the middle A M and the vertices A  are considered exact.Therefore, they can be used to evaluate the performance of the shape interpolation methods  S from Section 4. The fiber orientation tensors A  at the eight vertices of the cell are used as support points.The fiber orientation tensor A M , which is located in the center of the cell, is used as a comparison value for the interpolation Ā.Since the vertices   each have the same distance to the center  M , the weights at the supporting points can be chosen to be equal   = 0.125 for  = 1..8.The fiber orientation tensors calculated at the supporting points A  are given in the basis system { r ,   ,  z }.In order to perform the interpolation of the fiber orientation tensors, they must be given in the same basis system.Therefore, the fiber orientation tensors A  are transformed into the basis system { 1 ,  2 ,  3 }.The angle  M is chosen  M = ∕6.

Error metrics
In order to compare the interpolation results Ā with the numerical solution, scalar evaluation metrics are required.In total three rotation invariant scalars  (A) are used The analogous quantities for second-order tensors  have already been used by Krauß and Kärger [9] to evaluate their interpolation results.Therefore, ‖A‖ and FA 4 (A) will be used again for fourth-order fiber orientation tensors A. The invariant  ⟨4⟩ 2 (A) is used because it is part of the orthogonal invariants [22] that were very suitable for the shape interpolation second-order tensors .The fiber orientation tensors in the cell center A M , which are calculated by means of the numerical solution, are used as the reference solution.In each case, the relative errors  rel are considered 5.3.3.Results and discussion for fixed relative mesh sizes Fig. 7 shows the relative errors of the different interpolation methods for the problem described in Section 5.3.1.The closure-based interpolation method (CL) is used as a comparison interpolation method.The CL-method interpolates the according second-order fiber orientation tensors   by using the methods from Krauß and Kärger [9].Afterwards, the resulting interpolated tensor ĀCL is mapped to ĀCL using the IBOF-closure [56].To perform the IBOF-closure, fiberoripy [57] is used.
All evaluated error measures are defined in Section 5.3.2.The considered parameter set is given by  As an initial value A num 0 , the isotropic fiber orientation tensor A Iso is used since it guarantees orthotropic solutions A num for all radii .The trajectories of the relative errors in Fig. 7 partially overlap.In order to clearly distinguish the progressions of the different interpolation methods, the trajectories are marked by colored dots.In the limit  M ∕ 0 → ∞ (not shown in Fig. 7) the deviations of all interpolation methods disappear.
In general, Fig. 7 shows that the EU-method (red, Eq. ( 7)) produces the worst interpolation results for all considered errors.Similar to the DBC-method (purple, Section 4.3.1), the EU-method averages the tensor components in a Euclidean manner.However, the DBC-method rotates the tensors to their shape-defining coordinate system to realize the decomposition into orientation-and shape-defining features.This indicates that an interpolation method, where the tensor shape is separated from the tensor orientation, generally provides better results than Euclidean averaging of the tensor components in the global basis system.
The relative errors of the DBC-method, the DI-method (blue, Section 4.3.2) and the OI-method (turquoise, Section 4.3.3)hardly differ from each other.To interpolate the tensor shape using the OI-or the DI-methods, tensor invariants are calculated and interpolated.Hence, these interpolation methods are more complex than the DBC-method.The DBC-method does not require the calculation of tensor invariants (and an associated return path).Nevertheless, the DBC-method yields results that are not inferior to the results of the other decomposition-based interpolation methods.
The deviation of the results from the CL-method (green) is significantly smaller than that from the EU-method (red).Nevertheless, they clearly stand out from the deviations of the other interpolation methods.In particular, the deviations in the errors ‖ ⋅ ‖ and FA 4 (⋅) for large normalized radii ( M ∕ 0 > 3) should be emphasized here.The reason for these characteristic deviations in the range  M ∕ 0 > 3 is that the IBOF-closure cannot map the averaged second-order fiber orientation tensor ĀCL to A M .For this reason, the deviation of the CL-method appears like a systematic error in the interpolation.For large normalized radii  M ∕ 0 → ∞, the numerical solution A num is a unidirectional fiber orientation tensor, which can be obtained by the used IBOF-closure [56].Because of this, the deviation of the CL-method vanishes in the limit  M ∕ 0 → ∞.For the error measure  ⟨4⟩ 2 , the curve of the CL-method (green) lies on the curves of the decomposition-based interpolation methods.The reason for this is that the second-order fiber orientation tensors , which are interpolated in the CL-method, rely on the tensor invariants  ⟨4⟩  1 ,  ⟨4⟩ 2 and  ⟨4⟩ 3 in the same way as the OI-method (turquoise) and the DI-method (blue).

Results and discussion for variable relative mesh sizes
In Fig. 8, the relative deviations of the error measures defined in Section 5.3.2 are plotted for different normalized in-plane mesh sizes   ∕ 0 .Fig. 8 can be used to examine how a change of the in-plane mesh size   ∕ 0 affects the error measures considered.For the interpolation, the EU-method (red), the CL-method (green) and the shape interpolation using DI (blue) are used.The considered parameter set is Note that the aspect ratio of the cell  C is considered constant.Thus, the cell becomes thicker with an increasing in-plane mesh size   ∕ 0 .Fig. 8 reveals that the EU-method exhibits the largest errors for small normalized radii ( M ∕ 0 < 4).In contrast, the relative errors of the DI-method, for the same normalized in-plane mesh size   ∕ 0 , consistently remain smaller than those of the EU-method.For the error measure  ⟨4⟩  2 , the relative errors of the CL-method align with those of the DI-method.This alignment is attributed to both the DI-method and the CL-method, relying on the same orthogonal invariants  ⟨2⟩  for shape interpolation.For large normalized radii ( M ∕ 0 > 4) and small normalized in-plane mesh sizes (  ∕ 0 < 0.75), the relative errors for the norm ‖ ⋅ ‖ and the relative anisotropy FA 4 (⋅) of the EU-method are smaller than those of the CL-method.The reason for this is that the IBOF-closure cannot map the interpolated second-order fiber orientation tensor ĀCL to the numerical solution A M .For large normalized radii ( M ∕ 0 > 4) and small normalized in-plane mesh sizes (  ∕ 0 < 0.75), the error caused by the closure is larger than that caused by the interpolation.For small normalized radii ( M ∕ 0 < 4) and large normalized mesh sizes   ∕ 0 > 0.75, the CL-method (green) consistently outperforms the EU-method (red).The relative deviations of the EU-method and the DI-method become smaller with decreasing normalized in-plane mesh sizes   ∕ 0 .Considering FA 4 (⋅), the EU-method gives similar good results with   ∕ 0 = 0.5 as the DI-method with   ∕ 0 = 1.625.

Considering 𝐾 ⟨4⟩
2 , the EU-method gives similar good results for   ∕ 0 = 0.5 as the DI-method with   ∕ 0 = 1.25.This shows that the DI-method is able to give similar good interpolation results as the EU-method for much coarser meshes.In this specific example, the in-plane edge length of the hexahedral cells in Fig. 6(b) can be increased to about two to three times that of the EU-method.This allows the DI-method to achieve similar interpolation results without a significant increase in relative deviation.
Considering the error measure FA 4 (⋅) for large normalized radii ( M ∕ 0 > 3), the relative deviations of the CL-method are smaller for large normalized mesh sizes   ∕ 0 than for small normalized mesh sizes.This unexpected behavior is attributed to the employed closure approach.For small normalized radii ( M ∕ 0 < 4), the interpolated second-order fiber orientation tensor Ā does not map to A M due to the IBOF-closure.Therefore, for smaller normalized mesh sizes   ∕ 0 , the relative error of FA 4 (⋅) approaches the systematic error caused by closure.This error is problem-specific and may not solely reflect the quality of the interpolation method itself.

Summary of interpolation results
Overall, the EU-method consistently yields the poorest interpolation results across all scenarios.As the spatial discretization   ∕ 0 is refined, the errors decrease rapidly, and the interpolation results converge towards the exact solution.All results obtained using the EU-method are admissible fourth-order fiber orientation tensors.
The CL-method outperforms the EU-method significantly for small normalized radii ( M ∕ 0 < 3) and large normalized mesh sizes (  ∕ 0 > 0.75), but not always for small normalized mesh sizes (  ∕ 0 < 0.75).The slower convergence of the CL-method to the exact solution is attributed to the systematic error introduced by the closure.However, exclusively for small normalized mesh sizes and strongly aligned grid points (here  M ∕ 0 ≫ 5), the CL-method is inferior to the EU-method.Despite this, the CL-method yields admissible fourth-order fiber orientation tensors, and its computational effort is comparatively small, leveraging well-established methods from literature [9,57].In conclusion, the CL-method serves as a viable alternative to the EU-method.
The DBC-method consistently provides excellent interpolation results, yielding admissible fourth-order fiber orientation tensors.It is a comparatively simple interpolation method since it does not require the calculation of any tensor invariants.Its simplicity and superior performance, compared to the EU-method, highlight the effectiveness of separating tensor shape and orientation in spatial interpolation of fourth-order tensors.In conclusion, the DBC-method can be evaluated as a very good and efficient alternative to the EU-method.
The DI-method and OI-method are evaluated together for all errors since their relative deviations coincide.In each case, six tensor invariants are used for interpolation.Both interpolation methods utilize the three orthogonal invariants  ⟨2⟩  .Both methods provide excellent interpolation results.For the DI-and the OI-method, tensor invariants have to be calculated.Despite their algorithmic complexity, requiring the calculation of tensor invariants and the return path to the shape of the averaged fourth-order fiber orientation tensor, the DI-and OI-methods did not outperform the simpler DBC-method in the considered example.Nonetheless, the DI-and OI-methods are viable alternatives to the EU-method.
In Section 5.3.4 the influence of the normalized mesh size   ∕ 0 on the interpolation result of the DI-, the CL-and the EU-method were investigated.It was shown that the EU-method can only reach equally decent results (compared to the DI-method) with much finer meshes.

Conclusion and outlook
This work presents the development of interpolation methods for the decomposition-based interpolation of fourth-order fiber orientation tensors.Section 2 provides an overview of the related works.A literature review on spatial interpolation of fourthorder tensors indicates that the topic has not yet been considered in a mechanical context.Currently, the state of the art is limited to interpolation methods for averaging second-order tensors.
Section 4 describes the extraction of orientation-and shape-describing features from orthotropic fourth-order fiber orientation tensors.In contrast to second-order fiber orientation tensors, spectral decomposition did not provide an adequate decomposition of orientation-and shape-describing features for fourth-order fiber orientation tensors.However, known methods can be utilized for the interpolation of orientation-describing features.For the interpolation of shape-describing features, different approaches with varying complexity are described.Describing the shape of second-order tensors benefits from using tensor invariants whose gradients are orthogonal.It was shown that the gradients of the tensor invariants used for second-order tensors are still orthogonal when described by fourth-order tensors.Therefore, they serve as a proper choice to describe the shape of fourth-order orientation tensors.Section 5 derives the numerical solution for a continuous problem from fluid mechanics.This solution is utilized to generate reference data for the associated discrete problem.Based on the discrete problem, the interpolation methods developed in Section 4 are applied and evaluated.The results show a significant improvement in interpolation accuracy compared to conventional methods.
The new interpolation methods developed for the mechanical context of the CAE chain can significantly reduce information loss between simulation steps.This results in overall improved simulation results for the holistic component design process.Additionally, in design processes where a conventional mapping procedure requires a very fine mesh discretization, a coarser simulation mesh can be selected in upstream and downstream simulation steps.This can save computing time while maintaining the quality of the simulation results.
Within the scope of this work, only fiber orientation tensors with orthotropic material symmetry are considered.This is due to the fact that most fiber orientation tensors related to the virtual manufacturing process are at least orthotropic.For measured fourthorder fiber orientation tensors from CT-scans, which are triclinic rather than orthotropic, the developed interpolation techniques cannot be applied.This is due to the fact that it is impossible to find an orthonormal basis system that can distinctly describe the triclinic tensor's shape as outlined in this work.Therefore, future works will focus on the development of interpolation schemes for triclinic fourth-order fiber orientation tensors.To prove the condition, the identity  = A[] is used to apply the chain rule to the derivatives in Eq. (68).In index notation the condition is derived as follows ( = 1, 2, 3) The gradients of  ⟨2⟩  are calculated by Ennis and Kindlmann [22].They satisfy the condition

CRediT authorship contribution statement
The orthogonality condition in Eq. ( 68) can be evaluated by using the result from Eq.

Appendix B. Analytical solution for the point sprue problem
The following section presents the equations that analytically solve the continuous problem discussed in Section 5.1.Note that the equations are derived for the special case where the velocity vector  is independent of the plate thickness (e.g.,  max,0 =  min,0 ).Therefore, the analytical method is not valid for all cases.
In cylindrical coordinates the components are (in Mandel notation)

𝛼(𝑟) ,
where the components  ana 0 denote the components of the initial fiber orientation tensor A ana 0 .

Fig. 1 .
Fig. 1.Mapping of simulation data from source mesh to target mesh.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) The scalars  ⟨⟩  are called eigenvalues and are the roots of the characteristic polynomial () = det( −  ⟨⟩  ).The eigenvalues  ⟨⟩  of symmetric second-order tensors are always real numbers.If all eigenvalues  ⟨⟩  are unique, the associated normalized eigenvectors   are distinct and form an orthogonal basis.

Fig. 2 .
Fig. 2. Visualization of the loss of symmetry when the SP-method is used to average two fourth-order fiber orientation tensors that share the tensor shape.The orientations of the supporting points differ by 45 • around the  3 -axis.The weights are chosen to be  1 = 0.6 and  2 = 0.4.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Problem statement of the continuous point sprue problem.

Fig. 6 .
Fig. 6.Discretization of the equivalent point sprue problem in the plane and the thickness direction.

Fig. 7 .
Fig. 7. Relative errors for an isotropic initial value A num 0 and the parameterset from Eq. (65).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .
Fig. 8. Relative errors for an isotropic initial value A num 0 and variable mesh sizes   ∕ 0 ∈ [0.5..2].(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
[2]DI 3 are equal.Although the invariants  DI 1 ,  DI 2 and  DI 3 depend on A, their gradients are still pairwise orthogonal, as shown in Appendix A.The fourth-order deviator A Dev 4 is independent of .To describe A Dev 4 unequivocally, three invariants ( DI 4 , DI  5, DI  6) are required.IfA is given in the eigensystem of the according , the component matrix of A Dev 4 has the form (in Mandel notation)[2] Eq. (69) shows that the gradients of  ⟨4⟩  can be displayed in dependency of .The condition is reformulated in symbolic notation ⎣  z +  z ⊗   ), ⊗  z +  z ⊗  r ),  3 =  z ⊗  z ,To ensure a clear notation, the following abbreviation () for the denominator is defined as follows() = ( −2 0 + 1 − 2 0 −  0 −  0 )  0 4 e + 2 2 e  0 2 e  0 + 2 −2 e  0 6 e  0 +  −4 e  0 8 e  0 +  4 e  0 .(74)Finally, the components of the fourth-order fiber orientation tensor that solve the problem in Section 5.1 (for velocity vectors  that are independent of the plate thickness) can be calculated by r ⊗   +   ⊗  r ).