> restart: 
# 
# Code for the computation of equilibrium conformations of 
# fluid films with bending stiffness
# Scheme: Differences
# Structures: tube6, tube4, tube3, cylinder, sphere
# 
# Version 1.0,    6-2-2010
# 
# Author: Arnold Krawietz, Berlin, Germany,  krawietz@t-online.de  

> with(linalg):
> with(plots):
# Procedures:  symmetry condition on the boundaries
> boundv:=proc()
> local zeta,gxu,gyu,gzu,gxvtilde,gyvtilde,gzvtilde,guu,egv,
> gucrossex,gucrossey,gucrossez,wx,wy,wz,ww,psi,gx_oben_v,gy_oben_v,gz_oben_v,
> wgv,wm,zetav,gugv,gum,guv,gx_oben_u,gy_oben_u,gz_oben_u,
> dzeta,dgxu,dgyu,dgzu,dgxvtilde,dgyvtilde,dgzvtilde,degv,
> dgucrossex,dgucrossey,dgucrossez,dwx,dwy,dwz,
> dwgv,dwm,dzetav,dgugv,dgum,dguv,dguu,dpsi,hiii,
> dgx_oben_u,dgy_oben_u,dgz_oben_u,dgx_oben_v,dgy_oben_v,dgz_oben_v,
> hatzeta,hatgxu,hatgyu,hatgzu,hatgxvtilde,hatgyvtilde,hatgzvtilde,hategv,
> hatgucrossex,hatgucrossey,hatgucrossez,hatwx,hatwy,hatwz,hatzetav,hatgxv,hatgyv,hatgzv,
> hatgx_oben_u,hatgy_oben_u,hatgz_oben_u,hatgx_oben_v,hatgy_oben_v,hatgz_oben_v,
> hatwgv,hatwm,hatgugv,hatgum,hatguv,hatguu,hatpsi,hii,dhii,
> dhatwgv,dhatzetav,dhatgugv,dhatguv,dhatguu,dhatpsi,dhatwx,dhatwy,dhatwz,
> dhatgx_oben_u,dhatgy_oben_u,dhatgz_oben_u,dhatgx_oben_v,dhatgy_oben_v,dhatgz_oben_v;
> global wr,psir,obenu,obenv,dwr,dpsir,dobenu,dobenv,
> hatwr,hatpsir,hatobenu,hatobenv,
> dhatwr,dhatpsir,dhatobenu,dhatobenv;
> zeta:=0.5*a*(theta[iii]+theta[iii-1]);
> gxu:=xr[3*iii-2]-xr[3*iii-5];
> gyu:=xr[3*iii-1]-xr[3*iii-4];
> gzu:=xr[3*iii]-xr[3*iii-3];
> gxvtilde:=((x0m-x00)*(1-u/n)+(xnm-xn0)*u/n)/m+zeta*mmxv;
> gyvtilde:=((y0m-y00)*(1-u/n)+(ynm-yn0)*u/n)/m+zeta*mmyv;
> gzvtilde:=zeta*mmzv;
> egv:=e[1]*gxvtilde+e[2]*gyvtilde+e[3]*gzvtilde;
> gucrossex:=gyu*e[3]-gzu*e[2];
> gucrossey:=gzu*e[1]-gxu*e[3];
> gucrossez:=gxu*e[2]-gyu*e[1];
> wx:=egv*gucrossex;
> wy:=egv*gucrossey;
> wz:=egv*gucrossez;
> ww:=wx^2+wy^2+wz^2;
> psi:=1./sqrt(ww);
>   gx_oben_v:=e[1]/egv;
>   gy_oben_v:=e[2]/egv;
>   gz_oben_v:=e[3]/egv;
> wgv:=gucrossex*gxvtilde+gucrossey*gyvtilde+gucrossez*gzvtilde;
> wm:=gucrossex*mmx+gucrossey*mmy+gucrossez*mmz;
> zetav:= -wgv/wm;
> gugv:=gxu*gxvtilde+gyu*gyvtilde+gzu*gzvtilde;
> gum:=gxu*mmx+gyu*mmy+gzu*mmz;
> guv:=gugv+zetav*gum;
> guu:=gxu^2+gyu^2+gzu^2;
>   gx_oben_u:=(gxu-guv/egv*e[1])/guu;
>   gy_oben_u:=(gyu-guv/egv*e[2])/guu;
>   gz_oben_u:=(gzu-guv/egv*e[3])/guu;
> 
> wr[3*ii-2]:=wx; wr[3*ii-1]:=wy; wr[3*ii]:=wz;
> obenu[3*ii-2]:=gx_oben_u; obenu[3*ii-1]:=gy_oben_u; obenu[3*ii]:=gz_oben_u;
> obenv[3*ii-2]:=gx_oben_v; obenv[3*ii-1]:=gy_oben_v; obenv[3*ii]:=gz_oben_v;
> psir[ii]:=psi;
> 
> dzeta:=0.5*a*(dtheta[iii]+dtheta[iii-1]);
> dgxu:=dxr[3*iii-2]-dxr[3*iii-5];
> dgyu:=dxr[3*iii-1]-dxr[3*iii-4];
> dgzu:=dxr[3*iii]-dxr[3*iii-3];
> dgxvtilde:=dzeta*mmxv;
> dgyvtilde:=dzeta*mmyv;
> dgzvtilde:=dzeta*mmzv;
> degv:=e[1]*dgxvtilde+e[2]*dgyvtilde+e[3]*dgzvtilde;
> dgucrossex:=dgyu*e[3]-dgzu*e[2];
> dgucrossey:=dgzu*e[1]-dgxu*e[3];
> dgucrossez:=dgxu*e[2]-dgyu*e[1];
> dwx:=degv*gucrossex+egv*dgucrossex;
> dwy:=degv*gucrossey+egv*dgucrossey;
> dwz:=degv*gucrossez+egv*dgucrossez;
> 
> dpsi:=-psi^3*(wx*dwx+wy*dwy+wz*dwz);
>     dgx_oben_v:=-degv/egv^2*e[1];
>     dgy_oben_v:=-degv/egv^2*e[2];
>     dgz_oben_v:=-degv/egv^2*e[3];
> dwgv:=dgucrossex*gxvtilde+dgucrossey*gyvtilde+dgucrossez*gzvtilde+gucrossex*dgxvtilde+gucrossey*dgyvtilde+gucrossez*dgzvtilde;
> dwm:=dgucrossex*mmx+dgucrossey*mmy+dgucrossez*mmz;    
> dzetav:= -(dwgv+zetav*dwm)/wm;
> dgugv:=dgxu*gxvtilde+dgyu*gyvtilde+dgzu*gzvtilde+gxu*dgxvtilde+gyu*dgyvtilde+gzu*dgzvtilde;
> dgum:=dgxu*mmx+dgyu*mmy+dgzu*mmz;
> dguv:=dgugv+dzetav*gum+zetav*dgum;
> dguu:=2*(gxu*dgxu+gyu*dgyu+gzu*dgzu);
>     hiii:=(guv*degv-dguv*egv)/egv^2;
>     dgx_oben_u:=(-dguu*gx_oben_u +dgxu+hiii*e[1])/guu;
>     dgy_oben_u:=(-dguu*gy_oben_u +dgyu+hiii*e[2])/guu;
>     dgz_oben_u:=(-dguu*gz_oben_u +dgzu+hiii*e[3])/guu;
> 
> dwr[3*ii-2]:=dwx; dwr[3*ii-1]:=dwy; dwr[3*ii]:=dwz;   
> dobenu[3*ii-2]:=dgx_oben_u; dobenu[3*ii-1]:=dgy_oben_u; dobenu[3*ii]:=dgz_oben_u;
> dobenv[3*ii-2]:=dgx_oben_v; dobenv[3*ii-1]:=dgy_oben_v; dobenv[3*ii]:=dgz_oben_v;
> dpsir[ii]:=dpsi;
> 
> hatzeta:=0.5*a*(hattheta[iii]+hattheta[iii-1]);
> hatgxu:=hatxr[3*iii-2]-hatxr[3*iii-5];
> hatgyu:=hatxr[3*iii-1]-hatxr[3*iii-4];
> hatgzu:=hatxr[3*iii]-hatxr[3*iii-3];
> hatgxvtilde:=hatzeta*mmxv;
> hatgyvtilde:=hatzeta*mmyv;
> hatgzvtilde:=hatzeta*mmzv;
> hategv:=e[1]*hatgxvtilde+e[2]*hatgyvtilde+e[3]*hatgzvtilde;
> hatgucrossex:=hatgyu*e[3]-hatgzu*e[2];
> hatgucrossey:=hatgzu*e[1]-hatgxu*e[3];
> hatgucrossez:=hatgxu*e[2]-hatgyu*e[1];
> hatwx:=hategv*gucrossex+egv*hatgucrossex;
> hatwy:=hategv*gucrossey+egv*hatgucrossey;
> hatwz:=hategv*gucrossez+egv*hatgucrossez;
> hatpsi:=-psi^3*(wx*hatwx+wy*hatwy+wz*hatwz);
>     hatgx_oben_v:=-hategv/egv^2*e[1];
>     hatgy_oben_v:=-hategv/egv^2*e[2];
>     hatgz_oben_v:=-hategv/egv^2*e[3];
> hatwgv:=hatgucrossex*gxvtilde+hatgucrossey*gyvtilde+hatgucrossez*gzvtilde+gucrossex*hatgxvtilde+gucrossey*hatgyvtilde+gucrossez*hatgzvtilde;
> hatwm:=hatgucrossex*mmx+hatgucrossey*mmy+hatgucrossez*mmz;
> hatzetav:= -(hatwgv+zetav*hatwm)/wm;
> hatgugv:=hatgxu*gxvtilde+hatgyu*gyvtilde+hatgzu*gzvtilde+gxu*hatgxvtilde+gyu*hatgyvtilde+gzu*hatgzvtilde;
> hatgum:=hatgxu*mmx+hatgyu*mmy+hatgzu*mmz;
> hatguv:=hatgugv+hatzetav*gum+zetav*hatgum;
> hatguu:=2*(gxu*hatgxu+gyu*hatgyu+gzu*hatgzu);
>     hii:=(guv*hategv-hatguv*egv)/egv^2;
>     hatgx_oben_u:=(-hatguu*gx_oben_u +hatgxu+hii*e[1])/guu;
>     hatgy_oben_u:=(-hatguu*gy_oben_u +hatgyu+hii*e[2])/guu;
>     hatgz_oben_u:=(-hatguu*gz_oben_u +hatgzu+hii*e[3])/guu;    
> 
> hatwr[3*ii-2]:=hatwx; hatwr[3*ii-1]:=hatwy; hatwr[3*ii]:=hatwz;
> hatobenu[3*ii-2]:=hatgx_oben_u; hatobenu[3*ii-1]:=hatgy_oben_u; hatobenu[3*ii]:=hatgz_oben_u;
> hatobenv[3*ii-2]:=hatgx_oben_v; hatobenv[3*ii-1]:=hatgy_oben_v; hatobenv[3*ii]:=hatgz_oben_v;
> hatpsir[ii]:=hatpsi;
> 
> dhatwx:=hategv*dgucrossex+degv*hatgucrossex;
> dhatwy:=hategv*dgucrossey+degv*hatgucrossey;
> dhatwz:=hategv*dgucrossez+degv*hatgucrossez;
> dhatpsi:=3*dpsi*hatpsi/psi-psi^3*(dwx*hatwx+dwy*hatwy+dwz*hatwz+wx*dhatwx+wy*dhatwy+wz*dhatwz);
>     dhatgx_oben_v:=-2*degv*hatgx_oben_v/egv; 
>     dhatgy_oben_v:=-2*degv*hatgy_oben_v/egv;
>     dhatgz_oben_v:=-2*degv*hatgz_oben_v/egv;
> dhatwgv:=hatgucrossex*dgxvtilde+hatgucrossey*dgyvtilde+hatgucrossez*dgzvtilde+dgucrossex*hatgxvtilde+dgucrossey*hatgyvtilde+dgucrossez*hatgzvtilde;
> dhatzetav:= -(hatzetav*dwm+dzetav*hatwm+dhatwgv)/wm;
> dhatgugv:=hatgxu*dgxvtilde+hatgyu*dgyvtilde+hatgzu*dgzvtilde+dgxu*hatgxvtilde+dgyu*hatgyvtilde+dgzu*hatgzvtilde;
> dhatguv:=dhatgugv+dhatzetav*gum+hatzetav*dgum+dzetav*hatgum;
> dhatguu:=2*(dgxu*hatgxu+dgyu*hatgyu+dgzu*hatgzu);
>     dhii:=(dguv*hategv-dhatguv*egv-hatguv*degv)/egv^2-2*degv*hii/egv;
>     dhatgx_oben_u:=(-dhatguu*gx_oben_u -hatguu*dgx_oben_u+dhii*e[1] -dguu*hatgx_oben_u)/guu;
>     dhatgy_oben_u:=(-dhatguu*gy_oben_u -hatguu*dgy_oben_u +dhii*e[2] -dguu*hatgy_oben_u)/guu;
>     dhatgz_oben_u:=(-dhatguu*gz_oben_u -hatguu*dgz_oben_u+dhii*e[3] -dguu*hatgz_oben_u)/guu;
> 
> dhatwr[3*ii-2]:=dhatwx; dhatwr[3*ii-1]:=dhatwy; dhatwr[3*ii]:=dhatwz;
> dhatobenu[3*ii-2]:=dhatgx_oben_u; dhatobenu[3*ii-1]:=dhatgy_oben_u; dhatobenu[3*ii]:=dhatgz_oben_u;
> dhatobenv[3*ii-2]:=dhatgx_oben_v; dhatobenv[3*ii-1]:=dhatgy_oben_v; dhatobenv[3*ii]:=dhatgz_oben_v;
> dhatpsir[ii]:=dhatpsi;
> end:

