3-D phononic crystals with ultra-wide band gaps

In this paper gradient based topology optimization (TO) is used to discover 3-D phononic structures that exhibit ultra-wide normalized all-angle all-mode band gaps. The challenging computational task of repeated 3-D phononic band-structure evaluations is accomplished by a combination of a fast mixed variational eigenvalue solver and distributed Graphic Processing Unit (GPU) parallel computations. The TO algorithm utilizes the material distribution-based approach and a gradient-based optimizer. The design sensitivity for the mixed variational eigenvalue problem is derived using the adjoint method and is implemented through highly efficient vectorization techniques. We present optimized results for two-material simple cubic (SC), body centered cubic (BCC), and face centered cubic (FCC) crystal structures and show that in each of these cases different initial designs converge to single inclusion network topologies within their corresponding primitive cells. The optimized results show that large phononic stop bands for bulk wave propagation can be achieved at lower than close packed spherical configurations leading to lighter unit cells. For tungsten carbide - epoxy crystals we identify all angle all mode normalized stop bands exceeding 100%, which is larger than what is possible with only spherical inclusions.

design both phononic band-gap materials and structures. Gazonas et al. 30 and Bilal and Hussein 31 implemented a genetic algorithm based topology optimization method for the design of phononic band gap structures. Alternate structure types and materials exhibiting band-gap phenomenon have also been investigated. For example, Jensen 32 considered mass-spring structures, Diaz et al. 33 designed band-gap grid structures, Halkjaer et al. 34 maximized band gaps in plated structures, Olhoff et al. 35 and Halkjaer and Sigmund 36 optimized band-gap beam structures. Additionally, Vatanabe et al. 37 maximized phononic band gaps in piezocomposite materials, Liu et al. 38 explored the solid-solid phononic crystals for multiple separate band gaps with different polarizations, and Hedayatrasa et al. 39 optimized tunable phononic band gap plates under equibiaxial stretch.
Despite the considerable attention that topology optimization for 2-D phononic crystals has received, no work has been done on the topology optimization of 3-D phononic crystals. This is despite the potentially more useful nature of 3-D designs. Optimizing 2-D phononic crystals results in plate type designs which can have impacts on applications where wave propagation is constrained in 2-dimensions. However, wave propagation is inherently a 3-D phenomenon and optimization in 3-D can result in bulk materials which control wave propagation in all directions. This task is complicated by the challenging fact that phononic band structure evaluations are computationally expensive and that the computational complexity increases when band structure calculations are conducted repeatedly during the iterative process of optimization. The solution thus requires, first and foremost, an efficient phononic solver.
At this point there exist several numerical techniques for the evaluation of the phononic band structure. A good reference that discusses some of the most prominent techniques was published by Hussein 40 where the authors also presented a method of accelerating the existing algorithms through a secondary expansion. The Plane Wave Expansion 41 method (PWE) and the Finite Element method 42,43 are two of the most commonly used solvers owing to the ease of their implementation and their versatility. In this paper we have used a mixed variational method 44,45 to calculate phononic band structures. The mixed variational method is derived from the Hu-Washizu 46,47 variational theorem and it admits variations on both the stress and displacement fields. The mixed method has been known to converge faster than Rayleigh quotient which forms the basis of the traditional displacement based Finite Element method 48 . In a recently published comprehensive study we have shown that the mixed method also displays faster convergence than the PWE method 45 . In addition to using the mixed variational method as our solver we have achieved further computational accelerations by implementing it over distributed Graphical Processing Units 49 .
In this paper we have considered three main varieties of the cubic phononic crystal lattice (FCC, BCC, and SC). Our aim is to find 3-D topologies in a 2-material phononic crystal system that produce large all-angle, all-mode, normalized band gaps for each of the three symmetries considered. We evaluate the phononic band structures along the Irreducible Brillouin Zones (IBZ) 50,51 of the respective unit cells. The calculations are distributed over four compute nodes of a CPU-GPU hybrid cluster. We use a SIMP based topology optimization routine which is coupled with Heaviside projection for control over minimum feature sizes 52 . The sensitivity analyses for the eigenvalue problem required for TO are also calculated in parallel and through a vectorization process.

