Jan Grošelj, Marjeta Krajnc: Quartic splines on Powell-Sabin triangulations, notebook with examples

Functions

ClearAll[Norma]
(* Norm of a vector x *)
Norma[x_]:=(Sqrt[x.x])

ClearAll[bernstein];
(* Bernstein basis polynomials *)
bernstein[x_,n_,i_]:=Module[{},
If[x==0 && i==0,Return[1]];
If[x==1 && i==n,Return[1]];
Return[Binomial[n,i]*x^(i)*(1-x)^(n-i)]
]

ClearAll[bernsteinT]
(* Bernstein basis polynomials over a triangle *)
bernsteinT[U_,I_,n_]:=Module[{c1,c2,c3},
If[n==0,Return[1]];
If[I[[1]]==0,c1=1,c1=U[[1]]^I[[1]]];
If[I[[2]]==0,c2=1,c2=U[[2]]^I[[2]]];
If[I[[3]]==0,c3=1,c3=U[[3]]^I[[3]]];
Return[n!/(I[[1]]!I[[2]]!I[[3]]!)c1 c2 c3]
]

ClearAll[Baricentricne]
(* Function returns baricentric coordinates of a point T with respect to a triangle with vertices S1, S2, S3. *)
Baricentricne[S1_,S2_,S3_,T_]:=Module[{p0,p1,p2,p3},
p0=Det[{{1,1,1},{S1[[1]],S2[[1]],S3[[1]]},{S1[[2]],S2[[2]],S3[[2]]}}];
p1=Det[{{1,1,1},{T[[1]],S2[[1]],S3[[1]]},{T[[2]],S2[[2]],S3[[2]]}}];
p2=Det[{{1,1,1},{S1[[1]],T[[1]],S3[[1]]},{S1[[2]],T[[2]],S3[[2]]}}];
p3=Det[{{1,1,1},{S1[[1]],S2[[1]],T[[1]]},{S1[[2]],S2[[2]],T[[2]]}}];
Return[{p1,p2,p3}/p0]
]

ClearAll[DeCasteljau]
(* De Casteljau algorithm for control points given in KT. Function returns a matrix of all computed values.  *)
DeCasteljau[KT_,t_]:=Module[{M,n,r,i,j},
n=Length[KT]-1;
M=Table[Table[0,{i,0,n}],{j,0,n}];
For[i=0,i≤n,i++,
M[[i+1,1]]=KT[[i+1]];
];
For[r=1,r≤n,r++,
For[i=0,i≤n-r,i++,
M[[i+1,r+1]]=(1-t)M[[i+1,r]]+t M[[i+2,r]];
];
];
Return[M]
]

ClearAll[DomainPointsB]
(* Function computes domain points (baricentric) on the triangle for a polynomial of degree n. *)
DomainPointsB[n_]:=Module[{M,i,j},
M=Table[{0,0,0},{i,0,n},{j,0,n}];
For[i=0,i≤n,i++,
For[j=0,j≤n-i,j++,
M[[n-i+1,j+1]]={(n-i-j)/n,j/n,i/n};
];
];
Return[M]
]

ClearAll[DomainPoints]
(* Function computes domain points (cartesian) on the triangle with vertices V1, V2, V3 for a polynomial of degree n. *)
DomainPoints[V1_,V2_,V3_,n_]:=Module[{DP,TDP},
DP=DomainPointsB[n];
TDP=Table[DP[[i,j,1]]V1+DP[[i,j,2]]V2+DP[[i,j,3]]V3,{i,1,n+1},{j,1,n+1}];
Return[TDP];
]

ClearAll[BezierKontPoints]
(* Function computes Bezier control points on the triangle with vertices V1, V2, V3 for a polynomial of degree n with ordinates given in KT. *)
BezierKontPoints[V1_,V2_,V3_,KT_,n_]:=Module[{DP,TDP},
DP=DomainPointsB[n];
TDP=Table[Append[DP[[i,j,1]]V1+DP[[i,j,2]]V2+DP[[i,j,3]]V3,KT[[i,j]]],{i,1,n+1},{j,1,n+1}];
Return[TDP];
]

ClearAll[DomainInd]
(* Function computes indices of domain points on the triangle for a polynomial of degree n. *)
DomainInd[n_]:=Module[{M,i,j},
M=Table[{0,0,0},{i,0,n},{j,0,n}];
For[i=0,i≤n,i++,
For[j=0,j≤n-i,j++,
M[[n-i+1,j+1]]={(n-i-j),j,i};
];
];
Return[M]
]

ClearAll[NarisiTrikMrezo]
(* This function helps us to plot a triangle mesh. From a lower triangle matrix with vertices it computes a control mesh. More precisely, it returns a list of polygons that can be used in Graphics3D function. *)
NarisiTrikMrezo[MC_]:=Module[{m,tabPolygon,j,i},
m=Length[MC];
tabPolygon={};
For[j=1,j≤m-2,j++,
For[i=j,i≤m-1,i++,
tabPolygon=Append[tabPolygon,Polygon[{MC[[i,j]],MC[[i+1,j]],MC[[i+1,j+1]]}]];
];
For[i=j,i≤m-2,i++,
tabPolygon=Append[tabPolygon,Polygon[{MC[[i+1,j]],MC[[i+1,j+1]],MC[[i+2,j+1]]}]];
];
];
tabPolygon=Append[tabPolygon,Polygon[{MC[[m-1,m-1]],MC[[m,m-1]],MC[[m,m]]}]];
Return[tabPolygon]
]

ClearAll[NarisiKontPoligon]
(* This function helps us to plot a triangle mesh. From a lower triangle matrix with vertices it computes a control polygon. More precisely, it returns a list of lines of a polygon, that can be used in Graphics3D function. *)
NarisiKontPoligon[MC_]:=Module[{m,tabPolygon,j,i},
m=Length[MC];
tabPolygon={};
For[j=1,j≤m-2,j++,
For[i=j,i≤m-1,i++,
tabPolygon=Append[tabPolygon,Line[{MC[[i,j]],MC[[i+1,j]],MC[[i+1,j+1]],MC[[i,j]]}]];
];
For[i=j,i≤m-2,i++,
tabPolygon=Append[tabPolygon,Line[{MC[[i+1,j]],MC[[i+1,j+1]],MC[[i+2,j+1]],MC[[i+1,j]]}]];
];
];
tabPolygon=Append[tabPolygon,Line[{MC[[m-1,m-1]],MC[[m,m-1]],MC[[m,m]],MC[[m-1,m-1]]}]];
Return[tabPolygon]
]

ClearAll[ControlPointsToMatrix]
(* From a list of control points C function constitutes a matrix, so that the neighbouring control points are together. 
 We assume that in a list the basis functions are ordered as (n,0,0),(n-1,1,0),(n-1,0,1),(n-2,2,0),(n-2,1,1),(n-2,0,2), ..., (i,n-i,0),(i,n-i-1,i), (i,n-i-2,2),...,(i,0,n-i), ..., (0,0,n). 
In the returned lower triangular matrix the i-th line consist of control points with indexes (i,n-i,0),(i,n-i-1,i), (i,n-i-2,2),...,(i,0,n-i).
Here n is a degree of a surface. *)
ControlPointsToMatrix[C_,n_]:=Module[{M,z,i,j},
M=Table[{},{i,0,n},{j,0,n}];
M[[1,1]]=C[[1]];
z=2;
For[i=1,i≤n,i++,
For[j=z,j≤z+i,j++,
M[[i+1,j-z+1]]=C[[j]];
];
z=z+i+1;
];
Return[M]
]

ClearAll[ComputeValuesForGraph]
(* Function returns a (NN\[Cross]NN) matrix of values of the Bezier surface with given control points (in a list MC), where NN is the number of points on each edge.   *)
(* This function does not use de Casteljau algorithm! *)
ComputeValuesForGraph[MC_,NN_]:=Module[{ss,nn,i,j,u,v,w,V},
(* from the control points compute the surface ss *)
ss=0;
nn=Length[MC]-1;
For[i=0,i≤nn,i++,
For[j=0,j≤i,j++,
ss=ss+MC[[i+1,j+1]]nn!/((i-j)!(nn-i)! j!)v^(i-j) u^(nn-i)w^(j);
]
];
(* Compute the matrix of values *)
V=Table[{},{i,0,NN},{j,0,NN}];
For[j=0,j≤NN,j++,
For[i=j,i≤NN,i++,
V[[i+1,j+1]]=ss/.{u->(NN-i)/NN,v->(i-j)/NN,w->j/NN};
]
];
Return[V]
]

