# Direction of NeumannValue on boundary

When solving partial differential equations NeumannValue is used to specify the flux across the boundary. For normal fluxes the detailed notes state

They appear on the boundary ∂Ω of the region Ω and specify a flux across those edges in the direction of the outward normal.

This is nice and clear but what is the positive direction if the flux is tangential to the boundary? Tangential fluxes appear in stress calculations and are known as shear forces. Using the outward normal defines a local coordinate reference plane. What is the reference plane for shear forces?

Here is how I am trying to explore this question. First we need a stress operator.

Needs["NDSolveFEM"]

ClearAll[stressOperator, stressOperatorDynamic, u, v, w, x, y, z, Y, \
ν]; stressOperator[{Y_, ν_}] := {Inactive[
Div][{{0, 0, -((Y*ν)/((1 - 2*ν)*(1 + ν)))}, {0, 0,
0}, {-Y/(2*(1 + ν)), 0, 0}}.Inactive[Grad][
w[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, -((Y*ν)/((1 - 2*ν)*(1 + ν))),
0}, {-Y/(2*(1 + ν)), 0, 0}, {0, 0, 0}}.Inactive[Grad][
v[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0,
0}, {0, -Y/(2*(1 + ν)), 0}, {0,
u[x, y, z], {x, y, z}], {x, y, z}],
Inactive[Div][{{0, 0, 0}, {0,
0, -((Y*ν)/((1 -
2*ν)*(1 + ν)))}, {0, -Y/(2*(1 + ν)),
0}}.Inactive[Grad][w[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, -Y/(2*(1 + ν)),
0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0, 0}, {0, 0,
0}}.Inactive[Grad][u[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-Y/(2*(1 + ν)), 0,
0}, {0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0}, {0,
v[x, y, z], {x, y, z}], {x, y, z}],
Inactive[Div][{{0, 0, 0}, {0,
0, -Y/(2*(1 + ν))}, {0, -((Y*ν)/((1 -
v[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, 0, -Y/(2*(1 + ν))}, {0, 0,
0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0, 0}}.Inactive[
Grad][u[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-Y/(2*(1 + ν)), 0, 0}, {0, -Y/(2*(1 + ν)), 0}, {0,
0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν)))}}.Inactive[
Grad][w[x, y, z], {x, y, z}], {x, y, z}]};


Now we make a mesh of a simple cuboid region

ClearAll[Y, ν];
Len = 1;  (*length *)
ht = 0.125; (* height *)
wd = 0.5; (* width *)
materialParameters = {Y -> 10^3, ν -> 33/100};
ss = 1;(*Shear stress on beam*)
reg = Cuboid[{0, -wd/2, 0}, {Len, wd/2, ht}];
mesh = ToElementMesh[reg]; mesh["Wireframe"]


The first target is to put a shear force on the end at x = Len with the end at x = 0 clamped. To illustrate

Show[
Region[reg],
Graphics3D[{Black,
Table[Arrow[{{Len, y, 0}, {Len, y, ht}}], {y, -wd/2, wd/2, wd/8}]}]
]


Here the positive shear stress ss is put onto the end face with a NeumannValue and the other end is set to no deflection with a DirichletCondition

{uif, vif, wif} = NDSolveValue[{
stressOperator[{Y, ν} /. materialParameters] == {0, 0,
NeumannValue[ss, x == Len && -wd/2 <= y <= wd/2 && 0 <= z <= ht]},
DirichletCondition[u[x, y, z] == 0, x == 0],
DirichletCondition[v[x, y, z] == 0, x == 0],
DirichletCondition[w[x, y, z] == 0, x == 0]
},
{u, v, w}, {x, y, z} ∈ mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif, wif},
"ScalingFactor" -> 1];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]


As can be seen the deflection is in the direction of the shear force so this is as expected.

Now I change the target and put the shear force at x = 0 and clamp the end at x = Len.

{uif, vif, wif} = NDSolveValue[{
stressOperator[{Y, \[Nu]} /. materialParameters] == {0, 0,
NeumannValue[ss, x == 0]},
DirichletCondition[u[x, y, z] == 0, x == Len],
DirichletCondition[v[x, y, z] == 0, x == Len],
DirichletCondition[w[x, y, z] == 0, x == Len]
},
{u, v, w}, {x, y, z} \[Element] mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif, wif},
"ScalingFactor" -> 1];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]