> 
> boundu:=proc()
> local zeta,gxv,gyv,gzv,gxutilde,gyutilde,gzutilde,gvv,egu,
> gvcrossex,gvcrossey,gvcrossez,wx,wy,wz,ww,psi,gx_oben_u,gy_oben_u,gz_oben_u,
> wgu,wm,zetau,gvgu,gvm,guv,gx_oben_v,gy_oben_v,gz_oben_v,
> dzeta,dgxv,dgyv,dgzv,dgxutilde,dgyutilde,dgzutilde,degu,
> dgvcrossex,dgvcrossey,dgvcrossez,dwx,dwy,dwz,
> dwgu,dwm,dzetau,dgvgu,dgvm,dguv,dgvv,dpsi,hiii,
> dgx_oben_u,dgy_oben_u,dgz_oben_u,dgx_oben_v,dgy_oben_v,dgz_oben_v,
> hatzeta,hatgxv,hatgyv,hatgzv,hatgxutilde,hatgyutilde,hatgzutilde,hategu,
> hatgvcrossex,hatgvcrossey,hatgvcrossez,hatwx,hatwy,hatwz,
> hatwgu,hatwm,hatzetau,hatgvgu,hatgvm,hatguv,hatgvv,hatpsi,hii,
> hatgx_oben_u,hatgy_oben_u,hatgz_oben_u,hatgx_oben_v,hatgy_oben_v,hatgz_oben_v,
> dhatwx,dhatwy,dhatwz,dhatpsi,dhatwgu,dhatzetau,dhatgvgu,dhatguv,dhatgvv,dhii,
> dhatgx_oben_u,dhatgy_oben_u,dhatgz_oben_u,dhatgx_oben_v,dhatgy_oben_v,dhatgz_oben_v;
> global wr,psir,obenu,obenv,dwr,dpsir,dobenu,dobenv,
> hatwr,hatpsir,hatobenu,hatobenv,dhatwr,dhatpsir,dhatobenu,dhatobenv;
> zeta:=0.5*a*(theta[iii]+theta[iii-n1i]);
> gxv:=xr[3*iii-2]-xr[3*iii-3*n1i-2];
> gyv:=xr[3*iii-1]-xr[3*iii-3*n1i-1];
> gzv:=xr[3*iii]-xr[3*iii-3*n1i];
> gxutilde:=((xn0-x00)*(1-v/m)+(xnm-x0m)*v/m)/n+zeta*mmxu;
> gyutilde:=((yn0-y00)*(1-v/m)+(ynm-y0m)*v/m)/n+zeta*mmyu;
> gzutilde:=zeta*mmzu;
> egu:=e[1]*gxutilde+e[2]*gyutilde+e[3]*gzutilde;
> gvcrossex:=gyv*e[3]-gzv*e[2];
> gvcrossey:=gzv*e[1]-gxv*e[3];
> gvcrossez:=gxv*e[2]-gyv*e[1];
> wx:=-egu*gvcrossex;
> wy:=-egu*gvcrossey;
> wz:=-egu*gvcrossez;
> ww:=wx^2+wy^2+wz^2;
> psi:=1./sqrt(ww);
>   gx_oben_u:=e[1]/egu;
>   gy_oben_u:=e[2]/egu;
>   gz_oben_u:=e[3]/egu;
> wgu:=gvcrossex*gxutilde+gvcrossey*gyutilde+gvcrossez*gzutilde;
> wm:=gvcrossex*mmx+gvcrossey*mmy+gvcrossez*mmz;
> zetau:= -wgu/wm;
> gvgu:=gxv*gxutilde+gyv*gyutilde+gzv*gzutilde;
> gvm:=gxv*mmx+gyv*mmy+gzv*mmz;
> guv:=gvgu+zetau*gvm;
> gvv:=gxv^2+gyv^2+gzv^2;
> gx_oben_v:=(gxv-guv/egu*e[1])/gvv;
> gy_oben_v:=(gyv-guv/egu*e[2])/gvv;
> gz_oben_v:=(gzv-guv/egu*e[3])/gvv;
> 
> wr[3*ii-2]:=wx; wr[3*ii-1]:=wy; wr[3*ii]:=wz;
> obenu[3*ii-2]:=gx_oben_u; obenu[3*ii-1]:=gy_oben_u; obenu[3*ii]:=gz_oben_u;
> obenv[3*ii-2]:=gx_oben_v; obenv[3*ii-1]:=gy_oben_v; obenv[3*ii]:=gz_oben_v;
> psir[ii]:=psi;
> 
> dzeta:=0.5*a*(dtheta[iii]+dtheta[iii-n1i]);
> dgxv:=dxr[3*iii-2]-dxr[3*iii-3*n1i-2];
> dgyv:=dxr[3*iii-1]-dxr[3*iii-3*n1i-1];
> dgzv:=dxr[3*iii]-dxr[3*iii-3*n1i];
> dgxutilde:=dzeta*mmxu;
> dgyutilde:=dzeta*mmyu;
> dgzutilde:=dzeta*mmzu;
> degu:=e[1]*dgxutilde+e[2]*dgyutilde+e[3]*dgzutilde;
> dgvcrossex:=dgyv*e[3]-dgzv*e[2];
> dgvcrossey:=dgzv*e[1]-dgxv*e[3];
> dgvcrossez:=dgxv*e[2]-dgyv*e[1];
> dwx:=-degu*gvcrossex-egu*dgvcrossex;
> dwy:=-degu*gvcrossey-egu*dgvcrossey;
> dwz:=-degu*gvcrossez-egu*dgvcrossez;
> dpsi:=-psi^3*(wx*dwx+wy*dwy+wz*dwz);
>     dgx_oben_u:=-degu/egu^2*e[1];
>     dgy_oben_u:=-degu/egu^2*e[2];
>     dgz_oben_u:=-degu/egu^2*e[3];
> dwgu:=dgvcrossex*gxutilde+dgvcrossey*gyutilde+dgvcrossez*gzutilde+gvcrossex*dgxutilde+gvcrossey*dgyutilde+gvcrossez*dgzutilde;
> dwm:=dgvcrossex*mmx+dgvcrossey*mmy+dgvcrossez*mmz;
> dzetau:= -(dwgu+zetau*dwm)/wm;
> dgvgu:=dgxv*gxutilde+dgyv*gyutilde+dgzv*gzutilde+gxv*dgxutilde+gyv*dgyutilde+gzv*dgzutilde;
> dgvm:=dgxv*mmx+dgyv*mmy+dgzv*mmz;
> dguv:=dgvgu+dzetau*gvm+zetau*dgvm;
> dgvv:=2*(gxv*dgxv+gyv*dgyv+gzv*dgzv);
>     hiii:=(guv*degu-dguv*egu)/egu^2;
>     dgx_oben_v:=(-dgvv*gx_oben_v +dgxv+hiii*e[1])/gvv;
>     dgy_oben_v:=(-dgvv*gy_oben_v +dgyv+hiii*e[2])/gvv;
>     dgz_oben_v:=(-dgvv*gz_oben_v +dgzv+hiii*e[3])/gvv; 
> 
> dwr[3*ii-2]:=dwx; dwr[3*ii-1]:=dwy; dwr[3*ii]:=dwz;
> dobenu[3*ii-2]:=dgx_oben_u; dobenu[3*ii-1]:=dgy_oben_u; dobenu[3*ii]:=dgz_oben_u;
> dobenv[3*ii-2]:=dgx_oben_v; dobenv[3*ii-1]:=dgy_oben_v; dobenv[3*ii]:=dgz_oben_v;
> dpsir[ii]:=dpsi;
> 
> hatzeta:=0.5*a*(hattheta[iii]+hattheta[iii-n1i]);
> hatgxv:=hatxr[3*iii-2]-hatxr[3*iii-3*n1i-2];
> hatgyv:=hatxr[3*iii-1]-hatxr[3*iii-3*n1i-1];
> hatgzv:=hatxr[3*iii]-hatxr[3*iii-3*n1i];
> hatgxutilde:=hatzeta*mmxu;
> hatgyutilde:=hatzeta*mmyu;
> hatgzutilde:=hatzeta*mmzu;
> hategu:=e[1]*hatgxutilde+e[2]*hatgyutilde+e[3]*hatgzutilde;
> hatgvcrossex:=hatgyv*e[3]-hatgzv*e[2];
> hatgvcrossey:=hatgzv*e[1]-hatgxv*e[3];
> hatgvcrossez:=hatgxv*e[2]-hatgyv*e[1];
> hatwx:=-hategu*gvcrossex-egu*hatgvcrossex;
> hatwy:=-hategu*gvcrossey-egu*hatgvcrossey;
> hatwz:=-hategu*gvcrossez-egu*hatgvcrossez;
> hatpsi:=-psi^3*(wx*hatwx+wy*hatwy+wz*hatwz);
>     hatgx_oben_u:=-hategu/egu^2*e[1];
>     hatgy_oben_u:=-hategu/egu^2*e[2];
>     hatgz_oben_u:=-hategu/egu^2*e[3];
> hatwgu:=hatgvcrossex*gxutilde+hatgvcrossey*gyutilde+hatgvcrossez*gzutilde+gvcrossex*hatgxutilde+gvcrossey*hatgyutilde+gvcrossez*hatgzutilde;
> hatwm:=hatgvcrossex*mmx+hatgvcrossey*mmy+hatgvcrossez*mmz;
> hatzetau:= -(hatwgu+zetau*hatwm)/wm;
> hatgvgu:=hatgxv*gxutilde+hatgyv*gyutilde+hatgzv*gzutilde+gxv*hatgxutilde+gyv*hatgyutilde+gzv*hatgzutilde;
> hatgvm:=hatgxv*mmx+hatgyv*mmy+hatgzv*mmz;
> hatguv:=hatgvgu+hatzetau*gvm+zetau*hatgvm;
> hatgvv:=2*(gxv*hatgxv+gyv*hatgyv+gzv*hatgzv);
>     hii:=(guv*hategu-hatguv*egu)/egu^2;
>     hatgx_oben_v:=(-hatgvv*gx_oben_v +hatgxv+hii*e[1])/gvv;
>     hatgy_oben_v:=(-hatgvv*gy_oben_v +hatgyv+hii*e[2])/gvv;
>     hatgz_oben_v:=(-hatgvv*gz_oben_v +hatgzv+hii*e[3])/gvv;    
> 
> hatwr[3*ii-2]:=hatwx; hatwr[3*ii-1]:=hatwy; hatwr[3*ii]:=hatwz;
> hatobenu[3*ii-2]:=hatgx_oben_u; hatobenu[3*ii-1]:=hatgy_oben_u; hatobenu[3*ii]:=hatgz_oben_u;
> hatobenv[3*ii-2]:=hatgx_oben_v; hatobenv[3*ii-1]:=hatgy_oben_v; hatobenv[3*ii]:=hatgz_oben_v;
> hatpsir[ii]:=hatpsi;
> 
> 
> dhatwx:=-hategu*dgvcrossex-degu*hatgvcrossex;
> dhatwy:=-hategu*dgvcrossey-degu*hatgvcrossey;
> dhatwz:=-hategu*dgvcrossez-degu*hatgvcrossez;
> dhatpsi:=3*dpsi*hatpsi/psi-psi^3*(dwx*hatwx+dwy*hatwy+dwz*hatwz+wx*dhatwx+wy*dhatwy+wz*dhatwz);
>     dhatgx_oben_u:=-2*degu*hatgx_oben_u/egu; 
>     dhatgy_oben_u:=-2*degu*hatgy_oben_u/egu;
>     dhatgz_oben_u:=-2*degu*hatgz_oben_u/egu;
> dhatwgu:=hatgvcrossex*dgxutilde+hatgvcrossey*dgyutilde+hatgvcrossez*dgzutilde+dgvcrossex*hatgxutilde+dgvcrossey*hatgyutilde+dgvcrossez*hatgzutilde;
> dhatzetau:= -(hatzetau*dwm+dzetau*hatwm+dhatwgu)/wm;
> dhatgvgu:=hatgxv*dgxutilde+hatgyv*dgyutilde+hatgzv*dgzutilde+dgxv*hatgxutilde+dgyv*hatgyutilde+dgzv*hatgzutilde;
> dhatguv:=dhatgvgu+dhatzetau*gvm+hatzetau*dgvm+dzetau*hatgvm;
> dhatgvv:=2*(dgxv*hatgxv+dgyv*hatgyv+dgzv*hatgzv);
>     dhii:=(dguv*hategu-dhatguv*egu-hatguv*degu)/egu^2-2*degu*hii/egu;
>     dhatgx_oben_v:=(-dhatgvv*gx_oben_v -hatgvv*dgx_oben_v+dhii*e[1] -dgvv*hatgx_oben_v)/gvv;
>     dhatgy_oben_v:=(-dhatgvv*gy_oben_v -hatgvv*dgy_oben_v+dhii*e[2] -dgvv*hatgy_oben_v)/gvv;
>     dhatgz_oben_v:=(-dhatgvv*gz_oben_v -hatgvv*dgz_oben_v+dhii*e[3] -dgvv*hatgz_oben_v)/gvv;    
> 
> dhatwr[3*ii-2]:=dhatwx; dhatwr[3*ii-1]:=dhatwy; dhatwr[3*ii]:=dhatwz;
> dhatobenu[3*ii-2]:=dhatgx_oben_u; dhatobenu[3*ii-1]:=dhatgy_oben_u; dhatobenu[3*ii]:=dhatgz_oben_u;
> dhatobenv[3*ii-2]:=dhatgx_oben_v; dhatobenv[3*ii-1]:=dhatgy_oben_v; dhatobenv[3*ii]:=dhatgz_oben_v;
> dhatpsir[ii]:=dhatpsi;
> end:
> 
> 
# Structures:
# Please activate one of them as "true"
> sphere:=false:
> cylinder:=false:
> tube3:=false:
> tube4:=false:
> tube6:=true:

> 
> if sphere then
> typ:="sphere";
>  vop0:=a^3/3; 
>  vwp0:=100.*a^3
>  fi;

