path: root/Trivac/src/BIVCOL.f

*DECK BIVCOL
      SUBROUTINE BIVCOL (IPTRK,IMPX,IELEM,ICOL)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Selection of the unit matrices (mass, stiffness, etc.) for a finite
* element approximation.
*
*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
* IPTRK   L_TRACK pointer to the tracking information.
* IMPX    print parameter.
* IELEM   degree of the finite elements: =1: (linear polynomials);
*         =2: (parabolic polynomials);   =3: (cubic polynomials);
*         =4: (quartic polynomials).
* ICOL    type of quadrature used to integrate the mass matrices:
*         =1: (analytic integration);  =2: (Gauss-Lobatto quadrature)
*         =3: (Gauss-Legendre quadrature).
*         IELEM=1 with ICOL=2 is equivalent to finite differences.
*
*-----------------------------------------------------------------------
*
      USE GANLIB
*----
*  SUBROUTINE ARGUMENTS
*----
      TYPE(C_PTR) IPTRK
      INTEGER IMPX,IELEM,ICOL
*----
*  LOCAL VARIABLES
*----
      CHARACTER*40 HTYPE
      DOUBLE PRECISION DSUM
      REAL EL(2,2),TL(2),TSL(2),RL(2,2),RSL(2,2),QSL(2,2),TSL1(2),
     1 TSL2(2),RL2(2,2),RSL2(2,2)
      REAL RHA6(6,6),QHA6(6,6),RHL6(6,6),QHL6(6,6),RTA(3,3),QTA(3,3),
     1 RTL(3,3),QTL(3,3)
      REAL EP(3,3),TP(3),TSP(3),RP(3,3),VP(3,2),HP(3,2),RSP(3,3),
     1 QSP(3,3),EP1(3,3),TP1(3),TSP1(3),VP1(3,2),HP1(3,2),QSP1(3,3),
     2 EP2(3,3),TP2(3),TSP2(3),RP2(3,3),VP2(3,2),HP2(3,2),RSP2(3,3),
     3 QSP2(3,3)
      REAL EC(4,4),TC(4),TSC(4),RC(4,4),VC(4,3),HC(4,3),RSC(4,4),
     1 QSC(4,4),EC1(4,4),TC1(4),TSC1(4),VC1(4,3),HC1(4,3),QSC1(4,4),
     2 EC2(4,4),TC2(4),TSC2(4),RC2(4,4),VC2(4,3),HC2(4,3),RSC2(4,4),
     3 QSC2(4,4)
      REAL EQ(5,5),VQ(5,4),HQ(5,4),TQ(5),TSQ(5),QSQ(5,5)
      REAL RLQ(2,2),RLR(2,2),RL1Q(2,2),RL1R(2,2),RL2Q(2,2),RL2R(2,2)
*-----------------------------------------------------------------------
* THE BIVAC REFERENCE ELEMENT IS DEFINED BETWEEN POINTS -1/2 AND +1/2.
* THE COLLOCATION POLYNOMIALS CORRESPONDING TO APPROXIMATIONS LL2$, PL3,
* PL3$, PL3#, CL4, CL4$, CL4# AND QL5$ ARE PARTIALLY OR COMPLETELY
* ORTHONORMALIZED IN ORDER TO PRODUCE A SPARSE MASS MATRIX.
*-----------------------------------------------------------------------
*
******************* FINITE ELEMENT BASIC MATRICES **********************
*                                                                      *
      REAL T(5),TS(5),R(25),RS(25),Q(25),QS(25),V(20),H(20),E(25),
     1 RH(6,6),QH(6,6),RT(3,3),QT(3,3),R1DQ(4),R1DR(4)
*                                                                      *
* LC      : NUMBER OF POLYNOMIALS IN A COMPLETE 1-D BASIS.             *
* T       : CARTESIAN LINEAR PRODUCT VECTOR.                           *
* TS      : CYLINDRICAL LINEAR PRODUCT VECTOR.                         *
* R       : CARTESIAN MASS MATRIX.                                     *
* RS      : CYLINDRICAL MASS MATRIX.                                   *
* Q       : CARTESIAN STIFFNESS MATRIX.                                *
* QS      : CYLINDRICAL STIFFNESS MATRIX.                              *
* V       : NODAL COUPLING MATRIX.                                     *
* H       : PIOLAT (HEXAGONAL) COUPLING MATRIX.                        *
* E       : POLYNOMIAL COEFFICIENTS.                                   *
* RH      : HEXAGONAL MASS MATRIX.                                     *
* QH      : HEXAGONAL STIFFNESS MATRIX.                                *
* RT      : TRIANGULAR MASS MATRIX.                                    *
* QT      : TRIANGULAR STIFFNESS MATRIX.                               *
* R1DQ    : SPHERICAL MASS MATRIX.                                     *
* R1DR    : SPHERICAL MASS MATRIX.                                     *
*                                                                      *
************************************************************************
*
*----
*  LINEAR LAGRANGIAN POLYNOMIALS.
*----
      DATA EL/0.5,-1.0,0.5,1.0/
      DATA TL/0.5,0.5/
      DATA RL/
     $ 0.333333333333, 0.166666666667, 0.166666666667, 0.333333333333/
