[petsc-dev] broken DMGetInterpolation_DA !!!!

[Barry Smith]

On Jul 21, 2011, at 5:00 PM, Jed Brown wrote:

> On Thu, Jul 21, 2011 at 14:48, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
> 
>  Whoever wrote the DMGetInterpolation_DA_Q1() for when coordinates are provided in DA
> 
>   1) it is completely broken for the periodic case!
> 
> There is a serious semantic problem with periodic non-regular coordinates: we don't know how big the wrap element is. There is no way to determine this automatically, and there is no way to store it.

   So what? You don't use the coordinate information in computing the interpolation anyways so it seems to me you can do the same thing for the ghost elements you do for all the other elements in computing the interpolation matrix. 

   Barry

>  
> 
>   2) it doesn't actually use the coordinates provided by the user in the DA ; it uses uniform coordinates it cooks up on the fly
> 
> /* compute local coordinate arrays */
>    nxi  = ratioi + 1;
>    neta = ratioj + 1;
>    ierr = PetscMalloc(sizeof(PetscScalar)*nxi,&xi);CHKERRQ(ierr);
>    ierr = PetscMalloc(sizeof(PetscScalar)*neta,&eta);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));
>    }
> 
> These are reference coordinates. They get mapped using a Q1 local coordinate mapping. Dave wrote this code, but it was tested to give second order accuracy for non-uniform mapped grids.
> 
> Are you using non-uniform periodic grids? If you have regular periodic grids, but you still want to store coordinates explicitly, then we could write some special case to detect that you set coordinates that happened to be uniform, therefore we can guess that the wrap cell should also be uniform, and the code for no provided coordinates would be correct. What do you want to do?