Using OpenGL to plot a mesh

In the past I used Matlab or Octave for visualisation of the meshes [see 3]. However, this often is a waste of time since it takes ages to display the meshes (graphics) - my laptop at home is not up to date. Using OpenGL is amazingly (very) fast. Moreover, I can use it under Windows and Linux [see 4]. The example below shows a simple mesh.

Questions

1.) How can I make the mesh visible, i.e. a black outline for each element; 2.) assure that the aspect ratio is update automatically when the user resize.

      winapp
      PROGRAM Triangle
      INCLUDE <clearwin.ins>,nolist
      INCLUDE <opengl.ins>,nolist

      INTEGER i, ctrl, j
      type fpl_type
        integer  :: NELM,NNODE
        real,dimension(3,10)   :: ITYPE
        integer,dimension(3,10) :: NODE
        real,dimension(50) :: x,y
      end type fpl_type
      type(fpl_type) :: fpl
      double precision :: x,y
      i=winio@('%es%ca[Triagular meshes with OpenGL]&')
      i=winio@('%sp%ww[no_border]%og[static]%lw',
     &         0, 0, 500, 500, ctrl)
!
!        1----------2
!        |\\        /|
!        | \\      / |
!        |  \\    /  |
!        |   \\  /   |
!        |    \\/    |
!        3----4-----5
!        |    /\\    |
!        |   /  \\   |
!        |  /    \\  |
!        | /      \\ |
!        |/        \\|
!        6----------7
!
      fpl%nelm = 6
      fpl%nnode = 7
      fpl%node(1:3,1) = (/1,3,4/)
      fpl%node(1:3,2) = (/1,4,2/)
      fpl%node(1:3,3) = (/2,4,5/)
      fpl%node(1:3,4) = (/6,4,3/)
      fpl%node(1:3,5) = (/6,7,4/)
      fpl%node(1:3,6) = (/7,5,4/)
      fpl%itype(1:3,1) = (/1.0, 0.0, 0.0/)
      fpl%itype(1:3,2) = (/0.0, 1.0, 0.0/)
      fpl%itype(1:3,3) = (/0.0, 0.0, 1.0/)
      fpl%itype(1:3,4) = (/0.0, 0.0, 1.0/)
      fpl%itype(1:3,5) = (/0.0, 1.0, 0.0/)
      fpl%itype(1:3,6) = (/1.0, 0.0, 0.0/)
      fpl%x(1) = 0.0d0; fpl%y(1) = 1.0d0
      fpl%x(2) = 1.0d0; fpl%y(2) = 1.0d0
      fpl%x(3) = 0.0d0; fpl%y(3) = 0.5d0
      fpl%x(4) = 0.5d0; fpl%y(4) = 0.5d0
      fpl%x(5) = 1.0d0; fpl%y(5) = 0.5d0
      fpl%x(6) = 0.0d0; fpl%y(6) = 0.0d0
      fpl%x(7) = 1.0d0; fpl%y(7) = 0.0d0
      CALL glClearColor(0.0, 0.0, 0.0, 0.0)
      CALL glClear(GL_COLOR_BUFFER_BIT)
      CALL glColor3f(1.0, 1.0, 1.0)
      CALL glMatrixMode(GL_PROJECTION)
      CALL glLoadIdentity()
      CALL glOrtho(-0.5d0, 1.5d0, -0.5d0, 1.5d0, -0.5d0, 1.5d0)
      do i=1,fpl%nelm
        CALL glBegin(GL_TRIANGLES)
        CALL glColor3f(fpl%itype(1,i),fpl%itype(2,i),fpl%itype(3,i))
        CALL glVertex2f(fpl%x(fpl%node(1,i)),fpl%y(fpl%node(1,i)))
        CALL glVertex2f(fpl%x(fpl%node(2,i)),fpl%y(fpl%node(2,i)))
        CALL glVertex2f(fpl%x(fpl%node(3,i)),fpl%y(fpl%node(3,i)))
        CALL glEnd()
      enddo
      CALL glFlush()
      END

References

1.) Pie chart 2.) Examples 3.) Mesher 4.) F03GL

1. For visible mesh add one line (i think after CALL glMatrixMode(GL_PROJECTION) )

   call glPolygonMode( GL_FRONT_AND_BACK, GL_LINE )
   

2. to return back fill view from mesh view call

call glPolygonMode( GL_FRONT_AND_BACK, GL_FILL )

