Hyperelastic laws for 3d

Hello,

I have attempted to rewrite the hyperelastic macros given ElasticLaw2d.idp file into 3d version. But thus far I cant get any results from it in both static as dynamic problems. Can anyone give me some pointers if I made some mistakes in this macro file ? or is this approach not viable in 3d?

macro Strain3(d) 
[ 
dx(d[0]), 
dy(d[1]),
dz(d[2]), 
(dz(d[1])+dy(d[2]))/2.,
(dy(d[0])+dx(d[1]))/2.,
(dx(d[2])+dz(d[0]))/2.
  
] //EOM 

// --------------------------------------------------------------------------
// definition of C(d)  = $ {}^t F F$   and 2 differentiel
// with  , with $F = Id + \nabla d$ where $Id$ is the identity matrix 
macro C2(d)
[
1. + 2.*dx(d[0]) + dx(d[0])*dx(d[0]) + dx(d[1])*dx(d[1]) + dx(d[2])*dx(d[2]),
1. + 2.*dy(d[1]) + dy(d[0])*dy(d[0]) + dy(d[1])*dy(d[1]) + dy(d[2])*dy(d[2]),
1. + 2.*dz(d[2]) + dz(d[0])*dz(d[0]) + dz(d[1])*dz(d[1]) + dz(d[2])*dz(d[2]),
dz(d[1])+dy(d[2]) + dz(d[0])*dy(d[0]) + dz(d[1])*dy(d[1])+ dz(d[2])*dy(d[2]),
dy(d[0])+dx(d[1]) + dy(d[0])*dx(d[0]) + dy(d[1])*dx(d[1])+ dy(d[2])*dx(d[2]),
dz(d[0])+dx(d[2]) + dz(d[0])*dx(d[0]) + dz(d[1])*dx(d[1])+ dz(d[2])*dx(d[2])
 
] //
macro dC2(d,dd)
[
 2.*dx((dd)[0]) + 2.*dx((dd)[0])*dx(d[0]) + 2.*dx((dd)[1])*dx(d[1])+ 2.*dx((dd)[2])*dx(d[2]),
 2.*dy((dd)[1]) + 2.*dy((dd)[0])*dy(d[0]) + 2.*dy((dd)[1])*dy(d[1]) + 2.*dy((dd)[2])*dy(d[2]),
 2.*dz((dd)[2]) + 2.*dz((dd)[0])*dz(d[0]) + 2.*dz((dd)[1])*dz(d[1]) + 2.*dz((dd)[2])*dz(d[2]),
 dz((dd)[1]) + dy((dd)[2])  + dz((dd)[0])*dy(d[0]) + dz((dd)[1])*dy(d[1]) + dz((dd)[2])*dy(d[2]) 
                        + dz(d[0])*dy((dd)[0]) + dz(d[1])*dy((dd)[1]) + dz(d[2])*dy((dd)[2]),
 dy((dd)[0]) + dx((dd)[1])  + dy((dd)[0])*dx(d[0]) + dy((dd)[1])*dx(d[1]) + dy((dd)[2])*dx(d[2]) 
                        + dy(d[0])*dx((dd)[0]) + dy(d[1])*dx((dd)[1]) + dy(d[2])*dx((dd)[2]),						
 dz((dd)[0]) + dx((dd)[2])  + dz((dd)[0])*dx(d[0]) + dz((dd)[1])*dx(d[1]) + dz((dd)[2])*dx(d[2]) 
                        + dz(d[0])*dx((dd)[0]) + dz(d[1])*dx((dd)[1]) + dz(d[2])*dx((dd)[2]) 
] //
macro ddC2(dd,ddd)
[
 2.*dx((dd)[0])*dx(ddd[0]) + 2.*dx((dd)[1])*dx(ddd[1])+ 2.*dx((dd)[2])*dx(ddd[2]),
 2.*dy((dd)[0])*dy(ddd[0]) + 2.*dy((dd)[1])*dy(ddd[1]) + 2.*dy((dd)[2])*dy(ddd[2]),
 2.*dz((dd)[0])*dz(ddd[0]) + 2.*dz((dd)[1])*dz(ddd[1]) + 2.*dz((dd)[2])*dz(ddd[2]),
 dz((dd)[0])*dy(ddd[0]) + dz((dd)[1])*dy(ddd[1]) + dz((dd)[2])*dy(ddd[2]) 
                        + dz(ddd[0])*dy((dd)[0]) + dz(ddd[1])*dy((dd)[1]) + dz(ddd[2])*dy((dd)[2]),
dy((dd)[0])*dx(ddd[0]) + dy((dd)[1])*dx(ddd[1]) + dy((dd)[2])*dx(ddd[2]) 
                        + dy(ddd[0])*dx((dd)[0]) + dy(ddd[1])*dx((dd)[1]) + dy(ddd[2])*dx((dd)[2]),						
dz((dd)[0])*dx(ddd[0]) + dz((dd)[1])*dx(ddd[1]) + dz((dd)[2])*dx(ddd[2]) 
                        + dz(ddd[0])*dx((dd)[0]) + dz(ddd[1])*dx((dd)[1]) + dz(ddd[2])*dx((dd)[2]) 
] //