Results
Phononic band structure calculation. In the following calculations, the elastodynamic eigenvalue problem is formulated using the mixed variation method (Refer to Lu 45 for details). The propagation of waves in a three dimensional elastic medium is governed by j k jkmn mn ( , ) where λ = ω 2 , σ and u are the space and time dependent stress tensor and displacement vector respectively, ρ is the mass density and D is the compliance tensor. The Latin indices vary from 1 to 3 and subject to the summation conventions unless otherwise indicated. By varying independently on the stress and displacement field and enforcing Bloch periodic boundary conditions, (1) and (2) renders to the following functional stationary: and Q i are coordinates of the wave-vectors expressed in reciprocal lattice and T transforms the orthogonal coordinate system to the primary lattice of the unit cell. By substituting test functions into the mixed variation formulation (3) and setting the derivative of λ uσ with respect to the unknown displacement and stress coefficients equal to zero and then eliminating the stress coefficients through matrix manipulation, we obtain the following general matrix form of the eigenvalue problem V jkmn l l l 1 2 3 If M trigonometric expansion terms are used, i.e. α, β, γ vary from − M to M, then the size of the eigenvalue problem (4) will be 3(2M + 1) 3 × 3(2M + 1) 3 after considering the tensorial symmetries involved. Generally the mixed variation result will be neither upper nor lower bound of the eigenvalue solution 53 . The above method has been proven to converge faster to the real solution, in terms of the matrix size and corresponding relative error, than typical band structure algorithms such as the FE method, the Rayleigh quotient, and the PWE method 45,48 . Detailed convergence studies were presented by Lu and Srivastava 45 . Comparison was made in terms of the number of basis terms (trigonometric or real) needed in the approximate expansion. In summary, it was found that the mixed-variational method results in a higher accuracy for a given size of the eigenvalue problem.
Computation complexity and efficiency. Phononic band structure computation accuracy directly influences the objective function evaluation and its efficiency determines the tractability of topology optimization implementation. When the number of trigonometric terms, M, increases in the mixed variational formulation, the band structure shows decreasing relative error 45 . Therefore, we prefer to use a large M while keeping the matrix sizes manageable. In this study, 1029 trigonometric terms (M = 3) are used to compute the band structure. The related matrices Ω, Φ and H are of size 1029 × 1029, 2058 × 2058 and 1029 × 2058 respectively. Band structure evaluation has to be executed for 80 wave-vectors during each optimization iteration. The largest design domain in this study is a 48 3 mesh, which gives rise to 110,592 elements. Gradient-based optimization requires sensitivity analysis, or specifically the derivatives of Ω and Φ with respect to the elemental design variables indicating material concentration. This leads to 4.26 TB float type data which must be manipulated during each iteration. The computational cluster used for this work consists of 4 compute nodes, each of which has 4 NVIDIA GTX-780 graphic cards and 2 Intel(R) Xeon(R) E5-2630 v2 CPUs installed to form the mixed CPU-GPU architecture. Each GPU has 2304 CUDA cores and each CPU has 6 cores. In order to determine the parallel computation efficiency, we define an efficiency factor 49 which measures the performance improvement through the parallel computations over serial computations in terms of the time it takes to do the same problem through the two methods:  Fig. 1(b). As the number of trigonometric terms increases, the computation is 1000 times faster than the straightforward loop implementation. For a 2-D case, when M = 3, the matrix size is 98 × 98 and the parallel computation time is about 0.05 s when using 1000 elements. When using M = 3 in a 3-D case, the computation will be 263 times more complex, due to the complexity of eigenvalue problem being O(N 2.37 ), where N is the size of the matrix. Therefore, each eigenvalue problem takes about 13.2 s to be solved and it takes about 65 s to compute 80 eigenvalue problems in parallel on the 16GPUs to evaluate the necessary band structure.
A note on normalized band gaps. Band structure is a plot of the phononic eigenfrequencies for wave-vectors which span the boundaries of IBZ. On this plot the normalized band gap is calculated by taking the ratio between the band gap width and the mid gap frequency. This metric is independent of unit cell scaling and is indicative of the tendency of a phononic crystal to stop all waves in all directions (for 3-D). The simplest way to control normalized band gap sizes is to change the volume fractions of the material phases. Consider, for instance, simple phononic crystals made of tungsten carbide inclusions (ρ = 13800 kg/m 3 , E = 387.5559 GPa, ν = 0.3459) in epoxy matrix (ρ = 1180 kg/m 3 , E = 4.3438 GPa, ν = 0.3679). 1-D phononic crystals are layered composites (Fig. 2a). By changing the thickness of the material phases the maximal normalized band gap is seen to appear at a tungsten carbide volume fraction of 0.737, at which point the normalized band gap size is 155.5%. Since thickness is the only design variable in this 1-D case this maximal gap size is also the global optimum for the material choices. Some optimization work has been reported for obtaining specific phononic bandstructure of interest, such as design of phononic filter 54 and solving for transient wave propagation problem 55 . Circular inclusions inside a 1 × 1 square unit cell are considered for a 2-D phononic crystal in Fig. 2b. We observe that the maximum normalized band gap for this simple geometry has a value of 113.6% and it occurs at a volume fraction of 0.624. Several studies have been conducted on exploring 2-D phononic topologies for larger bandgaps. For example, Bilal and Hussein 31 and Liu et al. 38 have both achieved larger than 120% normalized band gaps in 2-D. Figure 2(c-e) show analogous 3-D phononic cystals with 1 × 1 × 1 cubic unit cells. The spherical inclusions lead to maximum normalized band gap values of 67.5%, 94.2% and 93.3% at volume fraction 0.412, 0.545 and 0.588 for SC, BCC and FCC lattices respectively. We further note here that it becomes progressively more difficult to obtain large band gaps as we consider crystals of higher space dimensions.
The 1-D phononic crystal case is simple enough to admit theoretical arguments on largest possible bandgap values 56 and the 2-D case received research interest in terms of topology optimization studies. However, the 3-D phononic crystal band gap optimization has not yet been reported in literature. It is clear that formal topology optimization should be able to reveal crystals with larger stop bands than those produced by the simple inclusions in Fig. (2). Topology optimized structures. For our study, topology optimization is performed over a 1 cm 3 cubic unit cell consisting of tungsten carbide and epoxy phases, allowing us to compare our results with previously published results 57 of ultra-large band gap phononic crystals. The topological variable ρ e indicates the relative volume fraction of epoxy in each element, with ρ e = 0 indicating the element contains only tungsten carbide and ρ e = 1 indicating only epoxy. During each iteration the stress and displacement fields are expanded using 1029 trigonometric terms (M = 3). After the optimized solution is obtained a larger number of expansion terms (M = 4) is used to calculate the final band structure at a higher accuracy. The objective function in this study is the normalized band gap between the 6th and 7th band, where a complete band gap opens naturally for a simple inclusion within a primitive cell. The longest single optimization iteration for the largest design domain, which consists of 48 3 elements without imposing symmetry, takes only 6.5 min with 22.5% of the time spent on band structure calculation and 45.5% on sensitivity analysis. Design step direction calculation and data storage takes 32% of the time in this case, however, if a smaller design domain is used, then this proportion will be smaller.
It should be noted that the considered topology optimization problem is nonconvex and thus the local minimum identified by the gradient-based optimizer is dependent on the topology used as the starting point for the optimization. Therefore, different material distributions are used as initial designs to help avoid getting trapped in low performance local minimum. Specifically, we have tested using spherical inclusions of various volume fractions and various homogeneous mixtures of tungsten carbide and epoxy as initial guesses. Of course, as in most topology optimization problems, there is no guarantee that we have been able to identify a global optimum.
Simple cubic lattice. As a first example, we consider the SC case defined on a relatively coarse mesh of 24 3 elements with basic symmetries employed. This means the optimization is performed over 1/8 of the total number of elements with reflection symmetries assumed to generate the rest of the SC primitive cell. These initial results are presented in Fig. 3 for different homogeneous initial distributions of material. When we take into account the periodicity of the unit cell, Fig. 3(a-c) show that for different homogeneous initial mixtures of tungsten carbide and epoxy, the optimization process converges to approximately the same topology after about 500 iterations. These solutions have similar normalized band gaps, which are 71.58%, 78.47% and 70.47% corresponding to tungsten carbide volume fractions 0.4374, 0.5015 and 0.4237, respectively.
The SC problem using symmetry and a uniform initial distribution of ρ e = 0.5 was re-solved on a finer mesh of 48 3 elements. The optimization is performed over only 1/8 of the total number of elements due to imposed symmetries. The size of the data which is needed to be manipulated at each iteration is 545 GB. The design evolution for this finer mesh case is shown in Fig. 4 and it is seen the algorithm evolves towards a single inclusion structure as in the previous coarse mesh cases. The optimization converges after about 1000 iterations and the optimized inclusion structure is a cross between a sphere and a cube with some variation on the surface (Fig. 5a). The structure has a normalized band gap ratio of 67.7% and tungsten carbide volume fraction 0.3907.
Since it is seen in both the coarse and fine mesh cases that the optimization process tends to evolve the topology into a structure that is equivalent to a single inclusion within a unit cell, we also used two variations of a centrally positioned stiff spherical inclusion as initial conditions for the optimization. The first case corresponds to a spherical inclusion with the largest volume fraction possible and the second corresponds to the one with largest normalized band gap (Fig. 2c). The latter case resulted in a final structure with larger normalized band gap after convergence (about 1450 iterations). Figure 5 show the results for the optimized SC design. The geometry of the inclusion, which has volume fraction 0.4771, resembles a cube with rounded corners. The corresponding stop band extends from 63.18 kHz to 143.17 kHz which is equivalent to a normalized band gap ratio of 80.1%. This is 18.67% larger than the maximum possible SC band gap with a fully dense spherical inclusion.
Finally, we note that we also attempted to relax the imposed symmetries and solve the SC case using an initial homogeneous distribution of ρ e = 0.5 in mesh of 48 3 elements. This requires manipulation of 4.26 TB float type data during each iteration and was thus significantly more computationally intensive than the previous case. The optimization tended towards highly asymmetric structures with non-existent or very small stop bands after several design iterations. This is likely due to the optimization process getting stuck in local minimas. The SC optimization problem without any additional symmetries is a very large computational problem with > 100,000 variables. To adequately explore it we require faster computational algorithms and we expect that relaxing the symmetries will indeed result in larger bandgaps. However, current computational resources do not permit us to study this problem adequately.
Body-centered cubic lattice. The BCC primitive cell is discretized into a 36 3 mesh giving rise to a total of 46,656 elements. Optimization is performed without any additional symmetry constraints. Figure 6 shows the design evolution, starting from a homogeneous initial design of ρ e = 0.5, progressing to a skewed cross-like structure, before converging after 186 iterations to a centrally-located large inclusion with smaller inclusions located near the unit cell corners as shown in Fig. 6(f). However, given the periodicity of the unit cell, the topology still corresponds to a single inclusion with the smaller inclusions joining with the large inclusion at appropriate locations in periodically repeated unit cells. The optimized structure features staggered inclusions which have the general appearance of a sphere but with some variations on the surface. The structure has normalized band gap ratio of 93.11% with a tungsten carbide volume fraction 0.4417.
We also considered two cases whose initial designs are spherical inclusions. As in the SC case, the first case corresponds to a spherical inclusion with the largest volume fraction possible and the second corresponds to the one with largest normalized band gap (Fig. 2d). The former case resulted in a larger normalized band gap, converging after 287 iterations. Figure 7 shows the optimized BCC design. The optimized BCC primitive cell contains the complete geometric information of the single inclusion. This topology is then assembled to result in its corresponding unit cell by repeating the primitive cells based on the built-in translation symmetry of the BCC lattice (Fig. 7(b)) Much like the solution reported in Fig. 6(f), the resulting topology resembles a staggered pattern of tungsten carbide inclusions which have the general appearance of spheres. However, large portions of the surfaces that face the body diagonals in these inclusions are flattened, as shown in Fig. 7(c). Although no internal symmetry constraint has been applied during the optimization process, the optimized inclusion structure has a centro-symmetric element arrangement. The optimized BCC structure has a band gap from 76.8 kHz to 230.49 kHz, which leads to a normalized band gap value of 100.03%. The corresponding volume fraction of tungsten carbide in the optimized structure is 0.5475. This is in contrast with spherical inclusion results where the maximum possible normalized band gap is 94.2% at a volume fraction of 0.5450.
Face-centered cubic lattice. The FCC primitive cell is discretized into a 36 3 mesh and optimization is performed without additional imposed symmetry constraints. The FCC case is solved using the same initial guesses as the BCC case: a homogeneous distribution with ρ e = 0.5 as well as a two spherical inclusion designs. As in the BCC case, the topology optimized FCC solution found using a uniform initial distribution of material converges to a solution featuring a single inclusion shown in Fig. 8(a), with the general shape being close to a sphere. It has a normalized band gap of 95.24% at a volume fraction 0.4861. The solution found using an initial spherical distribution corresponding to the largest band gap in Fig. 2(e) gave the best result, converging after 235 iterations to the optimized primitive and corresponding topologies shown in Fig. 8(b,c), respectively. Examining Fig. 8(c), it is clear that the shape of the tungsten carbide inclusion has the general appearance of a sphere with some variation on the surface. Although no internal symmetry constraint has been applied during the optimization process, it     is seen that the optimized inclusion structure has a centro-symmetric element arrangement. The optimized FCC structure has a stop band from 98.03 kHz to 292.55 kHz, which leads to a 99.61% normalized band gap. The corresponding volume fraction of tungsten carbide in the optimized structure is 0.5502.
According to Economou and Sigalas 58 the "cermet" topology, where inclusions consist of low-velocity materials surrounded by high-velocity matrix materials, is more favorable for the appearance of large elastic gaps. However, we could not find support for this in our studies. We found that beginning with cermet initial designs (Fig. 9a), where an epoxy sphere is embedded in tungsten carbide matrix, the optimization process invariably steered towards network topologies in search of large band gaps. By step 6, the optimization progression has already converted the cermet initial design into a network design where a tungsten carbide inclusion is surrounded by epoxy matrix (Fig. 9a). The optimization process converges after 134 iterations. The optimized structure has volume fraction 0.4752 and the stop band extends from 94.84 kHz to 265.48 kHz, which leads to 94.71% normalized band gap. Furthermore, in all our other optimization studies which involved beginning with a homogeneous initial design, the optimization process could have steered towards a cermet topology. However, it always resulted in network topologies.