instead. I do that using different mouse buttons (i use right mouse for mesh view and left mouse for fill view) or a %rb button with the key (see 'Mesh View' on the panel in MESH.FOR example). There also could be a way to plot both simultaneously, otherwise you will need to plot twice -first the fill, then the line mesh

3. For manipulation of view you better need to read and interpret mouse position and mouse clicks. Add callback function which is doing that like CBmouseGL01 in my mesh example

   i=winio@('%pv%^og[double,depth16]&',640,600,CBmouseGL01) 
   

4. For resizing you will need just to change z if i am not forgotten - in the demo it is done by pushing middle mouse button and simultaneously moving the mouse. see what is done there after this line

if(reasonMouseGL01.EQ.'MOUSE_middle_CLICK'.or.iflag_mouseGL01.eq.16) then

5. to visualize your own mesh you do not need even to change anything in MESH.FOR demo, you just need to read your x,y,z the way code understands it when you drag and drop the input file for it. I removed drag and drop option just for simplicity, but will add it. You can hard wire that with your own programming additions meantime since you've already done principal part of the job, everything else is just the evolutionary beautifications

The way i currently setup input in drag and drop is to define (x,y) for each node and z for its value x1,yA,zD x1,yB,zE x1,yC,zF ....... x2,yK,zP x2,yL,zQ x2,yM,zR ....... x3,yS,zT ..... That's because x in my case is always time slice which is equal for all y,z in a given moment

If you tell me how you do the output for your mesh i will probably make view arbitrary meshes possible when will have some free time in a month or so (reviews, reports, and conferences time! I will do that just for fun because i know the result will look a-m-a-z-i-n-g). For example if you define each (x,y) for the node of the triangle and one value z for all of 3 nodes of a triangle it may look like this x1a,y1a,x2a,y2a,x3a,y3a,z1a x1b,y1b,x2b,y2b,x3b,y3b,z1b .......

or some different way when you define x,y,z for each node

x1a,y1a,z1a,x2a,y2a,z2a,x3a,y3a,z3a ......

Hi Dan

Thanks for the tipps. I found the part in the MESH.FOR where the z-values are plotted (as polygons). This means I have to replace this part with my own data.

