c==============================================================================
c=== LIBRERIA DE CALCULO POR ELEMENTOS FINITOS PARA PROBLEMAS ELIPTICOS 1D ====
c==============================================================================

c______________________________________Libreria Lib1D__________________________
c
c AUTORES: F. Navarrina, I. Colominas, H. Gomez, M. Casteleiro
c                                                                 _
c        Dep. de Metodos Matematicos y de Rep., ETSICCP de La Coruna, UDC.
c
c FECHA ULTIMA REVISION: 20070524
c
c Copyright: Esta libreria es propiedad de su autor.
c            No puede utilizarse con fines comerciales sin conocimiento
c            y autorizacion de su autor.
c            Puede utilizarse con fines educativos sin ninguna restriccion.
c
c------------------------------------------------------------------------------
c===============================================================================
c================================== PRIMER NIVEL ===============================
c===============================================================================

c____________________________________________________________________CAL$Ini_Par
c     Inicializacion de parametros del programa

      Subroutine CAL$Ini_Par(ntype,ndime,npoin,ndofn,nnode,nstre,ngaue,
     .           nprop,ntloa,nevab,ntotv,posge,weige,shape,derie)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)
   
c.....Parametros que dependen del tipo de calculo
      if (ntype.eq.0) then   ! PROBLEMA UNIDIMENSIONAL ELASTICO
        ndime=1              ! 1 Coordenada por Punto
        ndofn=1              ! 1 Grado de Libertad por Nodo
        nstre=1              ! 1 Tension por Punto
        nprop=5              ! 5 Propiedades por Elemento:
                             !      {
                             !       1-> Modulo de Young           (E)    ,
                             !       2-> Modulo de Poisson         (nu)   ,
                             !       3-> Seccion Transversal       (A)    ,
                             !       4-> Densidad                  (rho)  ,
                             !       5-> Coeficiente de Dilatacion (alpha)
                             !      }
        ntloa=3              ! 3 Tipos de Carga:
                             !      {
                             !       1-> Peso Propio                      ,
                             !       2-> Centrifuga
                             !       3-> Puntual                          ,
                             !      }
      else
        call SRV$Error('CAL$Ini_Par: Opcion no implementada (NTYPE)')
      endif

c.....Parametros de almacenamiento generales
      nevab=nnode*ndofn      ! Grados de libertad de elemento
      ntotv=npoin*ndofn      ! Grados de libertad de la estructura

c.....Puntos de integracion, funciones de forma y sus derivadas
      call CAL$Ini_FF(nnode,ndime,ngaue,posge,weige,shape,derie)

      return
      end

C____________________________________________________________________CAL$Ini_AMC
c     Obtencion de los punteros de los grados de libertad en el almacenamiento
c     por columnas con renumeracion opcional


      Subroutine CAL$Ini_AMC(nelem,nnode,npoin,ndofn,ntotv,lnods,lpamc,
     .           lpntn,nrigi)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension lnods(nnode,nelem),lpamc(ntotv),lpntn(npoin)

      do ipoir=1,npoin           ! Inicializacion
        itot0=(ipoir-1)*ndofn
        lpamc(itot0+1)=0
      enddo

      do ielem=1,nelem           ! Altura de las columnas del perfil
        do inode=1,nnode-1
          itot0=(lpntn(lnods(inode,ielem))-1)*ndofn
          do jnode=inode+1,nnode
            jtot0=(lpntn(lnods(jnode,ielem))-1)*ndofn
            if (itot0.gt.jtot0) then
              kdifr=(itot0-jtot0)
              if (kdifr.gt.lpamc(itot0+1)) lpamc(itot0+1)=kdifr
            else
              kdifr=(jtot0-itot0)
              if (kdifr.gt.lpamc(jtot0+1)) lpamc(jtot0+1)=kdifr
            endif
          enddo
        enddo
      enddo

      ipoir=1                    ! Acumulacion del Perfil
        itot0=(ipoir-1)*ndofn
        kdifr=lpamc(itot0+1)
        idofn=1
          itotv=itot0+idofn
          lpamc(itotv)=               kdifr+idofn
        do idofn=2,ndofn
          itotv=itot0+idofn
          lpamc(itotv)=lpamc(itotv-1)+kdifr+idofn
        enddo
      do ipoir=2,npoin
        itot0=(ipoir-1)*ndofn
        kdifr=lpamc(itot0+1)
        do idofn=1,ndofn
          itotv=itot0+idofn
          lpamc(itotv)=lpamc(itotv-1)+kdifr+idofn
        enddo
      enddo

      nrigi=lpamc(ntotv)         ! Posiciones de memoria que ocupa la matriz
      return
      end

