# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020

################################Probelms###################
function inverter(nelx,nely,nelz)
    ndof = 3*(nelx+1)*(nely+1)*(nelz+1)
    # USER-DEFINED LOAD DOFs
    il = [0 nelx]; 
    jl = [nely nely]; 
    kl = [0 0];
    loadnid = kl .*(nelx+1) .*(nely+1) .+il .*(nely .+1) .+(nely+1 .-jl);
    loaddof = 3 .*loadnid[:] .- 2;

    loaddof=[getIndex3d(1,1,1,1,nelx+1,nely+1,nelz+1,3),
        getIndex3d(nelx+1,1,1,1,nelx+1,nely+1,nelz+1,3)]
    din = loaddof[1]; dout = loaddof[2];
    F = sparse(loaddof,[1;2],[1;-1],ndof,2);

    # USER-DEFINED SUPPORT FIXED DOFs
    # Top face
    m = Matlab.meshgrid(0:nelx,0:nelz)
    iif=m[1]
    kf=m[2]
    fixednid = kf .*(nelx+1) .*(nely+1) .+iif .*(nely+1) .+1;
    fixeddof_t = 3 .*fixednid[:] .-1;
    # Side face
    m = Matlab.meshgrid(0:nelx,0:nely)
    iif=m[1]
    jf=m[2]
    fixednid = iif .*(nely+1) .+(nely+1 .-jf);
    fixeddof_s = 3*fixednid[:];
    # Two pins
    iif = [0;0]; 
    jf = [0;1]; 
    kf = [nelz;nelz];
    fixednid = kf .*(nelx+1) .*(nely .+1) .+iif .*(nely .+1) .+(nely+1 .-jf)


    fixeddof_p = [3 .*fixednid[:]; 3 .*fixednid[:].-1; 3 .*fixednid[:].-2] # DOFs
    # Fixed DOFs
    fixeddof = union(fixeddof_t,fixeddof_s);
    fixeddof = union(fixeddof, fixeddof_p);
    fixeddof = sort(fixeddof);

    # Displacement vector
    U = zeros(ndof,2);
    freedofs = setdiff(1:ndof,fixeddof);

    return U,F,freedofs,din,dout
end

