rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
comsol_thermal_block_operators_output.m
1 function v = comsol_thermal_block_operators_output(model,model_data)
2 %function v = comsol_thermal_block_operators_output(model,model_data)
3 %
4 % This function returns the output operator v for the Comsol 2D stationary
5 % Thermal Block model.
6 %
7 % The output is defined as the integral over the bottom of the Thermal
8 % block:
9 % s(\mu) = l(u(\mu)) = \int_0^1 u(x,0) dx
10 % and can be comupted s(\mu) = v' * (rb_)sim_data.U;
11 %
12 % This function checks whether linear or quadratic finfite elements are
13 % used and uses a Simpson-rule (pdeg=2) or an trapezoidal rule (pdeg=1)
14 % If other Finite Elements are used, a trapezoidal rule is used (BUT THIS
15 % CASE HAS NOT BEEN TESTED, SO BE CAREFUL!!!)
16 %
17 % Input:
18 % - model
19 % - model_data
20 %
21 % Output:
22 % - v : output_operator: s=v'*sim-data.U;
23 %
24 % Oliver Zeeb, 2013/11/21
25 
26 
27 switch model.decomp_mode
28  case 0
29  error('not yet implemented, since unneeded!')
30  case 1 %components
31  if model.use_comsol
32  grid_info = model_data.grid_info;
33 
34  % get dof indices of the bottom
35  bnd_dof_ind = find(grid_info.dofs.coords(2,:) <= 10*eps);
36 
37  %get the according coords and values
38  bnd_coords = grid_info.dofs.coords(1,bnd_dof_ind);
39 
40  %resort for quadrature
41  [bnd_coord_sort,bnd_coord_sort_ind] = sort(bnd_coords);
42  %bnd_coord_sort is now:
43  %for pdeg=1: x1 x2 x3 x4 ...
44  %for pdeg=2: x1 (x1+x2)/2 x2 (x2+x3)/2 x3 ... --> coords in between of the mesh edges!
45 
46  if strcmp(model.element_type{1},'shlag') %check lagrange finite elements are used
47  v_dof = zeros(1,size(bnd_coords,2));
48  switch model.pdeg{1}
49  case 2 %fuer pdeg = 2: Simpson-Rule
50  v_dof(1) = 1/6*(bnd_coord_sort(3)-bnd_coord_sort(1));
51  v_dof(2:2:end-1) = diff(bnd_coord_sort(1:2:end))*2/3;
52  v_dof(3:2:end-2) = 1/6* (bnd_coord_sort(5:2:end) - bnd_coord_sort(1:2:end-3));
53  v_dof(end) = 1/6*(bnd_coord_sort(end)-bnd_coord_sort(end-2));
54  otherwise %pdeg = 1 or anything else: trapezoidal rule
55  v_dof(1) = 0.5*(bnd_coord_sort(2) - bnd_coord_sort(1));
56  v_dof(2:end-1) = 0.5 * (bnd_coord_sort(3:end) - bnd_coord_sort(1:end-2));
57  v_dof(end) = 0.5 * (bnd_coord_sort(end) - bnd_coord_sort(end-1));
58  end
59  else %no lagrangian elements!
60  warning('using trapezodial rule for output integration! Might be inaccurate since used element type is unknown!')
61  warning('THIS CASE HAS NOT BEEN TESTED EXTENSEVELY SO FAR...')
62  v_dof(1) = 0.5*(bnd_coord_sort(2) - bnd_coord_sort(1));
63  v_dof(2:end-1) = 0.5 * (bnd_coord_sort(3:end) - bnd_coord_sort(1:end-2));
64  v_dof(end) = 0.5 * (bnd_coord_sort(end) - bnd_coord_sort(end-1));
65  end
66 
67  %Blow up the Vector again to be able to calculate v' * sim_data.U
68  v = zeros(grid_info.ndofs,1);
69  bnd_dof_ind_sort = bnd_dof_ind(bnd_coord_sort_ind); %since v is sorted, the indices of the bnd_dofs have to be sorted accordingly
70  %now blow up the vector:
71  v(bnd_dof_ind_sort) = v_dof;
72 
73  %if eliminated data wanted: cut the vector!
74  if model.comsol_get_eliminated_data
75  v=v(model_data.ind_vectors.com_DOF_ind);
76  end
77 
78  v={v};
79 
80  else %~model.use_comsol
81  if model.comsol_get_eliminated_data
82  v = model_data.operators.output_comp_eliminated;
83  else %full data
84  v = model_data.operators.output_comp_full;
85  end
86  end
87 
88  case 2 %coefficients
89  v = 1;
90 end