ClearAll[PSSplineOrdinatesNew]
(* Function returns Bezier ordinates of a spline over one triangle with vertices V1,V2,V3. Input values R12,R13,R23 are vertices on lines betwen V1V2, V1V3, V2V3, respectively. Point Z is a midpoint from PS-split, that splits a triangle into six smaller triangles.     
A list 'PodT' gives the interpolation values for points (see fij in a paper). Similarly, lists 'PodW1','PodW2' give the interpolation values for edges (see gij in a paper). Lists 'omega' and 'BarWk1', 'BarWk2', 'BarWk3' represent baricentric coordinates as in the paper.  
*)
PSSplineOrdinatesNew[V1_,V2_,V3_,R12_,R23_,R13_,Z_,omega_,BarWk1_,BarWk2_,BarWk3_,PodT_,PodW1_,PodW2_,VredRob_]:=Module[{vv,mm,nn,zz,uu,oo,tt,qq,pp,rr,M1,M2,M3,M4,M5,M6,ω,WW1,WW2,WW3,TVV,resVi,resZi,resNi,BarLam,λ,i,resTi,resPi,eqo12,eqo31,eqo23,ν,resO1,resO2,bbb,BarBeta,β,vseVred,Nu1A,Nu1B,Bnu1,Nu2A,Nu2B,Bnu2,Nu3A,Nu3B,Bnu3,IzbNu,resRobA,RESALL},
TVV={V1,V2,V3};
M1=Table[0,{i,1,4},{j,1,4}];
M1[[1,1]]=zz;
M1[[2]]={(Subscript[nn, 1]),(Subscript[tt, 1,2]),0,0};
M1[[3]]={(Subscript[mm, 1]),(Subscript[oo, 1,2]),(Subscript[pp, 1,2]),0};
M1[[4]]={(Subscript[vv, 1]),(Subscript[uu, 1,2]),(Subscript[qq, 1,2]),(Subscript[rr, 1,2])};
M2=Table[0,{i,1,4},{j,1,4}];
M2[[1,1]]=zz;
M2[[2]]={(Subscript[nn, 2]),(Subscript[tt, 1,2]),0,0};
M2[[3]]={(Subscript[mm, 2]),(Subscript[oo, 2,1]),(Subscript[pp, 1,2]),0};
M2[[4]]={(Subscript[vv, 2]),(Subscript[uu, 2,1]),(Subscript[qq, 2,1]),(Subscript[rr, 1,2])};
M3=Table[0,{i,1,4},{j,1,4}];
M3[[1,1]]=zz;
M3[[2]]={(Subscript[nn, 2]),(Subscript[tt, 2,3]),0,0};
M3[[3]]={(Subscript[mm, 2]),(Subscript[oo, 2,3]),(Subscript[pp, 2,3]),0};
M3[[4]]={(Subscript[vv, 2]),(Subscript[uu, 2,3]),(Subscript[qq, 2,3]),(Subscript[rr, 2,3])};
M4=Table[0,{i,1,4},{j,1,4}];
M4[[1,1]]=zz;
M4[[2]]={(Subscript[nn, 3]),(Subscript[tt, 2,3]),0,0};
M4[[3]]={(Subscript[mm, 3]),(Subscript[oo, 3,2]),(Subscript[pp, 2,3]),0};
M4[[4]]={(Subscript[vv, 3]),(Subscript[uu, 3,2]),(Subscript[qq, 3,2]),(Subscript[rr, 2,3])};
M5=Table[0,{i,1,4},{j,1,4}];
M5[[1,1]]=zz;
M5[[2]]={(Subscript[nn, 3]),(Subscript[tt, 1,3]),0,0};
M5[[3]]={(Subscript[mm, 3]),(Subscript[oo, 3,1]),(Subscript[pp, 1,3]),0};
M5[[4]]={(Subscript[vv, 3]),(Subscript[uu, 3,1]),(Subscript[qq, 3,1]),(Subscript[rr, 1,3])};
M6=Table[0,{i,1,4},{j,1,4}];
M6[[1,1]]=zz;
M6[[2]]={(Subscript[nn, 1]),(Subscript[tt, 1,3]),0,0};
M6[[3]]={(Subscript[mm, 1]),(Subscript[oo, 1,3]),(Subscript[pp, 1,3]),0};
M6[[4]]={(Subscript[vv, 1]),(Subscript[uu, 1,3]),(Subscript[qq, 1,3]),(Subscript[rr, 1,3])};
(* Choose WWi *)
(Subscript[ω, 1])=omega[[1]];
(Subscript[ω, 2])=omega[[2]];
(Subscript[ω, 3])=omega[[3]];
WW1=(Subscript[ω, 1])V1 +(1-(Subscript[ω, 1]))Z;
WW2=(Subscript[ω, 2])V2 +(1-(Subscript[ω, 2]))Z;
WW3=(Subscript[ω, 3])V3 +(1-(Subscript[ω, 3]))Z;
(* Equations for interpolation in triangle points. *)
resVi=Join[Flatten[Table[{(Subscript[vv, i])->PodT[[i,1]],(Subscript[mm, i])->PodT[[i,1]]+1/3({PodT[[i,2]],PodT[[i,3]]}.(Z-TVV[[i]]))},{i,1,3}]],
{(Subscript[uu, 1,2])->PodT[[1,1]]+1/3({PodT[[1,2]],PodT[[1,3]]}.(R12-V1)),(Subscript[uu, 1,3])->PodT[[1,1]]+1/3({PodT[[1,2]],PodT[[1,3]]}.(R13-V1)),(Subscript[uu, 2,1])->PodT[[2,1]]+1/3({PodT[[2,2]],PodT[[2,3]]}.(R12-V2)),(Subscript[uu, 2,3])->PodT[[2,1]]+1/3({PodT[[2,2]],PodT[[2,3]]}.(R23-V2)),
(Subscript[uu, 3,2])->PodT[[3,1]]+1/3({PodT[[3,2]],PodT[[3,3]]}.(R23-V3)),(Subscript[uu, 3,1])->PodT[[3,1]]+1/3({PodT[[3,2]],PodT[[3,3]]}.(R13-V3))}];
(* Equations for interpolation in a middle point. *)
resZi={zz->PodW1[[1]]};
resNi={ (Subscript[nn, 1])->PodW1[[1]]+1/3({PodW1[[2]],PodW1[[3]]}.(V1-Z)),(Subscript[nn, 2])->PodW1[[1]]+1/3({PodW1[[2]],PodW1[[3]]}.(V2-Z)),(Subscript[nn, 3])->PodW1[[1]]+1/3({PodW1[[2]],PodW1[[3]]}.(V3-Z))};
(* Let R12= (Subscript[λ, 1,2]) V1 + (1-(Subscript[λ, 1,2])) V2, R23= (Subscript[λ, 2,3]) V2 + (1-(Subscript[λ, 2,3])) V3, R13= (Subscript[λ, 1,3]) V1 + (1-(Subscript[λ, 1,3])) V3. Then we can determine ordinates t12,t23,t13. *)
BarLam=Simplify[{(Subscript[λ, 1,2])->Baricentricne[V1,V2,Z,R12][[1]],(Subscript[λ, 2,3])->Baricentricne[V2,V3,Z,R23][[1]],(Subscript[λ, 1,3])->Baricentricne[V1,V3,Z,R13][[1]]}];
resTi=Simplify[{(Subscript[tt, 1,2])->(Subscript[λ, 1,2])(Subscript[nn, 1])+(1-(Subscript[λ, 1,2]))(Subscript[nn, 2]),(Subscript[tt, 2,3])->(Subscript[λ, 2,3])(Subscript[nn, 2])+(1-(Subscript[λ, 2,3]))(Subscript[nn, 3]),(Subscript[tt, 1,3])->(Subscript[λ, 1,3])(Subscript[nn, 1])+(1-(Subscript[λ, 1,3]))(Subscript[nn, 3])}];
resPi=Simplify[{(Subscript[pp, 1,2])->(Subscript[λ, 1,2])(Subscript[oo, 1,2])+(1-(Subscript[λ, 1,2]))(Subscript[oo, 2,1]),(Subscript[pp, 2,3])->(Subscript[λ, 2,3])(Subscript[oo, 2,3])+(1-(Subscript[λ, 2,3]))(Subscript[oo, 3,2]),(Subscript[pp, 1,3])->(Subscript[λ, 1,3])(Subscript[oo, 1,3])+(1-(Subscript[λ, 1,3]))(Subscript[oo, 3,1])}];
(* Equations for interpolation of directional derivatives. *)
i=1;
eqo12=Simplify[(BarWk1[[1]] (Subscript[vv, i])+ BarWk1[[2]]  (Subscript[mm, i])+ BarWk1[[3]] (Subscript[uu, i,i+1]))bernsteinT[{ (Subscript[ω, i]),1- (Subscript[ω, i]),0},{2,0,0},2]+(BarWk1[[1]] (Subscript[mm, i])+ BarWk1[[2]]  (Subscript[nn, i])+ BarWk1[[3]] (Subscript[oo, i,i+1]))bernsteinT[{ (Subscript[ω, i]),1- (Subscript[ω, i]),0},{1,1,0},2]+(BarWk1[[1]] (Subscript[nn, i])+ BarWk1[[2]] zz+ BarWk1[[3]] (Subscript[tt, i,i+1]))bernsteinT[{ (Subscript[ω, i]),1- (Subscript[ω, i]),0},{0,2,0},2]-1/3PodW2[[i]]];
i=2;
eqo23=Simplify[(BarWk2[[1]] (Subscript[vv, i])+ BarWk2[[2]]  (Subscript[mm, i])+ BarWk2[[3]] (Subscript[uu, i,i+1]))bernsteinT[{ (Subscript[ω, i]),1- (Subscript[ω, i]),0},{2,0,0},2]+(BarWk2[[1]] (Subscript[mm, i])+ BarWk2[[2]]  (Subscript[nn, i])+ BarWk2[[3]] (Subscript[oo, i,i+1]))bernsteinT[{ (Subscript[ω, i]),1- (Subscript[ω, i]),0},{1,1,0},2]+(BarWk2[[1]] (Subscript[nn, i])+ BarWk2[[2]] zz+ BarWk2[[3]] (Subscript[tt, i,i+1]))bernsteinT[{ (Subscript[ω, i]),1- (Subscript[ω, i]),0},{0,2,0},2]-1/3PodW2[[i]]];
i=3;
eqo31=Simplify[(BarWk3[[1]](Subscript[vv, i])+BarWk3[[2]] (Subscript[mm, i])+ BarWk3[[3]](Subscript[uu, i,Mod[i+1,3]]))bernsteinT[{ (Subscript[ω, i]),1- (Subscript[ω, i]),0},{2,0,0},2]+(BarWk3[[1]](Subscript[mm, i])+ BarWk3[[2]] (Subscript[nn, i])+ BarWk3[[3]](Subscript[oo, i,Mod[i+1,3]]))bernsteinT[{ (Subscript[ω, i]),1- (Subscript[ω, i]),0},{1,1,0},2]+(BarWk3[[1]](Subscript[nn, i])+ BarWk3[[2]]zz+ BarWk3[[3]](Subscript[tt, Mod[i+1,3],i]))bernsteinT[{ (Subscript[ω, i]),1- (Subscript[ω, i]),0},{0,2,0},2]-1/3PodW2[[i]]];
resO1=Simplify[Flatten[{Solve[eqo12==0,(Subscript[oo, 1,2])],Solve[eqo23==0,(Subscript[oo, 2,3])],Solve[eqo31==0,(Subscript[oo, 3,1])]}]];
bbb[1]=Baricentricne[V1,Z,R12,R13];
bbb[2]=Baricentricne[V2,Z,R23,R12];
bbb[3]=Baricentricne[V3,Z,R13,R23];
BarBeta=Simplify[Flatten[Table[Table[(Subscript[β, i,j])->bbb[i][[j]],{j,1,3}],{i,1,3}]]];
resO2={(Subscript[oo, 1,3])->(Subscript[β, 1,1]) (Subscript[mm, 1])+(Subscript[β, 1,2])(Subscript[nn, 1])+(Subscript[β, 1,3])(Subscript[oo, 1,2]),(Subscript[oo, 2,1])->(Subscript[β, 2,1]) (Subscript[mm, 2])+(Subscript[β, 2,2])(Subscript[nn, 2])+(Subscript[β, 2,3])(Subscript[oo, 2,3]),(Subscript[oo, 3,2])->(Subscript[β, 3,1]) (Subscript[mm, 3])+(Subscript[β, 3,2])(Subscript[nn, 3])+(Subscript[β, 3,3])(Subscript[oo, 3,1])}//.BarBeta;
resO1=Simplify[resO1/.resZi/.resTi/.resNi/.resVi/.BarLam];
resO2=Simplify[resO2/.resO1/.resZi/.resTi/.resNi/.resVi/.BarLam//.BarBeta];
vseVred=Flatten[{resVi,Simplify[resNi/.resVi],Simplify[resZi/.resVi/.BarLam],Simplify[resTi/.resNi/.resVi/.BarLam],resO1,resO2,resPi/.resO1/.resO2/.BarLam}];
(* Boundary values are set to zero. *)
resRobA={(Subscript[qq, 1,2])->VredRob[[1]],(Subscript[qq, 2,1])->VredRob[[2]],(Subscript[qq, 2,3])->VredRob[[3]],(Subscript[qq, 3,2])->VredRob[[4]],(Subscript[qq, 3,1])->VredRob[[5]],(Subscript[qq, 1,3])->VredRob[[6]],(Subscript[rr, 1,2])->VredRob[[7]],(Subscript[rr, 2,3])->VredRob[[8]],(Subscript[rr, 1,3])->VredRob[[9]]};
Return[{M1,M2,M3,M4,M5,M6}/.vseVred/.resRobA];
]


ClearAll[DeCasteljau3D]
(* De Casteljau algorithm with respect to 'v' for a triangle patch with given ordinates in a matrix M.*)
DeCasteljau3D[M_,v_]:=Module[{n,M1,M2,k,i,j,TT},
n=Length[M]-1;
TT={M};
M1=M;
For[k=1,k≤n,k++,
(* k-th step *)
M2=Table[0,{i,1,n+1-k},{j,1,n+1-k}];
For[i=0,i≤n-k,i++,
For[j=0,j≤i,j++,
M2[[i+1,j+1]]=v[[1]]M1[[i+2,j+1]]+v[[2]]M1[[i+2,j+2]]+v[[3]]M1[[i+1,j+1]];
];
] ;
TT=Append[TT,M2];
M1=M2;
];
Return[TT];
]

ClearAll[Blossom]
(* Function computes a blossom with respect to values 'TabV' with Bezier ordinates given in a matrix 'M' *)
Blossom[M_,TabV_]:=Module[{n,M1,M2,k,i,j,TT},
n=Length[M]-1;
TT={M};
M1=M;
For[k=1,k≤n,k++,
(* k-th step *)
M2=Table[0,{i,1,n+1-k},{j,1,n+1-k}];
For[i=0,i≤n-k,i++,
For[j=0,j≤i,j++,
M2[[i+1,j+1]]=TabV[[k,1]]M1[[i+2,j+1]]+TabV[[k,2]]M1[[i+2,j+2]]+TabV[[k,3]]M1[[i+1,j+1]];
];
] ;
TT=Append[TT,M2];
M1=M2;
];
Return[TT];
]

ClearAll[BezOrdinatesOneTRiangle]
(* Function returns Bezier ordinates of the spline over one triangle. Interpolation values are given in VF, DF, HF, VG.  *)
BezOrdinatesOneTRiangle[PSTriangle_,VF_,DF_,HF_,VG_,smeri_]:=Module[{b,M,V1,V2,V3,R12,R23,R13,Z,n,i,j,k,res1,v1,v2,TT,odv1a,odv1b,odv2a,odv2b,odv2c,res2,res3,smerB,e1,e2,e3,pom,res4,res5,res6,res7,λ,σ,σσ,BarLam,BarSigma,resAll},
n=4;
M=Table[Table[0,{i,1,n+1},{j,1,n+1}],{k,1,6}];
V1=PSTriangle[[1,1]];
V2=PSTriangle[[2,1]];
V3=PSTriangle[[4,1]];
R12=PSTriangle[[1,2]];
R23=PSTriangle[[3,2]];
R13=PSTriangle[[5,2]];
Z=PSTriangle[[1,3]];
For[k=1,k≤6,k++,
For[i=0,i≤n,i++,
For[j=0,j≤i,j++,
M[[k,i+1,j+1]]=(Subscript[b[k], i-j,j,n-i])
];
] ;
];
(* Interpolation at points V1,V2,V3.  *)
res1={};
res2={};
For[k=1,k≤6,k++,
res1=Append[res1,M[[k,5,1]]->VF[[k]]];
v1={-1,0,1};
v2={-1,1,0};
TT=Blossom[M[[k]],{v1,{1,0,0},{1,0,0},{1,0,0}}];
odv1a=TT[[2,4,1]];
TT=Blossom[M[[k]],{v2,{1,0,0},{1,0,0},{1,0,0}}];
odv1b=TT[[2,4,1]];
TT=Blossom[M[[k]],{v1,v1,{1,0,0},{1,0,0}}];
odv2a=TT[[3,3,1]];
TT=Blossom[M[[k]],{v2,v2,{1,0,0},{1,0,0}}];
odv2b=TT[[3,3,1]];
TT=Blossom[M[[k]],{v1,v2,{1,0,0},{1,0,0}}];
odv2c=TT[[3,3,1]];
res2=Append[res2,Flatten[{Flatten[Solve[DF[[k]].(PSTriangle[[k,3]]-PSTriangle[[k,1]])==4odv1a,M[[k,4,1]]]],Flatten[Solve[DF[[k]].(PSTriangle[[k,2]]-PSTriangle[[k,1]])==4odv1b,M[[k,5,2]]]],Flatten[Solve[(PSTriangle[[k,3]]-PSTriangle[[k,1]]).HF[[k]].(PSTriangle[[k,3]]-PSTriangle[[k,1]])==12odv2a,M[[k,3,1]]]],Flatten[Solve[(PSTriangle[[k,2]]-PSTriangle[[k,1]]).HF[[k]].(PSTriangle[[k,2]]-PSTriangle[[k,1]])==12odv2b,M[[k,5,3]]]],Flatten[Solve[(PSTriangle[[k,2]]-PSTriangle[[k,1]]).HF[[k]].(PSTriangle[[k,3]]-PSTriangle[[k,1]])==12odv2c,M[[k,4,2]]]]}]];
];
res2=Flatten[res2];
(* Interpolation at points Rij: *)
smerB={Baricentricne[V1,R12,Z,smeri[[1,2]]]-Baricentricne[V1,R12,Z,smeri[[1,1]]],Baricentricne[V2,R23,Z,smeri[[2,2]]]-Baricentricne[V2,R23,Z,smeri[[2,1]]],Baricentricne[V3,R13,Z,smeri[[3,2]]]-Baricentricne[V3,R13,Z,smeri[[3,1]]]};
res3={};
For[k=1,k≤3,k++,
e1=M[[2(k-1)+1,5,5]]==VG[[k,1]];
TT=Blossom[M[[2(k-1)+1]],{smerB[[k]],{0,1,0},{0,1,0},{0,1,0}}];
e2=4TT[[2,4,4]]==VG[[k,2]];
TT=Blossom[M[[2(k-1)+1]],{smerB[[k]],smerB[[k]],{0,1,0},{0,1,0}}];
e3=12TT[[3,3,3]]==VG[[k,3]];
res3=Join[res3,Flatten[Solve[{e1,e2,e3},{M[[2(k-1)+1,5,5]],M[[2(k-1)+1,4,4]],M[[2(k-1)+1,3,3]]}]]];
pom={M[[2(k-1)+2,5,5]]->(M[[2(k-1)+1,5,5]]/.res3),M[[2(k-1)+2,4,4]]->(M[[2(k-1)+1,4,4]]/.res3),M[[2(k-1)+2,3,3]]->(M[[2(k-1)+1,3,3]]/.res3)};
res3=Join[res3,pom];
];
(* C2 conditions over Z-Rij*)
res4={};
i=1;
j=2;
For[k=1,k≤3,k++,
res4=Append[res4,Solve[{M[[1+2(i-1),6-k,6-k]]==(Subscript[λ, i,j])M[[1+2(i-1),6-k,5-k]]+(Subscript[λ, j,i])M[[2+2(i-1),6-k,5-k]],
(Subscript[λ, i,j])^2(M[[1+2(i-1),6-k,6-k]]-2M[[1+2(i-1),6-k,5-k]]+M[[1+2(i-1),6-k,4-k]])==(Subscript[λ, j,i])^2(M[[2+2(i-1),6-k,6-k]]-2M[[2+2(i-1),6-k,5-k]]+M[[2+2(i-1),6-k,4-k]])},{M[[1+2(i-1),6-k,5-k]],M[[2+2(i-1),6-k,5-k]]}][[1]]]];
i=2;
j=3;
For[k=1,k≤3,k++,
res4=Append[res4,Solve[{M[[1+2(i-1),6-k,6-k]]==(Subscript[λ, i,j])M[[1+2(i-1),6-k,5-k]]+(Subscript[λ, j,i])M[[2+2(i-1),6-k,5-k]],
(Subscript[λ, i,j])^2(M[[1+2(i-1),6-k,6-k]]-2M[[1+2(i-1),6-k,5-k]]+M[[1+2(i-1),6-k,4-k]])==(Subscript[λ, j,i])^2(M[[2+2(i-1),6-k,6-k]]-2M[[2+2(i-1),6-k,5-k]]+M[[2+2(i-1),6-k,4-k]])},{M[[1+2(i-1),6-k,5-k]],M[[2+2(i-1),6-k,5-k]]}][[1]]]];
i=3;
j=1;
For[k=1,k≤3,k++,
res4=Append[res4,Solve[{M[[1+2(i-1),6-k,6-k]]==(Subscript[λ, i,j])M[[1+2(i-1),6-k,5-k]]+(Subscript[λ, j,i])M[[2+2(i-1),6-k,5-k]],
(Subscript[λ, i,j])^2(M[[1+2(i-1),6-k,6-k]]-2M[[1+2(i-1),6-k,5-k]]+M[[1+2(i-1),6-k,4-k]])==(Subscript[λ, j,i])^2(M[[2+2(i-1),6-k,6-k]]-2M[[2+2(i-1),6-k,5-k]]+M[[2+2(i-1),6-k,4-k]])},{M[[1+2(i-1),6-k,5-k]],M[[2+2(i-1),6-k,5-k]]}][[1]]]];
res4=Simplify[Flatten[res4]];
(* C1 conditions over Z-Vk:*)
res5={};
i=2;
j=3;
res5=Join[res5,{M[[i,2,1]]->1/(Subscript[σ, i,j,3])(M[[j,3,2]]-(Subscript[σ, i,j,1])M[[i,3,1]] -(Subscript[σ, i,j,2])M[[i,3,2]]),M[[j,2,1]]->1/(Subscript[σ, i,j,3])(M[[j,3,2]]-(Subscript[σ, i,j,1])M[[i,3,1]] -(Subscript[σ, i,j,2])M[[i,3,2]])}];
i=4;
j=5;
res5=Join[res5,{M[[i,2,1]]->1/(Subscript[σ, i,j,3])(M[[j,3,2]]-(Subscript[σ, i,j,1])M[[i,3,1]] -(Subscript[σ, i,j,2])M[[i,3,2]]),M[[j,2,1]]->1/(Subscript[σ, i,j,3])(M[[j,3,2]]-(Subscript[σ, i,j,1])M[[i,3,1]] -(Subscript[σ, i,j,2])M[[i,3,2]])}];
i=6;
j=1;
res5=Join[res5,{M[[i,2,1]]->1/(Subscript[σ, i,j,3])(M[[j,3,2]]-(Subscript[σ, i,j,1])M[[i,3,1]] -(Subscript[σ, i,j,2])M[[i,3,2]]),M[[j,2,1]]->1/(Subscript[σ, i,j,3])(M[[j,3,2]]-(Subscript[σ, i,j,1])M[[i,3,1]] -(Subscript[σ, i,j,2])M[[i,3,2]])}];
(* C2 conditions over Z-Rij*)
res6={};
k=4;
i=1;
j=2;
res6={M[[1+2(i-1),6-k,6-k]]->(Subscript[λ, i,j])M[[1+2(i-1),6-k,5-k]]+(Subscript[λ, j,i])M[[2+2(i-1),6-k,5-k]],M[[2+2(i-1),6-k,6-k]]->(Subscript[λ, i,j])M[[1+2(i-1),6-k,5-k]]+(Subscript[λ, j,i])M[[2+2(i-1),6-k,5-k]]};
i=2;
j=3;
res6=Join[res6,{M[[1+2(i-1),6-k,6-k]]->(Subscript[λ, i,j])M[[1+2(i-1),6-k,5-k]]+(Subscript[λ, j,i])M[[2+2(i-1),6-k,5-k]],M[[2+2(i-1),6-k,6-k]]->(Subscript[λ, i,j])M[[1+2(i-1),6-k,5-k]]+(Subscript[λ, j,i])M[[2+2(i-1),6-k,5-k]]}];
i=3;
j=1;
res6=Join[res6,{M[[1+2(i-1),6-k,6-k]]->(Subscript[λ, i,j])M[[1+2(i-1),6-k,5-k]]+(Subscript[λ, j,i])M[[2+2(i-1),6-k,5-k]],M[[2+2(i-1),6-k,6-k]]->(Subscript[λ, i,j])M[[1+2(i-1),6-k,5-k]]+(Subscript[λ, j,i])M[[2+2(i-1),6-k,5-k]]}];
res7=Table[M[[k,1,1]]->(Subscript[σσ, 1])M[[1,2,1]]+(Subscript[σσ, 2])M[[3,2,1]]+(Subscript[σσ, 3])M[[5,2,1]],{k,1,6}];
BarLam={(Subscript[λ, 1,2])->Baricentricne[V1,V2,V3,R12][[1]],(Subscript[λ, 2,1])->Baricentricne[V1,V2,V3,R12][[2]],(Subscript[λ, 2,3])->Baricentricne[V2,V3,V1,R23][[1]],(Subscript[λ, 3,2])->Baricentricne[V2,V3,V1,R23][[2]],(Subscript[λ, 1,3])->Baricentricne[V1,V3,V2,R13][[1]],(Subscript[λ, 3,1])->Baricentricne[V1,V3,V2,R13][[2]]};
BarSigma=Join[Table[(Subscript[σ, 2,3,k])->Baricentricne[V2,R12,Z,R23][[k]],{k,1,3}],Table[(Subscript[σ, 4,5,k])->Baricentricne[V3,R23,Z,R13][[k]],{k,1,3}],Table[(Subscript[σ, 6,1,k])->Baricentricne[V1,R13,Z,R12][[k]],{k,1,3}],Table[(Subscript[σσ, k])->Baricentricne[V1,V2,V3,Z][[k]],{k,1,3}]];
resAll=Simplify[Join[res1,res2,res3,res4,res5,res6,res7]/.BarLam/.BarSigma];
Return[M//.resAll];
]


(* ----------------- FUNCTIONS FOR A WORK WITH TRIANGULATIONS  --------------------- *)

ClearAll[DeterminePSTriangles]
(* Function returns a list of six triangles from a PS split. *)
DeterminePSTriangles[V1_,V2_,V3_,R12_,R23_,R13_,Z_]:={{V1,R12,Z},{V2,R12,Z},{V2,R23,Z},{V3,R23,Z},{V3,R13,Z},{V1,R13,Z}}

ClearAll[DolociParAlpha]
(* Function returns a matrix A with six columns that represent the values needed for a construction of  vertex basis functions.  *)
DolociParAlpha[Q1_,Q2_,Q3_,V1_]:=Module[{bb,MMM,smerX,smerY,a00,a10,a01,a20,a11,a02,i,A},
bb=Simplify[Baricentricne[Q1,Q2,Q3,V1]];
MMM={{{0,0,0},{0,0,0},{1,0,0}},{{0,0,0},{0,0,0},{0,1,0}},{{0,0,0},{1,0,0},{0,0,0}},{{0,0,0},{0,0,0},{0,0,1}},{{0,0,0},{0,1,0},{0,0,0}},{{1,0,0},{0,0,0},{0,0,0}}};
smerX=Baricentricne[Q1,Q2,Q3,{1,0}]-Baricentricne[Q1,Q2,Q3,{0,0}];
smerY=Baricentricne[Q1,Q2,Q3,{0,1}]-Baricentricne[Q1,Q2,Q3,{0,0}];
a00=Table[Last[Last[Last[Blossom[MMM[[i]],{bb,bb,bb}]]]],{i,1,6}];
a10=2Table[Last[Last[Last[Blossom[MMM[[i]],{smerX,bb,bb}]]]],{i,1,6}];
a01=2Table[Last[Last[Last[Blossom[MMM[[i]],{smerY,bb,bb}]]]],{i,1,6}];
a20=3/2 2Table[Last[Last[Last[Blossom[MMM[[i]],{smerX,smerX,bb}]]]],{i,1,6}];
a02=3/2 2Table[Last[Last[Last[Blossom[MMM[[i]],{smerY,smerY,bb}]]]],{i,1,6}];
a11=3/2 2Table[Last[Last[Last[Blossom[MMM[[i]],{smerX,smerY,bb}]]]],{i,1,6}];
A=Simplify[{a00,a10,a01,a20,a11,a02}];
Return[A]
]


ClearAll[Baricentricne1D]
Baricentricne1D[Q1_,Q2_,R_]:={((R-Q2).(Q1-Q2))/((Q1-Q2).(Q1-Q2)),1-((R-Q2).(Q1-Q2))/((Q1-Q2).(Q1-Q2))}

ClearAll[DolociParBeta]
(* Function returns a matrix B with three columns that represent the values needed for a construction of  edge basis functions.  *)
DolociParBeta[V1_,V2_,R12_,Q1_,Q2_]:=Module[{lam,pom,ksi1,ksi2,bb,smerW,smerB,MMM,b0,b1,b2,a,CC,B,a1},
{ksi1,ksi2}=Baricentricne1D[Q1,Q2,R12];
bb={0,ksi1,ksi2};
smerW={Q1,Q1+(Q2-Q1)/Norma[Q2-Q1]};
smerB=Prepend[Baricentricne1D[Q1,Q2,smerW[[2]]]-Baricentricne1D[Q1,Q2,smerW[[1]]],0];
MMM={{{0,0,0},{0,0,0},{0,0,1}},{{0,0,0},{0,1,0},{0,0,0}},{{1,0,0},{0,0,0},{0,0,0}}};
b0=Table[Last[Last[Last[Blossom[MMM[[i]],{bb,bb,bb}]]]],{i,1,3}];
b1=2Simplify[Table[Last[Last[Last[Blossom[MMM[[i]],{smerB,bb,bb}]]]],{i,1,3}]];
b2=2Simplify[Table[Last[Last[Last[Blossom[MMM[[i]],{smerB,smerB,bb}]]]],{i,1,3}]];
{a,a1}=Baricentricne1D[V1,V2,R12];
CC=2 a a1;
B={CC b0,CC b1,3/2CC b2};
Return[B];
]


ClearAll[BaznaFunZaTocko]
(* Function returns a list of six matrices with Bezier ordinates that determine a vertex basis function indexed by 'indT'. *)
BaznaFunZaTocko[PSTriangle_,podA_,indT_]:=Module[{n,M,V1,V2,V3,R12,R13,R23,Z,b,k,i,j,f,podF,ind1,ind2,VV,RR1,RR2,H,L1,L2,s,sig},
n=4;
M=Table[Table[0,{i,1,n+1},{j,1,n+1}],{k,1,6}];
V1=PSTriangle[[1,1]];
V2=PSTriangle[[2,1]];
V3=PSTriangle[[4,1]];
R12=PSTriangle[[1,2]];
R23=PSTriangle[[3,2]];
R13=PSTriangle[[5,2]];
Z=PSTriangle[[1,3]];
podF={(Subscript[f, 0,0])->podA[[1]],(Subscript[f, 1,0])->podA[[2]],(Subscript[f, 0,1])->podA[[3]],(Subscript[f, 2,0])->podA[[4]],(Subscript[f, 1,1])->podA[[5]],(Subscript[f, 0,2])->podA[[6]]};
H={{(Subscript[f, 2,0]),(Subscript[f, 1,1])},{(Subscript[f, 1,1]),(Subscript[f, 0,2])}}/.podF;
If[indT==1,ind1=1;ind2=6;VV=V1;RR1=R12;RR2=R13;L1=Baricentricne[V1,V2,V3,R12][[1]];L2=Baricentricne[V1,V2,V3,R13][[1]];sig=Baricentricne[V1,V2,V3,Z][[1]]];
If[indT==2,ind1=2;ind2=3;VV=V2;RR1=R12;RR2=R23;L1=Baricentricne[V1,V2,V3,R12][[2]];L2=Baricentricne[V1,V2,V3,R23][[2]];sig=Baricentricne[V1,V2,V3,Z][[2]]];
If[indT==3,ind1=4;ind2=5;VV=V3;RR1=R23;RR2=R13;L1=Baricentricne[V1,V2,V3,R23][[3]];L2=Baricentricne[V1,V2,V3,R13][[3]];sig=Baricentricne[V1,V2,V3,Z][[3]]];
(* ----------- *)
M[[ind1,5,1]]=(Subscript[f, 0,0]);
M[[ind1,4,1]]=(Subscript[f, 0,0])+1/4{(Subscript[f, 1,0]),(Subscript[f, 0,1])}.(Z-VV);
M[[ind1,5,2]]=(Subscript[f, 0,0])+1/4{(Subscript[f, 1,0]),(Subscript[f, 0,1])}.(RR1-VV);
M[[ind1,3,1]]=2((Subscript[f, 0,0])+1/4{(Subscript[f, 1,0]),(Subscript[f, 0,1])}.(Z-VV))-(Subscript[f, 0,0])+1/12(Z-VV).(H.(Z-VV));
M[[ind1,5,3]]=2((Subscript[f, 0,0])+1/4{(Subscript[f, 1,0]),(Subscript[f, 0,1])}.(RR1-VV))-(Subscript[f, 0,0])+1/12(RR1-VV).(H.(RR1-VV));
M[[ind1,4,2]]=((Subscript[f, 0,0])+1/4{(Subscript[f, 1,0]),(Subscript[f, 0,1])}.(Z-VV))+((Subscript[f, 0,0])+1/4{(Subscript[f, 1,0]),(Subscript[f, 0,1])}.(RR1-VV))-(Subscript[f, 0,0])+1/12(Z-VV).(H.(RR1-VV));
(* ----------- *)
M[[ind2,5,1]]=(Subscript[f, 0,0]);
M[[ind2,4,1]]=(Subscript[f, 0,0])+1/4{(Subscript[f, 1,0]),(Subscript[f, 0,1])}.(Z-VV);
M[[ind2,5,2]]=(Subscript[f, 0,0])+1/4{(Subscript[f, 1,0]),(Subscript[f, 0,1])}.(RR2-VV);
M[[ind2,3,1]]=2((Subscript[f, 0,0])+1/4{(Subscript[f, 1,0]),(Subscript[f, 0,1])}.(Z-VV))-(Subscript[f, 0,0])+1/12(Z-VV).(H.(Z-VV));
M[[ind2,5,3]]=2((Subscript[f, 0,0])+1/4{(Subscript[f, 1,0]),(Subscript[f, 0,1])}.(RR2-VV))-(Subscript[f, 0,0])+1/12(RR2-VV).(H.(RR2-VV));
M[[ind2,4,2]]=((Subscript[f, 0,0])+1/4{(Subscript[f, 1,0]),(Subscript[f, 0,1])}.(Z-VV))+((Subscript[f, 0,0])+1/4{(Subscript[f, 1,0]),(Subscript[f, 0,1])}.(RR2-VV))-(Subscript[f, 0,0])+1/12(Z-VV).(H.(RR2-VV));
(* ----------- *)
For[s=1,s≤3,s++,
M[[ind1,s+2,s+1]]=L1 M[[ind1,s+2,s]];
M[[ind1,s+2,s+2]]=L1^2 M[[ind1,s+2,s]];
M[[ind2,s+2,s+1]]=L2 M[[ind2,s+2,s]];
M[[ind2,s+2,s+2]]=L2^2 M[[ind2,s+2,s]];
];
(* ----------- *)
M[[ind1,2,1]]=sig M[[ind1,3,1]];
M[[ind2,2,1]]=sig M[[ind1,3,1]];
M[[ind1,2,2]]=L1 M[[ind1,2,1]];
M[[ind2,2,2]]=L2 M[[ind2,2,1]];
(* ----------- *)
M[[1,1,1]]=sig M[[ind1,2,1]];
For[s=2,s≤6,s++,
M[[s,1,1]]=M[[1,1,1]];
];
(* ----------- *)
If[indT==1,
For[s=2,s≤5,s++,
M[[2,s,s]]=M[[1,s,s]];
M[[5,s,s]]=M[[6,s,s]];
];
];
(* ----------- *)
If[indT==2,
For[s=2,s≤5,s++,
M[[1,s,s]]=M[[2,s,s]];
M[[4,s,s]]=M[[3,s,s]];
];
];
(* ----------- *)
If[indT==3,
For[s=2,s≤5,s++,
M[[3,s,s]]=M[[4,s,s]];
M[[6,s,s]]=M[[5,s,s]];
];
];
Return[M/.podF];
]


ClearAll[BaznaFunZaRob]
(* Function returns a list of six matrices with Bezier ordinates that determine an edge basis function indexed by 'indRL'. *)
BaznaFunZaRob[Triangle_,podB_,indRL_,smer_]:=Module[{M,V1,V2,V3,Z,R12,R13,R23,ind1,ind2,ksi,L1,L2,s,sig61,sig23,sig45,sig1,sig2,sig1Z,sig2Z,pom},
M=Table[Table[0,{i,1,5},{j,1,5}],{k,1,6}];
V1=Triangle[[1,1]];
V2=Triangle[[2,1]];
V3=Triangle[[4,1]];
R12=Triangle[[1,2]];
R23=Triangle[[3,2]];
R13=Triangle[[5,2]];
Z=Triangle[[1,3]];
If[indRL==1,ind1=1;ind2=2;ksi=1/(smer.(Z-R12));L1=Baricentricne[V1,V2,V3,R12][[1]];L2=Baricentricne[V1,V2,V3,R12][[2]];sig1=Baricentricne[V1,R13,Z,R12];sig2=Baricentricne[V2,R12,Z,R23];sig1Z=Baricentricne[V1,V2,V3,Z][[1]];sig2Z=Baricentricne[V1,V2,V3,Z][[2]]];
If[indRL==2,ind1=3;ind2=4;ksi=1/(smer.(Z-R23));L1=Baricentricne[V1,V2,V3,R23][[2]];L2=Baricentricne[V1,V2,V3,R23][[3]];sig1=Baricentricne[V2,R12,Z,R23];sig2=Baricentricne[V3,R23,Z,R13];sig1Z=Baricentricne[V1,V2,V3,Z][[2]];sig2Z=Baricentricne[V1,V2,V3,Z][[3]]];
If[indRL==3,ind1=5;ind2=6;ksi=1/(smer.(Z-R13));L1=Baricentricne[V1,V2,V3,R13][[3]];L2=Baricentricne[V1,V2,V3,R13][[1]];sig1=Baricentricne[V3,R23,Z,R13];sig2=Baricentricne[V1,R13,Z,R12];sig1Z=Baricentricne[V1,V2,V3,Z][[3]];sig2Z=Baricentricne[V1,V2,V3,Z][[1]]];
(* ----------- *)
M[[ind1,5,5]]=podB[[1]];
M[[ind2,5,5]]=M[[ind1,5,5]];
M[[ind1,4,4]]=podB[[1]]+1/(4 ksi)podB[[2]];
M[[ind2,4,4]]=M[[ind1,4,4]];
M[[ind1,3,3]]=podB[[1]]+1/(2 ksi)podB[[2]]+1/(12 ksi^2)podB[[3]];
M[[ind2,3,3]]=M[[ind1,3,3]];
(* ------------ *)
For[s=1,s≤3,s++,
M[[ind1,s+2,s+1]]=1/(2L1)M[[ind1,s+2,s+2]];
M[[ind2,s+2,s+1]]=1/(2L2)M[[ind1,s+2,s+2]];
];
(* ------------ *)
M[[ind1,2,1]]=1/(2 L1 sig1[[3]])M[[ind1,3,3]];
M[[ind2,2,1]]=-sig2[[2]]/(2 L2 sig2[[3]])M[[ind2,3,3]];
M[[ind1,2,2]]=L1 M[[ind1,2,1]]+L2 M[[ind2,2,1]];
M[[ind2,2,2]]=M[[ind1,2,2]];
M[[ind1,1,1]]=sig1Z M[[ind1,2,1]]+sig2Z M[[ind2,2,1]];
M[[ind2,1,1]]=M[[ind1,1,1]];
For[s=1,s≤6,s++,M[[s,1,1]]=M[[ind1,1,1]]];
(* ------------ *)
If[indRL==1,M[[6,2,1]]=M[[1,2,1]];M[[6,2,2]]=M[[1,2,1]]Baricentricne[V1,V2,V3,R13][[1]];M[[5,2,2]]=M[[6,2,2]];M[[3,2,1]]=M[[2,2,1]];M[[3,2,2]]=M[[2,2,1]]Baricentricne[V1,V2,V3,R23][[2]];M[[4,2,2]]=M[[3,2,2]]];
If[indRL==2,M[[2,2,1]]=M[[3,2,1]];M[[2,2,2]]=M[[3,2,1]]Baricentricne[V1,V2,V3,R12][[2]];M[[1,2,2]]=M[[2,2,2]];M[[5,2,1]]=M[[4,2,1]];M[[5,2,2]]=M[[4,2,1]]Baricentricne[V1,V2,V3,R13][[3]];M[[6,2,2]]=M[[5,2,2]]];
If[indRL==3,M[[4,2,1]]=M[[5,2,1]];M[[4,2,2]]=M[[5,2,1]]Baricentricne[V1,V2,V3,R23][[3]];M[[3,2,2]]=M[[4,2,2]];M[[1,2,1]]=M[[6,2,1]];M[[1,2,2]]=M[[6,2,1]]Baricentricne[V1,V2,V3,R12][[1]];M[[2,2,2]]=M[[1,2,2]]];
Return[M]
]

ClearAll[DolociBezOrdinateIzKoef]
(* Function computes Bezier ordinates of a spline determined by B-coefficients Cv, Ce. Lists 'SeznamVsehBezOrdinatVert' and 'SeznamVsehBezOrdinatEdge' represent Bezier ordinates of all basis functions. *)
DolociBezOrdinateIzKoef[SeznamCv_,SeznamCe_,SeznamVsehBezOrdinatVert_,SeznamVsehBezOrdinatEdge_,PST_]:=Module[{steviloTrik,sez,i,j,k,kk,i1,i2,SeznamBezOrdinatP,SeznamBezOrdinatE,pomM,SeznamKT,DD,ii,jj},
steviloTrik=Length[PST];
SeznamBezOrdinatP={};
For[i=1,i≤Length[SeznamVsehBezOrdinatVert],i++,
SeznamBezOrdinatP=Append[SeznamBezOrdinatP,SeznamCv[[i]]SeznamVsehBezOrdinatVert[[i]]];
];
SeznamBezOrdinatE={};
For[i=1,i≤Length[SeznamVsehBezOrdinatEdge],i++,
SeznamBezOrdinatE=Append[SeznamBezOrdinatE,SeznamCe[[i]]SeznamVsehBezOrdinatEdge[[i]]];
];
sez=Table[Table[Table[0,{i,1,5},{j,1,5}],{k,1,6}],{kk,1,steviloTrik}];
For[i1=1,i1≤steviloTrik,i1++,
For[i2=1,i2≤6,i2++,
pomM=Table[0,{i,1,5},{j,1,5}];
For[i=1,i≤Length[SeznamBezOrdinatP],i++,
pomM=pomM+SeznamBezOrdinatP[[i,i1,i2]];
] ;
For[i=1,i≤Length[SeznamBezOrdinatE],i++,
pomM=pomM+SeznamBezOrdinatE[[i,i1,i2]];
] ;
sez[[i1,i2]]=pomM;
]];
SeznamKT=Table[Table[Table[0,{i,1,5},{j,1,5}],{k,1,6}],{kk,1,steviloTrik}];
For[k=1,k≤steviloTrik,k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],4];
For[ii=0,ii≤4,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKT[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],sez[[k,kk,ii+1,jj+1]]];
];
];
];
];
Return[{sez,SeznamKT}]
]

