PageRenderTime 33ms CodeModel.GetById 9ms app.highlight 11ms RepoModel.GetById 1ms app.codeStats 1ms

/ClassicalField/TOOLS/MeshCreation/03_CreateField/cheart_prep3D.f

http://github.com/xyan075/examples
FORTRAN Legacy | 1030 lines | 719 code | 66 blank | 245 comment | 46 complexity | 8f6084f999e8344ea6b5e78d1a11f833 MD5 | raw file
   1!    CMHeart - Preparatory routines
   2!   =======================================
   3!     For use in generating general function spaces for nodal lagrange interpolations
   4!
   5!     Written by David A. Nordsletten (C) 2008
   6!     DPhil Student, Comlab, University of Oxford 2006-2008
   7!     PhD Student, Bioengineering Institute, University of Auckland 2005-2006
   8!
   9
  10      MODULE VARSKINDS
  11!       Module VARSKINDS - Written by David A. Nordsletten (C) 2008
  12!       DPhil Student, Comlab, University of Oxford 2006-2008
  13!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
  14!
  15!       This module defines the flags for single precision (SPR)
  16!       and double precision (DPR)
  17!
  18!          i.e.  1_DPR is a double precision variable 15 digits in length
  19!                1_SPR is a single precision variable 6 digits in length 
  20!
  21        INTEGER, PARAMETER :: SPR=SELECTED_REAL_KIND(6,15)
  22        INTEGER, PARAMETER :: DPR=SELECTED_REAL_KIND(15,307)
  23      END MODULE VARSKINDS
  24
  25
  26      MODULE VARS_PREP
  27!       Module VARS_PREP - Written by David A. Nordsletten (C) 2008
  28!       DPhil Student, Comlab, University of Oxford 2006-2008
  29!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
  30!
  31!       This module defines the common variables for the preparatory routines
  32!
  33        USE VARSKINDS
  34        INCLUDE 'cheart_prepvars.h'
  35        TYPE (ARRAY_MESH) Fl					! Input Spatial Tessellation
  36        TYPE (ARRAY_MESH), pointer :: Ml(:)			! Output Function Spaces
  37        TYPE (ARRAY_PROBLEM_BASE) E, E2				! Initial and final Basis sets
  38        DOUBLE PRECISION GLOBAL_TOL				! Global tolerance parameter
  39        INTEGER, pointer :: nL(:,:)				! Work array, (T)^-1 map
  40        INTEGER LC						! Common index for nL
  41        CHARACTER*60 MAPFILE, COFILE, iBASE, fBASE, FROOT, BNDRYFILE	! Input file names
  42        CHARACTER*90 NIMZ,EXTN					! Character variales
  43      END MODULE VARS_PREP
  44
  45
  46      PROGRAM CMHeart_preparatory
  47!       Program CMHeart_preparatory - Written by David A. Nordsletten (C) 2008
  48!       DPhil Student, Comlab, University of Oxford 2006-2008
  49!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
  50!
  51!       This program executes the preparatory routines for CMHeart
  52!
  53       USE VARS_PREP
  54       IMPLICIT NONE
  55       INTEGER I,J,K
  56       GLOBAL_TOL = 0.00000001
  57       EXTN = ''	! Set extension to basis root file
  58       CALL COMM_LINE_IO	! Obtain Command line inputs
  59       CALL INPUT_CORE		! Read command line input files and build basis
  60       CALL MAKE_NEIGH		! Construct neighboring array for spatial tessellation
  61       DO I = 1,E2%n_B
  62         Ml(I)%Lb = 0
  63         IF (INDEX(E2%B(I)%CL,'M') == 1) K = I
  64         CALL MAKE_FUNCTION_SPACE( Ml(I), E2%B(I) )	! Construct function space
  65         CALL PRINTZ(3,I)	! Output completion of function space
  66       END DO
  67       IF (Fl%Lb > 0) THEN
  68          Ml(K)%ID = K; CALL RECALCULATE_BOUNDARY( Ml( K ), E2, Fl, E )
  69       END IF
  70       CALL OUTPUT_FILE		! Output updated function spaces
  71       CALL PRINTZ(100,0)	! Output completion of routines
  72      END
  73
  74
  75      SUBROUTINE RECALCULATE_BOUNDARY( NW, Enw, OD, Eod )
  76!       Subprogram RECALCULATE_BOUNDARY - Written by David A. Nordsletten (C) 2008
  77!       DPhil Student, Comlab, University of Oxford 2006-2008
  78!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
  79!
  80        USE VARS_PREP
  81        TYPE (ARRAY_MESH) OD, NW
  82        TYPE (ARRAY_PROBLEM_BASE) Enw, Eod
  83        INTEGER I,J,K,m,n,n_l,P(3),S,n_o,LS(1000)
  84        DOUBLE PRECISION XInw(3), XIod(3), Ax, Bx
  85        NW%Lb = OD%Lb; ALLOCATE( NW%B(NW%Lb, 2 + Enw%FNODES) )
  86        NW%B(1:NW%Lb,1) = OD%B(1:OD%Lb,1); NW%B(1:NW%Lb,2 + Enw%FNODES) = OD%B(1:OD%Lb,2 + Eod%FNODES)
  87        DO I = 1,OD%Lb
  88          Ax = 0; XIod(1:3) = 0
  89          DO J = 1,Eod%B(OD%ID)%n
  90            DO K = 2, 1 + Eod%FNODES
  91             IF (OD%T(OD%B(I,1),J) == OD%B(I,K)) THEN
  92               XIod(1:3) = XIod(1:3) + Eod%B(OD%ID)%XI(J,1:3); Ax = Ax + 1; EXIT
  93             END IF
  94            END DO
  95          END DO
  96          XIod(1:3) = XIod(1:3) / Ax
  97          PRINT *, XIod(1:3)
  98          DO J = 1,Enw%FACES
  99            IF (OD%NBT(OD%B(I,1),J) == 0) THEN
 100              Ax = 0; XInw(1:3) = 0; n = 1
 101              DO K = 1,Enw%B(NW%ID)%n
 102                Bx = 0
 103                DO m = 1,Enw%n_pts_f
 104                   Bx = Bx + Enw%B(NW%ID)%I%Y_f(K,m,J)	! Calculate approximate integral on facet
 105                END DO
 106                IF (abs(Bx) > GLOBAL_TOL) THEN
 107                   n = n + 1; NW%B(I,n) = NW%T(NW%B(I,1),K)
 108                   XInw(1:3) = XInw(1:3) + Enw%B(NW%ID)%XI(K,1:3); Ax = Ax + 1;
 109                END IF
 110              END DO
 111              XInw(1:3) = XInw(1:3)/Ax - XIod(1:3)
 112              PRINT *, XInw(1:3)
 113              CALL P_NORM(XInw,3,2,Bx)
 114              IF (Bx < GLOBAL_TOL) EXIT
 115            END IF
 116          END DO
 117        END DO
 118      END SUBROUTINE RECALCULATE_BOUNDARY
 119
 120
 121      SUBROUTINE MAKE_NEIGH
 122!       Subprogram MAKE_NEIGH - Written by David A. Nordsletten (C) 2008
 123!       DPhil Student, Comlab, University of Oxford 2006-2008
 124!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 125!
 126!       This program builds a map of neighboring elements for each element of the
 127!       spatial tessellation.  This is accomplished by first constructing the inverse
 128!       map of Fl%T, identifying all elements for each node (nL).  The set of facet 
 129!       nodes is constructed for each facet (RW).  The neighbor of the facet J, is 
 130!       a unique element at the intersection of all I, nL(RW(I,J),:).
 131!
 132        USE VARS_PREP
 133        INTEGER I,J,K,m,n,n_l,P(3),S,n_o,LS(1000)
 134        INTEGER IRW(E%B(Fl%ID)%n,E%FACES)
 135        INTEGER RW(E%B(Fl%ID)%n,E%FACES),L_RW(E%FACES)
 136        LC = 100
 137        ALLOCATE(Fl%NBT(Fl%Lt,E%FACES),nL(Fl%Lx,LC))
 138        Fl%NBT(1:Fl%Lt,1:E%FACES) = 0		! Initialize neighboring array as zeros
 139        nL(1:Fl%Lx,LC) = 0			! nL(:,LC), no. of elements, initialized as zero
 140        DO I = 1,Fl%Lt
 141           DO J = 1,E%B(Fl%ID)%n
 142              nL(Fl%T(I,J),LC) = nL(Fl%T(I,J),LC) + 1	! increment no. elements for basis Fl%T(I,J)
 143              nL(Fl%T(I,J),nL(Fl%T(I,J),LC)) = I	! add element I to list for basis Fl%T(I,J)
 144           END DO
 145        END DO
 146        CALL BUILD_IRW(IRW,L_RW)	! Generate IRW
 147        DO I = 1,Fl%Lt
 148          CALL BUILD_RW(I,RW,IRW,L_RW)	! Generate RW and L_RW
 149          DO J = 1,E%FACES		! Loop over facets
 150             n_l = nL(RW(1,J),LC)	! set initial list length
 151             LS(1:n_l) = nL(RW(1,J),1:n_l)	! set initial list
 152             DO K = 1,n_l
 153               IF (LS(K) == I) LS(K) = 0	! delete current element from list
 154             END DO
 155             DO K = 2,L_RW(J)	! Loop over all basis on facet J
 156               DO m = 1,n_l	! Loop over all components from initial element list
 157                 IF (LS(m) > 0) THEN
 158                   S = 0
 159                   DO n = 1,nL(RW(K,J),LC)	! Loop over all elements associated with basis RW(K,J)
 160                     IF (nL(RW(K,J),n) == LS(m)) THEN ! If LS(m) is shared, keep it
 161                       S = 1; EXIT
 162                     END IF
 163                   END DO
 164                   IF (S == 0) LS(m) = 0	! LS(m) is not shared, delete it
 165                 END IF
 166               END DO
 167             END DO
 168             n = 0
 169             DO K = 1,n_l	! Loop over initial list
 170               IF (LS(K) > 0) THEN
 171                 n = n + 1; Fl%NBT(I,J) = LS(K)		! If component is at the intersection, set as neighbor
 172               END IF
 173             END DO
 174             IF (n > 1) CALL ERRORZ(6)		! If the number of interesection components is non-unique, quit
 175          END DO
 176        END DO
 177      END SUBROUTINE MAKE_NEIGH
 178
 179
 180      SUBROUTINE BUILD_IRW(IRW,L_RW)
 181!       Subprogram BUILD_IRW - Written by David A. Nordsletten (C) 2008
 182!       DPhil Student, Comlab, University of Oxford 2006-2008
 183!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 184!
 185!       This program builds arrays IRW and L_RW for spatial basis.  IRW (I, J) stores the
 186!       the Ith basis index on the Jth facet.  L_RW(J) is the number of basis on the Jth facet.
 187!       These arrays are computed by identifying facets based on those with possible values
 188!       on a facet.  In nodal lagrange interpolants, non-zeros indicate the basis is on
 189!       a given facet.
 190!
 191        USE VARS_PREP
 192        INTEGER I,J,K,ELEM,IRW(E%B(Fl%ID)%n,E%FACES),L_RW(E%FACES),FACE
 193        DOUBLE PRECISION Ax
 194        DO FACE = 1,E%FACES	! Loop over facets
 195           IRW(1:E%B(Fl%ID)%n,FACE) = 0		! Initialize RW as zeros
 196           L_RW(FACE) = 0			! Initialize L_RW as zeros
 197           DO I = 1,E%B(Fl%ID)%n		! Loop over all basis components of ELEM
 198              Ax = 0
 199              DO K = 1,E%n_pts_f
 200                 Ax = Ax + E%B(Fl%ID)%I%Y_f(I,K,FACE)	! Calculate approximate integral on facet
 201              END DO
 202              IF (abs(Ax) > GLOBAL_TOL) THEN
 203                  L_RW(FACE) = L_RW(FACE) + 1
 204                  IRW(L_RW(FACE),FACE) = I
 205              END IF
 206           END DO
 207        END DO
 208      END SUBROUTINE BUILD_IRW
 209
 210
 211      SUBROUTINE BUILD_RW(ELEM,RW,IRW,L_RW)
 212!       Subprogram BUILD_RW - Written by David A. Nordsletten (C) 2008
 213!       DPhil Student, Comlab, University of Oxford 2006-2008
 214!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 215!
 216!       This program builds arrays RW and L_RW for element ELEM.  RW (I, J) stores the
 217!       the Ith basis on the Jth facet.  L_RW(J) is the number of basis on the Jth facet.
 218!       These arrays are computed by identifying facets based on those with possible values
 219!       on a facet.  In nodal lagrange interpolants, non-zeros indicate the basis is on
 220!       a given facet.
 221!
 222        USE VARS_PREP
 223        INTEGER I,J,K,ELEM,RW(E%B(Fl%ID)%n,E%FACES),L_RW(E%FACES),FACE
 224        INTEGER IRW(E%B(Fl%ID)%n,E%FACES),ne
 225        DOUBLE PRECISION Ax
 226        DO FACE = 1,E%FACES	! Loop over facets
 227           DO I = 1,L_RW(FACE)	! Loop over number of basis for facet FACE
 228              RW(I,FACE) = Fl%T(ELEM,IRW(I,FACE))	! Set value of RW
 229           END DO
 230           K = 1; ne = 100000		! Reorder so that smallest element set is first
 231           DO I = 1,L_RW(FACE)
 232              IF (nL(RW(I,FACE),LC) < ne) THEN
 233                 ne = nL(RW(I,FACE),LC); K = I
 234              END IF
 235           END DO
 236           I = RW(1,FACE); J = RW(K,FACE)
 237           RW(1,FACE) = J; RW(K,FACE) = I
 238        END DO
 239      END SUBROUTINE BUILD_RW
 240
 241
 242      SUBROUTINE MAKE_FUNCTION_SPACE(AA,EE)
 243!       Subprogram MAKE_FUNCTION_SPACE - Written by David A. Nordsletten (C) 2008
 244!       DPhil Student, Comlab, University of Oxford 2006-2008
 245!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 246!
 247!       This program generates the function space (AA) on the spatial map (Fl) based
 248!       on the basis (EE).  This is process is straightforward for discontinuous fields.
 249!       In the case of continuous fields, facet basis will be shared.  Thus, as the 
 250!       space is computed for an element ELEM, a list of all previously defined elements 
 251!       ( I < ELEM) which share a corner, edge, or facet with ELEM is composed.  It is then
 252!       verified that a basis is unique before added to the coefficient list.
 253!
 254        USE VARS_PREP
 255        TYPE(ARRAY_MESH) AA
 256        TYPE(ARRAY_BASE) EE
 257        INTEGER I,J,K,m,n,P,S,ELEM,ELE_ARRAY(1000),L_EA
 258        DOUBLE PRECISION P1(3),P2(3),dL,WGT(8),TOL,ASH1,ASH2
 259        ALLOCATE(AA%T(Fl%Lt,EE%n),AA%X(EE%n*Fl%Lt,3))
 260        AA%Lt = Fl%Lt; AA%Lx = 0		! Initialize map size and coefficient size
 261        DO ELEM = 1,Fl%Lt		! Loop over elements of Fl
 262           IF (EE%DISCONT == 0) THEN	! If field is continuous
 263           L_EA = 0		! Set list length of previous elements to zero
 264           DO I = 1, E%B(Fl%ID)%n	! Loop over basis of ELEM
 265              DO J = 1,nL(Fl%T(ELEM,I),LC)	! Loop over elements sharing basis Fl%T(ELEM,I)
 266                 IF (nL(Fl%T(ELEM,I),J) < ELEM) THEN 	! If the element is previously encountered
 267                   m = 0
 268                   DO K = 1,L_EA
 269                    IF (nL(Fl%T(ELEM,I),J) == ELE_ARRAY(K)) THEN	! See if it exists on list ELE_ARRAY
 270                      m = 1; EXIT
 271                    END IF
 272                   END DO
 273                   IF (m == 0) THEN		! If not, add to list and update counter L_EA
 274                     L_EA = L_EA + 1;  ELE_ARRAY(L_EA) = nL(Fl%T(ELEM,I),J)
 275                   END IF
 276                 END IF
 277              END DO
 278           END DO
 279           END IF
 280           DO I = 1,EE%n	! Loop over the basis for Function Space
 281              CALL GET_SPATIAL_COORDINATE(E,ELEM,EE%XI(I,1:3),P1)	! Get a spatial coordinate
 282              m = 0
 283              IF (EE%DISCONT == 0) THEN		! If the field is continuous
 284                DO J = 1,L_EA		! Loop over previous element list
 285                   DO K = 1,EE%n	! Loop over the basis for Function Space
 286                     CALL P_NORM(P1(1:3) - AA%X(AA%T(ELE_ARRAY(J),K),1:3),3,2,dL)
 287                     IF (dL < GLOBAL_TOL) THEN		! Is the basis unique?
 288                       m = AA%T(ELE_ARRAY(J),K); EXIT		! If no, set m to its redundancy
 289                     END IF
 290                   END DO
 291                   IF (m > 0) EXIT	! If a redundancy has been found, exit
 292                END DO
 293              END IF
 294              DO J = 1,I-1	! Loop over previously defined basis within ELEM
 295                 CALL P_NORM(P1(1:3) - AA%X(AA%T(ELEM,J),1:3),3,2,dL)
 296                 IF (dL < GLOBAL_TOL) THEN 	! Is the basis unique?
 297                   m = AA%T(ELEM,J); EXIT		! If no, set m to its redundancy
 298                 END IF
 299              END DO
 300              IF (m == 0) THEN		! If no redundancy exists, add basis to coefficient map
 301                AA%Lx = AA%Lx + 1; AA%X(AA%Lx,1:3) = P1(1:3)
 302                m = AA%Lx
 303              END IF
 304              AA%T(ELEM,I) = m		! Add basis to function space
 305           END DO
 306        END DO
 307      END SUBROUTINE MAKE_FUNCTION_SPACE
 308
 309
 310      SUBROUTINE INPUT_CORE
 311!       Subprogram INPUT_CORE - Written by David A. Nordsletten (C) 2008
 312!       DPhil Student, Comlab, University of Oxford 2006-2008
 313!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 314!
 315!       This program reads in all input files.
 316!
 317        USE VARS_PREP
 318        INTEGER I,J,K,m,n
 319        DOUBLE PRECISION Ax,Bx(1000),Cx
 320        CHARACTER*60 NAMZ
 321        CHARACTER*1 ID
 322        NIMZ = TRIM(iBASE); CALL ADD_EXTENSION(NIMZ)	! Add root pathway to initial basis file name
 323        CALL LOAD_BASE_SET(E)		! Read initial basis
 324        NAMZ = TRIM(COFILE)
 325        OPEN(UNIT=1,FILE=NAMZ,STATUS='old',action='read')	! Read initial coefficient file
 326          READ(1,*) Ax
 327          Fl%Lx = INT(Ax); ALLOCATE(Fl%X(Fl%Lx,3))
 328          CALL READ_MESH_CORED(Fl%X,Fl%Lx,3)
 329        CLOSE(1)
 330        NAMZ = TRIM(MAPFILE); Fl%ID = 0
 331        OPEN(UNIT=1,FILE=NAMZ,STATUS='old',action='read')	! Find basis call letter
 332          READ(1,*,ERR=45) Ax, ID
 333 45       DO I = 1,E%n_B
 334             IF (INDEX(ID,E%B(I)%CL) == 1) Fl%ID = I	! If call letter found, set ID = I
 335          END DO
 336          IF (Fl%ID == 0) THEN		! If the call letter is not found, search basis for M
 337             DO I = 1,E%n_B
 338                IF (INDEX(E%B(I)%CL,'M') == 1) Fl%ID = I
 339             END DO
 340             IF (Fl%ID == 0) CALL ERRORZ(2)	! If no basis is found, exit
 341          END IF
 342        CLOSE(1)
 343        OPEN(UNIT=1,FILE=NAMZ,STATUS='old',action='read')	! Read initial spatial map
 344          READ(1,*) Ax
 345          Fl%Lt = INT(Ax); ALLOCATE(Fl%T(Fl%Lt,E%B(Fl%ID)%n))
 346          CALL READ_MESH_COREI(Fl%T,Fl%Lt,E%B(Fl%ID)%n)
 347        CLOSE(1)
 348        NIMZ = TRIM(fBASE); CALL ADD_EXTENSION(NIMZ)	! Add root pathway to final basis file name
 349        Fl%Lb = 0
 350        IF (BNDRYFILE .NE. '') THEN
 351           NAMZ = TRIM(BNDRYFILE)
 352           OPEN(UNIT=1,FILE=NAMZ,STATUS='old',action='read')	! Find basis call letter
 353             READ(1,*) Ax
 354             Fl%Lb = INT(Ax); ALLOCATE(Fl%B(Fl%Lb, E%FNODES + 2  ))
 355             CALL READ_MESH_COREI(Fl%B,Fl%Lb, E%FNODES + 2  )
 356           CLOSE(1)
 357        END IF
 358        CALL LOAD_BASE_SET(E2)		! Read initial basis
 359        ALLOCATE(Ml(E2%n_B))		! Allocate no. of Function spaces
 360        CALL PRINTZ(2,0)
 361      END SUBROUTINE INPUT_CORE
 362
 363
 364      SUBROUTINE READ_MESH_COREI(A,n,m)
 365!       Subprogram READ_MESH_COREI - Written by David A. Nordsletten (C) 2008
 366!       DPhil Student, Comlab, University of Oxford 2006-2008
 367!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 368!
 369!       This program reads in integer arrays for file 1.
 370!
 371        INTEGER I,n,m,A(n,m)
 372        DOUBLE PRECISION Bx(m)
 373        DO I = 1,n
 374          READ(1,*,END=30) Bx(1:m); A(I,1:m) = INT(Bx(1:m))
 375        END DO
 376        RETURN
 377 30     CALL ERRORZ(4)
 378      END SUBROUTINE READ_MESH_COREI
 379
 380
 381      SUBROUTINE READ_MESH_CORED(A,n,m)
 382!       Subprogram READ_MESH_CORED - Written by David A. Nordsletten (C) 2008
 383!       DPhil Student, Comlab, University of Oxford 2006-2008
 384!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 385!
 386!       This program reads in double arrays for file 1.
 387!
 388        INTEGER I,n,m
 389        DOUBLE PRECISION A(n,m)
 390        DO I = 1,n
 391          READ(1,*,END=35) A(I,1:m)
 392        END DO
 393        RETURN
 394 35     CALL ERRORZ(4)
 395      END SUBROUTINE READ_MESH_CORED
 396
 397
 398      SUBROUTINE OUTPUT_FILE
 399!       Subprogram OUTPUT_FILE - Written by David A. Nordsletten (C) 2008
 400!       DPhil Student, Comlab, University of Oxford 2006-2008
 401!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 402!
 403!       This program reads in all output files.
 404!
 405        USE VARS_PREP
 406        INTEGER I,J,K,m,n
 407        DOUBLE PRECISION Ax,Bx(100),Cx
 408        CHARACTER*30 NAMZ
 409        CALL PRINTZ(4,0)
 410        NAMZ = TRIM(FROOT)//'_FE.C'
 411        OPEN(UNIT=1,FILE=NAMZ,STATUS='unknown')		! Output updated coefficient file
 412          WRITE(1,*) Ml(1:E2%n_B)%Lx
 413          DO J = 1,E2%n_B
 414            DO I = 1,Ml(J)%Lx
 415              WRITE(1,*) Ml(J)%X(I,1:3)
 416            END DO
 417          END DO
 418        CLOSE(1)
 419        NAMZ = TRIM(FROOT)//'_FE.M'
 420        OPEN(UNIT=1,FILE=NAMZ,STATUS='unknown')		! Output updated map file
 421          WRITE(1,*) Ml(1:E2%n_B)%Lt, E2%B(1:E2%n_B)%CL
 422          DO J = 1,E2%n_B
 423            DO I = 1,Ml(J)%Lt
 424              WRITE(1,*) Ml(J)%T(I,1:E2%B(J)%n)
 425            END DO
 426          END DO
 427        CLOSE(1)
 428        NAMZ = TRIM(FROOT)//'_FE.NBT'
 429        OPEN(UNIT=1,FILE=NAMZ,STATUS='unknown')		! Output updated neighboring file
 430          WRITE(1,*) Fl%Lt
 431          DO I = 1,Fl%Lt
 432            WRITE(1,*) Fl%NBT(I,1:E%FACES)
 433          END DO
 434        CLOSE(1)
 435        DO K = 1,E2%n_B
 436          IF (Ml(K)%Lb > 0) THEN
 437            NAMZ = TRIM(FROOT)//'_FE.B'
 438            OPEN(UNIT=1,FILE=NAMZ,STATUS='unknown')		! Output updated neighboring file
 439              WRITE(1,*) Ml(K)%Lb
 440              DO I = 1,Ml(K)%Lb
 441                WRITE(1,*) Ml(K)%B(I,1:E2%FNODES+2)
 442              END DO
 443            CLOSE(1)
 444          END IF
 445        END DO
 446      END SUBROUTINE OUTPUT_FILE
 447
 448
 449      SUBROUTINE LOAD_BASE_SET(AA)
 450!       Subprogram LOAD_BASE_SET - Written by David A. Nordsletten (C) 2008
 451!       DPhil Student, Comlab, University of Oxford 2006-2008
 452!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 453!
 454!       This program reads in all input basis for base AA.
 455!
 456        USE VARS_PREP
 457        INTEGER I,STP,D,FLD
 458        CHARACTER*60 IN_CHAR
 459        TYPE(ARRAY_PROBLEM_BASE) AA
 460        NIMZ = TRIM(NIMZ); AA%n_B = 0; AA%HEXA = 0; AA%DM = 3
 461        OPEN(UNIT=1,FILE=NIMZ,STATUS='old',action='read')	! Read base file for initial parameters
 462          DO WHILE (0 < 1)
 463              READ(1,*,END=50) IN_CHAR
 464              IF (INDEX(IN_CHAR,'no_fields!') == 1)         READ(1,*) AA%n_B
 465              IF (INDEX(IN_CHAR,'no_gauss!') == 1)          READ(1,*) AA%n_pts
 466              IF (INDEX(IN_CHAR,'volume!') == 1)            READ(1,*) AA%VL
 467              IF (INDEX(IN_CHAR,'no_gauss_f!') == 1)        READ(1,*) AA%n_pts_f
 468              IF (INDEX(IN_CHAR,'no_ele_faces!') == 1)      READ(1,*) AA%FACES
 469              IF (INDEX(IN_CHAR,'no_ele_nodes_f!') == 1)    READ(1,*) AA%FNODES
 470              IF (INDEX(IN_CHAR,'hexa_basis!') == 1)        AA%HEXA = 1
 471              IF (INDEX(IN_CHAR,'domain_dimension!') == 1)  READ(1,*) AA%DM
 472              IF (INDEX(IN_CHAR,'TRI_BASIS!') == 1)  AA%TRI_BASIS = 1
 473              IF (INDEX(IN_CHAR,'TET_BASIS!') == 1)  AA%TET_BASIS = 1
 474              IF (INDEX(IN_CHAR,'QUAD_BASIS!') == 1) AA%QUAD_BASIS = 1
 475              IF (INDEX(IN_CHAR,'HEX_BASIS!') == 1)  AA%HEX_BASIS = 1
 476            END DO
 477 50     CLOSE(1)
 478        IF (AA%n_B == 0) CALL ERRORZ(3); ALLOCATE(AA%B(AA%n_B))		! Allocate base size
 479        OPEN(UNIT=1,FILE=NIMZ,STATUS='old',action='read')		! Read call letters
 480          DO WHILE (0 < 1)
 481              READ(1,*,END=55) IN_CHAR
 482              IF (INDEX(IN_CHAR,'call_letters!') == 1)       READ(1,*) AA%B(1:AA%n_B)%CL
 483            END DO
 484 55     CLOSE(1)
 485        AA%B(1:AA%n_B)%DISCONT = 0
 486        AA%n_ptsl = AA%n_pts
 487        AA%n_pts_fl = AA%n_pts_f
 488        IF (AA%HEXA == 0) THEN			! If tensor input is indicated
 489           D = 3
 490           ALLOCATE(AA%gpt_f(AA%n_pts_f,D,AA%FACES),AA%gw_f(AA%n_pts_f))
 491           ALLOCATE(AA%gpt(AA%n_pts,D),AA%gw(AA%n_pts))
 492        ELSE					! If full input is indicated
 493           D = 1
 494           AA%n_pts_f = AA%n_pts**REAL(AA%DM-1); AA%n_pts = AA%n_pts**REAL(AA%DM)
 495           ALLOCATE(AA%gpt(AA%n_pts,3),AA%gw(AA%n_pts))
 496           ALLOCATE(AA%gpt_f(AA%n_pts_f,3,AA%FACES),AA%gw_f(AA%n_pts_f))
 497        END IF
 498        OPEN(UNIT=1,FILE=NIMZ,STATUS='old',action='read')	! Read quadrature rule
 499          DO WHILE (0 < 1)
 500              READ(1,*,END=65) IN_CHAR
 501              IF (INDEX(IN_CHAR,'surface_area!') == 1)  READ(1,*) AA%VL_f(1:AA%FACES)
 502              IF (INDEX(IN_CHAR,'gauss_points!') == 1) THEN
 503                 DO I = 1,AA%n_ptsl
 504                   READ(1,*) AA%gpt(I,1:D)
 505                 END DO
 506              ELSE IF (INDEX(IN_CHAR,'gauss_weights!') == 1) THEN
 507                 DO I = 1,AA%n_ptsl
 508                   READ(1,*) AA%gw(I)
 509                 END DO
 510              ELSE IF (INDEX(IN_CHAR,'gauss_points_f!') == 1) THEN
 511                DO I = 1,AA%n_pts_fl
 512                  READ(1,*) AA%gpt_f(I,1:3,1)
 513                END DO
 514              ELSE IF (INDEX(IN_CHAR,'gauss_weights_f!') == 1) THEN
 515                DO I = 1,AA%n_pts_fl
 516                  READ(1,*) AA%gw_f(I)
 517                END DO
 518              END IF
 519            END DO
 520 65     CLOSE(1)
 521        IF (AA%HEXA == 1) CALL BUILD_HEXA_QUADRATURE(AA)	! Build Tensor quadrature rule
 522        CALL BUILD_FACE_GPT_ARRAY(AA)				! Build facet quadrature rule
 523        DO FLD = 1,AA%n_B		! Loop over number of basis
 524        OPEN(UNIT=1,FILE=NIMZ,STATUS='old',action='read')	! Read dimensions of basis CL
 525          DO WHILE (0 < 1)
 526              READ(1,*,END=70) IN_CHAR
 527              IF (INDEX(IN_CHAR,'no_basis_'//TRIM(AA%B(FLD)%CL)//'!') == 1) READ(1,*) AA%B(FLD)%n
 528              IF (INDEX(IN_CHAR,'dim_field_'//TRIM(AA%B(FLD)%CL)//'!') == 1) READ(1,*) AA%B(FLD)%DM
 529            END DO
 530 70     CLOSE(1)
 531        AA%B(FLD)%nl = AA%B(FLD)%n
 532        IF (AA%HEXA == 0) THEN		! If tensor input is indicated
 533           D = 3
 534           ALLOCATE(AA%B(FLD)%Q(AA%B(FLD)%n,AA%B(FLD)%n),AA%B(FLD)%P(3,AA%B(FLD)%n))
 535           ALLOCATE( AA%B(FLD)%XI(AA%B(FLD)%n,3) )
 536        ELSE				! If full input is indicated
 537           D = 1
 538           AA%B(FLD)%n  = AA%B(FLD)%n**REAL(AA%DM)
 539           ALLOCATE(AA%B(FLD)%Q(AA%B(FLD)%nl,AA%B(FLD)%nl),AA%B(FLD)%P(D,AA%B(FLD)%nl))
 540           ALLOCATE( AA%B(FLD)%B_ID(AA%B(FLD)%n,3), AA%B(FLD)%XI(AA%B(FLD)%n,3) )
 541        END IF
 542        OPEN(UNIT=1,FILE=NIMZ,STATUS='old',action='read')		! Read in Basis
 543          DO WHILE (0 < 1)
 544              READ(1,*,END=80) IN_CHAR
 545              IF (INDEX(IN_CHAR,TRIM(AA%B(FLD)%CL)//'_basis!') == 1) THEN
 546                 DO I = 1,AA%B(FLD)%nl
 547                 READ(1,*) AA%B(FLD)%Q(I,1:AA%B(FLD)%nl)
 548                 END DO
 549              ELSE IF (INDEX(IN_CHAR,TRIM(AA%B(FLD)%CL)//'_basis_discontinuous!') == 1) THEN
 550                 AA%B(FLD)%DISCONT = 1
 551              ELSE IF (INDEX(IN_CHAR,TRIM(AA%B(FLD)%CL)//'_xi_coordinates!') == 1) THEN
 552                 DO I = 1,AA%B(FLD)%nl
 553                 READ(1,*) AA%B(FLD)%XI(I,1:D)
 554                 END DO
 555              ELSE IF (INDEX(IN_CHAR,TRIM(AA%B(FLD)%CL)//'_P_basis!') == 1) THEN
 556                 DO I = 1,D
 557                 READ(1,*) AA%B(FLD)%P(I,1:AA%B(FLD)%nl)
 558                 END DO
 559              ELSE IF (INDEX(IN_CHAR,'basis_ordering_'//TRIM(AA%B(FLD)%CL)//'!') == 1) THEN
 560                DO I = 1,AA%B(FLD)%n
 561                  READ(1,*) AA%B(FLD)%B_ID(I,1:3)
 562                END DO
 563              END IF
 564            END DO
 565 80     CLOSE(1)
 566        IF (AA%HEXA == 1) CALL BUILD_HEXA_XI(AA,FLD)		! Build tensor xi points
 567        CALL LOAD_INTEGRALS(AA,FLD)				! Calculate basis integrals
 568        END DO
 569        RETURN
 570 90     CALL ERRORZ(5)
 571      END SUBROUTINE LOAD_BASE_SET
 572
 573
 574      SUBROUTINE LOAD_INTEGRALS(AA,FLD)
 575!       Subprogram LOAD_BASE_SET - Written by David A. Nordsletten (C) 2008
 576!       DPhil Student, Comlab, University of Oxford 2006-2008
 577!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 578!
 579!       This program loads integrals for base AA%B(FLD).  The basis is evaluated a given quadrature points.
 580!
 581        USE VARS_PREP
 582        TYPE(ARRAY_PROBLEM_BASE) AA
 583        INTEGER I,J,K,x,y,z,ID,ORD1(3),FLD,S
 584        DOUBLE PRECISION Ax(3),Bx,Cx(3),Dx
 585        ALLOCATE( AA%B(FLD)%I%Y(AA%B(FLD)%n,AA%n_pts) )
 586        ALLOCATE( AA%B(FLD)%I%Y_f(AA%B(FLD)%n,AA%n_pts_f,AA%FACES) )
 587        ORD1(1:3) = 0
 588        DO I = 1,AA%B(FLD)%n		! Loop over basis
 589           DO J = 1,AA%n_pts		! Loop over quadrature points
 590              CALL POLY(AA%B(FLD),AA%HEXA,I,AA%gpt(J,1:3),ORD1,AA%B(FLD)%I%Y(I,J))	! Calculate basis
 591           END DO
 592        END DO
 593        ORD1(1:3) = 0
 594        DO S = 1,AA%FACES	! Loop over facets
 595          DO I = 1,AA%B(FLD)%n		! Loop over basis
 596            DO J = 1,AA%n_pts_f			! Loop over facet quadrature points
 597              CALL POLY(AA%B(FLD),AA%HEXA,I,AA%gpt_f(J,1:3,S),ORD1,AA%B(FLD)%I%Y_f(I,J,S))	! Calculate basis
 598            END DO
 599          END DO
 600        END DO
 601      END SUBROUTINE LOAD_INTEGRALS
 602
 603
 604      SUBROUTINE BUILD_HEXA_XI(AA,FLD)
 605!       Subprogram BUILD_HEXA_XI - Written by David A. Nordsletten (C) 2008
 606!       DPhil Student, Comlab, University of Oxford 2006-2008
 607!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 608!
 609!       This routine builds the XI coordinate for the tensor basis AA%B(FLD)
 610!
 611        USE VARS_PREP
 612        TYPE(ARRAY_PROBLEM_BASE) AA
 613        INTEGER I,J,K,S,m,n,G,FLD
 614        DOUBLE PRECISION Bx(AA%B(FLD)%nl)
 615        Bx(1:AA%B(FLD)%nl) = AA%B(FLD)%XI(1:AA%B(FLD)%nl,1)
 616        DO I = 1,AA%B(FLD)%n		! Loop over basis
 617           AA%B(FLD)%XI(I,1:3) = Bx(AA%B(FLD)%B_ID(I,1:3))	! calculate xi point
 618        END DO
 619      END SUBROUTINE BUILD_HEXA_XI
 620
 621
 622      SUBROUTINE BUILD_HEXA_QUADRATURE(AA)
 623!       Subprogram BUILD_HEXA_QUADRATURE - Written by David A. Nordsletten (C) 2008
 624!       DPhil Student, Comlab, University of Oxford 2006-2008
 625!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 626!
 627!       This routine builds the volume/facet quadrature rule for the tensor basis AA
 628!
 629        USE VARS_PREP
 630        TYPE(ARRAY_PROBLEM_BASE) AA
 631        INTEGER I,J,K,S,m,n,G
 632        DOUBLE PRECISION Ax(AA%n_ptsl),Bx(AA%n_ptsl)
 633        Ax(1:AA%n_ptsl) = AA%gpt(1:AA%n_ptsl,1)
 634        Bx(1:AA%n_ptsl) = AA%gw(1:AA%n_ptsl)
 635        n = 0
 636        DO I = 1,AA%n_ptsl	! Loop over basis DIM 3
 637           IF ((AA%DM == 3).OR.(I == 1)) THEN	! If the dimension is 3 or I = 1
 638           DO J = 1,AA%n_ptsl	! Loop over basis DIM 2
 639              IF ((AA%DM >= 2).OR.(J == 1)) THEN	! If the dimension is > 2 or J = 1
 640              DO K = 1,AA%n_ptsl	! Loop over basis DIM 1
 641                 IF (AA%DM >= 1) THEN		! If the dimension is > 1 or J = 1
 642                   n = n + 1
 643                   AA%gpt(n,1:3) = 0; AA%gw(n) = 0	! Initialize volume quadrature rule, and build
 644                   IF (AA%DM >= 1) AA%gpt(n,1) = Ax(K)
 645                   IF (AA%DM >= 2) AA%gpt(n,2) = Ax(J)
 646                   IF (AA%DM >= 3) AA%gpt(n,3) = Ax(I)
 647                   IF (AA%DM == 1) AA%gw(n) = Bx(K)
 648                   IF (AA%DM == 2) AA%gw(n) = Bx(K)*Bx(J)
 649                   IF (AA%DM == 3) AA%gw(n) = Bx(K)*Bx(J)*Bx(I)
 650                 END IF
 651              END DO
 652              END IF
 653           END DO
 654           END IF
 655        END DO
 656        m = 0
 657        DO I = 1,AA%n_ptsl	! Loop over basis DIM 2
 658           IF ((AA%DM-1 == 2).OR.(I == 1)) THEN		! If the dimension is 2 or I = 1
 659           DO J = 1,AA%n_ptsl		! Loop over basis DIM 1
 660              IF ((AA%DM-1 >= 1).OR.(J == 1)) THEN	! If the dimension is > 1 or J = 1
 661                m = m + 1
 662                AA%gpt_f(m,1:3,1) = 0; AA%gw_f(m) = 0	! Initialize facet quadrature rule, and build
 663                IF (AA%DM-1 >= 1) AA%gpt_f(m,1,1) = Ax(J)
 664                IF (AA%DM-1 >= 2) AA%gpt_f(m,2,1) = Ax(I)
 665                IF (AA%DM-1 == 1) AA%gw_f(m) = Bx(J)
 666                IF (AA%DM-1 == 2) AA%gw_f(m) = Bx(J)*Bx(I)
 667              END IF
 668           END DO
 669           END IF
 670        END DO
 671        AA%VL_f(1:AA%FACES) = AA%VL**REAL(AA%DM-1)	! Initialize facet volume
 672        AA%VL = AA%VL**REAL(AA%DM)	! Initialize volume
 673      END SUBROUTINE BUILD_HEXA_QUADRATURE
 674
 675
 676      SUBROUTINE BUILD_FACE_GPT_ARRAY(AA)
 677!       Subprogram BUILD_FACE_GPT_ARRAY - Written by David A. Nordsletten (C) 2008
 678!       DPhil Student, Comlab, University of Oxford 2006-2008
 679!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 680!
 681!       This routine builds facet quadrature (must be preceeded by a call to BUILD_HEXA_QUADRATURE,
 682!       if necessary).
 683!
 684        USE VARS_PREP
 685        TYPE(ARRAY_PROBLEM_BASE) AA
 686        INTEGER I,J,K
 687        DOUBLE PRECISION Ax(3),Bx,Cx(3),Dx,ONN
 688        ONN = 1.0
 689        AA%gpt_f(1:AA%n_pts_f,AA%DM:3,1) = 0		! Initialize first quadrature rule
 690        DO I = 1,AA%FACES
 691           AA%nrm(1:3,I) = 0		! Initialize normal vectors
 692        END DO
 693        IF (AA%DM == 3) THEN		! If in R^3
 694         IF (AA%FACES == 4) THEN	! master element is a tetrahedron, build appropriate facet quadrature
 695          AA%nrm(3,1) = -ONN
 696
 697          AA%gpt_f(1:AA%n_pts_f,1,2) = AA%gpt_f(1:AA%n_pts_f,1,1)
 698          AA%gpt_f(1:AA%n_pts_f,2,2) = 0
 699          AA%gpt_f(1:AA%n_pts_f,3,2) = AA%gpt_f(1:AA%n_pts_f,2,1)
 700          AA%nrm(2,2) = -ONN
 701
 702
 703          AA%gpt_f(1:AA%n_pts_f,1,3) = 0
 704          AA%gpt_f(1:AA%n_pts_f,2,3) = AA%gpt_f(1:AA%n_pts_f,1,1)
 705          AA%gpt_f(1:AA%n_pts_f,3,3) = AA%gpt_f(1:AA%n_pts_f,2,1)
 706          AA%nrm(1,3) = -ONN
 707
 708          DO I = 1,AA%n_pts_f
 709          AA%gpt_f(I,1,4) = ONN-AA%gpt_f(I,1,1)-AA%gpt_f(I,2,1)
 710          END DO
 711          AA%gpt_f(1:AA%n_pts_f,2,4) = AA%gpt_f(1:AA%n_pts_f,1,1)
 712          AA%gpt_f(1:AA%n_pts_f,3,4) = AA%gpt_f(1:AA%n_pts_f,2,1)
 713          Ax(1) = ONN
 714          Ax(2) = 3.0
 715          Ax(2) = Ax(2)**0.5
 716          AA%nrm(1:3,4) = Ax(1)/Ax(2)
 717         ELSE IF (AA%FACES == 6) THEN	! master element is a cube, build appropriate facet quadrature
 718          AA%nrm(3,1) = -ONN
 719
 720          AA%gpt_f(1:AA%n_pts_f,1,2) = AA%gpt_f(1:AA%n_pts_f,1,1)
 721          AA%gpt_f(1:AA%n_pts_f,2,2) = 0
 722          AA%gpt_f(1:AA%n_pts_f,3,2) = AA%gpt_f(1:AA%n_pts_f,2,1)
 723          AA%nrm(2,2) = -ONN
 724
 725          AA%gpt_f(1:AA%n_pts_f,1,3) = 0
 726          AA%gpt_f(1:AA%n_pts_f,2,3) = AA%gpt_f(1:AA%n_pts_f,1,1)
 727          AA%gpt_f(1:AA%n_pts_f,3,3) = AA%gpt_f(1:AA%n_pts_f,2,1)
 728          AA%nrm(1,3) = -ONN
 729
 730          AA%gpt_f(1:AA%n_pts_f,1,4) = ONN
 731          AA%gpt_f(1:AA%n_pts_f,2,4) = AA%gpt_f(1:AA%n_pts_f,1,1)
 732          AA%gpt_f(1:AA%n_pts_f,3,4) = AA%gpt_f(1:AA%n_pts_f,2,1)
 733          AA%nrm(1,4) = ONN
 734
 735          AA%gpt_f(1:AA%n_pts_f,1,5) = AA%gpt_f(1:AA%n_pts_f,1,1)
 736          AA%gpt_f(1:AA%n_pts_f,2,5) = ONN
 737          AA%gpt_f(1:AA%n_pts_f,3,5) = AA%gpt_f(1:AA%n_pts_f,2,1)
 738          AA%nrm(2,5) = ONN
 739
 740          AA%gpt_f(1:AA%n_pts_f,1,6) = AA%gpt_f(1:AA%n_pts_f,1,1)
 741          AA%gpt_f(1:AA%n_pts_f,2,6) = AA%gpt_f(1:AA%n_pts_f,2,1)
 742          AA%gpt_f(1:AA%n_pts_f,3,6) = ONN
 743          AA%nrm(3,6) = ONN
 744         END IF
 745        ELSE IF (AA%DM == 2) THEN	! If in R^2
 746         IF (AA%FACES == 3) THEN	! master element is a triangle, build appropriate facet quadrature
 747          AA%nrm(2,1) = -ONN
 748
 749          AA%gpt_f(1:AA%n_pts_f,1,2) = 0
 750          AA%gpt_f(1:AA%n_pts_f,2,2) = AA%gpt_f(1:AA%n_pts_f,1,1)
 751          AA%gpt_f(1:AA%n_pts_f,3,2) = 0
 752          AA%nrm(1,2) = -ONN
 753
 754          DO I = 1,AA%n_pts_f
 755          AA%gpt_f(I,1,3) = ONN-AA%gpt_f(I,1,1)
 756          END DO
 757          AA%gpt_f(1:AA%n_pts_f,2,3) = AA%gpt_f(1:AA%n_pts_f,1,1)
 758          AA%gpt_f(1:AA%n_pts_f,3,3) = 0
 759          Ax(1) = ONN
 760          Ax(2) = 2.0
 761          Ax(2) = Ax(2)**0.5
 762          AA%nrm(1:2,3) = Ax(1)/Ax(2); AA%nrm(3,3) = 0
 763         ELSE IF (AA%FACES == 4) THEN	! master element is a square, build appropriate facet quadrature
 764          AA%nrm(2,1) = -ONN
 765
 766          AA%gpt_f(1:AA%n_pts_f,1,2) = AA%gpt_f(1:AA%n_pts_f,1,1)
 767          AA%gpt_f(1:AA%n_pts_f,2,2) = 0
 768          AA%gpt_f(1:AA%n_pts_f,3,2) = 0
 769          AA%nrm(1,2) = -ONN
 770
 771          AA%gpt_f(1:AA%n_pts_f,1,3) = AA%gpt_f(1:AA%n_pts_f,1,1)
 772          AA%gpt_f(1:AA%n_pts_f,2,3) = ONN
 773          AA%gpt_f(1:AA%n_pts_f,3,3) = 0
 774          AA%nrm(2,3) = ONN
 775
 776          AA%gpt_f(1:AA%n_pts_f,1,4) = ONN
 777          AA%gpt_f(1:AA%n_pts_f,2,4) = AA%gpt_f(1:AA%n_pts_f,1,1)
 778          AA%gpt_f(1:AA%n_pts_f,3,4) = 0
 779          AA%nrm(1,4) = ONN
 780         END IF
 781        END IF
 782      END SUBROUTINE BUILD_FACE_GPT_ARRAY
 783
 784
 785      SUBROUTINE GET_SPATIAL_COORDINATE(AA,ELEM,XI,Pt)
 786!       Subprogram GET_SPATIAL_COORDINATE - Written by David A. Nordsletten (C) 2008
 787!       DPhil Student, Comlab, University of Oxford 2006-2008
 788!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 789!
 790!       This routine evaluates spatial tessellation (Fl) with basis (AA), at xi point (XI)
 791!       in element Pt.
 792!
 793        USE VARS_PREP
 794        TYPE(ARRAY_PROBLEM_BASE) AA
 795        INTEGER I, J, K, ELEM, ORD1(3), ID
 796        DOUBLE PRECISION XI(3), Pt(3), P1(3), Ax
 797        Pt(1:3) = 0
 798        ORD1(1:3) = 0
 799        DO K = 1,AA%B(Fl%ID)%n	! Loop over basis
 800           CALL POLY(AA%B(Fl%ID),AA%HEXA,K,XI,ORD1,Ax)	! calculate basis K
 801           Pt(1:3) = Pt(1:3) + Fl%X(Fl%T(ELEM,K),1:3)*Ax	! Sum and add coefficient weighting
 802        END DO
 803      END SUBROUTINE GET_SPATIAL_COORDINATE
 804
 805
 806      SUBROUTINE POLY(AA,HEXA,BAS,X,ORD1,VAL)
 807!       Subprogram POLY - Written by David A. Nordsletten (C) 2008
 808!       DPhil Student, Comlab, University of Oxford 2006-2008
 809!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 810!
 811!       This routine evaluates basis (BAS) of (AA) with derivatives of order (ORD1) in each spatial
 812!       dimension.  The answer is provided in VAL
 813!
 814        USE VARS_PREP
 815        TYPE(ARRAY_BASE) AA
 816        INTEGER n,ORD1(3),I,J,P,BAS,K,HEXA
 817        DOUBLE PRECISION VAL,X(3)
 818        IF (HEXA == 1) THEN	! If tensor basis
 819           CALL POLY_HEX(AA,BAS,X,ORD1,VAL)
 820        ELSE			! If general basis
 821           CALL POLY_TET(AA,BAS,X,ORD1,VAL)
 822        END IF
 823      END SUBROUTINE POLY
 824
 825
 826      SUBROUTINE POLY_HEX(AA,BAS,X,ORD1,VAL)
 827!       Subprogram POLY_HEX - Written by David A. Nordsletten (C) 2008
 828!       DPhil Student, Comlab, University of Oxford 2006-2008
 829!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 830!
 831!       This routine evaluates basis (BAS) of (AA) with derivatives of order (ORD1) in each spatial
 832!       dimension.  The answer is provided in VAL
 833!
 834        USE VARS_PREP
 835        TYPE(ARRAY_BASE) AA
 836        INTEGER n,ORD1(3),I,J,P,BAS,K
 837        DOUBLE PRECISION VAL,Ax,Bx,Cx,Bo,X(3),Co(3)
 838        Co(1:3) = 0
 839        DO I = 1,AA%nl		! Loop over basis in a given coordinate dimension
 840           DO J = 1,AA%DM		! Loop over dimension
 841              IF (ORD1(J) <= AA%P(1,I)) THEN	! If power is great enough, that term is non-zero after diff
 842                 Ax = 1
 843                 DO K = 0,ORD1(J)-1	! Calculate the derivative of the term
 844                    Ax = Ax*(REAL(AA%P(1,I))-REAL(K))
 845                 END DO
 846                 Ax = Ax*(X(J)**(AA%P(1,I)-ORD1(J)))
 847                 Co(J) = Co(J) + AA%Q(AA%B_ID(BAS,J),I)*Ax	! add to current direction weighted by c.tensor
 848              END IF
 849           END DO
 850        END DO
 851        VAL = Co(1)
 852        DO I = 2,AA%DM		! Do tensor product
 853          VAL = VAL*Co(I)
 854        END DO
 855      END SUBROUTINE POLY_HEX
 856
 857
 858      SUBROUTINE POLY_TET(AA,BAS,X,ORD1,VAL)
 859!       Subprogram POLY_TET - Written by David A. Nordsletten (C) 2008
 860!       DPhil Student, Comlab, University of Oxford 2006-2008
 861!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 862!
 863!       This routine evaluates basis (BAS) of (AA) with derivatives of order (ORD1) in each spatial
 864!       dimension.  The answer is provided in VAL
 865!
 866        USE VARS_PREP
 867        TYPE(ARRAY_BASE) AA
 868        INTEGER n,ORD1(3),I,J,P,BAS,K
 869        DOUBLE PRECISION VAL,Ax,Bx,Cx,Bo,X(3),Co(3)
 870        VAL = 0
 871        DO I = 1,AA%n			! Loop over basis 
 872           Co(1:3) = 0; Bx = 1_DPR
 873           DO J = 1,AA%DM		! Loop over dimension
 874             IF (ORD1(J) <= AA%P(J,I)) THEN	! If power is great enough, that term is non-zero after diff
 875              Ax = 1
 876              DO K = 0,ORD1(J)-1	! Calculate the derivative of the term
 877                 Ax = Ax*(REAL(AA%P(J,I))-REAL(K))
 878              END DO
 879              Ax = Ax*(X(J)**(AA%P(J,I)-ORD1(J)))
 880              Co(J) = Ax
 881             END IF
 882             Bx = Bx*Co(J)	! calculate product of term
 883           END DO
 884           VAL = VAL + AA%Q(BAS,I)*Bx	! add to current value weighted by c.tensor
 885        END DO
 886      END SUBROUTINE POLY_TET
 887
 888
 889      SUBROUTINE P_NORM(A,n,p,dL)
 890!       Subprogram P_NORM - Written by David A. Nordsletten (C) 2008
 891!       DPhil Student, Comlab, University of Oxford 2006-2008
 892!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 893!
 894!       This routine evaluates dL = l_p (A(1:n)).  If 0, its the l_infty norm.
 895!
 896        USE VARSKINDS
 897        INTEGER n,I,p
 898        DOUBLE PRECISION A(n), dL
 899        dL = 0
 900        IF (p.NE.0) THEN	! Not zero, calculate norm
 901           DO I = 1,n
 902           dL = dL + (abs(A(I)))**REAL(p)
 903           END DO
 904           dL = dL**(1_DPR / REAL(p))
 905        ELSE			! Calculate l_infty norm
 906           DO I = 1,n
 907              dL = MAX(dL,abs(A(I)))
 908           END DO
 909        END IF
 910      END SUBROUTINE P_NORM
 911
 912
 913      SUBROUTINE ADD_EXTENSION(NAMZ)
 914!       Subprogram ADD_EXTENSION - Written by David A. Nordsletten (C) 2008
 915!       DPhil Student, Comlab, University of Oxford 2006-2008
 916!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 917!
 918!       This routine adds pathway to basis files to NAMZ
 919!
 920        USE VARS_PREP
 921        CHARACTER*90 NAMZ
 922        IF (TRIM(EXTN).NE.'') THEN
 923            NAMZ = TRIM(EXTN)//TRIM(NAMZ)
 924        END IF
 925      END SUBROUTINE ADD_EXTENSION
 926
 927
 928      SUBROUTINE COMM_LINE_IO
 929!       Subprogram COMM_LINE_IO - Written by David A. Nordsletten (C) 2008
 930!       DPhil Student, Comlab, University of Oxford 2006-2008
 931!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 932!
 933!       This routine gets the command line arguements for the preparatory routines
 934!
 935        USE VARS_PREP
 936        INTEGER I,J,K,CNT, n
 937        CHARACTER*60 NAMZ2
 938        CNT = COMMAND_ARGUMENT_COUNT(); I = 0; n = 0; BNDRYFILE = ''
 939        DO WHILE (I < CNT)
 940           I = I + 1
 941           CALL GET_COMMAND_ARGUMENT(I,NAMZ2,J,K)
 942           IF (INDEX(NAMZ2,'-iM') == 1) THEN		! If the initial map flag
 943               I = I + 1; n = n + 1
 944               CALL GET_COMMAND_ARGUMENT(I,MAPFILE,J,K)
 945               MAPFILE = TRIM(MAPFILE)
 946           ELSE IF (INDEX(NAMZ2,'-iC') == 1) THEN	! If the initial coefficient flag
 947               I = I + 1; n = n + 1
 948               CALL GET_COMMAND_ARGUMENT(I,COFILE,J,K)
 949               COFILE = TRIM(COFILE)
 950           ELSE IF (INDEX(NAMZ2,'-iGAMMA') == 1) THEN
 951               I = I + 1
 952               CALL GET_COMMAND_ARGUMENT(I,BNDRYFILE,J,K)
 953               BNDRYFILE = TRIM(BNDRYFILE)
 954           ELSE IF (INDEX(NAMZ2,'-iB') == 1) THEN	! If the initial basis flag
 955                I = I + 1; n = n + 1
 956               CALL GET_COMMAND_ARGUMENT(I,iBASE,J,K)
 957               iBASE = TRIM(iBASE)
 958           ELSE IF (INDEX(NAMZ2,'-fB') == 1) THEN	! If the final basis flag
 959                I = I + 1; n = n + 1
 960               CALL GET_COMMAND_ARGUMENT(I,fBASE,J,K)
 961               fBASE = TRIM(fBASE)
 962           ELSE IF (INDEX(NAMZ2,'-fO') == 1) THEN	! If the output flag
 963                I = I + 1; n = n + 1
 964               CALL GET_COMMAND_ARGUMENT(I,FROOT,J,K)
 965               FROOT = TRIM(FROOT)
 966           END IF
 967        END DO
 968        CALL PRINTZ(1,0)
 969        IF (n < 5) CALL ERRORZ(1)
 970      END SUBROUTINE COMM_LINE_IO
 971
 972
 973      SUBROUTINE PRINTZ(I,K)
 974!       Subprogram PRINTZ - Written by David A. Nordsletten (C) 2008
 975!       DPhil Student, Comlab, University of Oxford 2006-2008
 976!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
 977!
 978!       This routine prints messages and status for the preparatory routines
 979!
 980        USE VARS_PREP
 981        INTEGER I,K
 982        IF (I == 1) THEN
 983          PRINT *, ' ==>>  COMMAND LINE INPUTS COLLECTED  <<== '
 984          PRINT *, ' ==>>  MAP FILE  <<== '
 985          PRINT *, ' ----  '//TRIM(MAPFILE)
 986          PRINT *, ' ==>>  COEFFICIENT FILE  <<== '
 987          PRINT *, ' ----  '//TRIM(COFILE)
 988          PRINT *, ' ==>>  INITIAL BASIS  <<== '
 989          PRINT *, ' ----  '//TRIM(iBASE)
 990          PRINT *, ' ==>>  FINAL BASIS  <<== '
 991          PRINT *, ' ----  '//TRIM(fBASE)
 992          PRINT *, ' ==>>  ROOT OUTPUT FILE NAME  <<== '
 993          PRINT *, ' ----  '//TRIM(FROOT)
 994        END IF
 995        IF (I == 2) PRINT *, ' ==>>  STATUS :: INPUT FILES READ  <<=='
 996        IF (I == 3) THEN
 997          PRINT *, ' ==>>  STATUS :: SPACE '//TRIM(E2%B(K)%CL)//' COMPLETE  <<== '
 998        END IF
 999        IF (I == 4) THEN
1000          PRINT *, ' ==>>  OUTPUT FILES  <<== '
1001          PRINT *, ' ==>>  UPDATED MAP FILE  <<== '
1002          PRINT *, ' ----  '//TRIM(FROOT)//'_FE.T'
1003          PRINT *, ' ==>>  UPDATED COEFFICIENT FILE  <<== '
1004          PRINT *, ' ----  '//TRIM(FROOT)//'_FE.X'
1005          PRINT *, ' ==>>  NEIGHBORING FILE  <<== '
1006          PRINT *, ' ----  '//TRIM(FROOT)//'_FE.NBT'
1007        END IF
1008
1009        IF (I == 100) PRINT *, ' ==>>  STATUS :: DONE  <<=='
1010      END SUBROUTINE PRINTZ
1011
1012
1013      SUBROUTINE ERRORZ(I)
1014!       Subprogram ERRORZ - Written by David A. Nordsletten (C) 2008
1015!       DPhil Student, Comlab, University of Oxford 2006-2008
1016!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
1017!
1018!       This routine prints error messages and stops execution for the preparatory routines
1019!
1020        USE VARS_PREP
1021        INTEGER I,K
1022        IF (I == 1) PRINT *, ' ==>>  FORCE QUIT :: INCOMPLETE COMMAND LINE INPUT  <<== '
1023        IF (I == 2) PRINT *, ' ==>>  FORCE QUIT :: NO CALL LETTER or M BASIS  <<== '
1024        IF (I == 3) PRINT *, ' ==>>  FORCE QUIT :: NO FIELDS IN BASIS FILE <<== '
1025        IF (I == 4) PRINT *, ' ==>>  FORCE QUIT :: INCOMPLETE MAP or COEFFICIENT FILE(S) <<== '
1026        IF (I == 5) PRINT *, ' ==>>  FORCE QUIT :: PROBLEM ENCOUNTERED READING BASIS FILE <<== '
1027        IF (I == 6) PRINT *, ' ==>>  FORCE QUIT :: MULTIPLE ELEMENTS POINTING TO A SINGLE FACET <<== '
1028        PRINT *, ' ==>> QUITING <<== '
1029        STOP
1030      END SUBROUTINE ERRORZ