Topology Optimization - What does 2 and 1 mean in F(2,1)? (line 79 in the code)

This is a standard code for the topology optimization of a MBB beam as in the figure. The force has to be defined at the upper left corner. The design domain is descretized into finite elements, say 9 elements in the x-direction and 3 elements in the y-direction. Both node numbers and element numbers are numbered column wise from left to right. I am under the impression that the force is applied on the node. In that case the node number 1. But what does F(2,1) mean?

%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 %%% %%%% CODE MODIFIED FOR INCREASED SPEED, September 2002, BY OLE SIGMUND %%%
function top_mbb(nelx,nely,volfrac,penal,rmin);
% INITIALIZE
x(1:nely,1:nelx) = volfrac; 
loop = 0; 
change = 1.;
% START ITERATION
while change > 0.01  
  loop = loop + 1;
  xold = x;
% FE-ANALYSIS
  [U]=FE(nelx,nely,x,penal);         
% OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
  [KE] = lk;
  c = 0.;
  for ely = 1:nely
    for elx = 1:nelx
      n1 = (nely+1)*(elx-1)+ely; 
      n2 = (nely+1)* elx   +ely;
      Ue = U([2*n1-1;2*n1; 2*n2-1;2*n2; 2*n2+1;2*n2+2; 2*n1+1;2*n1+2],1);
      c = c + x(ely,elx)^penal*Ue'*KE*Ue;
      dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue;
    end
  end
% FILTERING OF SENSITIVITIES
  [dc]   = check(nelx,nely,rmin,x,dc);    
% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD
  [x]    = OC(nelx,nely,x,volfrac,dc); 
% PRINT RESULTS
  change = max(max(abs(x-xold)));
  disp([' It.: ' sprintf('%4i',loop) ' Obj.: ' sprintf('%10.4f',c) ...
       ' Vol.: ' sprintf('%6.3f',sum(sum(x))/(nelx*nely)) ...
        ' ch.: ' sprintf('%6.3f',change )])
% PLOT DENSITIES  
  colormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6);
end 
%%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xnew]=OC(nelx,nely,x,volfrac,dc)  
l1 = 0; l2 = 100000; move = 0.2;
while (l2-l1 > 1e-4)
  lmid = 0.5*(l2+l1);
  xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid)))));
  if sum(sum(xnew)) - volfrac*nelx*nely > 0;
    l1 = lmid;
  else
    l2 = lmid;
  end
end
%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dcn]=check(nelx,nely,rmin,x,dc)
dcn=zeros(nely,nelx);
for i = 1:nelx
  for j = 1:nely
    sum=0.0; 
    for k = max(i-floor(rmin),1):min(i+floor(rmin),nelx)
      for l = max(j-floor(rmin),1):min(j+floor(rmin),nely)
        fac = rmin-sqrt((i-k)^2+(j-l)^2);
        sum = sum+max(0,fac);
        dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k);
      end
    end
    dcn(j,i) = dcn(j,i)/(x(j,i)*sum);
  end
end
%%%%%%%%%% FE-ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [U]=FE(nelx,nely,x,penal)
[KE] = lk; 
K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));
F = sparse(2*(nely+1)*(nelx+1),1); U = zeros(2*(nely+1)*(nelx+1),1);
for elx = 1:nelx
  for ely = 1:nely
    n1 = (nely+1)*(elx-1)+ely; 
    n2 = (nely+1)* elx   +ely;
    edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2];
    K(edof,edof) = K(edof,edof) + x(ely,elx)^penal*KE;
  end
end
% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
F(2,1) = -1;
fixeddofs   = union([1:2:2*(nely+1)],[2*(nelx+1)*(nely+1)]);
alldofs     = [1:2*(nely+1)*(nelx+1)];
freedofs    = setdiff(alldofs,fixeddofs);
% SOLVING
U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:);      
U(fixeddofs,:)= 0;
%%%%%%%%%% ELEMENT STIFFNESS MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [KE]=lk
E = 210.; 
nu = 0.3;
k=[ 1/2-nu/6   1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 ... 
   -1/4+nu/12 -1/8-nu/8  nu/6       1/8-3*nu/8];
KE = E/(1-nu^2)*[ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8)
                  k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)
                  k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)
                  k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)
                  k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)
                  k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)
                  k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)
                  k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];

Any help would be appreciated.

Thank you!

[E]

Dear Anderson,

I would like to know how to set the design area to a circle

Answer 1

[Anderson Pereira]

Dear Patil,

Each node has two Degrees Of Freedom (DOF), one in the x-direction and other in the y-direction. The DOFs are labeled as

x-direction: (2·nn)-1

