code20.txt

! matrix2 maps from the 16 dof to the 20 pseudo-dof 
      double precision matrix2(20,16)
      matrix2(1,1)=1
      matrix2(1,2)=0
      matrix2(1,3)=0
      matrix2(1,4)=0
      matrix2(1,5)=0
      matrix2(1,6)=0
      matrix2(1,7)=0
      matrix2(1,8)=0
      matrix2(1,9)=0
      matrix2(1,10)=0
      matrix2(1,11)=0
      matrix2(1,12)=0
      matrix2(1,13)=0
      matrix2(1,14)=0
      matrix2(1,15)=0
      matrix2(1,16)=0
      matrix2(2,1)=0
      matrix2(2,2)=1
      matrix2(2,3)=0
      matrix2(2,4)=0
      matrix2(2,5)=0
      matrix2(2,6)=0
      matrix2(2,7)=0
      matrix2(2,8)=0
      matrix2(2,9)=0
      matrix2(2,10)=0
      matrix2(2,11)=0
      matrix2(2,12)=0
      matrix2(2,13)=0
      matrix2(2,14)=0
      matrix2(2,15)=0
      matrix2(2,16)=0
      matrix2(3,1)=0
      matrix2(3,2)=0
      matrix2(3,3)=1
      matrix2(3,4)=0
      matrix2(3,5)=0
      matrix2(3,6)=0
      matrix2(3,7)=0
      matrix2(3,8)=0
      matrix2(3,9)=0
      matrix2(3,10)=0
      matrix2(3,11)=0
      matrix2(3,12)=0
      matrix2(3,13)=0
      matrix2(3,14)=0
      matrix2(3,15)=0
      matrix2(3,16)=0
      matrix2(4,1)=0
      matrix2(4,2)=0
      matrix2(4,3)=0
      matrix2(4,4)=1
      matrix2(4,5)=0
      matrix2(4,6)=0
      matrix2(4,7)=0
      matrix2(4,8)=0
      matrix2(4,9)=0
      matrix2(4,10)=0
      matrix2(4,11)=0
      matrix2(4,12)=0
      matrix2(4,13)=0
      matrix2(4,14)=0
      matrix2(4,15)=0
      matrix2(4,16)=0
      matrix2(5,1)=0
      matrix2(5,2)=0
      matrix2(5,3)=0
      matrix2(5,4)=0
      matrix2(5,5)=1
      matrix2(5,6)=0
      matrix2(5,7)=0
      matrix2(5,8)=0
      matrix2(5,9)=0
      matrix2(5,10)=0
      matrix2(5,11)=0
      matrix2(5,12)=0
      matrix2(5,13)=0
      matrix2(5,14)=0
      matrix2(5,15)=0
      matrix2(5,16)=0
      matrix2(6,1)=0
      matrix2(6,2)=0
      matrix2(6,3)=0
      matrix2(6,4)=0
      matrix2(6,5)=0
      matrix2(6,6)=1
      matrix2(6,7)=0
      matrix2(6,8)=0
      matrix2(6,9)=0
      matrix2(6,10)=0
      matrix2(6,11)=0
      matrix2(6,12)=0
      matrix2(6,13)=0
      matrix2(6,14)=0
      matrix2(6,15)=0
      matrix2(6,16)=0
      matrix2(7,1)=0
      matrix2(7,2)=0
      matrix2(7,3)=0
      matrix2(7,4)=0
      matrix2(7,5)=0
      matrix2(7,6)=0
      matrix2(7,7)=1
      matrix2(7,8)=0
      matrix2(7,9)=0
      matrix2(7,10)=0
      matrix2(7,11)=0
      matrix2(7,12)=0
      matrix2(7,13)=0
      matrix2(7,14)=0
      matrix2(7,15)=0
      matrix2(7,16)=0
      matrix2(8,1)=0
      matrix2(8,2)=0
      matrix2(8,3)=0
      matrix2(8,4)=0
      matrix2(8,5)=0
      matrix2(8,6)=0
      matrix2(8,7)=0
      matrix2(8,8)=1
      matrix2(8,9)=0
      matrix2(8,10)=0
      matrix2(8,11)=0
      matrix2(8,12)=0
      matrix2(8,13)=0
      matrix2(8,14)=0
      matrix2(8,15)=0
      matrix2(8,16)=0
      matrix2(9,1)=0
      matrix2(9,2)=0
      matrix2(9,3)=0
      matrix2(9,4)=0
      matrix2(9,5)=0
      matrix2(9,6)=0
      matrix2(9,7)=0
      matrix2(9,8)=0
      matrix2(9,9)=1
      matrix2(9,10)=0
      matrix2(9,11)=0
      matrix2(9,12)=0
      matrix2(9,13)=0
      matrix2(9,14)=0
      matrix2(9,15)=0
      matrix2(9,16)=0
      matrix2(10,1)=0
      matrix2(10,2)=0
      matrix2(10,3)=0
      matrix2(10,4)=0
      matrix2(10,5)=0
      matrix2(10,6)=0
      matrix2(10,7)=0
      matrix2(10,8)=0
      matrix2(10,9)=0
      matrix2(10,10)=1
      matrix2(10,11)=0
      matrix2(10,12)=0
      matrix2(10,13)=0
      matrix2(10,14)=0
      matrix2(10,15)=0
      matrix2(10,16)=0
      matrix2(11,1)=0
      matrix2(11,2)=0
      matrix2(11,3)=0
      matrix2(11,4)=0
      matrix2(11,5)=0
      matrix2(11,6)=0
      matrix2(11,7)=0
      matrix2(11,8)=0
      matrix2(11,9)=0
      matrix2(11,10)=0
      matrix2(11,11)=1
      matrix2(11,12)=0
      matrix2(11,13)=0
      matrix2(11,14)=0
      matrix2(11,15)=0
      matrix2(11,16)=0
      matrix2(12,1)=0
      matrix2(12,2)=0
      matrix2(12,3)=0
      matrix2(12,4)=0
      matrix2(12,5)=0
      matrix2(12,6)=0
      matrix2(12,7)=0
      matrix2(12,8)=0
      matrix2(12,9)=0
      matrix2(12,10)=0
      matrix2(12,11)=0
      matrix2(12,12)=1
      matrix2(12,13)=0
      matrix2(12,14)=0
      matrix2(12,15)=0
      matrix2(12,16)=0
      matrix2(13,1)=0
      matrix2(13,2)=0
      matrix2(13,3)=0
      matrix2(13,4)=0
      matrix2(13,5)=0
      matrix2(13,6)=0
      matrix2(13,7)=0
      matrix2(13,8)=0
      matrix2(13,9)=0
      matrix2(13,10)=0
      matrix2(13,11)=0
      matrix2(13,12)=0
      matrix2(13,13)=1
      matrix2(13,14)=0
      matrix2(13,15)=0
      matrix2(13,16)=0
      matrix2(14,1)=0
      matrix2(14,2)=0
      matrix2(14,3)=0
      matrix2(14,4)=0
      matrix2(14,5)=0
      matrix2(14,6)=0
      matrix2(14,7)=0
      matrix2(14,8)=0
      matrix2(14,9)=0
      matrix2(14,10)=0
      matrix2(14,11)=0
      matrix2(14,12)=0
      matrix2(14,13)=0
      matrix2(14,14)=1
      matrix2(14,15)=0
      matrix2(14,16)=0
      matrix2(15,1)=0
      matrix2(15,2)=0
      matrix2(15,3)=0
      matrix2(15,4)=0
      matrix2(15,5)=0
      matrix2(15,6)=0
      matrix2(15,7)=0
      matrix2(15,8)=0
      matrix2(15,9)=0
      matrix2(15,10)=0
      matrix2(15,11)=0
      matrix2(15,12)=0
      matrix2(15,13)=0
      matrix2(15,14)=0
      matrix2(15,15)=1
      matrix2(15,16)=0
      matrix2(16,1)=0
      matrix2(16,2)=0
      matrix2(16,3)=0
      matrix2(16,4)=0
      matrix2(16,5)=0
      matrix2(16,6)=0
      matrix2(16,7)=0
      matrix2(16,8)=0
      matrix2(16,9)=0
      matrix2(16,10)=0
      matrix2(16,11)=0
      matrix2(16,12)=0
      matrix2(16,13)=0
      matrix2(16,14)=0
      matrix2(16,15)=0
      matrix2(16,16)=1
      matrix2(17,1)=0.3333333333333333
      matrix2(17,2)=(-2*x1 + x2 + x4)/18.
      matrix2(17,3)=(-2*y1 + y2 + y4)/18.
      matrix2(17,4)=(-2*z1 + z2 + z4)/18.
      matrix2(17,5)=0.3333333333333333
      matrix2(17,6)=(x1 - 2*x2 + x4)/18.
      matrix2(17,7)=(y1 - 2*y2 + y4)/18.
      matrix2(17,8)=(z1 - 2*z2 + z4)/18.
      matrix2(17,9)=0
      matrix2(17,10)=0
      matrix2(17,11)=0
      matrix2(17,12)=0
      matrix2(17,13)=0.3333333333333333
      matrix2(17,14)=(x1 + x2 - 2*x4)/18.
      matrix2(17,15)=(y1 + y2 - 2*y4)/18.
      matrix2(17,16)=(z1 + z2 - 2*z4)/18.
      matrix2(18,1)=0.3333333333333333
      matrix2(18,2)=(-2*x1 + x2 + x3)/18.
      matrix2(18,3)=(-2*y1 + y2 + y3)/18.
      matrix2(18,4)=(-2*z1 + z2 + z3)/18.
      matrix2(18,5)=0.3333333333333333
      matrix2(18,6)=(x1 - 2*x2 + x3)/18.
      matrix2(18,7)=(y1 - 2*y2 + y3)/18.
      matrix2(18,8)=(z1 - 2*z2 + z3)/18.
      matrix2(18,9)=0.3333333333333333
      matrix2(18,10)=(x1 + x2 - 2*x3)/18.
      matrix2(18,11)=(y1 + y2 - 2*y3)/18.
      matrix2(18,12)=(z1 + z2 - 2*z3)/18.
      matrix2(18,13)=0
      matrix2(18,14)=0
      matrix2(18,15)=0
      matrix2(18,16)=0
      matrix2(19,1)=0.3333333333333333
      matrix2(19,2)=(-2*x1 + x3 + x4)/18.
      matrix2(19,3)=(-2*y1 + y3 + y4)/18.
      matrix2(19,4)=(-2*z1 + z3 + z4)/18.
      matrix2(19,5)=0
      matrix2(19,6)=0
      matrix2(19,7)=0
      matrix2(19,8)=0
      matrix2(19,9)=0.3333333333333333
      matrix2(19,10)=(x1 - 2*x3 + x4)/18.
      matrix2(19,11)=(y1 - 2*y3 + y4)/18.
      matrix2(19,12)=(z1 - 2*z3 + z4)/18.
      matrix2(19,13)=0.3333333333333333
      matrix2(19,14)=(x1 + x3 - 2*x4)/18.
      matrix2(19,15)=(y1 + y3 - 2*y4)/18.
      matrix2(19,16)=(z1 + z3 - 2*z4)/18.
      matrix2(20,1)=0
      matrix2(20,2)=0
      matrix2(20,3)=0
      matrix2(20,4)=0
      matrix2(20,5)=0.3333333333333333
      matrix2(20,6)=(-2*x2 + x3 + x4)/18.
      matrix2(20,7)=(-2*y2 + y3 + y4)/18.
      matrix2(20,8)=(-2*z2 + z3 + z4)/18.
      matrix2(20,9)=0.3333333333333333
      matrix2(20,10)=(x2 - 2*x3 + x4)/18.
      matrix2(20,11)=(y2 - 2*y3 + y4)/18.
      matrix2(20,12)=(z2 - 2*z3 + z4)/18.
      matrix2(20,13)=0.3333333333333333
      matrix2(20,14)=(x2 + x3 - 2*x4)/18.
      matrix2(20,15)=(y2 + y3 - 2*y4)/18.
      matrix2(20,16)=(z2 + z3 - 2*z4)/18.