> if tube6 then
> typ:="tube6";
>  vop0:=2*a^3/3; 
>  vwp0:=2*a^3/3
> fi; 
> if cylinder then
> typ:="cylinder";
>  vop0:=a^3/2; 
>  vwp0:=100.*a^3
> fi; 
> if tube3 then
> typ:="tube3";
>  vop0:=a^3/54.; 
>  vwp0:=100.*a^3
> fi; 
> if tube4 then
> typ:="tube4";
>  vop0:=a^3/9/sqrt(3.); 
>  vwp0:=a^3/9/sqrt(3.)
> fi; 
> 
> if sphere or cylinder or tube6 then
>  barx00:=0.; barxn0:=1.; barx0m:=0.; barxnm:=1.; # important:floating point numbers
> bary00:=0.; baryn0:=0.; bary0m:=1.; barynm:=1.; # important:floating point numbers
> barzz0:=0.
> fi;
> if tube3 then
>  barx00:=-1/2/sqrt(6.); barxn0:=1/2/sqrt(6.); barx0m:=-1/2/sqrt(6.); barxnm:=1/2/sqrt(6.); 
> bary00:=-1/2.; baryn0:=-1/2.; bary0m:=-1/6.; barynm:=0.;
> barzz0:=-0.5/sqrt(6.)
> fi; 
> if tube4 then
>  barx00:=-1/3/sqrt(2.); barxn0:=1/3/sqrt(2.); barx0m:=-1/3/sqrt(2.); barxnm:=1/3/sqrt(2.);
> bary00:=-2/3.; baryn0:=-5/6.; bary0m:=-1/6.; barynm:=-1/3.;
> barzz0:=-1./sqrt(6.)
> fi; 
> 
# a denotes e with tube3 and tube 4 and e/4 with tube6; it is arbitrary but fixed with cylinders and spheres
> x00:=a*barx00; xn0:=a*barxn0; x0m:=a*barx0m; xnm:=a*barxnm;
> y00:=a*bary00; yn0:=a*baryn0; y0m:=a*bary0m; ynm:=a*barynm;
> zz0:=a*barzz0;
> u:='u': v:='v':
> if sphere or tube6 then
> 
>  mmz:=1./sqrt(1.+u^2/n^2+v^2/m^2);
>  mmx:=-u/n*mmz; 
> fi;  
> if sphere then
>  mmy:=-v/m*mmz
> fi;
> if tube6 then
>  mmy:=v/m*mmz
> fi; 
> 
> if cylinder then
>  mmz:=1./sqrt(1.+u^2/n^2); 
>  mmx:=-u/n*mmz; 
>  mmy:=0.
> fi;
> if tube3 then
> mmz:=1./sqrt(1.+(1.-2*u/n)^2+2/3.*(1-u/n)^2*(v/m)^2); 
> 
> mmx:=(1-2*u/n)*mmz; 
> mmy:=sqrt(2/3.)*(1-u/n)*v/m*mmz
> fi;
> 
> 
> if tube4 then
> mmz:=1./sqrt(1.+(1.-2*u/n)^2/3.+(u/n+3*v/m-2)^2/6.); 
> 
> mmx:=(1-2*u/n)/sqrt(3.)*mmz; 
> mmy:=(u/n+3*v/m-2)/sqrt(6.)*mmz
> fi; 
> 
> 
> mmxu:=simplify(diff(mmx,u)):
> mmyu:=simplify(diff(mmy,u)):
> mmzu:=simplify(diff(mmz,u)):
> mmxv:=simplify(diff(mmx,v)):
> mmyv:=simplify(diff(mmy,v)):
> mmzv:=simplify(diff(mmz,v)):

# Characterization of the mixture (change data if necessary):
> hh:=0.01; # thickness of the film
> 
> L_o:=10000.; #oil length of the mixture
> L_w:=10000.; #water length of the mixture
> kf:=1: # constraint on amphiphile mass
> ko:=0; # constraint on oil volume
> kw:=0; # constraint on water volume
> ka:=1; # allow variation of length a 
# Characterization of the film material (change data if necessary):
# phi denotes the energy density (per unit area) w, w_1 denotes lambda_1, w_2 denotes lambda_2, thetam denotes theta
> phi:=unapply(w0+w1*(trc/2-H0)^2+w2*(trc^2/4-dtc-D0^2)^2/
> (thetam*D0^2+trc^2/4-dtc),(trc,dtc));
> w0:=0.;
> w1:=1.;
> nu:='nu':
> w2:=(1-nu)/(1+nu)*w1;
> nu:=0.9;
> thetam:=1./3.;

> H0:=1.0;
> D0:=1.0;
> phi:=unapply(w0+w1*(trc/2-H0)^2+w2*(trc^2/4-dtc-D0^2)^2/
> (thetam*D0^2+trc^2/4-dtc),(trc,dtc));
> phitr:=D[1](phi);
> phidt:=D[2](phi);
> 
> phitrtr:=D[1](phitr);
> phitrdt:=D[2](phitr);
> phidtdt:=D[2](phidt);
# 
> 
# Grid (change data if necessary):
# n1i is r+1, n2i is s+1,  n is r,  m is s
> n1i:=16: n2i:=16: 
> if not type(n1i+n2i,integer) or n1i<4 or n2i<4 then print("Warning: wrong mesh") fi: 
> n1:=convert(n1i,float);
> n2:=convert(n2i,float);
> n:=n1-1.; m:=n2-1.; 
> ni:=n1i-1: mi:=n2i-1:
> nn:=n1i*n2i;
> 
# Generation of a starting vector
# Approach A:  Insert a starting vector (usually from some previous computation)
# The list has  (r+1)*(s+1)+5 entries
> ############## change if necessary
> thetali:= 
> [.6552413117, .6691173039, .7096582930, .7738786245, .857\
> 9398484, .9580037314, 1.070627666, .6151085476, .6286992420, .6683395981, .731\
> 3427697, .8140346834, .9126082400, 1.023899606, .5250033356, .5377159308, .575\
> 0608673, .6343877562, .7124900848, .8062090074, .9123418748, .4214517693, .433\
> 1052160, .4671537731, .5219278659, .5946476125, .6821836024, .7821469290, .315\
> 0332353, .3261599201, .3587782941, .4104343288, .4787794034, .5615990175, .656\
> 2662629, .2046024907, .2162410652, .2494603011, .3017014654, .3695216000, .450\
> 1596832, .5417412931, .8909337649e-1, .1015386527, .1375926246, .1933707284, .\
> 2645690207, .3473456869, .4393453079, 3.325636558, .2582630823, 0., 0., .69925\
> 61045];
> 
> theta:=vector(thetali):
> 
# Approach  B: Transition from one grid to another one (skip if not relevant):
> 
> n1iold:=26; n2iold:=26; # change if necessary
> 
>  
> thetapre:=NULL:
> for j to n2i do
> for i to n1i do
>   u:=(i-1)*(n1iold-1)/(n1-1); v:=(j-1)*(n2iold-1)/(n2-1);
>   flu:=floor(u); flv:=floor(v);
>   if flu=n1iold-1 then flu:=flu-1 fi;
>   if flv=n2iold-1 then flv:=flv-1 fi;
>   xi:=u-flu; eta:=v-flv;
>   kn00:=flu+1+flv*n1iold;
>   kn10:=kn00+1;
>   kn01:=kn00+n1iold;
>   kn11:=kn01+1;
> if flu>0 then axiu:=2*theta[kn01]-theta[kn01-1]-theta[kn11]  fi;
> if flu<n1iold-2 then axio:=2*theta[kn11]-theta[kn01]-theta[kn11+1] fi;
> if flv>0 then aetau:=2*theta[kn10]-theta[kn10-n1iold]-theta[kn11]  fi;
> if flv<n2iold-2 then aetao:=2*theta[kn11]-theta[kn10]-theta[kn11+n1iold] fi;
> if flu>0 and flu<n1iold-2 then axi:=(1-xi)*axiu+xi*axio fi;
> if flv>0 and flv<n2iold-2 then aeta:=(1-eta)*aetau+eta*aetao fi;
> if flu=0 then axi:=axio fi;
> if flu=n1iold-2 then axi:=axiu fi;
> if flv=0 then aeta:=aetao fi;
> if flv=n2iold-2 then aeta:=aetau fi;
>   tth:=theta[kn00]*(1-xi)*(1-eta)+theta[kn10]*xi*(1-eta)
>           +theta[kn01]*(1-xi)*eta+theta[kn11]*xi*eta+0.5*(axi*xi*(1-xi)+aeta*eta*(1-eta));
>   thetapre:=thetapre,tth
> od od;
> 
> nnold:=n1iold*n2iold;
> thetapre:=thetapre,theta[nnold+1],theta[nnold+2],theta[nnold+3],theta[nnold+4],theta[nnold+5]:
> thetali:=map(evalf,[thetapre]);
> theta:=vector(thetali):
> 
# Start of the computation:
# 
> ##############################################################################################################################
> wcalt:=1:
> acalt:=1:
> vocalt:=1:
> vwcalt:=1:
> 
> 
> mar:=NULL:
> mas:=NULL:
> 
> if ko=0 then theta[nn+3]:=0 fi:
> if kw=0 then theta[nn+4]:=0 fi:
# Label:   Come back to this point to start the next iteration step  
# XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
> 
# Checks at the four corners:
> 
> warning:=false:
> 
> 
> if theta[1]>1 then warning:=true fi:
> if tube6 and theta[1]<-1 then warning:=true fi:
> if tube4 and abs(theta[1])>1/sqrt(3.) then warning:=true fi:
> if tube3 and (theta[1]>0.5/sqrt(3.) or theta[1]<-1./sqrt(3.)) then warning:=true fi:

> 
> 
> if theta[n1i]>sqrt(2.) then warning:=true fi:
> if tube6 and theta[n1i]<-sqrt(2.) then warning:=true fi:
> if tube4 and abs(theta[n1i])>.5 then warning:=true fi:
> if tube3 and theta[n1i]>0.5/sqrt(3.) then warning:=true fi:
> 
> 
> if cylinder and theta[(n2i-1)*n1i+1]>1. then warning:=true fi:
> if theta[(n2i-1)*n1i+1]>sqrt(2.) then warning:=true fi:
> if tube6 and theta[(n2i-1)*n1i+1]<-sqrt(2.) then warning:=true fi:
> if tube4 and abs(theta[(n2i-1)*n1i+1])>.5 then warning:=true fi:
> if tube3 and (theta[(n2i-1)*n1i+1]>1./3. or theta[(n2i-1)*n1i+1]<-2./3.) then warning:=true fi:
> 
> if cylinder and theta[nn]>sqrt(2.) then warning:=true fi:
> if theta[nn]>sqrt(3.) then warning:=true fi:
> if tube6 and theta[nn]<-sqrt(3.) then warning:=true fi:
> if tube4 and abs(theta[nn])>1/sqrt(3.) then warning:=true fi:
> if tube3 and theta[nn]>0.5/sqrt(3.) then warning:=true fi:

> if warning then 
>  print("Error: Illegal displacement")
> fi:
> 
> 
# Diagrams of the distribution of theta values and their second differences over u and v:
# 
> 
> for j to n2i do
> 
> ztl[j]:=[0,theta[(j-1)*n1i+1]]:
> ddztl[j]:=NULL:

> for i to n1i-1 do
> ztl[j]:=ztl[j],[i,theta[(j-1)*n1i+i+1]]
> od;
> for i to n1i-2 do
> ddztl[j]:=ddztl[j],[i,theta[(j-1)*n1i+i]+theta[(j-1)*n1i+i+2]-2*theta[(j-1)*n1i+i+1]]
> od
> od:
> 
> 
> plot([seq([ztl[j]],j=1..n2i)]);
> plot([seq([ddztl[j]],j=1..n2i)]);
> 
> for i to n1i do
> 
> ztl[i]:=[0,theta[i]]:
> ddztl[i]:=NULL:

> for j to n2i-1 do
> ztl[i]:=ztl[i],[j,theta[j*n1i+i]]
> od;

> for j to n2i-2 do
> ddztl[i]:=ddztl[i],[j,theta[(j-1)*n1i+i] +theta[(j+1)*n1i+i]-2*theta[j*n1i+i] ]
> od
> 
> od:
> 
> plot([seq([ztl[i]],i=1..n1i)]);
> 
> 
> plot([seq([ddztl[i]],i=1..n1i)]);
> 
> 
> 
> 
# Initial values:
> 
> ns:=theta[nn+1]:
> f:=theta[nn+2]:
> po:=theta[nn+3]:
> pw:=theta[nn+4]:
> a:=theta[nn+5]:
# 
# Loci of the grid points:
> dtheta:=vector(nn):
> hattheta:=vector(nn):
> xseq:=NULL:
> dxseq:=NULL:
> hatxseq:=NULL:

> for j to n2i do 
> for i to n1i do
> 
> zeta:=a*theta[(j-1)*n1i+i];
> dzeta:=a*dtheta[(j-1)*n1i+i];
> hatzeta:=a*hattheta[(j-1)*n1i+i];
> 
> u:=i-1; v:=j-1;  
>        x:=x00*(1-u/n)*(1-v/m)+xn0*u/n*(1-v/m)
>        +x0m*(1-u/n)*v/m+xnm*u/n*v/m + zeta*mmx;
>       y:=y00*(1-u/n)*(1-v/m)+yn0*u/n*(1-v/m)
>        +y0m*(1-u/n)*v/m+ynm*u/n*v/m + zeta*mmy;
>       z:=zz0+zeta*mmz;
> xseq:=xseq,x,y,z;
>       dx:= dzeta*mmx;
>       dy:= dzeta*mmy;
>       dz:= dzeta*mmz;
> dxseq:=dxseq,dx,dy,dz;
>       hatx:= hatzeta*mmx;
>       haty:= hatzeta*mmy;
>       hatz:= hatzeta*mmz;
> hatxseq:=hatxseq,hatx,haty,hatz;
> od
> od;
> 
> xr:=vector([xseq]):
> dxr:=vector([dxseq]):
> hatxr:=vector([hatxseq]):

