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 maybe generate an error message??? Producing the incorrect answer is not
a good approach.
>
>
> 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.
I completely do not understand the situation. In the code
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);
}
but the variable coors and ccoors are NEVER used in the routine! How can this
routine supposedly support non-uniform mapped grids when the coordinates are
never even taken into account? And why is
ierr = DMDAVecGetArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr);
ierr = DMDAVecGetArray(ddc->da_coordinates,cvcoors,&ccoors);CHKERRQ(ierr);
even called if the result is not used?
What am I missing?
>
> Are you using non-uniform periodic grids?
No, they are uniform. But I want to store the coordinates explicitly.
> 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?