*     CYLINDRICAL OPTION MATRICES:
      DATA TSL/-0.083333333333, 0.083333333333/
      DATA RSL/
     $-0.083333333333, 0.000000000000, 0.000000000000, 0.083333333333/
      DATA QSL/
     $ 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000/
*     SPHERICAL OPTION MATRICES (ANALYTIC INTEGRATION):
      DATA RLQ/-.083333333333, 0.0, 0.0, 0.083333333333/
      DATA RLR/
     $ 0.033333333333, 0.008333333333, 0.008333333333, 0.033333333333/
*     GAUSS-LOBATTO (FINITE DIFFERENCES) MATRICES:
      DATA TSL1/-0.25,0.25/
*     SPHERICAL OPTION MATRICES (GAUSS-LOBATTO):
      DATA RL1Q/-0.166666666667, 0.0, 0.0, 0.166666666667/
      DATA RL1R/0.041666666667, 0.0, 0.0, 0.041666666667/
*     GAUSS-LEGENDRE (SUPERCONVERGENT) MATRICES:
      DATA TSL2/0.0,0.0/
      DATA RL2/0.25,0.25,0.25,0.25/
      DATA RSL2/0.0,0.0,0.0,0.0/
*     SPHERICAL OPTION MATRICES (SUPERCONVERGENT):
      DATA RL2Q/
     $ -0.03472222222,-0.0069444444444,-0.0069444444444, 0.048611111111/
      DATA RL2R/
     $  0.00925925926, 0.0185185185185, 0.0185185185185, 0.037037037037/
*
*     ANALYTIC INTEGRATION FOR HEXAGON AND TRIANGLE.
      DATA RHA6/
     > 0.158470 , 0.086580 , 0.036760 , 0.027870 , 0.036760 , 0.086580,
     > 0.086580 , 0.158470 , 0.086580 , 0.036760 , 0.027870 , 0.036760,
     > 0.036760 , 0.086580 , 0.158470 , 0.086580 , 0.036760 , 0.027870,
     > 0.027870 , 0.036760 , 0.086580 , 0.158470 , 0.086580 , 0.036760,
     > 0.036760 , 0.027870 , 0.036760 , 0.086580 , 0.158470 , 0.086580,
     > 0.086580 , 0.036760 , 0.027870 , 0.036760 , 0.086580 , 0.158470/
      DATA QHA6/
     > 0.760640 ,-0.161980 ,-0.169310 ,-0.098060 ,-0.169310 ,-0.161980,
     >-0.161980 , 0.760640 ,-0.161980 ,-0.169310 ,-0.098060 ,-0.169310,
     >-0.169310 ,-0.161980 , 0.760640 ,-0.161980 ,-0.169310 ,-0.098060,
     >-0.098060 ,-0.169310 ,-0.161980 , 0.760640 ,-0.161980 ,-0.169310,
     >-0.169310 ,-0.098060 ,-0.169310 ,-0.161980 , 0.760640 ,-0.161980,
     >-0.161980 ,-0.169310 ,-0.098060 ,-0.169310 ,-0.161980 , 0.760640/
      DATA RTA/
     > 1.0, 0.5, 0.5, 0.5, 1.0, 0.5, 0.5, 0.5, 1.0/
      DATA QTA/
     > 1.0,-0.5,-0.5,-0.5, 1.0,-0.5,-0.5,-0.5, 1.0/
