Skip to content

Commit 9c35ef4

Browse files
PetrKryslUCSDCopilot
PetrKryslUCSD
and
Copilot
committed
clean up
Co-authored-by: Copilot <copilot@github.com>
1 parent d56fa9e commit 9c35ef4

1 file changed

Lines changed: 68 additions & 147 deletions

File tree

‎examples/shells/statics/homogeneous/plates/simply_supp_square_plate_udl_examples.jl‎

Lines changed: 68 additions & 147 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,16 @@ using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_fie
1515
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
1616
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
1717

18+
const E = 30e6;
19+
const nu = 0.3;
20+
const L = 10.0;
21+
22+
loading(t) = 1.0e3*t^3
23+
24+
# From A Triangular Plate Bending Element Based on an Energy-Orthogonal Free,
25+
# by Felippa and Bergan, 1986, Table 3.
26+
w_analyt_sol(p, L, D) = -4.06235e-3 * p * L^4 / D
27+
1828
function _cond(A)
1929
emax = eigs(A, nev = 1, which = :LM, tol = 1e-3)[1]
2030
emin = eigs(A, nev = 1, which = :SM, tol = 1e-3)[1]
@@ -30,17 +40,10 @@ end
3040

3141
function _execute_t3ff_quarter_model(n = 2, tL_ratio = 0.01, visualize = true)
3242
formul = FEMMShellT3FFModule
33-
E = 30e6;
34-
nu = 0.3;
35-
L = 10.0;
3643
thickness = L * tL_ratio;
3744
D = E*thickness^3/12/(1-nu^2)
38-
p = 1.0*tL_ratio
39-
# analytical solution for the vertical deflection under the load
40-
# analytical solution for the vertical deflection under the load
41-
# From A Triangular Plate Bending Element Based on an Energy-Orthogonal Free,
42-
# by Felippa and Bergan, 1986, Table 3.
43-
analyt_sol=-4.06235e-3*p*L^4/D;
45+
p = loading(thickness)
46+
analyt_sol = w_analyt_sol(p, L, D)
4447

4548
tolerance = L/n/1000
4649
fens, fes = T3block(L/2,L/2,n,n);
@@ -110,28 +113,10 @@ function _execute_t3ff_quarter_model(n = 2, tL_ratio = 0.01, visualize = true)
110113
# Solve
111114
solve_blocked!(dchi, K, F)
112115
targetu = dchi.values[nl, 3][1]
113-
@info "Target: $(round(targetu, digits=8)), $(round(targetu/analyt_sol, digits = 4)*100)%"
116+
@info "w(center): $(round(targetu, digits=8)), $(round(targetu/analyt_sol, digits = 4)*100)%"
114117

115118
# Visualization
116119
if visualize
117-
118-
# Generate a graphical display of resultants
119-
# scalars = []
120-
# for nc in 1:3
121-
# fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
122-
# push!(scalars, ("m$nc", fld.values))
123-
# fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
124-
# push!(scalars, ("em$nc", fld.values))
125-
# end
126-
# vtkwrite("scordelis_lo_examples-$(n)-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
127-
# scalars = []
128-
# for nc in 1:3
129-
# fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
130-
# push!(scalars, ("n$nc", fld.values))
131-
# fld = elemfieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
132-
# push!(scalars, ("en$nc", fld.values))
133-
# end
134-
# vtkwrite("scordelis_lo_examples-$(n)-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
135120
scalars = []
136121
for nc in 1:2
137122
fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
@@ -154,13 +139,9 @@ end
154139

155140
function _execute_Q4RS_quarter_model(n = 2, tL_ratio = 0.01, visualize = true)
156141
formul = FEMMShellQ4RSModule
157-
E = 30e6;
158-
nu = 0.3;
159-
L = 10.0;
160-
thickness = L * tL_ratio;
142+
thickness = L * tL_ratio
161143
D = E*thickness^3/12/(1-nu^2)
162-
p = 1.0*tL_ratio
163-
# analytical solution for the vertical deflection under the load
144+
p = loading(thickness)
164145
# analytical solution for the vertical deflection under the load
165146
# From A Triangular Plate Bending Element Based on an Energy-Orthogonal Free,
166147
# by Felippa and Bergan, 1986, Table 3.
@@ -234,7 +215,7 @@ function _execute_Q4RS_quarter_model(n = 2, tL_ratio = 0.01, visualize = true)
234215
# Solve
235216
solve_blocked!(dchi, K, F)
236217
targetu = dchi.values[nl, 3][1]
237-
@info "Target: $(round(targetu, digits=8)), $(round(targetu/analyt_sol, digits = 4)*100)%"
218+
@info "w(center): $(round(targetu, digits=8)), $(round(targetu/analyt_sol, digits = 4)*100)%"
238219

