00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 static char const Id[] = "$Id: area.c 2698 2009-08-07 14:55:55Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "util.h"
00028
00029 int AreaCal(Mat AIpol, PetscReal *p_area)
00030 {
00031 PetscReal area;
00032
00033 int ngrows,ngcols;
00034 int nlrows,nlcols;
00035
00036 Vec test1,test2;
00037 PetscReal vmin,vmax;
00038 int imin,imax;
00039
00040 MagparFunctionLogBegin;
00041
00042 ierr = MatGetSize(AIpol,&ngrows,&ngcols);CHKERRQ(ierr);
00043 ierr = MatGetLocalSize(AIpol,&nlrows,&nlcols);CHKERRQ(ierr);
00044
00045 ierr = VecCreate(PETSC_COMM_WORLD,&test1);CHKERRQ(ierr);
00046 ierr = VecCreate(PETSC_COMM_WORLD,&test2);CHKERRQ(ierr);
00047 ierr = VecSetSizes(test1,nlcols,ngcols);CHKERRQ(ierr);
00048 ierr = VecSetSizes(test2,nlrows,ngrows);CHKERRQ(ierr);
00049 ierr = VecSetFromOptions(test1);CHKERRQ(ierr);
00050 ierr = VecSetFromOptions(test2);CHKERRQ(ierr);
00051
00052 ierr = VecSet(test1,1.0);CHKERRQ(ierr);
00053
00054 ierr = MatMult(AIpol,test1,test2);CHKERRQ(ierr);
00055
00056
00057 ierr = VecMin(test2,&imin,&vmin);CHKERRQ(ierr);
00058 ierr = VecMax(test2,&imax,&vmax);CHKERRQ(ierr);
00059
00060 if (vmin<0.0 || vmax>1.0+D_EPS) {
00061 PetscPrintf(PETSC_COMM_WORLD,
00062 "need to rescale matrix: %i: %e %i: %e\n",
00063 imin,vmin,imax,vmax
00064 );
00065
00066
00067
00068
00069 PetscReal *atest2;
00070 ierr = VecGetArray(test2,&atest2);CHKERRQ(ierr);
00071 for (int i=0;i<nlrows;i++) {
00072 if (atest2[i]>1.0) atest2[i]=1.0/atest2[i];
00073 }
00074 ierr = VecRestoreArray(test2,&atest2);CHKERRQ(ierr);
00075
00076 ierr = MatDiagonalScale(AIpol,test2,PETSC_NULL);CHKERRQ(ierr);
00077
00078 ierr = MatMult(AIpol,test1,test2);CHKERRQ(ierr);
00079 }
00080
00081
00082 ierr = VecMin(test2,&imin,&vmin);CHKERRQ(ierr);
00083 ierr = VecMax(test2,&imax,&vmax);CHKERRQ(ierr);
00084
00085 PetscPrintf(PETSC_COMM_WORLD,
00086 "min/max: %i: %e %i: %e\n",
00087 imin,vmin,imax,vmax
00088 );
00089
00090 assert(vmin>=0.0);
00091 assert(vmax<=1.0+D_EPS);
00092
00093
00094 ierr = VecSum(test2,&area);CHKERRQ(ierr);
00095 PetscPrintf(PETSC_COMM_WORLD,"area = %g/%i=%g%%\n",area,ngrows,area/ngrows*100);
00096
00097 *p_area=area;
00098
00099 ierr = VecDestroy(test1);
00100 ierr = VecDestroy(test2);
00101
00102 MagparFunctionLogReturn(0);
00103 }
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114 int PidSharedArea(GridData *gdata,int pid1,int pid2,PetscReal *sharedarea_area,PetscReal *sharedarea_point,PetscReal *sharedarea_norm)
00115 {
00116 MagparFunctionLogBegin;
00117
00118 int rank,size;
00119 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00120 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00121
00122 if (pid1==pid2)
00123 SETERRQ2(PETSC_ERR_ARG_INCOMP,"Pids have to be different to calculate shared area (%i==%i)!\n",pid1+1,pid2+1);
00124
00125 if (pid1>=gdata->n_prop)
00126 SETERRQ2(PETSC_ERR_ARG_INCOMP,"pid1 too large (%i>%i)!\n",pid1+1,gdata->n_prop);
00127 if (pid2>=gdata->n_prop)
00128 SETERRQ2(PETSC_ERR_ARG_INCOMP,"pid2 too large (%i>%i)!\n",pid1+1,gdata->n_prop);
00129 if (pid1<=-1 && pid2<=-1)
00130 SETERRQ2(PETSC_ERR_ARG_INCOMP,"pids too small (%i,%i)<1!\n",pid1+1,pid2+1);
00131
00132 PetscPrintf(PETSC_COMM_WORLD,"Search shared area for pids %i, %i\n",pid1+1,pid2+1);
00133
00134
00135 int low,high;
00136 ierr = VecGetOwnershipRange(gdata->elevol,&low,&high);CHKERRQ(ierr);
00137
00138 int *elevertall;
00139
00140 #ifdef HAVE_ELEVERTALL
00141 elevertall=gdata->elevertall;
00142 #else
00143 if (size>1) {
00144 ierr = PetscMalloc(NV*gdata->n_ele*sizeof(int),&elevertall);CHKERRQ(ierr);
00145 ierr = PetscMemcpy(elevertall+NV*low,gdata->elevert,NV*gdata->ln_ele*sizeof(int));CHKERRQ(ierr);
00146
00147 int *p;
00148 p=elevertall;
00149
00150 int recvcount[size];
00151 recvcount[rank]=NV*gdata->ln_ele;
00152 for (int i=0;i<size;i++) {
00153 ierr = MPI_Bcast(recvcount+i,1,MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
00154
00155 ierr = MPI_Bcast(p,recvcount[i],MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
00156 p+=recvcount[i];
00157 }
00158 assert(p==elevertall+NV*gdata->n_ele);
00159 }
00160 else {
00161 elevertall=gdata->elevert;
00162 }
00163 #endif
00164
00165 int *elepropall;
00166 if (size>1) {
00167 ierr = PetscMalloc(gdata->n_ele*sizeof(int),&elepropall);CHKERRQ(ierr);
00168 ierr = PetscMemcpy(elepropall+low,gdata->eleprop,gdata->ln_ele*sizeof(int));CHKERRQ(ierr);
00169
00170 int *p;
00171 p=elepropall;
00172
00173 int recvcount[size];
00174 recvcount[rank]=gdata->ln_ele;
00175 for (int i=0;i<size;i++) {
00176 ierr = MPI_Bcast(recvcount+i,1,MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
00177
00178 ierr = MPI_Bcast(p,recvcount[i],MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
00179 p+=recvcount[i];
00180 }
00181 assert(p==elepropall+gdata->n_ele);
00182 }
00183 else {
00184 elepropall=gdata->eleprop;
00185 }
00186
00187 int *elevertcoupled;
00188 ierr = PetscMalloc(NV*gdata->n_ele*sizeof(int),&elevertcoupled);CHKERRQ(ierr);
00189 ierr = PetscMemcpy(elevertcoupled,elevertall,NV*gdata->n_ele*sizeof(int));CHKERRQ(ierr);
00190
00191 #ifndef HAVE_ELEVERTALL
00192 if (size>1) {
00193 ierr = PetscFree(elevertall);CHKERRQ(ierr);
00194 }
00195 #endif
00196
00197 PetscReal sarea_area;
00198 sarea_area=0;
00199
00200 PetscReal sarea_point[ND];
00201 PetscReal sarea_norm[ND];
00202 for (int i=0;i<ND;i++) {
00203 sarea_point[i]=sarea_norm[i]=0.0;
00204 }
00205
00206
00207 if (!rank) {
00208
00209 #ifdef ADDONS
00210
00211 for (int i=0;i<NV*gdata->n_ele;i++) {
00212 if (gdata->vertcoupl[elevertcoupled[i]]>=0 &&
00213 gdata->vertcoupl[elevertcoupled[i]]<elevertcoupled[i]) {
00214 elevertcoupled[i]=gdata->vertcoupl[elevertcoupled[i]];
00215 }
00216 }
00217 #endif
00218
00219 PetscPrintf(PETSC_COMM_WORLD,"Calculating mesh connectivity\n");
00220
00221 int *ia,*ja;
00222 ierr = Mesh2Dual(gdata->n_ele,gdata->n_vert,elevertcoupled,&ia,&ja);CHKERRQ(ierr);
00223
00224 for (int i=1;i<gdata->n_ele+1;i++) {
00225 if (ia[i]-ia[i-1]>NF) {
00226 SETERRQ2(PETSC_ERR_MEMC,"Found element with %i neighbors (>NF)!\n",ia[i]-ia[i-1],NF);
00227 }
00228 }
00229
00230
00231 assert(ia[gdata->n_ele]<NF*gdata->n_ele);
00232
00233 PetscPrintf(PETSC_COMM_WORLD,"Searching shared faces\n");
00234 for (int i=0; i<gdata->n_ele; i++) {
00235 if (elepropall[i]!=pid1) continue;
00236
00237 if (pid2>=0) {
00238 for (int j2=ia[i];j2<ia[i+1];j2++) {
00239 int j;
00240 j=ja[j2];
00241
00242 if (elepropall[j]!=pid2) continue;
00243
00244 int s[NV-1];
00245 int k;
00246 k=TetSharedVertices(elevertcoupled+NV*i,elevertcoupled+NV*j,s);
00247 assert(k==NV-1);
00248
00249 PetscReal a;
00250 PetscReal n[ND];
00251 a=triarea(
00252 gdata->vertxyz+ND*s[0],
00253 gdata->vertxyz+ND*s[1],
00254 gdata->vertxyz+ND*s[2],
00255 n
00256 );
00257
00258 sarea_area+=a;
00259
00260 for (int p=0;p<ND;p++) {
00261 if (n[2]>=0.0) {
00262 sarea_norm[p]+=n[p]*a;
00263 }
00264 else {
00265 sarea_norm[p]-=n[p]*a;
00266 }
00267 for (int q=0;q<NV-1;q++) {
00268 sarea_point[p]+=gdata->vertxyz[ND*s[q]+p]*a/(NV-1.0);
00269 }
00270 }
00271 }
00272 }
00273 else if (pid2==-1) {
00274
00275 SETERRQ1(PETSC_ERR_ARG_INCOMP,"Program branch for pid2==%i not implemented!\n",pid2);
00276 }
00277 else if (pid2==-2) {
00278
00279 SETERRQ1(PETSC_ERR_ARG_INCOMP,"Program branch for pid2==%i not implemented!\n",pid2);
00280 }
00281 }
00282
00283 ierr = PetscFree(ia);CHKERRQ(ierr);
00284 ierr = PetscFree(ja);CHKERRQ(ierr);
00285
00286 }
00287
00288 ierr = MPI_Bcast(&sarea_area,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00289 ierr = MPI_Bcast(sarea_point,ND,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00290 ierr = MPI_Bcast(sarea_norm,ND,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00291
00292 if (size>1) {
00293 ierr = PetscFree(elepropall);CHKERRQ(ierr);
00294 }
00295
00296 *sharedarea_area=sarea_area;
00297
00298 if (sharedarea_point!=PETSC_NULL) {
00299 my_dscal(ND,1.0/sarea_area,sarea_point,1);
00300 my_dcopy(ND,sarea_point,1,sharedarea_point,1);
00301 }
00302
00303 if (sharedarea_norm!=PETSC_NULL) {
00304 my_dscal(ND,1.0/sarea_area,sarea_norm,1);
00305 my_dcopy(ND,sarea_norm,1,sharedarea_norm,1);
00306 }
00307
00308 MagparFunctionLogReturn(0);
00309 }
00310