The rotate and zoom functions works perfect - just what I was looking for. I also figured out where the beep came from - sound@ 😉

              call glBegin (GL_POLYGON)
              RelColor =(ZmatrGL01(ix,iy+1)-zminGL01)/
     *          (zmaxGL01+1.e-33-zminGL01)
              RelColLum=0.5 
              i = iColorSpectraLines2 ()
              Color_GL01(1)= min(1., ir111/255.)
              Color_GL01(2)= min(1., ig111/255.)
              Color_GL01(3)= min(1., ib111/255.)
              Color_GL01(4)= Transparency3D_GL01
              call glColor4fv(Color_GL01)
              call glVertex3f( XmatrGL01(ix,iy+1), YmatrGL01(iy+1),
     *          ZmatrGL01(ix,  iy+1))
              RelColor =(ZmatrGL01(ix,iy)-zminGL01)/
     *          (zmaxGL01+1.e-33-zminGL01)
              RelColLum=0.5
              i = iColorSpectraLines2 ()
              Color_GL01(1)= min(1., ir111/255.)
              Color_GL01(2)= min(1., ig111/255.)
              Color_GL01(3)= min(1., ib111/255.)
              Color_GL01(4)= Transparency3D_GL01
              call glColor4fv(Color_GL01)        
              call glVertex3f( XmatrGL01(ix,iy),   YmatrGL01(iy),  
     *         ZmatrGL01(ix,  iy))
              RelColor =(ZmatrGL01(ix+1,iy)-zminGL01)/
     *          (zmaxGL01+1.e-33-zminGL01)
              RelColLum=0.5
              i = iColorSpectraLines2 ()
              Color_GL01(1)= min(1., ir111/255.)
              Color_GL01(2)= min(1., ig111/255.)
              Color_GL01(3)= min(1., ib111/255.)
              Color_GL01(4)= Transparency3D_GL01 !  0.5
              call glColor4fv(Color_GL01)
              call glVertex3f( XmatrGL01(ix+1,iy),  YmatrGL01(iy),
     *          ZmatrGL01(ix+1,iy))
              RelColor =(ZmatrGL01(ix+1,iy+1)-zminGL01)/
     *          (zmaxGL01+1.e-33-zminGL01)
              RelColLum=0.5
              i = iColorSpectraLines2 ()
              Color_GL01(1)= min(1., ir111/255.)
              Color_GL01(2)= min(1., ig111/255.)
              Color_GL01(3)= min(1., ib111/255.)
              Color_GL01(4)= Transparency3D_GL01
              call glColor4fv(Color_GL01)
              call glVertex3f( XmatrGL01(ix+1,iy+1), YmatrGL01(iy+1),
     *          ZmatrGL01(ix+1,iy+1))
              call glEnd()

    ! Graphics window and menu: winio@('%mn[Item[subitem]]&', cb1)
    i = winio@('%es%ca[Geompack Meshplot]&')
    i = winio@('%sp%ww[no_border]&',200,200)
    i = winio@('%mn[&File[E&xit]]&', 'exit')
    i = winio@('%mi[PRG_ICON]&')
    i = WINIO@('%ff%bg[white]&')
    i = winio@('%nr%nd&')
    i = winio@('%1.3ob[invisible]&')
    i = winio@(' %ff&') 
    i = winio@(' %ff&') 
    !  Slide bars: win@('%Nsl[option]',value,min,max)
    i = winio@('   Spin Around X-axis%ta%^20sl[horizontal ]%ff&',      &
        RotAngleGL01_X, -180.D0, 180.D0,  cb_dummy)
    i = winio@('   Spin Around Y-axis%ta%^20sl[horizontal ]%ff&',      &
        RotAngleGL01_Y, -180.D0, 180.0D0, cb_dummy)
    i = winio@('   Spin Around Z-axis%ta%^20sl[horizontal ]%ff&',      &
        RotAngleGL01_Z, -180.D0, 180.0D0, cb_dummy)
    i = winio@('   Projection Position%ta%^20sl[horizontal ]%ff&',     &
        AltitudeProjecPlaneGL01_Z, -1.5D0, 1.5D0,cb_dummy)
    i = winio@('%cb&')
    i = winio@(' %ff&')
    i = winio@(' %ff&')
    ! Check boxes
    i = WINIO@('   %^`rb[ No Smooth Method]%ff&',                      &
        k_NoSmooth, '+', cb_dummy, cb_dummy)
    i = WINIO@('   %^`rb[ Quad Strip Smooth Method]%ff&',              &
        kQuadStripSmoothGL01, '+', cb_dummy, cb_dummy)
    i = WINIO@('   %^`rb[ Smooth Colors Inside Polygons]%ff&',         &
        k_smoothColorsInsidePolGL01,'+', cb_dummy, cb_dummy)
    ! Radio buttons
    i = WINIO@('   %^rb[ Mesh view]%ff&',                              &
        kMeshViewGL01, '+', cb_dummy, cb_dummy)
    i = WINIO@('   %^rb[ Print On Screen Window]%ff&',                 &
        kPrintOnScreenGL01, cb_dummy)
    i = WINIO@('%3ga&', kQuadStripSmoothGL01,                          &
        k_smoothColorsInsidePolGL01,k_NoSmooth)
    i = winio@('%cb&')
    i = winio@(' %ff&')
    i = winio@(' %ff&')
    i = winio@('   Frames per Sec  %ta%4rf%ff'//                       &
        '   Min, Max, Ave: %ta%4rf %4rf  %4rf&',FramesPerSecGL01,      &
        FramesPerSecGL01min,FramesPerSecGL01max,FramesPerSecGL01ave)
    ! Control definition
    i = winio@(' %ff %nl&')
    i = winio@('- Right mouse - mesh view of graph%ff&')
    i = winio@('- Ctrl + left mouse - move graph%ff&')
    i = winio@('- Middle mouse down + mouse move - zoom%ff&')
    i = winio@('- Default file for data read - piechart.dat%ff&')
    i = winio@('- Any file can be just drag-n-dropped%ff&')
    i = winio@('%cb&')
    i = winio@('%ob[depressed]&')
    i = winio@('%fn[Courier New]%bf%ts&',1.0d0)
    i = winio@('%pv%^og[double,depth16]&',640,600,cb_dummy)
    i = winio@('%pm[File[Exit]]&', 'exit')
    i = winio@('%cb&')
    i = winio@('%ac[Esc]&', '+','set',lw_ctrlGL01window,0, 'exit')
    i = winio@('%ac[Alt+T]&', 'set',k_StopSpinTestGL01,1)
    i = winio@('%lw', lw_ctrlGL01window)
    i = cb_dummy()

