1 2 #include "private/daimpl.h" /*I "petscdm.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMGetElements_DA_1D" 6 static PetscErrorCode DMGetElements_DA_1D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 7 { 8 PetscErrorCode ierr; 9 DM_DA *da = (DM_DA*)dm->data; 10 PetscInt i,xs,xe,Xs,Xe; 11 PetscInt cnt=0; 12 PetscFunctionBegin; 13 if (!da->e) { 14 ierr = DMDAGetCorners(dm,&xs,0,0,&xe,0,0);CHKERRQ(ierr); 15 ierr = DMDAGetGhostCorners(dm,&Xs,0,0,&Xe,0,0);CHKERRQ(ierr); 16 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 17 da->ne = 1*(xe - xs - 1); 18 ierr = PetscMalloc((1 + 2*da->ne)*sizeof(PetscInt),&da->e);CHKERRQ(ierr); 19 for (i=xs; i<xe-1; i++) { 20 da->e[cnt++] = (i-Xs ); 21 da->e[cnt++] = (i-Xs+1); 22 } 23 } 24 *nel = da->ne; 25 *nen = 2; 26 *e = da->e; 27 PetscFunctionReturn(0); 28 } 29 30 #undef __FUNCT__ 31 #define __FUNCT__ "DMGetElements_DA_2D" 32 static PetscErrorCode DMGetElements_DA_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 33 { 34 PetscErrorCode ierr; 35 DM_DA *da = (DM_DA*)dm->data; 36 PetscInt i,xs,xe,Xs,Xe; 37 PetscInt j,ys,ye,Ys,Ye; 38 PetscInt cnt=0, cell[4], ns=2, nn=3; 39 PetscInt c, split[] = {0,1,3, 40 2,3,1}; 41 PetscFunctionBegin; 42 if (!da->e) { 43 if (da->elementtype == DMDA_ELEMENT_P1) {ns=2; nn=3;} 44 if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=4;} 45 ierr = DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);CHKERRQ(ierr); 46 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);CHKERRQ(ierr); 47 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 48 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 49 da->ne = ns*(xe - xs - 1)*(ye - ys - 1); 50 ierr = PetscMalloc((1 + nn*da->ne)*sizeof(PetscInt),&da->e);CHKERRQ(ierr); 51 for (j=ys; j<ye-1; j++) { 52 for (i=xs; i<xe-1; i++) { 53 cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs); 54 cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs); 55 cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs); 56 cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs); 57 if (da->elementtype == DMDA_ELEMENT_P1) { 58 for (c=0; c<ns*nn; c++) 59 da->e[cnt++] = cell[split[c]]; 60 } 61 if (da->elementtype == DMDA_ELEMENT_Q1) { 62 for (c=0; c<ns*nn; c++) 63 da->e[cnt++] = cell[c]; 64 } 65 } 66 } 67 } 68 *nel = da->ne; 69 *nen = nn; 70 *e = da->e; 71 PetscFunctionReturn(0); 72 } 73 74 #undef __FUNCT__ 75 #define __FUNCT__ "DMGetElements_DA_3D" 76 static PetscErrorCode DMGetElements_DA_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 77 { 78 PetscErrorCode ierr; 79 DM_DA *da = (DM_DA*)dm->data; 80 PetscInt i,xs,xe,Xs,Xe; 81 PetscInt j,ys,ye,Ys,Ye; 82 PetscInt k,zs,ze,Zs,Ze; 83 PetscInt cnt=0, cell[8], ns=6, nn=4; 84 PetscInt c, split[] = {0,1,3,7, 85 0,1,7,4, 86 1,2,3,7, 87 1,2,7,6, 88 1,4,5,7, 89 1,5,6,7}; 90 PetscFunctionBegin; 91 if (!da->e) { 92 if (da->elementtype == DMDA_ELEMENT_P1) {ns=6; nn=4;} 93 if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=8;} 94 ierr = DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);CHKERRQ(ierr); 95 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);CHKERRQ(ierr); 96 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 97 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 98 ze += zs; Ze += Zs; if (zs != Zs) zs -= 1; 99 da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1); 100 ierr = PetscMalloc((1 + nn*da->ne)*sizeof(PetscInt),&da->e);CHKERRQ(ierr); 101 for (k=zs; k<ze-1; k++) { 102 for (j=ys; j<ye-1; j++) { 103 for (i=xs; i<xe-1; i++) { 104 cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 105 cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 106 cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 107 cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 108 cell[4] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 109 cell[5] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 110 cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 111 cell[7] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 112 if (da->elementtype == DMDA_ELEMENT_P1) { 113 for (c=0; c<ns*nn; c++) 114 da->e[cnt++] = cell[split[c]]; 115 } 116 if (da->elementtype == DMDA_ELEMENT_Q1) { 117 for (c=0; c<ns*nn; c++) 118 da->e[cnt++] = cell[c]; 119 } 120 } 121 } 122 } 123 } 124 *nel = da->ne; 125 *nen = nn; 126 *e = da->e; 127 PetscFunctionReturn(0); 128 } 129 130 #undef __FUNCT__ 131 #define __FUNCT__ "DMGetElements_DA" 132 PetscErrorCode DMGetElements_DA(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 133 { 134 DM_DA *da = (DM_DA*)dm->data; 135 PetscErrorCode ierr; 136 PetscFunctionBegin; 137 if (da->dim==-1) { 138 *nel = 0; *nen = 0; *e = PETSC_NULL; 139 } else if (da->dim==1) { 140 ierr = DMGetElements_DA_1D(dm,nel,nen,e);CHKERRQ(ierr); 141 } else if (da->dim==2) { 142 ierr = DMGetElements_DA_2D(dm,nel,nen,e);CHKERRQ(ierr); 143 } else if (da->dim==3) { 144 ierr = DMGetElements_DA_3D(dm,nel,nen,e);CHKERRQ(ierr); 145 } else { 146 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",da->dim); 147 } 148 149 PetscFunctionReturn(0); 150 } 151