ClearAll[IzracunTildeG]
(* Function computes \tilde{g} for each edge (see a paper for more precise explanation). Input data: koefG=values of g, koefV=coefficients at (vertex) basis functions, SeznamVsehBezOrdinat=Bezier ordinates of Bv, TT=triangles of initial triangulation, PST=triangles of splitted triangulation, DetEdge=list of coefficients that determine the edge points, IndeksiTrikRob= a list of indexes of triangles that contain a particular edge, Tsmeri=directions w.   *)
IzracunTildeG[koefG_,koefV_,SeznamVsehBezOrdinat_,TT_,PST_,TV_,DetEdge_,IndeksiTrikRob_,TsmeriV_]:=Module[{koefGb,indR,to1,to2,toR,indto1,indto2,i,j,kt,trik,indPST,pom,pom1,trik1,bb,smer,MM,v11,v12,v21,v22,v31,v32},
koefGb={};
For[indR=1,indR≤Length[TsmeriV],indR++,
{to1,to2,toR}=DetEdge[[indR]];
(* Determine indexes of points contained in a particular edge. *)
indto1=1;
indto2=1;
For[i=1,i≤Length[TV],i++,
If[to1==TV[[i]],indto1=i];
If[to2==TV[[i]],indto2=i];
]; 
(* Determine indexes of triangles that contain the edges. *)
kt=First[IndeksiTrikRob[[indR]]];
trik=TT[[kt]];
(* Determine PS triangles. *)
indPST={};
For[i=1,i≤6,i++,
For[j=1,j≤3,j++, 
If[toR==PST[[kt,i,j]],indPST=Append[indPST,i]];
];
];
(* Six basis functions with respect to a first point. *)
pom=Take[SeznamVsehBezOrdinat,{6(indto1-1)+1,6(indto1)}];
pom1=Take[koefV,{6(indto1-1)+1,6(indto1)}];
trik1=PST[[kt,indPST[[1]]]];
bb=Baricentricne[trik1[[1]],trik1[[2]],trik1[[3]],toR];
smer=Baricentricne[trik1[[1]],trik1[[2]],trik1[[3]],TsmeriV[[indR,2]]]-Baricentricne[trik1[[1]],trik1[[2]],trik1[[3]],TsmeriV[[indR,1]]];
MM=Table[pom[[i,kt,indPST[[1]]]],{i,1,6}];
v11=pom1.Table[Last[Last[Last[Blossom[MM[[i]],{bb,bb,bb,bb}]]]],{i,1,6}];
v21=pom1.(4Table[Last[Last[Last[Blossom[MM[[i]],{smer,bb,bb,bb}]]]],{i,1,6}]);
v31=pom1.(12Table[Last[Last[Last[Blossom[MM[[i]],{smer,smer,bb,bb}]]]],{i,1,6}]);
(* Six basis functions with respect to a second point. *)
pom=Take[SeznamVsehBezOrdinat,{6(indto2-1)+1,6(indto2)}];
pom1=Take[koefV,{6(indto2-1)+1,6(indto2)}];
trik1=PST[[kt,indPST[[1]]]];
bb=Baricentricne[trik1[[1]],trik1[[2]],trik1[[3]],toR];
smer=Baricentricne[trik1[[1]],trik1[[2]],trik1[[3]],TsmeriV[[indR,2]]]-Baricentricne[trik1[[1]],trik1[[2]],trik1[[3]],TsmeriV[[indR,1]]];
MM=Table[pom[[i,kt,indPST[[1]]]],{i,1,6}];
v12=pom1.Table[Last[Last[Last[Blossom[MM[[i]],{bb,bb,bb,bb}]]]],{i,1,6}];
v22=pom1.(4Table[Last[Last[Last[Blossom[MM[[i]],{smer,bb,bb,bb}]]]],{i,1,6}]);
v32=pom1.(12Table[Last[Last[Last[Blossom[MM[[i]],{smer,smer,bb,bb}]]]],{i,1,6}]);
koefGb=Append[koefGb,{v11+v12,v21+v22,v31+v32}];
];
Return[koefG-koefGb]
]

