/ClassicalField/TOOLS/MeshCreation/03_CreateField/cheart_prep3D.f
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