1SUBROUTINE LinearFormsAssembly( Model,Solver,dt,TransientSimulation ) 2!------------------------------------------------------------------------------ 3!****************************************************************************** 4! 5! Unit test for linear form / vectorized basis function computation 6! routines in Elmer 7! 8! ARGUMENTS: 9! 10! TYPE(Model_t) :: Model, 11! INPUT: All model information (mesh, materials, BCs, etc...) 12! 13! TYPE(Solver_t) :: Solver 14! INPUT: Linear & nonlinear equation solver options 15! 16! REAL(KIND=dp) :: dt, 17! INPUT: Timestep size for time dependent simulations 18! 19! LOGICAL :: TransientSimulation 20! INPUT: Steady state or transient simulation 21! 22!****************************************************************************** 23 USE DefUtils 24 USE LinearForms 25 USE ISO_C_BINDING 26#ifdef _OPENMP 27 USE omp_lib 28#endif 29!------------------------------------------------------------------------------ 30 IMPLICIT NONE 31!------------------------------------------------------------------------------ 32 TYPE(Model_t) :: Model 33 TYPE(Solver_t) :: Solver 34 35 REAL(KIND=dp) :: dt 36 LOGICAL :: TransientSimulation 37!------------------------------------------------------------------------------ 38! Local variables 39!------------------------------------------------------------------------------ 40 REAL(KIND=dp), PARAMETER :: tol1d = 1D-12, tol2d=1D-12, tol3d=1D-12 41 REAL(KIND=dp) :: float_P 42 INTEGER :: nerror, netest, P 43 LOGICAL :: Found 44 45 nerror = 0 46 float_P = ListGetCReal(GetSolverParams(), 'P', Found) 47 IF (.NOT. Found) float_P = 6.0_dp 48 P = NINT(float_P) 49 50 ! 1D tests 51 netest = TestLineElement(Solver, P, tol1d) 52 IF (netest /= 0) THEN 53 CALL Warn('LinearFormsAssembly','Line element contained errors') 54 END IF 55 nerror = nerror + netest 56 57 ! ! 2D tests 58 netest = TestTriangleElement(Solver, P, tol2d) 59 IF (netest /= 0) THEN 60 CALL Warn('LinearFormsAssembly','Triangle element contained errors') 61 END IF 62 nerror = nerror + netest 63 64 netest = TestQuadElement(Solver, P, tol2d) 65 IF (netest /= 0) THEN 66 CALL Warn('LinearFormsAssembly','Quad element contained errors') 67 END IF 68 nerror = nerror + netest 69 70 ! ! 3D tests 71 netest = TestTetraElement(Solver, P, tol3d) 72 IF (netest /= 0) THEN 73 CALL Warn('LinearFormsAssembly','Tetra element contained errors') 74 END IF 75 nerror = nerror + netest 76 77 netest = TestWedgeElement(Solver, P, tol3d) 78 IF (netest /= 0) THEN 79 CALL Warn('LinearFormsAssembly','Wedge element contained errors') 80 END IF 81 nerror = nerror + netest 82 83 netest = TestBrickElement(Solver, P, tol3d) 84 IF (netest /= 0) THEN 85 CALL Warn('LinearFormsAssembly','Brick element contained errors') 86 END IF 87 nerror = nerror + netest 88 89 ! Build solution norm for error checking 90 Solver % Variable % Norm = REAL(1+nerror,dp) 91 Solver % Variable % Values = REAL(1+nerror,dp) 92 93CONTAINS 94 95 FUNCTION TestLineElement(Solver, P, tol) RESULT(nerror) 96 IMPLICIT NONE 97 98 TYPE(Solver_t) :: Solver 99 INTEGER, INTENT(IN) :: P 100 REAL(kind=dp), INTENT(IN) :: tol 101 INTEGER :: nerror 102 103 nerror = TestElement(Solver, 202, P, tol) 104 END FUNCTION TestLineElement 105 106 FUNCTION TestTriangleElement(Solver, P, tol) RESULT(nerror) 107 IMPLICIT NONE 108 109 TYPE(Solver_t) :: Solver 110 INTEGER, INTENT(IN) :: P 111 REAL(kind=dp), INTENT(IN) :: tol 112 INTEGER :: nerror 113 114 nerror = TestElement(Solver, 303, P, tol) 115 END FUNCTION TestTriangleElement 116 117 FUNCTION TestQuadElement(Solver, P, tol) RESULT(nerror) 118 IMPLICIT NONE 119 120 TYPE(Solver_t) :: Solver 121 INTEGER, INTENT(IN) :: P 122 REAL(kind=dp), INTENT(IN) :: tol 123 INTEGER :: nerror 124 125 nerror = TestElement(Solver, 404, P, tol) 126 END FUNCTION TestQuadElement 127 128 FUNCTION TestTetraElement(Solver, P, tol) RESULT(nerror) 129 IMPLICIT NONE 130 131 TYPE(Solver_t) :: Solver 132 INTEGER, INTENT(IN) :: P 133 REAL(kind=dp), INTENT(IN) :: tol 134 INTEGER :: nerror 135 136 nerror = TestElement(Solver, 504, P, tol) 137 END FUNCTION TestTetraElement 138 139 FUNCTION TestWedgeElement(Solver, P, tol) RESULT(nerror) 140 IMPLICIT NONE 141 142 TYPE(Solver_t) :: Solver 143 INTEGER, INTENT(IN) :: P 144 REAL(kind=dp), INTENT(IN) :: tol 145 INTEGER :: nerror 146 147 nerror = TestElement(Solver, 706, P, tol) 148 END FUNCTION TestWedgeElement 149 150 FUNCTION TestBrickElement(Solver, P, tol) RESULT(nerror) 151 IMPLICIT NONE 152 153 TYPE(Solver_t) :: Solver 154 INTEGER, INTENT(IN) :: P 155 REAL(kind=dp), INTENT(IN) :: tol 156 INTEGER :: nerror 157 158 nerror = TestElement(Solver, 808, P, tol) 159 END FUNCTION TestBrickElement 160 161 FUNCTION TestElement(Solver, ecode, P, tol) RESULT(nerror) 162 IMPLICIT NONE 163 164 TYPE(Solver_t) :: Solver 165 INTEGER, INTENT(IN) :: ecode 166 INTEGER, INTENT(IN) :: P 167 REAL(kind=dp), INTENT(IN) :: tol 168 169 TYPE(Element_t), POINTER :: Element, SingleElement 170 TYPE(Mesh_t), POINTER :: NewMesh, OldMesh 171 REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), LOAD(:), FORCE(:), & 172 STIFFvec(:,:), FORCEvec(:) 173 174 INTEGER :: i, j, k, l, q, nerror, nbasis, nndof, allocstat, tag, nthr, & 175 nbasisvec, ndbasisdxvec, rep, dim, lm_eval, lm_eval_vec, NumGP 176 177 INTEGER, PARAMETER :: NREP = 100 178 REAL(kind=dp) :: t_start, t_end, t_tot, t_startvec, t_endvec, t_tot_vec 179 TYPE(GaussIntegrationPoints_t) :: Quadrature 180 181 nerror = 0 182 lm_eval = 0 183 lm_eval_vec = 0 184 t_tot = REAL(0,dp) 185 t_tot_vec = REAL(0,dp) 186 187 ! Create a mesh with a single element 188 OldMesh => Solver % Mesh 189 NewMesh => NULL() 190 CALL AllocateMeshAndPElement(NewMesh, ecode, P, SingleElement) 191 Solver % Mesh => NewMesh 192 193 ! Insert P element definitions to Solver mapping (sets P elements as "active") 194 IF (ALLOCATED(Solver % Def_Dofs)) THEN 195 tag = ecode / 100 196 Solver % Def_Dofs(tag,1,6) = P 197 END IF 198 IF (ecode==404) THEN 199 Quadrature = GaussPointsQuad((P+1)**2) 200 ELSE 201 Quadrature = GaussPoints(SingleElement, PReferenceElement=.TRUE.) 202 END IF 203 NumGP = Quadrature % N 204 205 !$OMP PARALLEL SHARED(Solver, SingleElement, ecode, tol) & 206 !$OMP PRIVATE(STIFF, FORCE, STIFFvec, FORCEvec, & 207 !$OMP LOAD, Element, nndof, nbasis, & 208 !$OMP allocstat, rep, t_start, t_end, & 209 !$OMP t_startvec, t_endvec) & 210 !$OMP REDUCTION(+:nerror,t_tot,t_tot_vec,lm_eval,lm_eval_vec) & 211 !$OMP DEFAULT(NONE) 212 213 ! Construct a temporary element based on the one created previously 214 Element => ClonePElement(SingleElement) 215 216 nndof = Element % Type % NumberOfNodes 217 nbasis = nndof + GetElementNOFDOFs( Element ) 218 219 ! Reserve workspace 220 ALLOCATE(STIFF(nbasis, nbasis), FORCE(nbasis), & 221 STIFFvec(nbasis, nbasis), FORCEvec(nbasis), & 222 LOAD(nndof), & 223 STAT=allocstat) 224 IF (allocstat /= 0) THEN 225 CALL Fatal('LinearForms',& 226 'Storage allocation for local matrices failed') 227 END IF 228 229 ! Initialize artificial load vector 230 LOAD = REAL(1,dp) 231 232 ! Warmup 233 CALL LocalMatrix( STIFF, FORCE, LOAD, Element, nndof, nbasis) 234 !$OMP BARRIER 235 t_start = ftimer() 236 DO rep=1,NREP 237 ! Construct local matrix 238!DIR$ NOINLINE 239 CALL LocalMatrix( STIFF, FORCE, LOAD, Element, nndof, nbasis) 240 END DO 241 t_end = ftimer() 242 lm_eval = NREP 243 244 ! Warmup 245 CALL LocalMatrixVec( STIFFvec, FORCEvec, LOAD, Element, nndof, nbasis) 246 !$OMP BARRIER 247 t_startvec = ftimer() 248 DO rep=1,NREP 249 ! Construct local matrix 250!DIR$ NOINLINE 251 CALL LocalMatrixVec( STIFFvec, FORCEvec, LOAD, Element, nndof, nbasis) 252 END DO 253 t_endvec = ftimer() 254 lm_eval_vec = NREP 255 256 nerror = TestLocalMatrix(nbasis, STIFF, STIFFvec, tol) 257 nerror = nerror + TestLocalForce(nbasis, FORCE, FORCEvec, tol) 258 259 t_tot = t_end - t_start 260 t_tot_vec = t_endvec - t_startvec 261 262 CALL DeallocatePElement(Element) 263 DEALLOCATE(STIFF, FORCE, STIFFvec, FORCEvec, LOAD) 264 265 !$OMP END PARALLEL 266 267 ! Normalize the times / thread 268 nthr = 1 269 !$ nthr = omp_get_max_threads() 270 t_tot = t_tot / nthr 271 t_tot_vec = t_tot_vec / nthr 272 CALL PrintTestData(SingleElement, NumGP, t_tot, lm_eval, & 273 t_tot_vec, lm_eval_vec) 274 275 CALL DeallocateTemporaryMesh(NewMesh) 276 Solver % Mesh => OldMesh 277 CALL DeallocatePElement(SingleElement) 278 279 IF (ALLOCATED(Solver % Def_Dofs)) THEN 280 tag = ecode / 100 281 Solver % Def_Dofs(tag,1,6) = 0 282 END IF 283 END FUNCTION TestElement 284 285 SUBROUTINE LocalMatrix( STIFF, FORCE, LOAD, Element, n, nd ) 286 IMPLICIT NONE 287 REAL(KIND=dp) CONTIG :: STIFF(:,:), FORCE(:) 288 REAL(KIND=dp) CONTIG, INTENT(IN) :: LOAD(:) 289 INTEGER :: n, nd 290 TYPE(Element_t), POINTER :: Element 291!------------------------------------------------------------------------------ 292 REAL(KIND=dp) :: Basis(nd),dBasisdx(nd,3),DetJ,LoadAtIP,Weight 293 LOGICAL :: Stat 294 INTEGER :: i,j,t,p,q,dim 295 TYPE(GaussIntegrationPoints_t) :: IP 296 297 TYPE(Nodes_t), SAVE :: Nodes 298 !$OMP THREADPRIVATE(Nodes) 299!------------------------------------------------------------------------------ 300 CALL GetReferenceElementNodes( Nodes, Element ) 301 STIFF = 0.0d0 302 FORCE = 0.0d0 303 304 dim = Element % Type % Dimension 305 !Numerical integration: 306 !---------------------- 307 IF (Element % TYPE % ElementCode / 100 == 4) THEN 308 IP = GaussPointsQuad((Element % PDefs % P+1)**2) 309 ELSE 310 IP = GaussPoints( Element, PReferenceElement=.TRUE. ) 311 END IF 312 313 DO t=1,IP % n 314 315 ! Basis function values & derivatives at the integration point: 316 !-------------------------------------------------------------- 317 stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), & 318 IP % W(t), detJ, Basis, dBasisdx) 319 320 ! The source term at the integration point: 321 !------------------------------------------ 322 LoadAtIP = SUM( Basis(1:n) * LOAD(1:n) ) 323 Weight = IP % s(t) * DetJ 324 ! STIFF=STIFF+(grad u, grad v) 325 STIFF(1:nd,1:nd) = STIFF(1:nd,1:nd) + Weight * & 326 MATMUL( dBasisdx(1:nd,1:dim), TRANSPOSE( dBasisdx(1:nd,1:dim) ) ) 327 328 DO p=1,nd 329 DO q=1,nd 330 ! STIFF=STIFF+(grad u,v) 331 ! ----------------------------------- 332 STIFF (p,q) = STIFF(p,q) + Weight * SUM(dBasisdx(q,1:dim)) * Basis(p) 333 334 ! STIFF=STIFF+(u,v) 335 STIFF(p,q) = STIFF(p,q) + Weight * Basis(q) * Basis(p) 336 END DO 337 END DO 338 339 FORCE(1:nd) = FORCE(1:nd) + IP % s(t) * DetJ * LoadAtIP * Basis(1:nd) 340 END DO 341 END SUBROUTINE LocalMatrix 342 343 SUBROUTINE LocalMatrixVec( STIFF, FORCE, LOAD, Element, n, nd ) 344!------------------------------------------------------------------------------ 345 IMPLICIT NONE 346 REAL(KIND=dp) CONTIG :: STIFF(:,:), FORCE(:) 347 REAL(KIND=dp) CONTIG, INTENT(IN) :: LOAD(:) 348 INTEGER, INTENT(IN) :: n, nd 349 TYPE(Element_t), POINTER :: Element 350!------------------------------------------------------------------------------ 351 LOGICAL :: Stat 352 INTEGER :: i, j, ngp, allocstat, gp 353 TYPE(GaussIntegrationPoints_t) :: IP 354 355 TYPE(Nodes_t), SAVE :: Nodes 356 REAL(KIND=dp), ALLOCATABLE, SAVE :: Basis(:,:), dBasisdx(:,:,:), & 357 DetJ(:), LoadAtIPs(:) 358 !$OMP THREADPRIVATE(Nodes, Basis, dBasisdx, DetJ, LoadAtIPs) 359!------------------------------------------------------------------------------ 360 CALL GetReferenceElementNodes( Nodes, Element ) 361 362 STIFF = REAL(0, dp) 363 FORCE = REAL(0, dp) 364 365 ! Get integration points 366 IF (Element % TYPE % ElementCode / 100 == 4) THEN 367 IP = GaussPointsQuad((Element % PDefs % P+1)**2) 368 ELSE 369 IP = GaussPoints( Element, PReferenceElement=.TRUE. ) 370 END IF 371 ngp = IP % n 372 373 ! Reserve workspace 374 IF (.NOT. ALLOCATED(Basis)) THEN 375 ALLOCATE(Basis(ngp,nd), dBasisdx(ngp,nd,3), & 376 DetJ(ngp), LoadAtIPs(ngp), STAT=allocstat) 377 IF (allocstat /= 0) THEN 378 CALL Fatal('LocalMatrixVec',& 379 'Storage allocation for local element basis failed') 380 END IF 381 ELSE IF (SIZE(Basis,1) /= ngp .OR. SIZE(Basis,2) /= nd) THEN 382 DEALLOCATE(Basis, dBasisdx, DetJ, LoadAtIPs) 383 ALLOCATE(Basis(ngp,nd), dBasisdx(ngp,nd,3), & 384 DetJ(ngp), LoadAtIPs(ngp), STAT=allocstat) 385 IF (allocstat /= 0) THEN 386 CALL Fatal('LocalMatrixVec',& 387 'Storage allocation for local element basis failed') 388 END IF 389 END IF 390 391 ! Compute values of all basis functions at all integration points 392 stat = ElementInfoVec( Element, Nodes, ngp, & 393 IP % U, IP % V, IP % W, DetJ, SIZE(Basis,2), Basis, dBasisdx ) 394 395 396 ! Compute actual integration weights (recycle memory space of DetJ) 397 DO i=1,ngp 398 DetJ(i) = Ip % s(i)*Detj(i) 399 END DO 400 401 ! STIFF=STIFF+(grad u, grad u) 402 CALL LinearForms_GradUdotGradU(ngp, nd, Element % TYPE % DIMENSION, & 403 dBasisdx, DetJ, STIFF) 404 ! STIFF=STIFF+(u,u) 405 CALL LinearForms_UdotU(ngp, nd, Element % TYPE % DIMENSION, Basis, DetJ, STIFF) 406 ! STIFF=STIFF+(grad u,v) 407 CALL LinearForms_GradUdotU(ngp, nd, Element % TYPE % DIMENSION, dBasisdx, Basis, DetJ, STIFF) 408 409 ! Source terms at IPs 410 !------------------------------------------ 411 ! LoadAtIPs(1:ngp) = MATMUL( Basis(1:ngp,1:n), LOAD(1:n) ) 412 CALL LinearForms_ProjectToU(ngp, n, Basis, LOAD, LoadAtIPs) 413 414 ! FORCE=FORCE+(u,f) 415 CALL LinearForms_UdotF(ngp, nd, Basis, DetJ, LoadAtIPs, FORCE) 416 END SUBROUTINE LocalMatrixVec 417 418 FUNCTION TestLocalMatrix(nbasis, STIFF1, STIFF2, tol) RESULT(nerror) 419 IMPLICIT NONE 420 421 INTEGER, INTENT(IN) :: nbasis 422 REAL(KIND=dp) CONTIG, INTENT(IN) :: STIFF1(:,:), STIFF2(:,:) 423 REAL(kind=dp), INTENT(IN) :: tol 424 INTEGER :: nerror 425 426 INTEGER :: i, j, dim 427 428 nerror = 0 429 ! Test element of local matrix 430 DO j=1,nbasis 431 DO i=1,nbasis 432 IF (ABS(STIFF1(i,j)-STIFF2(i,j)) >= tol) THEN 433 nerror = nerror + 1 434 WRITE (*,*) 'STIFF:', i,j,STIFF1(i,j), STIFF2(i,j) 435 END IF 436 END DO 437 END DO 438 END FUNCTION TestLocalMatrix 439 440 FUNCTION TestLocalForce(nbasis, FORCE1, FORCE2, tol) RESULT(nerror) 441 IMPLICIT NONE 442 443 INTEGER, INTENT(IN) :: nbasis 444 REAL(KIND=dp) CONTIG, INTENT(IN) :: FORCE1(:), FORCE2(:) 445 REAL(kind=dp), INTENT(IN) :: tol 446 INTEGER :: nerror 447 448 INTEGER :: i, dim 449 450 nerror = 0 451 ! Test element of local force vector 452 DO i=1,nbasis 453 IF (ABS(FORCE1(i)-FORCE2(i)) >= tol) THEN 454 nerror = nerror + 1 455 WRITE (*,*) 'FORCE:', i, FORCE1(i), FORCE2(i) 456 END IF 457 END DO 458 END FUNCTION TestLocalForce 459 460 SUBROUTINE GetReferenceElementNodes( ElementNodes, Element ) 461 TYPE(Nodes_t), TARGET :: ElementNodes 462 TYPE(Element_t) :: Element 463 464 INTEGER :: i, n, padn, sz, astat 465 466 n = Element % TYPE % NumberOfNodes 467 padn = n 468 469 ! TODO: Implement padding 470 IF (.NOT. ALLOCATED( ElementNodes % xyz)) THEN 471 ! Deallocate old storage 472 IF (ASSOCIATED(ElementNodes % x)) DEALLOCATE(ElementNodes % x) 473 IF (ASSOCIATED(ElementNodes % y)) DEALLOCATE(ElementNodes % y) 474 IF (ASSOCIATED(ElementNodes % z)) DEALLOCATE(ElementNodes % z) 475 476 ! Allocate new storage 477 ALLOCATE(ElementNodes % xyz(padn,3)) 478 ElementNodes % xyz = REAL(0,dp) 479 ElementNodes % x => ElementNodes % xyz(1:n,1) 480 ElementNodes % y => ElementNodes % xyz(1:n,2) 481 ElementNodes % z => ElementNodes % xyz(1:n,3) 482 ELSE IF (SIZE(ElementNodes % xyz, 1)<padn) THEN 483 DEALLOCATE(ElementNodes % xyz) 484 ALLOCATE(ElementNodes % xyz(padn,3)) 485 ElementNodes % xyz = REAL(0,dp) 486 ElementNodes % x => ElementNodes % xyz(1:n,1) 487 ElementNodes % y => ElementNodes % xyz(1:n,2) 488 ElementNodes % z => ElementNodes % xyz(1:n,3) 489 ELSE 490 ElementNodes % x => ElementNodes % xyz(1:n,1) 491 ElementNodes % y => ElementNodes % xyz(1:n,2) 492 ElementNodes % z => ElementNodes % xyz(1:n,3) 493 sz = SIZE(ElementNodes % xyz,1) 494 SELECT CASE(Element % TYPE % DIMENSION) 495 CASE(1) 496 DO i=n+1,sz 497 ElementNodes % xyz(i,1) = REAL(0,dp) 498 END DO 499 DO i=1,sz 500 ElementNodes % xyz(i,2) = REAL(0,dp) 501 ElementNodes % xyz(i,3) = REAL(0,dp) 502 END DO 503 CASE(2) 504 DO i=n+1,sz 505 ElementNodes % xyz(i,1) = REAL(0,dp) 506 ElementNodes % xyz(i,2) = REAL(0,dp) 507 END DO 508 DO i=1,sz 509 ElementNodes % xyz(i,3) = REAL(0,dp) 510 END DO 511 CASE(3) 512 DO i=n+1,sz 513 ElementNodes % xyz(i,1) = REAL(0,dp) 514 ElementNodes % xyz(i,2) = REAL(0,dp) 515 ElementNodes % xyz(i,3) = REAL(0,dp) 516 END DO 517 CASE DEFAULT 518 CALL Fatal('GetReferenceElementNodes','Unsupported element dimension') 519 END SELECT 520 END IF 521 522 IF (isPElement(Element)) THEN 523 CALL GetRefPElementNodes(Element % Type, & 524 ElementNodes % x, & 525 ElementNodes % y, & 526 ElementNodes % z) 527 ELSE 528 IF (ALLOCATED(Element % Type % NodeU)) THEN 529 !DIR$ IVDEP 530 DO i=1,n 531 ElementNodes % x(i) = Element % Type % NodeU(i) 532 END DO 533 END IF 534 IF (ALLOCATED(Element % Type % NodeV)) THEN 535 !DIR$ IVDEP 536 DO i=1,n 537 ElementNodes % y(i) = Element % Type % NodeV(i) 538 END DO 539 END IF 540 IF (ALLOCATED(Element % Type % NodeW)) THEN 541 !DIR$ IVDEP 542 DO i=1,n 543 ElementNodes % z(i) = Element % Type % NodeW(i) 544 END DO 545 END IF 546 END IF 547 END SUBROUTINE GetReferenceElementNodes 548 549 SUBROUTINE PrintTestData(Element, ngp, t_n1, evals1, t_n2, evals2) 550 551 IMPLICIT NONE 552 TYPE(Element_t) :: Element 553 INTEGER, INTENT(IN) :: ngp, evals1, evals2 554 REAL(kind=dp), INTENT(IN) :: t_n1, t_n2 555 556 WRITE (*,'(A,I0)') 'Element type=', Element % TYPE % ElementCode 557 IF (ASSOCIATED(Element % PDefs)) THEN 558 WRITE (*,'(A,I0)') 'Element polynomial degree=', Element % PDefs % P 559 END IF 560 WRITE (*,'(A,L1)') 'Active P element=', isActivePElement(Element) 561 WRITE (*,'(A,I0)') 'Element number of nodes=', Element % TYPE % NumberOfNodes 562 WRITE (*,'(A,I0)') 'Element number of nonnodal dofs=', GetElementNOFDOFs(Element) 563 WRITE (*,'(A,I0)') 'Element number of bubble dofs=', GetElementNOFBDOFs() 564 WRITE (*,'(A,I0)') 'Number of Gauss points=', ngp 565 WRITE (*,'(A,I0)') 'Nodal basis, number of local matrix evaluations=', evals1 566 WRITE (*,'(A,F12.9)') 'Nodal basis, local matrix assembly t(s):', t_n1 567 WRITE (*,'(A,F12.2)') 'Nodal basis, local matrix evaluations/sec:', evals1/t_n1 568 WRITE (*,'(A,I0)') 'Vector basis, number of local matrix evaluations=', evals2 569 WRITE (*,'(A,F12.9)') 'Vector basis, local matrix assembly t(s):', t_n2 570 WRITE (*,'(A,F12.2)') 'Vector basis, local matrix evaluations/sec:', evals2/t_n2 571 572 END SUBROUTINE PrintTestData 573 574 SUBROUTINE AllocateMeshAndPElement(Mesh, ElementCode, P, PElement) 575 IMPLICIT NONE 576 577 TYPE(Mesh_t), POINTER :: Mesh 578 INTEGER, INTENT(IN) :: ElementCode, P 579 TYPE(Element_t), POINTER :: PElement 580 INTEGER :: node, edge, face, astat 581 582 ! Allocate a mesh with a single element 583 Mesh => AllocateMesh() 584 ALLOCATE(Mesh % Elements(1), STAT=astat) 585 IF (astat /= 0) THEN 586 CALL Fatal('AllocateMeshAndElement','Allocation of mesh element failed') 587 END IF 588 Mesh % NumberOfBulkElements = 1 589 590 ! Construct P element 591 PElement => AllocateElement() 592 PElement % ElementIndex = 1 593 PElement % BodyId = 1 594 PElement % Type => GetElementType( ElementCode ) 595 CALL AllocatePDefinitions(PElement) 596 597 ! Add element node indexes to element 598 ALLOCATE(PElement % NodeIndexes(PElement % Type % NumberOfNodes), STAT=astat) 599 IF (astat /= 0) THEN 600 CALL Fatal('AllocateMeshAndElement','Allocation of node indices failed') 601 END IF 602 PElement % NodeIndexes(:) = [(node, node=1,PElement % Type % NumberOfNodes)] 603 604 IF (PElement % Type % Dimension > 1) THEN 605 ! Add element edge indexes to element 606 ALLOCATE(PElement % EdgeIndexes(PElement % Type % NumberOfEdges), STAT=astat) 607 IF (astat /= 0) THEN 608 CALL Fatal('AllocateMeshAndElement','Allocation of edge indices failed') 609 END IF 610 PElement % EdgeIndexes(:) = [(edge, edge=1,PElement % Type % NumberOfEdges)] 611 612 ! Add all element edges to mesh 613 ALLOCATE(Mesh % Edges(PElement % Type % NumberOfEdges), STAT=astat) 614 IF (astat /= 0) THEN 615 CALL Fatal('AllocateMeshAndElement','Allocation of mesh edges failed') 616 END IF 617 Mesh % NumberOfEdges = PElement % Type % NumberOfEdges 618 ! Mesh % MinEdgeDofs = HUGE(Mesh % MinEdgeDofs) 619 Mesh % MinEdgeDofs = 0 620 Mesh % MaxEdgeDofs = 0 621 DO edge=1,PElement % Type % NumberOfEdges 622 CALL InitializePElement(Mesh % Edges(edge)) 623 ! Allocate edge element (always 202, even in 3D) 624 Mesh % Edges(edge) % Type => GetElementType(202, .FALSE.) 625 CALL AllocatePDefinitions(Mesh % Edges(edge)) 626 Mesh % Edges(edge) % PDefs % P = P 627 Mesh % Edges(edge) % PDefs % isEdge = .TRUE. 628 Mesh % Edges(edge) % PDefs % GaussPoints = (P+1) ** Mesh % Edges(edge) % Type % DIMENSION 629 Mesh % Edges(edge) % PDefs % LocalNumber = edge 630 Mesh % Edges(edge) % BDOFs = GetBubbleDofs(Mesh % Edges(edge), P) 631 Mesh % MinEdgeDofs = MIN(Mesh % MinEdgeDofs, Mesh % Edges(edge) % BDOFs) 632 Mesh % MaxEdgeDofs = MAX(Mesh % MaxEdgeDofs, Mesh % Edges(edge) % BDOFs) 633 END DO 634 END IF 635 636 IF (PElement % Type % Dimension > 2) THEN 637 ! Add element face indexes to element 638 ALLOCATE(PElement % FaceIndexes(PElement % Type % NumberOfFaces)) 639 IF (astat /= 0) THEN 640 CALL Fatal('AllocateMeshAndElement','Allocation of face indices failed') 641 END IF 642 PElement % FaceIndexes(:) = [(face, face=1,PElement % Type % NumberOfFaces)] 643 644 ! Add element faces to mesh 645 ALLOCATE(Mesh % Faces(PElement % Type % NumberOfFaces), STAT=astat) 646 IF (astat /= 0) THEN 647 CALL Fatal('AllocateMeshAndElement','Allocation of mesh faces failed') 648 END IF 649 Mesh % NumberOfFaces = PElement % Type % NumberOfFaces 650 Mesh % MinFaceDofs = HUGE(Mesh % MinFaceDofs) 651 Mesh % MaxFaceDofs = 0 652 DO face=1,PElement % Type % NumberOfFaces 653 CALL InitializePElement(Mesh % Faces(face)) 654 ! Allocate edge element (always 303 or 404, depending on element) 655 SELECT CASE (ElementCode/100) 656 CASE(5) 657 Mesh % Faces(face) % Type => GetElementType(303, .FALSE.) 658 CASE(6) 659 IF ( face == 1 ) THEN 660 Mesh % Faces(Face) % Type => GetElementType( 404, .FALSE. ) 661 ELSE 662 Mesh % Faces(Face) % Type => GetElementType( 303, .FALSE. ) 663 END IF 664 CASE(7) 665 IF ( face <= 2 ) THEN 666 Mesh % Faces(Face) % Type => GetElementType( 303, .FALSE. ) 667 ELSE 668 Mesh % Faces(Face) % Type => GetElementType( 404, .FALSE. ) 669 END IF 670 CASE(8) 671 Mesh % Faces(Face) % Type => GetElementType( 404, .FALSE.) 672 CASE DEFAULT 673 CALL Fatal('AllocateMeshAndElement','Unknown element type') 674 END SELECT 675 CALL AllocatePDefinitions(Mesh % Faces(face)) 676 Mesh % Faces(face) % PDefs % P = P 677 Mesh % Faces(face) % PDefs % isEdge = .TRUE. 678 Mesh % Faces(face) % PDefs % GaussPoints = (P+1) ** Mesh % Faces(face) % Type % DIMENSION 679 Mesh % Faces(face) % PDefs % LocalNumber = face 680 Mesh % Faces(face) % BDOFs = GetBubbleDofs(Mesh % Faces(face), P) 681 Mesh % MinFaceDofs = MAX(Mesh % MinFaceDofs, Mesh % Faces(face) % BDOFs) 682 Mesh % MaxFaceDofs = MAX(Mesh % MaxFaceDofs, Mesh % Faces(face) % BDOFs) 683 END DO 684 END IF 685 686 PElement % BDofs = GetBubbleDofs( PElement, P ) 687 PElement % PDefs % P = P 688 IF (ElementCode == 504) THEN 689 PElement % PDefs % TetraType = 1 690 ELSE 691 PElement % PDefs % TetraType = 0 692 END IF 693 PElement % PDefs % isEdge = .FALSE. 694 ! Set max element dofs to mesh 695 Mesh % MaxElementDOFs = PElement % Type % NumberOfNodes + & 696 PElement % Type % NumberOfEdges * Mesh % MaxEdgeDOFs + & 697 PElement % Type % NumberOfFaces * Mesh % MaxFaceDOFs + & 698 PElement % BDOFs 699 700 PElement % PDefs % GaussPoints = (P+1) ** PElement % Type % DIMENSION 701 END SUBROUTINE AllocateMeshAndPElement 702 703 SUBROUTINE InitializePElement(Element) 704 IMPLICIT NONE 705 TYPE(Element_t) :: Element 706 707 Element % BDOFs = 0 708 Element % NDOFs = 0 709 Element % BodyId = 1 710 Element % Splitted = 0 711 Element % hK = 0 712 Element % ElementIndex = 0 713 Element % StabilizationMk = 0 714 NULLIFY( Element % TYPE ) 715 NULLIFY( Element % PDefs ) 716 NULLIFY( Element % BubbleIndexes ) 717 NULLIFY( Element % DGIndexes ) 718 NULLIFY( Element % NodeIndexes ) 719 NULLIFY( Element % EdgeIndexes ) 720 NULLIFY( Element % FaceIndexes ) 721 NULLIFY( Element % BoundaryInfo ) 722 END SUBROUTINE InitializePElement 723 724 FUNCTION ClonePElement(Element) RESULT(ClonedElement) 725 IMPLICIT NONE 726 TYPE(Element_t), POINTER :: Element, ClonedElement 727 728 ALLOCATE(ClonedElement) 729 CALL InitializePElement(ClonedElement) 730 ClonedElement % BDOFs = Element % BDOFs 731 ClonedElement % NDOFs = Element % NDOFs 732 ClonedElement % BodyId = Element % BodyId 733 ClonedElement % Splitted = Element % Splitted 734 ClonedElement % hK = Element % hK 735 ClonedElement % ElementIndex = Element % ElementIndex 736 ClonedElement % StabilizationMk = Element % StabilizationMk 737 ClonedElement % Type => Element % Type 738 IF (ASSOCIATED(Element % PDefs)) THEN 739 CALL AllocatePDefinitions(ClonedElement) 740 ClonedElement % PDefs % P = Element % PDefs % P 741 ClonedElement % PDefs % TetraType = Element % PDefs % TetraType 742 ClonedElement % PDefs % isEdge = Element % PDefs % isEdge 743 ClonedElement % PDefs % GaussPoints = Element % PDefs % GaussPoints 744 ClonedElement % PDefs % pyramidQuadEdge = Element % PDefs % pyramidQuadEdge 745 ClonedElement % PDefs % localNumber = Element % PDefs % localNumber 746 END IF 747 IF (ASSOCIATED( Element % NodeIndexes )) THEN 748 ALLOCATE(ClonedElement % NodeIndexes( SIZE(Element % NodeIndexes))) 749 ClonedElement % NodeIndexes(:) = Element % NodeIndexes(:) 750 END IF 751 IF (ASSOCIATED( Element % EdgeIndexes )) THEN 752 ALLOCATE(ClonedElement % EdgeIndexes( SIZE(Element % EdgeIndexes))) 753 ClonedElement % EdgeIndexes(:) = Element % EdgeIndexes(:) 754 END IF 755 IF (ASSOCIATED( Element % FaceIndexes )) THEN 756 ALLOCATE(ClonedElement % FaceIndexes( SIZE(Element % FaceIndexes))) 757 ClonedElement % FaceIndexes(:) = Element % FaceIndexes(:) 758 END IF 759 END FUNCTION ClonePElement 760 761 SUBROUTINE DeallocatePElement(PElement) 762 IMPLICIT NONE 763 764 TYPE(Element_t), POINTER :: PElement 765 766 IF (ASSOCIATED(PElement % PDefs)) DEALLOCATE(PElement % PDefs) 767 IF (ASSOCIATED(PElement % NodeIndexes)) DEALLOCATE(PElement % NodeIndexes) 768 IF (ASSOCIATED(PElement % EdgeIndexes)) DEALLOCATE(PElement % EdgeIndexes) 769 IF (ASSOCIATED(PElement % FaceIndexes)) DEALLOCATE(PElement % FaceIndexes) 770 DEALLOCATE(PElement) 771 END SUBROUTINE DeallocatePElement 772 773 SUBROUTINE DeallocateTemporaryMesh(Mesh) 774 IMPLICIT NONE 775 776 TYPE(Mesh_t), POINTER :: Mesh 777 INTEGER :: edge, face 778 779 IF (ASSOCIATED(Mesh % Elements)) DEALLOCATE(Mesh % Elements) 780 DO edge=1,Mesh % NumberOfEdges 781 IF (ASSOCIATED(Mesh % Edges(edge) % PDefs)) DEALLOCATE(Mesh % Edges(edge) % PDefs) 782 END DO 783 IF (ASSOCIATED(Mesh % Edges)) DEALLOCATE(Mesh % Edges) 784 DO face=1,Mesh % NumberOfFaces 785 IF (ASSOCIATED(Mesh % Faces(face) % PDefs)) DEALLOCATE(Mesh % Faces(face) % PDefs) 786 END DO 787 IF (ASSOCIATED(Mesh % Faces)) DEALLOCATE(Mesh % Faces) 788 789 DEALLOCATE( Mesh % Nodes ) 790 DEALLOCATE( Mesh ) 791 END SUBROUTINE DeallocateTemporaryMesh 792 793 ! Portable wall-clock timer 794 FUNCTION ftimer() RESULT(timerval) 795 IMPLICIT NONE 796 797 REAL(KIND=dp) :: timerval 798 INTEGER(KIND=8) :: t, rate 799 800#ifdef _OPENMP 801 timerval = OMP_GET_WTIME() 802#else 803 CALL SYSTEM_CLOCK(t,count_rate=rate) 804 timerval = REAL(t,dp)/rate 805#endif 806 END FUNCTION ftimer 807 808!------------------------------------------------------------------------------ 809END SUBROUTINE LinearFormsAssembly 810!------------------------------------------------------------------------------ 811