*
*     GAUSS-LOBATTO INTEGRATION FOR HEXAGON AND TRIANGLE.
      DATA RHL6/
     > 1.000000 , 0.000000 , 0.000000 , 0.000000 , 0.000000 , 0.000000,
     > 0.000000 , 1.000000 , 0.000000 , 0.000000 , 0.000000 , 0.000000,
     > 0.000000 , 0.000000 , 1.000000 , 0.000000 , 0.000000 , 0.000000,
     > 0.000000 , 0.000000 , 0.000000 , 1.000000 , 0.000000 , 0.000000,
     > 0.000000 , 0.000000 , 0.000000 , 0.000000 , 1.000000 , 0.000000,
     > 0.000000 , 0.000000 , 0.000000 , 0.000000 , 0.000000 , 1.000000/
      DATA QHL6/
     > 1.666667 ,-1.000000 , 0.166667 ,-1.000000 , 0.166667 , 0.000000,
     >-1.000000 , 1.666667 ,-1.000000 , 0.166667 , 0.000000 , 0.166667,
     > 0.166667 ,-1.000000 , 1.666667 , 0.000000 , 0.166667 ,-1.000000,
     >-1.000000 , 0.166667 , 0.000000 , 1.666667 ,-1.000000 , 0.166667,
     > 0.166667 , 0.000000 , 0.166667 ,-1.000000 , 1.666667 ,-1.000000,
     > 0.000000 , 0.166667 ,-1.000000 , 0.166667 ,-1.000000 , 1.666667/
      DATA RTL/
     > 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0/
      DATA QTL/
     > 1.0,-0.5,-0.5,-0.5, 1.0,-0.5,-0.5,-0.5, 1.0/
*----
*  PARABOLIC LAGRANGIAN POLYNOMIALS.
*----
      DATA EP/
     $-0.125000000000,-1.000000000000, 2.500000000000,
     $ 1.250000000000, 0.000000000000,-5.000000000000,
     $-0.125000000000, 1.000000000000, 2.500000000000/
      DATA TP/
     $ 0.083333333333, 0.833333333333, 0.083333333333/
      DATA RP/
     $ 0.125000000000, 0.000000000000,-0.041666666667,
     $ 0.000000000000, 0.833333333333, 0.000000000000,
     $-0.041666666667, 0.000000000000, 0.125000000000/
      DATA VP/
     $-1.000000000000, 0.000000000000, 1.000000000000,
     $ 1.443375672974,-2.886751345948, 1.443375672974/
      DATA HP/
     $ 0.083333333333, 0.833333333333, 0.083333333333,
     $-0.288675134595, 0.000000000000, 0.288675134595/
*     CYLINDRICAL OPTION MATRICES:
      DATA TSP/
     $-0.083333333333, 0.000000000000, 0.083333333333/
      DATA RSP/
     $-0.041666666667,-0.041666666667, 0.000000000000,
     $-0.041666666667, 0.000000000000, 0.041666666667,
     $ 0.000000000000, 0.041666666667, 0.041666666667/
      DATA QSP/
     $-0.833333333333, 0.833333333333, 0.000000000000,
     $ 0.833333333333, 0.000000000000,-0.833333333333,
     $ 0.000000000000,-0.833333333333, 0.833333333333/
*     GAUSS-LOBATTO (VARIATIONAL COLLOCATION METHOD) MATRICES:
      DATA EP1/0.0,-1.0,2.0,1.0,0.0,-4.0,0.0,1.0,2.0/
      DATA TP1/
     $ 0.166666666667, 0.666666666667, 0.166666666667/
      DATA VP1/
     $-1.000000000000, 0.0000000000000, 1.000000000000,
     $ 1.154700538379,-2.3094010767585, 1.154700538379/
      DATA HP1/
     $ 0.166666666667, 0.666666666667, 0.166666666667,
     $-0.288675134595, 0.000000000000, 0.288675134595/
*     CYLINDRICAL GAUSS-LOBATTO (VARIATIONAL COLLOCATION METHOD)
*     MATRICES.
      DATA TSP1/
     $-0.083333333333, 0.000000000000, 0.083333333333/
      DATA QSP1/
     $-0.666666666667, 0.666666666667, 0.000000000000,
     $ 0.666666666667, 0.000000000000,-0.666666666667,
     $ 0.000000000000,-0.666666666667, 0.666666666667/
*     GAUSS-LEGENDRE (SUPERCONVERGENT) MATRICES:
      DATA EP2/
     $-0.250000000000,-1.000000000000, 3.000000000000,
     $ 1.500000000000, 0.000000000000,-6.000000000000,
     $-0.250000000000, 1.000000000000, 3.000000000000/
      DATA TP2/
     $ 0.000000000000, 1.000000000000, 0.000000000000/
      DATA RP2/
     $ 0.083333333333, 0.000000000000,-0.083333333333,
     $ 0.000000000000, 1.000000000000, 0.000000000000,
     $-0.083333333333, 0.000000000000, 0.083333333333/
      DATA VP2/
     $-1.000000000000, 0.000000000000, 1.000000000000,
     $ 1.732050807569,-3.464101615138, 1.732050807569/
      DATA HP2/
     $ 0.000000000000, 1.000000000000, 0.000000000000,
     $-0.288675134595, 0.000000000000, 0.288675134595/