The deflection is again in the direction of the shear force. So although the normal is pointing in the opposite direction the shear force is still deflecting upwards. To me this suggests that the shear stress is following the global coordinates.

This is my guess at the rule. A shear stress applied on the boundary should be resolved in the direction of the global coordinates and takes its positive direction from the positive direction of the global coordinates.

A possible test would be to rotate the region through 45deg and then see what happens. Can I rely on my rule so that I get the stresses in the correct direction?

EDIT

I have just checked the simple case of a positive normal force applied on the surface of a cuboid which has an outward normal in the negative direction of the global coordinates. This is the configuration

Show[
Region[reg, Axes -> True],
Graphics3D[{Black,
Table[Arrow[{{0, y, ht/2}, {-ht, y, ht/2}}], {y, -wd/2, wd/2, wd/
8}]}]
]


The arrows are drawn on the surface at x = 0 with outward normal in the negative x direction.

This is what happens.

{uif, vif, wif} = NDSolveValue[{
stressOperator[{Y, ν} /.
materialParameters] == {NeumannValue[ss,
x == 0 && -wd/2 <= y <= wd/2 && 0 <= z <= ht], 0, 0},
DirichletCondition[u[x, y, z] == 0, x == Len],
DirichletCondition[v[x, y, z] == 0, x == Len],
DirichletCondition[w[x, y, z] == 0, x == Len]
},
{u, v, w}, {x, y, z} ∈ mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif, wif},
"ScalingFactor" -> 100];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]


If the rule in the documentation was correct then the deflection would be outward but the deflection is inwards. Thus my rule of resolving applied stresses along the global coordinates works for normal and shear forces. It is good to know that @user21 is working on this.

• In principal your thinking of the global reference frame is correct. However, there is not such rule implemented specifically, which means that this (correct) behavior needs to be explainable from the NeumannValue definition. I'd need to think about that.... Mar 15, 2019 at 10:13
• @user21 I agree it needs a lot of thinking! Good luck. While you are thinking about this you may also like to think about the positive direction of internal stress. Books often have a little cube with one normal and two shear stresses on each surface. However, taking the two surfaces for which the x-axis is perpendicular, for example, on which of these two surfaces is the stress presented in the final output?
– Hugh
Mar 15, 2019 at 10:24
• @user21 I have just added the simplest case of a normal stress. My rule still applies.
– Hugh
Mar 16, 2019 at 18:09
• @user21 (2) It might help if we had a contribution from someone who does electromagnetism because that is a vector potential needing tensors as well (I think). I may try and develop my thinking with further examples. Overall I think you have done a good job of adding finite elements to Mathematica.
– Hugh
Mar 18, 2019 at 9:51
• Thanks for the kind words. In the section The Relation between NeumannValue and Boundary Derivatives in the tutorial there is a somewhat related issue. It explains the relation of NeumannValue to the use of Derivative on the boundary. I think I'd start a more mathematical description of what your are seeing from that one. Mar 18, 2019 at 12:40

As I showed in this answer 225595, the appropriate momentum flux variable is traction or stress vector and not the stress tensor. Therefore, the normal must be dotted with the stress tensor to obtain the traction vector:

$${{\mathbf{T}}^{(\hat {\mathbf{n}})}} = \hat {\mathbf{n}} \cdot {\mathbf{\sigma }}$$

For scalar properties, such as energy and mass, the NeumannValue can simplify the formulation of a PDE on a complex geometry because we do not need to know the surface normals. We simply need to know whether the flux is into or out of the domain as shown in the following diagram.

Unfortunately, we require knowledge of the surface normal when considering vector properties such as momentum, so we lose that advantage with the NeumannValue. Just remember that we are using traction versus stress to define the NeumannValue.

