DENDRO: /net/www/csela/dendro/releases/Dendro/examples/include/stiffnessMatrix.h Source File

00001 
00015 #ifndef _STIFFNESSMATRIX_H_
00016 #define _STIFFNESSMATRIX_H_
00017 
00018 #include "feMatrix.h"
00019 
00033 class stiffnessMatrix : public ot::fem::feMatrix<stiffnessMatrix>
00034 {
00035   public:       
00036     stiffnessMatrix(daType da);
00047     inline bool ElementalMatVec(int i, int j, int k, PetscScalar ***in, PetscScalar ***out, double scale);
00048     inline bool ElementalMatVec(unsigned int idx, PetscScalar *in, PetscScalar *out, double scale);
00049 
00050     inline bool ElementalMatGetDiagonal(int i, int j, int k, PetscScalar ***diag, double scale);
00051     inline bool ElementalMatGetDiagonal(unsigned int idx, PetscScalar *diag, double scale);
00052     
00053     inline bool initStencils();
00054 
00055     bool preMatVec();
00056     bool postMatVec();
00057 
00058     void setNuVec(Vec nv) {
00059       nuvec = nv;
00060     }
00061 
00062    private:
00063     void*               m_nuarray; /* Diffusion coefficient array*/
00064     Vec                 nuvec;     
00065 
00066     double              m_dHx;
00067          double     m_nuval;
00068     double xFac, yFac, zFac;
00069     unsigned int maxD;
00070 
00071 };
00072 
00073 stiffnessMatrix::stiffnessMatrix(daType da) {
00074 #ifdef __DEBUG__
00075   assert ( ( da == PETSC ) || ( da == OCT ) );
00076 #endif
00077   m_daType = da;
00078   m_DA          = NULL;
00079   m_octDA       = NULL;
00080   m_stencil     = NULL;
00081 
00082   // initialize the stencils ...
00083   initStencils();
00084   if (da == OCT)
00085     initOctLut();
00086 }
00087 
00088 bool stiffnessMatrix::initStencils() {
00089   typedef int* int1Ptr;
00090   typedef int** int2Ptr;
00091 
00092   if ( m_daType == PETSC) {
00093     int Bjk[8][8] = {
00094       { 64,   0,   0, -16,   0, -16, -16, -16},
00095       {0,  64, -16,   0, -16,   0, -16, -16},
00096       {0, -16,  64,   0, -16, -16,   0, -16},
00097       { -16,   0,   0,  64, -16, -16, -16,   0},
00098       { 0, -16, -16, -16,  64,   0,   0, -16},
00099       { -16,   0, -16, -16,   0,  64, -16,   0},
00100       { -16, -16,   0, -16,   0, -16,  64,   0},
00101       { -16, -16, -16,   0, -16,   0,   0,  64}
00102     };
00103     
00104     int** Ajk = new int1Ptr[8];
00105     for (int j=0;j<8;j++) {
00106       Ajk[j] = new int[8];
00107       for (int k=0;k<8;k++) {
00108         Ajk[j][k] = Bjk[j][k];
00109       }//end k
00110     }//end j
00111     m_stencil = Ajk;
00112 
00113   } else {
00114     int Bijk[8][8][8] = {
00115       //Type-0: No Hanging
00116       {{ 64,   0,   0, -16,   0, -16, -16, -16},
00117         {0,  64, -16,   0, -16,   0, -16, -16},
00118         {0, -16,  64,   0, -16, -16,   0, -16},
00119         { -16,   0,   0,  64, -16, -16, -16,   0},
00120         { 0, -16, -16, -16,  64,   0,   0, -16},
00121         { -16,   0, -16, -16,   0,  64, -16,   0},
00122         { -16, -16,   0, -16,   0, -16,  64,   0},
00123         { -16, -16, -16,   0, -16,   0,   0,  64}},
00124       //Type-1: Y Hanging
00125       {{  80,  -8,  16, -16,  -8, -24, -16, -24},
00126         {  -8,  64,  -8,   0, -16,   0, -16, -16},
00127         {  16,  -8,  16,   0,  -8,  -8,   0,  -8},
00128         { -16,   0,   0,  64, -16, -16, -16,   0},
00129         {  -8, -16,  -8, -16,  64,   0,   0, -16},
00130         { -24,   0,  -8, -16,   0,  64, -16,   0},
00131         { -16, -16,   0, -16,   0, -16,  64,   0},
00132         { -24, -16,  -8,   0, -16,   0,   0,  64}},
00133       //Type-2: X and Y Hanging
00134       {{  88,  12,  12, -16, -16, -24, -24, -32},
00135         {  12,  16,  -4,   0,  -8,   0,  -8,  -8},
00136         {  12,  -4,  16,   0,  -8,  -8,   0,  -8},
00137         { -16,   0,   0,  64, -16, -16, -16,   0},
00138         { -16,  -8,  -8, -16,  64,   0,   0, -16},
00139         { -24,   0,  -8, -16,   0,  64, -16,   0},
00140         { -24,  -8,   0, -16,   0, -16,  64,   0},
00141         { -32,  -8,  -8,   0, -16,   0,   0,  64}},
00142       //Type-3: X, Y and Z Hanging
00143       {{  88,   8,   8, -24,   8, -24, -24, -40},
00144         {   8,  16,  -4,   0,  -4,   0,  -8,  -8},
00145         {   8,  -4,  16,   0,  -4,  -8,   0,  -8},
00146         { -24,   0,   0,  64,  -8, -16, -16,   0},
00147         {   8,  -4,  -4,  -8,  16,   0,   0,  -8},
00148         { -24,   0,  -8, -16,   0,  64, -16,   0},
00149         { -24,  -8,   0, -16,   0, -16,  64,   0},
00150         { -40,  -8,  -8,   0,  -8,   0,   0,  64}},
00151       //Type-4: XY and X and Y Hanging
00152       {{  84,  12,  12,   0, -20, -28, -28, -32},
00153         {  12,  20,   0,   4, -12,  -4, -12,  -8},
00154         {  12,   0,  20,   4, -12, -12,  -4,  -8},
00155         {   0,   4,   4,   4,  -4,  -4,  -4,   0},
00156         { -20, -12, -12,  -4,  64,   0,   0, -16},
00157         { -28,  -4, -12,  -4,   0,  64, -16,   0},
00158         { -28, -12,  -4,  -4,   0, -16,  64,   0},
00159         { -32,  -8,  -8,   0, -16,   0,   0,  64}},
00160       //Type-5: XY and X and Y and Z Hanging
00161       {{  80,   6,   6,  -2,   6, -28, -28, -40},
00162         {   6,  20,   0,   4,  -6,  -4, -12,  -8},
00163         {   6,   0,  20,   4,  -6, -12,  -4,  -8},
00164         {  -2,   4,   4,   4,  -2,  -4,  -4,   0},
00165         {   6,  -6,  -6,  -2,  16,   0,   0,  -8},
00166         { -28,  -4, -12,  -4,   0,  64, -16,   0},
00167         { -28, -12,  -4,  -4,   0, -16,  64,   0},
00168         { -40,  -8,  -8,   0,  -8,   0,   0,  64}},
00169       //Type-6:XY and YZ and X and Y and Z Hanging
00170       {{  70,   3,   2,  -3,   3, -32,  -3, -40},
00171         {   3,  20,  -3,   4,  -9,  -4,  -3,  -8},
00172         {   2,  -3,  22,   3,  -3, -16,   3,  -8},
00173         {  -3,   4,   3,   4,  -3,  -4,  -1,   0},
00174         {   3,  -9,  -3,  -3,  20,  -4,   4,  -8},
00175         { -32,  -4, -16,  -4,  -4,  64,  -4,   0},
00176         {  -3,  -3,   3,  -1,   4,  -4,   4,   0},
00177         { -40,  -8,  -8,   0,  -8,   0,   0,  64}},
00178       //Type-7:All 6 Hanging
00179       {{  58,  -2,  -2,  -4,  -2,  -4,  -4, -40},
00180         {  -2,  22,  -7,   3,  -7,   3,  -4,  -8},
00181         {  -2,  -7,  22,   3,  -7,  -4,   3,  -8},
00182         {  -4,   3,   3,   4,  -4,  -1,  -1,   0},
00183         {  -2,  -7,  -7,  -4,  22,   3,   3,  -8},
00184         {  -4,   3,  -4,  -1,   3,   4,  -1,   0},
00185         {  -4,  -4,   3,  -1,   3,  -1,   4,   0},
00186         { -40,  -8,  -8,   0,  -8,   0,   0,  64}}
00187     };
00188     int ***Aijk = new int2Ptr[8];
00189     for (int i=0;i<8;i++) {
00190       Aijk[i] = new int1Ptr[8];
00191       for (int j=0;j<8;j++) {
00192         Aijk[i][j] = new int[8];
00193         for (int k=0;k<8;k++) {
00194           Aijk[i][j][k] = Bijk[i][j][k];
00195         }//end k
00196       }//end j
00197     }//end i
00198     m_stencil = Aijk;
00199   }
00200   return true;
00201 }
00202 
00203 bool stiffnessMatrix::preMatVec() {
00204   // nuVec should be set directly into matrix outside the loop ...
00205   if (m_daType == PETSC) {
00206     PetscScalar ***nuarray; //
00207          // int ierr;
00208          // ierr = VecNorm(nuvec,NORM_INFINITY,&m_nuval); CHKERRQ(ierr);
00209 
00210     int ierr = DAVecGetArray(m_DA, nuvec, &nuarray);
00211 
00212     m_nuarray = nuarray;
00213     // compute Hx
00214     PetscInt mx,my,mz;
00215     ierr = DAGetInfo(m_DA,0, &mx, &my, &mz, 0,0,0,0,0,0,0); CHKERRQ(ierr); 
00216     CHKERRQ(ierr);
00217     
00218     m_dHx = m_dLx/(mx -1);
00219 
00220     m_dHx /= 192.0;
00221   } else {
00222     PetscScalar *nuarray; 
00223     // Get nuarray
00224     m_octDA->vecGetBuffer(nuvec, nuarray, false, false, true,m_uiDof);
00225     m_nuarray = nuarray;
00226 
00227     // compute Hx
00228     // For octree Hx values will change per element, so has to be 
00229     // computed inside the loop.
00230          maxD = m_octDA->getMaxDepth();
00231 
00232     
00233     // Get the  x,y,z factors 
00234     xFac = 1.0/((double)(1u << (maxD-1)));
00235     if (m_octDA->getDimension() > 1) {
00236       yFac = 1.0/((double)(1u << (maxD-1)));
00237       if (m_octDA->getDimension() > 2) {
00238         zFac = 1.0/((double)(1u << (maxD-1)));
00239       }
00240     }
00241   }
00242 
00243   return true;
00244 }
00245 
00246 bool stiffnessMatrix::ElementalMatVec(unsigned int i, PetscScalar *in, PetscScalar *out, double scale) {
00247   unsigned int lev = m_octDA->getLevel(i);
00248   double hx = xFac*(1u << (maxD - lev));
00249   double hy = yFac*(1u << (maxD - lev));
00250   double hz = zFac*(1u << (maxD - lev));
00251 
00252   double fac11 = -hx*scale/192.0;
00253   
00254   stdElemType elemType;
00255   unsigned int idx[8];
00256 
00257   unsigned char hangingMask = m_octDA->getHangingNodeIndex(i);    //  alignElementAndVertices(m_octDA, elemType, idx);       
00258   int ***Aijk = (int ***)m_stencil;
00259   
00260   alignElementAndVertices(m_octDA, elemType, idx);       
00261 
00262   PetscScalar *nuarray = (PetscScalar *)m_nuarray;
00263   for (int k = 0;k < 8;k++) {
00264     for (int j=0;j<8;j++) {
00265                 double fac1 = nuarray[idx[k]]*fac11;
00266       out[m_uiDof*idx[k]] +=  (fac1*(Aijk[elemType][k][j]))*in[m_uiDof*idx[j]];
00267                 //if(hangingMask) std::cout << fac1*Aijk[elemType][k][j] << " " ;
00268     }//end for j
00269          //      if(hangingMask) std::cout << std::endl;
00270   }//end for k
00271   //  if(hangingMask) std::cout << " /* ********************************************************************** */" << std::endl;
00272   return true;
00273 }
00274 
00275 bool stiffnessMatrix::ElementalMatVec(int i, int j, int k, PetscScalar ***in, PetscScalar ***out, double scale){
00276   int dof= m_uiDof;
00277   int idx[8][3]={
00278     {k, j, dof*i},
00279     {k,j,dof*(i+1)},
00280     {k,j+1,dof*i},
00281     {k,j+1,dof*(i+1)},
00282     {k+1, j, dof*i},
00283     {k+1,j,dof*(i+1)},
00284     {k+1,j+1,dof*i},
00285     {k+1,j+1,dof*(i+1)}               
00286   };             
00287 
00288   double stencilScale;
00289   int **Ajk = (int **)m_stencil;
00290 
00291   PetscScalar ***nuarray = (PetscScalar ***)m_nuarray;
00292   
00293   for (int q = 0; q < 8; q++) {
00294     for (int r = 0; r < 8; r++) {
00295       stencilScale = -(nuarray[k][j][i]*m_dHx)*scale;
00296                 // stencilScale = -(m_nuval*m_dHx)*scale;
00297                 //              stencilScale = -(nuarray[k][j][dof*i]*m_dHx)*scale;
00298                 //stencilScale = -(m_dHx)*scale;
00299       out[idx[q][0]][idx[q][1]][idx[q][2]] +=  stencilScale*Ajk[q][r]*in[idx[r][0]][idx[r][1]][idx[r][2]];
00300       if (dof > 1)      
00301         out[idx[q][0]][idx[q][1]][idx[q][2]+1] += 0.0;
00302     }
00303   }
00304 #ifdef __DEBUG__
00305   //            PetscPrintf(0,"stencil scale in stiffness matrix = %f\n",stencilScale);
00306 #endif
00307 
00308   return true;
00309 }
00310 
00311 bool stiffnessMatrix::postMatVec() {
00312   if ( m_daType == PETSC) {
00313          PetscScalar ***nuarray = (PetscScalar ***)m_nuarray;
00314     int ierr = DAVecRestoreArray(m_DA, nuvec, &nuarray);
00315     CHKERRQ(ierr);
00316   } else {
00317     PetscScalar *nuarray = (PetscScalar *)m_nuarray;
00318     m_octDA->vecRestoreBuffer(nuvec, nuarray, false, false, true,m_uiDof);
00319   }
00320 
00321   return true;
00322 }
00323 
00324 bool stiffnessMatrix::ElementalMatGetDiagonal(int i, int j, int k, PetscScalar ***diag, double scale) {
00325   int dof= m_uiDof;
00326   int idx[8][3]={
00327     {k, j, dof*i},
00328     {k,j,dof*(i+1)},
00329     {k,j+1,dof*i},
00330     {k,j+1,dof*(i+1)},
00331     {k+1, j, dof*i},
00332     {k+1,j,dof*(i+1)},
00333     {k+1,j+1,dof*i},
00334     {k+1,j+1,dof*(i+1)}               
00335   };             
00336 
00337   double stencilScale;
00338   int **Ajk = (int **)m_stencil;
00339 
00340   PetscScalar ***nuarray = (PetscScalar ***)m_nuarray;
00341   
00342   for (int q = 0; q < 8; q++) {
00343         // stencilScale = -(m_nuval*m_dHx)*scale;
00344     stencilScale = -(nuarray[k][j][i]*m_dHx)*scale;
00345     diag [idx[q][0]][idx[q][1]][idx[q][2]] +=  stencilScale*Ajk[q][q] ;
00346         if (dof > 1)    
00347       diag[idx[q][0]][idx[q][1]][idx[q][2]+1] += 0.0;
00348   }
00349 
00350   // std::cout << "Nuval is " << nuarray[k][j][i] << std::endl;
00351   return true;
00352 }
00353 
00354 bool stiffnessMatrix::ElementalMatGetDiagonal(unsigned int i, PetscScalar *diag, double scale) {
00355   unsigned int lev = m_octDA->getLevel(i);
00356   double hx = xFac*(1u << (maxD - lev));
00357   double hy = yFac*(1u << (maxD - lev));
00358   double hz = zFac*(1u << (maxD - lev));
00359 
00360   double fac11 = -hx*scale/192.0;
00361   
00362   stdElemType elemType;
00363   unsigned int idx[8];
00364 
00365   unsigned char hangingMask = m_octDA->getHangingNodeIndex(i);    //  alignElementAndVertices(m_octDA, elemType, idx);       
00366   int ***Aijk = (int ***)m_stencil;
00367   
00368   alignElementAndVertices(m_octDA, elemType, idx);       
00369 
00370   PetscScalar *nuarray = (PetscScalar *)m_nuarray;
00371   for (int k = 0;k < 8;k++) {
00372     double fac1 = nuarray[idx[k]]*fac11;
00373     diag[m_uiDof*idx[k]] +=  (fac1*(Aijk[elemType][k][k]));
00374   }//end for k
00375 
00376   return true;
00377 }
00378 
00379 
00380 #endif /*_STIFFNESSMATRIX_H_*/
00381