1 #include <petsc/private/dmmbimpl.h> /*I "petscdmmoab.h" I*/ 2 #include <petsc/private/vecimpl.h> 3 4 #include <petscdmmoab.h> 5 #include <MBTagConventions.hpp> 6 #include <moab/NestedRefine.hpp> 7 8 PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM, PetscInt*, PetscInt*, PetscInt*, PetscInt*, PetscBool); 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "DMCreateMatrix_Moab" 12 PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J) 13 { 14 PetscErrorCode ierr; 15 PetscInt innz = 0, ionz = 0, nlsiz; 16 DM_Moab *dmmoab = (DM_Moab*)dm->data; 17 PetscInt *nnz = 0, *onz = 0; 18 char *tmp = 0; 19 Mat A; 20 MatType mtype; 21 22 PetscFunctionBegin; 23 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 24 PetscValidPointer(J, 3); 25 26 /* next, need to allocate the non-zero arrays to enable pre-allocation */ 27 mtype = dm->mattype; 28 ierr = PetscStrstr(mtype, "baij", &tmp);CHKERRQ(ierr); 29 nlsiz = (tmp ? dmmoab->nloc : dmmoab->nloc * dmmoab->numFields); 30 31 /* allocate the nnz, onz arrays based on block size and local nodes */ 32 ierr = PetscCalloc2(nlsiz, &nnz, nlsiz, &onz);CHKERRQ(ierr); 33 34 /* compute the nonzero pattern based on MOAB connectivity data for local elements */ 35 ierr = DMMoab_Compute_NNZ_From_Connectivity(dm, &innz, nnz, &ionz, onz, (tmp ? PETSC_TRUE : PETSC_FALSE));CHKERRQ(ierr); 36 37 /* create the Matrix and set its type as specified by user */ 38 ierr = MatCreate((((PetscObject)dm)->comm), &A);CHKERRQ(ierr); 39 ierr = MatSetSizes(A, dmmoab->nloc * dmmoab->numFields, dmmoab->nloc * dmmoab->numFields, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 40 ierr = MatSetType(A, mtype);CHKERRQ(ierr); 41 ierr = MatSetBlockSize(A, dmmoab->bs);CHKERRQ(ierr); 42 ierr = MatSetDM(A, dm);CHKERRQ(ierr); /* set DM reference */ 43 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 44 45 if (!dmmoab->ltog_map) SETERRQ((((PetscObject)dm)->comm), PETSC_ERR_ORDER, "Cannot create a DMMoab Mat without calling DMSetUp first."); 46 ierr = MatSetLocalToGlobalMapping(A, dmmoab->ltog_map, dmmoab->ltog_map);CHKERRQ(ierr); 47 48 /* set preallocation based on different supported Mat types */ 49 ierr = MatSeqAIJSetPreallocation(A, innz, nnz);CHKERRQ(ierr); 50 ierr = MatMPIAIJSetPreallocation(A, innz, nnz, ionz, onz);CHKERRQ(ierr); 51 ierr = MatSeqBAIJSetPreallocation(A, dmmoab->bs, innz, nnz);CHKERRQ(ierr); 52 ierr = MatMPIBAIJSetPreallocation(A, dmmoab->bs, innz, nnz, ionz, onz);CHKERRQ(ierr); 53 54 /* clean up temporary memory */ 55 ierr = PetscFree2(nnz, onz);CHKERRQ(ierr); 56 57 /* set up internal matrix data-structures */ 58 ierr = MatSetUp(A);CHKERRQ(ierr); 59 60 MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); 61 62 *J = A; 63 PetscFunctionReturn(0); 64 } 65 66 67 #undef __FUNCT__ 68 #define __FUNCT__ "DMMoab_Compute_NNZ_From_Connectivity" 69 PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm, PetscInt* innz, PetscInt* nnz, PetscInt* ionz, PetscInt* onz, PetscBool isbaij) 70 { 71 PetscInt i, f, nloc, vpere, bs, n_nnz, n_onz, ivtx = 0; 72 PetscInt ibs, jbs, inbsize, iobsize, nfields, nlsiz; 73 DM_Moab *dmmoab = (DM_Moab*)dm->data; 74 moab::Range found; 75 std::vector<moab::EntityHandle> adjs, storage; 76 PetscBool isinterlaced; 77 moab::EntityHandle vtx; 78 moab::ErrorCode merr; 79 80 PetscFunctionBegin; 81 bs = dmmoab->bs; 82 nloc = dmmoab->nloc; 83 nfields = dmmoab->numFields; 84 isinterlaced = (isbaij || bs == nfields ? PETSC_TRUE : PETSC_FALSE); 85 nlsiz = (isinterlaced ? nloc : nloc * nfields); 86 87 /* loop over the locally owned vertices and figure out the NNZ pattern using connectivity information */ 88 for (moab::Range::const_iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, ivtx++) { 89 90 vtx = *iter; 91 /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects 92 to the current vertex. We can then decipher if a vertex is ghosted or not and compute the 93 non-zero pattern accordingly. */ 94 adjs.clear(); 95 if (dmmoab->hlevel) { 96 merr = dmmoab->hierarchy->get_adjacencies(vtx, dmmoab->dim, adjs); MBERRNM(merr); 97 } 98 else { 99 merr = dmmoab->mbiface->get_adjacencies(&vtx, 1, dmmoab->dim, true, adjs, moab::Interface::INTERSECT); MBERRNM(merr); 100 } 101 102 /* reset counters */ 103 n_nnz = n_onz = 0; 104 found.clear(); 105 106 /* loop over vertices and update the number of connectivity */ 107 for (unsigned jter = 0; jter < adjs.size(); ++jter) { 108 109 /* Get connectivity information in canonical ordering for the local element */ 110 const moab::EntityHandle *connect; 111 std::vector<moab::EntityHandle> cconnect; 112 if (dmmoab->hlevel) { 113 merr = dmmoab->hierarchy->get_connectivity(adjs[jter], dmmoab->hlevel, cconnect); MBERRNM(merr); 114 connect=cconnect.data(); 115 vpere=cconnect.size(); 116 } 117 else { 118 merr = dmmoab->mbiface->get_connectivity(adjs[jter], connect, vpere, false, &storage); MBERRNM(merr); 119 } 120 121 /* loop over each element connected to the adjacent vertex and update as needed */ 122 for (i = 0; i < vpere; ++i) { 123 /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */ 124 if (connect[i] == vtx || found.find(connect[i]) != found.end()) continue; /* make sure we don't double count shared vertices */ 125 if (dmmoab->vghost->find(connect[i]) != dmmoab->vghost->end()) n_onz++; /* update out-of-proc onz */ 126 else n_nnz++; /* else local vertex */ 127 found.insert(connect[i]); 128 } 129 } 130 storage.clear(); 131 132 if (isbaij) { 133 nnz[ivtx] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 134 if (onz) { 135 onz[ivtx] = n_onz; /* add ghost non-owned nodes */ 136 } 137 } 138 else { /* AIJ matrices */ 139 if (!isinterlaced) { 140 for (f = 0; f < nfields; f++) { 141 nnz[f * nloc + ivtx] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 142 if (onz) 143 onz[f * nloc + ivtx] = n_onz; /* add ghost non-owned nodes */ 144 } 145 } 146 else { 147 for (f = 0; f < nfields; f++) { 148 nnz[nfields * ivtx + f] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 149 if (onz) 150 onz[nfields * ivtx + f] = n_onz; /* add ghost non-owned nodes */ 151 } 152 } 153 } 154 } 155 156 for (i = 0; i < nlsiz; i++) 157 nnz[i] += 1; /* self count the node */ 158 159 for (ivtx = 0; ivtx < nloc; ivtx++) { 160 if (!isbaij) { 161 for (ibs = 0; ibs < nfields; ibs++) { 162 if (dmmoab->dfill) { /* first address the diagonal block */ 163 /* just add up the ints -- easier/faster rather than branching based on "1" */ 164 for (jbs = 0, inbsize = 0; jbs < nfields; jbs++) 165 inbsize += dmmoab->dfill[ibs * nfields + jbs]; 166 } 167 else inbsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */ 168 if (isinterlaced) nnz[ivtx * nfields + ibs] *= inbsize; 169 else nnz[ibs * nloc + ivtx] *= inbsize; 170 171 if (onz) { 172 if (dmmoab->ofill) { /* next address the off-diagonal block */ 173 /* just add up the ints -- easier/faster rather than branching based on "1" */ 174 for (jbs = 0, iobsize = 0; jbs < nfields; jbs++) 175 iobsize += dmmoab->dfill[ibs * nfields + jbs]; 176 } 177 else iobsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */ 178 if (isinterlaced) onz[ivtx * nfields + ibs] *= iobsize; 179 else onz[ibs * nloc + ivtx] *= iobsize; 180 } 181 } 182 } 183 else { 184 /* check if we got overzealous in our nnz and onz computations */ 185 nnz[ivtx] = (nnz[ivtx] > dmmoab->nloc ? dmmoab->nloc : nnz[ivtx]); 186 if (onz) onz[ivtx] = (onz[ivtx] > dmmoab->nloc ? dmmoab->nloc : onz[ivtx]); 187 } 188 } 189 /* update innz and ionz based on local maxima */ 190 if (innz || ionz) { 191 if (innz) *innz = 0; 192 if (ionz) *ionz = 0; 193 if (isbaij) { 194 for (i = 0; i < nlsiz; i++) { 195 if (innz && (nnz[i] > *innz)) *innz = nnz[i]; 196 if ((ionz && onz) && (onz[i] > *ionz)) *ionz = onz[i]; 197 } 198 } 199 else { 200 for (i = 0; i < nlsiz * nfields; i++) { 201 if (innz && (nnz[i] > *innz)) *innz = nnz[i]; 202 if ((ionz && onz) && (onz[i] > *ionz)) *ionz = onz[i]; 203 } 204 } 205 } 206 PetscFunctionReturn(0); 207 } 208 209 210 #undef __FUNCT__ 211 #define __FUNCT__ "DMMoabSetBlockFills_Private" 212 static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w, const PetscInt *fill, PetscInt **rfill) 213 { 214 PetscErrorCode ierr; 215 PetscInt i, j, *ifill; 216 217 PetscFunctionBegin; 218 if (!fill) PetscFunctionReturn(0); 219 ierr = PetscMalloc1(w * w, &ifill);CHKERRQ(ierr); 220 for (i = 0; i < w; i++) { 221 for (j = 0; j < w; j++) 222 ifill[i * w + j] = fill[i * w + j]; 223 } 224 225 *rfill = ifill; 226 PetscFunctionReturn(0); 227 } 228 229 230 /*@ 231 DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem 232 of the matrix returned by DMCreateMatrix(). 233 234 Logically Collective on DMDA 235 236 Input Parameter: 237 + dm - the DMMoab object 238 . dfill - the fill pattern in the diagonal block (may be NULL, means use dense block) 239 - ofill - the fill pattern in the off-diagonal blocks 240 241 Level: developer 242 243 Notes: This only makes sense when you are doing multicomponent problems but using the 244 MPIAIJ matrix format 245 246 The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries 247 representing coupling and 0 entries for missing coupling. For example 248 $ dfill[9] = {1, 0, 0, 249 $ 1, 1, 0, 250 $ 0, 1, 1} 251 means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with 252 itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the 253 diagonal block). 254 255 DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then 256 can be represented in the dfill, ofill format 257 258 Contributed by Glenn Hammond 259 260 .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly() 261 262 @*/ 263 PetscErrorCode DMMoabSetBlockFills(DM dm, const PetscInt *dfill, const PetscInt *ofill) 264 { 265 DM_Moab *dmmoab = (DM_Moab*)dm->data; 266 PetscErrorCode ierr; 267 268 PetscFunctionBegin; 269 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 270 ierr = DMMoabSetBlockFills_Private(dmmoab->numFields, dfill, &dmmoab->dfill);CHKERRQ(ierr); 271 ierr = DMMoabSetBlockFills_Private(dmmoab->numFields, ofill, &dmmoab->ofill);CHKERRQ(ierr); 272 PetscFunctionReturn(0); 273 } 274