A Set of Novel Procedures for Carbon Fiber Reinforcement on Complex Curved Surfaces Using Multi Axis Additive Manufacturing

There has been considerable research in recent years on the additive manufacturing (AM) of carbon fiber reinforced polymer (CFRP) parts based on the process of fused deposition modeling (FDM). The currently-applied steps within the manufacturing pipeline, such as slicing and path planning, consider only the planar case of filament deposition and mostly make no use of the possibility to place single pre-impregnated (prepreg) filaments. Classical methods such as tapelaying and laminating struggle with highly curved and complex geometries and require the costly production of molds, whereas when using AM, these geometries can be realized more easily and molds can be created using the same process. In this paper, a set of algorithms is presented that aims to resolve these problems. Criteria are formulated which enable the goal oriented development and evaluation of the presented methods and represent metrics for future methods. The developed algorithms enable the use of both continuous and discontinuous fiber patches in a much wider range of applications in designing and manufacturing of CFRPs. This opens up new possibilities in this promising field. The developed metrics and infrastructure further constitute progress in the field of multi-axis non-planar path planning for slicing algorithms in general and the conducted evaluation proves the formal applicability of the developed algorithms.


Introduction
Fused deposition modeling (FDM) as a form of 3D printing is the process of constructing a part by extruding a material through a nozzle, fusing extruded filament to the existing structure. Carbon fiber reinforced polymers (CFRPs) are composite materials achieving high strength-to-weight ratios and have seen increased application in the industry in domains such as aerospace, automotive and many others [1].
The process of being able to combine the benefits of FDM and CFRP has been a highly active topic in research over the last few years [1,2]. Dickson et al. state that the reason for this is the potential for CFRPs to match or exceed the mechanical performance of conventional composites [3]. With the development of a pre-impregnated (prepreg) fiber filament by the Markforged company and others, the possibility for printing carbon fiber was further enabled. A lot of research has gone into the evaluation of material properties [3,4], optimization of the printing process [5], printing hardware [2] and some investigations have been made into waste recycling [6,7]. Currently, a major bottleneck is the fact that classical 3D-Slicers lack adapted path planning methods [8,9]. Some publications have proposed planar strategies for utilizing carbon fibers effectively and industrial software exists only for some specific planar application scenarios [1].
Most of the methods currently used for manufacturing carbon fiber parts use molds where either sheets of carbon fiber cloth are layered or tape is placed onto the mold by use

Prior Work
The active field of research in CFRPs and FDM has resulted in a plethora of publications regarding this topic which are presented in the following. Prior to the availability of a continuous prepreg filament short fiber filaments were widely used. The history of short fiber FDM is described in detail by Such et al. in [12]. This process has given way to the use of continuous fibers in recent years as the research in fiber pre-impregnation has progressed. This change is due to continuous fibers resulting in parts with far greater load capacity [3]. In the area of printing and impregnating fibers, several methods for both ex situ and in situ impregnation have been proposed. Ex situ prepreg filament exhibits greater material quality with the drawback of having to be manufactured prior to printing [2]. Furthermore, the equipment and morphological properties have been analyzed and elaboration can be found in [5]. Process and parameter optimization has been investigated with a special focus on the fiber volume fraction and porosity [3].
Zhuo et al. point out that developed software needs to incorporate a stress analysis for CFRPs [9]. There have been multiple publications in which this is achieved in some form, although all were limited to the planar case. Chen et al., as well as Li et al., chose a field-based method, showing great potential [13,14]. Yao et al. developed a method for multi-degree-of-freedom path planning incorporating a stress analysis. However, the possible geometries are limited for this algorithm as it requires flattening of the part resulting in the inability to print parts with non-zero gaussian curvature, meaning local curvature in more than one direction [15]. Multi-axis printing also enables the printing process to function without supporting structures completely. Li et al. showed this by using a geodesic distance field-based approach [16].
Path planning has been an active field of research for the automated fiber placement process for a long time. Thus, methods for complex geometries have been established over the years with approximate algorithms for path planning on curved surfaces being developed as early as 2004 by Schueler et al. [17]. For closed surfaces, Han et al. proposed a solution taking the minimum turning radius into consideration [18]. Newer methods for robotic fiber placement include the application to multi-axis revolution surfaces developed by Hely et al. [19] and Xiao et al. using multi-guidelines on meshes to increase efficiency and applicability for complex geometries [20]. An extensive description and analysis of the state of the art of path planning methods used in automated fiber placement is performed in [21], where the methods using geodesic paths, fix angle paths, steering paths, seed curve projection, independent curves and parallel curves are presented. They have not been applied to the problem of 3D printing patches. Unwanted buckling and issues with highly curved surfaces remain a problem.
Another approach to planning paths for CFRPs is to apply topology optimization in order to simplify the load transmission paths. This was performed by Li et al. in [22] using scaffolding paths and Wang et al. in [23] using the medial axis transform as orientation. Both implemented path planning methods which aim to align the fibers with the principal directions of stress. This was applied to the use of short fibers by Kim et al. in [24]. The concept of using FDM-based molds for the manufacture of CFRPs was proposed in [25] eliminating the need for a separately built mold, reducing the cost and time needed for the production of CFRP parts. The first investigations into the multi-axis manufacture of CFRPs with FDM were made by Zhang et al. in [8] where the verification of the process proved its potential. Kallai et al. analyzed the multi-axis approach with a 12-axis machine in [11] identifying the necessities for enabling the flexible and dynamic robotic manufacturing process.
It can be concluded that the presented approaches establish the possibilities in the promising field of manufacturing CFRPs with FDM, but the lack of software and limitations of existing approaches result in a need for the development of path planning algorithms applied to the manufacturing process.

