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 #include <MBTagConventions.hpp> 6 7 static PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool); 8 9 #undef __FUNCT__ 10 #define __FUNCT__ "DMCreateMatrix_Moab" 11 PetscErrorCode DMCreateMatrix_Moab(DM dm, MatType mtype,Mat *J) 12 { 13 PetscErrorCode ierr; 14 ISLocalToGlobalMapping ltogb; 15 PetscInt innz,ionz,nlsiz; 16 DM_Moab *dmmoab=(DM_Moab*)dm->data; 17 PetscInt *nnz=0,*onz=0; 18 char *tmp=0; 19 moab::ErrorCode merr; 20 21 PetscFunctionBegin; 22 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 23 PetscValidPointer(J,3); 24 25 /* create the Matrix and set its type as specified by user */ 26 ierr = MatCreate(dmmoab->pcomm->comm(), J);CHKERRQ(ierr); 27 ierr = MatSetSizes(*J, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 28 ierr = MatSetType(*J, mtype);CHKERRQ(ierr); 29 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 30 31 /* next, need to allocate the non-zero arrays to enable pre-allocation */ 32 ierr = MatGetType(*J, &mltype);CHKERRQ(ierr); /* in case user overrode the default type from command-line, re-check the type */ 33 ierr = PetscStrstr(mltype, "baij", &tmp);CHKERRQ(ierr); 34 nsize = (tmp ? dmmoab->nloc:dmmoab->nloc*dmmoab->bs); 35 36 /* allocate the nnz, onz arrays based on block size and local nodes */ 37 ierr = PetscMalloc((nlsiz)*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 38 ierr = PetscMemzero(nnz,sizeof(PetscInt)*(nlsiz));CHKERRQ(ierr); 39 ierr = PetscMalloc(nlsiz*sizeof(PetscInt),&onz);CHKERRQ(ierr); 40 ierr = PetscMemzero(onz,sizeof(PetscInt)*nlsiz);CHKERRQ(ierr); 41 42 /* compute the nonzero pattern based on MOAB connectivity data for local elements */ 43 ierr = DMMoab_Compute_NNZ_From_Connectivity(dm,&innz,nnz,&ionz,onz,(tmp?PETSC_TRUE:PETSC_FALSE));CHKERRQ(ierr); 44 45 /* create the Matrix and set its type as specified by user */ 46 ierr = MatCreate(dmmoab->pcomm->comm(), J);CHKERRQ(ierr); 47 ierr = MatSetSizes(*J, dmmoab->nloc*dmmoab->nfields, dmmoab->nloc*dmmoab->nfields, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 48 ierr = MatSetBlockSize(*J, dmmoab->bs);CHKERRQ(ierr); 49 ierr = MatSetType(*J, mtype);CHKERRQ(ierr); 50 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 51 52 if (!dmmoab->ltog_map) SETERRQ(dmmoab->pcomm->comm(), PETSC_ERR_ORDER, "Cannot create a DMMoab Mat without calling DMSetUp first."); 53 ierr = MatSetLocalToGlobalMapping(*J,dmmoab->ltog_map,dmmoab->ltog_map);CHKERRQ(ierr); 54 ierr = ISLocalToGlobalMappingBlock(dmmoab->ltog_map,dmmoab->bs,<ogb); 55 ierr = MatSetLocalToGlobalMappingBlock(*J,ltogb,ltogb);CHKERRQ(ierr); 56 ierr = ISLocalToGlobalMappingDestroy(<ogb);CHKERRQ(ierr); 57 58 /* set preallocation based on different supported Mat types */ 59 ierr = MatSeqAIJSetPreallocation(*J, innz, nnz);CHKERRQ(ierr); 60 ierr = MatMPIAIJSetPreallocation(*J, innz, nnz, ionz, onz);CHKERRQ(ierr); 61 ierr = MatSeqBAIJSetPreallocation(*J, dmmoab->bs, innz, nnz);CHKERRQ(ierr); 62 ierr = MatMPIBAIJSetPreallocation(*J, dmmoab->bs, innz, nnz, ionz, onz);CHKERRQ(ierr); 63 64 ierr = MatSetDM(*J, dm);CHKERRQ(ierr); 65 ierr = PetscFree(nnz);CHKERRQ(ierr); 66 ierr = PetscFree(onz);CHKERRQ(ierr); 67 PetscFunctionReturn(0); 68 } 69 70 71 #undef __FUNCT__ 72 #define __FUNCT__ "DMMoab_Compute_NNZ_From_Connectivity" 73 PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm,PetscInt* innz,PetscInt* nnz,PetscInt* ionz,PetscInt* onz,PetscBool isbaij) 74 { 75 PetscInt i,f,nloc,vpere,bs,nsize,ivtx,n_nnz,n_onz; 76 DM_Moab *dmmoab = (DM_Moab*)dm->data; 77 const moab::EntityHandle *connect; 78 moab::Range adjs,found,allvlocal,allvghost; 79 moab::Range::iterator iter,jter; 80 std::vector<moab::EntityHandle> storage; 81 moab::EntityHandle vtx; 82 moab::ErrorCode merr; 83 84 PetscFunctionBegin; 85 bs = dmmoab->bs; 86 nloc = dmmoab->nloc; 87 nsize = (isbaij ? nloc:nloc*bs); 88 89 /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */ 90 merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,allvlocal,true);MBERRNM(merr); 91 merr = dmmoab->pcomm->filter_pstatus(allvlocal,PSTATUS_NOT_OWNED,PSTATUS_NOT,-1,&adjs);MBERRNM(merr); 92 allvghost = moab::subtract(allvlocal, adjs); 93 94 /* loop over vertices and update the number of connectivity */ 95 for (j=0;j<vpere;j++) { 96 97 vtx = *iter; 98 adjs.clear(); 99 /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects 100 to the current vertex. We can then decipher if a vertex is ghosted or not and compute the 101 non-zero pattern accordingly. */ 102 merr = dmmoab->mbiface->get_adjacencies(&vtx,1,dmmoab->dim,false,adjs,moab::Interface::INTERSECT); 103 104 /* loop over vertices and update the number of connectivity */ 105 for(jter = adjs.begin(); jter != adjs.end(); jter++) { 106 107 /* Get connectivity information in canonical ordering for the local element */ 108 merr = dmmoab->mbiface->get_connectivity(*jter,connect,vpere,false,&storage);MBERRNM(merr); 109 110 for (i=0; i<vpere; ++i) { 111 if (connect[i] == vtx || found.find(connect[i]) != found.end()) continue; 112 if (allvghost.find(connect[i]) != allvghost.end()) n_onz++; 113 else n_nnz++; 114 found.insert(connect[i]); 115 } 116 } 117 118 if (isbaij) { 119 nnz[ivtx]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 120 if (onz) onz[ivtx]=n_onz; /* add ghost non-owned nodes */ 121 } 122 else { 123 for (f=0;f<dmmoab->nfields;f++) { 124 nnz[dmmoab->nfields*ivtx+f]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 125 if (onz) onz[dmmoab->nfields*ivtx+f]=n_onz; /* add ghost non-owned nodes */ 126 } 127 } 128 } 129 130 if (innz) *innz=0; 131 if (ionz) *ionz=0; 132 for (i=0;i<nsize;i++) { 133 nnz[i]+=1; /* self count the node */ 134 /* check if we got overzealous */ 135 nnz[i]=(nnz[i]>dmmoab->nloc ? dmmoab->nloc:nnz[i]); 136 if (!isbaij) { 137 nnz[i]*=bs; 138 onz[i]*=bs; 139 } 140 141 if (innz && (nnz[i]>*innz)) *innz=nnz[i]; 142 if ((ionz && onz) && (onz[i]>*ionz)) *ionz=onz[i]; 143 } 144 PetscFunctionReturn(0); 145 } 146 147