# Topology Optimization Based Material Design for 3D Domains Using MATLAB
## 4 The MATLAB Code Implementation
### The Seven Components of the UCOpt3D Function
#### 4.1.4 Finite Element Analysis Section  
In the finite element analysis section the creation of the element stiffness matrix is added. The element stiffness matrix for a hexahedron element is created using function hexahedronElem, which is a modified version of the hexahedronElem function used in the homo3D script. The input arguments for this function are the element dimensions as well as the elasticity tensor required for creating the elemental stiffness matrix of a hexahedron element. Thus, the input arguments of the function are the following:  
```
functionke=hexahedronElem(a,b,c,C)
```  
where \(a,b,\) and \(c\) are the element dimensions and \(C\) is the elasticity tensor. The homogenized elasticity tensor provided by the homogenization formulation is provided as the input argument for the hexahedronElem function as follows:  
```
kE=hexahedronElem(1x/nelx/2,ly/nely/2,lz/nelz/2,CH);

#### 4.1.5 Objective Function--Sensitivity Analysis, Filtering, and Update Scheme Sections

The computation of the objective function \(c\) and its derivative with respect to the material densities \(dc\) relies on the process used for calculating the derivative of the stiffness matrix \(dKe\) with respect to the derivative of the homogenized elasticity tensor \(DCH\) for each finite element comprising the unit cell. To implement these computations, a MATLAB cell array has to be constructed first, where \(dKe\) will be computed for each unit cell element. Thus, function hexahedronElem has to be applied for every index of the \(DCH\) cell array. The MATLAB cellfun method is utilized to avoid using time consuming loops, as follows:
```
f=@(C)hexahedronElem(0.000001*1x/n1x/2,0.000001*1y/n1y/2,0.000001*1z/n1z/2,C); dKe=cellfun(f,DCH,'UniformOutput',false);
```
After the computation of the cell array \(dKe\), the derivative of the compliance is computed with respect to the macro elements. This computation is similar to that in \(Line\) 78 of the original top3d MATLAB code. It is implemented however for every finite element of the micro design domain using the \(dKe\) cell array. This is also performed using the cellfun MATLAB method as follows:
```
f=@(dKE)(reshape(sum((U(edofMat)*dKE). *U(edofMat),2),[nely, nelx,nelz])); ce=cellfun(f,dKe,'UniformOutput',false);
```
Thus, variable \(ce\) is a cell array whose dimensions coincide with those of the unit cell, and in every cell a tensor containing the product of \(Line\) 91 is assigned. For calculating the objective function \(c\) and its derivative \(dc\), the contents of the cell array \(ce\) are summed as follows:
```
c=sum(sum(sum(cell2 mat(cellfun(@sum,ce,'UniformOutput', false))))); dc=-cell2 mat(cellfun(@sum,cellfun(@sum,cellfun(@sum,ce,' UniformOutput',false),'UniformOutput',false),'UniformOutput',false));
```
While the last two components of the UCOpt3D function, i.e., the filtering technique as well as the update scheme, remain the same as implemented in the original top3d code, no information is therefore provided herein. Aiming to avoid the computation of the hexahedron element stiffness matrix from scratch on every iteration of the optimization process, as required both in the macro and micro domains, a new MATLAB class structure was developed. The class is called hexahedron and its purpose is to compute the element stiffness matrix of the hexahedron, utilizing variables such as weights and the strain-displacement matrix \(B\), which are computed once and then stored inside the class. Thus, the class constructor is as follows:
```
obj=hexahedron(a,b,c)
```
where \(a,b,\) and \(c\) are the half dimensions of the element in the \(x,y,\) and \(z\) dimensions.