c___________________________________________________________________CAL$Proc_Rig
c     Procesamiento de la matriz de rigidez del elemento IELEM
c     K_NIV = 0 -> Calculo de matrices D, B y S
c           = 1 -> Calculo de matrices D, B y S, Integracion y Ensamblaje

      Subroutine CAL$Proc_Rig(k_niv,ntype,ielem,nelem,nnode,npoin,ndime,
     .           ndofn,nstre,nmats,nprop,ngaue,nevab,ntotv,nrigi,lnods,
     .           coord,elcod,matno,props,weige,shape,derie,gcord,dvolu,
     .           xjace,djacb,xjaci,cartd,dmatx,bmatx,smatx,estif,lpamc,
     .           lpntn,rgmat)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension lnods(nnode,nelem),matno(nelem),props(nprop,nmats)

c.....Seleccion de coordenadas nodales del elemento
      call SRV$Sel_Subvector(nnode,lnods(1,ielem),ndime,npoin,coord,
     .     elcod)

c.....Selecion del tipo de material
      imats=matno(ielem)

c.....Generacion de matriz de rigidez del elemento
      call CAL$Gen_K(ntype,nnode,ndime,ndofn,nstre,nprop,ngaue,nevab,
     .     elcod,props(1,imats),weige,shape,derie,gcord,dvolu,xjace,
     .     djacb,xjaci,cartd,dmatx,bmatx,smatx)

      if (k_niv.eq.0) return

c.....Integracion de matriz de rigidez del elemento
      call CAL$IN_K(nstre,nevab,ngaue,dvolu,bmatx,smatx,estif)

c.....Ensamblaje de la matriz de rigidez de elemento en la general
      call CAL$Ensamb_K(nnode,ndofn,npoin,ntotv,nrigi,lnods(1,ielem),
     .     estif,lpamc,lpntn,rgmat)

      return
      end

c__________________________________________________________________CAL$Proc_FVol
c     Procesamiento de la fuerza volumetrica ILOAD del elemento IELEM
c     teniendo en cuenta los coeficientes de seguridad FLOAD(ncase,iload)
c     KLOAD(iload) = 1 -> peso propio [acel. gravedad(1:3) = RILOD(1:3,iload)]
c                    2 -> centrifuga  [vel. angular  (1:3) = RILOD(1:3,iload)]

      Subroutine CAL$Proc_FVol(ntype,ielem,iload,nelem,nnode,npoin,
     .           ndime,ndofn,nmats,nprop,ngaue,ntotv,nload,ncase,lnods,
     .           matno,props,shape,dvolu,bforc,gcord,kload,rilod,fload,
     .           eload,lpntn,gload)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension lnods(nnode,nelem),matno(nelem),props(nprop,nmats),
     .          kload(nload),rilod(3,nload),fload(ncase,nload),
     .          gload(ntotv,ncase)

      if (kload(iload).ne.1.and.kload(iload).ne.2) return

      cero=0.d+00

c.....Seleccion de tipo de material
      imats=matno(ielem)

c.....Generacion
      if      (kload(iload).eq.1) then                  ! Fuerzas Gravitatorias
        call CAL$Gen_FGra(ndofn,ndime,nprop,ngaue,
     .                    props(1,imats),rilod(1,iload),bforc)

      else if (kload(iload).eq.2) then                  ! Fuerzas Centrifugas
        call CAL$Gen_FCen(ntype,ndofn,ndime,nprop,ngaue,
     .                    props(1,imats),gcord,rilod(1,iload),bforc)
      endif

c.....Integracion
      call CAL$IN_FVol(nnode,ndofn,ngaue,dvolu,bforc,shape,eload)

c.....Ensamblaje
      do icase=1,ncase
        if (fload(icase,iload).ne.cero) then
          call CAL$Ensamb_F(nnode,ndofn,npoin,lnods(1,ielem),eload,
     .                      fload(icase,iload),lpntn,gload(1,icase))
        endif
      enddo

      return
      end

