rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
diffusivity_tensor_richards.m
1 function diff = diffusivity_tensor_richards(glob,U,params,callerid)
2 %function diffusivity = diffusivity_tensor_richards(glob,U,params,callerid)
3 %
4 % function computing the diffucivity tensor of a diffusion problem
5 % by pointwise evaluation in the point sequences indicated by X and
6 % Y and solution U.
7 %
8 % The callerid needs to be set in case the diffusivity tensor is cached. It
9 % adds a unique item to the argument list of a later inv_geo_trans_derivative
10 % call. That allows correct hashing.
11 %
12 % generated fields of diffusivity:
13 % epsilon: upper bound on diffusivity value
14 % K : diffusivity tensor, i.e. a sparse square matrix of size
15 % 2*dim x 2*dim
16 %
17 
18 % Bernard Haasdonk 11.4.2006
19 
20 % determine affine_decomposition_mode as integer
21 % glob column check
22 if params.debug
23  if ~isempty(glob) && size(glob,1) < size(glob,2)
24  warning('coordinates in variable glob are given row-wise, but expected them to be column-wise');
25  if params.debug > 2
26  keyboard;
27  end
28  end
29 end
30 decomp_mode = params.decomp_mode;
31 % flag indicating whether the computation respected the decomposition
32 
33 if decomp_mode~=0 && ~isempty(intersect(params.mu_names,'hill_height'))
34  error('affine decomposition not implemented!!');
35 end;
36 
37 diff = [];
38 diff.epsilon = 0;
39 
40 % [res1, res2] = inv_geo_trans_derivative(glob,{(1), (2)},{(1), (2)},params);
41 % row1 = [res1{1}, res1{2}];
42 % row2 = [res2{1}, res2{2}];
43 
44 % vlen = size(row1,1);
45 % temp0 = reshape([ sum(row1 .* row1, 2), sum(row2 .* row2, 2) ]',2*vlen,1);
46 % tempm1 = reshape([ sum(row2 .* row1, 2), zeros(vlen,1) ]', 2*vlen, 1);
47 % tempp1 = reshape([ zeros(vlen,1), sum(row1 .* row2, 2) ]', 2*vlen, 1);
48 % temp1 = [ sum(row1 .* row1, 2), sum(row1 .* row2, 2) ];
49 % temp2 = [ sum(row2 .* row1, 2), sum(row2 .* row2, 2) ];
50 % tempD = spdiags([tempm1,temp0,tempp1],-1:1,2*vlen,2*vlen);
51 tempD = real(diffusivity_cached(glob,params,callerid));
52 
53 % min(min(abs(atemp1 - temp1) < 1e-16))
54 % min(min(abs(atemp2 - temp2) < 1e-14))
55 % keyboard;
56 
57 % d.K1 = repmat(diff_k.K,1,2) .* temp1;
58 % d.K2 = repmat(diff_k.K,1,2) .* temp2;
59 
60 
61 % d.K = diff_k.K * tempD;
62 d.K = tempD;
63 
64 % d.K1 = params.k * temp1;
65 % d.K2 = params.k * temp2;
66 % keyboard
67 d.epsilon = max(max(d.K));
68 
69 diff = d;
70 
71 
72 %| \docupdate
function [ P1res , P2res ] = inv_geo_trans_derivative(model, glob, P1derivates, P2derivates, callerid)
computes entries of a geometry transformation function's inverse transposed jacobian ...