Wolfram Language & System Documentation Center

BoundaryUnitNormal

NDSolve`FEM`

BoundaryUnitNormal

BoundaryUnitNormal[x,y,…]

represents an outward-pointing unit normal vector on a region.

Details and Options

BoundaryUnitNormal can be used to construct partial differential equation boundary conditions that depend on the unit normal vector of the boundary.

BoundaryUnitNormal can be used with NeumannValue, DirichletCondition and NIntegrate.
BoundaryUnitNormal can be generated by boundary conditions like AcousticAbsorbingValue, HeatFluxValue or MassOutflowValue.
BoundaryUnitNormal can be used to specify a tangential on the boundary.
BoundaryUnitNormal will evaluate to a vector of length of the embedding dimension of the region when the boundary condition is discretized.
Components of the boundary unit normal can be accessed with Indexed.
For finite element approximations, the PDE is multiplied with a test function and integrated over . Integration by parts gives . The integrand in the boundary integral is replaced with the NeumannValue .
When a PDE specifies the Neumann value as , BoundaryUnitNormal can be used to model instead by specifying the NeumannValue as .
Conversely, when a PDE specifies the Neumann value as , BoundaryUnitNormal can be used to model instead by specifying the NeumannValue as .
At internal boundaries of a region, the boundary unit normal is not uniquely defined.
The value of the boundary unit normal will be computed by solving with a Dirichlet condition of on all boundaries including internal boundaries over the entire region . The boundary unit normal is then the gradient of normalized with .

Examples

open all close all

Basic Examples (1)

Load the package:

Wolfram Language code: <<NDSolve`FEM`

Solve a Poisson equation on a unit Disk:

Wolfram Language code: solution = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 1, DirichletCondition[u[x, y] == 0, True]}, u, {x, y}∈Disk[]]

Visualize the solution:

Wolfram Language code: Plot3D[solution[x, y], {x, y}∈Disk[]]

Compute the total flux through the boundary of the region through a second-order approximation of the boundary region:

Wolfram Language code:

NIntegrate[BoundaryUnitNormal[x, y].Grad[-solution[x, y], {x, y}], {x, y}∈ToBoundaryMesh[solution["ElementMesh"], "MeshOrder" -> 2]]

Compute the total flux through the boundary of a subregion:

Wolfram Language code:

NIntegrate[BoundaryUnitNormal[x, y].Grad[-solution[x, y], {x, y}], {x, y}∈RegionBoundary[Rectangle[{-1 / 2, -1 / 2}, {1 / 2, 1 / 2}]]]

Scope (6)

Solve a Poisson equation on a unit Disk:

Wolfram Language code: solution = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 1, DirichletCondition[u[x, y] == 0, True]}, u, {x, y}∈Disk[]]

Compute the total flux through the outer boundary:

Wolfram Language code: -NIntegrate[Grad[solution[x, y], {x, y}].BoundaryUnitNormal[x, y], {x, y}∈RegionBoundary[Disk[]]]

Compute the total flux through the boundary of the region through a second-order approximation of the boundary region:

Wolfram Language code:

-NIntegrate[Grad[solution[x, y], {x, y}].BoundaryUnitNormal[x, y], {x, y}∈ToBoundaryMesh[solution["ElementMesh"], "MeshOrder" -> 2]]

Specify a differential equation and a region:

Wolfram Language code:

op = Inactive[Div][-0.075*Inactive[Grad][c[x], {x}], {x}] + 1;
Ω = Line[{{0}, {1}}];

On the left, set up a NeumannValue with a boundary unit normal and solve the equation:

Wolfram Language code:

solution1 = NDSolveValue[{
	op == NeumannValue[BoundaryUnitNormal[x].{-c[x]}, x == 0]}, c, x∈Ω]

Visualize the solution:

Wolfram Language code: Plot[solution1[x], {x}∈Ω]

For this domain, the equivalent of the BoundaryUnitNormal at the left is :

Wolfram Language code:

solution2 = NDSolveValue[{
	op == NeumannValue[{-1}.{-c[x]}, x == 0]}, c, x∈Ω]

Show that the solutions are equal:

Wolfram Language code: Plot[solution1[x] - solution2[x], {x}∈Ω]

