fem_h10_boundary_inner_product_matrix.m

Go to the documentation of this file.

function G = fem_h10_boundary_inner_product_matrix(df_info, ...
                                                  elind,edgeind)
% function G =
% fem_h10_boundary_inner_product_matrix(df_info,elind,edgeind)
%
% function computing the matrix of the h10 inner product on the
% boundary specified through elind and edgeind. elind and edgeind
% have to be vectors of the same length.
%
% Return values:
%  G: ndofs x ndofs matrix with nonzero entries only for those
%  dofs, who lie on the specified boundary.

% I. Maier, 19.07.2011

qdeg = 2*df_info.pdeg;
ndofs = df_info.ndofs;
nlagrange_nodes = df_info.nlagrange_nodes;

j = sub2ind(size(df_info.grid.EL),elind,edgeind);
EL = df_info.grid.EL(j);

G  = spalloc(ndofs,ndofs,1);
for i = 1:nlagrange_nodes
  for j = 1:nlagrange_nodes    
    
    f = @(x) G_int_kernel(x,df_info,i,j,elind,edgeind);
    res = intervalquadrature(qdeg,f);
    res = res .* EL;

    gids_i = df_info.global_dof_index(elind,i);
    gids_j = df_info.global_dof_index(elind,j);

    G_tmp = sparse(gids_i,gids_j,res,ndofs,ndofs);
    G = G + G_tmp;
  end;
end;

function res = G_int_kernel(x,df_info,i,j,elind,edgeind)
% function res = G_int_kernel(x,df_info,i,j,elind,edgeind)

res = zeros(0,1);
for local_edge_id = 1:3
  lcoord = llocal2local(df_info.grid,local_edge_id,x);
  hat_phi_i = fem_evaluate_scalar_basis_function_derivative(df_info,lcoord,i);
  hat_phi_j = fem_evaluate_scalar_basis_function_derivative(df_info,lcoord,j);
  inds = find(edgeind == local_edge_id);
  switch (local_edge_id)
   case 1
    res = [res; ones(size(inds)) * hat_phi_i(1) * hat_phi_j(1)];
   case 3
    res = [res; ones(size(inds)) * hat_phi_i(2) * hat_phi_j(2)];
   case 2
    res = [res; ones(size(inds)) * (hat_phi_i(2)-hat_phi_i(1)) * ...
           (hat_phi_j(2) - hat_phi_j(1))];
   otherwise
  end;
end;