Convolutional approach to highly efficient wendland radial basis function level set structural optimization

Abstract

Systems, devices, methods and computer program products are described herein for a convolutional approach to highly efficient single parameter radial basis function level set structural optimization and inverse homogenization.

Description

CONVOLUTIONAL APPROACH TO HIGHLY EFFICIENT WENDLAND RADIAL BASIS FUNCTION LEVEL SET STRUCTURAL OPTIMIZATIONCROSS REFERENCE TO RELATED APPLICATIONS

[0001] This application claims priority to and benefit of U.S. provisional patent application 63 / 543,582 filed October 11, 2023; U.S. provisional patent application 63 / 668,024 filed July 5, 2024; and U.S. provisional patent application 63 / 668,590 filed July 8, 2024, each of which is fully incorporated by reference and made a part hereof.BACKGROUND

[0002] It is estimated that more than 200 million people worldwide are currently living with bone disease, and over 1.5 million Americans sustain bone fractures every year. Despite the prevalence of skeletal surrogates (see, Fig. 1), none have been able to capture its spatially variate 3D anisotropy (SV-3DA), inhibiting efforts to understand the effect of microstructural degeneration.

[0003] Current structural optimization methods suffer from a tradeoff between feature fidelity and computational cost. This is a result of mathematical simplifications associated with simplistic structural representations. The simplest, and therefore the most popular, method for structural optimization is the Single Isotropic Material with Penalization (SIMP) method. In this method, a region is discretized into cubic elements. The amount of material in these elements is modulated (uniformly within a given element) to adjust the respective element's contribution to a particular structural response. The gradient analysis here is straightforward, and as such, this approach has been widely adopted. However, this approach suffers from gradient vanishing issues, as well as design fidelity issues. As such, in commercial applications, it is often followed by a brief level set optimization. This level set optimization,referred to as the Level Set Method (LSM), is commonly highly intense computationally, and theoretically. A large part of the computational expense is that the nodes of the elements are moved in order to modulate shape, being that this is a shape optimization approach. It then yields a natural and printable / pre-manufacturable design. This is heavily used in the aerospace and automotive industries, among other industries. However, another issue is that the LSM is dependent on dense linear systems inversions, under classical theory. Because of this, the LSM has been limited to very small problems, dominantly 2D.

[0004] Therefore, systems, methods, devices and computer program products are desired that overcome challenges in the art, some of which are described above.SUMMARY

[0005] Systems, devices, methods, and computer program products are described herein that perform derivations on the LSM which make it so that the SIMP method is unnecessary, along with any intense shape fine-tuning. Departing from the status quote of linear systems and preliminary SIMP topology optimization, hexahedral element nodes are fixed in space, and material domains are adjusted on a fixed grid. The result is a structural optimization approach, powered by Graphical Processing Unit (GPU) accelerated Fast Fourier Transform (FFT) convolution. This methodology benefits from assumption based simplifications on the LSM, which actually make it almost as efficient as the SIMP method even for 3D problems. The algorithmic implementation approach disclosed herein addresses the issue of not being able to generate detailed 3D structures in optimization.

[0006] The disclosed systems, devices, methods, and computer program products provide a means of efficient next generation structural optimization utilizing gaming GPUs, without the need for supercomputers or quantum computing. The GPU is used to accelerate simulation and the FFT convolution of the compactly supported radial basis kernel with the design variables, the Heaviside of which defines the material domains. The disclosed approachhas been proven to design and print structures with effectively equivalent anisotropy to that of skeletal microstructures. This was done in material which had differing elastic properties to that of bone. This problem has neither been approached in literature nor commercially due to its high computational intensity, classically. As such, it is evident that the approach is capable of meeting next generation structural optimization demands. This will likely have an influence on many sectors which use computational optimization, including finance, medicine, military, transportation, civil engineering, logistics, etc. Any area where optimization is used will likely benefit from the invention.

[0007] In one aspect, a convolutional approach to sub-Gauss Point 3D inverse homogenization is disclosed and examples provided.

[0008] In another aspect, a 3D triclinic inverse homogenization assessment of design material elasticity and sub-femoral region is disclosed. This includes examples of high- throughput implementation on several cadaveric femoral scans, sampling multiple sub-femoral regions and designing in a range of materials. Further disclosed is GPU convolution accelerated high-resolution segmentation.

[0009] In yet another aspect, design and mechanical testing of bone with targeted triclinic spatially variate elasticity, and adjacent microstructure synergy is disclosed, which adapts the convolutional approach to sub-Gauss Point 3D inverse homogenization algorithms described herein. One disclosed implementation comprises performing the testing utilizing 3D with multiple triclinic elasticity targets, yielding up to 288 mechanical continuum constraint. This aspect further discloses and implements strategic concatenation of design variables and sensitivity, and de-concatenation for isolated design update.

[0010] Other systems, methods, features and / or advantages will be or may become apparent to one with skill in the art upon examination of the following drawings and detaileddescription. It is intended that all such additional systems, methods, features and / or advantages be included within this description and be protected by the accompanying claims.BRIEF DESCRIPTION OF THE DRAWINGS

[0011] The components in the drawings are not necessarily to scale relative to each other. Like reference numerals designate corresponding parts throughout the several views.

[0012] Fig. 1 illustrates an example of a biomechanical surrogate.

[0013] Fig. 2 illustrates an example method of bottom-up design used for spatially variate anisotropy mimicry.

[0014] Fig. 3 illustrates exemplary trabecular super-resolution (SR) performance.

[0015] Fig. 4 illustrates exemplary trabecular SR metrics.

[0016] Fig. 5 A illustrates a 80μm res. μCT; Fig. 5B illustrates a 6mm cube sample; Fig.5C illustrates a homogenized elasticity field; Fig. 5D illustrates anisotropy equivalent printable metamaterial; and Fig. 5 E illustrates exemplary fabrication.

[0017] Fig. 6A illustrates implicit masking in nTop; Fig. 6B illustrates elasticity mapping; Fig. 6C illustrates implicit geometry feature variation for isotropy mimicry; and Fig. 6D illustrates fabrication (IND 405, Carbon Ml).

[0018] Fig. 7 illustrates an exemplary triply periodic Gyroid, via nTop v4.25.2 where the shaded region represents the periodically repeating base unit cell.

[0019] Fig. 8A illustrates a 3D Level Set Function (LSF, φ(x)) with void domain (D), boundary domain (Γ), and solid domain (Ω) identified. These correspond to φ(x)<0, φ (x)=0, and φ(x)>0, respectively.

[0020] Fig. 8B illustrates that they do not overlap, and they ultimately define a solid area, which is expanded into 4D-LSF / 3D-domains as described herein.

[0021] Fig. 9A illustrates x is the coordinate being evaluated, r_i is the distance of x from the knot center (x_i) as a function of support radius (R), showing that the ith knot; and Fig. 9B illustrates evaluation of ith Wendland function at x, and (*)+ is described as a truncated power function, which allows knot i to contribute to the level set value of point x. As mathematically expressed, and pictorially depicted in Fig. 9B), the CS-RBF attenuates to zero at the compact support radius, R.

[0022] Fig. 10 illustrates spatially relative Wendland function values.

[0023] Fig. 11 illustrates a simplified illustration of the proposed approach to store DoF- by-DoF non-zero coefficients.

[0024] Fig. 12 illustrates a sigmoid adaptive move limit diagram.

[0025] Figs.l3A-13E illustrate FEA based computational homogenization has been used to reduce the complexity of large microstructural models, such as that shown in Fig. 13 A. Fig. 13A illustrates an 82 μm resolution HR-pQCT; Figs.13B and 13C show 4 mm trabecular patch ( 13B) and cubic sample ( 13C) ; FIG. 13D shows a homogenized elasticity field (via nTop 4.18.2 ). Note how the elasticity map is aligned with the path of principle compressive stress, supported by the GCM (Fig. 13E).

[0026] Figs.l4A-14F illustrate the outcomes of the statistical analysis described herein.All statistics presented heteroscedasticity, and as such, medians (short horizontal line within shaded area in notch plots) are those whose groupwise differences are statistically compared.

[0027] Figs.l5A-15G show MMA 3D IHP solution via CS-RBF cPLSM in overly rigid design material, where 15A shows a final design; 15B shows IPEI_ij matrix bars; 15C shows IHP iteration history; 15D shows XY projection of 15E, 3D elasticity field, derived from GCM; 15F shows XZ projection; and 15G shows YZ projection. ~18 seconds without visuals.

[0028] Figs. l6A-16G show SOMMA 3D IHP solution via CS-RBF cPLSM in overly rigid design material. Fig. 16A shows a final design; FIG. 16B shows IPEI_ij matrix bars; FIG.16C shows THP iteration history; FIG. 16D shows XY projection of Fig. 16E, 3D elasticity field, derived from GCM; FIG. 16F shows XZ projection; and FIG. 16G shows YZ projection. ~43 seconds without visuals.

[0029] Fig. 17 illustrates the solutions of Figs.15 and 16 displayed via nTop, serving as validation.

[0030] Fig. 18 illustrates that the design is periodic.

[0031] Fig. 19A shows a Volume slice; FIG. 19B shows FFT segmentation; Fig. 19C shows Femoral Head (FH), Neck Superior (NS), Neck Inferior (NI), Trochanter Superior (TS), Trochanter Middle (TM), and Trochanter Inferior (TI) sub-segmentation; Fig. 19D shows sample selection; and FIG. 19E shows the sample.

[0032] FIG. 20A illustrates a bone sample volume fraction, extreme outliers removed in visual. FIG. 20B shows a bone sample Hill’s Bulk Modulus, K_Hill, extreme outliers removed in visual. FIG. 20C shows IH symmetric Mean Absolute Percent Error (sMAPE). FIG. 20D shows percent error in K_Hill (PE(K_Hill )) of optimal design relative to skeletal target, extreme outliers removed in visual. FIG. 20E shows optimal design volume fraction percentage, f_(V%). FIG. 20F shows design iterations. * indicates p<0.05, ** indicates p<0.01, *** indicates p<0.001. Red marker indicates outliers. Rank based comparison of group medians used throughout.

