ellipt_compliant_detailed_simulation.m

function sim_data = ellipt_compliant_detailed_simulation(model,model_data)
%function sim_data = ellipt_compliant_detailed_simulation(model,model_data)
%
% computation of FEM solution for elliptic problem
% input: model_data, parameter independent contains grid, p,e,t,
%        model_data.dl
%        model.B1/B2
%        model.gammatop/bottom
%
% output: 
%      sim_data.bm: boundary condition matrix, dependent of mu
%      sim_data.u : vector of nodal values of solution in gridpoints
%      sim_data.s : real value, functional output of problem 
%       sim_data.plot_title: title of the plot
%           
%      following data are representation of the elliptic pde
%           -grad(c*grad u)+a*u=f
%      sim_data.a: 'reaction' term.
%      sim_data.f: 'source' term.
%      sim_data.c: 'diffusivity' term.
%      sim_data.A_matlab: left hand side, computed with assempde
%      sim_data.F_matlab: right hand side, computed with assempde
%      sim_data.M_matlab: Stiffness matrix
%      sim_data.K_matlab: Stiffness matrix
%      sim_data.F_matlab: Stiffness matrix
%      sim_data.Q_matlab: boundary condition (Neumann)
%      sim_data.G_matlab: boundary condition (Neumann)
%      sim_data.H_matlab: boundary condition (Dirichlet)
%      sim_data.B_matlab: boundary condition (Dirichlet)


% B. Haasdonk 21.7.2009
% S. Langhof 



sim_data.plot_title='Plot of a detailed simulation';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%create boundary matrix

%set boundary condition
%bottom+top: from left to right
%left+right: from bottom to top

dl=model_data.dl;

noc_bm=length(dl);                          %number of columns

check_left_region = dl(6,:);                %number of block left hand side
check_right_region = dl(7,:);               
check_region = dl(6:7,:);

zeros_left=find(check_left_region == 0);          %find 0's. this edge is an outer
zeros_right=find(check_right_region == 0);        %boundary

nonzeros_left=find(check_left_region);
nonzeros_right=find(check_right_region);

inner_boundary_index=intersect(nonzeros_left, nonzeros_right);  %index of non zero columns in dl

zero_index=find(check_region==0);
outer_boundary_index=ceil(zero_index/2);
outer_boundary_index=outer_boundary_index';

check_max_length = zeros(4,max(model.B1,model.B2));
for i=1:model.B1
    check_max_length(1,i)=length(model.gamma_bottom{i}{1})+length(model.gamma_bottom{i}{2})+4;          %bottom
    check_max_length(2,i)=length(model.gamma_top{i}{3})+length(model.gamma_top{i}{4})+8;                %top
end
for l=1:model.B2
    check_max_length(3,l)=length(model.gamma_left{l}{1})+length(model.gamma_left{l}{2})+4;              %left
    check_max_length(4,l)=length(model.gamma_right{l}{1})+length(model.gamma_right{l}{2})+4;            %right
end

nor_bm=max(max(check_max_length));              %number of rows BM


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Boundary matrix
BM=zeros(nor_bm, noc_bm);               %zeros BM Matrix


%%%%% x1/x2, y1/y2 values
check_y1=dl(4,:);
check_y2=dl(5,:);

check_x1=dl(2,:);
check_x2=dl(3,:);

%find indices of gamma top
one_y1_index=find(check_y1 == 1);
one_y2_index=find(check_y1 == 1);

one_y_index = intersect(one_y1_index, one_y2_index);

check_x1_one = check_x1(one_y_index);

top_tmp=0;
[tmp, top_tmp]=sort(check_x1_one);

k=0;
top_index=0;                            %
for l=top_tmp 
    k=k+1;
    top_index(k)=one_y_index(l);   %top_index: top indices from left to right
end

%find indices of gamma bottom

zero_y1_index=find(check_y1 == 0);
zero_y2_index=find(check_y2 == 0);

zero_y_index= intersect(zero_y1_index, zero_y2_index);