c__________________________________________________________________CAL$Proc_FPun
c     Procesamiento de la fuerza puntual ILOAD de valor PFORC(ndofn)
c     aplicada en el nodo IPOIN, teniendo en cuenta los coeficientes de
c     seguridad FLOAD(ncase,iload)

      Subroutine CAL$Proc_FPun(ipoin,iload,ndofn,npoin,ntotv,
     .           nload,ncase,pforc,fload,lpntn,gload)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension pforc(ndofn),fload(ncase,nload),gload(ntotv,ncase)

      cero=0.d+00

c.....ensamblaje
      do icase=1,ncase
        if (fload(icase,iload).ne.cero) then
          call CAL$Ensamb_F(    1,ndofn,npoin,ipoin,pforc,
     .                      fload(icase,iload),lpntn,gload(1,icase))
        endif
      enddo

      return
      end

c__________________________________________________________________CAL$Proc_Tens
c     Procesamiento de tensiones en los puntos de integracion del elemento IELEM

      Subroutine CAL$Proc_Tens(ielem,nelem,ncase,nnode,ndofn,npoin,
     .           nstre,nevab,ngaue,lnods,lpntn,xdisp,eldis,smatx,strsg)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension lnods(nnode,nelem),xdisp(ndofn,npoin,ncase),
     .          smatx(nstre,nevab,ngaue),strsg(nstre,ngaue,ncase)

      do icase=1,ncase

c.......selecciona desplazamientos del elemento
        call SRV$Sel_Subvector_2(nnode,lnods(1,ielem),lpntn,ndofn,npoin,
     .                           xdisp(1,1,icase),eldis)

c.......calculo de tensiones
        do igaue=1,ngaue
          call CAL$Cal_Tens(nstre,nevab,
     .                      smatx(1,1,igaue),eldis,strsg(1,igaue,icase))
        enddo

      enddo

      return
      end


c===============================================================================
c=============================== SEGUNDO NIVEL =================================
c===============================================================================

c______________________________________________________________________CAL$Gen_K
c     Generacion de la matriz de rigidez de elemento
   
      Subroutine CAL$Gen_K(ntype,nnode,ndime,ndofn,nstre,nprop,ngaue,
     .           nevab,elcod,props,weige,shape,derie,gcord,dvolu,xjace,
     .           djacb,xjaci,cartd,dmatx,bmatx,smatx)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension props(nprop),weige(ngaue),
     .          shape(nnode,ngaue),derie(ndime,nnode,ngaue),
     .          gcord(ndime,ngaue),dvolu(ngaue),
     .          xjace(ndime,ndime,ngaue),djacb(ngaue),
     .          xjaci(ndime,ndime,ngaue),cartd(ndime,nnode,ngaue),
     .          bmatx(nstre,nevab,ngaue),smatx(nstre,nevab,ngaue)

      young=props(1)
      poiss=props(2)
      thick=props(3)
      dense=props(4)

c.....calculo de la matriz D
      call CAL$Cal_D(ntype,young,poiss,nstre,dmatx)

c.....bucle sobre puntos de integracion gaussianos
      do igaue=1,ngaue

c.......calculo de posiciones de puntos de gauss
        call CAL$Interpol_FF(nnode,ndime,elcod,shape(1,igaue),
     .       gcord(1,igaue))

c.......calculo de jacobiano
        call CAL$Cal_Jac(nnode,ndime,elcod,derie(1,1,igaue),
     .       xjace(1,1,igaue))

c.......determinante e inversa del jacobiano
        call CAL$Cal_DIJac(ndime,xjace(1,1,igaue),djacb(igaue),
     .       xjaci(1,1,igaue))

c.......derivadas cartesianas
        call CAL$Cal_DerCart(nnode,ndime,xjaci(1,1,igaue),
     .       derie(1,1,igaue),cartd(1,1,igaue))

c.......generacion de la matriz B
        call CAL$Cal_B(ntype,nnode,ndime,ndofn,nstre,cartd(1,1,igaue),
     .       bmatx(1,1,igaue))

c.......calculo de la matriz S=DB
        call CAL$Cal_S(nstre,nevab,dmatx,bmatx(1,1,igaue),
     .       smatx(1,1,igaue))

c.......Calculo del volumen de integracion
        dvolu(igaue)=thick*djacb(igaue)*weige(igaue)

      enddo

      return
      end

