[petsc-users] query DMPlexCreateSection

[Dharmendar Reddy]

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()
>   DMPlexSetUp()
>   <loop over BC>
>     DMPlexSetConstraintIndices() and DMPlexSetFieldConstraintIndices()


Thanks
Reddy

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



-- 
-----------------------------------------------------
Dharmendar Reddy Palle
Graduate Student
Microelectronics Research center,
University of Texas at Austin,
10100 Burnet Road, Bldg. 160
MER 2.608F, TX 78758-4445
e-mail: dharmareddy84 at gmail.com
Phone: +1-512-350-9082
United States of America.
Homepage: https://webspace.utexas.edu/~dpr342