# Hermite main-interpolation functions for two-dimensional surface (update:24-07-07)

### Ali OZGUL (view profile)

24 Jul 2007 (Updated )

Hermiten

## Contents

% function [N_shape,N_shape2,N_shape3] = Hermiten(total_x,total_y,polynom_node_order,flag)
%Description:
%   Under bending or twisting effect a thin-plate element conforming
%   the shape function's derivative as basis displacement functions on
%   Pascal Triangle. This solution technique is non-conforming of
%   high-order-degree to plate element's shape functions and plate's
%   bending and twisting motion.
%
%  Syntax:
%      [N_shape,N_shape2,N_shape3]=Hermiten(totalx,totaly,nodevalue)
%
%       N_shape : Hermite base interpolation functions
%       N_shape2: n-order degree quadratic-surface shape function (bending motion)
%       N_shape3: n-order degree quadratic-surface shape function (bending+twisting motion)
%
%  Syntax.inputs
%      totalx   : two-dimensional interpolation area's "x" axis length
%      totaly   : two-dimensional interpolation area's "y" axis length
%      nodevalue: hermite basis functions nodevalues (interpolation node-degree)
%
%  Syntax.outputs
%      Nshape : n. order-degree Hermite basis interpolation functions
%      Shapefunctions :two-dimensional area plate shape functions for bending
%      motion
%
%
%  Example:
%  Example.inputs                        (flag==true)
%                 [A,B,C]=hermiten(1,1,3,1) <---|
%  Example.outputs
%
%A =
%
% 1-92*x^2+528*x^3-1088*x^4+768*x^5
%     x-12*x^2+52*x^3-96*x^4+64*x^5
%            64*x^2-256*x^3+256*x^4
%   -16*x^2+128*x^3-320*x^4+256*x^5
%    28*x^2-272*x^3+832*x^4-768*x^5
%       -2*x^2+20*x^3-64*x^4+64*x^5
%
%B =
%                     Node(1) ---->Node(n)
% w    [ (1-23*x^2+66*x^3-68*x^4+24*x^5)*(1-23*y^2+66*y^3-68*y^4+24*y^5),   (x-6*x^2+13*x^3-12*x^4+4*x^5)*(1-23*y^2+66*y^3-68*y^4+24*y^5),          (16*x^2-32*x^3+16*x^4)*(1-23*y^2+66*y^3-68*y^4+24*y^5),   (1-23*x^2+66*x^3-68*x^4+24*x^5)*(y-6*y^2+13*y^3-12*y^4+4*y^5),     (x-6*x^2+13*x^3-12*x^4+4*x^5)*(y-6*y^2+13*y^3-12*y^4+4*y^5),            (16*x^2-32*x^3+16*x^4)*(y-6*y^2+13*y^3-12*y^4+4*y^5),          (1-23*x^2+66*x^3-68*x^4+24*x^5)*(16*y^2-32*y^3+16*y^4),            (x-6*x^2+13*x^3-12*x^4+4*x^5)*(16*y^2-32*y^3+16*y^4),                   (16*x^2-32*x^3+16*x^4)*(16*y^2-32*y^3+16*y^4)]
% Qxx  [ (-46*x+198*x^2-272*x^3+120*x^4)*(1-23*y^2+66*y^3-68*y^4+24*y^5),   (1-12*x+39*x^2-48*x^3+20*x^4)*(1-23*y^2+66*y^3-68*y^4+24*y^5),            (32*x-96*x^2+64*x^3)*(1-23*y^2+66*y^3-68*y^4+24*y^5),   (-46*x+198*x^2-272*x^3+120*x^4)*(y-6*y^2+13*y^3-12*y^4+4*y^5),     (1-12*x+39*x^2-48*x^3+20*x^4)*(y-6*y^2+13*y^3-12*y^4+4*y^5),              (32*x-96*x^2+64*x^3)*(y-6*y^2+13*y^3-12*y^4+4*y^5),          (-46*x+198*x^2-272*x^3+120*x^4)*(16*y^2-32*y^3+16*y^4),            (1-12*x+39*x^2-48*x^3+20*x^4)*(16*y^2-32*y^3+16*y^4),                     (32*x-96*x^2+64*x^3)*(16*y^2-32*y^3+16*y^4)]
% Qyy  [ (1-23*x^2+66*x^3-68*x^4+24*x^5)*(-46*y+198*y^2-272*y^3+120*y^4),   (x-6*x^2+13*x^3-12*x^4+4*x^5)*(-46*y+198*y^2-272*y^3+120*y^4),          (16*x^2-32*x^3+16*x^4)*(-46*y+198*y^2-272*y^3+120*y^4),   (1-23*x^2+66*x^3-68*x^4+24*x^5)*(1-12*y+39*y^2-48*y^3+20*y^4),     (x-6*x^2+13*x^3-12*x^4+4*x^5)*(1-12*y+39*y^2-48*y^3+20*y^4),            (16*x^2-32*x^3+16*x^4)*(1-12*y+39*y^2-48*y^3+20*y^4),            (1-23*x^2+66*x^3-68*x^4+24*x^5)*(32*y-96*y^2+64*y^3),              (x-6*x^2+13*x^3-12*x^4+4*x^5)*(32*y-96*y^2+64*y^3),                     (16*x^2-32*x^3+16*x^4)*(32*y-96*y^2+64*y^3)]