*     CYLINDRICAL GAUSS-LEGENDRE (SUPERCONVERGENT) MATRICES:
      DATA TSP2/
     $-0.083333333333, 0.000000000000, 0.083333333333/
      DATA RSP2/
     $ 0.000000000000,-0.083333333333, 0.000000000000,
     $-0.083333333333, 0.000000000000, 0.083333333333,
     $ 0.000000000000, 0.083333333333, 0.000000000000/
      DATA QSP2/
     $-1.000000000000, 1.000000000000, 0.000000000000,
     $ 1.000000000000, 0.000000000000,-1.000000000000,
     $ 0.000000000000,-1.000000000000, 1.000000000000/
*----
*  CUBIC LAGRANGIAN POLYNOMIALS.
*----
      DATA EC/
     $-0.125000000000, 0.750000000000, 2.500000000000, -7.000000000000,
     $ 0.625000000000,-3.307189138831,-2.500000000000, 13.228756555323,
     $ 0.625000000000, 3.307189138831,-2.500000000000,-13.228756555323,
     $-0.125000000000,-0.750000000000, 2.500000000000,  7.000000000000/
      DATA TC/
     $ 0.083333333333, 0.416666666667, 0.416666666667, 0.083333333333/
      DATA RC/
     $ 0.066666666667, 0.000000000000, 0.000000000000, 0.016666666667,
     $ 0.000000000000, 0.416666666667, 0.000000000000, 0.000000000000,
     $ 0.000000000000, 0.000000000000, 0.416666666667, 0.000000000000,
     $ 0.016666666667, 0.000000000000, 0.000000000000, 0.066666666667/
      DATA VC/
     $-1.000000000000, 0.000000000000, 0.000000000000,1.000000000000,
     $ 1.443375672974,-1.443375672974,-1.443375672974,1.443375672974,
     $-1.565247584250, 2.958039891550,-2.958039891550,1.565247584250/
      DATA HC/
     $ 0.083333333333, 0.416666666667, 0.416666666667, 0.083333333333,
     $-0.086602540378,-0.381881307913, 0.381881307913, 0.086602540378,
     $ 0.186338998125,-0.186338998125,-0.186338998125, 0.186338998125/
*     CYLINDRICAL OPTION MATRICES:
      DATA TSC /
     $-0.025000000000,-0.110239637961, 0.110239637961, 0.025000000000/
      DATA RSC/
     $-0.025000000000,-0.015748519709, 0.015748519709, 0.000000000000,
     $-0.015748519709,-0.078742598544, 0.000000000000,-0.015748519709,
     $ 0.015748519709, 0.000000000000, 0.078742598544, 0.015748519709,
     $ 0.000000000000,-0.015748519709, 0.015748519709, 0.025000000000/
      DATA QSC/
     $-2.000000000000, 2.102396379610,-0.102396379610, 0.000000000000,
     $ 2.102396379610,-2.204792759220, 0.000000000000, 0.102396379610,
     $-0.102396379610, 0.000000000000, 2.204792759220,-2.102396379610,
     $ 0.000000000000, 0.102396379610,-2.102396379610, 2.000000000000/
*     GAUSS-LOBATTO (VARIATIONAL COLLOCATION METHOD) MATRICES:
      DATA EC1/
     $-0.125000000000, 0.250000000000, 2.500000000000, -5.000000000000,
     $ 0.625000000000,-2.795084971875,-2.500000000000, 11.180339887499,
     $ 0.625000000000, 2.795084971875,-2.500000000000,-11.180339887499,
     $-0.125000000000,-0.250000000000, 2.500000000000,  5.000000000000/
      DATA TC1/
     $ 0.083333333333, 0.416666666667, 0.416666666667, 0.083333333333/
      DATA VC1/
     $-1.000000000000, 0.000000000000, 0.000000000000,1.000000000000,
     $ 1.443375672974,-1.443375672974,-1.443375672974,1.443375672974,
     $-1.118033988750, 2.500000000000,-2.500000000000,1.118033988750/
      DATA HC1/
     $ 0.083333333333, 0.416666666667, 0.416666666667, 0.083333333333,
     $-0.144337567297,-0.322748612184, 0.322748612184, 0.144337567297,
     $ 0.186338998125,-0.186338998125,-0.186338998125, 0.186338998125/
