Skip to content

Commit 82decb0

Browse files
pjaappjaap
committed
Another 2D example
1 parent c7d5d4e commit 82decb0

2 files changed

Lines changed: 285 additions & 20 deletions

File tree

‎examples/Example212_PeriodicBoundary.jl‎

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

‎src/common_operators/combinedofs.jl‎

Lines changed: 83 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -84,8 +84,10 @@ function build_assembler!(CD::CombineDofs{UT, CT}, FE::Array{<:FEVectorBlock, 1}
8484
if CD.parameters[:verbosity] > 0
8585
@info ".... coupling $(length(coupling_matrix.nzval)) dofs"
8686
end
87-
function assemble(A::AbstractSparseArray{T}, b::AbstractVector{T}, assemble_matrix::Bool, assemble_rhs::Bool, kwargs...) where {T}
87+
function assemble!(A::AbstractSparseArray{T}, b::AbstractVector{T}, assemble_matrix::Bool, assemble_rhs::Bool, kwargs...) where {T}
8888
if assemble_matrix
89+
90+
# put FE adjacency information into this matrix row
8991
for dof_i in 1:size(coupling_matrix, 2)
9092
# this col-view is efficient
9193
coupling_i = @views coupling_matrix[:, dof_i]
@@ -94,23 +96,69 @@ function build_assembler!(CD::CombineDofs{UT, CT}, FE::Array{<:FEVectorBlock, 1}
9496
continue
9597
end
9698

97-
# get the coupled dofs of dof_i and the corresponding weights
98-
coupled_dofs, weights = findnz(coupling_i)
99+
targetrow = dof_i + offsetX
100+
101+
coupled_dofs_i, weights_i = findnz(coupling_i)
102+
103+
for dof_j in 1:size(coupling_matrix, 2)
104+
# this col-view is efficient
105+
coupling_j = @views coupling_matrix[:, dof_j]
106+
107+
if nnz(coupling_j) == 0
108+
# not both i and j on the boundary
109+
110+
targetcol = dof_j + offsetY
111+
for (dof_k, weight_k) in zip(coupled_dofs_i, weights_i)
112+
sourcerow = dof_k + offsetX
113+
sourcecol = targetcol
114+
val = A[sourcerow, sourcecol]
115+
_addnz(A, targetrow, targetcol, val, weight_k)
116+
end
117+
118+
119+
else
120+
# i and j on the same boundary
121+
122+
coupled_dofs_j, weights_j = findnz(coupling_j)
123+
targetcol = dof_j + offsetY
99124

100-
# transfer all assembly information of dof_i to the
101-
# coupled dofs: add the corresponding matrix row
102-
sourcerow = dof_i + offsetY
103-
for sourcecol in 1:size(A, 2)
104-
val = A[sourcerow, sourcecol]
105-
for dof_j in coupled_dofs
106-
targetrow = dof_j + offsetX
107-
targetcol = sourcecol
108-
_addnz(A, targetrow, targetcol, val, 1)
125+
for (dof_k, weight_k) in zip(coupled_dofs_i, weights_i)
126+
for (dof_l, weight_l) in zip(coupled_dofs_j, weights_j)
127+
128+
sourcerow = dof_k + offsetX
129+
sourcecol = dof_l + offsetY
130+
131+
val = A[sourcerow, sourcecol]
132+
# @show sourcerow, sourcecol, targetrow, targetcol, val, weight_k, weight_l
133+
_addnz(A, targetrow, targetcol, val, weight_k * weight_l)
134+
135+
end
136+
end
109137
end
110-
# eliminate the sourcerow
111-
A[sourcerow, sourcecol] = 0
112138
end
139+
end
140+
141+
# replace the geometric coupling rows
142+
for dof_i in 1:size(coupling_matrix, 2)
143+
# go through the rows of the matrix; use transposed matrix implicitly
144+
145+
# this row-view is not efficient :(
146+
coupling_i = coupling_matrix[dof_i, :]
147+
# do nothing if no coupling for dof_i
148+
if nnz(coupling_i) == 0
149+
continue
150+
end
151+
152+
153+
# get the coupled dofs of dof_i and the corresponding weights
154+
coupled_dofs, weights = findnz(coupling_i)
113155

156+
sourcerow = dof_i + offsetX
157+
158+
# eliminate the sourcerow
159+
for col in 1:size(A, 2)
160+
A[sourcerow, col] = 0
161+
end
114162

115163
# replace sourcerow with coupling linear combination
116164
_addnz(A, sourcerow, sourcerow, -1.0, 1)
@@ -119,11 +167,12 @@ function build_assembler!(CD::CombineDofs{UT, CT}, FE::Array{<:FEVectorBlock, 1}
119167
_addnz(A, sourcerow, dof_j + offsetY, weight, 1)
120168
end
121169

122-
flush!(A)
123170
end
171+
flush!(A)
124172
end
125173

126174
if assemble_rhs
175+
127176
for dof_i in 1:size(coupling_matrix, 2)
128177
# this col-view is efficient
129178
coupling_i = @views coupling_matrix[:, dof_i]
@@ -132,22 +181,36 @@ function build_assembler!(CD::CombineDofs{UT, CT}, FE::Array{<:FEVectorBlock, 1}
132181
continue
133182
end
134183

184+
# get the coupled dofs of dof_i and the corresponding weights
135185
coupled_dofs, weights = findnz(coupling_i)
136-
sourcerow = dof_i + offsetY
137186

138-
for dof_j in coupled_dofs
139-
targetrow = dof_j + offsetX
140-
b[targetrow] += b[sourcerow]
187+
# transfer all assembly information to dof_i
188+
targetrow = dof_i + offsetY
189+
for (dof_j, weight) in zip(coupled_dofs, weights)
190+
sourcerow = dof_j + offsetY
191+
b[targetrow] += weight * b[sourcerow]
192+
end
193+
end
194+
195+
196+
for dof_i in 1:size(coupling_matrix, 2)
197+
coupling_i = coupling_matrix[dof_i, :]
198+
# do nothing if no coupling for dof_i
199+
if nnz(coupling_i) == 0
200+
continue
141201
end
142202

203+
sourcerow = dof_i + offsetX
204+
143205
b[sourcerow] = 0.0
206+
144207
end
145208
end
146209

147210
return nothing
148211
end
149212

150-
CD.assembler = assemble
213+
CD.assembler = assemble!
151214
CD.FESX = FESX
152215
CD.FESY = FESY
153216
end

0 commit comments

Comments
 (0)