filternodes.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: filternodes.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 
00029 #define C_REQ -5   
00031 int FilterNodes(int *pn_vert, int *pn_ele, PetscReal **pvertxyz, int *elevert, Vec *pM, Vec *pVH1, Vec *pVH2)
00032 {
00033   Vec          M,VH1,VH2;
00034   PetscReal    *toM,*toVH1,*toVH2;
00035   PetscReal    *tnM,*tnVH1,*tnVH2;
00036 
00037   MagparFunctionLogBegin;
00038 
00039   int rank,size;
00040   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00041   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00042 
00043   if (rank) {
00044     int n_del;
00045     ierr = MPI_Bcast(&n_del,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00046 
00047     if (n_del!=0) {
00048       ierr = VecCreate(PETSC_COMM_WORLD,&M);CHKERRQ(ierr);
00049       ierr = VecSetSizes(M,0,PETSC_DECIDE);CHKERRQ(ierr);
00050       ierr = VecSetFromOptions(M);CHKERRQ(ierr);
00051       ierr = VecSetBlockSize(M,ND);CHKERRQ(ierr);
00052       ierr = VecDestroy(*pM);CHKERRQ(ierr);
00053       *pM=M;
00054 
00055       if (*pVH1!=PETSC_NULL) {
00056         ierr = VecDuplicate(M,&VH1);CHKERRQ(ierr);
00057         ierr = VecDestroy(*pVH1);CHKERRQ(ierr);
00058         *pVH1=VH1;
00059       }
00060       if (*pVH2!=PETSC_NULL) {
00061         ierr = VecDuplicate(M,&VH2);CHKERRQ(ierr);
00062         ierr = VecDestroy(*pVH2);CHKERRQ(ierr);
00063         *pVH2=VH2;
00064       }
00065     }
00066 
00067     MagparFunctionLogReturn(0);
00068   }
00069 
00070   int n_vert;         n_vert=*pn_vert;
00071   int n_ele;          n_ele=*pn_ele;
00072   PetscReal *vertxyz; vertxyz=*pvertxyz;
00073 
00074   int *nodesreq;
00075   ierr = PetscMalloc(n_vert*sizeof(int),&nodesreq);CHKERRQ(ierr);
00076   for (int i=0;i<n_vert;i++) nodesreq[i]=C_UNK;
00077 
00078   /* find all required nodes */
00079   for (int i=0;i<n_ele;i++) {
00080     for (int j=0;j<NV;j++) {
00081       nodesreq[elevert[i*NV+j]]=C_REQ;
00082     }
00083   }
00084 
00085   int n_req=0;
00086 
00087   /* create mapping */
00088   for (int i=0;i<n_vert;i++) {
00089     if (nodesreq[i]==C_REQ) {
00090       /* assign new id and update counter */
00091       nodesreq[i]=n_req++;
00092     }
00093   }
00094 
00095   /* get number of nodes to be deleted and inform other processors */
00096   int n_del;
00097   n_del=n_vert-n_req;
00098   ierr = MPI_Bcast(&n_del,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00099 
00100   if (n_del!=0) {
00101     PetscPrintf(PETSC_COMM_WORLD,
00102     "Unassigned vertices: %i out of %i\n",n_del,n_vert);
00103   }
00104   else {
00105     PetscPrintf(PETSC_COMM_WORLD,"No unassigned vertices\n");
00106     ierr = PetscFree(nodesreq);CHKERRQ(ierr);
00107     MagparFunctionLogReturn(0);
00108   }
00109 
00110   /* renumber nodes */
00111   for (int i=0;i<n_ele;i++) {
00112     for (int j=0;j<NV;j++) {
00113       elevert[i*NV+j]=nodesreq[elevert[i*NV+j]];
00114       assert(nodesreq[elevert[i*NV+j]]<n_req);
00115     }
00116     assert(
00117       elevert[i*NV+0] !=
00118       elevert[i*NV+1] &&
00119       elevert[i*NV+0] !=
00120       elevert[i*NV+2] &&
00121       elevert[i*NV+0] !=
00122       elevert[i*NV+3] &&
00123       elevert[i*NV+1] !=
00124       elevert[i*NV+2] &&
00125       elevert[i*NV+1] !=
00126       elevert[i*NV+3] &&
00127       elevert[i*NV+2] !=
00128       elevert[i*NV+3]
00129     );
00130   }
00131 
00132   PetscPrintf(PETSC_COMM_WORLD,"Updated elements\n");
00133 
00134   PetscReal *tvertxyz;
00135   ierr = PetscMalloc(ND*n_req*sizeof(PetscReal),&tvertxyz);CHKERRQ(ierr);
00136 
00137   ierr = VecCreate(PETSC_COMM_WORLD,&M);CHKERRQ(ierr);
00138   ierr = VecSetSizes(M,ND*n_req,PETSC_DECIDE);CHKERRQ(ierr);
00139   ierr = VecSetFromOptions(M);CHKERRQ(ierr);
00140   ierr = VecSetBlockSize(M,ND);CHKERRQ(ierr);
00141 
00142   ierr = VecGetArray(*pM,&toM);CHKERRQ(ierr);
00143   ierr = VecGetArray(M,&tnM);CHKERRQ(ierr);
00144 
00145   if (*pVH1!=PETSC_NULL) {
00146     ierr = VecDuplicate(M,&VH1);CHKERRQ(ierr);
00147     ierr = VecGetArray(*pVH1,&toVH1);CHKERRQ(ierr);
00148     ierr = VecGetArray(VH1,&tnVH1);CHKERRQ(ierr);
00149   }
00150   if (*pVH2!=PETSC_NULL) {
00151     ierr = VecDuplicate(M,&VH2);CHKERRQ(ierr);
00152     ierr = VecGetArray(*pVH2,&toVH2);CHKERRQ(ierr);
00153     ierr = VecGetArray(VH2,&tnVH2);CHKERRQ(ierr);
00154   }
00155 
00156   int k=0;
00157   /* update (and compress) vertxyz */
00158   for (int i=0;i<n_vert;i++) {
00159     if (nodesreq[i]!=C_UNK) {
00160       for (int j=0;j<ND;j++) {
00161         tvertxyz[ND*k+j]=vertxyz[ND*i+j];
00162         tnM[ND*k+j]=toM[ND*i+j];
00163         if (*pVH1!=PETSC_NULL) tnVH1[ND*k+j]=toVH1[ND*i+j];
00164         if (*pVH2!=PETSC_NULL) tnVH2[ND*k+j]=toVH2[ND*i+j];
00165       }
00166       k++;
00167     }
00168   }
00169   assert(k==n_req);
00170 
00171   ierr = VecRestoreArray(*pM,&toM);CHKERRQ(ierr);
00172   ierr = VecRestoreArray(M,&tnM);CHKERRQ(ierr);
00173   ierr = VecDestroy(*pM);CHKERRQ(ierr);
00174   *pM=M;
00175 
00176   if (*pVH1!=PETSC_NULL) {
00177     ierr = VecRestoreArray(*pVH1,&toVH1);CHKERRQ(ierr);
00178     ierr = VecRestoreArray(VH1,&tnVH1);CHKERRQ(ierr);
00179     ierr = VecDestroy(*pVH1);CHKERRQ(ierr);
00180     *pVH1=VH1;
00181   }
00182   if (*pVH2!=PETSC_NULL) {
00183     ierr = VecRestoreArray(*pVH2,&toVH2);CHKERRQ(ierr);
00184     ierr = VecRestoreArray(VH2,&tnVH2);CHKERRQ(ierr);
00185     ierr = VecDestroy(*pVH2);CHKERRQ(ierr);
00186     *pVH2=VH2;
00187   }
00188 
00189   PetscPrintf(PETSC_COMM_WORLD,
00190     "Updated Cartesian coordinates of %i vertices\n",
00191     n_req
00192   );
00193 
00194   /* update counters */
00195   *pn_vert=n_req;
00196   *pn_ele=n_ele;
00197 
00198   ierr = PetscFree(*pvertxyz);CHKERRQ(ierr);
00199   *pvertxyz=tvertxyz;
00200 
00201   ierr = PetscFree(nodesreq);CHKERRQ(ierr);
00202 
00203   MagparFunctionLogReturn(0);
00204 }

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