function model = exercise1_musterloesung_newton
% quadratic poisson equation on unit interval
%
% u on [0,1] such that
%
% - laplace (k u) + m u = 1
%    u(0) = u(1) = 0
%
% parameters k and m
%
% discretization with Finite differences
%
% inner product: l2-inner product, piecewise linear elements

% B. Haasdonk 6.12.2011

% for detailed_simulation
model.mu_names = {'k','m'};
model.mu_ranges = {[0.01,1],[0,10]};
model.xrange = [0,1];
model.xnumintervals = 1000;

model.dx = (model.xrange(2)-model.xrange(1))/model.xnumintervals;

model.k = 1;
model.m = 10;

% for detailed simulation
model.gen_model_data  = @my_gen_model_data;
model.detailed_simulation  = @my_detailed_simulation;
model.plot_sim_data = @my_plot_sim_data;
model.set_mu = @set_mu_default;
model.newton_residual = @my_newton_residual;
model.newton_residual_jacobian = @my_newton_residual_jacobian;
model.b_coefficients = @my_b_coefficients;
model.b_components = @my_b_components;
model.b_components_jacobian = @my_b_components;
model.f_coefficients = @my_f_coefficients;
model.f_components = @my_f_components;
model.a_coefficients = @my_a_coefficients;

model.a_components = @my_a_components;
model.a_components_jacobian = @my_a_components;
%model.a_components = @my_a_components2;
%model.a_components_jacobian = @my_a_components2_jacobian;

% for dummy demo_rb_gui
model.gen_detailed_data  = @my_gen_detailed_data;

%%% for real rb simulation:
model.rb_reconstruction = @my_rb_reconstruction;
model.gen_reduced_data  = @my_gen_reduced_data;
model.rb_simulation  = @my_rb_simulation;

% for demo rb gui
model.get_rb_size = @(model,detailed_data) size(detailed_data.RB,2);
model.get_mu = @get_mu_default;
model.is_stationary = 1;
model.reduced_data_subset = @my_reduced_data_subset;
model.get_dofs_from_sim_data = @(sim_data) sim_data.u;
model.gridtype = 'onedgrid';
model.no_lines = 0;
model.axis_equal = 0;

model.RB_Mtrain_size = 1000;
model.Newton_eps = 1e-10;

function model_data = my_gen_model_data(model)
params.xrange = model.xrange;
params.xnumintervals = model.xnumintervals;
model_data.grid = onedgrid(params);
model_data.xnelements = model.xnumintervals-1;
model_data.K = speye(model_data.xnelements) * model.dx;

function sim_data = my_detailed_simulation(model, model_data)

K = model_data.K;

% Newton loop
errsqr = 1e10;
k = 0;
uk = zeros(model_data.xnelements,1);
while errsqr>model.Newton_eps;
  F = model.newton_residual(model, model_data, uk);
  DF = model.newton_residual_jacobian(model, model_data, uk);

  min_hk = DF\F;
  uk = uk - min_hk;
  errsqr = min_hk'*K'*min_hk;
  k = k+1;
end;

sim_data.u = uk;
sim_data.iterations = k;
disp(['finished detailed_simulation with ',num2str(k),' iterations'])


%%%%%%%%%%%%%%%%%%%%%% data functions:

function F = my_newton_residual(model, model_data, uk)
a = model.a_coefficients(model) * model.a_components(model, model_data, uk);
b = model.b_coefficients(model) * model.b_components(model, model_data, uk);
f = model.f_coefficients(model) * model.f_components(model, model_data, uk);

F = (a + b)*uk - f;

function DF = my_newton_residual_jacobian(model, model_data, uk)
a = model.a_coefficients(model) * model.a_components_jacobian(model, model_data, uk);
b = model.b_coefficients(model) * model.b_components_jacobian(model, model_data, uk);

DF = (a + b);

function b_comp = my_b_components(model, model_data, uk)
% matrix b(xi_i,xi_j) of detailed model
% single component consisting of triangular matrix
n = model_data.xnelements;
dx = (model.dx)^(-2);
d = 2* dx * ones(n,1);
d1 = ones(n-1,1)*(-dx);
b_comp = spdiags([ [d1;0], d, [0;d1] ], [-1, 0, 1], n, n);

function b_coeff = my_b_coefficients(model)
% matrix b(xi_i,xi_j) of detailed model
b_coeff = [model.k];

function f_comp = my_f_components(model, model_data, uk)
% right hand side vector f(xi_i)
f_comp = ones(model_data.xnelements,1);

function f_coeff = my_f_coefficients(model)
% right hand side vector f(xi_i)
f_coeff = [1]; % always 1

function a_comp = my_a_components(model, model_data, uk)
% matrix a(xi_i,xi_j) decomposed
a_comp = speye(model_data.xnelements);

function a_coeff = my_a_coefficients(model)
% matrix a(xi_i,xi_j)
a_coeff = [model.m];

function a_comp = my_a_components2(model, model_data, uk)
% matrix a(xi_i,xi_j) decomposed
a_comp = spdiags(uk, 0, model_data.xnelements, model_data.xnelements);

function a_comp = my_a_components2_jacobian(model, model_data, uk)
% matrix a(xi_i,xi_j) decomposed
a_comp = 2*spdiags(uk, 0, model_data.xnelements, model_data.xnelements);


function a_coeff = my_a_coefficients2(model)
% matrix a(xi_i,xi_j)
a_coeff = [model.m];


%%%%%%% other help functions

function p = my_plot_sim_data(model, detailed_data, sim_data, plot_params)
uplot = [zeros(1,size(sim_data.u,2));sim_data.u;zeros(1,size(sim_data.u,2))];
if ~isfield(plot_params, 'color')
  plot_params.color = 'black';
end
p = plot(detailed_data.grid.X,uplot','Linewidth',2,'Color',plot_params.color);
set(gca,'xlim',model.xrange);
if ~isfield(plot_params, 'ylim')
  plot_params.ylim = [min(sim_data.u), max(sim_data.u)];
end
set(gca,'ylim', plot_params.ylim);
if isfield(plot_params,'title')
  title(plot_params.title);
end;
if isfield(plot_params,'legend')
  legend(plot_params.legend);
end;



