modifyprop_ser.c

Go to the documentation of this file.
00001 /*
00002     This file is part of magpar.
00003 
00004     Copyright (C) 2002-2009 Werner Scholz
00005 
00006     www:   http://www.magpar.net/
00007     email: magpar(at)magpar.net
00008 
00009     magpar is free software; you can redistribute it and/or modify
00010     it under the terms of the GNU General Public License as published by
00011     the Free Software Foundation; either version 2 of the License, or
00012     (at your option) any later version.
00013 
00014     magpar is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00017     GNU General Public License for more details.
00018 
00019     You should have received a copy of the GNU General Public License
00020     along with magpar; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 */
00023 
00024 static char const Id[] = "$Id: modifyprop_ser.c 2681 2009-07-31 04:30:53Z 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     /* set default values if option is missing */
00048     nsliceprop=0;
00049     PetscPrintf(PETSC_COMM_WORLD,
00050       "Option -nslicepropser not found, using default value %i!\n",
00051       nsliceprop
00052     );
00053   }
00054 
00055   /* exit if there is nothing to modify */
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   /* mark element if a whole face is on the boundary */
00072   for (int i=0;i<gdata->n_ele;i++) {
00073     elebnd[i]=C_INT;
00074     /* get number of neighbours for current element */
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   /* mark element if at least a node is on the boundary */
00083 /*
00084   for (int i=0;i<gdata->n_ele;i++) {
00085     elebnd[i]=C_INT;
00086     for (int j=0;j<NV;j++) {
00087       if (gdata->vertbndg2bnd[gdata->elevert[i*NV+j]]>=0) {
00088         if (elebnd[i]!=C_BND) {
00089           elebnd[i]=C_BND;
00090         }
00091       }
00092     }
00093   }
00094 */
00095 
00096   /* allocate array of new property ids */
00097   int *newpid;
00098   ierr = PetscMalloc(gdata->n_ele*sizeof(int),&newpid);CHKERRQ(ierr);
00099 
00100   int cnt1;
00101   cnt1=0;
00102 
00103   /* elements on surface get new (highest) property id */
00104   for (int i=0;i<gdata->n_ele;i++) {
00105     if (elebnd[i]==C_BND) {
00106 /*
00107       newpid[i]=gdata->eleprop[i]+gdata->n_prop;
00108       cnt1++;
00109 */
00110 
00111       /* special case: bottom surface */
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       /* special case: top surface */
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       /* assign new property id to other surface elements */
00132       else {
00133         newpid[i]=gdata->eleprop[i]+gdata->n_prop;
00134         cnt1++;
00135       }
00136     }
00137     else {
00138       /* retain old property id */
00139       newpid[i]=gdata->eleprop[i];
00140     }
00141   }
00142 
00143   /* modify property id of selected elements */
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   /* propagate new property ids into the interior */
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         /* get number of neighbours for current element */
00164         int t_nnb;
00165         t_nnb=ia[i+1]-ia[i];
00166 
00167         /* loop all neighbours */
00168         for (int k=0;k<t_nnb;k++) {
00169           /* get the id of our neighbour */
00170           int t_nbid;
00171           t_nbid=ja[ia[i]+k];
00172 
00173           /* set new property id unless we do already have a new one! */
00174           if (gdata->eleprop[t_nbid]<gdata->n_prop) {
00175             newpid[t_nbid]=gdata->eleprop[i];
00176           }
00177         }
00178       }
00179     }
00180 
00181     /* modify property id of selected elements */
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   /* update counter of property ids */
00205   gdata->n_prop = 2*gdata->n_prop+2;
00206 
00207   /* reread property ids */
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 

magpar - Parallel Finite Element Micromagnetics Package
Copyright (C) 2002-2009 Werner Scholz