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: regrefine.c 3037 2010-03-30 07:10:51Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "init/init.h"
00028 #include "util/util.h"
00029
00030 int RegularRefinement(GridData *gdata)
00031 {
00032 const int midnod[4][4]=
00033 {
00034 {C_UNK,4,6,7},
00035 {C_UNK,C_UNK,5,8},
00036 {C_UNK,C_UNK,C_UNK,9},
00037 {C_UNK,C_UNK,C_UNK,C_UNK}
00038 };
00039 const int child[8][4]=
00040 {
00041 {0,4,6,7},
00042 {4,1,5,8},
00043 {6,5,2,9},
00044 {7,8,9,3},
00045 {4,6,7,8},
00046 {4,5,6,8},
00047 {6,8,5,9},
00048 {6,7,8,9}
00049 };
00050
00051
00052 MagparFunctionLogBegin;
00053
00054 int rank,size;
00055 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00056 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00057
00058 int refine;
00059 PetscTruth flg;
00060 ierr = PetscOptionsGetInt(PETSC_NULL,"-refine",(PetscInt*)&refine,&flg);CHKERRQ(ierr);
00061 if (flg!=PETSC_TRUE) {
00062 refine=0;
00063 PetscPrintf(PETSC_COMM_WORLD,
00064 "Option -refine not found, using default value: %i\n",
00065 refine
00066 );
00067 }
00068
00069 PetscPrintf(PETSC_COMM_WORLD,
00070 "Requested number of refinement steps: %i\n",
00071 refine
00072 );
00073
00074 if (rank) {
00075 for (int r=0;r<refine;r++) {
00076 Vec Mnew;
00077 ierr = VecCreate(PETSC_COMM_WORLD,&Mnew);CHKERRQ(ierr);
00078 ierr = VecSetSizes(Mnew,ND*gdata->ln_vert,PETSC_DECIDE);CHKERRQ(ierr);
00079 ierr = VecSetFromOptions(Mnew);CHKERRQ(ierr);
00080 ierr = VecSetBlockSize(Mnew,ND);CHKERRQ(ierr);
00081
00082 ierr = VecDestroy(gdata->M);CHKERRQ(ierr);
00083 gdata->M=Mnew;
00084
00085 PetscReal devNorm;
00086 devNorm=RenormVec(gdata->M,1.0,0.0,ND);
00087 }
00088 MagparFunctionLogReturn(0);
00089 }
00090
00091
00092
00093
00094
00095
00096
00097 for (int r=0;r<refine;r++) {
00098 PetscPrintf(PETSC_COMM_WORLD,"Starting refinement step: %i\n",r+1);
00099
00100 int *nxadj,*nadjncy;
00101 ierr = Mesh2Nodal(gdata->n_ele,gdata->n_vert,gdata->elevert,&nxadj,&nadjncy);CHKERRQ(ierr);
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111 int n_edge;
00112 n_edge=nxadj[gdata->n_vert]/2;
00113 PetscPrintf(PETSC_COMM_WORLD,
00114 "Number of edges: %i\n",
00115 n_edge
00116 );
00117
00118 PetscPrintf(PETSC_COMM_WORLD,
00119 "Generating a FE mesh with %i nodes and %i elements.\n",
00120 gdata->n_vert+n_edge,
00121 8*gdata->n_ele
00122 );
00123
00124
00125
00126 int *nvert;
00127 ierr = PetscMalloc(2*n_edge*sizeof(int),&nvert);CHKERRQ(ierr);
00128
00129 for (int i=0;i<2*n_edge;i++) nvert[i]=C_UNK;
00130
00131
00132 int cnt;
00133 cnt=gdata->n_vert;
00134
00135
00136 PetscReal *nvertxyz;
00137 ierr = PetscMalloc(ND*(gdata->n_vert+n_edge)*sizeof(PetscReal),&nvertxyz);CHKERRQ(ierr);
00138
00139 ierr = PetscMemcpy(nvertxyz,gdata->vertxyz,ND*gdata->n_vert*sizeof(PetscReal));CHKERRQ(ierr);
00140
00141
00142 int *nelevert;
00143 ierr = PetscMalloc(NV*8*gdata->n_ele*sizeof(int),&nelevert);CHKERRQ(ierr);
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153 for (int i=0;i<gdata->n_ele;i++) {
00154
00155 int t_v[10];
00156 ierr = PetscMemcpy(t_v,gdata->elevert+i*NV,NV*sizeof(int));CHKERRQ(ierr);
00157
00158
00159 for (int m=0;m<NV-1;m++) {
00160 for (int n=m+1;n<NV;n++) {
00161
00162
00163
00164
00165 int lid,hid;
00166 lid=PetscMin(t_v[m],t_v[n]);
00167 hid=PetscMax(t_v[m],t_v[n]);
00168
00169 int nnb;
00170 nnb=nxadj[lid+1]-nxadj[lid];
00171
00172
00173
00174
00175 int j;
00176 for (j=0;j<nnb;j++) {
00177 if (nadjncy[nxadj[lid]+j]==hid) {
00178
00179
00180
00181
00182 if (nvert[nxadj[lid]+j] == C_UNK) {
00183 nvert[nxadj[lid]+j]=cnt;
00184
00185 for (int k=0;k<ND;k++) {
00186 nvertxyz[ND*cnt+k]=
00187 0.5*(gdata->vertxyz[ND*lid+k]+gdata->vertxyz[ND*hid+k]);
00188 }
00189
00190 cnt++;
00191 }
00192 break;
00193 }
00194 }
00195
00196
00197
00198 assert(j<nxadj[lid+1]-nxadj[lid]);
00199
00200 t_v[midnod[m][n]]=nvert[nxadj[lid]+j];
00201 }
00202 }
00203
00204 for (int j=0;j<8;j++) {
00205 for (int k=0;k<NV;k++) {
00206 nelevert[NV*8*i+NV*j+k]=t_v[child[j][k]];
00207 }
00208 }
00209 }
00210
00211 assert(cnt==gdata->n_vert+n_edge);
00212
00213 ierr = PetscFree(gdata->vertxyz);CHKERRQ(ierr);
00214 gdata->vertxyz=nvertxyz;
00215
00216 ierr = PetscFree(gdata->elevert);CHKERRQ(ierr);
00217 gdata->elevert=nelevert;
00218
00219
00220
00221 Vec Mnew;
00222 ierr = VecCreate(PETSC_COMM_WORLD,&Mnew);CHKERRQ(ierr);
00223 ierr = VecSetSizes(Mnew,ND*(gdata->n_vert+n_edge),PETSC_DECIDE);CHKERRQ(ierr);
00224 ierr = VecSetFromOptions(Mnew);CHKERRQ(ierr);
00225 ierr = VecSetBlockSize(Mnew,ND);CHKERRQ(ierr);
00226
00227 int cnt2;
00228 ierr = VecGetLocalSize(gdata->M,(PetscInt*)&cnt2);CHKERRQ(ierr);
00229 assert(cnt2==ND*gdata->n_vert);
00230 ierr = VecGetSize(gdata->M,(PetscInt*)&cnt2);CHKERRQ(ierr);
00231 assert(cnt2==ND*gdata->n_vert);
00232
00233 ierr = VecGetLocalSize(Mnew,(PetscInt*)&cnt2);CHKERRQ(ierr);
00234 assert(cnt2==ND*(gdata->n_vert+n_edge));
00235
00236 PetscReal *t_M, *t_Mnew;
00237 ierr = VecGetArray(gdata->M,&t_M);CHKERRQ(ierr);
00238 ierr = VecGetArray(Mnew,&t_Mnew);CHKERRQ(ierr);
00239
00240 PetscPrintf(PETSC_COMM_WORLD,"Updating M\n");
00241
00242 for (int i=0;i<gdata->n_vert;i++) {
00243 for (int k=0;k<ND;k++) {
00244 t_Mnew[ND*i+k]=t_M[ND*i+k];
00245 }
00246 for (int j=nxadj[i];j<nxadj[i+1];j++) {
00247 if (i<nadjncy[j]) {
00248 for (int k=0;k<ND;k++) {
00249 t_Mnew[ND*nvert[j]+k]=
00250 (t_M[ND*i+k]+t_M[ND*nadjncy[j]+k])/2.0;
00251 }
00252 }
00253 }
00254 }
00255
00256 ierr = VecRestoreArray(gdata->M,&t_M);CHKERRQ(ierr);
00257 ierr = VecRestoreArray(Mnew,&t_Mnew);CHKERRQ(ierr);
00258 ierr = VecDestroy(gdata->M);CHKERRQ(ierr);
00259 gdata->M=Mnew;
00260
00261 ierr = PetscFree(nvert);CHKERRQ(ierr);
00262
00263 PetscReal devNorm;
00264 devNorm=RenormVec(gdata->M,1.0,0.0,ND);
00265
00266
00267 int *neleprop;
00268 ierr = PetscMalloc(8*gdata->n_ele*sizeof(int),&neleprop);CHKERRQ(ierr);
00269 for (int i=0;i<gdata->n_ele;i++) {
00270 for (int j=0;j<8;j++) {
00271 neleprop[8*i+j]=gdata->eleprop[i];
00272 }
00273 }
00274 ierr = PetscFree(gdata->eleprop);CHKERRQ(ierr);
00275 gdata->eleprop=neleprop;
00276
00277 ierr = PetscFree(nxadj);CHKERRQ(ierr);
00278 ierr = PetscFree(nadjncy);CHKERRQ(ierr);
00279
00280
00281 gdata->n_vert += n_edge;
00282 gdata->ln_vert = gdata->n_vert;
00283
00284 gdata->n_ele *= 8;
00285 gdata->ln_ele = gdata->n_ele;
00286 }
00287
00288
00289 MagparFunctionLogReturn(0);
00290 }
00291