239220
# Visualization
240221
if visualize
@@ -296,24 +277,31 @@ function test_convergence_quarter()
296277
return true
297278
end
298279

299-
function _execute_t3ff_full_model(n = 2, tL_ratio = 0.01, simple_support = :hard, stab_alpha = 0.1, visualize = true)
300-
E = 30e6;
301-
nu = 0.3;
302-
L = 10.0;
280+
function _execute_t3ff_full_model(
281+
n=2,
282+
tL_ratio=0.01,
283+
simple_support=:hard,
284+
stab_alpha=0.1,
285+
tweak_mesh=false,
286+
visualize=true
287+
)
288+
formul = FEMMShellT3FFModule
303289
thickness = L * tL_ratio;
304290
D = E/12/(1-nu^2)*thickness^3
305-
p = 1.0*thickness^3
306-
# analytical solution for the vertical deflection under the load
307-
# analytical solution for the vertical deflection under the load
308-
# From A Triangular Plate Bending Element Based on an Energy-Orthogonal Free,
309-
# by Felippa and Bergan, 1986, Table 3.
310-
analyt_sol=-4.06235e-3*p*L^4/D;
311-
291+
p = loading(thickness)
292+
analyt_sol = w_analyt_sol(p, L, D)
293+
312294
tolerance = L/n/1000
313295
fens, fes = T3block(L,L,n,n);
314296
fens.xyz = xyz3(fens)
315297
fens.xyz[:, 1] .-= L/2
316298
fens.xyz[:, 2] .-= L/2
299+
nm = selectnode(fens; box = Float64[-L/2+L/n -L/2+L/n -L/2+L/n -L/2+L/n -Inf Inf], inflate = tolerance)
300+
if tweak_mesh
301+
shift = L/n/4
302+
fens.xyz[nm, 1] .+= shift
303+
fens.xyz[nm, 2] .+= shift
304+
end
317305

318306
mater = MatDeforElastIso(DeforModelRed3D, E, nu)
319307