# Tangential plane at the nodes:
> 
> psir:=vector((n1i+1)*(n2i+1)):
> dpsir:=vector((n1i+1)*(n2i+1)):
> hatpsir:=vector((n1i+1)*(n2i+1)):
> dhatpsir:=vector((n1i+1)*(n2i+1)):
> 
> wr:=vector(3*(n1i+1)*(n2i+1)):
> dwr:=vector(3*(n1i+1)*(n2i+1)):
> hatwr:=vector(3*(n1i+1)*(n2i+1)):
> dhatwr:=vector(3*(n1i+1)*(n2i+1)):
> 
> obenu:=vector(3*(n1i+1)*(n2i+1)):
> obenv:=vector(3*(n1i+1)*(n2i+1)):
> dobenu:=vector(3*(n1i+1)*(n2i+1)):
> dobenv:=vector(3*(n1i+1)*(n2i+1)):
> hatobenu:=vector(3*(n1i+1)*(n2i+1)):
> hatobenv:=vector(3*(n1i+1)*(n2i+1)):
> dhatobenu:=vector(3*(n1i+1)*(n2i+1)):
> dhatobenv:=vector(3*(n1i+1)*(n2i+1)):
> 
# Inner nodes:
> for j from 2 to n2i do
> for i from 2 to n1i do
> 
> ii:=(j-1)*(n1i+1)+i;
> iii:=(j-1)*n1i+i;
> 
> gxu:=0.5*(xr[3*iii-2] -xr[3*iii-5] +xr[3*(iii-n1i)-2] -xr[3*(iii-n1i)-5]);
> gyu:=0.5*(xr[3*iii-1] -xr[3*iii-4] +xr[3*(iii-n1i)-1] -xr[3*(iii-n1i)-4]);
> gzu:=0.5*(xr[3*iii] -xr[3*iii-3] +xr[3*(iii-n1i)] -xr[3*(iii-n1i)-3]);
> gxv:=0.5*(xr[3*iii-2] +xr[3*iii-5] -xr[3*(iii-n1i)-2] -xr[3*(iii-n1i)-5]);
> gyv:=0.5*(xr[3*iii-1] +xr[3*iii-4] -xr[3*(iii-n1i)-1] -xr[3*(iii-n1i)-4]);
> gzv:=0.5*(xr[3*iii] +xr[3*iii-3] -xr[3*(iii-n1i)] -xr[3*(iii-n1i)-3]);
> 
>     wx:=gyu*gzv-gyv*gzu;
>     wy:=gxv*gzu-gxu*gzv;
>     wz:=gxu*gyv-gxv*gyu;     
>     ww:=wx^2+wy^2+wz^2;
>     psi:=1./sqrt(ww);    
>     guu:=gxu^2+gyu^2+gzu^2;
>     guv:=gxu*gxv+gyu*gyv+gzu*gzv;
>     gvv:=gxv^2+gyv^2+gzv^2;
>     gx_oben_u:=psi^2*(gvv*gxu-guv*gxv);
>     gy_oben_u:=psi^2*(gvv*gyu-guv*gyv);
>     gz_oben_u:=psi^2*(gvv*gzu-guv*gzv);
>     gx_oben_v:=psi^2*(guu*gxv-guv*gxu);
>     gy_oben_v:=psi^2*(guu*gyv-guv*gyu);
>     gz_oben_v:=psi^2*(guu*gzv-guv*gzu);
> 
>     wr[3*ii-2]:=wx; wr[3*ii-1]:=wy; wr[3*ii]:=wz;
>     obenu[3*ii-2]:=gx_oben_u; obenu[3*ii-1]:=gy_oben_u; obenu[3*ii]:=gz_oben_u;
>     obenv[3*ii-2]:=gx_oben_v; obenv[3*ii-1]:=gy_oben_v; obenv[3*ii]:=gz_oben_v;
>     psir[ii]:=psi;
> 
> dgxu:=0.5*(dxr[3*iii-2] -dxr[3*iii-5] +dxr[3*(iii-n1i)-2] -dxr[3*(iii-n1i)-5]);
> dgyu:=0.5*(dxr[3*iii-1] -dxr[3*iii-4] +dxr[3*(iii-n1i)-1] -dxr[3*(iii-n1i)-4]);
> dgzu:=0.5*(dxr[3*iii] -dxr[3*iii-3] +dxr[3*(iii-n1i)] -dxr[3*(iii-n1i)-3]);
> dgxv:=0.5*(dxr[3*iii-2] +dxr[3*iii-5] -dxr[3*(iii-n1i)-2] -dxr[3*(iii-n1i)-5]);
> dgyv:=0.5*(dxr[3*iii-1] +dxr[3*iii-4] -dxr[3*(iii-n1i)-1] -dxr[3*(iii-n1i)-4]);
> dgzv:=0.5*(dxr[3*iii] +dxr[3*iii-3] -dxr[3*(iii-n1i)] -dxr[3*(iii-n1i)-3]);
>     dwx:=dgyu*gzv-gyv*dgzu+gyu*dgzv-dgyv*gzu;
>     dwy:=gxv*dgzu-dgxu*gzv+dgxv*gzu-gxu*dgzv;
>     dwz:=dgxu*gyv-gxv*dgyu+gxu*dgyv-dgxv*gyu;
> dpsi:=-psi^3*(wx*dwx+wy*dwy+wz*dwz);
> 
> dguu:=2*(gxu*dgxu+gyu*dgyu+gzu*dgzu);
> dgvv:=2*(gxv*dgxv+gyv*dgyv+gzv*dgzv);
> dguv:=gxu*dgxv+gyu*dgyv+gzu*dgzv +gxv*dgxu+gyv*dgyu+gzv*dgzu;
> 
>     dgx_oben_u:=2*dpsi/psi*gx_oben_u +psi^2*(dgvv*gxu-dguv*gxv+gvv*dgxu-guv*dgxv);
>     dgy_oben_u:=2*dpsi/psi*gy_oben_u +psi^2*(dgvv*gyu-dguv*gyv+gvv*dgyu-guv*dgyv);
>     dgz_oben_u:=2*dpsi/psi*gz_oben_u +psi^2*(dgvv*gzu-dguv*gzv+gvv*dgzu-guv*dgzv);
>     dgx_oben_v:=2*dpsi/psi*gx_oben_v +psi^2*(dguu*gxv-dguv*gxu+guu*dgxv-guv*dgxu);
>     dgy_oben_v:=2*dpsi/psi*gy_oben_v +psi^2*(dguu*gyv-dguv*gyu+guu*dgyv-guv*dgyu);
>     dgz_oben_v:=2*dpsi/psi*gz_oben_v +psi^2*(dguu*gzv-dguv*gzu+guu*dgzv-guv*dgzu);
> dwr[3*ii-2]:=dwx; dwr[3*ii-1]:=dwy; dwr[3*ii]:=dwz;
> dobenu[3*ii-2]:=dgx_oben_u; dobenu[3*ii-1]:=dgy_oben_u; dobenu[3*ii]:=dgz_oben_u;
> dobenv[3*ii-2]:=dgx_oben_v; dobenv[3*ii-1]:=dgy_oben_v; dobenv[3*ii]:=dgz_oben_v;
> dpsir[ii]:=dpsi;
> 
> hatgxu:=0.5*(hatxr[3*iii-2] -hatxr[3*iii-5] +hatxr[3*(iii-n1i)-2] -hatxr[3*(iii-n1i)-5]);
> hatgyu:=0.5*(hatxr[3*iii-1] -hatxr[3*iii-4] +hatxr[3*(iii-n1i)-1] -hatxr[3*(iii-n1i)-4]);
> hatgzu:=0.5*(hatxr[3*iii] -hatxr[3*iii-3] +hatxr[3*(iii-n1i)] -hatxr[3*(iii-n1i)-3]);
> hatgxv:=0.5*(hatxr[3*iii-2] +hatxr[3*iii-5] -hatxr[3*(iii-n1i)-2] -hatxr[3*(iii-n1i)-5]);
> hatgyv:=0.5*(hatxr[3*iii-1] +hatxr[3*iii-4] -hatxr[3*(iii-n1i)-1] -hatxr[3*(iii-n1i)-4]);
> hatgzv:=0.5*(hatxr[3*iii] +hatxr[3*iii-3] -hatxr[3*(iii-n1i)] -hatxr[3*(iii-n1i)-3]);
>     hatwx:=hatgyu*gzv-gyv*hatgzu+gyu*hatgzv-hatgyv*gzu;
>     hatwy:=gxv*hatgzu-hatgxu*gzv+hatgxv*gzu-gxu*hatgzv;
>     hatwz:=hatgxu*gyv-gxv*hatgyu+gxu*hatgyv-hatgxv*gyu;
> hatpsi:=-psi^3*(wx*hatwx+wy*hatwy+wz*hatwz);
> 
> hatguu:=2*(gxu*hatgxu+gyu*hatgyu+gzu*hatgzu);
> hatgvv:=2*(gxv*hatgxv+gyv*hatgyv+gzv*hatgzv);
> hatguv:=gxu*hatgxv+gyu*hatgyv+gzu*hatgzv +gxv*hatgxu+gyv*hatgyu+gzv*hatgzu;
> 
>     hatgx_oben_u:=2*hatpsi/psi*gx_oben_u +psi^2*(hatgvv*gxu-hatguv*gxv+gvv*hatgxu-guv*hatgxv);
>     hatgy_oben_u:=2*hatpsi/psi*gy_oben_u +psi^2*(hatgvv*gyu-hatguv*gyv+gvv*hatgyu-guv*hatgyv);
>     hatgz_oben_u:=2*hatpsi/psi*gz_oben_u +psi^2*(hatgvv*gzu-hatguv*gzv+gvv*hatgzu-guv*hatgzv);
>     hatgx_oben_v:=2*hatpsi/psi*gx_oben_v +psi^2*(hatguu*gxv-hatguv*gxu+guu*hatgxv-guv*hatgxu);
>     hatgy_oben_v:=2*hatpsi/psi*gy_oben_v +psi^2*(hatguu*gyv-hatguv*gyu+guu*hatgyv-guv*hatgyu);
>     hatgz_oben_v:=2*hatpsi/psi*gz_oben_v +psi^2*(hatguu*gzv-hatguv*gzu+guu*hatgzv-guv*hatgzu);
> hatwr[3*ii-2]:=hatwx; hatwr[3*ii-1]:=hatwy; hatwr[3*ii]:=hatwz;
> hatobenu[3*ii-2]:=hatgx_oben_u; hatobenu[3*ii-1]:=hatgy_oben_u; hatobenu[3*ii]:=hatgz_oben_u;
> hatobenv[3*ii-2]:=hatgx_oben_v; hatobenv[3*ii-1]:=hatgy_oben_v; hatobenv[3*ii]:=hatgz_oben_v;
> hatpsir[ii]:=hatpsi;
> 
>     dhatwx:=hatgyu*dgzv-dgyv*hatgzu+dgyu*hatgzv-hatgyv*dgzu;
>     dhatwy:=dgxv*hatgzu-hatgxu*dgzv+hatgxv*dgzu-dgxu*hatgzv;
>     dhatwz:=hatgxu*dgyv-dgxv*hatgyu+dgxu*hatgyv-hatgxv*dgyu;
> dhatpsi:=3*dpsi*hatpsi/psi-psi^3*(dwx*hatwx+dwy*hatwy+dwz*hatwz+wx*dhatwx+wy*dhatwy+wz*dhatwz);
> 
> dhatguu:=2*(dgxu*hatgxu+dgyu*hatgyu+dgzu*hatgzu);
> dhatgvv:=2*(dgxv*hatgxv+dgyv*hatgyv+dgzv*hatgzv);
> dhatguv:=dgxu*hatgxv+dgyu*hatgyv+dgzu*hatgzv +dgxv*hatgxu+dgyv*hatgyu+dgzv*hatgzu;
> 
>     dhatgx_oben_u:=2*(psi*dhatpsi-3*dpsi*hatpsi)/psi^2*gx_oben_u
> +2*(dpsi*hatgx_oben_u+hatpsi*dgx_oben_u)/psi
>  +psi^2*(dhatgvv*gxu-dhatguv*gxv+hatgvv*dgxu-hatguv*dgxv+dgvv*hatgxu-dguv*hatgxv);
>     dhatgy_oben_u:=2*(psi*dhatpsi-3*dpsi*hatpsi)/psi^2*gy_oben_u
>  +2*(dpsi*hatgy_oben_u+hatpsi*dgy_oben_u)/psi
>  +psi^2*(dhatgvv*gyu-dhatguv*gyv+hatgvv*dgyu-hatguv*dgyv+dgvv*hatgyu-dguv*hatgyv);
>     dhatgz_oben_u:=2*(psi*dhatpsi-3*dpsi*hatpsi)/psi^2*gz_oben_u
> +2*(dpsi*hatgz_oben_u+hatpsi*dgz_oben_u)/psi
>  +psi^2*(dhatgvv*gzu-dhatguv*gzv+hatgvv*dgzu-hatguv*dgzv+dgvv*hatgzu-dguv*hatgzv);
>     dhatgx_oben_v:=2*(psi*dhatpsi-3*dpsi*hatpsi)/psi^2*gx_oben_v
> +2*(dpsi*hatgx_oben_v+hatpsi*dgx_oben_v)/psi
>  +psi^2*(dhatguu*gxv-dhatguv*gxu+hatguu*dgxv-hatguv*dgxu+dguu*hatgxv-dguv*hatgxu);
>     dhatgy_oben_v:=2*(psi*dhatpsi-3*dpsi*hatpsi)/psi^2*gy_oben_v
> +2*(dpsi*hatgy_oben_v+hatpsi*dgy_oben_v)/psi
>  +psi^2*(dhatguu*gyv-dhatguv*gyu+hatguu*dgyv-hatguv*dgyu+dguu*hatgyv-dguv*hatgyu);
>     dhatgz_oben_v:=2*(psi*dhatpsi-3*dpsi*hatpsi)/psi^2*gz_oben_v 
> +2*(dpsi*hatgz_oben_v+hatpsi*dgz_oben_v)/psi
> +psi^2*(dhatguu*gzv-dhatguv*gzu+hatguu*dgzv-hatguv*dgzu+dguu*hatgzv-dguv*hatgzu);
> 
> dhatwr[3*ii-2]:=dhatwx; dhatwr[3*ii-1]:=dhatwy; dhatwr[3*ii]:=dhatwz;
> dhatobenu[3*ii-2]:=dhatgx_oben_u; dhatobenu[3*ii-1]:=dhatgy_oben_u; dhatobenu[3*ii]:=dhatgz_oben_u;
> dhatobenv[3*ii-2]:=dhatgx_oben_v; dhatobenv[3*ii-1]:=dhatgy_oben_v; dhatobenv[3*ii]:=dhatgz_oben_v;
> dhatpsir[ii]:=dhatpsi;
> od
> od:
> 
# Corner nodes:
> for j in [1,n2i+1] do 
> for i in [1,n1i+1] do
> if i=1 and j=1 then u:=0; v:=0; iii:=1 fi:
> if i=n1i+1 and j=1 then u:=n; v:=0; iii:=n1i fi:
> if i=1 and j=n2i+1 then u:=0; v:=m; iii:=(n2i-1)*n1i+1 fi:
> if i=n1i+1 and j=n2i+1 then u:=n; v:=m; iii:=n2i*n1i fi:
> ii:=(j-1)*(n1i+1)+i;
> zeta:=a*theta[iii]; 
> gxutilde:=((xn0-x00)*(1-v/m)+(xnm-x0m)*v/m)/n  +zeta*mmxu;
> gyutilde:= ((yn0-y00)*(1-v/m)+(ynm-y0m)*v/m)/n +zeta*mmyu;
> gzutilde:=zeta*mmzu;
> gxvtilde:=((x0m-x00)*(1-u/n)+(xnm-xn0)*u/n)/m+zeta*mmxv;
> gyvtilde:=((y0m-y00)*(1-u/n)+(ynm-yn0)*u/n)/m+zeta*mmyv;
> gzvtilde:=zeta*mmzv;
> 
> gxucrossm:=gyutilde*mmz-gzutilde*mmy;
> gyucrossm:=gzutilde*mmx-gxutilde*mmz;
> gzucrossm:=gxutilde*mmy-gyutilde*mmx;
> gxvcrossm:=gyvtilde*mmz-gzvtilde*mmy;
> gyvcrossm:=gzvtilde*mmx-gxvtilde*mmz;
> gzvcrossm:=gxvtilde*mmy-gyvtilde*mmx;
> wtildem:=gxutilde*gxvcrossm+gyutilde*gyvcrossm+gzutilde*gzvcrossm;
> psi:=1/wtildem;
> gx_oben_u:=psi*gxvcrossm;
> gy_oben_u:=psi*gyvcrossm;
> gz_oben_u:=psi*gzvcrossm;
> gx_oben_v:=-psi*gxucrossm;
> gy_oben_v:=-psi*gyucrossm;
> gz_oben_v:=-psi*gzucrossm;
> wx:=wtildem*mmx;
> wy:=wtildem*mmy;
> wz:=wtildem*mmz;
> wr[3*ii-2]:=wx; wr[3*ii-1]:=wy; wr[3*ii]:=wz;
> psir[ii]:=psi;
> obenu[3*ii-2]:=gx_oben_u; obenu[3*ii-1]:=gy_oben_u; obenu[3*ii]:=gz_oben_u;
> obenv[3*ii-2]:=gx_oben_v; obenv[3*ii-1]:=gy_oben_v; obenv[3*ii]:=gz_oben_v;
> 
> dgxu:=a*dtheta[iii]*mmxu;
> dgyu:= a*dtheta[iii]*mmyu;
> dgzu:=a*dtheta[iii]*mmzu;
> dgxv:=a*dtheta[iii]*mmxv;
> dgyv:=a*dtheta[iii]*mmyv;
> dgzv:=a*dtheta[iii]*mmzv;
> dgxucrossm:=dgyu*mmz-dgzu*mmy;
> dgyucrossm:=dgzu*mmx-dgxu*mmz;
> dgzucrossm:=dgxu*mmy-dgyu*mmx;
> dgxvcrossm:=dgyv*mmz-dgzv*mmy;
> dgyvcrossm:=dgzv*mmx-dgxv*mmz;
> dgzvcrossm:=dgxv*mmy-dgyv*mmx;
> dwtildem:=dgxu*gxvcrossm+dgyu*gyvcrossm+dgzu*gzvcrossm
> +gxutilde*dgxvcrossm+gyutilde*dgyvcrossm+gzutilde*dgzvcrossm;
> dpsi:=-psi^2*dwtildem;
> dwx:=dwtildem*mmx;
> dwy:=dwtildem*mmy;
> dwz:=dwtildem*mmz;
> dgx_oben_u:=dpsi*gxvcrossm+psi*dgxvcrossm;
> dgy_oben_u:=dpsi*gyvcrossm+psi*dgyvcrossm;
> dgz_oben_u:=dpsi*gzvcrossm+psi*dgzvcrossm;
> dgx_oben_v:=-dpsi*gxucrossm-psi*dgxucrossm;
> dgy_oben_v:=-dpsi*gyucrossm-psi*dgyucrossm;
> dgz_oben_v:=-dpsi*gzucrossm-psi*dgzucrossm;
> 
> dwr[3*ii-2]:=dwx; dwr[3*ii-1]:=dwy; dwr[3*ii]:=dwz;
> dobenu[3*ii-2]:=dgx_oben_u; dobenu[3*ii-1]:=dgy_oben_u; dobenu[3*ii]:=dgz_oben_u;
> dobenv[3*ii-2]:=dgx_oben_v; dobenv[3*ii-1]:=dgy_oben_v; dobenv[3*ii]:=dgz_oben_v;
> dpsir[ii]:=dpsi;
> 
> hatgxu:=a*hattheta[iii]*mmxu;
> hatgyu:= a*hattheta[iii]*mmyu;
> hatgzu:=a*hattheta[iii]*mmzu;
> hatgxv:=a*hattheta[iii]*mmxv;
> hatgyv:=a*hattheta[iii]*mmyv;
> hatgzv:=a*hattheta[iii]*mmzv;
> hatgxucrossm:=hatgyu*mmz-hatgzu*mmy;
> hatgyucrossm:=hatgzu*mmx-hatgxu*mmz;
> hatgzucrossm:=hatgxu*mmy-hatgyu*mmx;
> hatgxvcrossm:=hatgyv*mmz-hatgzv*mmy;
> hatgyvcrossm:=hatgzv*mmx-hatgxv*mmz;
> hatgzvcrossm:=hatgxv*mmy-hatgyv*mmx;
> hatwtildem:=hatgxu*gxvcrossm+hatgyu*gyvcrossm+hatgzu*gzvcrossm
> +gxutilde*hatgxvcrossm+gyutilde*hatgyvcrossm+gzutilde*hatgzvcrossm;
> hatpsi:=-psi^2*hatwtildem;
> hatwx:=hatwtildem*mmx;
> hatwy:=hatwtildem*mmy;
> hatwz:=hatwtildem*mmz;
> hatgx_oben_u:=hatpsi*gxvcrossm+psi*hatgxvcrossm;
> hatgy_oben_u:=hatpsi*gyvcrossm+psi*hatgyvcrossm;
> hatgz_oben_u:=hatpsi*gzvcrossm+psi*hatgzvcrossm;
> hatgx_oben_v:=-hatpsi*gxucrossm-psi*hatgxucrossm;
> hatgy_oben_v:=-hatpsi*gyucrossm-psi*hatgyucrossm;
> hatgz_oben_v:=-hatpsi*gzucrossm-psi*hatgzucrossm;
> 
> hatwr[3*ii-2]:=hatwx; hatwr[3*ii-1]:=hatwy; hatwr[3*ii]:=hatwz;
> hatobenu[3*ii-2]:=hatgx_oben_u; hatobenu[3*ii-1]:=hatgy_oben_u; hatobenu[3*ii]:=hatgz_oben_u;
> hatobenv[3*ii-2]:=hatgx_oben_v; hatobenv[3*ii-1]:=hatgy_oben_v; hatobenv[3*ii]:=hatgz_oben_v;
> hatpsir[ii]:=hatpsi;
> 
> dhatwtildem:=hatgxu*dgxvcrossm+hatgyu*dgyvcrossm+hatgzu*dgzvcrossm+dgxu*hatgxvcrossm+dgyu*hatgyvcrossm+dgzu*hatgzvcrossm;
> dhatpsi:=2*dpsi*hatpsi/psi-psi^2*dhatwtildem;
> dhatwx:=dhatwtildem*mmx;
> dhatwy:=dhatwtildem*mmy;
> dhatwz:=dhatwtildem*mmz;
> dhatgx_oben_u:=dhatpsi*gxvcrossm+hatpsi*dgxvcrossm+dpsi*hatgxvcrossm;
> dhatgy_oben_u:=dhatpsi*gyvcrossm+hatpsi*dgyvcrossm+dpsi*hatgyvcrossm;
> dhatgz_oben_u:=dhatpsi*gzvcrossm+hatpsi*dgzvcrossm+dpsi*hatgzvcrossm;
> dhatgx_oben_v:=-dhatpsi*gxucrossm-hatpsi*dgxucrossm-dpsi*hatgxucrossm;
> dhatgy_oben_v:=-dhatpsi*gyucrossm-hatpsi*dgyucrossm-dpsi*hatgyucrossm;
> dhatgz_oben_v:=-dhatpsi*gzucrossm-hatpsi*dgzucrossm-dpsi*hatgzucrossm;
> 
> dhatwr[3*ii-2]:=dhatwx; dhatwr[3*ii-1]:=dhatwy; dhatwr[3*ii]:=dhatwz;
> dhatobenu[3*ii-2]:=dhatgx_oben_u; dhatobenu[3*ii-1]:=dhatgy_oben_u; dhatobenu[3*ii]:=dhatgz_oben_u;
> dhatobenv[3*ii-2]:=dhatgx_oben_v; dhatobenv[3*ii-1]:=dhatgy_oben_v; dhatobenv[3*ii]:=dhatgz_oben_v;
> dhatpsir[ii]:=dhatpsi;
> 
> od
> od:
> 
> 
# Nodes within the boundaries:
> 
# Lower boundary:
> 
> if sphere or cylinder or tube6 or tube3 then
>  e:=vector([0,-1,0])
> fi;
> if tube4 then
>  e:=vector([-0.5/sqrt(3.),-sqrt(2./3.),-0.5])
> fi;
> for i from 2 to n1i do
> u:=i-1.5; v:=0;
> ii:=i;
> iii:=i;
> boundv():
> od:
> 
# Upper boundary:
> if sphere then
>  e:=vector([0,sqrt(.5),sqrt(.5)])
> fi;
> if cylinder then
>  e:=vector([0,1.,0])
> fi;
> if tube6 then
>  e:=vector([0,sqrt(.5),-sqrt(.5)])
> fi; 
> if tube3 then
>  e:=vector([-0.5*sqrt(0.5),0.5*sqrt(3.),-0.5*sqrt(0.5)])
> fi;
> if tube4 then
>  e:=vector([0.5/sqrt(3.),sqrt(2./3.),-0.5])
> fi;
> for i from 2 to n1i do
> u:=i-1.5; v:=m;
> ii:=(n1i+1)*n2i+i;
> iii:=n1i*(n2i-1)+i;
> boundv():
> od:
> 
# Left boundary:
> 
> if sphere or cylinder or tube6 then
>  e:=vector([-1,0,0])
> fi;
> if tube3 then
>  e:=vector([-sqrt(0.5),0,sqrt(0.5)])
> fi;
> if tube4 then
>  e:=vector([-0.5*sqrt(3.),0,0.5])
> fi;
> for j from 2 to n2i do
> u:=0; v:=j-1.5;
> ii:=(j-1)*(n1i+1)+1;
> iii:=(j-1)*n1i+1; 
> boundu():
> od:
# Right boundary:
> if sphere or cylinder or tube6 or tube3 then
>  e:=vector([sqrt(.5),0,sqrt(.5)])
> fi;
> if tube4 then
>  e:=vector([0.5*sqrt(3.),0,0.5])
> fi;
> for j from 2 to n2i do
> u:=n; v:=j-1.5;
> ii:=j*(n1i+1); 
> iii:=j*n1i;
> boundu()
> od:
> 
# Loop over the points of integration:
> u:='u': v:='v':
> wc:=0:
> ac:=0:
> a1c:=0:
> a2c:=0:
> intc:=0:
> volu:=0:

