diff --git a/fem/src/ElementDescription.F90 b/fem/src/ElementDescription.F90 index f85c50affb..406b4944cf 100644 --- a/fem/src/ElementDescription.F90 +++ b/fem/src/ElementDescription.F90 @@ -5387,7 +5387,7 @@ RECURSIVE FUNCTION FaceElementInfo( Element, Nodes, u, v, w, F, detF, & REAL(KIND=dp), OPTIONAL :: DivFBasis(:) !< The divergence of basis functions with respect to the local coordinates LOGICAL, OPTIONAL :: BDM !< If .TRUE., a basis for BDM space is constructed LOGICAL, OPTIONAL :: Dual !< If .TRUE., create an alternate dual basis - INTEGER, OPTIONAL :: BasisDegree(:) !< This a dummy parameter at the moment + INTEGER, OPTIONAL :: BasisDegree !< This has limited functionality at the moment LOGICAL, OPTIONAL :: ApplyPiolaTransform !< If .TRUE., perform the Piola transform so that, instead of b !< and Div b, return B(f(p)) and (div B)(f(p)) with B(x) the basis !< functions on the physical element and div the spatial divergence operator. @@ -5403,9 +5403,10 @@ RECURSIVE FUNCTION FaceElementInfo( Element, Nodes, u, v, w, F, detF, & INTEGER :: DOFs INTEGER :: n, dim, q, i, j, k, ni, nj, nk, I1, I2 INTEGER :: FDofMap(6,4), DofsPerFace, FaceIndices(4) + INTEGER :: Family, RTDegree REAL(KIND=dp) :: LF(3,3) REAL(KIND=dp) :: DivBasis(MaxDOFs) - REAL(KIND=dp) :: dLbasisdx(MAX(SIZE(Nodes % x),SIZE(Basis)),3), S, D1, D2, fun, dfun + REAL(KIND=dp) :: dLbasisdx(MAX(SIZE(Nodes % x),SIZE(Basis)),3), S, D1, D2, fun, dfun, wfun(2) REAL(KIND=dp) :: WorkBasis(24,3), WorkDivBasis(24) LOGICAL :: ReverseSign(6), CreateBDMBasis, Parallel @@ -5421,6 +5422,11 @@ RECURSIVE FUNCTION FaceElementInfo( Element, Nodes, u, v, w, F, detF, & !--------------------------------------------------------------------- CreateBDMBasis = .FALSE. IF ( PRESENT(BDM) ) CreateBDMBasis = BDM + RTDegree = 0 + IF (PRESENT(BasisDegree)) THEN + RTDegree = BasisDegree - 1 + IF (BasisDegree > 2) CALL Fatal('ElementDescription::FaceElementInfo', 'Unsupported element degree') + END IF CreateDualBasis = .FALSE. IF ( PRESENT(Dual) ) CreateDualBasis = Dual PerformPiolaTransform = .FALSE. @@ -5436,6 +5442,7 @@ RECURSIVE FUNCTION FaceElementInfo( Element, Nodes, u, v, w, F, detF, & dLbasisdx = 0.0d0 n = Element % TYPE % NumberOfNodes dim = Element % TYPE % DIMENSION + ! cdim = CoordinateSystemDimension() IF ( Element % TYPE % ElementCode == 101 ) THEN detF = 1.0d0 @@ -5451,7 +5458,20 @@ RECURSIVE FUNCTION FaceElementInfo( Element, Nodes, u, v, w, F, detF, & ! Remark: Using reference elements having the faces of the same area ! simplifies the implementation of element assembly procedures. !----------------------------------------------------------------------- - SELECT CASE(Element % TYPE % ElementCode / 100) + Family = Element % TYPE % ElementCode / 100 + SELECT CASE(Family) + CASE(2) + DO q=1,2 + Basis(q) = LineNodalPBasis(q, u) + dLBasisdx(q,1) = dLineNodalPBasis(q, u) + END DO + IF (RTDegree == 1) THEN + DOFs = 3 + ! Basis(3) = LineBubblePBasis(2, u) + ! dLBasisdx(q,1) = dLineBubblePBasis(2, u) + ELSE + DOFs = 2 + END IF CASE(3) DO q=1,n Basis(q) = TriangleNodalPBasis(q, u, v) @@ -5498,6 +5518,17 @@ RECURSIVE FUNCTION FaceElementInfo( Element, Nodes, u, v, w, F, detF, & !--------------------------------------------------------------------------- SELECT CASE(Element % TYPE % ElementCode / 100) + CASE(2) + ! TO DO: Implement possible sign reversions + FBasis(1,1) = -Basis(1) + DivBasis(1) = -dLBasisdx(q,1) + FBasis(2,1) = Basis(2) + DivBasis(2) = -dLBasisdx(q,2) + IF (RTDegree > 0) THEN + FBasis(3,1) = 4.0d0 * Basis(1) * Basis(2) + DivBasis(2) = 4.0d0 * dLBasisdx(1,1) * Basis(2) + 4.0d0 * Basis(1) * dLBasisdx(2,1) + END IF + CASE(3) !---------------------------------------------------------------- ! Note that the global orientation of face normal is taken to be @@ -5573,32 +5604,213 @@ RECURSIVE FUNCTION FaceElementInfo( Element, Nodes, u, v, w, F, detF, & END DO ELSE - DOFs = 3 + SELECT CASE (RTDegree) + CASE(0) + DOFs = 3 + + FBasis(1,1) = SQRT(3.0d0)/6.0d0 * u + FBasis(1,2) = -0.5d0 + SQRT(3.0d0)/6.0d0 * v + DivBasis(1) = SQRT(3.0d0)/3.0d0 + IF (ReverseSign(1)) THEN + FBasis(1,:) = -FBasis(1,:) + DivBasis(1) = -DivBasis(1) + END IF - FBasis(1,1) = SQRT(3.0d0)/6.0d0 * u - FBasis(1,2) = -0.5d0 + SQRT(3.0d0)/6.0d0 * v - DivBasis(1) = SQRT(3.0d0)/3.0d0 - IF (ReverseSign(1)) THEN - FBasis(1,:) = -FBasis(1,:) - DivBasis(1) = -DivBasis(1) - END IF + FBasis(2,1) = SQRT(3.0d0)/6.0d0 * (1.0d0 + u) + FBasis(2,2) = SQRT(3.0d0)/6.0d0 * v + DivBasis(2) = SQRT(3.0d0)/3.0d0 + IF (ReverseSign(2)) THEN + FBasis(2,:) = -FBasis(2,:) + DivBasis(2) = -DivBasis(2) + END IF - FBasis(2,1) = SQRT(3.0d0)/6.0d0 * (1.0d0 + u) - FBasis(2,2) = SQRT(3.0d0)/6.0d0 * v - DivBasis(2) = SQRT(3.0d0)/3.0d0 - IF (ReverseSign(2)) THEN - FBasis(2,:) = -FBasis(2,:) - DivBasis(2) = -DivBasis(2) - END IF + FBasis(3,1) = SQRT(3.0d0)/6.0d0 * (-1.0d0 + u) + FBasis(3,2) = SQRT(3.0d0)/6.0d0 * v + DivBasis(3) = SQRT(3.0d0)/3.0d0 + IF (ReverseSign(3)) THEN + FBasis(3,:) = -FBasis(3,:) + DivBasis(3) = -DivBasis(3) + END IF - FBasis(3,1) = SQRT(3.0d0)/6.0d0 * (-1.0d0 + u) - FBasis(3,2) = SQRT(3.0d0)/6.0d0 * v - DivBasis(3) = SQRT(3.0d0)/3.0d0 - IF (ReverseSign(3)) THEN - FBasis(3,:) = -FBasis(3,:) - DivBasis(3) = -DivBasis(3) - END IF + CASE(1) + ! + ! We use a non-hierarchic basis which is motivated by flux reconstruction. + ! The degrees of freedom associated with the element faces (edges) can be + ! integrated for a given flux q as (q.n,w) where the weights are the Lagrange + ! basis functions. + DOFs = 8 + !------------------------------------------------- + ! Two basis functions defined on the face 12. + !------------------------------------------------- + WorkBasis(3,1) = SQRT(3.0d0)/6.0d0 * u + WorkBasis(3,2) = -0.5d0 + SQRT(3.0d0)/6.0d0 * v + WorkDivBasis(3) = SQRT(3.0d0)/3.0d0 + IF (ReverseSign(1)) THEN + WorkBasis(3,:) = -WorkBasis(3,:) + WorkDivBasis(3) = -WorkDivBasis(3) + END IF + + wfun(1) = 4.0d0 * Basis(1) - 2.0d0 * Basis(2) + wfun(2) = 4.0d0 * Basis(2) - 2.0d0 * Basis(1) + WorkBasis(1,1:2) = wfun(1) * WorkBasis(3,1:2) + WorkBasis(2,1:2) = wfun(2) * WorkBasis(3,1:2) + WorkDivBasis(1) = wfun(1) * WorkDivBasis(3) + & + SUM(WorkBasis(3,1:2) * (4.0d0 * dLBasisdx(1,1:2) - 2.0d0 * dLBasisdx(2,1:2))) + WorkDivBasis(2) = wfun(2) * WorkDivBasis(3) + & + SUM(WorkBasis(3,1:2) * (4.0d0 * dLBasisdx(2,1:2) - 2.0d0 * dLBasisdx(1,1:2))) + + i = EdgeMap(1,1) + j = EdgeMap(1,2) + ni = Element % NodeIndexes(i) + IF (Parallel) ni=Mesh % ParallelInfo % GlobalDOFs(ni) + nj = Element % NodeIndexes(j) + IF (Parallel) nj=Mesh % ParallelInfo % GlobalDOFs(nj) + IF (nj FaceIndices(3) ) THEN + k = 3 + END IF + SELECT CASE(k) + CASE(1) + FBasis(7,1:2) = WorkBasis(1,1:2) + DivBasis(7) = WorkDivBasis(1) + FBasis(8,1:2) = WorkBasis(3,1:2) + DivBasis(8) = WorkDivBasis(3) + CASE(2) + FBasis(7,1:2) = WorkBasis(1,1:2) + DivBasis(7) = WorkDivBasis(1) + FBasis(8,1:2) = WorkBasis(2,1:2) + DivBasis(8) = WorkDivBasis(2) + CASE(3) + FBasis(7,1:2) = WorkBasis(2,1:2) + DivBasis(7) = WorkDivBasis(2) + FBasis(8,1:2) = WorkBasis(3,1:2) + DivBasis(8) = WorkDivBasis(3) + END SELECT + + END SELECT END IF CASE(4)