@@ -377,68 +365,27 @@ function _execute_t3ff_full_model(n = 2, tL_ratio = 0.01, simple_support = :hard
377365
# Solve
378366
solve_blocked!(dchi, K, F)
379367
targetu = dchi.values[nl, 3][1]
380-
@info "Target: $(round(targetu, digits=8)), $(round(targetu/analyt_sol, digits = 4)*100)%"
368+
@info "w(center): $(round(targetu, digits=8)), $(round(targetu/analyt_sol, digits = 4)*100)%"
381369
epsrel = abs(targetu - analyt_sol) / abs(analyt_sol)
382370
@info "Digits of accuracy: $(-log10(epsrel))"
383371
# Visualization
384372
if visualize
385-
U = gathersysvec(dchi)
386-
# Generate a graphical display of resultants
387-
# scalars = []
388-
# for nc in 1:3
389-
# fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
390-
# push!(scalars, ("m$nc", fld.values))
391-
# fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
392-
# push!(scalars, ("em$nc", fld.values))
393-
# end
394-
# vtkwrite("scordelis_lo_examples-$(n)-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
395-
# scalars = []
396-
# for nc in 1:3
397-
# fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
398-
# push!(scalars, ("n$nc", fld.values))
399-
# fld = elemfieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
400-
# push!(scalars, ("en$nc", fld.values))
401-
# end
402-
# vtkwrite("scordelis_lo_examples-$(n)-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
403-
# scalars = []
404-
# for nc in 1:2
405-
# fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
406-
# push!(scalars, ("q$nc", fld.values))
407-
# fld = elemfieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
408-
# push!(scalars, ("eq$nc", fld.values))
409-
# end
410-
# vtkwrite("simply_supp_square_plate_udl-full-n=$(n)-q.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
411-
412-
ocsys = CSys(3)
373+
vtkwrite("ss_sqpl_udl-t3ff-full-$(simple_support)-r=$(tL_ratio)-s=$(stab_alpha)-n=$(n)-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3]), ("ur", dchi.values[:, 4:6])])
374+
# ocsys = CSys(3)
413375
scalars = []
414376
for nc in 1:3
415377
fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys=ocsys)
416378
# fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
417379
push!(scalars, ("m$nc", fld.values))
418380
end
419-
vtkwrite("simply_supp_square_plate_udl-r=$(tL_ratio)-full-$(simple_support)-n=$(n)-m.vtu", fens, fes; scalars=scalars, vectors=[("u", dchi.values[:, 1:3])])
420-
scalars = []
421-
for nc in 1:3
422-
fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys=ocsys)
423-
# fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
424-
push!(scalars, ("n$nc", fld.values))
425-
end
426-
vtkwrite("simply_supp_square_plate_udl-r=$(tL_ratio)-full-$(simple_support)-n=$(n)-n.vtu", fens, fes; scalars=scalars, vectors=[("u", dchi.values[:, 1:3])])
381+
vtkwrite("ss_sqpl_udl-t3ff-full-$(simple_support)-r=$(tL_ratio)-s=$(stab_alpha)-m.vtu", fens, fes; scalars=scalars, vectors=[("u", dchi.values[:, 1:3])])
427382
scalars = []
428383
for nc in 1:2
429384
fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys=ocsys)
430385
# fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
431386
push!(scalars, ("q$nc", fld.values))
432387
end
433-
vtkwrite("simply_supp_square_plate_udl-r=$(tL_ratio)-full-$(simple_support)-n=$(n)-q.vtu", fens, fes; scalars=scalars, vectors=[("u", dchi.values[:, 1:3])])
434-
435-
# scattersysvec!(dchi, (L/4)/maximum(abs.(U)).*U)
436-
# update_rotation_field!(Rfield0, dchi)
437-
# plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
438-
# #plot_nodes(fens),
439-
# plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
440-
# dims = 1)
441-
# pl = render(plots)
388+
vtkwrite("ss_sqpl_udl-t3ff-full-$(simple_support)-r=$(tL_ratio)-s=$(stab_alpha)-n=$(n)-q.vtu", fens, fes; scalars=scalars, vectors=[("u", dchi.values[:, 1:3])])
442389
end
443390
return true
444391
end
@@ -452,17 +399,10 @@ function _execute_q4rs_full_model(
452399
visualize=true
453400
)
454401
formul = FEMMShellQ4RSModule
455-
E = 30e6;
456-
nu = 0.3;
457-
L = 10.0;
458402
thickness = L * tL_ratio;
459403
D = E/12/(1-nu^2)*thickness^3
460-
p = 1.0*thickness^3
461-
# analytical solution for the vertical deflection under the load
462-
# analytical solution for the vertical deflection under the load
463-
# From A Triangular Plate Bending Element Based on an Energy-Orthogonal Free,
464-
# by Felippa and Bergan, 1986, Table 3.
465-
analyt_sol=-4.06235e-3*p*L^4/D;
404+
p = loading(thickness)
405+
analyt_sol = w_analyt_sol(p, L, D)
466406