%C =
%                     Node(1)------->Node(n)
% w    [ (1-23*x^2+66*x^3-68*x^4+24*x^5)*(1-23*y^2+66*y^3-68*y^4+24*y^5),   (x-6*x^2+13*x^3-12*x^4+4*x^5)*(1-23*y^2+66*y^3-68*y^4+24*y^5),          (16*x^2-32*x^3+16*x^4)*(1-23*y^2+66*y^3-68*y^4+24*y^5),   (1-23*x^2+66*x^3-68*x^4+24*x^5)*(y-6*y^2+13*y^3-12*y^4+4*y^5),     (x-6*x^2+13*x^3-12*x^4+4*x^5)*(y-6*y^2+13*y^3-12*y^4+4*y^5),            (16*x^2-32*x^3+16*x^4)*(y-6*y^2+13*y^3-12*y^4+4*y^5),          (1-23*x^2+66*x^3-68*x^4+24*x^5)*(16*y^2-32*y^3+16*y^4),            (x-6*x^2+13*x^3-12*x^4+4*x^5)*(16*y^2-32*y^3+16*y^4),                   (16*x^2-32*x^3+16*x^4)*(16*y^2-32*y^3+16*y^4)]
% Qxx  [ (-46*x+198*x^2-272*x^3+120*x^4)*(1-23*y^2+66*y^3-68*y^4+24*y^5),   (1-12*x+39*x^2-48*x^3+20*x^4)*(1-23*y^2+66*y^3-68*y^4+24*y^5),            (32*x-96*x^2+64*x^3)*(1-23*y^2+66*y^3-68*y^4+24*y^5),   (-46*x+198*x^2-272*x^3+120*x^4)*(y-6*y^2+13*y^3-12*y^4+4*y^5),     (1-12*x+39*x^2-48*x^3+20*x^4)*(y-6*y^2+13*y^3-12*y^4+4*y^5),              (32*x-96*x^2+64*x^3)*(y-6*y^2+13*y^3-12*y^4+4*y^5),          (-46*x+198*x^2-272*x^3+120*x^4)*(16*y^2-32*y^3+16*y^4),            (1-12*x+39*x^2-48*x^3+20*x^4)*(16*y^2-32*y^3+16*y^4),                     (32*x-96*x^2+64*x^3)*(16*y^2-32*y^3+16*y^4)]
% Qyy  [ (1-23*x^2+66*x^3-68*x^4+24*x^5)*(-46*y+198*y^2-272*y^3+120*y^4),   (x-6*x^2+13*x^3-12*x^4+4*x^5)*(-46*y+198*y^2-272*y^3+120*y^4),          (16*x^2-32*x^3+16*x^4)*(-46*y+198*y^2-272*y^3+120*y^4),   (1-23*x^2+66*x^3-68*x^4+24*x^5)*(1-12*y+39*y^2-48*y^3+20*y^4),     (x-6*x^2+13*x^3-12*x^4+4*x^5)*(1-12*y+39*y^2-48*y^3+20*y^4),            (16*x^2-32*x^3+16*x^4)*(1-12*y+39*y^2-48*y^3+20*y^4),            (1-23*x^2+66*x^3-68*x^4+24*x^5)*(32*y-96*y^2+64*y^3),              (x-6*x^2+13*x^3-12*x^4+4*x^5)*(32*y-96*y^2+64*y^3),                     (16*x^2-32*x^3+16*x^4)*(32*y-96*y^2+64*y^3)]
% Qxy  [ (-46*x+198*x^2-272*x^3+120*x^4)*(-46*y+198*y^2-272*y^3+120*y^4),   (1-12*x+39*x^2-48*x^3+20*x^4)*(-46*y+198*y^2-272*y^3+120*y^4),            (32*x-96*x^2+64*x^3)*(-46*y+198*y^2-272*y^3+120*y^4),   (-46*x+198*x^2-272*x^3+120*x^4)*(1-12*y+39*y^2-48*y^3+20*y^4),     (1-12*x+39*x^2-48*x^3+20*x^4)*(1-12*y+39*y^2-48*y^3+20*y^4),              (32*x-96*x^2+64*x^3)*(1-12*y+39*y^2-48*y^3+20*y^4),            (-46*x+198*x^2-272*x^3+120*x^4)*(32*y-96*y^2+64*y^3),              (1-12*x+39*x^2-48*x^3+20*x^4)*(32*y-96*y^2+64*y^3),                       (32*x-96*x^2+64*x^3)*(32*y-96*y^2+64*y^3)]

