C INTERNATIONAL AVS CENTER C (This disclaimer must remain at the top of all files) C C WARRANTY DISCLAIMER C C This module and the files associated with it are distributed free of charge. C It is placed in the public domain and permission is granted for anyone to use, C duplicate, modify, and redistribute it unless otherwise noted. Some modules C may be copyrighted. You agree to abide by the conditions also included in C the AVS Licensing Agreement, version 1.0, located in the main module C directory located at the International AVS Center ftp site and to include C the AVS Licensing Agreement when you distribute any files downloaded from C that site. C C The International AVS Center, MCNC, the AVS Consortium and the individual C submitting the module and files associated with said module provide absolutely C NO WARRANTY OF ANY KIND with respect to this software. The entire risk as to C the quality and performance of this software is with the user. IN NO EVENT C WILL The International AVS Center, MCNC, the AVS Consortium and the individual C submitting the module and files associated with said module BE LIABLE TO C ANYONE FOR ANY DAMAGES ARISING FROM THE USE OF THIS SOFTWARE, INCLUDING, C WITHOUT LIMITATION, DAMAGES RESULTING FROM LOST DATA OR LOST PROFITS, OR ANY C SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES. C C This AVS module and associated files are public domain software unless C otherwise noted. Permission is hereby granted to do whatever you like with C it, subject to the conditions that may exist in copyrighted materials. Should C you wish to make a contribution toward the improvement, modification, or C general performance of this module, please send us your comments: why you C liked or disliked it, how you use it, and most important, how it helps your C work. We will receive your comments at avs@ncsc.org. C C Please send AVS module bug reports to avs@ncsc.org. C C NOTE: THIS MODULE AND SOURCE CODE IS FOR USE C WITH THE AVS SOFTWARE ENVIRONMENT ONLY C----------------------------------------------------------------------- C C (C) Copyright Evans & Sutherland Computer Limited 1992 C C EVANS & SUTHERLAND ASSUMES NO RESPONSIBILITY FOR ANY USE OR INABILITY C TO USE THIS SOFTWARE. THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY C OF ANY KIND WHATSOEVER. EVANS & SUTHERLAND DISCLAIMS ALL IMPLIED WARRANTIES C INCLUDING THOSE OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. C C PERMISSION IS GRANTED TO USE, COPY, MODIFY AND DISTRIBUTE THIS SOFTWARE C PROVIDED THAT: ALL SUBSEQUENT DISTRIBUTION IS FREE OF CHARGE; THE C COPYRIGHT NOTICE AND THIS PERMISSION NOTICE ARE PRESERVED IN ALL COPIES. C C----------------------------------------------------------------------- C C ISOLINE SLICE module C C Author: Mike French, Evans & Sutherland UK C C Date: 18 March 1992 C C----------------------------------------------------------------------- * SUBROUTINE AVSinit_modules * *-- include files INCLUDE '/usr/include/avs/avs.inc' * *-- external declarations EXTERNAL LINE_COMPUTE * *-- ports and parameters INTEGER OPORT1, IPORT1, IPORT2, IPORT3, IPARAM * *-- multi-module process flag INTEGER MULTI * *-- flags MULTI = IOR( COOPERATIVE, REENTRANT ) * *-- Set the module name * CALL AVSset_module_name( 'Isoline Slice', 'mapper' ) * *-- Create input port * *-- field data IPORT1 = AVScreate_input_port( 'input field', + 'field 2D scalar real 3-space', REQUIRED ) * *-- colormap IPORT2 = AVScreate_input_port( 'input colormap', + 'colormap', OPTIONAL ) * *-- contour level input IPORT3 = AVScreate_input_port( 'input levels', + 'field 1D scalar real uniform', OPTIONAL ) * *-- Create output port * *-- geometry OPORT1 = AVScreate_output_port( 'output geom', 'geom') * *-- Parameters * *-- line offset choice IPARAM = AVSadd_parameter( 'zflag', 'choice', 'no offset', + 'no offset!+!-!+ & -', '!' ) *-- offset factor IPARAM = AVSadd_parameter( 'offset', 'float', + 0.0035, 0.0001, 0.0050 ) *-- number of contours IPARAM = AVSadd_parameter( 'N', 'integer', 10, 2, 20 ) * *-- Automatically free space for the output data after each call * CALL AVSautofree_output( OPORT1 ) * *-- Set the compute process * CALL AVSset_compute_proc( LINE_COMPUTE ) * *-- Request single argument passing of fields and allow multi option * CALL AVSset_module_flags( IOR(SINGLE_ARG_DATA,MULTI) ) * * RETURN END * ************************************************************************ *-- *-- ===== Isoline Slice - compute function ===== *-- ************************************************************************ * INTEGER FUNCTION LINE_COMPUTE( + FPTR, CPTR, LPTR, OUTPUT, + ZFLAG, OFFSET, NCPAR ) * *-- include files INCLUDE '/usr/include/avs/avs.inc' INCLUDE '/usr/include/avs/geom.inc' * *-- declare arguments INTEGER FPTR, CPTR, LPTR, OUTPUT, NCPAR, SIZE REAL LO, HI, H(256), S(256), V(256), A(256) REAL OFFSET CHARACTER ZFLAG*(*) * *-- maximum grid dimensions INTEGER MXGRID PARAMETER (MXGRID=200) * *-- flag the interpolartion method *-- this could be made an AVS parameter, but vertex is always better CHARACTER INTERPOLATION_METHOD*6 PARAMETER (INTERPOLATION_METHOD='vertex') * *-- input connections INTEGER NOT_THERE, IS_THERE PARAMETER (NOT_THERE=0, IS_THERE=1) * *-- limit for real comparison REAL EPSILON PARAMETER (EPSILON=1E-5) * *-- axis indices INTEGER XAXIS, YAXIS, ZAXIS PARAMETER (XAXIS=1, YAXIS=2, ZAXIS=3) * *-- parameter widget display options INTEGER HIDE, SHOW PARAMETER (HIDE=0, SHOW=1) * *-- type of coordinate data for interpolation, switch for 3D LOGICAL DATA2D, DATA3D PARAMETER (DATA2D=.FALSE., DATA3D=.TRUE.) * *-- declare locals INTEGER OBJ, PERCNT, RESULT, AVS_ALL, NPT, N INTEGER FIND, FOFF, FTYPE, NSPACE, NDIM,DIM(2),NX,NY,NXY INTEGER IOFF(4), JOFF(4), I, J, K, II, JJ, IC, INDEX INTEGER PUNI, PIRR, POFF, IX, IY, IZ, LIND, LOFF INTEGER LDIM(1), NC, N_DIAL, O_DIAL REAL VERT(3,2), FIELD(1), L(1) REAL LEVEL, CLEVEL(20), R,G,B, COLOR(3,2), P(1) REAL X(4), Y(4),Z(4), F(4), XPT(2,2),YPT(2,2),ZPT(2,2) REAL XVAL, YVAL, ZVAL, XMAX, YMAX REAL OFFWID, OFFVAL, NO_Z, NO_ZPT REAL XMIN, YMIN, GRIDDX, GRIDDY REAL MINEXT(3), MAXEXT(3), DMIN, DMAX, E(3) REAL FNX, FNY, FNZ, VNX, VNY, VNZ, XOFF, YOFF, ZOFF REAL IRROFF(MXGRID,MXGRID,3) CHARACTER OBUFF*80, EBUFF*80 LOGICAL GRID_IS_UNIFORM, GRID_IS_IRREGULAR, FIRST * *-- Local data * *-- flag first pass DATA FIRST / .TRUE. / SAVE FIRST SAVE N_DIAL, O_DIAL * *-- cell traversal DATA IOFF / 0, 1, 1, 0 / DATA JOFF / 0, 0, 1, 1 / * *-- grid position function XVAL(I) = XMIN + FLOAT(I-1) * GRIDDX YVAL(J) = YMIN + FLOAT(J-1) * GRIDDY * *-- field data index function FIND(I,J) = FOFF + I + (J-1)*NX * *-- field coordinate index functions * *-- uniform field PUNI(I) = POFF + I * *-- irregular fields PIRR(I,J,N) = POFF + I + (J-1)*NX + (N-1)*NXY * *-- contour level index functions LIND(I) = LOFF + I * *-- constant definition AVS_ALL = IOR( IOR(AVS_MINVAL,AVS_MAXVAL), AVS_VALUE) * * *-- Initialize parameter dial visibility * IF (FIRST) THEN O_DIAL = HIDE RESULT = AVSparameter_visible( 'offset', O_DIAL ) N_DIAL = SHOW RESULT = AVSparameter_visible( 'N', N_DIAL ) FIRST = .FALSE. END IF * *-- Initialize input colourmap * IF (CPTR.NE.0) THEN RESULT = AVScolormap_get( CPTR,256,SIZE,LO,HI,H,S,V,A ) IF (RESULT.NE.1) THEN CALL AVSerror( 'Cannot get input colormap data' ) LINE_COMPUTE = 0 RETURN END IF END IF * *-- Create geometry information * *-- create frame object OBJ = GEOM_create_obj( GEOM_POLYTRI, GEOM_NULL ) * *-- Get field type * FTYPE = AVSfield_get_int( FPTR, AVS_FIELD_UNIFORM ) * IF (FTYPE.EQ.UNIFORM) THEN GRID_IS_UNIFORM = .TRUE. GRID_IS_IRREGULAR = .FALSE. ELSE IF (FTYPE.EQ.IRREGULAR) THEN GRID_IS_UNIFORM = .FALSE. GRID_IS_IRREGULAR = .TRUE. ELSE CALL AVSerror( 'Cannot contour rectilinear field (yet)' ) LINE_COMPUTE = 0 RETURN END IF * *-- Get dimensions, data pointer and coordinate space * NDIM = AVSfield_get_int( FPTR, AVS_FIELD_NDIM ) IF (NDIM.NE.2) THEN CALL AVSerror( 'Incorrect field dimension, must be 2D input' ) LINE_COMPUTE = 0 RETURN END IF * RESULT = AVSfield_get_dimensions( FPTR, DIM ) * IF ( (DIM(1).GE.MXGRID) .OR. + (DIM(2).GE.MXGRID) ) THEN CALL AVSerror( 'Field slice too large, increase downsize' ) LINE_COMPUTE = 0 RETURN END IF * NX = DIM(1) NY = DIM(2) NXY = NX*NY * NSPACE = AVSfield_get_int( FPTR, AVS_FIELD_NSPACE ) IF (NSPACE.NE.3) THEN CALL AVSerror( 'Incorrect field space dimension, must be 3D' ) LINE_COMPUTE = 0 RETURN END IF * RESULT = AVSfield_data_offset( FPTR, FIELD, FOFF ) * *-- Get data range and field extent * *-- get extent RESULT = AVSfield_get_extent( FPTR, MINEXT, MAXEXT ) E(1) = MAXEXT(1) - MINEXT(1) E(2) = MAXEXT(2) - MINEXT(2) E(3) = MAXEXT(3) - MINEXT(3) * *-- reset field minimum and maximum CALL AVSfield_reset_minmax( FPTR ) RESULT = AVSfield_get_minmax( FPTR, DMIN, DMAX ) * *-- find positions from points array *-- same as extents except for dimension with plane position RESULT = AVSfield_points_offset( FPTR, P, POFF ) * *-- Generate contour levels * IF (LPTR.NE.0) THEN * *-- contour level input exists * *-- remove 'N contour' parameter widget IF (N_DIAL.EQ.SHOW) THEN N_DIAL = HIDE RESULT = AVSparameter_visible( 'N', N_DIAL ) END IF * *-- get contour input data offset and dimension RESULT = AVSfield_data_offset( LPTR, L, LOFF ) RESULT = AVSfield_get_dimensions( LPTR, LDIM ) IF (RESULT.EQ.0) THEN CALL AVSerror( 'Failed to get contour level input dimensions') LINE_COMPUTE = 0 RETURN END IF NC = LDIM(1) * *-- copy contour input to the local array DO 210 IC=1,NC CLEVEL(IC) = L(LIND(IC)) 210 CONTINUE * ELSE * *-- contour level input does not exist *-- *-- replace 'N contour' parameter widget IF (N_DIAL.EQ.HIDE) THEN N_DIAL = SHOW RESULT = AVSparameter_visible( 'N', N_DIAL ) END IF * *-- levels to be generated locally NC = NCPAR DO 310 IC=1,NC CLEVEL(IC) = DMIN + (DMAX-DMIN)*FLOAT(IC-1)/FLOAT(NC-1) 310 CONTINUE * END IF * *-- Generate the grid parameters, widths and offsets * IF (GRID_IS_UNIFORM) THEN * *-- find plane of grid IF (P(PUNI(1)).EQ.P(PUNI(2))) THEN *-- x-plane IZ = XAXIS IY = YAXIS IX = ZAXIS ELSE IF (P(PUNI(3)).EQ.P(PUNI(4))) THEN *-- y-plane IX = XAXIS IZ = YAXIS IY = ZAXIS ELSE IF (P(PUNI(5)).EQ.P(PUNI(6))) THEN *-- z-plane IX = XAXIS IY = YAXIS IZ = ZAXIS ELSE *-- none of the planes CALL AVSerror( 'Cannot find plane of slice for uniform field') LINE_COMPUTE = 0 RETURN END IF * *-- allocate positions XMIN = P(PUNI(2*IX-1)) XMAX = P(PUNI(2*IX)) YMIN = P(PUNI(2*IY-1)) YMAX = P(PUNI(2*IY)) ZVAL = P(PUNI(2*IZ-1)) * *-- set uniform grid cell widths GRIDDX = (XMAX-XMIN) / FLOAT(NX-1) GRIDDY = (YMAX-YMIN) / FLOAT(NY-1) * END IF * *-- use total diagonal extent to scale offset OFFWID = SQRT( E(1)*E(1) + E(2)*E(2) + E(3)*E(3) ) * *-- Generate the offset value and manage the offset value widget * IF (ZFLAG.EQ.'no offset') THEN * *-- set zero offset OFFVAL = 0.0 * *-- remove offset dial widget IF (O_DIAL.EQ.SHOW) THEN O_DIAL = HIDE RESULT = AVSparameter_visible( 'offset', O_DIAL ) END IF * ELSE IF (ZFLAG.EQ.'-') THEN * *-- set negative offset value OFFVAL = -OFFSET * OFFWID * *-- replace offset dial widget IF (O_DIAL.EQ.HIDE) THEN O_DIAL = SHOW RESULT = AVSparameter_visible( 'offset', O_DIAL ) END IF * ELSE * *-- set positive offset value OFFVAL = OFFSET * OFFWID * *-- replace offset dial widget IF (O_DIAL.EQ.HIDE) THEN O_DIAL = SHOW RESULT = AVSparameter_visible( 'offset', O_DIAL ) END IF END IF * *-- Apply offsets * IF (GRID_IS_UNIFORM) THEN * *-- simple global offset for uniform grids ZVAL = ZVAL + OFFVAL * ELSE IF (GRID_IS_IRREGULAR) THEN * IF (INTERPOLATION_METHOD.EQ.'facet ') THEN * *-- initialise offsets XOFF = 0.0 YOFF = 0.0 ZOFF = 0.0 * ELSE IF ((INTERPOLATION_METHOD.EQ.'vertex') .AND. + (ZFLAG.NE.'no offset')) THEN * *-- calculate offsets based on vertex normals for irregular grids DO 500 I=1,NX DO 510 J=1,NY CALL VERTEX_NORMAL( NX, NY, P(POFF+1), I,J, VNX,VNY,VNZ ) IRROFF(I,J,1) = OFFVAL * VNX IRROFF(I,J,2) = OFFVAL * VNY IRROFF(I,J,3) = OFFVAL * VNZ 510 CONTINUE 500 CONTINUE * *-- end of branch on interpolation method END IF * *-- end of branch on grid type END IF * *-- Loop over contours * DO 600 IC=1,NC * *-- notify status PERCNT = 10 + INT( 80.0 * FLOAT(IC-1) / FLOAT(NC-1) ) CALL AVSmodule_status( "contouring", PERCNT ) * *-- set contour level LEVEL = CLEVEL(IC) * *-- set contour colour IF (CPTR.EQ.0) THEN * *-- make black contours R = 0.0 G = 0.0 B = 0.0 * ELSE * *-- set contour colour from colourmap IF (LEVEL.LE.LO) THEN INDEX = 1 ELSE IF (LEVEL.GE.HI) THEN INDEX = SIZE ELSE INDEX = 1 + INT( FLOAT(SIZE-1) * (LEVEL-LO) / (HI-LO) ) END IF * *-- convert colour from hsv to rgb colour model CALL HSVA_TO_RGB( H(INDEX), S(INDEX), V(INDEX), R, G, B ) * END IF * *-- load constant colour for vector COLOR(1,1) = R COLOR(2,1) = G COLOR(3,1) = B COLOR(1,2) = R COLOR(2,2) = G COLOR(3,2) = B * *-- Loop over grid cells * DO 700 J=1,NY-1 DO 710 I=1,NX-1 * *-- create contour lines for each type of cell IF (GRID_IS_UNIFORM) THEN * *-- generate uniform coordinates DO 810 K=1,4 II = I + IOFF(K) JJ = J + JOFF(K) F(K) = FIELD(FIND(II,JJ)) X(K) = XVAL(II) Y(K) = YVAL(JJ) 810 CONTINUE * *-- look-up the 2D line pattern for the uniform cell CALL CELL_LINES( LEVEL, DATA2D, X, Y, NO_Z, F, + NPT, XPT, YPT, NO_ZPT ) * *-- loop over number of line segments, either one or two DO 820 N=1,NPT * *-- map coordinates from plane to 3-space VERT(IX,1) = XPT(N,1) VERT(IY,1) = YPT(N,1) VERT(IZ,1) = ZVAL VERT(IX,2) = XPT(N,2) VERT(IY,2) = YPT(N,2) VERT(IZ,2) = ZVAL * *-- draw line segment(s) CALL GEOM_add_polyline(OBJ,VERT,COLOR,2,GEOM_COPY_DATA) IF (ZFLAG.EQ.'+ & -') THEN VERT(IZ,1) = ZVAL - 2.0*OFFVAL VERT(IZ,2) = ZVAL - 2.0*OFFVAL CALL GEOM_add_polyline( OBJ, VERT, COLOR, 2, + GEOM_COPY_DATA ) END IF * 820 CONTINUE * ELSE IF (GRID_IS_IRREGULAR) THEN * *-- facet normal displacement IF (INTERPOLATION_METHOD.EQ.'facet ') THEN * *-- load coordinates and field values DO 811 K=1,4 II = I + IOFF(K) JJ = J + JOFF(K) F(K) = FIELD(FIND(II,JJ)) X(K) = P(PIRR(II,JJ,1)) Y(K) = P(PIRR(II,JJ,2)) Z(K) = P(PIRR(II,JJ,3)) 811 CONTINUE * *-- look-up the 3D line pattern for the irregular cell CALL CELL_LINES( LEVEL,DATA3D,X,Y,Z,F, NPT,XPT,YPT,ZPT ) * *-- displace point along facet normal IF (ZFLAG.NE.'no offset') THEN CALL FACET_NORMAL( X, Y, Z, FNX, FNY, FNZ ) XOFF = OFFVAL * FNX YOFF = OFFVAL * FNY ZOFF = OFFVAL * FNZ END IF * *-- loop over number of line segments DO 821 N=1,NPT * *-- offset line vertices VERT(1,1) = XPT(N,1) + XOFF VERT(2,1) = YPT(N,1) + YOFF VERT(3,1) = ZPT(N,1) + ZOFF VERT(1,2) = XPT(N,2) + XOFF VERT(2,2) = YPT(N,2) + YOFF VERT(3,2) = ZPT(N,2) + ZOFF * *-- draw line segment CALL GEOM_add_polyline( OBJ, VERT, COLOR, + 2, GEOM_COPY_DATA ) * *-- draw lines on reverse side of slice IF (ZFLAG.EQ.'+ & -') THEN VERT(1,1) = XPT(N,1) - 2.0*XOFF VERT(2,1) = YPT(N,1) - 2.0*YOFF VERT(3,1) = ZPT(N,1) - 2.0*ZOFF VERT(1,2) = XPT(N,2) - 2.0*XOFF VERT(2,2) = YPT(N,2) - 2.0*YOFF VERT(3,2) = ZPT(N,2) - 2.0*ZOFF CALL GEOM_add_polyline( OBJ, VERT, COLOR, 2, + GEOM_COPY_DATA ) END IF * 821 CONTINUE * ELSE IF (INTERPOLATION_METHOD.EQ.'vertex') THEN * *-- add offsets before interpolation DO 812 K=1,4 II = I + IOFF(K) JJ = J + JOFF(K) F(K) = FIELD(FIND(II,JJ)) IF (ZFLAG.EQ.'no offset') THEN X(K) = P(PIRR(II,JJ,1)) Y(K) = P(PIRR(II,JJ,2)) Z(K) = P(PIRR(II,JJ,3)) ELSE X(K) = P(PIRR(II,JJ,1)) + IRROFF(II,JJ,1) Y(K) = P(PIRR(II,JJ,2)) + IRROFF(II,JJ,2) Z(K) = P(PIRR(II,JJ,3)) + IRROFF(II,JJ,3) END IF 812 CONTINUE * *-- look-up the 3D line pattern for the irregular cell CALL CELL_LINES( LEVEL,DATA3D,X,Y,Z,F, NPT,XPT,YPT,ZPT ) * *-- loop over number of line segments DO 822 N=1,NPT * *-- copy line vertices VERT(1,1) = XPT(N,1) VERT(2,1) = YPT(N,1) VERT(3,1) = ZPT(N,1) VERT(1,2) = XPT(N,2) VERT(2,2) = YPT(N,2) VERT(3,2) = ZPT(N,2) * *-- draw line segment CALL GEOM_add_polyline( OBJ, VERT, COLOR, + 2, GEOM_COPY_DATA ) * 822 CONTINUE * *-- draw on other side of slice if required IF (ZFLAG.EQ.'+ & -') THEN * *-- reverse offsets and redo interpolation DO 813 K=1,4 II = I + IOFF(K) JJ = J + JOFF(K) X(K) = X(K) - 2.0*IRROFF(II,JJ,1) Y(K) = Y(K) - 2.0*IRROFF(II,JJ,2) Z(K) = Z(K) - 2.0*IRROFF(II,JJ,3) 813 CONTINUE * *-- look-up the 3D line pattern for the irregular cell CALL CELL_LINES( LEVEL, DATA3D, X, Y, Z, F, + NPT, XPT, YPT, ZPT ) * *-- loop over number of line segments DO 823 N=1,NPT * *-- copy line vertices VERT(1,1) = XPT(N,1) VERT(2,1) = YPT(N,1) VERT(3,1) = ZPT(N,1) VERT(1,2) = XPT(N,2) VERT(2,2) = YPT(N,2) VERT(3,2) = ZPT(N,2) * *-- draw line segment CALL GEOM_add_polyline( OBJ, VERT, COLOR, + 2, GEOM_COPY_DATA ) * 823 CONTINUE * *-- end of branch on z flag END IF * *-- end of branch on interpolation method END IF * *-- end of branch on grid type END IF * *-- end of loop over grid cellls 710 CONTINUE 700 CONTINUE * *-- end of loop over contours 600 CONTINUE * *-- Copy objects to AVS * *-- initialize list of changes to null OUTPUT = GEOM_init_edit_list( OUTPUT ) * *-- add object geometry to output CALL GEOM_edit_geometry( OUTPUT, 'contour', OBJ ) * *-- fix contour lines relative to top level CALL GEOM_edit_transform_mode( OUTPUT,'contour','parent',0 ) * *-- push banded surfaces if necessary IF (ZFLAG.EQ.'no offset') THEN *-- use the new AVS4 geometry viewer option RESULT = AVScommand('kernel', + 'geom_set_render_mode -object top push_surface', + OBUFF, EBUFF ) END IF * *-- free the memory used while constructing the list CALL GEOM_destroy_obj( OBJ ) * *-- Successful completion * LINE_COMPUTE = 1 * * 999 RETURN END