<commit_before>c -*- Fortran -*-
c
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c
c                             <NAME>
c                       Rensselaer Polytechnic Institute
c                        (C) 2004  All Rights Reserved
c
c  Copyright 2004 Rensselaer Polytechnic Institute.
c  All worldwide rights reserved.  A license to use, copy, modify and
c  distribute this software for non-commercial research purposes only
c  is hereby granted, provided that this copyright notice and
c  accompanying disclaimer is not modified or removed from the software.
c
c  DISCLAIMER:  The software is distributed "AS IS" without any express
c  or implied warranty, including but not limited to, any implied
c  warranties of merchantability or fitness for a particular purpose
c  or any warranty of non-infringement of any current or pending patent
c  rights.  The authors of the software make no representations about
c  the suitability of this software for any particular purpose.  The
c  entire risk as to the quality and performance of the software is with
c  the user.  Should the software prove defective, the user assumes the
c  cost of all necessary servicing, repair or correction.  In
c  particular, neither Rensselaer Polytechnic Institute, nor the authors
c  of the software are liable for any indirect, special, consequential,
c  or incidental damages related to the software, to the maximum extent
c  the law permits.
c
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c
c
      subroutine pquadhex(sh,shj,gauss,infetype,intord)
c
c... Subroutine to compute shape functions in natural coordinates,
c    integration points, and weights for a quadratic (20-node)
c    hexahedron.
c
      include "implicit.inc"
c
c...  parameter definitions
c
      include "ndimens.inc"
      include "nshape.inc"
      include "nconsts.inc"
      include "rconsts.inc"
c
c...  subroutine arguments
c
      integer intord
      integer infetype(4)
      double precision sh(nsd+1,nenmax,ngaussmax)
      double precision shj(nsd+1,nenmax,ngaussmax)
      double precision gauss(nsd+1,ngaussmax)
c
c...  local constants
c
      double precision r(20),s(20),t(20)
      data r/-1d0, 1d0, 1d0,-1d0,-1d0, 1d0, 1d0,-1d0,
     &        0d0, 1d0, 0d0,-1d0, 0d0, 1d0, 0d0,-1d0,
     &       -1d0, 1d0, 1d0,-1d0/
      data s/-1d0,-1d0, 1d0, 1d0,-1d0,-1d0, 1d0, 1d0,
     &       -1d0, 0d0, 1d0, 0d0,-1d0, 0d0, 1d0, 0d0,
     &       -1d0,-1d0, 1d0, 1d0/
      data t/ 1d0, 1d0, 1d0, 1d0,-1d0,-1d0,-1d0,-1d0,
     &        1d0, 1d0, 1d0, 1d0,-1d0,-1d0,-1d0,-1d0,
     &        0d0, 0d0, 0d0, 0d0/
c
c...  intrinsic functions
c
      intrinsic abs,sqrt,dble
c
c...  user-defined functions
c
      double precision gquad,dgquad,hquad,dhquadr,dhquads,dhquadt
c
c...  local variables
c
      integer nen,ngauss,nec,nee,i,l,l1,l2,l3,nshsize,ngssize
      double precision g1,w1,w2,rr,ss,tt,rrw,ssw,ttw,drr,dss,dtt
      double precision uu,dur,dus,dut,v,vi,rp,ri,sp,si,tp,ti
c
c...  function definitions
c
      gquad(v,vi)=(half*(one-abs(vi))+half)*(one+v*vi+(abs(vi)-one)*v*v)
      dgquad(v,vi)=(half*(one-abs(vi))+half)*(vi+two*v*(abs(vi)-one))
      hquad(rp,sp,tp,ri,si,ti)=abs(ri*si*ti)*(ri*rp+si*sp+ti*tp-two)+
     & one-abs(ri*si*ti)
      dhquadr(ri,si,ti)=abs(ri*si*ti)*ri
      dhquads(ri,si,ti)=abs(ri*si*ti)*si
      dhquadt(ri,si,ti)=abs(ri*si*ti)*ti
c
c...  definitions
c
      nshsize=(nsd+1)*nenmax*ngaussmax
      ngssize=(nsd+1)*ngaussmax
c
c...  Quadratic hex definition
c
      nen=20
      ngauss=ieight
      nec=nsd*nen
      nee=ndof*nen
      do l=1,ngauss
        gauss(1,l)=r(l)*root3i
        gauss(2,l)=s(l)*root3i
        gauss(3,l)=t(l)*root3i
        gauss(4,l)=one
      end do
      if(intord.ne.2) then
        ngauss=27
        g1=sqrt(three/five)
        w1=five/nine
        w2=eight/nine
        l=0
        do l3=1,3
          tt=dble(l3-2)
          ttw=tt*w1+(abs(tt)-one)*w2
          do l2=1,3
            ss=dble(l2-2)
            ssw=ss*w1+(abs(ss)-one)*w2
            do l1=1,3
              rr=dble(l1-2)
              rrw=rr*w1+(abs(rr)-one)*w2
              l=l+1
              gauss(1,l)=rr*g1
              gauss(2,l)=ss*g1
              gauss(3,l)=tt*g1
              gauss(4,l)=rrw*ssw*ttw
            end do
          end do
        end do
      end if
