func_base.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of NS2DDV. %
% %
% Copyright(C) 2011-2018 C. Calgaro (caterina.calgaro@math.univ-lille1.fr) %
% E. Creusé (emmanuel.creuse@math.univ-lille1.fr) %
% T. Goudon (thierry.goudon@inria.fr) %
% A. Mouton (alexandre.mouton@math.univ-lille1.fr) %
% %
% NS2DDV is free software: you can redistribute it and/or modify it under the terms %
% of the GNU General Public License as published by the Free Software Foundation, %
% either version 3 of the License, or (at your option) any later version. %
% %
% NS2DDV is distributed in the hope that it will be useful, but WITHOUT ANY %
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A %
% PARTICULAR PURPOSE. See the GNU General Public License for more details. %
% %
% You should have received a copy of the GNU General Public License along with %
% NS2DDV. If not, see <http://www.gnu.org/licenses/>. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[base] = func_base(FE_approx, deriv_ordermax, int_points)
% INPUTS
% FE_approx Finite element method for which we want to compute the reference DOF
% deriv_ordermax Maximal order for the derivatives of the DOF to be computed
% int_points List of points where the DOF and their derivatives will be evaluated
%
% OUTPUT
% base The list of basis functions and their derivatives evaluated at points in int_points
% Allowed values for FE_approx : 'P1', 'P1B', 'P2'
d = (deriv_ordermax+1)*(deriv_ordermax+2)/2;
np = size(int_points,1);
switch FE_approx
case {'P1'}
% 3 DOF on the triangle
base = zeros(np, 3, d);
% Evaluate the DOF
if (deriv_ordermax >= 0)
base(:,1,1) = 1.0 - int_points(:,1) - int_points(:,2);
base(:,2,1) = int_points(:,1);
base(:,3,1) = int_points(:,2);
end
% Evaluate the first order derivatives of the DOF
if (deriv_ordermax >= 1)
% Apply d/dx
base(:,1,2) = -1.0;
base(:,2,2) = 1.0;
base(:,3,2) = 0.0;
% Apply d/dy
base(:,1,3) = -1.0;
base(:,2,3) = 0.0;
base(:,3,3) = 1.0;
end
% Higher order derivatives are equal to 0
if (deriv_ordermax >= 2)
base(:,:,4:d) = 0.;
end
case {'P1B'}
% 4 DOF on the triangle
base = zeros(np, 4, d);
% Evaluate the DOF
if (deriv_ordermax >= 0)
base(:,4,1) = 27.*(1. - int_points(:,1) - int_points(:,2)).*int_points(:,1).*int_points(:,2);
base(:,1,1) = 1. - int_points(:,1) - int_points(:,2) - base(:,4,1)/3.;
base(:,2,1) = int_points(:,1) - base(:,4,1)/3.;
base(:,3,1) = int_points(:,2) - base(:,4,1)/3.;
end
% Evaluate the first order derivatives of the DOF
if (deriv_ordermax >= 1)
% Apply d/dx
base(:,4,2) = -54.*int_points(:,1).*int_points(:,2) + 27.*int_points(:,2) - 27.*int_points(:,2).^2;
base(:,1,2) = -1. - base(:,4,2)/3.;
base(:,2,2) = 1. - base(:,4,2)/3.;
base(:,3,2) = -base(:,4,2)/3.;
% Apply d/dy
base(:,4,3) = -54.*int_points(:,2).*int_points(:,1) + 27.*int_points(:,1) - 27.*int_points(:,1).^2;
base(:,1,3) = -1. - base(:,4,3)/3.;
base(:,2,3) = -base(:,4,3)/3.;
base(:,3,3) = 1. - base(:,4,3)/3.;
end
% Evaluate the second order derivatives of the DOF
if (deriv_ordermax >= 2)
% Apply d^2/dx^2
base(:,4,4) = - 54.*int_points(:,2);
base(:,1,4) = - base(:,4,4)/3.;
base(:,2,4) = - base(:,4,4)/3.;
base(:,3,4) = - base(:,4,4)/3.;
% Apply d^2/dy^2
base(:,4,5) = - 54.*int_points(:,1);
base(:,1,5) = - base(:,4,5)/3.;
base(:,2,5) = - base(:,4,5)/3.;
base(:,3,5) = - base(:,4,5)/3.;
% Apply d^2/dxdy
base(:,4,6) = -54.*int_points(:,1) + 27. - 54.*int_points(:,2);
base(:,1,6) = - base(:,4,6)/3.;
base(:,2,6) = - base(:,4,6)/3.;
base(:,3,6) = - base(:,4,6)/3.;
end
% Evaluate the second order derivatives of the DOF
if (deriv_ordermax >= 3)
% Apply d^3/dx^3
base(:,4,7) = 0.;
base(:,1,7) = -base(:,4,7)/3.;
base(:,2,7) = -base(:,4,7)/3.;
base(:,3,7) = -base(:,4,7)/3.;
% Apply d^3/dy^3
base(:,4,8) = 0.;
base(:,1,8) = -base(:,4,8)/3.;
base(:,2,8) = -base(:,4,8)/3.;
base(:,3,8) = -base(:,4,8)/3.;
% Apply d^3/dx^2dy
base(:,4,9) = -54.;
base(:,1,9) = -base(:,4,9);
base(:,2,9) = -base(:,4,9);
base(:,3,9) = -base(:,4,9);
% Apply d^3/dxdy^2
base(:,4,10) = -54.;
base(:,1,10) = -base(:,4,10);
base(:,2,10) = -base(:,4,10);
base(:,3,10) = -base(:,4,10);
end
% Higher order derivatives are equal to 0
if (deriv_ordermax >= 4)
base(:,:,11:d) = 0.0;
end
case {'P2'}
base = zeros(np, 6, d);
% Evaluate the DOF
if (deriv_ordermax >= 0)
base(:,1,1) = 1. - 3.*int_points(:,1) - 3.*int_points(:,2) + 4.*int_points(:,1).*int_points(:,2) ...
+ 2.*int_points(:,1).^2 + 2.*int_points(:,2).^2;
base(:,2,1) = - int_points(:,1) + 2.*int_points(:,1).^2;
base(:,3,1) = - int_points(:,2) + 2.*int_points(:,2).^2;
base(:,4,1) = 4.*int_points(:,1).*int_points(:,2);
base(:,5,1) = 4.*int_points(:,2) - 4.*int_points(:,1).*int_points(:,2) - 4.*int_points(:,2).^2;
base(:,6,1) = 4.*int_points(:,1) - 4.*int_points(:,1).*int_points(:,2) - 4.*int_points(:,1).^2;
end
% Evaluate the first order derivatives of the DOF
if (deriv_ordermax >= 1)
% Apply d/dx
base(:,1,2) = -3. + 4.*int_points(:,2) + 4.*int_points(:,1);
base(:,2,2) = -1. + 4.*int_points(:,1);
base(:,3,2) = 0.;
base(:,4,2) = 4.*int_points(:,2);
base(:,5,2) = -4.*int_points(:,2);
base(:,6,2) = 4. - 4.*int_points(:,2) - 8.*int_points(:,1);
% Apply d/dy
base(:,1,3) = -3. + 4.*int_points(:,1) + 4.*int_points(:,2);
base(:,2,3) = 0.;
base(:,3,3) = -1. + 4.*int_points(:,2);
base(:,4,3) = 4.*int_points(:,1);
base(:,5,3) = 4. - 4.*int_points(:,1) - 8.*int_points(:,2);
base(:,6,3) = -4.*int_points(:,1);
end
% Evaluate the second order derivatives of the DOF
if (deriv_ordermax >= 2)
% Apply d^2/dx^2
base(:,1,4) = 4.;
base(:,2,4) = 4.;
base(:,3,4) = 0.;
base(:,4,4) = 0.;
base(:,5,4) = 0.;
base(:,6,4) = -8.;
% Apply d^2/dy^2
base(:,1,5) = 4.;
base(:,2,5) = 0.;
base(:,3,5) = 4.;
base(:,4,5) = 0.;
base(:,5,5) = -8.;
base(:,6,5) = 0.;
% Apply d^2/dxdy
base(:,1,6) = 4.;
base(:,2,6) = 0.;
base(:,3,6) = 0.;
base(:,4,6) = 4.;
base(:,5,6) = -4.;
base(:,6,6) = -4.;
end
% Higher order derivatives are equal to 0
if (deriv_ordermax >= 3)
base(:,:,7:d) = 0.0;
end
otherwise
error(sprintf('Degrees of freedom are implemented only for P1, P2, and P1B finite element methods.\n'));
end
end
