1 #include <stdio.h> 2 #include <stdlib.h> 3 #include <assert.h> 4 #include <math.h> 5 6 #include <mpi.h> 7 #include <petsc.h> 8 9 #include "common_c.h" 10 #include "FCMangle.h" 11 12 #define fillsparsecpetscs FortranCInterface_GLOBAL_(fillsparsecpetscs, FILLSPARSECPETSCS) 13 #define fillsparsecpetscc FortranCInterface_GLOBAL_(fillsparsecpetscc, FILLSPARSECPETSCC) 14 15 #define COLMAJ2D(row,col,numrow) (row-1)+(col-1)*numrow 16 #define COLMAJ3D(a,b,c,amax,bmax,cmax) (a-1)+amax*((b-1)+bmax*(c-1)) 17 #define ROWMAJ2D_ONE(row,col,numcol) (row-1)*numcol+(col-1) 18 typedef long long int gcorp_t; 19 20 void fillsparsecpetscs(gcorp_t* ieng, double* EGmass, Mat* lhsP) 21 { 22 int npro = propar.npro; 23 int nshl = shpdat.nshl; 24 double* mb = (double*) malloc(sizeof(double)*nshl*nshl); //block to insert 25 int e,i,j,aa; //following along with fillsparse.f 26 PetscInt* locat = (PetscInt*) malloc(sizeof(PetscInt)*nshl); 27 for(e=0;e<npro;e++) 28 { 29 for(aa=0;aa<nshl;aa++) locat[aa]=ieng[e+npro*aa]-1; 30 // for(aa=0;aa<nshl;aa++) assert(locat[aa]>=0); 31 for (i=0; i<nshl; i++) { // fill up Ke with respective egmass 32 for (j=0; j<nshl; j++) { 33 mb[nshl*i + j] = EGmass[e + npro*(i + nshl*j)]; 34 } 35 } 36 //MatSetValuesBlocked(*lhsP, nshl , locat, nshl, locat, mb, ADD_VALUES); 37 PetscInt petsc_nshl; 38 petsc_nshl = (PetscInt) nshl; 39 MatSetValues(*lhsP, petsc_nshl , locat, petsc_nshl, locat, mb, ADD_VALUES); 40 } 41 free(mb); 42 free(locat); 43 } 44 void fillsparsecpetscc(gcorp_t* ieng, double* EGmass, Mat* lhsP) 45 { 46 int npro = propar.npro; 47 int nshl = shpdat.nshl; 48 int nedof = conpar.nedof; 49 double* mb = (double*) malloc(sizeof(double)*nedof*nedof); //block to insert 50 int e,i,j,aa; //following along with fillsparse.f 51 //int* locat = (int*) malloc(sizeof(int)*nshl); 52 PetscInt* locat = (PetscInt*) malloc(sizeof(PetscInt)*nshl); 53 for(e=0;e<npro;e++) 54 { 55 for(aa=0;aa<nshl;aa++) locat[aa]=ieng[e+npro*aa]-1; 56 // for(aa=0;aa<nshl;aa++) assert(locat[aa]>=0); 57 for (i=0; i<nedof; i++) { /* fill up Ke with respective egmass */ 58 for (j=0; j<nedof; j++) { 59 mb[nedof*i + j] = EGmass[e + npro*(i + nedof*j)]; 60 } 61 } 62 //MatSetValuesBlocked(*lhsP, nshl , locat, nshl, locat, mb, ADD_VALUES); 63 PetscInt petsc_nshl; 64 petsc_nshl = (PetscInt) nshl; 65 MatSetValuesBlocked(*lhsP, petsc_nshl , locat, petsc_nshl, locat, mb, ADD_VALUES); 66 } 67 free(mb); 68 free(locat); 69 } 70