Code covered by the BSD License

# FELICITY

### Shawn Walker (view profile)

21 Apr 2011 (Updated )

Finite ELement Implementation and Computational Interface Tool for You

Get_BDM1_Interpolant_of_Function(Vtx,Tri,BDM1_DoFmap,Orient,FUNC)
```function uu_VEC = Get_BDM1_Interpolant_of_Function(Vtx,Tri,BDM1_DoFmap,Orient,FUNC)

Num_Tri = size(Tri,1);
Num_BDM1_DoFs = max(BDM1_DoFmap(:));

if (Num_Tri~=size(BDM1_DoFmap,1))
error('Number of Tri''s must equal number of rows of BDM1_DoFmap!');
end
if (Num_Tri~=size(Orient,1))
error('Number of Tri''s must equal number of rows of Orient!');
end

% init
uu_VEC = zeros(Num_BDM1_DoFs,1);
TAN = zeros(3,2);
NOR = zeros(3,2);

Orient_Dir = zeros(size(Orient,1),3);
Orient_Dir(Orient)  =  1;
Orient_Dir(~Orient) = -1;

% get quad rule on the unit interval

% loop through triangles
for index=1:Num_Tri
BDM1_i = BDM1_DoFmap(index,:);
% get x1,x2,x3 vertices
X1 = Vtx(Tri(index,1),:);
X2 = Vtx(Tri(index,2),:);
X3 = Vtx(Tri(index,3),:);
%     Node1 = X2;
%     Node2 = X3;
%     Node3 = X3;
%     Node4 = X1;
%     Node5 = X1;
%     Node6 = X2;
%     X_Mid_1 = 0.5 * (X2 + X3);
%     X_Mid_2 = 0.5 * (X3 + X1);
%     X_Mid_3 = 0.5 * (X1 + X2);
% compute t's and n's
TAN(1,:) = X3 - X2;
TAN(2,:) = X1 - X3;
TAN(3,:) = X2 - X1;
NOR(1,:) = [TAN(1,2), -TAN(1,1)];
NOR(2,:) = [TAN(2,2), -TAN(2,1)];
NOR(3,:) = [TAN(3,2), -TAN(3,1)];
% do not normalize!  that way the normal already has the edge length factor
% included.

% set DoFs
phi_tail = (1-Quad.Pt); % in local coordinates
uu_VEC(BDM1_i(1),1) = Orient_Dir(index,1) * eval_edge_integral(FUNC,phi_tail,X2,X3,NOR(1,:),Quad); % phi_1

uu_VEC(BDM1_i(3),1) = Orient_Dir(index,2) * eval_edge_integral(FUNC,phi_tail,X3,X1,NOR(2,:),Quad); % phi_3

uu_VEC(BDM1_i(5),1) = Orient_Dir(index,3) * eval_edge_integral(FUNC,phi_tail,X1,X2,NOR(3,:),Quad); % phi_5

%     % this is actually an integral!!! Not ``point-wise'' evaluation
%     % the half factor is due to integrating linear polys (dual basis) on each
%     % edge.
%     uu_VEC(BDM1_i(1),1) = Orient_Dir(index,1) * 0.5 * dot(NOR(1,:), FUNC(X2)); % phi_1
%     uu_VEC(BDM1_i(2),1) = Orient_Dir(index,1) * 0.5 * dot(NOR(1,:), FUNC(X3)); % phi_2
%
%     uu_VEC(BDM1_i(3),1) = Orient_Dir(index,2) * 0.5 * dot(NOR(2,:), FUNC(X3)); % phi_3
%     uu_VEC(BDM1_i(4),1) = Orient_Dir(index,2) * 0.5 * dot(NOR(2,:), FUNC(X1)); % phi_4
%
%     uu_VEC(BDM1_i(5),1) = Orient_Dir(index,3) * 0.5 * dot(NOR(3,:), FUNC(X1)); % phi_5
%     uu_VEC(BDM1_i(6),1) = Orient_Dir(index,3) * 0.5 * dot(NOR(3,:), FUNC(X2)); % phi_6

end

end

% note: N already has the edge length factor (Jacobian) included!