%  Example.inputs
%
%                  [A,B,C]=hermiten(1,1,3,0) <---|  or
%                  [A,B,C]=hermiten(1,1,3)   <---|
%  Example.outputs
%
%A =
%
% 1-23*x^2/a^2+66*x^3/a^3-68*x^4/a^4+24*x^5/a^5
%     x-6*x^2/a+13*x^3/a^2-12*x^4/a^3+4*x^5/a^4
%              16*x^2/a^2-32*x^3/a^3+16*x^4/a^4
%     -8*x^2/a+32*x^3/a^2-40*x^4/a^3+16*x^5/a^4
%    7*x^2/a^2-34*x^3/a^3+52*x^4/a^4-24*x^5/a^5
%          -x^2/a+5*x^3/a^2-8*x^4/a^3+4*x^5/a^4

%B =
%
%[ (1-23*x^2+66*x^3-68*x^4+24*x^5)*(1-23*y^2+66*y^3-68*y^4+24*y^5),   (x-6*x^2+13*x^3-12*x^4+4*x^5)*(1-23*y^2+66*y^3-68*y^4+24*y^5),          (16*x^2-32*x^3+16*x^4)*(1-23*y^2+66*y^3-68*y^4+24*y^5),   (1-23*x^2+66*x^3-68*x^4+24*x^5)*(y-6*y^2+13*y^3-12*y^4+4*y^5),     (x-6*x^2+13*x^3-12*x^4+4*x^5)*(y-6*y^2+13*y^3-12*y^4+4*y^5),            (16*x^2-32*x^3+16*x^4)*(y-6*y^2+13*y^3-12*y^4+4*y^5),          (1-23*x^2+66*x^3-68*x^4+24*x^5)*(16*y^2-32*y^3+16*y^4),            (x-6*x^2+13*x^3-12*x^4+4*x^5)*(16*y^2-32*y^3+16*y^4),                   (16*x^2-32*x^3+16*x^4)*(16*y^2-32*y^3+16*y^4)]
%[ (-46*x+198*x^2-272*x^3+120*x^4)*(1-23*y^2+66*y^3-68*y^4+24*y^5),   (1-12*x+39*x^2-48*x^3+20*x^4)*(1-23*y^2+66*y^3-68*y^4+24*y^5),            (32*x-96*x^2+64*x^3)*(1-23*y^2+66*y^3-68*y^4+24*y^5),   (-46*x+198*x^2-272*x^3+120*x^4)*(y-6*y^2+13*y^3-12*y^4+4*y^5),     (1-12*x+39*x^2-48*x^3+20*x^4)*(y-6*y^2+13*y^3-12*y^4+4*y^5),              (32*x-96*x^2+64*x^3)*(y-6*y^2+13*y^3-12*y^4+4*y^5),          (-46*x+198*x^2-272*x^3+120*x^4)*(16*y^2-32*y^3+16*y^4),            (1-12*x+39*x^2-48*x^3+20*x^4)*(16*y^2-32*y^3+16*y^4),                     (32*x-96*x^2+64*x^3)*(16*y^2-32*y^3+16*y^4)]
%[ (1-23*x^2+66*x^3-68*x^4+24*x^5)*(-46*y+198*y^2-272*y^3+120*y^4),   (x-6*x^2+13*x^3-12*x^4+4*x^5)*(-46*y+198*y^2-272*y^3+120*y^4),          (16*x^2-32*x^3+16*x^4)*(-46*y+198*y^2-272*y^3+120*y^4),   (1-23*x^2+66*x^3-68*x^4+24*x^5)*(1-12*y+39*y^2-48*y^3+20*y^4),     (x-6*x^2+13*x^3-12*x^4+4*x^5)*(1-12*y+39*y^2-48*y^3+20*y^4),            (16*x^2-32*x^3+16*x^4)*(1-12*y+39*y^2-48*y^3+20*y^4),            (1-23*x^2+66*x^3-68*x^4+24*x^5)*(32*y-96*y^2+64*y^3),              (x-6*x^2+13*x^3-12*x^4+4*x^5)*(32*y-96*y^2+64*y^3),                     (16*x^2-32*x^3+16*x^4)*(32*y-96*y^2+64*y^3)]