467407
tolerance = L/n/1000
468408
fens, fes = Q4block(L,L,n,n);
@@ -540,40 +480,12 @@ function _execute_q4rs_full_model(
540480
# Solve
541481
solve_blocked!(dchi, K, F)
542482
targetu = dchi.values[nl, 3][1]
543-
@info "Target: $(round(targetu, digits=8)), $(round(targetu/analyt_sol, digits = 4)*100)%"
483+
@info "w(center): $(round(targetu, digits=8)), $(round(targetu/analyt_sol, digits = 4)*100)%"
544484
epsrel = abs(targetu - analyt_sol) / abs(analyt_sol)
545485
@info "Digits of accuracy: $(-log10(epsrel))"
546486
# Visualization
547487
if visualize
548-
U = gathersysvec(dchi)
549-
# Generate a graphical display of resultants
550-
# scalars = []
551-
# for nc in 1:3
552-
# fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
553-
# push!(scalars, ("m$nc", fld.values))
554-
# fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
555-
# push!(scalars, ("em$nc", fld.values))
556-
# end
557-
# vtkwrite("scordelis_lo_examples-$(n)-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
558-
# scalars = []
559-
# for nc in 1:3
560-
# fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
561-
# push!(scalars, ("n$nc", fld.values))
562-
# fld = elemfieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
563-
# push!(scalars, ("en$nc", fld.values))
564-
# end
565-
# vtkwrite("scordelis_lo_examples-$(n)-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
566-
# scalars = []
567-
# for nc in 1:2
568-
# fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
569-
# push!(scalars, ("q$nc", fld.values))
570-
# fld = elemfieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
571-
# push!(scalars, ("eq$nc", fld.values))
572-
# end
573-
# vtkwrite("simply_supp_square_plate_udl-full-n=$(n)-q.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
574-
575488
vtkwrite("ss_sqpl_udl-q4rs-full-$(simple_support)-r=$(tL_ratio)-s=$(stab_alpha)-n=$(n)-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3]), ("ur", dchi.values[:, 4:6])])
576-
577489
# ocsys = CSys(3)
578490
scalars = []
579491
for nc in 1:3
@@ -582,28 +494,13 @@ function _execute_q4rs_full_model(
582494
push!(scalars, ("m$nc", fld.values))
583495
end
584496
vtkwrite("ss_sqpl_udl-q4rs-full-$(simple_support)-r=$(tL_ratio)-s=$(stab_alpha)-m.vtu", fens, fes; scalars=scalars, vectors=[("u", dchi.values[:, 1:3])])
585-
# scalars = []
586-
# for nc in 1:3
587-
# fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys=ocsys)
588-
# # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
589-
# push!(scalars, ("n$nc", fld.values))
590-
# end
591-
# vtkwrite("simply_supp_square_plate_udl-r=$(tL_ratio)-full-$(simple_support)-n=$(n)-n.vtu", fens, fes; scalars=scalars, vectors=[("u", dchi.values[:, 1:3])])
592497
scalars = []
593498
for nc in 1:2
594499
fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys=ocsys)
595500
# fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
596501
push!(scalars, ("q$nc", fld.values))
597502
end
598503
vtkwrite("ss_sqpl_udl-q4rs-full-$(simple_support)-r=$(tL_ratio)-s=$(stab_alpha)-n=$(n)-q.vtu", fens, fes; scalars=scalars, vectors=[("u", dchi.values[:, 1:3])])
599-
600-
# scattersysvec!(dchi, (L/4)/maximum(abs.(U)).*U)
601-
# update_rotation_field!(Rfield0, dchi)
602-
# plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
603-
# #plot_nodes(fens),
604-
# plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
605-
# dims = 1)
606-
# pl = render(plots)
607504
end
608505
return true
609506
end
@@ -663,7 +560,7 @@ function test_convergence_full_thickness(tL_ratios = [1.0e-1, 1.0e-2, 1.0e-4, 1.
663560
return true
664561
end
665562

666-
function test_stabilization(tL_ratios = [1.0e-2, 1.0e-4, 1.0e-8])
563+
function test_stabilization_q4rs(tL_ratios = [1.0e-2, 1.0e-4, 1.0e-8])
667564
stab_alpha = 0.001
668565
visualize = true
669566
support = :hard
@@ -681,13 +578,37 @@ function test_stabilization(tL_ratios = [1.0e-2, 1.0e-4, 1.0e-8])
681578
return true
682579
end
683580

581+
function test_stabilization_t3ff(tL_ratios = [1.0e-2, 1.0e-4, 1.0e-8])
582+
stab_alpha = 0.001
583+
visualize = true
584+
support = :hard
585+
tweak_mesh = true
586+
n = 16
587+
formul = FEMMShellT3FFModule
588+
@info "Simply supported square plate with uniform load, T3FF --------------------"
589+
for tL_ratio in tL_ratios
590+
@info ">>>>>>>>>>>>> thickness/length = $tL_ratio"
591+
for stab_alpha in [0.0, 0.000001, 0.001, 0.1]
592+
@info " stab_alpha = $stab_alpha"
593+
_execute_t3ff_full_model(n, tL_ratio, support, stab_alpha, tweak_mesh, visualize)
594+
end
595+
end
596+
return true
597+
end
598+
684599
function allrun()
600+
println("#####################################################")
601+
println("# test_stabilization_t3ff ")
602+
test_stabilization_t3ff()
603+
println("#####################################################")
604+
println("# test_stabilization_q4rs ")
605+
test_stabilization_q4rs()
685606
println("#####################################################")
686607
println("# test_convergence_full ")
687608
# test_convergence_full()
688609
println("#####################################################")
689610
println("# test_convergence_full_thickness ")
690-
test_convergence_full_thickness()
611+
# test_convergence_full_thickness()
691612
return true
692613
end # function allrun
693614

0 commit comments

Comments
 (0)