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,nsize,ivtx,n_nnz,n_onz; 80 DM_Moab *dmmoab = (DM_Moab*)dm->data; 81 const moab::EntityHandle *connect; 82 moab::Range adjs,found,allvlocal,allvghost; 83 moab::Range::iterator iter,jter; 84 std::vector<moab::EntityHandle> storage; 85 moab::EntityHandle vtx; 86 moab::ErrorCode merr; 87 88 PetscFunctionBegin; 89 bs = dmmoab->bs; 90 nloc = dmmoab->nloc; 91 nsize = (isbaij ? nloc:nloc*bs); 92 93 /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */ 94 merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,allvlocal,true);MBERRNM(merr); 95 merr = dmmoab->pcomm->filter_pstatus(allvlocal,PSTATUS_NOT_OWNED,PSTATUS_NOT,-1,&adjs);MBERRNM(merr); 96 allvghost = moab::subtract(allvlocal, adjs); 97 98 /* loop over the locally owned vertices and figure out the NNZ pattern using connectivity information */ 99 for(iter = dmmoab->vowned->begin(),ivtx=0; iter != dmmoab->vowned->end(); iter++,ivtx++) { 100 101 vtx = *iter; 102 adjs.clear(); 103 /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects 104 to the current vertex. We can then decipher if a vertex is ghosted or not and compute the 105 non-zero pattern accordingly. */ 106 merr = dmmoab->mbiface->get_adjacencies(&vtx,1,dmmoab->dim,false,adjs,moab::Interface::INTERSECT); 107 108 /* reset counters */ 109 n_nnz=n_onz=0; 110 found.clear(); 111 112 /* loop over vertices and update the number of connectivity */ 113 for(jter = adjs.begin(); jter != adjs.end(); jter++) { 114 115 /* Get connectivity information in canonical ordering for the local element */ 116 merr = dmmoab->mbiface->get_connectivity(*jter,connect,vpere,false,&storage);MBERRNM(merr); 117 118 /* loop over each element connected to the adjacent vertex and update as needed */ 119 for (i=0; i<vpere; ++i) { 120 if (connect[i] == vtx || found.find(connect[i]) != found.end()) continue; /* make sure we don't double count shared vertices */ 121 if (allvghost.find(connect[i]) != allvghost.end()) n_onz++; /* update out-of-proc onz */ 122 else n_nnz++; /* else local vertex */ 123 found.insert(connect[i]); 124 } 125 } 126 127 if (isbaij) { 128 nnz[ivtx]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 129 if (onz) onz[ivtx]=n_onz; /* add ghost non-owned nodes */ 130 } 131 else { /* AIJ matrices */ 132 for (f=0;f<dmmoab->numFields;f++) { 133 nnz[dmmoab->numFields*ivtx+f]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 134 if (onz) onz[dmmoab->numFields*ivtx+f]=n_onz; /* add ghost non-owned nodes */ 135 } 136 } 137 } 138 139 if (innz) *innz=0; 140 if (ionz) *ionz=0; 141 for (i=0;i<nsize;i++) { 142 nnz[i]+=1; /* self count the node */ 143 /* check if we got overzealous */ 144 nnz[i]=(nnz[i]>dmmoab->nloc ? dmmoab->nloc:nnz[i]); 145 if (!isbaij) { 146 nnz[i]*=bs; 147 if (onz) onz[i]*=bs; 148 } 149 150 /* update innz and ionz based on local maxima */ 151 if (innz && (nnz[i]>*innz)) *innz=nnz[i]; 152 if ((ionz && onz) && (onz[i]>*ionz)) *ionz=onz[i]; 153 } 154 PetscFunctionReturn(0); 155 } 156 157 158 #undef __FUNCT__ 159 #define __FUNCT__ "DMMoab_MatFillMatrixEntries_Private" 160 PetscErrorCode DMMoab_MatFillMatrixEntries_Private(DM dm, Mat A) 161 { 162 DM_Moab *dmmoab = (DM_Moab*)dm->data; 163 PetscInt nconn = 0,prev_nconn = 0; 164 const moab::EntityHandle *connect; 165 PetscScalar *locala=NULL; 166 PetscInt *dof_indices=NULL; 167 PetscErrorCode ierr; 168 169 PetscFunctionBegin; 170 /* loop over local elements */ 171 for(moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) { 172 const moab::EntityHandle ehandle = *iter; 173 174 /* Get connectivity information: */ 175 ierr = DMMoabGetElementConnectivity(dm, ehandle, &nconn, &connect);CHKERRQ(ierr); 176 177 /* if we have mixed elements or arrays have not been initialized - Allocate now */ 178 if (prev_nconn != nconn) { 179 if (locala) { 180 ierr = PetscFree(locala);CHKERRQ(ierr); 181 ierr = PetscFree(dof_indices);CHKERRQ(ierr); 182 } 183 ierr = PetscMalloc(sizeof(PetscScalar)*nconn*nconn*dmmoab->numFields*dmmoab->numFields,&locala);CHKERRQ(ierr); 184 ierr = PetscMemzero(locala,sizeof(PetscScalar)*nconn*nconn*dmmoab->numFields*dmmoab->numFields);CHKERRQ(ierr); 185 ierr = PetscMalloc(sizeof(PetscInt)*nconn,&dof_indices);CHKERRQ(ierr); 186 prev_nconn=nconn; 187 } 188 189 /* get the global DOF number to appropriately set the element contribution in the RHS vector */ 190 ierr = DMMoabGetDofsBlockedLocal(dm, nconn, connect, dof_indices);CHKERRQ(ierr); 191 192 /* set the values directly into appropriate locations. Can alternately use VecSetValues */ 193 ierr = MatSetValuesBlockedLocal(A, nconn, dof_indices, nconn, dof_indices, locala, INSERT_VALUES);CHKERRQ(ierr); 194 } 195 196 /* clean up memory */ 197 ierr = PetscFree(locala);CHKERRQ(ierr); 198 ierr = PetscFree(dof_indices);CHKERRQ(ierr); 199 200 /* finish assembly */ 201 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 202 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 203 PetscFunctionReturn(0); 204 } 205 206