[0033] Fig. 21 A shows symmetric Mean Absolute Percent Error across the constitutive matrices CTar9and C* (sMAPE, %) with respect to Elastic Modulus (E, GPa) in 14 materials (see Table 2, above. 90% 5,000 cycle bootstrap based Confidence Intervals (Cis) on median at each material, pooled across anatomic regions (material pooled and regionally stratified in Fig. 1C-21D. The 20% maximum sMAPE cutoff is measured, the solution spanning materials wherein maximum 90% CI is below or at the threshold (0.63 GPa ≤ E ≤ 30.70 GPa). Few notable materials are indicated with arrows pointing to respective point estimates (mean, splineinterpolated blue dotted dashed curve). 90% CI domains are given as spline interpolated transparent blue fills with the bounding approximated 90% limits as dashed blue curves. For all materials, the 90% Cis are given in orange with a dot at the computed mean and the caps terminating at the computed 90% CI ends. FIG. 21B shows PE(KHill) as a %. The transparent red shaded region represents the area covered by the elasticity bounds in which median sM APE is at most 20% with 90% confidence. FIG. 21C shows for the final designs. Note saturationat or near 30.70 GPa. FIG. 21D shows the number of design iterations prior to termination. Note the saturation, a result of intermediate material removal at iteration 60 and subsequent termination.

[0034] Figs.22A-22G show SIL30 IH design history, where FIG. 22A shows a final design; FIG. 22B shows a GCM component percent errors; FIGS.22C-22F show projections of the target and final axial elasticity surfaces; and FIG. 22E shows iteration history.

[0035] Figs. 23A-23G illustrate 18 GPa Elastic Modulus sample design history, whereFIG. 23 A shows a final design; FIG. 23B shows GCM component percent errors; FIGS. 23C- 23F show projections of the target and final axial elasticity surfaces; and FIG. 23E shows iteration history.

[0036] FIG. 24 shows target (left) and optimal (right) units cells at 18 GPa ElasticModulus. High elastic moduli are associated with less dense, high percent error, and high iteration count (Fig. 25). The given Ti6A14V example terminates with an sMAPE of 253.27%, PE(K_Hill ) of 92.86%, f_(V%) of 23.52%, and at 62 iterations. This is in line with the presented trends, as it shows how the design remains in an intermediate state between voluminous posts. Due to the difficulty associated with forming a lean, rigid connection, iteration proceeds for extended iterations until a hard cut in intermediate material occurs for termination at 60 iterations.

[0037] Fig. 25 shows Ti6A14V IH design history, where FIG. 25 A shows a final design;FIG. 25B shows GCM component percent errors; FIGS. 25C-25G show projections of the target and final axial elasticity surfaces; and FIG. 25E shows iteration history.

[0038] Fig. 26 illustrates the central skeletal microstructure in a 3x3x3 discretization and compounds relevant to the disclosed sequential cIH procedure

[0039] Fig. 27 illustrates a graphic procedure for the CS-RBF cPLSM cIH, which shows a workflow for the sequential CS-RBF cPLSM cIH procedure.

[0040] Fig. 28 illustrates the discretized femoral benchmark, with bone and discretization bounds shown in different shading.

[0041] Fig. 29 shows the central UC design and relevant compounds (columns). (Top) skeletal microstructures. (2nd row) inverse designs. Intermediate material is shown in blue, the methodologically relevant adjacent designs are in gray. The current design is in red. (3rd row) axial elasticity map. (Bottom) log-y scale function evaluation curves. The design variables driving the primary design, and all compounds are updated simultaneously.

[0042] Figs. 30A and 30B show the overall benchmark design (30B), split by cell, and the skeletal sample from Fig. 28 (30A) (duplicated for visual comparison). q=0.5 is used. Note similar spatially relative porosity and directionality, particularly in the mid-section. Figs. 30A and 30B show visualized via nTop, discretized benchmark sample (30A) and Ti-6A1-4V design (30B) via disclosed sequential CS-RBF cPLSM cIH procedure. Two-UC replication is used 01=0.5).

[0043] Figs. 31A-31F show the metric distributions for the benchmark across N.Statistical analysis is not performed as N=0 occurs once. Note that the evolving SOMMA approach yields consistent iteration count, hMAPE, sMAPE, f_(V%), and PTI. Duration increases seemingly linearly with N. In cIH, a total of 3.55 hours is spent in cIH, not considering forward homogenization evaluation of the skeletal microstructure(s) in thebeginning for target calculation. Fig. 31 shows cIH performance metrics on benchmark across N. Outliers not shown.

[0044] Fig. 32 gives the distribution of N across the cadaveric LI discretization. Fig.32 shows primary UC count vs. N.

[0045] Figs. 33A-33F show the metric distributions for the full vertebral body across N, as well as statistics.

[0046] Fig. 34A shows a cross-section of trabecular microstructure (yellow) and cortex(purple) of scanned vertebral body. FIG.34B shows a cross-section of designed model. FIG. 34C shows a supported model fabricated in Clear Loctite, on a build-plate.

[0047] Fig. 35 gives (35A) the processed force-displacement curves following uniaxial compression for the five surrogates and the cadaveric specimen, (35B) a boxplot of median surrogate stiffness and a line for specimen stiffness, and (35C) the percent error in surrogate stiffness relative to specimen stiffness, with a range in expected percent error given by mean hMAPE and sMAPE weighted by N (right). Per a one-sample two-tailed t-test, specimen stiffness is significantly lower than mean surrogate stiffness by 0.655 kN / mm (p « 0.05). This is anticipated as poor mechanical accuracy was yielded by the MMA approach used.

[0048] Fig. 36 illustrates an example computing device that may be used to practice aspects described herein.DETAILED DESCRIPTION

[0049] Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art. Methods and materials similar or equivalent to those described herein can be used in the practice or testing of the present disclosure. As used in the specification, and in the appended claims, the singular forms “a,” “an,” “the” include plural referents unless the context clearly dictatesotherwise. The term “comprising”, and variations thereof as used herein is used synonymously with the term “including” and variations thereof and are open, non-limiting terms. The terms “optional” or “optionally” used herein mean that the subsequently described feature, event or circumstance may or may not occur, and that the description includes instances where said feature, event or circumstance occurs and instances where it does not. As used herein, “exemplary” means an example of and is not intended to denote a preference or a preferred embodiment. Ranges may be expressed herein as from "about" one particular value, and / or to "about" another particular value. When such a range is expressed, an aspect includes from the one particular value and / or to the other particular value. Similarly, when values are expressed as approximations, by use of the antecedent "about," it will be understood that the particular value forms another aspect. It will be further understood that the endpoints of each of the ranges are significant both in relation to the other endpoint, and independently of the other endpoint.

[0050] The following acronyms may be used herein:

[0051] CT - Computed tomography

[0052] Al - Artificial intelligence

[0053] HR-pQCT - High resolution peripheral quantitative CT

[0054] NbN - Node by node

[0055] PCG - Preconditioned conjugate gradients

[0056] cIHP - compound inverse homogenization problem

[0057] cPLSM - convolutional parametric level set method

[0058] GPU - Graphical processing unit

[0059] SRGAN - Super resolution generative adversarial network

[0060] WDSR-B - Wide activation deep super resolution network type B

[0061] FEA - Finite element analysis

[0062] CS-RBF - Compactly supported radial basis function

[0063] CNN - Convolutional neural network

[0064] UC - Unit cell

[0065] FFT - Fast Fourier transform

[0066] cLSM - convolutional level set method

[0067] GCM - Generalized constitutive matrix

[0068] DbDPCG - DoF by DoF PCG

[0069] DoF - Degree of freedom

[0070] PTI - Percent total intermediate

[0071] SOMMA - Second order MMA

[0072] MMA - Method of moving asymptotes

[0073] sMAPE - symmetric mean absolute percent error

[0074] Systems, devices, methods and computer program products are described herein for a convolutional approach to highly efficient single parameter radial basis function level set structural optimization and inverse homogenization with degrees of freedom (DoF) by DoF finite element analysis. Derivations are performed on the LSM which make it so that the SIMP method is unnecessary, along with any intense shape fine-tuning. Departing from the status quote of linear systems and preliminary SIMP topology optimization, hexahedral element nodes are fixed in space, and material domains are adjusted on a fixed grid. The result is a structural optimization approach, powered by Graphical Processing Unit (GPU) accelerated Fast Fourier Transform (FFT) convolution. This methodology benefits from assumption based simplifications on the LSM, which actually make it almost as efficient as the SIMP method even for 3D problems. The algorithmic implementation approach disclosed herein addresses the issue of not being able to generate detailed 3D structures in optimization.

[0075] In particular, a bottom-up two-level structural optimization approach is disclosed that enables fabrication of anisotropy matched skeletal surrogates from patient data.

[0076] There are two main approaches to 2-level structural optimization: bottom-up and top-down. Top-down uses known macroscale loading to infer microstructural detail. Bottom- up uses spatially distributed microstructural target parameters to infer microstructural detail.

[0077] As shown in Fig. 2, a bottom-up design is used herein for spatially variate anisotropy mimicry.

[0078] The result of Step 1 from Fig. 2 is super-resolution. Here, 2D and 3D variants of the Wide Activation for Efficient and Accurate Image Super-Resolution (WDSR)-B Structural Recurrent Neural Network (SRNN) were trained on microstructural data for 4X SR. The SRGAN (Super Resolution Generative Adversarial Networks) perceptual loss framework was used to fine-tune SRNNs. Referring to Figs. 3 and 4, it can be seen that spectacular performance is achieved in peak signal-to-noise ratio (PSNR), trabecular thickness (Tb.Th.) % error, and anisotropic elasticity tensor accuracy.

[0079] Referring back to Fig. 2, the results of step 2 of the method shown in Fig. 2 is a microstructural inverse design. As shown in Fig. 5, the LSM is used to mimic the pure anisotropy of each unit cell, producing smooth periodic binary designs in high detail.

[0080] Referring again to Fig. 2, step 3 of Fig. 2 illustrates exemplary macrostructure assembly. As shown in Fig. 6, isotropy variable 2-level assembly has been performed

[0081] It should be appreciated that the logical operations described herein with respect to the various figures may be implemented (1) as a sequence of computer implemented acts or program modules (i.e., software) running on a computing device (e.g., the computing device described herein), (2) as interconnected machine logic circuits or circuit modules (i.e., hardware) within the computing device and / or (3) a combination of software and hardware of the computing device. Thus, the logical operations discussed herein are not limited to any specific combination of hardware and software. The implementation is a matter of choice dependent on the performance and other requirements of the computing device. Accordingly,the logical operations described herein are referred to variously as operations, structural devices, acts, or modules. These operations, structural devices, acts and modules may be implemented in software, in firmware, in special purpose digital logic, and any combination thereof. It should also be appreciated that more or fewer operations may be performed than shown in the figures and described herein. These operations may also be performed in a different order than those described herein.A Convolutional Approach To Sub-Gauss Point 3D Inverse Homogenization

[0082] Skeletal surrogates are used in many fields within medicine, such has human body modeling, medical device testing, and even as implants. Yet, the tools for rapid development of micromechanically matched anisotropic continua are few and far between. Described herein are tools for efficient, high-resolution, generalizable inverse homogenization in 3D. In some instances, the disclosed tools may be implemented on consumer grade hardware. For example, 30 inverse homogenizations are performed, each on two systems featuring NVIDIA graphics processing units (GPUs); one a hybrid tablet grade GTX 1060, and the other a desktop grade RTX 3090. Tests on different classes of anisotropy are performed; a skeletal triclinic sample and an artificial cubic Schoen Gyroid sample. Additionally, the influence that incorporation of the analytic Hessian has in the disclosed formulated structural optimization problem is assessed. Such modeling is done on a design material with a 1.8 GPa elastic modulus, DPR 10. Subsequently, validation testing is performed on the algorithms utilized with an overly stiff material having a 10 GPa elastic modulus, RigidlOk. Outcomes include: successful development of the disclosed convolution powered Level Set approach; integration of optimized pseudo-boundary adjustments for adaptive sensitivity; analytic Hessian formulation and deployment in the context of 3D LS IH; and effective generalization. Additionally, we succeed in validation on commercial modeling software, nTop. Implications of these findings include potential for continuum design with prescribed effective mechanicsin the computational modeling and additive manufacturing communities. Finally, application to bone may help yield the next generation of optimized skeletal surrogates.

[0083] In the early-90s, a landmark discovery in sub-discretization modeling was made in graphic design. This discovery has only recently been exploited to make sub-element structural optimization efficient. This was the discovery that convolution can be used to evaluate a B-Spline based Level Set Function (LSF), abstracting to the Radial Basis Function (RBF). In 2019, Wang et al. publicized their later patented [Rai R, Wang J. Periodic Cellular Structure Based Design for Additive Manufacturing Approach for Light Weighting and Optimizing Strong Functional Parts. US 20210216683A1, 2021, incorporated by reference] implementation of convolutional B-Spline collocation, which extends to the LS sensitivity analysis. This significantly accelerates typically tedious collocation (evaluation).

[0084] Without a convolutional approach to LS collocation, the Level Set Method(LSM) of Structural Optimization is intractable for large 3D problems. Examples of the standard dense linear system inversion procedure are presented in most literature on the topic of the LSM or LS modeling. Limitations of the standard approach include interpolation-based inaccuracy, and dependency on the assembly and inversion of a very large and dense linear system of equations.

[0085] A redeeming factor is the binary representation of the design domain, relative to the Single Isotropic Material with Penalization (SIMP) approach. By imposing a piecewise constant parametric component to the CS-RBF LSM, the CS-RBF Parametric LSM (CS-RBF PLSM) is attained. The LSM classically suffers from issues involving void definition and solution of the Hamilton- Jacobi PDE. To combat these issues, re-initialization was instituted, where the collapse or emergence of old or new voids can occur. Via piecewise separability of the LS, this method avoids solution of the Hamilton- Jacobi PDE, re-initialization, and even allows for voids to be opened and closed. This led to efficient LSM Structural Optimization.

[0086] A natural step forward is the disclosed CS-RBF convolutional called, PLSM(CS-RBF cPLSM). Via convolutional LS collocation, it avoids interpolation-based inaccuracy when computing the values of the LSF. It achieves enhanced evaluation accuracy without the need to form, simplify, or invert any dense matrices. The CS-RBF cPLSM approach we choose going forward avoids solution of the H-J PDE, re-initialization, and collocation; the primary sources of LSM computational cost.

[0087] Engineered continua are integral within structural design. The InverseHomogenization Problem (IHP) was first dubbed by Ole Sigmund in 1994. This was near the time in which homogenization was popularized in biomechanics, and with homogenized biological descriptors comes the question: if a material is had with constitutive properties differing from bone, how could this material be arranged to recreate effective mechanics? A gap exists in that this question has yet to be thoroughly explored in 3D with level set design accuracy for bone.

[0088] Typically, the Hessian is avoided when using a level set approach. This technical gap arises due to the computational intensity associated with LS evaluation / collocation, and typical re-initialization. The disclosed convolutional abstraction enables robust level-set Hessian evaluation even in 3D, providing further information for gradient based solvers.

[0089] The LS design boundary is moved spatially to evolve the design. A pseudo width of the boundary is imposed to enable sensitivity analysis spanning multiple knots. This pseudo width is set to be broad initially and may be reduced as the design nears optimality. However, the pseudo width has yet to be considered as a design variable in literature. A gap in literature presents itself in that the design boundary cannot grow if needed. This leads to prolonged evolution and excessive iterations. Disclosed is a continuous quantifier of intermediate material presence and a gradient based approach to pseudo boundary width control, duallyaccounting for constraints on homogenization targets. As such, design optimization is expedited.

[0090] Previous literature has demonstrated the sensitivity of a design with respect to extreme anisotropy. An important feature of such IHP solution approaches is the squaring of GCM residuals. This effectively enables a targeted reduction of large magnitude errors. However, to tightly mimic pure anisotropy, all terms must be closely matched. Yet, some may not need to be as tightly matched as others. For the disclosed constitutive constraints, a maximum target component normalization is utilized, and relative magnitude-based target percent error bounds. The efficacy of a single algorithm is demonstrated in two cases of anisotropy: a skeletal sample with triclinic elastic anisotropy, and a gyroid sample with cubic elastic anisotropy, both having equal bulk modulus (0.49 GPa). The cubic scenario features negligible coupling terms for dilatational (expansion-contraction) and distortional (shearing) deformation. In normalizing, the disclosed approach is robust to both small and large targets.

[0091] The following is a thorough introduction to the IHP and the steps from the LSM into the disclosed CS-RBF cPLSM. We see three essential representations in homogenization and inverse homogenization in order of decreasing storage size: 1) the fourth-rank elasticity tensor (81 components in 3D); 2) the more common symmetry based reduction to a square generalized constitutive matrix (GCM; 36 components in 3D); and 3) the further reduction to a set of unique GCM components (# depends on assumptions of anisotropy type, but at most 21 unique components for anisotropic elasticity). The source paper on inverse homogenization, [Sigmund O. Materials with prescribed constitutive parameters: An inverse homogenization problem. Int J Solids Struct 1994;31 :2313-29. https: / / doi.org / 10.1016 / 0020-7683(94)90154- 6, incorporated by reference], represents 1) as 2) as [E], and 3) individually as: forthe 7thunique tensor components’ ethelement mutual energy (via strain energy density, or SED), where^ represents a conversion from a fourth-ranktensor to a vector of unique ijkl combinations; and Dj for the target prescribed constitutive parameter for the unique tensor component. This is a fine approach, as one is a tensor, one is a matrix, and one is a single index vector. Additionally, the lettering is consistent. However, some literature in the field uses differing annotation, such as for the elasticity tensor, both homogenized and for the target. Meanwhile, other literature in the field represents the elasticity tensor as E^ , but represents the GCM as C, or by component as Qy. Other literature represents the tensor as clJklor c, and the matrix similarly to Sigmund as [c], Some literature only uses the matrix and refers to it as [C].

[0092] We will indicate: 1) the fourth-rank stiffness tensor as c, in Einstein notation asCijki', 2) the GCM as C, in Einstein notation as Qy; and 3) we will not be reducing it to unique components as we find, in practice, that the symmetric components are best accounted for byGCM occurrence in optimization. Note that (1):(1)

[0093] In matrix condensed and Einstein notation, know that Hooke’s law is as follows(2):(2)and that in tensor condensed and Einstein notation, it is (3):(3)

[0094] Mostly, Einstein notation will be used and will make it clear what is taking place.

[0095] Original Inverse Homogenization Problem Background

[0096] The IHP is meant to yield microstructures with prescribed effective properties.Given a prescribed GCM (4):

[0097] (4)

[0098] with d x d dimensionality (d=3 in 2D, 6 in 3D), the IHP is concerned with identifying the most optimal structure which achieves said mechanics subject to some additional constraint and objective. In Sigmund’s paper, a formal problem is posed in part by isolating and flattening unique components from (4) also shown in (5):

[0099] (5)

[0100] (5) minimizes weight, W, by changing beam / element, e, cross-sectional area / density, xe, across all elements, Ne, where a given node has d Degrees of Freedom (DoF)(see

[0044] for further information). Constraints include the range in the design variable, (i.e., cross- sectional area), the (cZ + (cZ2— cZ) / 2) independent unique flattened prescribed constitutive parameters, Df, and the residual of the sum of elemental energies, (6):

[0101]

[0102] where Yeis the volume of an element, Cpqrsis the constitutive property of the bulk material,is prescribed strain, and E* is the strain field from applying the prescribed strain to the microstructure in a periodic continuum fashion (7):

[0103] 1 s7' >){v: v is Y — periodic}

[0104] Wherein the displacement field is solved under ε°. ε* is computed from the displacement field using the strain-displacement matrix, as discussed in the following. The periodicity mentioned above is shown in 3D by Fig. 7, where FIG. 7 illustrates an exemplary triply periodic Gyroid, via nTop v4.25.2. The shaded region represents the periodically repeating base unit cell.

[0105] Reff ering to Figure 10 in Sigmund, above, note the portrayed solution to the IHP posed for the target 2D-GCM, portrayed in the same figure. The truss structures meet the constraints on effective properties. It is noteworthy that only terms resposible for orthotropic anisotropy (not exaustive however, as the shear term is absent) are targeted. The concept of periodicity is critical to homogenization. In practice, from a FEA perspective, periodic boundary conditions are imposed by setting nodal DoF displacements on one side of the FEM equal to those on the opposite side, doing so in all planes / axes, and tying a single node to stationarity while applying the prescribed macro-strain. The disclosedproblem specific implementation is described further herein, below.

[0106] A Level Set (LS) defines different values that a scalar field has as a contour(surface) in 2D (3D+). In their 2014 paper, Wang et. al. [Wang Y, Luo Z, Zhang N, Kang Z. Topological shape optimization of microstructural metamaterials using a level set method. Comput Mater Sci 2014;87:178-86, incorporated by reference] provides an exceptional explanation of the Parametric LS (PLS) and the PLS Method (PLSM). What follows is a summary. At a given value of the LS, there exists a boundary, Γ, of a 2D (3D) domain / structure, Ω , or on a 3D (4D) surface (8) (see Fig. 1):

[0108] Fig. 1 A illustrates a 3D Level Set Function (LSF, φ (x)) with void domain (D), boundary domain (F), and solid domain (Ω ) identified. These correspond to φ (x) < 0, φ (x) = 0,and φ(x) > 0, respectively. Fig. 8B illustrates that they do not overlap, and they ultimately define a solid area. This is expanded into 4D-LSF / 3D-domains within the work herein.

[0109] The teal (light shading) domain in Fig. IB is the solid region corresponding to the 1strow in (8) after assessing the region of x where the LSF, φ(x), is positive in Fig. 1 A. The boundary and void (D) follow. Given an arbitrarily high dimensional field, any lower dimensional domain can be defined analytically. This concept translates well into the space of structural optimization, where a parametric LS (PLS) can be tuned to define a domain that satisfies constraints or minimizes objectives on structural performance with fine boundary accuracy. First, the basis for the PLS, the Radial Basis Function (RBF), must be introduced. In particular, a Compactly Supported RBF (CS-RBF).

[0110] The CS-RBF employed by Wang et. al. is named the Wendlandn function afterHolger Wendland who defined the functions and explored their mathematical properties. He gave definition to the function with subscripts 3 and 1 which represent a single integration on the function raised to the 4th(i.e., (3+ 1)*11) power as follows in (9):

[0111] (9) ^

[0112] where is the Wendlandvi function located at the zthknot as depicted in Fig.2A, x is the coordinate being evaluated, is the distance of x from the knot center (xi) as a function of support radius (R) (10):

[0113] (10)

[0114] Fig. 2A illustrates the zthknot; and Fig. 9B illustrates evaluation of z'thWendland function at x.

[0115] and (*)+ is described as a truncated power function, which allows knot z to contribute to the level set value of point x. As mathematically expressed, and pictorially depicted in Fig. 2B), the CS-RBF attenuates to zero at the compact support radius, R.

[0116] Now that level sets and CS-RBFs are understood, the notion of compiling CS-RBFs to compute a level set will be described. Weighted superposition is employed in (1 1):

[0117] (11)

[0118] where cp is a scalar radial basis function, t is a pseudo time accounting for the pseudo-design iteration velocity of the design variable, i is the zthdesign knot index ranging from 1 to N, αi is the z-thdesign variable evolving across iterations where iterations are analogous to steps in the pseudo time, and is the spatially varying compactly supported Wendland-,.i function as defined in (9). Each blue dot in Fig. 2A and 9B represents a unique knot, and therefore a unique CS-RBF multiplier. A sample PLS surface is given in Fig. 1A and 8B. It is noteworthy that they are equally spaced. This is critical to the coming fixed grid expansion described herein.

[0119] The CS-RBF PLSM provides the method of evolving a CS-RBF to solve structural optimization problems. Through the multivariate calculus of variations, the terms α parameterizing φ enable perturbation towards minimization of functions of φ . The following derivation is adapted from

[0045] . Classically, the Hamilton- Jacobi Partial Differential Equation (H- J PDE) is used to express surface velocity as a function of φ, x, and t (12):

[0120] (12)

[0121] A common assumption is that only the normal velocity component, vn, contributes to boundary evolution. This allows the H-J PDE to be rewritten (13):

[0122] (13)

[0123] ω and α are functions of x and time t, respectively. This imposes separability upon φ , and hence the H-J PDE becomes (14):

[0124] (14)

[0125] Rearranging, vnis expressed as follows (15):

[0127] Wang et. al. goes on to explain how this indicates that the boundary velocity extends to the entire design domain as terms are evaluated at each knot. Because of this, they explain how “the original topological shape optimization problem is converted into a relatively easy size optimization problem, in which the expansion coefficients act as the design variables”. In the PLSM, additional schemes for velocity field extension are unnecessary. Furthermore, the H-J PDE can be used to cancel out time dependence rather than introduce a need to solve the PDE itself. Hence, the separability enabled by the RBF representation affords simplicity and reduction in implementation. With use of a fixed grid of knots, we describe a further simplification and reduction enabling 3D abstraction in the following.

