Skip to content

Commit 5aee3eb

Browse files
jmartfrutjmartfrut
committed
new DirichletBC tests
1 parent 9119a4e commit 5aee3eb

File tree

4 files changed

+34501
-0
lines changed

4 files changed

+34501
-0
lines changed

‎test/SimulationsTests/BoundaryConditions.jl‎

Lines changed: 144 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,144 @@
1+
using Gridap
2+
using Gridap.FESpaces
3+
4+
@testset "Dirichlet BC testing analytical mapping" begin
5+
6+
meshfile = "test_BC1.msh"
7+
geomodel = GmshDiscreteModel("./test/models/" * meshfile)
8+
9+
# Domains
10+
order = 1
11+
degree = 1 * order
12+
Ω = Interior(geomodel, tags=["Air"])
13+
dΩ = Measure(Ω, degree)
14+
15+
function Mapping(x, Λ)
16+
L = 0.1
17+
A = 0.98 * L
18+
xd = [x[1], x[2] + A * sin(2 * π * (x[1] + L / 2) / L)]
19+
u = (xd - [x[1], x[2]]) * Λ
20+
return VectorValue(u)
21+
end
22+
23+
evolu(Λ) = Λ
24+
dir_u_tags_air = ["uair_fixed", "Interface"]
25+
dir_u_values_air = [[0.0, 0.0], Λ -> (x -> Mapping(x, Λ))]
26+
dir_u_timesteps_air = [evolu, nothing]
27+
Du = DirichletBC(dir_u_tags_air, dir_u_values_air, dir_u_timesteps_air)
28+
29+
reffeu = ReferenceFE(lagrangian, VectorValue{2,Float64}, order)
30+
31+
Vu = TestFESpace(Ω, reffeu, Du, conformity=:H1)
32+
Uu = TrialFESpace(Vu, Du, 1.0)
33+
34+
@test norm(Uu.dirichlet_values) == 0.30990321069650995
35+
36+
TrialFESpace!(Uu, Du, 0.4)
37+
@test norm(Uu.dirichlet_values) == 0.12396128427860398
38+
TrialFESpace!(Uu, Du, 1.0)
39+
@test norm(Uu.dirichlet_values) == 0.30990321069650995
40+
41+
end
42+
43+
44+
@testset "Mesh movement stabilization" begin
45+
46+
meshfile = "test_BC2.msh"
47+
geomodel = GmshDiscreteModel("./test/models/" * meshfile)
48+
49+
μParams = [6456.9137547089595, 896.4633794151492,
50+
1.999999451256222,
51+
1.9999960497608036,
52+
11747.646562400318,
53+
0.7841068624959612, 1.5386288924587603]
54+
55+
model_vacuum_mech_ = NonlinearMooneyRivlin2D_CV(λ=1 * μParams[1], μ1=μParams[1], μ2=0.0, α1=6.0, α2=1.0, γ=6.0)
56+
model_vacuum_mech = HessianRegularization(mechano=model_vacuum_mech_, δ=1e-6 * μParams[1])
57+
58+
59+
# Domains
60+
order = 1
61+
degree = 1 * order
62+
bdegree = 1 * order
63+
Ωair = Interior(geomodel, tags=["Air"])
64+
dΩair = Measure(Ωair, degree)
65+
Ωsolid = Interior(geomodel, tags=["Solid"])
66+
67+
Γair_int = BoundaryTriangulation(Ωair, tags="Interface")
68+
nair_int = get_normal_vector(Γair_int)
69+
dΓair_int = Measure(Γair_int, bdegree)
70+
71+
Γsf = InterfaceTriangulation(Ωsolid, Ωair)
72+
nΓsf = get_normal_vector(Γsf)
73+
dΓsf = Measure(Γsf, bdegree)
74+
75+
76+
L = 0.1
77+
function Mapping(x, Λ)
78+
θmax = -1.0 * π / 2
79+
#θmax = -2.0*π/2*0.001
80+
A = 0.3 * L * Λ
81+
# xd = [x[1], x[2]+A*((x[1]+L/2)/L)^2]
82+
xd = [x[1], x[2] + A * sin(π * (x[1] + L / 2) / L)]
83+
84+
θ = θmax * Λ
85+
R = [[cos(θ) -sin(θ)]; [sin(θ) cos(θ)]]
86+
xd2 = R * (xd + [L / 2, 0]) - [L / 2, 0.0]
87+
u = (xd2 - [x[1], x[2]])
88+
return VectorValue(u)
89+
end
90+
91+
92+
# FE spaces
93+
reffeu = ReferenceFE(lagrangian, VectorValue{2,Float64}, order)
94+
reffeJ = ReferenceFE(lagrangian, Float64, order)
95+
reffeJx = ReferenceFE(lagrangian, Float64, order - 1)
96+
97+
# Test FE Spaces
98+
Vu_⁺ = TestFESpace(Ωair, reffeu, dirichlet_tags=["uair_fixed", "Interface"], conformity=:H1)
99+
Vu_⁻ = TestFESpace(Ωair, reffeu, dirichlet_tags=["uair_fixed", "Interface"], conformity=:H1)
100+
101+
uh⁺ = interpolate_everywhere(x -> Mapping(x, 1.0), Vu_⁺)
102+
uh⁻ = interpolate_everywhere(x -> Mapping(x, 0.0), Vu_⁻)
103+
104+
dir(Λ) = uh⁻ + (uh⁺ - uh⁻) * Λ
105+
106+
evolu(Λ) = Λ
107+
dir_u_tags_air = ["uair_fixed", "Interface"]
108+
dir_u_values_air = [[0.0, 0.0], Λ -> dir(Λ)]
109+
dir_u_timesteps_air = [evolu, nothing]
110+
Du_air = DirichletBC(dir_u_tags_air, dir_u_values_air, dir_u_timesteps_air)
111+
112+
Vu = TestFESpace(Ωair, reffeu, Du_air, conformity=:H1)
113+
Uu = TrialFESpace(Vu, Du_air, 1.0)
114+
TrialFESpace!(Uu, Du_air, 0.0)
115+
116+
DΨvacuum_mech = model_vacuum_mech(1.0)
117+
k = Kinematics(Mechano, Solid)
118+
F, H, J = get_Kinematics(k; Λ=1.0)
119+
120+
# Vacuum mechanics
121+
res_vacmech(Λ) = (u, v) -> ∫((∇(v)' ⊙ (DΨvacuum_mech[2] ∘ (F ∘ (∇(u)')))))dΩair
122+
123+
jac_vacmech(Λ) = (u, du, v) -> ∫(∇(v)' ⊙ ((DΨvacuum_mech[3] ∘ (F ∘ (∇(u)'))) ⊙ (∇(du)')))dΩair
124+
125+
α = CellState(1.0, dΩair)
126+
linesearch = Injectivity_Preserving_LS(α, Uu, Vu; maxiter=50, αmin=1e-16, ρ=0.5, c=0.95)
127+
nls_vacmech = Newton_RaphsonSolver(LUSolver(); maxiter=10, rtol=2, verbose=false, linesearch=linesearch)
128+
129+
xh = FEFunction(Uu, zero_free_values(Uu))
130+
comp_model_vacmech = StaticNonlinearModel(res_vacmech, jac_vacmech, Uu, Vu, Du_air; nls=nls_vacmech, xh=xh)
131+
args_vacmech = Dict(:stepping => (nsteps=1, maxbisec=5), :ProjectDirichlet =>true)
132+
133+
nsteps = 10
134+
flagconv = 1 # convergence flag 0 (max bisections) 1 (max steps)
135+
Δβ = 1.0 / nsteps
136+
nbisect = 0
137+
for t in 0:nsteps-1
138+
interpolate_everywhere!(x -> Mapping(x, Δβ * (1 + t)), get_free_dof_values(uh⁺), uh⁺.dirichlet_values, Vu_⁺)
139+
interpolate_everywhere!(x -> Mapping(x, Δβ * (t)), get_free_dof_values(uh⁻), uh⁻.dirichlet_values, Vu_⁻)
140+
solve!(comp_model_vacmech; args_vacmech...)
141+
end
142+
@test norm(xh.free_values) == 2.8132015601158087
143+
144+
end

‎test/SimulationsTests/runtests.jl‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,4 +10,7 @@ using Test
1010

1111
include("StaticMechanicalNeumannTest.jl")
1212

13+
include("BoundaryConditions.jl")
14+
15+
1316
end

0 commit comments

Comments
 (0)