/FluidMechanics/TOOLS/MeshCreation/02_CreateCMInput/cheart_visbasis3D.f
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