[0128] Note how for a given t, <p is a spatial function of a summation of a single unitaryWendland function that is spatially relocated to knot centers and multiplied with the corresponding knot value. The value of theWendland function at the knot furthermore corresponds to thevalue of the / Wendland function at the zthknot (see , which illustrates), which illustrates spatiallyrelative Wendland function values.

[0129] This observation leads to the following expansion on CS-RBF PLS evaluation / collocation (16):

[0130]

[0131] where a (aka, F1{A(t)}) is a matrix in 2D or a 3D tensor in 3D of equally spaced knot values (i.e., a fixed grid of values), * is the convolution operator, and ω (x) (aka,is a single unitary (occurring at origin) Wendland function evaluated at any set ofpoints relative to the origin / unitary Wendland center. The convolution operation can be performed in a windowed fashion via the template approach or spectrally via the Fast Fourier Transform (FFT, F{ } in (16)). FFTs are used in this work, along with the inverse FFT (iFFT, F-1{ }). They can be computed extremely efficiently on an inexpensive GPU for very large variable sets whenusing fixed grids [46-48]. As such, <p can be computed efficiently at all Gauss points, and with sub-Gauss point accuracy, achieved as follows:

[0132] Produce a compactly supported elemental list of Gauss points relative to the origin.

[0133] Evaluate ω at these points, yielding a Wendlandvi Gauss FE integration point volume.

[0134] Convolve this FE integration point volume with a periodic version of the α field.

[0135] Compute nD (multidimensional) FFT of convolutional operand, here α, and convolutional fdter, here ω .

[0136] Elementally multiply the two together.

[0137] Compute inverse nD FFT of the result and remove extra padding to arrive at the solution to the convolution.

[0138] Essentially, we convolve two FE integration point volumes : one representing the spatially defined unitary Wendland vi function, ω (x), and one representing the set of expansion coefficients, α, which belong to knots that happen to exist at nodal coordinates (note: this is the key observation for this approach to work). What results is a sub-Gauss point accurate φ field. Similar applications of convolutional LS structural optimization have been documented in recent literature, but none have been utilized to perform 3D IH, or by GPU.

[0139] What follows is a description of the objective, structural analysis method, constraints, and formal Optimization Problem definition for the disclosed CS-RBF cPLSM IHP.

[0140] In the formulation of this IHP, the goal is a printable structure. Towards this goal, IH must be rapid. The means to rapid computation, aside from sensitivity analysis, have already been described as pertinent to the CS-RBF cPLSM. The next consideration is that of printability, and specifically, whether the structure has floating objects. Digital Light Synthesis (DLS) based AM is considered here on. Floating objects do not contribute to the structuralresponse. Additionally, they have a good chance of impeding proper printing by obstructing the incident platform. In minimum compliance solutions, there are typically no floaters. This may also occur in minimal volume / weight designs subject to structural constraints. As the following must be devoid of floaters, the objective will be minimization of volume. As described herein regarding constraints, because a target homogenized property is had, the percent error relative to said target homogenized property can be used to constrain the design. Relevant to the present section on objective formulation is the fact that a percent is used in the constraint. To avoid a mismatch in the effect of optimization argument magnitude, the volume is also expressed as the overall percentage of the unit cell occupied by material. Volume is computed by integrating H(φ ), where (p is the radial basis field and H is the Heaviside function, over the microscale domain. As such, the objective is computed as follows slightly deviating from that in Wang, (17):

[0141] (17)

[0142] Here, fv%is the volume fraction percent, / is the objective function, |γ| is the overall volume of the design domain, D implies bounds on the design domain, dΩ is an infinitesimal volumetric step / sample in the design domain, e is the element index, g is the local Gauss point index, Neis the number of elements in the system, Ngis the number of Gauss points in an element (i.e., using 5-point Gauss Quadrature is a Gauss weight, and xgisa fixed grid of coordinates for the gthGauss point occurring within an element. Note that a fixed grid is achieved by using cubic hexahedral finite elements.

[0143] As previously discussed, and for reasons described in the ensuing discussion on constraints, volume fraction percent is used because the constraints will be expressed as percentages since GCM residuals can take on values of varying magnitude. Finally, it is essential to note that the above can be computed with relative ease even if there are millions of elementsdue to the convolutional approach. Similarly, 1stand 2ndorder sensitivity analyses are made possible through the approach.

[0144] If using the Interior Point Method (IPM) for optimization, then 2ndorder sensitivity (Hessian) would have to be computed or approximated. For this to be done, and without approximation, the Lagrangian must be twice differentiable. Additionally, it would need to be computationally feasible, which it is not for large compact support radii, R, and high discretization fidelity per our own investigation. Alternatively, if an approximation approach such as the Method of Moving Asymptotes (MMA) is used, which inherently has a built-in diagonalized approximation to the Hessian of the Lagrangian, then the Lagrangian must only be singly differentiable, and the entire Hessian does not need to be computationally tractable. Optionally, second order information may be incorporated into MMA, yielding the Second Order MMA (SOMMA)

[0062] . Here, the diagonal of the Hessian per constraint and the objective may be utilized. In the disclosed implementation, these components are nearly as computationally costly to compute as the Jacobian information and computational cost is less subject to R. Before computing the Jacobian, an expression for the velocity of level set dependent functions must be derived relative to the level set parameters, rather than time (18):100145] 618)

[0146] where y is an arbitrary function of the level set. Here, (14) is used to relate the temporal derivative of y to the product of sensitivity of y with respect to the mihdesign variable and the temporal derivative of said design variable.

[0147] What follows applies not only to the objective sensitivity, but also to the constraint sensitivity. Taking the temporal velocity of the boundary, via (13)-( 15) and (17), (19) is derived:

[0148] (19)

[0149] whereis the Dirac Delta on φ , which yields the boundary.

[0150] Rearranging (19) via (18) yields the sensitivity of the objective with respect to the expansion coefficient design variable, independent of time such that it is no longer considered (20):

[0151] (20)

[0152] where amrepresents theexpansion coefficient design variable,is theWendland function occurring at thedesign knot, and is the sensitivity of thevolume fraction percentage with respect to theexpansion coefficient design variable.

[0153] Typically,) and are approximated via differentiable piecewisefunctionsand [60,63]. This is done to control minima and maxima, as well asbeing computationally simplistic with respect to approximates such as the sigmoid function. However, when computed via GPU in MATLAB, exponential operations within sigmoid approximation are highly efficient. The sigmoid function is useful in the disclosed context since it is continuously differentiable beyond the first derivative. To account for the percent of intermediate material comprising the ultimate design, the boundary thickness parameter, A, is adjusted by a design variable

[0154]

[0155] (22)

[0156] (23)

[0157] where the following must always be true for A(24):

[0158] (24)

[0159] In this custom sigmoidal formulation, β is a scaling factor for the LS, isthe minimal Heaviside value to correspond with the boundary ends (here, 0.025), y is a scaling factor for the pseudo boundary width (here, 2.7), andis bounded such that the region of pseudomaterial is neither infinitesimally thin (sensitivity would vanish) nor excessively broad (optimization would be prolonged). It is initializedsuch that the boundary contains at least a Gauss point. Here, it is assumed that 3x3x3 Gauss Quadrature is used. In implementation, the approximates are not allowed to go to zero, preventing numerical instability. The sensitivity of the Heaviside with respect to the set a is computed via convolution and the product rule (25):

[0160]

[0161] where *xrepresents spatial convolution. Substituting the summation over all elements for convolution over the Gauss point fixed grids within FEMs of the Dirac Delta and a unitary Wendland function, the sensitivity for all design variables may be simultaneously computed using the convolution operation over the element dimensions (*e), as follows (26):

[0162]

[0163] where provides the total Jacobian of the objective. This exemplifiesthe power of the current approach. Via GPU accelerated FFT-based convolution, computation on a half a million elements with ample Gauss Quadrature and a generous compact support radius compute in under a tenth of a second (often less than a hundredth). Instead of performing aweighted summation of all Gauss points over all elements per each design variable, only a weighted summation of convolutions of all design points over all Gauss points must be performed.

[0164] The first-order sensitivity of (17) with respect to is (27):

[0165] (27)

[0166] For the total Hessian, a 3 -fold convolution is needed. In SOMMA, since only the Hessian diagonal is needed from the total Hessian, simplifications are used for further acceleration and memory reduction. This simplification avoids the 3 -fold convolution.

[0167] Per typical computational intensity, second derivative information is seldom included in LS optimization. There are situations in which second order information helps, and vice versa. As it has not been explored in 3D for the LSM, it is explored herein. In the SOMMA approach, second-order information is incorporated via the diagonal of the Hessian as per separability. The Hessian diagonal relating to knot expansion coefficients is computed as follows (derived below)(28):

[0168]

[0169] The total Hessian diagonal is computed as follows (29):

[0170] (29)

[0171] As noted herein, to compute a target effective stiffness tensor, , via theinverse homogenization problem, several direct homogenization analyses must be conducted. At each step of the optimization algorithm, a homogenized constitutive tensor, , must beevaluated for the current configuration of the micro-scale domain, Y. Such a task may be modeled as follows (30):

[0172]

[0173] where cpqrsis the constitutive tensor that varies along are prescribedeffective strains, and are unknown strain fluctuation fields that depict localized corrections todue to heterogeneity. Notice the superscript indices written above strain fields, like (ij) and (kl). These indices identify six different load cases that must be considered for the computation of all coefficients in . Namely, three cases of prescribed effective normal strain, and threeother cases of prescribed effective shear strain, covering all directions. For each of these scenarios, is found by solving the weak-form governing equations (31):

[0174] (31)

[0175] Notice that a weight field, wt, was introduced. It arises from the weighted residual method, is defined to be non-null, and can be associated with a virtual displacement field.

[0176] As is prescribed, the right-hand side of (31) is known and serves as forcingconditions. Further, following the assumption of linearity in displacement-strain relations, we can rewrite the problem to be solved for a displacement field, uL. as done in (32):

[0177] Ui .

[0178] By employing the Finite Element Method (FEM), it is assumed that, within each element (voxel, a discrete portion of Y), tq can be locally modeled as an interpolation of discrete displacement values. That is, within the domain of an element,we admit = Nede, where Neis a matrix of shape functions, and deis a vector of nodal displacements, both related to theelement. Derivatives of these shape functions can be employed to compute a matrix Besuch thatwithin Ye. By applying these considerations to (32), we arrive at(33), which yields an algebraic equation that involves large sparse symmetric matrix K, commonly called the stiffness matrix.

[0179]

[0180] where A is an assembly operator, 8 is the set of all elements, Ceis a 6x6 matrix representation of the fourth-rank constitutive tensor at each element, d is a global vector of nodal displacements, and f is a global vector of nodal forces. Notice that w is also a global vector, and it is known to be non-null, so the system of linear algebraic equations in (34) must hold true.

[0181] (34) .

[0182] After solving the system shown above for each load case, al can berecovered from each\ so that may be computed as in (30).

[0183] By means of FEA, one can think of the problem of numerical homogenization of elasticity as a computational procedure: For each of the six load cases, assemble a right-hand side vector and solve a large system of equations; then, compute the volume averaging integral in (30) to obtain all coefficients in c^ki. Among these steps, solving the systems can be, by far, the most significant task in terms of time expenditure, and it can also be particularly demanding in memory requirements. In conventional scenarios, the system of equations is the bottleneck of these analyses. When dealing with structured meshes, as images or regular grids, matrix-free implementations of Krylov subspace solvers in GPU are a promising approach to mitigate both the time and memory issues. Massive parallelization plays a key role in code acceleration, and matrix-free schemes allow for the global K to not be stored. This is done by substituting the sparse matrix-vector products within each Conjugate Gradient iteration, for example, by a convolutionlike element-by-element or node-by-node operation.

[0184] Tn this work, heterogeneity is assumed to be allowed within each voxel, and so the integrals that yield element-wise local stiffness matrices are numerically computed with 5d(53=125) Gauss points. In this regard, we cannot consider just a few possible configurations of local matrices, i.e. all of them must be either stored or recomputed on demand. To recompute the local integrals every time an element’s matrix is required would be possible, but impractical for simulations in our scope, as the computational cost of each matrix-vector product within the system solver would be significantly heightened. In that sense, we chose to store a list of local matrices. Notice, however, that storing all local matrices adds up to more memory requirements than that of storing only the non- zero coefficients of the global matrix, K. (Note: Every local matrix is 24 x 24 (576 coefficients), while there are 81 x 3 = 243 non-zero coefficients per node in the global matrix. In periodic meshes, # elements = # nodes.) As structured meshes are considered, we can improve on that. The idea is to store a list of non-zero local coefficients in terms of stencils instead of elements, as one would normally do in Finite Difference schemes, for example. Knowing that the disclosed meshes are periodic, that each node of the mesh has 27 neighboring nodes, and that each node has 3 degrees of freedom (DoFs), we know that there must be 81 nonzero coefficients per DoF (per line of the global matrix). Their positions within the global matrix can be quickly inferred via simple grid numbering rules. In a nutshell, we can store perfectly coalesced arrays of stencil coefficients per DoF, as illustrated in Fig. 3, which not only lowers memory allocation from conventional sparse data structures, but also speeds up the SIMT computations in GPU. Preliminary tests with an NVIDIA GeForce 3080 Ti graphics card indicate that matrix- vector products with this approach can be sped up by a factor of up to lOx in comparison to working with a list of element matrices.

[0185] Fig. 3 illustrates a simplified illustration of the proposed approach to store DoF- by-DoF non-zero coefficients. To the left is a plot of a 10x10 matrix that could represent a IDperiodic FEA problem. Each asterisk represents a non-zero coefficient. To the right, a depiction of how these coefficients would be stored in memory, as a contiguous array of coalesced data.

[0186] The effective homogenized elastic properties of the optimal solution do not need to exactly match the target values. If all the elastic components of the generated microstructures are off by a particular percentage such as being 5% lower than expected, the gross elastic response of the macrostructure is anticipated to only be off by at most 5% relative to theoretical response values. There will most certainly be a discrepancy between the mechanics of actual bone and what can be approximated by FEA homogenization and the mechanics of the fabrication material. However, the approach should in theory best recreate the anisotropy of the modeled structure. To that end, the constraints are as follows. Each unique component of the homogenized GCM is targeted. Such unique components are displayed in (35):

[0187]

[0188] The number of unique components is given by ):

[0189] (36):

[0190] (36)

[0191] Where gives the number of unique GCM components to analyze whenperforming full anisotropic 1H, which is also given by the total number of combinations possible in 2 samples out of (d + 1) objects, d is given as 6 in 3D. As such21. While there are 21 unique components of the target GCM, in practice as per the symmetry that produces this, there are 36 factors descriptive of the purely anisotropic homogenization. As such, all 36 components of the GCM will be constrained. Because 15 are repeated, only 21 need to be computed at each iteration with the 15 symmetric components being doubly weighted to account for recurrence.Because the GCM is a matrix, a two subscript convention will be used to denote row and column, respectively.

[0192] To compute the percent difference relative to the critical value, the followingPercent Error (PE) is computed (37):

[0193] (37)

[0194] where^ is the percent error of the GCM constituent in the row andcolumn (now on referred to as th component) as a function of the design variables. Whileis given as the component of the target GCM, (35), is computed as follows(38):

[0195] (38)

[0196] Bearing resemblance to (6), the component summed across e can be considered analogous to in (6). Like in (6) is the applied strain field (here a 6x1 vector). Thisis constant in all elements and at all Gauss points. Here, it is therow or column of the identitymatrix. represents the domain of the element. Cespans this domain. is the straindisplacement matrix, which is the same for all elements but varies at Gauss points. The approach to, and MATLAB code for, its computation for hexahedral elements is inspired by

[0056] . It is a 6x24 matrix in the integral notation and a Ngx6x24 tensor in the summation notation where the 1stand 2nddimensions are not indexed into, respectively. Again, there are 24 DoFs (d) in each hexahedral element.Xeis the displacement field in the element, which is a 24x1 vector. Thisis solved by imposing periodic boundary conditions and solving the response to the applied strain field as previously discussed and implemented in both in-house code. Herej is a placeholder for the residual strain -residual stress product. C represents the d x d (6x6) bulk materialelastic constitutive matrix. In this work, it is isotropic and has the same value at all elements andGauss points. It is modulated by varying a, which in turn modulates the presence of material spatially via the Heaviside scalar, While not explored here, it is possible to have afunction which interpolates between many materials, relevant to multi-material manufacturing.

[0197] Via a straightforward derivation, the sensitivity of thePE with respect to the zw* expansion coefficient is given by (39):

[0198] (39)

[0199] The sensitivity of theGCM constituent with respect to the m*11design variable, is given by (40):

[0200] (40)

[0201] which may be simplified via convolution and rearrangement, as done in (26), via the approach of convolving over all elements in the domain to obtain a direct expression for the Jacobian of the GCM constituent with respect to a (41):

[0202]

[0203] After convolution at the Gauss point domain-wide fixed-grid over allelements, numerical integration is performed via weighted sum over the Gaussian weights. The diagonal of the Hessian with respect to a is given by (42):

[0204] (42)

[0205] (41) and (42) are used to solve for the PE Jacobian and Hessian diagonal of the GCM component with respect to a efficiently, respectively (43) and (44):

[0206] (43)

[0207] (44)

[0208] Finally, constraints may be formulated for use within the optimization approach discussed in the following. Constraints on the component of PE are as follows (45):

[0209] (45)

[0210] Referring to (17)-(29), constraints on volume fraction percentage are as follows(46):

[0211] (46)

[0212] The minimum and maximum allowable design variable values at each iteration are as follows (47):

[0213] (47)

[0214] When the Heaviside approximation parameters are considered as optimization design variables, they must be accounted for in problem Jacobians. The sensitivity of PE with respect to Δ is computed as follows (48):

[0215] (48)

