Core Functions

dC[oC_,nC_,nCL_,i_]:=Module[{d,xs,x},
d=Length[nCL];
xs=Table[x[k],{k,d}];Function[f,Sum[D[f,nCL[[j]]]D[CoordinateTransformData[oC->nC,"Mapping",xs][[j]],x[i]],{j,d}]/.Thread[Rule[xs,CoordinateTransformData[nC->oC,"Mapping",nCL]]]]
]

derivativeC[oC_,nC_,nCL_][n__][f_][x__]:=Composition[ReplacePart[MapThread[(ff|->Nest[#1,ff,#2])&,{Table[dC[oC,nC,nCL,i],{i,Length[nCL]}],{n}}],0->Sequence]][f[nCL/.List->Sequence]]

Core Function Examples

$Assumptions={r>0,0<θ<π,0<ϕ<2π};

dC["Cartesian","Spherical",{r,θ,ϕ},1][f[r,θ,ϕ]]//Simplify//Expand

D[f[1][x,y,z],x]

%/.{Derivative->derivativeC["Cartesian","Spherical",{r,θ,ϕ}]}//Simplify//Expand

%/.{x__[f[1]][r,θ,ϕ]->x[f[#1]Cos[#2]Sin[#3]&][r,θ,ϕ],x__[f[2]][r,θ,ϕ]->x[f[#1]Sin[#2]Sin[#3]&][r,θ,ϕ],x__[f[3]][r,θ,ϕ]->x[f[#1]Cos[#3]&][r,θ,ϕ]}

The Elastic Equations

xt=Table[x[i],{i,3}];xs=(xt/.List->Sequence);

diuj=D[Table[u[j][xs],{j,3}],{xt}];

diuldjul=Table[Sum[D[u[l][xs],x[i]]D[u[l][xs],x[j]],{l,3}],{i,3},{j,3}];

ϵ=(1/2)(diuj+Transpose[diuj]+diuldjul);

U=Sum[λ ϵ[[i,i]]ϵ[[j,j]]+2μ ϵ[[i,j]]ϵ[[i,j]],{i,3},{j,3}];

dUddu=Table[D[U,{D[u[i][xs],{xt}]}],{i,3}];

cartesianEqs=Map[uu|->(DiracDelta[x[1],x[2],x[3]]==Total[MapIndexed[D[#1,x[#2[[1]]]]&,uu]]),dUddu]//Simplify;

sphericalEqs=cartesianEqs/.{Derivative->derivativeC["Cartesian","Spherical",{r,θ,ϕ}]};

radialEqs=sphericalEqs/.{x__[u[1]][r,θ,ϕ]->x[u[#1]Cos[#3]Sin[#2]&][r,θ,ϕ],x__[u[2]][r,θ,ϕ]->x[u[#1]Sin[#3]Sin[#2]&][r,θ,ϕ],x__[u[3]][r,θ,ϕ]->x[u[#1]Cos[#2]&][r,θ,ϕ]}//Simplify;

radialEqs=radialEqs/.DiracDelta[x[1],x[2],x[3]]->DiracDelta[r]/(2 π (r^2))//Simplify

numRadialEqs=radialeqs/.{x__[u]->x[U[#]/((#)^2)&],u[r]->U[r]/(r^2)};

polarBoundaries={U[1]==0,U'[1]==1};

sol=NDSolve[Join[{(numRadialEqs[[3]]/.{λ->(1/4),μ->(1/2),θ->π})},polarBoundaries],U,{r,(10^(-8)),1}];

DensityPlot[Log[Evaluate[-U[r]/(r^2)/.sol]]/.r->(Sqrt[(x^2)+(y^2)]),{x,-1/(Sqrt[2]),1/(Sqrt[2])},{y,-4/(6(Sqrt[2])),4/(6(Sqrt[2]))},AspectRatio->4/6,Frame->None,ColorFunction->"SouthwestColors",PlotPoints->100]