PageRenderTime 21ms CodeModel.GetById 16ms app.highlight 2ms RepoModel.GetById 1ms app.codeStats 0ms

/FluidMechanics/TOOLS/MeshCreation/02_CreateCMInput/cheart_visbasis3D.f

http://github.com/xyan075/examples
FORTRAN Legacy | 362 lines | 295 code | 29 blank | 38 comment | 18 complexity | 61a7294af502c6696daa3a613be748e6 MD5 | raw file
  1      SUBROUTINE LOAD_BASE_SET(AA)
  2!       Subprogram LOAD_BASE_SET - Written by David A. Nordsletten (C) 2008
  3!       DPhil Student, Comlab, University of Oxford 2006-2008
  4!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
  5!
  6!       This program reads in all input basis for base AA.
  7!
  8        USE CMGUI_VARS
  9        INTEGER I,STP,D,VFLD
 10        CHARACTER*60 IN_CHAR
 11        TYPE(ARRAY_PROBLEM_BASE) AA
 12        NIMZ = TRIM(NIMZ); AA%n_B = 0; AA%HEXA = 0; AA%DM = 3
 13        OPEN(UNIT=1,FILE=NIMZ,STATUS='old',action='read')	! Read base file for initial parameters
 14          DO WHILE (0 < 1)
 15              READ(1,*,END=50) IN_CHAR
 16              IF (INDEX(IN_CHAR,'no_fields!') == 1)         READ(1,*) AA%n_B
 17              IF (INDEX(IN_CHAR,'no_gauss!') == 1)          READ(1,*) AA%n_pts
 18              IF (INDEX(IN_CHAR,'volume!') == 1)            READ(1,*) AA%VL
 19              IF (INDEX(IN_CHAR,'no_gauss_f!') == 1)        READ(1,*) AA%n_pts_f
 20              IF (INDEX(IN_CHAR,'no_ele_faces!') == 1)      READ(1,*) AA%FACES
 21              IF (INDEX(IN_CHAR,'no_ele_nodes_f!') == 1)    READ(1,*) AA%FNODES
 22              IF (INDEX(IN_CHAR,'hexa_basis!') == 1)        AA%HEXA = 1
 23              IF (INDEX(IN_CHAR,'domain_dimension!') == 1)  READ(1,*) AA%DM
 24              IF (INDEX(IN_CHAR,'TRI_BASIS!') == 1)  AA%TRI_BASIS = 1
 25              IF (INDEX(IN_CHAR,'TET_BASIS!') == 1)  AA%TET_BASIS = 1
 26              IF (INDEX(IN_CHAR,'QUAD_BASIS!') == 1) AA%QUAD_BASIS = 1
 27              IF (INDEX(IN_CHAR,'HEX_BASIS!') == 1)  AA%HEX_BASIS = 1
 28            END DO
 29 50     CLOSE(1)
 30        IF (AA%n_B == 0) CALL ERRORZ(3); ALLOCATE(AA%B(AA%n_B))		! Allocate base size
 31        OPEN(UNIT=1,FILE=NIMZ,STATUS='old',action='read')		! Read call letters
 32          DO WHILE (0 < 1)
 33              READ(1,*,END=55) IN_CHAR
 34              IF (INDEX(IN_CHAR,'call_letters!') == 1)       READ(1,*) AA%B(1:AA%n_B)%CL
 35          END DO
 36 55     CLOSE(1)
 37        AA%B(1:AA%n_B)%DISCONT = 0
 38        AA%n_ptsl = AA%n_pts
 39        AA%n_pts_fl = AA%n_pts_f
 40        IF (AA%HEXA == 0) THEN			! If tensor input is indicated
 41           D = 3
 42           ALLOCATE(AA%gpt_f(AA%n_pts_f,D,AA%FACES),AA%gw_f(AA%n_pts_f))
 43           ALLOCATE(AA%gpt(AA%n_pts,D),AA%gw(AA%n_pts))
 44        ELSE					! If full input is indicated
 45           D = 1
 46           AA%n_pts_f = AA%n_pts**REAL(AA%DM-1); AA%n_pts = AA%n_pts**REAL(AA%DM)
 47           ALLOCATE(AA%gpt(AA%n_pts,3),AA%gw(AA%n_pts))
 48           ALLOCATE(AA%gpt_f(AA%n_pts_f,3,AA%FACES),AA%gw_f(AA%n_pts_f))
 49        END IF
 50        OPEN(UNIT=1,FILE=NIMZ,STATUS='old',action='read')	! Read quadrature rule
 51          DO WHILE (0 < 1)
 52              READ(1,*,END=65) IN_CHAR
 53              IF (INDEX(IN_CHAR,'surface_area!') == 1)  READ(1,*) AA%VL_f(1:AA%FACES)
 54              IF (INDEX(IN_CHAR,'gauss_points!') == 1) THEN
 55                 DO I = 1,AA%n_ptsl
 56                   READ(1,*) AA%gpt(I,1:D)
 57                 END DO
 58              ELSE IF (INDEX(IN_CHAR,'gauss_weights!') == 1) THEN
 59                 DO I = 1,AA%n_ptsl
 60                   READ(1,*) AA%gw(I)
 61                 END DO
 62              ELSE IF (INDEX(IN_CHAR,'gauss_points_f!') == 1) THEN
 63                DO I = 1,AA%n_pts_fl
 64                  READ(1,*) AA%gpt_f(I,1:3,1)
 65                END DO
 66              ELSE IF (INDEX(IN_CHAR,'gauss_weights_f!') == 1) THEN
 67                DO I = 1,AA%n_pts_fl
 68                  READ(1,*) AA%gw_f(I)
 69                END DO
 70              END IF
 71            END DO
 72 65     CLOSE(1)
 73        IF (AA%HEXA == 1) CALL BUILD_HEXA_QUADRATURE(AA)	! Build Tensor quadrature rule
 74        CALL BUILD_FACE_GPT_ARRAY(AA)				! Build facet quadrature rule
 75        DO VFLD = 1,AA%n_B		! Loop over number of basis
 76        IF (INDEX(AA%B(VFLD)%CL,'M') >= 1) AA%SPL = VFLD
 77        IF (INDEX(AA%B(VFLD)%CL,'V') >= 1) AA%VEL = VFLD
 78        IF (INDEX(AA%B(VFLD)%CL,'P') >= 1) AA%PRS = VFLD
 79        OPEN(UNIT=1,FILE=NIMZ,STATUS='old',action='read')	! Read dimensions of basis CL
 80          DO WHILE (0 < 1)
 81              READ(1,*,END=70) IN_CHAR
 82              IF (INDEX(IN_CHAR,'no_basis_'//TRIM(AA%B(VFLD)%CL)//'!') == 1) READ(1,*) AA%B(VFLD)%n
 83              IF (INDEX(IN_CHAR,'dim_field_'//TRIM(AA%B(VFLD)%CL)//'!') == 1) READ(1,*) AA%B(VFLD)%DM
 84            END DO
 85 70     CLOSE(1)
 86        AA%B(VFLD)%nl = AA%B(VFLD)%n
 87        IF (AA%HEXA == 0) THEN		! If tensor input is indicated
 88           D = 3
 89           ALLOCATE(AA%B(VFLD)%Q(AA%B(VFLD)%n,AA%B(VFLD)%n),AA%B(VFLD)%P(3,AA%B(VFLD)%n))
 90           ALLOCATE( AA%B(VFLD)%XI(AA%B(VFLD)%n,3) )
 91        ELSE				! If full input is indicated
 92           D = 1
 93           AA%B(VFLD)%n  = AA%B(VFLD)%n**REAL(AA%DM)
 94           ALLOCATE(AA%B(VFLD)%Q(AA%B(VFLD)%nl,AA%B(VFLD)%nl),AA%B(VFLD)%P(D,AA%B(VFLD)%nl))
 95           ALLOCATE( AA%B(VFLD)%B_ID(AA%B(VFLD)%n,3), AA%B(VFLD)%XI(AA%B(VFLD)%n,3) )
 96        END IF
 97        OPEN(UNIT=1,FILE=NIMZ,STATUS='old',action='read')		! Read in Basis
 98          DO WHILE (0 < 1)
 99              READ(1,*,END=80) IN_CHAR
100              IF (INDEX(IN_CHAR,TRIM(AA%B(VFLD)%CL)//'_basis!') == 1) THEN
101                 DO I = 1,AA%B(VFLD)%nl
102                 READ(1,*) AA%B(VFLD)%Q(I,1:AA%B(VFLD)%nl)
103                 END DO
104              ELSE IF (INDEX(IN_CHAR,TRIM(AA%B(VFLD)%CL)//'_basis_discontinuous!') == 1) THEN
105                 AA%B(VFLD)%DISCONT = 1
106              ELSE IF (INDEX(IN_CHAR,TRIM(AA%B(VFLD)%CL)//'_xi_coordinates!') == 1) THEN
107                 DO I = 1,AA%B(VFLD)%nl
108                 READ(1,*) AA%B(VFLD)%XI(I,1:D)
109                 END DO
110              ELSE IF (INDEX(IN_CHAR,TRIM(AA%B(VFLD)%CL)//'_P_basis!') == 1) THEN
111                 DO I = 1,D
112                 READ(1,*) AA%B(VFLD)%P(I,1:AA%B(VFLD)%nl)
113                 END DO
114              ELSE IF (INDEX(IN_CHAR,'basis_ordering_'//TRIM(AA%B(VFLD)%CL)//'!') == 1) THEN
115                DO I = 1,AA%B(VFLD)%n
116                  READ(1,*) AA%B(VFLD)%B_ID(I,1:3)
117                END DO
118              END IF
119            END DO
120 80     CLOSE(1)
121        IF (AA%HEXA == 1) CALL BUILD_HEXA_XI(AA,VFLD)		! Build tensor xi points
122        CALL LOAD_INTEGRALS(AA,VFLD)				! Calculate basis integrals
123        PRINT *, ' ==== >>>   BASIS FIELD  <<< ==== '
124        PRINT *, '   CALL LETTER     :: ', AA%B(VFLD)%CL
125        PRINT *, '   FIELD DIMENSION :: ', AA%B(VFLD)%DM
126        PRINT *, '   DIM  (S ( Tm )) :: ', AA%B(VFLD)%n
127
128        END DO
129        RETURN
130 90     CALL ERRORZ(5)
131      END SUBROUTINE LOAD_BASE_SET
132
133
134      SUBROUTINE LOAD_INTEGRALS(AA,VFLD)
135!       Subprogram LOAD_INTEGRALS - Written by David A. Nordsletten (C) 2008
136!       DPhil Student, Comlab, University of Oxford 2006-2008
137!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
138!
139!       This program loads integrals for base AA%B(VFLD).  The basis is evaluated a given quadrature points.
140!
141        USE CMGUI_VARS
142        TYPE(ARRAY_PROBLEM_BASE) AA
143        INTEGER I,J,K,x,y,z,ID,ORD1(3),VFLD,S
144        DOUBLE PRECISION Ax(3),Bx,Cx(3),Dx
145        ALLOCATE( AA%B(VFLD)%I%Y(AA%B(VFLD)%n,AA%n_pts) )
146        ALLOCATE( AA%B(VFLD)%I%Y_f(AA%B(VFLD)%n,AA%n_pts_f,AA%FACES) )
147        ALLOCATE( AA%B(VFLD)%I%dY(3,AA%B(VFLD)%n,AA%n_pts) )
148        ALLOCATE( AA%B(VFLD)%I%TdY(3,AA%B(VFLD)%n,AA%n_pts) )
149        ALLOCATE( AA%B(VFLD)%I%dY_f(3,AA%B(VFLD)%n,AA%n_pts_f,AA%FACES) )
150        DO I = 1,AA%B(VFLD)%n		! Loop over basis
151           ORD1(1:3) = 0
152           DO J = 1,AA%n_pts		! Loop over quadrature points
153              CALL POLY(AA,VFLD,I,AA%gpt(J,1:3),ORD1,AA%B(VFLD)%I%Y(I,J))	! Calculate basis
154           END DO
155           DO K = 1,3
156             ORD1(1:3) = 0; ORD1(K) = 1
157             DO J = 1,AA%n_pts		! Loop over quadrature points
158                CALL POLY(AA,VFLD,I,AA%gpt(J,1:3),ORD1,AA%B(VFLD)%I%dY(K,I,J))	! Calculate basis
159             END DO
160           END DO
161        END DO
162        DO S = 1,AA%FACES	! Loop over facets
163          DO I = 1,AA%B(VFLD)%n		! Loop over basis
164            ORD1(1:3) = 0
165            DO J = 1,AA%n_pts_f			! Loop over facet quadrature points
166              CALL POLY(AA,VFLD,I,AA%gpt_f(J,1:3,S),ORD1,AA%B(VFLD)%I%Y_f(I,J,S))	! Calculate basis
167            END DO
168            DO K = 1,3
169              ORD1(1:3) = 0; ORD1(K) = 1
170              DO J = 1,AA%n_pts_f		! Loop over quadrature points
171                CALL POLY(AA,VFLD,I,AA%gpt_f(J,1:3,S),ORD1,AA%B(VFLD)%I%dY_f(K,I,J,S))	! Calculate basis
172              END DO
173            END DO
174          END DO
175        END DO
176      END SUBROUTINE LOAD_INTEGRALS
177
178
179      SUBROUTINE BUILD_HEXA_XI(AA,VFLD)
180!       Subprogram BUILD_HEXA_XI - Written by David A. Nordsletten (C) 2008
181!       DPhil Student, Comlab, University of Oxford 2006-2008
182!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
183!
184!       This routine builds the XI coordinate for the tensor basis AA%B(VFLD)
185!
186        USE CMGUI_VARS
187        TYPE(ARRAY_PROBLEM_BASE) AA
188        INTEGER I,J,K,S,m,n,G,VFLD
189        DOUBLE PRECISION Bx(AA%B(VFLD)%nl)
190        Bx(1:AA%B(VFLD)%nl) = AA%B(VFLD)%XI(1:AA%B(VFLD)%nl,1)
191        DO I = 1,AA%B(VFLD)%n		! Loop over basis
192           AA%B(VFLD)%XI(I,1:3) = Bx(AA%B(VFLD)%B_ID(I,1:3))	! calculate xi point
193        END DO
194      END SUBROUTINE BUILD_HEXA_XI
195
196
197      SUBROUTINE BUILD_HEXA_QUADRATURE(AA)
198!       Subprogram BUILD_HEXA_QUADRATURE - Written by David A. Nordsletten (C) 2008
199!       DPhil Student, Comlab, University of Oxford 2006-2008
200!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
201!
202!       This routine builds the volume/facet quadrature rule for the tensor basis AA
203!
204        USE CMGUI_VARS
205        TYPE(ARRAY_PROBLEM_BASE) AA
206        INTEGER I,J,K,S,m,n,G
207        DOUBLE PRECISION Ax(AA%n_ptsl),Bx(AA%n_ptsl)
208        Ax(1:AA%n_ptsl) = AA%gpt(1:AA%n_ptsl,1)
209        Bx(1:AA%n_ptsl) = AA%gw(1:AA%n_ptsl)
210        n = 0
211        DO I = 1,AA%n_ptsl	! Loop over basis DIM 3
212           IF ((AA%DM == 3).OR.(I == 1)) THEN	! If the dimension is 3 or I = 1
213           DO J = 1,AA%n_ptsl	! Loop over basis DIM 2
214              IF ((AA%DM >= 2).OR.(J == 1)) THEN	! If the dimension is > 2 or J = 1
215              DO K = 1,AA%n_ptsl	! Loop over basis DIM 1
216                 IF (AA%DM >= 1) THEN		! If the dimension is > 1 or J = 1
217                   n = n + 1
218                   AA%gpt(n,1:3) = 0; AA%gw(n) = 0	! Initialize volume quadrature rule, and build
219                   IF (AA%DM >= 1) AA%gpt(n,1) = Ax(K)
220                   IF (AA%DM >= 2) AA%gpt(n,2) = Ax(J)
221                   IF (AA%DM >= 3) AA%gpt(n,3) = Ax(I)
222                   IF (AA%DM == 1) AA%gw(n) = Bx(K)
223                   IF (AA%DM == 2) AA%gw(n) = Bx(K)*Bx(J)
224                   IF (AA%DM == 3) AA%gw(n) = Bx(K)*Bx(J)*Bx(I)
225                 END IF
226              END DO
227              END IF
228           END DO
229           END IF
230        END DO
231        m = 0
232        DO I = 1,AA%n_ptsl	! Loop over basis DIM 2
233           IF ((AA%DM-1 == 2).OR.(I == 1)) THEN		! If the dimension is 2 or I = 1
234           DO J = 1,AA%n_ptsl		! Loop over basis DIM 1
235              IF ((AA%DM-1 >= 1).OR.(J == 1)) THEN	! If the dimension is > 1 or J = 1
236                m = m + 1
237                AA%gpt_f(m,1:3,1) = 0; AA%gw_f(m) = 0	! Initialize facet quadrature rule, and build
238                IF (AA%DM-1 >= 1) AA%gpt_f(m,1,1) = Ax(J)
239                IF (AA%DM-1 >= 2) AA%gpt_f(m,2,1) = Ax(I)
240                IF (AA%DM-1 == 1) AA%gw_f(m) = Bx(J)
241                IF (AA%DM-1 == 2) AA%gw_f(m) = Bx(J)*Bx(I)
242              END IF
243           END DO
244           END IF
245        END DO
246        AA%VL_f(1:AA%FACES) = AA%VL**REAL(AA%DM-1)	! Initialize facet volume
247        AA%VL = AA%VL**REAL(AA%DM)	! Initialize volume
248      END SUBROUTINE BUILD_HEXA_QUADRATURE
249
250
251      SUBROUTINE BUILD_FACE_GPT_ARRAY(AA)
252!       Subprogram BUILD_FACE_GPT_ARRAY - Written by David A. Nordsletten (C) 2008
253!       DPhil Student, Comlab, University of Oxford 2006-2008
254!       PhD Student, Bioengineering Institute, University of Auckland 2005-2006
255!
256!       This routine builds facet quadrature (must be preceeded by a call to BUILD_HEXA_QUADRATURE,
257!       if necessary).
258!
259        USE CMGUI_VARS
260        TYPE(ARRAY_PROBLEM_BASE) AA
261        INTEGER I,J,K
262        DOUBLE PRECISION Ax(3),Bx,Cx(3),Dx
263        AA%gpt_f(1:AA%n_pts_f,AA%DM:3,1) = 0_DPR	! Initialize first quadrature rule
264        DO I = 1,AA%FACES
265           AA%nrm(1:3,I) = 0_DPR			! Initialize normal vectors
266        END DO
267        IF (AA%DM == 3) THEN		! If in R^3
268         IF (AA%FACES == 4) THEN	! master element is a tetrahedron, build appropriate facet quadrature
269          AA%nrm(3,1) = -1_DPR
270
271          AA%gpt_f(1:AA%n_pts_f,1,2) = AA%gpt_f(1:AA%n_pts_f,1,1)
272          AA%gpt_f(1:AA%n_pts_f,2,2) = 0_DPR
273          AA%gpt_f(1:AA%n_pts_f,3,2) = AA%gpt_f(1:AA%n_pts_f,2,1)
274          AA%nrm(2,2) = -1_DPR
275
276
277          AA%gpt_f(1:AA%n_pts_f,1,3) = 0_DPR
278          AA%gpt_f(1:AA%n_pts_f,2,3) = AA%gpt_f(1:AA%n_pts_f,1,1)
279          AA%gpt_f(1:AA%n_pts_f,3,3) = AA%gpt_f(1:AA%n_pts_f,2,1)
280          AA%nrm(1,3) = -1_DPR
281
282          DO I = 1,AA%n_pts_f
283          AA%gpt_f(I,1,4) = 1_DPR-AA%gpt_f(I,1,1)-AA%gpt_f(I,2,1)
284          END DO
285          AA%gpt_f(1:AA%n_pts_f,2,4) = AA%gpt_f(1:AA%n_pts_f,1,1)
286          AA%gpt_f(1:AA%n_pts_f,3,4) = AA%gpt_f(1:AA%n_pts_f,2,1)
287          Ax(1) = 1_DPR
288          Ax(2) = 3_DPR
289          Ax(2) = Ax(2)**0.5
290          AA%nrm(1:3,4) = Ax(1)/Ax(2)
291         ELSE IF (AA%FACES == 6) THEN	! master element is a cube, build appropriate facet quadrature
292          AA%nrm(3,1) = -1_DPR
293
294          AA%gpt_f(1:AA%n_pts_f,1,2) = AA%gpt_f(1:AA%n_pts_f,1,1)
295          AA%gpt_f(1:AA%n_pts_f,2,2) = 0_DPR
296          AA%gpt_f(1:AA%n_pts_f,3,2) = AA%gpt_f(1:AA%n_pts_f,2,1)
297          AA%nrm(2,2) = -1_DPR
298
299          AA%gpt_f(1:AA%n_pts_f,1,3) = 0_DPR
300          AA%gpt_f(1:AA%n_pts_f,2,3) = AA%gpt_f(1:AA%n_pts_f,1,1)
301          AA%gpt_f(1:AA%n_pts_f,3,3) = AA%gpt_f(1:AA%n_pts_f,2,1)
302          AA%nrm(1,3) = -1_DPR
303
304          AA%gpt_f(1:AA%n_pts_f,1,4) = 1_DPR
305          AA%gpt_f(1:AA%n_pts_f,2,4) = AA%gpt_f(1:AA%n_pts_f,1,1)
306          AA%gpt_f(1:AA%n_pts_f,3,4) = AA%gpt_f(1:AA%n_pts_f,2,1)
307          AA%nrm(1,4) = 1_DPR
308
309          AA%gpt_f(1:AA%n_pts_f,1,5) = AA%gpt_f(1:AA%n_pts_f,1,1)
310          AA%gpt_f(1:AA%n_pts_f,2,5) = 1_DPR
311          AA%gpt_f(1:AA%n_pts_f,3,5) = AA%gpt_f(1:AA%n_pts_f,2,1)
312          AA%nrm(2,5) = 1_DPR
313
314          AA%gpt_f(1:AA%n_pts_f,1,6) = AA%gpt_f(1:AA%n_pts_f,1,1)
315          AA%gpt_f(1:AA%n_pts_f,2,6) = AA%gpt_f(1:AA%n_pts_f,2,1)
316          AA%gpt_f(1:AA%n_pts_f,3,6) = 1_DPR
317          AA%nrm(3,6) = 1_DPR
318         END IF
319        ELSE IF (AA%DM == 2) THEN	! If in R^2
320         IF (AA%FACES == 3) THEN	! master element is a triangle, build appropriate facet quadrature
321          AA%nrm(2,1) = -1_DPR
322
323          AA%gpt_f(1:AA%n_pts_f,1,2) = 0_DPR
324          AA%gpt_f(1:AA%n_pts_f,2,2) = AA%gpt_f(1:AA%n_pts_f,1,1)
325          AA%gpt_f(1:AA%n_pts_f,3,2) = 0_DPR
326          AA%nrm(1,2) = -1_DPR
327
328          DO I = 1,AA%n_pts_f
329          AA%gpt_f(I,1,3) = 1_DPR-AA%gpt_f(I,1,1)
330          END DO
331          AA%gpt_f(1:AA%n_pts_f,2,3) = AA%gpt_f(1:AA%n_pts_f,1,1)
332          AA%gpt_f(1:AA%n_pts_f,3,3) = 0_DPR
333          Ax(1) = 1_DPR
334          Ax(2) = 2_DPR
335          Ax(2) = Ax(2)**0.5
336          AA%nrm(1:2,3) = Ax(1)/Ax(2); AA%nrm(3,3) = 0_DPR
337         ELSE IF (AA%FACES == 4) THEN	! master element is a square, build appropriate facet quadrature
338          AA%nrm(2,1) = -1_DPR
339
340          AA%gpt_f(1:AA%n_pts_f,1,2) = AA%gpt_f(1:AA%n_pts_f,1,1)
341          AA%gpt_f(1:AA%n_pts_f,2,2) = 0_DPR
342          AA%gpt_f(1:AA%n_pts_f,3,2) = 0_DPR
343          AA%nrm(1,2) = -1_DPR
344
345          AA%gpt_f(1:AA%n_pts_f,1,3) = AA%gpt_f(1:AA%n_pts_f,1,1)
346          AA%gpt_f(1:AA%n_pts_f,2,3) = 1_DPR
347          AA%gpt_f(1:AA%n_pts_f,3,3) = 0_DPR
348          AA%nrm(2,3) = 1_DPR
349
350          AA%gpt_f(1:AA%n_pts_f,1,4) = 1_DPR
351          AA%gpt_f(1:AA%n_pts_f,2,4) = AA%gpt_f(1:AA%n_pts_f,1,1)
352          AA%gpt_f(1:AA%n_pts_f,3,4) = 0_DPR
353          AA%nrm(1,4) = 1_DPR
354         END IF
355        END IF
356      END SUBROUTINE BUILD_FACE_GPT_ARRAY
357
358
359
360
361
362