[0216] Via (23) and (40), the first and second order sensitivity of the constitutive components with respect toare (49) and (50):

[0217] (49)

[0218] (50)

[0219] The total Hessian diagonal for P is computed as follows (51):

[0220] (51)

[0221] The design variables driving the approximation to the Heaviside function are updated to minimize the Percent Total of Intermediate material, hence referred to as PTI. PTI may be formulated as (52):

[0222] (52)

[0223] where represents the spatially variate material density field. The totalvolume term cancels out from the integration.

[0224] (52) is discontinuous. Thus, it is not differentiable as given, and cannot be used as a constraint on its driving variables, α andTo circumvent this issue, a differentiable approximation is used, subject to a static sigmoid function (53):

[0225] (53)

[0226] where cr is the sigmoid function, p is the independent variable, a is a parameter for the width of the transition in the sigmoid function (here, 1.75), and c is a shift in the sigmoid function along p. (53) is used as follows to approximate (52) in (54):

[0227] (54)

[0228] where is the differentiable sigmoidal approximation to

[0229] In practice, (54) is constrained since it is differentiable, rather than (53) which is still used in early stop and tracking since it is most accurate. The following constraint is imposed on the design (55):

[0230] (55)

[0231] Differentiation is carried out via the quotient rule by breaking up the problem into multiple sub-problems (56):

[0232] (56)

[0233] where and represent the numerator and denominator,respectively, of the expression posed in (54). This allows for utilization of the quotient rule (57) (foregoing for brevity):

[0234] (57)

[0235]

[0236] (58)

[0237] <59)

[0238] is given as components of (25), and the lower and upper sigmoidfunctions, and(both relevant to (54)), are as follows (60)-(61):

[0239] (60)

[0240] (61)

[0241] and their sensitivity with respect to the Heaviside are as follows, via differentiation properties of the sigmoid function (62)-(63):

[0242] (62)

[0243] (63)

[0244] completing the set of equations required to compute (57).

[0245] A convolutional abstraction is implemented for efficient computation, derived as follows. First, the Jacobians can be solved (64)-(65):mnM71

[0248] where is given by (25), leading to the given convolutional expansions. (64)and (65) together enable the following approach to efficient computation of the Jacobian withrespect to the design variables α (66):

[0249] (66)

[0250] The most important aspect of PTI is the dependence onThe parameter directly controls the amount of intermediate material in the design. Updates in are dependent on havingthe gradient inwith respect to , as given by (67):

[0251] (67)

[0252] where

[0253] (68)

[0254] (69)

[0255] As such, (66) and (67) give the equations for the sensitivity of the intermediate material, (54), a differentiable approximation of (52), and its constraint, (55), with respect to the design variables α and , respectively. (26) and (27) give the sensitivity of the volume percentobjective, (17), and its constraints, (46), with respect to the design variables α and, respectively. Finally, (43) and (48) give the sensitivity of the constitutive percent error, (37), and its constraints, (45), with respect to the design variables α andrespectively.

[0256] If SOMMA is used, the diagonal of the totalHessian is needed. First, theMATLAB Symbolic Math Toolbox is used to take 2ndderivatives of (60) and (61) with respect to Then, the following are computed specifically under assumptions that they are for the Hessian diagonal with respect to α (i.e., n=m), otherwise further symbolic differentiation is needed to account for mixed derivatives (70)-(71):

[0257] (70)

[0258] (71 )

[0259] (70) and (71) can be evaluated for all α simultaneously via convolution (72)-(73):

[0260] (72)

[0261]

[0262] These equations are incorporated to find the Hessian diagonal of contributedby a (74):

[0263]

[0264] The second derivative of N and D with respect toare computed as follows(75)-

[0265] (76):

[0266] (75)

[0267] (76)

[0268] where the functional expression for is computed via the MATLAB SymbolicMath Toolbox.

[0269] The second derivative ofwith respect to is computed as follows (77):

[0270] (77)

[0271] The diagonal of the total Hessian of is given as (78):

[0272] (78)dA

[0273] With all objective and constraint evaluation and sensitivity expressions developed, prior to providing the disclosed structural optimization problem definition, we will briefly discuss the disclosed early termination criteria, accuracy metrics, and dynamic design variable move limits.

[0274] Herein, the preliminary derivations required for approximate intermediate material presence control in structural optimization are derived. Once intermediate material is totally removed from the system, so are all sensitivities (neglecting fictitious boundary values which serve to prevent numerical error). Once the boundary disappears, so does capacity for design domain evolution. Experimentation suggests that the optimization algorithm first best satisfies the constitutive targets, and then removes the intermediate material in the design. Hence, once a critical threshold in PTI is achieved, we find that any further iterations introduce great numerical instability. We propose a simple early termination criterion (79):

[0275] (79)

[0276] Since the target is less than the critical value, eventual satisfaction is typical even if constitutive constraints are not satisfied.

[0277] The following provides inverse homogenization metrics of harmonic MeanAbsolute Percent Error (hMAPE) (80) and symmetric Mean Absolute Percent Error (sMAPE) (81):

[0278] (80)

[0279] (81)

[0280] where d is the side length of the GCM as used earlier herein.

[0281] hMAPE skews the average error estimate towards small values of PE|;, which is important because large magnitude PE outliers present themselves. sMAPE provides a summary metric for PE that weighs residuals based on their magnitude. This is achieved by foregoing PE normalization and summing absolute residuals in the GCM. Normalization takes shape in the denominator, where the average magnitude of the GCM components is computed. Due to the sum of absolute residuals in the numerator, and the average magnitude in the denominator, symmetry exists, hence the term symmetric in sMAPE.

[0282] We propose that both hMAPE and sMAPE are essential in THP performance analysis. It is well known that skeletal microstructure homogenization often yields GCMs with small coupling terms, where coupling terms are those preserved in purely anisotropic / triclinic GCMs, not present in even orthotropic / orthorhombic GCMs

[0066] . These terms describe the coupling between extension and shear. As such, only triclinic materials can fully couple dilatational (expansion-contraction) and distortional (shearing) deformation. In a 6x6 GCM, these terms are those outside of the top left 3x3, and the lower right diagonal. Note how in (35), they are particularly small for the given bone specimen. However, C16and C36(hence, also C61and C63) are relatively large. As non-negligible coupling terms illustrate a departure from elastic symmetry, reasoning for the visual lack of symmetry is provided, which is the directional elasticity mapping for the GCM given in (35).