For the following, I am going to reduce the dimension to a 2D problem, just to save some memory and to simplify mesh construction. In 2D, we can represent the principal stresses with the follow diagram.

The traction vector is the vector sum of the principal stresses on a plane as shown in the figure taken from the Hooke's Law Wiki Article.

# Construction of Meshes

I will use ToBoundaryElementMesh so I can use ElementMarkers to refer to boundaries. Let's create a normal rectangular mesh and one that is rotated $$45^{\circ}$$

Needs["NDSolveFEM"]
(* Function to create bmesh rectangle with markers *)
bmeshrect[left_, bottom_, length_, height_, leftedge_, bottomedge_,
rightedge_, topedge_, angle_ : 0] :=
With[{right = left + length, top = bottom + height},
ToBoundaryMesh[
"Coordinates" -> RotationTransform[angle, First[#]] /@ # &@{{left,
bottom}(*1*), {right, bottom}(*2*), {right, top}(*3*), {left,
top}(*4*)},
"BoundaryElements" -> {LineElement[{{1, 2}(*bottom edge*)(*1*), {4,
1}(*left edge*)(*2*), {2, 3}(*3*)(*right edge*), {3,
4}(*4*)(*top edge*)}, {bottomedge, leftedge, rightedge,
topedge}]}]]
(* Unrotated *)
bmesh = bmeshrect[0, 0, 5, 1, 1, 2, 3, 2, 0 \[Degree]];
bmesh["Wireframe"["MeshElementMarkerStyle" -> Blue,
"MeshElementStyle" -> {Red, Green, Blue}, ImageSize -> Small]]
mesh = ToElementMesh[bmesh, MaxCellMeasure -> {"Length" -> .1}];
Show[mesh["Wireframe"],
mesh["Wireframe"["MeshElement" -> "PointElements",
"MeshElementMarkerStyle" -> Blue]]]
(* Rotated by 45 deg *)
bmesh = bmeshrect[0, 0, 5, 1, 1, 2, 3, 2, 45 \[Degree]];
bmesh["Wireframe"["MeshElementMarkerStyle" -> Blue,
"MeshElementStyle" -> {Red, Green, Blue}, ImageSize -> Small]]
mesh45 = ToElementMesh[bmesh, MaxCellMeasure -> {"Length" -> .1}];
Show[mesh45["Wireframe"],
mesh45["Wireframe"["MeshElement" -> "PointElements",
"MeshElementMarkerStyle" -> Blue]]]


# Set up Plane Stress Operator