c
      infetype(1)=ngauss
      infetype(2)=nen
      infetype(3)=nec
      infetype(4)=nee
c
      do l=1,ngauss
        do i=1,nen
          rr=gquad(gauss(1,l),r(i))
          ss=gquad(gauss(2,l),s(i))
          tt=gquad(gauss(3,l),t(i))
          uu=hquad(gauss(1,l),gauss(2,l),gauss(3,l),r(i),s(i),t(i))
          drr=dgquad(gauss(1,l),r(i))
          dss=dgquad(gauss(2,l),s(i))
          dtt=dgquad(gauss(3,l),t(i))
          dur=dhquadr(r(i),s(i),t(i))
          dus=dhquads(r(i),s(i),t(i))
          dut=dhquadt(r(i),s(i),t(i))
          sh(4,i,l)=rr*ss*tt*uu
          sh(1,i,l)=drr*ss*tt*uu+dur*rr*ss*tt
          sh(2,i,l)=dss*rr*tt*uu+dus*rr*ss*tt
          sh(3,i,l)=dtt*rr*ss*uu+dut*rr*ss*tt
        end do
      end do
      call dcopy(nshsize,sh,ione,shj,ione)
c
      return
      end
c
c version
c $Id: pquadhex.f,v 1.3 2004/07/07 19:40:22 willic3 Exp $
c
c Generated automatically by Fortran77Mill on Wed May 21 14:15:03 2003<commit_msg>Completely revised computation of shape functions, which was being done incorrectly before. Routine now contains commented out debugging statements.<commit_after>c -*- Fortran -*-
c
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c
c                             <NAME>
c                       Rensselaer Polytechnic Institute
c                        (C) 2004  All Rights Reserved
c
c  Copyright 2004 Rensselaer Polytechnic Institute.
c  All worldwide rights reserved.  A license to use, copy, modify and
c  distribute this software for non-commercial research purposes only
c  is hereby granted, provided that this copyright notice and
c  accompanying disclaimer is not modified or removed from the software.
c
c  DISCLAIMER:  The software is distributed "AS IS" without any express
c  or implied warranty, including but not limited to, any implied
c  warranties of merchantability or fitness for a particular purpose
c  or any warranty of non-infringement of any current or pending patent
c  rights.  The authors of the software make no representations about
c  the suitability of this software for any particular purpose.  The
c  entire risk as to the quality and performance of the software is with
c  the user.  Should the software prove defective, the user assumes the
c  cost of all necessary servicing, repair or correction.  In
c  particular, neither Rensselaer Polytechnic Institute, nor the authors
c  of the software are liable for any indirect, special, consequential,
c  or incidental damages related to the software, to the maximum extent
c  the law permits.
c
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c
c
      subroutine pquadhex(sh,shj,gauss,infetype,intord)
c
c... Subroutine to compute shape functions in natural coordinates,
c    integration points, and weights for a quadratic (20-node)
c    hexahedron.
c
      include "implicit.inc"
c
c...  parameter definitions
c
      include "ndimens.inc"
      include "nshape.inc"
      include "nconsts.inc"
      include "rconsts.inc"
c
c...  subroutine arguments
c
      integer intord
      integer infetype(4)
      double precision sh(nsd+1,nenmax,ngaussmax)
      double precision shj(nsd+1,nenmax,ngaussmax)
      double precision gauss(nsd+1,ngaussmax)
c
c...  local constants
c
      double precision r(20),s(20),t(20)
      data r/-1d0, 1d0, 1d0,-1d0,-1d0, 1d0, 1d0,-1d0,
     &        0d0, 1d0, 0d0,-1d0, 0d0, 1d0, 0d0,-1d0,
     &       -1d0, 1d0, 1d0,-1d0/
      data s/-1d0,-1d0, 1d0, 1d0,-1d0,-1d0, 1d0, 1d0,
     &       -1d0, 0d0, 1d0, 0d0,-1d0, 0d0, 1d0, 0d0,
     &       -1d0,-1d0, 1d0, 1d0/
      data t/ 1d0, 1d0, 1d0, 1d0,-1d0,-1d0,-1d0,-1d0,
     &        1d0, 1d0, 1d0, 1d0,-1d0,-1d0,-1d0,-1d0,
     &        0d0, 0d0, 0d0, 0d0/
c
c...  intrinsic functions
c
      intrinsic abs,sqrt,dble
c
c...  user-defined functions
c
      double precision gquad,dgquad,hquad,dhquadr,dhquads,dhquadt
c
c...  local variables
c
      integer nen,ngauss,nec,nee,i,l,l1,l2,l3,nshsize,ngssize
      double precision g1,rr,ss,tt,rrw,ssw,ttw,drr,dss,dtt
      double precision uu,dur,dus,dut,v,vi,rp,ri,sp,si,tp,ti
      double precision w(3)
