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 7 static PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool); 8 static PetscErrorCode DMMoab_MatFillMatrixEntries_Private(DM,Mat); 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "DMCreateMatrix_Moab" 12 PetscErrorCode DMCreateMatrix_Moab(DM dm,Mat *J) 13 { 14 PetscErrorCode ierr; 15 ISLocalToGlobalMapping ltogb; 16 PetscInt innz,ionz,nlsiz; 17 DM_Moab *dmmoab=(DM_Moab*)dm->data; 18 PetscInt *nnz=0,*onz=0; 19 char *tmp=0; 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->bs); 30 31 /* allocate the nnz, onz arrays based on block size and local nodes */ 32 ierr = PetscMalloc((nlsiz)*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 33 ierr = PetscMemzero(nnz,sizeof(PetscInt)*(nlsiz));CHKERRQ(ierr); 34 ierr = PetscMalloc(nlsiz*sizeof(PetscInt),&onz);CHKERRQ(ierr); 35 ierr = PetscMemzero(onz,sizeof(PetscInt)*nlsiz);CHKERRQ(ierr); 36 37 /* compute the nonzero pattern based on MOAB connectivity data for local elements */ 38 ierr = DMMoab_Compute_NNZ_From_Connectivity(dm,&innz,nnz,&ionz,onz,(tmp?PETSC_TRUE:PETSC_FALSE));CHKERRQ(ierr); 39 40 /* create the Matrix and set its type as specified by user */ 41 ierr = MatCreate(dmmoab->pcomm->comm(), J);CHKERRQ(ierr); 42 ierr = MatSetSizes(*J, dmmoab->nloc*dmmoab->numFields, dmmoab->nloc*dmmoab->numFields, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 43 ierr = MatSetBlockSize(*J, dmmoab->bs);CHKERRQ(ierr); 44 ierr = MatSetType(*J, mtype);CHKERRQ(ierr); 45 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 46 47 if (!dmmoab->ltog_map) SETERRQ(dmmoab->pcomm->comm(), PETSC_ERR_ORDER, "Cannot create a DMMoab Mat without calling DMSetUp first."); 48 ierr = MatSetLocalToGlobalMapping(*J,dmmoab->ltog_map,dmmoab->ltog_map);CHKERRQ(ierr); 49 ierr = ISLocalToGlobalMappingBlock(dmmoab->ltog_map,dmmoab->bs,<ogb); 50 ierr = MatSetLocalToGlobalMappingBlock(*J,ltogb,ltogb);CHKERRQ(ierr); 51 ierr = ISLocalToGlobalMappingDestroy(<ogb);CHKERRQ(ierr); 52 53 /* set preallocation based on different supported Mat types */ 54 ierr = MatSeqAIJSetPreallocation(*J, innz, nnz);CHKERRQ(ierr); 55 ierr = MatMPIAIJSetPreallocation(*J, innz, nnz, ionz, onz);CHKERRQ(ierr); 56 ierr = MatSeqBAIJSetPreallocation(*J, dmmoab->bs, innz, nnz);CHKERRQ(ierr); 57 ierr = MatMPIBAIJSetPreallocation(*J, dmmoab->bs, innz, nnz, ionz, onz);CHKERRQ(ierr); 58 59 /* clean up temporary memory */ 60 ierr = PetscFree(nnz);CHKERRQ(ierr); 61 ierr = PetscFree(onz);CHKERRQ(ierr); 62 63 /* set up internal matrix data-structures */ 64 ierr = MatSetUp(*J);CHKERRQ(ierr); 65 66 /* set DM reference */ 67 ierr = MatSetDM(*J, dm);CHKERRQ(ierr); 68 69 /* set the correct NNZ pattern by setting matrix entries - make the matrix ready to use */ 70 ierr = DMMoab_MatFillMatrixEntries_Private(dm,*J);CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 75 #undef __FUNCT__ 76 #define __FUNCT__ "DMMoab_Compute_NNZ_From_Connectivity" 77 PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm,PetscInt* innz,PetscInt* nnz,PetscInt* ionz,PetscInt* onz,PetscBool isbaij) 78 { 79 PetscInt i,f,nloc,vpere,bs,ivtx,n_nnz,n_onz; 80 PetscInt ibs,jbs,inbsize,iobsize,nfields; 81 DM_Moab *dmmoab = (DM_Moab*)dm->data; 82 const moab::EntityHandle *connect; 83 moab::Range adjs,found,allvlocal,allvghost; 84 moab::Range::iterator iter,jter; 85 std::vector<moab::EntityHandle> storage; 86 PetscBool isinterlaced; 87 moab::EntityHandle vtx; 88 moab::ErrorCode merr; 89 90 PetscFunctionBegin; 91 bs = dmmoab->bs; 92 nloc = dmmoab->nloc; 93 nfields = dmmoab->numFields; 94 isinterlaced=(isbaij || bs==nfields ? PETSC_FALSE : PETSC_TRUE); 95 96 /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */ 97 merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,allvlocal,true);MBERRNM(merr); 98 merr = dmmoab->pcomm->filter_pstatus(allvlocal,PSTATUS_NOT_OWNED,PSTATUS_NOT,-1,&adjs);MBERRNM(merr); 99 allvghost = moab::subtract(allvlocal, adjs); 100 101 /* loop over the locally owned vertices and figure out the NNZ pattern using connectivity information */ 102 for(iter = dmmoab->vowned->begin(),ivtx=0; iter != dmmoab->vowned->end(); iter++,ivtx++) { 103 104 vtx = *iter; 105 adjs.clear(); 106 /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects 107 to the current vertex. We can then decipher if a vertex is ghosted or not and compute the 108 non-zero pattern accordingly. */ 109 merr = dmmoab->mbiface->get_adjacencies(&vtx,1,dmmoab->dim,false,adjs,moab::Interface::INTERSECT); 110 111 /* reset counters */ 112 n_nnz=n_onz=0; 113 found.clear(); 114 115 /* loop over vertices and update the number of connectivity */ 116 for(jter = adjs.begin(); jter != adjs.end(); jter++) { 117 118 /* Get connectivity information in canonical ordering for the local element */ 119 merr = dmmoab->mbiface->get_connectivity(*jter,connect,vpere,false,&storage);MBERRNM(merr); 120 121 /* loop over each element connected to the adjacent vertex and update as needed */ 122 for (i=0; i<vpere; ++i) { 123 if (connect[i] == vtx || found.find(connect[i]) != found.end()) continue; /* make sure we don't double count shared vertices */ 124 if (allvghost.find(connect[i]) != allvghost.end()) n_onz++; /* update out-of-proc onz */ 125 else n_nnz++; /* else local vertex */ 126 found.insert(connect[i]); 127 } 128 } 129 130 if (isbaij) { 131 nnz[ivtx]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 132 if (onz) { 133 onz[ivtx]=n_onz; /* add ghost non-owned nodes */ 134 } 135 } 136 else { /* AIJ matrices */ 137 if (!isinterlaced) { 138 for (f=0;f<nfields;f++) { 139 nnz[f*nloc+ivtx]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 140 if (onz) 141 onz[f*nloc+ivtx]=n_onz; /* add ghost non-owned nodes */ 142 } 143 } 144 else { 145 for (f=0;f<nfields;f++) { 146 nnz[nfields*ivtx+f]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 147 if (onz) 148 onz[nfields*ivtx+f]=n_onz; /* add ghost non-owned nodes */ 149 } 150 } 151 } 152 } 153 154 for (i=0;i<nloc*nfields;i++) 155 nnz[i]+=1; /* self count the node */ 156 157 for (i=0;i<nloc;i++) { 158 if (!isbaij) { 159 for (ibs=0; ibs<nfields; ibs++) { 160 /* first address the diagonal block */ 161 if (dmmoab->dfill) { 162 /* just add up the ints -- easier/faster rather than branching based on "1" */ 163 for (jbs=0,inbsize=0; jbs<nfields; jbs++) 164 inbsize += dmmoab->dfill[ibs*nfields+jbs]; 165 } 166 else inbsize=bs; /* dense coupling since user didn't specify the component fill explicitly */ 167 if (isinterlaced) nnz[i*nfields+ibs]*=inbsize; 168 else nnz[ibs*nloc+i]*=inbsize; 169 170 if (onz) { 171 /* next address the off-diagonal block */ 172 if (dmmoab->ofill) { 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=bs; /* dense coupling since user didn't specify the component fill explicitly */ 178 if (isinterlaced) onz[i*nfields+ibs]*=iobsize; 179 else onz[ibs*nloc+i]*=iobsize; 180 } 181 } 182 } 183 else { 184 /* check if we got overzealous in our nnz computations */ 185 nnz[i]=(nnz[i]>dmmoab->nloc ? dmmoab->nloc:nnz[i]); 186 } 187 } 188 /* update innz and ionz based on local maxima */ 189 if (innz || ionz) { 190 if (innz) *innz=0; 191 if (ionz) *ionz=0; 192 for (i=0;i<nloc*nfields;i++) { 193 if (innz && (nnz[i]>*innz)) *innz=nnz[i]; 194 if ((ionz && onz) && (onz[i]>*ionz)) *ionz=onz[i]; 195 } 196 } 197 PetscFunctionReturn(0); 198 } 199 200 201 #undef __FUNCT__ 202 #define __FUNCT__ "DMMoab_MatFillMatrixEntries_Private" 203 PetscErrorCode DMMoab_MatFillMatrixEntries_Private(DM dm, Mat A) 204 { 205 DM_Moab *dmmoab = (DM_Moab*)dm->data; 206 PetscInt nconn = 0,prev_nconn = 0; 207 const moab::EntityHandle *connect; 208 PetscScalar *locala=NULL; 209 PetscInt *dof_indices=NULL; 210 PetscErrorCode ierr; 211 212 PetscFunctionBegin; 213 /* loop over local elements */ 214 for(moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) { 215 const moab::EntityHandle ehandle = *iter; 216 217 /* Get connectivity information: */ 218 ierr = DMMoabGetElementConnectivity(dm, ehandle, &nconn, &connect);CHKERRQ(ierr); 219 220 /* if we have mixed elements or arrays have not been initialized - Allocate now */ 221 if (prev_nconn != nconn) { 222 if (locala) { 223 ierr = PetscFree(locala);CHKERRQ(ierr); 224 ierr = PetscFree(dof_indices);CHKERRQ(ierr); 225 } 226 ierr = PetscMalloc(sizeof(PetscScalar)*nconn*nconn*dmmoab->numFields*dmmoab->numFields,&locala);CHKERRQ(ierr); 227 ierr = PetscMemzero(locala,sizeof(PetscScalar)*nconn*nconn*dmmoab->numFields*dmmoab->numFields);CHKERRQ(ierr); 228 ierr = PetscMalloc(sizeof(PetscInt)*nconn,&dof_indices);CHKERRQ(ierr); 229 prev_nconn=nconn; 230 } 231 232 /* get the global DOF number to appropriately set the element contribution in the RHS vector */ 233 ierr = DMMoabGetDofsBlockedLocal(dm, nconn, connect, dof_indices);CHKERRQ(ierr); 234 235 /* set the values directly into appropriate locations. Can alternately use VecSetValues */ 236 ierr = MatSetValuesBlockedLocal(A, nconn, dof_indices, nconn, dof_indices, locala, INSERT_VALUES);CHKERRQ(ierr); 237 } 238 239 /* clean up memory */ 240 ierr = PetscFree(locala);CHKERRQ(ierr); 241 ierr = PetscFree(dof_indices);CHKERRQ(ierr); 242 243 /* finish assembly */ 244 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 245 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 246 PetscFunctionReturn(0); 247 } 248 249 250 #undef __FUNCT__ 251 #define __FUNCT__ "DMMoabSetBlockFills_Private" 252 static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w,const PetscInt *fill,PetscInt **rfill) 253 { 254 PetscErrorCode ierr; 255 PetscInt i,j,*ifill; 256 257 PetscFunctionBegin; 258 if (!fill) PetscFunctionReturn(0); 259 ierr = PetscMalloc1(w*w,&ifill);CHKERRQ(ierr); 260 for (i=0;i<w;i++) { 261 for (j=0; j<w; j++) 262 ifill[i*w+j]=fill[i*w+j]; 263 } 264 265 *rfill = ifill; 266 PetscFunctionReturn(0); 267 } 268 269 270 #undef __FUNCT__ 271 #define __FUNCT__ "DMMoabSetBlockFills" 272 /*@ 273 DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem 274 of the matrix returned by DMCreateMatrix(). 275 276 Logically Collective on DMDA 277 278 Input Parameter: 279 + dm - the DMMoab object 280 . dfill - the fill pattern in the diagonal block (may be NULL, means use dense block) 281 - ofill - the fill pattern in the off-diagonal blocks 282 283 Level: developer 284 285 Notes: This only makes sense when you are doing multicomponent problems but using the 286 MPIAIJ matrix format 287 288 The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries 289 representing coupling and 0 entries for missing coupling. For example 290 $ dfill[9] = {1, 0, 0, 291 $ 1, 1, 0, 292 $ 0, 1, 1} 293 means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with 294 itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the 295 diagonal block). 296 297 DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then 298 can be represented in the dfill, ofill format 299 300 Contributed by Glenn Hammond 301 302 .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly() 303 304 @*/ 305 PetscErrorCode DMMoabSetBlockFills(DM dm,const PetscInt *dfill,const PetscInt *ofill) 306 { 307 DM_Moab *dmmoab=(DM_Moab*)dm->data; 308 PetscErrorCode ierr; 309 310 PetscFunctionBegin; 311 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 312 ierr = DMMoabSetBlockFills_Private(dmmoab->numFields,dfill,&dmmoab->dfill);CHKERRQ(ierr); 313 ierr = DMMoabSetBlockFills_Private(dmmoab->numFields,ofill,&dmmoab->ofill);CHKERRQ(ierr); 314 PetscFunctionReturn(0); 315 } 316