[petsc-users] related to DMPLExvecsetclsoure

[Dharmendar Reddy]

On Mon, Feb 3, 2014 at 5:49 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Mon, Feb 3, 2014 at 5:29 PM, Dharmendar Reddy <dharmareddy84 at gmail.com>
> wrote:
>>
>> Hello,
>>         I am trying to understand the layout of the values array
>> argument of the DMPLexVecSetClosure.
>>
>> Consider a triangular element with two fields at each node, say n, p
>> vertices : v1,v2,v3
>> If the value of fields in the local vector :  n1,p1, n2, p2, n3, p3
>>
>> If the equation set is: F(n,p) = 0 and G(n,p) = 0
>>
>> if the values of F and G evaluated at v1,v2, v3 are F1,F2,F3 and
>> G1,G2G3 respectively.
>>
>> Then the elemVec that goes into DMPLexVecSetClosure is:
>>
>> elemVec = { F1,F2,F3,G1G2,G3} is that correct ?
>>
>> Now, i am looking for some functionality which will let me access the
>> following mapping information, (local dof id)?
>>
>> what is the Local Dof Id , i.e, the position in elemVec given the
>> field id, component id and basis id (or node id for the case of P1
>> Lagrange ?)
>>
>>
>> The usage will be like this:
>
>
> If, as in this case, the fields have the same element, its just normal
> array ordering:
>
>   (f, b, c) = (f*numBasis + b)*numComponents + c

I used this approach in my current code which has same  number of
fields in all regions of the mesh.
I pre calculate and store the local dof ids in an array.

But now i am trying do a case where number of fields is not same in
all regions. I am a bit confused about the interface case.

Consider for example of interface region
 1----------2
 |   c1      |
 |             |
 3---------4
 |   c2      |
 |             |
 5---------6
the variables defined at nodes is shown below:
1: (phi),  2:(phi), [3,4,5,6]: (phi,n,p)

I think, c1 will have dof layout as :dofArray:
(phi3,phi4,phi2,phi1,n3,n4,p3,p4) assuming anti clockwise ordering of
nodes in ref cell. I can get dofArray using vecget closure on c1.



>
> In general, all field values are stacked together, and ordered by
> basis vec and then component.
>
>     Matt
>
>>
>> do fid=1,numField
>>    numBasis = getNumBasis(fid)
>>    numComponent = getNumComponent(fid)
>>    do bid=1,numBasis
>>       do cid =1,numCompoent
>>        localId = getLocalDofId( fid, bid, cid)
>>         elemResVec(localId) = residual value
>>      end do
>>    end do
>> enddo
>>
>> call
>> DMPlexVecSetClosure(dm,section,localResVec,cellId,elemResVec,ADD_VALUES,ierr)
>>
>> Thanks
>> Reddy
>
>
>
>
> --
> What most experimenters take for granted before they begin their experiments
> is infinitely more interesting than any results to which their experiments
> lead.
> -- Norbert Wiener