Our Method
In order to construct paths for reinforcing additively manufactured parts, evaluation criteria motivated by the process and material properties of carbon fiber prepreg filaments are first established. These serve as guidelines in the development of the algorithms and are employed as a way of measuring their performance.
A framework has been developed to embed the algorithms and enable further work: • The proposed procedure starts by analyzing the part to be manufactured in an external finite element method (FEM) software to extract the major principle stress vectors. An exemplary result can be found in Figure 1a. • The resulting information enables the partitioning of the surface to be reinforced into subsurfaces with the goal of achieving a high compliance of local fiber direction and directions of major principle stress. Figure 1b illustrates such a partitioning. This is similar to the zoning approach, used in the layup process. • The subsurfaces are triangulated with high precision to ensure minimal computation error. • For every subsurface a suitable algorithm is chosen. This is performed with the intent of maximizing the measured values of the evaluation criteria and relies on knowledge of the procedures. • The chosen algorithm is executed for every surface resulting in paths that include possible travel motions, an example of which can be seen in Figure 1c. • Post-processing of the paths yields a Gcode to be executed on the multi-axis printer. The focus of this paper will be the developed evaluation criteria and algorithms. The partitioning and the load cases for the finite element analysis (FEA) are subject to further optimization and will only be considered in an exemplary form. Some of the algorithms involve user input for choosing the most fitting fiber direction. This is motivated by the inability of current algorithms to deal with the high amount of divergence and turbulence observed in FEA results. This leads to the approaches discussed above that either only consider small areas of the entire dataset or the need to apply topology optimization, which can negatively effect the spectrum of possible applications. With user involvement at this step higher quality and consideration of the available data can be ensured, although automating this process may be considered for further research.

Methods and Evaluation Criteria
In this section, the process of extruding prepreg filament to form patches is analyzed to formulate the evaluation criteria that constitute the main analytical method. The manufacturing process begins by first printing the core volume with classical methods. On this core, the carbon fiber patches are applied increasing the overall structural integrity and interlaminar performance. To increase this effect the algorithms are designed to be able to align the filament to the main load transmission paths. This approach has been applied in [23,24]. Specifically, parts with open surfaces that cannot be manufactured by filament winding lend themselves to this approach.
In [23], it is stated that disordered stress distribution under different load conditions lead to the inability to follow stress distributions with a continuous path. They apply topological optimization which was not considered as an option in this work as explained above. This enables a higher quality surface finish with total coverage and eliminates the staircase effect often observed with planar FDM. Further elaboration on the process conditions and equipment such as the printing head and robotic setup for which this system was developed can be found in [11], although the algorithms presented in this paper are designed to work out of the box with any multi-axis system.
The software framework and the algorithms were developed in python with the help of the software library pyvista [26], a wrapper for vtk [27]. The user interface uses Qt5 [28] with pyqt [29] and numpy [30] was used extensively throughout. Distance computation was conducted by use of potpourri3d [31], networkx [32] was used for graph computations, scipy [33] for the savgol filter and trimesh [34] for subdiving, remeshing and intersecting the surface. The computer code for the algorithms will be made available upon request to the corresponding author.

