Re: [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 <knep...@gmail.com> wrote:
On Thu, Jan 2, 2014 at 4:11 AM, Dharmendar Reddy <dharmaredd...@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