00001
00005 #ifndef CIRCLE_PATTERN_CM
00006 #define CIRCLE_PATTERN_CM
00007 #define NOMINMAX
00008 #include <iostream>
00009 #define _USE_MATH_DEFINES
00010 #include <math.h>
00011
00012 class Mesh;
00013
00043 struct CirclePattern {
00044
00045 static double evaluateEnergy (Mesh * mesh, double * xx);
00046
00047 static void computeGradient (Mesh * mesh, double * xx, double * grad);
00048
00049 static void buildSparsityForHessian (Mesh * mesh, int * hessI, int * hessJ);
00050
00051 static void computeHessian (Mesh * mesh, double * xx, double * hessval);
00052
00053
00054 static double ImLi2Sum (const double& theta, const double& x,
00055 const double& tanHalfThetaStar, const double& clThetaStar);
00056
00058 static double Cl2 (const double & _x);
00059
00062 static double fTheta2(const double& cosTheta, const double& sinTheta, const double& x);
00063
00065 static double fThetaPrime2(const double & cosTheta, const double& sinTheta, const double & x) {
00066 double res = sinTheta/(cosh(x) - cosTheta);
00067 return res;
00068 }
00069
00071 static double IEEEremainder(const double & x, const double & div);
00072 };
00073
00074 inline double CirclePattern::fTheta2(const double& cosTheta, const double& sinTheta, const double& x){
00075 double ex = exp(x);
00076 double s1 = 1.0 - ex*cosTheta;
00077 double s2 = ex*sinTheta;
00078
00079 double res = M_PI - 2*atan2(s1, s2);
00080
00081 return res;
00082 }
00083
00084 inline double CirclePattern::IEEEremainder(const double & x, const double & div) {
00085
00086 double res;
00087 if (x >= 0) {
00088 if (x < div) res = x;
00089 else {
00090 double part = floor(x / div);
00091 res = x - part*div;
00092 }
00093 }
00094 else {
00095 double part = floor(x / div);
00096 res = x - part*div;
00097 }
00098 if (res < 0 || res > div) {
00099 std::cout << "x " << x << " div " << div << std::endl;
00100 std::cerr << "Wrong IEEEreminder!!!" << std::endl;
00101 }
00102 return res;
00103 }
00104
00105 inline double CirclePattern::Cl2 (const double & _x) {
00106 double TWO_PI = 2 * M_PI;
00107 double MINUS_LOG_2 = log(0.5);
00108 double BORDERLINE = 2.0944;
00109
00110 if (_x == 0.0) {
00111 return 0.0;
00112 }
00113
00114 double x = IEEEremainder(_x, TWO_PI);
00115
00116 if (x == 0.0) {
00117 return 0.0;
00118 }
00119
00120 if (fabs(x) <= BORDERLINE) {
00121 double xx = x * x;
00122 return ((((((((((((
00123 2.3257441143020875e-22 * xx
00124 + 1.0887357368300848e-20) * xx
00125 + 5.178258806090624e-19) * xx
00126 + 2.5105444608999545e-17) * xx
00127 + 1.2462059912950672e-15) * xx
00128 + 6.372636443183181e-14) * xx
00129 + 3.387301370953521e-12) * xx
00130 + 1.8978869988971e-10) * xx
00131 + 1.1482216343327455e-8) * xx
00132 + 7.873519778281683e-7) * xx
00133 + 0.00006944444444444444) * xx
00134 + 0.013888888888888888) * xx
00135 - log(fabs(x)) + 1.0) * x;
00136 }
00137 x += ((x > 0.0) ? - M_PI : M_PI);
00138 double xx = x * x;
00139 return ((((((((((((
00140 3.901950904063069e-15 * xx
00141 + 4.566487567193635e-14) * xx
00142 + 5.429792727596476e-13) * xx
00143 + 6.5812165661369675e-12) * xx
00144 + 8.167010963952222e-11) * xx
00145 + 1.0440290284867003e-9) * xx
00146 + 1.3870999114054669e-8) * xx
00147 + 1.941538399871733e-7) * xx
00148 + 2.927965167548501e-6) * xx
00149 + 0.0000496031746031746) * xx
00150 + 0.0010416666666666667) * xx
00151 + 0.041666666666666664) * xx
00152 + MINUS_LOG_2) * x;
00153 }
00154
00155 #endif
00156