Evaluation Criteria
The criteria on the path and properties of the surfaces to be reinforced are key for the goal-oriented formulation and development of the methods. They are chosen to minimize defects and allow the manufactured part to be of high quality. A lot of possible criteria such as isotropicity and manufacturability where addressed in the process of developing the algorithms, but the formal evaluation in this work will be limited to the four criteria presented in the following. As some of them are inversely correlated, trade-offs are discussed as well. The criteria are further illustrated by comparing them to similar aspects found in classical methods. Further elaboration on how the measured values for the criteria are obtained is included in Section 4.

Continuity
In general, discontinuities in the fiber are to be avoided as they lead to decreased tensile stress capacity. Here, the classical process of wrapping continuous fibers around a mandrel achieves the best results, while small patches of fabric for the layup process include a high number of discontinuities on the edges of the cut fabric. Additionally, the process of cutting the fiber in the printing head can lead to problems with the extrusion and could cause delays. This criterion is thus to be minimized to achieve a high efficiency and part quality.

Curvature
A major difference of the FDM process to classical methods is the high amount of local curvature in the path, as neither tape laying nor the layup process entail the placing of single bundles of fiber on open surfaces. It is still an active topic of research what the effect of high curvatures is on the mechanical properties of the resulting part, especially, for the case of fast changes in direction. Matsuzaki et al. identified the setting radius together with the fiber bundle size as the main parameters for printing quality [35]. Lower radii that are associated with high curvature lead to the inability to reach acceptable printing results. This motivates the assumption that high curvatures cause the brittle fiber bundles to fracture or twist, deteriorating the mechanical strength and lowering the printing accuracy, thus, this criterion is to be minimized.

Fiber Area Fraction
A common metric for the quality of fiber reinforced composites is the volume fraction. This measures the ratio of carbon fiber to the volume of a part. There are several methods of determining this important value and obtaining it is still an active topic of research, especially for parts manufactured with FDM. The presence of gaps or unfilled areas can serve as starting points for fractures and decrease the strength of the composite material. In this work, an equivalent measure is applied, but in a conceptual sense, meaning that the resulting index denotes the theoretical best value achieved by the path. The denomination as area fractions instead of volume fractions stems from the use of surfaces instead of volume bodies.

Load Adequacy
To enable a guided approach and fully utilize the possibility to regulate the direction and position of every fiber bundle, the load adequacy is taken as a criterion. To obtain an optimal direction of fibers the approach chosen in [23] is followed: the principle direction vectors are calculated from the FEA. As the fibers perform best under tensile stress, the direction of maximum principle stress is extracted and checked for alignment with the path in every point. This is sometimes called "principle of reinforcement" and the resulting fraction is to be maximized for best performance.

