Select Git revision
mainLevelset.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
OrthogonalPoly.cpp 10.39 KiB
#include "OrthogonalPoly.h"
double OrthogonalPoly::EvalLobatto(int order, double x)
{
double L=0;
double xsquare=pow(x,2);
switch(order) {
case(0):
L =0.5*(1-1*x);
return L;
case(1):
L =0.5*(1+x);
return L;
case(2):
L=(-1+pow(x,2));
L=L*0.5*pow(3./2.,0.5);
return L;
case(3):
L=x*(-1+xsquare);
L=L*0.5*pow(5./2.,0.5);
return L;
case(4):
L=1+xsquare*(-6+5.*xsquare);
L=L*1./8.*pow(7./2.,0.5);
return L;
case(5):
L=x*(3+xsquare*(-10+7.*xsquare));
L=L*3./8.*pow(2.,-0.5);
return L;
case(6):
L=-1+xsquare*(15+xsquare*(-35+21*xsquare));
L=L*1./16.*pow(11./2.,0.5);
return L;
case(7):
L=x*(-5+xsquare*(35+xsquare*(-63+33*xsquare)));
L=L*1./16.*pow(13./2.,0.5);
return L;
case(8):
L=5+xsquare*(-140+xsquare*(630+xsquare*(-924+429*xsquare)));
L=L*1./128.*pow(15./2.,0.5);
return L;
case(9):
L=x*(35+xsquare*(-420+xsquare*(1386+xsquare*(-1716+715*xsquare))));
L=L*1./128.*pow(17./2.,0.5);
return L;
case(10):
L=-7+xsquare*(315+xsquare*(-2310+xsquare*(6006+xsquare*(-6435+2431*xsquare))));
L=L*1./256.*pow(19./2.,0.5);
return L;
case(11):
L=x*(-63+xsquare*(1155+xsquare*(-6006+xsquare*(12870+xsquare*(-12155+4199*xsquare)))));
L=L*1./256.*pow(21./2.,0.5);
return L;
case(12):
L=21+xsquare*(-1386+xsquare*(15015+xsquare*(-60060+xsquare*(109395+xsquare*(-92378+29393*xsquare)))));
L=L*1./1024.*pow(23./2.,0.5);
return L;
case(13):
L=x*(231+xsquare*(-6006+xsquare*(45045+xsquare*(-145860+xsquare*(230945+xsquare*(-176358+52003*xsquare))))));
L=L*5./1024.*pow(2.,-0.5);
return L;
case(14):
L=-33+xsquare*(3003+xsquare*(-45045+xsquare*(255255+xsquare*(-692835+xsquare*(969969+xsquare*(-676039+185725*xsquare))))));
L=L*3./2048.*pow(3./2.,0.5);
return L;
case(15):
L=x*(-429 +xsquare*(15015+xsquare*(-153153+xsquare*(692835+xsquare*(-1616615+xsquare*(2028117+xsquare*(-1300075+334305*xsquare)))))));
L=L*1./2048.*pow(29./2.,0.5); return L;
default:
throw std::string("Lobatto functions are written for orders =< 15");
}
}
double OrthogonalPoly::EvalDLobatto(int order, double x){
double dL=0;
double xsquare=pow(x,2);
switch(order) {
case(0):
dL =-0.5;
return dL;
case(1):
dL =0.5;
return dL;
case(2):
dL=2*x;
dL=dL*0.5*pow(3./2.,0.5);
return dL;
case(3):
dL=-1+3*xsquare;
dL=dL*0.5*pow(5./2.,0.5);
return dL;
case(4):
dL=x*(-12+20*xsquare);
dL=dL*1./8.*pow(7./2.,0.5);
return dL;
case(5):
dL=3+xsquare*(-30+35.*xsquare);
dL=dL*3./8.*pow(2.,-0.5);
return dL;
case(6):
dL=x*(30+xsquare*(-140+126*xsquare));
dL=dL*1./16.*pow(11./2.,0.5);
return dL;
case(7):
dL=-5+xsquare*(105+xsquare*(-315+231*xsquare));
dL=dL*1./16.*pow(13./2.,0.5);
return dL;
case(8):
dL=x*(-280+xsquare*(2520+xsquare*(-5544+3432*xsquare)));
dL=dL*1./128.*pow(15./2.,0.5);
return dL;
case(9):
dL=35+xsquare*(-1260+xsquare*(6930+xsquare*(-12012+6435*xsquare)));
dL=dL*1./128.*pow(17./2.,0.5);
return dL;
case(10):
dL=x*(630+xsquare*(-9240+xsquare*(36036+xsquare*(-51480+24310 * xsquare))));
dL=dL*1./256.*pow(19./2.,0.5);
return dL;
case(11):
dL=-63+xsquare*(3465+xsquare*(-30030+xsquare*(90090+xsquare*(-109395+46189*xsquare))));
dL=dL*1./256.*pow(21./2.,0.5);
return dL;
case(12):
dL=x*(-2772+xsquare*(60060+xsquare*(-360360+xsquare*(875160+xsquare*(-923780+352716*xsquare)))));
dL=dL*1./1024.*pow(23./2.,0.5);
return dL;
case(13):
dL=231+xsquare*(-18018+xsquare*(225225+xsquare*(-1021020+xsquare*(2078505+xsquare*(-1939938+676039*xsquare)))));
dL=dL*5./1024.*pow(2.,-0.5);
return dL;
case(14):
dL=x*(6006+xsquare*(-180180+xsquare*(1531530+xsquare*(-5542680+xsquare*(9699690+xsquare*(-8112468+2600150*xsquare)))))); dL=dL*3./2048.*pow(3./2.,0.5);
return dL;
case(15):
dL=-429 +xsquare*(45045+xsquare*(-765765+xsquare*(4849845+xsquare*(-14549535+xsquare*(22309287+xsquare*(-16900975+5014575*xsquare))))));
dL=dL*1./2048.*pow(29./2.,0.5);
return dL;
default:
throw std::string("Lobatto functions are written for orders =< 15");
}
}
double OrthogonalPoly::EvalKernelFunction(int order, double x){
double phi=0;
double xsquare=pow(x,2);
switch(order){
case(0):
phi=-pow(6,0.5);
return phi;
case(1):
phi=-x*pow(10,0.5);
return phi;
case(2):
phi=1-5*xsquare;
phi=phi*0.5*pow(7./2.,0.5);
return phi;
case(3):
phi=x*(3-7*xsquare);
phi=phi*3./2.*pow(2,-0.5);
return phi;
case(4):
phi=-1+xsquare*(14-21*xsquare);
phi=phi*1./4.*pow(11./2.,0.5);
return phi;
case(5):
phi=x*(-5+xsquare*(30-33*xsquare));
phi=phi*1./4.*pow(13./2.,0.5);
return phi;
case(6):
phi=5+xsquare*(-135+xsquare*(495-429*xsquare));
phi=phi*1./32.*pow(15./2.,0.5);
return phi;
case(7):
phi=x*(35+xsquare*(-385+xsquare*(1001-715*xsquare)));
phi=phi*1./32.*pow(17./2.,0.5);
return phi;
case(8):
phi=-7+xsquare*(308+xsquare*(-2002+xsquare*(4004-2431*xsquare)));
phi=phi*1./64.*pow(19./2.,0.5);
return phi;
case(9):
phi=x*(-63+xsquare*(1092+xsquare*(-4914+xsquare*(7956-4199*xsquare))));
phi=phi*1./64.*pow(21./2.,0.5);
return phi;
case(10):
phi=21+xsquare*(-1365+xsquare*(13650+xsquare*(-46410+xsquare*(62985 -29393*xsquare))));
phi=phi*1./256.*pow(23./2.,0.5);
return phi;
case(11):
phi=x*(231+xsquare*(-5775+xsquare*(39270+xsquare*(-106590+xsquare*(124355-52003*xsquare)))));
phi=phi*5./256.*pow(2,-0.5);
return phi;
case(12):
phi=-33+xsquare*(2970+xsquare*(-42075+xsquare*(213180+xsquare*(-479655+xsquare*(490314-185725*xsquare)))));
phi=phi*3./512.*pow(3./2.,0.5); return phi;
case(13):
phi=x*(-429+xsquare*(14586+xsquare*(-138567+xsquare*(554268+xsquare*(-1062347+xsquare*(965770-334305*xsquare))))));
phi=phi*1./512.*pow(29./2.,0.5);
return phi;
default:
throw std::string("Lobatto functions are written for orders =< 15");
}
}
double OrthogonalPoly::EvalDKernelFunction(int order, double x){
double dphi=0;
double xsquare=pow(x,2);
switch(order){
case(0):
dphi=0;
return dphi;
case(1):
dphi=-pow(10,0.5);
return dphi;
case(2):
dphi=-10*x;;
dphi=dphi*0.5*pow(7./2.,0.5);
return dphi;
case(3):
dphi=3-21*xsquare;
dphi=dphi*3./2.*pow(2,-0.5);
return dphi;
case(4):
dphi=x*(28-84*xsquare);
dphi=dphi*1./4.*pow(11./2.,0.5);
return dphi;
case(5):
dphi=-5+xsquare*(90-165*xsquare);
dphi=dphi*1./4.*pow(13./2.,0.5);
return dphi;
case(6):
dphi= x*(-270+xsquare*(1980-2574*xsquare));
dphi=dphi*1./32.*pow(15./2.,0.5);
return dphi;
case(7):
dphi=35+xsquare*(-1155+xsquare*(5005-5005*xsquare));
dphi=dphi*1./32.*pow(17./2.,0.5);
return dphi;
case(8):
dphi=x*(616+xsquare*(-8008+xsquare*(24024-19448*xsquare)));
dphi=dphi*1./64.*pow(19./2.,0.5);
return dphi;
case(9):
dphi=-63+xsquare*(3276+xsquare*(-24570+xsquare*(55692-37791*xsquare)));
dphi=dphi*1./64.*pow(21./2.,0.5);
return dphi;
case(10):
dphi=x*(-2730+xsquare*(54600+xsquare*(-278460+xsquare*(503880-293930*xsquare))));
dphi=dphi*1./256.*pow(23./2.,0.5);
return dphi;
case(11):
dphi=231+xsquare*(-17325+xsquare*(196350+xsquare*(-746130+xsquare*(1119195-572033*xsquare))));
dphi=dphi*5./256.*pow(2,-0.5);
return dphi;
case(12):
dphi= x*(5940+xsquare*(-168300+xsquare*(1279080+xsquare*(-3837240+xsquare*(4903140-2228700*xsquare)))));
dphi=dphi*3./512.*pow(3./2.,0.5);
return dphi;
case(13):
dphi = -429+xsquare*(43758+xsquare*(-692835+xsquare*(3879876+xsquare*(-9561123+xsquare*(10623470-4345965*xsquare))))); dphi=dphi*1./512.*pow(29./2.,0.5);
return dphi;
default:
throw std::string("Lobatto functions are written for orders =< 15");
}
}
double OrthogonalPoly::EvalLegendre(int order, double x){
double L=0;
double xsquare=pow(x,2);
switch(order) {
case(0):
L =1;
return L;
case(1):
L =x;
return L;
case(2):
L=3./2.*xsquare-1./2.;
return L;
case(3):
L=0.5*x*(5*xsquare-3);
return L;
case(4):
L=(3+xsquare*(35*xsquare-30));
L=1./8.*L;
return L;
case(5):
L=x*(xsquare*(63*xsquare-70)+15);
L=1./8.*L;
return L;
case(6):
L=((231*xsquare-315)*xsquare+105)*xsquare-5;
L=1./16.*L;
return L;
case(7):
L=x*(((429*xsquare-693)*xsquare+315)*xsquare-35);
L=1./16.*L;
return L;
case(8):
L=(((6435*xsquare-12012)*xsquare+6930)*xsquare-1260)*xsquare+35;
L=1./128.*L;
return L;
case(9):
L=((((12155*xsquare-25740)*xsquare+18018)*xsquare-4620)*xsquare+315)*x;
L=1./128.*L;
return L;
case(10):
L=((((46189*xsquare-109395)*xsquare+90090)*xsquare-30030)*xsquare+3465)*xsquare- 63;
L=1./256.*L;
return L;
default:
throw std::string("Legendre functions are written for orders =< 10");
}
}
double OrthogonalPoly::EvalDLegendre(int order, double x){
double dL=0;
double xsquare=pow(x,2);
switch(order) {
case(0):
dL =0;
return dL;
case(1):
dL =1; return dL;
case(2):
dL=3*x;
return dL;
case(3):
dL=0.5*(15*xsquare-3);
return dL;
case(4):
dL=x*(140*xsquare-60);
dL=1./8.*dL;
return dL;
case(5):
dL=15+xsquare*(315*xsquare-210);
dL=1./8.*dL;
return dL;
case(6):
dL=x*(210+xsquare*(1386*xsquare-1260));
dL=1./16.*dL;
return dL;
case(7):
dL=((xsquare*3003-3465)*xsquare+945)*xsquare-35;
dL=1./16.*dL;
return dL;
case(8):
dL=x*(((51480*xsquare-72072)*xsquare+27720)*xsquare-2520);
dL=1./128.*dL;
return dL;
case(9):
dL=315+xsquare*(-13860+xsquare*(90090+xsquare*(-180180+109395*xsquare)));
dL=1./128.*dL;
return dL;
case(10):
dL=x*(6930+xsquare*(-120120+xsquare*(540540+xsquare*(-875160+461890*xsquare))));
dL=1./256.*dL;
return dL;
default:
throw std::string("Legendre functions are written for orders =< 10");
}
}