The original example I got from Dan without any functionality is given below. This offers the basic framework for the mesher (I am working on).

Many ClearWin functions needs a call-back function. As I unterstand it, one needs a function as parameter. To do so one could either: 1.) define the function as external (this I assume is Fortran 77 style); or 2.) as an interface (this seems to be 95 onwards).

Question

Is this the correct understanding of external/interface in the context of a function as a parameter?

winapp
program MeshGL01
    implicit none
    include <opengl.ins>,nolist
    include <clearwin.ins>,nolist

    integer :: i
    integer :: k_NoSmooth, kQuadStripSmoothGL01
    integer :: k_smoothColorsInsidePolGL01
    integer :: kMeshViewGL01
    integer :: kPrintOnScreenGL01
    integer :: k_StopSpinTestGL01, lw_ctrlGL01window

    real*8 :: RotAngleGL01_X, RotAngleGL01_Y, RotAngleGL01_Z
    real*8 :: AltitudeProjecPlaneGL01_Z
    real*8 :: FramesPerSecGL01, FramesPerSecGL01ave
    real*8 :: FramesPerSecGL01min, FramesPerSecGL01max

    interface
        integer*4 function cb_dummy()
        end function cb_dummy
    end interface

    kPrintOnScreenGL01          = 0
    k_StopSpinTestGL01          = 0
    k_NoSmooth                  = 0
    kQuadStripSmoothGL01        = 0
    k_smoothColorsInsidePolGL01 = 1
    kMeshViewGL01               = 0
    FramesPerSecGL01min         = 0
    FramesPerSecGL01max         = 0
    FramesPerSecGL01            = 0
    FramesPerSecGL01ave         = 0
!--------------------------
!
! Add previous post here
!
!--------------------------
end program MeshGL01

integer*4 function cb_dummy()
    implicit none
    cb_dummy = 1
    return
end function cb_dummy

RESOURCES
PRG_ICON ICON GEO.ICO

If use Fortran90 style I usually put all callback functions into the module after CONTAINS and USE it, then you do not need to declare them as external

For a test I plotted the same function in Matlab (or Octave). Changing the view requires some time and everthing is much slower than the same example using OpenGL.

Initially I thought the example is nice but nothing for me. In the meanwhile I know that this is exactly what I need. For each triangle I have the (x,y) coordinates and the vector potential, i.e. (x,y,z) for each node and the triangle node information. Between the glBegin() and glEnd() I only have to specify three vertices (i.e. call glVertex3f) instead of four as in the given example. Moreover, I can do the following: 1.) for the mesh only I can use glVertex2f; and 2.) for the vector potential solution glVertex3f.

iydim = 80;
ixdim = 80;

ZoomMesh01init   = 1.5;
ZoomMesh01       = ZoomMesh01init;
StretchXaxisGL01 = 1.2;
StretchYaxisGL01 = 1.2;

XmatrGL01 = zeros(ixdim,iydim);
YmatrGL01 = zeros(iydim);

for iy=1:iydim
	for ix=1:ixdim
		XmatrGL01(ix,iy) = ix/ixdim * ZoomMesh01 * StretchXaxisGL01;
      YmatrGL01(iy)    = iy/iydim * ZoomMesh01 * StretchYaxisGL01; 
    end
end

xCenterMesh = XmatrGL01(ixdim/2,iydim/2);
yCenterMesh = YmatrGL01(iydim/2);
zmaxGL01 = -1.e34;
zminGL01 =  1.e34;

for iy=1:iydim
    for ix=1:ixdim
		x = XmatrGL01 (ix,iy) - xCenterMesh;
		y = YmatrGL01 (iy) - yCenterMesh;
      radxx0 = sqrt( x.^2 + y.^2 );
      ZmatrGL01 (ix,iy)= 0.5*cos(15*radxx0) * exp(-5*radxx0^2)* ZoomMesh01 + 0.6;
    end
end

Dan, I tried you advice and it was easy to plot the mesh! Below is an image of the result.

You mentioned that you used examples you found in the FTN95 Help. Functions in the slab animation are similar. However, I tried this in the past and could not get it to work. I now realized that the example in the online help requires an include (that is not provided). In the FTN95 example directory one will find another slab example.

Question: What options do I have to print/save the image. I would prefer encapsulated postscript, but jpeg will do as well.