Algorithms
The algorithms described in the following aim to yield printing paths for reinforcement patches on additively manufactured structures. Partly adapted to fit this application from approaches in tape laying, classical 3D-slicing and CNC-machining, they are intended to serve as a starting point for non-planar path planning in manufacturing CFRP parts. In contrast to other approaches such as field-based methods, which struggle with divergence and turbulence, they are conceptualized to draw more from the geometry of the surface itself enabling a higher degree of control.
They are conceptually divided by their directionality, meaning the differentiation is performed in a classical way between direction parallel and contour parallel approaches. Generally, when working on non-planar manifold surfaces, the concept of direction in this sense is ambiguous. This is why several conceptually different methods for computing direction parallel paths are proposed. For the contour parallel approaches, the surface boundary is chosen as the outer contour, with the Geo-Spiral being an exception. The following algorithms are proposed: Direction parallel: • Fixed intrinsic angle (FInA); • Fixed extrinsic angle (FExA); • Seeded ortho-isogeodesics (Ortheo).
One essential concept for most of the algorithms is the notion of geodesic distance. It is computed directly on the mesh using the heat method proposed in [36]. The surface is referred to as T as it is a triangulated mesh and the distance function will be denoted as φ. For distance computation, the eikonal equation |∇φ| = 1 with boundary conditions φ| γ = 0 over some subset of the domain γ (point or curve) is stated. This expression defines the distance function as giving the distance along the domain, being zero at the boundary γ. Further elaboration can be found in [36].

Direction Parallel Algorithms
When constructing direction parallel paths, the notion of angle is of importance. In a lot of cases the tensile strength of fibers is aimed to be maximally utilized by laying the fibers along the principle directions of stress. To achieve this directionality, while finding a good compromise to the other criteria, the following strategies are employed: To have a fixed or constant intrinsic angle of the fiber, the lines used for construction of the path are set to be equidistant to neighboring lines in all points. To construct lines with this property, curves of equal value on φ are computed with the intersecting plane P being chosen by the user. These isogeodesics are connected at the ends to obtain a continuous path as displayed in Figure 2a for an example surface. This procedure ensures a high density of fiber while retaining a notion of directionality. Drawbacks are high curvatures when changing directions and during travel motions. The inputs for this algorithm are the triangular surface mesh T , the plane P, the spacing d and the travel height h r .
First, the mesh is split into two halves T 1 and T 2 by the plane P to retain a way of distinguishing the sides of the plane as φ is non-negative, which results in a multiplicity of curves for every distance. The seed curve γ for both meshes is calculated by intersecting the seed plane with the triangular surface mesh T . To obtain the geodesic distance, φ(x p ) is evaluated for every point in the mesh.
The isogeodesic curves {c i } i=1,...,n with c i = {x p |φ(x p ) = i · d} are represented as points in space resulting from the interpolated value of the distance function on T 1 and T 2 , respectively.
For every c i , the points in the curve are connected to an ordered path segment s i (k) with k denoting the point index. To connect the lines, they are sorted by distance from the seed curve and iteratively appended to the resulting path. The order of traversal of a curve is reversed every iteration to ensure the back and forth movement. If c i has n discontinuities, travel motions are incorporated between the corresponding segments s i,0 (k) to s i,n (k) by offsetting the respective boundary point by the travel height h r along the normal vector.
Upon computing and connecting both sides, the path is resampled to obtain an even point distribution and smoothed by a Savitzky-Golay filter to eliminate sharp turns.

Fixed Extrinsic Angle (FExA)
It may be of importance to achieve a certain direction of fibers in an extrinsic way, meaning considered as lying in a Euclidean space of higher dimension, in this case, the ambient space of three dimensions. An example for when this would be preferred could be that the load conditions dictate a clear, extrinsic direction, but usage of FInA would result in an intolerable deviation due to high local curvature of the surface. Fibers cannot be tightly packed in this case. Another effect of using the extrinsic perspective is the line being regarded as extended, meaning upon reaching a boundary, the curve will extend over to a possible continuation on the surface. This is not the case for the intrinsic perspective. The drawbacks to this approach are that some surfaces lead to a high number of travel motions and inconsistent spacing can be observed on highly curved surfaces as can be seen in Figure 2b.
Computation differs minimally from FInA. The only difference being the usage of minimum euclidean distance from the plane P instead of geodesic distance from the curve γ. All other computation is conducted equivalently.