Discussion
In this paper we have presented the first ever topology optimization results for 3-D phononic crystals. Specifically our objective was to reveal 3-D phononic unit cells comprised of two material phases which display large all angle, all mode phononic band gaps. Specifically, optimized results for simple cubic, body centered cubic, and face centered cubic crystal structures made up of tungsten carbide and epoxy phases are presented. We have shown that for all these cases large phononic stop bands for bulk wave propagation can be achieved at lower than close packed configurations for spherical inclusions. A summary of the band gap results can be found in Table 1, where we have compared our results with what is achievable through simple spherical inclusions. Specifically, it is possible to achieve normalized band gap of 67.5%, 94.2% and 93.3% using tungsten carbide spheres in SC, BCC and FCC configurations respectively. Topology optimization shows that the SC result can be significantly improved to 80.1% by modifying the shape of the inclusion, forming a cross between a sphere and a cube. The BCC and FCC optimization results also improve over their spherical inclusion counterparts. Specifically, the BCC optimized structure shows a normalized bang gap greater than 100%. Furthermore, it is interesting that the optimized structures shown in Figs 7(a) and 8(b) achieve large complete band gaps at relatively small volume fractions of the stiff and heavy inclusion. For instance the BCC optimized structure achieves a large normalized band gap at a tungsten-carbide volume fraction of 0.5475. This is in comparison with a volume fraction of 0.74 which corresponds to the case where spheres are tightly packed in a FCC configuration. This results in a unit cell which is 23.1% lighter than the close packed structure and still outperforms it in terms of its band gap size 57 . In  general, it is notable that largest normalized band gaps occur at less than close packed configuration both for spherical inclusions and optimized results. 3-D phononic band structure optimization is a highly computationally intensive task, enabled here through a GPU accelerated mixed variation method and mixed variation formulation based sensitivity analysis. For example, the 3-D SC lattice using a 48 3 mesh translates into an optimization problem with over 100,000 variables. Furthermore, with no assumed symmetry, the optimization process tends to get stuck in numerous local minima which correspond to highly asymmetric structures with small band gaps. For the SC case, therefore, we assume reflection symmetries over 1/8 of a unit cell. We studied various initial conditions with coarse (24 3 ) and fine (48 3 ) meshes. It was notable that almost all homogeneous initial conditions result in single contiguous inclusion topologies which are very similar to each other. All optimization runs which begin with spherical inclusions as their initial conditions also ended in simple contiguous inclusion topologies. This result was also noted for FCC and BCC runs. For both of these cases we used a 36 3 mesh configured along the appropriate unit vectors. No further symmetries were assumed. In general, it was noted that homogeneous initial condition resulted in inclusions with more pronounced asymmetrical features. On the other hand, the cases with spherical inclusions as initial conditions resulted in more symmetrical inclusions with higher normalized band gap values. In our studies we could not find cases which support that cermet topology is more favorable in generating elastic band gaps. In fact, beginning with cermet initial designs, the optimization process invariably steered towards network topologies in search of large band gaps.
In summary, we have presented topology optimization results which reveal the largest as yet reported all angle all mode normalized bandgaps in 3-D phononic crystals. These results pertain to the two material combination of tungsten-carbide and epoxy. Other material combinations and/or optimizing band gaps between other bands (other than the 6th and 7th) will likely result in different optimized topologies and different normalized bandgaps.

