<html><head><style type='text/css'>p { margin: 0; }</style></head><body><div style='font-family: Times New Roman; font-size: 12pt; color: #000000'>Hi Paul,<br><br>I am looking through the routine facind, and am wondering: Does this routine provide the loop limits for a particular face? Meaning, if I want to look at the points on face 1, the values it returns are what I loop through for r, s, and t?<br><br>Also, I can't seem to find the subroutine for facind2 to have a look at it. I have found many calls for it though...where is it located?<br><br>- Michael<br><br>----- Original Message -----<br>From: nek5000-users@lists.mcs.anl.gov<br>To: nek5000-users@lists.mcs.anl.gov<br>Sent: Tuesday, July 13, 2010 4:38:28 AM GMT -06:00 Guadalajara / Mexico City / Monterrey<br>Subject: Re: [Nek5000-users] r,s,t index for given point<br><br>
Hi Michael,
Points in the volume are stored in lexicographical order:
real u(lx1,ly1,lz1)
implies that u(i+1,j,k) follows u(i,j,k) in memory, so the
data layout (in memory) is:
u(1,1,1) u(2,1,1)...u(lx1,1,1) u(1,2,1) u(2,2,1)... u(lx1,ly1,lz1)
For u(i,j,k) the i index is associated with r, j with s, and k with t.
There are two numberings of the faces. The one that you're
most likely to interact with, and which is used for assigning
bcs, unit normals, surface Jacobians, etc. is what we refer to
as the "preprocessor notation". In 2D, preprocessor notation
uses the counter-clockwise ordering:
^ s
| Face 3
+-----------+
| | |
| | |
Face 4 | +-----|--> r
| |
| | Face 2
+-----------+
Face 1
In 3D, it's the same, save that you also have Face 5 assocated with
t=-1 and Face 6 associated with t=+1.
When you traverse a face, the indices describing that 2D manifold
advance in positive order and in that case the associated face
objects (e.g., unit normals and area) will have unit stride.
There are two ways to stride across a face f=1,...,6 (and there
are many routines implementing these).
Here is one:
integer e,f
call facind (i0,i1,j0,j1,k0,k1,f)
ia = 0
flux = 0
do k=k0,k1
do j=j0,j1
do i=i0,i1
ia = ia + 1
flux = flux + ( vx(i,j,k,e)*unx(ia,1,f,e)
$ vy(i,j,k,e)*uny(ia,1,f,e)
$ vz(i,j,k,e)*unz(ia,1,f,e) )*area(ia,1,f,e)
enddo
enddo
enddo
Here is another:
call facind2 (js1,jf1,jskip1,js2,jf2,jskip2,f)
ia = 0
flux = 0
do j2=js2,jf2,jskip2
do j1=js1,jf1,jskip1
ia = ia+1
flux = flux + ( vx(j1,j2,1,e)*unx(ia,1,f,e)
$ + vy(j1,j2,1,e)*uny(ia,1,f,e)
$ + vz(j1,j2,1,e)*unz(ia,1,f,e) )*area(ia,1,f,e)
enddo
enddo
I like the second approach because it is faster.
Note that unx,uny,unz are _face_ based structures, while vx
vy vz are volumetric.
I think if you study these loops for a bit you should be able
to see the connections between the two data layouts.
Paul
On Mon, 12 Jul 2010, nek5000-users@lists.mcs.anl.gov wrote:
Hi Paul,
To expand on this topic, I am wondering what the convention being used
to number the GL points for an element is? Could you explain this?
Thinking of it terms of the r-s-t diagram for a nekton element, I
know from previous examples that face 5 contains points 1 to
lx1*ly1. Also, face 6 should contain points (lz1-1)*lx1*ly1 to
lx1*ly1*lz1.
What I am trying to figure out, is what points are on the remaining
faces, and understanding the convention I think would be useful.
Thanks for the help!
>
> - Michael
>
>
> ----- Original Message -----
> From: nek5000-users@lists.mcs.anl.gov
> To: nek5000-users@lists.mcs.anl.gov
> Sent: Monday, July 12, 2010 5:13:47 AM GMT -06:00 US/Canada Central
> Subject: Re: [Nek5000-users] r,s,t index for given point
>
> fixed in r548.
> Stefan
>
> On Jul 12, 2010, at 10:49 AM, nek5000-users@lists.mcs.anl.gov wrote:
>
>>
>> Hi Markus,
>>
>> lglel reads "local to global element number" so input is
>> local element number and output is global element number.
>>
>> I'll be seeing Stefan this week so I'm sure we can straighten
>> out the interpolation routines.
>>
>> You're right about the other points - the ntotv wasn't flagged
>> by my compiler and the .his file was taken care of by my own
>> scripts, so I didn't catch either of them. Sorry about that.
>> I've updated the source so that those issues should be
>> resolved.
>>
>> Paul
>>
>>
>> On Mon, 12 Jul 2010, nek5000-users@lists.mcs.anl.gov wrote:
>>
>>> Hi,
>>>
>>> that was not it, I'm afraid (though I thought gllel transforms from local to global).
>>>
>>> But I also realized that the hpts routine does not deliver the correct values for me. To test this, I set up a simple 3d box with inflow, outflow and symmetry all around, Re very low, initial condition=boundary condition and 4 points.
>>> Please find all files attached, the output for the first point is OK, but then the values seem to shift positions and don't make sense.
>>>
>>> Also, when running it for 1 step or more, I get the error message:
>>> "
>>> dump history points
>>> done :: dump history points
>>> hisfile:
>>> /Users/m0s1978/nektonruns/rectest/rectest.his
>>> ERROR: .sch file already exists. ierr= 1
>>> "
>>>
>>> And lastly, the newest repo won't compile, here is the error message:
>>> "
>>> Error: /Users/m0s1978/nek5trial/trunk/nek/drive2.f, line 1551: This name has already been assigned a data type. [NVTOT]
>>> integer*8 ntot,nvtot,nptot
>>> ---------------------^
>>> "
>>>
>>> Thanks for any help,
>>> Markus
>>>
>>>
>>> nek5000-users@lists.mcs.anl.gov wrote:
>>>> Hi Markus,
>>>> It's possible that it's returning the local (to the processor)
>>>> element number. You would get the global (to the computational
>>>> domain, defined by the .rea file) element number via:
>>>>
>>>> integer e,eg
>>>>
>>>> eg = lglel(e)
>>>> Paul
>>>> On Sun, 11 Jul 2010, nek5000-users@lists.mcs.anl.gov wrote:
>>>>> Hi again,
>>>>> Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call:
>>>>>
>>>>> call findpts(inth_hpts,rcode,1,
>>>>> & proc,1,
>>>>> & elid,1,
>>>>> & rst,ndim,
>>>>> & dist,1,
>>>>> & pts(1,1),ndim,
>>>>> & pts(2,1),ndim,
>>>>> & pts(3,1),ndim,npoints)
>>>>> does not return the same element number (elid) as it can be found in the .rea file, for instance.
>>>>> For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point.
>>>>> Is there a way to transform the function return or am I missing something else?
>>>>> Thanks,
>>>>> Markus
>>>>> _______________________________________________
>>>>> Nek5000-users mailing list
>>>>> Nek5000-users@lists.mcs.anl.gov
>>>>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>>>> _______________________________________________
>>>> Nek5000-users mailing list
>>>> Nek5000-users@lists.mcs.anl.gov
>>>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>>>
>> _______________________________________________
>> Nek5000-users mailing list
>> Nek5000-users@lists.mcs.anl.gov
>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users@lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
><br>_______________________________________________
Nek5000-users mailing list
Nek5000-users@lists.mcs.anl.gov
https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
</div></body></html>