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