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: 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
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
00088 for (int i=0;i<n_vert;i++) {
00089 if (nodesreq[i]==C_REQ) {
00090
00091 nodesreq[i]=n_req++;
00092 }
00093 }
00094
00095
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
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
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
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 }