porsche_output_functional.m

function out_func = porsche_output_functional(model, model_data, uh)
%function out_func = porsche_output_functional(model, model_data, uh)
%
% function calculates the detailed output functional of the porsche model.
%
% The output-functional is:
% \int_{\Omega_{front,ref}} v \times t dx
% so the gradient of the solution (this is v) times the tangent on the
% front of the car.
%
% needed fields in model_data:
% model_data.front_ref_domain.linear_ind_front : linear indices of the
%                                                car-front
%
%
%
% Oliver Zeeb, 05.05.11


out_func=[];

%%%%%%%%%%%%%%%%%%%
% calculation of the gradient of the solution
%%%%%%%%%%%%%%%%%%%
% PASST DIE GRADIENTENBERECHNUNG????
%calculation of the gradient in the midpoints of the front-elements!
%since we need to calculate the gradient, which is not defined in the
%vertices of the grid, we need to get the midpoints of the elements
%and calculate the gradient in these points
% el_midpts_x=model_data.grid.CX(el_ind_front);
% el_midpts_y=model_data.grid.CY(el_ind_front);
%[uh_grad_front_x, uh_grad_front_y] = disc_gradient([el_midpts_x, el_midpts_y], uh);

%calculation of the gradient in the midpoints of the front-edges!
% edge_midpts_x=model_data.grid.ECX(linear_ind_front);
% edge_midpts_y=model_data.grid.ECY(linear_ind_front);
edge_length=model_data.front_ref_domain.edge_length_front;
edge_midpts=model_data.front_ref_domain.edge_midpts_front;
[uh_grad_front_x, uh_grad_front_y] = ...
    disc_gradient(edge_midpts, uh);
% PASST DIE GRADIENTENBERECHNUNG????
%%%%%%%%%%%%%%%%%%%

%Bei linearen finiten elementen ist der gradient auf dem gesamten element
%konstant!!! --> Ich darf in der Mitte auswerten und dann über den Rand
%integrieren.
% Daher die Ausgabe, dass nur pdeg=1 verwendet werden darf!
if(model.pdeg ~= 1)
    warning('Polynomial degree of Finite Elements is ~= 1!')
    warning('Output calculation is not correct!')
end

tx_front = model_data.front_ref_domain.tangent_front(:,1);
ty_front = model_data.front_ref_domain.tangent_front(:,2);
grad_times_t = tx_front.*uh_grad_front_x + ty_front.*uh_grad_front_y;
%integral as sum
out_func = edge_length' * grad_times_t;

% NOT NEEDED SO FAR:
%
%pressure=compute_pressure(model,uh_grad_front_x,uh_grad_front_y);
%
%X_front and Y_front are the vertices of the grid that belong to the
%front of the car:
% p_ind_front=model_data.grid.VI(linear_ind_front); %indices of the points belonging to the front
% X_front = model_data.grid.X(p_ind_front); 
% Y_front = model_data.grid.Y(p_ind_front); 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Tests zu dieser Funktion!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%TESTPLOT 1:
%FINDET ER DIE RICHTIGEN FRONT-ELEMENTE?!
%--> Frontelemente einzeln plotten!
%
% grid=model_data.grid;
% figure(2)
% hold on
% for k=1:length(el_ind_front)
%     x=grid.X(grid.VI(el_ind_front(k),:));
%     x=[x;grid.X(grid.VI(el_ind_front(k),1))]; %ersten punkt hinten hinhängen
%     
%     y=grid.Y(grid.VI(el_ind_front(k),:));
%     y=[y;grid.Y(grid.VI(el_ind_front(k),1))];
%     
%     plot(x,y)
% end
% axis equal
%ENDE: TESTPLOT 1

%TESTPLOT 2:
%FINDET ER DIE RICHTIGEN FRONT-PUNKTE?
%
% plot(X_front,Y_front,'.');
% axis equal
%ENDE: TESTPLOT 2

%TESTPLOT 3:
%PLOTTE DIE NORMALEN- UND TANGENTENVEKTOREN
%
% s = 1; %plot start
% e = length(nx_front); %plot end
% %e = 1;
% figure(1)
% quiver(model_data.grid.CX(el_ind_front(s:e)),model_data.grid.CY(el_ind_front(s:e)), nx_front(s:e), ny_front(s:e), 'blue')
% hold on
% quiver(model_data.grid.CX(el_ind_front(s:e)),model_data.grid.CY(el_ind_front(s:e)), tx_front(s:e), ty_front(s:e), 'red')
% axis equal
% legend('Normale','Tangente')
%ENDE: TESTPLOT 3


% TEST 1: Gradient in Normalenrichtung aufsummiert --> sollte eigentlich 0
% geben wegen Neumann-0-Randbedingungen
%
% grad_times_n = nx_front.*uh_grad_front_x + ny_front.*uh_grad_front_y;
% out_func_n = edge_length' * grad_times_n;
% NOCH SEHR HOCH?! Sollte doch eigentlich in der nähe von 0 liegen, wegen
% neumann randbedingungen...
% 645,7544 : bei originalgitter
% 340,7974 : einmal verfeinert in PDE-Toolbox
% 176,7265 : zwei mal verfeinert in PDE-Toolbox
%  91,1955 : drei mal verfeinert in PDE-Toolbox
% Fehler ist vor allem in knickpunkten sehr hoch!
%siehe dazu:
% x_plot=X_front(find(abs(grad_times_n)>5)); %wähle die x-punkte, in denen grad_times_n größer 10 ist
% y_plot=Y_front(find(abs(grad_times_n)>5)); % ... y-punkte ...
% plot(X_front,Y_front,'b.') 
% hold on
% plot(x_plot,y_plot,'rd')
%ENDE TEST 1