check_x1_zero=check_x1(zero_y_index);

bottom_tmp=0;
[tmp, bottom_tmp]=sort(check_x1_zero);
k=0;
bottom_index=0;                            %
for l=bottom_tmp 
    k=k+1;
    bottom_index(k)=zero_y_index(l);   %bottom index: bottom indices from left to right
end

%find left index

zero_x1_index=find(check_x1 == 0);
zero_x2_index=find(check_x2 == 0);

zero_x_index=intersect(zero_x1_index, zero_x2_index);

check_y1_zero=check_y1(zero_x_index);      %

left_tmp=0;
[tmp, left_tmp]=sort(check_y1_zero);
k=0;
left_index=0;                            %
for l=left_tmp 
    k=k+1;
    left_index(k)=zero_x_index(l);   %left index: left indices from left to right
end

%find right index

one_x1_index=find(check_x1 == 1);
one_x2_index=find(check_x2 == 1);

one_x_index=intersect(one_x1_index, one_x2_index);

check_y1_one=check_y1(one_x_index);      %

right_tmp=0;
[tmp, right_tmp]=sort(check_y1_one);
k=0;
right_index=0;                            %
for l=right_tmp 
    k=k+1;
    right_index(k)=one_x_index(l);   %right index: right indices from left to right
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%final creation of boundaries

%outer bound indices
outer_bound=[top_index, bottom_index, right_index, left_index];        
outer_bound=sort(outer_bound);

%inner bound indices
tmp_noc_vec=1:noc_bm;
for i=outer_bound
    tmp_noc_vec(i)=-1;
end
inner_bound=find(tmp_noc_vec >=0);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
%%%inner boundaries
std_inner_boundary_char=['0','0','1','0'];        %double(...) -> 48 48 49 48
std_inner_boundary_num=[0,repmat(1,1,5)];

std_inner_boundary=[std_inner_boundary_num,double(std_inner_boundary_char)];

if length(std_inner_boundary) < nor_bm              %fill up with zeros
    std_inner_boundary(nor_bm);
end

%positioning of inner boundary condition
for i=inner_bound
    BM(:,i)=std_inner_boundary(:);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%outer boundaries

%top boundary

%max length of vec
%l_top_bound=max(check_max_length(2,:));

%TOP=zeros(l_top_bound, model.B1);
clear fronttail
fronttail=repmat(1,1,4);        %definition of dim, diri, lq, lg

for i=1:model.B1
    lh=length(model.gamma_top{i}{3});        %number of chars for h
    lr=length(model.gamma_top{i}{4});        %number of chars for r
    outer_topboundary_num(i,:)=[fronttail,lh, lr, 48, 48];
end

outer_topboundary_num=outer_topboundary_num';        

tmp_char={};
for i=1:model.B1
    tmp_char3=model.gamma_top{i}{3};
    tmp_char4=model.gamma_top{i}{4};

    tmp_char{i}=[tmp_char3,tmp_char4];          %
end