Create a tangential for a NeumannValue:

Wolfram Language code:

Ω = RegionDifference[Rectangle[{0, 0}, {4, 4}], Disk[]];
solution = NDSolveValue[{Laplacian[u[x, y], {x, y}] == NeumannValue[Cross[BoundaryUnitNormal[x, y]].{1, 1}, x ^ 2 + y ^ 2 == 1], DirichletCondition[u[x, y] == 0, x == 4 && y == 0]}, u, {x, y}∈Ω];
StreamPlot[Evaluate[Grad[solution[x, y], {x, y}]], {x, y}∈Ω]

Make use of a Indexed component of a BoundaryUnitNormal to compute a NeumannValue:

Wolfram Language code:

Ω = RegionDifference[Rectangle[{0, 0}, {4, 4}], Disk[]];
solution = NDSolveValue[{Laplacian[u[x, y], {x, y}] == NeumannValue[{-Indexed[BoundaryUnitNormal[x, y], 2], 0}.{1, 1}, x ^ 2 + y ^ 2 == 1], DirichletCondition[u[x, y] == 0, x == 4 && y == 0]}, u, {x, y}∈Ω];
StreamPlot[Evaluate[Grad[solution[x, y], {x, y}]], {x, y}∈Ω]

Solve a Poisson equation on a geometry:

Wolfram Language code:

Ω = RegionDifference[RegionUnion[Disk[{0, 0}, 5], Rectangle[{0, 0}, {7, 5}]], Ellipsoid[{0, 0}, {3, 2}]];
solution = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 1, DirichletCondition[u[x, y] == 0, True]}, u, {x, y}∈Ω]

Plot the solution:

Wolfram Language code: Plot3D[solution[x, y], {x, y}∈Ω]

Compute the total flux through the boundary:

Wolfram Language code:

NIntegrate[BoundaryUnitNormal[x, y].Grad[-solution[x, y], {x, y}], {x, y}∈ToBoundaryMesh[solution["ElementMesh"], "MeshOrder" -> 2]]

Set up a NeumannValue with a boundary normal:

Wolfram Language code:

op = D[c[t, x], t] + Inactive[Div][-0.075*Inactive[Grad][c[t, x], {x}], {x}];
ifunNormal = NDSolveValue[{op == NeumannValue[-BoundaryUnitNormal[x].({-1} * c[t, x]), x == 0], c[0, x] == Exp[-100 * (x - 1 / 2) ^ 2]}, c, {t, 0, 2}, x∈Line[{{0}, {1}}]]

Use a Neumann 0 boundary condition and solve the equation again:

Wolfram Language code:

ifunNeumannZero = NDSolveValue[{op == NeumannValue[0, x == 0], c[0, x] == Exp[-100 * (x - 1 / 2) ^ 2]}, c, {t, 0, 2}, x∈Line[{{0}, {1}}]]

Inspect how the solutions start to differ over time:

Wolfram Language code:

Manipulate[Plot[{ifunNormal[t, x], ifunNeumannZero[t, x]}, {x, 0, 1}, PlotRange -> {0, 1}], {t, 0, 2}, SaveDefinitions -> True]

For this domain, the equivalent of the BoundaryUnitNormal at the left is :

Wolfram Language code:

ifunNormal2 = NDSolveValue[{op == NeumannValue[-{-1}.({-1} * c[t, x]), x == 0], c[0, x] == Exp[-100 * (x - 1 / 2) ^ 2]}, c, {t, 0, 2}, x∈Line[{{0}, {1}}]]

Check that the solutions are the same:

Wolfram Language code: Plot3D[ifunNormal[t, x] - ifunNormal2[t, x], {x, 0, 1}, {t, 0, 2}]

Applications (1)

The following example considers a Stokes flow with prescribed traction boundary conditions in an annulus region. On the outer edge, at , there are no-slip boundary conditions. These set the fluid velocity in the and directions to 0, . On the inner boundary, at , a traction is prescribed. The traction with the stress vector is given as:

Here, is the radial unit vector given as and is the tangent unit vector given as . and are the Cartesian unit vectors in the and directions, respectively. For this example, the traction is set in normal direction to and the tangential traction is set to . In other words, at the inner boundary the velocity is not specified, only the traction is specified.

