! 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 ! ! ********************************************************************** ! ! ** subroutine smooth ! ! ********************************************************************** ! SUBROUTINE smooth (npoin, x, y, ntri, nl, nm, nn, noutb, ninnb, & mx, mtri) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) ! ! *** this subroutine performs the required number of smoothing ! *** iterations on the mesh generated on every subdomain. ! DIMENSION x (mx), y (mx), nl (mtri), nm (mtri), nn (mtri), & xcoord (4999), ycoord (4999), nlt (4999) ! nfp = noutb + ninnb DO i = nfp + 1, npoin nlt (i) = 0 xcoord (i) = 0.d0 ycoord (i) = 0.d0 DO j = 1, ntri i1 = nl (j) i2 = nm (j) i3 = nn (j) IF (i == i1) then nlt (i) = nlt (i) + 1 xcoord (i) = xcoord (i) + x (i2) + x (i3) ycoord (i) = ycoord (i) + y (i2) + y (i3) END IF IF (i == i2) then nlt (i) = nlt (i) + 1 xcoord (i) = xcoord (i) + x (i1) + x (i3) ycoord (i) = ycoord (i) + y (i1) + y (i3) END IF IF (i == i3) then nlt (i) = nlt (i) + 1 xcoord (i) = xcoord (i) + x (i1) + x (i2) ycoord (i) = ycoord (i) + y (i1) + y (i2) END IF END DO END DO DO i = nfp + 1, npoin nd = nlt (i) * 2 x (i) = xcoord (i) / nd y (i) = ycoord (i) / nd END DO END SUBROUTINE smooth ! ********************************************************************** ! ! *** subroutine midnod ! ! ********************************************************************** ! SUBROUTINE midnod (npoin, x, y, ntri, nl, nm, nn, nlm, nmn, nnl, & nnpe, mxfn, mtfn) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) INTEGER, PARAMETER :: limit_1 = 15000, limit_2 = 25000 ! ??? ! = ntri = 2*mxfn ! ! *** this subroutine inserts midside nodes if six noded triangular ! *** elements are required. ! DIMENSION x (mxfn), y (mxfn), nl (mtfn), nm (mtfn), nn (mtfn), & nlm (mtfn), nmn (mtfn), nnl (mtfn), np (limit_1), nq (limit_1), & nr (limit_1), num (limit_2), l1 (limit_2), l2 (limit_2) ! nseg = 0 if ( 2*mxfn > limit_2 ) then print *,'midnod: expect limit_2 >= ', 2*mxfn if ( ntri > limit_1 ) then print *, 'midnod: increase limit_1 to at least ', ntri stop 'ERROR in midnod: increase limit_1' end if DO i = 1, ntri i1 = nl (i) i2 = nm (i) i3 = nn (i) np (i) = 0 nq (i) = 0 nr (i) = 0 IF (i == 1) then nseg = nseg + 1 !b xxx check against max if ( nseg > limit_2 ) then print *, 'midnod: increase limit_2 to at least ', nseg stop 'ERROR in midnod: increase limit_2' end if l1 (nseg) = i1 l2 (nseg) = i2 num (nseg) = 1 np (1) = nseg IF (nnpe /= 3) x (npoin + nseg) = (x (i1) + x (i2) ) / 2.d0 IF (nnpe /= 3) y (npoin + nseg) = (y (i1) + y (i2) ) / 2.d0 nseg = nseg + 1 !b xxx check against max if ( nseg > limit_2 ) then print *, 'midnod: increase limit_2 to at least ', nseg stop 'ERROR in midnod: increase limit_2' end if l1 (nseg) = i2 l2 (nseg) = i3 num (nseg) = 1 nq (1) = nseg IF (nnpe /= 3) x (npoin + nseg) = (x (i2) + x (i3) ) / 2.d0 IF (nnpe /= 3) y (npoin + nseg) = (y (i2) + y (i3) ) / 2.d0 nseg = nseg + 1 !b xxx check against max if ( nseg > limit_2 ) then print *, 'midnod: increase limit_2 to at least ', nseg stop 'ERROR in midnod: increase limit_2' end if l1 (nseg) = i3 l2 (nseg) = i1 nr (1) = nseg num (nseg) = 1 IF (nnpe /= 3) x (npoin + nseg) = (x (i3) + x (i1) ) / 2.d0 IF (nnpe /= 3) y (npoin + nseg) = (y (i3) + y (i1) ) / 2.d0 ELSE DO j = 1, nseg IF (num (j) == 2) GO TO 10 IF (i2 == l1 (j) .and. i1 == l2 (j) ) then num (j) = 2 np (i) = j ELSEIF (i3 == l1 (j) .and. i2 == l2 (j) ) then num (j) = 2 nq (i) = j ELSEIF (i1 == l1 (j) .and. i3 == l2 (j) ) then num (j) = 2 nr (i) = j END IF 10 CONTINUE END DO IF (np (i) == 0) then nseg = nseg + 1 !b xxx check against max if ( nseg > limit_2 ) then print *, 'midnod: increase limit_2 to at least ', nseg stop 'ERROR in midnod: increase limit_2' end if l1 (nseg) = i1 l2 (nseg) = i2 num (nseg) = 1 np (i) = nseg IF (nnpe /= 3) x (npoin + nseg) = (x (i1) + x (i2) ) / 2.d0 IF (nnpe /= 3) y (npoin + nseg) = (y (i1) + y (i2) ) / 2.d0 END IF IF (nq (i) == 0) then nseg = nseg + 1 !b xxx check against max if ( nseg > limit_2 ) then print *, 'midnod: increase limit_2 to at least ', nseg stop 'ERROR in midnod: increase limit_2' end if l1 (nseg) = i2 l2 (nseg) = i3 num (nseg) = 1 nq (i) = nseg IF (nnpe /= 3) x (npoin + nseg) = (x (i2) + x (i3) ) / 2.d0 IF (nnpe /= 3) y (npoin + nseg) = (y (i2) + y (i3) ) / 2.d0 END IF IF (nr (i) == 0) then nseg = nseg + 1 !b xxx check against max if ( nseg > limit_2 ) then print *, 'midnod: increase limit_2 to at least ', nseg stop 'ERROR in midnod: increase limit_2' end if l1 (nseg) = i3 l2 (nseg) = i1 num (nseg) = 1 nr (i) = nseg IF (nnpe /= 3) x (npoin + nseg) = (x (i3) + x (i1) ) / 2.d0 IF (nnpe /= 3) y (npoin + nseg) = (y (i3) + y (i1) ) / 2.d0 END IF END IF IF (nnpe /= 3) then nlm (i) = npoin + np (i) nmn (i) = npoin + nq (i) nnl (i) = npoin + nr (i) END IF END DO IF (nnpe /= 3) npoin = npoin + nseg END SUBROUTINE midnod ! ! ********************************************************************** ! ! ** subroutine bchek ! ! ********************************************************************** ! SUBROUTINE bchek (mtx, mx, mtri, mhol, mxin, npoin, ntri, nl, nm, & nn, noutb, inbnd, ninb, nodi, iflag, iflat, ibtri, & nbtri, nodtri, xb, yb, x, y) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) INTEGER, PARAMETER :: max_new = 99, max_k = 10 ! ! *** this subroutine rejects triangles straddling the boundary and ! *** enforces boundary integrity by replacing the rejected triangles ! *** appropriate ones. ! DIMENSION nl (mtri), nm (mtri), nn (mtri), iflag (mxin), & ibtri (mtri), ninb (mhol), nodi (mtx, mhol), nbtri (mtx), & nodtri (mtx, 20), x (mx), y (mx), xb (mtx), yb (mtx), & kpts (max_k), newnl (max_new), newnm (max_new), & newnn (max_new), newbt (max_new), iflat (mtri) ! ! *** find any boundary segment not part of the triangulation ! DO i = 1, ntri iflat (i) = 1 ibtri (i) = 0 END DO irej = 0 newtr = 0 ! ! *1* check the outer boundaries. ! ! *** make a list of triangles associated with each boundary node ! DO i = 1, noutb j = i + 1 IF (i == noutb) j = 1 iflag (i) = 0 nbtri (i) = 0 DO k = 1, ntri k1 = nl (k) k2 = nm (k) k3 = nn (k) IF (i == k1 .or. (i == k2 .or. i == k3) ) then nbtri (i) = nbtri (i) + 1 !b xxx check against max nodtri (i, nbtri (i) ) = k END IF IF (k1 == i .and. k2 == j) then ibtri (k) = 1 iflag (i) = k END IF IF (k2 == i .and. k3 == j) then IF (ibtri (k) == 0) ibtri (k) = 2 IF (ibtri (k) == 1) ibtri (k) = 4 iflag (i) = k END IF IF (k3 == i .and. k1 == j) then IF (ibtri (k) == 0) ibtri (k) = 3 IF (ibtri (k) == 2) ibtri (k) = 5 IF (ibtri (k) == 1) ibtri (k) = 6 iflag (i) = k END IF END DO END DO ! ! *** find the triangles intersecting missing segments ! iswit = 0 DO i = 1, noutb IF (iswit == 1) then iswit = 0 GO TO 5 END IF IF (iflag (i) == 0) then ipts = 0 j = i + 1 IF (j > noutb) j = 1 IF (nbtri (j) == 0) then j = i + 2 IF (j > noutb) j = 1 iswit = 1 END IF ax = xb (i) ay = yb (i) bx = xb (j) by = yb (j) ni = nbtri (i) DO k = 1, ni ktri = nodtri (i, k) IF (i == nl (ktri) ) then i1 = nm (ktri) i2 = nn (ktri) ELSEIF (i == nm (ktri) ) then i1 = nn (ktri) i2 = nl (ktri) ELSEIF (i == nn (ktri) ) then i1 = nl (ktri) i2 = nm (ktri) END IF cx = x (i1) cy = y (i1) dx = x (i2) dy = y (i2) CALL intsct (ax, ay, bx, by, cx, cy, dx, dy, iout) IF (iout == 1) then IF (iflat (ktri) == 1) then iflat (ktri) = 0 irej = irej + 1 END IF ipts = ipts + 1 ! check against max if ( ipts > max_k ) then print *, 'bchek: increase max_k to at least ', ipts stop 'ERROR in bchek: increase max_k' end if kpts (ipts) = i2 2 CALL adjcnt (mtri, ntri, nl, nm, nn, i1, i2, jtri, jn) IF (jtri /= 0) then IF (iflat (jtri) == 1) then iflat (jtri) = 0 irej = irej + 1 END IF IF (jn == j) GO TO 1 cx = x (i1) cy = y (i1) dx = x (jn) dy = y (jn) CALL intsct (ax, ay, bx, by, cx, cy, dx, dy, iout) IF (iout == 1) then ipts = ipts + 1 ! check against max if ( ipts > max_k ) then print *, 'bchek: increase max_k to at least ', ipts stop 'ERROR in bchek: increase max_k' end if kpts (ipts) = jn i2 = jn ELSE i1 = jn END IF GO TO 2 ELSE PRINT * , 'MESHING FAILURE HAS OCCURRED !' PRINT * , 'TRY AGAIN BY MODIFYING SOME PARAMETER.' PRINT * , 'IN CASE OF REPEATED FAILURES RE-ENTER THE' PRINT * , 'GEOMETRY DIVIDING IT INTO MORE SUBDOMAINS.' WRITE (13, * ) 'error in bchek' STOP END IF END IF END DO ! ! *** find two triangles one connecting i and the other connecting j ! *** with a common node. ! nj = nbtri (j) DO ii = 1, ni itri = nodtri (i, ii) IF (nl (itri) == i) then i1 = nm (itri) ELSEIF (nm (itri) == i) then i1 = nn (itri) ELSEIF (nn (itri) == i) then i1 = nl (itri) END IF DO ij = 1, nj jtri = nodtri (j, ij) IF (nl (jtri) == j) then j2 = nn (jtri) ELSEIF (nm (jtri) == j) then j2 = nl (jtri) ELSEIF (nn (jtri) == j) then j2 = nm (jtri) END IF IF (i1 == j2) then ipts = ipts + 1 ! check against max if ( ipts > max_k ) then print *, 'bchek: increase max_k to at least ', ipts stop 'ERROR in bchek: increase max_k' end if kpts (ipts) = i1 GO TO 1 END IF END DO END DO ! ! *** If no triangles are found connecting i and j assume that two ! *** triangles are needed to fill in the gap. The points to be used ! *** for this purpose are the first point after i and the first point ! *** before j (i1 and j2). ! IF (ipts == 0) then DO ii = 1, ni itri = nodtri (i, ii) IF (nl (itri) == i) then i1 = nm (itri) ELSEIF (nm (itri) == i) then i1 = nn (itri) ELSEIF (nn (itri) == i) then i1 = nl (itri) END IF ntri1 = 0 DO k = 1, ntri k1 = nl (k) k2 = nm (k) k3 = nn (k) IF (i1 == k1 .or. (i1 == k2 .or. i1 == k3) ) then DO ij = 1, nj jtri = nodtri (j, ij) IF (nl (jtri) == j) then j2 = nn (jtri) ELSEIF (nm (jtri) == j) then j2 = nl (jtri) ELSEIF (nn (jtri) == j) then j2 = nm (jtri) END IF IF (j2 == k1 .or. (j2 == k2 .or. j2 == k3) ) then ipts = ipts + 1 ! check against max if ( ipts > max_k ) then print *, 'bchek: increase max_k to >= ', ipts stop 'ERROR in bchek: increase max_k' end if kpts (ipts) = i1 ipts = ipts + 1 !b xxx check against max if ( ipts > max_k ) then print *, 'bchek: increase max_k to >= ', ipts stop 'ERROR in bchek: increase max_k' end if kpts (ipts) = j2 GO TO 1 END IF END DO END IF END DO END DO END IF IF (ipts == 0) then PRINT * , 'MESHING FAILURE HAS OCCURRED !' PRINT * , 'TRY AGAIN BY MODIFYING SOME PARAMETER.' PRINT * , 'IN CASE OF REPEATED FAILURES RE-ENTER THE' PRINT * , 'GEOMETRY DIVIDING IT INTO MORE SUBDOMAINS.' WRITE (13, * ) 'error in bchek' STOP END IF 1 newtr = newtr + 1 ! check against max if ( newtr > max_new ) then print *, 'bchek: increase max_new to at least ', newtr stop 'ERROR in bchek: increase max_new' end if newnl (newtr) = i newnm (newtr) = j newnn (newtr) = kpts (1) IF (iswit /= 1) newbt (newtr) = 1 IF (ipts >= 2) then DO ip = 2, ipts newtr = newtr + 1 ! check against max if ( newtr > max_new ) then print *, 'bchek: increase max_new to at least ', newtr stop 'ERROR in bchek: increase max_new' end if jp = ip - 1 newnl (newtr) = kpts (jp) newnm (newtr) = j newnn (newtr) = kpts (ip) END DO END IF END IF IF (iswit == 1) then newtr = newtr + 1 ! check against max if ( newtr > max_new ) then print *, 'bchek: increase max_new to at least ', newtr stop 'ERROR in bchek: increase max_new' end if newnl (newtr) = i newnm (newtr) = i + 1 newnn (newtr) = j newbt (newtr) = 4 END IF 5 CONTINUE END DO ! ! *2* check the inner boundaries. ! ! *** make a list of triangles associated with each boundary node ! DO inb = 1, inbnd nb = ninb (inb) DO in = 1, nb i = nodi (in, inb) jin = in + 1 !b xxx check against max IF (in == nb) jin = 1 j = nodi (jin, inb) iflag (in) = 0 nbtri (in) = 0 DO k = 1, ntri k1 = nl (k) k2 = nm (k) k3 = nn (k) IF (i == k1 .or. (i == k2 .or. i == k3) ) then nbtri (in) = nbtri (in) + 1 !b xxx check against max nodtri (in, nbtri (in) ) = k END IF IF (k1 == i .and. k2 == j) then ibtri (k) = 1 iflag (in) = k END IF IF (k2 == i .and. k3 == j) then IF (ibtri (k) == 0) ibtri (k) = 2 IF (ibtri (k) == 1) ibtri (k) = 4 iflag (in) = k END IF IF (k3 == i .and. k1 == j) then IF (ibtri (k) == 0) ibtri (k) = 3 IF (ibtri (k) == 2) ibtri (k) = 5 IF (ibtri (k) == 1) ibtri (k) = 6 iflag (in) = k END IF END DO END DO ! ! *** find the triangles intersecting missing segments ! DO in = 1, nb IF (iswit == 1) then iswit = 0 GO TO 6 END IF IF (iflag (in) == 0) then i = nodi (in, inb) jin = in + 1 IF (jin > nb) jin = 1 IF (nbtri (j) == 0) then jin = in + 2 IF (jin > nb) jin = 1 iswit = 1 END IF j = nodi (jin, inb) ni = nbtri (in) ipts = 0 ax = x (i) ay = y (i) bx = x (j) by = y (j) DO k = 1, ni ktri = nodtri (in, k) IF (i == nl (ktri) ) then i1 = nm (ktri) i2 = nn (ktri) ELSEIF (i == nm (ktri) ) then i1 = nn (ktri) i2 = nl (ktri) ELSEIF (i == nn (ktri) ) then i1 = nl (ktri) i2 = nm (ktri) END IF cx = x (i1) cy = y (i1) dx = x (i2) dy = y (i2) CALL intsct (ax, ay, bx, by, cx, cy, dx, dy, iout) IF (iout == 1) then IF (iflat (ktri) == 1) then iflat (ktri) = 0 irej = irej + 1 END IF ipts = ipts + 1 !b xxx check against max kpts (ipts) = i2 4 CALL adjcnt (mtri, ntri, nl, nm, nn, i1, i2, jtri, jn) IF (jtri /= 0) then IF (iflat (jtri) == 1) then iflat (jtri) = 0 irej = irej + 1 END IF IF (jn == j) GO TO 3 cx = x (i1) cy = y (i1) dx = x (jn) dy = y (jn) CALL intsct (ax, ay, bx, by, cx, cy, dx, dy, iout) IF (iout == 1) then ipts = ipts + 1 !b xxx check against max kpts (ipts) = jn i2 = jn ELSE i1 = jn END IF GO TO 4 ELSE PRINT * , 'MESHING FAILURE HAS OCCURRED !' PRINT * , 'TRY AGAIN BY MODIFYING SOME PARAMETER.' PRINT * , 'IN CASE OF REPEATED FAILURES RE-ENTER THE' PRINT * , 'GEOMETRY DIVIDING IT INTO MORE SUBDOMAINS.' WRITE (13, * ) 'error in bchek' STOP END IF END IF END DO ! ! *** find two triangles one connecting i and the other connecting j ! *** with a common node. ! nj = nbtri (jin) DO ii = 1, ni itri = nodtri (in, ii) IF (nl (itri) == i) then i1 = nm (itri) ELSEIF (nm (itri) == i) then i1 = nn (itri) ELSEIF (nn (itri) == i) then i1 = nl (itri) END IF DO ij = 1, nj jtri = nodtri (jin, ij) IF (nl (jtri) == j) then j2 = nn (jtri) ELSEIF (nm (jtri) == j) then j2 = nl (jtri) ELSEIF (nn (jtri) == j) then j2 = nm (jtri) END IF IF (i1 == j2) then ipts = ipts + 1 !b xxx check against max kpts (ipts) = i1 GO TO 3 END IF END DO END DO 3 newtr = newtr + 1 ! check against max if ( newtr > max_new ) then print *, 'bchek: increase max_new to at least ', newtr stop 'ERROR in bchek: increase max_new' end if newnl (newtr) = i newnm (newtr) = j newnn (newtr) = kpts (1) IF (iswit /= 1) newbt (newtr) = 1 IF (ipts >= 2) then DO ip = 2, ipts newtr = newtr + 1 ! check against max if ( newtr > max_new ) then print *, 'bchek: increase max_new to at least ', newtr stop 'ERROR in bchek: increase max_new' end if jp = ip - 1 newnl (newtr) = kpts (jp) newnm (newtr) = j newnn (newtr) = kpts (ip) END DO END IF END IF IF (iswit == 1) then newtr = newtr + 1 ! check against max if ( newtr > max_new ) then print *, 'bchek: increase max_new to at least ', newtr stop 'ERROR in bchek: increase max_new' end if newnl (newtr) = i newnm (newtr) = nodi (in + 1, inb) newnn (newtr) = j newbt (newtr) = 4 END IF 6 CONTINUE END DO END DO ! ! *** reject flagged triangles and incorporate new ones ! itri = 0 DO i = 1, ntri IF (iflat (i) == 1) then itri = itri + 1 !b xxx check against max nl (itri) = nl (i) nm (itri) = nm (i) nn (itri) = nn (i) ibtri (itri) = ibtri (i) END IF END DO DO i = 1, newtr itri = itri + 1 !b xxx check against max nl (itri) = newnl (i) nm (itri) = newnm (i) nn (itri) = newnn (i) ibtri (itri) = newbt (i) END DO ntri = itri END SUBROUTINE bchek ! ! ********************************************************************** ! ! ** subroutine intsct ! ! ********************************************************************** ! SUBROUTINE intsct (ax, ay, bx, by, cx, cy, dx, dy, iout) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) iout = 0 det = ( (bx - ax) * (dy - cy) ) - ( (by - ay) * (dx - cx) ) IF (det == 0) return s = ( (cy - dy) * (bx - dx) + (dx - cx) * (by - dy) ) / det t = ( (ay - by) * (bx - dx) + (bx - ax) * (by - dy) ) / det IF (s < - 1.d0 .or. 0.d0 < s) return IF (t < - 1.d0 .or. 0.d0 < t) return IF (s > - 1.d0 .and. 0.d0 > s) iout = 1 IF (t > - 1.d0 .and. 0.d0 > t) iout = 1 END SUBROUTINE intsct ! ! ********************************************************************** ! ! ** subroutine adjcnt ! ! ********************************************************************** ! SUBROUTINE adjcnt (mtri, ntri, nl, nm, nn, i1, i2, jtri, jn) DIMENSION nl (mtri), nm (mtri), nn (mtri) jtri = 0 DO j = 1, ntri j1 = nl (j) j2 = nm (j) j3 = nn (j) IF (i1 == j2 .and. i2 == j1) then jtri = j jn = j3 GO TO 1 ELSEIF (i1 == j3 .and. i2 == j2) then jtri = j jn = j1 GO TO 1 ELSEIF (i1 == j1 .and. i2 == j3) then jtri = j jn = j2 GO TO 1 END IF END DO 1 CONTINUE END SUBROUTINE adjcnt ! ! ********************************************************************** ! ! ** subroutine densit ! ! ********************************************************************** ! SUBROUTINE densit (xi, yi, xl, yl, val, ndpt, dens, tolrn, siz, & ord, mden, iadap) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) DIMENSION xl (mden), yl (mden), val (mden) ! param = 5.d0 IF (iadap /= 2) GO TO 999 dens = dsqrt (2.d0 * (xi**2 + yi**2 + param**2)) / 2.d0 / param GO TO 200 999 CONTINUE ! IF (iadap == 1) GO TO 99 tempn = 0. tempd = 0. IF (ndpt == 0) then dens = 1.d0 GO TO 20 END IF DO i = 1, ndpt temp = ( (xi - xl (i) ) **2 + (yi - yl (i) ) **2) IF (temp < tolrn) GO TO 10 temp = temp**ord tempd = tempd+1. / temp tempn = tempn + val (i) / temp END DO dens = tempn / tempd GO TO 20 10 dens = val (i) 20 CONTINUE GO TO 200 ! 99 CONTINUE IF (ndpt == 0) then dens = 1.d0 GO TO 40 END IF j = 0 DO i = 1, ndpt temp = ( (xi - xl (i) ) **2 + (yi - yl (i) ) **2) IF (temp < tolrn * 0.0000001d0) GO TO 30 j = j + 1 IF (j == 1) temp1 = temp IF (j == 2) then temp2 = temp IF (temp1 <= temp2) then smin1 = temp1 smin2 = temp2 imin1 = 1 imin2 = 2 ELSE smin1 = temp2 smin2 = temp1 imin1 = 2 imin2 = 1 END IF END IF IF (j == 3) then temp3 = temp IF (smin2 <= temp3) then smin3 = temp3 imin3 = 3 GO TO 29 ELSEIF (smin1 <= temp3) then smin3 = smin2 imin3 = imin2 smin2 = temp3 imin2 = 3 GO TO 29 ELSE smin3 = smin2 imin3 = imin2 smin2 = smin1 imin2 = imin1 smin1 = temp3 imin1 = 3 END IF END IF IF (j > 3) then IF (temp < smin1) then smin3 = smin2 imin3 = imin2 smin2 = smin1 imin2 = imin1 smin1 = temp imin1 = i GO TO 29 END IF IF (temp < smin2) then smin3 = smin2 imin3 = imin2 smin2 = temp imin2 = i GO TO 29 END IF IF (temp < smin3) then smin3 = temp imin3 = i END IF END IF 29 CONTINUE END DO ! *** interpolated value from the three nearest points tempd = (1.0d0 / smin1**ord) + (1.0d0 / smin1**ord) + (1.0d0 / & smin1**ord) tempn = (val (imin1) / smin1**ord) + (val (imin2) / smin2**ord) & + (val (imin3) / smin3**ord) dens = tempn / tempd ! *** minimum value ! dens=dmin1(val(imin1),val(imin2),val(imin3)) ! *** average of the two minimum values ! d1=dmax1(val(imin1),val(imin2),val(imin3)) ! if (d1 == val(imin1)) then ! dens=(val(imin2)+val(imin3))*0.5d0 ! else if (d1 == val(imin2)) then ! dens=(val(imin1)+val(imin3))*0.5d0 ! else if (d1 == val(imin3)) then ! dens=(val(imin2)+val(imin1))*0.5d0 ! end if GO TO 40 30 dens = val (i) 40 CONTINUE ! 200 CONTINUE END SUBROUTINE densit ! ! ! ********************************************************************** ! ! ** subroutine minsiz ! ! ********************************************************************** ! SUBROUTINE minsiz (xi, yi, ri, xl, yl, val, ndpt, tolrn, siz, ord,& szmin, mden, iadap) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) DIMENSION xl (mden), yl (mden), val (mden) ! ! *** this subroutine finds minimal mesh size requirement ! *** from the information on the density points within the square. ! ! *** basic mesh case ! IF (iadap == 0) then tempn = 0. tempd = 0. IF (ndpt == 0) then dens = 1.d0 GO TO 120 END IF DO i = 1, ndpt temp = ( (xi - xl (i) ) **2 + (yi - yl (i) ) **2) IF (temp < tolrn) GO TO 110 temp = temp**ord tempd = tempd+1. / temp tempn = tempn + val (i) / temp END DO dens = tempn / tempd GO TO 120 110 dens = val (i) 120 CONTINUE szmin = dens GO TO 100 END IF ! ! *** adaptive mesh case ! xtmax = xi + ri / 2.d0 xtmin = xi - ri / 2.d0 ytmax = yi + ri / 2.d0 ytmin = yi - ri / 2.d0 szmin = 1.0e+6 IF (ndpt == 0) then dens = 1.d0 GO TO 20 END IF DO i = 1, ndpt IF (xl (i) > xtmax .or. xl (i) < xtmin) GO TO 600 IF (yl (i) > ytmax .or. yl (i) < ytmin) GO TO 600 IF (val (i) < szmin) szmin = val (i) 600 CONTINUE END DO ! IF (szmin /= 1.0e+6) GO TO 100 tempn = 0. tempd = 0. rinf = siz 50 ninf = 0 DO i = 1, ndpt temp = ( (xi - xl (i) ) **2 + (yi - yl (i) ) **2) IF (temp < tolrn) GO TO 10 IF (temp > rinf) GO TO 60 temp = temp**ord tempd = tempd+1. / temp tempn = tempn + val (i) / temp ninf = ninf + 1 60 CONTINUE END DO IF (ninf == 0) then rinf = 2.d0 * rinf GO TO 50 END IF dens = tempn / tempd GO TO 20 10 dens = val (i) 20 CONTINUE szmin = dens ! 100 CONTINUE END SUBROUTINE minsiz !*************************************************************** ! ! ** subroutine refsiz ! ! ************************************************************** ! SUBROUTINE refsiz (xmax, xmin, ymax, ymin, xl, yl, val, ndpt, siz,& idpth, mden) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) DIMENSION xl (mden), yl (mden), val (mden) ! ! *** this subroutine finds reference mesh size for the current ! *** subdomain in order to determine inital cell division depth. ! root2 = 1.41421 szref = 0.d0 DO i = 1, ndpt IF (xl (i) > xmax .or. xl (i) < xmin) GO TO 60 IF (yl (i) > ymax .or. yl (i) < ymin) GO TO 60 IF (val (i) > szref) szref = val (i) 60 CONTINUE END DO ! IF (szref /= 0.d0) GO TO 100 DO i = 1, ndpt IF (val (i) > szref) szref = val (i) END DO IF (szref == 0.d0) szref = 1.d0 ! 100 CONTINUE xrang = xmax - xmin yrang = ymax - ymin range = dmax1 (xrang, yrang) ! write (13,*)'xmin=',xmin ! write (13,*)'xmax =',xmax ! write (13,*)'xrange=',xrang ! write (13,*)'yrange=',yrang ! write (13,*)'range=',range ! write (13,*)'siz=',siz ! write (13,*)'szref=',szref ! temp1 = dlog (range / root2 / siz / szref) IF (temp1 < 0.d0) temp1 = 0.d0 temp2 = dlog (2.d0) !b alog (2.d0) idpth = int (temp1 / temp2) ! write (13,*)'temp1=',temp1 ! write (13,*)'temp2=',temp2 ! write (13,*)'idpth=',idpth ! END SUBROUTINE refsiz ! !*********************************************************************** ! ! ** subroutine quadgen (to convert triangular elements to quads) ! ! ********************************************************************** ! SUBROUTINE quadgen (mpoin, melem, nel6, npoin, lmat, lnods, lnod4,& lma4, mbtri, mbqua, x, y) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) INTEGER, PARAMETER :: mseg = 14999, mfac = 29999, & mele4 = 9999, mpoi4 = 9999 DIMENSION x (mpoin), y (mpoin), lnods (melem, 6), lnod4 (melem, 4)& , mbtri (melem), iboun (mpoi4), lmat (melem), angle (mele4, 3), & quali (mseg), jqual (mseg), ielmf (mseg), mbqua (melem), jelmf ( & mseg), ielm (mseg), jelm (mseg), iflag (mele4), lma4 (melem), & nodf (mfac, 3), ncof (mfac), lelf (mfac), lelfn (mfac), lmatf ( & mfac) DO i = 1, mele4 iflag (i) = 0 END DO DO i = 1, nel6 i1 = lnods (i, 1) i2 = lnods (i, 3) i3 = lnods (i, 5) x1 = x (i1) x2 = x (i2) x3 = x (i3) y1 = y (i1) y2 = y (i2) y3 = y (i3) d1 = dsqrt ( (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) ) d2 = dsqrt ( (x2 - x3) * (x2 - x3) + (y2 - y3) * (y2 - y3) ) d3 = dsqrt ( (x3 - x1) * (x3 - x1) + (y3 - y1) * (y3 - y1) ) !b temp1=(d1*d1+d3*d3-d2*d2)/(2.0d0*d1*d3) !b temp2=(d1*d1+d2*d2-d3*d3)/(2.0d0*d1*d2) !b temp3=(d2*d2+d3*d3-d1*d1)/(2.0d0*d2*d3) temp1 = (d1 * d1 + d3 * d3 - d2 * d2) / (dabs (2.0d0 * d1 * d3) & + 1.d-8) temp2 = (d1 * d1 + d2 * d2 - d3 * d3) / (dabs (2.0d0 * d1 * d2) & + 1.d-8) temp3 = (d2 * d2 + d3 * d3 - d1 * d1) / (dabs (2.0d0 * d2 * d3) & + 1.d-8) angle (i, 1) = dacos (temp1) angle (i, 2) = dacos (temp2) angle (i, 3) = dacos (temp3) END DO CALL coinf6 (melem, mele4, mseg, mfac, nel6, ncoin, ncof, lelf, & lelfn, lmat, lmatf, lnods, nodf, ielmf, jelmf, ielm, jelm, angle, & quali) DO i = 1, ncoin jqual (i) = i END DO CALL shell (mseg, ncoin, jqual, quali) nel4 = 0 DO i = 1, ncoin j = ncoin - i + 1 IF (quali (j) == 0.d0) GO TO 100 iq = jqual (j) iel = ielm (iq) IF (iflag (iel) == 1) GO TO 99 jel = jelm (iq) IF (iflag (jel) == 1) GO TO 99 ielf = ielmf (iq) jelf = jelmf (iq) i1 = 2 * ielf - 1 i2 = i1 + 1 i3 = i2 + 1 IF (i3 > 6) i3 = 1 i4 = i3 + 1 i5 = i4 + 1 IF (i5 > 6) i5 = 1 i6 = i5 + 1 j1 = 2 * jelf - 1 j2 = j1 + 1 j3 = j2 + 1 IF (j3 > 6) j3 = 1 j4 = j3 + 1 j5 = j4 + 1 IF (j5 > 6) j5 = 1 j6 = j5 + 1 nel4 = nel4 + 1 lnod4 (nel4, 1) = lnods (iel, i1) lnod4 (nel4, 2) = lnods (jel, j4) lnod4 (nel4, 3) = lnods (iel, i2) lnod4 (nel4, 4) = lnods (iel, i6) lma4 (nel4) = lmat (iel) IF (mbtri (iel) /= 0) then IF (mbtri (iel) == 1) then IF (i6 == 2 .and. i1 == 3) then mbqua (nel4) = 1 nod1 = lnods (iel, i6) nod2 = lnods (iel, i1) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSEIF (mbtri (iel) == 2) then IF (i6 == 4 .and. i1 == 5) then mbqua (nel4) = 1 nod1 = lnods (iel, i6) nod2 = lnods (iel, i1) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSEIF (mbtri (iel) == 3) then IF (i6 == 6 .and. i1 == 1) then mbqua (nel4) = 1 nod1 = lnods (iel, i6) nod2 = lnods (iel, i1) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSE mbqua (nel4) = 1 nod1 = lnods (iel, i6) nod2 = lnods (iel, i1) iboun (nod1) = 1 iboun (nod2) = 1 END IF END IF IF (mbtri (jel) /= 0) then IF (mbtri (jel) == 1) then IF (j3 == 1 .and. j4 == 2) then mbqua (nel4) = 1 nod1 = lnods (jel, j3) nod2 = lnods (jel, j4) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSEIF (mbtri (jel) == 2) then IF (j3 == 3 .and. j4 == 4) then mbqua (nel4) = 1 nod1 = lnods (jel, j3) nod2 = lnods (jel, j4) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSEIF (mbtri (jel) == 3) then IF (j3 == 5 .and. j4 == 6) then mbqua (nel4) = 1 nod1 = lnods (jel, j3) nod2 = lnods (jel, j4) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSE mbqua (nel4) = 1 nod1 = lnods (jel, j3) nod2 = lnods (jel, j4) iboun (nod1) = 1 iboun (nod2) = 1 END IF END IF nel4 = nel4 + 1 !b xxx check against max lnod4 (nel4, 1) = lnods (jel, j4) lnod4 (nel4, 2) = lnods (jel, j5) lnod4 (nel4, 3) = lnods (jel, j6) lnod4 (nel4, 4) = lnods (jel, j2) lma4 (nel4) = lmat (iel) IF (mbtri (jel) /= 0) then IF (mbtri (jel) == 1) then IF (j4 == 2 .and. j5 == 3) then mbqua (nel4) = 1 nod1 = lnods (jel, j4) nod2 = lnods (jel, j5) iboun (nod1) = 1 iboun (nod2) = 1 END IF IF (j5 == 1 .and. j6 == 2) then mbqua (nel4) = 1 nod1 = lnods (jel, j5) nod2 = lnods (jel, j6) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSEIF (mbtri (jel) == 2) then IF (j4 == 4 .and. j5 == 5) then mbqua (nel4) = 1 nod1 = lnods (jel, j4) nod2 = lnods (jel, j5) iboun (nod1) = 1 iboun (nod2) = 1 END IF IF (j5 == 3 .and. j6 == 4) then mbqua (nel4) = 1 nod1 = lnods (jel, j5) nod2 = lnods (jel, j6) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSEIF (mbtri (jel) == 3) then IF (j4 == 6 .and. j5 == 1) then mbqua (nel4) = 1 nod1 = lnods (jel, j4) nod2 = lnods (jel, j5) iboun (nod1) = 1 iboun (nod2) = 1 END IF IF (j5 == 5 .and. j6 == 6) then mbqua (nel4) = 1 nod1 = lnods (jel, j5) nod2 = lnods (jel, j6) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSE mbqua (nel4) = 1 nod1 = lnods (jel, j4) nod2 = lnods (jel, j5) nod3 = lnods (jel, j6) iboun (nod1) = 1 iboun (nod2) = 1 iboun (nod3) = 1 END IF END IF nel4 = nel4 + 1 !b xxx check against max lnod4 (nel4, 1) = lnods (jel, j1) lnod4 (nel4, 2) = lnods (iel, i4) lnod4 (nel4, 3) = lnods (jel, j2) lnod4 (nel4, 4) = lnods (jel, j6) lma4 (nel4) = lmat (iel) IF (mbtri (jel) /= 0) then IF (mbtri (jel) == 1) then IF (j6 == 2 .and. j1 == 3) then mbqua (nel4) = 1 nod1 = lnods (jel, j6) nod2 = lnods (jel, j1) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSEIF (mbtri (jel) == 2) then IF (j6 == 4 .and. j1 == 5) then mbqua (nel4) = 1 nod1 = lnods (jel, j6) nod2 = lnods (jel, j1) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSEIF (mbtri (jel) == 3) then IF (j6 == 6 .and. j1 == 1) then mbqua (nel4) = 1 nod1 = lnods (jel, j6) nod2 = lnods (jel, j1) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSE mbqua (nel4) = 1 nod1 = lnods (jel, j6) nod2 = lnods (jel, j1) iboun (nod1) = 1 iboun (nod2) = 1 END IF END IF IF (mbtri (iel) /= 0) then IF (mbtri (iel) == 1) then IF (i3 == 1 .and. i4 == 2) then mbqua (nel4) = 1 nod1 = lnods (iel, i3) nod2 = lnods (iel, i4) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSEIF (mbtri (iel) == 2) then IF (i3 == 3 .and. i4 == 4) then mbqua (nel4) = 1 nod1 = lnods (iel, i3) nod2 = lnods (iel, i4) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSEIF (mbtri (iel) == 3) then IF (i3 == 5 .and. i4 == 6) then mbqua (nel4) = 1 nod1 = lnods (iel, i3) nod2 = lnods (iel, i4) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSE mbqua (nel4) = 1 nod1 = lnods (iel, i3) nod2 = lnods (iel, i4) iboun (nod1) = 1 iboun (nod2) = 1 END IF END IF nel4 = nel4 + 1 !b xxx check against max lnod4 (nel4, 1) = lnods (iel, i4) lnod4 (nel4, 2) = lnods (iel, i5) lnod4 (nel4, 3) = lnods (iel, i6) lnod4 (nel4, 4) = lnods (iel, i2) lma4 (nel4) = lmat (iel) IF (mbtri (iel) /= 0) then IF (mbtri (iel) == 1) then IF (i4 == 2 .and. i5 == 3) then mbqua (nel4) = 1 nod1 = lnods (iel, i4) nod2 = lnods (iel, i5) iboun (nod1) = 1 iboun (nod2) = 1 END IF IF (i5 == 1 .and. i6 == 2) then mbqua (nel4) = 1 nod1 = lnods (iel, i5) nod2 = lnods (iel, i6) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSEIF (mbtri (iel) == 2) then IF (i4 == 4 .and. i5 == 5) then mbqua (nel4) = 1 nod1 = lnods (iel, i4) nod2 = lnods (iel, i5) iboun (nod1) = 1 iboun (nod2) = 1 END IF IF (i5 == 3 .and. i6 == 4) then mbqua (nel4) = 1 nod1 = lnods (iel, i5) nod2 = lnods (iel, i6) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSEIF (mbtri (iel) == 3) then IF (i4 == 6 .and. i5 == 1) then mbqua (nel4) = 1 nod1 = lnods (iel, i4) nod2 = lnods (iel, i5) iboun (nod1) = 1 iboun (nod2) = 1 END IF IF (i5 == 5 .and. i6 == 6) then mbqua (nel4) = 1 nod1 = lnods (iel, i5) nod2 = lnods (iel, i6) iboun (nod1) = 1 iboun (nod2) = 1 END IF ELSE mbqua (nel4) = 1 nod1 = lnods (iel, i4) nod2 = lnods (iel, i5) nod3 = lnods (iel, i6) iboun (nod1) = 1 iboun (nod2) = 1 iboun (nod3) = 1 END IF END IF iflag (iel) = 1 iflag (jel) = 1 inod = lnods (iel, i2) in1 = lnods (iel, i1) in2 = lnods (jel, j1) in3 = lnods (iel, i5) in4 = lnods (jel, j5) x (inod) = (x (in1) + x (in2) + x (in3) + x (in4) ) * 0.25d0 y (inod) = (y (in1) + y (in2) + y (in3) + y (in4) ) * 0.25d0 99 CONTINUE END DO 100 CONTINUE npoi4 = npoin DO i = 1, nel6 IF (iflag (i) == 1) GO TO 999 npoi4 = npoi4 + 1 i1 = lnods (i, 1) i2 = lnods (i, 2) i3 = lnods (i, 3) i4 = lnods (i, 4) i5 = lnods (i, 5) i6 = lnods (i, 6) IF (mbtri (i) /= 0) then IF (mbtri (i) == 1) then iboun (i1) = 1 iboun (i2) = 1 iboun (i3) = 1 END IF IF (mbtri (i) == 2) then iboun (i3) = 1 iboun (i4) = 1 iboun (i5) = 1 END IF IF (mbtri (i) == 3) then iboun (i5) = 1 iboun (i6) = 1 iboun (i1) = 1 END IF IF (mbtri (i) == 4) then iboun (i1) = 1 iboun (i2) = 1 iboun (i3) = 1 iboun (i4) = 1 iboun (i5) = 1 END IF IF (mbtri (i) == 5) then iboun (i3) = 1 iboun (i4) = 1 iboun (i5) = 1 iboun (i6) = 1 iboun (i1) = 1 END IF IF (mbtri (i) == 6) then iboun (i5) = 1 iboun (i6) = 1 iboun (i1) = 1 iboun (i2) = 1 iboun (i3) = 1 END IF END IF x (npoi4) = (x (i1) + x (i3) + x (i5) ) / 3.d0 y (npoi4) = (y (i1) + y (i3) + y (i5) ) / 3.d0 nel4 = nel4 + 1 !b xxx check against max lnod4 (nel4, 1) = i1 lnod4 (nel4, 2) = i2 lnod4 (nel4, 3) = npoi4 lnod4 (nel4, 4) = i6 lma4 (nel4) = lmat (i) IF (mbtri (i) /= 0) mbqua (nel4) = 1 nel4 = nel4 + 1 !b xxx check against max lnod4 (nel4, 1) = i3 lnod4 (nel4, 2) = i4 lnod4 (nel4, 3) = npoi4 lnod4 (nel4, 4) = i2 lma4 (nel4) = lmat (i) IF (mbtri (i) /= 0) mbqua (nel4) = 1 nel4 = nel4 + 1 !b xxx check against max lnod4 (nel4, 1) = i5 lnod4 (nel4, 2) = i6 lnod4 (nel4, 3) = npoi4 lnod4 (nel4, 4) = i4 lma4 (nel4) = lmat (i) IF (mbtri (i) /= 0) mbqua (nel4) = 1 999 CONTINUE END DO npoin = npoi4 nel6 = nel4 END SUBROUTINE quadgen ! ! ******************************************************************** ! SUBROUTINE coinf6 (melem, mele4, mseg, mfac, nel6, ncoin, ncof, & lelf, lelfn, lmat, lmatf, lnods, nodf, ielmf, jelmf, ielm, jelm, & angle, quali) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) DIMENSION lnods (melem, 6), nodf (mfac, 3), ncof (mfac), & lelf (mfac), lelfn (mfac), angle (mele4, 3), quali (mseg), & ielmf (mseg), jelmf (mseg), ielm (mseg), jelm (mseg), & lmat (melem), lmatf (mfac) ! ! *** identify all the element faces. ! nfac = 0 DO i = 1, nel6 DO j = 1, 3 nfac = nfac + 1 !b xxx check against max lelf (nfac) = i ncof (nfac) = 0 lelf (nfac) = i lelfn (nfac) = j lmatf (nfac) = lmat (i) DO k = 1, 3 k1 = 2 * j + k - 2 IF (k1 > 6) k1 = 1 nodf (nfac, k) = lnods (i, k1) END DO END DO END DO ! ! *** identify coincident faces and establish quality criteria ! *** for triangle pairs. ! ncoin = 0 DO i = 1, nfac IF (ncof (i) /= 0) GO TO 10 DO j = 1, nfac IF (ncof (j) /= 0 .or. i == j) GO TO 20 nod1 = nodf (i, 2) nod2 = nodf (j, 2) mat1 = lmatf (i) mat2 = lmatf (j) IF ( (nod1 /= nod2) .or. (mat1 /= mat2) ) GO TO 20 ncof (i) = j ncof (j) = i ncoin = ncoin + 1 !b xxx check against max iel = lelf (i) jel = lelf (j) ieln = lelfn (i) iel1 = ieln + 1 IF (iel1 > 3) iel1 = 1 jeln = lelfn (j) jel1 = jeln + 1 IF (jel1 > 3) jel1 = 1 angl1 = angle (iel, ieln) + angle (jel, jel1) angl2 = angle (jel, jeln) + angle (iel, iel1) qtem1 = angl1 / 3.1415927d0 IF (qtem1 > 0.9d0) qtem1 = 0.0d0 IF (qtem1 > 0.5d0) qtem1 = 1.0d0 - qtem1 qtem2 = angl2 / 3.1415927d0 IF (qtem2 > 0.9d0) qtem2 = 0.0d0 IF (qtem2 > 0.5d0) qtem2 = 1.0d0 - qtem2 quali (ncoin) = qtem1 + qtem2 IF (qtem1 == 0.0d0 .or. qtem2 == 0.0d0) quali (ncoin) = 0.0d0 ielm (ncoin) = iel jelm (ncoin) = jel ielmf (ncoin) = ieln jelmf (ncoin) = jeln GO TO 10 20 CONTINUE END DO 10 CONTINUE END DO END SUBROUTINE coinf6 ! ! ********************************************************************** ! SUBROUTINE shell (m, n, ia, a) USE PRECISION_MODULE IMPLICIT real (DP)(a - h, o - z) DIMENSION a (m), ia (m) ! inc = 1 ! 10 IF (inc > n) GO TO 20 inc = 3 * inc + 1 IF (inc <= n) GO TO 10 ! 20 CONTINUE IF (inc <= 1) GO TO 100 inc = inc / 3 DO i = inc + 1, n v = a (i) iv = ia (i) j = i 40 IF (a (j - inc) <= v) GO TO 50 a (j) = a (j - inc) ia (j) = ia (j - inc) j = j - inc IF (j <= inc) GO TO 50 IF (a (j - inc) > v) GO TO 40 50 a (j) = v ia (j) = iv END DO IF (inc > 1) GO TO 20 ! 100 CONTINUE END SUBROUTINE shell ! !*********************************************************************** ! ! ** subroutine bquad (boundary nodes for 4 and 6 node elements) ! ! ********************************************************************** ! SUBROUTINE bquad (mtfn, ntri, nbpl, npinb, nbno, npin6, & nbno6, mbtri, lnods) ! mbpl mbpl, mbpl mbpl, DIMENSION npinb (189), nbno (189, 489), npin6 (189), & nbno6 (189, 979), mbtri (mtfn), lnods (mtfn, 6) ! mbpl, , mnpel DO i = 1, nbpl npin6 (i) = iabs (npinb (i) ) + iabs (npinb (i) ) - 1 IF (npinb (i) < 0) npin6 (i) = - npin6 (i) np = iabs (npinb (i) ) - 1 IF (np < 1) GO TO 100 DO j = 1, np jn1 = nbno (i, j) jn2 = nbno (i, j + 1) j1 = 2 * j - 1 j2 = 2 * j nbno6 (i, j1) = jn1 DO k = 1, ntri IF (mbtri (k) /= 0) then IF (mbtri (k) == 1) then kn1 = lnods (k, 1) kn2 = lnods (k, 3) IF (jn1 == kn1 .and. jn2 == kn2) jn = lnods (k, 2) IF (jn2 == kn1 .and. jn1 == kn2) jn = lnods (k, 2) END IF IF (mbtri (k) == 2) then kn1 = lnods (k, 3) kn2 = lnods (k, 5) IF (jn1 == kn1 .and. jn2 == kn2) jn = lnods (k, 4) IF (jn2 == kn1 .and. jn1 == kn2) jn = lnods (k, 4) END IF IF (mbtri (k) == 3) then kn1 = lnods (k, 5) kn2 = lnods (k, 1) IF (jn1 == kn1 .and. jn2 == kn2) jn = lnods (k, 6) IF (jn2 == kn1 .and. jn1 == kn2) jn = lnods (k, 6) END IF IF (mbtri (k) == 4) then kn1 = lnods (k, 1) kn2 = lnods (k, 3) IF (jn1 == kn1 .and. jn2 == kn2) jn = lnods (k, 2) IF (jn2 == kn1 .and. jn1 == kn2) jn = lnods (k, 2) kn1 = lnods (k, 3) kn2 = lnods (k, 5) IF (jn1 == kn1 .and. jn2 == kn2) jn = lnods (k, 4) IF (jn2 == kn1 .and. jn1 == kn2) jn = lnods (k, 4) END IF IF (mbtri (k) == 5) then kn1 = lnods (k, 3) kn2 = lnods (k, 5) IF (jn1 == kn1 .and. jn2 == kn2) jn = lnods (k, 4) IF (jn2 == kn1 .and. jn1 == kn2) jn = lnods (k, 4) kn1 = lnods (k, 5) kn2 = lnods (k, 1) IF (jn1 == kn1 .and. jn2 == kn2) jn = lnods (k, 6) IF (jn2 == kn1 .and. jn1 == kn2) jn = lnods (k, 6) END IF IF (mbtri (k) == 6) then kn1 = lnods (k, 1) kn2 = lnods (k, 3) IF (jn1 == kn1 .and. jn2 == kn2) jn = lnods (k, 2) IF (jn2 == kn1 .and. jn1 == kn2) jn = lnods (k, 2) kn1 = lnods (k, 5) kn2 = lnods (k, 1) IF (jn1 == kn1 .and. jn2 == kn2) jn = lnods (k, 6) IF (jn2 == kn1 .and. jn1 == kn2) jn = lnods (k, 6) END IF END IF END DO nbno6 (i, j2) = jn IF (j == np) nbno6 (i, j2 + 1) = jn2 END DO 100 CONTINUE END DO END SUBROUTINE bquad ! !*********************************************************************** ! ! ** subroutines for front optimisation ! ! ********************************************************************** ! SUBROUTINE FROMIN (MELEM, MNPEL, NPOIN, NELEM, NNODE, LNODS, & MATNO, NEN, nfold, nfnew) ! ! Reference parameter settings are as follows: ! ! melm = melem ! melm2 = 10*melem ! mpoin = mpoin ! mnpl = mnpel ! PARAMETER (melm = 9999, melm2 = 99990, mnpl = 9, mpoin = 9999) DIMENSION LNODS (MELEM, MNPEL), NODST (melm, mnpl), MATNO (MELEM) DIMENSION LIST (MPOIN), NDEG (MPOIN), NSTART (MPOIN), LEV (MPOIN),& NEWNUM (MPOIN), NADJ (melm2), IENM (melm), NEN (melem) ! DIMEN = (MELEM-1)*2+4 LOGICAL ALL MAXDEG = 9 MINMAX = 2000 ! Optimising scheme IALGOR = 0 ! ! Optimize front width INC = 1 ALL = .FALSE. IF (NNODE == 4) ALL = .TRUE. IF (NNODE >= 6) INC = 2 CALL MAXFRO (melem, mnpel, mpoin, melm, melm2, mnpl, nelem, npoin,& nnode, lnods, ienm, NFOLD) CALL CORNER (melem, mnpel, mpoin, nelem, nnode, lnods, list, INC, & NODCR) CALL RESEQ3 (melem, mnpel, mpoin, melm, mnpl, nelem, nnode, lnods,& nodst, list, NODCR) CALL SETUP (melem, mnpel, mpoin, melm2, nelem, nnode, lnods, ndeg,& nadj, NODCR, MAXDEG, INC, ALL) CALL DIAM (mpoin, melm2, nstart, ndeg, lev, nadj, NODCR, MAXDEG, & NS) CALL RESEQ1 (mpoin, melm2, ndeg, nstart, lev, newnum, nadj, NODCR,& MAXDEG, NS, MINMAX, IALGOR) CALL RESEQ2 (melem, mnpel, mpoin, melm, mnpl, nelem, nnode, lnods,& nodst, matno, newnum, ienm, nen, NODCR, INC, ALL) CALL MAXFRO (melem, mnpel, mpoin, melm, melm2, mnpl, nelem, npoin,& nnode, lnods, ienm, NFNEW) END SUBROUTINE FROMIN SUBROUTINE LEVEL (mpoin, melm2, ndeg, lev, nadj, LSD, MLW, NODCR, & NROOT, MAXDEG) !**************************************************************** ! Subprogram LEVEL - compute level structure rooted at NROOT !***************************************************************** DIMENSION NDEG (MPOIN), LEV (MPOIN), NADJ (melm2) ! Initialisation DO I = 1, NODCR LEV (I) = 0 END DO LEV (NROOT) = 1 KOUNT = 1 MLW = 1 ! Assign levels to vertices DO L = 2, NODCR LW = 0 DO 30 I = 1, NODCR IF (LEV (I) > 0) GO TO 30 NCS = NDEG (I) JSUB = (I - 1) * MAXDEG DO 20 JJ = 1, NCS NODE = NADJ (JSUB + JJ) IF (LEV (NODE) /= L - 1) GO TO 20 LSD = L LW = LW + 1 LEV (I) = L KOUNT = KOUNT + 1 IF (KOUNT == NODCR) GO TO 50 GO TO 30 20 END DO 30 END DO IF (LW > MLW) MLW = LW END DO 50 IF (LW > MLW) MLW = LW END SUBROUTINE LEVEL SUBROUTINE RESEQ2 (melem, mnpel, mpoin, melm, mnpl, nelem, nnode, & lnods, nodst, matno, newnum, ienm, nen, NODCR, INC, ALL) !********************************************************************* ! Subprogram RESEQ2 - Resquence element numbers to minimise ! the frontwidth ! - Reorder the elements in an ascending sequence ! of their lowest numbered nodes !********************************************************************** DIMENSION LNODS (MELEM, MNPEL), NODST (melm, mnpl), MATNO (MELEM),& NEWNUM (MPOIN), IENM (melm), NEN (melem) LOGICAL ALL DO I = 1, NELEM NEN (I) = 0 END DO KOUNT = 0 ! Loop over each new node number ! Loop only over corner nodes if ALL=.FALSE. DO I = 1, NODCR ! Loop over each element ! Skip to next element if already renumbered DO 30 J = 1, NELEM IF (NEN (J) > 0) GO TO 30 IF (ALL) INC = 1 ! Loop over each node in element ! Use only corner nodes if ALL=.FALSE. ! Assumed that corner nodes are listed first in nodal definition ! vectors if ALL=.FALSE. DO 20 K = 1, NNODE, INC N = LNODS (J, K) N = NEWNUM (N) IF (N /= I) GO TO 20 KOUNT = KOUNT + 1 NEN (J) = KOUNT IF (KOUNT == NELEM) GO TO 50 GO TO 30 20 END DO 30 END DO END DO 50 CONTINUE DO I = 1, NELEM IENM (NEN (I) ) = MATNO (I) DO J = 1, NNODE LNODS (NEN (I), J) = NODST (I, J) END DO END DO DO 65 I = 1, NELEM MATNO (I) = IENM (I) 65 END DO END SUBROUTINE RESEQ2 SUBROUTINE MAXFRO (melem, mnpel, mpoin, melm, melm2, mnpl, nelem, & npoin, nnode, lnods, ienm, KFRON) !***************************************************************** ! Subroutine MAXFRO - Computes the maximum frontwidth !****************************************************************** DIMENSION LNODS (MELEM, MNPEL), IENM (melm) NDOFN = 2 ! Zero NDFRO vector. DO IELEM = 1, NELEM IENM (IELEM) = 0 END DO DO IPOIN = 1, NPOIN KSTAR = 0 DO IELEM = 1, NELEM DO 90 INODE = 1, NNODE IF (LNODS (IELEM, INODE) /= IPOIN) GO TO 90 IF (KSTAR /= 0) GO TO 80 KSTAR = IELEM IENM (IELEM) = IENM (IELEM) + NDOFN 80 CONTINUE KLAST = IELEM NLAST = INODE 90 END DO END DO IF (KLAST < NELEM) IENM (KLAST + 1) = IENM (KLAST + 1) - NDOFN LNODS (KLAST, NLAST) = - IPOIN END DO NFRON = 0 KFRON = 0 DO IELEM = 1, NELEM NFRON = NFRON + IENM (IELEM) IF (NFRON > KFRON) KFRON = NFRON END DO DO IELEM = 1, NELEM DO INODE = 1, NNODE LNODS (IELEM, INODE) = IABS (LNODS (IELEM, INODE) ) END DO END DO END SUBROUTINE MAXFRO SUBROUTINE SETUP (melem, mnpel, mpoin, melm2, nelem, nnode, lnods,& ndeg, nadj, NODCR, MAXDEG, INC, ALL) !************************************************************ ! Subroutine SETUP - Compute adjancency list and degree node ! - Modified version of Collins routine ! - Use only corner nodes if ALL=.FALSE. !*************************************************************** DIMENSION LNODS (MELEM, MNPEL), NDEG (MPOIN), NADJ (melm2) LOGICAL ALL DO J = 1, NODCR NDEG (J) = 0 END DO DO J = 1, NELEM IF (ALL) INC = 1 DO I = 1, NNODE, INC JNTI = LNODS (J, I) JSUB = (JNTI - 1) * MAXDEG DO 40 II = 1, NNODE, INC IF (II == I) GO TO 40 JJT = LNODS (J, II) MEM1 = NDEG (JNTI) IF (MEM1 == 0) GO TO 30 DO III = 1, MEM1 IF (NADJ (JSUB + III) == JJT) GO TO 40 END DO 30 NDEG (JNTI) = NDEG (JNTI) + 1 NADJ (JSUB + NDEG (JNTI) ) = JJT 40 END DO END DO END DO END SUBROUTINE SETUP SUBROUTINE DIAM (mpoin, melm2, nstart, ndeg, lev, nadj, NODCR, & MAXDEG, NS) !****************************************************************** ! Subroutine DIAM - Compute set of psuedo-peripheral nodes !****************************************************************** DIMENSION NSTART (MPOIN), LEV (MPOIN), NDEG (MPOIN), NADJ (melm2) LOGICAL BETTER ! Begin iteration : select intial root node arbitarily ! and generate its level structure IROOT = 1 ITER = 0 10 ITER = ITER + 1 CALL LEVEL (mpoin, melm2, ndeg, lev, nadj, IDEPTH, IWIDTH, NODCR, & IROOT, MAXDEG) ! Create list of nodes which are at maxium distance from root node LHW = 0 DO 20 I = 1, NODCR IF (LEV (I) /= IDEPTH) GO TO 20 LHW = LHW + 1 !b xxx check against max NSTART (LHW) = I 20 END DO ! Store root on end of list of possible starting nodes NS = LHW + 1 !b xxx check against max NSTART (NS) = IROOT ! Loop over nodes at maxium distance from root node ! GENERATE level structure for each node ! Set switch if a level structure of greater depth occurs BETTER = .FALSE. DO 30 I = 1, LHW NEND = NSTART (I) CALL LEVEL (mpoin, melm2, ndeg, lev, nadj, NDEPTH, NWIDTH, & NODCR, NEND, MAXDEG) IF (NDEPTH < IDEPTH) GO TO 30 IF ( (NDEPTH == IDEPTH) .AND. (NWIDTH >= IWIDTH) ) GO TO 30 IROOT = NEND IDEPTH = NDEPTH IWIDTH = NWIDTH BETTER = .TRUE. 30 END DO IF (BETTER) GO TO 10 END SUBROUTINE DIAM SUBROUTINE CORNER (melem, mnpel, mpoin, nelem, nnode, lnods, list,& INC, NODCR) !********************************************************************** ! Subprogram CORNER - Assemble list of corner nodes for high order mesh ! with a single type of element ! - Also compute total number of corner nodes !********************************************************************** DIMENSION LNODS (MELEM, MNPEL), LIST (MPOIN) ! Initialise NODCR = 1 LIST (1) = LNODS (1, 1) ! Loop over corner nodes for each element DO IE = 1, NELEM DO 50 I = 1, NNODE, INC N = LNODS (IE, I) ! Check if corner node is already in list - if not, add it in DO J = 1, NODCR IF (N == LIST (J) ) GO TO 50 END DO NODCR = NODCR + 1 !b xxx check against max LIST (NODCR) = N 50 END DO END DO END SUBROUTINE CORNER SUBROUTINE RESEQ3 (melem, mnpel, mpoin, melm, mnpl, nelem, nnode, & lnods, nodst, list, NODCR) !******************************************************************* ! Subprogram RESEQ3 - Reorder corner nodes in ascending sequence !******************************************************************* DIMENSION LNODS (MELEM, MNPEL), NODST (melm, mnpl), LIST (MPOIN) DO J = 1, NELEM DO L = 1, NNODE NODST (J, L) = LNODS (J, L) END DO END DO ! Loop over each corner node DO 100 I = 1, NODCR ! Check if node I is already a corner node - if so,skip to next I DO J = 1, NODCR IF (LIST (J) == I) GO TO 100 END DO ! Scan through list of corner nodes ! If number of node K is > NODCR , swap node K and node I ! in nodal topology vectors for elements DO 60 K = 1, NODCR IF (LIST (K) <= NODCR) GO TO 60 DO IE = 1, NELEM DO 25 J = 1, NNODE N = LNODS (IE, J) IF (N == LIST (K) ) GO TO 20 IF (N == I) GO TO 21 GO TO 25 20 LNODS (IE, J) = I GO TO 25 21 LNODS (IE, J) = LIST (K) 25 END DO END DO ! Replace node K with node I in list of corner nodes LIST (K) = I GO TO 100 60 END DO 100 END DO END SUBROUTINE RESEQ3 SUBROUTINE RESEQ1 (mpoin, melm2, ndeg, nstart, lev, newnum, nadj, & NODCR, MAXDEG, NS, MINMAX, IALGOR) !******************************************************************* ! Subprogram RESEQ1 - Resequence nodes for minimum frontwidth !****************************************************************** DIMENSION NDEG (MPOIN), NSTART (MPOIN), LEV (MPOIN), NEWNUM ( & MPOIN), NADJ (melm2) LARGE = 5**5 ! Loop over set of starting nodes DO II = 1, NS I = NSTART (II) DO J = 1, NODCR LEV (J) = 0 END DO NIF = NDEG (I) MAXFRT = NIF LEV (I) = 1 ! Negate all NDEG entries for nodes which are adjacent to ! starting node I NCN = NDEG (I) JSUB = (I - 1) * MAXDEG DO J = 1, NCN N = NADJ (JSUB + J) NDEG (N) = - NDEG (N) END DO NDEG (I) = - NDEG (I) ! Loop over nodes to be renumbered DO 60 K = 2, NODCR MINNEW = LARGE LMIN = LARGE ! Loop over unnumbered nodes ! Skip to next node if old node is already renumbered ! Restrict search to active nodes for King scheme ! Ignore this restraint for the Levy scheme DO 40 J = 1, NODCR IF (LEV (J) > 0) GO TO 40 IF (IALGOR == 0 .AND. NDEG (J) > 0) GO TO 40 NEW = - 1 MIN = LARGE NCN = IABS (NDEG (J) ) LSUB = (J - 1) * MAXDEG ! Compute the increment in active nodes for each node J ! Compute when node was first active by checking for renumbered ! neighbours with lowest numbers DO 30 L = 1, NCN N = NADJ (LSUB + L) IF (NDEG (N) > 0) NEW = NEW + 1 IF (LEV (N) == 0) GO TO 30 IF (LEV (N) < MIN) MIN = LEV (N) 30 END DO ! Select node with smallest increment in active nodes ! In the case of a tie, select node which has been active the ! longest IF (IALGOR == 1 .AND. NDEG (J) > 0) NEW = NEW + 1 IF (NEW > MINNEW) GO TO 40 IF ( (NEW == MINNEW) .AND. (MIN >= LMIN) ) GO TO 40 MINNEW = NEW LMIN = MIN NEXT = J 40 END DO ! Renumber node and compute number of active nodes ! Abondon scheme if number of active nodes exeeds previous ! lowest maximum frontwidth LEV (NEXT) = K NIF = NIF + MINNEW IF (NIF > MAXFRT) MAXFRT = NIF IF (MAXFRT >= MINMAX) GO TO 80 ! Negate all NDEG entries for nodes which are ! adjacent to node just renumbered IF (IALGOR == 1 .AND. NDEG (NEXT) > 0) NDEG (NEXT) = - NDEG ( & NEXT) IF (MINNEW == - 1) GO TO 60 NCN = IABS (NDEG (NEXT) ) JSUB = (NEXT - 1) * MAXDEG DO 50 J = 1, NCN N = NADJ (JSUB + J) IF (NDEG (N) > 0) NDEG (N) = - NDEG (N) 50 END DO 60 END DO ! Store numbering scheme generated ! Reset NDEG to positive values DO J = 1, NODCR NEWNUM (J) = LEV (J) END DO MINMAX = MAXFRT 80 DO J = 1, NODCR NDEG (J) = IABS (NDEG (J) ) END DO END DO MINMAX = MINMAX + 1 END SUBROUTINE RESEQ1 ! !*********************************************************************** ! ! ** subroutines for profile optimisation ! ! ********************************************************************** ! SUBROUTINE promin (melem, mnpel, mpoin, npoin, nelem, nnode, lm2, & newno, nprof1, nprof2) ! ! Reference parameter settings are as follows: ! ! melm = melem ! mpoi = mpoin ! mpoi2 = mpoin*20 ! mnpl = mnpel ! PARAMETER (melm = 9999, mpoi = 9999, mpoi2 = 199999, mnpl = 9) ! ! *** PROGRAM FOR BANDWIDTH MINIMISATION ! DIMENSION lm (mnpl, melm), nadj (mpoi2), lev (mpoi), ns (mpoi), & nstart (mpoi), newnn (mpoi), NEWNO (mpoin), lm2 (melem, mnpel), & ncolm (mpoi), ndiag (mpoi) DO i = 1, nelem DO j = 1, nnode lm (j, i) = lm2 (i, j) END DO END DO isbw1 = 1 DO i = 1, nelem DO j = 1, nnode-1 jnod = lm (j, i) DO k = j + 1, nnode knod = lm (k, i) jsbw = iabs (jnod-knod) IF (isbw1 < jsbw + 1) isbw1 = jsbw + 1 END DO END DO END DO DO i = 1, npoin ncolm (i) = 1 END DO DO i = 1, nelem n1 = npoin DO j = 1, nnode k = lm (j, i) IF (k < n1) n1 = k END DO DO j = 1, nnode k = lm (j, i) n2 = k - n1 + 1 IF (ncolm (k) < n2) ncolm (k) = n2 END DO END DO ndiag (1) = 1 DO j = 2, npoin ndiag (j) = ndiag (j - 1) + ncolm (j) END DO nprof1 = ndiag (npoin) CALL COMPAD (melm, mnpl, mpoi2, LM, NADJ, NELEM, NPOIN, NNODE) CALL DIAMP (mpoi, mpoi2, NSTART, NADJ, LEV, NEWNN, NPOIN, & NS, NPASS) CALL RESEQP (mpoi, mpoi2, NADJ, NEWNN, NSTART, LEV, NPOIN, & NS, NPASS, MINMAX) DO iel = 1, nelem DO inod = 1, nnode nod = lm (inod, iel) newnod = newnn (nod) lm2 (iel, inod) = newnod END DO END DO isbw2 = 1 DO i = 1, nelem DO j = 1, nnode-1 jnod = lm2 (i, j) DO k = j + 1, nnode knod = lm2 (i, k) jsbw = iabs (jnod-knod) IF (isbw2 < jsbw + 1) isbw2 = jsbw + 1 END DO END DO END DO DO i = 1, npoin ncolm (i) = 1 newno (i) = newnn (i) END DO DO i = 1, nelem n1 = npoin DO j = 1, nnode k = lm2 (i, j) IF (k < n1) n1 = k END DO DO j = 1, nnode k = lm2 (i, j) n2 = k - n1 + 1 IF (ncolm (k) < n2) ncolm (k) = n2 END DO END DO ndiag (1) = 1 DO j = 2, npoin ndiag (j) = ndiag (j - 1) + ncolm (j) END DO nprof2 = ndiag (npoin) END SUBROUTINE promin ! SUBROUTINES FOR OPTIMUM BANDWIDTH RENUMBERING ! ! ! CALL SEQUENCE: ! ! MINMAX=NPOIN ! CALL COMPAD(LM,NADJ,NELEM,NPOIN,NNODE) ! CALL DIAM(NSTART,NADJ,LEV,NEWNN,NPOIN,NS,NPASS) ! CALL RESEQ1(NADJ,NEWNN,NSTART,LEV,NPOIN,NS,NPASS,MINMAX) ! ! ! VARIABLES: ! ! NPOIN --> NUMBER OF NODES IN THE MESH ! NELEM --> NUMBER OF ELEMENTS ! NNODE --> NUMBER OF NODES PER ELEMENT ! NPASS --> AUXILARY NUMBER ! ! ARRAYS: ! ! LM(NNODE,NELEM) --> NODAL CONNECTIVITIES ! NADJ(NMAX*NPOIN) --> AUXILIARY ARRAY (NMAX: MAXIMUM NUMBER OF ! NODES SURROUNDING ANY GIVEN NODE) ! LEV(NPOIN) --> AUXILIARY ARRAY ! NS(NPOIN) --> " " ! NSTART(NPOIN) --> " " ! NEWNN(NPOIN) --> NEW NODE NUMBER OF EACH NODE (OUTPUT) ! ! !----------------------------------------------------------------------- SUBROUTINE COMPAD (melm, mnpl, mpoi2, LM, NODAD, NELEM, & NPOIN, NNODE) !----------------------------------------------------------------------- DIMENSION LM (mnpl, melm), NODAD (mpoi2) ! ! FOR EACH DEGREE OF FREEDOM CREATES A LINKED LIST CONTAINING ! THOSE D.O.F. WHICH ARE CONNECTED WITH THE FORMER. THESE LISTS ! ARE STORED WITHIN NODAD AND THEIR BEGINING IS MARKED BY THE ! FIRST NPOIN POSITIONS OF NODAD. ! DO ILIB = 1, NPOIN NODAD (ILIB) = 0 END DO NFREE = NPOIN + 1 ! ! LOOP OVER ALL THE ELEMENTS TO EXAMINE ALL THE DOF-DOF CONNECTIONS, ! OBVIOUSLY SKIPS THOSE WHICH ARE ZERO ! DO IEL = 1, NELEM DO IDF = 1, NNODE ILIB = LM (IDF, IEL) IF (ILIB > 0) THEN DO JDF = 1, NNODE IF (JDF /= IDF) THEN JLIB = LM (JDF, IEL) IF (JLIB > 0) THEN ! ! THE DEGREE OF FREEDOM JLIB IS CONNECTED TO ILIB, THUS WE SHOULD ! INCLUDE IT IN ILIB'S LIST, HOWEVER FIRST WE MUST CHECK THAT IT IS ! ALREADY THERE (FROM THE ANALYSIS OF OTHER ELEMENTS) ! NADR = NODAD (ILIB) NADRO = ILIB 1 IF (NADR > 0) THEN KLIB = NODAD (NADR) IF (KLIB == JLIB) GO TO 10 NADRO = NADR + 1 NADR = NODAD (NADRO) GO TO 1 END IF ! ! IF JLIB WAS NOT FOUND IT IS INCLUDED IN THE LIST AND THE PREVIOUS ! LINK (POINTER) SET TO THE ADDRESS IN NODAD WERE JLIB IS PUT ! NODAD (NADRO) = NFREE NODAD (NFREE) = JLIB NODAD (NFREE+1) = 0 NFREE = NFREE+2 10 CONTINUE END IF END IF END DO END IF END DO END DO END SUBROUTINE COMPAD ! !----------------------------------------------------------------------- SUBROUTINE DIAMP (mpoi, mpoi2, NSTART, NADJ, LEV, NAUX, NPOIN, & NS, IPASS) !----------------------------------------------------------------------- ! ! COMPUTES A SET OF POSIBLES STARTINGS NPOIN ! !----------------------------------------------------------------------- DIMENSION NAUX (mpoi), NADJ (mpoi2), LEV (mpoi), NSTART (mpoi), & NS (mpoi) LOGICAL BETTER ! ! INITIALIZES "LEV" TO ZERO ! DO I = 1, NPOIN LEV (I) = 0 END DO IPASS = 0 ! ! SELECT INTIAL ROOT NODE ARBITARILY ! IROOT = 1 1 IF (IROOT > 0) THEN IPASS = IPASS + 1 BETTER = .TRUE. ! ! BEGIN ITERATION, CREATE THE LEVEL STRUCTURE OF THE ROOT AND ! CREATE A LIST OF NPOIN WHICH ARE AT MAXIMUM DISTANCE FROM ROOT ! 2 IF (BETTER) THEN CALL LEVELP (mpoi, mpoi2, NSTART, LEV, IDEPTH, NADJ, IWIDTH, & NPOIN, IROOT, LHW) ! ! STORE THE ROOT ! NS (IPASS) = IROOT ! ! LOOP OVER NPOIN AT MAXIUM DISTANCE FROM ROOT NODE ! GENERATE LEVEL STRUCTURE FOR EACH NODE ! SET SWITCH IF A LEVEL STRUCTURE OF GREATER DEPTH OCCURS ! BETTER = .FALSE. DO I = 1, LHW NEND = NSTART (I) CALL LEVELP (mpoi, mpoi2, NAUX, LEV, NDEPTH, NADJ, NWIDTH, & NPOIN, NEND, LR) IF (NDEPTH >= IDEPTH) THEN IF ( (NDEPTH /= IDEPTH) .OR. (NWIDTH < IWIDTH) ) THEN IROOT = NEND IDEPTH = NDEPTH IWIDTH = NWIDTH BETTER = .TRUE. END IF END IF END DO GO TO 2 END IF ! ! IF NOT ALL THE NPOIN HAVE BEEN ASSIGNED LEVELS, THE MESH MUST ! BE MADE UP OF UNCONNECTED SUBSTRUCTURES, THEREFORE A ROOT FOR ! EACH SUBSTRUCTURE MUST BE COMPUTED. WE START BY GIVING A RANDOM ! ROOT. ! IROOT = 0 DO I = 1, NPOIN IF (LEV (I) == 0) IROOT = I END DO GO TO 1 END IF END SUBROUTINE DIAMP ! !----------------------------------------------------------------------- SUBROUTINE LEVELP (mpoi, mpoi2, NDEG, LEV, LSD, NADJ, MLW, NPOIN, & NROOT, NLOC) !----------------------------------------------------------------------- ! ! SUBPROGRAM LEVEL-COMPUTE LEVEL STRUCTURE ROOTED AT NROOT ! !----------------------------------------------------------------------- LOGICAL BACK DIMENSION LEV (mpoi), NDEG (mpoi), NADJ (mpoi2) ! ! COMPUTES THE LEVELS OF A MESH OF POINTS (OR EQUIVALENTLY OF ! DEGREES OF FREEDOMS) GIVEN A ROOT POINT. THE LEVEL OF A POINT ! IS DEFINED BY INDUCTION AS FOLLOWS. THE ROOT IS SAID TO HAVE ! LEVEL=1, THEN THE GROUP OF POINTS WITH LEVEL K IS DEFINED AS ! THOSE POINTS NOT BELONGING TO GROUPS (1,2,...K-1) AND BEING ! IMMEDIATELY CONNECTED (I.E. WITHIN AN ELEMENT'S DISTANCE) WITH ! ANY NODE OF THE LEVEL K-1. ! DO I = 1, NPOIN LEV (I) = - ABS (LEV (I) ) END DO LEV (NROOT) = 1 LSD = 1 NLOC = 1 NLOCN = 0 MLW = 1 ! ! TO EVALUATE THE LEVEL OF EACH NODE THE PREVIOUS DEFINITION ! IS STRICTLY FOLLOWED. EACH LEVEL GROUP IS EVALUATED IN TURN ! BY MEANS OF EXAMINIG ALL THE CONNECTION OF THE PREVIOUS LEVEL ! MAKING USE OF THE NODE-NODE CONNECTIONS (STORED IN NADJ AS ! LINKED LISTS). THE PROCESS IS REPEATED UNTIL NO MORE LEVELS CAN ! BE FOUND (NLOC=0) ! NDEG (1) = NROOT BACK = .FALSE. 1 IF (NLOC > 0) THEN ! ! LOOP OVER ALL THE NPOIN OF THE PREVIOUS LEVEL, THESE ARE STORED ! IN NDEG TOGETHER WITH THE NPOIN OF PRESENT LEVEL, IN SUCH A WAY ! THAT ONE LEVEL STARTS FROM THE TOP AND THE OTHER FROM THE BOTTOM, ! THUS AVOIDING OVERLAPING. WHICH IS WHICH IS DETERMINE BY "BACK", ! THAT IS, IF BACK=TRUE THE OLD LIST STARTS FROM THE TOP AND RUNS ! DOWNWARDS WHILE THE CURRENT LIST STARTS AT THE BOTTOM AND RUNS ! UPWARDS AND VICEVERSE. ! DO IL = 1, NLOC IF (.NOT.BACK) IP = NDEG (IL) IF (BACK) IP = NDEG (NPOIN + 1 - IL) NAD = NADJ (IP) ! ! ALL NPOIN CONNECTED WITH NODE IP OF LEVEL K=LSD AND WITH LEVEL ! EQUAL ZERO SO FAR ARE ASSIGNED LEVEL LSD+1 AND STORED IN NDEG ! 2 IF (NAD > 0) THEN JP = NADJ (NAD) IF (LEV (JP) <= 0) THEN NLOCN = NLOCN + 1 LEV (JP) = LSD+1 IF (BACK) NDEG (NLOCN) = JP IF (.NOT.BACK) NDEG (NPOIN + 1 - NLOCN) = JP END IF NAD = NADJ (NAD+1) GO TO 2 END IF END DO ! ! NLOC IS THE NUMBER OF NPOIN IN THE PREVIOUS LEVEL, NLOCN ! IS THE NUMBER OF NPOIN IN THE CURRENT LEVEL AND MLW IS THE ! GLOBAL MAXIMIMUM OF NPOIN IN ANY LEVEL ! NLOCL = NLOC NLOC = NLOCN NLOCN = 0 IF (NLOC > MLW) MLW = NLOC LSD = LSD+1 BACK = .NOT.BACK GO TO 1 END IF ! ! THE NPOIN OF THE LAST LEVEL ARE RETURN AS OUTPUT OF THE SUBROUTINE ! SO IF THEY ARE UPSIDE DOWN THEY ARE COPIED PROPERLY ! NLOC = NLOCL LSD = LSD-1 IF (.NOT.BACK) THEN NHALF = NPOIN / 2 DO I = 1, NHALF NN = NDEG (I) NDEG (I) = NDEG (NPOIN + 1 - I) NDEG (NPOIN + 1 - I) = NN END DO END IF END SUBROUTINE LEVELP ! !----------------------------------------------------------------------- SUBROUTINE RESEQP (mpoi, mpoi2, NADJ, NEWNN, NACT, NODA, NPOIN, & NS, NPASS, MINMAX) !----------------------------------------------------------------------- ! ! SUBPROGRAM RESEQ1-RESEQUENCE NPOIN FOR MINIMUM FRONTWIDTH ! !----------------------------------------------------------------------- DIMENSION NADJ (mpoi2), NEWNN (mpoi), NACT (mpoi), NODA (mpoi), & NS (mpoi) ! ! KING'S SCHEME ! LARGE = 5**5 DO J = 1, NPOIN NEWNN (J) = 0 NODA (J) = 0 END DO KOUNT = 1 ! ! LOOP OVER ALL THE SUBSTRUCTERS AND RENUMBER INSIDE EACH OF THEM ! DO IPASS = 1, NPASS ISTART = NS (IPASS) NODA (ISTART) = LARGE NEWNN (ISTART) = KOUNT KOUNT = KOUNT + 1 NAC = 0 ! ! SET AS ACTIVE NPOIN THOSE THAT ARE ADJACENT TO THE STARTING ! NODE ! NAD = NADJ (ISTART) 1 IF (NAD > 0) THEN NPJ = NADJ (NAD) IF (NODA (NPJ) == 0) THEN NAC = NAC + 1 NODA (NPJ) = NAC NACT (NAC) = NPJ END IF NAD = NADJ (NAD+1) GO TO 1 END IF ! ! RENUMBER WHILE THERE ARE ACTIVES NPOIN TO SELECT ! 2 IF (NAC > 0) THEN MINNEW = LARGE LMIN = LARGE ! ! LOOP OVER ACTIVE NPOIN ! DO IACT = 1, NAC J = NACT (IACT) NEW = - 1 MIN = LARGE ! ! COMPUTE THE INCREMENT IN ACTIVE NPOIN FOR EACH NODE J ! COMPUTE WHEN THIS NODE WAS FIRST ACTIVE BY CHECKING FOR RENUMBERE ! NEIGHBOURS WITH LOWEST NUMBERS ! NAD = NADJ (J) 3 IF (NAD > 0) THEN N = NADJ (NAD) IF (NODA (N) == 0) NEW = NEW + 1 IF (NEWNN (N) /= 0) THEN IF (NEWNN (N) < MIN) MIN = NEWNN (N) END IF NAD = NADJ (NAD+1) GO TO 3 END IF ! ! SELECT NODE WITH SMALLEST INCREMENT IN ACTIVE NPOIN ! IN THE CASE OF A TIE, SELECT NODE WHICH HAS BEEN LONGEST ACTIVE ! IF (NEW <= MINNEW) THEN IF ( (NEW /= MINNEW) .OR. (MIN < LMIN) ) THEN MINNEW = NEW LMIN = MIN NEXT = J END IF END IF END DO ! ! RENUMBER NODE "NEXT" ! NEWNN (NEXT) = KOUNT KOUNT = KOUNT + 1 ! ! DEACTIVATE "NEXT" ! NPOS = NODA (NEXT) ILAST = NACT (NAC) NODA (ILAST) = NPOS NACT (NPOS) = ILAST NAC = NAC - 1 ! ! ACTIVATE THOSE NPOIN THAT ARE ADJACENT TO "NEXT" ! IF (MINNEW /= - 1) THEN NAD = IABS (NADJ (NEXT) ) 4 IF (NAD > 0) THEN N = NADJ (NAD) IF (NODA (N) == 0) THEN NAC = NAC + 1 NODA (N) = NAC NACT (NAC) = N END IF NAD = NADJ (NAD+1) GO TO 4 END IF END IF GO TO 2 END IF END DO END SUBROUTINE RESEQP ! ! ********************************************************************** ! ! ** subroutine dcheck ! ! ********************************************************************** ! SUBROUTINE dcheck (index, msize, string) CHARACTER string * 4 IF (index > msize) then PRINT * , string, ' DIMENSION EXCEEDED' STOP ELSEIF (index == 0) then PRINT * , 'MESHING FAILURE HAS OCCURRED !' PRINT * , 'TRY AGAIN BY MODIFYING SOME PARAMETER.' PRINT * , 'IN CASE OF REPEATED FAILURES RE-ENTER THE' PRINT * , 'GEOMETRY DIVIDING IT INTO MORE SUBDOMAINS.' STOP END IF END SUBROUTINE dcheck SUBROUTINE 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) !********************************************************************* ! ! *** THIS SUBROUTINE READS ALL PROBLEM DATA ! !********************************************************************* USE PRECISION_MODULE IMPLICIT REAL (DP)(A - H, O - Z) DIMENSION NBCTM (MBPLN), NBCFL (MBPLN), TEMBC (MCDPT), & FLUXP (MCDPT), COEFP (MCDPT), RADIP (MCDPT), AMBIP (MCDPT), & TEMIN (MMATR), CONDY (MMATR), CAPCY (MMATR), UVELM (MMATR), & VVELM (MMATR), CDVLU (MCDPT, MMATR), CPVLU (MCDPT, MMATR), & TVALU (MCDPT, MMATR) INTEGER :: I_O ! ! *** MATERIAL PROPERTIES ! !b need input NMATR !b xxx WRITE (7,*) 'Number of Materials defaults to ', NMATR !b IF ( NMATR > 0 ) THEN READ (10, * , IOSTAT = I_O) (CONDY (IM), CAPCY (IM), IM = 1, NMATR) IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR CONDY' WRITE (7, 912) 912 FORMAT(//' MAT.',1X,'CONDUCTIVITY',1X,'CAPACITY') DO IM = 1, NMATR WRITE (7, 906) IM, CONDY (IM), CAPCY (IM) 906 FORMAT(1X,I5,3F10.3) END DO END IF ! material input ! More non-linear data !b IF (ILINR == 1) THEN !b READ (10, * ) NCDPT !b IF (NCDPT > 0) READ (10, * ) ( (CDVLU (I, IM), CPVLU (I, IM), & !b TVALU (I, IM), I = 1, NCDPT), IM = 1, NMATR) !b ENDIF ! ! *** DIRICHLET BOUNDARY CONDITIONS (FIXED TEMPERATURE) ! READ (10, * , IOSTAT = I_O) NOBCD ! Number of EBC IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR NOBCD' WRITE (7, *) 'DIRICHLET BOUNDARY CONDITIONS = ', NOBCD IF (NOBCD > 0) THEN DO i = 1, nbpln ! Initalize all sides to no EBC nbctm (i) = 0 END DO DO I = 1, NOBCD READ (10, * , IOSTAT = I_O) TEMBC (I) ! All point or side temperatures IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR TEMBC ', I END DO READ (10, * , IOSTAT = I_O) num ! The side number of EBC groups = NOBCD ?? IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR EBC num' DO i = 1, num READ (10, * , IOSTAT = I_O) numb, numc ! side, - for 1st pt or + for all IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR numb, numc', i nbctm (numb) = numc ! 1 <= |numc| <= NOBCD END DO ENDIF ! ! *** NEUMANN BOUNDARY CONDITIONS (FIXED FLUX, CONVECTION, RADIATION) ! READ (10, * , IOSTAT = I_O) NOBCN ! Number of Neumann, Robin, or radiation IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR NOBCN' WRITE (7, *) 'NEUMANN BOUNDARY CONDITIONS = ', NOBCN IF (NOBCN > 0) THEN DO i = 1, nbpln ! Initalize all sides to no flux or mixed nbcfl (i) = 0 END DO DO I = 1, NOBCN ! Read pt or side values READ (10, * , IOSTAT = I_O) FLUXP (I), COEFP (I), RADIP (I), AMBIP (I) IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR FLUXP ', I END DO READ (10, * , IOSTAT = I_O) num ! The side number of RBC groups = NOBCN ?? IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR RBC num' DO i = 1, num READ (10, * , IOSTAT = I_O) numb, numc ! side, - for 1st edge or + for all IF ( I_O > 0 ) PRINT *, 'READ ERROR FOR 1st numb, numc ', i nbcfl (numb) = numc ! 1 <= |numc| <= NOBCN END DO ENDIF ! ! *** INITIAL TEMPERATURE FIELD FOR NONLINEAR OR TRANSIENT PROBLEMS ! !b IF (ILINR == 1.OR.ITRAN == 1) THEN !b READ (10, * ) (TEMIN (I), I = 1, NMATR) !b ENDIF ! ! *** VELOCITY FIELD FOR CONVECTION PROBLEMS ! !b IF (ICONV == 1) THEN !b READ (10, * ) (UVELM (I), VVELM (I), I = 1, NMATR) !b ENDIF END SUBROUTINE INDATA SUBROUTINE 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) !********************************************************************* ! !**** THIS SUBROUTINE APPLIES ALL INPUT DATA TO THE MESH ! !********************************************************************* USE PRECISION_MODULE IMPLICIT REAL (DP)(A - H, O - Z) DIMENSION NBCTM (MBPLN), NBCFL (MBPLN), TEMBC (MCDPT), & FLUXP (MCDPT), COEFP (MCDPT), RADIP (MCDPT), AMBIP (MCDPT), & TEMIN (MMATR), UVELM (MMATR), VVELM (MMATR) DIMENSION nbpoi (mbpln), nbref (mbpln, mpnod), nbeln (mbpln), & nbelf (mbpln, mpnod), nbele (mbpln, mpnod), inumb (MPOIN), & MTYPE (MELEM), NCONC (MELEM, MNPEL), COORD (MPOIN, MDIME), & NFIXD (MBOUN), FIXED (MBOUN), NFACB (MBOUN), NELMB (MBOUN), & FLUXE (MBOUN), COEFF (MBOUN), RADIA (MBOUN), AMBIT (MBOUN), & TEMPR (MPOIN), UVELO (MPOIN), VVELO (MPOIN), itemp (mpoin, mmatr) INTEGER :: COUNT, SIDES (NNPEL) INTEGER, INTENT (OUT) :: I_BC (MPOIN) !b ! ! *** DIRICHLET BOUNDARY CONDITIONS (FIXED TEMPERATURE) ! I_BC = 0 ; COUNT = 0 !b NELEM !b NFIXB = 0 NEUMN = 0 IF (NOBCD == 0) GOTO 201 ! no EBC anywhere !b print *,'at 5404' I_ONE = 1 !b DO 101 I = 1, NBPLN ! loop over sides NUM = NBCTM (I) ! the EBC flag for side i IF (NUM == 0) GOTO 101 ! no EBC here IF (NUM < 0) THEN ! EBC at first node only NOD = NBREF (I, NBPOI (I) ) ! the first node on side i DO K = 1, NFIXB ! loop over nodes with ECB IF (NFIXD (K) == NOD) THEN ! NOD has an EBC TEMPR (NOD) = TEMBC (IABS (NUM) ) ! set NOD BC value FIXED (K) = TEMBC (IABS (NUM) ) ! set k-th BC value GOTO 101 ! finished with search ENDIF END DO NFIXB = NFIXB + 1 ! update boundary condition count !b xxx check against max NFIXD (NFIXB) = NOD TEMPR (NOD) = TEMBC (IABS (NUM) ) FIXED (NFIXB) = TEMBC (IABS (NUM) ) I_BC (NOD) = 1 !b WRITE (23, '(2I6, 1PE13.5)') NOD, I_ONE, TEMPR (NOD) !b WRITE (33, '(2I6, 1PE13.5)') NOD, I_ONE, TEMPR (NOD) !b GOTO 101 ! finished with this EBC ENDIF DO 102 J = 1, NBPOI (I) ! loop on all pts on side i NOD = NBREF (I, J) ! j-th node number DO K = 1, NFIXB ! loop over nodes with ECB IF (NFIXD (K) == NOD) GOTO 102 ! single NOD value set, skip END DO NFIXB = NFIXB + 1 ! update boundary condition count !b xxx check against max NFIXD (NFIXB) = NOD TEMPR (NOD) = TEMBC (NUM) FIXED (NFIXB) = TEMBC (NUM) I_BC (NOD) = 1 !b WRITE (23, '(2I6, 1PE13.5)') NOD, I_ONE, TEMPR (NOD) !b WRITE (33, '(2I6, 1PE13.5)') NOD, I_ONE, TEMPR (NOD) !b 102 END DO 101 END DO PRINT *,'NOTE: NODE FLAG & BC_VALUE SAVED TO msh_bc_values.dat' !b print *,'at 5434' ! ! *** NEUMANN BOUNDARY CONDITIONS (FIXED FLUX, CONVECTION, RADIATION) ! !b print *,'at 5793 NOBCN, NBPLN =', NOBCN, NBPLN 201 IF (NOBCN == 0) GOTO 202 ! No NBC or RBC I_TWO = 2 ; SIDES = 0 !b type 2 element NNPF = 2 ; IF ( NNPEL == 6 ) NNPF = 3 !b nodes per face DO 105 I = 1, NBPLN ! loop over sides !b print *,'5516 side, # pts on side, NNPEL ', i, nbpoi (i), NNPEL !b !b print *,'all points ', (nbref (i, j), j = 1, nbpoi (i) ) !b NUM = NBCFL (I) ! the NBC flag for side i IF (NUM == 0) GOTO 105 ! no NBC !b need logic for NUM < 0 xxx else error here for first el bc DO J = 1, NBELN (I) ! loop over elements on side NEUMN = NEUMN + 1 !b xxx check against max NELMB (NEUMN) = NBELE (I, J) ! element number NFACB (NEUMN) = NBELF (I, J) ! face on this side FLUXE (NEUMN) = FLUXP (NUM) COEFF (NEUMN) = COEFP (NUM) RADIA (NEUMN) = RADIP (NUM) AMBIT (NEUMN) = AMBIP (NUM) j1 = (J-1) * (NNPF - 1) + 1 ; j2 = j1 + NNPF - 1 !b !b print *,'local elem ', j !b !b print *,'j1, j2 ', j1, j2 !b !b print *,'local elem nodes ', (nbref (i, jj) , jj = j1, j2) !b SIDES ( 1:NNPF ) = nbref (i, j1:j2 ) !b COUNT = COUNT + 1 !b WRITE (24, '(12I6)') COUNT, I_TWO, SIDES !b WRITE (34, '(12I6)') I_TWO, SIDES !b END DO 105 END DO !b print *,'at 5454' PRINT *,'NOTE: FLUX ELEMENT DATA SAVED TO msh_type_flux.dat' ! ! *** INITIAL TEMPERATURE FIELD FOR NONLINEAR OR TRANSIENT PROBLEMS ! 202 IF ( (ILINR == 1 .OR. ITRAN == 1) .and. nstep .eq. 0) THEN !b print *,'at 5459' DO I = 1, NPOIN TEMPR (I) = 0.0D0 DO j = 1, nmatr ITEMP (I, j) = 0 END DO inumb (I) = 0 END DO DO I = 1, NELEM !b print *,'at 5468 NNPEL', NNPEL MAT = MTYPE (I) DO J = 1, NNPEL !b print *,'at 5471' NOD = NCONC (I, J) IF (inumb (nod) .eq. 0) then !b print *,'at 5474' inumb (nod) = 1 itemp (nod, 1) = mat TEMPR (NOD) = TEMIN (MAT) ELSE !b print *,'at 5479' DO k = 1, inumb (nod) IF (itemp (nod, k) .eq.mat) goto 120 END DO inumb (nod) = inumb (nod) + 1 !b xxx check against max itemp (nod, inumb (nod) ) = mat TEMPR (NOD) = TEMPR (NOD) + TEMIN (MAT) 120 CONTINUE ENDIF END DO END DO DO I = 1, NPOIN !b TEMPR (I) = TEMPR (I) / DFLOAT (inumb (I) ) TEMPR (I) = TEMPR (I) / DBLE (inumb (I) ) !b END DO ! ! *** IMPOSE FIXED BOUNDARY CONDITIONS ON INITIAL TEMPERATURE FIELD ! DO I = 1, NFIXB TEMPR (NFIXD (I) ) = FIXED (I) END DO ENDIF ! INITIAL TEMPERATURE !b print *,'at 5501' ! ! *** VELOCITY FIELD FOR CONVECTION PROBLEMS ! IF (ICONV == 1) THEN !b print *,'at 5506 NPOIN, nmatr, NELEM ', NPOIN, nmatr, NELEM DO I = 1, NPOIN UVELO (I) = 0.0D0 VVELO (I) = 0.0D0 !b print *,'at 5510' DO j = 1, nmatr ITEMP (I, j) = 0 END DO !b print *,'at 5514'; if ( i > 29 ) stop '5515' !b inumb (I) = 0 END DO ! NPOIN !b print *,'at 5517 NELEM, NNPEL ', NELEM, NNPEL DO I = 1, NELEM MAT = MTYPE (I) ! Material type !b print *,'at 5520' DO J = 1, NNPEL NOD = NCONC (I, J) ! node in connectivity !b print *,'at 5523' if ( nod <= 0 ) stop '5523' IF (inumb (nod) .eq. 0) then inumb (nod) = 1 itemp (nod, 1) = mat uvelo (NOD) = uvelm (MAT) vvelo (NOD) = vvelm (MAT) ELSE !b print *,'at 5530 nod, inumb (nod) ', nod, inumb (nod) DO k = 1, inumb (nod) IF (itemp (nod, k) .eq. mat) goto 123 END DO inumb (nod) = inumb (nod) + 1 !b xxx check against max itemp (nod, inumb (nod) ) = mat uvelo (NOD) = uvelo (NOD) + uvelm (MAT) vvelo (NOD) = vvelo (NOD) + vvelm (MAT) 123 CONTINUE ENDIF END DO ! NNPEL !b print *,'at 5545' END DO ! NELEM DO I = 1, NPOIN !b print *,'at 5549 ,inumb (I) ' ,inumb (I) if ( inumb (I) <= 0 ) stop '5550' UVELO (I) = UVELO (I) / DBLE (inumb (I) ) VVELO (I) = VVELO (I) / DBLE (inumb (I) ) END DO ENDIF ! CONVECTION PROBLEMS !b print *,'at 5553' END SUBROUTINE TOMESH SUBROUTINE output_mesh (mdime, melem, mnpel, mpoin, ntot, ntr, nnpe, & coord, mtype, nconc, I_BC, NOBCN) USE PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: mdime, melem, mnpel, mpoin, & ntot, ntr, nnpe, mtype (melem), & nconc (melem, mnpel), I_BC (mpoin) INTEGER, INTENT (IN) :: NOBCN REAL (DP), INTENT (IN) :: coord (mpoin, mdime) INTEGER, PARAMETER :: I_ONE = 1 INTEGER :: j, k, MAX_TYPE ! *** ouput mesh ! 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, '(2i6, 2(1PE13.5))') j, I_BC(j), coord (j, 1:2) !b WRITE (31, '( i6, 2(1PE13.5))') I_BC(j), coord (j, 1:2) !b END DO PRINT *, 'NOTE: NODE FLAG & COORDINATES SAVED TO msh_bc_xyz.dat' ! check for number of element types MAX_TYPE = maxval (mtype(1:ntr)) !b IF ( MAX_TYPE > 1 .OR. NOBCN > 0 ) THEN !b 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)') j, mtype (j), (nconc (j, k), k = 1, nnpe)!b WRITE(32,'(12I6)') mtype (j), (nconc (j, k), k = 1, nnpe)!b ELSE WRITE(22,'(12I6)') j, mtype (j), nconc (j, 1), nconc (j, 3), & nconc (j, 5), nconc (j, 2), nconc (j, 4), nconc (j, 6) !b WRITE(32,'(12I6)') mtype (j), nconc (j, 1), nconc (j, 3), & nconc (j, 5), nconc (j, 2), nconc (j, 4), nconc (j, 6) !b END IF END DO ELSE ! only one type, omit it DO j = 1, ntr ! elem mat connectivity ! WRITE (17, * ) j, (nconc (j, k), k = 1, nnpe) IF ( nnpe <= 4 ) then !b WRITE (22, '(12I6)') j, (nconc (j, k), k = 1, nnpe)!b WRITE (32, '(12I6)') I_ONE, (nconc (j, k), k = 1, nnpe)!b ELSE WRITE (22, '(12I6)') j, nconc (j, 1), nconc (j, 3), & !b nconc (j, 5), nconc (j, 2), nconc (j, 4), nconc (j, 6) !b WRITE (32, '(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 END IF ! types PRINT *, 'NOTE: CONNECTIVITY SAVED TO msh_typ_nodes.dat' END SUBROUTINE output_mesh SUBROUTINE output_boundary (mbpln, mpnod, nbpln, nbele, nbelf, nbeln, & nbpoi, nbref) IMPLICIT NONE INTEGER, INTENT (IN) :: mbpln, mpnod, nbpln, & nbeln (mbpln), nbpoi (mbpln), nbref (mbpln, mpnod), & nbele (mbpln, mpnod), nbelf (mbpln, mpnod) INTEGER :: i, j, k ! *** 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 END SUBROUTINE output_boundary ! SUBROUTINE CONTRL (ITRAN, ILINR, IAXSY, IERST, ICONV, IPETR, & ! NITER, nitup, nitdn, TTIME, STIME, DTIME, dtmax, ALPHA, RELAX, & ! TOLER, PCENT, elmin, elmax, factr) !!********************************************************************* !! !! *** THIS SUBROUTINE READS ALL PROBLEM DATA !! !!********************************************************************* ! USE PRECISION_MODULE ! IMPLICIT REAL (DP)(A - H, O - Z) ! CHARACTER (len=72) :: title !b ! ! PCENT = 1.0D0 ! NITER = 1 !! !! *** READ AND WRITE TITLE !! !print *, 'at 5547' ! READ (10, 920) TITLE ; WRITE (7, 920) TITLE ! 920 FORMAT(12A6) !! !! *** READ AND WRITE CONTROL PARAMETERS !! ! !b READ (10, * ) ITRAN, ILINR, IAXSY, IERST, ICONV, IPETR ! WRITE (7, 901) ITRAN, ILINR, IAXSY, IERST, ICONV, IPETR ! 901 FORMAT (/, 'ITRAN =',I4, 'ILINR =',I4, 'IAXSY =',I4, & ! & 'IERST =',I4, 'ICONV =',I4, 'IPETR =',I4) ! !! Transient data !!b IF (ITRAN == 1) READ (10, * ) TTIME, STIME, DTIME, ALPHA !!b IF (ITRAN == 1) WRITE (7, 902) TTIME, STIME, DTIME, ALPHA !!b902 FORMAT(//8H TTIME =,F5.3,3X,8H STIME =,F5.3,3X, & !!b &8H DTIME =,F5.3,3X,8H ALPHA =,F5.3) ! !! Non-linear data !!b IF (ILINR == 1) READ (10, * ) NITER, TOLER, RELAX !!b IF (ILINR == 1) WRITE (7, 903) NITER, TOLER, RELAX !!b903 FORMAT(//8H NITER =,I4,4X,8H TOLER =,F5.3,3X, & !!b &8H RELAX =,F5.3) !!b IF (ILINR == 1.and.itran.eq.1) READ (10, * ) NITUP, NITDN, & !!b DTMAX, factr !!b IF (ILINR == 1.and.itran.eq.1) WRITE (7, 905) NITUP, NITDN, & !!b DTMAX, factr !!b905 FORMAT(//8H NITUP =,I4,4X,8H NITDN =,I4,4X,8H DTMAX =,F5.3,3x,& !!b &8H FACTR =,F5.3) ! !! Error estimate data !!b IF (ierst.eq.1) read (10, * ) pcent, elmin, elmax !!b IF (ierst.eq.1) WRITE (7, 904) pcent, elmin, elmax !!b904 FORMAT(/8H PCENT =,F5.3,3X,8H ELMIN =,F5.3,3X,3X,8H ELMAX =,F5.3) ! end SUBROUTINE CONTRL