Seeded Ortho-Isogeodesics (Ortheo)
The core idea of this algorithm is to generate a path with subcurves that are orthogonal to the isogeodesics computed from a seed curve. This approach stems from the observation that the resulting path can be utilized in combination with the FInA algorithm to form a grid-like structure similar to the carbon fiber fabric used in layer stacking, which can be seen in Figure 3. This composition, where the fibers meet at right angles is unachievable by the classical method of draping a fabric, as the process of draping around a curvature necessarily introduces gaps in the cloth, which are filled with extra segments by the algorithm inducing the need for discontinuities. Additionally, it offers a new notion of directionality: at every point in the path the subcurves point along the gradient ∇φ. It was hypothesized that the maximum principle directions would align with ∇φ if the load is applied in the direction of this gradient at the boundary curves. The algorithm in its current form is in need of optimization for manufacturability as a lot of short segments are included in the path that, when regarding the distance between cutter and nozzle, can impede the ability to extrude the material. As this algorithm has not even been described for the planar case as of yet, the elaboration will be performed in more detail in the following.
The inputs for this algorithm are the triangular surface mesh T , the isocontour resolution d iso , the spacing d, the travel height h r and a triangular mesh T s regarded as the "source mesh" from which the seed curve γ can be calculated. This is performed to enable the use of other existing components as the source, where the load could be applied.
Computation starts by finding γ from the intersection of T and T s . Like before, the isocontours c i = {x p |φ(x p ) = i · d iso } are extracted after applying the heat method, note the usage of d iso for the distance between the isocontours instead of d. The set of path segments are initialized as 1-tupels S path = {s j with j ∈ 1, . . . , l arc d } by extracting equidistant points from γ with spacing d and arc length l arc of γ. The path segments are then grown by iterating over the c i . For every c i+1 , the path segments s j are extended by finding the closest points in c i+1 to the points s j (i) and appending it to the tupel. If the spacing between two chosen points s j (i + 1) and s j+1 (i + 1) on c i+1 exceeds d, a new path segment is started by inserting an additional s j into S path and adjusting the indices for the every s l with l > j. This way, gaps between path segments, which emerge due to the curvature of the surface are filled with new path segments, decreasing the unfilled area and increasing the area fraction while maintaining the directionality of the paths. Similarly, when the spacing decreases enough, segments are ended and removed from S path . The resulting path segments are resampled and smoothed by a Savitzky-Golay filter to eliminate sharp turns. When the iteration reaches the last isocontour the path segments are connected with travel motions by offsetting start-and endpoints of the path segments by h r along the normal direction which yields the final path. A result of this procedure is displayed in Figure 2c for an exemplary surface. Because of the high amount of travel motions, a lot of cutting can occur when fiber paths are short. A possible improvement would be to significantly decrease the number of travel motions by implementing a back and forth motion in a weaving manner, as it cannot be achieved continuously due to the curvatureinduced possibility that startpoints might be unreachable from previous endpoints without crossing already traversed regions.

Contour Parallel Algorithms
As was mentioned in the paragraphs above, the high maximum curvature of direction parallel paths is suspected to impede part quality due to fractures and twists of fibers. Additionally, surfaces with holes inevitably lead to travel motions in most situations. To avoid high curvatures and travel motions, the notion of contour parallel algorithm is introduced. These generally follow the boundary or contour and produce long segments with low curvature. They can also open the opportunity of load adequate path planning from a manufacturing perspective, similar to direction parallel paths, especially if the computed vectors of maximum principle stress closely follow a certain pattern in a subarea of the surface. This pattern can be extracted by the user and will then be followed by the path as the contour is traced and repeated inwards.
One general observation is to be made when dealing with contour parallel algorithms: φ can have several local maxima, making a spiral-like traversal nontrivial. Further description on why the local maxima lead to the impossibility for a trivial solution to the continuous path planning problem can be found in [37].

