OK, here is the entire code.
I created a simple mesh to test DMPlex.
As you can see, the GetCone is not working.
As I said, I followed one of the examples for DMPlex using Fortran, I'm not
an exper in F90.
Best regards,
Bernardo M. Rocha
On Fri, Apr 6, 2018 at 11:48 PM, Matthew Knepley <[email protected]> wrote:
> If you send the whole code, we can try it here. I am assuming that the
> examples run for you.
>
> Matt
>
> On Fri, Apr 6, 2018 at 5:47 PM, Bernardo Rocha <
> [email protected]> wrote:
>
>> Hello,
>>
>> I can't get a a simple code to work for DMPlexGetCone.
>> I have tried different ways of handling the array/pointer which is the
>> third argument of DMPlexGetCone,
>> but neither of them worked for me.
>>
>> The last version that I got compiling is this one, which is very much
>> based on the
>> following example: .../petsc-3.7.6/src/dm/impls/p
>> lex/examples/tutorials/ex1f90.F
>>
>> PetscInt, target, dimension(6) :: icone
>> PetscInt, pointer :: pIcone(:)
>>
>> call DMPlexCreateFromCellList(PETSC_COMM_WORLD,
>> & ndim,numel,numnp,nen,PETSC_TRUE,
>> & ielement,ndim,xnodes,dm,ierr)
>> ...
>> call DMPlexGetCone(dm,0,pIcone,ierr)
>> write(*,*) icone(1),icone(2),icone(3),icone(4),icone(5),icone(6)
>> call DMPlexRestoreCone(dm,0,pIcone,ierr)
>>
>> I always get a segmentation fault.
>>
>> I've created a simple cuboid mesh with 2x2x2 hexaedral elements.
>> Many other routines are working and give reasonable results.
>>
>> Output of DMView
>>
>> DM Object: 1 MPI processes
>> type: plex
>> DM_0x84000000_0 in 3 dimensions:
>> 0-cells: 27
>> 1-cells: 66
>> 2-cells: 44
>> 3-cells: 8
>> Labels:
>> depth: 4 strata of sizes (27, 66, 44, 8)
>>
>>
>>
>>
>> On Thu, Apr 5, 2018 at 10:57 PM, Bernardo Rocha <
>> [email protected]> wrote:
>>
>>> OK, thanks a lot.
>>> I had to rename the files and add -ffixed-form as argument to the
>>> compiler.
>>> Best regards,
>>> Bernardo
>>>
>>> On Thu, Apr 5, 2018 at 8:21 PM, Matthew Knepley <[email protected]>
>>> wrote:
>>>
>>>> On Thu, Apr 5, 2018 at 7:16 PM, Bernardo Rocha <
>>>> [email protected]> wrote:
>>>>
>>>>> Hello everyone,
>>>>>
>>>>> I've been trying to use DMPlex in a Fortran 77 application
>>>>> (unfortunately this is legacy code and I can't move to F90).
>>>>>
>>>>> Is there a way to use DMPlexGetCone on it?
>>>>> Although the documentation online says it is only available for F90,
>>>>> I'm wondering if there is a trick like the ones in the users manual
>>>>> for routines
>>>>> that return an array. I tried, but I'm not sure it it is right.
>>>>>
>>>>
>>>> We have standardized our array handling with F90, and I don't think we
>>>> will go back
>>>> to supporting new things for F77-style arrays.
>>>>
>>>> That being said, I think all you have to do is compile that source with
>>>> an F90 compiler. You
>>>> should not have to change anything. If you want the F90 array to look
>>>> like an F77 array,
>>>> just pass it as a function argument (I think this works).
>>>>
>>>> Thanks,
>>>>
>>>> Matt
>>>>
>>>>
>>>>> Best regards,
>>>>> Bernardo M. Rocha
>>>>>
>>>> --
>>>> 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
>>>>
>>>> https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
>>>>
>>>
>>>
>>
>
>
> --
> 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
>
> https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
>
c
program testdmplex
c
implicit real*8 (a-h,o-z)
c
#include <petsc/finclude/petsc.h>
#include <petsc/finclude/petscdmplex.h>
c
DM m
PetscErrorCode ierr
PetscInt icell
PetscInt ipstart, ipend
PetscInt, target, dimension(6) :: icone
PetscInt, pointer :: pIcone(:)
PetscMPIInt size,rank
c
integer ielem(64)
real*8 xx(3),yy(3),zz(3)
real*8 xcoords(27,3)
real*8 xnodes(3*27)
c
c definitions
c
ndim = 3 ! number of dimensions
nen = 8 ! number of element nodes
numel = 8 ! number of elements (2x2x2)
numnp = 27 ! number of nodes (3x3x3)
c
c coordinates
c
xdx = 0.5d0
do i=1,3
xx(i) = (i-1)*xdx
yy(i) = (i-1)*xdx
zz(i) = (i-1)*xdx
end do
c
id = 1
do k=1,3
z = zz(k)
do j=1,3
y = yy(j)
do i=1,3
x = xx(i)
xcoords(id,1) = x
xcoords(id,2) = y
xcoords(id,3) = z
id = id + 1
end do
end do
end do
c
id = 0
do i=1,numnp
id = id + 1
xnodes(id) = xcoords(i,1)
id = id + 1
xnodes(id) = xcoords(i,2)
id = id + 1
xnodes(id) = xcoords(i,3)
end do
c
write(*,*) "coordinates"
do i=1,3*numnp,3
write(*,*) xnodes(i), xnodes(i+1), xnodes(i+2)
end do
c
c element connectivity
c
ielem(1) = 1
ielem(2) = 2
ielem(3) = 5
ielem(4) = 4
ielem(5) = 10
ielem(6) = 11
ielem(7) = 14
ielem(8) = 13
c
ielem(9) = 2
ielem(10) = 3
ielem(11) = 6
ielem(12) = 5
ielem(13) = 11
ielem(14) = 12
ielem(15) = 15
ielem(16) = 14
c
ielem(17) = 4
ielem(18) = 5
ielem(19) = 8
ielem(20) = 7
ielem(21) = 13
ielem(22) = 14
ielem(23) = 17
ielem(24) = 16
c
ielem(25) = 5
ielem(26) = 6
ielem(27) = 9
ielem(28) = 8
ielem(29) = 14
ielem(30) = 15
ielem(31) = 18
ielem(32) = 17
c
ielem(33) = 10
ielem(34) = 11
ielem(35) = 14
ielem(36) = 13
ielem(37) = 19
ielem(38) = 20
ielem(39) = 23
ielem(40) = 22
c
ielem(41) = 11
ielem(42) = 12
ielem(43) = 15
ielem(44) = 14
ielem(45) = 20
ielem(46) = 21
ielem(47) = 24
ielem(48) = 23
c
ielem(49) = 13
ielem(50) = 14
ielem(51) = 17
ielem(52) = 16
ielem(53) = 22
ielem(54) = 23
ielem(55) = 26
ielem(56) = 25
c
ielem(57) = 14
ielem(58) = 15
ielem(59) = 18
ielem(60) = 17
ielem(61) = 23
ielem(62) = 24
ielem(63) = 27
ielem(64) = 26
c
do i=1,64
ielem(i) = ielem(i) - 1
end do
c
write(*,*) "elements"
do i=1,8
k = (i-1)*8
write(*,*) ielem(k+1),ielem(k+2),ielem(k+3),ielem(k+4),
& ielem(k+5),ielem(k+6),ielem(k+7),ielem(k+8)
end do
c
c testing DMPlex in Fortran
c
call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
if (size .ne. 1) then
call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
if (rank .eq. 0) then
write(6,*) 'This is a uniprocessor example only!'
endif
SETERRQ(PETSC_COMM_WORLD,1,' ',ierr)
endif
c
c DMPlex
c
call DMPlexCreateFromCellList(PETSC_COMM_WORLD,
& ndim,numel,numnp,nen,PETSC_TRUE,
& ielem,ndim,xnodes,dm,ierr)
CHKERRQ(ierr)
call DMView(dm,PETSC_VIEWER_STDOUT_WORLD,ierr)
CHKERRQ(ierr)
call DMPlexGetChart(dm,ipstart,ipend,ierr)
CHKERRQ(ierr)
write(*,*) "DMPlex chart start=",ipstart
write(*,*) "DMPlex chart end=",ipend
c
c trying to get the cone of element 0
c
icsz = 0
icell = 0
call DMPlexGetConeSize(dm, icell, icsz, ierr)
CHKERRQ(ierr)
write(*,*) "cone for cell", icell, " has size",icsz
c
c prepare for testing GetCone for cell 0
c
do i=1,6
icone(i) = 0
end do
pIcone => icone
call DMPlexGetCone(dm,0,pIcone,ierr)
CHKERRQ(ierr)
write(*,*)
write(*,*) "results of DMPlexGetCone"
write(*,*) icone(1),icone(2),icone(3),icone(4),icone(5),icone(6)
write(*,*)
call DMPlexRestoreCone(dm,0,pIcone,ierr)
CHKERRQ(ierr)
c
c end
c
call PetscFinalize(ierr)
CHKERRQ(ierr)
end program testdmplex
