
% The theoretical part of this example can be found in the OpenFEM manual
% "Application examples : Heat equation"

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Mesh the model - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

femesh('reset')
FEnode=[1 0 0 0  0 0 0; 2 0 0 0  0 5 0;
        3 0 0 0  0 5 4; 4 0 0 0  0 0 4];
FEel0=[Inf abs('quad4');1 2 3 4 1 1];
femesh(';divide 4 5;addsel;');
femesh(';selgroup1;extrude 10 1 0 0;orientel0'); % group 1: volume elements
FEelt=FEel0;
femesh(';selelt selface&withoutnode{x~=0};addsel'); % group 2: surface 
                                                    % elements (x=0)
femesh(';selelt selface&withoutnode{x~=10};addsel'); % group 3: surface 
                                                     % elements (x=10)
femesh(';selelt selface&withnode{x~=0}&withnode{x~=10};addsel'); % group 4: 
                                       % surface elements (all the other faces)
femesh('set group1 mat 1 pro 1')
femesh('set group2 mat 1 pro 2 name q4p')
femesh('set group3 mat 2 pro 2 name q4p')
femesh('set group4 mat 1 pro 2 name q4p')
model=femesh('model');
                                   
model.pl=[1 fe_mat('m_heat','SI',1) 400 1 1 1   % k rho C alpha  
          2 fe_mat('m_heat','SI',1) 400 1 1 1e10  % k rho C (alpha_0: not used)
          ];
model.il=[1 fe_mat('p_heat','SI',1) 0 2 3;   % volume integration
          2 fe_mat('p_heat','SI',2) 0 2 3];  % exchange surface integration

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Build the loading conditions on nodes of the interface
% This contains g+\alpha \theta_ext which will be integrated on the surface

% **** data ****
g1=25;
alpha=1; theta_ext=20;
f=40;
% **************

% boundary conditions on face x==0
data1 = struct('sel','group2',...
            'eltsel','group2',...
            'def',g1,'DOF',.20);
model = fe_case(model,'Fvol','Surface_load1',data1);
        % the surface integrals correspond to volume integrals on
        % the surface elements representing the boundary

% boundary conditions on face x==10
mdof0=feutil('findnode group3',model)+.20;
data2 = struct('DOF',mdof0,'def',ones(size(mdof0))*theta_ext);
model = fe_case(model,'DofSet','Surface_load2',data2);

% boundary conditions on other faces
node=feutil('getnode group4',model);
data3 = struct('sel','group4','def',25 - node(:,5).^2/20,'DOF',node(:,1)+.20);
model = fe_case(model,'Fvol','Surface_load3',data3);

% constant internal heat source
data4 = struct('sel','group1','def',[f],'DOF',.20);
model = fe_case(model,'Fvol','Internal source',data4);

% Assemble matrices and RHS - - - - - - - - - - - - - - - - - - - - - - - - -

DOF=feutil('getdof',model);unique(round(rem(DOF,1)*100))
[Case,model.DOF]=fe_mknl('init -gstate',model);

% try; if RunOpt.debug
%  sdtdef('diag',12);  % very slow elem0.m file implementation, just check
%  % the right-hand side is in dc.def(:,2) which is modified even though it is
%  % now returned by fe_mknl
%  dc2=struct('def',zeros(length(model.DOF),2),'DOF',model.DOF);
%  k=fe_mknl('assemble not',model,Case,dc2,1);
%  sdtdef('diag',0); % modified here
% end;end
sdtdef('diag',0); 

dc=struct('def',zeros(length(model.DOF),2),'DOF',model.DOF);
k=fe_mknl('assemble not',model,Case,dc,1);

% if (norm(dc2.def(:,2)-dc.def(:,2))>1e-6) error('on RHS'); 
% else disp('RHS OK'); end

% computation of the response to internal forces and enforced displacement
ind=fe_c(model.DOF,Case.DofIn,'ind');cind=setdiff(1:length(model.DOF),ind);
resp=struct('def',zeros(size(model.DOF)),'DOF',model.DOF);
resp.def(ind,1)=Case.TIn(ind,:);
F=dc.def(cind,2)-k(cind,ind)*Case.TIn(ind,:);
resp.def(cind,1)=ofact(k(cind,cind),F);

% View results  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

feplot(model,resp); try;fecom('colordata20');end

[i1,n1]=feutil('findnode z==2 & y==2',model);
figure(10); % (figure 4.1 of the manual)
plot(n1(:,5),fe_c(resp.DOF,i1)*resp.def);
%title('Temperature on the y=z=0 edge');
xlabel('x coordinate'); ylabel('temperature')

% This was used to generate the manual figure
% set(10,'position',[520 911 560 199]);
% convert('-pautops2pdf','../tex/plots/heat_eq.eps')

[i1,n1]=feutil('findnode x==5 & y==2',model);
figure(11);subplot(211);
plot(n1(:,7),fe_c(resp.DOF,i1)*resp.def)
title('Temperature at points x=5, y=2');

[i1,n1]=feutil('findnode x==5 & z==2',model);
subplot(212);plot(n1(:,6),fe_c(resp.DOF,i1)*resp.def)
title('Temperature at points x=5, z=2');

feutil('infoelt',model)


%	Etienne Balmes, Frederic Bourquin, A. Nassiopoulos
%       Copyright (c) 2001-2005 by INRIA and SDTools, All Rights Reserved.
%       Use under OpenFEM trademark.html license and LGPL.txt library license
%       $Revision: 1.11 $  $Date: 2006/01/25 15:53:41 $
