! Reference: H.C. Huang, A.S. Usmani, Finite Element Analysis for Heat ! Transfer, London, Springer-Verlag, 1994. ! ********************************************************************** ! ** Program mesh2d automatically generates 2-d meshes given a set ** ! ** of outer and inner boundary segments enclosing the domain to be ** ! ** meshed, and a reference element size. It generates node points ** ! ** on the boundary segments, and within the domain according to a ** ! ** given mesh density. These points are connected to produce the ** ! ** mesh using Watson's algorithm for Delaunay triangulation. ** ! ** A. S. Usmani sep. 1987 ** ! ********************************************************************** ! Translated from FORTRAN 66 to Fortran 90, corrected, generalized. ! Copyright 1999, J. E. Akin, all rights reserved. MODULE PRECISION_MODULE IMPLICIT NONE ! Use hardware defaults for single and double precision INTEGER, PARAMETER :: SP = KIND (1.0) ! Single precision INTEGER, PARAMETER :: DP = KIND (1.d0) ! Double precision ! Note: Use SELECTED_INT_KIND or SELECTED_REAL_KIND for user ! defined portable precision. However, user choice may not be ! on all hardware so check it with functionality below. CONTAINS SUBROUTINE CHECK_HARDWARE_PRECISION ! Verify that the hardware supports the request IF ( KIND (SP) < 0 ) THEN PRINT *, 'ERROR: SP PRECISION NOT AVAILABLE' STOP 'SET SP = KIND (1.0) IN MODULE PRECISION_MODULE' END IF IF ( KIND (DP) < 0 ) THEN PRINT *, 'ERROR: DP PRECISION NOT AVAILABLE' STOP 'SET DP = KIND (1.d0) IN MODULE PRECISION_MODULE' END IF END SUBROUTINE CHECK_HARDWARE_PRECISION END MODULE PRECISION_MODULE MODULE keyword_buffer ! For nonadvancing keyword, data input ! Copyright J. E. Akin, 1998 ! Input lines (of up to LINE_SIZE characters) have a keyword (of up ! to MAX_KEY characters) and may be followed by a mixture of ! integers, reals, or strings separated by a gap. ! Strings with more than one word should be in quotes ("), or use ! the underscore (_) instead of gaps. ! A gap is a blank space, comma, or horizontal tab. implicit none public :: get_int, get_key, get_real, get_string ! subroutines public :: MAX_KEY, LINE_SIZE, WORD_SIZE ! size parameters public :: ECHO_UNIT, ECHO_KEY, ECHO_OPEN ! for debug private ! everything not listed as public integer, parameter :: LINE_SIZE = 80, MAX_KEY = 16 integer, parameter :: WORD_SIZE = LINE_SIZE - MAX_KEY integer, parameter :: ECHO_UNIT = 37 character(len=5), parameter :: FMT_LINE = "(a80)" character(len=1), parameter :: BLA = ' ', COM = ',' character(len=1), parameter :: TAB = achar(9) character(len=LINE_SIZE) :: buffer ! input buffer line logical :: ECHO_KEY = .FALSE. ! true for debug logical :: ECHO_OPEN = .FALSE. ! echo file status integer :: LINES_READ = 0 ! line information logical :: LINE_FINISHED = .FALSE. ! line information integer :: ioResult ! 0=no error, <0=EOF integer :: NEXT_START, START, STOP, LENGTH, N_CHARS contains ! encapsulated functionality subroutine open_echo_unit !--------------------------------------------------------- ! Open the file to echo keywords that are parsed !--------------------------------------------------------- implicit none integer :: ok ! file status if ( ECHO_KEY ) then open (ECHO_UNIT, file = 'keyword_input.echo', & action ='write', status = 'unknown', & iostat = ok) if ( ok == 0 ) ECHO_OPEN = .TRUE. else print *,'Note: set keyword_buffer ECHO_KEY = .true.' ECHO_OPEN = .FALSE. end if end subroutine open_echo_unit function is_a_gap (new) result (gap) !--------------------------------------------------------- ! Test if a character is a blank, comma, or tab !--------------------------------------------------------- implicit none character(len=1), intent(in) :: new logical :: gap gap = .false. ! initialize if (new == BLA .or. new == TAB .or. new == COM) gap = .true. end function is_a_gap subroutine locate_next_word !--------------------------------------------------------- ! Locate start and stop characters of the next word !--------------------------------------------------------- implicit none integer :: j ! loops do j = NEXT_START, LENGTH ! find non gap if ( .not. is_a_gap (buffer(j:j)) ) then ! word starts START = j ; exit ! this loop end if end do ! for START STOP = START ! initialize NEXT_START = 1 ! initialize LINE_FINISHED = .TRUE. ! initialize do j = START, LENGTH ! find next gap if ( is_a_gap (buffer(j:j)) ) then ! found gap if ( j /= LENGTH ) LINE_FINISHED = .FALSE. NEXT_START = j ! passed word exit ! this loop else ; STOP = j ! update stop end if end do ! for STOP N_CHARS = STOP - START + 1 ! number of characters end subroutine locate_next_word subroutine get_key (key_word) !--------------------------------------------------------- ! get the first keyword on line !--------------------------------------------------------- implicit none character(len=MAX_KEY), intent(out) :: key_word character(len=13) :: S_FORMAT ! for debug ! read a line of input into the buffer key_word = ' ' read (5, FMT_LINE, iostat = ioResult) buffer ! input buffer line if ( ioResult < 0 ) then ! test for EOF key_word = 'EOF' print *, 'get_key: EOF found after line ', LINES_READ ! debug else LINES_READ = LINES_READ + 1 ! increment count LENGTH = MAX (1, len_trim (buffer)) ! true length START = 1 ; STOP = LENGTH ; NEXT_START = 1 ! initialize word call locate_next_word ! find keyword key_word (1:N_CHARS) = buffer(START:STOP) ! copy keyword if ( ECHO_KEY ) then ! begin debug echo ! create format via internal file write (S_FORMAT, '("(A", I2, ",1X)")' ) N_CHARS write (ECHO_UNIT, S_FORMAT, advance='no') & key_word (1:N_CHARS) ! append key_word end if ! echo end if ! EOF end subroutine get_key subroutine get_string (string) !--------------------------------------------------------- ! get the next string on the line !--------------------------------------------------------- implicit none character(len=WORD_SIZE), intent(out) :: string character(len=13) :: S_FORMAT ! for debug integer :: j ! loops string = ' ' ; N_CHARS = 1 ! initialize if ( .not. LINE_FINISHED ) then ! a string should be present call locate_next_word ! find the string size if ( N_CHARS > 1 .and. buffer(START:START) == '"' ) then do j = START + 1, LENGTH ! locate end of quotes STOP = j if ( buffer(j:j) == '"' ) exit ! this loop end do ! for next quote N_CHARS = STOP - START + 1 ! revise number of characters end if ! string starts with a quote string (1:N_CHARS) = buffer(START:STOP) ! copy from buffer if ( ECHO_KEY ) then ! begin debug echo ! create format via internal file write (S_FORMAT, '("(A", I2, ",1X)")' ) N_CHARS write (ECHO_UNIT, S_FORMAT, advance='no') & string (1:N_CHARS) ! append string end if ! echo end if ! LINE_FINISHED end subroutine get_string subroutine get_int (fixed) !--------------------------------------------------------- ! get the next integer !--------------------------------------------------------- implicit none integer, intent(out) :: fixed ! the integer if ( LINE_FINISHED ) then ! assign default of zero fixed = 0 ; N_CHARS = 1 else ! should be present call locate_next_word ! locate integer read (buffer(START:STOP), *) fixed ! read integer from buffer if ( ECHO_KEY ) write (ECHO_UNIT, '(I9)', advance='no') fixed end if ! LINE_FINISHED end subroutine get_int subroutine get_real (float) ! get the next real number !--------------------------------------------------------- ! get the next real number !--------------------------------------------------------- Use precision_module implicit none real(dp), intent(out) :: float ! the real number if ( LINE_FINISHED ) then ! assign default of zero float = 0.d0 ; N_CHARS = 1 else ! should be present call locate_next_word ! locate real number read (buffer(START:STOP), *) float ! read float from buffer if ( ECHO_KEY ) write (ECHO_UNIT, '(1PE14.6, 1X)', & advance='no') float ! append float end if end subroutine get_real END MODULE keyword_buffer ! copyright 1999, J. E. Akin, all rights reserved. Module H_Mesh_Constants Use Precision_Module IMPLICIT NONE Integer :: MELEM, MPOIN, MNPEL, MMATR, MBOUN, MCDPT, MDIME, & MPROF, MDENS, MBPLN, MPNOD, & mhol, mreg, matt, melm, melm2, mnpl, mpoi, mpoi2, mtx, mx, & mxfn, mtri, mtfn, mxin, mele4, mpoi4, mfac Integer :: NELEM, NPOIN, NNPEL, NMATR, NBOUN, NCDPT, NDIME, & NGAUS, NPROF, NDENS, NFRON, NBPLN, NPNOD, & nhol, nreg, natt Integer :: M_SEG, M_USER, N_WARN !b ! dummies for debuging key words CHARACTER (len=72) :: title Integer :: T_SETS, T_STEPS Real(dp) :: TS_A, TS_B Logical :: AXISYMMETRIC CONTAINS SUBROUTINE SET_H_MESH_DEFAULTS ! WARNING: Multiple limits are set but NEVER checked in the original ! code. Count checks are being added. MBPLN = 99 ! max number of boundary segments for a geometry MBOUN = 999 ! max number of input boundaries MCDPT = 4 ! max number of points defining segment properties MDIME = 2 ! max number of spatial dimensions MELEM = 9999 ! max number of elements to generate MNPEL = 9 ! max number of nodes per element MPNOD = 299 ! max number of nodes or el's on a boundary segment MPOIN = 9999 ! max number of nodes to generate MPROF = 16999 ! max number of profile (skyline) coefficients mfac = 29999 ! max number of element faces in mesh mhol = 9 ! max number of hole regions in mesh mreg = 59 ! max number of total regions in mesh M_SEG = 100 ! max number of segments (el) on a boundary M_USER = 69 ! max number of points on user defined boundary matt = mreg - mhol ! max number of material regions MMATR = matt ! max number of material regions MDENS = MPOIN ! max number of density (size) points mele4 = MELEM ! max number of Q4 elements allowed mpoi4 = MPOIN ! max number of nodes in a Q4 mesh ! Alternate names used for above data mtx = mpoin/2 ; mx = mpoin mxfn = mpoin ; mtri = melem mtfn = melem ; mxin = mpoin*4 melm = melem ; melm2 = 10*melem ; mnpl = mnpel mpoi = mpoin ; mpoi2 = mpoin*20 ! in promin END SUBROUTINE SET_H_MESH_DEFAULTS End Module H_Mesh_Constants program mesher !b USE PRECISION_MODULE IMPLICIT NONE ! Reference parameter settings are as follows: ! mtx = mpoin/2, mx = mpoin, mxfn = mpoin, mtri = melem, ! mtfn = melem, mxin = mpoin*4 integer, parameter :: mtx = 4999, mx = 9999, mxfn = 9999, & mtri = 9999, mtfn = mtri, mhol = 9, mreg = 59, matt = 50, & mxin = mx*4, mtem = mtri, mpoin = mx, mdime = 2, mden = mx, & melem = mtfn, mnpel = 6 integer, parameter :: mbpln = 199, mpnod = 299 !b !b mden = mxfn ?? = mx !b mpoin = mxfn ?? = mx INTEGER, PARAMETER :: MMATR = 2, MBOUN = 999, MCDPT = 4, & MPROF = 16999, MDENS = mx !b MDENS = MPOIN = mden = melem ?? ! For no transient, or non-linear INTEGER, PARAMETER :: ITRAN = 0, ILINR = 0, IPETR = 0, & NITUP = 0, NITDN = 0, IAXSY = 0, IERST = 0, & ICONV = 1 REAL (DP), PARAMETER :: TTIME = 0.d0, STIME =0.d0, DTIME = 0.d0, & DTMAX = 0.d0, ALPHA = 0.d0, RELAX =0.d0, & TOLER =0.d0, FACTR = 0.d0 ! For control INTEGER :: IPASS, IOPTN, NSTEP, ITRAD, NITRA, ICOAR, NITER REAL (DP) :: PCENT, ELMIN, ELMAX, TIME, RCENT ! For MESH2D INTEGER :: IADAP, IFRNT, NPOIN, NELEM, NNPEL, NMATR, NDENS !b , & !b NBPLN, NBPOI, NBREF, NBELN, NBELF, NBELE, & !b NBELT, NBRET !b INTEGER :: MTYPE (MELEM), NCONC (MELEM, MNPEL) INTEGER :: itype INTEGER :: ntot, ntr, nnpe, nmat, ndpt, nbpln INTEGER :: ntri !b delete INTEGER :: npinb (189), nbno (189, 489), npin6 (189), & nbno6 (189, 979), nbpoi (mbpln), nbref (mbpln, mpnod), & nbeln (mbpln), nbret (mbpln, mpnod), nbelt (mbpln, mpnod), & NEWEL (melem), nbelf (mbpln, mpnod), nbele (mbpln, mpnod), & mtype (melem), nconc (melem, mnpel), NEWNO (mpoin) REAL (DP) :: coord (mpoin, mdime), xd (mden), yd (mden), vald (mden) ! For INDATA INTEGER :: NCDPT, NOBCD, NOBCN INTEGER :: NFIXD (MBOUN), NFACB (MBOUN), NELMB (MBOUN) REAL(DP) :: CONDY (MMATR), CAPCY (MMATR), FIXED (MBOUN), & FLUXE (MBOUN), COEFF (MBOUN), RADIA (MBOUN), & AMBIT (MBOUN), TEMPR (MPOIN), CDVLU (MCDPT, MMATR), & CPVLU (MCDPT, MMATR), TVALU (MCDPT, MMATR) INTEGER :: NBCTM (MBPLN), NBCFL (MBPLN) REAL(DP) :: TEMBC (MCDPT), FLUXP (MCDPT), COEFP (MCDPT), & RADIP (MCDPT), AMBIP (MCDPT), TEMIN (MMATR), & UVELM (MMATR), VVELM (MMATR) !For TOMESH INTEGER :: inumb (MPOIN), itemp (mpoin, mmatr),& NFIXB, NEUMN, I_O INTEGER :: I_BC (MPOIN) !b REAL(DP) :: UVELO (MPOIN), VVELO (MPOIN) CHARACTER (len=72) :: title !b OPEN (UNIT = 7, STATUS = 'UNKNOWN', FILE = 'in_put.out') OPEN (UNIT = 10, STATUS = 'OLD', FILE = 'h_mesh.dat') OPEN (unit = 11, status = 'old', file = 'h_geom.dat') OPEN (unit = 12, status = 'old', file = 'densit.dat') OPEN (UNIT = 13, STATUS = 'UNKNOWN', FILE = 'history.log') OPEN (unit = 16, status = 'unknown', file = 'bound.res') OPEN (unit = 17, status = 'unknown', file = 'mesh.res') OPEN (unit = 21, status = 'unknown', file = 'msh_bc_xyz.dat') OPEN (unit = 22, status = 'unknown', file = 'msh_typ_nodes.dat') OPEN (unit = 23, status = 'unknown', file = 'msh_bc_values.dat') OPEN (unit = 24, status = 'unknown', file = 'msh_type_flux.dat') OPEN (unit = 31, status = 'unknown', file = 'msh_bc_xyz.tmp') OPEN (unit = 32, status = 'unknown', file = 'msh_typ_nodes.tmp') OPEN (unit = 33, status = 'unknown', file = 'msh_bc_values.tmp') OPEN (unit = 34, status = 'unknown', file = 'msh_type_flux.tmp') WRITE (6, 403) ; 403 FORMAT ( & /10x,'please input element type.', & /10x,'1 -- triangle with three nodes;', & /10x,'2 -- triangle with six nodes;', & /10x,'3 -- quadrilateral with four nodes;') !b , & !b /10x,'4 -- quadrilateral with eight nodes;', & !b /10x,'5 -- quadrilateral with nine nodes;') READ (5, * ) itype IF (itype == 1) nnpe = 3 IF (itype == 2) nnpe = 6 IF (itype == 3) nnpe = 4 !b IF (itype == 4) nnpe = 8 !b IF (itype == 5) nnpe = 9 NNPEL = nnpe !b notation duplicated ! *** READ CONTROL DATA ! IPASS = 0 ; IOPTN = 1 ; NSTEP = 0 ; IFRNT = 0 ; ITRAD = 0 NITRA = 10 ; ICOAR = 1 ; PCENT = 1.d0 PCENT = 1.0D0 ; NITER = 1 ! ! *** READ AND WRITE TITLE ! READ (10, 920, IOSTAT = I_O) TITLE ; WRITE (7, 920) TITLE IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR TITLE' 920 FORMAT(12A6) WRITE (7, *) TITLE ! ! *** READ AND WRITE CONTROL PARAMETERS TIME = STIME ; RCENT = PCENT ifrnt = 0 ; nmat = 1 ; NMATR = 1 !b WRITE (7, 901) ITRAN, ILINR, IAXSY, IERST, ICONV, IPETR 901 FORMAT (/, 'ITRAN =',I4, ' ILINR =',I4, ' IAXSY =',I4, & & ' IERST =',I4, ' ICONV =',I4, ' IPETR =',I4) ! ! *** READ IN GEOMETRY DATA AND CREATE MESH ! iadap = 0 ! basic case, 1 = adaptive, 2 = special CALL mesh2d (mpoin, melem, mnpel, mbpln, mpnod, mdime, & mden, iadap, ifrnt, ntot, ntr, nnpe, nmat, ndpt, nbpln, & nbpoi, nbref, nbeln, nbelf, nbele, nconc, mtype, nbelt, & nbret, newel, newno, coord, xd, yd, vald, itype, npoin, ntri ) !b !b nbret, newel, newno, coord, xd, yd, vald) nelem = ntr !b ! *** READ PROBLEM DATA ! CALL INDATA (MMATR, MBPLN, MCDPT, NMATR, NCDPT, NOBCD, NOBCN, & NBPLN, ITRAN, ILINR, ICONV, NBCTM, NBCFL, TEMBC, FLUXP, COEFP, & RADIP, AMBIP, TEMIN, CONDY, CAPCY, CDVLU, CPVLU, TVALU, UVELM, & VVELM) ! ! *** TRANSFER PROBLEM DATA TO MESH ! 888 CALL TOMESH (MPOIN, MELEM, MNPEL, MMATR, MBPLN, MCDPT, MPNOD, & MDIME, MBOUN, ITRAN, ILINR, ICONV, NPOIN, NELEM, NNPEL, NMATR, & NOBCD, NOBCN, NFIXB, NEUMN, NBPLN, NCONC, MTYPE, NBPOI, NBREF, & NBELN, NBELF, NBELE, NBCTM, NBCFL, NFIXD, NELMB, NFACB, FIXED, & TEMBC, FLUXE, COEFF, RADIA, AMBIT, FLUXP, COEFP, RADIP, AMBIP, & TEMPR, TEMIN, UVELO, VVELO, UVELM, VVELM, COORD, ITEMP, INUMB, & NSTEP, I_BC) ! OUTPUT TO MODEL PROGRAM CALL output_mesh (mdime, melem, mnpel, mpoin, ntot, ntr, nnpe, & coord, mtype, nconc, I_BC, NOBCN) CALL output_boundary (mbpln, mpnod, nbpln, nbele, nbelf, nbeln, & nbpoi, nbref) end program mesher !b subroutine h_mesh_key (key) !--------------------------------------------------------- ! free format input of application data !--------------------------------------------------------- use H_Mesh_Constants ! for data definable here use keyword_buffer ! for MAX_KEY, WORD_SIZE implicit none character(len=MAX_KEY), intent (in) :: key character(len=WORD_SIZE) :: on_off character(len=WORD_SIZE) :: word integer, parameter :: limit = 10 ! max allowed bad words integer, save :: count = 0 ! number of bad words word = 'null' on_off = 'off' select case ( key ) ! *** add any application specific control keyword case here *** ! (they must be previously defined in module System_Constants) ! Case (' ! *** end of application dependent case actions *** ! Check allowed dummy keys of lines to skip Case (' ' ) ! ignore blank line Case ('#' ) ! ignore comment line Case ('!' ) ! ignore comment line Case ('?' ) ! ignore comment line ! problem title Case ('title' ) ; call get_string (title) ! program logical controls Case ('axisymmetric') ; AXISYMMETRIC = .TRUE. ! program control integers Case ('data_set') ; call get_int (T_SETS) ! test multiple inputs per line (two integers) !b Case ('test_2_i') ; call get_int (T_SETS) !b print *,'T_SETS ', T_SETS !b call get_int (T_STEPS) !b print *,'T_STEPS ', T_STEPS !b ! program control reals Case ('TS_B') ; call get_real (TS_B) ! time step ratio TS_A = 1.d0 - TS_B ! time step ratio ! test multiple inputs per line (two reals) !b Case ('test_2_r') ; call get_real (TS_A) !b print *, 'TS_A ', TS_A !b call get_real (TS_B) !b print *, 'TS_B ', TS_B !b ! test multiple inputs per line (integer then a real) !b Case ('test_i_r') ; call get_int (T_SETS) !b print *, 'T_SETS ', T_SETS !b call get_real (TS_B) !b print *, 'TS_B ', TS_B !b ! all other keywords Case default ! all other words Select Case (key(1:1)) ! on first character only Case (' ' ) ! ignore blank line Case ('#' ) ! ignore comment line Case ('!' ) ! ignore comment line Case ('?' ) ! ignore comment line Case default ! all other words print *, 'WARNING, h_mesh_key: unknown keyword ', key n_warn = n_warn + 1 count = count + 1 if ( count >= limit ) then print *, 'STOP: reached limit on unknown words' print *, 'Expecting to find "end" or "quit"' stop 'No end or quit keyword found in control' end if ! likely user error end select ! on first character only end select ! from key end subroutine h_mesh_key SUBROUTINE mesh2d (mpoin, melem, mnpel, mbpln, mpnod, mdime, & mden, iadap, ifrnt, ntot, ntr, nnpe, nmat, ndpt, nbpln, & nbpoi, nbref, nbeln, nbelf, nbele, nconc, mtype, nbelt, & nbret, newel, newno, coord, xd, yd, vald, itype, npoin, ntri) !b !b nbret, newel, newno, coord, xd, yd, vald) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) INTEGER, INTENT (IN) :: mpoin, melem, mnpel, mbpln, mpnod, & mden, iadap, ifrnt, itype INTEGER, INTENT (OUT) :: ntot, ntr, nnpe, nmat, ndpt, nbpln, & nbpoi, nbref, nbeln, nbelf, nbele, nconc, mtype, nbelt, & nbret, newel, newno, ntri ! ! ********************************************************************** ! ** ** ! ** This program automatically generates 2-d meshes given a set of ** ! ** outer and inner boundary segments enclosing the domain to be ** ! ** meshed, and a reference element size. It generates node points ** ! ** on the boundary segments, and within the domain according to a ** ! ** given mesh density. These points are connected to produce the ** ! ** mesh using Watson's algorithm for delaunay triangulation. ** ! ** ** ! ** A. S. Usmani sep. 1987 ** ! ** ** ! ********************************************************************** ! ! Reference parameter settings are as follows: ! ! mtx = mpoin/2 ! mx = mpoin ! mxfn = mpoin ! mtri = melem ! mtfn = melem ! mxin = mpoin*4 ! ! For stand-alone program ! parameter (mtx=4999,mx=9999,mxfn=9999,mtri=9999,mtfn=mtri, ! . mhol=9,mreg=59,matt=50,mxin=mx*4,mtem=mtri, ! . mpoin=mxfn,mdime=2, ! . mden=mxfn,melem=mtfn,mnpel=6,mbpln=99,mpnod=299) ! For subroutine PARAMETER (mtx = 4999, mx = 9999, mxfn = 9999, mtri = 9999, & mtfn = mtri, mhol = 9, mreg = 59, matt = 50, & mxin = mx * 4, mtem = mtri) DIMENSION x (mx), y (mx), xc (mtri), yc (mtri), r2 (mtri), & xt (mtx), yt (mtx), szt (mtx), xnuwh (mxin), ynuwh (mxin), & rnuwh (mxin), xnu (mxin), ynu (mxin), rnu (mxin), & xnukp (mx), ynukp (mx), xx (3), xy (3), & xorig (mhol), yorig (mhol), xd (mden), yd (mden), & vald (mden), xm (mtem), ym (mtem), rm (mtem), xtemp (mtx), & ytemp (mtx), stemp (mtx), xb (mtx), yb (mtx), szb (mtx), & xfin (mxfn), yfin (mxfn), xquad (mxin), yquad (mxin), & rquad (mxin), coord (mpoin, mdime) INTEGER :: nl (mtri), nm (mtri), nn (mtri), iflwh (mxin), & iflag (mxin), icrsp (mxin), lnods (mtfn, 6), & !b lnod4 (mtfn, 4), l (199), m (199), n (199), nref (199), & !b l1 (199), l2 (199), ninb (mhol), nodi (mtx, mhol), & iatt (matt), lma4 (mtfn), ibtri (mtri), nnb (mreg), & nrfnnb (mtx, mreg), nchk (mhol), nodref (mx), & nodchk (mx), nlm (mtfn), nmn (mtfn), nnl (mtfn), & mbtri (mtfn), mbqua (mtfn), iflat (mtri), nbtri (mtx), & nodtri (mtx, 20), ncref (mx), ml (mtfn), mm (mtfn), & mn (mtfn), nat (mtfn), iflin (mxin), iornt (mxin) DIMENSION npinb (189), nbno (189, 489), npin6 (189), & !b nbno6 (189, 979), nbpoi (mbpln), nbref (mbpln, mpnod), & !b nbeln (mbpln), nbret (mbpln, mpnod), nbelt (mbpln, mpnod), & NEWEL (melem), nbelf (mbpln, mpnod), nbele (mbpln, mpnod), & mtype (melem), nconc (melem, mnpel), NEWNO (mpoin) !b print *, 'entered mesh2d ' !b iprog = 1 !b 0 pi = 3.1415926536d0 !b search elsewhere beta0 = pi / 2.d0 ipert = 1 nitsm = 3 ! number of smoothing cycles CALL input (nt, xt, yt, szt, s, nbnds, natrib, iatt, incr, ipol, & nnb, nrfnnb, tolrn, ndpt, xd, yd, vald, nnpe, ord, nbpl, nbnp, & npinb, nbno, mtx, matt, mreg, mden, iadap, ifrnt, nitsm, iprog) ntri = 0 ndum = 0 ! ! *** for triangle quadtree node insertion, set `ityqd'=2 ! ! if (nnpe == 3) ityqd=2 ! if (nnpe == 6) ityqd=2 IF (nnpe == 3) ityqd = 1 IF (nnpe == 6) ityqd = 1 IF (nnpe == 4) ityqd = 1 IF (nnpe == 8) ityqd = 1 IF (nnpe == 9) ityqd = 1 !b print *, 'iprog, itype, nnpe, ityqd ', iprog, itype, nnpe, ityqd ! ! *** generate points for all domains in an imaginary square ! *** based on the given density distribution ! CALL sort (nt, xt, yt, xmax, xmin, ymax, ymin, mtx) write (13,*)'xmin=',xmin write (13,*)'xmax =',xmax CALL refsiz (xmax, xmin, ymax, ymin, xd, yd, vald, ndpt, & s, idpth, mden) write (13,*)'idpth=',idpth WRITE (13, * ) 'start nodins' CALL nodins (nt, xt, yt, tolrn, xd, yd, vald, ndpt, xnuwh, & ynuwh, rnuwh, npwh, s, ord, xmax, xmin, ymax, & ymin, mtx, mden, ityqd, iadap, idpth, xquad, & yquad, rquad, iflin, iornt, mxin) !b DO i = 1, npwh !b iflwh (i) = 0 !b END DO iflwh (1:npwh) = 0 !b ! ! *** ! WRITE (13, * ) 'end nodins' WRITE (13, * ) 'nt=', nt WRITE (13, * ) 'npwh=', npwh WRITE (13, * ) 'xnuwh(1)=', xnuwh (1) WRITE (13, * ) 'ynuwh(1)=', ynuwh (1) ! ! *** ! !b DO i = 1, nt !b nodchk (i) = 0 !b END DO nodchk (1:nt) = 0 !b IF (nbnds == 0) GO TO 50 DO i = 1, nbnds nchk (i) = 0 k = nrfnnb (1, i) xorig (i) = xt (k) yorig (i) = yt (k) END DO 50 nbeg = nbnds + 1 nend = natrib + nbnds ntot = nt ntr = 0 DO im = nbeg, nend npoin = 0 iatm = im - nbnds iat = iatt (iatm) noutb = nnb (im) npoin = npoin + noutb DO i = 1, noutb nod = nrfnnb (i, im) xb (i) = xt (nod) yb (i) = yt (nod) x (i) = xb (i) y (i) = yb (i) szb (i) = szt (nod) nodref (i) = nod END DO ! CALL sort (noutb, xb, yb, xmax, xmin, ymax, ymin, mtx) newp = 0 DO i = 1, npwh IF (iflwh (i) == 0) then newp = newp + 1 !b xxx check against max xnu (newp) = xnuwh (i) ynu (newp) = ynuwh (i) rnu (newp) = rnuwh (i) icrsp (newp) = i END IF END DO ! ! *** ! WRITE (13, * ) 'start to delete external points', nbnds ! inbnd = 0 ninnb = 0 IF (nbnds == 0) GO TO 70 DO i = 1, nbnds IF (nchk (i) == 1) GO TO 60 x0 = xorig (i) y0 = yorig (i) CALL inorout (x0, y0, noutb, xb, yb, ichk, mtx, iangl) IF (ichk == 0) GO TO 60 nchk (i) = 1 inbnd = inbnd+1 ninb (inbnd) = nnb (i) ntemp = ninb (inbnd) DO j = 1, ntemp k = npoin + j itemp = nrfnnb (j, i) nodi (j, inbnd) = k xtemp (j) = xt (itemp) ytemp (j) = yt (itemp) x (k) = xtemp (j) y (k) = ytemp (j) stemp (j) = szt (itemp) nodref (k) = itemp END DO ninnb = ninnb + ntemp npoin = npoin + ntemp CALL reject (i, ntemp, iflag, xtemp, ytemp, stemp, iflwh, & newp, xnu, ynu, rnu, icrsp, xm, ym, rm, xmax, & xmin, ymax, ymin, 1, mtx, mxin, mtem) 60 CONTINUE END DO 70 CALL reject (im, noutb, iflag, xb, yb, szb, iflwh, newp, & xnu, ynu, rnu, icrsp, xm, ym, rm, xmax, xmin, & ymax, ymin, 2, mtx, mxin, mtem) DO i = 1, newp iflwh (icrsp (i) ) = im END DO WRITE (13, * ) 'end deleting of points' ! ! *** introduce perturbation of inserted points ! *** to avoid the degeneration of triangulation ! WRITE (13, * ) 'npoin=', npoin DO i = 1, newp xnukp (i) = xnu (i) ynukp (i) = ynu (i) END DO nptmp = npoin iiptb = 0 newpt = newp IF (newpt == 0) newpt = 1 alpha = pi / 2.d0 * dsqrt (5.0d0) / 2.d0 IF (ipert /= 0) then DO i = 1, newp xnu (i) = xnukp (i) & + rnu (i) / 100000.d0 * dcos (float (i) * alpha) ynu (i) = ynukp (i) & + rnu (i) / 100000.d0 * dsin (float (i) * alpha) END DO END IF incr = 1 ipol = 2 ! GO TO 76 75 CONTINUE iiptb = iiptb + 1 beta = float (iiptb) * pi / 2.d0 * dsqrt (5.0d0) + pi / 4.d0 i = ipoit - nptmp xnu (i) = xnukp (i) + rnu (i) / 100000.d0 * dcos (beta) ynu (i) = ynukp (i) + rnu (i) / 100000.d0 * dsin (beta) 76 CONTINUE ! DO i = 1, newp k = nptmp + i x (k) = xnu (i) y (k) = ynu (i) END DO npoin = nptmp + newp gama = 1.1d0 const = 2.d0 * s ilarg = 0 ! rotat = - 1 * dsqrt (2.0d0) * pi / 4.d0 xcbox = (xmin + xmax) / 2 ycbox = (ymin + ymax) / 2 xsbox = xmax - xmin ysbox = ymax - ymin rbox = dsqrt (2.0d0) / 2 * dmax1 (xsbox, ysbox) ! 77 CONTINUE ilarg = ilarg + 1 flarg = float (ilarg) rotat = rotat + dsqrt (2.0d0) * pi / 4 const = gama * const ! write (13,*)'const=',const ! rbox=rbox+const rbox = rbox * (0.5d0 + flarg / (flarg + 1) ) rotri = 2 * rbox roqua = dsqrt (2.0d0) * rbox ! ! *** define the first fictitious triangles encompassing the set of data ! *** points given. ! IF (ipol == 1) then ntri = 1 nl (1) = npoin + 1 nm (1) = npoin + 2 nn (1) = npoin + 3 i1 = npoin + 1 i2 = npoin + 2 i3 = npoin + 3 !b xxx check against max ! x(i1)=xmin-const ! y(i1)=ymin-const ! x(i2)=xmax+ymax+2.*const ! y(i2)=y(i1) ! x(i3)=x(i1) ! y(i3)=x(i2) x (i1) = xcbox + rotri * dcos (rotat) y (i1) = ycbox + rotri * dsin (rotat) x (i2) = xcbox + 1.1 * rotri * dcos (2 * pi / 3 + rotat) y (i2) = ycbox + 1.1 * rotri * dsin (2 * pi / 3 + rotat) x (i3) = xcbox + rotri * dcos (4 * pi / 3 + rotat) y (i3) = ycbox + rotri * dsin (4 * pi / 3 + rotat) ! write (13,*)'initial polygon (triangle)' ! write (13,*)' x(i1)=',x(i1),' x(i2)=',x(i2),' x(i3)=',x(i3) ! write (13,*)' y(i1)=',y(i1),' y(i2)=',y(i2),' y(3)=',y(i3) ELSEIF (ipol == 2) then ntri = 2 nl (1) = npoin + 1 nm (1) = npoin + 2 nn (1) = npoin + 3 nl (2) = nl (1) nm (2) = nn (1) nn (2) = npoin + 4 !b xxx check against max i1 = npoin + 1 i2 = npoin + 2 i3 = npoin + 3 i4 = npoin + 4 ! xhf=(xmax-xmin)/2 ! yhf=(ymax-ymin)/2 ! xav=(xmax+xmin)/2 ! yav=(ymax+ymin)/2 ! x(i1)=xmin-const-yhf ! y(i1)=yav ! x(i2)=xav ! y(i2)=ymin-const-xhf ! x(i3)=xmax+const+yhf ! y(i3)=y(i1) ! x(i4)=x(i2) ! y(i4)=ymax+const+xhf x (i1) = xcbox + roqua * dcos (rotat) y (i1) = ycbox + roqua * dsin (rotat) x (i2) = xcbox + 1.1 * roqua * dcos (pi / 2 + rotat) y (i2) = ycbox + 1.1 * roqua * dsin (pi / 2 + rotat) x (i3) = xcbox + roqua * dcos (pi + rotat) y (i3) = ycbox + roqua * dsin (pi + rotat) x (i4) = xcbox + 1.2 * roqua * dcos (3 * pi / 2 + rotat) y (i4) = ycbox + 1.2 * roqua * dsin (3 * pi / 2 + rotat) ! write (13,*)' xcbox=',xcbox,' ycbox=',ycbox ! write (13,*)' rotri=',rotri,' roqua=',roqua ! write (13,*)' rotat=',rotat ! write (13,*)'initial polygon (quad)' ! write (13,*)' x(i1)=',x(i1),' x(i2)=',x(i2) ! write (13,*)' x(i3)=',x(i3),' x(i4)=',x(i4) ! write (13,*)' y(i1)=',y(i1),' y(i2)=',y(i2) ! write (13,*)' y(i3)=',y(i3),' y(i4)=',y(i4) END IF DO i = 1, ntri xx (1) = x (nl (i) ) xy (1) = y (nl (i) ) xx (2) = x (nm (i) ) xy (2) = y (nm (i) ) xx (3) = x (nn (i) ) xy (3) = y (nn (i) ) CALL circum (xx, xy, cx, cy, rsq, icirc) xc (i) = cx yc (i) = cy r2 (i) = rsq END DO ! ! *** loop over the no. of node points (npoin), one by one, to obtain ! *** the final triangulaion. ! add = 1.d0 - 0.5d0 * (1.d0 / incr) tnpoin = npoin tempo = tnpoin / incr + add mt = tempo DO ii = 1, incr i = 0 DO jj = 1, mt IF (jj == 1) i = ii IF (jj > 1) i = i + incr IF (i > npoin) GO TO 20 nrej = 0 ! ! *** loop over the no. of triangles so far to identify the ones whose ! *** circumcircle contains 'i', i.e. the current node point. ! DO j = 1, ntri xi = x (i) yi = y (i) xcj = xc (j) ycj = yc (j) r2j = r2 (j) CALL check (xi, yi, xcj, ycj, r2j, ichk) IF (ichk == 1) then nrej = nrej + 1 !b xxx check against max l (nrej) = nl (j) m (nrej) = nm (j) n (nrej) = nn (j) nref (nrej) = j END IF END DO CALL inserp (nrej, l, m, n, l1, l2, npols) nrej2 = nrej + 2 IF (npols /= nrej2) then ipoit = i ! if (iiptb < 4) then WRITE (13, 10) npols, nrej 10 FORMAT (//10x,'triangulation failed as the no. of ', & & 'new triangles',2x,i3/10x,'is not equal to the ', & & 'no. of rejected triangles'/10x,i3,2x,'plus 2.'//) IF (ipoit <= nptmp) then WRITE (13, * ) 'inserp failure for point', i WRITE (13, * ) 'the current point is on boundary' WRITE (13, * ) 'recovering' incr = incr + 2 IF (ipol == 1) ipol = 2 IF (ipol == 2) ipol = 1 GO TO 77 ELSE WRITE (13, * ) 'inserp failure for point', i WRITE (13, * ) 'the current point is inside subdomain' WRITE (13, * ) 'perturbing point' GO TO 75 END IF ! else ! write(13,15)npols,nrej ! 15 format(10x,'program stopped as the no. of new ', ! 1 'triangles',2x,i3//10x,'is not equal to the no. ', ! 2 'of rejected triangles'//10x,i3,2x,'plus 2,', ! 3 'and 4 times of perturbation have been done.') ! stop ! end if END IF !xxxx CALL update (x, y, nl, nm, nn, ntri, xc, yc, r2, npols, & l1, l2, nref, i, mx, mtri, icirc) IF (icirc == 1) then ipoit = i IF (ipoit <= nptmp) then WRITE (13, * ) 'update failure for point', i WRITE (13, * ) 'the current point is on boundary' WRITE (13, * ) 'recovering' incr = incr + 2 IF (ipol == 1) ipol = 2 IF (ipol == 2) ipol = 1 GO TO 77 ELSE WRITE (13, * ) 'update failure for point', i WRITE (13, * ) 'the current point is inside subdomain' WRITE (13, * ) 'perturbing point' GO TO 75 END IF END IF 20 CONTINUE END DO END DO ! ! *** reject from the final triangulation the external triangles ! CALL final (npoin, ntri, nl, nm, nn, noutb, inbnd, ninb, & nodi, iflag, mxin, mx, mtx, mtri, mhol, xb, & yb, x, y, xtemp, ytemp) CALL bchek (mtx, mx, mtri, mhol, mxin, npoin, ntri, nl, nm, & nn, noutb, inbnd, ninb, nodi, iflag, iflat, & ibtri, nbtri, nodtri, xb, yb, x, y) IF (iprog == 1) then WRITE (6, 501) im - nbnds 501 FORMAT(//10x,'for subdomain ',i2, & & /10x,'enter 0 for no smoothing' & & /10x,'enter x for x smoothing iterations'//' !') !b READ (5, * ) nitsm END IF nitsm = 3 !b IF (nitsm /= 0) then DO jism = 1, nitsm CALL smooth (npoin, x, y, ntri, nl, nm, nn, noutb, & ninnb, mx, mtri) END DO END IF CALL collect (npoin, x, y, ntri, nl, nm, nn, iat, ntot, & xfin, yfin, ntr, ml, mm, mn, nat, noutb, & ninnb, nodref, nodchk, ncref, ibtri, mbtri, & mx, mtri, mxfn, mtfn) !b print *,'at 692 mpoin, npoin ', mpoin, npoin !b !b print *,'at 692 melem, ntr, ntri ', melem, ntr, ntri !b if ( npoin > mpoin ) stop '694 npoin > mpoin' !b if ( ntri > melem ) stop '694 ntri > melem' !b if ( ntr > melem ) stop '694 ntr > melem' !b END DO ! ! *** interpolate the final triangulation for midside nodes and output ! CALL conver (ntot, xfin, yfin, ntr, ml, mm, mn, nlm, nmn, nnl, & nat, nnpe, natrib, mxfn, mtfn, mbtri, mbqua, iatt, & matt, nmat, lnods, lnod4, lma4, nbpl, npinb, nbno, & npin6, nbno6) ! ! *** copy the nodal connectivity and coordinates into HEAT2D arrays ! DO i = 1, ntr IF (nnpe == 3) then nconc (i, 1) = ml (i) nconc (i, 2) = mm (i) nconc (i, 3) = mn (i) mtype (i) = nat (i) ELSEIF (nnpe == 6) then nconc (i, 1) = ml (i) nconc (i, 2) = nlm (i) nconc (i, 3) = mm (i) nconc (i, 4) = nmn (i) nconc (i, 5) = mn (i) nconc (i, 6) = nnl (i) mtype (i) = nat (i) ELSEIF (nnpe == 4) then nconc (i, 1) = lnod4 (i, 1) nconc (i, 2) = lnod4 (i, 2) nconc (i, 3) = lnod4 (i, 3) nconc (i, 4) = lnod4 (i, 4) mtype (i) = lma4 (i) !b need 8 & 9 here too xxx END IF !b print *,'729 i, nconc ', i, nconc (i, 1:nnpe) !b END DO ! ! *** find the nodes and elements on the boundary ! nbpln = 0 DO 301 ibpl = 1, nbpl ! if (npinb(ibpl) <= 0) GO TO 301 nbpln = nbpln + 1 !b xxx check mbpln here IF (nnpe == 3) then nb = iabs (npinb (ibpl) ) nbpoi (nbpln) = nb DO i = 1, nb nbret (nbpln, i) = nbno (ibpl, i) END DO ELSE nb = iabs (npin6 (ibpl) ) nbpoi (nbpln) = nb DO i = 1, nb nod = nbno6 (ibpl, i) nbret (nbpln, i) = nod END DO END IF nbeln (nbpln) = 0 IF (nnpe /= 4) then DO ii = 1, npinb (ibpl) - 1 jj = ii + 1 i = nbno (ibpl, ii) j = nbno (ibpl, jj) DO k = 1, ntr IF (mbtri (k) /= 0) then k1 = ml (k) k2 = mm (k) k3 = mn (k) IF (i == k2 .and. j == k1) then nbeln (nbpln) = nbeln (nbpln) + 1 !b xxx check against max nbelf (nbpln, nbeln (nbpln) ) = 1 nbelt (nbpln, nbeln (nbpln) ) = k ELSEIF (i == k3 .and. j == k2) then nbeln (nbpln) = nbeln (nbpln) + 1 !b xxx check against max nbelf (nbpln, nbeln (nbpln) ) = 2 nbelt (nbpln, nbeln (nbpln) ) = k ELSEIF (i == k1 .and. j == k3) then nbeln (nbpln) = nbeln (nbpln) + 1 !b xxx check against max nbelf (nbpln, nbeln (nbpln) ) = 3 nbelt (nbpln, nbeln (nbpln) ) = k END IF END IF END DO END DO ELSE DO ii = 1, npin6 (ibpl) - 1 jj = ii + 1 !b xxx check against max i = nbno6 (ibpl, ii) j = nbno6 (ibpl, jj) DO k = 1, ntr IF (mbqua (k) /= 0) then k1 = lnod4 (k, 1) k2 = lnod4 (k, 2) k3 = lnod4 (k, 3) k4 = lnod4 (k, 4) IF (i == k2 .and. j == k1) then nbeln (nbpln) = nbeln (nbpln) + 1 !b xxx check against max nbelf (nbpln, nbeln (nbpln) ) = 1 nbelt (nbpln, nbeln (nbpln) ) = k ELSEIF (i == k3 .and. j == k2) then nbeln (nbpln) = nbeln (nbpln) + 1 !b xxx check against max nbelf (nbpln, nbeln (nbpln) ) = 2 nbelt (nbpln, nbeln (nbpln) ) = k ELSEIF (i == k4 .and. j == k3) then nbeln (nbpln) = nbeln (nbpln) + 1 !b xxx check against max nbelf (nbpln, nbeln (nbpln) ) = 3 nbelt (nbpln, nbeln (nbpln) ) = k ELSEIF (i == k1 .and. j == k4) then nbeln (nbpln) = nbeln (nbpln) + 1 !b xxx check against max nbelf (nbpln, nbeln (nbpln) ) = 4 nbelt (nbpln, nbeln (nbpln) ) = k END IF END IF END DO END DO END IF 301 END DO ! ! *** minimize the bandwidth/profile ! IF (ifrnt == 0) then ! renumber the nodes only CALL promin (melem, MNPEL, mpoin, Ntot, Ntr, NNpe, & nconc, NEWNO, npold, npnew) !b print *, 'Re-ordered nodes' !b do i = 1, Ntr !b !b print *,'819 i, nconc ', i, nconc (i, 1:nnpe) !b !b end do !b DO ipoin = 1, ntot newnod = newno (ipoin) coord (newnod, 1) = xfin (ipoin) coord (newnod, 2) = yfin (ipoin) END DO DO ibpl = 1, nbpln kelm = nbeln (ibpl) DO j = 1, kelm nbele (ibpl, j) = nbelt (ibpl, j) END DO knod = nbpoi (ibpl) DO j = 1, knod jnod = nbret (ibpl, j) nbref (ibpl, j) = NEWNO (jnod) END DO END DO END IF ! ! *** minimize the frontwidth ! IF (ifrnt == 1) then ! renumber the elements only CALL fromin (melem, MNPEL, Ntot, Ntr, NNpe, nconc, & Mtype, NEWEL, nfold, nfnew) DO ibpl = 1, nbpln kelm = nbeln (ibpl) DO j = 1, kelm lelm = nbelt (ibpl, j) nbele (ibpl, j) = NEWEL (lelm) END DO knod = nbpoi (ibpl) DO j = 1, knod nbref (ibpl, j) = nbret (ibpl, j) END DO END DO DO i = 1, ntot coord (i, 1) = xfin (i) coord (i, 2) = yfin (i) END DO END IF ! ! *** minimize the frontwidth and the bandwidth/profile ! IF (ifrnt == 2) then ! renumber both nodes and elements CALL fromin (melem, MNPEL, Ntot, Ntr, NNpe, nconc, & Mtype, NEWEL, nfold, nfnew) DO ibpl = 1, nbpln kelm = nbeln (ibpl) DO j = 1, kelm lelm = nbelt (ibpl, j) nbele (ibpl, j) = NEWEL (lelm) END DO END DO CALL promin (melem, MNPEL, mpoin, Ntot, Ntr, NNpe, & nconc, NEWNO, npold, npnew) DO ipoin = 1, ntot newnod = newno (ipoin) coord (newnod, 1) = xfin (ipoin) coord (newnod, 2) = yfin (ipoin) END DO DO ibpl = 1, nbpln knod = nbpoi (ibpl) DO j = 1, knod jnod = nbret (ibpl, j) nbref (ibpl, j) = NEWNO (jnod) END DO END DO END IF !!b IF (iprog == 0) return !! *** output boundary nodes, elements, and faces ! WRITE (16, * ) nbpln ! number of boundary sides ! DO i = 1, nbpln ! WRITE (16, * ) nbpoi (i), nbeln (i) ! # pts, elem on side i ! WRITE (16, * ) (nbref (i, j), j = 1, nbpoi (i) ) ! nodes on side ! DO k = 1, nbeln (i) ! WRITE (16, * ) nbele (i, k), nbelf (i, k) ! el & face on side ! END DO ! END DO ! !! *** ouput mesh ! I_ZERO = 0 ; I_ONE = 1 !b ! WRITE (17, * ) 'a_comment_for_finding_the_start' ! WRITE (17, '(3i6)') ntot, ntr, nnpe ! nodes, elems, node/el ! DO j = 1, ntot ! node x, y ! WRITE (17, '(i6,f15.6,2x,f15.6)') j, coord (j, 1:2) ! WRITE (21, '(i6, 2(1PE13.5))') I_ZERO, coord (j, 1:2) !b ! END DO ! PRINT *, 'NOTE: NODE FLAG & COORDINATES SAVED TO msh_bc_xyz.dat' ! DO j = 1, ntr ! elem mat connectivity ! WRITE (17, * ) j, mtype (j), (nconc (j, k), k = 1, nnpe) ! IF ( nnpe <= 4 ) then !b ! WRITE (22, '(12I6)') I_ONE, (nconc (j, k), k = 1, nnpe)!b ! ELSE ! WRITE (22, '(12I6)') I_ONE, nconc (j, 1), nconc (j, 3), & !b ! nconc (j, 5), nconc (j, 2), nconc (j, 4), nconc (j, 6) !b ! END IF ! END DO ! PRINT *, 'NOTE: ELEMENT CONNECTIVITY SAVED TO msh_typ_nodes.dat' END SUBROUTINE mesh2d ! ! ******************************************************************* ! ! ** subroutine input ! ! ******************************************************************* ! SUBROUTINE input (n, x, y, sz, s, nbnds, natrib, iatt, incr, & ipol, nnb, nrfnnb, tolrn, ndpt, xd, yd, & vald, nnpe, ord, nbpl, nbnp, npinb, nbno, & mtx, matt, mreg, mden, iadap, ifrnt, & nitsm, iprog) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine generates node points from given boundary info ! DIMENSION x (mtx), y (mtx), sz (mtx), nnb (mreg), iatt (matt), & nrfnnb (mtx, mreg) !b M_SEG , M_SEG DIMENSION nrfb (989), nspt (189), nrfs (489, 189), & xd (mden), yd (mden), vald (mden) !b M_SEG M_SEG mreg, M_SEG DIMENSION npinb (189), nbno (189, 489), nsdm (108, 100), & nsdms (108), nstyp (189) !b mreg M_SEG s = 1.0d0 incr = 1 ipol = 2 ord = 1.0d0 IF (iprog == 1) then !b print *,'11 872 nbnds, natrib, tolrn, ifrnt, nitsm, iadap' !b READ (11, * ) nbnds, natrib, tolrn, ifrnt, nitsm, iadap ifrnt = 0 ; iadap = 0 !b add nnpe here !b READ (11, * , IOSTAT = I_O) nbnds, natrib, tolrn, nitsm IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR nbnds, natrib, tolrn, nits' WRITE (7, *) 'nbnds, natrib, tolrn, nitsm = ', nbnds, natrib, tolrn, nitsm !b print *,nbnds, natrib, tolrn, ifrnt, nitsm, iadap !b print *,'nbnds, natrib, nnpe, tolrn', nbnds, natrib, nnpe, tolrn !b print *,'12 874' !b READ (12, * ) ndpt, (xd (i), yd (i), vald (i), i = 1, ndpt) READ (12, * ) ndpt !b WRITE (7, *) 'NUMBER OF DENSITY POINTS = ', ndpt !b READ (12, * ) (xd (i), yd (i), vald (i), i = 1, ndpt) !b DO I = 1, ndpt !b READ (12, * ) xd (i), yd (i), vald (i) !b END DO !b print *,'ndpt ', xd(ndpt), yd(ndpt), vald(ndpt) IF (nnpe == 4) s = s * 2.0d0 sb = s ELSE !b print *,'11 878' READ (11, * , IOSTAT = I_O) nbnds, natrib, nnpe, tolrn IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR nbnds, natrib, nnpe, tol' s = 1.0d0 IF (nnpe == 4) s = 2.0d0 sb = s IF (nnpe == 3 .or. nnpe == 6) sb = 0.9d0 * s IF (iadap == 0) then !b print *,'14 855 density' !b READ (14, * ) ndpt, (xd (i), yd (i), vald (i), i = 1, ndpt) READ (12, * ) ndpt, (xd (i), yd (i), vald (i), i = 1, ndpt) END IF END IF ! n = 0 nseg1 = 0 nseg = 0 nstage = 1 IF (nbnds == 0) GO TO 10 ! No holes DO i = 1, nbnds ! Define each hole isdm = i isegn = 0 nbnod = 0 !b print *,'11 897' READ (11, * , IOSTAT = I_O) noseg ! Hole region number IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR noseg IN HOLE ', i CALL read (nstage, noseg, sb, n, x, y, sz, nseg, nspt, nrfs, & nbnod, nrfb, nseg1, ndpt, xd, yd, vald, tolrn, ord, mtx, mden,& iadap, isdm, isegn, nsdm, nstyp) nsdms (i) = isegn nnb (i) = nbnod DO j = 1, nbnod nrfnnb (j, i) = nrfb (j) END DO END DO nseg1 = nseg 10 nstage = 2 ntotb = nbnds + natrib ! natrib >= 1 always nstg2 = nbnds + 1 !b print *,'11 912' READ (11, * , IOSTAT = I_O) (iatt (j), j = 1, natrib) ! material of each region IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR iatt', j WRITE (7, *) 'Region material codes = ',iatt (1:natrib) DO i = nstg2, ntotb isdm = i isegn = 0 nbnod = 0 !b print *,'11 919' READ (11, * , IOSTAT = I_O) noseg ! number of sides on region i IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR noseg' WRITE (7, *) 'number of sides on region ', i, ' is ', noseg CALL read (nstage, noseg, sb, n, x, y, sz, nseg, nspt, nrfs, & nbnod, nrfb, nseg1, ndpt, xd, yd, vald, tolrn, ord, mtx, mden, & iadap, isdm, isegn, nsdm, nstyp) nsdms (i) = isegn nnb (i) = nbnod DO j = 1, nbnod nrfnnb (j, i) = nrfb (j) END DO END DO ! ! *** copy segment information into arrays for boundary conditions ! nbpl = nseg DO i = 1, nseg npinb (i) = nspt (i) DO j = 1, npinb (i) nbno (i, j) = nrfs (j, i) END DO IF (nstyp (i) == 4) npinb (i) = - npinb (i) END DO END SUBROUTINE input ! ! ******************************************************************** ! ! ** subroutine read ! ! ******************************************************************** ! SUBROUTINE read (nst, noseg, s, n, x, y, sz, nseg, nspt, nrfs, & nbnod, nrfb, nseg1, ndpt, xd, yd, vald, tolrn, ord, mtx, mden, & iadap, isdm, isegn, nsdm, nstyp) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine reads domain and attribute boundaries ! DIMENSION x (mtx), y (mtx), nrfb (989), nspt (189), & nrfs (489, 189), sz (mtx), xd (mden), yd (mden), vald (mden) !b DIMENSION xx (69), yy (69), nsdm (108, 100), nstyp (189) !b M_USER M_USER lastyp = 0 ; ifstyp = 0 ; ifstsg = 0 ; lastsg = 0 ! initialize DO j = 1, noseg ! number of sides for this region npt = 0 !b print *,'11 965' READ (11, * , IOSTAT = I_O) ityp ! type 1 line 2 arc+- 3 user 4 other side IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR ityp ', j, noseg WRITE (7, *) 'Segment ', j, ' line type ', ityp IF ( ityp < 1 .OR. ityp > 4 ) STOP 'h_geom invalid line type' IF (j == 1) then ! first side IF (ityp == 1) then ! first line READ (11, * , IOSTAT = I_O) x1, y1, x2, y2 ! first line IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR x1, y1, x2, y2 ', j firstx = x1 ; firsty = y1 prevx = x2 ; prevy = y2 END IF IF (ityp == 2) then ! first arc READ (11, * , IOSTAT = I_O) ax, ay, bx, by, rad ! first arc IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR ax, ay, bx, by, rad ', j firstx = ax ; firsty = ay prevx = bx ; prevy = by END IF IF (ityp == 3) then ! first user curve READ (11, * ) numpt, (xx (k), yy (k), k = 1, numpt) firstx = xx (1) ; firsty = yy (1) prevx = xx (numpt) ; prevy = yy (numpt) END IF IF (ityp == 4) then ! previously defined side READ (11, * ) iseg npts = nspt (iseg) nod1 = nrfs (1, iseg) nod2 = nrfs (npts, iseg) firstx = x (nod1) ; firsty = y (nod1) prevx = x (nod2) ; prevy = y (nod2) END IF ELSEIF (j == noseg) then ! last side IF (ityp == 1) then ! last line x1 = prevx ; y1 = prevy x2 = firstx ; y2 = firsty END IF IF (ityp == 2) then ! last arc ax = prevx ; ay = prevy bx = firstx ; by = firsty READ (11, * , IOSTAT = I_O) rad IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR rad ', j END IF IF (ityp == 3) then ! last user READ (11, * ) numpt, (xx (k), yy (k), k = 2, numpt - 1) xx (1) = prevx ; yy (1) = prevy xx (numpt) = firstx ; yy (numpt) = firsty END IF IF (ityp == 4) then ! last side READ (11, * ) iseg END IF ELSE IF (ityp == 1) then ! intermediate line x1 = prevx ; y1 = prevy READ (11, * , IOSTAT = I_O) x2, y2 IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR x2, y2 ', j prevx = x2 ; prevy = y2 END IF IF (ityp == 2) then ! intermediate arc ax = prevx ; ay = prevy READ (11, * , IOSTAT = I_O) bx, by, rad IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR bx, by, rad ', j prevx = bx ; prevy = by END IF IF (ityp == 3) then ! user defined group of points xx (1) = prevx ; yy (1) = prevy !b xxx must have numpt <= M_USER READ (11, * ) numpt, (xx (k), yy (k), k = 2, numpt) prevx = xx (numpt) ; prevy = yy (numpt) END IF IF (ityp == 4) then ! intermediate user READ (11, * ) iseg npts = nspt (iseg) nod2 = nrfs (npts, iseg) prevx = x (nod2) ; prevy = y (nod2) END IF END IF IF (ityp == 1) then CALL segdiv1 (x1, y1, x2, y2, s, nbnod, nrfb, j, noseg, n, & x, y, sz, nst, nseg, nseg1, nspt, nrfs, & ifstsg, lastsg, ifstyp, ndpt, xd, yd, vald, & tolrn, lastyp, ord, mtx, mden, iadap) isegn = isegn + 1 !b xxx check against max nsdm (isdm, isegn) = nseg nstyp (nseg) = ityp END IF IF (ityp == 2) then CALL segdiv2 (ax, ay, bx, by, rad, s, nbnod, nrfb, j, noseg, & n, x, y, sz, nst, nseg, nseg1, nspt, nrfs, & ifstsg, lastsg, ifstyp, ndpt, xd, yd, vald, & tolrn, lastyp, ord, mtx, mden, iadap) isegn = isegn + 1 !b xxx check against max nsdm (isdm, isegn) = nseg nstyp (nseg) = ityp END IF IF (ityp == 3) then ! user defined set of points isegn = isegn + 1 !b xxx check against max nsdm (isdm, isegn) = nseg nseg = nseg + 1 !b xxx check against max nspt (nseg) = numpt nstyp (nseg) = ityp kst = 1 numpt1 = numpt - 1 !b xxx check against max IF (j > 1) then kst = 2 N = n + 1 !b xxx check against max nbnod = nbnod+1 !b xxx check against max x (n) = xx (1) y (n) = yy (1) nrfb (nbnod) = n nrfs (1, nseg) = n END IF DO k = kst, numpt1 nbnod = nbnod+1 !b xxx check against max n = n + 1 !b xxx check against max nrfb (nbnod) = n nrfs (k, nseg) = n x (n) = xx (k) y (n) = yy (k) xn = x (n) yn = y (n) CALL densit (xn, yn, xd, yd, vald, ndpt, dens, & tolrn, s, ord, mden, iadap) sz (n) = s * dens END DO IF (j /= noseg) then nrfs (numpt, nseg) = n + 1 ELSE nod = nrfs (1, ifstsg) IF (ifstyp /= 4) then npts = nspt (ifstsg) nod = nrfs (npts, ifstsg) END IF nrfs (numpt, nseg) = nod END IF IF (lastyp == 4) then npts = nspt (lastsg) nod = nrfs (npts, lastsg) l = n - nspt (nseg) + 3 nrfs (1, nseg) = nod nbnod = nbnod-nspt (nseg) + 1 ncoun = 1 DO k = l, n ncoun = ncoun + 1 !b xxx check against max nrfs (ncoun, nseg) = k - 1 nbnod = nbnod+1 nrfb (nbnod) = k - 1 x (k - 1) = x (k) y (k - 1) = y (k) sz (k - 1) = sz (k) END DO IF (j /= noseg) nrfs (nspt (nseg), nseg) = n n = n - 1 END IF IF (nst == 2 .and. j == 1) then t = tolrn * 1.0d-2 DO k = 1, n IF (k /= nrfb (1) ) then IF (dabs(x (k) - x1) < t .and. & dabs(y (k) - y1) < t) then nod = k GO TO 191 END IF END IF END DO GO TO 192 191 l = n - nspt (nseg) + 3 nrfs (1, nseg) = nod nbnod = 1 nrfb (nbnod) = nod ncoun = 1 DO k = l, n ncoun = ncoun + 1 !b xxx check against max nrfs (ncoun, nseg) = k - 1 nbnod = nbnod+1 nrfb (nbnod) = k - 1 x (k - 1) = x (k) y (k - 1) = y (k) sz (k - 1) = sz (k) END DO nrfs (nspt (nseg), nseg) = n n = n - 1 END IF 192 CALL invseg (nseg, nspt, nrfs) END IF IF (ityp == 4) then isegn = isegn + 1 !b xxx check against max nsdm (isdm, isegn) = iseg nstyp (iseg) = ityp kst = 1 npts = nspt (iseg) nod = nrfs (1, iseg) IF (lastyp /= 0 .and. lastyp /= 4) then nrfs (1, lastsg) = nod nbnod = nbnod+1 nrfb (nbnod) = nod END IF IF (j > 1) kst = 2 kpt = npts IF (j == noseg) kpt = npts - 1 DO k = kst, kpt nbnod = nbnod+1 nod = nrfs (k, iseg) nrfb (nbnod) = nod END DO END IF IF (j == 1) then ifstsg = nseg IF (ityp == 4) ifstsg = iseg ifstyp = ityp END IF lastyp = ityp lastsg = nseg IF (ityp == 4) lastsg = iseg END DO END SUBROUTINE read ! ! ******************************************************************** ! ! ** subroutine segdiv1 ! ! ******************************************************************** ! SUBROUTINE segdiv1 (x1, y1, x2, y2, s, nbnod, nrfb, j, noseg, & n, x, y, sz, nst, nseg, nseg1, nspt, nrfs, & ifstsg, lastsg, ifstyp, ndpt, xd, yd, vald, & tolrn, lastyp, ord, mtx, mden, iadap) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine divides a st. line segment of boundary into ! *** smaller ones according to the given element size and the local ! *** mesh density. ! DIMENSION x (mtx), y (mtx), nrfb (989), nspt (189), & nrfs (489, 189), sz (mtx), dist (489), xd (mden), & yd (mden), vald (mden) d = 0.d0 nseg = nseg + 1 !b xxx check against max nspt (nseg) = 0 xlen = x2 - x1 ylen = y2 - y1 dl = dsqrt (xlen**2 + ylen**2) temp = huge (temp) !b 99999. IF (xlen /= 0.) temp = ylen / xlen n = n + 1 !b xxx check against max nbnod = nbnod+1 nrfb (nbnod) = n x (n) = x1 y (n) = y1 nspt (nseg) = nspt (nseg) + 1 CALL densit (x1, y1, xd, yd, vald, ndpt, dens, tolrn, s, ord, & mden, iadap) div1 = dens * s sz (n) = div1 CALL divis (xlen, ylen, temp, div1, tdx, tdy) txn = x1 + tdx tyn = y1 + tdy CALL densit (txn, tyn, xd, yd, vald, ndpt, dens, tolrn, s, ord, & mden, iadap) div2 = dens * s IF (div2 < div1) then div = (div1 * div1) / (2.d0 * div1 - div2) ELSE div = div2 END IF d = d+div IF (d >= dl) GO TO 30 10 dist (nspt (nseg) ) = div CALL divis (xlen, ylen, temp, div, xdiv, ydiv) n = n + 1 !b xxx check against max x (n) = x (n - 1) + xdiv y (n) = y (n - 1) + ydiv nspt (nseg) = nspt (nseg) + 1 xn = x (n) yn = y (n) CALL densit (xn, yn, xd, yd, vald, ndpt, dens, tolrn, s, ord, & mden, iadap) div1 = dens * s CALL divis (xlen, ylen, temp, div1, tdx, tdy) txn = xn + tdx tyn = yn + tdy CALL densit (txn, tyn, xd, yd, vald, ndpt, dens, tolrn, s, ord, & mden, iadap) div2 = dens * s IF (div2 < div1) then div = (div1 * div1) / (2.d0 * div1 - div2) ELSE div = div2 END IF d = d+div IF (d >= dl) GO TO 20 GO TO 10 20 CONTINUE dist (nspt (nseg) ) = dl - d+div CALL densit (x2, y2, xd, yd, vald, ndpt, dens, tolrn, s, ord, & mden, iadap) divj = dens * s endiv = dist (nspt (nseg) ) IF (endiv > divj / 2.) then err = divj - endiv ELSE err = (divj - dist (nspt (nseg) - 1) ) - endiv n = n - 1 nspt (nseg) = nspt (nseg) - 1 END IF dist (nspt (nseg) ) = divj DO k = 1, nspt (nseg) - 1 l = n - nspt (nseg) + k m = l + 1 !b xxx check against max dist (k) = dist (k) - (dist (k) / (dl + err) ) * err div = dist (k) CALL divis (xlen, ylen, temp, div, xdv, ydv) x (m) = x (l) + xdv y (m) = y (l) + ydv xm = x (m) ym = y (m) CALL densit (xm, ym, xd, yd, vald, ndpt, dens, tolrn, s, ord, & mden, iadap) sz (m) = s * dens nrfs (k, nseg) = l nbnod = nbnod+1 nrfb (nbnod) = m END DO 30 nrfs (nspt (nseg), nseg) = n nspt (nseg) = nspt (nseg) + 1 IF (j /= noseg) then nrfs (nspt (nseg), nseg) = n + 1 ELSE nod = nrfs (1, ifstsg) IF (ifstyp /= 4) then npts = nspt (ifstsg) nod = nrfs (npts, ifstsg) END IF nrfs (nspt (nseg), nseg) = nod END IF IF (lastyp == 4) then npts = nspt (lastsg) nod = nrfs (npts, lastsg) l = n - nspt (nseg) + 3 nrfs (1, nseg) = nod nbnod = nbnod-nspt (nseg) + 1 ncoun = 1 DO k = l, n ncoun = ncoun + 1 !b xxx check against max nrfs (ncoun, nseg) = k - 1 nbnod = nbnod+1 nrfb (nbnod) = k - 1 x (k - 1) = x (k) y (k - 1) = y (k) sz (k - 1) = sz (k) END DO IF (j /= noseg) nrfs (nspt (nseg), nseg) = n n = n - 1 END IF IF (nst == 2 .and. j == 1) then t = tolrn * 1.0d-2 DO k = 1, n IF (k /= nrfb (1) ) then IF (dabs (x (k) - x1) < t .and. & dabs (y (k) - y1) < t) then nod = k GO TO 191 END IF END IF END DO GO TO 192 191 l = n - nspt (nseg) + 3 nrfs (1, nseg) = nod nbnod = 1 nrfb (nbnod) = nod ncoun = 1 DO k = l, n ncoun = ncoun + 1 !b xxx check against max nrfs (ncoun, nseg) = k - 1 nbnod = nbnod+1 nrfb (nbnod) = k - 1 x (k - 1) = x (k) y (k - 1) = y (k) sz (k - 1) = sz (k) END DO nrfs (nspt (nseg), nseg) = n n = n - 1 END IF 192 CALL invseg (nseg, nspt, nrfs) END SUBROUTINE segdiv1 ! ! ******************************************************************** ! ! ** subroutine segdiv2 ! ! ******************************************************************** ! SUBROUTINE segdiv2 (ax, ay, bx, by, r, s, nbnod, nrfb, j, noseg,& n, x, y, sz, nst, nseg, nseg1, nspt, nrfs, & ifstsg, lastsg, ifstyp, ndpt, xd, yd, vald, & tolrn, lastyp, ord, mtx, mden, iadap) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine devides a circular segment of boundary into ! *** smaller ones according to the given element size and local ! *** mesh density. ! DIMENSION x (mtx), y (mtx), nrfb (989), sz (mtx), ang (489), & xd (mden), yd (mden), vald (mden) DIMENSION nspt (189), nrfs (489, 189) nseg = nseg + 1 !b xxx check against max nspt (nseg) = 0 angl = 0.d0 d = dsqrt ( (ax - bx) **2 + (ay - by) **2) cx1 = (bx - ax) / d cy1 = (by - ay) / d CALL newpt (ax, ay, cx, cy, cx1, cy1, d, r, r) IF (r < 0.d0) r = - r cox = (ax - cx) / r six = (ay - cy) / r IF ( (d / r) >= 2.d0) then WRITE (13, 501) STOP ELSE theta = dacos (1.d0 - 0.5d0 * (d / r) **2) END IF a2 = ax * (by - cy) - bx * (ay - cy) + cx * (ay - by) nbnod = nbnod+1 n = n + 1 !b xxx check against max nrfb (nbnod) = n x (n) = ax y (n) = ay nspt (nseg) = nspt (nseg) + 1 CALL densit (ax, ay, xd, yd, vald, ndpt, dens, tolrn, s, & ord, mden, iadap) chrd1 = dens * s sz (n) = chrd1 IF (chrd1 >= d) chrd1 = d div1 = dacos (1.d0 - 0.5d0 * (chrd1 / r) **2) divn = angl + div1 IF (a2 < 0.d0) divn = - divn xp = r * dcos (divn) yp = r * dsin (divn) txn = (cox * xp - six * yp) + cx tyn = (six * xp + cox * yp) + cy CALL densit (txn, tyn, xd, yd, vald, ndpt, dens, tolrn, s, & ord, mden, iadap) chrd2 = dens * s IF (chrd2 < chrd1) then chrd = (chrd1 * chrd1) / (2.d0 * chrd1 - chrd2) ELSE chrd = chrd2 END IF IF (chrd >= d) chrd = d div = dacos (1.d0 - 0.5 * (chrd / r) **2) angl = angl + div IF (angl >= theta) GO TO 30 10 ang (nspt (nseg) ) = div n = n + 1 !b xxx check against max divn = angl IF (a2 < 0.d0) divn = - divn xp = r * dcos (divn) yp = r * dsin (divn) x (n) = (cox * xp - six * yp) + cx y (n) = (six * xp + cox * yp) + cy nspt (nseg) = nspt (nseg) + 1 xn = x (n) yn = y (n) CALL densit (xn, yn, xd, yd, vald, ndpt, dens, tolrn, s, & ord, mden, iadap) chrd1 = dens * s IF (chrd1 >= d) chrd1 = d div1 = dacos (1.d0 - 0.5 * (chrd1 / r) **2) divn = angl + div1 IF (a2 < 0.d0) divn = - divn xp = r * dcos (divn) yp = r * dsin (divn) txn = (cox * xp - six * yp) + cx tyn = (six * xp + cox * yp) + cy CALL densit (txn, tyn, xd, yd, vald, ndpt, dens, tolrn, s, & ord, mden, iadap) chrd2 = dens * s IF (chrd2 < chrd1) then chrd = (chrd1 * chrd1) / (2.d0 * chrd1 - chrd2) ELSE chrd = chrd2 END IF IF (chrd >= d) chrd = d div = dacos (1.d0 - 0.5 * (chrd / r) **2) angl = angl + div IF (angl >= theta) GO TO 20 GO TO 10 20 CONTINUE ang (nspt (nseg) ) = theta - angl + div CALL densit (bx, by, xd, yd, vald, ndpt, dens, tolrn, s, & ord, mden, iadap) chrd = dens * s IF (chrd >= d) chrd = d divj = dacos (1.d0 - 0.5 * (chrd / r) **2) endiv = ang (nspt (nseg) ) IF (endiv >= divj / 2.) then err = divj - endiv ELSE err = divj - ang (nspt (nseg) - 1) - endiv n = n - 1 nspt (nseg) = nspt (nseg) - 1 END IF ang (nspt (nseg) ) = divj angl = 0.d0 DO k = 1, nspt (nseg) - 1 l = n - nspt (nseg) + k m = l + 1 !b xxx check against max ang (k) = ang (k) - (ang (k) / (theta + err) ) * err angl = angl + ang (k) divn = angl IF (a2 < 0.d0) divn = - divn xp = r * dcos (divn) yp = r * dsin (divn) x (m) = (cox * xp - six * yp) + cx y (m) = (six * xp + cox * yp) + cy xm = x (m) ym = y (m) CALL densit (xm, ym, xd, yd, vald, ndpt, dens, tolrn, s, & ord, mden, iadap) sz (m) = s * dens nrfs (k, nseg) = l nbnod = nbnod+1 nrfb (nbnod) = m END DO 30 nrfs (nspt (nseg), nseg) = n nspt (nseg) = nspt (nseg) + 1 IF (j /= noseg) then nrfs (nspt (nseg), nseg) = n + 1 ELSE nod = nrfs (1, ifstsg) IF (ifstyp /= 4) then npts = nspt (ifstsg) nod = nrfs (npts, ifstsg) END IF nrfs (nspt (nseg), nseg) = nod END IF IF (lastyp == 4) then npts = nspt (lastsg) nod = nrfs (npts, lastsg) l = n - nspt (nseg) + 3 nrfs (1, nseg) = nod nbnod = nbnod-nspt (nseg) + 1 ncoun = 1 DO k = l, n ncoun = ncoun + 1 !b xxx check against max nrfs (ncoun, nseg) = k - 1 nbnod = nbnod+1 !b xxx check against max nrfb (nbnod) = k - 1 x (k - 1) = x (k) y (k - 1) = y (k) sz (k - 1) = sz (k) END DO IF (j /= noseg) nrfs (nspt (nseg), nseg) = n n = n - 1 END IF IF (nst == 2 .and. j == 1) then t = tolrn * 1.0d-2 DO k = 1, n IF (k /= nrfb (1) ) then print *, 'WARNING: segdiv2, unverified use of ax,ay occurred' !b IF (dabs (x (k) - x1) < t .and. & !b xxx y1 UNDEFINED !b dabs (y (k) - y1) < t) then !b xxx y1 UNDEFINED IF (dabs (x (k) - ax) < t .and. & !b xxx y1 UNDEFINED dabs (y (k) - ay) < t) then !b xxx y1 UNDEFINED nod = k GO TO 191 END IF END IF END DO GO TO 192 191 l = n - nspt (nseg) + 3 nrfs (1, nseg) = nod nbnod = 1 nrfb (nbnod) = nod ncoun = 1 DO k = l, n ncoun = ncoun + 1 !b xxx check against max nrfs (ncoun, nseg) = k - 1 nbnod = nbnod+1 nrfb (nbnod) = k - 1 x (k - 1) = x (k) y (k - 1) = y (k) sz (k - 1) = sz (k) END DO nrfs (nspt (nseg), nseg) = n n = n - 1 END IF 192 CALL invseg (nseg, nspt, nrfs) 501 FORMAT(//10x,'***************program stopped****************' & & /10x,'arc segments of 180 degrees or more prohibited'//) PRINT *,'Exiting SUBROUTINE segdiv2' END SUBROUTINE segdiv2 ! !********************************************************************* ! ! ** subroutine sort ! ! ******************************************************************** ! SUBROUTINE sort (npoin, x, y, xmax, xmin, ymax, ymin, mtx) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine extracts from the given data points ! *** the maximum and minimum coordinate values. ! DIMENSION x (mtx), y (mtx) xmin = x (1) xmax = x (1) ymin = y (1) ymax = y (1) DO i = 2, npoin IF (x (i) > xmax) xmax = x (i) IF (x (i) < xmin) xmin = x (i) IF (y (i) > ymax) ymax = y (i) IF (y (i) < ymin) ymin = y (i) END DO END SUBROUTINE sort ! ! ******************************************************************** ! ! ** subroutine invseg ! ! ******************************************************************** ! SUBROUTINE invseg (nseg, nspt, nrfs) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine inverts the order of the segments recorded ! *** for attribute boundaries ! !b M_SEG M_SEG DIMENSION nspt (189), nrfs (489, 189), ntemp (489) npts = nspt (nseg) !b xxx check against max DO i = 1, npts ntemp (i) = nrfs (i, nseg) END DO j = npts + 1 DO i = 1, npts k = j - i nrfs (k, nseg) = ntemp (i) END DO END SUBROUTINE invseg ! ! ******************************************************************** ! ! ** subroutine reject ! ! ******************************************************************** ! SUBROUTINE reject (ir, n, iflag, x, y, sz, iflwh, np, px, py, & pr, icrsp, xm, ym, rm, xmax, xmin, ymax, & ymin, iopt, mtx, mxin, mtem) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine rejects the points external to domain boundaries ! DIMENSION xm (mtem), ym (mtem), rm (mtem), x (mtx), y (mtx), & sz (mtx), iflwh (mxin), px (mxin), py (mxin), & pr (mxin), iflag (mxin), icrsp (mxin) DO i = 1, np iflag (i) = 1 END DO DO i = 1, np xi = px (i) yi = py (i) IF (xi > xmax .or. xi < xmin) GO TO 10 IF (yi > ymax .or. yi < ymin) GO TO 10 DO j = 1, n xj = x (j) yj = y (j) sj = sz (j) dist1 = dabs (xi - xj) dist2 = dabs (yi - yj) dist = dsqrt ( (xi - xj) **2 + (yi - yj) **2) ! if (dist1 <= 0.3161*sj .and. dist2 <= 0.3161*sj) GO TO 10 ! if (dist < 0.707*sj) GO TO 10 IF (dist < 0.3 * sj) GO TO 10 END DO CALL inorout (xi, yi, n, x, y, ichk, mtx, iangl) IF (iopt == 1 .and. ichk == 0) iflag (i) = 0 IF (iopt == 1 .and. ichk == 1) iflwh (icrsp (i) ) = ir IF (iopt == 2 .and. ichk == 1) iflag (i) = 0 IF (iangl == 1) iflag (i) = 1 10 CONTINUE END DO m = 0 DO i = 1, np IF (iflag (i) == 1) GO TO 20 m = m + 1 !b xxx check against max xm (m) = px (i) ym (m) = py (i) rm (m) = pr (i) icrsp (m) = icrsp (i) 20 CONTINUE END DO np = m DO i = 1, np px (i) = xm (i) py (i) = ym (i) pr (i) = rm (i) END DO END SUBROUTINE reject ! ! ******************************************************************** ! ! ** subroutine inorout ! ! ******************************************************************** ! SUBROUTINE inorout (xi, yi, n, x, y, ichk, mtx, iangl) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine checks whether a given point lies inside or ! *** outside a polygon. ! DIMENSION x (mtx), y (mtx), d (999), dk (999), a2 (999) ichk = 0 iangl = 0 anglm = - 0.71d0 tupi = 2.d0 * 3.1415927d0 sum = 0.d0 DO j = 1, n k = j + 1 !b xxx check against max IF (j == n) k = 1 d (j) = dsqrt ( (x (j) - xi) **2 + (y (j) - yi) **2) dk (j) = dsqrt ( (x (j) - x (k) ) **2 + (y (j) - y (k) ) **2) a2 (j) = x (j) * (y (k) - yi) - x (k) * (y (j) - yi) & + xi * (y (j) - y (k) ) END DO DO j = 1, n k = j + 1 IF (j == n) k = 1 temp1 = (d (j) * d (j) + d (k) * d (k) - dk (j) * dk (j) ) & / (2.d0 * d (j) * d (k) ) IF (temp1 > 1.d0) temp2 = 0.d0 IF (temp1 < - 1.d0) temp2 = 3.1415927d0 IF (temp1 < 1.d0 .and. temp1 > - 1.d0) temp2 = dacos (temp1) IF (temp1 < anglm) iangl = 1 IF (a2 (j) < 0.) temp2 = - temp2 ! if(dabs(a2(j)) < 1.0d-9) iangl=1 ! if (iangl == 1) GO TO 30 sum = sum + temp2 END DO temp = dabs (sum) - dabs (tupi) IF (dabs (temp) > 0.0001) GO TO 10 ichk = 1 GO TO 30 10 IF (dabs (sum) > 0.0001) then WRITE (13, 200) sum 200 FORMAT (//10x,'second subroutine entered as - sum =',f5.3, & & /10x,'does not equal either 2*pi or zero'//) CALL inorou2 (xi, yi, n, x, y, ichk, mtx) END IF 30 CONTINUE END SUBROUTINE inorout ! ! ********************************************************************** ! ! ** subroutine inorou2 ! ! ********************************************************************** ! SUBROUTINE inorou2 (xi, yi, n, x, y, ichk, mtx) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine checks whether a given point lies inside or ! *** outside a polygon. ! DIMENSION x (mtx), y (mtx) ichk = 0 inum = 0 xj = huge (xj) ! 1.0d12 yj = huge (yj) / 10 ! 1.0d11 DO i = 1, n j = i + 1 IF (i == n) j = 1 xk = x (i) yk = y (i) xl = x (j) yl = y (j) CALL intsct (xi, yi, xj, yj, xk, yk, xl, yl, iout) IF (iout == 1) inum = inum + 1 END DO rinum = float (inum) / 2.d0 inum = inum / 2 rem = rinum - float (inum) IF (rem > 0.49) ichk = 1 END SUBROUTINE inorou2 ! ! ********************************************************************** ! ! ** subroutine circum ! ! ********************************************************************** ! SUBROUTINE circum (x, y, xc, yc, r2, icirc) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine calculates the circumcentre and square of the ! *** radius of the circumcirle of a triangle, given the coordinates ! *** of its node points. ! DIMENSION x (3), y (3) icirc = 0 a = (x (1) * x (1) - x (2) * x (2) ) + (y (1) * y (1) & - y (2) * y (2) ) b = (x (1) * x (1) - x (3) * x (3) ) + (y (1) * y (1) & - y (3) * y (3) ) d = (4. * (x (1) - x (2) ) * (y (1) - y (3) ) ) & - (4. * (y (1) - y (2) ) * (x (1) - x (3) ) ) IF (d == 0.0d0) then WRITE (13, * ) 'warning: d = 0.0 in circum' ! write (13,*)' x(1)=',x(1),' x(2)=',x(2),' x(3)=',x(3) ! write (13,*)' y(1)=',y(1),' y(2)=',y(2),' y(3)=',y(3) icirc = 1 RETURN END IF c11 = 2. * (y (1) - y (3) ) / d c12 = 2. * (y (1) - y (2) ) / d c21 = 2. * (x (1) - x (3) ) / d c22 = 2. * (x (1) - x (2) ) / d xc = c11 * a - c12 * b yc = c22 * b - c21 * a r2 = (x (1) - xc) * (x (1) - xc) + (y (1) - yc) * (y (1) - yc) END SUBROUTINE circum ! ! ! ********************************************************************** ! ! ** subroutine check ! ! ********************************************************************** ! SUBROUTINE check (xi, yi, xcj, ycj, r2j, ichk) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine checks whether a new node point lies within the ! *** circumcircle of a triangle or not ! ichk = 0 a = (xi - xcj) * (xi - xcj) b = (yi - ycj) * (yi - ycj) c = a + b ! r2j=0.999999999*r2j r2j = 0.999999999999d0 * r2j IF (c < r2j) ichk = 1 END SUBROUTINE check ! ! ********************************************************************** ! ! ** subroutine inserp ! ! ********************************************************************** ! SUBROUTINE inserp (nrej, l, m, n, l1, l2, npols) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine finds the insertion polynomial for the currently ! *** inserted node. ! DIMENSION l (199), m (199), n (199), a (199), b (199), c (199), & d (199), e (199), f (199), l1 (199), l2 (199) DO i = 1, nrej a (i) = l (i) b (i) = m (i) c (i) = m (i) d (i) = n (i) e (i) = n (i) f (i) = l (i) END DO npols = 0 nrej1 = nrej - 1 IF (nrej1 == 0) GO TO 50 DO i = 1, nrej1 ii = i + 1 DO j = ii, nrej IF (a (i) /= 0 .and. b (i) /= 0) then IF (a (i) == b (j) .and. b (i) == a (j) ) then a (i) = 0 b (i) = 0 a (j) = 0 b (j) = 0 ELSEIF (a (i) == d (j) .and. b (i) == c (j) ) then a (i) = 0 b (i) = 0 c (j) = 0 d (j) = 0 ELSEIF (a (i) == f (j) .and. b (i) == e (j) ) then a (i) = 0 b (i) = 0 e (j) = 0 f (j) = 0 END IF END IF IF (c (i) /= 0 .and. d (i) /= 0) then IF (c (i) == b (j) .and. d (i) == a (j) ) then c (i) = 0 d (i) = 0 a (j) = 0 b (j) = 0 ELSEIF (c (i) == d (j) .and. d (i) == c (j) ) then c (i) = 0 d (i) = 0 c (j) = 0 d (j) = 0 ELSEIF (c (i) == f (j) .and. d (i) == e (j) ) then c (i) = 0 d (i) = 0 e (j) = 0 f (j) = 0 END IF END IF IF (e (i) /= 0 .and. f (i) /= 0) then IF (e (i) == b (j) .and. f (i) == a (j) ) then e (i) = 0 f (i) = 0 a (j) = 0 b (j) = 0 ELSEIF (e (i) == d (j) .and. f (i) == c (j) ) then e (i) = 0 f (i) = 0 c (j) = 0 d (j) = 0 ELSEIF (e (i) == f (j) .and. f (i) == e (j) ) then e (i) = 0 f (i) = 0 e (j) = 0 f (j) = 0 END IF END IF END DO END DO DO i = 1, nrej IF (a (i) /= 0 .and. b (i) /= 0) then npols = npols + 1 !b xxx check against max l1 (npols) = a (i) l2 (npols) = b (i) END IF IF (c (i) /= 0 .and. d (i) /= 0) then npols = npols + 1 !b xxx check against max l1 (npols) = c (i) l2 (npols) = d (i) END IF IF (e (i) /= 0 .and. f (i) /= 0) then npols = npols + 1 !b xxx check against max l1 (npols) = e (i) l2 (npols) = f (i) END IF END DO GO TO 100 50 l1 (1) = a (1) l2 (1) = b (1) l1 (2) = c (1) l2 (2) = d (1) l1 (3) = e (1) l2 (3) = f (1) npols = 3 100 CONTINUE END SUBROUTINE inserp ! ! ********************************************************************** ! ! ** subroutine update ! ! ********************************************************************** ! SUBROUTINE update (x, y, nl, nm, nn, ntri, xc, yc, r2, npols, l1, & l2, nref, i, mx, mtri, icirc) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine updates the triangulation by rejecting the ! *** tringles whose circumcircle contains the currently inserted ! *** node point and joining the node point with the vertices of the ! *** insertion polygon defined by subroutine inserp. it also ! *** calculates the circumcentre and radius(sq.) of the new triangles. ! DIMENSION x (mx), y (mx), nl (mtri), nm (mtri), nn (mtri), & xc (mtri), yc (mtri), r2 (mtri), ax (3), ay (3), nref (199), & l1 (199), l2 (199) CHARACTER string * 4 k = 0 n = npols - 2 string = 'mtri' DO j = 1, npols IF (j <= n) then nt = nref (j) ELSE k = k + 1 nt = ntri + k !b xxx check against max CALL dcheck (nt, mtri, string) END IF nl (nt) = l1 (j) nm (nt) = l2 (j) nn (nt) = i i1 = nl (nt) i2 = nm (nt) i3 = nn (nt) ax (1) = x (i1) ay (1) = y (i1) ax (2) = x (i2) ay (2) = y (i2) ax (3) = x (i3) ay (3) = y (i3) CALL circum (ax, ay, cx, cy, rsq, icirc) ! IF (icirc == 1) then WRITE (13, * ) 'warning: emergency exit from update' RETURN END IF ! xc (nt) = cx yc (nt) = cy r2 (nt) = rsq END DO ntri = ntri + 2 END SUBROUTINE update ! ! ********************************************************************** ! ! ** subroutine final ! ! ********************************************************************** ! SUBROUTINE final (npoin, ntri, nl, nm, nn, noutb, inbnd, ninb, & nodi, iflag, mxin, mx, mtx, mtri, mhol, xb, yb, x, y, xtemp, & ytemp) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine rejects from the final triangulation the external ! *** triangles. ! DIMENSION nl (mtri), nm (mtri), nn (mtri), iflag (mxin), nodi ( & mtx, mhol), ninb (mhol), xb (mtx), yb (mtx), xtemp (mtx), ytemp ( & mtx), x (mx), y (mx) nin = 0 nout = 0 DO 100 i = 1, ntri IF (nl (i) > npoin) GO TO 40 IF (nm (i) > npoin) GO TO 40 IF (nn (i) > npoin) GO TO 40 GO TO 45 40 nout = nout + 1 GO TO 100 45 CONTINUE DO 10 k1 = 1, noutb IF (nl (i) /= k1) GO TO 10 GO TO 11 10 END DO GO TO 50 11 DO 12 k2 = 1, noutb IF (nm (i) /= k2) GO TO 12 GO TO 13 12 END DO GO TO 50 13 DO 14 k3 = 1, noutb IF (nn (i) /= k3) GO TO 14 IF (k1 < k2 .and. k2 < k3) GO TO 101 IF (k2 < k3 .and. k3 < k1) GO TO 101 IF (k1 < k2 .and. k3 < k1) GO TO 101 nout = nout + 1 GO TO 100 14 END DO GO TO 50 101 x1 = xb (k1) y1 = yb (k1) x2 = xb (k2) y2 = yb (k2) x3 = xb (k3) y3 = yb (k3) x12 = 0.5 * (x1 + x2) y12 = 0.5 * (y1 + y2) x23 = 0.5 * (x3 + x2) y23 = 0.5 * (y3 + y2) x31 = 0.5 * (x1 + x3) y31 = 0.5 * (y1 + y3) xp = x1 + (x23 - x1) / 9.d0 yp = y1 + (y23 - y1) / 9.d0 CALL inorout (xp, yp, noutb, xb, yb, ichk, mtx, iangl) IF (ichk == 1) GO TO 102 nout = nout + 1 GO TO 100 102 xq = x2 + (x31 - x2) / 9.d0 yq = y2 + (y31 - y2) / 9.d0 CALL inorout (xq, yq, noutb, xb, yb, ichk, mtx, iangl) IF (ichk == 1) GO TO 103 nout = nout + 1 GO TO 100 103 xr = x3 + (x12 - x3) / 9.d0 yr = y3 + (y12 - y3) / 9.d0 CALL inorout (xr, yr, noutb, xb, yb, ichk, mtx, iangl) IF (ichk == 1) GO TO 90 nout = nout + 1 GO TO 100 50 IF (inbnd == 0) GO TO 90 DO 30 j = 1, inbnd nj = ninb (j) DO 31 j1 = 1, nj IF (nl (i) /= nodi (j1, j) ) GO TO 31 GO TO 32 31 END DO GO TO 30 32 DO 33 j2 = 1, nj IF (nm (i) /= nodi (j2, j) ) GO TO 33 GO TO 34 33 END DO GO TO 30 34 DO 35 j3 = 1, nj IF (nn (i) /= nodi (j3, j) ) GO TO 35 IF (j1 < j2 .and. j2 < j3) GO TO 111 IF (j2 < j3 .and. j3 < j1) GO TO 111 IF (j1 < j2 .and. j3 < j1) GO TO 111 nout = nout + 1 GO TO 100 35 END DO GO TO 30 111 DO j4 = 1, nj nod = nodi (j4, j) xtemp (j4) = x (nod) ytemp (j4) = y (nod) END DO x1 = x (nl (i) ) y1 = y (nl (i) ) x2 = x (nm (i) ) y2 = y (nm (i) ) x3 = x (nn (i) ) y3 = y (nn (i) ) x12 = 0.5 * (x1 + x2) y12 = 0.5 * (y1 + y2) x23 = 0.5 * (x3 + x2) y23 = 0.5 * (y3 + y2) x31 = 0.5 * (x1 + x3) y31 = 0.5 * (y1 + y3) xp = x1 + (x23 - x1) / 9.d0 yp = y1 + (y23 - y1) / 9.d0 CALL inorout (xp, yp, nj, xtemp, ytemp, ichk, mtx, iangl) IF (ichk == 0) GO TO 112 nout = nout + 1 GO TO 100 112 xq = x2 + (x31 - x2) / 9.d0 yq = y2 + (y31 - y2) / 9.d0 CALL inorout (xq, yq, nj, xtemp, ytemp, ichk, mtx, iangl) IF (ichk == 0) GO TO 113 nout = nout + 1 GO TO 100 113 xr = x3 + (x12 - x3) / 9.d0 yr = y3 + (y12 - y3) / 9.d0 CALL inorout (xr, yr, nj, xtemp, ytemp, ichk, mtx, iangl) IF (ichk == 0) GO TO 90 nout = nout + 1 GO TO 100 30 END DO 90 nin = nin + 1 !b xxx check against max iflag (nin) = i 100 END DO IF (ntri /= (nin + nout) ) then WRITE (13, 200) nin, nout 200 FORMAT (//10x,'program stopped as the no. of rejected triangles ' & & /10x,i3,' plus the no. of included triangles ',i3/10x & & ,'is not equal to the total ',i3//) STOP END IF ntri = nin DO i = 1, ntri j = iflag (i) nl (i) = nl (j) nm (i) = nm (j) nn (i) = nn (j) END DO END SUBROUTINE final ! ********************************************************** ! ! ** subroutine conver ! ! ********************************************************** ! SUBROUTINE conver (npoin, x, y, ntri, nl, nm, nn, nlm, nmn, nnl, & nat, nnpe, natrib, mxfn, mtfn, mbtri, mbqua, iatt, matt, nmat, & lnods, lnod4, lma4, nbpl, npinb, nbno, npin6, nbno6) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) DIMENSION x (mxfn), y (mxfn), nl (mtfn), nm (mtfn), nn (mtfn), & nlm (mtfn), nmn (mtfn), nnl (mtfn), nat (mtfn), iatt (matt), & mbtri (mtfn), mbqua (mtfn), lnods (mtfn, 6), lnod4 (mtfn, 4), & lma4 (mtfn) !b M_SEG M_SEG M_SEG DIMENSION npinb (189), nbno (189, 489), npin6 (189), & nbno6 (189, 979) IF (nnpe /= 3) then CALL midnod (npoin, x, y, ntri, nl, nm, nn, nlm, nmn, & nnl, nnpe, mxfn, mtfn) DO i = 1, ntri lnods (i, 1) = nl (i) lnods (i, 2) = nlm (i) lnods (i, 3) = nm (i) lnods (i, 4) = nmn (i) lnods (i, 5) = nn (i) lnods (i, 6) = nnl (i) END DO END IF IF (nnpe /= 3) CALL bquad (mtfn, ntri, nbpl, npinb, nbno, & npin6, nbno6, mbtri, lnods) IF (nnpe /= 3 .and. nnpe /= 6) then CALL quadgen (mxfn, mtfn, ntri, npoin, nat, lnods, lnod4, & lma4, mbtri, mbqua, x, y) END IF nmat = 0 DO i = 1, natrib IF (nmat < iatt (i) ) nmat = iatt (i) END DO END SUBROUTINE conver ! ! ********************************************************************** ! ! ** subroutine collect ! ! ********************************************************************** ! SUBROUTINE collect (npoin, x, y, ntri, nl, nm, nn, iat, ntot, & xfin, yfin, ntr, ml, mm, mn, nat, noutb, ninnb, nodref, nodchk, & ncref, ibtri, mbtri, mx, mtri, mxfn, mtfn) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! DIMENSION x (mx), y (mx), nl (mtri), nm (mtri), nn (mtri), & xfin (mxfn), yfin (mxfn), ml (mtfn), mm (mtfn), mn (mtfn), & nat (mtfn), ncref (mx), nodref (mx), nodchk (mx), ibtri (mtri), & mbtri (mtfn) CHARACTER string * 4 !b print *,'collect ntot ', ntot !b string = 'mtfn' nb = noutb + ninnb DO i = 1, npoin IF (i <= nb) then j = nodref (i) IF (nodchk (j) == 0) then xfin (j) = x (i) yfin (j) = y (i) nodchk (j) = 1 END IF ncref (i) = j ELSE j = ntot + i - nb xfin (j) = x (i) yfin (j) = y (i) END IF END DO DO i = 1, ntri j = ntr + i CALL dcheck (j, mtfn, string) IF (nl (i) <= nb) then k = ncref (nl (i) ) ml (j) = k ELSE ml (j) = nl (i) + ntot - nb END IF IF (nm (i) <= nb) then k = ncref (nm (i) ) mm (j) = k ELSE mm (j) = nm (i) + ntot - nb END IF IF (nn (i) <= nb) then k = ncref (nn (i) ) mn (j) = k ELSE mn (j) = nn (i) + ntot - nb END IF nat (j) = iat mbtri (j) = ibtri (i) END DO ntot = ntot + npoin - nb ntr = ntr + ntri !b print *,'collect ntot, npoin, nb, ntri, ntr ',ntot, npoin, nb, ntri, ntr END SUBROUTINE collect ! ! ********************************************************************** ! ! ** subroutine divis ! ! ********************************************************************** ! SUBROUTINE divis (xl, yl, temp, div, xdiv, ydiv) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) IF (xl == 0.) then xdiv = 0. IF (yl > 0.) then ydiv = div ELSE ydiv = - div END IF ELSE xdiv = div / dsqrt (1. + temp**2) IF (xl < 0.) xdiv = - xdiv ydiv = xdiv * temp END IF END SUBROUTINE divis ! ! ********************************************************************** ! ! *** subroutine nodins (written by YAO ZHENG) ! ! ********************************************************************** ! SUBROUTINE nodins (n, x, y, tolrn, xd, yd, val, ndpt, xnu, ynu, & rnu, newp, siz, ord, xmax, xmin, ymax, ymin, mtx, mden, ityqd, & iadap, idpth, xquad, yquad, rquad, iflag, iornt, mxin) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine performs quadtree decomposition ! *** to generate points in an imagary square ! *** based on the given density distribution ! DIMENSION x (mtx), y (mtx), val (mden), xd (mden), yd (mden), & xnu (mxin), ynu (mxin), rnu (mxin) DIMENSION xquad (mxin), yquad (mxin), rquad (mxin), iflag (mxin), & iornt (mxin) ! ! mnew -- maximal number of quadtree cells mnew = 6999 ! iflag=-1 -- to be subdivided ! iflag= 0 -- not yet tested ! iflag= 1 -- finished subdivision ! ityqd= 1 -- square quadtree ! ityqd= 2 -- triangle quadtree ! iornt= 1 -- upwards (triangle cells) ! iornt=-1 -- downwards (triangle cells) root2 = 1.41421 root3 = 1.73205 ! IF (ityqd == 2) GO TO 100 ! square quadtree (ityqd = 1) jdpth = 0 newp = 0 xrang = xmax - xmin yrang = ymax - ymin range = dmax1 (xrang, yrang) xcnt = (xmax + xmin) / 2.d0 ycnt = (ymax + ymin) / 2.d0 ! ! write (13,*) 'xcnt=',xcnt ! write (13,*) 'ycnt=',ycnt ! xquad (1) = xcnt yquad (1) = ycnt rquad (1) = range itest = 1 xtest = xquad (1) ytest = yquad (1) rtest = rquad (1) IF (jdpth < idpth) then iflag (1) = - 1 GO TO 7 END IF CALL minsiz (xtest, ytest, rtest, xd, yd, val, ndpt, tolrn, siz, & ord, szmin, mden, iadap) ! sgrid=0.9306d0*siz*szmin sgrid = root2 * siz * szmin IF (rquad (1) > sgrid) then iflag (1) = - 1 ELSE iflag (1) = 1 END IF 7 CONTINUE nquad = 1 nfnsh = 0 nnew = 1 ! IF (iflag (1) == 1) then xnu (newp + 1) = xquad (itest) ynu (newp + 1) = yquad (itest) rnu (newp + 1) = rquad (itest) newp = newp + 1 !b xxx check against max GO TO 30 END IF ! nfnsh = nquad xquad (nquad+1) = xquad (1) - rquad (1) / 4.d0 yquad (nquad+1) = yquad (1) - rquad (1) / 4.d0 rquad (nquad+1) = rquad (1) / 2.d0 iflag (nquad+1) = 0 xquad (nquad+2) = xquad (1) + rquad (1) / 4.d0 yquad (nquad+2) = yquad (1) - rquad (1) / 4.d0 rquad (nquad+2) = rquad (1) / 2.d0 iflag (nquad+2) = 0 xquad (nquad+3) = xquad (1) + rquad (1) / 4.d0 yquad (nquad+3) = yquad (1) + rquad (1) / 4.d0 rquad (nquad+3) = rquad (1) / 2.d0 iflag (nquad+3) = 0 xquad (nquad+4) = xquad (1) - rquad (1) / 4.d0 yquad (nquad+4) = yquad (1) + rquad (1) / 4.d0 rquad (nquad+4) = rquad (1) / 2.d0 iflag (nquad+4) = 0 nquad = nquad+4 ! ! write (13,*)'sgrid=',sgrid ! write (13,*)'xquad1=',xquad(1) ! write (13,*)'yquad1=',yquad(1) ! write (13,*)'rquad1=',rquad(1) ! write (13,*)'iflag1=',iflag(1) ! 10 IF (nnew == 1) GO TO 15 DO i = nfnsh + 1, nnew IF (iflag (i) == - 1) then xquad (nquad+1) = xquad (i) - rquad (i) / 4.d0 yquad (nquad+1) = yquad (i) - rquad (i) / 4.d0 rquad (nquad+1) = rquad (i) / 2.d0 iflag (nquad+1) = 0 xquad (nquad+2) = xquad (i) + rquad (i) / 4.d0 yquad (nquad+2) = yquad (i) - rquad (i) / 4.d0 rquad (nquad+2) = rquad (i) / 2.d0 iflag (nquad+2) = 0 xquad (nquad+3) = xquad (i) + rquad (i) / 4.d0 yquad (nquad+3) = yquad (i) + rquad (i) / 4.d0 rquad (nquad+3) = rquad (i) / 2.d0 iflag (nquad+3) = 0 xquad (nquad+4) = xquad (i) - rquad (i) / 4.d0 yquad (nquad+4) = yquad (i) + rquad (i) / 4.d0 rquad (nquad+4) = rquad (i) / 2.d0 iflag (nquad+4) = 0 nquad = nquad+4 END IF END DO 15 nfnsh = nnew nnew = nquad jdpth = jdpth + 1 ! DO itest = nfnsh + 1, nnew xtest = xquad (itest) ytest = yquad (itest) rtest = rquad (itest) ! CALL boxcheck (xtest, ytest, rtest, iout, xmax, xmin, ymax, ymin) IF (iout == 1) then iflag (itest) = 1 GO TO 17 END IF ! IF (jdpth < idpth) then iflag (itest) = - 1 GO TO 17 END IF CALL minsiz (xtest, ytest, rtest, xd, yd, val, ndpt, tolrn, siz, & ord, szmin, mden, iadap) ! sgrid=0.9306d0*siz*szmin sgrid = root2 * siz * szmin IF (rquad (itest) > sgrid) then iflag (itest) = - 1 ELSE iflag (itest) = 1 xnu (newp + 1) = xquad (itest) ynu (newp + 1) = yquad (itest) rnu (newp + 1) = rquad (itest) newp = newp + 1 !b xxx check against max END IF 17 CONTINUE END DO ! nmove = 0 DO itest = nfnsh + 1, nnew IF (iflag (itest) == - 1) then GO TO 20 ELSE nmove = nmove+1 END IF END DO 20 nfnsh = nfnsh + nmove ! ! write (13,*)'nfnsh=',nfnsh ! ! if (nnew > mnew) write (13,*)'***' IF (nnew > mnew) GO TO 30 IF (nfnsh == nnew) GO TO 30 ! if (nfnsh /= nnew) go to 10 ! if (nmove /= 0) go to 10 GO TO 10 30 CONTINUE ! ! write (13,*)'iflag1=',iflag(1) ! write (13,*)'iflag2=',iflag(2) ! write (13,*)'iflag3=',iflag(3) ! write (13,*)'iflag4=',iflag(4) ! write (13,*)'iflag5=',iflag(5) ! write (13,*)'iflag6=',iflag(6) ! write (13,*)'iflag7=',iflag(7) ! write (13,*)'iflag8=',iflag(8) ! write (13,*)'iflag9=',iflag(9) ! write (13,*)'iflag10=',iflag(10) ! write (13,*)'nmove=',nmove ! write (13,*)'nquad=',nquad ! write (13,*)'nnew=',nnew ! write (13,*)'newp=',newp ! write (13,*)'xnu(1)=',xnu(1) ! write (13,*)'ynu(1)=',ynu(1) ! GO TO 200 ! ! triangle quadtree (ityqd = 2) 100 CONTINUE jdpth = 0 newp = 0 xrang = xmax - xmin yrang = ymax - ymin range = root2 * root3 * dmax1 (xrang, yrang) xcnt = (xmax + xmin) / 2.d0 ycnt = (ymax + ymin) / 2.d0 ! ! write (13,*) 'xcnt=',xcnt ! write (13,*) 'ycnt=',ycnt ! xquad (1) = xcnt yquad (1) = ycnt rquad (1) = range iornt (1) = 1 itest = 1 xtest = xquad (1) ytest = yquad (1) rtest = rquad (1) IF (jdpth < idpth) then iflag (1) = - 1 GO TO 107 END IF CALL minsiz (xtest, ytest, rtest, xd, yd, val, ndpt, tolrn, siz, & ord, szmin, mden, iadap) ! sgrid=0.9306d0*siz*szmin sgrid = 2.d0 * root2 * siz * szmin IF (rquad (1) > sgrid) then iflag (1) = - 1 ELSE iflag (1) = 1 END IF 107 CONTINUE nquad = 1 nfnsh = 0 nnew = 1 ! IF (iflag (1) == 1) then xnu (newp + 1) = xquad (itest) ynu (newp + 1) = yquad (itest) + rquad (itest) / 2.d0 / root3 rnu (newp + 1) = rquad (itest) xnu (newp + 2) = xquad (itest) - rquad (itest) / 4.d0 ynu (newp + 2) = yquad (itest) - rquad (itest) / 4.d0 / root3 rnu (newp + 2) = rquad (itest) xnu (newp + 3) = xquad (itest) + rquad (itest) / 4.d0 ynu (newp + 3) = yquad (itest) - rquad (itest) / 4.d0 / root3 rnu (newp + 3) = rquad (itest) newp = newp + 3 GO TO 130 END IF ! nfnsh = nquad xquad (nquad+1) = xquad (1) yquad (nquad+1) = yquad (1) + rquad (1) / 2.d0 / root3 rquad (nquad+1) = rquad (1) / 2.d0 iornt (nquad+1) = iornt (1) iflag (nquad+1) = 0 xquad (nquad+2) = xquad (1) - rquad (1) / 4.d0 yquad (nquad+2) = yquad (1) - rquad (1) / 4.d0 / root3 rquad (nquad+2) = rquad (1) / 2.d0 iornt (nquad+2) = iornt (1) iflag (nquad+2) = 0 xquad (nquad+3) = xquad (1) yquad (nquad+3) = yquad (1) rquad (nquad+3) = rquad (1) / 2.d0 iornt (nquad+3) = - iornt (1) iflag (nquad+3) = 0 xquad (nquad+4) = xquad (1) + rquad (1) / 4.d0 yquad (nquad+4) = yquad (1) - rquad (1) / 4.d0 / root3 rquad (nquad+4) = rquad (1) / 2.d0 iornt (nquad+4) = iornt (1) iflag (nquad+4) = 0 nquad = nquad+4 ! ! write (13,*)'sgrid=',sgrid ! write (13,*)'xquad1=',xquad(1) ! write (13,*)'yquad1=',yquad(1) ! write (13,*)'rquad1=',rquad(1) ! write (13,*)'iflag1=',iflag(1) ! 110 IF (nnew == 1) GO TO 115 DO i = nfnsh + 1, nnew IF (iflag (i) == - 1) then xquad (nquad+1) = xquad (i) yquad (nquad+1) = yquad (i) + iornt (i) * rquad (i) / 2.d0 / & root3 rquad (nquad+1) = rquad (i) / 2.d0 iornt (nquad+1) = iornt (i) iflag (nquad+1) = 0 xquad (nquad+2) = xquad (i) - iornt (i) * rquad (i) / 4.d0 yquad (nquad+2) = yquad (i) - iornt (i) * rquad (i) / 4.d0 / & root3 rquad (nquad+2) = rquad (i) / 2.d0 iornt (nquad+2) = iornt (i) iflag (nquad+2) = 0 xquad (nquad+3) = xquad (i) yquad (nquad+3) = yquad (i) rquad (nquad+3) = rquad (i) / 2.d0 iornt (nquad+3) = - iornt (i) iflag (nquad+3) = 0 xquad (nquad+4) = xquad (i) + iornt (i) * rquad (i) / 4.d0 yquad (nquad+4) = yquad (i) - iornt (i) * rquad (i) / 4.d0 / & root3 rquad (nquad+4) = rquad (i) / 2.d0 iornt (nquad+4) = iornt (i) iflag (nquad+4) = 0 nquad = nquad+4 END IF END DO 115 nfnsh = nnew nnew = nquad jdpth = jdpth + 1 ! DO itest = nfnsh + 1, nnew xtest = xquad (itest) ytest = yquad (itest) rtest = rquad (itest) ! CALL boxcheck (xtest, ytest, rtest, iout, xmax, xmin, ymax, ymin) IF (iout == 1) then iflag (itest) = 1 GO TO 117 END IF ! IF (jdpth < idpth) then iflag (itest) = - 1 GO TO 117 END IF CALL minsiz (xtest, ytest, rtest, xd, yd, val, ndpt, tolrn, siz, & ord, szmin, mden, iadap) ! sgrid=0.9306d0*siz*szmin sgrid = 2.d0 * root2 * siz * szmin IF (rquad (itest) > sgrid) then iflag (itest) = - 1 ELSE iflag (itest) = 1 IF (iornt (itest) == 1) then xnu (newp + 1) = xquad (itest) ynu (newp + 1) = yquad (itest) + rquad (itest) / 2.d0 / root3 rnu (newp + 1) = rquad (itest) xnu (newp + 2) = xquad (itest) - rquad (itest) / 4.d0 ynu (newp + 2) = yquad (itest) - rquad (itest) / 4.d0 / root3 rnu (newp + 2) = rquad (itest) xnu (newp + 3) = xquad (itest) + rquad (itest) / 4.d0 ynu (newp + 3) = yquad (itest) - rquad (itest) / 4.d0 / root3 rnu (newp + 3) = rquad (itest) newp = newp + 3 !b xxx check against max ELSE xnu (newp + 1) = xquad (itest) ynu (newp + 1) = yquad (itest) rnu (newp + 1) = rquad (itest) newp = newp + 1 !b xxx check against max END IF END IF 117 CONTINUE END DO ! nmove = 0 DO itest = nfnsh + 1, nnew IF (iflag (itest) == - 1) then GO TO 120 ELSE nmove = nmove+1 END IF END DO 120 nfnsh = nfnsh + nmove ! ! write (13,*)'nfnsh=',nfnsh ! ! if (nnew > mnew) write (13,*)'***' IF (nnew > mnew) GO TO 130 IF (nfnsh == nnew) GO TO 130 ! if (nfnsh /= nnew) go to 110 ! if (nmove /= 0) go to 110 GO TO 110 130 CONTINUE ! ! write (13,*)'iflag1=',iflag(1) ! write (13,*)'iflag2=',iflag(2) ! write (13,*)'iflag3=',iflag(3) ! write (13,*)'iflag4=',iflag(4) ! write (13,*)'iflag5=',iflag(5) ! write (13,*)'iflag6=',iflag(6) ! write (13,*)'iflag7=',iflag(7) ! write (13,*)'iflag8=',iflag(8) ! write (13,*)'iflag9=',iflag(9) ! write (13,*)'iflag10=',iflag(10) ! write (13,*)'nmove=',nmove ! write (13,*)'nquad=',nquad ! write (13,*)'nnew=',nnew ! write (13,*)'newp=',newp ! write (13,*)'xnu(1)=',xnu(1) ! write (13,*)'ynu(1)=',ynu(1) ! 200 CONTINUE irm = 0 DO i = 1, newp xtest = xnu (i) ytest = ynu (i) rtest = 0.d0 CALL boxcheck (xtest, ytest, rtest, iout, xmax, xmin, & ymax, ymin) IF (iout == 1) then irm = irm + 1 ELSE xnu (i - irm) = xnu (i) ynu (i - irm) = ynu (i) rnu (i - irm) = rnu (i) END IF END DO newp = newp - irm END SUBROUTINE nodins ! ! ********************************************************************** ! ! ** subroutine boxcheck ! ! ********************************************************************** ! SUBROUTINE boxcheck (x, y, r, iout, xmax, xmin, ymax, ymin) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) iout = 1 IF ( (x + r / 2.d0) < xmin) RETURN !b GO TO 10 IF ( (x - r / 2.d0) > xmax) RETURN !b GO TO 10 IF ( (y + r / 2.d0) < ymin) RETURN !b GO TO 10 IF ( (y - r / 2.d0) > ymax) RETURN !b GO TO 10 iout = 0 !b 10 CONTINUE END SUBROUTINE boxcheck ! ! ********************************************************************** ! ! ** subroutine newpt ! ! ********************************************************************** ! SUBROUTINE newpt (xin, yin, x, y, cx, cy, s, sz1, sz2) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) fctr = 1 IF (sz1 < 0.) fctr = - fctr IF (sz1 < 0.) sz1 = sz1 * fctr cs = (sz1**2 + s**2 - sz2**2) / (2. * sz1 * s) const = 1.d0 - cs**2 IF (const > 0.d0) then sn = dsqrt (const) xn = sz1 * cs yn = sz1 * sn * fctr ELSE xn = s / 2.d0 yn = 0.d0 WRITE (13, 501) 501 FORMAT(//10x,'*************** warning ***************',& & /10x,'the dist b/w adjacent points is greater' & & /10x,'than the sum of the elem sizes of these'//) END IF x = (cx * xn - cy * yn) + xin y = (cy * xn + cx * yn) + yin END SUBROUTINE newpt ! ! ********************************************************************** ! ! ** subroutine area ! ! ********************************************************************** ! SUBROUTINE area (x, y, n, result) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! DIMENSION x (1000), y (1000) temp1 = y (1) * (x (2) - x (n) ) + y (n) * (x (1) - x (n - 1) ) temp = 0.d0 DO i = 2, n - 1 temp = temp + y (i) * (x (i + 1) - x (i - 1) ) END DO !b result=(temp1+temp)/-2. result = - (temp1 + temp) / 2. END SUBROUTINE area