Finite Element Method

How to talk to it and make it your friend.

Paritosh Mokhasi,
Kernel Developer, WRI.

proxying for Oliver.

Please send all complaints, protests and grievances to Oliver.
Please send all compliments, acclaim and monetary rewards to Nina and Paritosh. 

We have been trained by the Numerics Jedi Master Rob Knapp.

Join the Conversation #WolframTechConf

NDSolve capability

Where to get detailed FEM information?

Tutorials - in order of difficulty

Solving PDE with FEM

Finite Element Method Usage Tips

Element Mesh Generation / Visualization

Finite Element Programming

Wolfram Technology Conferences

PDEs and Finite Elements

Workshop: PDEs and Finite Elements

Applications of Numerical Computation

Writing Your Own PDE Solvers - The Finite Element Method Interface

Videos (Virtual Conference)

Solving PDEs Over Regions with Symbolics and Numerics

Solving Engineering Problems with Mathematica’s PDE Tools

Stackexchange Mathematica

FEM

Differential Equations

How to use FEM in NDSolve?

({{-(∇^2)u=1}, {Ω:Disk}, {Γ:u=0|x≥0}})

Setup: Specify equation, domain and boundary conditions.

eqn=-(InactiveDiv[(InactiveGrad[u[x,y], {x,y}, Inactive]), {x,y}, Inactive])==1;
Ω=Disk[];
(Subscript[Γ, D])=DirichletCondition[u[x,y]==0,x≥0 || x^2+y^2==1/4];

Solve

ufun=NDSolveValue[{eqn,(Subscript[Γ, D])},u,{x,y}∈Ω]

Visualize

Plot3D[ufun[x,y],{x,y}∈Ω]

What can NDSolve-FEM solve?

d(∂/(∂t))u+∇·(-c∇u-α u+γ)+β·∇u+a u=f

c		Matrix
α, β, γ	Vector
d,a,f		Scalar

Laplace:                                                 ∇·(-c∇u-α u+γ)+β·∇u+a u=0

Poisson:                                                ∇·(-c∇u-α u+γ)+β·∇u+a u=f

Helmholtz:                                       ∇·(-c∇u-α u+γ)+β·∇u+a u=f

Convection/Reaction:              ∇·(-c∇u-α u+γ)+β·∇u+a u=f

Heat:                                         (∂/(∂t))u+∇·(-c∇u-α u+γ) +β·∇u+a u=f

Wave:                        ((∂^2)/(∂(t^2)))u+(∂/(∂t))u+∇·(-c∇u-α u+γ)+β·∇u+a u=0

Systems of PDEs

Boundary Conditions: Dirichlet and Neumann

h u= r  on  ∂Ω
(OverVector[n])·(c ∇u + α u - γ)= g - q u  on  ∂Ω
g - q u= 0  default

Dirichlet conditions

DirichletCondition[u[x,y,..], pred]

({{-(∇^2)u=1}, {Ω:Disk}, {Γ:u=0|x≥0}})

Setup:

eqn=-(InactiveDiv[(InactiveGrad[u[x,y], {x,y}, Inactive]), {x,y}, Inactive])==1;
Ω=Disk[];
(Subscript[Γ, D])=DirichletCondition[u[x,y]==0,x≥0];

Solve and Visualize:

ufun=NDSolveValue[{eqn,(Subscript[Γ, D])},u,{x,y}∈Ω];
Plot3D[ufun[x,y],{x,y}∈Ω]

Dirichlet conditions: Multiple materials