[URL=http://imageshack.us/photo/my-images/88/meshgv.jpg/][/URL]

Uploaded with [URL=http://imageshack.us]ImageShack.us[/URL]

Hi Dan, below is the subroutine I use for reading the mesh. The basic format is as follows: 1.) line 1: nde 2.) line 2: nelm, nnode, testa, nlines 3.) nelm lines: (node(i,j),j=1,nde), ITYPE(i), (NEIGH(i,j),j=1,NDE) nodes of triangle, element type, neighboring triangles
4.) nlines lines: (LINE(i,j),j=1,2) begin en end nodes of each line segment 5.) nnode lines: x(i), y(i) node coordinates

Definitions: nde: number of nodes (always thee for triangular mesh) nelm: number of elements nnode: number of nodes testa: testa=0 → no vector potential at node testa=1 → vector potential at node node: triangle nodes itype: type of material (is an integer) neigh: neighboring nodes of a triangle lines: line segments x: x-coordinates y: y-coordinates

SUBROUTINE read_fpl(fls,fpl)
    use FEMPtype
    implicit none
    type(file_type),intent(in) :: fls
    type(fpl_type),intent(out) :: fpl
!
    character(len=256) :: fileToRead
    integer,PARAMETER :: LU = 11

    INTEGER :: I,j,NDE
    
    fileToRead = TRIM(fls%dir_fpl)//fls%SEP//fls%fpl_file    
    OPEN(UNIT=LU,FILE=fileToRead)
    REWIND LU
    read (LU,*,ERR=80) fpl%NDE
    NDE = fpl%NDE
    read (LU,*,ERR=80) fpl%NELM,fpl%NNODE,fpl%TESTA,fpl%NLINES
    DO I=1,fpl%NELM
        READ (LU,*,ERR=80) (fpl%node(i,j),j=1,NDE),fpl%ITYPE(I),(fpl%NEIGH(i,j),j=1,NDE)
    END DO
    DO I=1,fpl%NLINES
        READ (LU,*,ERR=80) (fpl%LINE(i,j),j=1,2)
    END DO
    DO I=1,fpl%NNODE
        READ (LU,*,ERR=80) fpl%x(i),fpl%y(i)
    END DO
    CLOSE(UNIT=LU,STATUS='KEEP')
    return
    80 print *,'read_fpl: Error in while reading the file...'
    STOP
    RETURN
END subroutine

Why do you read NEIGH ? You should find this automatically; For each edge of a triangle; find the other triangle that shares this edge. ( there should be only one) if it exists, then select the other triangle's other node else nominate this edge as an edge to the region. Edges could be defined as ordered nodes, n1:n2 where n1 < n2

An outcome of this analysis would be to identify: all external edges in the model all edges at a boundary between materials ( other triangle has different itype) ( else internal edge, else duplicte connection edge )

This would be a useful subset of edges to draw, confirming some basic geometry. You could draw the 3 (or 4) types of edges with a different colour to validate the model.

John

Do you make it in 3D or just 2D? Do you color it? Post here, it should be looking very nice

Does animation in my mesh example work? With Salford own OpenGL examples i do not need to add anything, i only commented last paragraph in opengl.ins which always caused my OpenGL not to work at all (including i think up to right today)

As to saving image -- i just copy it from the screen. I also have saving function somewhere programmed into the code (John Horspool suggested here in the forum how to do that exactly, as i remember) but could not find it

Hi John and Dan,

thanks for your comments!

Indeed, I actually never (directly) use the neigh variable - it was part of the format I use to come into when I started with the FEM code. I do have somewhere a subroutine outline, which determine the boundaries as you mentioned (I need to search at home for the subroutine - will post it). Since I work with rotational machines this allows me to automatically detect the inner/outer radius as well as the periodic boundaries.

Dan, your example works perfect - including the animation. The actual mesh is 2D - since I only work with 2D field problems (in the screenshot above you can see the function in your original code). The vector potential solution is then visualised in 3D, but here most people prefer the field lines in 2D. I prefer the 3D view since can then visually check if the solution is fine (actually some time ago a found a bug in the solution that was only clear by viewing the 3D plot!).

At present I also do a screenshot. It would however be nice to have an menu with save as...

Hi John,

in a previous thread I posted the code for the very simple mesher I had. Once the triangles are known, the following code calculates the outline and the neigbors. The fpl is only a type declaration for all the mesh/solution variables - this makes handling the data much easier.