c___________________________________________________________________CAL$Gen_FGra
c     Generacion, en los puntos de integracion, del vector de fuerzas volume-
c     tricas originadas por la aceleracion gravitatoria VGRAV

      Subroutine CAL$Gen_FGra(ndofn,ndime,nprop,ngaue,props,vgrav,bforc)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension props(nprop),vgrav(3),bforc(ndofn,ngaue)

      dense=props(4)
      do igaue=1,ngaue
        do idofn=1,ndofn
          bforc(idofn,igaue)=vgrav(idofn)*dense
        enddo
      enddo

      return
      end

c___________________________________________________________________CAL$Gen_FCen
c     Generacion, en los puntos de integracion, del vector de fuerzas volume-
c     tricas originadas por la rotacion en torno al origen con velocidad VGIRO.

      Subroutine CAL$Gen_FCen(ntype,ndofn,ndime,nprop,ngaue,
     .                        props,gcord,vgiro,bforc)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)
   
      dimension props(nprop),gcord(ndime,ngaue),bforc(ndofn,ngaue),
     .          vgiro(3)


      if (ntype.eq.0) then
        dense=props(4)
        w2=vgiro(2)*vgiro(2)+vgiro(3)*vgiro(3)
        do igaue=1,ngaue
          bforc(1,igaue)=w2*gcord(1,igaue)*dense
        enddo

      else
        call SRV$Error('CAL$Gen_FCen: Opcion no implementada (NTYPE)')
      endif

      return
      end


c===============================================================================
c=========================== INTEGRACION NUMERICA ==============================
c===============================================================================

c_______________________________________________________________________CAL$IN_K
c                                                                      T
c     Integracion numerica de la matriz de rigidez del elemento: K  = B  . S
c                                                                ~e   ~    ~
c     (se calcula solo la triangular inferior, por simetria)

      Subroutine CAL$IN_K(nstre,nevab,ngaue,dvolu,bmatx,smatx,estif)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension dvolu(ngaue),bmatx(nstre,nevab,ngaue),
     .          smatx(nstre,nevab,ngaue),estif(nevab,nevab)

      cero=0.d+00

      do ievab=1,nevab
        do jevab=ievab,nevab
          estif(jevab,ievab)=.0d+00
          do igaue=1,ngaue
            p=cero
            do istre=1,nstre
              p=p+bmatx(istre,ievab,igaue)*smatx(istre,jevab,igaue)
            enddo
            estif(jevab,ievab)=estif(jevab,ievab)+dvolu(igaue)*p
          enddo
        enddo
      enddo

      return
      end

c____________________________________________________________________CAL$IN_FVol
c                                                        _      T   _
c     Integracion numerica de las fuerzas volumetricas:  f =   N  . b
c                                                              ~

      Subroutine CAL$IN_FVol(nnode,ndofn,ngaue,dvolu,bforc,shape,eload)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension dvolu(ngaue),bforc(ndofn,ngaue),shape(nnode,ngaue),
     .          eload(ndofn,nnode)

      cero=0.d+00

      do inode=1,nnode
        do idofn=1,ndofn
          eload(idofn,inode)=cero
          do igaue=1,ngaue
            eload(idofn,inode)=eload(idofn,inode)+dvolu(igaue)*
     .                         bforc(idofn,igaue)*shape(inode,igaue)
          enddo
        enddo
      enddo

      return
      end

c===============================================================================
c============================ FUNCIONES DE FORMA ===============================
c===============================================================================

c_____________________________________________________________________CAL$Ini_FF
c     Inicializacion de los puntos (POSGE) y pesos (WEIGE) de integracion,
c     y precalculo de las funciones de forma (SHAPE) y sus derivadas (DERIE)

      Subroutine CAL$Ini_FF(nnode,ndime,ngaue,posge,weige,shape,derie)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension posge(ndime,ngaue),weige(ngaue),
     .          shape(nnode,ngaue),derie(ndime,nnode,ngaue)