The fluid stress tensor is given by:

Stokes equation is given by:

and the continuity equation:

Set up parameters, geometry and refined mesh:

Wolfram Language code:

pars = <|a -> 1 / 5, b -> 1, μ -> 1|>;
Ω = Annulus[{0, 0}, {a, b}] /. pars;
mesh = ToElementMesh[Ω, "RegionHoles" -> {0, 0}, MeshRefinementFunction -> Function[{vertices, area}, area > 0.000025 (0.1 + 20 Norm[Mean[vertices]])]];

Set up a Stokes flow operator:

Wolfram Language code:

op = {Inactive[Div][{{-2 * μ, 0}, {0, -μ}}.Inactive[Grad][u[x, y], {x, y}], {x, y}] + 
	Inactive[Div][{{0, 0}, {-μ, 0}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] + Inactive[Div][{1, 0} * p[x, y], {x, y}], 
	Inactive[Div][{{0, -μ}, {0, 0}}.Inactive[Grad][u[x, y], {x, y}], {x, y}] + 
	Inactive[Div][{{-μ, 0}, {0, -2 * μ}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] + 
	Inactive[Div][{0, 1} * p[x, y], {x, y}], 
	Derivative[1, 0][u][x, y] + Derivative[0, 1][v][x, y]
	} /. pars;

Set up a no-slip condition on the outer boundary:

Wolfram Language code: ΓWall = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, ElementMarker == 2]

On the inner surface, you can compute the tangent unit normal from the outward-pointing unit normal with a cross product. Use a cross product to compute the unit tangent from the unit normal vector :

Wolfram Language code: Cross[{Subscript[n, x], Subscript[n, y]}]

To extract the first and second components of the cross product for the and directions, respectively, Indexed is used.

Set up the traction boundary conditions:

Wolfram Language code:

Γu = NeumannValue[Indexed[Cross[BoundaryUnitNormal[x, y]], 1] * ((x / a) ^ 2 - (y / a) ^ 2), ElementMarker == 1];
Γv = NeumannValue[Indexed[Cross[BoundaryUnitNormal[x, y]], 2] * ((x / a) ^ 2 - (y / a) ^ 2), ElementMarker == 1];

Equivalently, note that the tangent unit normal can also be given as because the unit normal vector on the inner surface is and the unit tangent then is :

Wolfram Language code:

NeumannValue[(y / a) * ((x / a) ^ 2 - (y / a) ^ 2), ElementMarker == 1];
NeumannValue[-(x / a) * ((x / a) ^ 2 - (y / a) ^ 2), ElementMarker == 1];

Without specifying a boundary condition for the pressure, the pressure value will be floating and NDSolve will give a warning that not enough boundary conditions are specified:

Wolfram Language code:

{ufun, vfun, pfun} = NDSolveValue[{op == {Γu, Γv, 0} /. pars, ΓWall}, {u, v, p}, {x, y}∈mesh, Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}, DependentVariables -> {u, v, p}]

Visualize the pressure solution and the velocity field:

Wolfram Language code: Show[{DensityPlot[pfun[x, y], {x, y}∈mesh, ...], StreamPlot[{ufun[x, y], vfun[x, y]}, {x, y}∈mesh, ...]}, ...]

To verify that the traction boundary condition works, you can compute radial and tangential components of the velocities and plot them at for all . For the radial velocity, a solution proportional to is expected, and for the tangential velocity, a solution proportional to is expected.

Define a function to compute the radial and tangential components of the velocities:

Wolfram Language code:

vr[x_, y_] := x / Sqrt[x ^ 2 + y ^ 2] * ufun[x, y] + y / Sqrt[x ^ 2 + y ^ 2] * vfun[x, y]
vt[x_, y_] := -(y / Sqrt[x ^ 2 + y ^ 2]) * ufun[x, y] + x / Sqrt[x ^ 2 + y ^ 2] * vfun[x, y]

Plot the radial component of the velocities and a function proportional to :

Wolfram Language code:

Plot[Evaluate[{vr[a * Cos[θ], a * Sin[θ]], -Sin[2θ] * a / 20} /. pars], {θ, 0, 2π}, PlotLegends -> {"vr", "∝ sin(2θ)"}]

