1 function v = comsol_thermal_block_operators_output(model,model_data)
2 %
function v = comsol_thermal_block_operators_output(model,model_data)
4 % This
function returns the output
operator v
for the Comsol 2D stationary
7 % The output is defined as the integral over the bottom of the Thermal
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;
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!!!)
22 % - v : output_operator: s=v'*sim-data.U;
24 % Oliver Zeeb, 2013/11/21
27 switch model.decomp_mode
29 error('not yet implemented, since unneeded!')
32 grid_info = model_data.grid_info;
34 % get dof indices of the bottom
35 bnd_dof_ind = find(grid_info.dofs.coords(2,:) <= 10*eps);
37 %get the according coords and values
38 bnd_coords = grid_info.dofs.coords(1,bnd_dof_ind);
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!
46 if strcmp(model.element_type{1},
'shlag') %check lagrange finite elements are used
47 v_dof = zeros(1,size(bnd_coords,2));
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));
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));
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;
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);
80 else %~model.use_comsol
81 if model.comsol_get_eliminated_data
82 v = model_data.operators.output_comp_eliminated;
84 v = model_data.operators.output_comp_full;