Geodesic First in Spiral Out (Geo-FISO)
The FISO Algorithm developed in [37] for the planar case offers a possibility for complete continuity with low curvature. The strategy of traversing isocontours by first directly moving into the maxima of the distance transform and then spiraling out along the contours allows for low curvature on inward and outward movement. This enables a continuous path to reach multiple local maxima of the distance transform without travel motions, which would not be possible with classical contour parallel algorithms. To demonstrate the algorithm on a complex and curved surface, the Stanford bunny with an open base was chosen as an exemplary surface as can be seen in Figure 4a. The inputs for this algorithm are the triangular surface mesh T , the spacing d and the travel height h r .
γ is taken as the boundary of the surface. If the surface has holes, γ will have a multiplicity of subcurves from which φ is computed. The isocontours are again given by {c i } i=1,...,n with c i = {x p |φ(x p ) = i · d} with n discontinuities leading to continuous subcurves s i,j with j ∈ {0, . . . , n}.
A weighted digraph G = {V, E} is then constructed with nodes v i,j ∈ V corresponding to the subcurves s i,j and edges E = {e i,j;i+1,k | ∀i, j, k}, i.e., connected for all nodes belonging to c i and c i+1 . The weight function is given by f : E → R + with f (e i,j;i+1,k ) = d min (s i,j , s i+1,k ), with d min denoting the minimum euclidean distance over all points in s i,j and s i+1,k , respectively. A minimum spanning tree T is then computed over this graph, yielding the connections that have to be made when linking the subcurves to a continuous path. This minimum spanning tree is traversed recursively starting with the root node v 0,0 which is always set as the outermost subcontour of γ. In each step of the recursion, the child nodes of the current node v i,j are continually checked for distance, while traversing the current contour and checking for possibilities to reroute from the current position to one of them and back. Once a possibility for rerouting is detected, meaning a segment where the distance between the curve and the current path stay below a threshold t dist for a n thresh path points, the connection is created and a recursive call is made. Important for a correct connection is to detect the winding direction and connect the curves in the orientation that avoids crossing or overlaps. From experience, t dist is chosen as 1.2 · d to allow for some numerical buffer when rerouting the contours. The choice of n thresh depends on the resampling resolution r resample with n thresh = d r resample . Further description of the whole process for the planar case can be found in [37].
Generally, the usage of the geodesic distance from the boundary can lead to voids at the aforementioned local maxima of φ. Unwanted gaps can be created by contouring this function at the fixed values i · d if a local maximum at p fulfills (j + 0.5)d < φ(p) < (j + 1)d for some j. This is addressed by the next algorithm.

Geodesic Contour Parallel Spirals (Geo-Contour)
In some cases, the gaps created by missing contours in the local maxima are of concern as they are suspected to degrade material quality by introducing points where fractures or cracks form and decrease the fiber area fraction. To mitigate this, the global continuity is sacrificed in favor of spirals that travel into the maxima, filling them with a final path segment. This would not be achievable with the with Geo-FISO as two segments, created by splitting a closed curve, are needed to retain the inward-outward movement necessary for the algorithm. The approach described here is displayed in Figure 4b.
The inputs and first steps are identical to those in Geo-FISO. Upon obtaining the minimum spanning tree, all nodes v i,j with degree(v i,j ) > 2 are denoted as rerouting contours and removed from the tree. The remaining connected components represent the "spiralable" areas. These are the areas around local maxima in φ, which can be attributed as belonging to this local maximum. The corresponding curves are then rerouted to obtain a continuous path. For each of these spirals, the points that lie at the local maximum are computed with special caution for numerical errors and an additional curve is appended to the path which traverses these points. After obtaining all spirals, they are connected by travel motions as before and the rerouting contours are traversed.

Geodesic Spiral Around Seed-Point (Geo-Spiral)
The core focus in this algorithm is the continuous curvature of a circular path which was studied in [35], and thus, is better understood than the rapid and irregular direction changes in more classical path planning methods. By selecting a single seed point and creating a spiral that continually increases with geodesic distance, the observed high curvatures when changing directions incorporated in the preceding algorithms are prevented. The big drawback hereby being the inability to completely fill the surface with fiber, taking a hard toll on the fiber area fraction. Nevertheless, this algorithm serves to establish a good control for further testing and experimentation, when checking the material quality, as the induced gaps and the curvature are expected to be minimal.
The input to the algorithm is a point p seed ∈ T on the surface mesh, which is to be selected by a user, the triangular surface mesh T , the spacing d and the travel height h r . The selected point p seed serves as the seed γ for calculating the geodesic distance function φ via the heat method. As before, the isocontours c i are extracted from φ. For this algorithm, the procedure is terminated once a multiplicity is detected, resulting in the advantage of no travel motions and discontinuities, while sacrificing fiber area fraction. To this end, only the c i up to the first i are kept for which |{s i,j } j=1,...,m | > 1. These path segments are then reindexed to have their starting point with minimal distance to the boundary. This way, once the first contour is discontinuous due to hitting the boundary, it can be rerouted at the intersection point, including the open contours with no multiplicity in the path. If there are such open contours they are continued to be traversed with reversing directions until no further connections can be made.

