Skip to content

Commit 8513f34

Browse files
pjaappjaap
committed
Fix inconsistency in the CombineDofs operator treating the coupling matrix
Add two examples to demonstate the feature.
1 parent 601a7ee commit 8513f34

4 files changed

Lines changed: 453 additions & 25 deletions

File tree

‎examples/Example212_PeriodicBoundary2D.jl‎

Lines changed: 180 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,180 @@
1+
#=
2+
3+
# 212 : Periodic Boundary 2D
4+
([source code](@__SOURCE_URL__))
5+
6+
This is a simple demonstration and validation of the generic periodic boundary operator.
7+
8+
We construct an unstructured periodic 2D grid and solve a simple linear elastic problem
9+
with periodic boundary coupling along the x-axis.
10+
11+
![](example212.png)
12+
=#
13+
14+
module Example212_PeriodicBoundary2D
15+
16+
using ExtendableFEM
17+
using ExtendableGrids
18+
using SimplexGridFactory
19+
using GridVisualize
20+
using Triangulate
21+
using UnicodePlots
22+
using StaticArrays
23+
using LinearAlgebra
24+
using Test #hide
25+
using SparseArrays
26+
27+
## enumerate the boundary regions
28+
const reg_left = 4
29+
const reg_right = 2
30+
const reg_dirichlet = 1
31+
const reg_default = 3
32+
33+
## 2D reduction of the material used in Example 312
34+
## in Voigt notation
35+
function material_tensor()
36+
c11 = 396.0
37+
c12 = 137.0
38+
c44 = 116.0
39+
40+
return @SArray [
41+
c11 c12 0
42+
c12 c11 0
43+
0 0 c44
44+
]
45+
end
46+
47+
## generate the kernels for the linear problem
48+
## 𝐂: Hooke tensor, 𝑓: body force
49+
function make_kernels(𝐂, 𝑓)
50+
51+
## linear stress-strain mapping
52+
bilinear_kernel!(σ, εv, qpinfo) = mul!(σ, 𝐂, εv)
53+
54+
## plain body force
55+
linear_kernel!(result, qpinfo) = (result .= 𝑓)
56+
57+
return bilinear_kernel!, linear_kernel!
58+
end
59+
60+
61+
"""
62+
create a 2D grid with Dirichlet boundary region at the bottom center
63+
"""
64+
function create_grid(; h, height, width)
65+
builder = SimplexGridBuilder(; Generator = Triangulate)
66+
67+
## bottom points
68+
b1 = point!(builder, 0, 0)
69+
b2 = point!(builder, 0.45 * width, 0)
70+
b3 = point!(builder, 0.5 * width, 0)
71+
b4 = point!(builder, 0.55 * width, 0)
72+
b5 = point!(builder, width, 0)
73+
74+
## top points
75+
t1 = point!(builder, 0, height)
76+
t2 = point!(builder, width / 2, height)
77+
t3 = point!(builder, width, height)
78+
79+
## default faces
80+
facetregion!(builder, reg_default)
81+
facet!(builder, b1, b2)
82+
facet!(builder, b4, b5)
83+
facet!(builder, t1, t2)
84+
facet!(builder, t2, t3)
85+
86+
## left face
87+
facetregion!(builder, reg_left)
88+
facet!(builder, b1, t1)
89+
90+
## right face
91+
facetregion!(builder, reg_right)
92+
facet!(builder, b5, t3)
93+
94+
## Dirichlet face
95+
facetregion!(builder, reg_dirichlet)
96+
facet!(builder, b3, b4)
97+
facet!(builder, b2, b3)
98+
99+
## divider
100+
facetregion!(builder, reg_default)
101+
facet!(builder, t2, b3)
102+
103+
104+
cellregion!(builder, 1)
105+
maxvolume!(builder, h)
106+
regionpoint!(builder, width / 3, height / 2)
107+
108+
cellregion!(builder, 2)
109+
# much finer grid on the right half to make periodic coupling non-trivial
110+
maxvolume!(builder, 0.1 * h)
111+
regionpoint!(builder, 2width / 3, height / 2)
112+
113+
return simplexgrid(builder)
114+
end
115+
116+
function main(;
117+
order = 1,
118+
periodic = true,
119+
Plotter = nothing,
120+
force = 10.0,
121+
h = 1.0e-2,
122+
width = 6.0,
123+
height = 1.0,
124+
kwargs...
125+
)
126+
xgrid = create_grid(; h, width, height)
127+
128+
## create finite element space and solution vector
129+
if order == 1
130+
FES = FESpace{H1P1{2}}(xgrid)
131+
elseif order == 2
132+
FES = FESpace{H1P2{2, 2}}(xgrid)
133+
end
134+
135+
## problem description
136+
PD = ProblemDescription()
137+
u = Unknown("u"; name = "displacement")
138+
assign_unknown!(PD, u)
139+
140+
𝐂 = material_tensor()
141+
𝑓 = force * [0, 1]
142+
143+
bilinear_kernel!, linear_kernel! = make_kernels(𝐂, 𝑓)
144+
assign_operator!(PD, BilinearOperator(bilinear_kernel!, [εV(u, 1.0)]; kwargs...))
145+
assign_operator!(PD, LinearOperator(linear_kernel!, [id(u)]; kwargs...))
146+
147+
assign_operator!(PD, HomogeneousBoundaryData(u; regions = [reg_dirichlet]))
148+
149+
if periodic
150+
function give_opposite!(y, x)
151+
y .= x
152+
y[1] = width - x[1]
153+
return nothing
154+
end
155+
156+
coupling_matrix = get_periodic_coupling_matrix(FES, xgrid, reg_left, reg_right, give_opposite!)
157+
display(coupling_matrix)
158+
assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
159+
end
160+
161+
sol = solve(PD, FES)
162+
163+
plt = GridVisualizer(; Plotter, size = (1300, 800))
164+
165+
magnification = 1
166+
displaced_grid = deepcopy(xgrid)
167+
displace_mesh!(displaced_grid, sol[1], magnify = magnification)
168+
gridplot!(plt, displaced_grid, linewidth = 1, title = "displaced mesh, $(magnification)x magnified", scene3d = :LScene)
169+
170+
return sol, plt
171+
end
172+
173+
generateplots = ExtendableFEM.default_generateplots(Example212_PeriodicBoundary2D, "example212.png") #hide
174+
function runtests() #hide
175+
sol, plt = main() #hide
176+
@test maximum(view(sol[1])) ≈ 1.3447465095618172 #hide
177+
return nothing #hide
178+
end #hide
179+
180+
end # module

