How to implement Einstein summation convention with differential operators

Aha, simpler than I thought. Assuming all I guessed in the comments is correct:

BeginPackage["einstein`"]

allowtensor;
$tensordimension = 3;

Begin["`private`"]

expand[func_, {}] := # &
expand[func_, var_] := 
 Function[s, func[s, ##] &[Sequence @@ ({#, $tensordimension} & /@ var)], HoldAll]

tensor[index_List] := 
 Function[{expr}, 
  With[{count = Count[expr // Unevaluated, #, Infinity, Heads -> True] & /@ index}, 
   expand[Table, Pick[index, OddQ[#] && # > 0 & /@ count]]@
    expand[Sum, Pick[index, EvenQ[#] && # > 0 & /@ count]]@expr], HoldAll]

flat[expr_List] := Flatten@expr
flat[expr_] := expr

SetAttributes[allowtensor, HoldFirst]
allowtensor[a_ + b_, index_List] := allowtensor[a, index] + allowtensor[b, index]
allowtensor[c_ a_, index_List] /; FreeQ[Unevaluated@c, Alternatives @@ index] := 
 c allowtensor[a, index]
allowtensor[a_ == b_, index_List] := 
 flat@allowtensor[a, index] == flat@allowtensor[b, index] // Thread
allowtensor[expr_, index_List] := tensor[index][expr]

End[]
EndPackage[]

The following is not necessary, but will make the output pretty:

rule[var_] := var[i__] :> Subscript[var, Sequence @@ x /@ {i}]

drule = Derivative[id__][f_][args__] :> 
   TraditionalForm[
    HoldForm@D[f, ##] &[
     Sequence @@ (DeleteCases[
         Transpose[{{args}, {id}}], {_, 0}] /. {x_, 1} :> x)]];

Then let's check. Some preparation:

inde = {x, y, z};

Clear@x; x[i_] := inde[[i]];

Oh, I've used x both for function definition and independent variable, which isn't a good practice, but this is just a toy example and we know what we're doing, so let it be.

Now check the first example:

allowtensor[D[σ[i, j][x, y, z], x[j]] + F[i] == 0, {i, j}] /. 
  rule /@ {σ, F} /. drule

enter image description here

The second:

ϵ = LeviCivitaTensor[3];

allowtensor[ϵ[[m, i, k]] ϵ[[n, l, j]] D[e[i, j][x, y, z], x[k], x[l]] == 0,
            {i, j, k, l, m, n}] /. 
    e[i_, j_] /; i > j -> e[j, i] /. rule[e] /. 
  drule // DeleteDuplicates

enter image description here

Notice the output is eliminated to 6 equations because of the symmetry, which should have been clarified in the body of question.


Let me try to partially answer. Partially for the following reason: I know how to implement index vector and tensor notations and how to work with them. I also wanted to implement the Einstein convention and failed. However, even without it one can successfully use the index notations.

Let us first introduce the Kronecker, \[Delta] and Levi-Civita, ee tensors:

Subscript[δ, i_, j_] := KroneckerDelta[i, j];
Subscript[ee, i_, j_, k_] := Signature[{i, j, k}];

Let us try them. This looks as

enter image description here

on your screen. I mean that on the screen it looks as we traditionally used to denote vectors and tensors in the index notations, but in the StackExchange is is clumsy. Therefore, in the following I include the screenshots.

Subscript[ee, 1, 2, 3]
Subscript[ee, 1, 1, 3]

(* 1

0 *)

This is the contraction of the Levi-Civita with the Kronecker tensor

enter image description here

    Sum[Subscript[ee, i, j, k]*Subscript[δ, i, k], {i, 1, 3}, {k, 
       1, 3}] /. j -> 3
(* 0 *)

This is the example of a vector product:

enter image description here

Subscript[s, i_] := 
  Sum[Subscript[ee, i, j, k] Subscript[a, j] Subscript[b, k], {j, 1, 
    3}, {k, 1, 3}];
Subscript[s, 1]

-Subscript[a, 3] Subscript[b, 2] + Subscript[a, 2] Subscript[b, 3]

Here is an example of an electrodynamics calculation of the magnetic field enter image description here as a part of a Fresnel problem within this technique

enter image description here

This is a simple example from the elasticity theory (since you seem to be interested precisely in this area):

enter image description here

Subscript[ϵ, 1, 1] = 
  1/Ε*(Subscript[σ, 1, 
     1] - ν*(Subscript[σ, 2, 2] + Subscript[σ, 3, 
        3]));
Subscript[ϵ, 2, 2] = 
  1/Ε*(Subscript[σ, 2, 
     2] - ν*(Subscript[σ, 1, 1] + Subscript[σ, 3, 
        3]));
Subscript[ϵ, 3, 3] = 
  1/Ε*(Subscript[σ, 3, 
     3] - ν*(Subscript[σ, 1, 1] + Subscript[σ, 2, 
        2]));
expr = (Sum[Subscript[ϵ, i, i], {i, 1, 3}]) /. 
   Subscript[σ, 3, 
    3] -> ν*(Subscript[σ, 1, 1] + Subscript[σ, 2, 
       2]) // Factor

(* -(((1 + ν) (-1 + 2 ν) (Subscript[σ, 1, 1] + 
    Subscript[σ, 2, 2]))/Ε)  *)

I have more examples from elasticity theory including operating with derivatives and Green functions. However, I feel that this answer is already too long.

Have fun!