c.....Puntos y Pesos de Integracion
      if (ndime.eq.1) then                                    ! PROBLEMAS 1D

        if (ngaue.eq.2) then
          posge(1,1)=-.577350269189626d+00
          weige(  1)= .100000000000000d+01
          posge(1,2)= .577350269189626d+00
          weige(  2)= .100000000000000d+01

        else if (ngaue.eq.3) then
          posge(1,1)=-.774596669241483d+00
          weige(  1)= .555555555555556d+00
          posge(1,2)= .000000000000000d+00
          weige(  2)= .888888888888889d+00
          posge(1,3)= .774596669241483d+00
          weige(  3)= .555555555555556d+00

        else
          call SRV$Error('CAL$Ini_FF: Opcion no implementada (NGAUE)')
        endif

      else
        call SRV$Error('CAL$Ini_FF: Opcion no implementada (NDIME)')
      endif

c.....valoracion de funciones de forma en puntos de integracion
      do igaue=1,ngaue
        call CAL$Val_FF(nnode,ndime,posge(1,igaue),
     .                  shape(  1,igaue),derie(1,1,igaue))
      enddo

      return
      end

c_____________________________________________________________________CAL$Val_FF
c     Valoracion de las funciones de forma (SHAPE) y sus derivadas (DERIE)
c     en un punto de coordenadas conocidas (CORIS) en el sistema local.

      Subroutine CAL$Val_FF(nnode,ndime,coris,shape,derie)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension coris(ndime),shape(nnode),derie(ndime,nnode)

c-----PROBLEMAS UNIDIMENSIONALES
      if (ndime.eq.1) then
        s=coris(1)

        if (nnode.eq.2) then                            ! Elemento 1D2N/LSEG1
          shape(  1)=  (1.d+00-s)/2.d+00
          derie(1,1)= -.5d+00
          shape(  2)=  (1.d+00+s)/2.d+00
          derie(1,2)=  .5d+00

        else if (nnode.eq.3) then                       ! Elemento 1D3N/LSEG2
          shape(  1)=s*(s-1.d+00)/2.d+00
          derie(1,1)=s-.5d+00
          shape(  2)=s*(s+1.d+00)/2.d+00
          derie(1,2)=s+.5d+00
          shape(  3)=1.d+00-s*s
          derie(1,3)=-2.d+00*s

        else
          call SRV$Error('CAL$Val_FF: Opcion no implementada (NNODE)')
        endif

      else
        call SRV$Error('CAL$Val_FF: Opcion no implementada (NDIME)')
      endif

      return
      end

c________________________________________________________________CAL$Interpol_FF
c     Calculo de valores interpolados (PINTR) a partir de valores
c     nodales (PNODE) mediante las funciones de forma (SHAPE)

      Subroutine CAL$Interpol_FF(nnode,nd,pnode,shape,pintr)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

c.....Producto: PINTR(nd)=PNODE(nd,nnode)*SHAPE(nnode)
      call SRV$Mult_M_M(nd,nnode,1,pnode,shape,pintr)

      return
      end

c===============================================================================
c============================ FORMULACION MEF ==================================
c===============================================================================

c____________________________________________________________________CAL$Cal_Jac
c     Calculo de matriz jacobiana

      Subroutine CAL$Cal_Jac(nnode,ndime,elcod,derie,xjace)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

c.....Producto: XJACE(ndime,ndime)=ELCOD(ndime,nnode)*DERIE(ndime,nnode)
      call SRV$Mult_M_Mt(ndime,nnode,ndime,elcod,derie,xjace)

      return
      end

c__________________________________________________________________CAL$Cal_DIJac
c     Calculo del determinante y la inversa de la matriz jacobiana

      Subroutine CAL$Cal_DIJac(ndime,xjace,djacb,xjaci)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension xjace(ndime,ndime),xjaci(ndime,ndime)

      if(ndime.eq.1) then
        djacb=xjace(1,1)
        xjaci(1,1)=1.d+00/djacb
        djacb=dabs(djacb)
      else
        call SRV$Error('CAL$Cal_DIJac: Opcion no implementada (NDIME)')
      endif

      return
      end

c________________________________________________________________CAL$Cal_DerCart
c     Calculo de las derivadas cartesianas

      Subroutine CAL$Cal_DerCart(nnode,ndime,xjaci,derie,cartd)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

c.....Producto: CARTD(ndime,nnode)=XJACI(ndime,ndime)*DERIE(ndime,nnode)
      call SRV$Mult_Mt_M(ndime,ndime,nnode,xjaci,derie,cartd)

      return
      end