// --------------------------------------------------------------------------
// definition of the 3 invariant I[0], I[1], I[2]in 3d 
// so the value is an array size 3 . 
// I2C (I (C) 2d ..) 
// 
//  I[0] = trace(C) 
//  I[1] = (trace(C)^2 - trace(C^2) )/2
//  I[2] = det(C) 
//   ************** BIZARRE FH ***********************

macro I2C(C) 
[
  C[0] + C[1] + C[2] , 
  C[0]*C[1] + C[1]*C[2] + C[0]*C[2] - C[4]*C[4] - C[3]*C[3] - C[5]*C[5],
  C[0]*C[1]*C[2] +C[4]*C[3]*C[5] +C[5]*C[4]*C[3] - C[5]*C[1]*C[5]- C[4]*C[4]*C[2]-C[0]*C[3]*C[3]
] //EOM

macro dI2C(C,dC) 
 [ 
   dC[0] + dC[1] + dC[2] , 
   dC[0]*C[1]+C[0]*dC[1]+dC[1]*C[2]+C[1]*dC[2]+ dC[0]*C[2]+C[0]*dC[2]-2*C[4]*dC[4] -2*C[3]*dC[3]-2*C[5]*dC[5], 
   dC[0]*C[1]*C[2]+C[0]*dC[1]*C[2]+C[0]*C[1]*dC[2]
   +dC[4]*C[3]*C[5]+C[4]*dC[3]*C[5]+C[4]*C[3]*dC[5]
   +dC[5]*C[4]*C[3]+C[5]*dC[4]*C[3]+C[5]*C[4]*dC[3]
   -2*dC[5]*C[1]*C[5]-C[5]*dC[1]*C[5]
   -2*dC[4]*C[4]*C[2]- C[4]*C[4]*dC[2]
   -dC[0]*C[3]*C[3]-2*C[0]*dC[3]*C[3]
 ] // 

macro  ddI2C(C,dC,ddC) 
[ 
   0. , 
  dC[0]*ddC[1]+ddC[0]*dC[1]+dC[1]*ddC[2]+ddC[1]*dC[2]+ dC[0]*ddC[2]+ddC[0]*dC[2]-2*ddC[4]*dC[4] -2*ddC[3]*dC[3]-2*ddC[5]*dC[5],
    dC[0]*ddC[1]*C[2]+dC[0]*C[1]*ddC[2]
    +ddC[0]*dC[1]*C[2]+C[0]*dC[1]*ddC[2]
	+ddC[0]*C[1]*dC[2]+C[0]*ddC[1]*dC[2]
	+dC[4]*ddC[3]*C[5]+dC[4]*C[3]*ddC[5]
	+ddC[5]*dC[4]*C[3]+C[5]*dC[4]*ddC[3]
	+ddC[4]*C[3]*dC[5]+C[4]*ddC[3]*dC[5]
	+dC[5]*ddC[4]*C[3]+dC[5]*C[4]*ddC[3]
	+ddC[5]*dC[4]*C[3]+C[5]*dC[4]*ddC[3]
	+ddC[5]*C[4]*dC[3]+C[5]*ddC[4]*dC[3]
	-2*dC[5]*ddC[1]*C[5]-2*dC[5]*C[1]*ddC[5]
	-2*ddC[5]*dC[1]*C[5]
	 -2*dC[4]*ddC[4]*C[2] -2*dC[4]*C[4]*ddC[2]
	- 2*ddC[4]*C[4]*dC[2]
	-2*dC[0]*ddC[3]*C[3]
	-2*ddC[0]*dC[3]*C[3]-2*C[0]*dC[3]*ddC[3]
 ] // 
 

macro I2d(d) I2C(C2(d))  //
macro dI2d(d,dd)  dI2C( C2(d) , dC2(d,(dd)) ) //
macro ddI2d(d,dd,ddd)  ( ddI2C(C2(d), dC2(d,(dd)),dC2(d,(ddd))  ) + dI2C( C2(d) , ddC2((dd),(ddd)) ) ) //
macro W2d(d) W(I2d(d)) //
macro dW2d(d,dd)  dW( I2d(d) , dI2d(d,(dd))  ) //
macro ddW2d(d,dd,ddd)  
(  ddW( I2d(d) , dI2d(d,(dd)) ,dI2d(d,(ddd)) )  + dW( I2d(d) , ddI2d(d,(dd),(ddd)) ) ) //