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: calAsq2fe.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 CalAsq2fe(GridData *gdata,PetscReal *vertxyz2,PetscReal *bbox,int *pix,
00030 int slice_id,Mat *pAIpol)
00031 {
00032 int i,j,k,l;
00033 Mat AIpol;
00034
00035
00036 MagparFunctionLogBegin;
00037
00038 int rank,size;
00039 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00040 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00041
00042 PetscPrintf(PETSC_COMM_WORLD,
00043 "Calculating interpolation matrix sq->fe for volume with pid %i\n"
00044 "in bbox (%g,%g,%g) (%g,%g,%g) with %i x %i x %i pixels\n",
00045 slice_id+1,
00046 bbox[0],bbox[1],bbox[2],bbox[ND+0],bbox[ND+1],bbox[ND+2],
00047 pix[0],pix[1],pix[2]
00048 );
00049
00050 if (!rank) {
00051 ierr = MatCreateMPIAIJ(
00052 PETSC_COMM_WORLD,
00053 ND*gdata->ln_vert,ND*pix[0]*pix[1],
00054 ND*gdata->n_vert,ND*pix[0]*pix[1],
00055 NV,PETSC_NULL,NV,PETSC_NULL,
00056 &AIpol
00057 );CHKERRQ(ierr);
00058 }
00059 else {
00060 ierr = MatCreateMPIAIJ(
00061 PETSC_COMM_WORLD,
00062 ND*gdata->ln_vert,0,
00063 ND*gdata->n_vert,ND*pix[0]*pix[1],
00064 NV,PETSC_NULL,NV,PETSC_NULL,
00065 &AIpol
00066 );CHKERRQ(ierr);
00067 }
00068 ierr = MatSetFromOptions(AIpol);CHKERRQ(ierr);
00069
00070 for (l=0;l<gdata->ln_vert;l++) {
00071 PetscReal *v, p[ND];
00072
00073 if (slice_id>=0 && gdata->vertprop[l] != slice_id) continue;
00074
00075
00076
00077
00078
00079 v=vertxyz2+ND*gdata->vertl2g[l];
00080 i=(int) floor(pix[0]*(v[0]-bbox[0])/(bbox[ND+0]-bbox[0]));
00081 j=(int) floor(pix[1]*(v[1]-bbox[1])/(bbox[ND+1]-bbox[1]));
00082
00083 PetscReal a,b;
00084 a=b=-2;
00085 if (i<0) {i=0;a=-1;}
00086 if (j<0) {j=0;b=-1;}
00087 if (i>=pix[0]) {i=pix[0]-1;a=+1;}
00088 if (j>=pix[1]) {j=pix[1]-1;b=+1;}
00089
00090
00091
00092
00093
00094 p[0]=bbox[0]+i*(bbox[ND+0]-bbox[0])/pix[0];
00095 p[1]=bbox[1]+j*(bbox[ND+1]-bbox[1])/pix[1];
00096 p[2]=0;
00097
00098
00099 if (a<-1.5) a=pix[0]*(v[0]-p[0])/(bbox[ND+0]-bbox[0])*2.0-1.0;
00100 if (b<-1.5) b=pix[1]*(v[1]-p[1])/(bbox[ND+1]-bbox[1])*2.0-1.0;
00101
00102 for (k=0;k<ND;k++) {
00103 PetscReal matele;
00104
00105 matele=(1-a)*(1-b)/4.0;
00106 ierr = MatSetValue(AIpol,
00107 ND*gdata->vertl2g[l]+k,
00108 ND*((i+0)*pix[1]+(j+0))+k,
00109 matele,ADD_VALUES);CHKERRQ(ierr);
00110
00111 matele=(1-a)*(1+b)/4.0;
00112 ierr = MatSetValue(AIpol,
00113 ND*gdata->vertl2g[l]+k,
00114 ND*((i+0)*pix[1]+(j+1))+k,
00115 matele,ADD_VALUES);CHKERRQ(ierr);
00116
00117 matele=(1+a)*(1-b)/4.0;
00118 ierr = MatSetValue(AIpol,
00119 ND*gdata->vertl2g[l]+k,
00120 ND*((i+1)*pix[1]+(j+0))+k,
00121 matele,ADD_VALUES);CHKERRQ(ierr);
00122
00123 matele=(1+a)*(1+b)/4.0;
00124 ierr = MatSetValue(AIpol,
00125 ND*gdata->vertl2g[l]+k,
00126 ND*((i+1)*pix[1]+(j+1))+k,
00127 matele,ADD_VALUES);CHKERRQ(ierr);
00128 }
00129 }
00130
00131 PetscInfo(0,"AIpol matrix assembly complete\n");
00132 ierr = MatAssemblyBegin(AIpol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00133 ierr = MatAssemblyEnd(AIpol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00134
00135
00136
00137
00138
00139 PetscPrintf(PETSC_COMM_WORLD,"Checking for missing nodes\n");
00140
00141 Vec test1,test2;
00142 int cnt=0;
00143
00144 ierr = VecDuplicate(gdata->M,&test2);
00145 ierr = VecCreate(PETSC_COMM_WORLD,&test1);CHKERRQ(ierr);
00146 if (!rank) {
00147 ierr = VecSetSizes(test1,ND*pix[0]*pix[1],ND*pix[0]*pix[1]);CHKERRQ(ierr);
00148 }
00149 else{
00150 ierr = VecSetSizes(test1,0,ND*pix[0]*pix[1]);CHKERRQ(ierr);
00151 }
00152 ierr = VecSetFromOptions(test1);CHKERRQ(ierr);
00153
00154 ierr = VecSet(test1,1.0);CHKERRQ(ierr);
00155
00156 ierr = MatMult(AIpol,test1,test2);CHKERRQ(ierr);
00157
00158 PetscReal *atest2;
00159 ierr = VecGetArray(test2,&atest2);CHKERRQ(ierr);
00160
00161 for (l=0;l<gdata->ln_vert;l++) {
00162 PetscReal *v;
00163
00164 if (slice_id>=0 && gdata->vertprop[l] != slice_id) continue;
00165 if (atest2[ND*l]>0.5) continue;
00166
00167
00168
00169
00170
00171 v=vertxyz2+ND*gdata->vertl2g[l];
00172
00173 i=(int) floor(pix[0]*(v[0]-bbox[0])/(bbox[ND+0]-bbox[0]));
00174 j=(int) floor(pix[1]*(v[1]-bbox[1])/(bbox[ND+1]-bbox[1]));
00175
00176
00177
00178
00179
00180
00181
00182
00183 if (i<0) i=0;
00184 if (j<0) j=0;
00185 if (i>=pix[0]) i=pix[0]-1;
00186 if (j>=pix[1]) j=pix[1]-1;
00187
00188 for (k=0;k<ND;k++) {
00189 ierr = MatSetValue(AIpol,
00190 ND*gdata->vertl2g[l]+k,
00191 ND*((i+1)*pix[1]+(j+1))+k,
00192 1.0,ADD_VALUES);CHKERRQ(ierr);
00193 }
00194 cnt++;
00195 }
00196
00197 ierr = VecRestoreArray(test2,&atest2);CHKERRQ(ierr);
00198
00199 PetscInfo(0,"AIpol matrix assembly complete\n");
00200 ierr = MatAssemblyBegin(AIpol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00201 ierr = MatAssemblyEnd(AIpol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00202 ierr = PrintMatInfoAll(AIpol);CHKERRQ(ierr);
00203
00204 PetscSynchronizedPrintf(PETSC_COMM_WORLD,
00205 "<%i>Added %i nodes\n",
00206 rank,cnt
00207 );
00208 PetscSynchronizedFlush(PETSC_COMM_WORLD);
00209
00210 assert(cnt==0);
00211
00212
00213
00214 PetscReal area;
00215 ierr = AreaCal(AIpol,&area);CHKERRQ(ierr);
00216
00217 *pAIpol=AIpol;
00218
00219 MagparFunctionLogReturn(0);
00220 }
00221