src/include/vec3.h

   1 #ifndef _vec3_h_
   2 #define _vec3_h_
   3 #include <math.h>
   4 #include "petsc.h"
   5
   6 //--------------------------------------------------------------------
   7 //this is the basic operate of vec and mat. speed is the most important thing
   8 typedef struct
   9 {
  10         PetscScalar v[3];
  11 }Vec3;
  12
  13 typedef struct
  14 {
  15         PetscScalar m[9];
  16 }Mat3;
  17
  18 inline void Set_Vec3_zero(Vec3 &vec)
  19 {
  20         vec.v[0]=0;
  21         vec.v[1]=0;
  22         vec.v[2]=0;
  23 }
  24
  25 inline void Set_Vec3(Vec3 &vec,PetscScalar value)
  26 {
  27         vec.v[0]=value;
  28         vec.v[1]=value;
  29         vec.v[2]=value;
  30 }
  31
  32 inline void Set_Vec3(Vec3 &vec,PetscScalar v1, PetscScalar v2,PetscScalar v3)
  33 {
  34         vec.v[0]=v1;
  35         vec.v[1]=v2;
  36         vec.v[2]=v3;
  37 }
  38
  39 inline Vec3 operator +(const Vec3 &vec1,const Vec3 &vec2)
  40 {
  41         Vec3 sum;
  42         sum.v[0]=vec1.v[0]+vec2.v[0];
  43         sum.v[1]=vec1.v[1]+vec2.v[1];
  44         sum.v[2]=vec1.v[2]+vec2.v[2];
  45         return sum;
  46 }
  47
  48 inline Vec3 operator +(const Vec3 &vec1,const PetscScalar value)
  49 {
  50         Vec3 sum;
  51         sum.v[0]=vec1.v[0]+value;
  52         sum.v[1]=vec1.v[1]+value;
  53         sum.v[2]=vec1.v[2]+value;
  54         return sum;
  55 }
  56
  57 inline Vec3 operator +(const PetscScalar value,const Vec3 &vec1)
  58 {
  59         Vec3 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         return sum;
  64 }
  65
  66 inline Vec3 operator -(const Vec3 &vec1,const Vec3 &vec2)
  67 {
  68         Vec3 sub;
  69         sub.v[0]=vec1.v[0]-vec2.v[0];
  70         sub.v[1]=vec1.v[1]-vec2.v[1];
  71         sub.v[2]=vec1.v[2]-vec2.v[2];
  72         return sub;
  73 }
  74
  75 inline Vec3 operator -(const Vec3 &vec1,const PetscScalar value)
  76 {
  77         Vec3 sum;
  78         sum.v[0]=vec1.v[0]-value;
  79         sum.v[1]=vec1.v[1]-value;
  80         sum.v[2]=vec1.v[2]-value;
  81         return sum;
  82 }
  83
  84 inline Vec3 operator -(const PetscScalar value,const Vec3 &vec1)
  85 {
  86         Vec3 sum;
  87         sum.v[0]=value-vec1.v[0];
  88         sum.v[1]=value-vec1.v[1];
  89         sum.v[2]=value-vec1.v[2];
  90         return sum;
  91 }
  92
  93 inline Vec3 operator -(const Vec3 &vec1)
  94 {
  95         Vec3 sub;
  96         sub.v[0]=-vec1.v[0];
  97         sub.v[1]=-vec1.v[1];
  98         sub.v[2]=-vec1.v[2];
  99         return sub;
 100 }
 101
 102 inline Vec3 operator *(const Vec3 &vec,const PetscScalar &value)
 103 {
 104         Vec3 result;
 105         result.v[0]=vec.v[0]*value;
 106         result.v[1]=vec.v[1]*value;
 107         result.v[2]=vec.v[2]*value;
 108         return result;
 109 }
 110
 111 inline Vec3 operator *(const PetscScalar &value,const Vec3 &vec)
 112 {
 113         Vec3 result;
 114         result.v[0]=vec.v[0]*value;
 115         result.v[1]=vec.v[1]*value;
 116         result.v[2]=vec.v[2]*value;
 117         return result;
 118 }
 119
 120 inline Vec3 operator /(const Vec3 &vec,const PetscScalar &value)
 121 {
 122         Vec3 result;
 123         result.v[0]=vec.v[0]/value;
 124         result.v[1]=vec.v[1]/value;
 125         result.v[2]=vec.v[2]/value;
 126         return result;
 127 }
 128
 129 inline Vec3 operator /(const PetscScalar &value,const Vec3 &vec)
 130 {
 131         Vec3 result;
 132         result.v[0]=value/vec.v[0];
 133         result.v[1]=value/vec.v[1];
 134         result.v[2]=value/vec.v[2];
 135         return result;
 136 }
 137
 138 inline Vec3 operator *(const Vec3 &vec1,const Vec3 &vec2)
 139 {
 140         Vec3 result;
 141         result.v[0]=vec1.v[0]*vec2.v[0];
 142         result.v[1]=vec1.v[1]*vec2.v[1];
 143         result.v[2]=vec1.v[2]*vec2.v[2];
 144         return result;
 145 }
 146
 147 inline Vec3 operator /(const Vec3 &vec1,const Vec3 &vec2)
 148 {
 149         Vec3 result;
 150         result.v[0]=vec1.v[0]/vec2.v[0];
 151         result.v[1]=vec1.v[1]/vec2.v[1];
 152         result.v[2]=vec1.v[2]/vec2.v[2];
 153         return result;
 154 }
 155
 156
 157 //----------------------------------------------------
 158 inline void  Set_Mat3_zero(Mat3 &A)
 159 {
 160          A.m[0]=0;       A.m[1]=0;       A.m[2]=0;
 161          A.m[3]=0;       A.m[4]=0;       A.m[5]=0;
 162          A.m[6]=0;       A.m[7]=0;       A.m[8]=0;
 163 }
 164
 165 inline void Set_Mat3_I(Mat3 &A)
 166 {
 167          A.m[0]=1;       A.m[1]=0;       A.m[2]=0;
 168          A.m[3]=0;       A.m[4]=1;       A.m[5]=0;
 169          A.m[6]=0;       A.m[7]=0;       A.m[8]=1;
 170
 171 }
 172
 173 inline Mat3 Set_Mat3_Scale(const Mat3 &A, const Vec3 &a)
 174 {
 175         Mat3 C;
 176         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];
 177         C.m[3]=A.m[3]*a.v[1];    C.m[4]=A.m[4]*a.v[1];    C.m[5]=A.m[5]*a.v[1];
 178         C.m[6]=A.m[6]*a.v[2];    C.m[7]=A.m[7]*a.v[2];    C.m[8]=A.m[8]*a.v[2];
 179         return C;
 180
 181 }
 182
 183 inline Mat3 operator +(const Mat3 &A, const Mat3 &B)
 184 {
 185         Mat3 C;
 186         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];
 187         C.m[3]=A.m[3]+B.m[3];    C.m[4]=A.m[4]+B.m[4];    C.m[5]=A.m[5]+B.m[5];
 188         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];
 189         return C;
 190 }
 191
 192 inline Mat3 operator -(const Mat3 &A, const Mat3 &B)
 193 {
 194         Mat3 C;
 195         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];
 196         C.m[3]=A.m[3]-B.m[3];    C.m[4]=A.m[4]-B.m[4];    C.m[5]=A.m[5]-B.m[5];
 197         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];
 198         return C;
 199 }
 200
 201 inline Mat3  operator -(const Mat3 &A)
 202 {
 203         Mat3 C;
 204         C.m[0]=-A.m[0];    C.m[1]=-A.m[1];    C.m[2]=-A.m[2];
 205         C.m[3]=-A.m[3];    C.m[4]=-A.m[4];    C.m[5]=-A.m[5];
 206         C.m[6]=-A.m[6];    C.m[7]=-A.m[7];    C.m[8]=-A.m[8];
 207         return C;
 208 }
 209
 210 inline Mat3  operator *(const Mat3 &A,const PetscScalar &a)
 211 {
 212         Mat3 C;
 213         C.m[0]=a*A.m[0];    C.m[1]=a*A.m[1];    C.m[2]=a*A.m[2];
 214         C.m[3]=a*A.m[3];    C.m[4]=a*A.m[4];    C.m[5]=a*A.m[5];
 215         C.m[6]=a*A.m[6];    C.m[7]=a*A.m[7];    C.m[8]=a*A.m[8];
 216         return C;
 217 }
 218
 219 inline Mat3  operator *(const PetscScalar &a,const Mat3 &A)
 220 {
 221         Mat3 C;
 222         C.m[0]=a*A.m[0];    C.m[1]=a*A.m[1];    C.m[2]=a*A.m[2];
 223         C.m[3]=a*A.m[3];    C.m[4]=a*A.m[4];    C.m[5]=a*A.m[5];
 224         C.m[6]=a*A.m[6];    C.m[7]=a*A.m[7];    C.m[8]=a*A.m[8];
 225         return C;
 226 }
 227
 228 inline Mat3  operator /(const Mat3 &A, const PetscScalar &a)
 229 {
 230         Mat3 C;
 231         C.m[0]=A.m[0]/a;    C.m[1]=A.m[1]/a;    C.m[2]=A.m[2]/a;
 232         C.m[3]=A.m[3]/a;    C.m[4]=A.m[4]/a;    C.m[5]=A.m[5]/a;
 233         C.m[6]=A.m[6]/a;    C.m[7]=A.m[7]/a;    C.m[8]=A.m[8]/a;
 234         return C;
 235 }
 236
 237 inline Mat3  operator /(const Mat3 &A, const Vec3 &a)
 238 {
 239         Mat3 C;
 240         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];
 241         C.m[3]=A.m[3]/a.v[1];    C.m[4]=A.m[4]/a.v[1];    C.m[5]=A.m[5]/a.v[1];
 242         C.m[6]=A.m[6]/a.v[2];    C.m[7]=A.m[7]/a.v[2];    C.m[8]=A.m[8]/a.v[2];
 243         return C;
 244 }
 245
 246 inline Vec3 operator *(const Mat3 &A, const Vec3 &x)
 247 {
 248         Vec3 b;
 249         b.v[0]=A.m[0]*x.v[0] + A.m[1]*x.v[1] + A.m[2]*x.v[2];
 250         b.v[1]=A.m[3]*x.v[0] + A.m[4]*x.v[1] + A.m[5]*x.v[2];
 251         b.v[2]=A.m[6]*x.v[0] + A.m[7]*x.v[1] + A.m[8]*x.v[2];
 252         return b;
 253 }
 254
 255 inline Mat3 operator *(const Mat3 &A, const Mat3 &B)
 256 {
 257         Mat3 C;
 258         C.m[0]=A.m[0]*B.m[0] + A.m[1]*B.m[3] + A.m[2]*B.m[6];
 259         C.m[1]=A.m[0]*B.m[1] + A.m[1]*B.m[4] + A.m[2]*B.m[7];
 260         C.m[2]=A.m[0]*B.m[2] + A.m[1]*B.m[5] + A.m[2]*B.m[8];
 261
 262         C.m[3]=A.m[3]*B.m[0] + A.m[4]*B.m[3] + A.m[5]*B.m[6];
 263         C.m[4]=A.m[3]*B.m[1] + A.m[4]*B.m[4] + A.m[5]*B.m[7];
 264         C.m[5]=A.m[3]*B.m[2] + A.m[4]*B.m[5] + A.m[5]*B.m[8];
 265
 266         C.m[6]=A.m[6]*B.m[0] + A.m[7]*B.m[3] + A.m[8]*B.m[6];
 267         C.m[7]=A.m[6]*B.m[1] + A.m[7]*B.m[4] + A.m[8]*B.m[7];
 268         C.m[8]=A.m[6]*B.m[2] + A.m[7]*B.m[5] + A.m[8]*B.m[8];
 269         return C;
 270 }
 271
 272
 273 inline void mat3_transport(Mat3 &A)
 274 {
 275
 276         PetscScalar swap;
 277         for (int i=0;i<3;i++)
 278           for (int j=0;j<=i;j++)
 279             {
 280              swap=A.m[i*3+j];
 281              A.m[i*3+j]=A.m[j*3+i];
 282              A.m[j*3+i]=swap;
 283             }
 284 }
 285
 286 inline ostream& operator << ( ostream& out, const Mat3 & A)
 287 {
 288     for (int i=0;i<3;i++)
 289     {
 290       for (int j=0;j<3;j++)
 291         out<<A.m[i*3+j]<<" ";
 292       out<<endl;
 293     }
 294     return out;
 295 }
 296
 297 #endif
 298