Initial commit from Polytechnique Montreal

1 files changed, 95 insertions, 0 deletions

diff --git a/Trivac/src/FLDBH2.f b/Trivac/src/FLDBH2.f
new file mode 100755
index 0000000..7dd7db4
--- /dev/null
+++ b/Trivac/src/FLDBH2.f
@@ -0,0 +1,95 @@
+*DECK FLDBH2
+      SUBROUTINE FLDBH2 (ISPLH,NEL,NUN,NELEM,EVECT,VOL,IDL,KN,QFR,RH,RT)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Calculation of the averaged flux with a linear Lagrangian finite
+* element or mesh corner finite difference method in hexagonal geometry.
+*
+*Copyright:
+* Copyright (C) 2002 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version
+*
+*Author(s): A. Hebert
+*
+*Parameters: input
+* ISPLH   type of hexagonal mesh-splitting: =1 for complete hexagons;
+*         >1 for triangular mesh-splitting.
+* NEL     number of hexagons.
+* NUN     number of unknowns per energy group.
+* NELEM   number of finite elements (hexagons or triangles) excluding
+*         the virtual elements.
+* EVECT   variational coefficients of the flux. The information is
+*         contained in position EVECT(1) to EVECT(LL4) where LL4 is the
+*         order of the system matrices.
+* VOL     volume of each hexagon.
+* IDL     position of the average flux component associated with each
+*         hexagon.
+* KN      element-ordered unknown list. The dimension of KN is equal
+*         to (LC+1)*NELEM where LC=6 (hexagons) or 3 (triangles).
+* QFR     element-ordered albedo information. The dimension of QFR is
+*         equal to (LC+1)*NELEM.
+* RH      unit matrix
+* RT      unit matrix
+*
+*Parameters: output
+* EVECT   averaged fluxes. The information is contained in positions
+*         EVECT(IDL(I)).
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER ISPLH,NEL,NUN,NELEM,IDL(NEL),KN(*)
+      REAL EVECT(NUN),VOL(NEL),QFR(*),RH(6,6),RT(3,3)
+*----
+*  LOCAL VARIABLES
+*----
+      REAL T(6)
+*----
+*  COMPUTE THE LINEAR PRODUCT VECTOR T
+*----
+      IF(ISPLH.EQ.1) THEN
+*        HEXAGONAL BASIS.
+         LC=6
+         DO 15 I=1,6
+         T(I)=0.0
+         DO 10 J=1,6
+         T(I)=T(I)+RH(I,J)
+   10    CONTINUE
+   15    CONTINUE
+         CONST=1.5*SQRT(3.0)
+      ELSE
+*        TRIANGULAR BASIS.
+         LC=3
+         DO 25 I=1,3
+         T(I)=0.0
+         DO 20 J=1,3
+         T(I)=T(I)+RT(I,J)
+   20    CONTINUE
+   25    CONTINUE
+         CONST=0.25*SQRT(3.0)
+      ENDIF
+*
+      DO 30 KHEX=1,NEL
+      IF(IDL(KHEX).NE.0) EVECT(IDL(KHEX))=0.0
+   30 CONTINUE
+      NUM1=0
+      DO 60 K=1,NELEM
+      KHEX=KN(NUM1+LC+1)
+      IF(VOL(KHEX).EQ.0.0) GO TO 50
+      DO 40 I=1,LC
+      IND1=KN(NUM1+I)
+      IF(IND1.EQ.0) GO TO 40
+      SS=T(I)*QFR(NUM1+LC+1)/(CONST*VOL(KHEX))
+      EVECT(IDL(KHEX))=EVECT(IDL(KHEX))+SS*EVECT(IND1)
+   40 CONTINUE
+   50 NUM1=NUM1+LC+1
+   60 CONTINUE
+      RETURN
+      END