ClearAll[KoefKvaziinterpolanta1]
(* Function returns B-coefficients of a Hermite type quasi-interpolant for a function f. *)
KoefKvaziinterpolanta1[f_,x_,y_,SeznamVsehBezOrdinat_,TT_,PST_,PSTrik_,TV_,DetEdge_,IndeksiTrikRob_,Tsmeri_,TsmeriV_,PSLines_]:=Module[{Df,Hf,koefCV,k,c1,c2,c3,c4,c5,c6,koefGa,i,koefG,koefCE,indR,zeta1,zeta2,zeta3,lam1,lam2,lam3,CC,norm,pom},
Df={D[f,x],D[f,y]};
Hf={{D[f,{x,2}],D[D[f,x],y]},{D[D[f,x],y],D[f,{y,2}]}};
koefCV={};
For[k=1,k≤Length[TV],k++,
c1=f+Df.(PSTrik[[k,1]]-TV[[k]])+1/3(PSTrik[[k,1]]-TV[[k]]).Hf.(PSTrik[[k,1]]-TV[[k]])/.{x->TV[[k,1]],y->TV[[k,2]]};
c2=f+1/2Df.(PSTrik[[k,1]]+PSTrik[[k,2]]-2TV[[k]])+1/3(PSTrik[[k,1]]-TV[[k]]).Hf.(PSTrik[[k,2]]-TV[[k]])/.{x->TV[[k,1]],y->TV[[k,2]]};
c3=f+1/2Df.(PSTrik[[k,1]]+PSTrik[[k,3]]-2TV[[k]])+1/3(PSTrik[[k,1]]-TV[[k]]).Hf.(PSTrik[[k,3]]-TV[[k]])/.{x->TV[[k,1]],y->TV[[k,2]]};
c4=f+Df.(PSTrik[[k,2]]-TV[[k]])+1/3(PSTrik[[k,2]]-TV[[k]]).Hf.(PSTrik[[k,2]]-TV[[k]])/.{x->TV[[k,1]],y->TV[[k,2]]};
c5=f+1/2Df.(PSTrik[[k,2]]+PSTrik[[k,3]]-2TV[[k]])+1/3(PSTrik[[k,2]]-TV[[k]]).Hf.(PSTrik[[k,3]]-TV[[k]])/.{x->TV[[k,1]],y->TV[[k,2]]};
c6=f+Df.(PSTrik[[k,3]]-TV[[k]])+1/3(PSTrik[[k,3]]-TV[[k]]).Hf.(PSTrik[[k,3]]-TV[[k]])/.{x->TV[[k,1]],y->TV[[k,2]]};
koefCV=Append[koefCV,{c1,c2,c3,c4,c5,c6}];
];
koefCV=Flatten[koefCV];
(* ---- *)
koefGa={};
For[i=1,i≤Length[Tsmeri],i++,
koefGa=Append[koefGa,{f,Df.Tsmeri[[i]],Tsmeri[[i]].Hf.Tsmeri[[i]]}/.{x->DetEdge[[i,3,1]],y->DetEdge[[i,3,2]]}]
];
koefG=IzracunTildeG[koefGa,koefCV,SeznamVsehBezOrdinat,TT,PST,TV,DetEdge,IndeksiTrikRob,TsmeriV];
koefCE={};
For[indR=1,indR≤Length[DetEdge],indR++,{zeta1,zeta2,zeta3}=Baricentricne[PSLines[[indR,1]],PSLines[[indR,2]],First[DetEdge[[indR]]],Last[DetEdge[[indR]]]];
{lam1,lam2,lam3}=Baricentricne[DetEdge[[indR,1]],DetEdge[[indR,2]],PSLines[[indR,2]],DetEdge[[indR,3]]];
CC=2lam1 lam2;
norm=Norm[PSLines[[indR,2]]-PSLines[[indR,1]]];
pom=1/CC{koefG[[indR,1]]-koefG[[indR,2]]zeta2 norm+1/3koefG[[indR,3]]zeta2^2norm^2,koefG[[indR,1]]+1/2koefG[[indR,2]](zeta1-zeta2) norm-1/3koefG[[indR,3]]zeta1 zeta2 norm^2,koefG[[indR,1]]+koefG[[indR,2]]zeta1 norm+1/3koefG[[indR,3]]zeta1^2norm^2};
koefCE=Append[koefCE,pom]
];
koefCE=Flatten[koefCE];
Return[{koefCV,koefCE}];
]

ClearAll[LagPol]
(* Function returns Lagrange interpolating polynomial that interpolates function 'ff' at 'toc'. *)
LagPol[toc_,ff_,x_,y_]:=Module[{dd,mat,tabLag,dmat,res,mat1,kk,k,vv,pol},
dd=IdentityMatrix[6];
mat=Table[{1,x,y,x^2,x y,y^2}/.{x->toc[[i,1]],y->toc[[i,2]]},{i,1,6}];
tabLag={};
For[kk=1,kk≤6,kk++,
vv=dd[[kk]];
dmat=Det[mat];
res={};
For[k=1,k≤6,k++,
mat1=Transpose[mat];
mat1[[k]]=vv;
res=Append[res,Det[Transpose[mat1]]/dmat];
];
tabLag=Append[tabLag,res.{1,x,y,x^2,x y,y^2}];
];
pol=Simplify[Table[ff/.{x->toc[[i,1]],y->toc[[i,2]]},{i,1,6}].tabLag];
Return[pol]
]

ClearAll[LagPolDegN]
(* Function returns Lagrange interpolating polynomial of degree n, that interpolates function 'ff' at a uniform lattice of a triangle <V1,V2,V3>. *)
LagPolDegN[V1_,V2_,V3_,n_,ff_,x_,y_]:=Module[{toc,vred1,vred2,koef,a,baza, res, IntPol,i,j,pol},
toc=Flatten[Table[(n-i-j)/n V1+ j/n V2+i/n V3,{i,0,n},{j,0,n-i}],1];
vred1=Table[ff/.{x->toc[[i,1]],y->toc[[i,2]]},{i,1,Length[toc]}];
baza=Flatten[Table[x^i y^j,{i,0,n},{j,0,n-i}]];
koef=Flatten[Table[(Subscript[a, i,j]),{i,0,n},{j,0,n-i}]];
pol=koef.baza;
vred2=Table[pol/.{x->toc[[i,1]],y->toc[[i,2]]},{i,1,Length[toc]}];
res=Solve[vred1==vred2,koef];
IntPol=pol/.res[[1]];
Return[IntPol];
]

ClearAll[KoefKvaziinterpolanta2]
(* Function returns B-coefficients of a Lagrange type quasi-interpolant for a function f. Local polynomials of degree 'stopnja' are used. *)
KoefKvaziinterpolanta2[f_,stopnja_,x_,y_,SeznamVsehBezOrdinat_,TT_,PST_,PSTrik_,TV_,DetEdge_,IndeksiTrikRob_,Tsmeri_,TsmeriV_,PSLines_]:=Module[{IntPolinomi,fp,Dfp,Hfp,koefCV,k,c1,c2,c3,c4,c5,c6,koefGa,i,koefG,koefCE,indR,zeta1,zeta2,zeta3,lam1,lam2,lam3,CC,norm,pom,kat,V1,V2,V3,vec,vec1,intp},
IntPolinomi={};
For[i=1,i≤Length[TV],i++,
IntPolinomi=Append[IntPolinomi,LagPolDegN[PSTrik[[i,1]],PSTrik[[i,2]],PSTrik[[i,3]],stopnja,f,x,y]];
];
koefCV={};
For[k=1,k≤Length[TV],k++,
fp=IntPolinomi[[k]];
Dfp={D[fp,x],D[fp,y]};
Hfp={{D[fp,{x,2}],D[D[fp,x],y]},{D[D[fp,x],y],D[fp,{y,2}]}};
c1=fp+Dfp.(PSTrik[[k,1]]-TV[[k]])+1/3(PSTrik[[k,1]]-TV[[k]]).Hfp.(PSTrik[[k,1]]-TV[[k]])/.{x->TV[[k,1]],y->TV[[k,2]]};
c2=fp+1/2Dfp.(PSTrik[[k,1]]+PSTrik[[k,2]]-2TV[[k]])+1/3(PSTrik[[k,1]]-TV[[k]]).Hfp.(PSTrik[[k,2]]-TV[[k]])/.{x->TV[[k,1]],y->TV[[k,2]]};
c3=fp+1/2Dfp.(PSTrik[[k,1]]+PSTrik[[k,3]]-2TV[[k]])+1/3(PSTrik[[k,1]]-TV[[k]]).Hfp.(PSTrik[[k,3]]-TV[[k]])/.{x->TV[[k,1]],y->TV[[k,2]]};
c4=fp+Dfp.(PSTrik[[k,2]]-TV[[k]])+1/3(PSTrik[[k,2]]-TV[[k]]).Hfp.(PSTrik[[k,2]]-TV[[k]])/.{x->TV[[k,1]],y->TV[[k,2]]};
c5=fp+1/2Dfp.(PSTrik[[k,2]]+PSTrik[[k,3]]-2TV[[k]])+1/3(PSTrik[[k,2]]-TV[[k]]).Hfp.(PSTrik[[k,3]]-TV[[k]])/.{x->TV[[k,1]],y->TV[[k,2]]};
c6=fp+Dfp.(PSTrik[[k,3]]-TV[[k]])+1/3(PSTrik[[k,3]]-TV[[k]]).Hfp.(PSTrik[[k,3]]-TV[[k]])/.{x->TV[[k,1]],y->TV[[k,2]]};
koefCV=Append[koefCV,{c1,c2,c3,c4,c5,c6}];
];
koefCV=Flatten[koefCV];
(* ---- *)
(* edge coefficients *)
koefGa={};
For[kat=1,kat≤Length[Tsmeri],kat++,
V1=PSLines[[kat,1]];
V2=PSLines[[kat,2]];
vec=V2-V1;
vec1={vec[[2]],-vec[[1]]};
V3=V1+vec1;
intp=LagPolDegN[V1,V2,V3,stopnja,f,x,y];
fp=intp;
Dfp={D[fp,x],D[fp,y]};
Hfp={{D[fp,{x,2}],D[D[fp,x],y]},{D[D[fp,x],y],D[fp,{y,2}]}};
koefGa=Append[koefGa,{intp,Dfp.Tsmeri[[kat]],Tsmeri[[kat]].Hfp.Tsmeri[[kat]]}/.{x->DetEdge[[kat,3,1]],y->DetEdge[[kat,3,2]]}];
];
(* coefficients tilde{g}:*)
koefG=IzracunTildeG[koefGa,koefCV,SeznamVsehBezOrdinat,TT,PST,TV,DetEdge,IndeksiTrikRob,TsmeriV];
koefCE={};
For[indR=1,indR≤Length[DetEdge],indR++,{zeta1,zeta2,zeta3}=Baricentricne[PSLines[[indR,1]],PSLines[[indR,2]],First[DetEdge[[indR]]],Last[DetEdge[[indR]]]];
{lam1,lam2,lam3}=Baricentricne[DetEdge[[indR,1]],DetEdge[[indR,2]],PSLines[[indR,2]],DetEdge[[indR,3]]];
CC=2lam1 lam2;
norm=Norm[PSLines[[indR,2]]-PSLines[[indR,1]]];
pom=1/CC{koefG[[indR,1]]-koefG[[indR,2]]zeta2 norm+1/3koefG[[indR,3]]zeta2^2norm^2,koefG[[indR,1]]+1/2koefG[[indR,2]](zeta1-zeta2) norm-1/3koefG[[indR,3]]zeta1 zeta2 norm^2,koefG[[indR,1]]+koefG[[indR,2]]zeta1 norm+1/3koefG[[indR,3]]zeta1^2norm^2};
koefCE=Append[koefCE,pom];
];
koefCE=Flatten[koefCE];
Return[{koefCV,koefCE}];
]

ClearAll[Incenters]
(* Function returns incenters of a triangle with vertices V1,V2,V3.  *)
Incenters[V1_,V2_,V3_]:=Module[{bar,a,b,c},
a=Norm[V2-V3];
b=Norm[V1-V3];
c=Norm[V2-V1];
Return[(a V1+b V2+c V3)/(a+b+c)]
]

ClearAll[VsebujeRob]
(* Function returns 1, if a triangle 'Trik' contains an edge 'Edge', else it returns 0. *)
VsebujeRob[Trik_,Edge_]:=Module[{},
If[({Trik[[1]],Trik[[2]]}==Edge)||({Trik[[2]],Trik[[1]]}==Edge)||({Trik[[1]],Trik[[3]]}==Edge)||({Trik[[3]],Trik[[1]]}==Edge)||({Trik[[2]],Trik[[3]]}==Edge)||({Trik[[3]],Trik[[2]]}==Edge),Return[1],Return[0]];
]
ClearAll[VsebujeTocko]
(* Function returns 1, if a triangle 'Trik' contains a point 'tocka', else it returns 0. *)
VsebujeTocko[Trik_,tocka_]:=Module[{},
If[(Trik[[1]]==tocka)||(Trik[[2]]==tocka)||(Trik[[3]]==tocka),Return[1],Return[0]];
]