Plot the tangential component of the velocities and a function proportional to :

Wolfram Language code:

Plot[Evaluate[{vt[a * Cos[θ], a * Sin[θ]], -Cos[2θ] * a / 5} /. pars], {θ, 0, 2π}, PlotLegends -> {"vt", "∝ cos(2θ)"}]

Next, compute the normal stress and the shear at the inner boundary and check that they match the prescribed traction boundary conditions. The stress tensor is given by:

The normal and tangential components of the traction vector can be computed by

and

Create a function to compute the stress tensor :

Wolfram Language code:

gradU = Grad[{ufun[x, y], vfun[x, y]}, {x, y}];
sigma[x_, y_] := Evaluate[(μ * (gradU + Transpose[gradU]) - pfun[x, y] * IdentityMatrix[2]) /. pars];

Compute the unit normal vector:

Wolfram Language code: normal = ElementMeshBoundaryUnitNormal[mesh]

Define a function for the tangent vector:

Wolfram Language code: tangent[x_, y_] := Cross[normal[x, y]]

Define a function to compute the normal stress:

Wolfram Language code: normalStress[x_ ? NumberQ, y_ ? NumberQ] := normal[x, y].(sigma[x, y].normal[x, y]);

Define a function to compute the shear stress:

Wolfram Language code: shearStress[x_ ? NumberQ, y_ ? NumberQ] := tangent[x, y].(sigma[x, y].normal[x, y]);

Visualize the computed normal and shear stress at :

Wolfram Language code:

stresses = {normalStress[a * Cos[θ], a * Sin[θ]], shearStress[a * Cos[θ], a * Sin[θ]]} /. pars;
Plot[stresses, {θ, 0, 2 Pi}, ...]

Visualize the difference of the normal and shear stress at from the boundary conditions:

Wolfram Language code: Plot[Evaluate[stresses - {0, Cos[2θ]}], {θ, 0, 2 Pi}, ...]

Properties & Relations (2)

Compute the boundary unit normal field:

Wolfram Language code: boundaryUnitNormal = ElementMeshBoundaryUnitNormal[ToElementMesh[Disk[]]]

Visualize the boundary unit normal:

Wolfram Language code: VectorPlot[boundaryUnitNormal[x, y], {x, y}∈Disk[]]

Visualize the boundary unit normal on the boundary only:

Wolfram Language code: VectorPlot[boundaryUnitNormal[x, y], {x, y}∈Disk[], VectorPoints -> ToBoundaryMesh[Disk[]]["Coordinates"]]

The boundary unit normal is computed by solving a Poisson equation over the region and specifying 0 Dirichlet conditions. Compute a Poisson equation over a unit Disk:

Wolfram Language code:

potential = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 1, DirichletCondition[u[x, y] == 0, True]}, u[x, y], {x, y}∈Disk[], "ExtrapolationHandler" -> {Automatic, "WarningMessage" -> False}];

Compute the normalized gradient of the potential:

Wolfram Language code: unitNormal = Normalize[Grad[potential, {x, y}], Sqrt[Total[# ^ 2]]&];

Visualize the unit normal:

Wolfram Language code: VectorPlot[unitNormal, {x, y}∈Disk[]]

This is the same as computing the boundary unit normal of an ElementMesh:

Wolfram Language code:

boundaryUnitNormal = ElementMeshBoundaryUnitNormal[ToElementMesh[Disk[]]];
VectorPlot[boundaryUnitNormal[x, y], {x, y}∈Disk[]]

Top

More Learning

Tech Support

Wolfram Solutions

Wolfram Solutions For Education

Get Started

Grow Your Skills

Work with Us

Educational Programs for Adults

Educational Programs for Youth

Read

BoundaryUnitNormal

Details and Options

Examples

Basic Examples (1)

Scope (6)

Applications (1)

Properties & Relations (2)

Text

CMS

APA

BibTeX

BibLaTeX

BoundaryUnitNormal

Details and Options

Examples

Basic Examples (1)

Scope (6)

Applications (1)

Properties & Relations (2)

See Also

Text

CMS

APA

BibTeX

BibLaTeX