# 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
During the creation of the class through its constructor, i.e., the class method called initialize is called, and its results is then stored inside the class as follows:
```
10 functionobj=hexahedron(a,b,c)
11%hexahedronConstructaninstanceofthisclass
12%computeweightsandBforlateruse
13[weights,B]=obj.initialize(a,b,c);
14 obj.weights=weights;
15 obj.Bv=B;
16 end
```
where variables _weights_ and \(B\) are the weight and strain-displacement matrix of each Gauss point as computed in _Lines_ 102 to 139 of the hexahedron function of the original homo3D code. Thus, class method initialize contains the implementation of these _Lines_. Two instances of the hexahedron class are utilized in this formulation, one for the topology optimization part and one for the homogenization one. The creation of these instances is performed in the initialization stage and they are replacing _Line_ 24 of the UCOpt3D function as follows:
```
12%KE=lk_H8(nu);
13h1=hexahedron(1x/nelx/2,ly/nely/2,lz/nelz/2);
14h2=hexahedron(0.000001*lx/nlx/2,0.000001*ly/nly/2,0.000001*lz/2);
```
By utilizing the call variables _weights_ and \(B\) of the class hexahedron, we can compute the element stiffness matrix _ke_ as well as the force vector required by the homogenization procedure into two different class methods. _The first class method requires the elasticity tensor for the computation of the stiffness matrix as follows:
```
18 function[ke,fe]=elemStiffness(obj,C)
19 ke=zeros(24,24);
20 fe=zeros(24,6);
21 fori=1:size(obj.weights,1)
22 ke=ke+obj.weights(i)*obj.Bv{i}'*C*obj.Bv{i};
23 fe=fe+obj.weights(i)*obj.Bv{i}'*C;
24 end
```
The procedure performed in the elemStiffness class method is the same as the one implemented by the hexahedronElem function. To integrate the implementation of the class method into the homo3D function, two different _Lines_ of the homo3D function have to be modified. The first one is denoted in _Line_ 84 of the homo3D function where the element stiffness matrix is computed. The specific _Line_ is modified as follows:
```
18%KE=hexahedronElem(a,b,c,CH);
19 KE=h1.elemStiffness(CH);
```
The second one is shown in _Line_ 89 where the derivative of the element stiffness matrix is computed. The specific _Line_ is modified as follows:
```
19%f=@(c)hexahedronElem(0.000001*lx/nlx/2,0.000001*ly/nly/2,0.000001*lz/nlz/2,C);
20%dKe=cellfun(f,DCH,'UniformOutput',false);
31 dKe=cellfun(@h1.elemStiffness,DCH,'UniformOutput',false);Moving to the second class method, this method requires [PERSON]'s modulus \(E\) and [PERSON]'s ratio \(nu\) for the computation of the element stiffness matrix as follows:  
```
2 function[ke,fe]=elemStiffnessC(obj,E,nu)
3C=obj.elasticityTensor(E,nu);
4[ke,fe]=elemStiffness(obj,C);
5 end
```  
This class method utilizes two other class methods to perform the computation. The first is the elasticityTensor class method, which creates the elasticity tensor using [PERSON]'s modulus \(E\) and [PERSON]'s ratio \(nu\), as follows:  
```
8 function[C]=elasticityTensor(E,nu)
9C=E/((1+nu)*(1-2*nu))*[1-nununu000;nu1-nununu000;...
10 nunu1-nu000;000(1-2*nu)/20;...
110000(1-2*nu)/2];
110 end
```  
The second class method refers to the elemStiffness described above. In order to integrate the implementation of class hexahedron into the homo3D script, an extra input argument has to be added to the homo3D function, as follows:  
```
1 function[CH,DCH]=homo3D(1x,ly,lz,E,dE,nu,Emin,hk)
```  
where \(hk\) is the instance of the hexahedron class. To compute the element stiffness matrix using the \(hk\) class instance inside the homo3D function, \(Line\) 12 has to be replaced, as follows:  
```
12%[ke1,fe1]=hexahedronY(dx/2,dy/2,dz/2,nu);
13[ke,fe]=hk.elemStiffnessC(1,nu);
```  
The sensitivity and density filtering techniques are available in the UCOpt3D function in order to deal with the checkerboard effect; both techniques are implemented in accordance to the original top3d and top88 MATLAB codes. Similar to the UCOpt function presented in Ref. [17] the variable \(ft\) is used to denote which filtering technique is to be applied.