! the q-matrix maps from the 20 polynomial coeffs. to the 20 pseudo-dof
      double precision q(20,20)
      q(1,1)=1
      q(1,2)=z1
      q(1,3)=z1**2
      q(1,4)=z1**3
      q(1,5)=y1
      q(1,6)=y1*z1
      q(1,7)=y1*z1**2
      q(1,8)=y1**2
      q(1,9)=y1**2*z1
      q(1,10)=y1**3
      q(1,11)=x1
      q(1,12)=x1*z1
      q(1,13)=x1*z1**2
      q(1,14)=x1*y1
      q(1,15)=x1*y1*z1
      q(1,16)=x1*y1**2
      q(1,17)=x1**2
      q(1,18)=x1**2*z1
      q(1,19)=x1**2*y1
      q(1,20)=x1**3
      q(2,1)=0
      q(2,2)=0
      q(2,3)=0
      q(2,4)=0
      q(2,5)=0
      q(2,6)=0
      q(2,7)=0
      q(2,8)=0
      q(2,9)=0
      q(2,10)=0
      q(2,11)=1
      q(2,12)=z1
      q(2,13)=z1**2
      q(2,14)=y1
      q(2,15)=y1*z1
      q(2,16)=y1**2
      q(2,17)=2*x1
      q(2,18)=2*x1*z1
      q(2,19)=2*x1*y1
      q(2,20)=3*x1**2
      q(3,1)=0
      q(3,2)=0
      q(3,3)=0
      q(3,4)=0
      q(3,5)=1
      q(3,6)=z1
      q(3,7)=z1**2
      q(3,8)=2*y1
      q(3,9)=2*y1*z1
      q(3,10)=3*y1**2
      q(3,11)=0
      q(3,12)=0
      q(3,13)=0
      q(3,14)=x1
      q(3,15)=x1*z1
      q(3,16)=2*x1*y1
      q(3,17)=0
      q(3,18)=0
      q(3,19)=x1**2
      q(3,20)=0
      q(4,1)=0
      q(4,2)=1
      q(4,3)=2*z1
      q(4,4)=3*z1**2
      q(4,5)=0
      q(4,6)=y1
      q(4,7)=2*y1*z1
      q(4,8)=0
      q(4,9)=y1**2
      q(4,10)=0
      q(4,11)=0
      q(4,12)=x1
      q(4,13)=2*x1*z1
      q(4,14)=0
      q(4,15)=x1*y1
      q(4,16)=0
      q(4,17)=0
      q(4,18)=x1**2
      q(4,19)=0
      q(4,20)=0
      q(5,1)=1
      q(5,2)=z2
      q(5,3)=z2**2
      q(5,4)=z2**3
      q(5,5)=y2
      q(5,6)=y2*z2
      q(5,7)=y2*z2**2
      q(5,8)=y2**2
      q(5,9)=y2**2*z2
      q(5,10)=y2**3
      q(5,11)=x2
      q(5,12)=x2*z2
      q(5,13)=x2*z2**2
      q(5,14)=x2*y2
      q(5,15)=x2*y2*z2
      q(5,16)=x2*y2**2
      q(5,17)=x2**2
      q(5,18)=x2**2*z2
      q(5,19)=x2**2*y2
      q(5,20)=x2**3
      q(6,1)=0
      q(6,2)=0
      q(6,3)=0
      q(6,4)=0
      q(6,5)=0
      q(6,6)=0
      q(6,7)=0
      q(6,8)=0
      q(6,9)=0
      q(6,10)=0
      q(6,11)=1
      q(6,12)=z2
      q(6,13)=z2**2
      q(6,14)=y2
      q(6,15)=y2*z2
      q(6,16)=y2**2
      q(6,17)=2*x2
      q(6,18)=2*x2*z2
      q(6,19)=2*x2*y2
      q(6,20)=3*x2**2
      q(7,1)=0
      q(7,2)=0
      q(7,3)=0
      q(7,4)=0
      q(7,5)=1
      q(7,6)=z2
      q(7,7)=z2**2
      q(7,8)=2*y2
      q(7,9)=2*y2*z2
      q(7,10)=3*y2**2
      q(7,11)=0
      q(7,12)=0
      q(7,13)=0
      q(7,14)=x2
      q(7,15)=x2*z2
      q(7,16)=2*x2*y2
      q(7,17)=0
      q(7,18)=0
      q(7,19)=x2**2
      q(7,20)=0
      q(8,1)=0
      q(8,2)=1
      q(8,3)=2*z2
      q(8,4)=3*z2**2
      q(8,5)=0
      q(8,6)=y2
      q(8,7)=2*y2*z2
      q(8,8)=0
      q(8,9)=y2**2
      q(8,10)=0
      q(8,11)=0
      q(8,12)=x2
      q(8,13)=2*x2*z2
      q(8,14)=0
      q(8,15)=x2*y2
      q(8,16)=0
      q(8,17)=0
      q(8,18)=x2**2
      q(8,19)=0
      q(8,20)=0
      q(9,1)=1
      q(9,2)=z3
      q(9,3)=z3**2
      q(9,4)=z3**3
      q(9,5)=y3
      q(9,6)=y3*z3
      q(9,7)=y3*z3**2
      q(9,8)=y3**2
      q(9,9)=y3**2*z3
      q(9,10)=y3**3
      q(9,11)=x3
      q(9,12)=x3*z3
      q(9,13)=x3*z3**2
      q(9,14)=x3*y3
      q(9,15)=x3*y3*z3
      q(9,16)=x3*y3**2
      q(9,17)=x3**2
      q(9,18)=x3**2*z3
      q(9,19)=x3**2*y3
      q(9,20)=x3**3
      q(10,1)=0
      q(10,2)=0
      q(10,3)=0
      q(10,4)=0
      q(10,5)=0
      q(10,6)=0
      q(10,7)=0
      q(10,8)=0
      q(10,9)=0
      q(10,10)=0
      q(10,11)=1
      q(10,12)=z3
      q(10,13)=z3**2
      q(10,14)=y3
      q(10,15)=y3*z3
      q(10,16)=y3**2
      q(10,17)=2*x3
      q(10,18)=2*x3*z3
      q(10,19)=2*x3*y3
      q(10,20)=3*x3**2
      q(11,1)=0
      q(11,2)=0
      q(11,3)=0
      q(11,4)=0
      q(11,5)=1
      q(11,6)=z3
      q(11,7)=z3**2
      q(11,8)=2*y3
      q(11,9)=2*y3*z3
      q(11,10)=3*y3**2
      q(11,11)=0
      q(11,12)=0
      q(11,13)=0
      q(11,14)=x3
      q(11,15)=x3*z3
      q(11,16)=2*x3*y3
      q(11,17)=0
      q(11,18)=0
      q(11,19)=x3**2
      q(11,20)=0
      q(12,1)=0
      q(12,2)=1
      q(12,3)=2*z3
      q(12,4)=3*z3**2
      q(12,5)=0
      q(12,6)=y3
      q(12,7)=2*y3*z3
      q(12,8)=0
      q(12,9)=y3**2
      q(12,10)=0
      q(12,11)=0
      q(12,12)=x3
      q(12,13)=2*x3*z3
      q(12,14)=0
      q(12,15)=x3*y3
      q(12,16)=0
      q(12,17)=0
      q(12,18)=x3**2
      q(12,19)=0
      q(12,20)=0
      q(13,1)=1
      q(13,2)=z4
      q(13,3)=z4**2
      q(13,4)=z4**3
      q(13,5)=y4
      q(13,6)=y4*z4
      q(13,7)=y4*z4**2
      q(13,8)=y4**2
      q(13,9)=y4**2*z4
      q(13,10)=y4**3
      q(13,11)=x4
      q(13,12)=x4*z4
      q(13,13)=x4*z4**2
      q(13,14)=x4*y4
      q(13,15)=x4*y4*z4
      q(13,16)=x4*y4**2
      q(13,17)=x4**2
      q(13,18)=x4**2*z4
      q(13,19)=x4**2*y4
      q(13,20)=x4**3
      q(14,1)=0
      q(14,2)=0
      q(14,3)=0
      q(14,4)=0
      q(14,5)=0
      q(14,6)=0
      q(14,7)=0
      q(14,8)=0
      q(14,9)=0
      q(14,10)=0
      q(14,11)=1
      q(14,12)=z4
      q(14,13)=z4**2
      q(14,14)=y4
      q(14,15)=y4*z4
      q(14,16)=y4**2
      q(14,17)=2*x4
      q(14,18)=2*x4*z4
      q(14,19)=2*x4*y4
      q(14,20)=3*x4**2
      q(15,1)=0
      q(15,2)=0
      q(15,3)=0
      q(15,4)=0
      q(15,5)=1
      q(15,6)=z4
      q(15,7)=z4**2
      q(15,8)=2*y4
      q(15,9)=2*y4*z4
      q(15,10)=3*y4**2
      q(15,11)=0
      q(15,12)=0
      q(15,13)=0
      q(15,14)=x4
      q(15,15)=x4*z4
      q(15,16)=2*x4*y4
      q(15,17)=0
      q(15,18)=0
      q(15,19)=x4**2
      q(15,20)=0
      q(16,1)=0
      q(16,2)=1
      q(16,3)=2*z4
      q(16,4)=3*z4**2
      q(16,5)=0
      q(16,6)=y4
      q(16,7)=2*y4*z4
      q(16,8)=0
      q(16,9)=y4**2
      q(16,10)=0
      q(16,11)=0
      q(16,12)=x4
      q(16,13)=2*x4*z4
      q(16,14)=0
      q(16,15)=x4*y4
      q(16,16)=0
      q(16,17)=0
      q(16,18)=x4**2
      q(16,19)=0
      q(16,20)=0
      q(17,1)=1
      q(17,2)=(z1 + z2 + z4)/3.
      q(17,3)=(z1 + z2 + z4)**2/9.
      q(17,4)=(z1 + z2 + z4)**3/27.
      q(17,5)=(y1 + y2 + y4)/3.
      q(17,6)=((y1 + y2 + y4)*(z1 + z2 + z4))/9.
      q(17,7)=((y1 + y2 + y4)*(z1 + z2 + z4)**2)/27.
      q(17,8)=(y1 + y2 + y4)**2/9.
      q(17,9)=((y1 + y2 + y4)**2*(z1 + z2 + z4))/27.
      q(17,10)=(y1 + y2 + y4)**3/27.
      q(17,11)=(x1 + x2 + x4)/3.
      q(17,12)=((x1 + x2 + x4)*(z1 + z2 + z4))/9.
      q(17,13)=((x1 + x2 + x4)*(z1 + z2 + z4)**2)/27.
      q(17,14)=((x1 + x2 + x4)*(y1 + y2 + y4))/9.
      q(17,15)=((x1 + x2 + x4)*(y1 + y2 + y4)*(z1 + z2 + z4))/27.
      q(17,16)=((x1 + x2 + x4)*(y1 + y2 + y4)**2)/27.
      q(17,17)=(x1 + x2 + x4)**2/9.
      q(17,18)=((x1 + x2 + x4)**2*(z1 + z2 + z4))/27.
      q(17,19)=((x1 + x2 + x4)**2*(y1 + y2 + y4))/27.
      q(17,20)=(x1 + x2 + x4)**3/27.
      q(18,1)=1
      q(18,2)=(z1 + z2 + z3)/3.
      q(18,3)=(z1 + z2 + z3)**2/9.
      q(18,4)=(z1 + z2 + z3)**3/27.
      q(18,5)=(y1 + y2 + y3)/3.
      q(18,6)=((y1 + y2 + y3)*(z1 + z2 + z3))/9.
      q(18,7)=((y1 + y2 + y3)*(z1 + z2 + z3)**2)/27.
      q(18,8)=(y1 + y2 + y3)**2/9.
      q(18,9)=((y1 + y2 + y3)**2*(z1 + z2 + z3))/27.
      q(18,10)=(y1 + y2 + y3)**3/27.
      q(18,11)=(x1 + x2 + x3)/3.
      q(18,12)=((x1 + x2 + x3)*(z1 + z2 + z3))/9.
      q(18,13)=((x1 + x2 + x3)*(z1 + z2 + z3)**2)/27.
      q(18,14)=((x1 + x2 + x3)*(y1 + y2 + y3))/9.
      q(18,15)=((x1 + x2 + x3)*(y1 + y2 + y3)*(z1 + z2 + z3))/27.
      q(18,16)=((x1 + x2 + x3)*(y1 + y2 + y3)**2)/27.
      q(18,17)=(x1 + x2 + x3)**2/9.
      q(18,18)=((x1 + x2 + x3)**2*(z1 + z2 + z3))/27.
      q(18,19)=((x1 + x2 + x3)**2*(y1 + y2 + y3))/27.
      q(18,20)=(x1 + x2 + x3)**3/27.
      q(19,1)=1
      q(19,2)=(z1 + z3 + z4)/3.
      q(19,3)=(z1 + z3 + z4)**2/9.
      q(19,4)=(z1 + z3 + z4)**3/27.
      q(19,5)=(y1 + y3 + y4)/3.
      q(19,6)=((y1 + y3 + y4)*(z1 + z3 + z4))/9.
      q(19,7)=((y1 + y3 + y4)*(z1 + z3 + z4)**2)/27.
      q(19,8)=(y1 + y3 + y4)**2/9.
      q(19,9)=((y1 + y3 + y4)**2*(z1 + z3 + z4))/27.
      q(19,10)=(y1 + y3 + y4)**3/27.
      q(19,11)=(x1 + x3 + x4)/3.
      q(19,12)=((x1 + x3 + x4)*(z1 + z3 + z4))/9.
      q(19,13)=((x1 + x3 + x4)*(z1 + z3 + z4)**2)/27.
      q(19,14)=((x1 + x3 + x4)*(y1 + y3 + y4))/9.
      q(19,15)=((x1 + x3 + x4)*(y1 + y3 + y4)*(z1 + z3 + z4))/27.
      q(19,16)=((x1 + x3 + x4)*(y1 + y3 + y4)**2)/27.
      q(19,17)=(x1 + x3 + x4)**2/9.
      q(19,18)=((x1 + x3 + x4)**2*(z1 + z3 + z4))/27.
      q(19,19)=((x1 + x3 + x4)**2*(y1 + y3 + y4))/27.
      q(19,20)=(x1 + x3 + x4)**3/27.
      q(20,1)=1
      q(20,2)=(z2 + z3 + z4)/3.
      q(20,3)=(z2 + z3 + z4)**2/9.
      q(20,4)=(z2 + z3 + z4)**3/27.
      q(20,5)=(y2 + y3 + y4)/3.
      q(20,6)=((y2 + y3 + y4)*(z2 + z3 + z4))/9.
      q(20,7)=((y2 + y3 + y4)*(z2 + z3 + z4)**2)/27.
      q(20,8)=(y2 + y3 + y4)**2/9.
      q(20,9)=((y2 + y3 + y4)**2*(z2 + z3 + z4))/27.
      q(20,10)=(y2 + y3 + y4)**3/27.
      q(20,11)=(x2 + x3 + x4)/3.
      q(20,12)=((x2 + x3 + x4)*(z2 + z3 + z4))/9.
      q(20,13)=((x2 + x3 + x4)*(z2 + z3 + z4)**2)/27.
      q(20,14)=((x2 + x3 + x4)*(y2 + y3 + y4))/9.
      q(20,15)=((x2 + x3 + x4)*(y2 + y3 + y4)*(z2 + z3 + z4))/27.
      q(20,16)=((x2 + x3 + x4)*(y2 + y3 + y4)**2)/27.
      q(20,17)=(x2 + x3 + x4)**2/9.
      q(20,18)=((x2 + x3 + x4)**2*(z2 + z3 + z4))/27.
      q(20,19)=((x2 + x3 + x4)**2*(y2 + y3 + y4))/27.
      q(20,20)=(x2 + x3 + x4)**3/27.

