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: calAfe2sq.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "util.h"
00028
00029 int findeleslice(int ln_ele,int *elevert,PetscReal *vertxyz2,PetscReal zcut,int *nele,int *sliceele)
00030 {
00031
00032 MagparFunctionLogBegin;
00033
00034 int rank,size;
00035 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00036 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00037
00038 int cnt=0;
00039 for (int i=0;i<ln_ele;i++) {
00040 PetscReal sgn;
00041
00042 sgn=vertxyz2[ND*elevert[NV*i+0]+2]-zcut;
00043
00044 for (int j=1;j<NV;j++) {
00045
00046
00047
00048
00049 if (sgn*(vertxyz2[ND*elevert[NV*i+j]+2]-zcut) <= 0.0) {
00050
00051 sliceele[cnt++]=i;
00052 break;
00053 }
00054 }
00055 }
00056
00057 PetscSynchronizedPrintf(PETSC_COMM_WORLD,
00058 "<%i>Found %i elements cut by slice plane\n",
00059 rank,cnt
00060 );
00061 PetscSynchronizedFlush(PETSC_COMM_WORLD);
00062
00063 *nele=cnt;
00064
00065 MagparFunctionLogReturn(0);
00066 }
00067
00068 int CalAfe2sq(GridData *gdata,PetscReal *vertxyz2,PetscReal *bbox,int *pix,int slice_id,Mat *pAIpol)
00069 {
00070 MagparFunctionLogBegin;
00071
00072 int rank,size;
00073 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00074 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00075
00076 PetscPrintf(PETSC_COMM_WORLD,
00077 "Calculating interpolation matrix fe->sq for volume with pid %i\n"
00078 "in bbox (%g,%g,%g) (%g,%g,%g) with %i x %i x %i pixels\n",
00079 slice_id+1,
00080 bbox[0],bbox[1],bbox[2],bbox[ND+0],bbox[ND+1],bbox[ND+2],
00081 pix[0],pix[1],pix[2]
00082 );
00083
00084
00085
00086
00087 int *sliceele;
00088 ierr = PetscMalloc(gdata->ln_ele*sizeof(int),&sliceele);CHKERRQ(ierr);
00089 int cnt;
00090 ierr = findeleslice(gdata->ln_ele,gdata->elevert,vertxyz2,(bbox[2]+bbox[ND+2])/2.0,&cnt,sliceele);CHKERRQ(ierr);
00091
00092 int sum;
00093 ierr=MPI_Allreduce((int*)&cnt,(int*)&sum,1,MPI_INT,MPI_SUM,PETSC_COMM_WORLD);CHKERRQ(ierr);
00094 PetscPrintf(PETSC_COMM_WORLD,
00095 "Found total %i elements cut by slice plane.\n",
00096 sum
00097 );
00098
00099 if (sum==0) {
00100 PetscPrintf(PETSC_COMM_WORLD,"Slice plane outside the FE mesh.\n");
00101
00102 ierr = PetscFree(sliceele);CHKERRQ(ierr);
00103 *pAIpol=PETSC_NULL;
00104
00105 MagparFunctionLogReturn(0);
00106 }
00107
00108 Mat AIpol;
00109 if (!rank) {
00110 ierr = MatCreateMPIAIJ(
00111 PETSC_COMM_WORLD,
00112 ND*pix[0]*pix[1],ND*gdata->ln_vert,
00113 ND*pix[0]*pix[1],ND*gdata->n_vert,
00114 2*NV,PETSC_NULL,2*NV,PETSC_NULL,
00115 &AIpol
00116 );CHKERRQ(ierr);
00117 }
00118 else {
00119 ierr = MatCreateMPIAIJ(
00120 PETSC_COMM_WORLD,
00121 0,ND*gdata->ln_vert,
00122 ND*pix[0]*pix[1],ND*gdata->n_vert,
00123 NV,PETSC_NULL,NV,PETSC_NULL,
00124 &AIpol
00125 );CHKERRQ(ierr);
00126 }
00127 ierr = MatSetFromOptions(AIpol);CHKERRQ(ierr);
00128
00129 PetscPrintf(PETSC_COMM_WORLD,"Calculating interpolation matrix...\n");
00130
00131 for (int k=0; k<cnt; k++) {
00132 PetscReal elebox[2*ND];
00133 int pixbox[2*ND];
00134 int t_ele;
00135
00136 ierr = ProgressBar(k,cnt,10);
00137
00138
00139 t_ele=sliceele[k];
00140
00141
00142
00143
00144
00145
00146 if (slice_id>=0 && gdata->eleprop[t_ele]!=slice_id) continue;
00147 if (slice_id<=-2 && gdata->eleprop[t_ele]==-slice_id-1) continue;
00148
00149
00150 for (int i=0;i<ND;i++) {
00151 elebox[i]=PETSC_MAX;
00152 elebox[ND+i]=PETSC_MIN;
00153 }
00154 for (int l=0; l<NV; l++) {
00155 for (int i=0;i<ND;i++) {
00156 elebox[i]= PetscMin(elebox[i], vertxyz2[ND*gdata->elevert[NV*t_ele+l]+i]);
00157 elebox[ND+i]=PetscMax(elebox[ND+i],vertxyz2[ND*gdata->elevert[NV*t_ele+l]+i]);
00158 }
00159 }
00160
00161 for (int i=0;i<ND;i++) {
00162 assert(elebox[ND+i]>elebox[i]);
00163 }
00164
00165
00166 for (int i=0;i<ND;i++) {
00167 pixbox[i]= (int) floor(pix[i]*(elebox[ i]-bbox[i])/(bbox[ND+i]-bbox[i]));
00168 pixbox[ND+i]=(int) ceil(pix[i]*(elebox[ND+i]-bbox[i])/(bbox[ND+i]-bbox[i]));
00169 }
00170
00171
00172 for (int i=pixbox[0];i<=pixbox[ND+0];i++) {
00173
00174 if (i<0 || i>=pix[0]) continue;
00175
00176 for (int j=pixbox[1];j<=pixbox[ND+1];j++) {
00177
00178 if (j<0 || j>=pix[1]) continue;
00179
00180 PetscReal p[ND],rhs[ND];
00181
00182
00183
00184
00185
00186 p[0]=bbox[0]+i*(bbox[ND+0]-bbox[0])/pix[0];
00187 p[1]=bbox[1]+j*(bbox[ND+1]-bbox[1])/pix[1];
00188 p[2]=(bbox[2]+bbox[ND+2])/2.0;
00189
00190 if (
00191 barycent(
00192 p,
00193 vertxyz2+ND*gdata->elevert[NV*t_ele+0],
00194 vertxyz2+ND*gdata->elevert[NV*t_ele+1],
00195 vertxyz2+ND*gdata->elevert[NV*t_ele+2],
00196 vertxyz2+ND*gdata->elevert[NV*t_ele+3],
00197 rhs
00198 )
00199 ) {
00200 for (int l=0;l<ND;l++) {
00201 PetscReal matele;
00202
00203 matele=(1-rhs[0]-rhs[1]-rhs[2]);
00204 ierr = MatSetValue(AIpol,
00205 ND*(i*pix[1]+j)+l,
00206 ND*gdata->elevert[NV*t_ele+3]+l,
00207 matele,ADD_VALUES);CHKERRQ(ierr);
00208
00209 matele=rhs[0];
00210 ierr = MatSetValue(AIpol,
00211 ND*(i*pix[1]+j)+l,
00212 ND*gdata->elevert[NV*t_ele+0]+l,
00213 matele,ADD_VALUES);CHKERRQ(ierr);
00214
00215 matele=rhs[1];
00216 ierr = MatSetValue(AIpol,
00217 ND*(i*pix[1]+j)+l,
00218 ND*gdata->elevert[NV*t_ele+1]+l,
00219 matele,ADD_VALUES);CHKERRQ(ierr);
00220
00221 matele=rhs[2];
00222 ierr = MatSetValue(AIpol,
00223 ND*(i*pix[1]+j)+l,
00224 ND*gdata->elevert[NV*t_ele+2]+l,
00225 matele,ADD_VALUES);CHKERRQ(ierr);
00226 }
00227 }
00228 }
00229 }
00230 }
00231 ierr = ProgressBar(1,1,10);
00232
00233 ierr = PetscFree(sliceele);CHKERRQ(ierr);
00234
00235 PetscInfo(0,"AIpol matrix assembly complete\n");
00236 ierr = MatAssemblyBegin(AIpol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00237 ierr = MatAssemblyEnd(AIpol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00238 ierr = PrintMatInfoAll(AIpol);CHKERRQ(ierr);
00239
00240
00241 PetscReal area;
00242 ierr = AreaCal(AIpol,&area);CHKERRQ(ierr);
00243
00244 *pAIpol=AIpol;
00245
00246 MagparFunctionLogReturn(0);
00247 }
00248