function wing(nelx,nely,nelz)
    ndof = 3*(nelx+1)*(nely+1)*(nelz+1)
    # USER-DEFINED LOAD DOFs
    # loaddof=[getIndex3d1(1,nely+1,1,1,nelx+1,nely+1,nelz+1,3),
    #     getIndex3d1(nelx+1,nely+1,1,1,nelx+1,nely+1,nelz+1,3)]
    # din = loaddof[1]; 
    # dout = loaddof[2];

    # F = sparse(loaddof,[1;2],[1;-1],ndof,2);

    loaddof=[getIndex3d(1,nely+1,1,1,nelx+1,nely+1,nelz+1,3),
        getIndex3d(nelx+1,nely+1,1,2,nelx+1,nely+1,nelz+1,3)]
    din = loaddof[1]; 
    dout = loaddof[2];

    F = sparse(loaddof,[1;2],[1;1],ndof,2);

    # F = fill(0.0,ndof,2);

    # F[loaddof[1],1]=1
    # F[loaddof[2],2]=-1
    # F[getIndex3d(1,1,1,2,nelx+1,nely+1,nelz+1,3):3:getIndex3d(nelx+1,nely+1,nelz+1,2,nelx+1,nely+1,nelz+1,3),1] .+= -0.01
    # F[getIndex3d(1,1,1,2,nelx+1,nely+1,nelz+1,3),1] = 1
    # F[loaddof[1],1].=F[loaddof[1],1].+1.0
    # F = sparse(loaddof,[1;2],[1;1],ndof,2);
    # F[:]


    # USER-DEFINED SUPPORT FIXED DOFs
    # m= Matlab.meshgrid(0:0,0:nely,0:nelz)# Coordinates
    # iif=m[1]
    # jf=m[2]
    # kf=m[3]
    # fixednid = kf.*(nelx+1).*(nely+1) .+iif .*(nely .+1) .+(nely .+1 .-jf) # Node IDs
    # fixeddof = [3 .*fixednid[:]; 3 .*fixednid[:].-1; 3 .*fixednid[:].-2] # DOFs

    # Side face
    m = Matlab.meshgrid(0:nelx,0:nely)
    iif=m[1]
    jf=m[2]
    fixednid = iif .*(nely+1) .+(nely+1 .-jf);
    fixeddof_s = 3*fixednid[:];

    fixeddof=union(fixeddof_s,getIndex3d(1,1,(nelz+1),1:3,nelx+1,nely+1,nelz+1,3))
    fixeddof=union(fixeddof,getIndex3d(1,2,(nelz+1),1:3,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(getIndex3d(1,1,1:(nelz+1),1,nelx+1,nely+1,nelz+1,3),getIndex3d(1,1,1:(nelz+1),2,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(fixeddof,getIndex3d(1,1,1:(nelz+1),3,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(fixeddof,getIndex3d(1,2,1:(nelz+1),1,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(fixeddof,getIndex3d(1,2,1:(nelz+1),2,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(fixeddof,getIndex3d(1,2,1:(nelz+1),3,nelx+1,nely+1,nelz+1,3))
    fixeddof = sort(fixeddof);


    # Displacement vector
    U = zeros(ndof,2);
    freedofs = setdiff(1:ndof,fixeddof);

    return U,F,freedofs,din,dout
end

function wing1(nelx,nely,nelz)
    ndof = 3*(nelx+1)*(nely+1)*(nelz+1)
    # USER-DEFINED LOAD DOFs
    # loaddof=[getIndex3d1(1,nely+1,1,1,nelx+1,nely+1,nelz+1,3),
    #     getIndex3d1(nelx+1,nely+1,1,1,nelx+1,nely+1,nelz+1,3)]
    # din = loaddof[1]; 
    # dout = loaddof[2];

    # F = sparse(loaddof,[1;2],[1;-1],ndof,2);

    loaddof=[getIndex3d(nelx/2,nely+1,1,1,nelx+1,nely+1,nelz+1,3),
        getIndex3d(nelx+1,nely+1,1,2,nelx+1,nely+1,nelz+1,3)]
    din = loaddof[1]; 
    dout = loaddof[2];

    F = sparse(loaddof,[1;2],[-1;-1],ndof,2);


    # Side face
    m = Matlab.meshgrid(0:nelx,0:nely)
    iif=m[1]
    jf=m[2]
    fixednid = iif .*(nely+1) .+(nely+1 .-jf);
    fixeddof_s = 3*fixednid[:];

    fixeddof=union(fixeddof_s,getIndex3d(1,1,(nelz+1),1:3,nelx+1,nely+1,nelz+1,3))
    fixeddof=union(fixeddof,getIndex3d(1,2,(nelz+1),1:3,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(getIndex3d(1,1,1:(nelz+1),1,nelx+1,nely+1,nelz+1,3),getIndex3d(1,1,1:(nelz+1),2,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(fixeddof,getIndex3d(1,1,1:(nelz+1),3,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(fixeddof,getIndex3d(1,2,1:(nelz+1),1,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(fixeddof,getIndex3d(1,2,1:(nelz+1),2,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(fixeddof,getIndex3d(1,2,1:(nelz+1),3,nelx+1,nely+1,nelz+1,3))
    fixeddof = sort(fixeddof);


    # Displacement vector
    U = zeros(ndof,2);
    freedofs = setdiff(1:ndof,fixeddof);

    return U,F,freedofs,din,dout
end

function wing2(nelx,nely,nelz)
    ndof = 3*(nelx+1)*(nely+1)*(nelz+1)
    # USER-DEFINED LOAD DOFs
    # loaddof=[getIndex3d1(1,nely+1,1,1,nelx+1,nely+1,nelz+1,3),
    #     getIndex3d1(nelx+1,nely+1,1,1,nelx+1,nely+1,nelz+1,3)]
    # din = loaddof[1]; 
    # dout = loaddof[2];

    # F = sparse(loaddof,[1;2],[1;-1],ndof,2);

    loaddof=[getIndex3d(nelx+1,nely+1,1,2,nelx+1,nely+1,nelz+1,3),getIndex3d(1,nely+1,1,1,nelx+1,nely+1,nelz+1,3)]
    din = loaddof[1]; 
    dout = loaddof[2];

    F = sparse(loaddof,[1;2],[2;1],ndof,2);

    # F = fill(0.0,ndof,2);

    # F[loaddof[1],1]=1
    # F[loaddof[2],2]=-1
    # F[getIndex3d(1,1,1,2,nelx+1,nely+1,nelz+1,3):3:getIndex3d(nelx+1,nely+1,nelz+1,2,nelx+1,nely+1,nelz+1,3),1] .+= -0.01
    # F[getIndex3d(1,1,1,2,nelx+1,nely+1,nelz+1,3),1] = 1
    # F[loaddof[1],1].=F[loaddof[1],1].+1.0
    # F = sparse(loaddof,[1;2],[1;1],ndof,2);
    # F[:]


    # USER-DEFINED SUPPORT FIXED DOFs
    # m= Matlab.meshgrid(0:0,0:nely,0:nelz)# Coordinates
    # iif=m[1]
    # jf=m[2]
    # kf=m[3]
    # fixednid = kf.*(nelx+1).*(nely+1) .+iif .*(nely .+1) .+(nely .+1 .-jf) # Node IDs
    # fixeddof = [3 .*fixednid[:]; 3 .*fixednid[:].-1; 3 .*fixednid[:].-2] # DOFs

    # Side face
    m = Matlab.meshgrid(0:nelx,0:nely)
    iif=m[1]
    jf=m[2]
    fixednid = iif .*(nely+1) .+(nely+1 .-jf);
    fixeddof_s = 3*fixednid[:];

    fixeddof=union(fixeddof_s,getIndex3d(1,1,(nelz+1),1:3,nelx+1,nely+1,nelz+1,3))
    fixeddof=union(fixeddof,getIndex3d(1,2,(nelz+1),1:3,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(getIndex3d(1,1,1:(nelz+1),1,nelx+1,nely+1,nelz+1,3),getIndex3d(1,1,1:(nelz+1),2,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(fixeddof,getIndex3d(1,1,1:(nelz+1),3,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(fixeddof,getIndex3d(1,2,1:(nelz+1),1,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(fixeddof,getIndex3d(1,2,1:(nelz+1),2,nelx+1,nely+1,nelz+1,3))
    # fixeddof=union(fixeddof,getIndex3d(1,2,1:(nelz+1),3,nelx+1,nely+1,nelz+1,3))
    fixeddof = sort(fixeddof);


    # Displacement vector
    U = zeros(ndof,2);
    freedofs = setdiff(1:ndof,fixeddof);

    return U,F,freedofs,din,dout
end