Justification of a certain algorithm for shape optimization in 3D elasticity

The trabecular bone adapts its form to mechanical loads and is able to form structures that are both lightweight but very stiff. In this sense, it is a problem (for the Nature or living entities) similar to the structural optimization, especially the topology optimization. The structural optimization system developed by the authors previously was based directly on the bio-mechanical models. The paper aims to clarify the approach of using the bio-mechanical models of the trabecular bone remodelling phenomenon as a base for the structural optimization. The justified algorithm for shape optimization mimics the natural phenomenon to satisfy the condition of constant value of the strain energy density on the structural surface, issuing from the considerations in the area of mechanics. The main body of the paper contains the analysis of the stiffness optimization problem in the framework of speed method approach to shape optimization. The given numerical example (cantilever beam in bending) includes the influence of the structural surface curvature on the optimization procedure. The presented in the paper approach can be used as a method of structural optimization unrelated already to trabecular bone remodelling phenomenon.


Introduction
The field of structural optimization is still a relatively new field undergoing rapid changes in methods and focus. There are exciting applications of the methods of structural optimizations in the automotive, aerospace, civil engineering, machine design and other engineering fields. As a result of the growing pace of applications, research into structural optimization methods is increasingly driven by real-life problems. The major challenge faced by researchers in structural optimization is to develop methods that are suitable for numerical solution of shape-topology optimization problems. Another major challenge is the high computational cost associated with the analysis of many complex real-life problems. In the paper we combine the method inspired by bone evolution with the shape gradient framework for minimization of the energy type functionals in 3D elasticity. In this way the computational cost of the proposed method becomes reasonable for 3D elastic structures.
The classical approach to structural optimization can be found e.g., in monographs (Arora 2015;Haftka and Gürdal 1992). The mathematical analysis of the so-called velocity method in shape optimization is performed, e.g., in Sokołowski and Zolesio (1992).
It is widely believed that the trabecular bone adapts its form to mechanical loads, in accordance with the Wolff's law (Wolff 1892). Healthy tissue of trabecular bone has a very sophisticated structure. The tissue forms a network of beams called trabeculae and this structure is capable of handling a wide range of loads. The length of the trabecula can be 100 or 200 μm, whereas its diameter is approximately 50 μm. This structure is continually rebuilt so that the whole bone tissue is replaced during the course of few years. The process is called trabecular bone adaptation or remodelling.
Wolff's law is not indeed a mechanical law but rather a hypothesis, that the trabecular bone is a self optimizing material. Also others researchers (Cowin et al. 1976;Huiskes et al. 2000) suggested, that the trabecular tissue is able to adapt its structural form in reaction to the external load. In this sense, it is a problem similar to the structural optimization, especially the topology optimization. The model of the bone remodelling process discussed in the paper is based on the regulatory model of Huiskes. The main idea of the model consists of a regulatory mechanism (on the bone surface only) between bone resorption and formation. Bone mass increases above a certain level of mechanical stimulation (measured by the strain energy density on the tissue surface) and decreases below the other certain energy level. When the strain energy density on the structural surface is between these two levels the bone mass is maintained and this rate of bone mechanical stimulation is called adaptation or lazy zone, because then the bone does not react to the changes in mechanical stimulation level. The processes of resorption and formation are coupled in basic multicellular units, regions smaller than the single trabecula where the sequence of successive bone tissue resorption and formation occurs and this is a local change. However, bone formation depends on global mechanical stimulation of the whole bone.
The strain energy density on the structural surface appears also in the area of optimization research, distant from bio-mechanical studies (Wasiutyński 1960;Pedersen 2003) where one can find the theorem, that for the stiffest design the energy density along the shape to be designed must be constant. The similarity between trabecular bone remodelling and topology optimization is used in two opposite directions.
Some researchers use the optimization methods to recreate the real biological structural forms of bone. They use a variety of approaches and methods. Sigmund (1999) discussed the problem if the trabecular bone structure has optimal stiffness. He used SIMP method to recreate bone structures and concluded, that the trabecular bone structure does not have optimal stiffness. Wu et al. (2017) used also SIMP method to recreate bone like structures. The modified optimization method has been used and the modification were local volume constraints. Also Jang and Kim (2008) use topology optimization for the trabecular bone structure prediction. Lee et al. (2015) expanded this model and developed homeostasis-based aging model for trabecular changes. Marzban et al. (2015) uses the trajectorial architecture theory of optimization to predict the remodelling of material microstructure and structural organisation under mechanical loading. Goda et al. (2016) presented combined bone internal and external remodelling model relying on the viewpoint of configurational mechanics with unified treatment of volumetric and surface growth, based on Eshelby stress. Brampton et al. (2012) presented application of level set topology optimization to predict the trabecular bone architecture. There are also the variety of reasearch concerning the multiscale approach to trabecular bone remodelling. Coelho et al. (2011) presented two-scale topology optimization procedure for simulation of trabecular bone remodelling phenomenon.  and Fernandes et al. (1999) developed (using homogenisation theory) a global analytical parametric microstructural model for the remodelling phenomenon. Makowski and Kuś (2016) proposed the methodology of optimization patientspecific trabecular bone scaffolds. The optimization is performed on the basis of homogenized orthotropic elastic properties of trabecular bone tissue and three-scale numerical model.
The opposite direction of using the similarity between trabecular bone remodelling and topology optimization is to apply the bio-mechanical observations and models to the structural optimization issues. On such an assumption the structural optimization system based on trabecular bone surface adaptation was developed (Nowak 2006). The domain independence, functional configurations during the process of optimization and possibility to solve multiple load problems are the unique features of the biomimetic method, useful in mechanical design. Additionally, such an approach allows to comprise the optimization of size, shape and topology in one numerical procedure, where the lazy zone concept is an important component of the algorithm. Klarbring and Torstenfelt (2012) discusses the connection between apparent density-type bone remodelling theories and density formulations of topology optimization and involves a lazy zone concept from bone remodelling theory to this optimization approach. Nutu (2015) gives an interpretation of bone remodelling models parameters, like the homeostatic equilibrium threshold, during the process of topology optimization. Florio (2015) developed the gradient-free shape optimization method based on bone adaptation model. The paper contains discussion of the different stress measures as a base for the remodelling scenarios. Author concludes, that use of different stress meassures does not influence the trends in shape changes predicted by these models. The similar discussion can be found even earlier in the paper by Huiskes (2000). He used also the strain energy density as a remodelling signal by a choice only and suggested, that other mechanical variables could be used as well. In our work we will try to show, that Huiskes choice although intuitive was appropriate.
The paper aims to clarify the approach of using the bio-mechanical models of the trabecular bone remodelling phenomenon directly as a base for the structural optimization (Nowak 2010). In general the problem of stiffest design (compliance minimization) has no solution. If the volume of the object is increasing, the compliance is decreasing. Thus, in the standard approach to the energy based topology optimization the additional constraint has to be added. Usually the volume of the material is limited. When using the volume constraints the additional stress constraints should be introduced. As a result the topology optimization has to be separated from the shape and size optimization within the numerical procedure (Krog et al. 2004). In case of the biomimetic approach presented here, the role of an additional constraint plays the strain energy density on the structural surface (which is between two assumed levels forming the lazy zone) and the volume or structural mass results from the optimization procedure. In the standard SIMP topology optimization method the volume constraint is used. Also in the stress-based topology optimization approach (Bruggi 2008;Le et al. 2010) the volume constraint is present. To start the optimization procedure the volume constraint has to be imposed. In biomimetic approach, instead of imposing volume constraint we parameterize shapes by the assumed energy density, which may be quite accurately predicted from the yield criteria. In the proposed numerical procedure the shape and topology optimization is performed without volume constraints.
The plan of the paper is as follows. First we formulate the problem and briefly describe the biomimetic algorithm used in Nowak (2010). Then the very short description of the speed method used in shape optimization is given, together with demonstration that all ingredients of the resulting formula for the shape derivative are significant. The shape derivative allows to obtain a straightforward justification of the well known necessary optimality condition for our problem. Then we reformulate the biomimetic algorithm as the shape optimization procedure and use speed method in order to obtain the shape gradient of an appropriate domain functional.
Finally, we use the part of the shape gradient containing the curvature of the variable part of the boundary improving the descent direction used in the gradient type method of optimization.
It is known, that in order to assure the convergence of optimization algorithm, the exact value of the steepest descend direction is not required, but only the direction of improvement. This is the case of purely biomimetic approach, which works. However, we may expect that using better approximation of the exact gradient should improve the performance of the optimization procedure, resulting in smoother trajectory in the design space and faster convergence.