Evaluation
The developed algorithms have been applied to a number of problems in the process of development to study their behavior and adherence to the aforementioned criteria. The following exemplary part and resulting surfaces have been constructed to cover a wide variety of edge cases to serve a better understanding of the algorithm's behavior. An exemplary planned path can be seen in Figure 1c. The part's geometry and by extension the analysis is not exhaustive. This evaluation is meant to serve as a starting point and explanation of the general procedure.
The part used in this evaluation resembles an angle bracket with continuously curved surfaces as can be seen in Figure 5a. The FEA, design and partitioning as explained in Section 1.2 has been performed in FreeCAD and Autodesk ® Fusion360. The load cases examined here are displayed in Figure 6. These are chosen to model a weight exerting force on the angle bracket in a position where the bracket is mounted on a wall or ceiling. As only the maximum principle directions are of interest, the vectors are displayed in normalized form. Upon analysis of the part, taking the shown directions into account, a partitioning into five different surfaces is chosen according to the maximum principle directions, as seen in Figure 5b. The process of partitioning is performed by feature and visually based manual user input. Here, knowledge of the algorithms resulting directionality is required to produce a partitioning that delivers a good result. How many subsurfaces are to be constructed depends on the specific load cases and the part geometry. This is similar to the industry practice of zone-based design in the production of carbon fiber composites in the layup process.

Evaluation Criteria
Every criterion is evaluated without taking segments of travel motions into account, meaning that only the path segments lying on the surface are considered. As the means of acquiring the indices for the criteria are pivotal they are described in the following.

Continuity
Measuring continuity is simply performed by counting the number of rerouting procedures needed to follow the path which corresponds exactly to the number of discontinuities.

Curvature
This criterion is meant to give an index on whether the path contains high curvatures which might induce fractures in the fibers and twists in the fiber bundles, diminishing material quality. For every point in the path, the angle is computed as α(i) = arccos( p i − p i−1 , p i+1 − p i ) with the angles for start-and endpoints of each continuous segment being set to zero. To the end of extracting a meaningful indication on how high the curvature is when changing direction, only the top 1% of computed angles are selected and the mean is calculated yielding the value for this criterion.

Fiber Area Fraction
To compute the fiber area fraction, the area A s of the triangulated surface is first extracted by summing the area of its triangles. Then, the geodesic distance field φ path on the mesh is calculated with the path points serving as the source γ. The area for which φ path (x p ) > d is the uncovered area A u . This is calculated by summing the triangle areas of a subdivided mesh, where all vertices fulfill φ path (x p ) > d. Lastly 1−A u A s gives us the index for the area fraction. This value indicates the amount of used area, with 0 representing no fiber deposited and 1 meaning complete coverage.

Load Adequacy
The index for load adequacy is calculated for every point excluding the last point in the path by finding the closest maximum principle direction vector v from the data set and checking for alignment by computing the dot product v, u with the vector u = p i − p i+1 connecting the current point with the next one. Before calculation u is normalized, as the sign does not matter in this case, the absolute value of the resulting dot product is taken. If several load cases are to be examined at this point, the highest alignment is chosen. Finally, the index for load adequacy is calculated as the mean of the computed dot products for all points. This value will be between 0 and 1 representing the amount of alignment with the load transmission paths.