c______________________________________________________________________CAL$Cal_D
c     Calculo de la matriz D del material
c                          ~

      Subroutine CAL$Cal_D(ntype,young,poiss,nstre,dmatx)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension dmatx(nstre,nstre)

      if (ntype.eq.0) then
        dmatx(1,1)=young
      else
        call SRV$Error('CAL$Cal_D: Opcion no implementada (NTYPE)')
      endif

      return
      end

c______________________________________________________________________CAL$Cal_B
c     Calculo de la matriz B = L N
c                          ~   ~ ~

      Subroutine CAL$Cal_B(ntype,nnode,ndime,ndofn,nstre,cartd,bmatx)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension cartd(ndime,nnode),bmatx(nstre,ndofn,nnode)

      if (ntype.eq.0) then
        do inode=1,nnode
          bmatx(1,1,inode)=cartd(1,inode)
        enddo
      else
        call SRV$Error('CAL$CAL$Cal_B: Opcion no implementada (NTYPE)')
      endif

      return
      end

c______________________________________________________________________CAL$Cal_S
c     Calculo de la matriz S = D . B
c                          ~   ~   ~

      Subroutine CAL$Cal_S(nstre,nevab,dmatx,bmatx,smatx)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

c.....Producto: SMATX(nstre,nevab)=DMATX(nstre,nstre)*BMATX(nstre,nevab)
      call SRV$Mult_M_M(nstre,nstre,nevab,dmatx,bmatx,smatx)

      return
      end

c___________________________________________________________________CAL$Cal_Tens
c                              _     _
c     Calculo de las tensiones s = S a
c                                  ~

      Subroutine CAL$Cal_Tens(nstre,nevab,smatx,eldis,strsg)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension smatx(nstre,nevab),strsg(nstre)

c.....Producto: STRSG(nstre)=SMATX(nstre,nevab)*ELDIS(nevab)
      call SRV$Mult_M_M(nstre,nevab,1,smatx,eldis,strsg)

      return
      end


c===============================================================================
c============================= GESTION DE MEMORIA ==============================
c===============================================================================

c___________________________________________________________________CAL$Ensamb_K
c     Ensamblaje de la matriz de rigidez de elemento (ESTIF)
c     en la matriz de rigidez global (RGMAT) en perfil (LPAMC)
c     considerando la renumeracion nodal (LPNTN)

      Subroutine CAL$Ensamb_K(nnode,ndofn,npoin,ntotv,nrigi,
     .                        lnods,estif,lpamc,lpntn,rgmat)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)
   
      dimension lnods(nnode),estif(ndofn,nnode,ndofn,nnode),
     .          lpamc(ntotv),lpntn(npoin),rgmat(nrigi)

      do inode=1,nnode
        ipoir=lpntn(lnods(inode))
        do idofn=1,ndofn
          itotv=(ipoir-1)*ndofn+idofn

          jnode=inode
            jpoir=ipoir
            do jdofn=idofn,ndofn
              jtotv=(jpoir-1)*ndofn+jdofn
                lpjj=iabs(lpamc(jtotv))
                lpj0=lpjj-jtotv
                lp=lpj0+itotv
              rgmat(lp)=rgmat(lp)+estif(jdofn,jnode,idofn,inode)
            enddo

          do jnode=inode+1,nnode
            jpoir=lpntn(lnods(jnode))
            do jdofn=1,ndofn
              jtotv=(jpoir-1)*ndofn+jdofn
              if (jtotv.le.itotv) then
                lpii=iabs(lpamc(itotv))
                lpi0=lpii-itotv
                lp=lpi0+jtotv
              else
                lpjj=iabs(lpamc(jtotv))
                lpj0=lpjj-jtotv
                lp=lpj0+itotv
              endif
              rgmat(lp)=rgmat(lp)+estif(jdofn,jnode,idofn,inode)
            enddo
          enddo

        enddo
      enddo

      return
      end


c___________________________________________________________________CAL$Ensamb_F
c     Ensamblaje del vector de fuerzas de elemento (ELOAD)
c     en el vector de fuerzas global (GLOAD)
c     considerando la renumeracion nodal (LPNTN)
c     y teniendo en cuenta el coeficiente de seguridad (FLOAD)

      Subroutine CAL$Ensamb_F(nnode,ndofn,npoin,
     .                        lnods,eload,fload,lpntn,gload)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)
   
      dimension lnods(nnode),eload(ndofn,nnode),lpntn(npoin),
     .          gload(ndofn,npoin)

      do inode=1,nnode
        ipoir=lpntn(lnods(inode))
        do idofn=1,ndofn
          gload(idofn,ipoir)=gload(idofn,ipoir)+
     .                       eload(idofn,inode)*fload
        enddo
      enddo

      return
      end

