praktikum_ellipt_model.m

function model = praktikum_ellipt_model(params);
%function model = praktikum_ellipt_model(params);
%
% small example of a model, i.e. a structure describing the data
% functions and geometry information of a general elliptic equation consisting 
% of diffusion, convection, reaction equation:
%
% - div ( a(x) grad u(x)) + div (b(x) u(x)) + c(x) u(x) = f(x)    on Omega
%                                                 u(x)) = g_D(x)  on Gamma_D
%                                       a(x) grad u(x)) = g_N(x)  on Gamma_N
% 
%  Here, we denote the functions as
%                   u: solution (if known, useful for validation purpose)
%                   f: source
%                   a: diffusivity
%                   b: velocity
%                   c: reaction
%                 g_D: Dirichlet boundary values
%                 g_N: Neumann boundary values
%
% Each function allows the evaluation in many points
% simultaneuously by
%
%        model.source(glob)
% or     model.source(glob,params)
%
% where glob is a n times 2 matrix of row-wise points. The result
% is a n times 1 vector of resulting values of f.
%
% additionally, the model has a function, which determines, whether
% a point lies on a Dirichlet or Neumann boundary:
%        
%           model.boundary_type(glob)
%                0 no boundary (inner edge or point)
%               -1 indicates Dirichlet-boundary
%               -2 indicates Neumann-boundary
%
% The data functions given in this model are given on the exercise sheet
% Omega = [-1,1]^2 \ [0,1]^2
% with Gamma_D = ecke oben rechts, Gamma_N = rand vom ursprünglichen 4eck
%
%    -div ( a div(u)) = f
%                   u = g_D   on   Gamma_D
%        (a div(u))*n = g_N  on Gamma_N
%
% with exact solution u(x) = ||x||^(alpha), 
%
% params is an optional parameter, perhaps useful later

% B. Haasdonk 23.11.2010

if nargin<1
  params = [];
end;

% check if user wants to have pure dirichlet boundary:
if ~isfield(params,'all_dirichlet_boundary')
  params.all_dirichlet_boundary = 1;
end;


%n=10;
%m=20;
%k=3*n+6*m;

%nor=zeros(k,2);
%for i=(3*n+1):(3*n+m)
%    nor(i,2)=1;
%end
%for i=(3*n+m+1):(3*n+2*m)
%    nor(i,1)=-1;
%end
%for i=(3*n+2*m+1):(3*n+3*m)
%    nor(i,2)=-1;
%end
%for i=(3*n+3*m+1):(3*n+4*m)
%    nor(i,1)=1;
%end
%for i=(3*n+4*m+1):(3*n+5*m)
%    nor(i,2)=1;
%end
%for i=(3*n+5*m+1):k
%    nor(i,1)=1;
%end

model = [];
alpha=0.5;
% params is an optional parameter, perhaps useful later
model.source = @(glob,params) ...
           -alpha*sqrt(glob(:,1).^2+glob(:,2).^2).^(alpha-2).*(3*glob(:,1)+4)-alpha*(alpha-2)*(glob(:,1)+2).*sqrt(glob(:,1).^2+glob(:,2).^2).^(alpha-3);
model.reaction = @(glob,params) zeros(size(glob,1),1);
model.velocity = @(glob,params) zeros(size(glob,1),1);
model.diffusivity = @(glob,params) glob(:,1)+2;
%model.diffusivity_gradient = @(glob,params) [1,0];

%  Dirichlet everywhere, see function below.
if params.all_dirichlet_boundary
  model.boundary_type = @all_dirichlet_boundary_type;
else
  model.boundary_type = @mixed_boundary_type;
end;
model.dirichlet_values = @(glob,params) sqrt(glob(:,1).^2+glob(:,2).^2).^alpha;
model.normals = @my_normals;
model.neumann_values = @(glob,params) alpha*sqrt(glob(:,1).^2+glob(:,2).^2).^(alpha-2).*(glob(:,1)+2).*sum(glob.*model.normals(glob),2);

% solution is known:
model.solution = @(glob,params) ...
             sqrt(glob(:,1).^2+glob(:,2).^2).^alpha;

% all edges of unit square are dirichlet, other inner
function res = all_dirichlet_boundary_type(glob,params)
res = zeros(size(glob,1),1);
i  = find(glob(:,1)<-1+1e-10);
res(i)= -1.0;
i  = find(glob(:,1)>1-1e-10);
res(i)= -1.0;
i  = find(glob(:,2)>1-1e-10);
res(i)= -1.0;
i  = find(glob(:,2)<-1+1e-10);
res(i)= -1.0;
i  = find((glob(:,1)>0-1e-10) & (glob(:,2)>1e-10));
res(i)= -1.0;
i  = find((glob(:,1)>1e-10) & (glob(:,2)>-1e-10));
res(i)= -1.0;

function res = mixed_boundary_type(glob,params)
res = zeros(size(glob,1),1);
i  = find(glob(:,1)<-1+1e-10);
res(i)= -2.0;
i  = find(glob(:,1)>1-1e-10);
res(i)= -2.0;
i  = find(glob(:,2)>1-1e-10);
res(i)= -2.0;
i  = find(glob(:,2)<-1+1e-10);
res(i)= -2.0;
i  = find((glob(:,1)>0-1e-10) & (glob(:,2)>1e-10));
res(i)= -1.0;
i  = find((glob(:,1)>1e-10) & (glob(:,2)>-1e-10));
res(i)= -1.0;

function res = my_normals(glob,params);
res = zeros(size(glob,1),2); % each row one normal
i  = find((glob(:,1)>0-1e-10) & (glob(:,2)>1e-10));
res(i,1)= 1.0;
i  = find((glob(:,1)>1e-10) & (glob(:,2)>-1e-10));
res(i,2)= 1.0;
i  = find(glob(:,1)<-1+1e-10);
res(i,1)= -1.0;
i  = find(glob(:,1)>1-1e-10);
res(i,1)= 1.0;
i  = find(glob(:,2)>1-1e-10);
res(i,2)= 1.0;
i  = find(glob(:,2)<-1+1e-10);
res(i,2)= -1.0;