ClearAll[DolociRij]
(* Function returns points Rij for each edge of a triangulation. *)
DolociRij[TE_,TZ_,IndeksiTrikRob_]:=Module[{s,t,k,Rij,TR},
TR={};
For[k=1,k≤Length[TE],k++,
If[Length[IndeksiTrikRob[[k]]]==1,Rij=(TE[[k,1]]+TE[[k,2]])/2,
Rij=((TE[[k,1]]+s(TE[[k,2]]-TE[[k,1]]))/.Flatten[Solve[TZ[[IndeksiTrikRob[[k,1]]]]+t(TZ[[IndeksiTrikRob[[k,2]]]]-TZ[[IndeksiTrikRob[[k,1]]]])==TE[[k,1]]+s(TE[[k,2]]-TE[[k,1]]),{s,t}]]);
];
TR=Append[TR,Rij];
];
Return[TR]
]


ClearAll[DolociPSTrik1]
(* Function returns PS triangles for each point of a triangulation. NOTE THAT THESE TRIANGLES ARE NOT OPTIMAL!!! *)
DolociPSTrik1[TV1_,Toc_]:=Module[{n,d,ind,j,d1,toc},
d=0;
ind=0;
n=Length[Toc];
For[j=1,j≤n,j++,
d1=Norm[Toc[[j]]-TV1];
If[N[d1,100]≥N[d,100],ind=j;d=d1];
];
toc={{3TV1[[1]]-2Toc[[ind,1]],3TV1[[2]]-2Toc[[ind,2]]},{Toc[[ind,1]]+Sqrt[3](TV1[[2]]-Toc[[ind,2]]),Toc[[ind,2]]-Sqrt[3](TV1[[1]]-Toc[[ind,1]])},{Toc[[ind,1]]-Sqrt[3](TV1[[2]]-Toc[[ind,2]]),Toc[[ind,2]]+Sqrt[3](TV1[[1]]-Toc[[ind,1]])}};
Return[toc]
]

ClearAll[DolociPSTrikTip1]
(*  Function returns PS triangles for each point of a Type 1 triangulation. *)
DolociPSTrikTip1[TV1_,h_]:=Module[{toc},
toc ={
{TV1[[1]]-3h/(2(2+Sqrt[2])),TV1[[2]]+(1-Sqrt[2])h/(2(2+Sqrt[2]))},
{TV1[[1]]+5h/(2(2+Sqrt[2])),TV1[[2]]+(1+3Sqrt[2])h/(2(2+Sqrt[2]))},
{TV1[[1]]+(2Sqrt[2]-3)h/(2(2+Sqrt[2])),TV1[[2]]+(1-3Sqrt[2])h/(2(2+Sqrt[2]))}
};
Return[toc]
]


ClearAll[BazneFunkcijeZaTocke]
(* Function returns basis functions associated to points (list of B-coefficients) for a triangulation TT with vertices TT and edges TE.  *)
BazneFunkcijeZaTocke[TV_,TT_,TE_]:=Module[{TZ,IndeksiTrik,pom,i,j,IndeksiTrikRob,IndRobTrik,IndeksiRobV,TRE,TR,PST,TPS,PSTrik,SeznamAlf,SeznamVsehBezOrdinat,SeznamVsehKT,indT,indBF,SeznamBezOrdinat,ss,ind,indTLok,pod,kt,MF,SeznamKT,k,kk,DD,ii,jj},
TZ={};
For[i=1,i≤Length[TT],i++,
TZ=Append[TZ,Incenters[TT[[i,1]],TT[[i,2]],TT[[i,3]]]]
];
IndeksiTrik={};
For[i=1,i≤Length[TV],i++,
pom={};
For[j=1,j≤Length[TT],j++,
If[VsebujeTocko[TT[[j]],TV[[i]]]==1,pom=Append[pom,j]];
];
IndeksiTrik=Append[IndeksiTrik,pom];
];
IndeksiTrikRob={};
For[i=1,i≤Length[TE],i++,
pom={};
For[j=1,j≤Length[TT],j++,
If[VsebujeRob[TT[[j]],TE[[i]]]==1,pom=Append[pom,j]];
];
IndeksiTrikRob=Append[IndeksiTrikRob,pom];
];
IndRobTrik={};
For[i=1,i≤Length[TT],i++,
pom={0,0,0};
For[j=1,j≤Length[TE],j++,
If[({TT[[i,1]],TT[[i,2]]}==TE[[j]])||({TT[[i,2]],TT[[i,1]]}==TE[[j]]),pom[[1]]=j];
If[({TT[[i,2]],TT[[i,3]]}==TE[[j]])||({TT[[i,3]],TT[[i,2]]}==TE[[j]]),pom[[2]]=j];
If[({TT[[i,1]],TT[[i,3]]}==TE[[j]])||({TT[[i,3]],TT[[i,1]]}==TE[[j]]),pom[[3]]=j];
];
IndRobTrik=Append[IndRobTrik,pom];
];
IndeksiRobV={};
For[i=1,i≤Length[TV],i++,
pom={};
For[j=1,j≤Length[TE],j++,
If[(TV[[i]]==TE[[j,1]])||(TV[[i]]==TE[[j,2]]),pom=Append[pom,j]];
];
IndeksiRobV=Append[IndeksiRobV,pom];
];
TRE=DolociRij[TE,TZ,IndeksiTrikRob];
TR={};
For[i=1,i≤Length[TT],i++,
TR=Append[TR,{TRE[[IndRobTrik[[i,1]]]],TRE[[IndRobTrik[[i,2]]]],TRE[[IndRobTrik[[i,3]]]]}];
];
PST=Table[DeterminePSTriangles[TT[[k,1]],TT[[k,2]],TT[[k,3]],TR[[k,1]],TR[[k,2]],TR[[k,3]],TZ[[k]]],{k,1,Length[TT]}];
TPS=TV;
For[i=1,i≤Length[TV],i++,
TPS[[i]]=Join[Table[1/2(TV[[i]]+TRE[[IndeksiRobV[[i,j]]]]),{j,1,Length[IndeksiRobV[[i]]]}],Table[1/2(TV[[i]]+TZ[[IndeksiTrik[[i,j]]]]),{j,1,Length[IndeksiTrik[[i]]]}]];
];
(*If[PSTrik1=={},PSTrik=Table[DolociPSTrik1[TV[[i]],TPS[[i]]],{i,1,Length[TV]}],PSTrik=PSTrik1];*)
PSTrik=Table[DolociPSTrik1[TV[[i]],TPS[[i]]],{i,1,Length[TV]}];
SeznamAlf=Table[Transpose[DolociParAlpha[PSTrik[[indT,1]],PSTrik[[indT,2]],PSTrik[[indT,3]],TV[[indT]]]],{indT,1,Length[TV]}];
(* ------- *)
SeznamVsehBezOrdinat={};
SeznamVsehKT={};
For[indT=1,indT<=Length[TV],indT++,
For[indBF=1,indBF≤6,indBF++,
SeznamBezOrdinat=Table[Table[Table[0,{i,1,5},{j,1,5}],{k,1,6}],{kk,1,Length[TT]}];
For[ss=1,ss≤Length[IndeksiTrik[[indT]]],ss++,
pod=SeznamAlf[[indT,indBF]];
(* which triangle *)
kt=IndeksiTrik[[indT,ss]];
(* Index of a point in this triangle (localy). *)
If[TT[[kt,1]]==TV[[indT]],indTLok=1];
If[TT[[kt,2]]==TV[[indT]],indTLok=2];
If[TT[[kt,3]]==TV[[indT]],indTLok=3];
(* Bezier ordinates: *)
MF=BaznaFunZaTocko[PST[[kt]],pod,indTLok];
SeznamBezOrdinat[[kt]]=MF;
];
SeznamKT=Table[Table[Table[0,{i,1,5},{j,1,5}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],4];
For[ii=0,ii≤4,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKT[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],SeznamBezOrdinat[[k,kk,ii+1,jj+1]]];
];
];
];
];
SeznamVsehBezOrdinat=Append[SeznamVsehBezOrdinat,SeznamBezOrdinat];
SeznamVsehKT=Append[SeznamVsehKT,SeznamKT];
];
];
Return[{SeznamVsehBezOrdinat,SeznamVsehKT}]
]

ClearAll[BazneFunkcijeZaRob]
(* Function returns basis functions associated to edges (list of B-coefficients) for a triangulation TT with vertices TT and edges TE.  *)
BazneFunkcijeZaRob[TV_,TT_,TE_]:=Module[{TZ,IndeksiTrik,pom,i,j,IndeksiTrikRob,IndRobTrik,IndeksiRobV,TRE,TR,PST,TPS,SeznamBet,SeznamVsehBezOrdinatRob,SeznamVsehKTRob,indR,indBF,DetEdge,PSLines,SeznamBezOrdinat,SeznamBezOrdinatRob,toc1,Tsmeri,toc2,trik,s,indRL,ind,ind1,ind2,ss,indTLok,pod,kt,MF,SeznamKT,k,kk,DD,ii,jj},
TZ={};
For[i=1,i≤Length[TT],i++,
TZ=Append[TZ,Incenters[TT[[i,1]],TT[[i,2]],TT[[i,3]]]]
];
IndeksiTrik={};
For[i=1,i≤Length[TV],i++,
pom={};
For[j=1,j≤Length[TT],j++,
If[VsebujeTocko[TT[[j]],TV[[i]]]==1,pom=Append[pom,j]];
];
IndeksiTrik=Append[IndeksiTrik,pom];
];
IndeksiTrikRob={};
For[i=1,i≤Length[TE],i++,
pom={};
For[j=1,j≤Length[TT],j++,
If[VsebujeRob[TT[[j]],TE[[i]]]==1,pom=Append[pom,j]];
];
IndeksiTrikRob=Append[IndeksiTrikRob,pom];
];
IndRobTrik={};
For[i=1,i≤Length[TT],i++,
pom={0,0,0};
For[j=1,j≤Length[TE],j++,
If[({TT[[i,1]],TT[[i,2]]}==TE[[j]])||({TT[[i,2]],TT[[i,1]]}==TE[[j]]),pom[[1]]=j];
If[({TT[[i,2]],TT[[i,3]]}==TE[[j]])||({TT[[i,3]],TT[[i,2]]}==TE[[j]]),pom[[2]]=j];
If[({TT[[i,1]],TT[[i,3]]}==TE[[j]])||({TT[[i,3]],TT[[i,1]]}==TE[[j]]),pom[[3]]=j];
];
IndRobTrik=Append[IndRobTrik,pom];
];
IndeksiRobV={};
For[i=1,i≤Length[TV],i++,
pom={};
For[j=1,j≤Length[TE],j++,
If[(TV[[i]]==TE[[j,1]])||(TV[[i]]==TE[[j,2]]),pom=Append[pom,j]];
];
IndeksiRobV=Append[IndeksiRobV,pom];
];
TRE=DolociRij[TE,TZ,IndeksiTrikRob];
TR={};
For[i=1,i≤Length[TT],i++,
TR=Append[TR,{TRE[[IndRobTrik[[i,1]]]],TRE[[IndRobTrik[[i,2]]]],TRE[[IndRobTrik[[i,3]]]]}];
];
PST=Table[DeterminePSTriangles[TT[[k,1]],TT[[k,2]],TT[[k,3]],TR[[k,1]],TR[[k,2]],TR[[k,3]],TZ[[k]]],{k,1,Length[TT]}];
DetEdge=TE;
For[i=1,i≤Length[TE],i++,
DetEdge[[i]]=Append[TE[[i]],TRE[[i]]]
];
PSLines=TE;
For[i=1,i≤Length[TE],i++,
If[Length[IndeksiTrikRob[[i]]]==1,PSLines[[i]]={TRE[[i]],1/2(TRE[[i]] + TZ[[IndeksiTrikRob[[i,1]]]])},
PSLines[[i]]={1/2(TRE[[i]] + TZ[[IndeksiTrikRob[[i,1]]]]),1/2(TRE[[i]] + TZ[[IndeksiTrikRob[[i,2]]]])}
];
];
SeznamBet=Table[Transpose[DolociParBeta[DetEdge[[indR,1]],DetEdge[[indR,2]],DetEdge[[indR,3]],PSLines[[indR,1]],PSLines[[indR,2]]]],{indR,1,Length[DetEdge]}];
Tsmeri=Table[(PSLines[[k,2]]-PSLines[[k,1]])/Norm[(PSLines[[k,2]]-PSLines[[k,1]])],{k,1,Length[PSLines]}];
(* ------- *)
SeznamVsehBezOrdinatRob={};
SeznamVsehKTRob={};
For[indR=1,indR<=Length[DetEdge],indR++,
For[indBF=1,indBF≤3,indBF++,
SeznamBezOrdinatRob=Table[Table[Table[0,{i,1,5},{j,1,5}],{k,1,6}],{kk,1,Length[TT]}];
For[ss=1,ss≤Length[IndeksiTrikRob[[indR]]],ss++,
pod=SeznamBet[[indR,indBF]];
(* which triangle *)
kt=IndeksiTrikRob[[indR,ss]];
(* Index of an edge in this triangle (localy). *)
trik=TT[[kt]];
toc1=DetEdge[[indR]][[1]];
toc2=DetEdge[[indR]][[2]];
For[s=1,s≤3,s++,
If[toc1==trik[[s]],ind1=s];
If[toc2==trik[[s]],ind2=s];
];
ind=Sort[{ind1,ind2}];
If[ind=={1,2},indRL=1];
If[ind=={2,3},indRL=2];
If[ind=={1,3},indRL=3];
(* Bezier ordinates:*)
MF=BaznaFunZaRob[PST[[kt]],pod,indRL,Tsmeri[[indR]]];
SeznamBezOrdinatRob[[kt]]=MF;
];
SeznamKT=Table[Table[Table[0,{i,1,5},{j,1,5}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],4];
For[ii=0,ii≤4,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKT[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],SeznamBezOrdinatRob[[k,kk,ii+1,jj+1]]];
];
];
];
];
SeznamVsehBezOrdinatRob=Append[SeznamVsehBezOrdinatRob,SeznamBezOrdinatRob];
SeznamVsehKTRob=Append[SeznamVsehKTRob,SeznamKT];
]
];
Return[{SeznamVsehBezOrdinatRob,SeznamVsehKTRob}]
]

ClearAll[IzrKolicin]
(* Function returns a list of data needed to construct basis functions from a list of data points, edges and triangles. *)IzrKolicin[TV_,TT_,TE_]:=Module[{TZ,i,j,k,IndeksiTrik,pom,IndeksiTrikRob,IndRobTrik,IndeksiRobV,TRE,TR,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV},
TZ={};
For[i=1,i≤Length[TT],i++,
TZ=Append[TZ,Incenters[TT[[i,1]],TT[[i,2]],TT[[i,3]]]]
];
IndeksiTrik={};
For[i=1,i≤Length[TV],i++,
pom={};
For[j=1,j≤Length[TT],j++,
If[VsebujeTocko[TT[[j]],TV[[i]]]==1,pom=Append[pom,j]];
];
IndeksiTrik=Append[IndeksiTrik,pom];
];
IndeksiTrikRob={};
For[i=1,i≤Length[TE],i++,
pom={};
For[j=1,j≤Length[TT],j++,
If[VsebujeRob[TT[[j]],TE[[i]]]==1,pom=Append[pom,j]];
];
IndeksiTrikRob=Append[IndeksiTrikRob,pom];
];
IndRobTrik={};
For[i=1,i≤Length[TT],i++,
pom={0,0,0};
For[j=1,j≤Length[TE],j++,
If[({TT[[i,1]],TT[[i,2]]}==TE[[j]])||({TT[[i,2]],TT[[i,1]]}==TE[[j]]),pom[[1]]=j];
If[({TT[[i,2]],TT[[i,3]]}==TE[[j]])||({TT[[i,3]],TT[[i,2]]}==TE[[j]]),pom[[2]]=j];
If[({TT[[i,1]],TT[[i,3]]}==TE[[j]])||({TT[[i,3]],TT[[i,1]]}==TE[[j]]),pom[[3]]=j];
];
IndRobTrik=Append[IndRobTrik,pom];
];
IndeksiRobV={};
For[i=1,i≤Length[TV],i++,
pom={};
For[j=1,j≤Length[TE],j++,
If[(TV[[i]]==TE[[j,1]])||(TV[[i]]==TE[[j,2]]),pom=Append[pom,j]];
];
IndeksiRobV=Append[IndeksiRobV,pom];
];
TRE=DolociRij[TE,TZ,IndeksiTrikRob];
TR={};
For[i=1,i≤Length[TT],i++,
TR=Append[TR,{TRE[[IndRobTrik[[i,1]]]],TRE[[IndRobTrik[[i,2]]]],TRE[[IndRobTrik[[i,3]]]]}];
];
PST=Table[DeterminePSTriangles[TT[[k,1]],TT[[k,2]],TT[[k,3]],TR[[k,1]],TR[[k,2]],TR[[k,3]],TZ[[k]]],{k,1,Length[TT]}];
TPS=TV;
For[i=1,i≤Length[TV],i++,
TPS[[i]]=Join[Table[1/2(TV[[i]]+TRE[[IndeksiRobV[[i,j]]]]),{j,1,Length[IndeksiRobV[[i]]]}],Table[1/2(TV[[i]]+TZ[[IndeksiTrik[[i,j]]]]),{j,1,Length[IndeksiTrik[[i]]]}]];
];
PSTrik=Table[DolociPSTrik1[TV[[i]],TPS[[i]]],{i,1,Length[TV]}];
DetEdge=TE;
For[i=1,i≤Length[TE],i++,
DetEdge[[i]]=Append[TE[[i]],TRE[[i]]]
];
PSLines=TE;
For[i=1,i≤Length[TE],i++,
If[Length[IndeksiTrikRob[[i]]]==1,PSLines[[i]]={TRE[[i]],1/2(TRE[[i]] + TZ[[IndeksiTrikRob[[i,1]]]])},
PSLines[[i]]={1/2(TRE[[i]] + TZ[[IndeksiTrikRob[[i,1]]]]),1/2(TRE[[i]] + TZ[[IndeksiTrikRob[[i,2]]]])}
];
];
Tsmeri=Table[(PSLines[[k,2]]-PSLines[[k,1]])/Norm[(PSLines[[k,2]]-PSLines[[k,1]])],{k,1,Length[PSLines]}];
TsmeriV=Table[{PSLines[[k,1]],PSLines[[k,1]]+(PSLines[[k,2]]-PSLines[[k,1]])/Norm[(PSLines[[k,2]]-PSLines[[k,1]])]},{k,1,Length[PSLines]}];
Return[{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}]
]