%C =
%
%[                 (1-23*x^2+66*x^3-68*x^4+24*x^5)*(1-23*y^2+66*y^3-68*y^4+24*y^5),                   (x-6*x^2+13*x^3-12*x^4+4*x^5)*(1-23*y^2+66*y^3-68*y^4+24*y^5),                          (16*x^2-32*x^3+16*x^4)*(1-23*y^2+66*y^3-68*y^4+24*y^5),                   (1-23*x^2+66*x^3-68*x^4+24*x^5)*(y-6*y^2+13*y^3-12*y^4+4*y^5),                     (x-6*x^2+13*x^3-12*x^4+4*x^5)*(y-6*y^2+13*y^3-12*y^4+4*y^5),                            (16*x^2-32*x^3+16*x^4)*(y-6*y^2+13*y^3-12*y^4+4*y^5),                          (1-23*x^2+66*x^3-68*x^4+24*x^5)*(16*y^2-32*y^3+16*y^4),                            (x-6*x^2+13*x^3-12*x^4+4*x^5)*(16*y^2-32*y^3+16*y^4),                                   (16*x^2-32*x^3+16*x^4)*(16*y^2-32*y^3+16*y^4)]
%[                 (-46*x+198*x^2-272*x^3+120*x^4)*(1-23*y^2+66*y^3-68*y^4+24*y^5),                   (1-12*x+39*x^2-48*x^3+20*x^4)*(1-23*y^2+66*y^3-68*y^4+24*y^5),                            (32*x-96*x^2+64*x^3)*(1-23*y^2+66*y^3-68*y^4+24*y^5),                   (-46*x+198*x^2-272*x^3+120*x^4)*(y-6*y^2+13*y^3-12*y^4+4*y^5),                     (1-12*x+39*x^2-48*x^3+20*x^4)*(y-6*y^2+13*y^3-12*y^4+4*y^5),                              (32*x-96*x^2+64*x^3)*(y-6*y^2+13*y^3-12*y^4+4*y^5),                          (-46*x+198*x^2-272*x^3+120*x^4)*(16*y^2-32*y^3+16*y^4),                            (1-12*x+39*x^2-48*x^3+20*x^4)*(16*y^2-32*y^3+16*y^4),                                     (32*x-96*x^2+64*x^3)*(16*y^2-32*y^3+16*y^4)]
%[                 (1-23*x^2+66*x^3-68*x^4+24*x^5)*(-46*y+198*y^2-272*y^3+120*y^4),                   (x-6*x^2+13*x^3-12*x^4+4*x^5)*(-46*y+198*y^2-272*y^3+120*y^4),                          (16*x^2-32*x^3+16*x^4)*(-46*y+198*y^2-272*y^3+120*y^4),                   (1-23*x^2+66*x^3-68*x^4+24*x^5)*(1-12*y+39*y^2-48*y^3+20*y^4),                     (x-6*x^2+13*x^3-12*x^4+4*x^5)*(1-12*y+39*y^2-48*y^3+20*y^4),                            (16*x^2-32*x^3+16*x^4)*(1-12*y+39*y^2-48*y^3+20*y^4),                            (1-23*x^2+66*x^3-68*x^4+24*x^5)*(32*y-96*y^2+64*y^3),                              (x-6*x^2+13*x^3-12*x^4+4*x^5)*(32*y-96*y^2+64*y^3),                                     (16*x^2-32*x^3+16*x^4)*(32*y-96*y^2+64*y^3)]
%[ (-46*x/a^2+198*x^2/a^3-272*x^3/a^4+120*x^4/a^5)*(-46*y+198*y^2-272*y^3+120*y^4),     (1-12*x/a+39*x^2/a^2-48*x^3/a^3+20*x^4/a^4)*(-46*y+198*y^2-272*y^3+120*y^4),                (32*x/a^2-96*x^2/a^3+64*x^3/a^4)*(-46*y+198*y^2-272*y^3+120*y^4),   (-46*x/a^2+198*x^2/a^3-272*x^3/a^4+120*x^4/a^5)*(1-12*y+39*y^2-48*y^3+20*y^4),       (1-12*x/a+39*x^2/a^2-48*x^3/a^3+20*x^4/a^4)*(1-12*y+39*y^2-48*y^3+20*y^4),                  (32*x/a^2-96*x^2/a^3+64*x^3/a^4)*(1-12*y+39*y^2-48*y^3+20*y^4),            (-46*x/a^2+198*x^2/a^3-272*x^3/a^4+120*x^4/a^5)*(32*y-96*y^2+64*y^3),                (1-12*x/a+39*x^2/a^2-48*x^3/a^3+20*x^4/a^4)*(32*y-96*y^2+64*y^3),                           (32*x/a^2-96*x^2/a^3+64*x^3/a^4)*(32*y-96*y^2+64*y^3)]