(* set material parameters *)
materialParameters = {Y -> 10^3, ν -> 33/100};
(* set up factor matrix to be used in subsequent stress calcs *)
pfac = Y/(1 - ν^2)*{{1, ν, 0}, {ν, 1, 0}, {0,
0, (1 - ν)/2}};
fac = pfac /. materialParameters;
ClearAll[ν, Y]
parmop = {Inactive[
Div][({{0, -((Y ν)/(1 - ν^2))}, {-((Y (1 - ν))/(
2 (1 - ν^2))), 0}}.Inactive[Grad][v[x, y], {x, y}]), {x,
y}] + Inactive[
Div][({{-(Y/(1 - ν^2)),
0}, {0, -((Y (1 - ν))/(2 (1 - ν^2)))}}.Inactive[
Grad][u[x, y], {x, y}]), {x, y}],
Inactive[
Div][({{0, -((Y (1 - ν))/(2 (1 - ν^2)))}, {-((Y ν)/(
1 - ν^2)), 0}}.Inactive[Grad][u[x, y], {x, y}]), {x,
y}] + Inactive[
Div][({{-((Y (1 - ν))/(2 (1 - ν^2))),
v[x, y], {x, y}]), {x, y}]};
(* pde plane stress operator *)
op = parmop /. materialParameters;


# Solve PDE System for Unrotated and Rotated Meshes

## Fix Left Boundary

(* Fix left boundary *)
dcx = DirichletCondition[u[x, y] == 0., ElementMarker == 1];
dcy = DirichletCondition[v[x, y] == 0., ElementMarker == 1];
(* Create parametric functions of both meshes *)
pfun = ParametricNDSolveValue[{op == {NeumannValue[sx,
ElementMarker == 3], NeumannValue[sy, ElementMarker == 3]},
dcx, dcy}, {u, v}, {x, y} \[Element] mesh, {sx, sy}];
pfun45 = ParametricNDSolveValue[{op == {NeumannValue[sx,
ElementMarker == 3], NeumannValue[sy, ElementMarker == 3]},
dcx, dcy}, {u, v}, {x, y} \[Element] mesh45, {sx, sy}];
(* Set normal, stress, and traction *)
n = {1, 0};
stress = {{0, 1}, {1, 0}};
traction = n.stress;
(* Solve and display unrotated mesh *)
{ufun, vfun} = pfun @@ traction;
msh = ufun["ElementMesh"];
Show[{msh["Wireframe"["MeshElement" -> "BoundaryElements"]],
ElementMeshDeformation[msh, {ufun, vfun}][
"Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]
(* Solve and display rotated mesh *)
{ufun, vfun} = pfun45 @@ RotationTransform[45 \[Degree]][traction];
msh = ufun["ElementMesh"];
Show[{msh["Wireframe"["MeshElement" -> "BoundaryElements"]],
ElementMeshDeformation[msh, {ufun, vfun}][
"Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]


## Fix Right Boundary

(* Fix right boundary *)
dcx = DirichletCondition[u[x, y] == 0., ElementMarker == 3];
dcy = DirichletCondition[v[x, y] == 0., ElementMarker == 3];
(* Create parametric functions of both meshes *)
pfun = ParametricNDSolveValue[{op == {NeumannValue[sx,
ElementMarker == 1], NeumannValue[sy, ElementMarker == 1]},
dcx, dcy}, {u, v}, {x, y} \[Element] mesh, {sx, sy}];
pfun45 = ParametricNDSolveValue[{op == {NeumannValue[sx,
ElementMarker == 1], NeumannValue[sy, ElementMarker == 1]},
dcx, dcy}, {u, v}, {x, y} \[Element] mesh45, {sx, sy}];
(* Set normal, stress, and traction *)
n = {-1, 0};
stress = {{0, 1}, {1, 0}};
traction = n.stress;
(* Solve and display unrotated mesh *)
{ufun, vfun} = pfun @@ traction;
msh = ufun["ElementMesh"];
Show[{msh["Wireframe"["MeshElement" -> "BoundaryElements"]],
ElementMeshDeformation[msh, {ufun, vfun}][
"Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]
(* Solve and display rotated mesh *)
{ufun, vfun} = pfun45 @@ RotationTransform[45 \[Degree]][traction];
msh = ufun["ElementMesh"];
Show[{msh["Wireframe"["MeshElement" -> "BoundaryElements"]],
ElementMeshDeformation[msh, {ufun, vfun}][
"Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]


These results are consistent with the principal stress diagram and what I believe Hugh expected. I also showed that without too much effort, we can apply the approach to a rotated mesh. The key is to recognize that we need to convert stress into a traction vector.

# Extension to 3D

For completeness, I will extend this methodology to the 3D case.

## Create 3D Meshes

In 2D, Boundary Elements for NeumannValue and Point Elements for DirichletCondition are handled automatically if one first creates a bmesh with ToBoundaryMesh. In 3D, I did not find an analogous approach. However, one can use PointMarkerFunction and BoundaryMarkerFunction to assign nodal and elemental markers. The following workflow will create a mesh where we mark the left and right face:

(* Create 3D Mesh *)
(* 3D Parameters *)
Len = 1;  (*length *)
ht = 0.125; (* height *)
wd = 0.5; (* width *)
(* Material Parameters *)
materialParameters = {Y -> 10^3, \[Nu] -> 33/100};
ss = 1;(*Shear stress on beam*)
reg = Cuboid[{0, -wd/2, 0}, {Len, wd/2, ht}];
{{xmn, xmx}, {ymn, ymx}, {zmn, zmx}} = RegionBounds@reg;
(* Point Marker Function for Dirichlet Conditions *)
pointMarkerFunction =
Compile[{{coords, _Real, 2}, {pMarker, _Integer, 1}},
Block[{x = #[[1]], y = #[[2]], z = #[[3]], epsilon},
epsilon = 10^-6.;
Which[
Abs[x - xmn] <= epsilon, 1,
Abs[x - xmx] <= epsilon, 3,
Abs[y - ymn] <= epsilon, 2,
Abs[y - ymx] <= epsilon, 2,
Abs[z - zmn] <= epsilon, 2,
Abs[z - zmx] <= epsilon, 2,
True, 0
]] & /@ coords];
(* Boundary Marker Function for NeumannValues *)
boundaryMarkerFunction =
Compile[{{boundaryElementCoords, _Real,
3}, {boundaryElementPointMarkers, _Integer, 2}}, Which[
Union[#] == {1}, 1,
Union[#] == {3}, 3,
True, 2 ] & /@ boundaryElementPointMarkers];
mesh3D = ToElementMesh[reg,
"PointMarkerFunction" -> pointMarkerFunction,
"BoundaryMarkerFunction" -> boundaryMarkerFunction];
groups = mesh3D["BoundaryElementMarkerUnion"]
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp
mesh3D["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]
(* Rotate Mesh 45 deg *)
mesh3D45 =
ToElementMesh[
"Coordinates" ->
RotationTransform[-Pi/4, {0, 1, 0}, First[#]] /@ # &@
mesh3D["Coordinates"], "MeshElements" -> mesh3D["MeshElements"],
"BoundaryElements" -> mesh3D["BoundaryElements"],
"PointElements" -> mesh3D["PointElements"]];
groups = mesh3D45["BoundaryElementMarkerUnion"]
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp
mesh3D45["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]


## Fix Left Edge and Positive $$\sigma_{xz}$$ on Right Edge

(* Hugh's 3D Stress Operator *)
ClearAll[stressOperator, stressOperatorDynamic, u, v, w, x, y, z, Y, \
\

ν]; stressOperator[{Y_, ν_}] := {Inactive[

Div][{{0, 0, -((Y*ν)/((1 - 2*ν)*(1 + ν)))}, {0, 0,
0}, {-Y/(2*(1 + ν)), 0, 0}}.Inactive[Grad][
w[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, -((Y*ν)/((1 - 2*ν)*(1 + ν))),
0}, {-Y/(2*(1 + ν)), 0, 0}, {0, 0, 0}}.Inactive[
v[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0,
0}, {0, -Y/(2*(1 + ν)), 0}, {0,
u[x, y, z], {x, y, z}], {x, y, z}],
Inactive[Div][{{0, 0, 0}, {0,
0, -((Y*ν)/((1 -

2*ν)*(1 + ν)))}, {0, -Y/(2*(1 + ν)),
0}}.Inactive[Grad][w[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, -Y/(2*(1 + ν)),
0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0,
0}, {0, 0,
0}}.Inactive[Grad][u[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-Y/(2*(1 + ν)), 0,
0}, {0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))),
0}, {0,
v[x, y, z], {x, y, z}], {x, y, z}],
Inactive[Div][{{0, 0, 0}, {0,
0, -Y/(2*(1 + ν))}, {0, -((Y*ν)/((1 -
v[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, 0, -Y/(2*(1 + ν))}, {0, 0,
0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0,
0}}.Inactive[
Grad][u[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[

Div][{{-Y/(2*(1 + ν)), 0, 0}, {0, -Y/(2*(1 + ν)), 0}, {0,

0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν)))}}.Inactive[
Grad][w[x, y, z], {x, y, z}], {x, y, z}]};
(* Fix Left Edge/NeumannValue on Right Edge *)
parmop = stressOperator[{Y, ν}];
op = parmop /. materialParameters;
dcx = DirichletCondition[u[x, y, z] == 0., ElementMarker == 1];
dcy = DirichletCondition[v[x, y, z] == 0., ElementMarker == 1];
dcz = DirichletCondition[w[x, y, z] == 0., ElementMarker == 1];
nv = {NeumannValue[sx, ElementMarker == 3],
NeumannValue[sy, ElementMarker == 3],
NeumannValue[sz, ElementMarker == 3]};
pfun = ParametricNDSolveValue[{op == nv, dcx, dcy, dcz}, {u, v,
w}, {x, y, z} ∈ mesh3D, {sx, sy, sz}];
pfun45 = ParametricNDSolveValue[{op == nv, dcx, dcy, dcz}, {u, v,
w}, {x, y, z} ∈ mesh3D45, {sx, sy, sz}];


Now solve and view solutions of unrotated and rotated meshes:

(* Set normal, z-stress, and traction for right edge *)
n = {1, 0, 0};
stress = {{0, 0, 1}, {0, 0, 0}, {1, 0, 0}};
traction = n.stress;
(* Solve and display unrotated mesh *)
{ufun, vfun, wfun} = pfun @@ traction;
msh = ufun["ElementMesh"];
Show[{msh["Wireframe"["MeshElement" -> "BoundaryElements"]],
ElementMeshDeformation[msh, {ufun, vfun, wfun}][
"Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]},
ViewPoint -> {0, -Infinity, 0}]
(* Solve and display rotated mesh *)
{ufun, vfun, wfun} =
pfun45 @@ RotationTransform[-Pi/4, {0, 1, 0}][traction];
msh45 = ufun["ElementMesh"];
Show[{msh45["Wireframe"["MeshElement" -> "BoundaryElements"]],
ElementMeshDeformation[msh45, {ufun, vfun, wfun}][
"Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]},
ViewPoint -> {0, -Infinity, 0}]


## Fix Right Edge and Positive $$\sigma_{xz}$$ on Left Edge

(* Fix Left Edge/NeumannValue on Right Edge *)
dcx = DirichletCondition[u[x, y, z] == 0., ElementMarker == 3];
dcy = DirichletCondition[v[x, y, z] == 0., ElementMarker == 3];
dcz = DirichletCondition[w[x, y, z] == 0., ElementMarker == 3];
nv = {NeumannValue[sx, ElementMarker == 1],
NeumannValue[sy, ElementMarker == 1],
NeumannValue[sz, ElementMarker == 1]};
pfun = ParametricNDSolveValue[{op == nv, dcx, dcy, dcz}, {u, v,
w}, {x, y, z} \[Element] mesh3D, {sx, sy, sz}];
pfun45 = ParametricNDSolveValue[{op == nv, dcx, dcy, dcz}, {u, v,
w}, {x, y, z} \[Element] mesh3D45, {sx, sy, sz}];
(* Set normal, z-stress, and traction for left edge *)
n = {-1, 0, 0};
stress = {{0, 0, 1}, {0, 0, 0}, {1, 0, 0}};
traction = n.stress;
(* Solve and display unrotated mesh *)
{ufun, vfun, wfun} = pfun @@ traction;
msh = ufun["ElementMesh"];
Show[{msh["Wireframe"["MeshElement" -> "BoundaryElements"]],
ElementMeshDeformation[msh, {ufun, vfun, wfun}][
"Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]},
ViewPoint -> {0, -Infinity, 0}]
(* Solve and display rotated mesh *)
{ufun, vfun, wfun} =
pfun45 @@ RotationTransform[-Pi/4, {0, 1, 0}][traction];
msh45 = ufun["ElementMesh"];
Show[{msh45["Wireframe"["MeshElement" -> "BoundaryElements"]],
ElementMeshDeformation[msh45, {ufun, vfun, wfun}][
"Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]},
ViewPoint -> {0, -Infinity, 0}]


The approach works for 3D as well, but setting up boundary markers is more involved.

• Beautifully compiled answer. Jul 13, 2020 at 3:55