ClearAll[IzrKolicinTip1]
(* Function returns a list of data needed to construct basis functions from a list of data points, edges and triangles. *)
IzrKolicinTip1[TV_,TT_,TE_,h_]:=Module[{TZ,i,j,k,IndeksiTrik,pom,IndeksiTrikRob,IndRobTrik,IndeksiRobV,TRE,TR,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV},
TZ={};
For[i=1,i≤Length[TT],i++,
TZ=Append[TZ,Incenters[TT[[i,1]],TT[[i,2]],TT[[i,3]]]]
];
IndeksiTrik={};
For[i=1,i≤Length[TV],i++,
pom={};
For[j=1,j≤Length[TT],j++,
If[VsebujeTocko[TT[[j]],TV[[i]]]==1,pom=Append[pom,j]];
];
IndeksiTrik=Append[IndeksiTrik,pom];
];
IndeksiTrikRob={};
For[i=1,i≤Length[TE],i++,
pom={};
For[j=1,j≤Length[TT],j++,
If[VsebujeRob[TT[[j]],TE[[i]]]==1,pom=Append[pom,j]];
];
IndeksiTrikRob=Append[IndeksiTrikRob,pom];
];
IndRobTrik={};
For[i=1,i≤Length[TT],i++,
pom={0,0,0};
For[j=1,j≤Length[TE],j++,
If[({TT[[i,1]],TT[[i,2]]}==TE[[j]])||({TT[[i,2]],TT[[i,1]]}==TE[[j]]),pom[[1]]=j];
If[({TT[[i,2]],TT[[i,3]]}==TE[[j]])||({TT[[i,3]],TT[[i,2]]}==TE[[j]]),pom[[2]]=j];
If[({TT[[i,1]],TT[[i,3]]}==TE[[j]])||({TT[[i,3]],TT[[i,1]]}==TE[[j]]),pom[[3]]=j];
];
IndRobTrik=Append[IndRobTrik,pom];
];
IndeksiRobV={};
For[i=1,i≤Length[TV],i++,
pom={};
For[j=1,j≤Length[TE],j++,
If[(TV[[i]]==TE[[j,1]])||(TV[[i]]==TE[[j,2]]),pom=Append[pom,j]];
];
IndeksiRobV=Append[IndeksiRobV,pom];
];
TRE=DolociRij[TE,TZ,IndeksiTrikRob];
TR={};
For[i=1,i≤Length[TT],i++,
TR=Append[TR,{TRE[[IndRobTrik[[i,1]]]],TRE[[IndRobTrik[[i,2]]]],TRE[[IndRobTrik[[i,3]]]]}];
];
PST=Table[DeterminePSTriangles[TT[[k,1]],TT[[k,2]],TT[[k,3]],TR[[k,1]],TR[[k,2]],TR[[k,3]],TZ[[k]]],{k,1,Length[TT]}];
TPS=TV;
For[i=1,i≤Length[TV],i++,
TPS[[i]]=Join[Table[1/2(TV[[i]]+TRE[[IndeksiRobV[[i,j]]]]),{j,1,Length[IndeksiRobV[[i]]]}],Table[1/2(TV[[i]]+TZ[[IndeksiTrik[[i,j]]]]),{j,1,Length[IndeksiTrik[[i]]]}]];
];
PSTrik=Table[DolociPSTrikTip1[TV[[i]],h],{i,1,Length[TV]}];
DetEdge=TE;
For[i=1,i≤Length[TE],i++,
DetEdge[[i]]=Append[TE[[i]],TRE[[i]]]
];
PSLines=TE;
For[i=1,i≤Length[TE],i++,
If[Length[IndeksiTrikRob[[i]]]==1,PSLines[[i]]={TRE[[i]],1/2(TRE[[i]] + TZ[[IndeksiTrikRob[[i,1]]]])},
PSLines[[i]]={1/2(TRE[[i]] + TZ[[IndeksiTrikRob[[i,1]]]]),1/2(TRE[[i]] + TZ[[IndeksiTrikRob[[i,2]]]])}
];
];
Tsmeri=Table[(PSLines[[k,2]]-PSLines[[k,1]])/Norm[(PSLines[[k,2]]-PSLines[[k,1]])],{k,1,Length[PSLines]}];
TsmeriV=Table[{PSLines[[k,1]],PSLines[[k,1]]+(PSLines[[k,2]]-PSLines[[k,1]])/Norm[(PSLines[[k,2]]-PSLines[[k,1]])]},{k,1,Length[PSLines]}];
Return[{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}]
]


ClearAll[VseBazneFunkcijeZaQuarticSpline]
(* Function returns basis B splines of degree 4 associated to points and edges.  *)
VseBazneFunkcijeZaQuarticSpline[TV_,TT_,TE_,PSTrik1_]:=Module[{TZ,IndeksiTrik,pom,i,j,IndeksiTrikRob,IndRobTrik,IndeksiRobV,TRE,TR,PST,TPS,PSTrik,TsmeriV,SeznamBet,SeznamVsehBezOrdinatRob,SeznamVsehKTRob,indR,indBF,DetEdge,PSLines,SeznamBezOrdinat,SeznamBezOrdinatRob,toc1,Tsmeri,toc2,trik,s,indRL,ind,ind1,ind2,ss,indTLok,pod,kt,MF,SeznamKT,k,kk,DD,ii,jj,SeznamAlf,SeznamVsehBezOrdinat,SeznamVsehKT,indT},
{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}=IzrKolicin[TV,TT,TE];
(* --------------- basis functions associated to points ------------------ *)
If[PSTrik1=={},PSTrik=Table[DolociPSTrik1[TV[[i]],TPS[[i]]],{i,1,Length[TV]}],PSTrik=PSTrik1];
(*PSTrik=Table[DolociPSTrik1[TV[[i]],TPS[[i]]],{i,1,Length[TV]}];*)
SeznamAlf=Table[Transpose[DolociParAlpha[PSTrik[[indT,1]],PSTrik[[indT,2]],PSTrik[[indT,3]],TV[[indT]]]],{indT,1,Length[TV]}];
(* ------- *)
SeznamVsehBezOrdinat={};
SeznamVsehKT={};
For[indT=1,indT<=Length[TV],indT++,
For[indBF=1,indBF≤6,indBF++,
SeznamBezOrdinat=Table[Table[Table[0,{i,1,5},{j,1,5}],{k,1,6}],{kk,1,Length[TT]}];
For[ss=1,ss≤Length[IndeksiTrik[[indT]]],ss++,
pod=SeznamAlf[[indT,indBF]];
(* which triangle *)
kt=IndeksiTrik[[indT,ss]];
(* Index of a point in this triangle (localy). *)
If[TT[[kt,1]]==TV[[indT]],indTLok=1];
If[TT[[kt,2]]==TV[[indT]],indTLok=2];
If[TT[[kt,3]]==TV[[indT]],indTLok=3];
(* Bezier ordinates: *)
MF=BaznaFunZaTocko[PST[[kt]],pod,indTLok];
SeznamBezOrdinat[[kt]]=MF;
];
SeznamKT=Table[Table[Table[0,{i,1,5},{j,1,5}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],4];
For[ii=0,ii≤4,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKT[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],SeznamBezOrdinat[[k,kk,ii+1,jj+1]]];
];
];
];
];
SeznamVsehBezOrdinat=Append[SeznamVsehBezOrdinat,SeznamBezOrdinat];
SeznamVsehKT=Append[SeznamVsehKT,SeznamKT];
];
];
(* --------------- basis functions associated to edges ------------------ *)
SeznamBet=Table[Transpose[DolociParBeta[DetEdge[[indR,1]],DetEdge[[indR,2]],DetEdge[[indR,3]],PSLines[[indR,1]],PSLines[[indR,2]]]],{indR,1,Length[DetEdge]}];
Tsmeri=Table[(PSLines[[k,2]]-PSLines[[k,1]])/Norm[(PSLines[[k,2]]-PSLines[[k,1]])],{k,1,Length[PSLines]}];
(* ------- *)
SeznamVsehBezOrdinatRob={};
SeznamVsehKTRob={};
For[indR=1,indR<=Length[DetEdge],indR++,
For[indBF=1,indBF≤3,indBF++,
SeznamBezOrdinatRob=Table[Table[Table[0,{i,1,5},{j,1,5}],{k,1,6}],{kk,1,Length[TT]}];
For[ss=1,ss≤Length[IndeksiTrikRob[[indR]]],ss++,
pod=SeznamBet[[indR,indBF]];
(* which triangle *)
kt=IndeksiTrikRob[[indR,ss]];
(* Index of an edge in this triangle (localy). *)
trik=TT[[kt]];
toc1=DetEdge[[indR]][[1]];
toc2=DetEdge[[indR]][[2]];
For[s=1,s≤3,s++,
If[toc1==trik[[s]],ind1=s];
If[toc2==trik[[s]],ind2=s];
];
ind=Sort[{ind1,ind2}];
If[ind=={1,2},indRL=1];
If[ind=={2,3},indRL=2];
If[ind=={1,3},indRL=3];
(* Beziero ordinates: *)
MF=BaznaFunZaRob[PST[[kt]],pod,indRL,Tsmeri[[indR]]];
SeznamBezOrdinatRob[[kt]]=MF;
];
SeznamKT=Table[Table[Table[0,{i,1,5},{j,1,5}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],4];
For[ii=0,ii≤4,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKT[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],SeznamBezOrdinatRob[[k,kk,ii+1,jj+1]]];
];
];
];
];
SeznamVsehBezOrdinatRob=Append[SeznamVsehBezOrdinatRob,SeznamBezOrdinatRob];
SeznamVsehKTRob=Append[SeznamVsehKTRob,SeznamKT];
]
];
Return[{SeznamVsehBezOrdinat,SeznamVsehKT,SeznamVsehBezOrdinatRob,SeznamVsehKTRob}];
]


ClearAll[RegTocke]
(* Function returns regular triangulation on a unit square. *)
RegTocke[NN_]:=Module[{TV,TE,TT,i,j},
TV=Flatten[Table[{i/NN,j/NN},{i,0,NN},{j,0,NN}],1];
TE={};
For[i=1,i≤NN,i++,
For[j=0,j≤NN,j++,
TE=Append[TE,{{(i-1)/NN,j/NN},{(i)/NN,j/NN}}];
TE=Append[TE,{{j/NN,(i-1)/NN},{j/NN,(i)/NN}}];
]
];
For[j=0,j≤NN,j++,
For[i=1,i≤NN-j,i++,
TE=Append[TE,{{(i-1+j)/NN,(i-1)/NN},{(i+j)/NN,(i)/NN}}];
];
];
For[j=1,j≤NN,j++,
For[i=j+1,i≤NN,i++,
TE=Append[TE,{{(i-1-j)/NN,(i-1)/NN},{(i-j)/NN,(i)/NN}}];
];
];
TT={};
For[i=1,i≤NN,i++,
For[j=1,j≤NN,j++,
TT=Append[TT,{{(i-1)/NN,(j-1)/NN},{(i)/NN,(j-1)/NN},{(i)/NN,(j)/NN}}];
TT=Append[TT,{{(i-1)/NN,(j-1)/NN},{(i-1)/NN,(j)/NN},{(i)/NN,(j)/NN}}];
]
];
Return[{TV,TE,TT}]
]

ClearAll[Red1];
(* Function computes the approximation order of two consecutive measurements. *)
Red1[napaka1_,h_,napaka2_,k_]:=(Log[napaka1]-Log[napaka2])/(Log[h]-Log[k])

ClearAll[RedAproksimacije];
(* Function computes an approximation order. *)
RedAproksimacije[TabH_,TabN_]:=Module[{tabela,i},
If[Length[TabH]≠Length[TabN],Return[{}]];
tabela={{"h","napaka","eksponent padanja"}};
tabela=Append[tabela,{TabH[[1]],TabN[[1]],"/"}];
For[i=2,i≤Length[TabH],i++,
tabela=Append[tabela,{TabH[[i]],TabN[[i]],Red1[TabN[[i-1]],TabH[[i-1]],TabN[[i]],TabH[[i]]]}];
];
Return[tabela];
]

Examples of B splines

Print["Points of a triangulation:"]
TV={{-5,0},{2,3},{8,0},{4,-3},{0,-4},{0,0},{4,0}}//N
Print["Triangles:"]
TT={{TV[[1]],TV[[6]],TV[[2]]},{TV[[6]],TV[[7]],TV[[2]]},{TV[[7]],TV[[3]],TV[[2]]},{TV[[7]],TV[[4]],TV[[3]]},{TV[[7]],TV[[5]],TV[[4]]},{TV[[6]],TV[[5]],TV[[7]]},{TV[[1]],TV[[5]],TV[[6]]}}
Print["Edges:"]
TE={{TV[[1]],TV[[6]]},{TV[[2]],TV[[6]]},{TV[[1]],TV[[2]]},{TV[[6]],TV[[7]]},{TV[[2]],TV[[7]]},{TV[[2]],TV[[3]]},{TV[[3]],TV[[7]]},{TV[[3]],TV[[4]]},{TV[[4]],TV[[7]]},{TV[[4]],TV[[5]]},{TV[[5]],TV[[7]]},{TV[[5]],TV[[6]]},{TV[[1]],TV[[5]]}}

Print["Computation of all data that are needed. These data are shown on a graph. "]
{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}=IzrKolicin[TV,TT,TE];

Print["Optimal PS triangles (imported): "]
PSTrik={{{-3.382327,-1.294138},{-5.000000,0.000000},{-2.532291,1.057590}},{{-0.375224,1.982047},{2.000000,3.000000},{4.195565,1.902217}},{{6.058997,0.970502},{8.000000,-0.000000},{6.653020,-1.010235}},{{4.203223,-1.851611},{5.300806,-2.400403},{2.699194,-3.599597}},{{2.067180,-2.509721},{0.464887,-4.371910},{-1.625806,-2.699355}},{{-1.476559,0.390594},{1.000000,0.535184},{0.178121,-1.689032}},{{5.693728,0.279952},{2.592063,-1.746597},{2.579517,0.575025}}}

Print["Graphs:"]
G1=Table[Line[{TT[[k,1]],TT[[k,2]],TT[[k,3]],TT[[k,1]]}],{k,1,Length[TT]}];
G2=Table[Table[NarisiKontPoligon[DomainPoints[PST[[k,l,1]],PST[[k,l,2]],PST[[k,l,3]],1]],{l,1,6}],{k,1,Length[TT]}];
G3=Table[Table[NarisiKontPoligon[DomainPoints[PST[[k,l,1]],PST[[k,l,2]],PST[[k,l,3]],4]],{l,1,6}],{k,1,Length[TT]}];
G4=Table[Polygon[PSTrik[[k]]],{k,1,Length[TV]}];
G5=Table[Line[PSLines[[k]]],{k,1,Length[PSLines]}];

Graphics[{GrayLevel[0.9],G3,GrayLevel[0.6],G2,GrayLevel[0],G1}]

G13D=Table[Line[{Append[TT[[k,1]],0],Append[TT[[k,2]],0],Append[TT[[k,3]],0],Append[TT[[k,1]],0]}],{k,1,Length[TT]}];
PST3D=Table[Table[Append[#,0]&/@PST[[k,ll]],{ll,1,6}],{k,1,Length[PST]}];
G23D=Table[Table[NarisiKontPoligon[DomainPoints[PST3D[[k,l,1]],PST3D[[k,l,2]],PST3D[[k,l,3]],1]],{l,1,6}],{k,1,Length[TT]}];
G33D=Table[Table[NarisiKontPoligon[DomainPoints[PST3D[[k,l,1]],PST3D[[k,l,2]],PST3D[[k,l,3]],4]],{l,1,6}],{k,1,Length[TT]}];
grafTriangulacije3D=Graphics3D[{GrayLevel[0.9],G33D,GrayLevel[0.6],G23D,GrayLevel[0],G13D}]
grafTriangulacije3Da=Graphics3D[{GrayLevel[0.6],G23D,GrayLevel[0],G13D}];

Graphics[{GrayLevel[0.9],G3,GrayLevel[0.6],G2,GrayLevel[0],G1,GrayLevel[0.9],G4,Thickness[0.003],RGBColor[0,0,1],G5,RGBColor[1,0,0],PointSize[0.005],Point[TV],Point[Flatten[TPS,1]],PointSize[0.008],RGBColor[0,1,0]},Axes->True]

Print["Computation of all B splines:"]
{SeznamVsehBezOrdinat,SeznamVsehKT,SeznamVsehBezOrdinatRob,SeznamVsehKTRob}=VseBazneFunkcijeZaQuarticSpline[TV,TT,TE,PSTrik];

Print["With a parameter 'izb' we choose a basis function (associated to a vertex) and we draw its graph. "]
kolTrik=10; (* resolution parameter *)
izb=31;
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamVsehKT[[izb,k,kk]],kolTrik]],{kk,1,6}],{k,1,Length[TT]}]];
graf1a=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom},AspectRatio->1];
izb=32;
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamVsehKT[[izb,k,kk]],kolTrik]],{kk,1,6}],{k,1,Length[TT]}]];
graf1b=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom},AspectRatio->1];
izb=33;
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamVsehKT[[izb,k,kk]],kolTrik]],{kk,1,6}],{k,1,Length[TT]}]];
graf1c=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom},AspectRatio->1];
izb=34;
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamVsehKT[[izb,k,kk]],kolTrik]],{kk,1,6}],{k,1,Length[TT]}]];
graf1d=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom},AspectRatio->1];
izb=35;
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamVsehKT[[izb,k,kk]],kolTrik]],{kk,1,6}],{k,1,Length[TT]}]];
graf1e=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom},AspectRatio->1];
izb=36;
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamVsehKT[[izb,k,kk]],kolTrik]],{kk,1,6}],{k,1,Length[TT]}]];
graf1f=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom},AspectRatio->1];