! the shape function monomials 
      double precision n(20)
      n(1)=1
      n(2)=z
      n(3)=z**2
      n(4)=z**3
      n(5)=y
      n(6)=y*z
      n(7)=y*z**2
      n(8)=y**2
      n(9)=y**2*z
      n(10)=y**3
      n(11)=x
      n(12)=x*z
      n(13)=x*z**2
      n(14)=x*y
      n(15)=x*y*z
      n(16)=x*y**2
      n(17)=x**2
      n(18)=x**2*z
      n(19)=x**2*y
      n(20)=x**3

! the shape function monomials first gradients 
      double precision gn(20,3)
      gn(1,1)=0
      gn(1,2)=0
      gn(1,3)=0
      gn(2,1)=0
      gn(2,2)=0
      gn(2,3)=1
      gn(3,1)=0
      gn(3,2)=0
      gn(3,3)=2*z
      gn(4,1)=0
      gn(4,2)=0
      gn(4,3)=3*z**2
      gn(5,1)=0
      gn(5,2)=1
      gn(5,3)=0
      gn(6,1)=0
      gn(6,2)=z
      gn(6,3)=y
      gn(7,1)=0
      gn(7,2)=z**2
      gn(7,3)=2*y*z
      gn(8,1)=0
      gn(8,2)=2*y
      gn(8,3)=0
      gn(9,1)=0
      gn(9,2)=2*y*z
      gn(9,3)=y**2
      gn(10,1)=0
      gn(10,2)=3*y**2
      gn(10,3)=0
      gn(11,1)=1
      gn(11,2)=0
      gn(11,3)=0
      gn(12,1)=z
      gn(12,2)=0
      gn(12,3)=x
      gn(13,1)=z**2
      gn(13,2)=0
      gn(13,3)=2*x*z
      gn(14,1)=y
      gn(14,2)=x
      gn(14,3)=0
      gn(15,1)=y*z
      gn(15,2)=x*z
      gn(15,3)=x*y
      gn(16,1)=y**2
      gn(16,2)=2*x*y
      gn(16,3)=0
      gn(17,1)=2*x
      gn(17,2)=0
      gn(17,3)=0
      gn(18,1)=2*x*z
      gn(18,2)=0
      gn(18,3)=x**2
      gn(19,1)=2*x*y
      gn(19,2)=x**2
      gn(19,3)=0
      gn(20,1)=3*x**2
      gn(20,2)=0
      gn(20,3)=0

