Jed,
You will need to explain to me all this code or I am going to have to rip
it all out since Matt wants a release.
That is, all the DMDA interpolation stuff that was carefully put in under
your supervision will be removed unless I understand why it is there. I asked
you once before and you didn't really answer.
For example in the 3d interpolation it has a
if (!vcoors) {
then two blocks of code NEITHER of which actually use the coordinate
information AT ALL. So why the two blocks of code? Is one wrong and one right
and if so why not just get rid of the wrong one. BTW the second version does
not work for periodic conditions which is bad! Is your intention to use the
coordinate information on each level to generate an interpolation matrix that
depends on the coordinates or not? Why is the vcoors flag used as a check here?
Are you assuming that the vcoors existing has some special meaning? I
understand that the code
Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
comes from the trilinear finite element basic functions to define an
interpolation but since you are evaluating them on a uniform mesh
for (li=0; li<nxi; li++) {
xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
}
for (lj=0; lj<neta; lj++) {
eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
}
for (lk=0; lk<nzeta; lk++) {
zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1));
}
why bother? Isn't this always going to give you a uniform grid interpolation?
Also src/dm/examples/tests/ex36.c seems to be suppose to test the interpolation
for nonuniform grids but it doesn't seem to work
arry-smiths-macbook-pro:tests barrysmith$ ./ex36 -dim 2 -nl 1 -cmap 1
[2 x 2]=>[4 x 4], interp err = 1.4873e+00
and doesn't have any "nightly" output files so isn't a test at all. Is there a
test that is actually run in the nightly for this?
I have been totally confused about this stuff for a long time and it needs
to be resolved.
Barry
ierr = DMDAGetCoordinates(daf,&vcoors);CHKERRQ(ierr);
if (vcoors) {
ierr = DMDAGetGhostedCoordinates(dac,&cvcoors);CHKERRQ(ierr);
ierr = DMDAVecGetArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr);
ierr = DMDAVecGetArray(ddc->da_coordinates,cvcoors,&ccoors);CHKERRQ(ierr);
}
/* loop over local fine grid nodes setting interpolation for those*/
if (!vcoors) {
for (l=l_start; l<l_start+p_f; l++) {
for (j=j_start; j<j_start+n_f; j++) {
for (i=i_start; i<i_start+m_f; i++) {
/* convert to local "natural" numbering and then to PETSc global
numbering */
row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) +
m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
i_c = (i/ratioi);
j_c = (j/ratioj);
l_c = (l/ratiok);
/*
Only include those interpolation points that are truly
nonzero. Note this is very important for final grid lines
in x and y directions; since they have no right/top neighbors
*/
x = ((double)(i - i_c*ratioi))/((double)ratioi);
y = ((double)(j - j_c*ratioj))/((double)ratioj);
z = ((double)(l - l_c*ratiok))/((double)ratiok);
nc = 0;
/* one left and below; or we are right on it */
col =
dof*(m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c)+m_ghost_c*(j_c-j_start_ghost_c)+(i_c-i_start_ghost_c));
cols[nc] = idx_c[col]/dof;
v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
if (i_c*ratioi != i) {
cols[nc] = idx_c[col+dof]/dof;
v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. -
(2.0*z-1.));
}
if (j_c*ratioj != j) {
cols[nc] = idx_c[col+m_ghost_c*dof]/dof;
v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. -
(2.0*z-1.));
}
if (l_c*ratiok != l) {
cols[nc] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof;
v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. +
(2.0*z-1.));
}
if (j_c*ratioj != j && i_c*ratioi != i) {
cols[nc] = idx_c[col+(m_ghost_c+1)*dof]/dof;
v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. -
(2.0*z-1.));
}
if (j_c*ratioj != j && l_c*ratiok != l) {
cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof;
v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. +
(2.0*z-1.));
}
if (i_c*ratioi != i && l_c*ratiok != l) {
cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof;
v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. +
(2.0*z-1.));
}
if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof;
v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. +
(2.0*z-1.));
}
ierr =
MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr);
}
}
}
} else {
PetscScalar *xi,*eta,*zeta;
PetscInt li,nxi,lj,neta,lk,nzeta,n;
PetscScalar Ni[8];
/* compute local coordinate arrays */
nxi = ratioi + 1;
neta = ratioj + 1;
nzeta = ratiok + 1;
ierr = PetscMalloc(sizeof(PetscScalar)*nxi,&xi);CHKERRQ(ierr);
ierr = PetscMalloc(sizeof(PetscScalar)*neta,&eta);CHKERRQ(ierr);
ierr = PetscMalloc(sizeof(PetscScalar)*nzeta,&zeta);CHKERRQ(ierr);
for (li=0; li<nxi; li++) {
xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
}
for (lj=0; lj<neta; lj++) {
eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
}
for (lk=0; lk<nzeta; lk++) {
zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1));
}
for (l=l_start; l<l_start+p_f; l++) {
for (j=j_start; j<j_start+n_f; j++) {
for (i=i_start; i<i_start+m_f; i++) {
/* convert to local "natural" numbering and then to PETSc global
numbering */
row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) +
m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
i_c = (i/ratioi);
j_c = (j/ratioj);
l_c = (l/ratiok);
/* remainders */
li = i - ratioi * (i/ratioi);
if (i==mx-1){ li = nxi-1; }
lj = j - ratioj * (j/ratioj);
if (j==my-1){ lj = neta-1; }
lk = l - ratiok * (l/ratiok);
if (l==mz-1){ lk = nzeta-1; }
/* corners */
col =
dof*(m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c)+m_ghost_c*(j_c-j_start_ghost_c)+(i_c-i_start_ghost_c));
cols[0] = idx_c[col]/dof;
Ni[0] = 1.0;
if ( (li==0) || (li==nxi-1) ) {
if ( (lj==0) || (lj==neta-1) ) {
if ( (lk==0) || (lk==nzeta-1) ) {
ierr =
MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);CHKERRQ(ierr);
continue;
}
}
}
/* edges + interior */
/* remainders */
if (i==mx-1){ i_c--; }
if (j==my-1){ j_c--; }
if (l==mz-1){ l_c--; }
col = dof*(m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) +
m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
cols[0] = idx_c[col]/dof; /* one left and below; or we are right on
it */
cols[1] = idx_c[col+dof]/dof; /* one right and below */
cols[2] = idx_c[col+m_ghost_c*dof]/dof; /* one left and above */
cols[3] = idx_c[col+(m_ghost_c+1)*dof]/dof; /* one right and above */
cols[4] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof; /* one left and
below and front; or we are right on it */
cols[5] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof; /* one right
and below, and front */
cols[6] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof;/* one
left and above and front*/
cols[7] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof; /*
one right and above and front */
Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
for (n=0; n<8; n++) {
if( PetscAbsScalar(Ni[n])<1.0e-32) { cols[n]=-1; }
}
ierr =
MatSetValues(mat,1,&row,8,cols,Ni,INSERT_VALUES);CHKERRQ(ierr);
}
}
}
ierr = PetscFree(xi);CHKERRQ(ierr);
ierr = PetscFree(eta);CHKERRQ(ierr);
ierr = PetscFree(zeta);CHKERRQ(ierr);
}