Results
Of all possible combinations of the surfaces depicted in Figure 5b and all the algorithms, ten were chosen and analyzed. These have been selected to cover all algorithms in several scenarios, the only exception being Ortheo and Geo-Spiral as the performance is evident from one test. The results from this evaluation are shown in Table 1 and a selection of the chosen combinations can be seen in Figure 7a-e. What can be observed is that the Geo-Contour and Geo-FISO algorithms show a low curvature while maintaining a level of fiber area fraction and load adequacy. The direction parallel algorithms excluding Ortheo show higher curvature because of the sharper angles when changing directions but have a consistently high fiber area fraction. They can be observed to incorporate discontinuities for parts with holes or complex and large surfaces. Ortheo yields especially high load adequacy and the lowest curvature while introducing a lot of discontinuities as expected. Geo-Spiral performs poorly in every aspect with even curvature being high, which is suspected to have its cause in the rerouting, which could be improved by smoother transitions. In total, the results show the trade-offs and differences supporting the claim that having multiple algorithms is necessary to adapt to the varying conditions when planning the reinforcement patches.  To the end of analyzing the variety of applications, a qualitative feature-based matching is made in the following. All algorithms are numerically stable and deliver results in edge cases, as well as highly curved surfaces. Additionally, computation times were observed to be in the lower minute range. The primitive cartesian smoothing with the Savitzky-Golay filter can lead to collisions with the surface for highly curved surfaces if its degree is set too high.

Direction Parallel Algorithms
The direction parallel algorithms exhibit a higher degree of control outside the partitioning as the direction can be set. The user-defined seed curve has a direct influence on the results and is as relevant as the subsurface boundary and thus underlies the same principle of user experience being necessary in application. In general, the performance of direction parallel algorithms correlated with the relative dimension in the chosen direction. Especially for long and thin parts this choice has an effect. Additionally, highly concave structures lead to a sharp increase in travel motions for the direction parallel algorithms as such gaps prevent the movements from being achieved without discontinuities.

Contour Parallel Algorithms
The major influencing factor for the performance of the contour parallel algorithms is the boundary geometry and the resulting distance field. If φ has a high number of maxima Geo-Contour will show a lot of travel motions. In this case, the choice of Geo-FISO is favorable. Geo-Spiral has a major dependency on boundary geometry. On convex and highly curved surfaces it can terminate due to the isocontours self intersecting. Thus, complex geometries with protruding features lead to a poor performance. In summary, Geo-FISO offers the best performance for the widest applicability.

Conclusions
A procedure is proposed in this paper including six algorithms for planning paths of reinforcement patches on complex and curved surfaces. This represents an important step into field planning paths for the multi-axis manufacture of CFRPs with FDM which is currently one of the major bottlenecks for this process. Evaluation criteria were formulated and motivated by the process to enable the guided development and serve as a metric for methods developed in the future. The algorithms were described and evaluated regarding the established criteria to assess the applicability and the spectrum of possible surfaces to be covered. The concept of isocontours on geodesic fields is employed to yield paths with low curvature, high fiber density, few discontinuities and high load adequacy. Methods to obtain these values were developed and applied to measure the performance of the algorithms. The presented results show promising capabilities for this application and others such as planning paths of slices when constructing a volume body. To develop and evaluate such an approach will be highly beneficial for load adequacy as no multi-axis slicer software for CFRPs has yet been made available.
Further work includes the possibility for some of the proposed methods such as Geo-Spiral to be optimized to enable termination under less restrictive conditions and further increase performance. Domain-based smoothing could eliminate defects introduced by currently used smoothing methods. Additional analysis of the process of partitioning and the load cases for FEA is needed to increase the applicability of the described methods and quality of results. An automation of partitioning and planning would be favorable, as user experience and knowledge of algorithms is still required in the process. Moreover, a comprehensive evaluation and experimental validation is needed to increase the effectiveness of the deployment and application.
This paper precedes any complete analysis of the multi-axis manufacturing process for CFRPs including slicing, construction and planning and gives a starting point for the promising field of planning paths for CFRPs manufactured by use of FDM.  Data Availability Statement: Computer code will be made available upon request to the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: Seeded ortho-isogeodesics Geo-FISO Geodesic first in spiral out Geo-Contour Geodesic contour parallel spirals Geo-Spiral Geodesic spiral around seed-point