Show[{grafTriangulacije3Da,graf1a},BoxRatios->{1,1,1/4},Boxed->False,ViewPoint->{0,0,2}]
Show[{grafTriangulacije3Da,graf1b},BoxRatios->{1,1,1/4},Boxed->False,ViewPoint->{0,0,2}]
Show[{grafTriangulacije3Da,graf1c},BoxRatios->{1,1,1/4},Boxed->False,ViewPoint->{0,0,2}]
Show[{grafTriangulacije3Da,graf1d},BoxRatios->{1,1,1/4},Boxed->False,ViewPoint->{0,0,2}]
Show[{grafTriangulacije3Da,graf1e},BoxRatios->{1,1,1/4},Boxed->False,ViewPoint->{0,0,2}]
Show[{grafTriangulacije3Da,graf1f},BoxRatios->{1,1,1/4},Boxed->False,ViewPoint->{0,0,2}]

Print["With a parameter 'izb' we choose a basis function (associated to an edge) and we draw its graph. "]
kolTrik=10;  (* resolution parameter *)
izb=10;
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamVsehKTRob[[izb,k,kk]],kolTrik]],{kk,1,6}],{k,1,Length[TT]}]];
graf2a=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom},AspectRatio->1];
izb=11;
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamVsehKTRob[[izb,k,kk]],kolTrik]],{kk,1,6}],{k,1,Length[TT]}]];
graf2b=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom},AspectRatio->1];
izb=12;
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamVsehKTRob[[izb,k,kk]],kolTrik]],{kk,1,6}],{k,1,Length[TT]}]];
graf2c=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom},AspectRatio->1];

Show[{grafTriangulacije3Da,graf2a},BoxRatios->{1,1,1/4},Boxed->False,ViewPoint->{0,0,2}]
Show[{grafTriangulacije3Da,graf2b},BoxRatios->{1,1,1/4},Boxed->False,ViewPoint->{0,0,2}]
Show[{grafTriangulacije3Da,graf2c},BoxRatios->{1,1,1/4},Boxed->False,ViewPoint->{0,0,2}]

Quasi-interpolants

First example

Print["Function that is approximated: "]
ClearAll[x,y]
f=Sin[x]Cos[4/3y]+2

Print["Points of a triangulation:"]
TV={{-5,0},{2,3},{8,0},{4,-3},{0,-4},{0,0},{4,0}}//N
Print["Triangles:"]
TT={{TV[[1]],TV[[6]],TV[[2]]},{TV[[6]],TV[[7]],TV[[2]]},{TV[[7]],TV[[3]],TV[[2]]},{TV[[7]],TV[[4]],TV[[3]]},{TV[[7]],TV[[5]],TV[[4]]},{TV[[6]],TV[[5]],TV[[7]]},{TV[[1]],TV[[5]],TV[[6]]}}
Print["Edges:"]
TE={{TV[[1]],TV[[6]]},{TV[[2]],TV[[6]]},{TV[[1]],TV[[2]]},{TV[[6]],TV[[7]]},{TV[[2]],TV[[7]]},{TV[[2]],TV[[3]]},{TV[[3]],TV[[7]]},{TV[[3]],TV[[4]]},{TV[[4]],TV[[7]]},{TV[[4]],TV[[5]]},{TV[[5]],TV[[7]]},{TV[[5]],TV[[6]]},{TV[[1]],TV[[5]]}}
{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}=IzrKolicin[TV,TT,TE];

Print["Graph of a triangulation: "]
G1=Table[Line[{TT[[k,1]],TT[[k,2]],TT[[k,3]],TT[[k,1]]}],{k,1,Length[TT]}];
G2=Table[Table[NarisiKontPoligon[DomainPoints[PST[[k,l,1]],PST[[k,l,2]],PST[[k,l,3]],1]],{l,1,6}],{k,1,Length[TT]}];
G3=Table[Table[NarisiKontPoligon[DomainPoints[PST[[k,l,1]],PST[[k,l,2]],PST[[k,l,3]],4]],{l,1,6}],{k,1,Length[TT]}];
Graphics[{GrayLevel[0.9],G3,GrayLevel[0.6],G2,GrayLevel[0],G1}]

G13D=Table[Line[{Append[TT[[k,1]],0],Append[TT[[k,2]],0],Append[TT[[k,3]],0],Append[TT[[k,1]],0]}],{k,1,Length[TT]}];
PST3D=Table[Table[Append[#,0]&/@PST[[k,ll]],{ll,1,6}],{k,1,Length[PST]}];
G23D=Table[Table[NarisiKontPoligon[DomainPoints[PST3D[[k,l,1]],PST3D[[k,l,2]],PST3D[[k,l,3]],1]],{l,1,6}],{k,1,Length[TT]}];
G33D=Table[Table[NarisiKontPoligon[DomainPoints[PST3D[[k,l,1]],PST3D[[k,l,2]],PST3D[[k,l,3]],4]],{l,1,6}],{k,1,Length[TT]}];
grafTriangulacije3D=Graphics3D[{GrayLevel[0.9],G33D,GrayLevel[0.6],G23D,GrayLevel[0],G13D}]
grafTriangulacije3Da=Graphics3D[{GrayLevel[0.6],G23D,GrayLevel[0],G13D}];

Print["Graph of a function f(x,y) over a triangulation: "]
kol=10;   (* resolution parameter *)
SeznamKTF=Table[Table[Table[0,{i,1,kol+1},{j,1,kol+1}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],kol];
For[ii=0,ii≤kol,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKTF[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],f/.{x->DD[[ii+1,jj+1,1]],y->DD[[ii+1,jj+1,2]]}];
];
];
];
];
pom1=Flatten[Table[Table[NarisiTrikMrezo[SeznamKTF[[k,kk]]],{kk,1,6}],{k,1,Length[TT]}]];
grafF=Graphics3D[{Opacity[0.5],PointSize[Large],RGBColor[0.1,0.4,0.4],EdgeForm[None],Glow[RGBColor[0.6,0.6,0.6]],Lighting->"Neutral",pom1},BoxRatios->{1,1,1}];
grafF=Graphics3D[{PointSize[Large],RGBColor[0.7,0.1,0.1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom1},BoxRatios->{1,1,1}];
Show[{grafTriangulacije3Da,grafF},BoxRatios->{1,1,1}]

Print["Optimal PS triangles (imported): "]
PSTrik={{{-3.382327,-1.294138},{-5.000000,0.000000},{-2.532291,1.057590}},{{-0.375224,1.982047},{2.000000,3.000000},{4.195565,1.902217}},{{6.058997,0.970502},{8.000000,-0.000000},{6.653020,-1.010235}},{{4.203223,-1.851611},{5.300806,-2.400403},{2.699194,-3.599597}},{{2.067180,-2.509721},{0.464887,-4.371910},{-1.625806,-2.699355}},{{-1.476559,0.390594},{1.000000,0.535184},{0.178121,-1.689032}},{{5.693728,0.279952},{2.592063,-1.746597},{2.579517,0.575025}}}

Print["Computation of all basis functions:"]
{SeznamVsehBezOrdinat,SeznamVsehKT,SeznamVsehBezOrdinatRob,SeznamVsehKTRob}=N[VseBazneFunkcijeZaQuarticSpline[TV,TT,TE,PSTrik]];

Print["B-coefficients of a Hermite type quasi-interpolant: "]
{KoefCV,KoefCE}=KoefKvaziinterpolanta1[f,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines]

Print["Graph of a Hermite type quasi-interpolant:"]
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]],{kk,1,6}],{k,1,Length[TT]}]];
grafS=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom},BoxRatios->{1,1,1}];
grafS=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Lighting->"Neutral",pom},BoxRatios->{1,1,1}];
Show[{grafTriangulacije3Da,grafS,grafF},BoxRatios->{1,1,1},Axes->True]
Print["Error:"]
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]]

Print["B-coefficients of a Lagrange type quasi-interpolant that reproduce polynomials of degree 2: "]
{KoefCV,KoefCE}=KoefKvaziinterpolanta2[f,2,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines]
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]],{kk,1,6}],{k,1,Length[TT]}]];
grafS=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Lighting->"Neutral",pom},BoxRatios->{1,1,1}];
Show[{grafTriangulacije3Da,grafS,grafF},BoxRatios->{1,1,1},Axes->True]
Print["Error:"]
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]]

Print["B-coefficients of a Lagrange type quasi-interpolant that reproduce polynomials of degree 4: "]
{KoefCV,KoefCE}=KoefKvaziinterpolanta2[f,4,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines]
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]],{kk,1,6}],{k,1,Length[TT]}]];
grafS=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Lighting->"Neutral",pom},BoxRatios->{1,1,1}];
Show[{grafTriangulacije3Da,grafS,grafF},BoxRatios->{1,1,1},Axes->True]
Print["Error:"]
nap=Max[Flatten[Table[Table[Abs[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]]],{kk,1,6}],{k,1,Length[TT]}]]]

Second example

Print["Function that is approximated: "]
f=N[Sin[x]Cos[y]Exp[(x+y)/10]+1,20]

Print["Points of a triangulation:"]
TV={{-5,0},{2,3},{8,0},{4,-3},{0,-4},{0,0},{4,0}}//N
Print["Triangles:"]
TT={{TV[[1]],TV[[6]],TV[[2]]},{TV[[6]],TV[[7]],TV[[2]]},{TV[[7]],TV[[3]],TV[[2]]},{TV[[7]],TV[[4]],TV[[3]]},{TV[[7]],TV[[5]],TV[[4]]},{TV[[6]],TV[[5]],TV[[7]]},{TV[[1]],TV[[5]],TV[[6]]}}
Print["Edges:"]
TE={{TV[[1]],TV[[6]]},{TV[[2]],TV[[6]]},{TV[[1]],TV[[2]]},{TV[[6]],TV[[7]]},{TV[[2]],TV[[7]]},{TV[[2]],TV[[3]]},{TV[[3]],TV[[7]]},{TV[[3]],TV[[4]]},{TV[[4]],TV[[7]]},{TV[[4]],TV[[5]]},{TV[[5]],TV[[7]]},{TV[[5]],TV[[6]]},{TV[[1]],TV[[5]]}}
{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}=IzrKolicin[TV,TT,TE];

Print["Graph of a triangulation: "]
G1=Table[Line[{TT[[k,1]],TT[[k,2]],TT[[k,3]],TT[[k,1]]}],{k,1,Length[TT]}];
G2=Table[Table[NarisiKontPoligon[DomainPoints[PST[[k,l,1]],PST[[k,l,2]],PST[[k,l,3]],1]],{l,1,6}],{k,1,Length[TT]}];
G3=Table[Table[NarisiKontPoligon[DomainPoints[PST[[k,l,1]],PST[[k,l,2]],PST[[k,l,3]],4]],{l,1,6}],{k,1,Length[TT]}];
Graphics[{GrayLevel[0.9],G3,GrayLevel[0.6],G2,GrayLevel[0],G1}]

G13D=Table[Line[{Append[TT[[k,1]],0],Append[TT[[k,2]],0],Append[TT[[k,3]],0],Append[TT[[k,1]],0]}],{k,1,Length[TT]}];
PST3D=Table[Table[Append[#,0]&/@PST[[k,ll]],{ll,1,6}],{k,1,Length[PST]}];
G23D=Table[Table[NarisiKontPoligon[DomainPoints[PST3D[[k,l,1]],PST3D[[k,l,2]],PST3D[[k,l,3]],1]],{l,1,6}],{k,1,Length[TT]}];
G33D=Table[Table[NarisiKontPoligon[DomainPoints[PST3D[[k,l,1]],PST3D[[k,l,2]],PST3D[[k,l,3]],4]],{l,1,6}],{k,1,Length[TT]}];
grafTriangulacije3D=Graphics3D[{GrayLevel[0.9],G33D,GrayLevel[0.6],G23D,GrayLevel[0],G13D}]
grafTriangulacije3Da=Graphics3D[{GrayLevel[0.6],G23D,GrayLevel[0],G13D}];

Print["Graph of a function f(x,y) over a triangulation: "]
kol=10;    (* resolution parameter *)
SeznamKTF=Table[Table[Table[0,{i,1,kol+1},{j,1,kol+1}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],kol];
For[ii=0,ii≤kol,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKTF[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],f/.{x->DD[[ii+1,jj+1,1]],y->DD[[ii+1,jj+1,2]]}];
];
];
];
];
pom1=Flatten[Table[Table[NarisiTrikMrezo[SeznamKTF[[k,kk]]],{kk,1,6}],{k,1,Length[TT]}]];
grafF=Graphics3D[{Opacity[0.5],PointSize[Large],RGBColor[0.1,0.4,0.4],EdgeForm[None],Glow[RGBColor[0.6,0.6,0.6]],Lighting->"Neutral",pom1},BoxRatios->{1,1,1}];
grafF=Graphics3D[{PointSize[Large],RGBColor[0.7,0.1,0.1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom1},BoxRatios->{1,1,1}];
Show[{grafTriangulacije3Da,grafF},BoxRatios->{1,1,1}]

Print["Optimal PS triangles (imported): "]
PSTrik={{{-3.382327,-1.294138},{-5.000000,0.000000},{-2.532291,1.057590}},{{-0.375224,1.982047},{2.000000,3.000000},{4.195565,1.902217}},{{6.058997,0.970502},{8.000000,-0.000000},{6.653020,-1.010235}},{{4.203223,-1.851611},{5.300806,-2.400403},{2.699194,-3.599597}},{{2.067180,-2.509721},{0.464887,-4.371910},{-1.625806,-2.699355}},{{-1.476559,0.390594},{1.000000,0.535184},{0.178121,-1.689032}},{{5.693728,0.279952},{2.592063,-1.746597},{2.579517,0.575025}}}

Print["Computation of all basis functions:"]
{SeznamVsehBezOrdinat,SeznamVsehKT,SeznamVsehBezOrdinatRob,SeznamVsehKTRob}=N[VseBazneFunkcijeZaQuarticSpline[TV,TT,TE,PSTrik]];

Print["B-coefficients of a Hermite type quasi-interpolant: "]
{KoefCV,KoefCE}=KoefKvaziinterpolanta1[f,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines]
Print["Graph of a Hermite type quasi-interpolant:"]
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]],{kk,1,6}],{k,1,Length[TT]}]];
grafS=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom},BoxRatios->{1,1,1}];
grafS=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Lighting->"Neutral",pom},BoxRatios->{1,1,1}];
Show[{grafTriangulacije3Da,grafS,grafF},BoxRatios->{1,1,1},Axes->True]
Print["Error:"]
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]]

Print["B-coefficients of a Lagrange type quasi-interpolant that reproduce polynomials of degree 2: "]
{KoefCV,KoefCE}=KoefKvaziinterpolanta2[f,2,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines]
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]],{kk,1,6}],{k,1,Length[TT]}]];
grafS=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Lighting->"Neutral",pom},BoxRatios->{1,1,1}];
Show[{grafTriangulacije3Da,grafS,grafF},BoxRatios->{1,1,1},Axes->True]
Print["Error:"]
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]]

Print["B-coefficients of a Lagrange type quasi-interpolant that reproduce polynomials of degree 4: "]
{KoefCV,KoefCE}=KoefKvaziinterpolanta2[f,4,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines]
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]],{kk,1,6}],{k,1,Length[TT]}]];
grafS=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Lighting->"Neutral",pom},BoxRatios->{1,1,1}];
Show[{grafTriangulacije3Da,grafS,grafF},BoxRatios->{1,1,1},Axes->True]
Print["Error:"]
nap=Max[Flatten[Table[Table[Abs[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]]],{kk,1,6}],{k,1,Length[TT]}]]]

Quasi-interpolant on a regular triangulation

Print["Regular type-I triangulation (Subscript[Δ, 8]):"]
{TV,TE,TT}=N[RegTocke[8]];
{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}=IzrKolicinTip1[TV,TT,TE,1/8];
G1=Table[Line[{TT[[k,1]],TT[[k,2]],TT[[k,3]],TT[[k,1]]}],{k,1,Length[TT]}];
G2=Table[Table[NarisiKontPoligon[DomainPoints[PST[[k,l,1]],PST[[k,l,2]],PST[[k,l,3]],1]],{l,1,6}],{k,1,Length[TT]}];
G3=Table[Table[NarisiKontPoligon[DomainPoints[PST[[k,l,1]],PST[[k,l,2]],PST[[k,l,3]],4]],{l,1,6}],{k,1,Length[TT]}];
Graphics[{GrayLevel[0],G1}]