> dwc:=0:
> dac:=0:
> da1c:=0:
> da2c:=0:
> dvolu:=0:

> dhatwc:=0:
> dhatac:=0:
> dhata1c:=0:
> dhata2c:=0:
> dhatvolu:=0:
> 
> a_dwc:=0:
> a_dac:=0:
> a_dvoc:=0:
> a_dvwc:=0:
> atimesa_intc:=0:
> 
> zi1:=0:
> ziq1:=0:
> zi2:=0:
> ziq2:=0:
> zi3:=0:
> ziq3:=0:
> zi4:=0:
> ziq4:=0:
> zi5:=0:
> ziq5:=0:
> zi6:=0:
> ziq6:=0:
> zi7:=0:
> ziq7:=0:
> zi8:=0:
> ziq8:=0:
> zi9:=0:
> ziq9:=0:
> 
> tbetrseq:=NULL:
> mbetrseq:=NULL:
# 
> for j to n2i do
> v:=j-1.;
> wt2:=1.;
> if j=1 then v:=0.25; wt2:=0.5 fi;
> if j=n2i then v:=n2i-1.25; wt2:=0.5 fi;
> out1[j]:=NULL:
> out2[j]:=NULL:
> out3[j]:=NULL:
> out4[j]:=NULL:
> out5[j]:=NULL:
> out6[j]:=NULL:
> out7[j]:=NULL:
> out8[j]:=NULL:
> out9[j]:=NULL:
> 
> for i to n1i do
> u:=i-1.;
> wt1:=1.;
> if i=1 then u:=0.25; wt1:=0.5 fi;
> if i=n1i then u:=n1i-1.25; wt1:=0.5 fi; 
> 
> if i=n1i then print(i,j) fi;
> 
> ii:=(j-1)*(n1i+1)+i;
>      
> gx_oben_u:=0.25*(obenu[3*ii+3*n1i+4] +obenu[3*ii+1] +obenu[3*ii+3*n1i+1] +obenu[3*ii-2]);
> gy_oben_u:=0.25*(obenu[3*ii+3*n1i+5] +obenu[3*ii+2] +obenu[3*ii+3*n1i+2] +obenu[3*ii-1]);
> gz_oben_u:=0.25*(obenu[3*ii+3*n1i+6] +obenu[3*ii+3] +obenu[3*ii+3*n1i+3] +obenu[3*ii]);
> gx_oben_v:=0.25*(obenv[3*ii+3*n1i+4] +obenv[3*ii+1] +obenv[3*ii+3*n1i+1] +obenv[3*ii-2]);
> gy_oben_v:=0.25*(obenv[3*ii+3*n1i+5] +obenv[3*ii+2] +obenv[3*ii+3*n1i+2] +obenv[3*ii-1]);
> gz_oben_v:=0.25*(obenv[3*ii+3*n1i+6] +obenv[3*ii+3] +obenv[3*ii+3*n1i+3] +obenv[3*ii]);
>   psi:=0.25*(psir[ii+n1i+2] +psir[ii+1] +psir[ii+n1i+1] +psir[ii]);
> wxu:=0.5*(wr[3*ii+3*n1i+4] +wr[3*ii+1] -wr[3*ii+3*n1i+1] -wr[3*ii-2])/wt1;
> wyu:=0.5*(wr[3*ii+3*n1i+5] +wr[3*ii+2] -wr[3*ii+3*n1i+2] -wr[3*ii-1])/wt1;
> wzu:=0.5*(wr[3*ii+3*n1i+6] +wr[3*ii+3] -wr[3*ii+3*n1i+3] -wr[3*ii])/wt1;
> wxv:=0.5*(wr[3*ii+3*n1i+4] -wr[3*ii+1] +wr[3*ii+3*n1i+1] -wr[3*ii-2])/wt2;
> wyv:=0.5*(wr[3*ii+3*n1i+5] -wr[3*ii+2] +wr[3*ii+3*n1i+2] -wr[3*ii-1])/wt2;
> wzv:=0.5*(wr[3*ii+3*n1i+6] -wr[3*ii+3] +wr[3*ii+3*n1i+3] -wr[3*ii])/wt2;
>      cuu:=-psi*(gx_oben_u*wxu+gy_oben_u*wyu+gz_oben_u*wzu);
>      cuv:=-psi*(gx_oben_u*wxv+gy_oben_u*wyv+gz_oben_u*wzv);
>      cvu:=-psi*(gx_oben_v*wxu+gy_oben_v*wyu+gz_oben_v*wzu);
>      cvv:=-psi*(gx_oben_v*wxv+gy_oben_v*wyv+gz_oben_v*wzv);
>      inv1:=cuu+cvv;
>      inv2:=cuu*cvv-cuv*cvu;
>      hi10:=max(0.,inv1^2/4-inv2);
>      c1:=inv1/2+sqrt(hi10);
>      c2:=inv1-c1;
> 
> accontr:=1./psi*wt1*wt2;
> wc:=wc+phi(inv1,inv2)*accontr;
> ac:=ac+accontr;
> a1c:=a1c+inv1*accontr;
> a2c:=a2c+inv2*accontr;
> hi5:=phitr(inv1,inv2)*inv1+2*phidt(inv1,inv2)*inv2;
> intc:=intc+hi5*accontr;
> 
> zeta:=a*theta[(j-1)*n1i+i];
> dzeta:=a*dtheta[(j-1)*n1i+i];
> hatzeta:=a*hattheta[(j-1)*n1i+i];
> 
> if i=1 and j>1 and j<n1i then 
>   zeta:=a*(0.75*theta[(j-1)*n1i+1]+0.25*theta[(j-1)*n1i+2]);
>   dzeta:=a*(0.75*dtheta[(j-1)*n1i+1]+0.25*dtheta[(j-1)*n1i+2]);
>   hatzeta:=a*(0.75*hattheta[(j-1)*n1i+1]+0.25*hattheta[(j-1)*n1i+2]);
> fi;
> 
> if i=n1i and j>1 and j<n2i then 
>   zeta:=a*(0.75*theta[j*n1i]+0.25*theta[j*n1i-1]); 
>   dzeta:=a*(0.75*dtheta[j*n1i]+0.25*dtheta[j*n1i-1]);
>   hatzeta:=a*(0.75*hattheta[j*n1i]+0.25*hattheta[j*n1i-1]);
> fi;
> 
> if j=1 and i>1 and i<n1i then 
>   zeta:=a*(0.75*theta[i]+0.25*theta[n1i+i]);
>   dzeta:=a*(0.75*dtheta[i]+0.25*dtheta[n1i+i]);
>   hatzeta:=a*(0.75*hattheta[i]+0.25*hattheta[n1i+i]);
> fi;
> 
> if j=n2i and i>1 and i<n1i then 
>   zeta:=a*(0.75*theta[(n2i-1)*n1i+i]+0.25*theta[(n2i-2)*n1i+i]);
>   dzeta:=a*(0.75*dtheta[(n2i-1)*n1i+i]+0.25*dtheta[(n2i-2)*n1i+i]);
>   hatzeta:=a*(0.75*hattheta[(n2i-1)*n1i+i]+0.25*hattheta[(n2i-2)*n1i+i]);
> fi;
> 
> if i=1 and j=1 then
>  zeta:=a*(0.5625*theta[1]+0.1875*(theta[2]+theta[n1i+1])+0.0625*theta[n1i+2]);
>  dzeta:=a*(0.5625*dtheta[1]+0.1875*(dtheta[2]+dtheta[n1i+1])+0.0625*dtheta[n1i+2]);
>  hatzeta:=a*(0.5625*hattheta[1]+0.1875*(hattheta[2]+hattheta[n1i+1])+0.0625*hattheta[n1i+2]);
> fi;
> 
> if i=n1i and j=1 then
>  zeta:=a*(0.5625*theta[n1i]+0.1875*(theta[n1i-1]+theta[2*n1i])+0.0625*theta[2*n1i-1]);
>  dzeta:=a*(0.5625*dtheta[n1i]+0.1875*(dtheta[n1i-1]+dtheta[2*n1i])+0.0625*dtheta[2*n1i-1]);
>  hatzeta:=a*(0.5625*hattheta[n1i]+0.1875*(hattheta[n1i-1]+hattheta[2*n1i])+0.0625*hattheta[2*n1i-1]);
> fi;
> 
> if i=1 and j=n2i then
>  zeta:=a*(0.5625*theta[(n2i-1)*n1i+1]+0.1875*(theta[(n2i-1)*n1i+2]+theta[(n2i-2)*n1i+1])+0.0625*theta[(n2i-2)*n1i+2]);
>  dzeta:=a*(0.5625*dtheta[(n2i-1)*n1i+1]+0.1875*(dtheta[(n2i-1)*n1i+2]+dtheta[(n2i-2)*n1i+1])+0.0625*dtheta[(n2i-2)*n1i+2]);
>  hatzeta:=a*(0.5625*hattheta[(n2i-1)*n1i+1]+0.1875*(hattheta[(n2i-1)*n1i+2]+hattheta[(n2i-2)*n1i+1])+0.0625*hattheta[(n2i-2)*n1i+2]);
> fi;
> 
> if i=n1i and j=n2i then
>  zeta:=a*(0.5625*theta[nn]+0.1875*(theta[nn-1]+theta[(n2i-1)*n1i])+0.0625*theta[(n2i-1)*n1i-1]);
>  dzeta:=a*(0.5625*dtheta[nn]+0.1875*(dtheta[nn-1]+dtheta[(n2i-1)*n1i])+0.0625*dtheta[(n2i-1)*n1i-1]);
>  hatzeta:=a*(0.5625*hattheta[nn]+0.1875*(hattheta[nn-1]+hattheta[(n2i-1)*n1i])+0.0625*hattheta[(n2i-1)*n1i-1]);
> fi;
> 
>    higxu:=(xn0-x00)/n*(1-v/m)+(xnm-x0m)/n*v/m +zeta*mmxu/2.;
>    higyu:=(yn0-y00)/n*(1-v/m)+(ynm-y0m)/n*v/m +zeta*mmyu/2.;    
>    higzu:=zeta*mmzu/2.;
>    higxv:=(x0m-x00)*(1-u/n)/m+(xnm-xn0)*u/n/m+zeta*mmxv/2.;
>    higyv:=(y0m-y00)*(1-u/n)/m+(ynm-yn0)*u/n/m+zeta*mmyv/2.;    
>    higzv:=zeta*mmzv/2.;
>     hiwx:=(higyu*higzv-higyv*higzu)*zeta
>             +(mmyu*mmzv-mmyv*mmzu)*zeta^3/12.;
>     hiwy:=(higxv*higzu-higxu*higzv)*zeta
>             +(mmzu*mmxv-mmzv*mmxu)*zeta^3/12.;
>     hiwz:=(higxu*higyv-higxv*higyu)*zeta
>            +(mmxu*mmyv-mmxv*mmyu)*zeta^3/12.;
> 
> volu:=volu+(hiwx*mmx+hiwy*mmy+hiwz*mmz)*wt1*wt2;
> 
>  
> 
>  mtensuu:=-phitr(inv1,inv2)-phidt(inv1,inv2)*cvv;
>  mtensuv:=phidt(inv1,inv2)*cuv;
>  mtensvu:=phidt(inv1,inv2)*cvu;
>  mtensvv:=-phitr(inv1,inv2)-phidt(inv1,inv2)*cuu;
>  
>  hi11:=phi(inv1,inv2)-phidt(inv1,inv2)*inv2-f*kf;
>  ttensuu:=hi11 -phitr(inv1,inv2)*cuu;
>  ttensuv:=-phitr(inv1,inv2)*cuv;
>  ttensvu:=-phitr(inv1,inv2)*cvu;
>  ttensvv:=hi11 -phitr(inv1,inv2)*cvv;
> 
> memb:=ttensuu*cuu+ttensuv*cvu+ttensvu*cuv+ttensvv*cvv;
> memb2:=(-mtensuu*cuu-mtensuv*cvu)*cuu+(-mtensuu*cuv-mtensuv*cvv)*cvu
>   +(-mtensvu*cuu-mtensvv*cvu)*cuv+(-mtensvu*cuv-mtensvv*cvv)*cvv;
> memba:=phi(inv1,inv2);
> pn:=-po*(1-inv1*hh/2+inv2*hh^2/4)+pw*(1+inv1*hh/2+inv2*hh^2/4);
> memd:=memb+pn; 
> memx:=(memba-wcalt/acalt)*inv1;
> memy:=(po*vocalt+pw*vwcalt)/acalt*inv1;
> 
> ####################
> 
> gupuu:=gx_oben_u*gx_oben_u+gy_oben_u*gy_oben_u+gz_oben_u*gz_oben_u;
> gupuv:=gx_oben_u*gx_oben_v+gy_oben_u*gy_oben_v+gz_oben_u*gz_oben_v;
> gupvv:=gx_oben_v*gx_oben_v+gy_oben_v*gy_oben_v+gz_oben_v*gz_oben_v;
> 
> ttensupuv:=gupuu*ttensvu+gupuv*ttensvv;
> ttensupvu:=gupuv*ttensuu+gupvv*ttensuv;
> ttensupuu:=gupuu*ttensuu+gupuv*ttensuv;
> ttensupvv:=gupuv*ttensvu+gupvv*ttensvv;
> intttens:=max(0.,ttensuu^2+2*ttensuv*ttensvu+ttensvv^2);
> tbetrseq:=tbetrseq,sqrt(intttens);
> 
> mtensupuv:=gupuu*mtensvu+gupuv*mtensvv;
> mtensupvu:=gupuv*mtensuu+gupvv*mtensuv;
> mtensupuu:=gupuu*mtensuu+gupuv*mtensuv;
> mtensupvv:=gupuv*mtensvu+gupvv*mtensvv;
> intmtens:=max(0.,mtensuu^2+2*mtensuv*mtensvu+mtensvv^2);
> mbetrseq:=mbetrseq,sqrt(intmtens);
> 
> 
> ########################
> 
> out1[j]:=out1[j],[u/n,(c1+c2)*0.5];
> out2[j]:=out2[j],[u/n,(c1-c2)*0.5];
> out3[j]:=out3[j],[u/n,memx];
> out4[j]:=out4[j],[u/n,memy];
> out5[j]:=out5[j],[u/n,-memb2];
> out6[j]:=out6[j],[u/n,memb]:
> out7[j]:=out7[j],[u/n,pn];
> out8[j]:=out8[j],[u/n,-memd];
> out9[j]:=out9[j],[u/n,memba];
> 
> zi1:=zi1+(c1+c2)*0.5/psi*wt1*wt2;
> zi2:=zi2+(c1-c2)*0.5/psi*wt1*wt2;
> zi3:=zi3+memx/psi*wt1*wt2;
> zi4:=zi4+memy/psi*wt1*wt2;
> zi5:=zi5-memb2/psi*wt1*wt2;
> zi6:=zi6+memb/psi*wt1*wt2;
> zi7:=zi7+pn/psi*wt1*wt2;
> zi8:=zi8-memd/psi*wt1*wt2;
> zi9:=zi9+memba/psi*wt1*wt2;
> ziq1:=ziq1+(c1+c2)^2*0.25/psi*wt1*wt2;
> ziq2:=ziq2+(c1-c2)^2*0.25/psi*wt1*wt2;
> ziq3:=ziq3+memx^2/psi*wt1*wt2;
> ziq4:=ziq4+memy^2/psi*wt1*wt2;
> ziq5:=ziq5+memb2^2/psi*wt1*wt2;
> ziq6:=ziq6+memb^2/psi*wt1*wt2;
> ziq7:=ziq7+pn^2/psi*wt1*wt2;
> ziq8:=ziq8+memd^2/psi*wt1*wt2;
> ziq9:=ziq9+memba^2/psi*wt1*wt2;
> 
> dgx_oben_u:=0.25*(dobenu[3*ii+3*n1i+4] +dobenu[3*ii+1] +dobenu[3*ii+3*n1i+1] +dobenu[3*ii-2]);
> dgy_oben_u:=0.25*(dobenu[3*ii+3*n1i+5] +dobenu[3*ii+2] +dobenu[3*ii+3*n1i+2] +dobenu[3*ii-1]);
> dgz_oben_u:=0.25*(dobenu[3*ii+3*n1i+6] +dobenu[3*ii+3] +dobenu[3*ii+3*n1i+3] +dobenu[3*ii]);
> dgx_oben_v:=0.25*(dobenv[3*ii+3*n1i+4] +dobenv[3*ii+1] +dobenv[3*ii+3*n1i+1] +dobenv[3*ii-2]);
> dgy_oben_v:=0.25*(dobenv[3*ii+3*n1i+5] +dobenv[3*ii+2] +dobenv[3*ii+3*n1i+2] +dobenv[3*ii-1]);
> dgz_oben_v:=0.25*(dobenv[3*ii+3*n1i+6] +dobenv[3*ii+3] +dobenv[3*ii+3*n1i+3] +dobenv[3*ii]);
>   dpsi:=0.25*(dpsir[ii+n1i+2] +dpsir[ii+1] +dpsir[ii+n1i+1] +dpsir[ii]);
> dwxu:=0.5*(dwr[3*ii+3*n1i+4] +dwr[3*ii+1] -dwr[3*ii+3*n1i+1] -dwr[3*ii-2])/wt1;
> dwyu:=0.5*(dwr[3*ii+3*n1i+5] +dwr[3*ii+2] -dwr[3*ii+3*n1i+2] -dwr[3*ii-1])/wt1;
> dwzu:=0.5*(dwr[3*ii+3*n1i+6] +dwr[3*ii+3] -dwr[3*ii+3*n1i+3] -dwr[3*ii])/wt1;
> dwxv:=0.5*(dwr[3*ii+3*n1i+4] -dwr[3*ii+1] +dwr[3*ii+3*n1i+1] -dwr[3*ii-2])/wt2;
> dwyv:=0.5*(dwr[3*ii+3*n1i+5] -dwr[3*ii+2] +dwr[3*ii+3*n1i+2] -dwr[3*ii-1])/wt2;
> dwzv:=0.5*(dwr[3*ii+3*n1i+6] -dwr[3*ii+3] +dwr[3*ii+3*n1i+3] -dwr[3*ii])/wt2;
>      dcuu:=dpsi/psi*cuu-psi*(dgx_oben_u*wxu+dgy_oben_u*wyu+dgz_oben_u*wzu+gx_oben_u*dwxu+gy_oben_u*dwyu+gz_oben_u*dwzu);
>      dcuv:=dpsi/psi*cuv-psi*(dgx_oben_u*wxv+dgy_oben_u*wyv+dgz_oben_u*wzv+gx_oben_u*dwxv+gy_oben_u*dwyv+gz_oben_u*dwzv);
>      dcvu:=dpsi/psi*cvu-psi*(dgx_oben_v*wxu+dgy_oben_v*wyu+dgz_oben_v*wzu+gx_oben_v*dwxu+gy_oben_v*dwyu+gz_oben_v*dwzu);
>      dcvv:=dpsi/psi*cvv-psi*(dgx_oben_v*wxv+dgy_oben_v*wyv+dgz_oben_v*wzv+gx_oben_v*dwxv+gy_oben_v*dwyv+gz_oben_v*dwzv);
> dinv1:=dcuu+dcvv;
> dinv2:=dcuu*cvv+cuu*dcvv-dcuv*cvu-cuv*dcvu;
> dphi:=phitr(inv1,inv2)*dinv1+phidt(inv1,inv2)*dinv2;
> daccontr:=-dpsi/psi^2*wt1*wt2;
> dac:=dac+daccontr;
> dwccontr:=phi(inv1,inv2)*daccontr+dphi*accontr;
> dwc:=dwc+dwccontr;
> da1ccontr:=inv1*daccontr+dinv1*accontr;
> da1c:=da1c+da1ccontr;
> da2ccontr:=inv2*daccontr+dinv2*accontr;
> da2c:=da2c+da2ccontr;
> 
>    dhigxu:=dzeta*mmxu/2.;
>    dhigyu:=dzeta*mmyu/2.;    
>    dhigzu:=dzeta*mmzu/2.;
>    dhigxv:=dzeta*mmxv/2.;
>    dhigyv:=dzeta*mmyv/2.;    
>    dhigzv:=dzeta*mmzv/2.;    
>     dhiwx:=(dhigyu*higzv-higyv*dhigzu+higyu*dhigzv-dhigyv*higzu)*zeta
>           +(higyu*higzv-higyv*higzu+(mmyu*mmzv-mmyv*mmzu)*zeta^2/4.)*dzeta;
>     dhiwy:=(higxv*dhigzu-dhigxu*higzv+dhigxv*higzu-higxu*dhigzv)*zeta
>           +(higxv*higzu-higxu*higzv+(mmzu*mmxv-mmzv*mmxu)*zeta^2/4.)*dzeta;
>     dhiwz:=(dhigxu*higyv-higxv*dhigyu+higxu*dhigyv-dhigxv*higyu)*zeta
>           +(higxu*higyv-higxv*higyu+(mmxu*mmyv-mmxv*mmyu)*zeta^2/4.)*dzeta;
> dvolucontr:=(dhiwx*mmx+dhiwy*mmy+dhiwz*mmz)*wt1*wt2;
> dvolu:=dvolu+dvolucontr;  
> 
> hatgx_oben_u:=0.25*(hatobenu[3*ii+3*n1i+4] +hatobenu[3*ii+1] +hatobenu[3*ii+3*n1i+1] +hatobenu[3*ii-2]);
> hatgy_oben_u:=0.25*(hatobenu[3*ii+3*n1i+5] +hatobenu[3*ii+2] +hatobenu[3*ii+3*n1i+2] +hatobenu[3*ii-1]);
> hatgz_oben_u:=0.25*(hatobenu[3*ii+3*n1i+6] +hatobenu[3*ii+3] +hatobenu[3*ii+3*n1i+3] +hatobenu[3*ii]);
> hatgx_oben_v:=0.25*(hatobenv[3*ii+3*n1i+4] +hatobenv[3*ii+1] +hatobenv[3*ii+3*n1i+1] +hatobenv[3*ii-2]);
> hatgy_oben_v:=0.25*(hatobenv[3*ii+3*n1i+5] +hatobenv[3*ii+2] +hatobenv[3*ii+3*n1i+2] +hatobenv[3*ii-1]);
> hatgz_oben_v:=0.25*(hatobenv[3*ii+3*n1i+6] +hatobenv[3*ii+3] +hatobenv[3*ii+3*n1i+3] +hatobenv[3*ii]);
> 
> hatpsi:=0.25*(hatpsir[ii+n1i+2] +hatpsir[ii+1] +hatpsir[ii+n1i+1] +hatpsir[ii]);
> 
> hatwxu:=0.5*(hatwr[3*ii+3*n1i+4] +hatwr[3*ii+1] -hatwr[3*ii+3*n1i+1] -hatwr[3*ii-2])/wt1;
> hatwyu:=0.5*(hatwr[3*ii+3*n1i+5] +hatwr[3*ii+2] -hatwr[3*ii+3*n1i+2] -hatwr[3*ii-1])/wt1;
> hatwzu:=0.5*(hatwr[3*ii+3*n1i+6] +hatwr[3*ii+3] -hatwr[3*ii+3*n1i+3] -hatwr[3*ii])/wt1;
> hatwxv:=0.5*(hatwr[3*ii+3*n1i+4] -hatwr[3*ii+1] +hatwr[3*ii+3*n1i+1] -hatwr[3*ii-2])/wt2;
> hatwyv:=0.5*(hatwr[3*ii+3*n1i+5] -hatwr[3*ii+2] +hatwr[3*ii+3*n1i+2] -hatwr[3*ii-1])/wt2;
> hatwzv:=0.5*(hatwr[3*ii+3*n1i+6] -hatwr[3*ii+3] +hatwr[3*ii+3*n1i+3] -hatwr[3*ii])/wt2;
>      hatcuu:=hatpsi/psi*cuu-psi*(hatgx_oben_u*wxu+hatgy_oben_u*wyu+hatgz_oben_u*wzu
>                 +gx_oben_u*hatwxu+gy_oben_u*hatwyu+gz_oben_u*hatwzu);
>      hatcuv:=hatpsi/psi*cuv-psi*(hatgx_oben_u*wxv+hatgy_oben_u*wyv+hatgz_oben_u*wzv
>                 +gx_oben_u*hatwxv+gy_oben_u*hatwyv+gz_oben_u*hatwzv);
>      hatcvu:=hatpsi/psi*cvu-psi*(hatgx_oben_v*wxu+hatgy_oben_v*wyu+hatgz_oben_v*wzu
>                  +gx_oben_v*hatwxu+gy_oben_v*hatwyu+gz_oben_v*hatwzu);
>      hatcvv:=hatpsi/psi*cvv-psi*(hatgx_oben_v*wxv+hatgy_oben_v*wyv+hatgz_oben_v*wzv
>                  +gx_oben_v*hatwxv+gy_oben_v*hatwyv+gz_oben_v*hatwzv);
> hatinv1:=hatcuu+hatcvv;
> hatinv2:=hatcuu*cvv+cuu*hatcvv-hatcuv*cvu-cuv*hatcvu;
> hatphi:=phitr(inv1,inv2)*hatinv1+phidt(inv1,inv2)*hatinv2;
> 
> dhatgx_oben_u:=0.25*(dhatobenu[3*ii+3*n1i+4] +dhatobenu[3*ii+1] +dhatobenu[3*ii+3*n1i+1] +dhatobenu[3*ii-2]);
> dhatgy_oben_u:=0.25*(dhatobenu[3*ii+3*n1i+5] +dhatobenu[3*ii+2] +dhatobenu[3*ii+3*n1i+2] +dhatobenu[3*ii-1]);
> dhatgz_oben_u:=0.25*(dhatobenu[3*ii+3*n1i+6] +dhatobenu[3*ii+3] +dhatobenu[3*ii+3*n1i+3] +dhatobenu[3*ii]);
> dhatgx_oben_v:=0.25*(dhatobenv[3*ii+3*n1i+4] +dhatobenv[3*ii+1] +dhatobenv[3*ii+3*n1i+1] +dhatobenv[3*ii-2]);
> dhatgy_oben_v:=0.25*(dhatobenv[3*ii+3*n1i+5] +dhatobenv[3*ii+2] +dhatobenv[3*ii+3*n1i+2] +dhatobenv[3*ii-1]);
> dhatgz_oben_v:=0.25*(dhatobenv[3*ii+3*n1i+6] +dhatobenv[3*ii+3] +dhatobenv[3*ii+3*n1i+3] +dhatobenv[3*ii]);
> 
> dhatpsi:=0.25*(dhatpsir[ii+n1i+2] +dhatpsir[ii+1] +dhatpsir[ii+n1i+1] +dhatpsir[ii]);
> 
> dhatwxu:=0.5*(dhatwr[3*ii+3*n1i+4] +dhatwr[3*ii+1] -dhatwr[3*ii+3*n1i+1] -dhatwr[3*ii-2])/wt1;
> dhatwyu:=0.5*(dhatwr[3*ii+3*n1i+5] +dhatwr[3*ii+2] -dhatwr[3*ii+3*n1i+2] -dhatwr[3*ii-1])/wt1;
> dhatwzu:=0.5*(dhatwr[3*ii+3*n1i+6] +dhatwr[3*ii+3] -dhatwr[3*ii+3*n1i+3] -dhatwr[3*ii])/wt1;
> dhatwxv:=0.5*(dhatwr[3*ii+3*n1i+4] -dhatwr[3*ii+1] +dhatwr[3*ii+3*n1i+1] -dhatwr[3*ii-2])/wt2;
> dhatwyv:=0.5*(dhatwr[3*ii+3*n1i+5] -dhatwr[3*ii+2] +dhatwr[3*ii+3*n1i+2] -dhatwr[3*ii-1])/wt2;
> dhatwzv:=0.5*(dhatwr[3*ii+3*n1i+6] -dhatwr[3*ii+3] +dhatwr[3*ii+3*n1i+3] -dhatwr[3*ii])/wt2;
>      
> dhatcuu:=(psi*dhatpsi-2*dpsi*hatpsi)/psi^2*cuu
> -psi*(dhatgx_oben_u*wxu+dhatgy_oben_u*wyu+dhatgz_oben_u*wzu+gx_oben_u*dhatwxu+gy_oben_u*dhatwyu+gz_oben_u*dhatwzu
> +hatgx_oben_u*dwxu+hatgy_oben_u*dwyu+hatgz_oben_u*dwzu+dgx_oben_u*hatwxu+dgy_oben_u*hatwyu+dgz_oben_u*hatwzu)
> +(dpsi*hatcuu+hatpsi*dcuu)/psi;
> dhatcuv:=(psi*dhatpsi-2*dpsi*hatpsi)/psi^2*cuv
> -psi*(dhatgx_oben_u*wxv+dhatgy_oben_u*wyv+dhatgz_oben_u*wzv+gx_oben_u*dhatwxv+gy_oben_u*dhatwyv+gz_oben_u*dhatwzv
> +hatgx_oben_u*dwxv+hatgy_oben_u*dwyv+hatgz_oben_u*dwzv+dgx_oben_u*hatwxv+dgy_oben_u*hatwyv+dgz_oben_u*hatwzv)
> +(dpsi*hatcuv+hatpsi*dcuv)/psi;
> dhatcvu:=(psi*dhatpsi-2*dpsi*hatpsi)/psi^2*cvu
> -psi*(dhatgx_oben_v*wxu+dhatgy_oben_v*wyu+dhatgz_oben_v*wzu+gx_oben_v*dhatwxu+gy_oben_v*dhatwyu+gz_oben_v*dhatwzu
> +hatgx_oben_v*dwxu+hatgy_oben_v*dwyu+hatgz_oben_v*dwzu+dgx_oben_v*hatwxu+dgy_oben_v*hatwyu+dgz_oben_v*hatwzu)
> +(dpsi*hatcvu+hatpsi*dcvu)/psi;
> dhatcvv:=(psi*dhatpsi-2*dpsi*hatpsi)/psi^2*cvv
> -psi*(dhatgx_oben_v*wxv+dhatgy_oben_v*wyv+dhatgz_oben_v*wzv+gx_oben_v*dhatwxv+gy_oben_v*dhatwyv+gz_oben_v*dhatwzv
> +hatgx_oben_v*dwxv+hatgy_oben_v*dwyv+hatgz_oben_v*dwzv+dgx_oben_v*hatwxv+dgy_oben_v*hatwyv+dgz_oben_v*hatwzv)
> +(dpsi*hatcvv+hatpsi*dcvv)/psi;
> dhatinv1:=dhatcuu+dhatcvv;
> dhatinv2:=hatcuu*dcvv+dcuu*hatcvv-hatcuv*dcvu-dcuv*hatcvu+dhatcuu*cvv+cuu*dhatcvv-dhatcuv*cvu-cuv*dhatcvu;
> dhatphi:=(phitrtr(inv1,inv2)*dinv1+phitrdt(inv1,inv2)*dinv2)*hatinv1
>          +(phitrdt(inv1,inv2)*dinv1+phidtdt(inv1,inv2)*dinv2)*hatinv2
>          +phitr(inv1,inv2)*dhatinv1+phidt(inv1,inv2)*dhatinv2;
> 
> dhataccontr:= (2*dpsi*hatpsi-psi*dhatpsi)/psi^3*wt1*wt2;
> dhatac:=dhatac+dhataccontr;
> hataccontr:=-hatpsi/psi^2*wt1*wt2;
> dhatwc:=dhatwc+phi(inv1,inv2)*dhataccontr+dhatphi*accontr+hatphi*daccontr+dphi*hataccontr;
> dhata1c:=dhata1c+inv1*dhataccontr+dhatinv1*accontr+hatinv1*daccontr+dinv1*hataccontr;
> dhata2c:=dhata2c+inv2*dhataccontr+dhatinv2*accontr+hatinv2*daccontr+dinv2*hataccontr;
> 
> gxutilde:=((xn0-x00)*(1-v/m)+(xnm-x0m)*v/m)/n  +zeta*mmxu;
> gyutilde:= ((yn0-y00)*(1-v/m)+(ynm-y0m)*v/m)/n +zeta*mmyu;
> gzutilde:=zeta*mmzu;
> gxvtilde:=((x0m-x00)*(1-u/n)+(xnm-xn0)*u/n)/m+zeta*mmxv;
> gyvtilde:=((y0m-y00)*(1-u/n)+(ynm-yn0)*u/n)/m+zeta*mmyv;
> gzvtilde:=zeta*mmzv;
> dhathiwx:=gyutilde*mmzv-gzutilde*mmyv+mmyu*gzvtilde-mmzu*gyvtilde;
> dhathiwy:=gzutilde*mmxv-gxutilde*mmzv+mmzu*gxvtilde-mmxu*gzvtilde;
> dhathiwz:=gxutilde*mmyv-gyutilde*mmxv+mmxu*gyvtilde-mmyu*gxvtilde;
> dhatvolu:=dhatvolu+(dhathiwx*mmx+dhathiwy*mmy+dhathiwz*mmz)*dzeta*hatzeta*wt1*wt2;
> 
>  hi1:=phitrtr(inv1,inv2)*inv1+phitrdt(inv1,inv2)*2*inv2+phitr(inv1,inv2);
>  hi2:=phitrdt(inv1,inv2)*inv1+phidtdt(inv1,inv2)*2*inv2+2*phidt(inv1,inv2);
> a_dwc:=a_dwc+2/a*dwccontr-1/a*(hi5*daccontr+(hi1*dinv1+hi2*dinv2)*accontr);
> a_dac:=a_dac+2/a*daccontr;
> a_dvoc:=a_dvoc+(-3*dvolucontr-hh*daccontr+hh^2/8*da1ccontr)/a;
> a_dvwc:=a_dvwc+(3*dvolucontr-hh*daccontr-hh^2/8*da1ccontr)/a;
> 
> hi6:=phitr(inv1,inv2)*inv1-phitrtr(inv1,inv2)*inv1^2-4*phitrdt(inv1,inv2)*inv1*inv2
>      -4*phidtdt(inv1,inv2)*inv2^2;
> atimesa_intc:=atimesa_intc+hi6*accontr;
> od od:
> voc:=vop0-volu-hh/2*ac+hh^2/8*a1c-hh^3/24*a2c:
> vwc:=vwp0+volu-hh/2*ac-hh^2/8*a1c-hh^3/24*a2c:
> dvoc:=-dvolu-hh/2*dac+hh^2/8*da1c-hh^3/24*da2c:
> dvwc:= dvolu-hh/2*dac-hh^2/8*da1c-hh^3/24*da2c:
> dhatvoc:=-dhatvolu-hh/2*dhatac+hh^2/8*dhata1c-hh^3/24*dhata2c:
> dhatvwc:= dhatvolu-hh/2*dhatac-hh^2/8*dhata1c-hh^3/24*dhata2c:

> denergy:=dwc-f*kf*dac-po*ko*dvoc-pw*kw*dvwc:
> dhatenergy:=dhatwc-f*kf*dhatac-po*ko*dhatvoc-pw*kw*dhatvwc:
> a_denergy:=a_dwc-f*kf*a_dac-po*ko*a_dvoc-pw*kw*a_dvwc:

> rs:=vector(nn+5):
> mit:=matrix(nn+5,nn+5,0):
> for jj to n2i do
> for ij to n1i do
>   kc:=(jj-1)*n1i+ij:  
>   rs[kc]:=-ns*coeff(denergy,dtheta[kc]);
> 
>   mit[kc,nn+1]:=kf*coeff(denergy,dtheta[kc]);
>   mit[nn+1,kc]:=kf*coeff(denergy,dtheta[kc]);
>   mit[kc,nn+2]:=-kf*ns*coeff(dac,dtheta[kc]);
>   mit[nn+2,kc]:=-kf*ns*coeff(dac,dtheta[kc]);
>   mit[kc,nn+3]:=-ko*ns*coeff(dvoc,dtheta[kc]);
>   mit[nn+3,kc]:=-ko*ns*coeff(dvoc,dtheta[kc]);
>   mit[kc,nn+4]:=-kw*ns*coeff(dvwc,dtheta[kc]);
>   mit[nn+4,kc]:=-kw*ns*coeff(dvwc,dtheta[kc]);
>   mit[kc,nn+5]:=ka*ns*coeff(a_denergy,dtheta[kc]);
>   mit[nn+5,kc]:=ka*ns*coeff(a_denergy,dtheta[kc]);
>   rowkc:=ns*coeff(dhatenergy,dtheta[kc]);
> 
> for mm to n2i do
> for ll to n1i do
>   kb:=(mm-1)*n1i+ll:
>   mit[kc,kb]:=coeff(rowkc,hattheta[kb])  
> od
> od
> od;
> print(jj)
> od:
> 
# Supplementing the last five rows and colums of the matrix:
> 
> rs[nn+1]:=-(wc-f*kf*ac-po*ko*voc-pw*kw*vwc)*kf:
> rs[nn+2]:=-(1-ns*ac)*kf:
> bed_o:=L_o-ns*voc:
> bed_w:=L_w-ns*vwc:
> rs[nn+3]:=-bed_o*ko:
> rs[nn+4]:=-bed_w*kw:
> rhi9:=kf*2/a*ac:
> rhi10:=ko/a*(3*(vop0-volu)-hh*ac+hh^2/8*a1c):
> rhi11:=kw/a*(3*(vwp0+volu)-hh*ac-hh^2/8*a1c):
> rhi8:=(2*wc-intc)/a-f*rhi9-po*rhi10-pw*rhi11:
> rs[nn+5]:=-ka*ns*rhi8:
> 
> mit[nn+1,nn+3]:=-kf*ko*voc:
> mit[nn+3,nn+1]:=-kf*ko*voc:
> mit[nn+1,nn+4]:=-kf*kw*vwc:
> mit[nn+4,nn+1]:=-kf*kw*vwc:
> mit[nn+1,nn+2]:=-kf*ac:
> mit[nn+2,nn+1]:=-kf*ac:
> mit[nn+1,nn+1]:=1-kf:
> mit[nn+2,nn+2]:=1-kf:
> mit[nn+3,nn+3]:=1-ko:
> mit[nn+4,nn+4]:=1-kw:
> mit[nn+1,nn+5]:=ka*rhi8:
> mit[nn+5,nn+1]:=ka*rhi8:
> mit[nn+2,nn+5]:=-ka*ns*rhi9:
> mit[nn+5,nn+2]:=-ka*ns*rhi9:
> mit[nn+3,nn+5]:=-ka*ns*rhi10:
> mit[nn+5,nn+3]:=-ka*ns*rhi10:
> mit[nn+4,nn+5]:=-ka*ns*rhi11:
> mit[nn+5,nn+4]:=-ka*ns*rhi11:
> mit[nn+5,nn+5]:=ka*ns/a^2*( 2*wc-intc-atimesa_intc
> -2*f*kf*ac-po*ko*(6*(vop0-volu)-hh*ac)-pw*kw*(6*(vwp0+volu)-hh*ac)) +1-ka:
> 
# 
# Residuum vector:
> 
> evalm(rs);
> 
> 
# 
# Solution of the Newton-Raphson equations
> 
> ds:=linsolve(mit,rs):
> evalm(ds);
> 
> 
# Discussion of the results of the last iteration step:
> ser:=seq(abs(rs[i]),i=1..nn):
> ses:=seq(abs(ds[i]),i=1..nn):

