1 function grad = gradient_approx(model,model_data,U, NU_ind, edge)
2 %
function grad = gradient_approx(model,model_data,U, [NU_ind], edge)
4 %
function that computes an approximation of the gradient over the edge
7 % Martin Drohmann 13.05.2008
9 grid = model_data.grid;
10 nels = grid.nelements;
15 %
get a 6-environment of each cell (see get_enbi
for details)
16 ENBI = get_enbi(grid, edge, model.tstep);
17 % initialise UE, a 4 environment that stores interpolated values at the
18 % edges tips and exact values at the two adjacent cell centers. From UE
19 % the gradient is calculated by lineare interpolation (lsq)
22 %%%% Handle dirichlet cells ("ghost" cells with index -1) %%%%
23 % if a cell adjacent to the current edge is a _dirichlet_ "ghost" cell,
24 % entries in the matrix UE have to be replaced by exact evaluations at
25 % the dirichlet boundary. The following cases have to be separated:
26 % case 1: one of the adjacent cells is a ghost
27 % => replace the cell's entry in UE with the dirichlet value at
29 % case 2/3: one cell adjacent to the tip is a ghost
30 % => replace the tip's entry in UE with the dirichlet value at the
33 % find the _row_ index (edge inedex) of the three cases
34 case1_dir_ind = find(ENBI(:,4) == -1);
35 case2_dir_ind = find(ENBI(:,5) == -1);
36 case3_dir_ind = find(ENBI(:,6) == -1);
38 % the case23_mapper maps to the correct edge tip (correct column index)
41 case23_mapper = [ mod(edge, 4) + 1, edge ];
43 case23_mapper = [ edge, mod(edge, 4) + 1 ];
46 % neuman boundaries!!! For all other ghost cells the entries of
47 % adjacent non-ghost cells are copied. Ghost cells in the corners of
48 % the
rectgrid are assigned in the end by copying one of the already
49 % assigned adjacent ghost-cell. It is unspecified which ghost cell is
51 if(any(any(ENBI(NU_ind, :) == -10)))
52 disp('Attention in gradient_approx! This line should never be executed');
54 nb_ind = find(ENBI(:,2) <= 0);
55 ENBI(nb_ind, 2) = ENBI(nb_ind, 1);
56 nb_ind = find(ENBI(:,3) <= 0);
57 ENBI(nb_ind, 3) = ENBI(nb_ind, 1);
58 nb_ind = find(ENBI(:,4) <= 0);
59 ENBI(nb_ind, 4) = ENBI(nb_ind, 1);
60 nb_ind = find(ENBI(:,5) == -2);
61 ENBI(nb_ind, 5) = ENBI(nb_ind, 2);
62 nb_ind = find(ENBI(:,5) <= 0);
63 ENBI(nb_ind, 5) = ENBI(nb_ind, 4);
64 nb_ind = find(ENBI(:,6) == -2);
65 ENBI(nb_ind, 6) = ENBI(nb_ind, 3);
66 nb_ind = find(ENBI(:,6) <= 0);
67 ENBI(nb_ind, 6) = ENBI(nb_ind, 4);
69 % now every ghost cell should have a positive index
70 % assert(max(max(ENBI <= 0)) == 0);
73 % + ytrans(grid.CX(ENBI(non_zero_ind)), grid.CY(ENBI(non_zero_ind)));
75 % construct the 4-environment UE
76 UE(NU_ind,[1,2]) = UENBI(NU_ind,[1,4]);
77 UE(NU_ind,3) = sum(UENBI(NU_ind,[1,2,4,5]), 2) / 4;
78 UE(NU_ind,4) = sum(UENBI(NU_ind,[1,3,4,6]), 2) / 4;
79 % update the 4-environment concerning dirichlet cells
80 Xdir = grid.ESX(case1_dir_ind, edge);
81 Ydir = grid.ESY(case1_dir_ind, edge);
82 UE(case1_dir_ind, 2) = model.dirichlet_values_ptr([Xdir,Ydir], model);
83 Xdir = grid.X(grid.VI(case2_dir_ind, case23_mapper(1) ));
84 Ydir = grid.Y(grid.VI(case2_dir_ind, case23_mapper(1) ));
85 UE(case2_dir_ind, 3) = model.dirichlet_values_ptr([Xdir,Ydir], model);
86 Xdir = grid.X(grid.VI(case3_dir_ind, case23_mapper(2) ));
87 Ydir = grid.Y(grid.VI(case3_dir_ind, case23_mapper(2) ));
88 UE(case3_dir_ind, 4) = model.dirichlet_values_ptr([Xdir,Ydir], model);
90 horiz_length = grid.EL(1,1);
91 vert_length = grid.EL(1,2);
94 grad = [(UE(NU_ind,2)-UE(NU_ind,1))./horiz_length, (UE(NU_ind,4)-UE(NU_ind,3))./vert_length];
96 grad = [(UE(NU_ind,3)-UE(NU_ind,4))./horiz_length, (UE(NU_ind,2)-UE(NU_ind,1))./vert_length];
98 grad = [(UE(NU_ind,1)-UE(NU_ind,2))./horiz_length, (UE(NU_ind,4)-UE(NU_ind,3))./vert_length];
100 grad = [(UE(NU_ind,3)-UE(NU_ind,4))./horiz_length, (UE(NU_ind,1)-UE(NU_ind,2))./vert_length];
A cartesian rectangular grid in two dimensions with axis parallel elements.