Geometry

00001 #ifndef NMWR_GB_ALGEBRAIC_PRIMITIVES_3D_H
00002 #define NMWR_GB_ALGEBRAIC_PRIMITIVES_3D_H
00003 
00004 
00005 // $LICENSE
00006 
00007 
00008 
00017 template<class POINT>
00018 struct dimension_dependent_primitives_3d 
00019   : public virtual basic_algebraic_primitives<POINT> {
00020 
00023   static scalar det3(const POINT& p1, const POINT& p2, const POINT& p3)
00024     {
00025       int l = pt::LowerIndex(p1);
00026       return(
00027              p1[l]*(p2[l+1]*p3[l+2] - p2[l+2]*p3[l+1])
00028              -p1[l+1]*(p2[l]*p3[l+2] - p2[l+2]*p3[l])
00029              +p1[l+2]*(p2[l]*p3[l+1] - p2[l+1]*p3[l])
00030              );
00031     }
00032 
00035   static scalar triangle_area(const POINT& p1, const POINT& p2, const POINT& p3)
00036     { return 0.5*norm_2(vectorproduct(p1-p2,p1-p3));}
00037 
00040   static POINT vectorproduct(const POINT& p, const POINT& q)   // p x q 
00041     { 
00042       int lp = LowerIndex(p)-1;
00043       int lq = LowerIndex(q)-1; 
00044       return (POINT(p[lp+2]*q[lq+3]-p[lp+3]*q[lq+2], 
00045                     p[lp+3]*q[lq+1]-p[lp+1]*q[lq+3], 
00046                     p[lp+1]*q[lq+2]-p[lp+2]*q[lq+1]));
00047     }
00048 
00054   static void inverse(POINT const& r1, POINT const& r2, POINT const& r3,
00055                       POINT      & i1, POINT      & i2, POINT      & i3)
00056     {
00057       // inversion by Cramer's rule 
00058       scalar d = det3(r1,r2,r3);
00059       REQUIRE_ALWAYS( (d != 0.0), "Matrix singular!\n",1);
00060       scalar d_inv = 1.0/d;
00061       int li = pt::LowerIndex(r1);
00062       POINT e[3]; // unit vectors
00063       for(int i = 0; i < 3; ++i) {
00064         e[i] = pt::Origin(3);
00065         e[i][i+ li] = 1;
00066       }
00067       i1[li+0] = d_inv * det3(e[0],r2,r3); 
00068       i2[li+0] = d_inv * det3(e[1],r2,r3); 
00069       i3[li+0] = d_inv * det3(e[2],r2,r3); 
00070 
00071       i1[li+1] = d_inv * det3(r1,e[0],r3); 
00072       i2[li+1] = d_inv * det3(r1,e[1],r3); 
00073       i3[li+1] = d_inv * det3(r1,e[2],r3); 
00074 
00075       i1[li+2] = d_inv * det3(r1,r2,e[0]); 
00076       i2[li+2] = d_inv * det3(r1,r2,e[1]); 
00077       i3[li+2] = d_inv * det3(r1,r2,e[2]); 
00078     }
00079 
00084   static void solve(POINT const& A1, POINT const& A2, POINT const& A3,
00085                     POINT      & x,  POINT const& b)
00086     {
00087       // solve by Cramer's rule
00088       scalar d = det3(A1,A2,A3);
00089       REQUIRE_ALWAYS( (d != 0.0), "Matrix singular!\n",1);
00090       scalar d_inv = 1.0/d;
00091       int li = pt::LowerIndex(A1);
00092       x[li+0] = d_inv * det3(b, A2,A3);
00093       x[li+1] = d_inv * det3(A1,b ,A3);
00094       x[li+2] = d_inv * det3(A1,A2,b );
00095     }
00096 
00097   static real matrixnorm_infinity(POINT const& A1, POINT const& A2, POINT const& A3)
00098     {
00099       real max_sum = 0;
00100       for(int i = pt::LowerIndex(A1); i <= pt::UpperIndex(A1); ++i) {
00101         real sum = fabs(A1[i]) + fabs(A2[i]) + fabs(A3[i]);
00102         max_sum = (max_sum < sum ? sum : max_sum);
00103       }
00104       return max_sum;
00105     }
00106 
00110   static real matrixnorm_1(POINT const& A1, POINT const& A2, POINT const& A3)
00111     {
00112       real max_sum = norm_1(A1);
00113       real sum     = norm_1(A2);
00114       max_sum = (max_sum < sum ? sum : max_sum);
00115       sum          = norm_1(A3);
00116       max_sum = (max_sum < sum ? sum : max_sum);
00117       return max_sum;
00118     }
00119 
00120   // the 2-norm is the magnitude of the larges eigenvalue.
00121   // This is not yet implemented.
00122   static real matrixnorm_frobenius(POINT const& A1, POINT const& A2, POINT const& A3)
00123     {
00124       return sqrt(  squared_norm_2(A1) 
00125                   + squared_norm_2(A2)
00126                   + squared_norm_2(A3));
00127     }
00128 
00129   struct Norm_frobenius {
00130     real operator()(POINT const& A1, POINT const& A2, POINT const& A3) const
00131       { return matrixnorm_frobenius(A1,A2,A3);}
00132   };
00133   typedef Norm_frobenius Norm_F;
00134 
00135   struct Norm_1 {
00136     real operator()(POINT const& A1, POINT const& A2, POINT const& A3) const
00137       { return matrixnorm_1(A1,A2,A3);}
00138   };
00139   typedef Norm_1 Norm_columnsum;
00140 
00141   struct Norm_infinity {
00142     real operator()(POINT const& A1, POINT const& A2, POINT const& A3) const
00143       { return matrixnorm_infinity(A1,A2,A3);}
00144   };
00145   typedef Norm_infinity Norm_oo;
00146   typedef Norm_infinity Norm_rowsum;
00147 
00148   template<class NORM>
00149   static real condition(POINT const& A1, POINT const& A2, POINT const& A3,
00150                         NORM const& N)
00151     {
00152       POINT inv[3];
00153       for(int i = 0; i<3; ++i)
00154         pt::ConstructWithDim(3,inv[i]);
00155       // FIXME: should be
00156       // if (det3(inv[0],inv[1],inv[2]) <= numeric_traits<real>::epsilon()*max(|Ai|)
00157       //   return numeric_traits<real>::infinity();
00158       double eps = 1.0/1.e10;
00159       if(fabs(det3(A1,A2,A3)) <= (eps*(norm_2(A1) + norm_2(A2) + norm_2(A3))))
00160         return 1.0e100;
00161      
00162       inverse(A1,A2,A3,inv[0],inv[1],inv[2]);
00163       return N(A1,A2,A3) * N(inv[0],inv[1],inv[2]);
00164     }
00165 
00166 
00167 };
00168 
00169 
00170 
00171 #endif
00172 
00173 
00174 
00175 
00176
Geometry/primitives3d.h