Problem setting
The goal is to maximize the stiffness of a structure, that is minimization of the functional Here, Ω represents domain of the elasticity system, u the displacement, V 0 = |Ω 0 | a given volume, Γ 0 part of the boundary with Dirichlet condition, Γ 1 part of the boundary loaded by traction forces.
The heuristic algorithm (Nowak 2010) reads as follows: -it is assumed that the energy density σ (u) : ε(u) has a constant value λ on Γ v ; -if at a given point on Γ v this density is bigger than λ + s then the boundary in moved outside; -if at a given point on Γ v this density is smaller than λ−s then the boundary in moved inside; -these steps are repeated until equilibrium is achieved; -the value of λ is modified if the final design is unsatisfactory.
The developed finite element mesh generator was originally dedicated to mesh creation for biological entities. The basic input to the optimization system is the geometrical model of the structure. Since the visualization for the biological entities is based on the digital images e.g. Computed Tomography, the input to the system is also based on the collection of the 2-dimensional images. The images are obtained automatically by cutting a 3D CAD model, in a predefined direction, into a flat intersections. After preliminary graphical operations the images of slices are directly used for building of the 3-dimensional finite element mesh. The 2-dimensional image is first translated into the bitmap, where "0" represents the void and "1" the object. On the bitmap the initial step of discretization is executed. The aim of this first step is to describe the areas with material. The discretization procedure produces a 2-dimensional network of tetragonal elements, according to the object image shape.
The discretized 2-dimensional image is projected onto the subsequent one. If there are areas containing material on both images, the boxes are created. Each box is in turn translated into 6 tetrahedral volume elements. We divide hexahedral boxes into the 6 tetrahedral elements to be able to perform smoothing of the finite element mesh. We introduce additional tetrahedral elements using existing nodes in such a way that the number of equations remains the same, but the stress and strain distribution are more smooth. The volumetric mesh is ready to be used in computations of the distribution of the strain energy density. After the finite element method computations, for each box the information about the strain energy distribution is in turn associated with the description of the structural surface position on the 2dimensional images. The details of the mesh generation are depicted in Fig. 1.
In the biomimetic algorithm, if the strain energy density is less than the lower limit for the lazy zone, the small amount of material is removed from the surface. If the strain energy density is greater than the upper limit for the lazy zone, the small amount of material is added on the surface. The structural evolution is separated from the finite element method computations and modeled with use of 2dimensional images. If there is the need to add material on the surface, the additional layer of "1"'s, representing the material is placed instead of "0"'s (case A' in Fig. 2). If there is the need to remove material from the surface, the additional layer of "0"'s, representing the void is placed instead of "1"'s (case A" in Fig. 2). In this way the structural surface is moved in the virtual space and it is possible to simulate structural evolution including separation or merging of the structural elements, without problems with finite element computations. The range for lazy zone is needed if the algorithm is to give satisfying results in case of multiple objectives (e.g. several loadings (Nowak 2010)) and consecutive "alternate cases" modification steps, as well as for ensuring the stopping of the procedure. As we know, the exact solution without tolerance zone usually does not exists. The heuristic biomimetic procedure is based directly on the biomechanical model of the trabecular bone remodelling phenomenon. To compute the strain energy density on the structural surface the Finite Element Method is used. To compare the procedure to the standard topology optimization method, the typical example i.e. the bending cantilever beam was chosen. The biomimetic optimization system like it is in the case of living entities is 3-dimensional, and the definition of the design domain is not needed. Moreover, also the definition of the boundary conditions could be different for each step of the simulation. It is assumed for the presented example of the cantilever beam in bending, that all nodes of the one edge could be fixed (clamped wall) and the bending force is applied to the middle of the opposite edge. Material parameters were assumed as follows: Young's modulus 2 · 10 11 Pa, Poisson's ratio 0.3. The obtained results with the different starting configuration i.e. full domain and just a stick connecting the clamped wall (the potential surface, where the no displacement boundary conditions could be defined) are presented in Figs. 3 Fig. 3 Optimization results of the cantilever beam bending using of the biomimetic procedure, based on the trabecular bone surface adaptation -standard topology optimization domain full of material and 4 respectively. Obtained results are close to the well known solution. The final configuration are not unique due to existence of the insensitivity zone, different starting configurations and the possibility of choosing the boundary conditions position during the simulation based on the prescribed area. If the position of the boundary condition is fixed for the whole starting configuration the results have the same form (Nowak 2006).

Speed method for 3D elasticity
In this section we shall sketch the idea of the speed method, as used in shape optimization in 3D elasticity. The regularity conditions concerning the objects appearing in the Fig. 4 Optimization results of the cantilever beam bending using of the biomimetic procedure, based on the trabecular bone surface adaptation -empty domain presentation will not be discussed, we refer the reader to the monograph by Sokołowski and Zolesio (1992). We refer also to Delfour and Zolésio (2011) for further developments, and to Plotnikov and Sokołowski (2012) for the case of nonlinear problems governed by compressible Navier-Stokes equations. There are other approches like traction method (Shimoda et al. 1998(Shimoda et al. , 2012. However, the idea of traction method is related to a specific approach to the parametrization of the variable domain by means of elastic deformation. Such a parametrization does not allow for topology changes. Let V(x) be a smooth vector field defined on IR 3 and T t (x) a transformation defined by the formula Let also Ω 0 ⊂ IR 3 be a certain domain and Ω t its image by means of T t (•), that is Ω t = T t (Ω 0 ). Observe that T 0 (Ω 0 ) = Ω 0 . Next we take a function u t (x) defined on Ω t (here u t (x) does not have to be a displacement) and dependent on Ω t explicitly or for example by means of the boundary value problem posed in Ω t .
Using this notation we define in Ω 0 the shape derivative of the function u t with respect to the (fixed) velocity field V(x): The shape derivative can be evaluated by using an appropriate integral functional and it should be distinguished from material derivative, not used in further presentation, which is defined aṡ There is a relation between these two objects, namely which is not employed it the paper. The shape derivative is used for analysing the behaviour of domain functionals depending on both Ω t and u t . Let us consider the shape functionals in the form of volume or surface integrals and, denoting by Γ t = ∂Ω t , If we define shape derivatives of these functionals in the direction of the field V(x) as then it can be proved (Delfour and Zolésio 2011;Plotnikov and Sokołowski 2012), that Here κ(x) denotes an average curvature of the surface Γ 0 at the point x, namely hence Ω t = B(0, R + t) and V(x).n(x) = 1 on Γ 0 . The variable r denotes polar coordinate around 0. As function depending on Ω t we take and as F 1 and F 2 the function F (u) = u 3 . Then it is easy to compute Now we shall use the derived formulas. After obvious calculations As a result the a 1 + b 1 = J 1 , as expected. Similarly, Again the result a 2 + b 2 + c 2 = J 2 is in agreement with appropriate formulas. What is important, we see that no part of them can be neglected without causing discrepancy.
The effect of curvature may be illustrated in the even more straightforward way. Let us assume, that we have the integral over the sphere S Then, assuming that we want to increase radius, V.n = δR, On the other hand, and it is indeed the average curvature, if 1/R 1 , 1/R 2 denote principal curvatures at the point on the sphere.

Justification of assumption about constant energy density
Let us define the Lagrangian for the problem under considerations and rewrite the state equation in the weak form Here u t , ϕ ∈ H 1 Γ 0 (Ω t ) and Ω t = T t (Ω). For fixed field V(x) which vanishes on all the boundary of ∂Ω with the exception of Γ v , the Lagrange function depends only on two scalar variables (t, λ). Then we take shape derivative of both Lagrange function and weak state equation using speed method. At the local minimum this first derivative should vanish. The shape derivative is denoted here by (•) .
After substituting ϕ := u in the weak state equation and ϕ := u in it's shape derivative we get Using this result in the derivative of Lagrange function gives Since at the stationary point this should hold for any vector field V(x) on Γ v , then Of course we do not know the value of λ. Similar results were obtained previously by others researchers (Pedersen 2003). Presented here result is based on the shape derivative concept and does not need any additional assumptions.
The assumed value of the strain energy density on the part of the boundary subject to modification could be related to the material properties. Change in the assumed value of the strain energy density results in change of the structural form -topology and volume. In this way, the final structural volume results from the optimization procedure. Instead of imposing volume constraint we parameterize shapes by the assumed energy density, which may be quite accurately predicted from the yield criteria. Although the presented result shows the similarity between trabecular bone remodelling and structural optimization, the results has no biological background and has been rigidly proven mechanically and mathematically.

Shape modification using shape derivative
Let us define the function F (z) penalizing the deviation of z from λ and taking into account the insensitivity zone Using this function we replace the heuristic biomimetic strategy by equivalent minimization of the functional After taking shape derivative of this functional (observe that F = dF /dz) we have Here κ is an average curvature at the given point on Γ v and In order to get rid of u we introduce an adjoint equation for p in the weak form where p, ϕ ∈ H 1 Γ 0 (Ω). Then substituting again ϕ := p in the shape derivative of the weak state equation and ϕ := u in the adjoint equation we obtain As a result When the expression in square brackets is negative, then we should add material (V.n > 0), otherwise remove material (V.n < 0). This formula is in complete agreement with intuition: the first term in brackets represents nonlocal influence of boundary modification, the second term describes the change of integral due to spatial variability of F , and the third takes into account the increase or decrease of the area of the surface itself. It is easy to notice, that if energy density lies in the insensitivity zone everywhere, then p ≡ 0 and J λ (Ω) = 0, as expected. The notable difference to the heuristic algorithm is contained in the non-local term: changing the boundary at a certain point changes the energy density everywhere.
The tendency to decrease the area of the variable boundary Γ v may be deleterious, since we only demand, that energy density on it should be constant. Therefore it is advisable to modify the functional by averaging it over the variable surface. The modified functional has the form It's shape derivative is obtained immediately from the same adjoint equation. Let us denote byF the average value of F on Γ v , which again is in perfect agreement with engineering intuition for the specific shape optimization problem. Such derivative may be used for modification of the variable boundary. Namely, if we take where t denotes the step length, then which means that the step was taken in the descent direction.
Implementation remarks Sometimes, due to the method used in parametrization of the domain, it is impossible to move the nodes in the direction perpendicular to the surface Γ v . In such a situation helps the observation that we do not need the steepest descend of the functional. It is enough to ensure only it's decrease. Let us assume for example that the nodes of discretization can move only on certain parallel planes (slices of the domain), and that m 1 , m 2 are basis vectors on these planes, m 1 = 1, m 2 = 1. Furthermore we assume that Γ v is never parallel to these planes, which means that never m 1 .n(x) = 0 and m 2 .n(x) = 0 hold simultaneously. Then we take the vector field in the form V = t 1 m 1 + t 2 m 2 and denote the integrand in (30) as The following possibilities for t 1 and t 2 guarantee the decrease of the functional.
-The first choice: which gives for the step-length t > 0 -The second choice: which gives again negative result When is the construction valid? Let us assume that ∂Ω, Γ 0 , Γ v , Γ 1 consist of smooth patches and data are smooth on Γ 1 . Let us denote also by M the space curve consisting of all the common edges between these patches, and byB(0, r) a closed ball around origin. Then all the constructions described above have meaning if the variable part of the boundary (subset of Γ v ) is contained in Γ v \ (M + B(0, r)) for some fixed r > 0. This requirement is needed in order to give meaning to the derivatives of energy density. In addition, higher Sobolev norms of u over Ω \ (M + B(0, r)) may be bounded by appropriate norms of data over Γ 1 , see Nazarov and Plamenevsky (1994).

Approximation of external normal and curvature
Computation of the normal Let us assume, that ∂Ω is approximated by triangular patches with vertices in points Denote also by S 1 (x) = x/ x the normalization operation.
Next we select certain vertex p 0 and find all other vertices p k , k = 1 . . . K connected to p 0 by edges. In this way we obtain also the triangles T k surrounding p 0 . Each triangle has an internal normal n(T k ).
We approximate the internal normal to the surface ∂Ω by the normalized sum n 0 = S 1 k=K k=1 n(T k ) .
In this way we also get the approximation of the plane tangent to the surface ∂Ω at p 0 : Close to p 0 the surface ∂Ω may be represented by the function of two variables connected with H (p 0 ). We introduce the local coordinate system with origin at p 0 and axes z along n 0 , ξ 1 along b 1 , ξ 2 along b 2 where b 1 is any unit vector perpendicular to n 0 and b 2 = n 0 × b 1 . Then such function has at p k values Computation of curvatures It well known, that if the surface is in the neighbourhood of origin represented by the graph of the function z = f (x 1 , x 2 ) and we denote p = f /1 (0, 0), q = f /2 (0, 0), r = f /11 (0, 0), s = f /12 (0, 0), t = f /22 (0, 0), then the average curvature may be computed from formula If the plane {x 1 , x 2 } is tangent to this graph, that is p = q = 0, then the principal curvatures coincide with eigenval- Under this tangency assumptions, the function z(ξ k 1 , ξ k 2 ) may be approximated by The problem is thus reduced to finding Hessian at the origin. It may be done using least square method. Let us definê Then we have This gives directly the approximation of κ: Testing the approximation In order to test the robustness of the above approximation, we generate two 2 × 2 matrices A 1 and A 2 with entries distributed according to uniform density over interval [-1,1] and a vector [c 1 , c 2 ] with the same property. Then we use the symmetric matrix On the left the view of the surface approximated on the grid (axis labels -coordinate system x, y, z). On the right the contour plot of the typical surface (axis labels -coordinate system x1, x2). The polygon marks a projection of the triangular surface mesh on the plane, δ = 0.31 over the plane {x 1 , x 2 }. The graph of this function is not tangent to the plane at the origin and it may have positive or negative average curvature. Here K denotes the range of the curvatures used in numerical experiment. The next step is to introduce a 3-dimensional grid with spacing δ in all directions. Then we create a typical neighbourhood of triangles around origin with x 1 , x 2 coordinates of vertices on the sides of the square [−1, 1] × [−1, 1]. Finally, the x i values are snapped to the grid, resulting in x g i = δ · round(x i /δ). The third coordinate of vertices is similarly snapped to the vertical grid, z g (x g ) = δ · round(z(x g )/δ). The size of these triangles is in typical relation to the variability of z, as encountered in real examples of surface approximation. In Fig. 5 we see two types of surfaces and the triangles. During the experiment real values of curvature were obtained using exact values and formula (33). In the numerical computation of curvature there where two sources of errors. The values of z g (x g ) are only approximate, due to the snapping to the grid, and their accuracy depends on the grid size δ. The second source of error comes from the approximation of n(0), so the rotated coordinates b 1 , b 2 are only approximately tangent to the surface. Thus the formula (34) is only approximately true.
We have performed two experiments with different grid sizes, corresponding to different ranges of grid points used for identifying the curvature. In Fig. 6 we see scatter plots of real κ against its values computed according to the procedure outlined above. In order to depict the situation more intuitively, we may imagine the surface supported on the pile of cubes with edges of length δ. It turns out that the method is quite robust and insensitive to the parameter δ, if it is smaller than 1/3. The vector n(0) is also approximated quite accurately.

Numerical example
In this section, numerical examples of structural optimization using both the formulation with and without surface curvature measuring term are shown. The typical topology optimization problem -the cantilever beam in bending with the one edge fixed and bending force applied to the middle of the opposite edge was analysed. Material parameters were assumed as follows: Young's modulus 2 · 10 11 Pa, Poisson's ratio 0.3.
To compare optimization results with, and without curvature measuring term results for three different sizes of the lazy zone was computed, and are depicted: -wide size of the insensitivity zone 40-100 MPa in Fig. 7, -middle size of the insensitivity zone 45-100 MPa in Fig. 8, -narrow size of the insensitivity zone 50-100 MPa in Fig. 9.
The curvature is incorporated into the procedure in the following way. The surface value of the strain energy density is treated in the algorithm as a mean value for the thin layer (one layer of the solid finite elements) forming the structural surface. The computed averaged strain energy density for surface elements in this area are identified with the first two terms in brackets in the formula (30), including the impact of both of them. The last term is identified with the influence of the surface curvature. The improvement of descend direction should come from the partially heuristic modification.
-If κ > 0 and F >F , then after biomimetic modification the boundary is additionally moved inside Ω by 50% of the biomimetic step. -If κ < 0 and F >F , then after biomimetic modification the boundary is additionally moved outside Ω by 50% of the biomimetic step.
Thus the curvature term locally enhances or diminishes the purely biomimetic modification of the structure. The numerical experiments seem to confirm this claim and speed method furnishes better descent direction compared to the previous results already published.
The optimal configurations are presented in Fig. 10. The left column presents the optimization results without the curvature measuring term, when the right column presents the optimization results when taking into account the curvature of the evolving structural surface. The rows presents results for different insensitivity zone size. The pictures presented from the top to the bottom illustrate results for the narrowing of the insensitivity zone. The upper level of the insensitivity zone remains the same, while the lower level approaches the upper one. Results still tend to the Fig. 7 Bending of the cantilever beam -the optimization steps and compliance change for the wide size of the insensitivity zone. Left: optimization results without the curvature measuring term. Right: optimization results with the curvature measuring term Fig. 8 Bending of the cantilever beam -the optimization steps and compliance change for the middle size of the insensitivity zone. Left: optimization results without the curvature measuring term. Right: optimization results with the curvature measuring term Fig. 9 Bending of the cantilever beam -the optimization steps and compliance change for the narrow size of the insensitivity zone. Left: optimization results without the curvature measuring term. Right: optimization results with the curvature measuring term Fig. 10 Bending of the cantilever beam -final optimal configurations for the different insensitivity zone size (the upper level remains the same while the lower level approaches the upper one). Left: optimization results without the curvature measuring term. Right: optimization results with the curvature measuring term well known solution. Including the surface curvature measure term helps the numerical procedure to find better the geometrical features during the structural optimization and decrease the number of necessary modification steps.

Conclusions
In the paper the new formulation of the biomimetic optimization based on the trabecular bone phenomenon was presented. The stiffest design is obtained by adding or removal material on the structural surface in the virtual space. The three terms completed the formula to steer the material adding or removal: non-local influence of boundary modification, the spatial variability of the function penalizing the deviation from the assumed strain energy density on the structural surface and measuring of the increase or decrease of the area of the surface itself, by measuring the surface curvature. The notable difference to the heuristic algorithm is contained in the non-local term. Thus, changing the boundary at a certain point changes the energy density everywhere. The structural evolution is based on the shape gradient approximation by the speed method and is separated from the finite element method computations and modeled with use of 2-dimensional images. Finite element method is used to compute the distribution of the strain energy density only.
The assumed value of the strain energy density on the part of the boundary subject to modification could be related to the material properties. Change in the assumed value of the strain energy density results in change of the structural form -topology and volume. In this way, the final structural volume results from the optimization procedure. Instead of imposing volume constraint we parameterize shapes by the assumed energy density, which may be quite accurately predicted from the yield criteria. The remodelling phenomenon with the lazy zone and strain energy density equalization on the trabecular bone surface assumptions could be described from the stiffness point of view. We proved with use of shape derivative (21), that the maximization of a structure stiffness needs the structural form, having on the part of the boundary subject to modification constant value of the strain energy density.
The functional configurations during the process of optimization allows to include size, shape and topology optimization in one numerical procedure. Presented results are the justification that it is possible to use a different approach to the problem of structural optimization. Instead of the assumption of a constant volume, the assumption of a constant strain energy density on the surface of the structure is used and the volume or structural mass results from the optimization procedure. Therefore, presented in the paper approach can be used as a method of structural optimization unrelated already to trabecular bone remodelling phenomenon. This does not change the fact, that studies of the phenomenon of bone remodelling may contribute to achieving better results also in the field of structural optimization. A "lazy zone" concept is here a good example. This interdisciplinary approach can be very fruitful for all disciplines, enabling the flow of concepts from biology to engineering and vice versa.