rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
init_values_transformed_blobs_richards.m
1 function U0 = init_values_transformed_blobs_richards(glob,params)
2 %function U0 = init_values_transformed_blobs_richards(glob,params)
3 %
4 % function constructing the initial values of the convection diffusion
5 % problem in the specified global points glob and parameters.
6 % It returns an initial data function that is mosly homogeneous with several
7 % (seven) blob like structures of higher concentration that drops exponentially
8 % away from their centres.
9 %
10 % required fields in params
11 % c_init: constant for homogeneous initial data to be returned
12 % blob_height: addend to c_init for the maximum concentration of the blobs
13 % gravity: physical gravity
14 % hill_height: geometry parameter for hill of upper gemeotries' boundary
15 %
16 % in 'coefficient' mode the model_data structure is empty
17 
18 % Martin Drohmann 23.9.2009
19 
20 % glob column check
21 if params.debug
22  if ~isempty(glob) && size(glob,1) < size(glob,2)
23  warning('coordinates in variable glob are given row-wise, but expected them to be column-wise');
24  if params.debug > 2
25  keyboard;
26  end
27  end
28 end
29 decomp_mode = params.decomp_mode;
30 
31 if any(ismember(fieldnames(params), 'blob_height'))
32  h = params.blob_height;
33 else
34  h = 0.15;
35 end
36 
37 if decomp_mode == 2
38  U0 = [ h, h, h, h, h, h, h, params.c_init];
39  U0 = [ U0, params.gravity, params.gravity * params.hill_height ];
40 else
41  X = glob(:,1);
42  Y = glob(:,2);
43  B{1} = ((X(:) - 0.5).^2 +(Y(:)-0.7).^2);
44  B{2} = ((X(:) - 0.8).^2 +(Y(:)-0.9).^2);
45  B{3} = ((X(:) - 0.8).^2 +(Y(:)-0.2).^2);
46  B{4} = ((X(:) - 0.2).^2 +(Y(:)-0.2).^2);
47  B{5} = ((X(:) - 0.2).^2 +(Y(:)-0.9).^2);
48  B{6} = ((X(:) - 0.4).^2 +(Y(:)-0.4).^2);
49  B{7} = ((X(:) - 0.6).^2 +(Y(:)-0.6).^2);
50 
51  I=cell(1,7);
52  % p(mu) = mu*(1-x)
53  % p_mu = spline_select(params);
54  for i=1:7
55  I{i} = (B{i} < 0.05.^2);
56  end
57  if decomp_mode == 0
58  U0 = h .* I{1} + params.c_init * ones(length(X),1);
59  % TODO: replace ztrans with something derivated from myspline
60 
61  ztrans = Y.*(1-X);
62  U0 = U0 + params.gravity * (Y + params.hill_height*ztrans);
63 
64  for i = 2:7
65  U0 = U0 + h .* exp(-B{i}) .* I{i};
66  end
67  elseif decomp_mode == 1
68  ztrans = Y.*(1-X);
69  U0 = cell(1,8);
70  U0{1} = I{1};
71  for i = 2:7
72  U0{i} = I{i};
73  end
74  U0{8} = ones(length(X),1);
75  U0{9} = Y;
76  U0{10} = ztrans;
77  else
78  error(['decomp_mode number ' params.decomp_mode, ' is unknown.']);
79  end
80 end
81 %| \docupdate