1!/*****************************************************************************/ 2! * 3! * Elmer, A Finite Element Software for Multiphysical Problems 4! * 5! * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland 6! * 7! * This library is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU Lesser General Public 9! * License as published by the Free Software Foundation; either 10! * version 2.1 of the License, or (at your option) any later version. 11! * 12! * This library is distributed in the hope that it will be useful, 13! * but WITHOUT ANY WARRANTY; without even the implied warranty of 14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 15! * Lesser General Public License for more details. 16! * 17! * You should have received a copy of the GNU Lesser General Public 18! * License along with this library (in file ../LGPL-2.1); if not, write 19! * to the Free Software Foundation, Inc., 51 Franklin Street, 20! * Fifth Floor, Boston, MA 02110-1301 USA 21! * 22! *****************************************************************************/ 23! 24!/****************************************************************************** 25! * 26! * Authors: Mikko Byckling, Juha Ruokolainen 27! * Email: Juha.Ruokolainen@csc.fi 28! * Web: http://www.csc.fi/elmer 29! * Address: CSC - IT Center for Science Ltd. 30! * Keilaranta 14 31! * 02101 Espoo, Finland 32! * 33! * Original Date: 23 Aug 2004 34! * 35! *****************************************************************************/ 36 37!> \ingroup ElmerLib 38!> \{ 39 40!------------------------------------------------------------------------------- 41!> Module defining mappings for p elements. These include nodal points 42!> contained by faces and edges, element boundary maps (edges for 2d elements, 43!> faces for 3d) and mappings from faces to edge numbers. Mappings defined in 44!> this module are compatible with basis functions defined in module 45!> PElementBase. 46!------------------------------------------------------------------------------- 47 48MODULE PElementMaps 49 USE Types 50 51 IMPLICIT NONE 52 53 ! Private mappings. For access use get[Element][Type]Map(i) 54 PRIVATE QuadEdgeMap, TriangleEdgeMap, & 55 TetraEdgeMap1, TetraFaceMap1, TetraFaceEdgeMap1, & 56 TetraEdgeMap2, TetraFaceMap2, TetraFaceEdgeMap2, & 57 BrickEdgeMap, BrickFaceMap, BrickFaceEdgeMap, & 58 WedgeEdgeMap, WedgeFaceMap, WedgeFaceEdgeMap, & 59 PyramidEdgeMap, PyramidFaceMap, PyramidFaceEdgeMap, & 60 MInit 61 ! Mappings 62 INTEGER, TARGET, SAVE :: QuadEdgeMap(4,2), TriangleEdgeMap(3,2), & 63 TetraEdgeMap1(6,2), TetraFaceMap1(4,3), TetraFaceEdgeMap1(4,3), & 64 TetraEdgeMap2(6,2), TetraFaceMap2(4,3), TetraFaceEdgeMap2(4,3),& 65 BrickEdgeMap(12,2), BrickFaceMap(6,4), BrickFaceEdgeMap(6,4), & 66 WedgeEdgeMap(9,2), WedgeFaceMap(5,4), WedgeFaceEdgeMap(5,4), & 67 PyramidEdgeMap(8,2), PyramidFaceMap(5,4), PyramidFaceEdgeMap(5,4) 68 69 LOGICAL, SAVE :: MInit = .FALSE. 70 !$OMP THREADPRIVATE(MInit, QuadEdgeMap, TriangleEdgeMap, & 71 !$OMP& TetraEdgeMap1, TetraFaceMap1, TetraFaceEdgeMap1, & 72 !$OMP& TetraEdgeMap2, TetraFaceMap2, TetraFaceEdgeMap2,& 73 !$OMP& BrickEdgeMap, BrickFaceMap, BrickFaceEdgeMap, & 74 !$OMP& WedgeEdgeMap, WedgeFaceMap, WedgeFaceEdgeMap, & 75 !$OMP& PyramidEdgeMap, PyramidFaceMap, PyramidFaceEdgeMap) 76CONTAINS 77 78 ! MAPPINGS 79 80 ! First some direct mappings to elements. These should not be used directly 81 ! unless element type is implicitly known from context. Better way is to use 82 ! getElement[Boundary,Edge,Face]Map -routines. 83 84 ! Call: localEdge = getQuadEdge(i) 85 ! 86 ! Function returns mapping from edge number to edge endpoints 87 88 FUNCTION getQuadEdgeMap(i) RESULT(localEdge) 89 IMPLICIT NONE 90 91 INTEGER, INTENT(IN) :: i 92 INTEGER, DIMENSION(2) :: localEdge 93 94 IF (.NOT. MInit) CALL InitializeMappings() 95 96 localEdge(:) = QuadEdgeMap(i,:) 97 END FUNCTION getQuadEdgeMap 98 99 ! Call: localEdge = getTriangleEdge(i) 100 ! 101 ! Function returns mapping from edge number to edge endpoints 102 103 FUNCTION getTriangleEdgeMap(i) RESULT(localEdge) 104 IMPLICIT NONE 105 106 INTEGER, INTENT(IN) :: i 107 INTEGER, DIMENSION(2) :: localEdge 108 109 IF (.NOT. MInit) CALL InitializeMappings() 110 111 localEdge(:)=TriangleEdgeMap(i,:) 112 END FUNCTION getTriangleEdgeMap 113 114 ! Call: localEdge = getBrickEdgeMap(i) 115 ! 116 ! Function returns mapping from edge number to edge endpoints 117 118 FUNCTION getBrickEdgeMap(i) RESULT(localEdge) 119 IMPLICIT NONE 120 121 INTEGER, INTENT(IN) :: i 122 INTEGER, DIMENSION(2) :: localEdge 123 124 IF (.NOT. MInit) CALL InitializeMappings() 125 126 localEdge(:) = BrickEdgeMap(i,:) 127 END FUNCTION getBrickEdgeMap 128 129 ! Call: localFace = getBrickFaceMap(i) 130 ! 131 ! Function returns mapping from face number to face nodes 132 133 FUNCTION getBrickFaceMap(i) RESULT(localFace) 134 IMPLICIT NONE 135 136 INTEGER, INTENT(IN) :: i 137 INTEGER, DIMENSION(4) :: localFace 138 139 IF (.NOT. MInit) CALL InitializeMappings() 140 141 localFace(:) = BrickFaceMap(i,:) 142 END FUNCTION getBrickFaceMap 143 144 ! Call: localEdge = getFaceEdgeMap(face, localNode) 145 ! 146 ! getFaceEdgeMap returns number of local edge when given face and 147 ! its local node number. Node number is treated as edges beginning point 148 149 FUNCTION getBrickFaceEdgeMap(face, localNode) RESULT(localEdge) 150 IMPLICIT NONE 151 CHARACTER(LEN=MAX_NAME_LEN) :: msg 152 153 ! Parameters 154 INTEGER, INTENT(IN) :: face, localNode 155 ! Variables 156 INTEGER :: localEdge 157 158 IF (.NOT. MInit) CALL InitializeMappings() 159 160 localEdge = BrickFaceEdgeMap(face,localNode) 161 162 IF (localEdge == 0) THEN 163 WRITE (msg,'(A,I2,I3)') 'Unknown combination node for (face,node)', face,localNode 164 CALL Fatal('getBrickFaceEdgeMap', msg) 165 END IF 166 END FUNCTION getBrickFaceEdgeMap 167 168 169 FUNCTION getTetraEdgeMap(i,TYPE) RESULT(edge) 170 IMPLICIT NONE 171 172 INTEGER, INTENT(IN) :: i 173 INTEGER, INTENT(IN), OPTIONAL :: TYPE 174 INTEGER :: t 175 INTEGER, DIMENSION(2) :: edge 176 177 IF (.NOT. MInit) CALL InitializeMappings() 178 179 ! If type not present use default (1) 180 t = 1 181 IF (PRESENT(TYPE)) t = TYPE 182 183 ! Select edge map by tetra type 184 SELECT CASE (t) 185 CASE (1) 186 edge(:) = TetraEdgeMap1(i,:) 187 CASE (2) 188 edge(:) = TetraEdgeMap2(i,:) 189 CASE DEFAULT 190 CALL Fatal('PElementMaps::getTetraEdgeMap','Unknown tetra type') 191 END SELECT 192 END FUNCTION getTetraEdgeMap 193 194 195 FUNCTION getTetraFaceMap(i,TYPE) RESULT(face) 196 IMPLICIT NONE 197 198 INTEGER, INTENT(IN) :: i 199 INTEGER, INTENT(IN), OPTIONAL :: TYPE 200 INTEGER :: t 201 INTEGER, DIMENSION(3) :: face 202 203 IF (.NOT. MInit) CALL InitializeMappings() 204 205 ! If type not present use default (1) 206 t = 1 207 IF (PRESENT(TYPE)) t = TYPE 208 209 ! Select face map by tetra type 210 SELECT CASE(t) 211 CASE (1) 212 face(:) = TetraFaceMap1(i,:) 213 CASE (2) 214 face(:) = TetraFaceMap2(i,:) 215 CASE DEFAULT 216 CALL Fatal('PElementMaps::getTetraFaceMap','Unknown tetra type') 217 END SELECT 218 END FUNCTION getTetraFaceMap 219 220 FUNCTION getWedgeEdgeMap(i) RESULT(edge) 221 IMPLICIT NONE 222 223 INTEGER, INTENT(IN) :: i 224 INTEGER, DIMENSION(2) :: edge 225 226 IF (.NOT. MInit) CALL InitializeMappings() 227 228 edge(:) = WedgeEdgeMap(i,:) 229 END FUNCTION getWedgeEdgeMap 230 231 232 FUNCTION getWedgeFaceMap(i) RESULT(face) 233 IMPLICIT NONE 234 235 INTEGER, INTENT(IN) :: i 236 INTEGER, DIMENSION(4) :: face 237 238 IF (.NOT. MInit) CALL InitializeMappings() 239 240 face(:) = WedgeFaceMap(i,:) 241 END FUNCTION getWedgeFaceMap 242 243 244 FUNCTION getPyramidEdgeMap(i) RESULT(edge) 245 IMPLICIT NONE 246 247 INTEGER, INTENT(IN) :: i 248 INTEGER, DIMENSION(2) :: edge 249 250 IF (.NOT. MInit) CALL InitializeMappings() 251 252 edge(:) = PyramidEdgeMap(i,:) 253 END FUNCTION getPyramidEdgeMap 254 255 256 FUNCTION getPyramidFaceMap(i) RESULT(face) 257 IMPLICIT NONE 258 259 INTEGER, INTENT(IN) :: i 260 INTEGER, DIMENSION(4) :: face 261 262 IF (.NOT. MInit) CALL InitializeMappings() 263 264 face(:) = PyramidFaceMap(i,:) 265 END FUNCTION getPyramidFaceMap 266 267 268!------------------------------------------------------------------------------ 269!> Mapping from element local edge or face number to nodes contained in 270!> that edge or face. 271!------------------------------------------------------------------------------ 272 FUNCTION getElementBoundaryMap(Element, i) RESULT(map) 273!------------------------------------------------------------------------------ 274! 275! ARGUMENTS: 276! Type(Element_t) :: Element 277! INPUT: Element to get map for 278! 279! INTEGER, INTENT(IN) :: i 280! INPUT: Local number of elements edge or face 281! 282! FUNCTION VALUE: 283! INTEGER :: map(4) 284! Map containing local node numbers of given local edge or face 285! 286!------------------------------------------------------------------------------ 287 IMPLICIT NONE 288 289 TYPE(Element_t) :: Element 290 INTEGER, INTENT(IN) :: i 291 292 INTEGER :: map(4) 293 294 IF (.NOT. MInit) CALL InitializeMappings() 295 map = 0 296 297 ! Function is not defined for non p elements 298 !IF (.NOT. ASSOCIATED(Element % PDefs)) THEN 299 ! CALL Warn('PElementMaps::getElementBoundaryMap','Element not p element') 300 ! RETURN 301 !END IF 302 303 SELECT CASE(Element % TYPE % ElementCode / 100) 304 CASE (3) 305 map(1:2) = getTriangleEdgeMap(i) 306 CASE (4) 307 map(1:2) = getQuadEdgeMap(i) 308 CASE (5) 309 map(1:3) = getTetraFaceMap(i,Element % PDefs % TetraType) 310 CASE (6) 311 map(1:4) = getPyramidFaceMap(i) 312 CASE (7) 313 map(1:4) = getWedgeFaceMap(i) 314 CASE (8) 315 map(1:4) = getBrickFaceMap(i) 316 CASE DEFAULT 317 CALL Fatal('PElementMaps::getElementBoundaryMap','Unsupported element type') 318 END SELECT 319 END FUNCTION getElementBoundaryMap 320 321 322!------------------------------------------------------------------------------ 323!> Mapping from element local face to local edges contained in face. Given 324!> element and local face number this routine returns numbers of local edges 325!> on face. 326!------------------------------------------------------------------------------ 327 FUNCTION getFaceEdgeMap( Element, i) RESULT(map) 328!------------------------------------------------------------------------------ 329! 330! ARGUMENTS: 331! Type(Element_t) :: Element 332! INPUT: Element to get map for 333! 334! INTEGER, INTENT(IN) :: i 335! INPUT: Local number of element face 336! 337! FUNCTION VALUE: 338! INTEGER :: map(4) 339! Map containing local numbers of edges on face 340! 341!------------------------------------------------------------------------------ 342 IMPLICIT NONE 343 344 TYPE(Element_t) :: Element 345 INTEGER, INTENT(IN) :: i 346 347 INTEGER :: elementCode, map(4) 348 349 elementCode = Element % TYPE % ElementCode 350 351 IF (.NOT. MInit) CALL InitializeMappings() 352 353 ! Function is not defined for non p elements 354 !IF (.NOT. ASSOCIATED(Element % PDefs)) THEN 355 ! CALL Warn('PElementMaps::getFaceEdgeMap','Element not p element') 356 ! map = 0 357 ! RETURN 358 !END IF 359 360 SELECT CASE(elementCode / 100) 361 CASE (5) 362 map = 0 363 SELECT CASE (Element % PDefs % TetraType) 364 CASE (1) 365 map(1:3) = TetraFaceEdgeMap1(i,:) 366 CASE (2) 367 map(1:3) = TetraFaceEdgeMap2(i,:) 368 CASE DEFAULT 369 CALL Fatal('PElementMaps::getFaceEdgeMap','Unknown tetra type') 370 END SELECT 371 CASE (6) 372 map(1:4) = PyramidFaceEdgeMap(i,:) 373 CASE (7) 374 map(1:4) = WedgeFaceEdgeMap(i,:) 375 CASE (8) 376 map(1:4) = BrickFaceEdgeMap(i,:) 377 CASE DEFAULT 378 CALL Fatal('PElementMaps::getFaceEdgeMap','Unsupported element type') 379 END SELECT 380 END FUNCTION getFaceEdgeMap 381 382 383!------------------------------------------------------------------------------ 384!> Get mappings for given element to element edges and their nodes. Given 385!> element, this routine returns a map containing nodes (endpoints) of 386!> elements edges. 387!------------------------------------------------------------------------------ 388 SUBROUTINE GetElementEdgeMap( Element, map ) 389!------------------------------------------------------------------------------ 390! 391! ARGUMENTS: 392! Type(Element_t) :: Element 393! INPUT: Element to get map for 394! 395! INTEGER :: map(:,:) 396! OUTPUT: Map containing local node numbers of local edges 397! 398!------------------------------------------------------------------------------ 399 IMPLICIT NONE 400 TYPE(Element_t) :: Element 401 INTEGER, POINTER :: map(:,:) 402 403 IF (.NOT. MInit) CALL InitializeMappings() 404 405 ! Function is not defined for non p elements 406 IF (.NOT. ASSOCIATED(Element % PDefs)) THEN 407 CALL Warn('PElementMaps::GetElementEdgeMap','Element not p element') 408 map = 0 409 RETURN 410 END IF 411 412 SELECT CASE (Element % TYPE % ElementCode / 100) 413 CASE (3) 414 map => TriangleEdgeMap 415 CASE (4) 416 map => QuadEdgeMap 417 CASE (5) 418 SELECT CASE( Element % PDefs % TetraType ) 419 CASE (1) 420 map => TetraEdgeMap1 421 CASE (2) 422 map => TetraEdgeMap2 423 CASE DEFAULT 424 CALL Fatal('PElementMaps::GetElementEdgeMap','Unknown tetra type for p element') 425 END SELECT 426 CASE (6) 427 map => PyramidEdgeMap 428 CASE (7) 429 map => WedgeEdgeMap 430 CASE (8) 431 map => BrickEdgeMap 432 CASE DEFAULT 433 CALL Fatal('PElementMaps::GetElementEdgeMap','Unsupported element type') 434 END SELECT 435 END SUBROUTINE GetElementEdgeMap 436 437 438!------------------------------------------------------------------------------ 439!> Get mappings for given element to element faces and their nodes. Given 440!> element, this routine returns a map containing nodes (endpoints) of 441!> elements face. 442!------------------------------------------------------------------------------ 443 SUBROUTINE GetElementFaceMap( Element, faceMap ) 444!------------------------------------------------------------------------------ 445! 446! ARGUMENTS: 447! Type(Element_t) :: Element 448! INPUT: Element to get map for 449! 450! INTEGER :: map(:,:) 451! OUTPUT: Map containing local node numbers of local faces 452! 453!------------------------------------------------------------------------------ 454 IMPLICIT NONE 455 456 TYPE(Element_t) :: Element 457 INTEGER, POINTER :: faceMap(:,:) 458 459 IF (.NOT. MInit) CALL InitializeMappings() 460 461 ! Function is not defined for non p elements 462 IF (.NOT. ASSOCIATED(Element % PDefs)) THEN 463 CALL Warn('PElementMaps::GetElementFaceMap','Element not p element') 464 NULLIFY(faceMap) 465 RETURN 466 END IF 467 468 SELECT CASE (Element % TYPE % ElementCode / 100) 469! CASE (3) 470! facemap => TriangleEdgeMap 471! CASE (4) 472! facemap => QuadEdgeMap 473 CASE (5) 474 SELECT CASE( Element % PDefs % TetraType ) 475 CASE (1) 476 faceMap => TetraFaceMap1 477 CASE (2) 478 faceMap => TetraFaceMap2 479 CASE DEFAULT 480 CALL Fatal('PElementMaps::GetElementFaceMap','Unknown tetra type for p element') 481 END SELECT 482 CASE (6) 483 faceMap => PyramidFaceMap 484 CASE (7) 485 faceMap => WedgeFaceMap 486 CASE (8) 487 faceMap => BrickFaceMap 488 CASE DEFAULT 489 CALL Fatal('PElementMaps::GetElementFaceMap','Unsupported element type') 490 END SELECT 491 END SUBROUTINE GetElementFaceMap 492 493 494!------------------------------------------------------------------------------ 495!> Get mappings for given element to elements faces and their edge. Given 496!> element, this routine returns a map containing local edge numbers of 497!> elements faces. 498!------------------------------------------------------------------------------ 499 SUBROUTINE GetElementFaceEdgeMap( Element, faceEdgeMap ) 500!------------------------------------------------------------------------------ 501! 502! ARGUMENTS: 503! Type(Element_t) :: Element 504! INPUT: Element to get map for 505! 506! INTEGER :: map(:,:) 507! OUTPUT: Map containing local edge numbers of local faces 508! 509!------------------------------------------------------------------------------ 510 IMPLICIT NONE 511 512 TYPE(Element_t) :: Element 513 INTEGER, POINTER :: faceEdgeMap(:,:) 514 515 IF (.NOT. MInit) CALL InitializeMappings() 516 517 ! Function is not defined for non p elements 518 IF (.NOT. ASSOCIATED(Element % PDefs)) THEN 519 CALL Warn('PElementMaps::GetElementFaceEdgeMap','Element not p element') 520 NULLIFY(faceEdgeMap) 521 RETURN 522 END IF 523 524 SELECT CASE (Element % TYPE % ElementCode / 100) 525 CASE (5) 526 SELECT CASE( Element % PDefs % TetraType ) 527 CASE (1) 528 faceEdgeMap => TetraFaceEdgeMap1 529 CASE (2) 530 faceEdgeMap => TetraFaceEdgeMap2 531 CASE DEFAULT 532 CALL Fatal('PElementMaps::GetElementFaceEdgeMap','Unknown tetra type for p element') 533 END SELECT 534 CASE (6) 535 faceEdgeMap => PyramidFaceEdgeMap 536 CASE (7) 537 faceEdgeMap => WedgeFaceEdgeMap 538 CASE (8) 539 faceEdgeMap => BrickFaceEdgeMap 540 CASE DEFAULT 541 CALL Fatal('PElementMaps::GetElementFaceEdgeMap','Unsupported element type') 542 END SELECT 543 END SUBROUTINE GetElementFaceEdgeMap 544 545 546!------------------------------------------------------------------------------ 547!> This subroutine initializes element mappings. 548!------------------------------------------------------------------------------ 549 SUBROUTINE InitializeMappings() 550!------------------------------------------------------------------------------ 551 IMPLICIT NONE 552 553 CALL Info('PElementMaps::InitializeMappings','Initializing mappings for elements',Level=10) 554 555 ! Quad edge mappings 556 QuadEdgeMap(1,:) = (/ 1,2 /) 557 QuadEdgeMap(2,:) = (/ 2,3 /) 558 QuadEdgeMap(3,:) = (/ 4,3 /) 559 QuadEdgeMap(4,:) = (/ 1,4 /) 560 561 ! Triangle edge mappings 562 TriangleEdgeMap(1,:) = (/ 1,2 /) 563 TriangleEdgeMap(2,:) = (/ 2,3 /) 564 TriangleEdgeMap(3,:) = (/ 3,1 /) 565 566 ! Brick edge mappings 567 BrickEdgeMap(1,:) = (/ 1,2 /) 568 BrickEdgeMap(2,:) = (/ 2,3 /) 569 BrickEdgeMap(3,:) = (/ 4,3 /) 570 BrickEdgeMap(4,:) = (/ 1,4 /) 571 BrickEdgeMap(5,:) = (/ 5,6 /) 572 BrickEdgeMap(6,:) = (/ 6,7 /) 573 BrickEdgeMap(7,:) = (/ 8,7 /) 574 BrickEdgeMap(8,:) = (/ 5,8 /) 575 BrickEdgeMap(9,:) = (/ 1,5 /) 576 BrickEdgeMap(10,:) = (/ 2,6 /) 577 BrickEdgeMap(11,:) = (/ 3,7 /) 578 BrickEdgeMap(12,:) = (/ 4,8 /) 579 580 ! Brick face mappings 581 BrickFaceMap(1,:) = (/ 1,2,3,4 /) ! xi,eta 582 BrickFaceMap(2,:) = (/ 5,6,7,8 /) ! xi,eta 583 BrickFaceMap(3,:) = (/ 1,2,6,5 /) ! xi,zeta 584 BrickFaceMap(4,:) = (/ 2,3,7,6 /) ! eta,zeta 585 ! BrickFaceMap(5,:) = (/ 3,4,8,7 /) 586 BrickFaceMap(5,:) = (/ 4,3,7,8 /) 587 ! BrickFaceMap(6,:) = (/ 4,1,5,8 /) 588 BrickFaceMap(6,:) = (/ 1,4,8,5 /) 589 590 BrickFaceEdgeMap(1,:) = (/ 1,2,3,4 /) 591 BrickFaceEdgeMap(2,:) = (/ 5,6,7,8 /) 592 BrickFaceEdgeMap(3,:) = (/ 1,10,5,9 /) 593 BrickFaceEdgeMap(4,:) = (/ 2,11,6,10 /) 594 ! BrickFaceEdgeMap(5,:) = (/ 3,12,7,11 /) 595 BrickFaceEdgeMap(5,:) = (/ 3,11,7,12 /) 596 ! BrickFaceEdgeMap(6,:) = (/ 4,9,8,12 /) 597 BrickFaceEdgeMap(6,:) = (/ 4,12,8,9 /) 598 599 ! Tetra edge mappings (not needed for enforcing parity!) 600 ! Type 1 601 TetraEdgeMap1(1,:) = (/ 1,2 /) 602 TetraEdgeMap1(2,:) = (/ 2,3 /) 603 TetraEdgeMap1(3,:) = (/ 1,3 /) 604 TetraEdgeMap1(4,:) = (/ 1,4 /) 605 TetraEdgeMap1(5,:) = (/ 2,4 /) 606 TetraEdgeMap1(6,:) = (/ 3,4 /) 607 ! Type 2 608 TetraEdgeMap2(1,:) = (/ 1,2 /) 609 TetraEdgeMap2(2,:) = (/ 3,2 /) 610 TetraEdgeMap2(3,:) = (/ 1,3 /) 611 TetraEdgeMap2(4,:) = (/ 1,4 /) 612 TetraEdgeMap2(5,:) = (/ 2,4 /) 613 TetraEdgeMap2(6,:) = (/ 3,4 /) 614 615 ! Tetra face mappings (not needed for enforcing parity!) 616 ! Type 1 617 TetraFaceMap1(1,:) = (/ 1,2,3 /) 618 TetraFaceMap1(2,:) = (/ 1,2,4 /) 619 TetraFaceMap1(3,:) = (/ 2,3,4 /) 620 TetraFaceMap1(4,:) = (/ 1,3,4 /) 621 ! Type 2 622 TetraFaceMap2(1,:) = (/ 1,3,2 /) 623 TetraFaceMap2(2,:) = (/ 1,2,4 /) 624 TetraFaceMap2(3,:) = (/ 3,2,4 /) 625 TetraFaceMap2(4,:) = (/ 1,3,4 /) 626 627 ! Type 1 628 TetraFaceEdgeMap1(1,:) = (/ 1,2,3 /) 629 TetraFaceEdgeMap1(2,:) = (/ 1,5,4 /) 630 TetraFaceEdgeMap1(3,:) = (/ 2,6,5 /) 631 TetraFaceEdgeMap1(4,:) = (/ 3,6,4 /) 632 ! Type 2 633 TetraFaceEdgeMap2(1,:) = (/ 3,2,1 /) 634 TetraFaceEdgeMap2(2,:) = (/ 1,5,4 /) 635 TetraFaceEdgeMap2(3,:) = (/ 2,5,6 /) 636 TetraFaceEdgeMap2(4,:) = (/ 3,6,4 /) 637 638 ! Wedge edge mappings 639 WedgeEdgeMap(1,:) = (/ 1,2 /) 640 WedgeEdgeMap(2,:) = (/ 2,3 /) 641 WedgeEdgeMap(3,:) = (/ 3,1 /) 642 WedgeEdgeMap(4,:) = (/ 4,5 /) 643 WedgeEdgeMap(5,:) = (/ 5,6 /) 644 WedgeEdgeMap(6,:) = (/ 6,4 /) 645 WedgeEdgeMap(7,:) = (/ 1,4 /) 646 WedgeEdgeMap(8,:) = (/ 2,5 /) 647 WedgeEdgeMap(9,:) = (/ 3,6 /) 648 649 ! Wedge face mappings 650 WedgeFaceMap(1,:) = (/ 1,2,3,0 /) 651 WedgeFaceMap(2,:) = (/ 4,5,6,0 /) 652 WedgeFaceMap(3,:) = (/ 1,2,5,4 /) 653 WedgeFaceMap(4,:) = (/ 2,3,6,5 /) 654 WedgeFaceMap(5,:) = (/ 3,1,4,6 /) 655 656 WedgeFaceEdgeMap(1,:) = (/ 1,2,3,0 /) 657 WedgeFaceEdgeMap(2,:) = (/ 4,5,6,0 /) 658 WedgeFaceEdgeMap(3,:) = (/ 1,8,4,7 /) 659 WedgeFaceEdgeMap(4,:) = (/ 2,9,5,8 /) 660 WedgeFaceEdgeMap(5,:) = (/ 3,7,6,9 /) 661 662 ! Pyramid edge mappings 663 PyramidEdgeMap(1,:) = (/ 1,2 /) 664 PyramidEdgeMap(2,:) = (/ 2,3 /) 665 PyramidEdgeMap(3,:) = (/ 4,3 /) 666 PyramidEdgeMap(4,:) = (/ 1,4 /) 667 PyramidEdgeMap(5,:) = (/ 1,5 /) 668 PyramidEdgeMap(6,:) = (/ 2,5 /) 669 PyramidEdgeMap(7,:) = (/ 3,5 /) 670 PyramidEdgeMap(8,:) = (/ 4,5 /) 671 672 ! Pyramid face mappings 673 PyramidFaceMap(1,:) = (/ 1,2,3,4 /) 674 PyramidFaceMap(2,:) = (/ 1,2,5,0 /) 675 PyramidFaceMap(3,:) = (/ 2,3,5,0 /) 676 PyramidFaceMap(4,:) = (/ 3,4,5,0 /) 677 PyramidFaceMap(5,:) = (/ 4,1,5,0 /) 678 679 PyramidFaceEdgeMap(1,:) = (/ 1,2,3,4 /) 680 PyramidFaceEdgeMap(2,:) = (/ 1,6,5,0 /) 681 PyramidFaceEdgeMap(3,:) = (/ 2,7,6,0 /) 682 PyramidFaceEdgeMap(4,:) = (/ 3,8,7,0 /) 683 PyramidFaceEdgeMap(5,:) = (/ 4,5,8,0 /) 684 685 MInit = .TRUE. 686 END SUBROUTINE InitializeMappings 687 688!------------------------------------------------------------------------------ 689 FUNCTION getEdgeDOFs( Element, p ) RESULT(EdgeDOFs) 690!------------------------------------------------------------------------------ 691 IMPLICIT NONE 692!------------------------------------------------------------------------------ 693 TYPE(Element_t) :: Element 694 INTEGER :: EdgeDOFs 695 INTEGER, INTENT(IN) :: p 696 697 IF (.NOT. ASSOCIATED(Element % PDefs) ) THEN 698 EdgeDOFs = 0 699 RETURN 700 END IF 701 702 EdgeDOFs = MAX(0, p-1) 703!------------------------------------------------------------------------------ 704 END FUNCTION getEdgeDOFs 705!------------------------------------------------------------------------------ 706 707 708!------------------------------------------------------------------------------ 709!> Based on element face polynomial degree p, return degrees of freedom for 710!> given face. 711!------------------------------------------------------------------------------ 712 FUNCTION getFaceDOFs( Element, p, faceNumber ) RESULT(faceDOFs) 713!------------------------------------------------------------------------------ 714! 715! ARGUMENTS: 716! Type(Element_t), POINTER :: Element 717! INPUT: Element to get face dofs to 718! 719! INTEGER :: p 720! INPUT: Face polynomial degree p 721! 722! INTEGER :: faceNumber 723! INPUT: Local number of face for element (important for wedges and 724! pyramids). 725! 726! FUNCTION VALUE: 727! REAL(KIND=dp) :: faceDOFs 728! number of face dofs for Element 729! 730!------------------------------------------------------------------------------ 731 IMPLICIT NONE 732 733 TYPE(Element_t) :: Element 734 INTEGER, INTENT(IN) :: p 735 INTEGER, INTENT(IN), OPTIONAL :: faceNumber 736 INTEGER :: faceDOFs 737 738 ! This function is not defined for non p elements 739 IF (.NOT. ASSOCIATED(Element % PDefs) ) THEN 740 faceDOFs = 0 741 RETURN 742 END IF 743 744 faceDOFs = 0 745 SELECT CASE(Element % TYPE % ElementCode / 100) 746 ! Tetrahedron 747 CASE (5) 748 IF (p >= 3) faceDOFs = (p-1)*(p-2)/2 749 ! Pyramid 750 CASE (6) 751 SELECT CASE(faceNumber) 752 CASE (1) 753 IF (p >= 4) faceDOFs = (p-2)*(p-3)/2 754 CASE (2:5) 755 IF (p >= 3) faceDOFs = (p-1)*(p-2)/2 756 END SELECT 757 ! Wedge 758 CASE (7) 759 SELECT CASE(faceNumber) 760 CASE (1,2) 761 IF (p >= 3) faceDOFs = (p-1)*(p-2)/2 762 CASE (3:5) 763 IF (p >= 4) faceDOFs = (p-2)*(p-3)/2 764 END SELECT 765 ! Brick 766 CASE (8) 767 IF (p >= 4) faceDOFs = (p-2)*(p-3)/2 768 CASE DEFAULT 769 CALL Warn('MeshUtils::getFaceDOFs','Unsupported p element type') 770 faceDOFs = p 771 END SELECT 772 773 faceDOFs = MAX(0, faceDOFs) 774 END FUNCTION getFaceDOFs 775 776 777!------------------------------------------------------------------------------ 778!> Based on the polynomial degree p of the element, return the number of 779!> bubble functions (the count of bubble DOFs). 780!> NOTE: The returned value is not the bubble count for an approximation 781!> based on the space Q_p of polynomials of degree at most p in each variable 782!> separately. 783!------------------------------------------------------------------------------ 784 FUNCTION getBubbleDOFs( Element, p) RESULT(bubbleDOFs) 785!------------------------------------------------------------------------------ 786! 787! ARGUMENTS: 788! Type(Element_t), POINTER :: Element 789! INPUT: Element to get bubble dofs to 790! 791! INTEGER :: p 792! INPUT: Element polynomial degree p 793! 794! FUNCTION VALUE: 795! REAL(KIND=dp) :: bubbleDOFs 796! number of bubble dofs for Element 797! 798!------------------------------------------------------------------------------ 799 IMPLICIT NONE 800 801 TYPE(Element_t) :: Element 802 INTEGER, INTENT(IN) :: p 803 INTEGER :: bubbleDOFs 804 805 ! This function is not defined for non p elements 806 IF (.NOT. ASSOCIATED(Element % PDefs) ) THEN 807 bubbleDOFs = 0 808 RETURN 809 END IF 810 811 ! Select by element type 812 bubbleDOFs = 0 813 SELECT CASE (Element % TYPE % ElementCode / 100) 814 ! Line 815 CASE (2) 816 IF (p >= 2) bubbleDOFs = p - 1 817 ! Triangle 818 CASE (3) 819 IF (p >= 3) bubbleDOFs = (p-1)*(p-2)/2 820 ! Quad 821 CASE (4) 822 IF (p >= 4) bubbleDOFs = (p-2)*(p-3)/2 823 ! Tetrahedron 824 CASE (5) 825 IF (p >= 4) bubbleDOFs = (p-1)*(p-2)*(p-3)/6 826 ! Pyramid 827 CASE (6) 828 IF (p >= 4) bubbleDOFs = (p-1)*(p-2)*(p-3)/6 829 ! Wedge 830 CASE (7) 831 IF (p >= 5) bubbleDOFs = (p-2)*(p-3)*(p-4)/6 832 ! Brick 833 CASE (8) 834 IF (p >= 6) bubbleDOFs = (p-3)*(p-4)*(p-5)/6 835 CASE DEFAULT 836 CALL Warn('MeshUtils::getBubbleDOFs','Unsupported p element type') 837 bubbleDOFs = p 838 END SELECT 839 840 bubbleDOFs = MAX(0, bubbleDOFs) 841!------------------------------------------------------------------------------ 842 END FUNCTION getBubbleDOFs 843!------------------------------------------------------------------------------ 844 845 846!------------------------------------------------------------------------------ 847 FUNCTION isActivePElement(Element) RESULT(retVal) 848!------------------------------------------------------------------------------ 849 IMPLICIT NONE 850 851 TYPE(Element_t), INTENT(IN) :: Element 852 INTEGER :: c 853 LOGICAL :: retVal 854 855 retVal = isPelement(Element) 856 857 IF(ASSOCIATED(CurrentModel % Solver)) THEN 858 IF(ALLOCATED(CurrentModel % Solver % Def_Dofs)) THEN 859 c = Element % Type % ElementCode / 100 860 retVal = retVal.AND.ANY(CurrentModel % Solver % Def_Dofs(c,:,6)>0) 861 END IF 862 END IF 863 864!------------------------------------------------------------------------------ 865 END FUNCTION isActivePElement 866!------------------------------------------------------------------------------ 867 868 869!------------------------------------------------------------------------------ 870!> Checks if given element is a p element. 871!------------------------------------------------------------------------------ 872 FUNCTION isPElement( Element ) RESULT(retVal) 873!------------------------------------------------------------------------------ 874! 875! ARGUMENTS: 876! Type(Element_t) :: Element 877! INPUT: Element to check 878! 879! FUNCTION VALUE: 880! LOGICAL :: retVal 881! .TRUE. if given element is a p triangle, .FALSE. otherwise 882! 883!------------------------------------------------------------------------------ 884 IMPLICIT NONE 885 886 TYPE(Element_t), INTENT(IN) :: Element 887 LOGICAL :: retVal 888 889 retVal = ASSOCIATED(Element % PDefs) 890!------------------------------------------------------------------------------ 891 END FUNCTION isPElement 892!------------------------------------------------------------------------------ 893 894 895!------------------------------------------------------------------------------ 896!> Function checks if given element is p element triangle 897!------------------------------------------------------------------------------ 898 FUNCTION isPTriangle( Element ) RESULT(retVal) 899!------------------------------------------------------------------------------ 900! 901! ARGUMENTS: 902! Type(Element_t) :: Element 903! INPUT: Element to check 904! 905! FUNCTION VALUE: 906! LOGICAL :: retVal 907! .TRUE. if given element is a p triangle, .FALSE. otherwise 908! 909!------------------------------------------------------------------------------ 910 IMPLICIT NONE 911 912 TYPE(Element_t), INTENT(IN) :: Element 913 LOGICAL :: retVal 914 915 ! Check elementcode and p element flag 916 retVal = Element % TYPE % ElementCode/100==3 .AND. isPElement(Element) 917!------------------------------------------------------------------------------ 918 END FUNCTION isPTriangle 919!------------------------------------------------------------------------------ 920 921 922!------------------------------------------------------------------------------ 923!> Function checks if given element is p element quad 924!------------------------------------------------------------------------------ 925 FUNCTION isPQuad( Element ) RESULT(retVal) 926!------------------------------------------------------------------------------ 927! 928! ARGUMENTS: 929! Type(Element_t) :: Element 930! INPUT: Element to check 931! 932! FUNCTION VALUE: 933! LOGICAL :: retVal 934! .TRUE. if given element is a p quad, .FALSE. otherwise 935! 936!------------------------------------------------------------------------------ 937 IMPLICIT NONE 938 939 TYPE(Element_t), INTENT(IN) :: Element 940 LOGICAL :: retVal 941 942 ! Check elementcode and p element flag 943 retVal = Element % TYPE % ElementCode/100==4 .AND. isPElement(Element) 944!------------------------------------------------------------------------------ 945 END FUNCTION isPQuad 946!------------------------------------------------------------------------------ 947 948 949!------------------------------------------------------------------------------ 950!> Function checks if given element is p element tetra 951!------------------------------------------------------------------------------ 952 FUNCTION isPTetra( Element ) RESULT(retVal) 953!------------------------------------------------------------------------------ 954! 955! ARGUMENTS: 956! Type(Element_t) :: Element 957! INPUT: Element to check 958! 959! FUNCTION VALUE: 960! LOGICAL :: retVal 961! .TRUE. if given element is a p tetra, .FALSE. otherwise 962! 963!------------------------------------------------------------------------------ 964 IMPLICIT NONE 965 966 TYPE(Element_t), INTENT(IN) :: Element 967 LOGICAL :: retVal 968 969 retVal = Element % TYPE % ElementCode/100==5 .AND. isPElement(Element) 970!------------------------------------------------------------------------------ 971 END FUNCTION isPTetra 972!------------------------------------------------------------------------------ 973 974 975!------------------------------------------------------------------------------ 976!> Function checks if given element is p element wedge 977!------------------------------------------------------------------------------ 978 FUNCTION isPWedge( Element ) RESULT(retVal) 979!------------------------------------------------------------------------------ 980! 981! ARGUMENTS: 982! Type(Element_t) :: Element 983! INPUT: Element to check 984! 985! FUNCTION VALUE: 986! LOGICAL :: retVal 987! .TRUE. if given element is a p wedge, .FALSE. otherwise 988! 989!------------------------------------------------------------------------------ 990 IMPLICIT NONE 991 992 TYPE(Element_t), INTENT(IN) :: Element 993 LOGICAL :: retVal 994 995 retVal = Element % TYPE % ElementCode/100==7 .AND. isPElement(Element) 996!------------------------------------------------------------------------------ 997 END FUNCTION isPWedge 998!------------------------------------------------------------------------------ 999 1000 1001!------------------------------------------------------------------------------ 1002!> Function checks if given element is p element pyramid 1003!------------------------------------------------------------------------------ 1004 FUNCTION isPPyramid( Element ) RESULT(retVal) 1005!------------------------------------------------------------------------------ 1006! 1007! ARGUMENTS: 1008! Type(Element_t) :: Element 1009! INPUT: Element to check 1010! 1011! FUNCTION VALUE: 1012! LOGICAL :: retVal 1013! .TRUE. if given element is a p pyramid, .FALSE. otherwise 1014! 1015!------------------------------------------------------------------------------ 1016 IMPLICIT NONE 1017 1018 TYPE(Element_t), INTENT(IN) :: ELement 1019 LOGICAL :: retVal 1020 1021 retVal = Element % TYPE % ElementCode/100==6 .AND. isPElement(Element) 1022!------------------------------------------------------------------------------ 1023 END FUNCTION isPPyramid 1024!------------------------------------------------------------------------------ 1025 1026!------------------------------------------------------------------------------ 1027!> Function checks if given element is p element brick 1028!------------------------------------------------------------------------------ 1029 FUNCTION isPBrick( Element ) RESULT(retVal) 1030!------------------------------------------------------------------------------ 1031! 1032! ARGUMENTS: 1033! Type(Element_t) :: Element 1034! INPUT: Element to check 1035! 1036! FUNCTION VALUE: 1037! LOGICAL :: retVal 1038! .TRUE. if given element is a p brick, .FALSE. otherwise 1039! 1040!------------------------------------------------------------------------------ 1041 IMPLICIT NONE 1042 1043 TYPE(Element_t), INTENT(IN) :: ELement 1044 LOGICAL :: retVal 1045 1046 retVal = Element % TYPE % ElementCode/100==8 .AND. isPElement(Element) 1047!------------------------------------------------------------------------------ 1048 END FUNCTION isPBrick 1049!------------------------------------------------------------------------------ 1050 1051!------------------------------------------------------------------------------ 1052!> Get the number of Gauss points for P-elements. 1053!------------------------------------------------------------------------------ 1054 FUNCTION getNumberOfGaussPoints( Element, Mesh ) RESULT(ngp) 1055!------------------------------------------------------------------------------ 1056 IMPLICIT NONE 1057 TYPE(Mesh_t) :: Mesh 1058 TYPE(Element_t) :: Element 1059 INTEGER :: ngp 1060!------------------------------------------------------------------------------ 1061 INTEGER :: edgeP, faceP, bubbleP, TrueBubbleP, nb, maxp 1062 1063 IF (.NOT. ASSOCIATED(Element % PDefs)) THEN 1064 CALL Warn('PElementBase::getNumberOfGaussPoints','Element not p element') 1065 ngp = 0 1066 RETURN 1067 END IF 1068 1069 ! Max p of edges 1070 edgeP = 0 1071 IF ( Element % TYPE % DIMENSION == 2 .OR. & 1072 Element % TYPE % DIMENSION == 3) THEN 1073 edgeP = getEdgeP( Element, Mesh ) 1074 END IF 1075 1076 ! Max p of faces 1077 faceP = 0 1078 IF ( Element % TYPE % DIMENSION == 3 ) THEN 1079 faceP = getFaceP(Element, Mesh ) 1080 END IF 1081 1082 ! Element bubble p 1083 bubbleP = 0 1084 TrueBubbleP = 0 1085 IF (Element % BDOFs > 0) THEN 1086 bubbleP = Element % PDefs % P 1087 1088 SELECT CASE( Element % TYPE % ElementCode / 100 ) 1089 CASE(3) 1090 nb = MAX( GetBubbleDOFs( Element, bubbleP ), Element % BDOFs ) 1091 bubbleP = CEILING( ( 3.0d0+SQRT(1.0d0+8.0d0*nb) ) / 2.0d0 - AEPS) 1092 1093 CASE(4) 1094 nb = MAX( GetBubbleDOFs( Element, bubbleP ), Element % BDOFs ) 1095 TrueBubbleP = CEILING( ( 5.0d0+SQRT(1.0d0+8.0d0*nb) ) / 2.0d0 - AEPS ) 1096 bubbleP = TrueBubbleP - 2 1097 1098 CASE(5) 1099 nb = MAX( GetBubbleDOFs(Element, bubbleP ), Element % BDOFs ) 1100 bubbleP = CEILING(1/3d0*(81*nb+3*SQRT(-3d0+729*nb**2))**(1/3d0)+1d0 / & 1101 (81*nb+3*SQRT(-3d0+729*nb**2))**(1/3d0)+2 - AEPS) 1102 1103 CASE(6) 1104 nb = MAX( GetBubbleDOFs(Element, bubbleP ), Element % BDOFs ) 1105 bubbleP = CEILING(1/3d0*(81*nb+3*SQRT(-3d0+729*nb**2))**(1/3d0)+1d0 / & 1106 (81*nb+3*SQRT(-3d0+729*nb**2))**(1/3d0)+2 - AEPS) - 1 1107 1108 CASE(7) 1109 nb = MAX( GetBubbleDOFs( Element, bubbleP ), Element % BDOFs ) 1110 bubbleP = CEILING(1/3d0*(81*nb+3*SQRT(-3d0+729*nb**2))**(1/3d0)+1d0 / & 1111 (81*nb+3*SQRT(-3d0+729*nb**2))**(1/3d0)+3 - AEPS) - 2 1112 1113 CASE(8) 1114 nb = MAX( GetBubbleDOFs(Element, bubbleP ), Element % BDOFs ) 1115 bubbleP = CEILING(1/3d0*(81*nb+3*SQRT(-3d0+729*nb**2))**(1/3d0)+1d0 / & 1116 (81*nb+3*SQRT(-3d0+729*nb**2))**(1/3d0)+4 - AEPS) - 4 1117 END SELECT 1118 END IF 1119 1120 ! Special quadrature may be available: 1121 IF (Element % TYPE % ElementCode / 100 == 4) THEN 1122 ! The true polynomial degree is as follows 1123 maxp = MAX(1, edgeP, faceP, TrueBubbleP) 1124 ! but this would replace the true bubble degree by a tampered value: 1125 !maxp = MAX(1, edgeP, faceP, BubbleP) 1126 1127 ! Economic quadratures cannot be used if an explicit bubble augmentation is used 1128 ! with lower-order finite elements: 1129 IF ( .NOT.(Element % PDefs % P < 4 .AND. Element % BDOFs>0) ) THEN 1130 IF (maxp > 1 .AND. maxp <= 8) THEN 1131 !PRINT *, 'SETTING SPECIAL NGP' 1132 !PRINT *, 'MAXP=',MAXP 1133 SELECT CASE(maxp) 1134 CASE(2) 1135 ngp = 8 1136 CASE(3) 1137 ngp = 12 1138 CASE(4) 1139 ngp = 20 1140 CASE(5) 1141 ngp = 25 1142 CASE(6) 1143 ngp = 36 1144 CASE(7) 1145 ngp = 45 1146 CASE(8) 1147 ngp = 60 1148 END SELECT 1149 RETURN 1150 END IF 1151 END IF 1152 END IF 1153 ! Get the number r of Gauss points for the product of two basis functions: 1154 ! r = (2*max(p)+1)/2 1155 maxp = MAX(1, edgeP, faceP, bubbleP) + 1 1156 ! The number of Gauss points based on the Cartesian product (more efficient 1157 ! rules could be defined): 1158 ngp = maxp ** Element % TYPE % DIMENSION 1159!------------------------------------------------------------------------------ 1160 END FUNCTION getNumberOfGaussPoints 1161!------------------------------------------------------------------------------ 1162 1163 1164!------------------------------------------------------------------------------ 1165 FUNCTION getEdgeP( Element, Mesh ) RESULT(edgeP) 1166!------------------------------------------------------------------------------ 1167 IMPLICIT NONE 1168 1169 TYPE(Mesh_t) :: Mesh 1170 TYPE(Element_t) :: Element 1171 TYPE(Element_t), POINTER :: Edge 1172 1173 INTEGER :: edgeP, i 1174 1175 IF (.NOT. ASSOCIATED(Element % PDefs)) THEN 1176 CALL Warn('PElementBase::getEdgeP','Element not p element') 1177 edgeP = 0 1178 RETURN 1179 END IF 1180 1181 ! Get max p of edges of element if any 1182 edgeP = 0 1183 IF (ASSOCIATED(Element % EdgeIndexes)) THEN 1184 DO i=1, Element % TYPE % NumberOfEdges 1185 Edge => Mesh % Edges(Element % EdgeIndexes(i)) 1186 ! Here if edge has no dofs it effectively has degree of 0 1187 IF (Edge % BDOFs <= 0) CYCLE 1188 edgeP = MAX(edgeP,Edge % PDefs % P) 1189 END DO 1190 END IF 1191!------------------------------------------------------------------------------ 1192 END FUNCTION getEdgeP 1193!------------------------------------------------------------------------------ 1194 1195!------------------------------------------------------------------------------ 1196 FUNCTION getFaceP( Element, Mesh ) RESULT(faceP) 1197!------------------------------------------------------------------------------ 1198 IMPLICIT NONE 1199 1200 TYPE(Element_t) :: Element 1201 TYPE(Element_t), POINTER :: Face 1202 INTEGER :: faceP, i 1203 TYPE(Mesh_t) :: Mesh 1204 1205 IF (.NOT. ASSOCIATED(Element % PDefs)) THEN 1206 CALL Warn('PElementBase::getFaceP','Element not p element') 1207 faceP = 0 1208 RETURN 1209 END IF 1210 1211 faceP = 0 1212 IF (ASSOCIATED(Element % FaceIndexes)) THEN 1213 DO i=1,Element % TYPE % NumberOfFaces 1214 Face => Mesh % Faces( Element % FaceIndexes(i) ) 1215 ! Here if face has no dofs it effectively has degree of 0 1216 IF (Face % BDOFs <= 0) CYCLE 1217 faceP = MAX(faceP, Face % PDefs % P) 1218 END DO 1219 END IF 1220!------------------------------------------------------------------------------ 1221 END FUNCTION getFaceP 1222!------------------------------------------------------------------------------ 1223 1224!------------------------------------------------------------------------------ 1225!> Get the number of Gauss points for the faces of p-elements 1226!------------------------------------------------------------------------------ 1227 FUNCTION getNumberOfGaussPointsFace( Face, Mesh ) RESULT(ngp) 1228!------------------------------------------------------------------------------ 1229 IMPLICIT NONE 1230 1231 TYPE(Element_t), POINTER :: Face 1232 TYPE(Mesh_t) :: Mesh 1233 INTEGER :: ngp 1234!------------------------------------------------------------------------------ 1235 TYPE(Element_t), POINTER :: Edge 1236 INTEGER :: edgeP, faceP, i, maxp 1237!------------------------------------------------------------------------------ 1238 ! Get the max p of edges contained in the face 1239 edgeP = 0 1240 DO i=1, Face % TYPE % NumberOfEdges 1241 Edge => Mesh % Edges ( Face % EdgeIndexes(i) ) 1242 edgeP = MAX(edgeP, Edge % PDefs % P) 1243 END DO 1244 1245 IF (Face % BDOFs <= 0) THEN 1246 ! If no face dofs, the max p is defined by the edges 1247 maxp = edgeP 1248 ELSE 1249 maxp = MAX(edgeP, Face % PDefs % P) 1250 END IF 1251 1252 ! An economic quadrature may be available: 1253 IF (Face % TYPE % ElementCode / 100 == 4) THEN 1254 !IF ( .NOT.(maxp < 4 .AND. Face % BDOFs>0) ) THEN 1255 IF (maxp > 1 .AND. maxp <= 8) THEN 1256 SELECT CASE(maxp) 1257 CASE(2) 1258 ngp = 8 1259 CASE(3) 1260 ngp = 12 1261 CASE(4) 1262 ngp = 20 1263 CASE(5) 1264 ngp = 25 1265 CASE(6) 1266 ngp = 36 1267 CASE(7) 1268 ngp = 45 1269 CASE(8) 1270 ngp = 60 1271 END SELECT 1272 RETURN 1273 END IF 1274 END IF 1275 1276 ! Get the standard number of Gauss points: 1277 i = maxp + 1 1278 ngp = i ** 2 1279!------------------------------------------------------------------------------ 1280 END FUNCTION getNumberOfGaussPointsFace 1281!------------------------------------------------------------------------------ 1282 1283 1284!------------------------------------------------------------------------------ 1285!> Subroutine for getting reference p element nodes (because these are NOT 1286!> yet defined in element description files) 1287!------------------------------------------------------------------------------ 1288 SUBROUTINE GetRefPElementNodes(Element, U, V, W) 1289!------------------------------------------------------------------------------ 1290 IMPLICIT NONE 1291 TYPE(ElementType_t) :: Element 1292 REAL(KIND=dp) :: U(:), V(:), W(:) 1293!-------------------------------------------------------------------------------- 1294 INTEGER :: n 1295!-------------------------------------------------------------------------------- 1296 ! Reserve space for element nodes 1297 n = Element % NumberOfNodes 1298 1299 IF(.NOT.ALLOCATED(element % N_NodeU).AND.ALLOCATED(Element % NodeU) ) THEN 1300 ALLOCATE(Element % N_NodeU(n), Element % N_NodeV(n), Element % N_NodeW(n)) 1301 element % N_NodeU = element % NodeU 1302 element % N_NodeV = element % NodeV 1303 element % N_NodeW = element % NOdeW 1304 END IF 1305 1306 ! Select by element type given 1307 SELECT CASE(Element % ElementCode / 100) 1308 ! Line 1309 CASE(2) 1310 U(1:n) = (/ -1d0,1d0 /) 1311 ! Triangle 1312 CASE(3) 1313 U(1:n) = (/ -1d0,1d0,0d0 /) 1314 V(1:n) = (/ 0d0,0d0,SQRT(3.0d0) /) 1315 ! Quad 1316 CASE(4) 1317 U(1:n) = (/ -1d0,1d0,1d0,-1d0 /) 1318 V(1:n) = (/ -1d0,-1d0,1d0,1d0 /) 1319 ! Tetrahedron 1320 CASE(5) 1321 U(1:n) = (/ -1d0,1d0,0d0,0d0 /) 1322 V(1:n) = (/ 0d0,0d0,SQRT(3.0d0),1.0d0/SQRT(3.0d0) /) 1323 W(1:n) = (/ 0d0,0d0,0d0,2*SQRT(2.0d0/3.0d0) /) 1324 ! Pyramid 1325 CASE(6) 1326 U(1:n) = (/ -1d0,1d0,1d0,-1d0,0d0 /) 1327 V(1:n) = (/ -1d0,-1d0,1d0,1d0,0d0 /) 1328 W(1:n) = (/ 0d0,0d0,0d0,0d0,SQRT(2.0d0) /) 1329 ! Wedge 1330 CASE(7) 1331 U(1:n) = (/ -1d0,1d0,0d0,-1d0,1d0,0d0 /) 1332 V(1:n) = (/ 0d0,0d0,SQRT(3.0d0),0d0,0d0,SQRT(3.0d0) /) 1333 W(1:n) = (/ -1d0,-1d0,-1d0,1d0,1d0,1d0 /) 1334 ! Brick 1335 CASE(8) 1336 U(1:n) = (/ -1d0,1d0,1d0,-1d0,-1d0,1d0,1d0,-1d0 /) 1337 V(1:n) = (/ -1d0,-1d0,1d0,1d0,-1d0,-1d0,1d0,1d0 /) 1338 W(1:n) = (/ -1d0,-1d0,-1d0,-1d0,1d0,1d0,1d0,1d0 /) 1339 CASE DEFAULT 1340 WRITE(Message,'(A,I0)') 'Unknown element type: ',Element % ElementCode 1341 CALL Warn('PElementMaps::GetRefPElementNodes',Message) 1342 END SELECT 1343!------------------------------------------------------------------------------ 1344 END SUBROUTINE GetRefPElementNodes 1345!------------------------------------------------------------------------------ 1346 1347 1348END MODULE 1349 1350!> \} 1351