cdebug      integer idb,jdb,kdb
c
c...  function definitions
c
      gquad(v,vi)=(half*(one-abs(vi))+half)*(one+v*vi+(abs(vi)-one)*v*v)
      dgquad(v,vi)=(half*(one-abs(vi))+half)*(vi+two*v*(abs(vi)-one))
      hquad(rp,sp,tp,ri,si,ti)=abs(ri*si*ti)*(ri*rp+si*sp+ti*tp-two)+
     & one-abs(ri*si*ti)
      dhquadr(ri,si,ti)=abs(ri*si*ti)*ri
      dhquads(ri,si,ti)=abs(ri*si*ti)*si
      dhquadt(ri,si,ti)=abs(ri*si*ti)*ti
c
cdebug      write(6,*) "Hello from pquadhex_f!"
c
c
c...  definitions
c
      nshsize=(nsd+1)*nenmax*ngaussmax
      ngssize=(nsd+1)*ngaussmax
c
c...  Quadratic hex definition
c
      nen=20
      ngauss=ieight
      nec=nsd*nen
      nee=ndof*nen
      do l=1,ngauss
        gauss(1,l)=r(l)*root3i
        gauss(2,l)=s(l)*root3i
        gauss(3,l)=t(l)*root3i
        gauss(4,l)=one
      end do
      if(intord.ne.2) then
        ngauss=27
        g1=sqrt(three/five)
        w(1)=five/nine
        w(2)=eight/nine
        w(3)=five/nine
        l=0
        do l3=1,3
          tt=dble(l3-itwo)
          ttw=w(l3)
          do l2=1,3
            ss=dble(l2-itwo)
            ssw=w(l2)
            do l1=1,3
              rr=dble(l1-itwo)
              rrw=w(l1)
              l=l+1
              gauss(1,l)=rr*g1
              gauss(2,l)=ss*g1
              gauss(3,l)=tt*g1
              gauss(4,l)=rrw*ssw*ttw
            end do
          end do
        end do
      end if
c*      if(intord.ne.2) then
c*        ngauss=27
c*        g1=sqrt(three/five)
c*        w1=five/nine
c*        w2=eight/nine
c*        l=0
c*        do l3=1,3
c*          tt=dble(l3-2)
c*          ttw=tt*w1+(abs(tt)-one)*w2
c*          do l2=1,3
c*            ss=dble(l2-2)
c*            ssw=ss*w1+(abs(ss)-one)*w2
c*            do l1=1,3
c*              rr=dble(l1-2)
c*              rrw=rr*w1+(abs(rr)-one)*w2
c*              l=l+1
c*              gauss(1,l)=rr*g1
c*              gauss(2,l)=ss*g1
c*              gauss(3,l)=tt*g1
c*              gauss(4,l)=rrw*ssw*ttw
c*            end do
c*          end do
c*        end do
c*      end if
cdebug      do idb=1,ngauss
cdebug        write(6,*) "gauss:",(gauss(jdb,idb),jdb=1,4)
cdebug      end do
c
      infetype(1)=ngauss
      infetype(2)=nen
      infetype(3)=nec
      infetype(4)=nee
c
      do l=1,ngauss
        do i=1,nen
          rr=gquad(gauss(1,l),r(i))
          ss=gquad(gauss(2,l),s(i))
          tt=gquad(gauss(3,l),t(i))
          uu=hquad(gauss(1,l),gauss(2,l),gauss(3,l),r(i),s(i),t(i))
          drr=dgquad(gauss(1,l),r(i))
          dss=dgquad(gauss(2,l),s(i))
          dtt=dgquad(gauss(3,l),t(i))
          dur=dhquadr(r(i),s(i),t(i))
          dus=dhquads(r(i),s(i),t(i))
          dut=dhquadt(r(i),s(i),t(i))
          sh(4,i,l)=rr*ss*tt*uu
          sh(1,i,l)=drr*ss*tt*uu+dur*rr*ss*tt
          sh(2,i,l)=dss*rr*tt*uu+dus*rr*ss*tt
          sh(3,i,l)=dtt*rr*ss*uu+dut*rr*ss*tt
        end do
      end do
cdebug      do idb=1,ngauss
cdebug        do jdb=1,nen
cdebug          write(6,*) "sh:",(sh(kdb,jdb,idb),kdb=1,4)
cdebug        end do
cdebug      end do
      call dcopy(nshsize,sh,ione,shj,ione)
c
      return
      end
c
c version
c $Id: pquadhex.f,v 1.4 2004/08/02 21:20:44 willic3 Exp $
c
c Generated automatically by Fortran77Mill on Wed May 21 14:15:03 2003
{'fraction_non_alphanumeric': 0.15194949236413274, 'fraction_numerical': 0.04504735090862554, 'mean_word_length': 2.4804038004750595, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2746011', 'n_tokens_mistral': 4884, 'n_tokens_neox': 4614, 'n_words': 1053}