[求助] matlab 热力耦合分析

function [freeDOF] = boundaryCond(DOF)
% This function defines the boundary conditions to be applied.              %  case 1              

global BC  CONNEC DOMAIN NODES

nXElem   = DOMAIN(1);                                                       % Number of elements in the x-direction
nYElem   = DOMAIN(2);                                                       % Number of elements in the y-direction
fixedTDOF = NaN(1,(DOMAIN(1)+1)*(DOMAIN(2)+1));                              % Initialize vector of fixed DOFs
fixedDOF = NaN(1,(DOMAIN(1)+1)*(DOMAIN(2)+1));   
TDOF     = 1*max(max(NODES));
index = 1;
switch BC
    case 1                                                                  % Fix the temperature in bottom and top edge of the domain
        
      
              for i = 1:nXElem+1;
                if i == 1;
                fixedTDOF(index) = node;                                        % Fix temperature in bottom edge of the domain
                fixedTDOFvalue=0;                                                  %    bcval2 = 0;
                index = index+1;
                 node = node+1;                                             
               
                 elseif node == (nYElem)*(nXElem+1)+1;
                fixedTDOF(index) = node;                                        %  bcval1 = -10;
                fixedTDOFvalue=-10;
                index = index+1;
                 node = node+1;
                                                            
                end  
               
             end
        
   
         
         
         
    case 2                                                                  % Boundary conditions for an edge crack
        fixedDOF = [2*(nXElem+1)-1 2*(nXElem+1) 2*(nXElem+1)*(nYElem+1)-1];
   
    case 3                                                                  % Fix the bottom edge of the domain
        for i = 1:nXElem+1
            if i == 1
                fixedDOF(index:(index+1)) = [2*i-1 2*i];                    % Fix displacement in x and y-direction
                index = index+2;
            else
                fixedDOF(index) = 2*i;                                      % Fix displacement in y-direction
                index = index+1;
            end
        end            
    case 4                                                                  % Boundary conditions for first quadrant
        node = 1;
        for i = 1:nYElem+1                                                  % Rollers left edge
            fixedDOF(index) = 2*node-1;
            index = index+1;
            node = node+nXElem+1;
        end               
        for i = 1:nXElem+1
            fixedDOF(index) = 2*i;                                          % Fix displacement in y-direction
            index = index+1;
        end      
end
                                         
fixedTDOF(:,isnan(fixedTDOF(1,:))) = [];  
freeTDOF = setdiff(1:TDOF,fixedTDOF);
fixedDOF(:,isnan(fixedDOF(1,:))) = [];
freeDOF = setdiff(1:DOF,fixedDOF);                                          % Solve for the free DOF