Methods
Topology optimization formulation. In topology optimization of periodic composite materials, the goal is to optimize the distribution of two (or more) base material phases across the unit cell, which for finite element-based approaches, reduces to determining whether each element is to contain base material 1 or material 2 (see e.g. refs [14][15][16][17]. This fundamentally is a binary (or integer) programming problem of extremely high dimension, motivating relaxation of the binary condition and representation of each element's material properties as a continuous combination of the two base materials. In order to obtain a binary design, penalization methods such as the Solid Isotropic Material with Penalization (SIMP) method are used to make mixtures of the two materials at a location inefficient. Interestingly, Sigmund and Jensen 29 found that the use of penalization is not required in the design of band-gap structures as sharp contrasts in stiffness (Young's modulus) are desirable to produce large band-gaps. This allows use of a simple linear interpolation model for Young's modulus, given as where E 1 and E 2 denote the Young's modulus corresponding to the two base materials. We note this is equivalent to the SIMP interpolation for composites 59 with exponent penalty term set to one. The goal of the optimization is to maximize the gap between the t-th mode and (t + 1)th mode, and thus we use the objective function given as e e 1 Note that we do not impose a volume constraint and allow the algorithm to freely distribute the two base materials, although it would straightforward to restrict the design problem in that manner.
We want to emphasize that Eq. (11) generally indicates an asymmetric eigenvalue problem. As a result the eigenvalues and corresponding eigenvectors are complex. Sensitivity calculations of complex eigenvalues require the use of normalized left and right eigenvectors, ψ and θ, in the sense that for each eigenvalue the corresponding pair of left and right eigenvectors should satisfy ψ * θ = 1, such that Although the closed form expressions for Ω −1 and Φ −1 are implicit, based on the invertibility their derivatives with respect to design variable ρ e can be written as, The sensitivity of the natural frequencies 29 can be calculated by differentiating (12) with respect to the design variables ρ e , as follows The sensitivity of the objective function can now be calculated by differentiating (10). Finally we note that the Heaviside Projection Method (HPM) 52 is used within these formulations to control the minimum length scale of designed features. In particular, a reduced design variable field 60 is adopted with design variables spaced at two times the finite element size. In HPM, the continuum design ρ e variables are expressed as a closed-form function of an independent design variable field without any other changes to the above equations. For brevity, the details are omitted here; however, the reader is referred to Guest 61,62 for full algorithmic details. Here, a continuation strategy is used on the Heaviside parameter, beginning at β HPM = 1. The projection radius is 1.4% of the unit cell length. It should be noted that HPM is capable of controlling the minimum length scale of stiff and/or compliant designed features and thus preventing solution mesh dependency, although such dependencies have not been observed for phononic band-gap materials 29 .
Vectorization and parallel computations. In order to calculate the entire band structure of a unit cell the matrices have to be assembled and the eigenvalues have to be calculated at multiple wave-vector points along the edge of the IBZ. This results in considerable computational complexity. However, since the assembly and eigenvalue solving processes are independent of each other they can be executed in parallel if the formulation is properly recast. The most basic computational unit in the formulation is the following integral in (5): . The integrand can be expanded as It is a constant matrix independent from the wave-vector coordinates Q i . The above equation can be viewed as the outer product of f and its own complex conjugate, where is a vector of size (2 M + 1) 3 . Now the integral (16) can be rewritten as where p, q = 1, 2, … , (2 M + 1) 3 . It is a constant global matrix which only needs to be calculated once. To compute these matrices using Graphical Processing Units we need to pass the vectors, volumes and centroids from the CPU to the GPU. On the GPU, the computation kernels are executed by a grid of thread blocks, where each thread has a unique id which corresponds to a set of indices e, p, q. Since the actual computation on each thread is relatively simple and many threads are operating in parallel, the method shows significantly reduced computation times compared to serial computations over a CPU. Furthermore, the band structure computations for different wave-vectors along the IBZ can be distributed over multiple GPUs in a distributed GPU cluster. In this case each compute node will solve only a part of the band structure thus decreasing the computation time further. Sensitivity is calculated element-wise and therefore, parallel computation can significantly accelerate the process. Substituting (19) into (6) and (7), the derivative of Ω and Φ with respect to the design variable of the e-th element, ρ e , can be calculated by e pq e e 1 2 ( ) ( )  64 . The discretization of the design domain should allow the formation of all possible variations during the optimization process, however, it is not very economic to use the entire cube as the design domain. In order to maintain the basic crystal structure during the optimization process, the geometry information has to be exactly the same on each lattice point. This can be realized by enforcing translation symmetries along their symmetry axes. As shown in Fig. 10(a-c), SC lattice's translation symmetry is along the lattice edges, whereas BCC and FCC lattice symmetry axes are along the body and face diagonals respectively. These translation symmetry axes are the vectors of the primitive cells which contains the geometry information of exactly one lattice point. Band structure calculations and sensitivity analyses are implemented on primitive cells. As demonstrated by Dong et al. 65,66 , reduction of symmetry is favorable in generating ultra-wide bandgaps, therefore, BCC and FCC primitive cells are discretized into 36 3 elements and no further symmetries are assumed for them (Fig. 10e,f). The primitive cell of SC lattice is the cubic unit cell itself and it is discretized into 48 3 elements. We had difficulty converging to a meaningful solution without imposing symmetry in the SC case where homogeneous material distribution is used as the initial design. Following Bilal and Hussein 31 , who implemented C 4v rotational symmetry to their 2-D square lattice optimization, we reduced the SC design domain to 1/8 of the primitive cell as shown in Fig. 10(d) and reconstruct the lattice by taking reflection symmetries, in which the reflecting mirrors are the orthogonal planes. Figure 10(g-i) show the IBZ boundaries of SC, BCC, and FCC lattices respectively in the reciprocal space. During topology optimization the band structure is calculated for wave-vectors spanning these boundaries and then band gap location and size are extracted.