add reference comments and cleanup code in diffmtc

OpenFAST · Sep 28, 2017 · abb0d0d · abb0d0d
1 parent d8dddc2
commit abb0d0d
Showing 1 changed file with 76 additions and 56 deletions.
diff --git a/modules-local/beamdyn/src/BeamDyn.f90 b/modules-local/beamdyn/src/BeamDyn.f90
@@ -2660,81 +2660,101 @@ SUBROUTINE BD_GyroForce( nelem, idx_qp, m )
 END SUBROUTINE BD_GyroForce
 
 
-
-
 !-----------------------------------------------------------------------------------------------------------------------------------
 !> calculate Lagrangian interpolant tensor at ns points where basis
 !! functions are assumed to be associated with (np+1) GLL points on [-1,1]
 SUBROUTINE BD_diffmtc( p,GLL_nodes,Shp,ShpDer )
 
+   ! See Bauchau equations 17.1 - 17.5
+
    TYPE(BD_ParameterType), INTENT(IN   )  :: p              !< Parameters
    REAL(BDKi),             INTENT(IN   )  :: GLL_nodes(:)   !< GLL_nodes(p%nodes_per_elem): location of the (p%nodes_per_elem) p%GLL points
    REAL(BDKi),             INTENT(INOUT)  :: Shp(:,:)       !< p%Shp    (or another Shp array for when we add outputs at arbitrary locations)
    REAL(BDKi),             INTENT(INOUT)  :: ShpDer(:,:)    !< p%ShpDer (or another Shp array for when we add outputs at arbitrary locations)
 
-   REAL(BDKi)                  :: dnum
-   REAL(BDKi)                  :: den
-   REAL(BDKi),        PARAMETER:: eps = SQRT(EPSILON(eps)) !1.0D-08
-   INTEGER(IntKi)              :: l
-   INTEGER(IntKi)              :: j
-   INTEGER(IntKi)              :: i
-   INTEGER(IntKi)              :: k
+   REAL(BDKi)                  :: dnum, den
+   REAL(BDKi),       PARAMETER :: eps = SQRT(EPSILON(eps)) !1.0D-08
+   INTEGER(IntKi)              :: i, j, k, l
 
+   ! initialize shape functions to 0
    Shp(:,:)     = 0.0_BDKi
    ShpDer(:,:)  = 0.0_BDKi
-
-
+   
+   ! shape function derivative
    do j = 1,p%nqp
       do l = 1,p%nodes_per_elem
-
-       if ((abs(p%QPtN(j)-1.).LE.eps).AND.(l.EQ.p%nodes_per_elem)) then           !adp: FIXME: do we want to compare to eps, or EqualRealNos???
-         ShpDer(l,j) = REAL((p%nodes_per_elem)*(p%nodes_per_elem-1), BDKi)/4.0_BDKi
-       elseif ((abs(p%QPtN(j)+1.).LE.eps).AND.(l.EQ.1)) then
-         ShpDer(l,j) = -REAL((p%nodes_per_elem)*(p%nodes_per_elem-1), BDKi)/4.0_BDKi
-       elseif (abs(p%QPtN(j)-GLL_nodes(l)).LE.eps) then
-         ShpDer(l,j) = 0.0_BDKi
-       else
-         ShpDer(l,j) = 0.0_BDKi
-         den = 1.0_BDKi
-         do i = 1,p%nodes_per_elem
-           if (i.NE.l) then
-             den = den*(GLL_nodes(l)-GLL_nodes(i))
-           endif
-           dnum = 1.0_BDKi
-           do k = 1,p%nodes_per_elem
-             if ((k.NE.l).AND.(k.NE.i).AND.(i.NE.l)) then
-               dnum = dnum*(p%QPtN(j)-GLL_nodes(k))
-             elseif (i.EQ.l) then
-               dnum = 0.0_BDKi
-             endif
-           enddo
-           ShpDer(l,j) = ShpDer(l,j) + dnum
-         enddo
-         ShpDer(l,j) = ShpDer(l,j)/den
-       endif
-     enddo
+
+         ! if this quadrature point is located at 1/-1 AND there is exactly 1 node per element
+         if ( (abs(p%QPtN(j)-1.).le.eps).and.(p%nodes_per_elem.eq.1) ) then           !adp: FIXME: do we want to compare to eps, or EqualRealNos???
+            ShpDer(l,j) = REAL( p%nodes_per_elem * (p%nodes_per_elem-1), BDKi )/4.0_BDKi
+
+         ! if this quadrature point is located at 1/-1 AND this is the first node of a few
+         elseif ( (abs(p%QPtN(j)+1.).le.eps).and.(l.eq.1) ) then
+            ShpDer(l,j) = -REAL((p%nodes_per_elem)*(p%nodes_per_elem-1), BDKi)/4.0_BDKi
+
+         ! if this quadrature point is located at the GLL node
+         elseif ( abs(p%QPtN(j)-GLL_nodes(l)).le.eps ) then
+            ShpDer(l,j) = 0.0_BDKi
+
+         ! otherwise, ...
+         else  
+            ! initialize denominator to 1
+            den = 1.0_BDKi
+
+            ! loop over all nodes in this element
+            do i = 1,p%nodes_per_elem
+
+               if (i.NE.l) then
+                  den = den*(GLL_nodes(l)-GLL_nodes(i))
+               end if
+
+               ! initialize the numerator to 1
+               dnum = 1.0_BDKi
+
+               ! loop over all nodes in this element
+               do k = 1,p%nodes_per_elem
+
+                  if ( (k.NE.l).AND.(k.NE.i).AND.(i.NE.l) ) then
+                     dnum = dnum*(p%QPtN(j)-GLL_nodes(k))
+                  elseif (i.EQ.l) then
+                     dnum = 0.0_BDKi
+                  endif
+
+               enddo
+
+               ShpDer(l,j) = ShpDer(l,j) + dnum
+            enddo
+            ShpDer(l,j) = ShpDer(l,j)/den
+         endif
+      enddo
    enddo
-
+
+   ! shape function
    do j = 1,p%nqp
       do l = 1,p%nodes_per_elem
-
-       if(abs(p%QPtN(j)-GLL_nodes(l)).LE.eps) then
-         Shp(l,j) = 1.0_BDKi
-       else
-         dnum = 1.0_BDKi
-         den  = 1.0_BDKi
-         do k = 1,p%nodes_per_elem
-           if (k.NE.l) then
-             den  = den *(GLL_nodes(l) - GLL_nodes(k))
-             dnum = dnum*(p%QPtN(j) - GLL_nodes(k))
-           endif
-         enddo
-         Shp(l,j) = dnum/den
-       endif
-     enddo
+
+         ! if this quadrature point is located at the GLL node
+         if ( abs(p%QPtN(j)-GLL_nodes(l)).LE.eps ) then
+            Shp(l,j) = 1.0_BDKi
+
+         else
+            dnum = 1.0_BDKi
+            den  = 1.0_BDKi
+
+            do k = 1,p%nodes_per_elem
+
+               if (k.NE.l) then
+                  den  = den *(GLL_nodes(l) - GLL_nodes(k))
+                  dnum = dnum*(p%QPtN(j) - GLL_nodes(k))
+               endif
+
+            enddo
+
+            Shp(l,j) = dnum/den
+
+         endif
+      enddo
    enddo
-
-
  END SUBROUTINE BD_diffmtc