C_________________________________________________________________CAL$Imp_Vincul
c     Formacion del vector de grados de libertad prescritos (NGFIX)
c     a partir de las condiciones de vinculacion (NOFIX)
c     considerando la renumeracion nodal (LPNTN)

      Subroutine CAL$Imp_Vincul(nvfix,ndofn,npoin,nofix,ngfix,lpntn)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension nofix(2,nvfix),ngfix(nvfix),lpntn(npoin)

      do ivfix=1,nvfix
        ipoin=nofix(1,ivfix)
        idofn=nofix(2,ivfix)
        ipoir=lpntn(ipoin)
        itotv=(ipoir-1)*ndofn+idofn
        ngfix(ivfix)=itotv
      enddo

      return
      end


c===============================================================================
c======================== Subrutinas Generales de Servicio =====================
c===============================================================================

c_________________________________________________________________SRV$Ini_Vector
c     Inicializacion de un vector a cero

      Subroutine SRV$Ini_Vector(n,v)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension v(n)

      cero=0.d+00

      do i=1,n
        v(i)=cero
      enddo

      return
      end

c______________________________________________________________SRV$Sel_Subvector
c     Seleccion de un subvector mediante un redireccionamiento

      Subroutine SRV$Sel_Subvector(nn,ln,nd,np,co,el)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension ln(nn),co(nd,np),el(nd,nn)

      do in=1,nn
        ip=ln(in)
        do id=1,nd
          el(id,in)=co(id,ip)
        enddo
      enddo

      return
      end

c____________________________________________________________SRV$Sel_Subvector_2
c     Seleccion de un subvector mediante dos redireccionamientos consecutivos

      Subroutine SRV$Sel_Subvector_2(nn,ln,ln2,nd,np,co,el)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension ln(nn),ln2(np),co(nd,np),el(nd,nn)

      do in=1,nn
        ip=ln2(ln(in))
        do id=1,nd
          el(id,in)=co(id,ip)
        enddo
      enddo

      return
      end


c___________________________________________________________________SRV$Mult_M_M
c     Producto Matricial C         = A         . B
c                        ~ (n x m)   ~ (n x l)   ~ (l x m)

      Subroutine SRV$Mult_M_M(n,l,m,a,b,c)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension a(n,l),b(l,m),c(n,m)

      cero=0.d+00

      do j=1,m
        do i=1,n
          c(i,j)=cero
          do k=1,l
            c(i,j)=c(i,j)+a(i,k)*b(k,j)
          enddo  
        enddo
      enddo

      return
      end

c__________________________________________________________________SRV$Mult_M_Mt
c                                                  T
c     Producto Matricial C         = A         . B
c                        ~ (n x m)   ~ (n x l)   ~ (m x l)

      Subroutine SRV$Mult_M_Mt(n,l,m,a,b,c)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension a(n,l),b(m,l),c(n,m)

      cero=0.d+00

      do j=1,m
        do i=1,n
          c(i,j)=cero
          do k=1,l
            c(i,j)=c(i,j)+a(i,k)*b(j,k)
          enddo  
        enddo
      enddo

      return
      end

c__________________________________________________________________SRV$Mult_Mt_M
c                                      T
c     Producto Matricial C         = A         . B
c                        ~ (n x m)   ~ (l x n)   ~ (l x m)

      Subroutine SRV$Mult_Mt_M(n,l,m,a,b,c)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension a(l,n),b(l,m),c(n,m)

      cero=0.d+00

      do j=1,m
        do i=1,n
          c(i,j)=cero
          do k=1,l
            c(i,j)=c(i,j)+a(k,i)*b(k,j)
          enddo  
        enddo
      enddo

      return
      end

c______________________________________________________________________SRV$Error
c     Impresion de mensaje de error y detencion del programa

      Subroutine SRV$Error(mensaje)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      character*(*) mensaje

      write(6,100) mensaje
 100  format(/' ERROR > ',a/)

      stop
      end
