Hyperelastic flux definitions  for 3D

Definition of field variables  q=(v,F), and their flux  f=(P,v⊗δ) for hyperelastic media described by conservation laws
(1)	(Subscript[(ρ v), t])=div P(F)
(Subscript[F, t])=div (v⊗δ), 
rewritten as
(2)	(Subscript[q, t])+div f(q)=0,
with ρ - mass density, v - velocity, P first Piola-Kirchhoff stress, F - deformation gradient, v⊗δ - exterior product of velocity with identity matrix.  The linearized version of (2) is
(3)	(Subscript[q, t])+A·div(q)=0. 
with the flux Jacobian A=∂f/∂q.

Notation

Tensors are generally denoted by double strike Latin letters such as: 𝕩, 𝔽 , with keyboard shortcuts EscdsxEsc, EscdsFEsc. The scalar components are denoted by the indexed, corresponding normal font Latin letter, 𝕩={(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}. Tensors with standard notation by Greek letters are denoted by the corresponding script letter ℯ, with keyboard shortcut EscdseEsc, and scalar components denoted by the indexed Greek letters, ℯ={{(Subscript[ϵ, 11]),(Subscript[ϵ, 12]),(Subscript[ϵ, 13])},{(Subscript[ϵ, 21]),(Subscript[ϵ, 22]),(Subscript[ϵ, 23])},{(Subscript[ϵ, 31]),(Subscript[ϵ, 32]),(Subscript[ϵ, 33])}}. When the scalar components of the tensor are not used, a tensor can also be denoted by a normal face Greek letter.

Define number of space dimensions

d=3

Define identity matrix, Dirac delta

𝕀=IdentityMatrix[d]; δ=𝕀; MatrixForm[δ]

Define some utility functions for Einstein summation notation

idx2[i_,j_]:=ToString[i]<>ToString[j];
idx3[i_,j_,k_]:=ToString[i]<>ToString[j]<>ToString[k];
idx4[i_,j_,k_,l_]:=ToString[i]<>ToString[j]<>ToString[k]<>ToString[l];

Flux Jacobians

Define q, f(q)

Define q=(v,F). First, the velocity vector

𝕧=Table[(Subscript[v, i]),{i,d}]

Next, the deformation gradient tensor

𝔽=Table[(Subscript[F, idx2[i,j]]),{i,d},{j,d}]

The deformation gradient tensor represented as a vector

LinearizeFTensor = Flatten[Table[
(Subscript[F, idx2[i,j]])->(Subscript[F, Position[Flatten[𝔽],(Subscript[F, idx2[i,j]])][[1,1]]]),{i,d},{j,d}
]]

ℒ𝔽 =Flatten[𝔽 /. LinearizeFTensor]

𝕢=Flatten[{𝕧,ℒ𝔽}]

Define f(q)=(P(F),v⊗δ). First the first Piola-Kirchhoff (PK1) tensor

ℙ=Table[(Subscript[P, idx2[i,j]]),{i,d},{j,d}]; TableForm[ℙ]

The vector representation of the PK1 tensor

LinearizePTensor = Flatten[Table[
(Subscript[P, idx2[i,j]])->(Subscript[P, Position[Flatten[ℙ],(Subscript[P, idx2[i,j]])][[1,1]]]),{i,d},{j,d}
]]

ℒℙ =Flatten[ℙ /. LinearizePTensor]

The flux of the deformation gradient 𝕧⊗δ, a rank-3 tensor, represented as a rank-2 array 𝕎=𝔓(𝕧⊗δ).

𝕎=Flatten[Outer[Times,𝕧,δ],1]; TableForm[𝕎]

The overall flux 𝕗(𝕢) for quantities within 𝕢. This is an array of dimensions (d+(d^2))×d, i.e., 6×4 for d=2.

𝕗=Join[ℙ /. LinearizePTensor,𝕎,1]; TableForm[𝕗]

𝕢

Define flux Jacobian ∂f/∂q

The flux Jacobian, using the notation (Subscript[P, l,m])=∂(Subscript[P, l])/∂(Subscript[F, m]), with l,m indices over the vector representations of the ℙ, 𝔽 tensors.

dPdvRules = Flatten[Table[dq[(Subscript[P, l]),(Subscript[v, j])]->0,{l,(d^2)},{j,d}]]

dPdFRules = Flatten[Table[dq[(Subscript[P, l]),(Subscript[F, m])]->(Subscript[P, l,m]),{l,(d^2)},{m,(d^2)}]]

dvdvRules = Flatten[Table[dq[(Subscript[v, i]),(Subscript[v, j])]->δ[[i,j]],{i,d},{j,d}]]

zerodvRules = Flatten[Table[dq[0,(Subscript[v, j])]->0,{j,d}]]

zerodFRules = Flatten[Table[dq[0,(Subscript[F, l])]->0,{l,(d^2)}]]

dvdFRules = Flatten[Table[dq[(Subscript[v, i]),(Subscript[F, m])]->0,{i,d},{m,(d^2)}]]

dqdfRules = Flatten[Table[dq[(Subscript[P, l]),(Subscript[F, m])]->(Subscript[P, l,m]),{l,(d^2)},{m,(d^2)}]]

A[k_]:=Table[dq[𝕗[[α,k]],𝕢[[β]]],{α,d+(d^2)},{β,d+(d^2)}] /. dPdFRules /. dPdvRules /. dvdFRules /. dvdvRules /. zerodFRules /. zerodvRules

Display flux Jacobians

A[1] // TableForm

A[2] // TableForm

A[3] // TableForm

Isolate blocks within flux Jacobian

H = Map[#[[1;;d,d+1;;d+(d^2)]]&,{A[1],A[2],A[3]}];
Map[TableForm[#]&,H]

U = Map[#[[d+1;;d+(d^2),1;;d]]&,{A[1],A[2],A[3]}];
Map[TableForm[#]&,U]