‎examples/Example312_PeriodicBoundary3D.jl‎

Lines changed: 196 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,196 @@
1+
#=
2+
3+
# 312 : Periodic Boundary 3D
4+
([source code](@__SOURCE_URL__))
5+
6+
This is a simple demonstration and validation of the generic periodic boundary operator.
7+
8+
We construct an unstructured periodic 3D grid and solve a simple linear elastic problem
9+
with periodic coupling along the x-axis.
10+
11+
![](example312.png)
12+
=#
13+
14+
module Example312_PeriodicBoundary3D
15+
16+
using ExtendableFEM
17+
using ExtendableGrids
18+
using SimplexGridFactory
19+
using GridVisualize
20+
using TetGen
21+
using UnicodePlots
22+
using StaticArrays
23+
using LinearAlgebra
24+
using Test #hide
25+
26+
const reg_left = 1
27+
const reg_right = 2
28+
const reg_dirichlet = 3
29+
const reg_default = 4
30+
31+
32+
# define the Hooke tensor for AlN (from DOI 10.1063/1.1368156)
33+
function material_tensor()
34+
c11 = 396.0
35+
c12 = 137.0
36+
c13 = 108.0
37+
c33 = 373.0
38+
c44 = 116.0
39+
40+
return @SArray [
41+
c11 c12 c13 0 0 0
42+
c12 c11 c13 0 0 0
43+
c13 c13 c33 0 0 0
44+
0 0 0 c44 0 0
45+
0 0 0 0 c44 0
46+
0 0 0 0 0 c44
47+
]
48+
end
49+
50+
## generate the kernels for the linear problem
51+
## 𝐂: Hooke tensor, 𝑓: body force
52+
function make_kernels(𝐂, 𝑓)
53+
54+
## linear stress-strain mapping
55+
bilinear_kernel!(σ, εv, qpinfo) = mul!(σ, 𝐂, εv)
56+
57+
## body force
58+
linear_kernel!(result, qpinfo) = (result .= 𝑓)
59+
60+
return bilinear_kernel!, linear_kernel!
61+
end
62+
63+
64+
"""
65+
create 3D grid with Dirichlet boundary region at the bottom center
66+
"""
67+
function create_grid(; h, height, width, depth)
68+
builder = SimplexGridBuilder(; Generator = TetGen)
69+
70+
## bottom points
71+
b01 = point!(builder, 0, 0, 0)
72+
b02 = point!(builder, 0.45 * width, 0, 0)
73+
b03 = point!(builder, 0.55 * width, 0, 0)
74+
b04 = point!(builder, width, 0, 0)
75+
76+
b11 = point!(builder, 0, depth, 0)
77+
b12 = point!(builder, 0.45 * width, depth, 0)
78+
b13 = point!(builder, 0.55 * width, depth, 0)
79+
b14 = point!(builder, width, depth, 0)
80+
81+
## top points
82+
t01 = point!(builder, 0, 0, height)
83+
t02 = point!(builder, width, 0, height)
84+
85+
t11 = point!(builder, 0, depth, height)
86+
t12 = point!(builder, width, depth, height)
87+
88+
## center points
89+
c01 = point!(builder, 0.5 * width, 0, 0)
90+
c02 = point!(builder, 0.5 * width, 0, height)
91+
c11 = point!(builder, 0.5 * width, depth, 0)
92+
c12 = point!(builder, 0.5 * width, depth, height)
93+
94+
## default faces
95+
facetregion!(builder, reg_default)
96+
facet!(builder, b01, b02, b12, b11)
97+
facet!(builder, b03, b04, b14, b13)
98+
facet!(builder, [t01, c02, c12, t11])
99+
facet!(builder, [c02, t02, t12, c12])
100+
facet!(builder, [b01, b02, c01, c02, t01])
101+
facet!(builder, [c01, b03, b04, t02, c02])
102+
facet!(builder, [b11, b12, c11, c12, t11])
103+
facet!(builder, [c11, b13, b14, t12, c12])
104+
facet!(builder, c01, c02, c12, c11)
105+
106+
## left face
107+
facetregion!(builder, reg_left)
108+
facet!(builder, b01, t01, t11, b11)
109+
110+
## right face
111+
facetregion!(builder, reg_right)
112+
facet!(builder, b04, t02, t12, b14)
113+
114+
## Dirichlet face
115+
facetregion!(builder, reg_dirichlet)
116+
facet!(builder, [b02, c01, c11, b12])
117+
facet!(builder, [c01, b03, b13, c11])
118+
119+
cellregion!(builder, 1)
120+
maxvolume!(builder, h)
121+
regionpoint!(builder, width / 3, depth / 2, height / 2)
122+
cellregion!(builder, 2)
123+
## finer grid on the right half to make the periodic coupling non-trivial
124+
maxvolume!(builder, 0.3 * h)
125+
regionpoint!(builder, 2 * width / 3, depth / 2, height / 2)
126+
127+
return simplexgrid(builder)
128+
end
129+
130+
function main(;
131+
order = 1,
132+
periodic = true,
133+
Plotter = nothing,
134+
force = 1.0,
135+
h = 1.0e-4,
136+
width = 6.0,
137+
height = 0.2,
138+
depth = 1,
139+
kwargs...
140+
)
141+
142+
xgrid = create_grid(; h, width, height, depth)
143+
144+
## create finite element space and solution vector
145+
if order == 1
146+
FES = FESpace{H1P1{3}}(xgrid)
147+
elseif order == 2
148+
FES = FESpace{H1P2{3, 3}}(xgrid)
149+
end
150+
151+
## problem description
152+
PD = ProblemDescription()
153+
u = Unknown("u"; name = "displacement")
154+
assign_unknown!(PD, u)
155+
156+
𝐂 = material_tensor()
157+
𝑓 = force * [0, 0, 1]
158+
159+
bilinear_kernel!, linear_kernel! = make_kernels(𝐂, 𝑓)
160+
assign_operator!(PD, BilinearOperator(bilinear_kernel!, [εV(u, 1.0)]; kwargs...))
161+
assign_operator!(PD, LinearOperator(linear_kernel!, [id(u)]; kwargs...))
162+
163+
assign_operator!(PD, HomogeneousBoundaryData(u; regions = [reg_dirichlet]))
164+
165+
if periodic
166+
function give_opposite!(y, x)
167+
y .= x
168+
y[1] = width - x[1]
169+
return nothing
170+
end
171+
coupling_matrix = get_periodic_coupling_matrix(FES, xgrid, reg_left, reg_right, give_opposite!)
172+
display(coupling_matrix)
173+
assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
174+
end
175+
176+
## solve
177+
sol = solve(PD, FES; kwargs...)
178+
179+
plt = GridVisualizer(; Plotter, size = (1200, 800))
180+
magnification = 1
181+
displaced_grid = deepcopy(xgrid)
182+
displace_mesh!(displaced_grid, sol[1], magnify = magnification)
183+
gridplot!(plt, displaced_grid, linewidth = 1, title = "displaced mesh, $(magnification)x magnified", scene3d = :LScene)
184+
185+
return sol, plt
186+
187+
end
188+
189+
generateplots = ExtendableFEM.default_generateplots(Example312_PeriodicBoundary3D, "example312.png") #hide
190+
function runtests() #hide
191+
sol, plt = main() #hide
192+
@test maximum(view(sol[1])) ≈ 1.8038107003663026 #hide
193+
return nothing #hide
194+
end #hide
195+
196+
end # module

0 commit comments

Comments
 (0)