1 #include <petsc-private/dmmbimpl.h> /*I "petscdm.h" I*/ 2 #include <petsc-private/vecimpl.h> /*I "petscdm.h" I*/ 3 4 #include <petscdmmoab.h> 5 6 static PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool); 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "DMCreateMatrix_Moab" 10 PetscErrorCode DMCreateMatrix_Moab(DM dm, MatType mtype,Mat *J) 11 { 12 PetscErrorCode ierr; 13 ISLocalToGlobalMapping ltog; 14 MatType mltype; 15 PetscInt i,nloc,count,dof,innz,ionz,nsize; 16 DM_Moab *dmmoab = (DM_Moab*)dm->data; 17 moab::Range *range=dmmoab->vowned; 18 PetscInt *nnz,*onz; 19 char *tmp=0; 20 moab::ErrorCode merr; 21 22 PetscFunctionBegin; 23 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 24 PetscValidPointer(J,3); 25 nloc = dmmoab->vowned->size() * dmmoab->bs; 26 27 /* create the Matrix and set its type as specified by user */ 28 ierr = MatCreate(dmmoab->pcomm->comm(), J);CHKERRQ(ierr); 29 ierr = MatSetSizes(*J, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 30 ierr = MatSetType(*J, mtype);CHKERRQ(ierr); 31 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 32 33 /* next, need to allocate the non-zero arrays to enable pre-allocation */ 34 ierr = MatGetType(*J, &mltype);CHKERRQ(ierr); /* in case user overrode the default type from command-line, re-check the type */ 35 ierr = PetscStrstr(mltype, "baij", &tmp);CHKERRQ(ierr); 36 nsize = (tmp ? dmmoab->nloc:dmmoab->nloc*dmmoab->bs); 37 38 /* allocate the nnz, onz arrays based on block size and local nodes */ 39 ierr = PetscMalloc(nsize*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 40 ierr = PetscMemzero(nnz,sizeof(PetscInt)*nsize);CHKERRQ(ierr); 41 ierr = PetscMalloc(nsize*sizeof(PetscInt),&onz);CHKERRQ(ierr); 42 ierr = PetscMemzero(onz,sizeof(PetscInt)*nsize);CHKERRQ(ierr); 43 44 /* compute the nonzero pattern based on MOAB connectivity data for local elements */ 45 ierr = DMMoab_Compute_NNZ_From_Connectivity(dm,&innz,nnz,&ionz,onz,(tmp?PETSC_TRUE:PETSC_FALSE));CHKERRQ(ierr); 46 47 if (dmmoab->bs > 1 && tmp) { 48 /* Block matrix created, now set local to global mapping */ 49 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,range->size(),dmmoab->gsindices,PETSC_COPY_VALUES,<og);CHKERRQ(ierr); 50 ierr = MatSetLocalToGlobalMappingBlock(*J,ltog,ltog);CHKERRQ(ierr); 51 52 /* Clean up */ 53 ierr = ISLocalToGlobalMappingDestroy(<og);CHKERRQ(ierr); 54 55 } else { 56 ierr = MatSetBlockSize(*J, dmmoab->bs);CHKERRQ(ierr); 57 } 58 59 /* set preallocation based on different supported Mat types */ 60 ierr = MatSeqAIJSetPreallocation(*J, innz, nnz);CHKERRQ(ierr); 61 ierr = MatMPIAIJSetPreallocation(*J, innz, nnz, ionz, onz);CHKERRQ(ierr); 62 ierr = MatSeqBAIJSetPreallocation(*J, dmmoab->bs, innz, nnz);CHKERRQ(ierr); 63 ierr = MatMPIBAIJSetPreallocation(*J, dmmoab->bs, innz, nnz, ionz, onz);CHKERRQ(ierr); 64 65 ierr = MatSetDM(*J, dm);CHKERRQ(ierr); 66 ierr = PetscFree(nnz);CHKERRQ(ierr); 67 ierr = PetscFree(onz);CHKERRQ(ierr); 68 PetscFunctionReturn(0); 69 } 70 71 72 #undef __FUNCT__ 73 #define __FUNCT__ "DMMoab_Compute_NNZ_From_Connectivity" 74 PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm,PetscInt* innz,PetscInt* nnz,PetscInt* ionz,PetscInt* onz,PetscBool isbaij) 75 { 76 PetscErrorCode ierr; 77 PetscInt i,j,f,nloc,vpere,bs,nsize,nghost_found; 78 DM_Moab *dmmoab = (DM_Moab*)dm->data; 79 const moab::EntityHandle *connect; 80 moab::Range *vowned=dmmoab->vowned,*elocal=dmmoab->elocal,*eghost=dmmoab->eghost,adjs,visited; 81 moab::Range::iterator iter; 82 PetscInt dof,doff,ndofs; 83 moab::Tag id_tag; 84 moab::ErrorCode merr; 85 PetscSection section; 86 87 PetscFunctionBegin; 88 bs = dmmoab->bs; 89 nloc = dmmoab->nloc; 90 nsize = (isbaij ? nloc:nloc*bs); 91 92 ierr = DMMoabGetLocalToGlobalTag(dm,&id_tag);CHKERRQ(ierr); 93 94 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 95 96 for(iter = elocal->begin(),i=0; iter != elocal->end(); iter++,i++) { 97 /* Get connectivity information in canonical ordering for the local element */ 98 merr = dmmoab->mbiface->get_connectivity((*iter),connect,vpere);MBERRNM(merr); 99 100 nghost_found=0; 101 /* loop over vertices and update the number of connectivity */ 102 for (j=0;j<vpere;j++) { 103 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 104 if (giter != dmmoab->vghost->end()) nghost_found++; 105 } 106 107 /* loop over vertices and update the number of connectivity */ 108 for (j=0;j<vpere;j++) { 109 110 /* 111 moab::Range::const_iterator giter = visited.find(connect[j]); 112 if (giter != visited.end()) continue; 113 114 visited.insert(connect[j]); 115 */ 116 117 /* loop over vertices and update the number of connectivity */ 118 for(jter = adjs.begin(); jter != adjs.end(); jter++) { 119 120 /* Get connectivity information in canonical ordering for the local element */ 121 merr = dmmoab->mbiface->get_connectivity(*jter,connect,vpere,false,&storage);MBERRNM(merr); 122 123 for (i=0; i<vpere; ++i) { 124 if (connect[i] == vtx || found.find(connect[i]) != found.end()) continue; 125 if (dmmoab->vghost->find(connect[i]) != dmmoab->vghost->end()) n_onz++; 126 else n_nnz++; 127 found.insert(connect[i]); 128 } 129 } 130 131 merr = dmmoab->mbiface->tag_get_data(id_tag,&vtx,1,&gid);MBERRNM(merr); 132 133 if (isbaij) { 134 nnz[ivtx]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 135 if (onz) onz[ivtx]=n_onz; /* add ghost non-owned nodes */ 136 } 137 else { 138 for (f=0;f<dmmoab->nfields;f++) { 139 nnz[dmmoab->nfields*ivtx+f]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 140 if (onz) onz[dmmoab->nfields*ivtx+f]=n_onz; /* add ghost non-owned nodes */ 141 } 142 } 143 } 144 145 visited.clear(); 146 147 if (innz) *innz=0; 148 if (ionz) *ionz=0; 149 for (i=0;i<nsize;i++) { 150 nnz[i]+=1; /* self count the node */ 151 if (!isbaij) { 152 nnz[i]*=bs; 153 onz[i]*=bs; 154 } 155 156 if (innz && nnz[i]>*innz) *innz=nnz[i]; 157 if (ionz && onz[i]>*ionz) *ionz=onz[i]; 158 } 159 PetscFunctionReturn(0); 160 } 161 162