[petsc-users] query DMPlexCreateSection

[Brad Aagaard]

On 01/10/2014 11:19 AM, Dharmendar Reddy wrote:
> Hello,
>          I was able to create the petsc section as per your advice.
>>    DMPlexSetChart()
>>    <loop over region 1>
>>      DMPlexSetDof() and DMPlexSetFieldDof()
>>    <loop over region 2>
>>      DMPlexSetDof() and DMPlexSetFieldDof()
>
> I could program the code upto this step.
>
>    ! Create Section
>    call PetscSectionCreate(comm, section, ierr)
>    ! This call sets the number of fields in the section. Default number
> of compents per field is 1
>    call PetscSectionSetNumFields(section, numField, ierr)
>    ! Set the number of componets per field
>    do fid=0,numField-1
>      call PetscSectionSetFieldComponents(section, fid, numComp(fid+1), ierr)
>    end do
>    ! Set the chart for section
>    call DMPlexGetChart(meshDM, cellIdStart, cellIdEnd, ierr)
>    call PetscSectionSetChart(section, cellIdStart, cellIdEnd, ierr)
>    print*,'Sucessfully set chart',cellIdStart,cellIdEnd
>    rgnLabel = "region"
>    call DMPlexGetLabelIdIS(meshDM,rgnLabel,regionIS,ierr)
>    call ISView(regionIS,PETSC_VIEWER_STDOUT_SELF,ierr)
>    call ISGetIndicesF90(regionIS, regionId, ierr)
>    print*, 'region', regionId
>    do jc=1,size(regionId) ! Loop Over region1
>      rgnId = regionId(jc)
>      call DMPlexGetStratumIS(meshDM, "region",rgnId, regionCellIS, ierr)
>      !print*,'Number of cells in region',rgnId
>      call ISGetSize(regionCellIS, numCell, ierr)
>      !print*,'Number of cells in region',rgnId, 'is',numCell
>      call ISGetIndicesF90(regionCellIS, regionCell, ierr)
>      do ic=1,numCell ! Get the cells of this region and loop over Cells
>        cellId = regionCell(ic)
>        ! Get the dimension of the cell
>        call DMPlexGetLabelValue(meshDM, "depth", cellId, cellDim, ierr)
>        do fid=0,numField-1 ! Loop over fields
>           ! Set field only if the field is defined for this region
>           !print*, 'settinf Field
> for',cellId,cellDim,fid,rgnId,numDof(fid+1,cellDim+1)
>           if(fieldInRegion(fid+1,rgnId)) then
>              call PetscSectionSetFieldDof(section, cellId, fid,
> numDof(fid+1,cellDim+1), ierr)
>           end if
>        end do
>        call  PetscSectionSetDof(section, cellId,
> numDofTot(cellDim+1,rgnId), ierr)
>      end do
>      call ISRestoreIndicesF90(regionCellIS, regionCell, ierr)
>    end do
>    call ISRestoreIndicesF90(regionIS, regionId, ierr)
>
>
> Now I need to work on the code for BC. Can you help me with a bit more
> detail here.
>>    <loop over BC>
>>      DMPlexSetConstraintDof() and DMPlexSetFieldConstraintDof()

Use DMPlexSetConstraintDof() to set the number of DOF constrained at 
each point on the BC. DMPlexSetFieldConstraintDof() is used if you have 
multiple fields in a section.

>>    DMPlexSetUp()
>>    <loop over BC>
>>      DMPlexSetConstraintIndices() and DMPlexSetFieldConstraintIndices()

Use DMPlexSetConstraintIndices() to set the indices of the DOF that are 
constrained.

For example, one might constrain the DOF normal and/or tangential to a 
boundary in an elasticity problem. You first set the size for all the 
points, allocate via DMPlexSetUp, and then set the indices for the 
constrained DOF at each point.

Regards,
Brad

> On Thu, Jan 2, 2014 at 9:56 AM, Matthew Knepley <knepley at gmail.com> wrote:
>> On Thu, Jan 2, 2014 at 4:11 AM, Dharmendar Reddy <dharmareddy84 at gmail.com>
>> wrote:
>>>
>>> Hello,
>>>           I am trying to use DMPlexCreateSection from fortran. I was
>>> able to solve a Poisson equation earlier using DMPlex for handling
>>> mesh data.
>>>
>>> Now i need to solve a system of equations: (poisson + continuity
>>> equations)
>>>
>>> - div (grad phi) = (C + n)  --- (1)
>>>    div (J ) = 0  ---------------  (2)
>>>    J = n grad(phi)
>>> ( C is constant)
>>> Simulation domain is defined as, rectangular regions 1 and 2 shown below.
>>>
>>> -------------
>>> |     1        |
>>> -------------
>>> |               |
>>> |       2      |
>>> |               |
>>> --------------
>>>
>>> Now, equation 1 is defined in region 1 and 2
>>> and equation 2 is defined only for region 2.
>>>
>>> How do i setup the section  ? DMPlexCreateSection applies the given
>>> DOF layout  per cell to all cells in the mesh.
>>>
>>> I am not sure if all the function calls inside
>>> DMPlexCreateSectionIntial and DMPlexCreateSectionBCDof  are accessible
>>> from Fotran.
>>
>>
>> You don't want them anyway since they apply to the whole domain. The control
>> flow could be:
>>
>>    DMPlexSetChart()
>>    <loop over region 1>
>>      DMPlexSetDof() and DMPlexSetFieldDof()
>>    <loop over region 2>
>>      DMPlexSetDof() and DMPlexSetFieldDof()
>>    <loop over BC>
>>      DMPlexSetConstraintDof() and DMPlexSetFieldConstraintDof()
>>    DMPlexSetUp()
>>    <loop over BC>
>>      DMPlexSetConstraintIndices() and DMPlexSetFieldConstraintIndices()
>>
>> We could try and package some of this up if it looks generic.
>>
>>    Thanks,
>>
>>       Matt
>>
>>>
>>> Thanks
>>> Reddy
>>>
>>> --
>>> -----------------------------------------------------
>>> Dharmendar Reddy Palle
>>
>>
>>
>>
>> --
>> 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
>
>
>