*     CYLINDRICAL GAUSS-LOBATTO (VARIATIONAL COLLOCATION METHOD)
*     MATRICES:
      DATA TSC1/
     $-0.041666666667,-0.093169499062, 0.093169499062, 0.041666666667/
      DATA QSC1/
     $-1.666666666667, 1.765028323958,-0.098361657292, 0.000000000000,
     $ 1.765028323958,-1.863389981250, 0.000000000000, 0.098361657292,
     $-0.098361657292, 0.000000000000, 1.863389981250,-1.765028323958,
     $ 0.000000000000, 0.098361657292,-1.765028323958, 1.666666666667/
*     GAUSS-LEGENDRE (SUPERCONVERGENT) MATRICES:
      DATA EC2/
     $-0.125000000000, 1.500000000000, 2.500000000000,-10.000000000000,
     $ 0.625000000000,-3.952847075210,-2.500000000000, 15.811388300842,
     $ 0.625000000000, 3.952847075210,-2.500000000000,-15.811388300842,
     $-0.125000000000,-1.500000000000, 2.500000000000, 10.000000000000/
      DATA TC2/
     $ 0.083333333333, 0.416666666667, 0.416666666667, 0.083333333333/
      DATA RC2/
     $ 0.041666666667, 0.000000000000, 0.000000000000, 0.041666666667,
     $ 0.000000000000, 0.416666666667, 0.000000000000, 0.000000000000,
     $ 0.000000000000, 0.000000000000, 0.416666666667, 0.000000000000,
     $ 0.041666666667, 0.000000000000, 0.000000000000, 0.041666666667/
      DATA VC2/
     $-1.000000000000, 0.000000000000, 0.000000000000,1.000000000000,
     $ 1.443375672974,-1.443375672974,-1.443375672974,1.443375672974,
     $-2.236067977500, 3.535533905933,-3.535533905933,2.236067977500/
      DATA HC2/
     $ 0.083333333333, 0.416666666667, 0.416666666667, 0.083333333333,
     $ 0.000000000000,-0.456435464588, 0.456435464588, 0.000000000000,
     $ 0.186338998125,-0.186338998125,-0.186338998125, 0.186338998125/
*     CYLINDRICAL GAUSS-LEGENDRE (SUPERCONVERGENT) MATRICES:
      DATA TSC2/
     $ 0.000000000000,-0.131761569174, 0.131761569174, 0.000000000000/
      DATA RSC2/
     $ 0.000000000000,-0.032940392293, 0.032940392293, 0.000000000000,
     $-0.032940392293,-0.065880784587, 0.000000000000,-0.032940392293,
     $ 0.032940392293, 0.000000000000, 0.065880784587, 0.032940392293,
     $ 0.000000000000,-0.032940392293, 0.032940392293, 0.000000000000/
      DATA QSC2/
     $-2.500000000000, 2.567615691737,-0.067615691737, 0.000000000000,
     $ 2.567615691737,-2.635231383474, 0.000000000000, 0.067615691737,
     $-0.067615691737, 0.000000000000, 2.635231383474,-2.567615691737,
     $ 0.000000000000, 0.067615691737,-2.567615691737, 2.500000000000/
*----
*  QUARTIC LAGRANGIAN POLYNOMIALS.
*----
      DATA EQ/
     $  0.000000000000, 0.750000000000,-1.500000000000,-7.000000000000,
     $ 14.000000000000, 0.000000000000,-2.673169155391, 8.166666666667,
     $ 10.692676621564,-32.666666666667, 1.000000000000, 0.000000000000,
     $-13.333333333333, 0.000000000000,37.333333333333, 0.000000000000,
     $  2.673169155391,8.166666666667,-10.692676621564,-32.666666666667,
     $  0.000000000000,-0.750000000000,-1.500000000000, 7.000000000000,
     $ 14.000000000000/
      DATA TQ/
     $ 0.050000000000, 0.272222222222, 0.355555555556, 0.272222222222,
     $ 0.050000000000/
      DATA VQ/
     $-1.000000000000, 0.00000000000, 0.000000000000, 0.00000000000,
     $ 1.000000000000, 1.55884572681,-0.943005439677,-1.23168057427,
     $-0.943005439677, 1.55884572681,-1.565247584250, 2.39095517873,
     $ 0.000000000000,-2.39095517873, 1.565247584250, 1.058300524426,
     $-2.46936789032 , 2.8221347318 ,-2.46936789032 , 1.058300524426/
      DATA HQ/
     $ 0.050000000000, 0.272222222222, 0.355555555556, 0.272222222222,
     $ 0.050000000000,-0.086602540378,-0.308670986291, 0.000000000000,
     $ 0.308670986291, 0.086602540378, 0.111803398875, 0.086958199125,
     $-0.397523196000, 0.086958199125, 0.111803398875,-0.132287565553,
     $ 0.202072594216, 0.000000000000,-0.202072594216, 0.132287565553/
