#### /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
324        NAMZ = TRIM(COFILE)
327          Fl%Lx = INT(Ax); ALLOCATE(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
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)
345          Fl%Lt = INT(Ax); ALLOCATE(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
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
359        ALLOCATE(Ml(E2%n_B))		! Allocate no. of Function spaces
360        CALL PRINTZ(2,0)
361      END SUBROUTINE INPUT_CORE
362
363
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)
379
380
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
392        END DO
393        RETURN
394 35     CALL ERRORZ(4)
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
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
462          DO WHILE (0 < 1)
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
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
480          DO WHILE (0 < 1)
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
499          DO WHILE (0 < 1)
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
505                 END DO
506              ELSE IF (INDEX(IN_CHAR,'gauss_weights!') == 1) THEN
507                 DO I = 1,AA%n_ptsl
509                 END DO
510              ELSE IF (INDEX(IN_CHAR,'gauss_points_f!') == 1) THEN
511                DO I = 1,AA%n_pts_fl
513                END DO
514              ELSE IF (INDEX(IN_CHAR,'gauss_weights_f!') == 1) THEN
515                DO I = 1,AA%n_pts_fl
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
525          DO WHILE (0 < 1)
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
543          DO WHILE (0 < 1)
545              IF (INDEX(IN_CHAR,TRIM(AA%B(FLD)%CL)//'_basis!') == 1) THEN
546                 DO 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
554                 END DO
555              ELSE IF (INDEX(IN_CHAR,TRIM(AA%B(FLD)%CL)//'_P_basis!') == 1) THEN
556                 DO I = 1,D
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
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)
572
573
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
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
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
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
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
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
```