/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. END DO
  128. RETURN
  129. 90 CALL ERRORZ(5)
  130. END SUBROUTINE LOAD_BASE_SET
  131. SUBROUTINE LOAD_INTEGRALS(AA,VFLD)
  132. ! Subprogram LOAD_INTEGRALS - Written by David A. Nordsletten (C) 2008
  133. ! DPhil Student, Comlab, University of Oxford 2006-2008
  134. ! PhD Student, Bioengineering Institute, University of Auckland 2005-2006
  135. !
  136. ! This program loads integrals for base AA%B(VFLD). The basis is evaluated a given quadrature points.
  137. !
  138. USE CMGUI_VARS
  139. TYPE(ARRAY_PROBLEM_BASE) AA
  140. INTEGER I,J,K,x,y,z,ID,ORD1(3),VFLD,S
  141. DOUBLE PRECISION Ax(3),Bx,Cx(3),Dx
  142. ALLOCATE( AA%B(VFLD)%I%Y(AA%B(VFLD)%n,AA%n_pts) )
  143. ALLOCATE( AA%B(VFLD)%I%Y_f(AA%B(VFLD)%n,AA%n_pts_f,AA%FACES) )
  144. ALLOCATE( AA%B(VFLD)%I%dY(3,AA%B(VFLD)%n,AA%n_pts) )
  145. ALLOCATE( AA%B(VFLD)%I%TdY(3,AA%B(VFLD)%n,AA%n_pts) )
  146. ALLOCATE( AA%B(VFLD)%I%dY_f(3,AA%B(VFLD)%n,AA%n_pts_f,AA%FACES) )
  147. DO I = 1,AA%B(VFLD)%n ! Loop over basis
  148. ORD1(1:3) = 0
  149. DO J = 1,AA%n_pts ! Loop over quadrature points
  150. CALL POLY(AA,VFLD,I,AA%gpt(J,1:3),ORD1,AA%B(VFLD)%I%Y(I,J)) ! Calculate basis
  151. END DO
  152. DO K = 1,3
  153. ORD1(1:3) = 0; ORD1(K) = 1
  154. DO J = 1,AA%n_pts ! Loop over quadrature points
  155. CALL POLY(AA,VFLD,I,AA%gpt(J,1:3),ORD1,AA%B(VFLD)%I%dY(K,I,J)) ! Calculate basis
  156. END DO
  157. END DO
  158. END DO
  159. DO S = 1,AA%FACES ! Loop over facets
  160. DO I = 1,AA%B(VFLD)%n ! Loop over basis
  161. ORD1(1:3) = 0
  162. DO J = 1,AA%n_pts_f ! Loop over facet quadrature points
  163. CALL POLY(AA,VFLD,I,AA%gpt_f(J,1:3,S),ORD1,AA%B(VFLD)%I%Y_f(I,J,S)) ! Calculate basis
  164. END DO
  165. DO K = 1,3
  166. ORD1(1:3) = 0; ORD1(K) = 1
  167. DO J = 1,AA%n_pts_f ! Loop over quadrature points
  168. 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
  169. END DO
  170. END DO
  171. END DO
  172. END DO
  173. END SUBROUTINE LOAD_INTEGRALS
  174. SUBROUTINE BUILD_HEXA_XI(AA,VFLD)
  175. ! Subprogram BUILD_HEXA_XI - Written by David A. Nordsletten (C) 2008
  176. ! DPhil Student, Comlab, University of Oxford 2006-2008
  177. ! PhD Student, Bioengineering Institute, University of Auckland 2005-2006
  178. !
  179. ! This routine builds the XI coordinate for the tensor basis AA%B(VFLD)
  180. !
  181. USE CMGUI_VARS
  182. TYPE(ARRAY_PROBLEM_BASE) AA
  183. INTEGER I,J,K,S,m,n,G,VFLD
  184. DOUBLE PRECISION Bx(AA%B(VFLD)%nl)
  185. Bx(1:AA%B(VFLD)%nl) = AA%B(VFLD)%XI(1:AA%B(VFLD)%nl,1)
  186. DO I = 1,AA%B(VFLD)%n ! Loop over basis
  187. AA%B(VFLD)%XI(I,1:3) = Bx(AA%B(VFLD)%B_ID(I,1:3)) ! calculate xi point
  188. END DO
  189. END SUBROUTINE BUILD_HEXA_XI
  190. SUBROUTINE BUILD_HEXA_QUADRATURE(AA)
  191. ! Subprogram BUILD_HEXA_QUADRATURE - Written by David A. Nordsletten (C) 2008
  192. ! DPhil Student, Comlab, University of Oxford 2006-2008
  193. ! PhD Student, Bioengineering Institute, University of Auckland 2005-2006
  194. !
  195. ! This routine builds the volume/facet quadrature rule for the tensor basis AA
  196. !
  197. USE CMGUI_VARS
  198. TYPE(ARRAY_PROBLEM_BASE) AA
  199. INTEGER I,J,K,S,m,n,G
  200. DOUBLE PRECISION Ax(AA%n_ptsl),Bx(AA%n_ptsl)
  201. Ax(1:AA%n_ptsl) = AA%gpt(1:AA%n_ptsl,1)
  202. Bx(1:AA%n_ptsl) = AA%gw(1:AA%n_ptsl)
  203. n = 0
  204. DO I = 1,AA%n_ptsl ! Loop over basis DIM 3
  205. IF ((AA%DM == 3).OR.(I == 1)) THEN ! If the dimension is 3 or I = 1
  206. DO J = 1,AA%n_ptsl ! Loop over basis DIM 2
  207. IF ((AA%DM >= 2).OR.(J == 1)) THEN ! If the dimension is > 2 or J = 1
  208. DO K = 1,AA%n_ptsl ! Loop over basis DIM 1
  209. IF (AA%DM >= 1) THEN ! If the dimension is > 1 or J = 1
  210. n = n + 1
  211. AA%gpt(n,1:3) = 0; AA%gw(n) = 0 ! Initialize volume quadrature rule, and build
  212. IF (AA%DM >= 1) AA%gpt(n,1) = Ax(K)
  213. IF (AA%DM >= 2) AA%gpt(n,2) = Ax(J)
  214. IF (AA%DM >= 3) AA%gpt(n,3) = Ax(I)
  215. IF (AA%DM == 1) AA%gw(n) = Bx(K)
  216. IF (AA%DM == 2) AA%gw(n) = Bx(K)*Bx(J)
  217. IF (AA%DM == 3) AA%gw(n) = Bx(K)*Bx(J)*Bx(I)
  218. END IF
  219. END DO
  220. END IF
  221. END DO
  222. END IF
  223. END DO
  224. m = 0
  225. DO I = 1,AA%n_ptsl ! Loop over basis DIM 2
  226. IF ((AA%DM-1 == 2).OR.(I == 1)) THEN ! If the dimension is 2 or I = 1
  227. DO J = 1,AA%n_ptsl ! Loop over basis DIM 1
  228. IF ((AA%DM-1 >= 1).OR.(J == 1)) THEN ! If the dimension is > 1 or J = 1
  229. m = m + 1
  230. AA%gpt_f(m,1:3,1) = 0; AA%gw_f(m) = 0 ! Initialize facet quadrature rule, and build
  231. IF (AA%DM-1 >= 1) AA%gpt_f(m,1,1) = Ax(J)
  232. IF (AA%DM-1 >= 2) AA%gpt_f(m,2,1) = Ax(I)
  233. IF (AA%DM-1 == 1) AA%gw_f(m) = Bx(J)
  234. IF (AA%DM-1 == 2) AA%gw_f(m) = Bx(J)*Bx(I)
  235. END IF
  236. END DO
  237. END IF
  238. END DO
  239. AA%VL_f(1:AA%FACES) = AA%VL**REAL(AA%DM-1) ! Initialize facet volume
  240. AA%VL = AA%VL**REAL(AA%DM) ! Initialize volume
  241. END SUBROUTINE BUILD_HEXA_QUADRATURE
  242. SUBROUTINE BUILD_FACE_GPT_ARRAY(AA)
  243. ! Subprogram BUILD_FACE_GPT_ARRAY - 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 routine builds facet quadrature (must be preceeded by a call to BUILD_HEXA_QUADRATURE,
  248. ! if necessary).
  249. !
  250. USE CMGUI_VARS
  251. TYPE(ARRAY_PROBLEM_BASE) AA
  252. INTEGER I,J,K
  253. DOUBLE PRECISION Ax(3),Bx,Cx(3),Dx
  254. AA%gpt_f(1:AA%n_pts_f,AA%DM:3,1) = 0_DPR ! Initialize first quadrature rule
  255. DO I = 1,AA%FACES
  256. AA%nrm(1:3,I) = 0_DPR ! Initialize normal vectors
  257. END DO
  258. IF (AA%DM == 3) THEN ! If in R^3
  259. IF (AA%FACES == 4) THEN ! master element is a tetrahedron, build appropriate facet quadrature
  260. AA%nrm(3,1) = -1_DPR
  261. AA%gpt_f(1:AA%n_pts_f,1,2) = AA%gpt_f(1:AA%n_pts_f,1,1)
  262. AA%gpt_f(1:AA%n_pts_f,2,2) = 0_DPR
  263. AA%gpt_f(1:AA%n_pts_f,3,2) = AA%gpt_f(1:AA%n_pts_f,2,1)
  264. AA%nrm(2,2) = -1_DPR
  265. AA%gpt_f(1:AA%n_pts_f,1,3) = 0_DPR
  266. AA%gpt_f(1:AA%n_pts_f,2,3) = AA%gpt_f(1:AA%n_pts_f,1,1)
  267. AA%gpt_f(1:AA%n_pts_f,3,3) = AA%gpt_f(1:AA%n_pts_f,2,1)
  268. AA%nrm(1,3) = -1_DPR
  269. DO I = 1,AA%n_pts_f
  270. AA%gpt_f(I,1,4) = 1_DPR-AA%gpt_f(I,1,1)-AA%gpt_f(I,2,1)
  271. END DO
  272. AA%gpt_f(1:AA%n_pts_f,2,4) = AA%gpt_f(1:AA%n_pts_f,1,1)
  273. AA%gpt_f(1:AA%n_pts_f,3,4) = AA%gpt_f(1:AA%n_pts_f,2,1)
  274. Ax(1) = 1_DPR
  275. Ax(2) = 3_DPR
  276. Ax(2) = Ax(2)**0.5
  277. AA%nrm(1:3,4) = Ax(1)/Ax(2)
  278. ELSE IF (AA%FACES == 6) THEN ! master element is a cube, build appropriate facet quadrature
  279. AA%nrm(3,1) = -1_DPR
  280. AA%gpt_f(1:AA%n_pts_f,1,2) = AA%gpt_f(1:AA%n_pts_f,1,1)
  281. AA%gpt_f(1:AA%n_pts_f,2,2) = 0_DPR
  282. AA%gpt_f(1:AA%n_pts_f,3,2) = AA%gpt_f(1:AA%n_pts_f,2,1)
  283. AA%nrm(2,2) = -1_DPR
  284. AA%gpt_f(1:AA%n_pts_f,1,3) = 0_DPR
  285. AA%gpt_f(1:AA%n_pts_f,2,3) = AA%gpt_f(1:AA%n_pts_f,1,1)
  286. AA%gpt_f(1:AA%n_pts_f,3,3) = AA%gpt_f(1:AA%n_pts_f,2,1)
  287. AA%nrm(1,3) = -1_DPR
  288. AA%gpt_f(1:AA%n_pts_f,1,4) = 1_DPR
  289. AA%gpt_f(1:AA%n_pts_f,2,4) = AA%gpt_f(1:AA%n_pts_f,1,1)
  290. AA%gpt_f(1:AA%n_pts_f,3,4) = AA%gpt_f(1:AA%n_pts_f,2,1)
  291. AA%nrm(1,4) = 1_DPR
  292. AA%gpt_f(1:AA%n_pts_f,1,5) = AA%gpt_f(1:AA%n_pts_f,1,1)
  293. AA%gpt_f(1:AA%n_pts_f,2,5) = 1_DPR
  294. AA%gpt_f(1:AA%n_pts_f,3,5) = AA%gpt_f(1:AA%n_pts_f,2,1)
  295. AA%nrm(2,5) = 1_DPR
  296. AA%gpt_f(1:AA%n_pts_f,1,6) = AA%gpt_f(1:AA%n_pts_f,1,1)
  297. AA%gpt_f(1:AA%n_pts_f,2,6) = AA%gpt_f(1:AA%n_pts_f,2,1)
  298. AA%gpt_f(1:AA%n_pts_f,3,6) = 1_DPR
  299. AA%nrm(3,6) = 1_DPR
  300. END IF
  301. ELSE IF (AA%DM == 2) THEN ! If in R^2
  302. IF (AA%FACES == 3) THEN ! master element is a triangle, build appropriate facet quadrature
  303. AA%nrm(2,1) = -1_DPR
  304. AA%gpt_f(1:AA%n_pts_f,1,2) = 0_DPR
  305. AA%gpt_f(1:AA%n_pts_f,2,2) = AA%gpt_f(1:AA%n_pts_f,1,1)
  306. AA%gpt_f(1:AA%n_pts_f,3,2) = 0_DPR
  307. AA%nrm(1,2) = -1_DPR
  308. DO I = 1,AA%n_pts_f
  309. AA%gpt_f(I,1,3) = 1_DPR-AA%gpt_f(I,1,1)
  310. END DO
  311. AA%gpt_f(1:AA%n_pts_f,2,3) = AA%gpt_f(1:AA%n_pts_f,1,1)
  312. AA%gpt_f(1:AA%n_pts_f,3,3) = 0_DPR
  313. Ax(1) = 1_DPR
  314. Ax(2) = 2_DPR
  315. Ax(2) = Ax(2)**0.5
  316. AA%nrm(1:2,3) = Ax(1)/Ax(2); AA%nrm(3,3) = 0_DPR
  317. ELSE IF (AA%FACES == 4) THEN ! master element is a square, build appropriate facet quadrature
  318. AA%nrm(2,1) = -1_DPR
  319. AA%gpt_f(1:AA%n_pts_f,1,2) = AA%gpt_f(1:AA%n_pts_f,1,1)
  320. AA%gpt_f(1:AA%n_pts_f,2,2) = 0_DPR
  321. AA%gpt_f(1:AA%n_pts_f,3,2) = 0_DPR
  322. AA%nrm(1,2) = -1_DPR
  323. AA%gpt_f(1:AA%n_pts_f,1,3) = AA%gpt_f(1:AA%n_pts_f,1,1)
  324. AA%gpt_f(1:AA%n_pts_f,2,3) = 1_DPR
  325. AA%gpt_f(1:AA%n_pts_f,3,3) = 0_DPR
  326. AA%nrm(2,3) = 1_DPR
  327. AA%gpt_f(1:AA%n_pts_f,1,4) = 1_DPR
  328. AA%gpt_f(1:AA%n_pts_f,2,4) = AA%gpt_f(1:AA%n_pts_f,1,1)
  329. AA%gpt_f(1:AA%n_pts_f,3,4) = 0_DPR
  330. AA%nrm(1,4) = 1_DPR
  331. END IF
  332. END IF
  333. END SUBROUTINE BUILD_FACE_GPT_ARRAY