[0283] Typically, bone is assumed to abide by orthotropic elasticity. However, this example clearly illustrates that bone can present non-negligible mechanical coupling between dilatational and distortional deformation. sMAPE foregoes the effects of low magnitude GCM residuals. This ensures that the most dominant components of the GCM, whether they be coupled or not, are accurately targeted. hMAPE removes outliers from the estimate of GCM difference, ensuring that there are dominantly low percent errors between GCMs. The previously described PEj;- constraints ensure reductions in all GCM residuals proportional to the maximum target magnitude. Hence, by the point in which sMAPE and hMAPE are sufficiently low, it is probable that PE(7is generally reduced as well. To ensure that the design has undergone significant optimization towards its practical limitations, PTI must be critically low as well. As such, we tend to observe that hMAPE and sMAPE are reduced, and then PTI is reduced. We explore this in the disclosed results. We also assess performance by the number of design iterations, overall duration, and duration per iteration.

[0284] We chose to employ dynamic design variable move limits. The caption of Fig.4 provides a description of the approach. Fig. 4 illustrates a ssigmoid adaptive move limit diagram. The upper move limit is 0.25 and the lower move limit is 0.10. At ~6 iterations, it starts a sigmoid decline from the maximum and reaches the minimum around 26 iterations smoothly. The design variables are always bounded to ±1.00, adding perspective to the given move limits within the context of structural optimization herein.

[0285] It is a sigmoidal decline from the initial maximum move limit to a reduced move limit. It smoothly declines from 0.25 at 6 iterations to 0.10 at 26 iterations. The curve follows (82):

[0286] (82)

[0287] where |d(d. v. ) |maxrepresents the maximum allowable absolute change in the design variables, a and A if the Heaviside parameter is optimized, and iter, is the iteration number.We find that this helps promote initial rapid evolution of the design, and subsequent gentle evolution, avoiding unsteady designs alterations (please see associated video of design evolution).

[0288] Note thatwas limited to the range (—1, 1] as presented in (21) to ensure that the move limit applied to a is valid foras well, in addition to ensuring change in value is similar in order of magnitude across design variables.

[0289] With the objective, constraints, sensitivity analysis, and termination criteria defined, we may define the optimization problem.

[0290] The following is defined to appropriately constrain constitutive components following normalization

[0291] (83):

[0292] (83)

[0293] Combining equations, iterative design optimization is performed by applying theMMA approximation

[0067] , and separately by applying the SOMMA approximation

[0062] , to the following constrained minimization problem (84):

[0294] (84)

[0295] Note the absence of KKT metrics (i.e., They are standardoutput in the MATLAB implementation of such MM A algorithms

[0067] . We do not use KKT based metrics because the problem is established in a general sense for various kinds of IHPs, featuring problem specific metrics.

[0296] The MMA algorithm serves to yield a separable approximate to a nonlinear constrained optimization problem, here the IHP (84). This explains the capacity to use the diagonal of the Hessian in SOMMA. MMA and SOMMA use the same coefficient sets in approximation. The disclosed values are as follows, (85):

[0297]

[0298] Svanberg’s MATLAB implementation of the MM A approximation comes with a solver for the approximate subproblem, subsolv, in addition to formulating the approximate subproblem, mmasub. We developed a GPU accelerated variant of the codes, still fully in MATLAB. As there are more design variables than constraints, the disclosed subsolv algorithm employs the direct method rather than the adjoint method. We removed the extremely cumbersome matrix-(spdiags matrix) products and replaced them with equivalent matrix-vector element-wise multiplication. This greatly improved speed and efficiency and yields GPU compatibility.

[0299] We opted for less refinement of the solution. Within the subsolv algorithm, instead of enabling 200 max iterations and 50 max refinements per iteration of the approximate solution, we used 14 and 10, respectively. This allowed us to quickly find solutions involving the analytic gradients. We likely had to call mmasub more often than if we enabled fine solution prior to incorporation. However, this equalized the time spent in the MMA update and the time spent in simulation. Additionally, designs still converged rapidly. With this, we proceed to discuss the computing systems used and experimentation.

[0300] As shown in (40), above, the sensitivity with respect to an expansion coefficient, am, of a function that is linear with respect to the Heaviside of the LSF,, reduces to a multiplication of the non- -dependent parts with the volumetric Dirac Delta of the LSF,, and the unitary Wendland function, relocated to the center point of the':knot.It is shown in (41 ) that to can be convolved over the volume yielded by the product of the derivativeof and a linear field (in (41), the strain energy density) to simultaneously solve (40)for all m.

[0301] Herein, as in (38), use the following form for such a function (86):

[0302] (86)

[0303] where f is an arbitrary linearly integrated metric dependent on a spatial field of values, y, and

[0304] The sensitivity ofwith respect to amis (87):

[0305] (87)

[0306] and using (23), the sensitivity of / with respect to A is (88):

[0307]

[0308] Above is detailed howis simply an elementally relocated variant of m, since the spatial element distribution coincides with the spatial distribution of knots. This is used to perform the simultaneous solution to (87) for all m as follows

[0309] (89)

[0310]

[0311] where [ ] indicates a column vector.

[0312] (89) provides the part of the Jacobian not contributed by A, which has its own gradient formulation separate from that for α.

[0313] The Hessian components with respect to α may be considered as providing the sensitivity with respect to theexpansion coefficient of the sensitivity with respect to theexpansion coefficient. It is represented as such (90):

[0314] (90)

[0315] where { } indicates a row vector and implies that the contents are to be transposed, symbolizes the derivative of the approximated Dirac Delta function.For flexibility in testing alternate Heaviside approximations, the derivative of the Dirac Delta function is computed symbolically in MATLAB via the Symbolic Math Toolbox, along with the first derivative and derivatives with respect toHence, the lengthy equation is omitted herein.

[0316] represents the A-fold convolution. Here, N=3. This is not just two sequential convolutions, as is easily conceptualized. It is the very rare simultaneous convolution of two kernels independently on top of each other (contained right of the operator) and the convolved volume / signal (left of the operator). We have developed two GPU accelerated implementations in MATLAB (one template based, the other FFT based). The result is the influence that one knot has on the domain, and then the overlap in sensitivity attributed by the surrounding compactly supported knots. It is a tensor of rank-(N-1).

[0317] SOMMA requires the diagonal of the Hessian. This is represented in red within(90). It does not require N-fold convolution. Since m=n, the diagonal of the Hessian with respect to a reduces to (91):

[0318] (91)

[0319] which is approximately as easy to compute as any of the Jacobians presented like

[0320] (89).

[0321] The total Hessian includes components with respect toThis is computed as follows (92):

[0322] (92), . . . .. . . . . , , . ,

[0323] where depends on symbolic derivatives with respect towhich areprimarily carried out in MATLAB’s Symbolic Math Toolbox. Note that this is not the case for all metrics, such a where the second derivative computations are more involved, as discussedherein.

[0324] In the diagonalized Hessian, this is appended to the bottom of (91 ) as follows to create a column vector with length m+1 (93):

[0325] Evaluation Of The Above-Described Approach On Different Systems, Target Categories, And IHP Approximations (First Vs. Second Order MMA), And Validation Of The Approach

[0326] The program is ubiquitous in that it has optional mixed-second derivative utilization and can run on a system with a sufficiently modern NVIDIA GPU and a corresponding MATLAB version. To assess generality, the algorithms were tested on two systems: a custom desktop PC featuring a 24 GB vRAM NVIDIA GeForce RTX 3090, to be referred to as 3090; and a hybrid tablet, the 15” Microsoft Surface Book 2 with the default NVIDIA GeForce GTX 1060 with 6 GB vRAM, to be referred to as 1060. 3090 is representative of a relatively high-end computing system, while 1060 is among the most affordable computing options available. Note that GPU memory sharing was disabled on both systems for efficiency.

[0327] The goal of the disclosed statistical comparisons will be to assess the impacts of anisotropy type (Triclinic vs. Cubic), system (3090 vs. 1060), and use of mixed-second derivatives (MMA vs. SOMMA) on the performance metrics presented herein. This fall into two categories: percent statistics (PTI, hMAPE, and sMAPE) and time (iterations, duration, and duration per iteration). Tests involve 30 replicates per each of the 8 possible scenarios. The a field is initialized based on the Schoen Gyroid, a triply periodic minimal surface which partially inspired this work. Note the random noise applied to the a field. All IHP solutions in this performance test are done with DPR 10 being the design material, having an elastic modulus of 1.8 GPa and an idealized Poisson’ s ratio of 0.35. Materials in this work have a 0.35 idealized Poisson’ s ratio, per algorithm sensitivity.

[0328] For each of the metrics, the following approach to statistical analysis is performed in MATLAB :

[0329] Heteroscedasticity test on all groups

[0330] If heteroscedastic, perform a tied ranks transformation across the whole set

[0331] 3-way ANOVA for variable primary effects

[0332] If variable(s) do not significantly influence the metric, pool metric data per the variable(s), avoiding effect assessment

[0333] With respect to variables that significantly influence the metric, perform an N- way ANOVA with interactions enabled, followed by Tukey HSD multiple comparisons and Bonferroni correction

[0334] If the metric was transformed in Step 1 , it still is in this step

[0335] A significance level of 0.05 is used in all tests. Exemplary iteration histories and topologies are to be provided as result.

[0336] FEA based computational homogenization has been used to reduce the complexity of large microstructural models, such as that shown in Fig. 13A. Fig. 13A illustrates an 82 μm resolution HR-pQCT; Figs. 13B and 13 C show 4 mm trabecular patch (13B) and cubic sample (13C); FIG. 13D shows a homogenized elasticity field (via nTop 4.18.2 ). Note how the elasticity map is aligned with the path of principle compressive stress, supported by the GCM (Fig. 13E). A large model can be broken down into constituent samples (as shown in Figs. 13B-13C, which are then reconstituted into a drastically simplified model. To achieve this, the key is an effective approach to continuum averaging, which is the homogenization process. The mixture of materials, here marrow and bone, can be represented as a single finite element via the homogenization approach. In the case of the microstructure presented in Fig. 13C, the generalized constitutive matrix (GCM) in Fig. 13E yielding the elasticity map in Fig. 13D defines the material properties of an equivalent single finite element. The structure embodied by these reducedelements enables tractable large-scale simulation, at the cost of possible load path mismatches and associated inaccuracies.

[0337] We ran a low target IH, with a target of 0.5% instead of 7.5%. Forvalidation, we assigned a 10 GPa elastic modulus to the design material, while maintaining the 5 GPa elastic modulus in the target from Fig. 13. This is done using SOMMA and comparatively MMA on the 3090 system. We validated efficacy by comparing both skeletal sample and designed topology homogenized elasticity using the commercial software, nTop (v5.4.2).

[0338] Figs. 14A-14F present the outcomes of the statistical analysis described herein.All statistics presented heteroscedasticity, and as such, medians (short horizontal line within shaded area in notch plots) are those whose groupwise differences are statistically compared. Figs, 14A-14F show notched box plots and indicators of statistically significant trends in assessment metrics. Red dot-dashed lines indicate that the respective group has a median statistically different from that in all other groups depicted, p < 0.05. Many significant trends are effectively depicted in this way. Black dashed lines act to separate MMA plots from SOMMA plots. Where orange and blue notch plots are both depicted, comparisons between the bone and the gyroid influence the metric. Otherwise (see 14C), anisotropy type did not have statistically significant effects, and both are pooled under blue. Bars with asterisks (see 14D) act to display group specific statistically insightful comparisons. ***, the only asterisks set depicted, represents

[0339] Fig. 14A, depicting IH duration, indicates that all groups present significantly different medians with respect to one another, except for the 3090-SOMMA Cubic and Triclinic anisotropy groupings with a pooled median duration of 32.21 seconds The 3090-MMA-Triclinic group had the lowest median duration, at 14.89 seconds.

[0340] System did not influence number of iterations, PTI, hMAPE, or sMAPE (p =0.751, 0.163, 0.543, and 0.868, respectively).

[0341] Tn iterations (Fig. 14B), all groups had significantly distinct medians (p < 0.05).At 7 iterations, the Triclinic-MMA group had the lowest median iteration count (p < 0.05). The Triclinic-SOMMA group had the most at 14 iterations (p < 0.05).

[0342] The SOMMA-Triclinic group had significantly higher median PTI than any of the other groupings presented at 0.735% (p < 0.05, Fig. 14D). All groups had median PTI less than 1%, indicating that the termination criteria were consistently satisfied.

[0343] The Triclinic groups had significantly lower hMAPE than either of the Cubic groups (p < 0.05, Fig. 14E). For the Triclinic materials, the SOMMA group had a significantly lower hMAPE than the MMA group (p < 0.05), with the SOMMA and MMA groups at 3.50% and 5.43%, respectively. The two Cubic groups had similar medians, with a pooled value of 13.31% (P « 1).

[0344] The MMA algorithm was subject to higher groupwise median sMAPE than that of SOMMA (p < 0.05, Fig. 14F), with the Cubic and Triclinic groups having respective median values of 9.27% and 15.86%. With SOMMA, the Cubic and Triclinic groups had similar medians (p = 0.307), pooled at 4.75%, which is within the desired 7.5%.

[0345] In computing the average duration per iteration, target class timing variability is mostly removed. Fig. 14C displays how anisotropy type does not influence duration per iteration (p = 0.818). As such, Triclinic and Cubic results per this metric are pooled. Stratified by system and algorithm, all groups presented statistically distinct medians (p < 0.05). The 1060-SOMMA group presents the greatest median duration per iteration, at 9.28 seconds / iteration (p < 0.05), while the 3090-MMA group presented the lowest median duration per iteration at 2.07 seconds / iteration (p < 0.05). Interestingly, the 3090-SOMMA group presented a significantly higher median duration per iteration than the 3090-MMA group, but only by 0.32 seconds per iteration (p < 0.05).

[0346] Meanwhile, the second order adaptation on the 1060 yielded a 2.42 second per iteration increase (p < 0.05). The 3090 appears to more effectively distribute the same workload,which is likely a matter of memory bandwidth differences (936 GB / s on the 3090 vs. 192 GB / s on the 1060) and the number of CUDA cores (10,496 on the 3090 vs. 1,280 on the 1060).

[0347] Per the above, the SOMMA approach preserves mechanical accuracy despite differences in anisotropy class. The 3090 system preserves mechanical accuracy while reducing overall duration. We proceed to validation, as described herein. Figs. 15A-15G and Figs. 16A-16G display exemplary runs respectively with MMA and SOMMA. The same settings and targets are used as previously described herein. However, to further test the two approaches in their ability to significantly alter the initial design, we used an overly rigid material inspired by RigidlOk (10 GPa elastic modulus). This contrasts with the 1.8 GPa material modeled herein. Figs.l5A-15G show MMA 3D IHP solution via CS-RBF cPLSM in overly rigid design material, where 15A shows a final design; 15B shows |PE| matrix bars; 15C shows IHP iteration history; 15D shows XY projection of 15E, 3D elasticity field, derived from GCM; 15F shows XZ projection; and 15G shows YZ projection. ~18 seconds without visuals.

[0348] Figs. 16A-16G show SOMMA 3D IHP solution via CS-RBF cPLSM in overly rigid design material. Fig. 16A shows a final design; FIG. 16B shows |PE|,-ymatrix bars; FIG. 16C shows IHP iteration history; FIG. 16D shows XY projection of Fig. 16E, 3D elasticity field, derived from GCM; FIG. 16F shows XZ projection; and FIG. 16G shows YZ projection. ~43 seconds without visuals.

[0349] The two final designs following MMA and SOMMA (Fig. 15A and Fig. 16A, respectively) appear similar. However, |PE| matrix bars (Fig. 15B and Fig. 16B depict stark differences in mechanical accuracy. The SOMMA algorithm yielded a more mechanically accurate design. This is further supported when analyzing the objective and constraint iteration histories, Fig. 15C and Fig. 16C. MMA appears to quickly reduce PTI, whereas SOMMA first reduces the mechanical constraints, increasing PTI, and then once mechanical constraints are satisfied, PTI is reduced. This follows the trends inferred from the results of batched analysis on the 1.8 GPa elasticmodulus material. It appears that SOMMA requires more iterations than MMA, but that it better meets the desired constraints, and in the fashion desired. Figs. 15D-15G and Fig. 16D-16G display corresponding results, where SOMMA better fits the desired continuum mechanics in this scenario. Note that the gridded and solid elasticity surfaces represent the directional elasticity for the design and target, respectively. Fig. 17 displays this via nTop, serving as validation. sMAPE was 5.70%, hMAPE on normalized PE was 0.04%, standard hMAPE was 0.96%, MAPE on normalized PE was 0.91%, and standard MAPE was 37.88%. Clearly, the disclosed PE normalization approach is effective in removing negligible factors which influence standard PE without visually altering the mechanics. Fig. 18 shows the design is periodic.

[0350] Based on modem GPU accelerated homogenization, a 3D LS IH on consumer grade hardware is disclosed herein. We then provide background on the potential for 3D LS IH in skeletal microstructure inverse design and discloses that it is a component for generalized inverse macrostructural design. Subsequently, we highlight the efficacy of convolutional approaches in structural optimization, proposing the disclosed CS-RBF cPLSM approach. We state that this avoids the primary sources of LSM computational cost.

[0351] Analysis identifies anisotropy class and second order information as key factors influencing mechanical accuracy, measures of duration, and final PTI. Validation suggests that the MMA algorithm exploits 10x weighting of PTI in (85) such that the PTI constraint is satisfied before the constitutive constraints, leading to early termination as intermediate material presence diminishes (Fig. Figs.15 and Fig. 16). While the batched analysis is not as difficult as the validation, since a lower elastic modulus material is used, these results reappear. Fig. 14F displays how second order information yields lower median sMAPE, consistent across anisotropy class. Fig. 14D displays a more conservative reduction in PTI when SOMMA is used, typically being closer to the stopping criteria of 1 than in MMA. Fewer iterations, per Fig. 14B, indicate a relatively simplistic trend in MMA IH.

[0352] This trend is best visualized in the rigorous validation for MMA, Fig. Figs. 15C, where PTI is promptly reduced. Meanwhile, SOMMA, as visualized in Fig. 16C, promptly increases PTI to increase the sensitivity domain. During this, the constitutive targets are satisfied. Once satisfied, PTI is reduced. SOMMA, while taking ~2x the iterations and 0.32 seconds longer per iteration, strikes a balance between PTI reduction and constitutive targeting not made in MMA when the initial mechanics are far from the target (Fig. Figs.15 and Fig. 16).

[0353] This study investigated the efficacy of the disclosed LS approach on different anisotropy types, systems, and if second order derivative information altered performance. A key finding is that the inclusion of the Hessian in MMA increases intermediate material presence for bolstered sensitivity, satisfies constitutive constraints, and then reduced the boundary until it was negligible. The standard MMA algorithm was not as effective. Additionally, the disclosed approach successfully extends to more than one anisotropy class. We also find that the disclosed3D LSM is effective in skeletal inverse homogenization.

[0354] We achieved 3D LS IH which benefitted from an FFT-based GPU-acceleratedLSM. Additionally, we incorporated the Hessian into the LSM. This study resolves the issue where the LSM had only truly been viable in 2D. The disclosed convolutional abstraction enabled derivations regarding the LSM that would otherwise be intractable to implement. The disclosed approach achieves anisotropy mimicry and yields smooth binary designs quickly and generally on consumer grade hardware. This study resolves the presented gaps in the field, as introduced primarily herein. These include high-level gaps such as 3D skeletal anisotropy mimicry and low- level gaps like Hessian inclusion in LS MMA implementation. Fig. 18 shows the design is periodic.

[0355] The approach resulting from this work would be a valuable addition to the repertoire of the modern computational modeler. Particularly, it could benefit those who utilize nTop. The modeling program has underpinnings like those of the presented approach, both beingbased on signed distance functions (SDFs). Advancements in engineered continua are at the cutting edge of automotive lightweighting. The approach presented here may be modified such that auxetic properties may be targeted, providing new value to those in the automotive engineering industry and in government such as at the National Highway Traffic Safety Administration. Finally, the approach may yield great value in skeletal biomechanics. The developed algorithm may be utilized to build up a gross scale skeletal surrogate from image data. This has the potential to influence human body model design, medical device testing, and implant design.

[0356] To recapitulate key findings, we have successfully developed the disclosed CSRBF cPLSM such that skeletal and generalized IH are feasible on consumer grade hardware. We integrated optimized pseudo-boundary adjustments to meet the problem defined, and the analytic Hessian to effectively generalize the problem to difficult cases. As negligible coupling terms may or may not present themselves, we formulated dually bounded maximum-absolute-target- component normalized constitutive residual constraints. This made the approach effective on different classes of anisotropic elasticity.3D Triclinic Inverse Homogenization Assessment Of Design Material Elasticity And SubFemoral Region

[0357] In another aspect, a 3D triclinic inverse homogenization assessment of design material elasticity and sub-femoral region is disclosed. This includes examples of high- throughput implementation on several cadaveric femoral scans, sampling multiple subfemoral regions and designing in a range of materials. Further disclosed is GPU convolution accelerated high-resolution segmentation.

[0358] Owing to its dynamic load-driven architecture, the skeletal microstructure is a rich source of mechanical variability. Recent literature has attempted to mimic anisotropy characteristics of the skeletal microstructure, for example via re-orientation of implicit derived structures and application of additive manufacturing (AM) considerations. Theseapproaches are generally limited to orthotropic anisotropy. However, without extra simplification, bone abides by pure anisotropy (aka, triclinic). For clinical translation, the development of a validated pure anisotropy 3D microstructure mimicry approach is essential. Considerations for successful generalized skeletal microstructure Inverse Homogenization (IH) includes design domain fidelity; sub-considerations of fabrication approach (i.e., intermediate material in optimized design); fabrication material properties; and skeletally relevant microstructural anisotropy variability. Described herein are developed curves that provide insight on material selection to those in the skeletal biomimicry community. This is accomplished by assessing the effects of elasticity on metrics describing the performance of the above pre-validated IH algorithm. Herein, the algorithm is used on the largest set of skeletal microstructures to be inversely homogenized in 3D. Additionally, triclinic properties are targeted.

[0359] Described herein are methods providing insight on: dataset acquisition; robust microstructure sample selection / preparation (GPU FFT convolutional masking approach, large scan semi-non-rigid atlas landmark / mesh-based scan registration, and distance field selection of nearest 4 mm cubes by trabecular volume fraction via convolution, for IH); implementation of, and adaptations on, the custom 3D IH approach; and statistical analysis (sub-femoral region analysis and design material elasticity analysis). Results are presented, which highlights the outcomes of IH on select materials as well as broad results and statistical outcomes. A discussion is included, which provides an explanation of algorithm and statistical outcomes. Finally, conclusions are presented.

[0360] The femur is chosen for its abundance of microstructural variability. This work aims to assess mimicry characteristics across different possible printing material elasticities and within anatomy. We conducted this work in two phases: (1) digital microstructure sampling; (2) microstructure inverse homogenization. Herein, a GPU accelerated approach toboth landmark registration and analytically solved constrained optimal microstructure sampling is disclosed. Also, modifications to the CS-RBF cPLSM 3D IH approach are described. The formal optimization problem definition is presented along with statistical comparisons to be made. To test the disclosed invention, performance variability across anatomic region and material elasticity are assessed, respectively.

[0361] All computations described herein were performed on a custom Windows 11 machine. Relevant hardware includes dual NVIDIA 3090 GPUs with 24 GB VRAM each, an AMD Ryzen 9 5950X CPU, and 128 GB RAM. The gaming series GPUs are not Tensor Compute Core (TCC) enabled, and do not provide gains via cross-GPU communication. Instead, sequential IH was performed concurrently across both GPUs. It is to be appreciated that other hardware can be used to practice aspects described herein

[0362] Described herein is the formation of a diverse skeletal microstructure dataset for subsequent inverse homogenization and ensuing analyses. Samples vary broadly in volume fraction such that the probable range in structural content is accounted for. These samples are from distinct skeletal regions, and a diverse (i.e., male+f emale) set of hosts.

[0363] A cadaveric femoral HR-pQCT dataset made public by Hazrati et al is used herein for its abundance of microstructural data and open availability [Hazrati-Marangalou J, Ito K, Cataldi M, Taddei F, Van Rietbergen B. A novel approach to estimate trabecular bone anisotropy using a database approach. J Biomech 2013;46:2356-62. https: / / doi.Org / 10.1016 / i.ibiomech.2013.07.042, incorporated by reference.]. The dataset presents a mix of HR-pQCT and pCT data. Because HR-pQCT data can be acquired from living human subjects, and to maintain consistency, the scans performed on the 1stGeneration XtremeCT HR-pQCT scanner (XtremeCT, Scanco Medical AG, Bruttisellen, Switzerland) are considered for this study. 24 scans with 82 μm isotropic resolution are considered, with donor data summarized by Table 1, below.Table 1. Donor metadata summary.

[0364] Fig. 19A shows a Volume slice; FIG. 19B shows FFT segmentation; Fig. 19C shows Femoral Head (FH), Neck Superior (NS), Neck Inferior (NI), Trochanter Superior (TS), Trochanter Middle (TM), and Trochanter Inferior (TI) sub-segmentation; Fig. 19D shows sample selection; and FIG. 19E shows the sample.

[0365] Sub-femoral region may influence microstructural morphology. Skeletal microstructure samples are selected from six regions of the femur: Femoral Head (FH), Neck Superior (NS), Neck Inferior (NI), Trochanteric Superior (TS), Trochanteric Middle (TM), and Trochanteric Inferior (TI). The following 3-stage process is used to automatically isolate specimens satisfying criteria described herein: (1) GPU based FFT convolutional cortical, trabecular, marrow, and exterior masking following the automated CPU based approach described in [Buie HR, Campbell GM, Klinck RJ, MacNeil JA, Boyd SK. Automatic segmentation of cortical and trabecular compartments based on a dual threshold technique for in vivo micro-CT bone analysis. Bone 2007;41 :505-15. https: / / doi.Org / 10.1016 / i.bone.2007.07.007, incorporated by reference], (2) Semi-non-rigid atlas landmark / mesh-based scan registration based on public MATLAB repositories, (3) Distance field selection of nearest 4 mm (~49 voxel side length) cubes with convolution based trabecular volume fraction over 20% and less than the subject regional 95thpercentile.

[0366] Step 1 is necessary because the cortical and trabecular regions have distinct properties. The trabecular-cortical interface is a-periodic orthogonal to the cortex. Because of this, the cortical region is typically rejected when homogenizing the trabecular region, and vice versa. There may even be volume slices consisting of material / air outside the bone if thesample is taken too close to the boundary of the cortical region. Tn contrast, valuable information presents itself when mechanics are within reason, even if air / cortical bone are present. We appropriately account for this in the disclosed statistical analysis by building on the segmentation workflow by Buie and colleagues, on GPU. Their approach makes use of segmentation and morphologic erosion, dilation, closing, smoothing, connectivity operations, etc. to automatically segment marrow, trabecular bone, cortical bone, and the exterior from an image given minimal user input threshold(s). We introduce convolutional morphologic operations, building on the step-and-save approach. Exploiting this property, a rapid GPUbased FFT segmentation algorithm was implemented using MATLAB R2022b. The result is presented in Table 1. With the rapid and efficient material segmentation algorithm, we proceed to Stage 2.

[0367] For each scan, independent of left vs. right, similar sub-femoral segmentations were needed. For consistency, an automated semi-non-rigid atlas landmark / mesh-based scan registration algorithm was developed, again featuring GPU acceleration in MATLAB R2022b. To achieve this, a single scan was segmented via the approach in Stage 1, and then the periosteal surface was meshed using available algorithms in MATLAB. A set of planes and lines were then defined via linear operations to partition the sub-femoral regions (FH, NS, NT, TS, TM, and TI) based on landmarks to represent physical partitioning. Once this atlas is created, it is possible to proceed to the main part of this Stage. For a given femoral scan, the algorithm from Stage 1 is utilized to generate a femoral segmentation. A rough mesh is then generated for the periosteal surface. Rough to fine rigid translations and rotations are utilized to bring the atlas mesh closer to the subject mesh through a grided k- nearest neighbors search. A mirroring and subsequent second rigid registration are then performed, as to avoid dependence on left vs. right annotation, which is not always provided though it is for this dataset. Following the final rigid registration, a non-rigid Iterative Closest Point (TCP)algorithm, originally sourced from the MATLAB community, is modified to benefit from GPU acceleration. Both meshes having the same number of vertices, this is deployed on the atlas mesh to nearly match the vertex coordinates of the subject mesh, which are relocated minimally. Particular atlas vertex indices correspond to landmarks which are then used to define the sub-femoral segmenting planes and lines. A nearest neighbor analysis yields the corresponding scan mesh indices for the landmarks, and as such, the landmarks are automatically identified. These landmarks are then used to define the lines and planes which separate the six sub-femoral regions. Upon visual inspection, this worked flawlessly on all scans, independent of mesh quality. The result is presented in Table 1.

[0368] With Stage 2 algorithms complete, Stage 3 algorithms can be developed. The purpose of this step is to find cubic samples appropriate for IH. These samples would be closest to the center of each of a given scan’s sub-regions. They should present sufficient trabecular content for FEA homogenization. Additionally, while samples should be selected randomly, it should be guided by a combination of sub-region centrality, sufficient trabecular volume fraction, and tolerable trabecular volume fraction. Trabecular volume fraction is considered here to be a compactly supported average of the trabecular mask over a 4 mm (49x49x49 voxel) window. Since convolution can be considered as a tensorial moving dotproduct, the trabecular volume fraction for every possible cubic specimen in a given volume can be simultaneously and efficiently computed via the GPU FFT accelerated step-and-save convolution algorithm described in Stage 1 (94):

[0369] (94)

[0370] where fv%,trab.gives the volume fraction of trabecular bone in a scan for the set of cubic specimens,andgive the discrete FFT and iFFT, respectively is thetrabecular mask within the given scan sub-region, and h is the convolutional averaging kernel. Crucially, the kernel is padded with zeros prior to the FFT to match the volume segmentdimensionality, denoted viap. Bold indicates the whole set in a scan, since multipleare computed simultaneously, where i indicates the z-thunit cell where there are NCeilsunit cells (95):

[0371]

[0372] Following computation of the trabecular volume fraction field, regional analysis is performed. First, the sub-femoral region is segmented. The 3D Euclidean distance field is then computed for the sub-region via a built-in fast approach in MATLAB, summarized here as (96):

[0373]

[0374] where dregis the regional distance field, DEuciideanis the Euclidean distance function, is the center point position of the ithunit cell, and cregis the position of the point at the center of the sub-femoral region under consideration. This distance field is multiplied with the logical masked region, and to maintain proximity to the center of the region, the innermost 4 cm region is considered. Each 4 mm specimen must have a trabecular volume fraction of at least 20% and at most the upper 95thpercentile in scan and sub-region- specific trabecular sample volume fraction within 4 cm of the region center (97):

[0375]

[0376] where gives the 95thpercentile for the given set. Thisequation indicates that specimens are desired which have typically sufficient material content, but not excessive bone content. We allow cortical specimens as TI is mostly comprised of them at the given minimum tolerable trabecular volume fraction.

[0377] The task of finding the specimen with optimal unit cell index, z*, while satisfying constraints on acceptable and centrality may be considered a minimization problem(98):

[0378]

[0379] where is the sub-femoral region domain and isafunction yielding auniformly distributed random value between a and b. In practice we impose the similarity constraint,withthe following inequality constraints (99):

[0380]

[0381] Because the distance transform is had, the optimal unit cell in each set can easily be found, solving the above constrained minimization problem as follows (100):

[0382] (100)

[0383] whereis the Hadamard product, best utilized when the operation is vectorized, done in practice. () is shown as a function of individual i to emphasize its’ dependent variables. The expression within is evaluated for all nonzero values (minimization doesnot accept cells not satisfying the constraints). The right-hand side ofis a binary rank-three tensor. As exemplified in Figs. 19D-19E, optimal unit cells are isolated across the domains. This completes Stage 3 of sampling. Next, we briefly describe the methodology for sample inverse homogenization, performance metrics, and comparisons to be made.

[0384] Demonstrated herein is the capacity of the CS-RBF cPLSM 3D IH algorithm in different kinds of materials and microstructural regions. The results are validated againstcommercial 3D homogenization in nTop. Hence, in addition to efficiency and speed, it is dependable for the task.

[0385] To assess the importance of design material elasticity, a broad spectrum of material is considered, some real and some hypothetical for the purpose of filling trends.Please see Table 2, below, for the list of material names and typical Elastic Moduli.

[0386] Table 2. Inverse design material elastic moduli.

[0387] The 14 materials range in Elastic Modulus from Carbon® SIL30 at 1 MPa to Ti- 6A1-4V at 110 GPa. Linear elasticity and a Poisson’s Ratio of 0.35 are assumed for all forward and inverse homogenization materials.

[0388] The IHP is defined in (101 ):

[0389] (101)-

[0390] To emphasize mechanical accuracy, we set P to 0.5% instead of a functionwith a minimum of 0.5% and a maximum of 7.5% depending on ij. Per the demonstrated mechanical accuracy, we employ the same GPU accelerated SOMMA algorithm. This is the same as the validation problem definition used herein.

[0391] Our only changes are that we reduce the emphasis on andconstraints relative to To do this, we simply change their c coefficients from 104to 102, and their d coefficients from 101to 10-1. The c and d coefficients for arekept at 103and 10°, respectively. Hence, depending on the initial relative residuals for the constraints, we place a larger emphasis on mechanical accuracy than on intermediate material reduction and swift termination, as well as the volume fraction constraints. Hence, theSOMMA coefficients are as follows (102):

[0392] (102)

[0393] The goal of this work is to assess the influence of anatomic region and design material elasticity on performance metrics for IH. Given its capacity to compare triclinic materials as demonstrated herein, for both regional and material assessment, we compare a custom formulation of symmetric Mean Absolute Percent Error (sMAPE), given by (103):

[0394] (103)

[0395] We also compare percent error (PE) in Kniib design volume fraction { theobjective), and the number of required design iterations. KHillcomputation is described by [Zhou YC, Zhao C, Wang F, Sun YJ, Zheng LY, Wang XH. Theoretical prediction and experimental investigation on the thermal and mechanical properties of bulk β-Yb2Si2O7. Journal of the American Ceramic Society 2013;96:3891-900. https: / / doi.Org / 10.l 11 l / jace.12618, incorporated by reference]. In region specifically, we compare skeletal sample volume fraction percentage as well as KHill•

[0396] Regarding statistical comparisons, region is categorical and material, represented by Elastic Modulus, is continuous. One-way comparisons are utilized on region to compare categories. If the groups are homoscedastic and normal, one-way ANOVA is used, and when heteroscedastic and / or non-normal, the Kruskal-Wallis test is used. If significant trends between groups are found at a significance level of 0.05, it is proceeded by Tukey HSDmultiple comparisons with Bonferroni correction. Tf the groups are heteroscedastic and / or non-normal, tied rank transformation is used first.

[0397] Regarding material comparisons, a different approach is required. As will be described herein, material has a nonlinear effect on the assessed metrics. Because it is continuous, nonlinear trend curves are developed. At each material, we approximate the median and generate 90% confidence intervals on the approximate. Regional effects are used to inform trend development.

[0398] Figs. 20A and 20B give regional comparisons on the sampled skeletal microstructures, and Fig. 20C-20F give regional comparisons on outcomes of IH. The sampled skeletal microstructures vary in bone volume over total volume, and Hill’s Bulk Modulus (KHill). This is indicated by several mixed significant trends in group medians via Figs. 2Fig. 0A and Fig. 20B, respectively. In contrast, regarding IH performance metrics, the only region to significantly influence sMAPE, and iterations is TI (p<0.05). In every test,rank based comparisons on group median were required. PE(KHill) was the only metric where no significant trends were observed in regional medians (p=0.1181).

[0399] Fig. 20 illustrates regional notched boxplots, where FIG. 20 A illustrates a bone sample volume fraction, extreme outliers removed in visual. FIG. 20B shows a bone sample Hill’s Bulk Modulus, KHill, extreme outliers removed in visual. FIG. 20C shows IH symmetric Mean Absolute Percent Error (sMAPE). FIG. 20D shows percent error in / cHill(PE( / cHill) ) of optimal design relative to skeletal target, extreme outliers removed in visual. FIG. 20E shows optimal design volume fraction percentage FIG. 20F shows design iterations. * indicatesp<0.05, ** indicates p<0.01, *** indicates p<0.001. Red marker indicates outliers. Rank based comparison of group medians used throughout.

[0400] Table 3, below, gives the total number of collected samples by region, and the percentage of which contained cortical bone. TI dominantly contains samples with corticalbone, and TI is the only region in which IH performance metrics significantly differ with respect to other regions (Figs. 20C-20F. Because of this, we exclude TI from IH performance analysis with respect to material, allowing us to merge the other regions as they do not present region-region differences with respect to one another following Bonferroni corrected multiple comparisons without TI.

[0401] Table 3. Percentage of regional samples containing cortical bone. Region with highest cortical presence highlighted in the last row.

[0402] Figs. 21A-21D presents the IH performance metrics pooled across all subfemoral regions. Median curves are given, as well as 5,000 cycle bootstrap based 90% Cis on the metric median at each material. As with regional analysis, no outlier removal is performed. The critical / viable region is given as that which corresponds to an sMAPE upper CI of at least 20%. This is found to exist between 0.63 GPa and 30.70 GPa (Fig. A) and is represented by a transparent red patch in Figs. 21B and 21D.

[0403] Fig. 21 A shows symmetric Mean Absolute Percent Error across the constitutive matrices CTargand C* (sMAPE, %) with respect to Elastic Modulus (E, GPa) in 14 materials (see Table 2, above. 90% 5,000 cycle bootstrap based Confidence Intervals (Cis) on median at each material, pooled across anatomic regions (material pooled and regionally stratified inFigs. 2Fig. 1C-21D. The 20% maximum sMAPE cutoff is measured, the solution spanning materials wherein maximum 90% CI is helow or at the threshold (0.63 GPa < E < 30.70 GPa). Few notable materials are indicated with arrows pointing to respective point estimates (mean, spline interpolated blue dotted dashed curve). 90% CI domains are given as spline interpolated transparent blue fills with the bounding approximated 90% limits as dashed blue curves. For all materials, the 90% Cis are given in orange with a dot at the computed mean and the caps terminating at the computed 90% CI ends. FIG. 21B shows as aThe transparent red shaded region represents the area covered by theelasticity bounds in which median sMAPE is at most 20% with 90% confidence. FIG. 21C shows for the final designs. Note saturation at or near 30.70 GPa. FIG. 21D shows the number of design iterations prior to termination. Note the saturation, a result of intermediate material removal at iteration 60 and subsequent termination.

[0404] Referring to FIGS. 22, 23, 24 and 25, which show that for a FH sample that sMAPE, PE(KHill), and iterations follow a saturating U-shaped trend with respect to material elasticity. Low elastic moduli are associated with more dense, high percent error, and tumultuous iteration count (Fig. 22). The given SIL30 example terminates with an sMAPE of 292.55%, of of 75.13%, and at 7 iterations. This is in line with thepresented trends, as it shows how material fills in the space, and because constitutive targets are so difficult to satisfy, the SOMMA algorithm quickly satisfies the PTI constraint.

[0405] Figs. 22A-22G show SIL30 IH design history, where FIG. 22A shows a final design; FIG. 22B shows a GCM component percent errors; FIGS. 22C-22F show projections of the target and final axial elasticity surfaces; and FIG. 22E shows iteration history.

[0406] Mid-range elastic moduli are associated with low density, low percent error, and medium iteration count (Fig. 23). The given 18 GPa example terminates with an sMAPE of5.10%,of 3.06%, , and at 29 iterations. This is in line with thepresented trends, as it shows how such a material has a fvo / onear the saturation limit. Asdescribed herein, is increased while P are reduced. Once aresatisfactorily reduced, is then reduced, here slowly due to the disclosed changes inthe SOMMA parameters. Not part of the material testing, we inferred such performance based on elasticity. Validation via nTop is provided in Fig. 24. Note the similarity in the GCMs.

[0407] Fig. 23 illustrates 18 GPa Elastic Modulus sample design history, where FIG. 23A shows a final design; FIG. 23B shows GCM component percent errors; FIGS. 23C-23F show projections of the target and final axial elasticity surfaces; and FIG. 23E shows iteration history.

[0408] FIG. 24 shows target (left) and optimal (right) units cells at 18 GPa Elastic Modulus. High elastic moduli are associated with less dense, high percent error, and high iteration count (Fig. 25). The given Ti6A14V example terminates with an sMAPE of 253.27%, PE( / cHill) of 92.86%, of 23.52%, and at 62 iterations. This is in line with the presentedtrends, as it shows how the design remains in an intermediate state between voluminous posts. Due to the difficulty associated with forming a lean, rigid connection, iteration proceeds for extended iterations until a hard cut in intermediate material occurs for termination at 60 iterations.

[0409] Fig. 25 shows Ti6A14V IH design history, where FIG. 25 A shows a final design; FIG. 25B shows GCM component percent errors; FIGS. 25C-25G show projections of the target and final axial elasticity surfaces; and FIG. 25E shows iteration history.

[0410] In summary, disclosed is a test of the disclosed CS-RBF cPLSM 3D IH algorithm on triclinic targets from microstructurally diverse anatomic regions, demonstrate the importance of design material on IH outcomes, and yield an efficient and rapid automated rule-abiding microstructural sampling approach. We achieved these aims by applyingconvolution in unique ways, here in efficiently computing volume fraction and in performing morphologic operations.Design And Mechanical Testing Of Bone With Targeted Triclinic Spatially Variate Elasticity, And Adjacent Microstructure Synergy

[0411] In yet another aspect, design and mechanical testing of bone with targeted triclinic spatially variate elasticity, and adjacent microstructure synergy is disclosed, which adapts the convolutional approach to sub-Gaussian Point 3D inverse homogenization algorithms described herein. One disclosed implementation comprises performing the testing utilizing 3D with multiple triclinic elasticity targets, yielding up to 288 mechanical continuum constraint. This aspect further discloses and implements strategic concatenation of design variables and sensitivity, and de-concatenation for isolated design update.

[0412] Biomedical engineering innovation and epidemiological health disparities have jointly advanced in recent decades. Between 2006 and 2014, total hip replacement (THA) and revision (rTHA) increased in prevalence by 69.5% and 28.5%, respectively. In contrast, the US population had increased by approximately 6.68% within that timeframe. The rate of change in THA prevalence greatly exceeds that of population growth.

[0413] Using THA as a correlative to fracture risk, an enhanced understanding of skeletal fracture risk prediction and mitigation is essential. 1-year post operative mortality following hip fracture is about 27.3%. To combat this issue, research efforts in skeletal biomechanics and aging have accelerated tremendously.

[0414] Elasticity mismatch-based stress shielding is a primary cause for revision surgery. To improve microstructural synergy, recent work has focused on the “bottom-up” design of implants. For example, in bone, this may mean that individual scan data / models are discretized spatially. The mechanics (i.e., bulk modulus) of each functional Unit Cell, or UC, are then computed. Subsequently, the microstructures may be designed in a way which bringseffective properties closer to that of the native spatially located bone, while preserving connectivity with designed neighboring UCs. In computational biomechanics, there is ample support for this multiscale modeling approach.

[0415] Since the 1990s’, effective mechanics of the detailed skeletal microstructure have been investigated via applied continuum mechanics. Today, we use similar classic approaches, scaled with modem hardware such as the graphical processing unit (GPU). In skeletal biomechanics, the modelled macroscale domain is quite finite, which theoretically violates continuum homogenization assumptions necessary in representing a detailed microstructure with a single effective material. Yet, despite this, the skeletal biomechanics community has made extensive use of computational homogenization to simplify otherwise intractably large, detailed, models for simulation. Other than bone, this approach has also been utilized in different materials such as calcifying cartilage. Crucially, the approach of representing several million element skeletal structures with a few thousand elements, and still obtaining accurate mechanics, has been validated.

[0416] Hence, under the assumption that homogenized UC mechanics can be mimicked via Inverse Homogenization (IH) continuum optimization, a macroscale structure may be designed bottom-up to effectively mimic macros tructural mechanics. This; however, has been a topic of much research in recent years.

[0417] There are both old and new challenges in modern bottom-up inverse design. Old challenges include those present when considering fabrication of continuum structures which theoretically repeat indefinitely. It has been well understood that Saint-Venant’s principle is non- negligible in homogenized continuum structures. A thorough analysis on the influence of the ratio of Representative Volume Element (RVE) size to UC size (p = LRVE / LUC) relative to effective macroscale properties is presented with a focus on stress concentrations.

[0418] New challenges include adjacent microstructure synergy. This compound unit cell IH approach (herein, cIH) attempts to support appropriate communication of load paths across UCs for accurate gross mechanics. This approach resolves UC connectivity issues, which do not naturally occur in many top-down approaches. Though in some instances designs featuring= 0.5, as two UCs yield over 80% efficacy in full scale analysis relative to homogenized material analysis,= 0.25 yields 90% efficacy. However, for practicality in micromechanical mimicry and fabrication on similar length scales, large, detailed designs cannot easily achieve = 0.25 per computational memory and printer resolution, at least inour experience with skeletal micro-mimicry. This lends great insight into the tradeoff between representational accuracy and practicality, influencing the course of the ensuing work.

[0419] The disclosed development of the Compactly Supported Radial Basis Function convolutional Parametric Level Set Method (CS-RBF cPLSM) has enabled the design of 3D microstructures with commercial solver validated mimicry of linear elastic FE-homogenized anisotropy on a triclinic skeletal target and a cubic Gyroid target. We then applied the algorithm towards the evaluation of IH performance metrics on a spectrum of materials, sorted by Elastic Modulus, as well as several sub- femoral targets in 24 cadaveric specimens scanned at 82 μm isotropic resolution via High Resolution peripheral Quantitative Computed Tomography (HR-pQCT).

[0420] In contrast with previously presented bottom-up methods, via derivation, we observe inherent convolutional smoothing within CS-RBF cPLSM collocation and sensitivity analysis. We utilized the Wendland function for enhanced efficiency. This approach was used to solve a maximum target normalized Percent Error (PE) variation on the classic IHP. We incorporated derived intermediate material control for automated sensitivity adjustments and rapid termination. Furthermore, we implemented the Second Order Method of MovingAsymptotes (SOMMA), using the Hessian diagonal, to emphasize constitutive targeting. Thishas not been done in IH literature otherwise. This approach, unlike other IH approaches, provides extensive computational efficiency while yielding 3D sub-Gauss point accurate LS designs.

[0421] Because of the derived presence of convolution, expansion coefficient fields, a, yield spatial LSFs, which are smooth and continuous and may be repeated in a triplyperiodic fashion. Further, this yields triply periodic binary continua upon application of the approximated Heaviside function, Per derivations presented herein, thisperiodic smoothing applies to sensitivity analysis as well. We hypothesize that these properties make the CS-RBF cPLSM an ideal candidate for a sequential compound cell TH procedure. Our rationale is that individual cells can be designed one at a time, and that the effects of neighboring cells may be accounted for concurrently in each design iteration such thaadapts at neighbor boundaries. Depending on compact support radius, R, and the spatial distribution of design knots, = 0.5 may be sufficient.

[0422] We disclose that if spatially variate homogenized mechanics can be synergistically replicated, then macrostructural mechanics can be generally mimicked, which is then tested. There are several gaps in conventional technology we resolve in the disclosed approach: (1) 3D compound IH with an optional second order approximate optimization approach; (2) bottom-up multiscale structural optimization which inherently adapts the current design to neighboring designs; (3) preserved targeting of current native bone properties in the current design while also targeting the compound properties of the “rectangular box”, resulting from concatenation of the current native bone with the neighboring native bone and similarly for the spatially located designs, without altering the finalized design variables for the neighboring design; implementing all on a custom Windows 11 PC on an NVIDIA RTX 3090 GPU in MATLAB and partially in CUDA C; (4) A 3D cIH benchmark for algorithm comparison on triclinic targets ;(5) mechanical testing of the result;(6) Benchmark testing on Ti-6A1-4V; (7) A neighborhood planning algorithm for UC compounds and sequential progression; and Adaptive MMA parameters based on sMAPE.

[0423] Presented below is a sequential compound UC planning algorithm, approximate optimization algorithm settings, performance metrics to collect, benchmark development, anatomic data acquisition, full surrogate design and fabrication, statistical comparisons, and finally mechanical testing for validation. Subsequently, collected performance metrics are presented herein, compared across compound count, and compare to findings from prior literature.

[0424] All computations were performed on a custom Windows 11 machine. Relevant hardware includes an NVIDIA 3090 GPU with 24 GB VRAM, an AMD Ryzen 9 5950X CPU, and 128 GB RAM. It is to be appreciated; however, that other machines and components may be used in implementations of the disclosure.

[0425] The disclosed compound Inverse Homogenization Problem (cIHP) is defined in (104):

[0426]

[0427] where r represents the current UC number, f represents the adjacent previous UC number which is in the θ axial direction, represents the concatenation of theT and ( UC design variable sets in the 9 axial direction,is the critical maximum absolute target value normalized percent error for theGCM component in the θ-direclionUC compound, is the respective normalized percent error,is the volume fraction percentage for the r UC design, and the rest of thecomponents such as and Stop( are acting on the current design,

[0428] Please see Fig. 26, which displays the central skeletal microstructure in a 3x3x3 macrostructural discretization, r and ( follow a planning algorithm. In the following, the planning algorithm is disclosed for the sequential cIH process. Note that the 3x3x3 example later becomes our benchmark. Fig. 26 shows the central skeletal microstructure in a 3x3x3 discretization and compounds relevant to the disclosed sequential cIH procedure.

[0429] A bottom-up approach is developed herein, not requiring cells solely interpolating designs. Provided is a reduced approach to programmatic multi compound-cell constrained optimization: (1) Label the set of UCs with their index sorted axially by X, Y, and Z; (2) Establish a face-to-face prior UC adjacent connectivity list; (2.1) i.e., first cell has nothing behind it, to the left, or above it.; (2.2) Current UC index is denoted as r, the prior adjacent UC index is denoted asis the axis index (0 for no connection, 1 for -X adjacency, 2 for -Y adjacency, and 3 for -Z adjacency) of adjacency, used respectively for bivariate concatenation / compounding function and univariate de-concatenation function(3)For the first UC: (3.1) Solve IHP; (3.2) Store: (3.2.1) The current optimal expansion coefficient set to the system storage, as many are anticipated in large structures; (3.2.2)The current optimal and targeted homogenized elasticity matrixand ; (4)Progress to the next UC in the order of +X, +Y, +Z (assuming that there may be no more UCs across θ): (4.1) If N or θ #= 0, where N is the number of current -X, -Y, and -Z connections, perform a compound cell inverse homogenization (cIH) as follows: (4.1.1) From systemstorage, obtain the adjacent optimal cell design (ctf ), in correspondence with the connection list; (4.1.2) Holdfixed while modifying and as per N+l (adjacent + current)homogenization / sensitivity analyses: In 3D, the number of constraints on the multi-compound IHP becomes 36(N + 1) + 3. This accounts for the N cIH and 1 immediate IH constraints as well as 2 volume fraction constraints and the PTI constraint; (4.1.3) Pass the new set of objective / constraint evaluations and sensitivities to the MMA / SOMMA algorithm to iterate aT; (4.1.4) Repeat previous step until termination; (4.2) Otherwise, perform regular IH (5) Repeat the previous step until all Unit Cells have been optimized to best satisfy their compound representation(s).

[0430] Alternately, see , which shows a7, which shows a graphic procedure for the CS- RBF cPLSM cIH, which shows a workflow for the sequential CS-RBF cPLSM cIH procedure. A volume is discretized. UCs are labelled in order of +X, +Y, and +Z. The connectivity table is computed. On the first row, 0 = 0 and N = 0. So, and there areneither concatenations nor de-concatenations. is instantiated and the target,jscomputed. Perform cIH (here, only a standard IH). Store needed information. Check if the following row exists. If it does not, exit. If it does, progress to the next table row. Then, cycle through the first four columns, storing initial andwhich is computed percompound by column, for 0 0 are created by compounding the table based prior adjacent UCs with the current design variables, . cIH is performed by appending allconstraint evaluations and sensitivities for all designs T and f, and updating onlyandaThas an influence on the total design via compact support and load paths. The optimal immediate UC expansion coefficient field, is stored, along with the resulting homogenizedelasticity matrices for the designs and their targets,andCheck that the next row, theUC, exists. As for the 1stUC, if it does not, exit. If it does, repeat the loop untilit does not. In this way, cIH is sequentially performed to generate the cell designs required for macrostructural assembly.

[0431] A benchmark is disclosed which exemplifies anisotropy mimicry in 3D. Bone, typically approximated as orthotropic with historically acceptable error, is more accurately triclinic / non-orthotropic. Specifically, the assumption of orthotropy holds to within 6.1% in normalized residual differences for trabecular microstructure cubes of 4 mm side length. Hence, such sized microstructures are typical. Note that in the disclosed compounding approach, we exceed 4 mm side length, with certain dimensions having 8 mm side length.

[0432] The benchmark is developed using publicly available elderly cadaveric femur HR-pQCT data at 82 μm resolution. A cube of side length 12 mm is extracted from a scan. Subsequently, the cube is discretized into 4 mm sub-cubes, shown in Fig. 28, which shows the discretized femoral benchmark, with bone and discretization bounds shown in different shading.

[0433] We apply the developed approach to this assembly of microstructures. The design material is Ti-6A1-4V, a popular material in implant design. This material has an Elastic Modulus of 110 GPa. We employ SOMMA using the following settings (105):

[0434] (105)

[0435] For m c and d adapt at each iteration if symmetric Mean AbsolutePercent Error (sMAPE, discussed herein) is less than or equal to 20%, as follows (106):

[0436] (106)

[0437] Further disclosed is an inverse design of a vertebral body targeting accurate Spatially Variate Triclinic Elasticity (SVTE). Per the computing time associated with such a task, we limit the disclosed approach to the design of a single vertebral body. We were provided an LI specimen, frozen past seven years and thawed at the time of scanning and testing. The donor was an 88-year-old male, passed from chronic kidney disease. Using a second-generation HR-pQCT scanner (XtremeCT II, Scanco Medical AG, Briitisellen, Switzerland), we scanned the specimens at an isometric voxel size of 60 μm.

[0438] Due to the porosity of such specimens, and the materials available to us, we proceed to design the vertebral body in IND407, Clear Loctite. This material has an Elastic Modulus between 1.385 GPa and 1.5 GPa and allows us to observe the internal architecture. For speed in design, we employ the MMA algorithm rather than SOMMA using the same first order approximate optimization settings . Fabrication took place on a Carbon DLS™ (Digital Light Synthesis) Ml 3D printer (Ml, Carbon Inc., California, USA). The printer has an inplane resolution of 75 μm and an out of plane resolution between 25 μm and 100 μm depending on the desired level of detail and print duration.

[0439] Five ideally identical surrogate LI replicates were printed. We developed custom loading platforms conformal to the vertebral body end plates. These would ensure that samples and specimens would have the same orientation, and that the microstructure was not irregularly constrained. Additionally, platforms prevent sample / specimen translation. We employ quasi-static compression at 1.5 mm / s until a 10 mm displacement or 75% reduction in load magnitude relative to peak. This is done on an MTS uniaxial mechanical testing system.

[0440] In both the benchmark and full vertebral design, we assess the following metrics.Given capacity to compare triclinic materials, we compare hMAPE, iteration count, duration,volume fraction percentage PTI, and a cIH adaptation on the disclosed formulation ofsymmetric Mean Absolute Percent Error (sMAPE), given by (107):

[0441] (107)

[0442] sMAPE and hMAPE are computed using all compounds. Regarding statistical comparisons, the number of compounds, N, that a UC has is ordinal / categorical. One-way comparisons are utilized herein. The Kruskal -Wallis test is used per low sample size and heteroscedasticity. If significant trends between groups are found at a significance level of0.05, it is proceeded by Tukey HSD multiple comparisons with Bonferroni correction following tied rank transformation. Statistical tests are not performed on the benchmark due to the single N = 0 case.

[0443] Fig. 29 displays the central UC design and the relevant compounds. Note the conformal designs achieved using consistent r design variables across compounds. While not necessarily unique to the disclosed approach with respect to standard structural optimization problems, inherent convolution bridges designs with boundary smoothness not seen in prior cIH literature. Note how elasticity fields are not all perfectly matched. This is because Ti-6A1- 4V is a challenging material for femoral inverse design. Despite this, the disclosed sMAPE adaptive SOMMA approach terminates with a fine binary design, especially given the challenge of satisfying triclinic constraints of four UCs simultaneously. Fig. 29 shows the central UC design and relevant compounds (columns). (Top) skeletal microstructures. (2ndrow) inverse designs. Intermediate material is shown in blue, the methodologically relevant adjacent designs are in gray. The current design is in red. (3rdrow) axial elasticity map. (Bottom) log-y scale function evaluation curves. The design variables driving the primary design, and all compounds are updated simultaneously.

[0444] Figs. 30A and 30B show the overall benchmark design (30B), split by cell, and the skeletal sample from Fig. 28 (30A) (duplicated for visual comparison), r) = 0.5 is used. Note similar spatially relative porosity and directionality, particularly in the mid-section. Figs. 30A and 30B show visualized via nTop, discretized benchmark sample (30A) and Ti- 6A1-4V design (30B) via disclosed sequential CS-RBF cPLSM cIH procedure. Two-UC replication is used (77 = 0.5).

[0445] Figs. 31A-31F show the metric distributions for the benchmark across N. Statistical analysis is not performed as N = 0 occurs once. Note that the evolving SOMMA approach yields consistent iteration count, hM APE, sMAPE, fv%, and PTI. Duration increases seemingly linearly with N. In cIH, a total of 3.55 hours is spent in cIH, not considering forward homogenization evaluation of the skeletal microstructure(s) in the beginning for target calculation. Fig. 31 shows cIH performance metrics on benchmark across N. Outliers not shown.

[0446] Fig. 32 gives the distribution of N across the cadaveric LI discretization. Fig. 32 shows primary UC count vs. N.

[0447] Figs. 33A-33F show the metric distributions for the full vertebral body across N, as well as statistics. Relative to the benchmark, poor mechanical accuracy is observed as anticipated. This is because second order information and sMAPE based adaptation are foregone. Additionally, c and d that emphasize PTI reduction over mechanical accuracy are utilized for speed. So, short iteration count is observed. Additionally, Clear Loctite is used. However, considering evaluation on a porous vertebral body, it is still a challenge despite use of a material with relatively lower elasticity. The proposed algorithm executed cIH over 4.07 days. The Kruskal- Wallis test implies, with p < 0.05, that there are significant trends in all metrics with respect to N. Via the post-hoc procedure, it is evident that increasing N yields increasing duration. In the other metrics, there exists a multitude of significant trends. Themost notable include the increase in sMAPE for N 0 (p < 0.05), the increase in fvo / ofor N = 2, 3 (p < 0.05), and the apparent reduction in PTI for N > 0. These strongly imply a focus on PTI reduction rather than mechanics. Fig. 33 shows notched boxplots for cIH performance metrics on the vertebral body across N. * indicates p < 0.05, ** indicates p < 0.01, *** indicates p < 0.001. Red / shaded markers indicate outliers. Rank based comparison of group medians used throughout.

[0448] Fig. 34 shows the scan, inverse design, and fabrication per the vertebral body with t] = 0.5. Note the apparent preservation of the cortex. The algorithm effectively produced rigid designs in correspondence with increased material Elastic Modulus. The cortex is given in purple while the trabecular region is given in yellow. The same automated convolutional segmentation procedure we developed for the femur is utilized herein. A coarse variant of the model presented middle in Fig. 34 is fabricated, as shown in Fig. 34B, right. Fig. 34A shows a cross-section of trabecular microstructure (yellow) and cortex (purple) of scanned vertebral body. FIG.34B shows a cross-section of designed model. FIG. 34C shows a supported model fabricated in Clear Loctite, on a build-plate.

[0449] Fig. 35 gives (35A) the processed force-displacement curves following uniaxial compression for the five surrogates and the cadaveric specimen, (35B) a boxplot of median surrogate stiffness and a line for specimen stiffness, and (35C) the percent error in surrogate stiffness relative to specimen stiffness, with a range in expected percent error given by mean hMAPE and sMAPE weighted by N (right). Per a one-sample two-tailed t-test, specimen stiffness is significantly lower than mean surrogate stiffness by 0.655 kN / mm (p « 0.05). This is anticipated as poor mechanical accuracy was yielded by the MMA approach used.

[0450] hMAPE gives a lower bound on anticipated percent error, and sMAPE gives an estimate of percent error based on large magnitude constitutive components and their relatively large magnitude residuals. A one-sample right-tailed t-test indicates that meansurrogate stiffness percent error is greater than N-weighted mean hMAPE of 24.44% by 13.45% (p < 0.05). Per Fig. (c), it is close to the minima of surrogate stiffness percent error. By a one-sample two-tailed t-test indicates that mean surrogate stiffness percent error is not significantly different than N-weighted mean sMAPE of 50.04% (p = 0.095). 95% confidence intervals on surrogate stiffness percent error range from 22.43% to 53.35%. These results indicate a potential to estimate percent error in stiffness based on variations of percent error in computational modeling.

[0451] While the computational metrics do not explain variability in stiffness across surrogates, they do provide bounds on expected percent error relative to the target, presuming accurate biomechanical modeling following 7-year frozen storage. This itself is expected to negatively impact indentation testing results, and presumably Elastic Modulus.

[0452] Fig. 35A shows post-processed force -displacement curves following uniaxial compression for the five surrogates and the cadaveric specimen. The approximate linear stiffness region is given in solid green. The intercept with the x-axis is used to zero displacement, and the extending line is given in dashed green per curve. FIG. 35B is a boxplot of median surrogate stiffness and a line for specimen stiffness. FIG. 35C shows percent error in surrogate stiffness relative to specimen stiffness. The range in expected percent error is given by mean hMAPE and sMAPE weighted by N.

[0453] The results were abstracted to 3D using the convolutional level set method.Computing Environment

[0454] Referring to Fig. 36, an example computing device 3000 upon which embodiments of the invention may be implemented is illustrated. In some instances, the computing device 3000 of Fig. 36 comprises a Graphical Processing Unit (GPU), as referred to herein. It should be understood that the example computing device 3000 is only one example of a suitable computing environment upon which embodiments of the invention may beimplemented. Optionally, the computing device 3000 can be a well-known computing system including, but not limited to, personal computers, servers, handheld or laptop devices, multiprocessor systems, microprocessor-based systems, network personal computers (PCs), minicomputers, mainframe computers, embedded systems, and / or distributed computing environments including a plurality of any of the above systems or devices. Distributed computing environments enable remote computing devices, which are connected to a communication network or other data transmission medium, to perform various tasks. In the distributed computing environment, the program modules, applications, and other data may be stored on local and / or remote computer storage media. The disclosed computing device 3000 may comprise all or a part of the controller used to implement the above-described methods of inverter control.

[0455] In its most basic configuration, computing device 3000 typically includes at least one processing unit 3006 and system memory 3004. Depending on the exact configuration and type of computing device, system memory 3004 may be volatile (such as random-access memory (RAM)), non-volatile (such as read-only memory (ROM), flash memory, etc.), or some combination of the two. This most basic configuration is illustrated in FIG. 36 by dashed line 3002. The processing unit 3006 may be a standard programmable processor that performs arithmetic and logic operations necessary for operation of the computing device 3000. The computing device 3000 may also include a bus or other communication mechanism for communicating information among various components of the computing device 3000.

[0456] Computing device 3000 may have additional features / functionality. For example, computing device 3000 may include additional storage such as removable storage 3008 and non-removable storage 3010 including, but not limited to, magnetic or optical disks or tapes. Computing device 3000 may also contain network connection(s) 3016 that allow the device to communicate with other devices. Computing device 3000 may also have inputdevice(s) 3014 such as a keyboard, mouse, touch screen, etc. Output device(s) 3012 such as a display, speakers, printer, etc. may also be included. The additional devices may be connected to the bus in order to facilitate communication of data among the components of the computing device 3000. All these devices are well known in the art and need not be discussed at length here.

[0457] The processing unit 3006 may be configured to execute program code encoded in tangible, computer-readable media. Tangible, computer-readable media refers to any media that is capable of providing data that causes the computing device 3000 (i.e., a machine) to operate in a particular fashion. Various computer-readable media may be utilized to provide instructions to the processing unit 3006 for execution. Example tangible, computer-readable media may include, but is not limited to, volatile media, non-volatile media, removable media and non-removable media implemented in any method or technology for storage of information such as computer readable instructions, data structures, program modules or other data. System memory 3004, removable storage 3008, and non-removable storage 3010 are all examples of tangible, computer storage media. Example tangible, computer-readable recording media include, but are not limited to, an integrated circuit (e.g., field-programmable gate array or application-specific IC), a hard disk, an optical disk, a magneto-optical disk, a floppy disk, a magnetic tape, a holographic storage medium, a solid-state device, RAM, ROM, electrically erasable program read-only memory (EEPROM), flash memory or other memory technology, CD-ROM, digital versatile disks (DVD) or other optical storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices.

[0458] In an example implementation, the processing unit 3006 may execute program code stored in the system memory 3004. For example, the bus may carry data to the system memory 3004, from which the processing unit 3006 receives and executes instructions. Thedata received by the system memory 3004 may optionally be stored on the removable storage 3008 or the non-removable storage 3010 before or after execution by the processing unit 3006.

[0459] It should be understood that the various techniques described herein may be implemented in connection with hardware or software or, where appropriate, with a combination thereof. Thus, the methods and apparatuses of the presently disclosed subject matter, or certain aspects or portions thereof, may take the form of program code (i.e., instructions) embodied in tangible media, such as floppy diskettes, CD-ROMs, hard drives, or any other machine-readable storage medium wherein, when the program code is loaded into and executed by a machine, such as a computing device, the machine becomes an apparatus for practicing the presently disclosed subject matter. In the case of program code execution on programmable computers, the computing device generally includes a processor, a storage medium readable by the processor (including volatile and non-volatile memory and / or storage elements), at least one input device, and at least one output device. One or more programs may implement or utilize the processes described in connection with the presently disclosed subject matter, e.g., through the use of an application programming interface (API), reusable controls, or the like. Such programs may be implemented in a high-level procedural or object- oriented programming language to communicate with a computer system. However, the program(s) can be implemented in assembly or machine language, if desired. In any case, the language may be a compiled or interpreted language and it may be combined with hardware implementations .

[0460] Commercial applications for the disclosed systems, methods, devices and computer program products include: a. skeletal surrogates for surgical training and practice, skeletal surrogates for biomechanics testing and medical device performance validation (this could include groupings of patient populations by age, gender, ethnicity, disease state, etc.), physiologically accurate injurydummies for use in dynamic injury (i.e., automotive, military, etc.). The technology and production services would be scalable. Customers could include clinicians, biomechanics researchers, medical device manufacturers, and injury testing labs. b. Biomechanically equivalent implants, where risk of stress shielding based fractures would be minimized due to analogous mechanics.

[0461] Alternate uses of the disclosed systems, methods, devices and computer program products include: a. the mimicry of porous media phenomena such as anisotropic permeability, permisivity, heat transport, etc. This would make it useful in a whole host of engineering domains. b. The use of this continuous optimization approach in the finance domain. c. General structural optimization tasks, such as jet wing optimization. d. It can be abstracted for use with any type of material and can be used for dynamic problems. e. Repeatable biomechanical testing, medical device performance evaluation, and surgical training or practice.

[0462] Although the subject matter has been described in language specific to structural features and / or methodological acts, it is to be understood that the subject matter defined in the appended claims is not necessarily limited to the specific features or acts described above. Rather, the specific features and acts described above are disclosed as example forms of implementing the claims.

[0463] Inventive concepts disclosed and described herein include systems, methods, devices and computer program products for:[1] a novel and highly efficient DoF by DoF algorithm for solution of indeterminate mechanical problems wherein each element has a unique local stiffness matrix, critical in this endeavor as each element has varying material distributions continuously.[2] a GPU accelerated nD FFT convolution algorithm for MATLAB, critical to efficiency and low memory usage.[3] a novel GPU accelerated nD FFT convolution approach to finite element integration point stress and strain evaluation.[4] a novel GPU accelerated approach to compactly supported-radial basis function (CS- RBF) evaluation for the parametric level set method (CS-RBF PLSM) for sub-element accuracy, ELIMINATING the need for a dense linear systems solution and ENABLING expansion to very large 3D problems, where previously only small 2D problems were feasible.[5] a novel convolutional approach to sensitivity analysis and visualization:[5a] a novel approach to computation of the Jacobian of the CS-RBF PLSM design variables.[5b] a novel approach to visualizing the Jacobian of the CS-RBF PLSM design variables; this innovation enables the user to visualize what parts of a structure a responsive to, and responsible for resistance to, loading.[5c] a novel approach to computation of the Hessian of the CS-RBF PLSM design variables (never before done, due to what was previously extremely computationally intensive).[5ca] a rapid FFT approach to 3-way convolution. This extremely rare process is where two convolutional kernels are had for a given signal, and for each pass of the first kernel, the second kernel is completely convolved over the domain, which ultimately yields an output where all possibly effects of design variables with respect to each other are computed. This process has yet to be elucidated in literature to the best of the authors knowledge. An extremely efficient FFT convolution approach was derived,and an implementation executed in MATLAB. However, for the same reason as to why the Hessian has never been needed in CS-RBF PLSM structural optimization, there are approximation approaches with their own efficient approximations of the Hessian based on the Jacobian and other information. As such, in the current implementation of the overall invention described here, this algorithm is optional. Still, it is a novelty and likely has a large impact to be had either in the future or in other problems / situations.[5d] a novel approach to visualization of the Hessian of the CS-RBF PLSM design variables (never before computed, and thus, never before visualized).[6] a GPU accelerated variant of the famous Method of Moving Asymptotes (MMA) algorithm, and a GPU accelerated variant of the Globally Convergent MMA (GC-MMA) algorithm, both featuring updates yielding improved speed and computational efficiency.[7] a novel implementation in solving the inverse homogenization problem in 3D. This problem is where you are given only the effective properties of a metamaterial / metastructure / microstructure, and you need to generate a structure with the same or approximate properties in a given material. This problem has never been solved simultaneously for all 21 unique components of the generalized elasticity matrix in 3D. It has not been possible with the SIMP method due to numerical instabilities such as checkerboarding and the occurrence of significantly prevalent intermediate densities. In this approach, material is defined continuously, and there are tight limits on intermediate densities to the point that they are negligible.[8] A direct to STL algorithm for structures defined by the CS-RBF design variable solutions, with sub-Gauss Point accuracy.

Claims

CLAIMSWhat is claimed is:

1. A method of bottom-up inverse design for micromechanical spatially variate anisotropy mimicry, comprising: obtaining an image of a microstructure; spatially mapping linear elastic anisotropy of the microstructure; and creating a printable macrostructure with synergistic microstructures.

2. The method of claim 1 , wherein obtaining an image of the microstructure comprises obtaining an image of a microstructure of bone.

3. The method of claim 2, wherein the bone is cadaveric, and the image is obtained using Micro-computed tomography (micro-CT).

4. The method of claim 2, wherein the bone is of a living being, and the image is obtained using small region high-resolution peripheral computed tomography (HR-pQCT).

5. The method of claim 2, wherein the bone is of a living being, and the image is obtained using large region 3D Al Super Resolution.

6. The method of claim 5, wherein large region 3D Al super resolution comprises: gathering low resolution clinical CT imaging data; employing a pretrained 3D super-resolution generative adversarial network (SRGAN) framework for skeletal microstructure super resolution, comprising:a wide activation for efficient and accurate image super-resolution-B (WDSR-B) generator; a standard convolutional neural network (CNN) discriminator; and a voxel shuffle algorithm for sub-voxel up-sampling; and gather generated super resolution inferred data into macro- volume and continue.

7. The method of any one of the above claims, wherein spatially mapping the linear elastic anisotropy comprises NbN-PCG FEA Homogenization.

8. The method of claim 7, wherein NbN-PCG FEA Homogenization comprises:Step 1 - compute a list of local finite element matrices, as in:Step la - evaluating cp via convolution;Step lb. Perform numerical integration of local stiffness matrices and load vectors ( unitary residual strain) via Gaussian Quadrature, wherein several Gauss points should be considered, so that heterogeneity due to phi can be captured;Step 2 - From the list of local stiffness matrices, assemble a list of stencil coefficients, wherein these stencils are stored in a column-major fashion, for coalesced access in massively parallel matrix-vector products;Step 3 - From the list of local force vectors, assemble global force vectors;Step 4 - Copy the list of stencils and global force vectors from CPU to GPU ( host to device);Step 5 - Solve the systems of equations for each load vector using Krylov subspace solvers, wherein the Jacobi-preconditioned PCG method is a straightforward and efficient option of solver, which depends on:Step 5a - Implementing sparse matrix-vector products in GPU as convolution- like operations that use the list of stencils;Step 5b - Implementing daxpy and reduce operations in GPU;Step 6 - Copy the resulting displacement vectors back, from GPU to CPU ( device to host);Step 7 - Compute resulting stress fields from the displacement vectors, and integrate these stress fields to obtain effective stress at the macroscale; andStep 8. Compute the effective elasticity tensor from the obtained effective stress fields.

9. The method of claim 8, wherein evaluating cp via convolution comprises: via convolution, analytically evaluating cp; solving for loading and stiffness matrices in each element by direct Gaussian-Integration on approximate Heaviside of cp; store coalesced arrays of DoF stencil coefficients, taking advantage of gridded regularity: and solve the DoF displacements using a matrix-free implementation of Krylov subspace solvers in GPU.

10. The method of any one of the above claims, wherein creating the printable macrostructure with synergistic microstructures comprises using a sequential compound inverse homogenization in 3D with convolutional level set method procedure.

11. The method of claim 10, wherein the sequential compound inverse homogenization in 3D with convolutional level set method procedure comprises:(1) discretizing a volume of the microstructure volume into a set of unit cells (UCs);(2) labeling the set of UCs in +X, +Y, +Z order;(3) establishing a prior UC 6-connectivity / compounding table comprising a list of UCs where there are adjacent UCs to attain synergistic mechanics with one another, wherein the prior UC 6-connectivity / compounding table is planned by:Step 1 - importing a volume, wherein the volume is determined using HR- pQCT or microCT volumetric data;Step 2 - discretizing the volume;Step 3 - labeling the discrete units as the Unit Cells (UC) increasing in order of X+, Y+, Z+ to create a UC label volume;Step 4 - padding the UC label volume with a layer of zeros;Step 5 - create three variants respectively shifted in the +1 X, +1 Y, and +1 Z directions;Step 6 - trimming off the shifted cells outside of the original padding domain;Step 7 - flatten the sets and only keep the labels that are located where valid UCs are in the initial label volume, wherein this directly gives the -X, -Y, and -Z adjacent UCs (else, 0) for a given valid UC;Step 8 - store the valid UC indices in a column vector to the table (Column 1);Step 9 - storing column vectors for each directional adjacency as given by Step 7, zeros allowed to reflect a non-adjacency (Columns 2 - 4); andStep 10 - by row, computing adjacent cell count, N, by counting the number of nonzeros in the row across Columns 2 - 4 (Column 5).(4) entering a sequential design loop, comprising:(4.1) with insight from Columns 1-4 of the UC 6-connectivity / compounding table, in all compounding axes, θ, obtain:current UC design variables, a set of adjacent prior UC design variables andcompound UC homogenized elasticity matrices, wherein atT and there is neither concatenation nor de-concatenation of bone nor design;(4.2) perform compound inverse homogenization (cIH )on all N+l designs simultaneously, wherein N is the number of adjacent prior UCs per the UC 6- connectivity / compounding table, Columns 2-4, given in Column 5 directly;(4.3) store and, where * indicates final optimized design;(4.4) determine is the r+1 UC exist from a next row of the UC 6- connectivity / compounding table, Col. 1 , next row;If it does not exist, then all UCs have been designed, exit, marking termination of algorithm;If it does exist, then progress to next row, progressing UC index, T, in order of +X, +Y, +Z; and repeat steps 4.1-4.4.

12. The method of claim 11, wherein discretizing the volume of the microstructure volume comprises forming approximately 4mm side length Unit Cells (UC).

13. The method of claim 11, wherein performing compound inverse homogenization (cIH )on all N+l designs simultaneously comprises:Step 1 - initializing design variables, andStep 2 - via a convolutional compactly supported radial basis function (CS-RBF) collocation, yield an LS function (LSF), with sub-Gauss point accuracy;Step 3 - homogenize the Heaviside evaluation (pseudo binarization) of , solving fordisplacements at each Degree of Freedom (DoF);Step 4 - Tf Percent Total Intermediate material (PTI) is less than a critical value (i.e., 1 %), Exit, wherein PTI is largely governed by A which is initially large;Step 5 - compute a Jacobian (gradient based first order design sensitivity) and optionally a Hessian diagonal (gradient based second order design sensitivity);Step 6 - Update a and A, optionally using a Second Order Method of Moving Asymptotes (SOMMA), or a standard first order MM A algorithm or similar gradient based structural optimization / optimization algorithm(s), wherein if a symmetric Mean Absolute Percent Error (sMAPE) is less than a critical value (i.e., 20%), adjust SOMMA / MM A parameters to shift focus from reducing constitutive residuals to reducing PTI; andStep 7- repeat Steps 2-7.

14. The method of claim 13, wherein the convolutional compactly supported radial basis function (CS-RBF) collocation comprises: producing a compactly supported elemental list of Gauss points relative to the origin; evaluating a 3D Wendland3,l function, co, at these points, yielding a Wendland3,l Gauss FE integration point volume; convolve the FE integration point volume with a periodic version of the a field, comprising: computing nD (multidimensional) FFT of convolutional operand, here a, and convolutional filter, here co; elementally multiplying the two together; and computing an inverse nD FFT of the result and remove extra padding to arrive at the solution to the convolution.

15. The method of claim 13, wherein homogenizing the Heaviside evaluation (pseudo binarization) of cp, solving for displacements at each Degree of Freedom (DoF), comprises: via convolution, analytically evaluating cp; solving for loading and stiffness matrices in each element by direct Gaussian-Integration on approximate Heaviside of cp; store coalesced arrays of DoF stencil coefficients, taking advantage of gridded regularity; and solve the DoF displacements using a matrix-free implementation of Krylov subspace solvers in GPU.

16. The method of claim 15, analytically evaluating cp comprises: producing a compactly supported elemental list of Gauss points relative to the origin; evaluating a 3D Wendland3,l function, co, at these points, yielding a Wendland3,l Gauss FE integration point volume; convolve the FE integration point volume with a periodic version of the a field, comprising: computing nD (multidimensional) FFT of convolutional operand, here a, and convolutional filter, here co; elementally multiplying the two together; and computing an inverse nD FFT of the result and remove extra padding to arrive at the solution to the convolution.

17. The method of any of the above claims, wherein spatially mapping linear elastic anisotropy of the microstructure comprises:Step 1 - compute a list of local finite element matrices, as in:Step 1.1 - evaluate cp via convolution;Step 1.2 - perform numerical integration of local stiffness matrices and load vectors (unitary residual strain) via Gaussian Quadrature, wherein several Gauss points are considered, so that heterogeneity due to cp can be captured;Step 2 - from the list of local force vectors, assemble global force vectors;Step 3 - copy the list of local stiffness matrices and global force vectors from CPU to GPU (host to device);Step 4 - solve the systems of equations for each load vector using Krylov subspace solvers, wherein the Jacobi-preconditioned PCG method is a straightforward and efficient option of solver, which depends on:Step 4. 1 - Implementing sparse matrix-vector products in GPU as Node-by-node operators that recompute global stiffness coefficients on-the-fly, from the neighboring local matrices around each node;Step 4.2 - Implementing daxpy and reduce operations in GPU;Step 5 - copy the resulting displacement vectors back, from GPU to CPU (device to host);Step 6 - computing resulting stress fields from the displacement vectors, and integrate these stress fields to obtain effective stress at the macroscale; andStep 7 - computing the effective elasticity tensor from the obtained effective stress fields.

18. The method of claim 17, wherein evaluating cp via convolution comprises: producing a compactly supported elemental list of Gauss points relative to the origin; evaluating a 3D Wendland3,l function, co, at these points, yielding a Wendland3,l GaussFE integration point volume;convolve the FE integration point volume with a periodic version of the a field, comprising: computing nD (multidimensional) FFT of convolutional operand, here a, and convolutional filter, here co; elementally multiplying the two together; and computing an inverse nD FFT of the result and remove extra padding to arrive at the solution to the convolution.

19. A system comprising one or more processors in communication with a memory, wherein the one or more processors execute computer-executable instructions stored in the memory, said computer-executable instructions stored in the memory causing the one or more processors to implement the steps of any of the above method claims.

20. A non-transitory computer-program product comprising computer-executable instructions stored on the computer-program product, said computer-executable instructions for causing one or more processors to implement the steps of any of the above method claims.

Images