Print["Regular type-I triangulation (Subscript[Δ, 4]):"]
{TV,TE,TT}=N[RegTocke[4]];
{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}=IzrKolicinTip1[TV,TT,TE,1/4];
G1=Table[Line[{TT[[k,1]],TT[[k,2]],TT[[k,3]],TT[[k,1]]}],{k,1,Length[TT]}];
Graphics[{GrayLevel[0],G1}]

Print["Graph of a triangulation, PS triangles, PS lines, ..."]
G1=Table[Line[{TT[[k,1]],TT[[k,2]],TT[[k,3]],TT[[k,1]]}],{k,1,Length[TT]}];
G2=Table[Table[NarisiKontPoligon[DomainPoints[PST[[k,l,1]],PST[[k,l,2]],PST[[k,l,3]],1]],{l,1,6}],{k,1,Length[TT]}];
G3=Table[Table[NarisiKontPoligon[DomainPoints[PST[[k,l,1]],PST[[k,l,2]],PST[[k,l,3]],4]],{l,1,6}],{k,1,Length[TT]}];
Graphics[{GrayLevel[0.9],G3,GrayLevel[0.6],G2,GrayLevel[0],G1}]

G13D=Table[Line[{Append[TT[[k,1]],0],Append[TT[[k,2]],0],Append[TT[[k,3]],0],Append[TT[[k,1]],0]}],{k,1,Length[TT]}];
PST3D=Table[Table[Append[#,0]&/@PST[[k,ll]],{ll,1,6}],{k,1,Length[PST]}];
G23D=Table[Table[NarisiKontPoligon[DomainPoints[PST3D[[k,l,1]],PST3D[[k,l,2]],PST3D[[k,l,3]],1]],{l,1,6}],{k,1,Length[TT]}];
G33D=Table[Table[NarisiKontPoligon[DomainPoints[PST3D[[k,l,1]],PST3D[[k,l,2]],PST3D[[k,l,3]],4]],{l,1,6}],{k,1,Length[TT]}];
G4=Table[Polygon[PSTrik[[k]]],{k,1,Length[TV]}];
G5=Table[Line[PSLines[[k]]],{k,1,Length[PSLines]}];
grafTriangulacije3D=Graphics3D[{GrayLevel[0.9],G33D,GrayLevel[0.6],G23D,GrayLevel[0],G13D}];
grafTriangulacije3Da=Graphics3D[{GrayLevel[0.6],G23D,GrayLevel[0],G13D}];
Graphics[{GrayLevel[0.9],G4,GrayLevel[0.9],G3,GrayLevel[0.6],G2,GrayLevel[0],G1,Thickness[0.003],RGBColor[0,0,1],G5,RGBColor[1,0,0],PointSize[0.005],Point[TV]},Axes->True]
Graphics[{EdgeForm[{GrayLevel[0.8],Thickness[0.003]}],GrayLevel[0.85],G4,GrayLevel[0.8],Thickness[0.0001],G3,GrayLevel[0.4],Thickness[0.001],G2,GrayLevel[0],Thickness[0.0015],G1,Thickness[0.0035],RGBColor[0.3,0.3,0.3],G5,RGBColor[0.3,0.3,0.3],PointSize[0.008],Point[TV],Point[Flatten[TPS,1]],PointSize[0.008],RGBColor[1,1,1]}]

Print["Function that is approximated: "]
f=N[Sin[5x]Cos[5y]Log[1+x+y],20]
(*f=N[3/4Exp[-((9x-2)^2+(9y-2)^2)]+3/4Exp[-(9x+1)^2/49-(9y+1)/10]+1/2Exp[-((9x-7)^2+(9y-3)^2/4)]-1/5Exp[-(9x-4)^2-(9y-7)^2],20]*)

Print["Graph of a function f(x,y) over a triangulation: "]
kol=10;    (* resolution parameter *)
SeznamKTF=Table[Table[Table[0,{i,1,kol+1},{j,1,kol+1}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],kol];
For[ii=0,ii≤kol,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKTF[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],f/.{x->DD[[ii+1,jj+1,1]],y->DD[[ii+1,jj+1,2]]}];
];
];
];
];
pom1=Flatten[Table[Table[NarisiTrikMrezo[SeznamKTF[[k,kk]]],{kk,1,6}],{k,1,Length[TT]}]];
grafF=Graphics3D[{PointSize[Large],RGBColor[0.7,0.1,0.1],EdgeForm[None],Glow[RGBColor[0.3,0.3,0.3]],Lighting->"Neutral",pom1},BoxRatios->{1,1,1}];
Show[{grafF},BoxRatios->{1,1,1},Axes->True]

Print["Computation of all basis functions:"]
{SeznamVsehBezOrdinat,SeznamVsehKT,SeznamVsehBezOrdinatRob,SeznamVsehKTRob}=N[VseBazneFunkcijeZaQuarticSpline[TV,TT,TE,PSTrik]];

Print["B-coefficients of a Hermite type quasi-interpolant: "]
{KoefCV,KoefCE}=KoefKvaziinterpolanta1[f,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines]

Print["Graph of a Hermite type quasi-interpolant:"]
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]],{kk,1,6}],{k,1,Length[TT]}]];
grafS=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Lighting->"Neutral",pom},BoxRatios->{1,1,1}];
Show[{grafS,grafF,Graphics3D[{G13D}]},BoxRatios->{1,1,1},Axes->True]
Print["Error:"]
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]]

Print["B-coefficients of a Lagrange type quasi-interpolant that reproduce polynomials of degree 2: "]
{KoefCV,KoefCE}=KoefKvaziinterpolanta2[f,2,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines]

Print["Graph of a Lagrange type quasi-interpolant:"]
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]],{kk,1,6}],{k,1,Length[TT]}]];
grafS=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Lighting->"Neutral",pom},BoxRatios->{1,1,1}];
Show[{grafS,grafF,Graphics3D[{G13D}]},BoxRatios->{1,1,1},Axes->True]
Print["Error:"]
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]]

Print["B-coefficients of a Lagrange type quasi-interpolant that reproduce polynomials of degree 4: "]
{KoefCV,KoefCE}=KoefKvaziinterpolanta2[f,4,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines];

Print["Graph of a Lagrange type quasi-interpolant:"]
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
pom=Flatten[Table[Table[NarisiTrikMrezo[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]],{kk,1,6}],{k,1,Length[TT]}]];
grafS=Graphics3D[{PointSize[Large],RGBColor[0.3,0.4,1],EdgeForm[None],Lighting->"Neutral",pom},BoxRatios->{1,1,1}];
Show[{grafS,grafF,Graphics3D[{G13D}]},BoxRatios->{1,1,1},Axes->True]
Print["Error:"]
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]]

Approximation order

First example

Print["Function that is approximated: "]
f=N[3/4Exp[-((9x-2)^2+(9y-2)^2)]+3/4Exp[-(9x+1)^2/49-(9y+1)/10]+1/2Exp[-((9x-7)^2+(9y-3)^2/4)]-1/5Exp[-(9x-4)^2-(9y-7)^2],20]

Print["Approximation order for a Hermite type quasi-interpolant. Note that PS triangles are not optimal! "]

TabNap={};
For[NN=2,NN≤8,NN++,
(* Points of a triangulation: *)
{TV,TE,TT}=N[RegTocke[NN]];
{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}=IzrKolicin[TV,TT,TE];
{SeznamVsehBezOrdinat,SeznamVsehKT}=N[BazneFunkcijeZaTocke[TV,TT,TE]];
{SeznamVsehBezOrdinatRob,SeznamVsehKTRob}=N[BazneFunkcijeZaRob[TV,TT,TE]];
{KoefCV,KoefCE}=KoefKvaziinterpolanta1[f,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines];
(* coefficients of the approximant *)
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
(* values of a function f *)
kol=25;
SeznamKTF=Table[Table[Table[0,{i,1,kol+1},{j,1,kol+1}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],kol];
For[ii=0,ii≤kol,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKTF[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],f/.{x->DD[[ii+1,jj+1,1]],y->DD[[ii+1,jj+1,2]]}];
];
];
];
];
(* error *)
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]];
TabNap=Append[TabNap,nap];
Print[nap];
]

Print["Table of the approximation error: "]
TabH=Table[1/NN,{NN,2,8}];
RedAproksimacije[TabH,TabNap]//MatrixForm

Print["Approximation order for a Lagrange type quasi-interpolant that reproduce polynomials of degree 4:"]

TabNap={};
For[NN=2,NN≤8,NN++,
(* Points of a triangulation: *)
{TV,TE,TT}=N[RegTocke[NN]];
{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}=IzrKolicin[TV,TT,TE];
{SeznamVsehBezOrdinat,SeznamVsehKT}=N[BazneFunkcijeZaTocke[TV,TT,TE]];
{SeznamVsehBezOrdinatRob,SeznamVsehKTRob}=N[BazneFunkcijeZaRob[TV,TT,TE]];
{KoefCV,KoefCE}=KoefKvaziinterpolanta2[f,4,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines];
(* coefficients of the approximant *)
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
(* values of a function f *)
kol=25;
SeznamKTF=Table[Table[Table[0,{i,1,kol+1},{j,1,kol+1}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],kol];
For[ii=0,ii≤kol,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKTF[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],f/.{x->DD[[ii+1,jj+1,1]],y->DD[[ii+1,jj+1,2]]}];
];
];
];
];
(* error *)
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]];
TabNap=Append[TabNap,nap];
Print[nap];
]

Print["Table of the approximation error: "]
TabH=Table[1/NN,{NN,2,8}];
RedAproksimacije[TabH,TabNap]//MatrixForm

Second example

Print["Function that is approximated: "]
f=N[Sin[5x]Cos[5y],20]

Print["Approximation order for a Hermite type quasi-interpolant. Note that PS triangles are not optimal! "]

TabNap={};
For[NN=2,NN≤8,NN++,
(* Points of a triangulation: *)
{TV,TE,TT}=N[RegTocke[NN]];
{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}=IzrKolicin[TV,TT,TE];
{SeznamVsehBezOrdinat,SeznamVsehKT}=N[BazneFunkcijeZaTocke[TV,TT,TE]];
{SeznamVsehBezOrdinatRob,SeznamVsehKTRob}=N[BazneFunkcijeZaRob[TV,TT,TE]];
{KoefCV,KoefCE}=KoefKvaziinterpolanta1[f,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines];
(* coefficients of the approximant *)
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
(* values of a function f *)
kol=25;
SeznamKTF=Table[Table[Table[0,{i,1,kol+1},{j,1,kol+1}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],kol];
For[ii=0,ii≤kol,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKTF[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],f/.{x->DD[[ii+1,jj+1,1]],y->DD[[ii+1,jj+1,2]]}];
];
];
];
];
(* error *)
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]];
TabNap=Append[TabNap,nap];
Print[nap];
]

Print["Table of the approximation error: "]
TabH=Table[1/NN,{NN,2,8}];
RedAproksimacije[TabH,TabNap]//MatrixForm

Print["Approximation order for a Lagrange type quasi-interpolant that reproduce polynomials of degree 4:"]

TabNap={};
For[NN=2,NN≤8,NN++,
(* Points of a triangulation: *)
{TV,TE,TT}=N[RegTocke[NN]];
{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}=IzrKolicin[TV,TT,TE];
{SeznamVsehBezOrdinat,SeznamVsehKT}=N[BazneFunkcijeZaTocke[TV,TT,TE]];
{SeznamVsehBezOrdinatRob,SeznamVsehKTRob}=N[BazneFunkcijeZaRob[TV,TT,TE]];
{KoefCV,KoefCE}=KoefKvaziinterpolanta2[f,4,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines];
(* coefficients of the approximant *)
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
(* values of a function f *)
kol=25;
SeznamKTF=Table[Table[Table[0,{i,1,kol+1},{j,1,kol+1}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],kol];
For[ii=0,ii≤kol,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKTF[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],f/.{x->DD[[ii+1,jj+1,1]],y->DD[[ii+1,jj+1,2]]}];
];
];
];
];
(* error *)
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]];
TabNap=Append[TabNap,nap];
Print[nap];
]

Print["Table of the approximation error: "]
TabH=Table[1/NN,{NN,2,8}];
RedAproksimacije[TabH,TabNap]//MatrixForm

Third example

Print["Function that is approximated: "]
f=N[Sin[5x]Cos[5y]Log[1+x+y],20]

Print["Approximation order for a Hermite type quasi-interpolant. Note that PS triangles are not optimal! "]

TabNap={};
For[NN=2,NN≤8,NN++,
(* Points of a triangulation: *)
{TV,TE,TT}=N[RegTocke[NN]];
{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}=IzrKolicinTip1[TV,TT,TE,1/NN];
{SeznamVsehBezOrdinat,SeznamVsehKT,SeznamVsehBezOrdinatRob,SeznamVsehKTRob}=N[VseBazneFunkcijeZaQuarticSpline[TV,TT,TE,PSTrik]];
{KoefCV,KoefCE}=KoefKvaziinterpolanta1[f,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines];
(* coefficients of the approximant *)
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
(* values of a function f *)
kol=25;
SeznamKTF=Table[Table[Table[0,{i,1,kol+1},{j,1,kol+1}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],kol];
For[ii=0,ii≤kol,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKTF[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],f/.{x->DD[[ii+1,jj+1,1]],y->DD[[ii+1,jj+1,2]]}];
];
];
];
];
(* error *)
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]];
TabNap=Append[TabNap,nap];
Print[nap];
]

Print["Table of the approximation error: "]
TabH=Table[1/NN,{NN,2,8}];
RedAproksimacije[TabH,TabNap]//MatrixForm

Print["Approximation order for a Lagrange type quasi-interpolant that reproduce polynomials of degree 4:"]

TabNap={};
For[NN=2,NN≤8,NN++,
(* Points of a triangulation: *)
{TV,TE,TT}=N[RegTocke[NN]];
{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}=IzrKolicinTip1[TV,TT,TE,1/NN];
{SeznamVsehBezOrdinat,SeznamVsehKT,SeznamVsehBezOrdinatRob,SeznamVsehKTRob}=N[VseBazneFunkcijeZaQuarticSpline[TV,TT,TE,PSTrik]];
{KoefCV,KoefCE}=KoefKvaziinterpolanta2[f,4,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines];
(* coefficients of the approximant *)
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
(* values of a function f *)
kol=25;
SeznamKTF=Table[Table[Table[0,{i,1,kol+1},{j,1,kol+1}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],kol];
For[ii=0,ii≤kol,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKTF[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],f/.{x->DD[[ii+1,jj+1,1]],y->DD[[ii+1,jj+1,2]]}];
];
];
];
];
(* error *)
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]];
TabNap=Append[TabNap,nap];
Print[nap];
]

Print["Table of the approximation error: "]
TabH=Table[1/NN,{NN,2,8}];
RedAproksimacije[TabH,TabNap]//MatrixForm

Print["Approximation order for a Lagrange type quasi-interpolant that reproduce polynomials of degree 2:"]

TabNap={};
For[NN=2,NN≤8,NN++,
(* Points of a triangulation: *)
{TV,TE,TT}=N[RegTocke[NN]];
{TZ,IndeksiTrik,IndeksiTrikRob,IndRobTrik,IndeksiRobV,PST,TPS,PSTrik,DetEdge,PSLines,Tsmeri,TsmeriV}=IzrKolicinTip1[TV,TT,TE,1/NN];
{SeznamVsehBezOrdinat,SeznamVsehKT,SeznamVsehBezOrdinatRob,SeznamVsehKTRob}=N[VseBazneFunkcijeZaQuarticSpline[TV,TT,TE,PSTrik]];
(*{SeznamVsehBezOrdinat,SeznamVsehKT}=N[BazneFunkcijeZaTocke[TV,TT,TE]];*)
(*{SeznamVsehBezOrdinatRob,SeznamVsehKTRob}=N[BazneFunkcijeZaRob[TV,TT,TE]];*)
{KoefCV,KoefCE}=KoefKvaziinterpolanta2[f,2,x,y,SeznamVsehBezOrdinat,TT,PST,PSTrik,TV,DetEdge,IndeksiTrikRob,Tsmeri,TsmeriV,PSLines];
(* coefficients of the approximant *)
{sez,SeznamKT}=DolociBezOrdinateIzKoef[KoefCV,KoefCE,N[SeznamVsehBezOrdinat],N[SeznamVsehBezOrdinatRob],PST];
(* values of a function f *)
kol=25;
SeznamKTF=Table[Table[Table[0,{i,1,kol+1},{j,1,kol+1}],{k,1,6}],{kk,1,Length[TT]}];
For[k=1,k≤Length[TT],k++,
For[kk=1,kk≤6,kk++,
DD=DomainPoints[PST[[k,kk,1]],PST[[k,kk,2]],PST[[k,kk,3]],kol];
For[ii=0,ii≤kol,ii++,
For[jj=0,jj≤ii,jj++,
SeznamKTF[[k,kk,ii+1,jj+1]]=Append[DD[[ii+1,jj+1]],f/.{x->DD[[ii+1,jj+1,1]],y->DD[[ii+1,jj+1,2]]}];
];
];
];
];
(* error *)
nap=Max[Abs[Flatten[Table[Table[ComputeValuesForGraph[SeznamKT[[k,kk]],kol]-SeznamKTF[[k,kk]],{kk,1,6}],{k,1,Length[TT]}]]]];
TabNap=Append[TabNap,nap];
Print[nap];
]

Print["Table of the approximation error: "]
TabH=Table[1/NN,{NN,2,8}];
RedAproksimacije[TabH,TabNap]//MatrixForm