*     CYLINDRICAL GAUSS-LOBATTO (VARIATIONAL COLLOCATION METHOD)
*     MATRICES:
      DATA TSQ/
     $-0.025000000000,-0.089105638513, 0.000000000000, 0.089105638513,
     $ 0.025000000000/
      DATA QSQ/
     $-3.000000000000, 3.237234826568,-0.266666666667, 0.029431840099,
     $ 0.000000000000, 3.237234826568,-4.158263130608, 0.950460144139,
     $ 0.000000000000,-0.029431840099,-0.266666666667, 0.950460144139,
     $ 0.000000000000,-0.950460144139, 0.266666666667, 0.029431840099,
     $ 0.000000000000,-0.950460144139, 4.158263130608,-3.237234826568,
     $ 0.000000000000,-0.029431840099, 0.266666666667,-3.237234826568,
     $ 3.000000000000/
*
      LC=IELEM+1
      IF((IELEM.EQ.1).AND.(ICOL.EQ.1)) THEN
*        LL2
         HTYPE='LINEAR LAGRANGIAN POLYNOMIALS'
         DO 20 I=1,LC
         T(I)=TL(I)
         TS(I)=TSL(I)
         DO 10 J=1,LC
         R((J-1)*LC+I)=RL(I,J)
         RS((J-1)*LC+I)=RSL(I,J)
         QS((J-1)*LC+I)=QSL(I,J)
         R1DQ((J-1)*LC+I)=RLQ(I,J)
         R1DR((J-1)*LC+I)=RLR(I,J)
         E((J-1)*LC+I)=EL(I,J)
   10    CONTINUE
   20    CONTINUE
         V(1)=-1.0
         V(2)=1.0
         H(1)=0.5
         H(2)=0.5
         DO 40 I=1,6
         DO 30 J=1,6
         RH(I,J)=RHA6(I,J)
         QH(I,J)=QHA6(I,J)
   30    CONTINUE
   40    CONTINUE
         DO 60 I=1,3
         DO 50 J=1,3
         RT(I,J)=RTA(I,J)*SQRT(3.0)/24.0
         QT(I,J)=QTA(I,J)/SQRT(3.0)
   50    CONTINUE
   60    CONTINUE
      ELSE IF((IELEM.EQ.1).AND.(ICOL.EQ.2)) THEN
*        LL2$
         HTYPE='FINITE DIFFERENCES'
         DO 80 I=1,LC
         T(I)=TL(I)
         TS(I)=TSL1(I)
         DO 70 J=1,LC
         R((J-1)*LC+I)=0.0
         RS((J-1)*LC+I)=0.0
         QS((J-1)*LC+I)=QSL(I,J)
         R1DQ((J-1)*LC+I)=RL1Q(I,J)
         R1DR((J-1)*LC+I)=RL1R(I,J)
         E((J-1)*LC+I)=EL(I,J)
   70    CONTINUE
         R((I-1)*LC+I)=TL(I)
         RS((I-1)*LC+I)=TSL1(I)
   80    CONTINUE
         V(1)=-1.0
         V(2)=1.0
         H(1)=0.5
         H(2)=0.5
         DO 100 I=1,6
         DO 90 J=1,6
         RH(I,J)=RHL6(I,J)*SQRT(3.0)/4.0
         QH(I,J)=QHL6(I,J)*SQRT(3.0)
   90    CONTINUE
  100    CONTINUE
         DO 120 I=1,3
         DO 110 J=1,3
         RT(I,J)=RTL(I,J)*SQRT(3.0)/12.0
         QT(I,J)=QTL(I,J)/SQRT(3.0)
  110    CONTINUE
  120    CONTINUE
      ELSE IF((IELEM.EQ.1).AND.(ICOL.EQ.3)) THEN
