/
PoissonAgFEM.jl
111 lines (81 loc) · 2.37 KB
/
PoissonAgFEM.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
module PoisonAgFEMTests
using IterativeSolvers: cg
using Preconditioners: AMGPreconditioner, SmoothedAggregation
using STLCutters
using Gridap
using GridapEmbedded
using Test
using STLCutters: compute_stl_model
using STLCutters: read_stl, merge_nodes, get_bounding_box
# Manufactured solution
u(x) = x[1] + x[2] - x[3]
f(x) = - Δ(u)(x)
ud(x) = u(x)
stlpath = joinpath(@__DIR__,"data/Bunny-LowPoly.stl")
@time geo = STLGeometry( stlpath )
n = 28
δ = 0.2
pmin,pmax = get_bounding_box(geo)
diagonal = pmax-pmin
pmin = pmin - diagonal*δ
pmax = pmax + diagonal*δ
partition = (n,n,n)
bgmodel = CartesianDiscreteModel(pmin,pmax,partition)
# Cut the background model
cutter = STLCutter()
@time cutgeo,facet_to_inoutcut = cut(cutter,bgmodel,geo)
strategy = AggregateAllCutCells()
aggregates = aggregate(strategy,cutgeo,facet_to_inoutcut)
# Setup integration meshes
Ω_bg = Triangulation(bgmodel)
Ω = Triangulation(cutgeo)
Γd = EmbeddedBoundary(cutgeo)
# Setup normal vectors
n_Γd = get_normal_vector(Γd)
#writevtk(Ω,"trian_O")
#writevtk(Γd,"trian_Gd",cellfields=["normal"=>n_Γd])
#writevtk(Triangulation(bgmodel),"bgtrian")
# Setup Lebesgue measures
order = 1
degree = 2*order
dΩ = Measure(Ω,degree)
dΓd = Measure(Γd,degree)
#dΓg = Measure(Γg,degree)
vol = sum( ∫(1)*dΩ )
surf = sum( ∫(1)*dΓd )
@show vol,surf
# Setup FESpace
model = DiscreteModel(cutgeo)
Vstd = FESpace(model,ReferenceFE(lagrangian,Float64,order),conformity=:H1)
@time V = AgFEMSpace(Vstd,aggregates) #,Vser)
U = TrialFESpace(V)
# Weak form
γd = 10.0
h = (pmax - pmin)[1] / partition[1]
a(u,v) =
∫( ∇(v)⋅∇(u) ) * dΩ +
∫( (γd/h)*v*u - v*(n_Γd⋅∇(u)) - (n_Γd⋅∇(v))*u ) * dΓd
l(v) =
∫( v*f ) * dΩ +
∫( (γd/h)*v*ud - (n_Γd⋅∇(v))*ud ) * dΓd
# FE problem
@time op = AffineFEOperator(a,l,U,V)
A = get_matrix(op)
b = get_vector(op)
@time p = AMGPreconditioner{SmoothedAggregation}(A)
@time x = cg(A,b,verbose=true,Pl=p,reltol=1e-10)
@time uh = FEFunction(U,x)
@time e = u - uh
# Postprocess
l2(u) = sqrt(sum( ∫( u*u )*dΩ ))
h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ ))
el2 = l2(e)
eh1 = h1(e)
ul2 = l2(uh)
uh1 = h1(uh)
#colors = color_aggregates(aggregates,bgmodel)
#writevtk(Ω_bg,"trian",celldata=["aggregate"=>aggregates,"color"=>colors],cellfields=["uh"=>uh])
#writevtk(Ω,"results",cellfields=["uh"=>uh])
@test el2/ul2 < 1.e-9
@test eh1/uh1 < 1.e-9
end # module