> mar:=[fakt,max(ser),seq(rs[i],i=nn+1..nn+5)],mar:
> convert([mar],matrix);
> mas:=[fakt,max(ses),seq(ds[i],i=nn+1..nn+5)],mas:
> convert([mas],matrix);
> 
> mean1:=zi1/ac;
> dev1:=sqrt(max(0.0,ziq1/ac-mean1^2));
> mean2:=zi2/ac;
> dev2:=sqrt(max(0.0,ziq2/ac-mean2^2));
> mean3:=zi3/ac;
> dev3:=sqrt(max(0.0,ziq3/ac-mean3^2));
> mean4:=zi4/ac;
> dev4:=sqrt(max(0.0,ziq4/ac-mean4^2));
> mean5:=zi5/ac;
> dev5:=sqrt(max(0.0,ziq5/ac-mean5^2));
> mean6:=zi6/ac;
> dev6:=sqrt(max(0.0,ziq6/ac-mean6^2));
> mean7:=zi7/ac;
> dev7:=sqrt(max(0.0,ziq7/ac-mean7^2));
> mean8:=zi8/ac;
> dev8:=sqrt(max(0.0,ziq8/ac-mean8^2));
> mean9:=zi9/ac;
> dev9:=sqrt(max(0.0,ziq9/ac-mean9^2));

