%~Mouliné par MaN_auto v.0.37.1 (6604bb9b) 2026-02-02 09:41:00
\documentclass[CRMECA,Unicode,biblatex,published]{cedram}

\addbibresource{CRMECA_Placidi_20260010.bib}

\usepackage{amssymb}
\usepackage{siunitx}
\usepackage{booktabs}

\newcommand{\dif}{\mathop{}\!{\operatorfont{d}}}

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}
\let\oldhat\hat
\renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}}
\let\oldforall\forall
\renewcommand*{\forall}{\mathrel{\oldforall}}

%%%---------------------------------------------------------------------------

%%% ARROWS

%%% Implications
\newcommand*{\myimplies}{\ensuremath{\mathchoice{\quad \Longrightarrow \quad}{\ \Rightarrow \ }{\Rightarrow}{\Rightarrow}}}

%%%---------------------------------------------------------------------------

%%% BAR & OVERLINE

%%% Better widebar macro
\makeatletter
\let\save@mathaccent\mathaccent
\newcommand*\if@single[3]{%
  \setbox0\hbox{${\mathaccent"0362{#1}}^H$}%
  \setbox2\hbox{${\mathaccent"0362{\kern0pt#1}}^H$}%
  \ifdim\ht0=\ht2 #3\else #2\fi
  }
%The bar will be moved to the right by a half of \macc@kerna, which is computed by amsmath:
\newcommand*\rel@kern[1]{\kern#1\dimexpr\macc@kerna}
%If there's a superscript following the bar, then no negative kern may follow the bar;
%an additional {} makes sure that the superscript is high enough in this case:
\newcommand*\widebar[1]{\@ifnextchar^{{\wide@bar{#1}{0}}}{\wide@bar{#1}{1}}}
%Use a separate algorithm for single symbols:
\newcommand*\wide@bar[2]{\if@single{#1}{\wide@bar@{#1}{#2}{1}}{\wide@bar@{#1}{#2}{2}}}
\newcommand*\wide@bar@[3]{%
  \begingroup
  \def\mathaccent##1##2{%
%Enable nesting of accents:
    \let\mathaccent\save@mathaccent
%If there's more than a single symbol, use the first character instead (see below):
    \if#32 \let\macc@nucleus\first@char \fi
%Determine the italic correction:
    \setbox\z@\hbox{$\macc@style{\macc@nucleus}_{}$}%
    \setbox\tw@\hbox{$\macc@style{\macc@nucleus}{}_{}$}%
    \dimen@\wd\tw@
    \advance\dimen@-\wd\z@
%Now \dimen@ is the italic correction of the symbol.
    \divide\dimen@ 3
    \@tempdima\wd\tw@
    \advance\@tempdima-\scriptspace
%Now \@tempdima is the width of the symbol.
    \divide\@tempdima 10
    \advance\dimen@-\@tempdima
%Now \dimen@ = (italic correction / 3) - (Breite / 10)
    \ifdim\dimen@>\z@ \dimen@0pt\fi
%The bar will be shortened in the case \dimen@<0 !
    \rel@kern{0.6}\kern-\dimen@
    \if#31
      \overline{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@}%
      \advance\dimen@0.4\dimexpr\macc@kerna
%Place the combined final kern (-\dimen@) if it is >0 or if a superscript follows:
      \let\final@kern#2%
      \ifdim\dimen@<\z@ \let\final@kern1\fi
      \if\final@kern1 \kern-\dimen@\fi
    \else
      \overline{\rel@kern{-0.6}\kern\dimen@#1}%
    \fi
  }%
  \macc@depth\@ne
  \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar
  \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax}%
  \macc@set@skewchar\relax
  \let\mathaccentV\macc@nested@a
%The following initialises \macc@kerna and calls \mathaccent:
  \if#31
    \macc@nested@a\relax111{#1}%
  \else
%If the argument consists of more than one symbol, and if the first token is
%a letter, use that letter for the computations:
    \def\gobble@till@marker##1\endmarker{}%
    \futurelet\first@char\gobble@till@marker#1\endmarker
    \ifcat\noexpand\first@char A\else
      \def\first@char{}%
    \fi
    \macc@nested@a\relax111{\first@char}%
  \fi
  \endgroup
}
\makeatother

\let\oldbar\bar
\renewcommand*{\bar}[1]{{\mathchoice{\widebar{#1}}{\widebar{#1}}{\widebar{#1}}{\oldbar{#1}}}}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\braces}{\{}{\}}
\DeclarePairedDelimiter{\bracks}{[}{]}

%%%---------------------------------------------------------------------------

\title{Novel material anisotropy optimization based upon granular micromechanics}
\alttitle{Optimisation de l'anisotropie des nouveaux matériaux basée sur la micromécanique granulaire}

\author{\firstname{Luca} \lastname{Placidi}\IsCorresp}
\email{luca.placidi@unipegaso.it}
\address{Department of Engineering, Telematic University Pegaso, Centro Direzionale ISOLA F2, Napoli, 80143, Italy}
\address{International Research Center M\&MoCS, Universit\`a degli Studi dell'Aquila, L'Aquila, Italy}
 
\author{\firstname{Anil} \lastname{Misra}}
\address[2]{International Research Center M\&MoCS, Universit\`a degli Studi dell'Aquila, L'Aquila, Italy}
\address{Department of Civil and Environmental Engineering, Florida International University, Miami, FL 33174-1630, USA}

\author{\firstname{Emilio} \lastname{Barchiesi}}
\address[2]{International Research Center M\&MoCS, Universit\`a degli Studi dell'Aquila, L'Aquila, Italy}
\address{Department of Architecture, design, and urban planning, Universit\`a degli Studi di Sassari, Alghero, Italy}

\author{\firstname{Raimondo} \lastname{Luciano}}
\address{Department of Engineering, Universit degli Studi di Napoli ``Parthenope'', Napoli, 80133, Italy}

\author{\firstname{Francesco} \lastname{Fabbrocino}}
\address[1]{Department of Engineering, Telematic University Pegaso, Centro Direzionale ISOLA F2, Napoli, 80143, Italy}

\thanks{The authors gratefully acknowledge the financial support of the Italian Ministry of University and Research (MUR), Research Grant PRIN 2022 PNRR No.~P202278LFC\_01 on ``DRASTIC -- Digitalized smaRt And SusTaInable Concrete''. Besides, the research reported in the present contribution was also carried out as part of the project ``Metamaterials design and synthesis with applications to infrastructure engineering'' funded by the MUR Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2022, grant~20228CPHN5}
\CDRGrant[Italian Ministry of University and Research (MUR)]{PRIN 2022 PNRR No.~P202278LFC\_01}
\CDRGrant[Italian Ministry of University and Research (MUR)]{20228CPHN5}

\editornote{Article submitted by invitation}
\alteditornote{Article soumis sur invitation}

\begin{abstract}
We propose an evolutionary algorithm that seeks to determine the optimal anisotropy of a deforming body in response to a given applied mechanical load, under the constraint of assigned mass. The algorithm is based for the first time upon a granular micromechanics approach to determine the effective material behavior, making use of an orientation-dependent distribution of normal and tangential elastic grain-grain interactions, whose associated stiffnesses are assumed to depend on an orientation-dependent angular mass density. This novel idea is intrinsically simple and takes advantage of both those penalization techniques, that are generally used in topological optimization, and on those basic concepts of continuum granular micromechanics that are particularly prone to be used in this field. The algorithm is initialized with an isotropic distribution of mass such that the total mass exceeds the desired one. Grain-grain interaction stiffnesses along orientations that are only lightly stressed by the applied load are penalized and the angular mass density is accordingly reduced along these orientations. This yields a non-uniform orientation-dependent mass density and, in turn, an anisotropic constitutive law. The proposed algorithm is numerically evaluated for two load cases implying homogeneous deformations and the response to loading produced by the optimal effective fourth-rank elasticity tensor is compared with that obtained by isotropic angular mass density reduction. The proposed algorithm can be employed for engineering microstructures and as a building block for topology optimization algorithms.
\end{abstract}

\begin{altabstract}
Nous proposons un algorithme évolutionnaire visant à déterminer l'anisotropie optimale d'un corps déformable sous l'effet d'une charge mécanique appliquée, sous la contrainte d'une masse imposée. Cet algorithme repose, pour la première fois, sur une approche de micromécanique granulaire pour déterminer le comportement effectif du matériau. Il exploite une distribution, dépendante de l'orientation, des interactions élastiques normales et tangentielles entre les grains, dont les rigidités associées sont supposées dépendre d'une densité de masse angulaire également dépendante de l'orientation. Cette idée novatrice est intrinsèquement simple et tire parti à la fois des techniques de pénalisation généralement utilisées en optimisation topologique et des concepts fondamentaux de la micromécanique granulaire continue, particulièrement pertinents dans ce domaine. L'algorithme est initialisé avec une distribution isotrope de la masse telle que la masse totale excède celle souhaitée. Les rigidités des interactions entre les grains selon les orientations faiblement sollicitées par la charge appliquée sont pénalisées, et la densité de masse angulaire est réduite en conséquence le long de ces orientations. Il en résulte une densité de masse non uniforme dépendant de l'orientation et, par conséquent, une loi constitutive anisotrope. L'algorithme proposé est évalué numériquement pour deux cas de charge impliquant des déformations homogènes, et la réponse à la charge obtenue avec le tenseur d'élasticité effectif optimal de rang quatre est comparée à celle obtenue par réduction isotrope de la densité de masse angulaire. Cet algorithme peut être utilisé pour l'ingénierie des microstructures et comme élément de base pour les algorithmes d'optimisation topologique.
\end{altabstract}

\keywords{\kwd{Granular micromechanics} \kwd{anisotropy optimization} \kwd{orientation-dependent stiffness} \kwd{metamaterials} \kwd{multiscale materials design}}

\altkeywords{\kwd{Micromécanique granulaire} \kwd{optimisation de l’anisotropie} \kwd{rigidité dépendante de l’orientation} \kwd{métamatériaux} \kwd{conception de matériaux multi-échelles}}

\dateposted{2026-03-17}
\begin{document}
%\input{CR-pagedemetas}
%\end{document}
\maketitle

\section{Introduction}

For many structural mechanics applications it is meaningful to utilize materials with specific anisotropy. Indeed, many natural structural systems, such as bone tissue and wood, develop during their growth a strong material anisotropy that allows them to resist applied loads optimally~\cite{Collins_et_al_2004,Mattheck_Tesari_2004,Wainwright_2020,Zioupos_2004}. Topology optimization approaches in the context of linear elasticity have been widely investigated in the literature~\cite{Banichuk_1981,li_et_al_2020,Mirzendehdel_et_al_2018}. Nevertheless, research efforts in material property or so-called multi-scale optimization are relatively rare, see the recent publication~\cite{key-4} and the seminal papers~\cite{key-1,key-2,key-3}, which pioneered the use of the homogenization method in topology and shape optimization, and contributed to the development of level-set methods and shape derivative techniques, enabling smooth handling of topological changes such as the creation or merging of holes and connections during optimization. Indeed, topology optimization typically aims at finding the best geometry at a single length scale that, for a given mass, is able to guarantee the optimal elastic behavior for given loading/boundary conditions. Relevant applications include, as an instance, composite cable-stayed bridges design~\cite{fabbrocino_et_al_2017}. In this context, the problem is complexified by the fact that the resulting mass distribution may result in the appearance of rotary and micro-inertia, see the review paper~\cite{Eremeyev_et_al_2024}, which require a refinement of the structural analyses for the accurate determination of the functional response within the optimization process. Similarly, refinement of structural analyses within topological optimization processes may be needed in cases in which, explicitly or tacitly, a relationship between Young's modulus and mass density is assumed. This is particularly relevant to the cases of optimized porous materials, such as foams and bones~\cite{Giorgio_et_al_2017_EML,Giorgio_et_al_2019}, for functionally graded plates~\cite{Altenbach_eremeyev_2008,Altenbach_eremeyev_2009}, and when surface phenomena are significant~\cite{Eremeyev_2016}. Refined modeling is also critical when dealing with bone remodeling~\cite{allena_et_al_2005,allena_et_al_2024,scerrato_et_al_2022}. It is worth noting that optimization algorithms studied in the literature are typically of evolutionary type and are bio-inspired, often mimicking bone tissue growth or other biological growth mechanisms~\cite{Giorgio_et_al_2017_MMS}. Among these approaches we find the classical Solid Isotropic Material with Penalization (SIMP) or Evolutionary Structural Optimization (ESO)-type approaches. We also note that optimization methods are applied not only in purely mechanical context~\cite{alessandroni_et_al_2001,rosi_et_al_2010}. In addition to topology, the anisotropy of the material's mechanical response can be designed in an optimal way, too~\cite{desmorat_B_2013,Ranaivomiarana_et_al_2019,Yan_et_al_2019}. This kind of optimization is particularly suitable for the design of optimal microstructures in metamaterials~\cite{Auffray_et_al_2023,Catapano_et_al_2015,desmorat_B_et_al_2020} and composite materials~\cite{barretta_et_al_2015,luciano_willis_2003,luciano_willis_2005}, notwithstanding the need to properly take into account size effects~\cite{Ciallella_et_al_2023,Yang_H_et_al_2023}. Metamaterials~\cite{dell_Isola_et_al_2024,dell_Isola_et_al_2022} are indeed typically designed with an engineering purpose~\cite{giorgio_Et_al_2024}, to optimize one or more desired features~\cite{pranno_et_al_2025,Stilz_et_al_2022}. The identification of a continuum model for metamaterials~\cite{Giorgio_2016,Shekarchizadeh_Et_al_2021} is generally convenient from a computational point of view. However, this poses some peculiar challenges~\cite{dell_Isola_et_al_2023,Eremeyev_et_al_2020,Solyaev_et_al_2022} and optimization concepts need to be adapted to alternative computational strategies~\cite{Reiher_Et_al_2017,Saitta_et_al_2022_eabe}. In the context of fibered metamaterials the optimization of fiber reinforcement is also relevant~\cite{fraternali_et_al_2015}. For granular materials external loading and the features of grain-grain mechanical interactions lead to the emergence and evolution of anisotropy~\cite{Bryant_et_al_2022,Jia_et_al_2017,Wang_et_al_2017}.

Microstructure optimization can be pursued through homogenization procedures~\cite{Barchiesi_et_al_2023,Rahali_et_al_2015} in both static and dynamic cases. Wave propagation features~\cite{Rosi_et_al_2013} are notable dynamic features and can be analyzed to create optimal waveguides. The micro-macro identification performed within the granular micromechanics approach~\cite{Misra_et_al_2022} and the application of the principle of virtual work provide the framework for devising an optimization algorithm that can lead to the design of metamaterials with pre-determined anisotropy~\cite{dellisola_misra_2023}. Higher order gradient~\cite{Eugster_et_al_2022} or nonlocal~\cite{Barretta_et_al_2020,barretta_et_al_2024} effects may also be considered, as these effects are significant whenever long range interactions are relevant. While such higher-order effects are not addressed in this paper, the method that we introduce can be generalized for these purposes; similarly, while this paper addresses only the static case, the proposed method can also target dynamic features such as eigenmodes.

In this paper we integrate the granular micromechanics approach (GMA) into an optimization algorithm aiming at computing the fourth-rank elasticity tensor needed to support optimally a given load with fixed mass. In GMA, a continuous orientation-dependent distribution of elastic grain-grain interactions is utilized to determine material behavior, whose anisotropy can be modulated by varying the stiffness of elastic interactions in the various orientations. The proposed optimization algorithm exploits this aspect of GMA to weight the angular density of elastic energy with an orientation-dependent existence index. The algorithm applies standard techniques used in topological optimization field and reduces the angular mass density for those orientations that do not contribute to the elastic strain energy significantly. Thus, the proposed method essentially adapts the classical principle of penalizing underutilized materials, common in Evolutionary Structural Optimization (ESO), to an orientation-dependent setting within GMA. Such an adaptation is particularly advantageous because it drastically simplifies existing anisotropy optimization frameworks achieving similar results.

In Section~\ref{sec:A-recap-of} we give a recap of granular micromechanics for the sake of self consistency. In this section we also include the necessary relation between mass and stiffnesses. In Section~\ref{sec:An-isotropic-standard} we lay the foundation of an isotropic optimization algorithm based on an evolutionary scheme, which is preliminary to the presentation of the anisotropic optimization algorithm presented in Section~\ref{sec:The-anisotropic-optimization}. In Section~\ref{sec:Numerical-example} we give two simple numerical examples to prove the effectiveness of the anisotropic algorithm. For the sake of simplicity we address only homogeneous deformations, since they do not require to solve a boundary value problem. The presented cases include uniaxial and simple shear loading. The results are given in terms of the newly introduced orientation-dependent mass density. Further, the reaction forces produced by the evolved optimal anisotropic stiffness tensor are compared with respect to those produced by isotropic mass reduction to show, in the same way as established anisotropy optimization methods (e.g., SIMP with orientation variables, homogenization-based microstructure optimization) do, the better performance of the former in terms of lower stiffness reduction for the same mass reduction. In the two case studies (uniaxial strain and simple shear), the optimal anisotropy directions are intuitive, meaning not only that the algorithm confirms physical intuition but it demonstrates capability for realistic, heterogeneous stress states.

\section{Recap of the granular micromechanics approach}\label{sec:A-recap-of}

\subsection{General case}

In the GMA, the elastic internal energy per unit area (henceforth called \emph{elastic internal energy density}) $U$ in the 2D case is given as,
\begin{equation}\label{eq:density elastic energy-1}
U=\int_{\mathcal{S}^1}\mathcal{U} \dif S,
\end{equation}
where the integral in~\eqref{eq:density elastic energy-1} is conceived over all the orientations $\hat{c}$ belonging to the unit circle $S^1$ and $\dif S$ is the angular element over $S^1$. The elastic energy density per orientation, $\mathcal{U}$, is expressed in the following quadratic form:
\begin{equation}\label{eq:density density}
\mathcal{U}=\frac{1}{2}k_{\eta}u_{\eta}^{2}+\frac{1}{2}k_{\tau}u_{\tau}^{2},
\end{equation}
where $k_{\eta}$ is the normal stiffness related to the orientation characterized by the unit vector $\hat{c}$,
\[
\hat{c}=\cos\theta\hat{e}_{1}+\sin\theta\hat{e}_{2},
\]
or by the corresponding angle $\theta$ with respect to the horizontal orientation characterized by the unit vector $\hat{e}_{1}$ ($\hat{e}_{2}$ is the vertical unit vector), $k_{\tau}$ is the tangent stiffness related to the same orientation $\hat{c}$. Further, $u_{\eta}$ is the normal displacement, which is a scalar, given as,
\begin{equation}\label{eq:ueta1}
u_{\eta}=\frac{1}{2}u^{np}\cdot\hat{c}=\frac{1}{2}u_{i}^{np}\hat{c}_{i},
\end{equation}
and $u_{\tau}$ is the tangent displacement, which is a vector,
\begin{equation}\label{eq:utau1}
u_{\tau}=u^{np}-\parens{u^{np}\cdot\hat{c}}\hat{c}=u^{np}-2u_{\eta}\hat{c}.
\end{equation}
In index notation, the tangent displacement reads as
\[
u_{\tau i}=u_{i}^{np}-u_{j}^{np}u_{j}^{np}\hat{c}_{i},
\]
where a sum over repeated indices is considered. In~\eqref{eq:ueta1}--\eqref{eq:utau1} we use the following definition of objective relative displacement:
\begin{equation}\label{eq:objreldispl}
u^{np}=2G\hat{c}L,
\qquad u_{i}^{np}=2G_{ij}\hat{c}_{j}L,
\end{equation}
where $G$ is the Green--Saint-Venant strain tensor,
\begin{equation}\label{eq:green saint venant}
G=\frac{1}{2}\parens{F^{T}F-I},
\qquad G_{ij}=\frac{1}{2}\parens{F_{ki}F_{kj}-\delta_{ij}},
\end{equation}
and $L$ is an average distance between neighbouring grain-pairs. In the last equation, the symbol~$F$ denotes the deformation gradient,
\begin{equation}\label{eq:deformation gradiet}
F=\nabla\chi,\qquad F_{ij}=\chi_{i,j},
\end{equation}
where $\chi$ is the placement function, $\nabla$ is the gradient operator, and subscript commas stand for differentiation with respect to the components of the position $X$ (in the reference configuration) indicated by the following index.

The normal displacement $u_{\eta}$ is therefore given by insertion of~\eqref{eq:objreldispl} into~\eqref{eq:ueta1},
\begin{equation}\label{eq:normal displacement}
u_{\eta}=G_{ij}\hat{c}_{i}\hat{c}_{j}L,
\end{equation}
and the squared tangent displacement by insertion of~\eqref{eq:objreldispl} into~\eqref{eq:utau1},
\begin{equation}\label{eq: squared tangent displacement final form}
u_{\tau}^{2}=L^{2}G_{ab}G_{ij}\parens*{\delta_{ia}\hat{c}_{b}\hat{c}_{j}+\delta_{ja}\hat{c}_{b}\hat{c}_{i}+\delta_{ib}\hat{c}_{a}\hat{c}_{j}+\delta_{jb}\hat{c}_{a}\hat{c}_{i}-4\hat{c}_{i}\hat{c}_{j}\hat{c}_{a}\hat{c}_{b}}.
\end{equation}
Insertion of~\eqref{eq:normal displacement} and~\eqref{eq: squared tangent displacement final form} into the elastic internal energy density~\eqref{eq:density elastic energy-1} yields
\begin{multline*}
U=\int_{\mathcal{S}^1}\braces[\Bigg]{\frac{1}{2}k_{\eta}L^{2}G_{ij}G_{ab}\hat{c}_{i}\hat{c}_{j}\hat{c}_{a}\hat{c}_{b}
\\	\frac{1}{2}k_{\tau}L^{2}G_{ab}G_{ij}\parens*{\delta_{ia}\hat{c}_{b}\hat{c}_{j}+\delta_{ja}\hat{c}_{b}\hat{c}_{i}+\delta_{ib}\hat{c}_{a}\hat{c}_{j}+\delta_{jb}\hat{c}_{a}\hat{c}_{i}-4\hat{c}_{i}\hat{c}_{j}\hat{c}_{a}\hat{c}_{b}}} \dif S,
\end{multline*}
or, taking into account that strain $G$ and the average distance between neighbouring grain-pairs~$L$ do not depend upon the orientation $\hat{c}$,
\[
U=\frac{1}{2}G_{ab}G_{ij}L^{2}\int_{\mathcal{S}^1}\braces[\big]{\parens{k_{\eta}-4k_{\tau}}\hat{c}_{i}\hat{c}_{j}\hat{c}_{a}\hat{c}_{b}
k_{\tau}\parens*{\delta_{ia}\hat{c}_{b}\hat{c}_{j}+\delta_{ja}\hat{c}_{b}\hat{c}_{i}+\delta_{ib}\hat{c}_{a}\hat{c}_{j}+\delta_{jb}\hat{c}_{a}\hat{c}_{i}}} \dif S.
\]
The expression above can be written more concisely as,
\begin{equation}\label{eq:elastic energy general}
U=\frac{1}{2}\mathbb{C}_{\mathit{abij}}G_{ij}G_{ab},
\end{equation}
where the stiffness tensor $\mathbb{C}$ has been identified,
\begin{equation}\label{eq:stiffness identification C}
\mathbb{C}_{\mathit{abij}}=L^{2}\int_{\mathcal{S}^1}\braces[\big]{\parens{k_{\eta}-4k_{\tau}}\hat{c}_{i}\hat{c}_{j}\hat{c}_{a}\hat{c}_{b}
k_{\tau}\parens*{\delta_{ia}\hat{c}_{b}\hat{c}_{j}+\delta_{ja}\hat{c}_{b}\hat{c}_{i}+\delta_{ib}\hat{c}_{a}\hat{c}_{j}+\delta_{jb}\hat{c}_{a}\hat{c}_{i}}} \dif S.
\end{equation}
A more synthetic way for expressing the identification of the stiffness tensor is the following one:
\begin{equation}\label{eq:stiffness identification C shynthetic}
\mathbb{C}_{\mathit{abij}}=L^{2}\int_{\mathcal{S}^1}\bracks*{\parens{k_{\eta}-4k_{\tau}}\hat{c}_{i}\hat{c}_{j}\hat{c}_{a}\hat{c}_{b}+4k_{\tau}\delta_{(i(a}\hat{c}_{b)}\hat{c}_{j)}} \dif S,
\end{equation}
where subscript parentheses stand for the standard single symmetrization rule,
\[
A_{\parens*{ij}}=\frac{1}{2}\parens{A_{ij}+A_{ji}}.
\]
This notation is extended to nested parentheses,
\[
A_{(a(bc)d)}=\frac{1}{2}\parens*{A_{(abcd)}+A_{(acbd)}}=\frac{1}{4}\parens*{A_{abcd}+A_{acbd}+A_{dbca}+A_{dcba}}.
\]

\subsection{Isotropic identification}

In the general anisotropic case, the normal $k_{\eta}$ and tangential $k_{\tau}$ stiffnesses are functions of the orientation $\hat{c}$. In the isotropic case there is no such functional dependence. In the 2D case, where the integration of~\eqref{eq:stiffness identification C} or~\eqref{eq:stiffness identification C shynthetic} is performed over the unit circle $S^1$, we have,
\begin{equation}\label{eq:iso2D}
k_{\eta}=\frac{\bar{k}_{\eta}}{2\pi},
\qquad k_{\tau}=\frac{\bar{k}_{\tau}}{2\pi},
\end{equation}
where $\bar{k}_{\eta}$ and $\bar{k}_{\tau}$ are the integrated stiffnesses over the set of orientations. An analytical evaluation of~\eqref{eq:stiffness identification C shynthetic} taking into account~\eqref{eq:iso2D} yields the following identification of the non null components of the fourth-rank stiffness tensor, see also~\cite{barchiesietlZAMM2021}:
\begin{align}
	& \mathbb{C}_{1111}=\mathbb{C}_{2222}=\frac{1}{8}L^{2}\parens{3\bar{k}_{\eta}+4\bar{k}_{\tau}}, \label{eq:C1111C2222}
\\	& \mathbb{C}_{1122}=\mathbb{C}_{2211}=\frac{1}{8}L^{2}\parens{\bar{k}_{\eta}-4\bar{k}_{\tau}}, \label{eq:C1122C2211}
\\	& \mathbb{C}_{1212}=\mathbb{C}_{1221}=\mathbb{C}_{2112}=\mathbb{C}_{2121}=\frac{1}{8}L^{2}\parens{\bar{k}_{\eta}+4\bar{k}_{\tau}}. \label{eq:C1212}
\end{align}
The remaining components of the stiffness tensor are null in the isotropic case. Considering the plain strain case, the identification of the 2D Lam parameters is given by
\[
\lambda=\mathbb{C}_{1122}=\frac{1}{8}L^{2}\parens{\bar{k}_{\eta}-4\bar{k}_{\tau}},
\qquad \mu=\mathbb{C}_{1212}=\frac{L^{2}}{8}\parens{\bar{k}_{\eta}+4\bar{k}_{\tau}}=\frac{1}{2}\parens*{\mathbb{C}_{1111}-\mathbb{C}_{1122}}.
\]
The Young's modulus $Y_{2D}$ and Poisson's ratio $\nu_{2D}$ can then be computed as (cf.~\cite[(A.12)]{barchiesietlZAMM2021}),
\begin{equation}\label{eq:Y2dnu2d}
Y_{2D}=4\mu\frac{\lambda+\mu}{\lambda+2\mu}=\frac{L^{2}\bar{k}_{\eta}\parens{\bar{k}_{\eta}+4\bar{k}_{\tau}}}{3\bar{k}_{\eta}+4\bar{k}_{\tau}},
\qquad \nu_{2D}=\frac{\lambda}{\lambda+2\mu}=\frac{\bar{k}_{\eta}-4\bar{k}_{\tau}}{3\bar{k}_{\eta}+4\bar{k}_{\tau}}.
\end{equation}
The inversion of~\eqref{eq:Y2dnu2d} gives the normal and tangential integrated stiffnesses in terms of both the Young's modulus $Y_{2D}$ and the Poisson's ratio $\nu_{2D}$,
\begin{equation}\label{eq:ketaktau2Dinterms of youngnu}
\bar{k}_{\eta}=Y_{2D}\frac{2}{L^{2}\parens*{1-\nu_{2D}}},
\qquad \bar{k}_{\tau}=Y_{2D}\frac{\parens*{3\nu_{2D}-1}}{2L^{2}\parens{\nu_{2D}^{2}-1}}.
\end{equation}
Finally, insertion of~\eqref{eq:ketaktau2Dinterms of youngnu} into~\eqref{eq:C1111C2222}, \eqref{eq:C1122C2211} and~\eqref{eq:C1212} yields
\begin{equation}\label{eq:CYrelations}
\mathbb{C}_{1111}=\frac{Y_{2D}}{\parens{1-\nu_{2D}^{2}}},
\qquad \mathbb{C}_{1122}=\frac{Y_{2D}\nu_{2D}}{\parens{1-\nu_{2D}^{2}}}=\nu_{2D}\mathbb{C}_{1111},
\qquad \mathbb{C}_{1212}=\frac{Y_{2D}}{2\parens*{1+\nu_{2D}}}.
\end{equation}

\subsection{Mass density and stiffness}

An important quantity in the optimization problem statement is the mass $M$ of the domain $\mathcal{B}$, defined as
\[
M=\int_{\mathcal{B}}\rho \dif A,
\]
where $\rho$ is the mass per unit area, also called the \emph{mass density}. In the literature it is generally assumed that the Young's modulus~\cite{Giorgio_et_al_2017_EML,Giorgio_et_al_2017_MMS,Giorgio_et_al_2019} is related to the mass density in the following form:
\begin{equation}\label{eq:youngdensity}
Y_{2D}=A\rho^{n},
\end{equation}
where $A$ and $n$ are two constitutive parameters. While this assumption is meaningful for the isotropic case, for the anisotropic case we need to extend it. Such a generalization is pursued assuming the following relationship linking the mass density $\rho$ and the normal $k_{\eta}$ and tangential~$k_{\tau}$ stiffnesses:
\begin{equation}\label{eq:youngdensity-1}
k_{\eta}=A_{\eta}\rho_{c}^{m},
\qquad k_{\tau}=A_{\tau}\rho_{c}^{m},
\qquad\rho=\int_{\mathcal{S}^1}\rho_{c} \dif S,
\end{equation}
where $A_{\eta}$, $A_{\tau}$, and $m$ are three constitutive parameters and an orientation-dependent mass density $\rho_{c}$ has been introduced, see also~\cite{Faria_2001}. It is worth to note that, in the isotropic case of an orientation-independent mass density, the angular density $\rho_{c}$ is given as
\begin{equation}\label{eq: isotropic identification of roc}
\rho_{c}=\frac{\bar{\rho}_{c}}{2\pi} \myimplies \rho=\bar{\rho}_{c},
\end{equation}
and, if $m$ is set to be equal to $n$, equation~\eqref{eq:youngdensity} is retrieved. Indeed, when integrated over the set of possible orientations, the integrated orientation-dependent mass density $\bar{\rho}_{c}$ is identified with the mass density $\rho$. Insertion of~\eqref{eq:iso2D} and~\eqref{eq: isotropic identification of roc} into~\eqref{eq:youngdensity-1} yields
\begin{equation}\label{eq: isotropic representation1}
k_{\eta}=\frac{\bar{k}_{\eta}}{2\pi}=A_{\eta}\parens*{\frac{\bar{\rho}_{c}}{2\pi}}^{m},
\qquad k_{\tau}=\frac{\bar{k}_{\tau}}{2\pi}=A_{\tau}\parens*{\frac{\bar{\rho}_{c}}{2\pi}}^{m}.
\end{equation}
This yields the following explicit expressions of the integrated stiffnesses:
\begin{equation}\label{eq: isotropic representation2}
\bar{k}_{\eta}=A_{\eta}\parens*{2\pi}^{1-m}\bar{\rho}_{c}^{m},
\qquad\bar{k}_{\tau}=A_{\tau}\parens*{2\pi}^{1-m}\bar{\rho}_{c}^{m},
\end{equation}
allowing to write, by their insertion into~\eqref{eq:Y2dnu2d}, the Young's modulus $Y_{2D}$ and the Poisson's ratio $\nu_{2D}$~as
\begin{equation}\label{eq:young and poisson isotropic deriviation}
Y_{2D}=A_{\eta}\parens*{2\pi}^{1-m}L^{2}\frac{\parens{A_{\eta}+4A_{\tau}}}{3A_{\eta}+4A_{\tau}}\rho^{m},
\qquad\nu_{2D}=\frac{A_{\eta}-4A_{\tau}}{3A_{\eta}+4A_{\tau}},
\end{equation}
where~\eqref{eq: isotropic identification of roc} has been used. From the comparison of~\eqref{eq:youngdensity} and~\eqref{eq:young and poisson isotropic deriviation} we have the identification of the two constitutive parameters in~\eqref{eq:youngdensity}
\[
n=m,
\qquad A=A_{\eta}\parens*{2\pi}^{1-m}L^{2}\frac{\parens{A_{\eta}+4A_{\tau}}}{3A_{\eta}+4A_{\tau}},
\]
in terms of the three constitutive parameters in~\eqref{eq:youngdensity-1}.

\subsection{Stress-strain relationships}

The second Piola stress tensor can be derived from the density of the strain energy defined in~\eqref{eq:elastic energy general},
\begin{equation}\label{eq:stress_general}
S_{ij}=\frac{\partial U}{\partial G_{ij}}=\mathbb{C}_{\mathit{ijab}}G_{ab}.
\end{equation}
As a consequence, in the isotropic case we have the standard representation,
\[
S_{ij}=2\mu G_{ij}+\lambda\delta_{ij}G_{aa}.
\]
In particular, for the uniaxial strain case,
\begin{equation}\label{eq:uniaxialstraindefinition}
G_{11}=\epsilon,
\qquad G_{12}=G_{21}=G_{22}=0,
\end{equation}
we have the following anisotropic stress response:
\begin{equation}\label{eq:anisotropic_stress_response_for_uniaxial_horizontal_strain}
S_{11}=\mathbb{C}_{1111}\epsilon,
\qquad S_{12}=S_{21}=\mathbb{C}_{1112}\epsilon,
\qquad S_{22}=\mathbb{C}_{1122}\epsilon,
\end{equation}
which simplifies in the isotropic case to
\begin{equation}\label{eq:isotropic_stress_response_for_uniaxial_horizontal_strain}
S_{11}=\parens{2\mu+\lambda}\epsilon,
\qquad S_{12}=S_{21}=0,
\qquad S_{22}=\lambda\epsilon.
\end{equation}
In terms of Young's modulus and Poisson's ratio, from~\eqref{eq:CYrelations}, the components of the second Piola stress tensor read as
\[
S_{11}=\frac{Y_{2D}}{\parens{1-\nu_{2D}^{2}}}\epsilon,
\qquad S_{12}=S_{21}=0,
\qquad S_{22}=\frac{Y_{2D}\nu_{2D}}{\parens{1-\nu_{2D}^{2}}}\epsilon,
\]
and, in terms of mass density, from~\eqref{eq:youngdensity},
\begin{equation}\label{eq:usreactions}
S_{11}=\frac{A\rho^{n}}{\parens{1-\nu_{2D}^{2}}}\epsilon,
\qquad S_{12}=S_{21}=0,
\qquad S_{22}=\frac{A\rho^{n}\nu_{2D}}{\parens{1-\nu_{2D}^{2}}}\epsilon.
\end{equation}
Finally, for the simple shear strain case,
\begin{equation}\label{eq:shearstraindefinition}
G_{12}=G_{21}=\epsilon,
\qquad G_{11}=G_{22}=0,
\end{equation}
we have the following anisotropic stress response:
\begin{equation}\label{eq:anisotropic_stress_response_for_shear_strain}
S_{12}=S_{12}=2\mathbb{C}_{1212}\epsilon,
\qquad S_{11}=2\mathbb{C}_{1112}\epsilon,
\qquad S_{22}=2\mathbb{C}_{1222}\epsilon,
\end{equation}
which is simplified in the isotropic case to
\begin{equation}\label{eq:isotropic_stress_response_for_shear_strain}
S_{12}=S_{12}=2\mu\epsilon,
\qquad S_{11}=S_{22}=0.
\end{equation}
In terms of Young's modulus and Poisson's ratio, from~\eqref{eq:CYrelations}, we have
\[
S_{12}=S_{12}=\frac{Y_{2D}}{1+\nu_{2D}}\epsilon,
\qquad S_{11}=S_{22}=0,
\]
and, in terms of mass density, from~\eqref{eq:youngdensity},
\begin{equation}\label{eq:ssreactions}
S_{12}=S_{12}=\frac{A\rho^{n}}{1+\nu_{2D}}\epsilon,
\qquad S_{11}=S_{22}=0.
\end{equation}

\section{Isotropic optimization algorithm}\label{sec:An-isotropic-standard}

In this section we introduce an \emph{isotropic optimization algorithm}. The purpose is to better explain the new anisotropic algorithm introduced in the next section.

We introduce, for every step $h$ of the optimization algorithm, an existence index field $s_{h}$ such that
\begin{equation}\label{eq: density hth step}
\rho_{h}=s_{h}\rho,
\quad s_{h}\in\bracks*{0,1},
\quad \forall X\in\mathcal{B},
\quad h=0,1,2,\dots,N,
\end{equation}
where $\rho_{h}$ is the density of the material of the $h$-th iteration of the algorithm and the condition
\begin{equation}\label{eq: initial index isotropic case}
s_{0}=1,
\quad \forall X\in\mathcal{B},
\end{equation}
guarantees that, at the beginning of the iterative scheme ($h=0$), the material has initial mass density $\rho$ and, therefore, from~\eqref{eq:youngdensity}, the initial Young's modulus $Y_{2D}$. It is worth to note that if, for a given point of the domain $\mathcal{B}$, the existence index field is null, then no material is there and, therefore, the corresponding Young's modulus is null from~\eqref{eq:youngdensity}. The optimization algorithm specifies the rules governing the evolution of the existing index field for each point of the domain.

From~\eqref{eq:elastic energy general}, for the $h$-th step of the optimization algorithm, the elastic energy density, per unit area, reads as,
\begin{equation}\label{eq: interna elastic energy density step h}
U_{h}=\frac{1}{2}\mathbb{C}_{\mathit{abij}}^{h}G_{ij}G_{ab},
\end{equation}
where, because of~\eqref{eq:stiffness identification C shynthetic}, \eqref{eq:iso2D}, \eqref{eq:ketaktau2Dinterms of youngnu}, \eqref{eq:youngdensity} and~\eqref{eq: density hth step}, we have
\begin{equation}\label{eq: stiffness isotropic step h}
\mathbb{C}_{\mathit{abij}}^{h}=\frac{A}{\pi\parens{1-\nu_{2D}^{2}}}\rho^{n}s_{h}^{n}\int_{\mathcal{S}^1}\bracks*{4\nu_{2D}\hat{c}_{i}\hat{c}_{j}\hat{c}_{a}\hat{c}_{b}+\parens*{1-3\nu_{2D}}\delta_{(i(a}\hat{c}_{b)}\hat{c}_{j)}} \dif S.
\end{equation}
For each point, at the $(h+1)$-th step, the isotropic optimization algorithm assigns the following value to the existence index,
\begin{equation}\label{eq:isotropic optimization algotithm}
s_{h+1}=s_{h}-\varepsilon\frac{U_{\max}-U_{h}}{U_{\max}},
\quad \forall X\in\mathcal{B},
\end{equation}
where $\varepsilon$ is a small parameter and $U_{\max}>0$ is a strain energy density threshold to be chosen. Herein, we choose $U_{\max}$ as the maximum value of the strain energy density that is attained at the $h$-th step. The idea at the basis of~\eqref{eq:isotropic optimization algotithm} is that the mass density is decreased for those positions $X$ that have low strain energy per unit area. In other words, the lower the strain energy for that position at the $h$-th step, the lower the mass density (and hence stiffness) prescribed by the algorithm for the same position at the $(h+1)$-th step. Besides, when $U_{h}$ equals $U_{\max}$ at the $h$-step, then the mass density at the $(h+1)$-th step remains unchanged.

At each step of the algorithm the total mass is also calculated,
\[
M_{h}=\int_{\mathcal{B}}\rho_{h} \dif A.
\]
This mass is, by the definition of the reduction of the existence index~\eqref{eq:isotropic optimization algotithm} and from~\eqref{eq: density hth step}, reduced. The same occurs to the stiffness, because of~\eqref{eq:youngdensity-1}. The end of the process is prescribed \emph{a priori}. Indeed, if we intend to reduce the initial mass $M_{0}$ by a pre-assigned reduction coefficient $\alpha$ belonging to the range $\bracks*{0,1}$, the algorithm ends at $h=\bar{h}$ when the following condition is reached:
\begin{equation}\label{eq: desired reduction}
M_{\bar{h}}\leq\alpha M_{0},
\quad \alpha\in\bracks*{0,1}.
\end{equation}
Besides, at each step the following conditions must be guaranteed for each point,
\begin{equation}\label{eq:stop condition}
U_{h}\leq U_{\max},
\qquad s_{h}\geq0.
\end{equation}
If the first condition of~\eqref{eq:stop condition} is not respected, then the algorithm has to be stopped even if the target expressed in~\eqref{eq: desired reduction} has not yet been reached. The second condition of~\eqref{eq:stop condition} can be kept true through the following \emph{ad hoc} assumption,
\[
s_{h}\leq\bar{s} \myimplies s_{h+1}=\bar{s}
\]
where $\bar{s}$ prescribes, because of~\eqref{eq: density hth step}, a minimum admissible mass density and, therefore, because of~\eqref{eq:youngdensity}, a minimum admissible value of the 2D Young's modulus and for a given load, of the deformation energy density, per unit area.

It is worth noting that the optimization algorithm reported in this section can be modified taking into account specific further restrictions. For example, it is possible to specify a region for which the reduction rule~\eqref{eq:isotropic optimization algotithm} is not employed and reduction does not take place. Besides, one could initialize the algorithm from a configuration to be reinforced, where we seek to add mass (and thus increase stiffness) in an optimal way. In that case, it suffices to change the sign in front of $\varepsilon$ in~\eqref{eq:isotropic optimization algotithm}, ending the algorithm once the final mass is reached from below and not from above like in~\eqref{eq: desired reduction}.

\section{Anisotropic optimization algorithm}\label{sec:The-anisotropic-optimization}

In this section we propose a new anisotropic optimization algorithm.

Substitution of the normal and tangential stiffness representations in terms of the orientation-dependent mass density~\eqref{eq:youngdensity-1}\textsubscript{1,2} into the general anisotropic form of the stiffness tensor~\eqref{eq:stiffness identification C shynthetic} yields
\begin{equation}\label{eq:stiffness representation anisotropic alg1}
\mathbb{C}_{\mathit{abij}}=L^{2}\int_{\mathcal{S}^1}\rho_{c}^{m}\bracks*{\parens{A_{\eta}-4A_{\tau}}\hat{c}_{i}\hat{c}_{j}\hat{c}_{a}\hat{c}_{b}+4A_{\tau}\delta_{(i(a}\hat{c}_{b)}\hat{c}_{j)}} \dif S.
\end{equation}
In the isotropic optimization algorithm, we introduced in~\eqref{eq: density hth step} an existence index field $s_{h}$, which was initialized from unity in~\eqref{eq: initial index isotropic case}. The algorithm makes this index evolve in such a way that a domain with non homogeneous mass density is produced. In the anisotropic optimization algorithm we generalize this approach by introducing an orientation-dependent existence index~$s_{h}^{c}$ such that
\begin{equation}\label{eq:definition orientation dependent mass density for opt alg}
\rho_{c,h}=s_{h}^{c}\rho_{c},
\quad s_{h}^{c}\in\bracks*{0,1},
\quad\forall X\in\mathcal{B},
\quad\forall\hat{c}\in\mathcal{S}^1,
\quad h=0,1,2,\dots,N.
\end{equation}
The condition
\begin{equation}\label{eq:initial condition odei}
s_{0}^{c}=1,
\quad \forall X\in\mathcal{B},
\quad \forall\hat{c}\in\mathcal{S}^1
\end{equation}
guarantees that
\[
\rho_{c,0}=\rho_{c}
\]
and hence that, at the beginning of the iterative scheme, when $h=0$, the material possesses an initial mass density, according to~\eqref{eq:youngdensity-1}\textsubscript{3}, equal to
\begin{equation}\label{eq:ro e roc0}
\rho=\int_{\mathcal{S}^1}\rho_{c} \dif S=2\pi\rho_{c,0}.
\end{equation}
The anisotropic optimization algorithm prescribes an evolution of the orientation-dependent existence index $s_{h}^{c}$ for each point of the domain and for each orientation. Therefore, for the $h$-th step, the elastic energy density $U_{h}$ can still be represented with~\eqref{eq: interna elastic energy density step h}, except that the representation of the stiffness matrix in~\eqref{eq: stiffness isotropic step h} is not valid anymore and we have from~\eqref{eq:stiffness representation anisotropic alg1} and~\eqref{eq:definition orientation dependent mass density for opt alg}
\begin{equation}\label{eq:Ch}
\mathbb{C}_{\mathit{abij}}^{h}=L^{2}\int_{\mathcal{S}^1}\parens{s_{h}^{c}\rho_{c}}^{m}\bracks*{\parens{A_{\eta}-4A_{\tau}}\hat{c}_{i}\hat{c}_{j}\hat{c}_{a}\hat{c}_{b}+4A_{\tau}\delta_{(i(a}\hat{c}_{b)}\hat{c}_{j)}} \dif S.
\end{equation}

A natural anisotropic generalization of the optimization algorithm in~\eqref{eq:isotropic optimization algotithm} assigns the new orientation-dependent existence index for each point of the domain and for each orientation as follows:
\begin{equation}\label{eq:anisotropic optimization algotithm}
s_{h+1}^{c}=s_{h}^{c}-\varepsilon\frac{\mathcal{U}_{\max}-\mathcal{U}_{h}}{\mathcal{U}_{\max}},
\quad \forall X\in\mathcal{B},
\quad \forall \hat{c}\in\mathcal{S}^1,
\end{equation}
where $\varepsilon$ is, again, a small parameter. The elastic internal energy per unit area and per unit orientation $\mathcal{U}$ is defined in~\eqref{eq:density density} and, for a given step of the optimization algorithm, is $\mathcal{U}_{h}$. From~\eqref{eq:normal displacement}--\eqref{eq: squared tangent displacement final form} and~\eqref{eq:youngdensity-1}\textsubscript{1,2}, we have
\[
\mathcal{U}_{h}=\frac{1}{2}L^{2}\parens{s_{h}^{c}\rho_{c}}^{m}\bracks*{\parens{A_{\eta}-4A_{\tau}}\hat{c}_{i}\hat{c}_{j}\hat{c}_{a}\hat{c}_{b}+4A_{\tau}\delta_{(i(a}\hat{c}_{b)}\hat{c}_{j)}}G_{ij}G_{ab}.
\]
The quantity $\mathcal{U}_{\max}$ is a strain energy density per unit orientation threshold. The idea on which~\eqref{eq:anisotropic optimization algotithm} is based is the following. The orientation-dependent mass density of the next step is prescribed in such a way that it is decreased for those positions $X$ and orientations $\hat{c}$ that have low strain energy density per unit orientation. In other words, the lower the strain energy per unit orientation for that position and for that orientation at the $h$-step, the lower the orientation-dependent mass density for the same position and the same orientation at the $(h+1)$-th step. Besides, if the orientation-dependent density of the elastic strain energy $\mathcal{U}_{h}$ reaches the threshold value $\mathcal{U}_{\max}$ at the $h$-step, then the orientation-dependent mass density at the $(h+1)$-step remains unchanged with respect to the $h$-th step.

At each step of the algorithm, the total mass is calculated as an integral of the orientation-dependent mass density over the domain $\mathcal{B}$ and over the orientation space $\mathcal{S}^1$,
\[
M_{h}=\int_{\mathcal{B}}\bracks*{\int_{\mathcal{S}^1}\rho_{c,h} \dif S} \dif A.
\]
Thus, at each step, by the definition of the reduction of the orientation-dependent existence index defined by~\eqref{eq:anisotropic optimization algotithm} and from~\eqref{eq:definition orientation dependent mass density for opt alg}, the total mass is reduced along with the stiffness of the body, because of~\eqref{eq:youngdensity-1}. The end of the process is prescribed \emph{a priori}. If we intend to reduce the initial mass $M_{0}$ with a pre-assigned reduction coefficient $\alpha$ belonging to the range $\bracks*{0,1}$, the algorithm ends at $h=\bar{h}$ when the following condition is reached:
\begin{equation}\label{eq: desired reduction-1}
M_{\bar{h}}\leq\alpha M_{0}.
\end{equation}
At each step, the following further conditions must be guaranteed for each point and for each orientation:
\begin{equation}\label{eq:stop condition-1}
\mathcal{U}_{h}\leq\mathcal{U}_{\max},
\quad s_{h}^{c}\geq0,
\quad \forall X\in\mathcal{B},
\quad \forall\hat{c}\in\mathcal{S}^1.
\end{equation}
If the first condition in~\eqref{eq:stop condition-1} is not respected, then the algorithm has to be stopped even if the target expressed in~\eqref{eq: desired reduction} has not yet been reached. The second condition of~\eqref{eq:stop condition-1} can be kept true with the following \emph{ad hoc} assumption,
\begin{equation}\label{eq:threshold orient dep ei}
s_{h}^{c}\leq\bar{s}^c \myimplies s_{h+1}^{c}=\bar{s}^c,
\end{equation}
where $\bar{s}^c$ prescribes with~\eqref{eq: density hth step} an arbitrary minimum admissible value of the orientation-dependent existence index, which guarantees, at each point and for each orientation, not only its positive semi-definiteness but also, because of~\eqref{eq:youngdensity-1}, the positive definiteness of the elastic energy density per unit area.

As for the isotropic algorithm, it is worth noting that the optimization algorithm reported in this section can be modified taking into account specific further restrictions. For example, it is possible to specify a region for which the reduction rule~\eqref{eq:anisotropic optimization algotithm} is not employed and reduction does not take place. Besides, one could initialize the algorithm from a configuration to be reinforced, where we seek to add mass (and thus increase stiffness) in an optimal way. In that case, it suffices to change the sign in front of $\varepsilon$ in~\eqref{eq:anisotropic optimization algotithm}, ending the algorithm once the final mass is reached from below and not from above like in~\eqref{eq: desired reduction-1}. The anisotropic algorithm may be particularly suitable for those applications where adding fibers in a given orientation is easier than adding fibers with orientations distributed randomly. We finally mention that, as for the isotropic algorithm, if the material behaves differently in tension and in compression, which is a quite standard case for cementitious and fibered materials, then the proposed algorithm can be adapted without any substantial changes of the main concepts.

\section{Numerical examples}\label{sec:Numerical-example}

\subsection{General setting}

We take into account only study cases with homogeneously deforming body, so that the optimization process acts only on the anisotropy evolution while keeping unchanged the topology of the domain $\mathcal{B}$. The employed material parameters are listed in Table~\ref{tab:Material-parameters}. In the initial configuration the material is isotropic (with Young's modulus $Y_{2D}$ and Poisson's ratio $\nu_{2D}$) and the body has a unit volume, so that the value of the initial mass $M_{0}$ coincides with that of the initial mass density $\rho_{0}$. Because of~\eqref{eq:ro e roc0}, we have an initial orientation-independent mass density $\rho_{c,0}$. Once fixed the value of the intergranular distance $L$, the initial normal and tangential stiffnesses are deduced through equations~\eqref{eq:iso2D} and~\eqref{eq:ketaktau2Dinterms of youngnu}. Further, owing to~\eqref{eq: isotropic identification of roc} and~\eqref{eq: isotropic representation1}, once the value of the exponent $m$ is known, the values of the constitutive parameters $A_{\eta}$ and $A_{\tau}$ can be determined.

\begin{table}[h!]
\caption{Material parameters of the initial isotropic configuration.}\label{tab:Material-parameters}
\begin{tabular}{@{}cccccc@{}}
\toprule
\boldmath{$Y_{2D}$} {[}N/m{]} & \boldmath{$\nu_{2D}$} {[}1{]} & \boldmath{$M_{0}$} {[}\unit{Kg}{]} & \boldmath{$\rho_{0}$} {[}\unit{Kg/m^2}{]} & \boldmath{$L$} {[}m{]} & \boldmath{$m$} {[}1{]}\\
\midrule
$10^{9}$ & $0.3$ & 1 & 1 & $0.1$ & 3\\
\bottomrule
\end{tabular}
\end{table}

The algorithm further requires the specification of the applied strain and other parameters, which are reported in Table~\ref{tab:Algorithmic-parameters}. The strain $\epsilon$ is a datum needed in the definitions~\eqref{eq:uniaxialstraindefinition} and~\eqref{eq:shearstraindefinition} of the loading cases. This datum has to be integrated with the parameter $\varepsilon$ introduced in~\eqref{eq:anisotropic optimization algotithm}, the minimum admissible value $\bar{s}^c$ of the orientation-dependent existence index introduced in~\eqref{eq:threshold orient dep ei}, and the mass fraction in~\eqref{eq: desired reduction-1} that we want to achieve at the end of the application of the algorithm.

\begin{table}
\caption{Strain value considered for the problems explored in Sections~\ref{subsec:Horizontal-strain} and~\ref{subsec:Homogeneous-simple-shear} and values of further parameters needed to run the presented algorithm.}\label{tab:Algorithmic-parameters}
\begin{tabular}{@{}cccc@{}}
\toprule
\boldmath{$\epsilon$} {[}1{]} & \boldmath{$\varepsilon$} {[}1{]} & \boldmath{$\bar{s}^{c}$} {[}1{]} & \boldmath{$\alpha$} {[}1{]}\\
\midrule
$10^{-5}$ & $10^{-2}$ & $10^{-2}$ & $0.3$\\
\bottomrule
\end{tabular}
\end{table}

The results presented in the next sections will include firstly the orientation-dependent existence index against the orientation. The objective is to identify those orientations which undergo a density reduction and, because of~\eqref{eq:youngdensity-1}, also a stiffness reduction. The reduction of total mass will be then shown against the algorithm step number. For a given final reduction of mass, we then show the corresponding reduction in the stress response. Finally, such a reduction is compared with that obtained through the isotropic algorithm. It will be proved that the isotropic algorithm yields a higher stiffness reduction with respect to the anisotropic one.

\subsection{Uniaxial strain case}\label{subsec:Horizontal-strain}

The first numerical experiment concerns the uniaxial strain problem expressed by~\eqref{eq:uniaxialstraindefinition} and graphically represented on the left-hand side of Figure~\ref{fig:Evolution-of-the orientation dependet mass density}. The use of the anisotropic algorithm yields the evolution of the orientation-dependent existence index $s_{h}^{c}$ shown in Figure~\ref{fig:Evolution-of-the orientation dependet mass density} for some representative steps. It can be seen that at the beginning of the algorithm, when $h=0$, the index is uniform and equal to unity because of condition~\eqref{eq:initial condition odei}. However, since the vertical orientations, with $\theta=\pi/2$ and with $\theta=3\pi/2$, do not contribute to the strain energy, the corresponding orientation-dependent mass density (and hence existence index) are reduced. When, for a given orientation, the existence index reaches the minimum threshold $\bar{s}^c$, then mass density reduction is stopped according to~\eqref{eq:threshold orient dep ei}. Besides, the horizontal orientations, with $\theta=0$, $\theta=\pi$ and $\theta=2\pi$, give the maximum contribution to the strain energy density and, therefore, are not reduced by the algorithm.

\begin{figure}[h!]
\includegraphics[scale=0.4]{fig_1_odmduhs}
\caption{Evolution of the orientation-dependent existence index for selected iteration steps for the uniaxial strain problem schematized on the left-hand side. At the top, the iteration numbers range from the initial isotropic configuration ($h=0$) with $M_{0}=1~\unit{Kg}$ to an intermediate configuration ($h=110$) with $M_{110}\backsimeq0.5~\unit{Kg}$. At the bottom, the iteration numbers range from an intermediate configuration ($h=120$) with $M_{120}\backsimeq0.5~\unit{Kg}$ to the final step ($h=240$) with $M_{240}\backsimeq0.3~\unit{Kg}$. At the bottom, we make explicit not only that the horizontal orientations are so important that the algorithm maintains their associated density at the initial level, but also that the vertical orientations are so unimportant that the optimization algorithm penalizes them so much that their associated density vanishes.}\label{fig:Evolution-of-the orientation dependet mass density}
\end{figure}

The reduction of the mass is also shown in Figure~\ref{fig:Mass-evolution-for_us}, which shows a decreasing trend that stops, according to~\eqref{eq: desired reduction-1}, when the prescribed fraction $\alpha$ is reached. In Figure~\ref{fig:Mass-evolution-for_us} we also plot the line corresponding to $0.5~\unit{Kg}$, which delimits in terms of total mass the plots at the top and those at the bottom of Figure~\ref{fig:Evolution-of-the orientation dependet mass density}. It is interesting to note that, after that value of the mass is reached, mass reduction is slower with respect to the iteration number.

\begin{figure}
\includegraphics[scale=0.5]{fig_2_massevolution_us}
\caption{Mass evolution for the uniaxial strain problem. Two lines corresponding to $M_{112}\backsimeq0.5~\unit{Kg}$ and $M_{269}\backsimeq0.3~\unit{Kg}$ are also shown.}\label{fig:Mass-evolution-for_us}
\end{figure}

Another important index is related to the optimized fourth-order elasticity tensor $\mathbb{C}$. In particular a quantitative measure of the resulting macroscopic anisotropy is provided in Figure~\ref{fig:anisotropyindex} by the ratio $r$,
\begin{equation}\label{eq: r}
r=\frac{\mathbb{C}_{1111}}{\mathbb{C}_{2222}},
\end{equation}
between the horizontal $\mathbb{C}_{1111}$ and the vertical $\mathbb{C}_{2222}$ stiffness. It is worth to note that such a ratio increases from the isotropic value at $r=1$ to a very large value at $r=r_{\max}\approx35$. However, after a certain value, the occurrence of the threshold $\bar{s}^c$, defined in~\eqref{eq:threshold orient dep ei}, induces an opposite effect. The reason is that the deformed orientations continue to reduce their weight and the non-deformed orientations are prescribed to be constant over the algorithm index.

\begin{figure}[h!]
\includegraphics[scale=0.5]{fig_3_macroscopic_anisotropy}
\caption{Evolution of the anisotropic index defined in~\eqref{eq: r}.}\label{fig:anisotropyindex}
\end{figure}

Finally, the evolution of the response $S_{11}$ produced by the evolution of the fourth-rank stiffness tensor, computed using~\eqref{eq:anisotropic_stress_response_for_uniaxial_horizontal_strain}\textsubscript{1}, is given in Figure~\ref{fig:Evolution-of-the orientation dependet mass density-1}, where also the evolution of the same horizontal reaction obtained with the isotropic algorithm~\eqref{eq:usreactions}\textsubscript{1} is plotted, evaluated at each step as follows:
\[
S_{11h}=\mathbb{C}_{1111}^{h}\epsilon.
\]
We recall that the stiffness at the $h$-th step in the equation above is computed by means of~\eqref{eq: stiffness isotropic step h}. From Figure~\ref{fig:Evolution-of-the orientation dependet mass density-1} it is evident that the reduction of the stress reaction is higher for the isotropic algorithm, indicating that the reduction of mass suggested by the anisotropic algorithm is more efficient. In Figure~\ref{fig:Evolution-of-the orientation dependet mass density-1} we also represent the line corresponding to total mass equal to $0.5~\unit{Kg}$. We observe that, on the one hand, the isotropic algorithm attains a close-to-zero stiffness when total mass reaches $0.5~\unit{Kg}$ and that further mass reduction yields a limited stiffness reduction. On the other hand, the stiffness prescribed by the anisotropic algorithm for the mass threshold at $0.5~\unit{Kg}$ is not negligible and further mass reduction yields non-negligible stiffness reduction.

\begin{figure}
\includegraphics[scale=0.5]{fig_4_comparisonstiffness_us}
\caption{Evolution of the horizontal reaction for the uniaxial strain case according to the anisotropic and isotropic algorithms.}\label{fig:Evolution-of-the orientation dependet mass density-1}
\end{figure}

\begin{figure}[h!]
\includegraphics[scale=0.4]{fig_5_odmduss}
\caption{Evolution of the orientation-dependent existence index for selected iteration steps for the simple shear strain problem schematized on the left-hand side. At the top, the iteration numbers range from the initial isotropic configuration ($h=0$) with $M_{0}=1~\unit{Kg}$ to an intermediate configuration ($h=140$) with $M_{140}\backsimeq0.5~\unit{Kg}$. At the bottom, the iteration numbers range from an intermediate configuration ($h=150$) with $M_{150}\backsimeq0.5~\unit{Kg}$ to the final step ($h=330$) with $M_{330}\backsimeq0.3~\unit{Kg}$. At the bottom, we make explicit not only that the oblique orientations are so important that the algorithm maintains their associated density at the initial level.}\label{fig:Evolution-of-the orientationeiss}
\end{figure}

\subsection{Simple shear strain case}\label{subsec:Homogeneous-simple-shear}

The second numerical experiment concerns the simple shear strain problem expressed by~\eqref{eq:shearstraindefinition} and graphically represented on the left-hand side of Figure~\ref{fig:Evolution-of-the orientationeiss}. The use of the anisotropic algorithm yields the evolution of the orientation-dependent existence index $s_{h}^{c}$ shown in Figure~\ref{fig:Evolution-of-the orientationeiss} for some selected steps. Initial steps are accounted for in the plot at the top, while final ones are accounted for at the bottom. We notice that, at the beginning, when $h=0$, the index is uniformly equal to unity by definition~\eqref{eq:initial condition odei}. However, both the vertical orientations, with $\theta=\pi/2$ and $\theta=3\pi/2$, and the horizontal ones, with $\theta=0$, $\theta=\pi$, and $\theta=2\pi$, do not contribute to the strain energy density and, therefore, the associated contributions to the mass density (and existence index field) are reduced. It is seen that, in accordance with~\eqref{eq:threshold orient dep ei}, when the existence index reaches the minimum threshold $\bar{s}^c$, there is no further reduction. Besides, the oblique orientations, with $\theta=\pi/4$, $\theta=3\pi/4$, $\theta=5\pi/4$, and $\theta=7\pi/4$, give the maximum contribution to the strain energy density and, therefore, the associated contributions to the mass density are not reduced by the anisotropic algorithm.

The mass reduction is plotted in Figure~\ref{fig:Mass-evolution-for_us} against the step number.

\begin{figure}
\includegraphics[scale=0.5]{fig_6_massevolution_ss}
\caption{Mass evolution for the simple shear strain problem. Two lines corresponding to $M_{143}\backsimeq0.5~\unit{Kg}$ and $M_{335}\backsimeq0.3~\unit{Kg}$ are also shown.}\label{fig:Mass-evolution-for_us-1}
\end{figure}

Finally, the evolution of the response $S_{12}$ produced by the evolution of the fourth-rank stiffness tensor, computed using~\eqref{eq:anisotropic_stress_response_for_uniaxial_horizontal_strain}\textsubscript{1}, is given in Figure~\ref{fig:Evolution-of-thess}, where also the evolution of the same shear reaction obtained with the isotropic algorithm~\eqref{eq:usreactions}\textsubscript{1} is plotted, evaluated at each step as follows,
\[
S_{12h}=2\mathbb{C}_{1212}^{h}\epsilon.
\]
From Figure~\ref{fig:Evolution-of-thess} it is evident that the reduction of the stress reaction is higher for the isotropic algorithm, indicating that the reduction of mass suggested by the anisotropic algorithm is more efficient.

\begin{figure}[h!]
\includegraphics[scale=0.5]{fig_7_comparisonstiffness_ss}
\caption{Evolution of the horizontal reaction for the simple shear strain case according to the anisotropic and isotropic algorithms.}\label{fig:Evolution-of-thess}
\end{figure}

\section{Conclusions and outlooks}

We have, for the first time, developed an optimization algorithm that exploits the micro-macro identification provided by the Granular Micromechanics Approach (GMA) to optimize material anisotropy under general loading conditions. The presented numerical examples demonstrate that the granular micromechanics framework is ideally suited for anisotropy optimization. The main advantage of GMA lies in its use of continuous distributions of elastic interactions over the space of orientations to characterize material behavior. Accordingly, in our optimization strategy, we introduce an orientation-dependent existence index governed by each direction's contribution to the strain energy density, thereby promoting material removal primarily along directions that contribute minimally to load transfer. As shown by the numerical tests, this leads to a smaller stiffness loss for a prescribed mass reduction when compared to isotropic mass-reduction benchmarks. Practically, the algorithm naturally aligns material anisotropy with the dominant kinematic and stress-resultant paths.

This work deliberately focuses on homogeneous strain states to clearly isolate the effects of the anisotropic update rule. While this focus clarifies the underlying mechanism, it also represents a limitation: realistic boundary-value problems generally involve heterogeneous fields and geometric constraints that can influence the optimal orientation distribution. For this reason, a short-term outlook involves applying the same algorithm to problems combining anisotropic and topological optimization, which require solving boundary-value problems. Moreover, manufacturability constraints (e.g., admissible fiber angles or feasible print paths in additive manufacturing) have not been enforced and should be incorporated in future engineering-oriented applications.

From a methodological standpoint, several aspects merit further investigation, including: (i)~the sensitivity of results to the chosen stiffness-density relationship and to the discretization of the orientation space; (ii)~the use of regularization or filtering techniques to suppress spurious oscillations in the orientation field; (iii)~convergence and stopping criteria specifically tailored to anisotropy-driven updates must still be implemented; and (iv)~robustness with respect to uncertainties in loading and material parameters must still be analyzed. Benchmarking against established topology-optimization schemes (e.g., SIMP- or ESO-type baselines) would also help quantify performance gains in stiffness retention relative to mass savings.

On the application side, the short-term path we envisage, coupling anisotropic and topological optimization, can be pursued by embedding the current orientation-update rule within a standard finite element loop, thereby addressing nonuniform strain fields and general boundary conditions. Mid-term directions include: (i)~extending the formulation to dynamic objectives (e.g., wave steering, band-gap widening, or mode shaping); (ii)~targeting stability-related metrics (e.g., buckling load factors, imperfection sensitivity, or selected eigenmodes); and (iii)~incorporating process-aware constraints such as curvilinear fiber placement and layer-wise printability.

Finally, further outlooks of the present work are experimental or high-fidelity numerical validation on representative metamaterial cells would strengthen the connection between the granular micromechanics rationale and continuum-level effective responses. Parallel implementations and GPU acceleration could make the approach competitive for large-scale problems, while multi-objective formulations (e.g., stiffness-to-weight optimization combined with directional compliance bounds) could improve the flexibility of the anisotropic update paradigm.

\section*{CRediT authorship contribution statement}

\begin{description}
\item[Luca Placidi] conceptualization; methodology; formal analysis; visualization; writing ---~original draft; supervision; project administration.
\item[Anil Misra] conceptualization; methodology; validation; writing ---~review \& editing; supervision.
\item[Emilio Barchiesi] methodology; software; formal analysis; visualization; writing ---~review \& editing.
\item[Raimondo Luciano] validation; resources; writing ---~review \& editing.
\item[Francesco Fabbrocino] supervision; resources; writing ---~review \& editing; project administration.
\end{description}

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\printbibliography

\end{document}