subroutine makeneigbors(fpl)
    use FEMPtype
    implicit none
    type(fpl_type),intent(inout) :: fpl

    integer com,m,no1(5),no2(5)
    integer p(5),ct1,ct2,no,i,j,k,l
    integer mar1,mar2,mar3
    integer :: nde,nelm

    fpl%nlines = 0
    fpl%neigh = 0
    nde = fpl%nde
    nelm = fpl%nelm
    do i=1,nelm
        m = 1
        do j=1,nelm
            if ((i == j).and.(i == nelm)) goto 99
            if (j == i) goto 98
            com = 0
            do k=1,nde
                do l=1,nde
                    if (fpl%node(i,k).eq.fpl%node(j,l)) then
                        com = com+1
                        if (com.eq.1) no1(m) = fpl%node(i,k)
                        if (com.eq.1) mar1 = k
                    end if
                    if (com.eq.2) then
                        no2(m) = fpl%node(i,k)
                        mar2   = k
                        if ((mar1+mar2).eq.3) mar3 = 1
                        if ((mar1+mar2).eq.5) mar3 = 2
                        if ((mar1+mar2).eq.4) mar3 = 3
                        fpl%neigh(i,mar3)=j
                        if ((fpl%itype(i).ne.fpl%itype(j)).and.(fpl%neigh(j,1).eq.0)) then
                            fpl%nlines         = fpl%nlines+1
                            fpl%line(fpl%nlines,1) = no1(m)
                            fpl%line(fpl%nlines,2) = no2(m)
                        end if
                        m   = m+1
                        com = 0
                        goto 99
                    end if
                end do
            end do
            99 continue

            if ((m.lt.(nde+1)).and.(m.gt.1).and.(j.eq.nelm)) then
                do while (m.lt.(nde+1))
                    do ct1=1,(nde-1)
                        p(ct1)=fpl%node(i,ct1)*fpl%node(i,ct1+1)
                    end do
                    p(nde)=fpl%node(i,1)*fpl%node(i,nde)
                    do ct1=1,nde
                        no = 0
                        do ct2=1,(m-1)
                            if (p(ct1).eq.(no1(ct2)*no2(ct2))) no = 1
                        end do
                        if (no.eq.0) then
                            no1(m) = fpl%node(i,ct1)
                            no2(m) = p(ct1)/fpl%node(i,ct1)
                        end if
                    end do
                    fpl%nlines     = fpl%nlines+1
                    fpl%line(fpl%nlines,1) = no1(m)
                    fpl%line(fpl%nlines,2) = no2(m)
                    m              = m+1
                end do
            end if
            98 continue
        end do
    end do
    return
end subroutine makeneigbors

here is a link about saving opengl screen.

https://forums.silverfrost.com/Forum/Topic/1822&highlight=opengl

Nice would be to see saving OpenGL screen the way it is done in regular graphics

Don't forget to make small code examples while you learn OpenGL so that more people will come after your steps. Hopefully there will be the whole community of OpenGL programmers from which we will learn quicker then from all manuals. May be some day someone will create Visual OpenGL in Visual Clearwin so that we will stop programming and will just assemble our programs with mouse clicks or finger taps like from lego pieces.

Hi Dan,

good idea - I will do that.

