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: calAfe2fe.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 "util.h"
00028
00029 int CalAfe2fe(GridData *gdata,PetscReal *shift,int *slice_id,Mat *pAIpol)
00030 {
00031 MagparFunctionLogBegin;
00032
00033 PetscPrintf(PETSC_COMM_WORLD,"Calculating interpolation matrix fe->fe\n");
00034
00035 Mat AIpol;
00036 ierr = MatCreateMPIAIJ(
00037 PETSC_COMM_WORLD,
00038 ND*gdata->ln_vert,ND*gdata->ln_vert,
00039 ND*gdata->n_vert,ND*gdata->n_vert,
00040 NV,PETSC_NULL,NV,PETSC_NULL,
00041 &AIpol
00042 );CHKERRQ(ierr);
00043 ierr = MatSetFromOptions(AIpol);CHKERRQ(ierr);
00044
00045
00046
00047
00048
00049
00050 for (int i=0;i<gdata->ln_vert;i++) {
00051 for (int k=0;k<ND;k++) {
00052 ierr = MatSetValue(
00053 AIpol,
00054 ND*gdata->vertl2g[i]+k,
00055 ND*gdata->vertl2g[i]+k,
00056 1.0,
00057 INSERT_VALUES
00058 );CHKERRQ(ierr);
00059 }
00060 }
00061
00062 #define DIST 1e-12
00063 PetscReal tshift[ND];
00064 my_dcopy(ND,shift,1,tshift,1);
00065 my_dscal(ND,1+DIST,tshift,1);
00066 tshift[0]+=DIST;
00067 tshift[1]+=DIST;
00068 tshift[2]+=DIST;
00069 PetscPrintf(PETSC_COMM_WORLD,"Shifting nodes by (%g,%g,%g)\n",tshift[0],tshift[1],tshift[2]);
00070
00071 PetscPrintf(PETSC_COMM_WORLD,"Calculating interpolation matrix...\n");
00072 for (int i=0;i<gdata->ln_ele;i++) {
00073 ierr = ProgressBar(i,gdata->ln_ele,10);
00074
00075 int j;
00076 for (j=0;j<gdata->n_prop;j++) {
00077 if (gdata->eleprop[i]==slice_id[j]) break;
00078 }
00079 if (j==gdata->n_prop) continue;
00080
00081 for (int k=0;k<NV;k++) {
00082 int v;
00083 v=gdata->elevert[NV*i+k];
00084
00085 PetscReal p[ND],rhs[ND];
00086
00087
00088 for (int m=0;m<ND;m++) {
00089 p[m]=gdata->vertxyz[ND*v+m];
00090 p[m]+=tshift[m];
00091 }
00092
00093 if (
00094 barycent(
00095 p,
00096 gdata->vertxyz+ND*gdata->elevert[NV*i+0],
00097 gdata->vertxyz+ND*gdata->elevert[NV*i+1],
00098 gdata->vertxyz+ND*gdata->elevert[NV*i+2],
00099 gdata->vertxyz+ND*gdata->elevert[NV*i+3],
00100 rhs
00101 )
00102 ) {
00103 for (int l=0;l<ND;l++) {
00104 PetscReal matele;
00105
00106 matele=(1-rhs[0]-rhs[1]-rhs[2]);
00107 ierr = MatSetValue(AIpol,
00108 ND*v+l,
00109 ND*gdata->elevert[NV*i+3]+l,
00110 matele,INSERT_VALUES);CHKERRQ(ierr);
00111
00112 matele=rhs[0];
00113 ierr = MatSetValue(AIpol,
00114 ND*v+l,
00115 ND*gdata->elevert[NV*i+0]+l,
00116 matele,INSERT_VALUES);CHKERRQ(ierr);
00117
00118 matele=rhs[1];
00119 ierr = MatSetValue(AIpol,
00120 ND*v+l,
00121 ND*gdata->elevert[NV*i+1]+l,
00122 matele,INSERT_VALUES);CHKERRQ(ierr);
00123
00124 matele=rhs[2];
00125 ierr = MatSetValue(AIpol,
00126 ND*v+l,
00127 ND*gdata->elevert[NV*i+2]+l,
00128 matele,INSERT_VALUES);CHKERRQ(ierr);
00129 }
00130 }
00131 }
00132 }
00133
00134 ierr = ProgressBar(1,1,10);
00135
00136 PetscInfo(0,"AIpol matrix assembly complete\n");
00137 ierr = MatAssemblyBegin(AIpol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00138 ierr = MatAssemblyEnd(AIpol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00139
00140 #if 0
00141
00142
00143
00144
00145
00146 Vec Vdiag;
00147 ierr = VecDuplicate(gdata->M,&Vdiag);CHKERRQ(ierr);
00148
00149
00150 ierr = MatGetDiagonal(AIpol,Vdiag);CHKERRQ(ierr);
00151
00152 PetscReal *d;
00153 ierr = VecGetArray(Vdiag,&d);CHKERRQ(ierr);
00154 for (int i=0;i<ND*gdata->ln_vert;i++) {
00155 if (d[gdata->vertl2g[i]]==0.0) d[gdata->vertl2g[i]]=1.0;
00156 }
00157 ierr = VecRestoreArray(Vdiag,&d);CHKERRQ(ierr);
00158
00159 ierr = MatDiagonalSet(AIpol,Vdiag,INSERT_VALUES);CHKERRQ(ierr);
00160 #endif
00161
00162 ierr = PrintMatInfoAll(AIpol);CHKERRQ(ierr);
00163
00164
00165
00166
00167
00168
00169 *pAIpol=AIpol;
00170
00171 MagparFunctionLogReturn(0);
00172 }