*        LL2#
         HTYPE='SUPERCONVERGENT LINEAR POLYNOMIALS'
         DO 140 I=1,LC
         T(I)=TL(I)
         TS(I)=TSL2(I)
         DO 130 J=1,LC
         R((J-1)*LC+I)=RL2(I,J)
         RS((J-1)*LC+I)=RSL2(I,J)
         QS((J-1)*LC+I)=QSL(I,J)
         R1DQ((J-1)*LC+I)=RL2Q(I,J)
         R1DR((J-1)*LC+I)=RL2R(I,J)
         E((J-1)*LC+I)=EL(I,J)
  130    CONTINUE
  140    CONTINUE
         V(1)=-1.0
         V(2)=1.0
         H(1)=0.5
         H(2)=0.5
      ELSE IF((IELEM.EQ.2).AND.(ICOL.EQ.1)) THEN
*        PL3
         HTYPE='PARABOLIC LAGRANGIAN POLYNOMIALS'
         DO 170 I=1,LC
         T(I)=TP(I)
         TS(I)=TSP(I)
         DO 150 J=1,LC-1
         V((J-1)*LC+I)=VP(I,J)
         H((J-1)*LC+I)=HP(I,J)
  150    CONTINUE
         DO 160 J=1,LC
         R((J-1)*LC+I)=RP(I,J)
         RS((J-1)*LC+I)=RSP(I,J)
         QS((J-1)*LC+I)=QSP(I,J)
         E((J-1)*LC+I)=EP(I,J)
  160    CONTINUE
  170    CONTINUE
      ELSE IF((IELEM.EQ.2).AND.(ICOL.EQ.2)) THEN
*        PL3$
         HTYPE='PARABOLIC COLLOCATION METHOD'
         DO 200 I=1,LC
         T(I)=TP1(I)
         TS(I)=TSP1(I)
         DO 180 J=1,LC-1
         V((J-1)*LC+I)=VP1(I,J)
         H((J-1)*LC+I)=HP1(I,J)
  180    CONTINUE
         DO 190 J=1,LC
         R((J-1)*LC+I)=0.0
         RS((J-1)*LC+I)=0.0
         QS((J-1)*LC+I)=QSP1(I,J)
         E((J-1)*LC+I)=EP1(I,J)
  190    CONTINUE
         R((I-1)*LC+I)=TP1(I)
         RS((I-1)*LC+I)=TSP1(I)
  200    CONTINUE
      ELSE IF((IELEM.EQ.2).AND.(ICOL.EQ.3)) THEN
*        PL3#
         HTYPE='PARABOLIC SUPERCONVERGENT POLYNOMIALS'
         DO 230 I=1,LC
         T(I)=TP2(I)
         TS(I)=TSP2(I)
         DO 210 J=1,LC-1
         V((J-1)*LC+I)=VP2(I,J)
         H((J-1)*LC+I)=HP2(I,J)
  210    CONTINUE
         DO 220 J=1,LC
         R((J-1)*LC+I)=RP2(I,J)
         RS((J-1)*LC+I)=RSP2(I,J)
         QS((J-1)*LC+I)=QSP2(I,J)
         E((J-1)*LC+I)=EP2(I,J)
  220    CONTINUE
  230    CONTINUE
      ELSE IF((IELEM.EQ.3).AND.(ICOL.EQ.1)) THEN
*        CL4
         HTYPE='CUBIC LAGRANGIAN POLYNOMIALS'
         DO 260 I=1,LC
         T(I)=TC(I)
         TS(I)=TSC(I)
         DO 240 J=1,LC-1
         V((J-1)*LC+I)=VC(I,J)
         H((J-1)*LC+I)=HC(I,J)
  240    CONTINUE
         DO 250 J=1,LC
         R((J-1)*LC+I)=RC(I,J)
         RS((J-1)*LC+I)=RSC(I,J)
         QS((J-1)*LC+I)=QSC(I,J)
         E((J-1)*LC+I)=EC(I,J)
  250    CONTINUE
  260    CONTINUE
      ELSE IF((IELEM.EQ.3).AND.(ICOL.EQ.2)) THEN