clc

clear
total_x = 1
total_y = 1
polynom_node_order = 2
flag=1

%Function error control
if (nargin < 3) || (total_x*total_y) == 0
error('sub-function runing error for inputs variables')
end

%flag true or false
if nargin==3
flag=false;
end

%Graphic drawing sensibility
%try ,plot_sensibility = 5 and 15;
plot_sensibility = 8;

%Total degree of freedom
syms dispfunc x y a b real
polynom_degree = 2*polynom_node_order; %[w, Qyy]

%Hermite interpolation's boundary order-degrees
if polynom_node_order < 2 || polynom_node_order > 9
error('This hermite polynomal system is non-conform');
end

total_x =

1

total_y =

1

polynom_node_order =

2

flag =

1



## Produce the Hermite initiative function

    for polynom_order_degree = 1:polynom_degree
if polynom_order_degree == 1;
disp_func = x^(polynom_order_degree-1); %order=1 for first term
else
disp_func = [disp_func  x^(polynom_order_degree-1)];
end
end

step_length_x = 0 : (a)/(polynom_node_order-1) : a ;
%     (1)node---------(2)node--------(3)node-------n(node)
%      |               |              |            |
%      <--steplengthx-->

% Boundary Conditions
for i=1: polynom_node_order;
A(2*(i-1)+1,:) = subs(disp_func , {x},{step_length_x(i)}) ; %   w-freedom
A(2*(i-1)+2,:) = subs(diff(disp_func,x,1) ...
,{x},{step_length_x(i)}) ; % Qyy-freedom
% Rotation initiative function
%    Qyy = diff(disp_func,x,1);
end

% Produce Hermite interpolation functions for interpolation-order-degree
if rank(A)==size(A,1)
N_shape=(disp_func*A^-1)'; %Conforming-shape functions
else
error('This solution not-possible for selected boundary conditions')
end