({{∇({{{, ({{(x^2)+(y^2)≤1/4, -10}, {True, -1}})}})∇u=1}, {Ω:Disk}, {Γ:u=0|x≥0}})

Setup:

c=If[x^2+y^2≤1/4,{{10,0},{0,10}},{{1,0},{0,1}}];
eqn=(InactiveDiv[(-c.(InactiveGrad[u[x,y], {x,y}, Inactive])), {x,y}, Inactive])==1;
Ω=Disk[];
(Subscript[Γ, D])=DirichletCondition[u[x,y]==0,x≥0];

Solve & Visualize:

ufun=NDSolveValue[{eqn,(Subscript[Γ, D])},u,{x,y}∈Ω];
Plot3D[ufun[x,y],{x,y}∈Ω]

Dirichlet conditions: Internal domains.

({{-(∇^2)u=1}, {Ω:Disk}, {Γ:u=0|(x≥0|| (x^2)+(y^2)==1/4)}})

Setup:

eqn = Laplacian[u[x,y],{x,y}]==-1;
Ω = NDSolve`FEM`ToElementMesh[ImplicitRegion[x^2+y^2==1 ||x^2+y^2==1/4,{x,y}],"RegionHoles"->None,"MeshOrder"->1];
(Subscript[Γ, D])=DirichletCondition[u[x,y]==0,x≥0 ||x^2+y^2==1/4];

Solve & Visualize:

ufun=NDSolveValue[{eqn,(Subscript[Γ, D])},u,{x,y}∈Ω];
Plot3D[ufun[x,y],{x,y}∈Ω]

What should we expect for this case?

ufun=NDSolveValue[{eqn,DirichletCondition[u[x,y]==0,x≥0]},u,{x,y}∈Ω];
Plot3D[ufun[x,y],{x,y}∈Ω]

Neumann conditions

NeumannValue[g - q*u[x,y,..], pred]

FEM derivation:

(∫_Ω)∇·(-c ∇u-α u+γ) ϕ+β·∇u ϕ+a u ϕ-f ϕdx=0

Integration by parts gives weak form:

(∫_Ω)(c ∇u+α u-γ) ·∇ϕ+β·∇u ϕ+a u ϕ-f ϕ dx - (∫_(∂Ω))(OverVector[n])·(c ∇u+α u-γ)ϕdx=0

Substitute the definition of Neumann value and take is on the RHS:

(∫_Ω)(c ∇u+α u-γ) ·∇ϕ+β·∇u ϕ+a u ϕ-f ϕ dx = (∫_(∂Ω))(g-q u) ϕdx

See also: Finite Element Method Usage Tips.

NeumannValue

({{3u''(x)+15u'(x)+200u(x)=0,}, {(Subscript[Γ, D]): u(1) = 1,
(Subscript[Γ, N]): u(2) - 4u'(2)=5}})

Recall that the Neumann conditions are enforced in the system as:

(∫_Ω)(c ∇u+α u-γ) ·∇ϕ+β·∇u ϕ+a u ϕ-f ϕ dx = (∫_(∂Ω))(OverVector[n])·(c ∇u+α u-γ)ϕdx,
(∫_Ω)(c ∇u+α u-γ) ·∇ϕ+β·∇u ϕ+a u ϕ-f ϕ dx = (∫_(∂Ω))(g-q u) ϕdx

dsol = DSolveValue[{3u''[x]+15u'[x]+200u[x]==0,u[1]==1,u[2]-4u'[2]==5},u,x];
Plot[dsol[x],{x,1,2}, PlotRange->All]

Ω=Line[{{1},{2}}];

op1 = 3u''[x]+15u'[x]+200 u[x];

Equation: ∇·(3 ∇u ) +15·∇u + 200 u=0,
(OverVector[n])·(-3 ∇u )= g - q u,
-3 u'(2) = g - q u,

To Match (Subscript[Γ, N]): - 4u'(2)=5-u(2)  with the Neumann condition in the system we need g = ((3\[Cross]5)/4), q = (3/4)

femsol1 = NDSolveValue[{op1==NeumannValue[(3/4)*(5-u[x]),x==2], u[1]==1},u,{x}∈Ω];
Plot[(femsol1[x]-dsol[x]),{x,1,2}, PlotRange->All]

NeumannValue

({{3u''(x)+15u'(x)+200u(x)=0,}, {(Subscript[Γ, D]): u(1) = 1,
(Subscript[Γ, N]): u(2) - 4u'(2)=5}})

op2=Inactive[Div][{{3}}. Inactive[Grad][u[x],{x}]+Inactive[Times][{15},u[x]],{x}]+200*u[x];

Equation: ∇·(3 ∇u  + 15 u) + 200 u=0,
(OverVector[n])·(-3 ∇u -15u)= g - q u,
-3 u'(2)-15u(2) = g - q u,

To Match (Subscript[Γ, N]): - 4u'(2)=5-u(2)  with the Neumann condition in the system we need g = ((3\[Cross]5)/4), q = (15+(3/4))

femsol2 =NDSolveValue[{op2 == NeumannValue[(3/4)(5-21u[x]),x==2],u[1]==1},u,x∈Ω];
Plot[femsol2[x]-dsol[x],{x,1,2}]

NeumannValue: Multiple boundary conditions.

({{-232.6(∇^2)u=0}, {Ω:Mesh}, {(Subscript[Γ, N1]):(x/100+2)*58.15*(273.15+20-u)|EM1}, {(Subscript[Γ, N2]):1000|EM2}})

bmesh=;

Visualize the location of each element marker:

bmesh1 = NDSolve`FEM`ToBoundaryMesh[bmesh];
bmesh1["BoundaryElementMarkerUnion"]
Show[bmesh1["Wireframe"[ElementMarker==1,"MeshElementStyle"->EdgeForm[Black]]],bmesh1["Wireframe"[ElementMarker==2,"MeshElementStyle"->EdgeForm[Red]]]]

Solve and visualize:

uif=NDSolveValue[-232.6 Laplacian[u[x,y,z],{x,y,z}]==NeumannValue[1000,ElementMarker==2]+NeumannValue[(x/100+2) 58.15*(273.15+20-u[x,y,z]),ElementMarker==1],u,{x,y,z}∈bmesh,Method->{"PDEDiscretization"->{"FiniteElement","MeshOptions"->{"MaxCellMeasure"->10}}}];
NDSolve`FEM`ElementMeshSurfacePlot3D[uif,Boxed->False,Lighting->"Neutral",ViewPoint->{-2,-2,-1.6}]

NeumannValue is part of the system because...

2 PDEs with dependent variables u and v:

∇·(-(Subscript[c, 11]) ∇u -(Subscript[c, 12])∇v -(Subscript[α, 11]) u +(Subscript[α, 12]) v+ (Subscript[γ, 1]) ) +(Subscript[β, 11])·∇u +(Subscript[β, 12])·∇v + (Subscript[a, 11]) u+ (Subscript[a, 12]) v ==(Subscript[f, 1])
∇·(-(Subscript[c, 21]) ∇u -(Subscript[c, 22])∇v -(Subscript[α, 21]) u +(Subscript[α, 22]) v+ (Subscript[γ, 2]) ) +(Subscript[β, 21])·∇u +(Subscript[β, 22])·∇v + (Subscript[a, 21]) u+ (Subscript[a, 22]) v==(Subscript[f, 2])

In this case, the generalized Neumann conditions are

(OverVector[n])·((Subscript[c, 11]) ∇u +(Subscript[c, 12])∇v)= (Subscript[g, 1]) - (Subscript[q, 11]) u  - (Subscript[q, 12]) v 
(OverVector[n])·((Subscript[c, 21]) ∇u +(Subscript[c, 22])∇v)= (Subscript[g, 2]) - (Subscript[q, 21]) u  - (Subscript[q, 22]) v

(* plane stress *)
Thread[({(InactiveDiv[({{-(Y/(1-(ν^2))),0},{0,-((Y (1-ν))/(2 (1-(ν^2))))}}.(InactiveGrad[u[x,y], {x,y}, Inactive])), {x,y}, Inactive])+(InactiveDiv[({{0,-((Y ν)/(1-(ν^2)))},{-((Y (1-ν))/(2 (1-(ν^2)))),0}}.(InactiveGrad[v[x,y], {x,y}, Inactive])), {x,y}, Inactive]),(InactiveDiv[({{0,-((Y (1-ν))/(2 (1-(ν^2))))},{-((Y ν)/(1-(ν^2))),0}}.(InactiveGrad[u[x,y], {x,y}, Inactive])), {x,y}, Inactive])+(InactiveDiv[({{-((Y (1-ν))/(2 (1-(ν^2)))),0},{0,-(Y/(1-(ν^2)))}}.(InactiveGrad[v[x,y], {x,y}, Inactive])), {x,y}, Inactive])}/.{Y->1,ν->1/3})=={0,NeumannValue[-1,x==5]}//Activate]//TableForm

Since NeumannValue[-1,x==5] does not depend on u or v  we want to associate it with an equation!

Error Analysis: Better Accuracy using Curved Elements.

reg = Disk[];
(π - NIntegrate[1,{x,y}∈#])&/@{reg, DiscretizeRegion[reg], ToElementMesh[reg]}

reg = ImplicitRegion[x^2+y^2≤1,{x,y}];
(π - NIntegrate[1,{x,y}∈#])&/@{reg, DiscretizeRegion[reg], ToElementMesh[reg]}

Why is there a variation in the results?

mr=DiscretizeRegion[Disk[],MaxCellMeasure->{"Area"->Infinity},AccuracyGoal->1];
Show[Graphics[Circle[]],mr]

mesh=ToElementMesh[Disk[],MaxCellMeasure->{"Area"->Infinity},AccuracyGoal->1];
Show[Graphics[Circle[]],ElementMeshCurvedWireframe[mesh],mesh["Wireframe"["MeshElement"->"PointElements","MeshElementStyle"->Directive[{PointSize[0.02],Red}]]]]

Error Analysis: Effect of discretization on PDEs

(∇^2)u = f(x,y), 
u((x,y)∈Γ) = 0.
Γ := (x^2)+(y^2)≤2 && (((x-(1/2)))^2)+(y^2) ≥  1/2,
(Subscript[u, exact])(x,y)=(2-((x^2)+(y^2)))((((x-(1/2)))^2) + (y^2) - (1/2))(1+Sin(π x)Sin(π y))

exactSol[x_,y_]:= (2-(x^2+y^2))*((x-1/2)^2 + y^2 - 1/2)*(1+Sin[Pi*x]*Sin[Pi*y]);
rhs[x_,y_]=Laplacian[exactSol[x,y],{x,y}];

reg = ImplicitRegion[x^2+y^2≤2 && (x-1/2)^2+y^2>=1/2,{x,y}];
Plot3D[exactSol[x,y],{x,y}∈reg, Axes->True,PlotRange->All]

Case 1: region is discretized -> boundary is first order accurate.

errd = Table[
mr = DiscretizeRegion[reg, MaxCellMeasure-> {"Length"-> 0.2/i}];
usol = NDSolveValue[{Laplacian[u[x,y],{x,y}]==rhs[x,y], DirichletCondition[u[x,y]==0,True]},u,{x,y}∈mr];
emesh = usol["Coordinates"][[1]];
h = Mean[PropertyValue[{mr,1},MeshCellMeasure]];
erri =  Mean[Abs[usol["ValuesOnGrid"] - exactSol@@Transpose[emesh["Coordinates"]]]];
{h,erri}
,{i,{1,2,4,8,16}}];
refErr = NestList[#/{2,4}&,errd[[1]],4];
ListLogLogPlot[{errd,refErr}, Joined->True, PlotStyle->{Red,{Blue,Dashed}},PlotMarkers->Automatic, PlotLegends->{"FEM","O((h^2))"}, Axes->False, Frame->True, FrameLabel -> {"Average vertex sep. dist (h)","Mean Abs. Error"}, FrameStyle->Directive[Black,14], ImageSize->500]

Case 2: region is symbolic -> boundary is second order accurate.

err = Table[
usol = NDSolveValue[{Laplacian[u[x,y],{x,y}]==rhs[x,y], DirichletCondition[u[x,y]==0,True]},u,{x,y}∈reg, Method->{"PDEDiscretization"-> {"FiniteElement", "MeshOptions"-> {"MaxCellMeasure"-> {"Length"-> 0.2/i}}}}];
emesh = usol["Coordinates"][[1]];
h = Mean[PropertyValue[{MeshRegion[emesh],1},MeshCellMeasure]];
erri =  Mean[Abs[usol["ValuesOnGrid"] - exactSol@@Transpose[emesh["Coordinates"]]]];
{h,erri}
,{i,{1,2,4,8,16}}];
refErr2ndOrder = NestList[#/{2,2^2}&,errd[[1]]+{0,0.005},4];
refErr4thOrder = NestList[#/{2,2^3}&,err[[1]] + {0, 0.0001},4];
ListLogLogPlot[{err,errd,refErr2ndOrder, refErr4thOrder}, Joined->True, PlotStyle->{{Red,Thick},{Brown,Thick},{Blue,Dashed}, {Black,Dashed}},PlotMarkers->Automatic, PlotLegends->{"FEM - Symbolic Region","FEM - Discretized Region","O((h^2))","O((h^4))"},Axes->False, Frame->True, FrameLabel -> {"Average vertex sep. dist (h)","Mean Abs. Error"}, FrameStyle->Directive[Black,14], ImageSize->500]

Error Analysis: Effect of discretization on PDEs

(∇^2)u = f(x,y), 
u((x,y)∈Γ) = 0.
Γ := -1 ≤ x≤1 && -1 ≤ y≤1,
(Subscript[u, exact])(x,y)=(2-((x^2)+(y^2)))((((x-(1/2)))^2) + (y^2) - (1/2))(Sin(π x)Sin(π y))

exactSol[x_,y_]:= (2-(x^2+y^2))*((x-1/2)^2 + y^2 - 1/2)*(Sin[Pi*x]*Sin[Pi*y]);
rhs[x_,y_]=Laplacian[exactSol[x,y],{x,y}];
reg = Rectangle[{-1,-1},{1,1}];

errd = Table[
mr = DiscretizeRegion[reg, MaxCellMeasure-> {"Length"-> 0.2/i}];
usol = NDSolveValue[{Laplacian[u[x,y],{x,y}]==rhs[x,y], DirichletCondition[u[x,y]==0,True]},u,{x,y}∈mr];
emesh = usol["Coordinates"][[1]];
h = Mean[PropertyValue[{mr,1},MeshCellMeasure]];
erri =  Mean[Abs[usol["ValuesOnGrid"] - exactSol@@Transpose[emesh["Coordinates"]]]];
{h,erri}
,{i,{1,2,4,8,16}}];
refErr2ndOrder = NestList[#/{2,2^2}&,errd[[1]] ,4];
refErr4thOrder = NestList[#/{2,2^3}&,errd[[1]],4];
p1 = Plot3D[exactSol[x,y],{x,y}∈reg, Axes->True,PlotRange->All, ImageSize->400];
{p1, ListLogLogPlot[{errd,refErr2ndOrder, refErr4thOrder}, Joined->True, PlotStyle->{Red,{Blue,Dashed},{Black,Dashed}},PlotMarkers->Automatic, PlotLegends->{"FEM","O((h^2))", "O((h^4))"}, Axes->False, Frame->True, FrameLabel -> {"Average vertex sep. dist (h)","Mean Abs. Error"}, FrameStyle->Directive[Black,14], ImageSize->500]}

Error Analysis: Residual Checking

When checking solutions by reinserting into equations, be careful: The interpolation is second order and the PDEs are too! Recovery of derivatives is difficult.

exact[x_,y_]:=(x^3+5y^2)
op=-Laplacian[u[x,y],{x,y}]+10u[x,y];
f=op/.u->exact;

eqn = op - f;
Ω=Disk[];

sol=NDSolveValue[{eqn==0,DirichletCondition[u[x,y]==exact[x,y],True]},u,{x,y}∈Ω];
Plot3D[sol[x,y],{x,y}∈Ω]

check=eqn/.u->sol;
GraphicsRow[{Plot3D[check,{x,y}∈Ω,Axes->True, ColorFunction->Automatic], Plot3D[sol[x,y]-exact[x,y],{x,y}∈Ω,Axes->True,PlotRange->All]}, ImageSize->800]

Thanks

Message from Oliver:

Init

Needs["NDSolve`FEM`"]
$HistoryLength=0;

Needs["NDSolve`FEM`"]
$HistoryLength=0;
SetOptions[#,{ColorFunction->"TemperatureMap",AspectRatio->Automatic}]& /@{ContourPlot,DensityPlot};
SetOptions[RegionPlot,AspectRatio->Automatic];
SetOptions[#,{ColorFunction->"TemperatureMap",AspectRatio->Automatic,Boxed->False,Axes->False}]&/@{Plot3D};

planeStress[Y_,ν_]:={(InactiveDiv[({{0,-((Y ν)/(1-(ν^2)))},{-((Y (1-ν))/(2 (1-(ν^2)))),0}}.(InactiveGrad[v[x,y], {x,y}, Inactive])), {x,y}, Inactive])+(InactiveDiv[({{-(Y/(1-(ν^2))),0},{0,-((Y (1-ν))/(2 (1-(ν^2))))}}.(InactiveGrad[u[x,y], {x,y}, Inactive])), {x,y}, Inactive]),(InactiveDiv[({{0,-((Y (1-ν))/(2 (1-(ν^2))))},{-((Y ν)/(1-(ν^2))),0}}.(InactiveGrad[u[x,y], {x,y}, Inactive])), {x,y}, Inactive])+(InactiveDiv[({{-((Y (1-ν))/(2 (1-(ν^2)))),0},{0,-(Y/(1-(ν^2)))}}.(InactiveGrad[v[x,y], {x,y}, Inactive])), {x,y}, Inactive])}


ElementMeshCurvedWireframe[mesh_ElementMesh]:=Block[{triIncidentsOrder,quadIncidentsOrder,interpolatingQuadBezierCurve,interpolatingQuadBezierCurveComplex,triEdges,triIncidents},

triIncidentsOrder={1,4,2,5,3,6};
quadIncidentsOrder={1,5,2,6,3,7,4,8};

interpolatingQuadBezierCurve[pts_List]/;Length[pts]==3:=BezierCurve[{pts[[1]],(1/2) (-pts[[1]]+4 pts[[2]]-pts[[3]]),pts[[3]]}];
interpolatingQuadBezierCurve[pts_List,"ControlPoints"]/;(Length[pts]==3):={pts[[1]],(1/2) (-pts[[1]]+4 pts[[2]]-pts[[3]]),pts[[3]]};

interpolatingQuadBezierCurve[ptslist_List]:=interpolatingQuadBezierCurve/@ptslist;

interpolatingQuadBezierCurveComplex[coords_,indices_]:=interpolatingQuadBezierCurve[Map[coords[[#]]&,indices]];

triEdges=Partition[triIncidentsOrder , 3, 2, {1,1} ];
triIncidents=Join@@ElementIncidents[mesh["MeshElements"]];
coords=mesh["Coordinates"];
Graphics[{
interpolatingQuadBezierCurveComplex[coords,triIncidents[[All, #]]]&/@triEdges
}]
]

NeumannValue: Extra slide.

({{3u''(x)+15u'(x)+200u(x)=0,}, {(Subscript[Γ, D]): u(1) = 1,
(Subscript[Γ, N]): u(2) - 4u'(2)=5}})

op3=Inactive[Div][{{3}}. Inactive[Grad][u[x],{x}]+Inactive[Times][{5},u[x]],{x}]+10*u'[x]+200*u[x];

Equation: ∇·(3 ∇u  + 5 u) + 10 ∇u + 200 u=0,
(OverVector[n])·(-3 ∇u -5u)= g - q u,
-3 u'(2)-5u(2) = g - q u,

To Match (Subscript[Γ, N]): - 4u'(2)=5-u(2)  with the Neumann condition in the system we need g = ((3\[Cross]5)/4), q = (5+(3/4))

femsol3 =NDSolveValue[{op3 == NeumannValue[(1/4)(15-23u[x]),x==2],u[1]==1},u,x∈Ω];
Plot[femsol3[x]-dsol[x],{x,1,2}]