*        CL4$
         HTYPE='CUBIC COLLOCATION METHOD'
         DO 290 I=1,LC
         T(I)=TC1(I)
         TS(I)=TSC1(I)
         DO 270 J=1,LC-1
         V((J-1)*LC+I)=VC1(I,J)
         H((J-1)*LC+I)=HC1(I,J)
  270    CONTINUE
         DO 280 J=1,LC
         R((J-1)*LC+I)=0.0
         RS((J-1)*LC+I)=0.0
         QS((J-1)*LC+I)=QSC1(I,J)
         E((J-1)*LC+I)=EC1(I,J)
  280    CONTINUE
         R((I-1)*LC+I)=TC1(I)
         RS((I-1)*LC+I)=TSC1(I)
  290    CONTINUE
      ELSE IF((IELEM.EQ.3).AND.(ICOL.EQ.3)) THEN
*        CL4#
         HTYPE='SUPERCONVERGENT CUBIC POLYNOMIALS'
         DO 320 I=1,LC
         T(I)=TC2(I)
         TS(I)=TSC2(I)
         DO 300 J=1,LC-1
         V((J-1)*LC+I)=VC2(I,J)
         H((J-1)*LC+I)=HC2(I,J)
  300    CONTINUE
         DO 310 J=1,LC
         R((J-1)*LC+I)=RC2(I,J)
         RS((J-1)*LC+I)=RSC2(I,J)
         QS((J-1)*LC+I)=QSC2(I,J)
         E((J-1)*LC+I)=EC2(I,J)
  310    CONTINUE
  320    CONTINUE
      ELSE IF((IELEM.EQ.4).AND.(ICOL.EQ.2)) THEN
*        QL5$
         HTYPE='QUARTIC COLLOCATION METHOD'
         DO 350 I=1,LC
         T(I)=TQ(I)
         TS(I)=TSQ(I)
         DO 330 J=1,LC-1
         V((J-1)*LC+I)=VQ(I,J)
         H((J-1)*LC+I)=HQ(I,J)
  330    CONTINUE
         DO 340 J=1,LC
         R((J-1)*LC+I)=0.0
         RS((J-1)*LC+I)=0.0
         QS((J-1)*LC+I)=QSQ(I,J)
         E((J-1)*LC+I)=EQ(I,J)
  340    CONTINUE
         R((I-1)*LC+I)=TQ(I)
         RS((I-1)*LC+I)=TSQ(I)
  350    CONTINUE
      ELSE
         CALL XABORT('BIVCOL: TYPE OF FINITE ELEMENT NOT AVAILABLE.')
      ENDIF
*----
*  COMPUTE THE CARTESIAN STIFFNESS MATRIX FROM THE TENSORIAL PRODUCT OF
*  TWO NODAL COUPLING MATRICES.
*----
      DO 380 I=1,LC
      DO 370 J=1,LC
      DSUM=0.0D0
      DO 360 K=1,LC-1
      DSUM=DSUM+V((K-1)*LC+I)*V((K-1)*LC+J)
  360 CONTINUE
      Q((J-1)*LC+I)=REAL(DSUM)
  370 CONTINUE
  380 CONTINUE
      IF(IMPX.GT.0) WRITE (6,'(/9H BIVCOL: ,A40)') HTYPE
*----
*  SAVE THE UNIT MATRICES ON LCM.
*----
      CALL LCMSIX(IPTRK,'BIVCOL',1)
      CALL LCMPUT(IPTRK,'T',LC,2,T)
      CALL LCMPUT(IPTRK,'TS',LC,2,TS)
      CALL LCMPUT(IPTRK,'R',LC*LC,2,R)
      CALL LCMPUT(IPTRK,'RS',LC*LC,2,RS)
      IF(IELEM.EQ.1) THEN
         CALL LCMPUT(IPTRK,'RSH1',LC*LC,2,R1DQ)
         CALL LCMPUT(IPTRK,'RSH2',LC*LC,2,R1DR)
      ENDIF
      CALL LCMPUT(IPTRK,'Q',LC*LC,2,Q)
      CALL LCMPUT(IPTRK,'QS',LC*LC,2,QS)
      CALL LCMPUT(IPTRK,'V',LC*(LC-1),2,V)
      CALL LCMPUT(IPTRK,'H',LC*(LC-1),2,H)
      CALL LCMPUT(IPTRK,'E',LC*LC,2,E)
      IF((IELEM.EQ.1).AND.(ICOL.LE.2)) THEN
         CALL LCMPUT(IPTRK,'RH',36,2,RH)
         CALL LCMPUT(IPTRK,'QH',36,2,QH)
         CALL LCMPUT(IPTRK,'RT',9,2,RT)
         CALL LCMPUT(IPTRK,'QT',9,2,QT)
      ENDIF
      CALL LCMSIX(IPTRK,' ',2)
      RETURN
      END