src/include/vec5.h

   1 #ifndef _vec5_h_
   2 #define _vec5_h_
   3 #include <math.h>
   4 #include "petsc.h"
   5
   6 //this is the basic operate of vec and mat. speed is the most important thing
   7 typedef struct
   8 {
   9   PetscScalar v[5];
  10 }
  11 Vec5;
  12
  13 typedef struct
  14 {
  15   PetscScalar m[25];
  16 }
  17 Mat5;
  18
  19 inline void Set_Vec5_zero(Vec5 &vec)
  20 {
  21   vec.v[0]=0;
  22   vec.v[1]=0;
  23   vec.v[2]=0;
  24   vec.v[3]=0;
  25   vec.v[4]=0;
  26 }
  27
  28 inline void Set_Vec5(Vec5 &vec,PetscScalar value)
  29 {
  30   vec.v[0]=value;
  31   vec.v[1]=value;
  32   vec.v[2]=value;
  33   vec.v[3]=value;
  34   vec.v[4]=value;
  35 }
  36
  37 inline void Set_Vec5(Vec5 &vec,PetscScalar v1,PetscScalar v2,PetscScalar v3,PetscScalar v4, PetscScalar v5)
  38 {
  39   vec.v[0]=v1;
  40   vec.v[1]=v2;
  41   vec.v[2]=v3;
  42   vec.v[3]=v4;
  43   vec.v[4]=v5;
  44 }
  45
  46 inline Vec5 operator +(const Vec5 &vec1,const Vec5 &vec2)
  47 {
  48   Vec5 sum;
  49   sum.v[0]=vec1.v[0]+vec2.v[0];
  50   sum.v[1]=vec1.v[1]+vec2.v[1];
  51   sum.v[2]=vec1.v[2]+vec2.v[2];
  52   sum.v[3]=vec1.v[3]+vec2.v[3];
  53   sum.v[4]=vec1.v[4]+vec2.v[4];
  54   return sum;
  55 }
  56
  57 inline Vec5 operator +(const Vec5 &vec1,const PetscScalar value)
  58 {
  59   Vec5 sum;
  60   sum.v[0]=vec1.v[0]+value;
  61   sum.v[1]=vec1.v[1]+value;
  62   sum.v[2]=vec1.v[2]+value;
  63   sum.v[3]=vec1.v[3]+value;
  64   sum.v[4]=vec1.v[4]+value;
  65   return sum;
  66 }
  67
  68 inline Vec5 operator +(const PetscScalar value,const Vec5 &vec1)
  69 {
  70   Vec5 sum;
  71   sum.v[0]=vec1.v[0]+value;
  72   sum.v[1]=vec1.v[1]+value;
  73   sum.v[2]=vec1.v[2]+value;
  74   sum.v[3]=vec1.v[3]+value;
  75   sum.v[4]=vec1.v[4]+value;
  76   return sum;
  77 }
  78
  79 inline Vec5 operator -(const Vec5 &vec1,const Vec5 &vec2)
  80 {
  81   Vec5 sub;
  82   sub.v[0]=vec1.v[0]-vec2.v[0];
  83   sub.v[1]=vec1.v[1]-vec2.v[1];
  84   sub.v[2]=vec1.v[2]-vec2.v[2];
  85   sub.v[3]=vec1.v[3]-vec2.v[3];
  86   sub.v[4]=vec1.v[4]-vec2.v[4];
  87   return sub;
  88 }
  89
  90 inline Vec5 operator -(const Vec5 &vec1,const PetscScalar value)
  91 {
  92   Vec5 sum;
  93   sum.v[0]=vec1.v[0]-value;
  94   sum.v[1]=vec1.v[1]-value;
  95   sum.v[2]=vec1.v[2]-value;
  96   sum.v[3]=vec1.v[3]-value;
  97   sum.v[4]=vec1.v[4]-value;
  98   return sum;
  99 }
 100
 101 inline Vec5 operator -(const PetscScalar value,const Vec5 &vec1)
 102 {
 103   Vec5 sum;
 104   sum.v[0]=value-vec1.v[0];
 105   sum.v[1]=value-vec1.v[1];
 106   sum.v[2]=value-vec1.v[2];
 107   sum.v[3]=value-vec1.v[3];
 108   sum.v[4]=value-vec1.v[4];
 109   return sum;
 110 }
 111
 112 inline Vec5 operator *(const Vec5 &vec,const PetscScalar &value)
 113 {
 114   Vec5 result;
 115   result.v[0]=vec.v[0]*value;
 116   result.v[1]=vec.v[1]*value;
 117   result.v[2]=vec.v[2]*value;
 118   result.v[3]=vec.v[3]*value;
 119   result.v[4]=vec.v[4]*value;
 120   return result;
 121 }
 122
 123 inline Vec5 operator *(const PetscScalar &value,const Vec5 &vec)
 124 {
 125   Vec5 result;
 126   result.v[0]=vec.v[0]*value;
 127   result.v[1]=vec.v[1]*value;
 128   result.v[2]=vec.v[2]*value;
 129   result.v[3]=vec.v[3]*value;
 130   result.v[4]=vec.v[4]*value;
 131   return result;
 132 }
 133
 134 inline Vec5 operator /(const Vec5 &vec,const PetscScalar &value)
 135 {
 136   Vec5 result;
 137   result.v[0]=vec.v[0]/value;
 138   result.v[1]=vec.v[1]/value;
 139   result.v[2]=vec.v[2]/value;
 140   result.v[3]=vec.v[3]/value;
 141   result.v[4]=vec.v[4]/value;
 142   return result;
 143 }
 144
 145 inline Vec5 operator /(const PetscScalar &value,const Vec5 &vec)
 146 {
 147   Vec5 result;
 148   result.v[0]=value/vec.v[0];
 149   result.v[1]=value/vec.v[1];
 150   result.v[2]=value/vec.v[2];
 151   result.v[3]=value/vec.v[3];
 152   result.v[4]=value/vec.v[4];
 153   return result;
 154 }
 155
 156 inline Vec5 operator *(const Vec5 &vec1,const Vec5 &vec2)
 157 {
 158   Vec5 result;
 159   result.v[0]=vec1.v[0]*vec2.v[0];
 160   result.v[1]=vec1.v[1]*vec2.v[1];
 161   result.v[2]=vec1.v[2]*vec2.v[2];
 162   result.v[3]=vec1.v[3]*vec2.v[3];
 163   result.v[4]=vec1.v[4]*vec2.v[4];
 164   return result;
 165 }
 166
 167 inline Vec5 operator /(const Vec5 &vec1,const Vec5 &vec2)
 168 {
 169   Vec5 result;
 170   result.v[0]=vec1.v[0]/vec2.v[0];
 171   result.v[1]=vec1.v[1]/vec2.v[1];
 172   result.v[2]=vec1.v[2]/vec2.v[2];
 173   result.v[3]=vec1.v[3]/vec2.v[3];
 174   result.v[4]=vec1.v[4]/vec2.v[4];
 175   return result;
 176 }
 177
 178 //--------------------------------------------------------------------
 179
 180 inline PetscScalar vmaxvalue(const Vec5 &a)
 181 {
 182   PetscScalar v=a.v[0];
 183   if(a.v[1]>v)  v=a.v[1];
 184   if(a.v[2]>v)  v=a.v[2];
 185   if(a.v[3]>v)  v=a.v[3];
 186   if(a.v[4]>v)  v=a.v[4];
 187   return v;
 188 }
 189
 190 inline Vec5 vmax(const Vec5 &a,const Vec5 &b)
 191 {
 192   Vec5 result;
 193   result.v[0]=a.v[0]>b.v[0] ? a.v[0]:b.v[0];
 194   result.v[1]=a.v[1]>b.v[1] ? a.v[1]:b.v[1];
 195   result.v[2]=a.v[2]>b.v[2] ? a.v[2]:b.v[2];
 196   result.v[3]=a.v[3]>b.v[3] ? a.v[3]:b.v[3];
 197   result.v[4]=a.v[4]>b.v[4] ? a.v[4]:b.v[4];
 198   return result;
 199 }
 200
 201 inline Vec5 vmin(const Vec5 &a,const Vec5 &b)
 202 {
 203   Vec5 result;
 204   result.v[0]=a.v[0]<b.v[0] ? a.v[0]:b.v[0];
 205   result.v[1]=a.v[1]<b.v[1] ? a.v[1]:b.v[1];
 206   result.v[2]=a.v[2]<b.v[2] ? a.v[2]:b.v[2];
 207   result.v[3]=a.v[3]<b.v[3] ? a.v[3]:b.v[3];
 208   result.v[4]=a.v[4]<b.v[4] ? a.v[4]:b.v[4];
 209   return result;
 210 }
 211
 212 inline Vec5 vabs(const Vec5 &a)
 213 {
 214   Vec5 result;
 215   result.v[0]=fabs(a.v[0]);
 216   result.v[1]=fabs(a.v[1]);
 217   result.v[2]=fabs(a.v[2]);
 218   result.v[3]=fabs(a.v[3]);
 219   result.v[4]=fabs(a.v[4]);
 220   return result;
 221 }
 222
 223 inline Vec5 vsign(const Vec5 &a)
 224 {
 225   Vec5 result;
 226   result.v[0]=a.v[0]>0.0? 1.0:-1.0;
 227   result.v[1]=a.v[1]>0.0? 1.0:-1.0;
 228   result.v[2]=a.v[2]>0.0? 1.0:-1.0;
 229   result.v[3]=a.v[3]>0.0? 1.0:-1.0;
 230   result.v[4]=a.v[4]>0.0? 1.0:-1.0;
 231   return result;
 232 }
 233
 234 //-------------------------------------------------------------------------
 235
 236 inline void  Set_Mat5_zero(Mat5 &A)
 237 {
 238   A.m[0]=0;      A.m[1]=0;       A.m[2]=0;       A.m[3]=0;       A.m[4]=0;
 239   A.m[5]=0;      A.m[6]=0;       A.m[7]=0;       A.m[8]=0;       A.m[9]=0;
 240   A.m[10]=0;     A.m[11]=0;      A.m[12]=0;      A.m[13]=0;      A.m[14]=0;
 241   A.m[15]=0;     A.m[16]=0;      A.m[17]=0;      A.m[18]=0;      A.m[19]=0;
 242   A.m[20]=0;     A.m[21]=0;      A.m[22]=0;      A.m[23]=0;      A.m[24]=0;
 243 }
 244
 245 inline void Set_Mat5_I(Mat5 &A)
 246 {
 247   A.m[0]=1;      A.m[1]=0;       A.m[2]=0;       A.m[3]=0;       A.m[4]=0;
 248   A.m[5]=0;      A.m[6]=1;       A.m[7]=0;       A.m[8]=0;       A.m[9]=0;
 249   A.m[10]=0;     A.m[11]=0;      A.m[12]=1;      A.m[13]=0;      A.m[14]=0;
 250   A.m[15]=0;     A.m[16]=0;      A.m[17]=0;      A.m[18]=1;      A.m[19]=0;
 251   A.m[20]=0;     A.m[21]=0;      A.m[22]=0;      A.m[23]=0;      A.m[24]=1;
 252
 253 }
 254
 255 inline Mat5 operator +(const Mat5 &A, const Mat5 &B)
 256 {
 257   Mat5 C;
 258   C.m[0]=A.m[0]+B.m[0];    C.m[1]=A.m[1]+B.m[1];    C.m[2]=A.m[2]+B.m[2];    C.m[3]=A.m[3]+B.m[3];    C.m[4]=A.m[4]+B.m[4];
 259   C.m[5]=A.m[5]+B.m[5];    C.m[6]=A.m[6]+B.m[6];    C.m[7]=A.m[7]+B.m[7];    C.m[8]=A.m[8]+B.m[8];    C.m[9]=A.m[9]+B.m[9];
 260   C.m[10]=A.m[10]+B.m[10]; C.m[11]=A.m[11]+B.m[11]; C.m[12]=A.m[12]+B.m[12]; C.m[13]=A.m[13]+B.m[13]; C.m[14]=A.m[14]+B.m[14];
 261   C.m[15]=A.m[15]+B.m[15]; C.m[16]=A.m[16]+B.m[16]; C.m[17]=A.m[17]+B.m[17]; C.m[18]=A.m[18]+B.m[18]; C.m[19]=A.m[19]+B.m[19];
 262   C.m[20]=A.m[20]+B.m[20]; C.m[21]=A.m[21]+B.m[21]; C.m[22]=A.m[22]+B.m[22]; C.m[23]=A.m[23]+B.m[23]; C.m[24]=A.m[24]+B.m[24];
 263   return C;
 264 }
 265
 266 inline Mat5 operator -(const Mat5 &A, const Mat5 &B)
 267 {
 268   Mat5 C;
 269   C.m[0]=A.m[0]-B.m[0];    C.m[1]=A.m[1]-B.m[1];    C.m[2]=A.m[2]-B.m[2];    C.m[3]=A.m[3]-B.m[3];    C.m[4]=A.m[4]-B.m[4];
 270   C.m[5]=A.m[5]-B.m[5];    C.m[6]=A.m[6]-B.m[6];    C.m[7]=A.m[7]-B.m[7];    C.m[8]=A.m[8]-B.m[8];    C.m[9]=A.m[9]-B.m[9];
 271   C.m[10]=A.m[10]-B.m[10]; C.m[11]=A.m[11]-B.m[11]; C.m[12]=A.m[12]-B.m[12]; C.m[13]=A.m[13]-B.m[13]; C.m[14]=A.m[14]-B.m[14];
 272   C.m[15]=A.m[15]-B.m[15]; C.m[16]=A.m[16]-B.m[16]; C.m[17]=A.m[17]-B.m[17]; C.m[18]=A.m[18]-B.m[18]; C.m[19]=A.m[19]-B.m[19];
 273   C.m[20]=A.m[20]-B.m[20]; C.m[21]=A.m[21]-B.m[21]; C.m[22]=A.m[22]-B.m[22]; C.m[23]=A.m[23]-B.m[23]; C.m[24]=A.m[24]-B.m[24];
 274   return C;
 275 }
 276
 277 inline Mat5  operator *(const Mat5 &A,const PetscScalar &a)
 278 {
 279   Mat5 C;
 280   C.m[0]=a*A.m[0];    C.m[1]=a*A.m[1];    C.m[2]=a*A.m[2];    C.m[3]=a*A.m[3];    C.m[4]=a*A.m[4];
 281   C.m[5]=a*A.m[5];    C.m[6]=a*A.m[6];    C.m[7]=a*A.m[7];    C.m[8]=a*A.m[8];    C.m[9]=a*A.m[9];
 282   C.m[10]=a*A.m[10];  C.m[11]=a*A.m[11];  C.m[12]=a*A.m[12];  C.m[13]=a*A.m[13];  C.m[14]=a*A.m[14];
 283   C.m[15]=a*A.m[15];  C.m[16]=a*A.m[16];  C.m[17]=a*A.m[17];  C.m[18]=a*A.m[18];  C.m[19]=a*A.m[19];
 284   C.m[20]=a*A.m[20];  C.m[21]=a*A.m[21];  C.m[22]=a*A.m[22];  C.m[23]=a*A.m[23];  C.m[24]=a*A.m[24];
 285   return C;
 286 }
 287
 288 inline Mat5  operator *(const PetscScalar &a,const Mat5 &A)
 289 {
 290   Mat5 C;
 291   C.m[0]=a*A.m[0];    C.m[1]=a*A.m[1];    C.m[2]=a*A.m[2];    C.m[3]=a*A.m[3];    C.m[4]=a*A.m[4];
 292   C.m[5]=a*A.m[5];    C.m[6]=a*A.m[6];    C.m[7]=a*A.m[7];    C.m[8]=a*A.m[8];    C.m[9]=a*A.m[9];
 293   C.m[10]=a*A.m[10];  C.m[11]=a*A.m[11];  C.m[12]=a*A.m[12];  C.m[13]=a*A.m[13];  C.m[14]=a*A.m[14];
 294   C.m[15]=a*A.m[15];  C.m[16]=a*A.m[16];  C.m[17]=a*A.m[17];  C.m[18]=a*A.m[18];  C.m[19]=a*A.m[19];
 295   C.m[20]=a*A.m[20];  C.m[21]=a*A.m[21];  C.m[22]=a*A.m[22];  C.m[23]=a*A.m[23];  C.m[24]=a*A.m[24];
 296   return C;
 297 }
 298
 299 inline Mat5  operator /(const Mat5 &A, const PetscScalar &a)
 300 {
 301   Mat5 C;
 302   C.m[0]=A.m[0]/a;    C.m[1]=A.m[1]/a;    C.m[2]=A.m[2]/a;    C.m[3]=A.m[3]/a;    C.m[4]=A.m[4]/a;
 303   C.m[5]=A.m[5]/a;    C.m[6]=A.m[6]/a;    C.m[7]=A.m[7]/a;    C.m[8]=A.m[8]/a;    C.m[9]=A.m[9]/a;
 304   C.m[10]=A.m[10]/a;  C.m[11]=A.m[11]/a;  C.m[12]=A.m[12]/a;  C.m[13]=A.m[13]/a;  C.m[14]=A.m[14]/a;
 305   C.m[15]=A.m[15]/a;  C.m[16]=A.m[16]/a;  C.m[17]=A.m[17]/a;  C.m[18]=A.m[18]/a;  C.m[19]=A.m[19]/a;
 306   C.m[20]=A.m[20]/a;  C.m[21]=A.m[21]/a;  C.m[22]=A.m[22]/a;  C.m[23]=A.m[23]/a;  C.m[24]=A.m[24]/a;
 307   return C;
 308 }
 309
 310 inline Mat5  operator /(const Mat5 &A, const Vec5 &a)
 311 {
 312   Mat5 C;
 313   C.m[0]=A.m[0]/a.v[0];    C.m[1]=A.m[1]/a.v[0];    C.m[2]=A.m[2]/a.v[0];    C.m[3]=A.m[3]/a.v[0];    C.m[4]=A.m[4]/a.v[0];
 314   C.m[5]=A.m[5]/a.v[1];    C.m[6]=A.m[6]/a.v[1];    C.m[7]=A.m[7]/a.v[1];    C.m[8]=A.m[8]/a.v[1];    C.m[9]=A.m[9]/a.v[1];
 315   C.m[10]=A.m[10]/a.v[2];  C.m[11]=A.m[11]/a.v[2];  C.m[12]=A.m[12]/a.v[2];  C.m[13]=A.m[13]/a.v[2];  C.m[14]=A.m[14]/a.v[2];
 316   C.m[15]=A.m[15]/a.v[3];  C.m[16]=A.m[16]/a.v[3];  C.m[17]=A.m[17]/a.v[3];  C.m[18]=A.m[18]/a.v[3];  C.m[19]=A.m[19]/a.v[3];
 317   C.m[20]=A.m[20]/a.v[4];  C.m[21]=A.m[21]/a.v[4];  C.m[22]=A.m[22]/a.v[4];  C.m[23]=A.m[23]/a.v[4];  C.m[24]=A.m[24]/a.v[4];
 318   return C;
 319 }
 320
 321 inline Vec5 operator *(const Mat5 &A, const Vec5 &x)
 322 {
 323   Vec5 b;
 324   b.v[0]=A.m[0]*x.v[0]+A.m[1]*x.v[1]+A.m[2]*x.v[2]+A.m[3]*x.v[3]+A.m[4]*x.v[4];
 325   b.v[1]=A.m[5]*x.v[0]+A.m[6]*x.v[1]+A.m[7]*x.v[2]+A.m[8]*x.v[3]+A.m[9]*x.v[4];
 326   b.v[2]=A.m[10]*x.v[0]+A.m[11]*x.v[1]+A.m[12]*x.v[2]+A.m[13]*x.v[3]+A.m[14]*x.v[4];
 327   b.v[3]=A.m[15]*x.v[0]+A.m[16]*x.v[1]+A.m[17]*x.v[2]+A.m[18]*x.v[3]+A.m[19]*x.v[4];
 328   b.v[4]=A.m[20]*x.v[0]+A.m[21]*x.v[1]+A.m[22]*x.v[2]+A.m[23]*x.v[3]+A.m[24]*x.v[4];
 329   return b;
 330 }
 331
 332
 333 inline void mat5_transport(Mat5 &A)
 334 {
 335   PetscScalar swap;
 336   for (int i=0;i<5;i++)
 337     for (int j=0;j<=i;j++)
 338     {
 339       swap=A.m[i*5+j];
 340       A.m[i*5+j]=A.m[j*5+i];
 341       A.m[j*5+i]=swap;
 342     }
 343 }
 344
 345 inline ostream& operator << ( ostream& out, const Mat5 & A)
 346 {
 347     for (int i=0;i<5;i++)
 348     {
 349       for (int j=0;j<5;j++)
 350         out<<A.m[i*5+j]<<" ";
 351       out<<endl;
 352     }
 353     return out;
 354 }
 355
 356 #endif
 357