One thing I try to figure is the mapping form the 'real world' to the OpenGL window. The basic idea is to plot data to fit the available screen area. In the example below you can see that only a part of the mesh is shown. Can OpenGL manage this automatically or do I have to take care of that → OpenGL uses a left-handed-system and right-handed-system - the internet is sometimes confusing:(

In a previous thread a few tipps when using only ClearWin was discussed.

[URL=http://imageshack.us/photo/my-images/196/mesh2r.jpg/][/URL]

It is easy to imagine that system can manage that automatically in 2D if you set the rules. In 3D this was realized in physics packages of modern games - they precisely find all limits, view angles and depth of view, boundaries, they ray trace light, do hydrodynamics simulations etc, it's almost real virtual world. We need some higher level language similar to what game developed created for their game models but this one will be for scientific functions and methods. Ones such 'world' will be created users will come adding their contributions and expanding it.

I can get the (book) examples to work. However, when I try my own one it does not work properly - somehow I miss something. Below is an example code to read the coordinates of Lake Superior. The data file is online available: test vector.

winapp
program lsup
  implicit none
  INCLUDE <clearwin.ins>,nolist
  INCLUDE <opengl.ins>,nolist
! Constants
  integer, parameter :: inunit = 11, maxnc = 30, maxnv = 400, maxvc = 500
  integer, parameter :: dp=kind(0.0d0)
! Character
  character ( len = 256 ) filename, rgname
  real ( kind = dp ) :: tolin, angspc, angtol
  real ( kind = dp ) :: kappa, dmin, vcl(2,maxvc)
  real ( kind = dp ) :: xmin,xmax,ymin,ymax
  integer :: case, i, nvc, nmin, ntrid, msglvl, ncur, nv
  integer :: icur(maxnc), ivrt(maxnv), nvbc(maxnc), ctrl
! The data file is available online at: 
! http://people.sc.fsu.edu/~jburkardt/f_src/geompack2/lsup.in
  filename = '..\\test\\geo-03\\lsup'
  open ( unit = inunit, file = TRIM(filename)//'.in')
  read ( inunit, * ) rgname
  read ( inunit, * ) tolin, angspc, angtol, kappa, dmin, nmin, ntrid
  read ( inunit, * ) case, nvc, ncur, msglvl
  read ( inunit, * ) nvbc(1:ncur)
  if ( abs(case) == 2 ) then
    read ( inunit, * ) icur(1:ncur)
  end if
  read ( inunit, * ) ( vcl(1,i), vcl(2,i), i = 1, nvc )
  if ( abs(case) == 2 ) then
    nv = sum ( nvbc(1:ncur) )
    read ( inunit, * ) ivrt(1:nv)
  end if
  close ( unit = inunit )
  xmin =  1d6
  xmax = -1d6
  ymin =  1d6
  ymax = -1d6
  do i=1,nvbc(1)
     xmin = min(xmin,vcl(1,i))
     xmax = max(xmax,vcl(1,i))
     ymin = min(ymin,vcl(2,i))
     ymax = max(ymax,vcl(2,i))
  enddo
  i=winio@('%es%ca[Triagular meshes with OpenGL]&')
  i=winio@('%sp%ww[no_border]%og[static]%lw',0, 0, 800, 600, ctrl)
  CALL glClearColor(0.0, 0.0, 0.0, 0.0)
  CALL glClear(GL_COLOR_BUFFER_BIT)
  CALL glColor3f(1.0, 1.0, 1.0)
  CALL glMatrixMode(GL_PROJECTION)
  call glPolygonMode( GL_FRONT_AND_BACK, GL_LINE )
  CALL glOrtho(xmin, xmax, ymin, ymax, -10.0, 10.0)
  CALL glBegin(GL_POLYGON)
  CALL glColor3f(0.0,1.0,0.0)
  do i=1,nvbc(1)
      CALL glVertex2f(vcl(1,i),vcl(2,i))
  enddo
  CALL glEnd()
  CALL glFlush()
end

From what I figure out one must shift (and scale) the data to fit it to the current window. The following code helped - a screen shot is shown below (on a Linux machine). Anyway, plotting a mesh is such a trivial problem that is not to mention. Many OpenGL examples are looking at 3D examples that are way to complicated for my application (at the moment).

xmin =  5d6
xmax = -5d6
ymin =  5d6
ymax = -5d6
do i=1,fpl%nnodes
    xmin = min(xmin,fpl%x(i))
    xmax = max(xmax,fpl%x(i))
    ymin = min(ymin,fpl%y(i))
    ymax = max(ymax,fpl%y(i))
enddo
w = xmax-xmin
h = ymax-ymin
x_0 = xmin+w/2.0d0
y_0 = ymin+h/2.0d0

[URL=http://imageshack.us/photo/my-images/689/trivials.jpg/][/URL]

Hi Dan

I (re)discovered the viewer of the Triangle mesher, i.e. ShowMe - below is a screenshot. However, it only works under Linux (at home).

Anyway, what I can suggest for the viewer you sent me is the 'control pannel' at the lower left. With OpenGL and such a menu one should be able to write a viewer that can be compiled with FTN95 and whatever compiler one prefer.

Regarding the eps-format: the viewer has a EPS 'button'. This exports the figure in eps-format which I can include this in my documents. Very useful!

At present I use the GEOMPACK mesher. The output I write to the Triangle file format and use the ShowMe viewer.

[URL=http://imageshack.us/photo/my-images/215/showme.png/][/URL]