%%Creating new figure on visible-figure
%    if isempty(get(0,'CurrentFigure'))
%       figure(1)
%    else
%       figure(get(0,'CurrentFigure')+1)
%    end

%if sub-program running secondly than closed hold on
%    clf reset;


## Ploting Hermite Main-interpolation functions

hold on
for i=1: size(N_shape,1)

%       Note:this application not running directly
%       plotdata = subs(subs(Nshape(i),{a,x},{totalx, 0:totalx/(8*polynom_node_order):totalx});
plotdata = subs( ...
subs(N_shape(i),{a},{total_x}) , ...
{x},{0:total_x/(plot_sensibility*polynom_node_order):total_x});

line([i*(total_x)/(2*polynom_node_order)       ...
i*(total_x)/(2*polynom_node_order)],     ...
[min(plotdata) max(plotdata)],'LineStyle','-.')    ;

text ( i*(total_x)/(2*polynom_node_order), ...
max(plotdata), ...
['N(' num2str(i) ')max'], 'FontSize',7,'FontWeight','bold' );

text ( i*(total_x)/(2*polynom_node_order), ...
min(plotdata), ...
['N(' num2str(i) ')min'], 'FontSize',7,'FontWeight','bold' );

%Graphical (x)-and-(y) values
positionx = 0:total_x/(plot_sensibility*polynom_node_order):total_x;
positiony = plotdata;

for j=1:(size(plotdata,2)-2)
if polynom_node_order < 3
text ( positionx(j),     ...
positiony(j)+0.02,...
['N(' num2str(i) ')= ' num2str(plotdata(j))],'FontSize',6);
else
text (positionx(j),      ...
positiony(j)+0.02, ...
['N' num2str(i) ] ,'FontSize',5);
end

end%-----(j)---|

plot( 0:total_x/(plot_sensibility*polynom_node_order):total_x, ...
plotdata, '-ro',                       ...
'LineWidth',1.4,                       ...
'Color',[rand rand rand],              ...
'MarkerEdgeColor','k',                 ...
'MarkerFaceColor',[rand rand rand],    ...
'MarkerSize',5);

end %--------(i)--|

% Ploting texts
xlabel(['Hermite interpolation degree =', num2str(polynom_node_order)]);
ylabel(' Shape functions N(i) ');
title('Shape functions (\psi) ');


## Open dimension 2d-surface hermite polynoms

    Nw   = [];  %vertical displacement freedom.
NQyy = [];  %rotation freedom on y-axis (right hand rule).
NQxx = [];  %rotation freedom on x-axis (right hand rule).
NQxy = [];  %twisting freedom on xy surface.

%N_shapex = N_shape;
%    N_shapey = subs ( ...
%                    subs(N_shape,{x,a},{y,b}) , {b},{totaly});
N_shapey = subs(N_shape,{x,a},{y,total_y});
N_shapex = subs(N_shape,{a},{total_x});

% Produce two-dimensional Hermite interpolation functions
for j=1:polynom_node_order;      % y value;
for i=1:polynom_node_order;  % x value;

Nw   = [Nw    N_shapex(i)*N_shapey(j) ];                    %Nx(i)*Ny(j)         , force effect   - z axis
NQyy = [NQyy  diff(N_shapex(i),x,1)*N_shapey(j) ];          %Ny(j)*[d/dx].Nx(i)  , bending moment - yy axis
NQxx = [NQxx  N_shapex(i)*diff(N_shapey(j),y,1) ];          %Nx(i)*[d/dy].Ny(j)  , bending moment - xx axis
NQxy = [NQxy  diff(N_shape(i),x,1)*diff(N_shapey(j),y,1) ]; %[d/dx].Nx(i)*[d/dy].Ny(j)  , twisting moment - xy surface

end
end

N_shape2=[ Nw
NQyy
NQxx  ];

N_shape3=[ N_shape2
NQxy];


## Hermite interpolation Function on 2D-Surface

%for i=1:1 %polynom_node_order^2 -->all function's graphical interface
for i=1:1
%create new figure-per unit node with three degree of freedom
figure(get(0,'CurrentFigure')+1)

%subplot[3,2,1]
subplot(3,2,1);
plot_fh =vectorize(subs(Nw(i),{a,b},{total_x,total_y}));
ezmesh(plot_fh,[0 total_x 0 total_y],20);colormap([0 0 1]);
if polynom_node_order < 5
title(plot_fh,'FontSize',8);
end
%Ezxontour[3,2,2]
subplot(3,2,2); ezcontour(plot_fh,[0 total_x 0 total_y]);
title('Ezcontour','FontSize',8);

%subplot[3,2,3]
subplot(3,2,3);
plot_fh =vectorize(subs(NQyy(i),{a,b},{total_x,total_y}));
ezmesh(plot_fh,[0 total_x 0 total_y],20);colormap([0 0 1]);
if polynom_node_order < 5
title(plot_fh,'FontSize',8);
end
%Ezxontour[3,2,4]
subplot(3,2,4); ezcontour(plot_fh,[0 total_x 0 total_y]);
title('Ezcontour','FontSize',8);

%subplot[3,2,5]
subplot(3,2,5);
plot_fh =vectorize(subs(NQxx(i),{a,b},{total_x,total_y}));
ezmesh(plot_fh,[0 total_x 0 total_y],20);colormap([0 0 1]);
if polynom_node_order < 5
title(plot_fh,'FontSize',8);
end
%Ezxontour[3,2,6]
subplot(3,2,6); ezcontour(plot_fh,[0 total_x 0 total_y]);
title('Ezcontour','FontSize',8);

end

%Numeric functions.
if flag==true
N_shape  = subs(N_shape ,{a},{total_x/(polynom_node_order-1)});
N_shape2 = subs(N_shape2,{a,b},{total_x,total_y});
N_shape3 = subs(N_shape3,{a,b},{total_x,total_y});
end

%}
% Several examples;
%                             (9node-27Dof plate element)
%Rank 27 (non-conforming)
%w=[                   1 ...
%                   x     y ...
%               x^2   x*y   y^2 ...
%            x^3 x^2*y    x*y^2  y^3 ...
%         x^4+x^3*y  x^2*y^2   x*y^3+y^4 ...        Pascal tree
%     x^5 x^4*y x^3*y^2  x^2*y^3 x*y^4 y^5 ...
%       x^5*y+x^4*y^2  x^3*y^3  x^2*y^4+x*y^5 ...
%          x^5*y^2+x^4*y^3  x^3*y^4+x^2*y^5 ...
%                x^5*y^3 x^4*y^4 x^3*y^5]
%
%Rank 26
%w=[                   1 ...
%                   x     y ...
%               x^2   x*y   y^2 ...
%            x^3 x^2*y    x*y^2  y^3 ...
%         x^4+x^3*y    x^2*y^2  x*y^3+y^4 ...
%     x^5 x^4*y x^3*y^2  x^2*y^3 x*y^4 y^5 ...
%       x^5*y x^4*y^2  x^3*y^3  x^2*y^4 x*y^5 ...
%          x^5*y^2+x^4*y^3  x^3*y^4+x^2*y^5 ...
%                      x^4*y^4 ];
%
% Rank 26
%w=[                   1 ...
%                   x     y ...
%               x^2   x*y   y^2 ...
%            x^3 x^2*y    x*y^2  y^3 ...
%         x^4 x^3*y   x^2*y^2  x*y^3 y^4 ...
%     x^5 x^4*y x^3*y^2  x^2*y^3 x*y^4 y^5 ...
%       x^5*y+x^4*y^2  x^3*y^3  x^2*y^4+x*y^5 ...
%          x^5*y^2+x^4*y^3  x^3*y^4+x^2*y^5 ...
%                    x^4*y^4 ];
%
% Rank 25
%w=[                   1 ...
%                   x     y ...
%               x^2   x*y   y^2 ...
%            x^3 x^2*y    x*y^2  y^3 ...
%         x^4+x^3*y          x*y^3+y^4 ...
%     x^5 x^4*y x^3*y^2  x^2*y^3 x*y^4 y^5 ...
%       x^5*y x^4*y^2  x^3*y^3  x^2*y^4 x*y^5 ...
%          x^5*y^2 x^4*y^3  x^3*y^4 x^2*y^5 ];
%
% Rank 26
%w=[                   1 ...
%                   x     y ...
%               x^2   x*y   y^2 ...
%            x^3 x^2*y      x*y^2  y^3 ...
%         x^4 x^3*y            x*y^3 y^4 ...
%     x^5 x^4*y x^3*y^2  x^2*y^3 x*y^4 y^5 ...
%       x^5*y x^4*y^2  x^3*y^3  x^2*y^4 x*y^5 ...
%          x^5*y^2+x^4*y^3  x^3*y^4+x^2*y^5 ];
%
% Rank 24
%w=[                   1 ...
%                   x     y ...
%               x^2   x*y   y^2 ...
%            x^3 x^2*y      x*y^2  y^3 ...
%         x^4 x^3*y            x*y^3 y^4 ...
%     x^5 x^4*y x^3*y^2  x^2*y^3 x*y^4 y^5 ...
%       x^5*y+x^4*y^2  x^3*y^3  x^2*y^4+x*y^5 ...
%          x^5*y^2 x^4*y^3  x^3*y^4 x^2*y^5 ];
%
% Rank 24
%w=[                   1 ...
%                   x     y ...
%               x^2   x*y   y^2 ...
%            x^3 x^2*y      x*y^2  y^3 ...
%         x^4 x^3*y            x*y^3 y^4 ...
%     x^5 x^4*y x^3*y^2  x^2*y^3 x*y^4 y^5 ...
%       x^5*y        x^3*y^3      x*y^5 ...
%          x^5*y^2 x^4*y^3  x^3*y^4 x^2*y^5 ];
%Rank 27 (non-conforming)
%w=[                   1 ...
%                   x     y ...
%               x^2   x*y   y^2 ...
%            x^3 x^2*y      x*y^2  y^3 ...
%         x^4+x^3*y    x^2*y^2    x*y^3+y^4 ...
%     x^5+x^4*y   x^3*y^2+x^2*y^3 x*y^4+y^5 ...
%       x^5*y+x^4*y^2  x^3*y^3  x^2*y^4+x*y^5 ...
%          x^5*y^2 x^4*y^3  x^3*y^4 x^2*y^5 ...
%               x^5*y^3  x^4*y^4   x^3*y^5 ]
%Rank 27 (non-conforming)
%w=[                   1 ...
%                   x     y ...
%               x^2   x*y   y^2 ...
%            x^3 x^2*y      x*y^2  y^3 ...
%         x^4      x^3*y+x*y^3      y^4 ...
%     x^5+x^4*y   x^3*y^2 x^2*y^3 x*y^4+y^5 ...
%       x^5*y+x^4*y^2  x^3*y^3  x^2*y^4+x*y^5 ...
%          x^5*y^2 x^4*y^3  x^3*y^4 x^2*y^5 ...
%               x^5*y^3  x^4*y^4   x^3*y^5 ]
%Rank 27 (non-conforming)
%w=[                   1 ...
%                   x     y ...
%               x^2   x*y   y^2 ...
%            x^3 x^2*y      x*y^2  y^3 ...
%         x^4      x^3*y+x*y^3      y^4 ...
%     x^5 x^4*y   x^3*y^2+x^2*y^3 x*y^4 y^5 ...
%       x^5*y+x^4*y^2  x^3*y^3  x^2*y^4+x*y^5 ...
%          x^5*y^2 x^4*y^3+x^3*y^4 x^2*y^5 ...
%               x^5*y^3  x^4*y^4   x^3*y^5 ]
%Rank 27 (non-conforming)
%w=[                   1 ...
%                   x     y ...
%               x^2   x*y   y^2 ...
%            x^3 x^2*y      x*y^2  y^3 ...
%         x^4      x^3*y+x*y^3      y^4 ...
%     x^5 x^4*y   x^3*y^2+x^2*y^3 x*y^4 y^5 ...
%       x^5*y+x^4*y^2  x^3*y^3  x^2*y^4+x*y^5 ...
%          x^5*y^2 x^4*y^3+x^3*y^4 x^2*y^5 ...
%               x^5*y^3  x^4*y^4   x^3*y^5 ]
%