for i=1:model.B1
    outer_topboundary(:,i)=[outer_topboundary_num(:,i);double(tmp_char{i})'];
    if length(outer_topboundary(:,i)) < nor_bm
        outer_topboundary(nor_bm,i) = 0;
    end
end


%positioning of top boundary conditions
outer_topbound_loop=outer_topboundary;
for i=top_index
        BM(:,i)=outer_topbound_loop(:,1);
        outer_topbound_loop(:,1)=[];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%bottom boundaries

clear fronttail
fronttail=[1,0];        %definition of dim, diri, lq, lg

for i=1:model.B1
    lq=length(model.gamma_bottom{i}{1});        %number of chars for 
    lg=length(model.gamma_bottom{i}{2});        %number of chars for 
    outer_bottomboundary_num(i,:)=[fronttail,lq, lg];
end

outer_bottomboundary_num=outer_bottomboundary_num';        

tmp_char={};
for i=1:model.B1
    tmp_char1=model.gamma_bottom{i}{1};
    tmp_char2=model.gamma_bottom{i}{2};

    tmp_char{i}=[tmp_char1,tmp_char2];          %
end

a=[];
for i=1:model.B1
    a=[];
    a=[outer_bottomboundary_num(:,i);double(tmp_char{i})'];
    if length(a) < nor_bm
        a(nor_bm);
    end
    outer_bottomboundary(:,i)=a;
    
%     outer_bottomboundary(:,i)=[outer_bottomboundary_num(:,i);double(tmp_char{i})'];
%     if length(outer_bottomboundary(:,i)) < nor_bm
%         outer_bottomboundary(nor_bm,i) = 0;
%     end
end

%positioning of bottom boundary conditions
outer_bottombound_loop=outer_bottomboundary;
for i=bottom_index
        BM(:,i)=outer_bottombound_loop(:,1);
        outer_bottombound_loop(:,1)=[];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%left boundaries
clear fronttail
fronttail=[1,0];        %definition of dim, diri, lq, lg

for i=1:model.B2
    lq=length(model.gamma_left{i}{1});        %number of chars for 
    lg=length(model.gamma_left{i}{2});        %number of chars for 
    outer_leftboundary_num(i,:)=[fronttail,lq, lg];
end

outer_leftboundary_num=outer_leftboundary_num';        

tmp_char={};
for i=1:model.B2
    tmp_char1=model.gamma_left{i}{1};
    tmp_char2=model.gamma_left{i}{2};

    tmp_char{i}=[tmp_char1,tmp_char2];          %
end

a=[];
for i=1:model.B2
    a=[];
    a=[outer_leftboundary_num(:,i);double(tmp_char{i})'];
    if length(a) < nor_bm
        a(nor_bm);
    end
    outer_leftboundary(:,i)=a;
    
%     outer_bottomboundary(:,i)=[outer_bottomboundary_num(:,i);double(tmp_char{i})'];
%     if length(outer_bottomboundary(:,i)) < nor_bm
%         outer_bottomboundary(nor_bm,i) = 0;
%     end
end

%positioning of left boundary conditions
outer_leftbound_loop=outer_leftboundary;
for i=left_index
        BM(:,i)=outer_leftbound_loop(:,1);
        outer_leftbound_loop(:,1)=[];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%right boundaries
clear fronttail
fronttail=[1,0];        %definition of dim, diri, lq, lg

for i=1:model.B2
    lq=length(model.gamma_right{i}{1});        %number of chars for 
    lg=length(model.gamma_right{i}{2});        %number of chars for 
    outer_rightboundary_num(i,:)=[fronttail,lq, lg];
end

outer_rightboundary_num=outer_rightboundary_num';        

tmp_char={};
for i=1:model.B2
    tmp_char1=model.gamma_right{i}{1};
    tmp_char2=model.gamma_right{i}{2};

    tmp_char{i}=[tmp_char1,tmp_char2];          %
end

a=[];
for i=1:model.B2
    a=[];
    a=[outer_rightboundary_num(:,i);double(tmp_char{i})'];
    if length(a) < nor_bm
        a(nor_bm);
    end
    outer_rightboundary(:,i)=a;
    
%     outer_bottomboundary(:,i)=[outer_bottomboundary_num(:,i);double(tmp_char{i})'];
%     if length(outer_bottomboundary(:,i)) < nor_bm
%         outer_bottomboundary(nor_bm,i) = 0;
%     end
end

%positioning of left boundary conditions
outer_rightbound_loop=outer_rightboundary;
for i=right_index
        BM(:,i)=outer_rightbound_loop(:,1);
        outer_rightbound_loop(:,1)=[];
end

%
sim_data.bm=BM;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% falls stueckweise konstante Diffusivitaet:
% modulo sortierung:
%
% c = ''
% for p=1:params.B1*params.B2
%    mu_p = getfield(model,['mu',num2str(p)])
%    c = [c,'!',numstr(mu_p)]


% allgemeimer:
%  c = richtige_kombination_und_parameter von model.diffusitivy_ptr(...)


%%%%COM Berechnung
%for i=1:length(t)
%    tmp_t=t(1:3,i);
%    tmp_p=p(:,tmp_t);
%    tmp_sum=sum(tmp_p,2);        
%    com_t{i}=tmp_sum/3;         %CELL ARRAY!!
%end

% glob = aus com_t konstruieren;


%%%%%simplify p,e,t call

p=model_data.grid.p;
e=model_data.grid.e;
t=model_data.grid.t;
% Mittelwert von allen Dreiecken:
ntria = size(t,2);
tp = t(1:3,:);
tp(:);
px= p(1,tp);
px = reshape(px,3,ntria);
cx = mean(px)';
% selbes spielchen fuer py
py=p(2,tp);
py = reshape(py,3,ntria);
cy = mean(py)';
glob = [cx,cy];
% Diffusiviy:
c = model.diffusivity(glob,model);
sim_data.c=c';



% Source-Term: f = model.source(glob,model);

% ....

f=zeros(1,length(c));
sim_data.f=f;

% Reaction: a = model.reaction(glob,model);

a=zeros(1,length(c));
sim_data.a=a;

% possible to use dir_vals = model.dirichlet_values(dglob,model)

% .........
% dglob = ...
% dir_vals = model.dirichlet_values(dglob,model)
% BM = dir_vals ..;

% possible to use model.neumann_values(nglob,model)   % g_N(x)
% .........
% nglob = ...
% neu_vals = model.neumann_values(bglob,model)
% BM = neu_vals ..;

% assempde and solve

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%ASSEMBLING+solving%%%%%%%%%%%%%%
% assemblieren:
A=[];
F=[];
p=model_data.grid.p;
e=model_data.grid.e;
t=model_data.grid.t;
a=sim_data.a;
c=sim_data.c;
f=sim_data.f;
BM=sim_data.bm;

%%assembling system with 'stiff springs' (dirichlet conditions)
% [A,F]=assempde(BM,p,e,t,c,a,f);               %A stiffness matrix, F right-hand side
% u=A\F;
% sim_data.u=u;
% 
% sim_data.A_matlab_assempde=A;
% sim_data.F_matlab_assempde=F;
%%%assembling system, first eliminating diri condi
B=[];
ud=[];
u_tmp=[];
[A,F,B,ud]=assempde(BM,p,e,t,c,a,f);
u_tmp=A\F;
u=B*u_tmp+ud;

Afull = speye(size(p,2));
non_diri_values=find(p(2,:)~=1);
Afull(non_diri_values,non_diri_values) = A;

Bfull = ud;
Bfull(non_diri_values)= F;
ufull = Afull\Bfull;

%Afull = A * B';
%Bfull = F + A * B' * ud;


% A2 = B * A
% B 

%keyboard;

sim_data.u=ufull;
sim_data.B = Bfull;
sim_data.A = Afull;
% eventuell ueberfluessig:
sim_data.ud=ud;
sim_data.u_tmp=u_tmp;
sim_data.A_matlab=A;
sim_data.F_matlab=F;
sim_data.B_matlab=B;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %split into assembly of area int contribution + boundary condi contribution
K=[];           %stiffness matrix
M=[];           %mass matrix
F=[];           %right-hand side
[K,M,F]=assema(p,t,c,a,f);
%
Q=[];
G=[];
H=[];
B=[];
[Q,G,H,B]=assemb(BM,p,e);

sim_data.K_matlab = K;
sim_data.M_matlab = M;
sim_data.F_matlab = F;

sim_data.Q_matlab = Q;
sim_data.G_matlab = G;
sim_data.H_matlab = H;
sim_data.B_matlab = B;
% %
% %%%%%%% 3 different types
% %u=assempde(K,M,F,Q,G,H,B);
% %sim_data.u=u;
% %%%%%%%%%%%%%%%%%%%%%%
% K1=[];
% F1=[];
% [K1,F1]=assempde(K,M,F,Q,G,H,B);
% u_split_matlab=K1\F1;
% sim_data.A_split_assempde = K1;
% sim_data.F_split_assempde = F1;
% sim_data.u_split_matlab=u_split_matlab;
% %%%%%%%%%%%%%%%%%%%%%%%
% K2=[];
% F2=[];
% B=[];
% ud=[];
% u_tmp=[];
% [K2,F2,B,ud]=assempde(K,M,F,Q,G,H,B);
% u_tmp=K2\F2;
% u=B*u_tmp+ud;
% sim_data.K2_no_stiff = K2;
% sim_data.F2_no_stiff = F2;
% sim_data.B_no_stiff = B;
% sim_data.ud_no_stiff = ud;
% sim_data.u_no_stiff_springs=u;



%old
%% A = ...;
%% f = ...;
%u=assempde(BM,p,e,t,c,a,f);
%% u = A\f;
%sim_data.u=u;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output computation: Integral over boundary segment

%sim_data.s=thermalblock_output_function(model.plot_output,model_data)

%sim_data.s=thermalblock_output_function(glob,model)

%%shortcuts
p=model_data.grid.p;
e=model_data.grid.e;
t=model_data.grid.t;

% ausarbeiten des folgenden:

calc_mode = model.output_calc_mode;

%calc_mode = 2; % 2 only lower edge, 1 all edges

if calc_mode == 1
ei6 = find(e(6,:)==0);  %indices of outer boundary
ei7 = find(e(7,:)==0);  %edges
ei = [ei6,ei7];
integral = 0;
% glob = Mittelpunkt von Kante / midpoint of edge

oe=e(:,ei);     %outer boundary edges
oestart=oe(1,:);        %index of starting point
oeend=oe(2,:);  %index of ending point

poestart=p(:,oestart);  %x,y coordinates of correspondign points
poeend=p(:,oeend);              %      ''               ''

glob=0.5.*(poestart+poeend)';

weight = model.output_weight_function(glob,model);
% %%  eventuell indexmenge von nichtnull-weights bestimmen
u_start = sim_data.u(oestart);
u_end = sim_data.u(oeend);
u_values_in_midpoints = 0.5*(u_start + u_end);          %column
edge_lengths = sqrt(sum((poeend-poestart).^2)); %row 
edge_lengths=edge_lengths(:);           %column
integral = edge_lengths.*u_values_in_midpoints.*weight;

sim_data.s=sum(integral);

else % calc_mode ==2                    %only lower edge

zero_y_ind=find(p(2,:)==0);             %find indices with y=0

x_values=p(1,:);
x_value=x_values(zero_y_ind);   %corresponding x values with y=0

[x_value_sort,i]=sort(x_value); %sort vector
glob = [x_value_sort(:),zeros(length(x_value),1)];


u_zero_y=u(zero_y_ind);         %u values at y=0
u_zero_y_sort=u_zero_y(i);              %corresponding u values sorted

%the following section tests, if the vector x_value contains 0,1, because
%we want to compute the integral in the intervall [0,1]. 
if min(x_value) ~= 0
        warning('x=0 not found!!! extrapolating values!!')
        x_value=[0,x_value];
        u_zero=interp1(x_value_sort, u_zero_y_sort,0, 'cubic');
        u_zero_y_sort=[u_zero,u_zero_y_sort];
end

if max(x_value) ~= 1
        warning('x=1 not found!!! extrapolating values!!')      
        x_value=[x_value ,1];
        u_one = interp1(x_value_sort,u_zero_y_sort,1,'cubic');
        u_zero_y_sort=[u_zero_y_sort,u_one];
end

weight_values = model.output_weight_function(glob,model); 
integrand_values = u_zero_y_sort.* weight_values;
sim_data.s = trapz(x_value_sort, integrand_values);

end; 

%if model.plot_output ==1
%       figure;
%       plot(x_value_sort, u_zero_y_sort)
%end