rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
fem_matrix_diff_integral_kernel.m
1 function res = fem_matrix_diff_integral_kernel(x,model,df_info,i,j)
2 %function res = fem_matrix_diff_integral_kernel(x,model,df_info,i,j)
3 %
4 % auxiliary function for integral kernel for A_diff
5 % integral kernel on all elements simultaneously
6 % here x is in reference triangle!
7 % ``f(x) = (\nabla \hat\phi_j)^T * JIT^T * A^T * (JIT * \nabla \hat\phi_i)``
8 % multiplication with `|det(DF)|` is realized in caller after quadrature.
9 % function can integrate cell-array valued function
10 
11 % B. Haasdonk 22.2.2011
12 % I. Maier 24.03.2011
13 
14 gradient_hat_phi_i = fem_evaluate_scalar_basis_function_derivative(df_info,x,i)';
15 gradient_hat_phi_j = fem_evaluate_scalar_basis_function_derivative(df_info,x,j)';
16 a = model.diffusivity_tensor(df_info.grid,1:df_info.grid.nelements,x,model); % local model!
17 JIT = df_info.grid.JIT;
18 
19 if ~iscell(a)
20  % the following are vectorial of size nelements!
21  % AT = Atransposed = (a1,a2; a3,a4)
22  ATJ_11 = a(:,1) .* JIT(:,1,1) + a(:,2) .* JIT(:,2,1);
23  ATJ_21 = a(:,3) .* JIT(:,1,1) + a(:,4) .* JIT(:,2,1);
24  ATJ_12 = a(:,1) .* JIT(:,1,2) + a(:,2) .* JIT(:,2,2);
25  ATJ_22 = a(:,3) .* JIT(:,1,2) + a(:,4) .* JIT(:,2,2);
26  %JTATJ_11 = ...
27  % ( a(:,1) .* JIT(:,1,1) + a(:,3) .* JIT(:,2,1)) .* JIT(:,1,1) + ...
28  % ( a(:,1) .* JIT(:,1,2) + a(:,3) .* JIT(:,2,2)) .* JIT(:,2,1);
29  %JTATJ_12 = ...
30  % ( a(:,1) .* JIT(:,1,1) + a(:,3) .* JIT(:,2,1)) .* JIT(:,1,2) + ...
31  % ( a(:,1) .* JIT(:,1,2) + a(:,3) .* JIT(:,2,2)) .* JIT(:,2,2);
32  %JTATJ_21 = ...
33  % ( a(:,2) .* JIT(:,1,1) + a(:,4) .* JIT(:,2,1)) .* JIT(:,1,1) + ...
34  % ( a(:,2) .* JIT(:,1,2) + a(:,4) .* JIT(:,2,2)) .* JIT(:,2,1);
35  %JTATJ_22 = ...
36  %% ( a(:,2) .* JIT(:,1,1) + a(:,4) .* JIT(:,2,1)) .* JIT(:,1,2) + ...
37  % ( a(:,2) .* JIT(:,1,2) + a(:,4) .* JIT(:,2,2)) .* JIT(:,2,2);
38  JTATJ_11 = JIT(:,1,1).* ATJ_11 + JIT(:,2,1).* ATJ_21;
39  JTATJ_12 = JIT(:,1,1).* ATJ_12 + JIT(:,2,1).* ATJ_22;
40  JTATJ_21 = JIT(:,1,2).* ATJ_11 + JIT(:,2,2).* ATJ_21;
41  JTATJ_22 = JIT(:,1,2).* ATJ_12 + JIT(:,2,2).* ATJ_22;
42  res = ...
43  JTATJ_11 * gradient_hat_phi_j(1)* gradient_hat_phi_i(1) + ...
44  JTATJ_12 * gradient_hat_phi_j(1)* gradient_hat_phi_i(2)+ ...
45  JTATJ_21 * gradient_hat_phi_j(2)* gradient_hat_phi_i(1)+ ...
46  JTATJ_22 * gradient_hat_phi_j(2)* gradient_hat_phi_i(2);
47 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48 else % iscell!
49  res = cell(1,length(a));
50  for q= 1:length(a);
51  % the following are vectorial of size nelements!
52  % AT = Atransposed = (a1,a2; a3,a4)
53  ATJ_11 = a{q}(:,1) .* JIT(:,1,1) + a{q}(:,2) .* JIT(:,2,1);
54  ATJ_21 = a{q}(:,3) .* JIT(:,1,1) + a{q}(:,4) .* JIT(:,2,1);
55  ATJ_12 = a{q}(:,1) .* JIT(:,1,2) + a{q}(:,2) .* JIT(:,2,2);
56  ATJ_22 = a{q}(:,3) .* JIT(:,1,2) + a{q}(:,4) .* JIT(:,2,2);
57  %JTATJ_11 = ...
58  % ( a(:,1) .* JIT(:,1,1) + a(:,3) .* JIT(:,2,1)) .* JIT(:,1,1) + ...
59  % ( a(:,1) .* JIT(:,1,2) + a(:,3) .* JIT(:,2,2)) .* JIT(:,2,1);
60  %JTATJ_12 = ...
61  % ( a(:,1) .* JIT(:,1,1) + a(:,3) .* JIT(:,2,1)) .* JIT(:,1,2) + ...
62  % ( a(:,1) .* JIT(:,1,2) + a(:,3) .* JIT(:,2,2)) .* JIT(:,2,2);
63  %JTATJ_21 = ...
64  % ( a(:,2) .* JIT(:,1,1) + a(:,4) .* JIT(:,2,1)) .* JIT(:,1,1) + ...
65  % ( a(:,2) .* JIT(:,1,2) + a(:,4) .* JIT(:,2,2)) .* JIT(:,2,1);
66  %JTATJ_22 = ...
67  %% ( a(:,2) .* JIT(:,1,1) + a(:,4) .* JIT(:,2,1)) .* JIT(:,1,2) + ...
68  % ( a(:,2) .* JIT(:,1,2) + a(:,4) .* JIT(:,2,2)) .* JIT(:,2,2);
69  JTATJ_11 = JIT(:,1,1).* ATJ_11 + JIT(:,2,1).* ATJ_21;
70  JTATJ_12 = JIT(:,1,1).* ATJ_12 + JIT(:,2,1).* ATJ_22;
71  JTATJ_21 = JIT(:,1,2).* ATJ_11 + JIT(:,2,2).* ATJ_21;
72  JTATJ_22 = JIT(:,1,2).* ATJ_12 + JIT(:,2,2).* ATJ_22;
73  res{q} = ...
74  JTATJ_11 * gradient_hat_phi_j(1)* gradient_hat_phi_i(1) + ...
75  JTATJ_12 * gradient_hat_phi_j(1)* gradient_hat_phi_i(2)+ ...
76  JTATJ_21 * gradient_hat_phi_j(2)* gradient_hat_phi_i(1)+ ...
77  JTATJ_22 * gradient_hat_phi_j(2)* gradient_hat_phi_i(2);
78  end;
79 end;
80 
81