modifyprop_ser.c
Go to the documentation of this file.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: modifyprop_ser.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 "init.h"
00028 #include "io/magpario.h"
00029 #include "util/util.h"
00030
00031 int ModifyPropSer(GridData *gdata)
00032 {
00033 MagparFunctionLogBegin;
00034
00035 int rank,size;
00036 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00037 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00038
00039 if (rank) {
00040 MagparFunctionLogReturn(0);
00041 }
00042
00043 PetscTruth flg;
00044 int nsliceprop;
00045 ierr = PetscOptionsGetInt(PETSC_NULL,"-nslicepropser",(PetscInt*)&nsliceprop,&flg);CHKERRQ(ierr);
00046 if (flg!=PETSC_TRUE) {
00047
00048 nsliceprop=0;
00049 PetscPrintf(PETSC_COMM_WORLD,
00050 "Option -nslicepropser not found, using default value %i!\n",
00051 nsliceprop
00052 );
00053 }
00054
00055
00056 if (nsliceprop<=0) {
00057 PetscPrintf(PETSC_COMM_WORLD,
00058 "nslicepropser==%i, nothing modified!\n",
00059 nsliceprop
00060 );
00061 MagparFunctionLogReturn(0);
00062 }
00063
00064
00065 int *ia,*ja;
00066 ierr = Mesh2Dual(gdata->n_ele,gdata->n_vert,gdata->elevert,&ia,&ja);CHKERRQ(ierr);
00067
00068 int *elebnd;
00069 ierr = PetscMalloc(gdata->n_ele*sizeof(int),&elebnd);CHKERRQ(ierr);
00070
00071
00072 for (int i=0;i<gdata->n_ele;i++) {
00073 elebnd[i]=C_INT;
00074
00075 int t_nnb;
00076 t_nnb=ia[i+1]-ia[i];
00077 if (t_nnb<NF) {
00078 elebnd[i]=C_BND;
00079 }
00080 }
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097 int *newpid;
00098 ierr = PetscMalloc(gdata->n_ele*sizeof(int),&newpid);CHKERRQ(ierr);
00099
00100 int cnt1;
00101 cnt1=0;
00102
00103
00104 for (int i=0;i<gdata->n_ele;i++) {
00105 if (elebnd[i]==C_BND) {
00106
00107
00108
00109
00110
00111
00112 if (
00113 gdata->vertxyz[ND*gdata->elevert[NV*i+0]+2]+
00114 gdata->vertxyz[ND*gdata->elevert[NV*i+1]+2]+
00115 gdata->vertxyz[ND*gdata->elevert[NV*i+2]+2]+
00116 gdata->vertxyz[ND*gdata->elevert[NV*i+3]+2]<3
00117 ) {
00118 newpid[i]=2*gdata->n_prop;
00119 cnt1++;
00120 }
00121
00122 else if (
00123 gdata->vertxyz[ND*gdata->elevert[NV*i+0]+2]+
00124 gdata->vertxyz[ND*gdata->elevert[NV*i+1]+2]+
00125 gdata->vertxyz[ND*gdata->elevert[NV*i+2]+2]+
00126 gdata->vertxyz[ND*gdata->elevert[NV*i+3]+2]>77
00127 ) {
00128 newpid[i]=2*gdata->n_prop+1;
00129 cnt1++;
00130 }
00131
00132 else {
00133 newpid[i]=gdata->eleprop[i]+gdata->n_prop;
00134 cnt1++;
00135 }
00136 }
00137 else {
00138
00139 newpid[i]=gdata->eleprop[i];
00140 }
00141 }
00142
00143
00144 for (int i=0;i<gdata->n_ele;i++) {
00145 if (newpid[i]>=gdata->n_prop) {
00146 gdata->eleprop[i]=newpid[i];
00147 newpid[i]=C_UNK;
00148 }
00149 }
00150 PetscPrintf(PETSC_COMM_WORLD,
00151 "Modified property id of %i surface elements\n",
00152 cnt1
00153 );
00154
00155
00156 int cnt2;
00157 cnt2=0;
00158
00159
00160 for (int j=1;j<nsliceprop;j++) {
00161 for (int i=0;i<gdata->n_ele;i++) {
00162 if (gdata->eleprop[i]>=gdata->n_prop) {
00163
00164 int t_nnb;
00165 t_nnb=ia[i+1]-ia[i];
00166
00167
00168 for (int k=0;k<t_nnb;k++) {
00169
00170 int t_nbid;
00171 t_nbid=ja[ia[i]+k];
00172
00173
00174 if (gdata->eleprop[t_nbid]<gdata->n_prop) {
00175 newpid[t_nbid]=gdata->eleprop[i];
00176 }
00177 }
00178 }
00179 }
00180
00181
00182 for (int i=0;i<gdata->n_ele;i++) {
00183 if (newpid[i]>=gdata->n_prop) {
00184 gdata->eleprop[i]=newpid[i];
00185 cnt2++;
00186 newpid[i]=C_UNK;
00187 }
00188 }
00189 }
00190
00191 PetscPrintf(PETSC_COMM_WORLD,
00192 "Modified property id of %i additional elements\n",
00193 cnt2
00194 );
00195
00196 PetscPrintf(PETSC_COMM_WORLD,
00197 "Modified property id of %i elements in total\n",
00198 cnt1+cnt2
00199 );
00200
00201 assert(cnt1+cnt2<=gdata->n_ele);
00202
00203
00204
00205 gdata->n_prop = 2*gdata->n_prop+2;
00206
00207
00208 ierr = PetscFree(gdata->propdat);CHKERRQ(ierr);
00209 ierr = ReadKrn(gdata);CHKERRQ(ierr);
00210
00211 ierr = PetscFree(newpid);CHKERRQ(ierr);
00212 ierr = PetscFree(elebnd);CHKERRQ(ierr);
00213
00214 ierr = PetscFree(ia);CHKERRQ(ierr);
00215 ierr = PetscFree(ja);CHKERRQ(ierr);
00216
00217 MagparFunctionLogReturn(0);
00218 }
00219