y-direction: (2·nn)

where nn is the node number.

As the load is applied in the y-direction of node 1, the correspondent DOF is 2*1 = 2. So, answering your question, line 79 in the code: F(2,1) = -1;, is a force of -1 applied in the y direction of node 1.

Best,

Anderson

[MATHEW KURIAN]

How can i change the angle or orientation of the force?

Say i want it to act at 45 degrees and in the +ve x & y plane.

[Anderson Pereira]

Let the node number where you want to apply the load be nn. Moreover, let theta be the angle from the x-axis measured anti-clockwise and P be the magnutude of the load. Thus, the load vector F my be written

F(2*nn-1,1) = P*cos(theta);

F(2*nn,1) = P*sin(theta);

[MATHEW KURIAN]

Thank You.

[Sagar Gupta Satrasala]

Hello,

I am also working on this project. Can you please explain what is happening in this part,

K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1)),

F = sparse(2*(nely+1)*(nelx+1),1), U = zeros(2*(nely+1)*(nelx+1),1) in FE-Analysis.

[Anderson Pereira]

Hi Sagar,

This is just the creation of the MATLAB sparse arrays.

K is the global stiffness matrix

F is the global force vector

2*(nelx+1)*(nely+1) is the total number of degrees of freedom (DOF)

Best,

Anderson

[Sagar Gupta Satrasala]

Hello,

I know it is a global stiffness matrix but i dont understand how this (2*(nelx+1)*(nely+1)) works.

Can you elaborate the explanation.

Thank you

[Anderson Pereira]

Hi,

There is no much to add.

(2*(nelx+1)*(nely+1)) is a scalar, the size of the arrays.

Anderson

[Sagar Gupta Satrasala]

Hi,

Can you tell me where the dimentions(length/width/thickness) of th beam are given in this code.

I need to change the thickness values for better stresses.

[Ahmed Adel]

so how can i change the boundary conditions to solve those ploblems?

[Majeed]

Sugar,

Regarding your question about the number of DOFs. there is nelx * nely element in the body, and in each row there is nelx+1 node. Also, the same in each column nely+1 node. Therefore, the total number of the DOF is equal to the number of nodes x2.

Regards,,

[Majeed]

Sugar,

Also there is no thickness in this problem since it is 2D.

[Anas Hadni]

@Anderson Pereira Dear Aderson,

Why are DOFs defined as such? Since each node has 2 DOF, what happens for node 5 for example? following what you wrote that would mean than in the x direction there would be a 9th DOF. Obviously I am misunderstanding something...

Also, I feel like with only two inputs for F, we do not have enough information to fully model the force. If the first input is the DOF, and the second one is the node, how do we know if the first input is a DOF in the x or y direction?

I have real trouble understanding the way this code is written, I don't know if these are conventions that I don't know about or if it was the author's choice...

Thank you in advance!

[Anderson Pereira]

Hi Anas,

Some comments bellow:

Q#1: Why are DOFs defined as such?

A: They are defined as such to avoid storing them explicitly.

Q#2: Since each node has 2 DOF, what happens for node 5 for example? following what you wrote that would mean than in the x direction there would be a 9th DOF. Obviously I am misunderstanding something...

A:You understand right, node 5 will have the dofs 9 and 10, x and y direction, respectively.

Q#3:Also, I feel like with only two inputs for F, we do not have enough information to fully model the force. If the first input is the DOF, and the second one is the node, how do we know if the first input is a DOF in the x or y direction?

A: F is a vector with size 2*(number of nodes) x one (column vector). Each "i" position of "F" corresponds to a prescribed load. If the load is applied in the x-direction of node 5 we need to use the following:

F(9) = load.

Q#4:I have real trouble understanding the way this code is written, I don't know if these are conventions that I don't know about or if it was the author's choice...

A: Many matlab based FEM codes use the same idea.

Best,

Anderson

Answer 2

[Anderson Pereira]

Hi Sagar,

The finite element model adopted is independent of inplane dimension scalings and have unit thickness. Moreover, an isotropic material with with elastic modulus E and Poisson’s ratio nu was considered.

Best,

Anderson

[Sagar Gupta Satrasala]

Is there a way to introduce plane dimensions in this code?

[Anderson Pereira]

You can only define the aspect ratio of the structure which is given by the ratio of elements in the horizontal (nelx) and the vertical direction (nely).

[Rajnandini Das]

Hello, Can you please explain me how did they get the equation for n1 and n2 and also for ue?

Answer 3

[Piyali Dey]

Hello, I have been working on a project to obtain a topology of a compliant mechanism. Can you help me change the objective to maximum displacement.