# 
> wcalt:=wc:
> acalt:=ac:
> vocalt:=voc:
> vwcalt:=vwc:

> fakt:=1.0;
> #fakt:=0.3; # reduce the step by an appropriate factor if there are problems with the convergence
> 
> theta:=evalm(theta+fakt*ds):
> 
# If , in your opinion, the results are not accurate enough, go back to the label and start the next iteration !!!! -> -> -> -> 
# Otherwise proceed and save the results 
> 
# Note: mean3, dev3, mean4, dev4, outp3, and outp4 are only correct if at least two iterations have been performed!

> #########################################################################
> tnum:=63180020; # change the identification number if necessary
# Graphical representation:
> 
> outp1||tnum:=[seq([out1[j]],j=1..n2i)]:
> outp2||tnum:=[seq([out2[j]],j=1..n2i)]:
> outp3||tnum:=[seq([out3[j]],j=1..n2i)]:
> outp4||tnum:=[seq([out4[j]],j=1..n2i)]:
> outp5||tnum:=[seq([out5[j]],j=1..n2i)]:
> outp6||tnum:=[seq([out6[j]],j=1..n2i)]:
> outp7||tnum:=[seq([out7[j]],j=1..n2i)]:
> outp8||tnum:=[seq([out8[j]],j=1..n2i)]:
> outp9||tnum:=[seq([out9[j]],j=1..n2i)]:

> plot(outp1||tnum);
> plot(outp2||tnum);
> plot(outp3||tnum);
> plot(outp4||tnum);
> plot(outp5||tnum);
> plot(outp6||tnum);
> plot(outp7||tnum);
> plot(outp8||tnum);
> plot(outp9||tnum);
# Storage of the final results:
> 
> daten||tnum:=[tnum,hh,n1i,n2i,L_o,nu,D0,H0,a*theta[1],ns*wc,a,po,ac,voc,vwc,wc,ns,f,pw,bed_o,bed_w,
> ns*voc,ns*vwc,mean1,dev1,mean2,dev2,mean3,dev3,mean4,dev4,mean5,dev5,mean6,dev6,mean7,dev7,mean8,dev8,mean9,dev9,a1c,a2c,
> evalf(-6/Pi*a2c)];
> thetalist||tnum:=convert(theta,list);
> 
> save daten||tnum,thetalist||tnum,diffdata||tnum:
> save outp1||tnum,outp2||tnum,outp3||tnum,outp4||tnum,outp5||tnum,
> outp6||tnum,outp7||tnum,outp8||tnum,outp9||tnum,diffplots||tnum:
# 
> 
#  Spatial representation of the structure:
> 
> 
> form:=proc(u,v)
> global theta;
> local kn00,kn10,kn01,kn11,xi,eta,flu,flv;
>   flu:=floor(0.9999*u);
>   flv:=floor(0.9999*v);
>   kn00:=flu+1+flv*(ni+1);
>   kn10:=kn00+1;
>   kn01:=kn00+ni+1;
>   kn11:=kn10+ni+1;
> subs(xi=u-flu,eta=v-flv,
>       a*(theta[kn00]*(1-xi)*(1-eta)+theta[kn10]*xi*(1-eta)
>           +theta[kn01]*(1-xi)*eta+theta[kn11]*xi*eta))          
>   end:

> u:='u': v:='v': 
>        xxp:=x00*(1-u/n)*(1-v/m)+xn0*u/n*(1-v/m)
>           +x0m*(1-u/n)*v/m+xnm*u/n*v/m + form(u,v)*mmx:
>        
> yyp:=y00*(1-u/n)*(1-v/m)+yn0*u/n*(1-v/m)
>           +y0m*(1-u/n)*v/m+ynm*u/n*v/m + form(u,v)*mmy:
>       
>  zzp:=zz0+form(u,v)*mmz:
> 
> if sphere or cylinder or tube6 then
>  spc4:=spacecurve([[0,0,0],[a,0,0],[a,a,0],[0,a,0],[0,0,0]],thickness=3)
> fi:
> 
> if cylinder then
>  spc5:=spacecurve({[[0,2*a,-a],[0,0,-a],[2*a,0,-a],[0,0,a],[0,2*a,a],[2*a,2*a,-a],[2*a,0,-a],[0,0,-a],
> [0,0,a],[0,2*a,a],[0,2*a,-a],[2*a,2*a,-a]]},linestyle=3)
> fi:

> if sphere then 
>  spc5:=spacecurve({[[0,2*a,-a],[0,0,-a],[2*a,0,-a],[0,0,a],[2*a,2*a,-a],[0,2*a,-a],[0,0,a],
>  [2*a,0,-a],[2*a,2*a,-a],[0,0,a],[0,0,-a]]},linestyle=3)
> fi:

> if tube6 then
>  spc5:=spacecurve({[[0,0,a],[0,0,-a],[2*a,0,-a],[0,0,a]],
>  [[0,0,-a],[0,2*a,a],[0,0,a]],
>  [[0,0,a],[2*a,0,-a],[0,2*a,a],[0,a,a]],
>  [[0,0,-a],[0,2*a,a],[2*a,0,-a],[a,0,-a]]},linestyle=3)
> fi:
> 
> if tube4 then
> spc4:=spacecurve([[-a/sqrt(18.),-2*a/3.,-a/sqrt(6.)],[a/sqrt(18.),-5*a/6.,-a/sqrt(6.)],[a/sqrt(18.),-a/3.,-a/sqrt(6.)],
> [-a/sqrt(18.),-a/6.,-a/sqrt(6.)],[-a/sqrt(18.),-2*a/3.,-a/sqrt(6.)]],thickness=3):
> spc5:=spacecurve({[[0,-a,0],[0,0,0],[-sqrt(2.)/3.*a,-a/3.,-sqrt(2./3.)*a],[0,-a,0],[sqrt(2.)/3.*a,-2*a/3.,-sqrt(2./3.)*a],
> [0,0,0],[sqrt(2.)/3.*a,-2*a/3.,-sqrt(2./3.)*a],[-sqrt(2.)/3.*a,-a/3.,-sqrt(2./3.)*a]]},linestyle=3)
> fi:

> if tube3 then
> spc4:=spacecurve([[-a/2/sqrt(6.),-a/2.,-a/2/sqrt(6.)],[a/2/sqrt(6.),-a/2.,-a/2/sqrt(6.)],[a/2/sqrt(6.),0,-a/2/sqrt(6.)],[-a/2/sqrt(6.),-a/6.,-a/2/sqrt(6.)],[-a/2/sqrt(6.),-a/2.,-a/2/sqrt(6.)]],thickness=3):
> spc5:=spacecurve({[[-a/sqrt(6.)*1.5,-a/2.,-a/sqrt(6.)*1.5],[0,-a/2.,0],[a/sqrt(6.)*1.5,-a/2.,-a/sqrt(6.)*1.5]],[[0,-a/2.,0],[0,0,0],   [a/sqrt(6.)*1.5,0,-a/sqrt(6.)*1.5]],[[0,0,0],[-a/sqrt(6.)*1.5,-a/2,-a/sqrt(6.)*1.5]]},linestyle=3):
> #spc5:=spacecurve({[[-a/sqrt(6.),-a/3.,-a/sqrt(6.)],[-a/sqrt(6.),-a/2.,-a/sqrt(6.)],[a/sqrt(6.),-a/2.,-a/sqrt(6.)],[0,-a/2.,0],[-a/sqrt(6.),-a/2.,-#a/sqrt(6.)]],[[0,-a/2.,0],[0,0,0],   [a/sqrt(6.),0,-a/sqrt(6.)],[a/sqrt(6.),-a/2.,-a/sqrt(6.)] ],
> # [[0,0,0],[-a/sqrt(6.),-a/3.,-a/sqrt(6.)], [a/sqrt(6.),0,-a/sqrt(6.)]]},linestyle=3)
> fi:
> 
> ku:=plot3d([xxp,yyp,zzp],u=0..n,v=0..m):
> 
> 
> display({spc5,spc4,ku},axes=framed,color=black,style=hidden,scaling=constrained,
> orientation=[-145,60],labels=["x","y","z"],tickmarks=[5,5,5]);
> 
>