! the shape function monomials second gradients 
      double precision ggn(20,3,3)
      ggn(1,1,1)=0
      ggn(1,1,2)=0
      ggn(1,1,3)=0
      ggn(1,2,1)=0
      ggn(1,2,2)=0
      ggn(1,2,3)=0
      ggn(1,3,1)=0
      ggn(1,3,2)=0
      ggn(1,3,3)=0
      ggn(2,1,1)=0
      ggn(2,1,2)=0
      ggn(2,1,3)=0
      ggn(2,2,1)=0
      ggn(2,2,2)=0
      ggn(2,2,3)=0
      ggn(2,3,1)=0
      ggn(2,3,2)=0
      ggn(2,3,3)=0
      ggn(3,1,1)=0
      ggn(3,1,2)=0
      ggn(3,1,3)=0
      ggn(3,2,1)=0
      ggn(3,2,2)=0
      ggn(3,2,3)=0
      ggn(3,3,1)=0
      ggn(3,3,2)=0
      ggn(3,3,3)=2
      ggn(4,1,1)=0
      ggn(4,1,2)=0
      ggn(4,1,3)=0
      ggn(4,2,1)=0
      ggn(4,2,2)=0
      ggn(4,2,3)=0
      ggn(4,3,1)=0
      ggn(4,3,2)=0
      ggn(4,3,3)=6*z
      ggn(5,1,1)=0
      ggn(5,1,2)=0
      ggn(5,1,3)=0
      ggn(5,2,1)=0
      ggn(5,2,2)=0
      ggn(5,2,3)=0
      ggn(5,3,1)=0
      ggn(5,3,2)=0
      ggn(5,3,3)=0
      ggn(6,1,1)=0
      ggn(6,1,2)=0
      ggn(6,1,3)=0
      ggn(6,2,1)=0
      ggn(6,2,2)=0
      ggn(6,2,3)=1
      ggn(6,3,1)=0
      ggn(6,3,2)=1
      ggn(6,3,3)=0
      ggn(7,1,1)=0
      ggn(7,1,2)=0
      ggn(7,1,3)=0
      ggn(7,2,1)=0
      ggn(7,2,2)=0
      ggn(7,2,3)=2*z
      ggn(7,3,1)=0
      ggn(7,3,2)=2*z
      ggn(7,3,3)=2*y
      ggn(8,1,1)=0
      ggn(8,1,2)=0
      ggn(8,1,3)=0
      ggn(8,2,1)=0
      ggn(8,2,2)=2
      ggn(8,2,3)=0
      ggn(8,3,1)=0
      ggn(8,3,2)=0
      ggn(8,3,3)=0
      ggn(9,1,1)=0
      ggn(9,1,2)=0
      ggn(9,1,3)=0
      ggn(9,2,1)=0
      ggn(9,2,2)=2*z
      ggn(9,2,3)=2*y
      ggn(9,3,1)=0
      ggn(9,3,2)=2*y
      ggn(9,3,3)=0
      ggn(10,1,1)=0
      ggn(10,1,2)=0
      ggn(10,1,3)=0
      ggn(10,2,1)=0
      ggn(10,2,2)=6*y
      ggn(10,2,3)=0
      ggn(10,3,1)=0
      ggn(10,3,2)=0
      ggn(10,3,3)=0
      ggn(11,1,1)=0
      ggn(11,1,2)=0
      ggn(11,1,3)=0
      ggn(11,2,1)=0
      ggn(11,2,2)=0
      ggn(11,2,3)=0
      ggn(11,3,1)=0
      ggn(11,3,2)=0
      ggn(11,3,3)=0
      ggn(12,1,1)=0
      ggn(12,1,2)=0
      ggn(12,1,3)=1
      ggn(12,2,1)=0
      ggn(12,2,2)=0
      ggn(12,2,3)=0
      ggn(12,3,1)=1
      ggn(12,3,2)=0
      ggn(12,3,3)=0
      ggn(13,1,1)=0
      ggn(13,1,2)=0
      ggn(13,1,3)=2*z
      ggn(13,2,1)=0
      ggn(13,2,2)=0
      ggn(13,2,3)=0
      ggn(13,3,1)=2*z
      ggn(13,3,2)=0
      ggn(13,3,3)=2*x
      ggn(14,1,1)=0
      ggn(14,1,2)=1
      ggn(14,1,3)=0
      ggn(14,2,1)=1
      ggn(14,2,2)=0
      ggn(14,2,3)=0
      ggn(14,3,1)=0
      ggn(14,3,2)=0
      ggn(14,3,3)=0
      ggn(15,1,1)=0
      ggn(15,1,2)=z
      ggn(15,1,3)=y
      ggn(15,2,1)=z
      ggn(15,2,2)=0
      ggn(15,2,3)=x
      ggn(15,3,1)=y
      ggn(15,3,2)=x
      ggn(15,3,3)=0
      ggn(16,1,1)=0
      ggn(16,1,2)=2*y
      ggn(16,1,3)=0
      ggn(16,2,1)=2*y
      ggn(16,2,2)=2*x
      ggn(16,2,3)=0
      ggn(16,3,1)=0
      ggn(16,3,2)=0
      ggn(16,3,3)=0
      ggn(17,1,1)=2
      ggn(17,1,2)=0
      ggn(17,1,3)=0
      ggn(17,2,1)=0
      ggn(17,2,2)=0
      ggn(17,2,3)=0
      ggn(17,3,1)=0
      ggn(17,3,2)=0
      ggn(17,3,3)=0
      ggn(18,1,1)=2*z
      ggn(18,1,2)=0
      ggn(18,1,3)=2*x
      ggn(18,2,1)=0
      ggn(18,2,2)=0
      ggn(18,2,3)=0
      ggn(18,3,1)=2*x
      ggn(18,3,2)=0
      ggn(18,3,3)=0
      ggn(19,1,1)=2*y
      ggn(19,1,2)=2*x
      ggn(19,1,3)=0
      ggn(19,2,1)=2*x
      ggn(19,2,2)=0
      ggn(19,2,3)=0
      ggn(19,3,1)=0
      ggn(19,3,2)=0
      ggn(19,3,3)=0
      ggn(20,1,1)=6*x
      ggn(20,1,2)=0
      ggn(20,1,3)=0
      ggn(20,2,1)=0
      ggn(20,2,2)=0
      ggn(20,2,3)=0
      ggn(20,3,1)=0
      ggn(20,3,2)=0
      ggn(20,3,3)=0