From db995aaebeff3b0fbc364f6ece2c3c97b8467dd5 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Mon, 9 Mar 2020 19:36:22 +0100 Subject: [PATCH 1/4] Added IsotropicDamageTest --- test/GridapTests/IsotropicDamageTests.jl | 134 +++++++++++++++++++++++ test/GridapTests/runtests.jl | 2 + 2 files changed, 136 insertions(+) create mode 100644 test/GridapTests/IsotropicDamageTests.jl diff --git a/test/GridapTests/IsotropicDamageTests.jl b/test/GridapTests/IsotropicDamageTests.jl new file mode 100644 index 000000000..2f562efb3 --- /dev/null +++ b/test/GridapTests/IsotropicDamageTests.jl @@ -0,0 +1,134 @@ +module IsotropicDamageTests + +using Gridap +using Gridap.FESpaces +using LinearAlgebra + +# max imposed displacement +const udx_max = 0.05 + +# Elastic model +const E = 2.1e4 # Pa +const ν = 0.3 # dim-less +const λ = (E*ν)/((1+ν)*(1-2*ν)) +const μ = E/(2*(1+ν)) +σe(ε) = λ*tr(ε)*one(ε) + 2*μ*ε # Pa +τ(ε) = sqrt(inner(ε,σe(ε))) # Pa^(1/2) + +# Damage model +const σ_u = 5.5 # Pa +const r_0 = σ_u / sqrt(E) # Pa^(1/2) +const H = 0.1 # dim-less + +function d(r) + 1 - q(r)/r +end + +function q(r) + r_0 + H*(r-r_0) +end + +function update(ε_in,r_in,d_in) + τ_in = τ(ε_in) + if τ_in <= r_in + r_out = r_in + d_out = d_in + damaged = false + else + r_out = τ_in + d_out = d(r_out) + damaged = true + end + damaged, r_out, d_out +end + +@law function σ(ε_in,r_in,d_in) + + _, _, d_out = update(ε_in,r_in,d_in) + + (1-d_out)*σe(ε_in) + +end + +@law function dσ(dε_in,ε_in,r_in,d_in) + + damaged, r_out, d_out = update(ε_in,r_in,d_in) + + if ! damaged + return (1-d_out)*σe(dε_in) + + else + c_inc = ((q(r_out) - H*r_out)*inner(σe(ε_in),dε_in))/(r_out^3) + return (1-d_out)*σe(dε_in) - c_inc*σe(ε_in) + + end +end + +function main(;n,nsteps) + + domain = (0,1,0,1,0,1) + partition = (n,n,n) + model = CartesianDiscreteModel(domain,partition) + + labeling = get_face_labeling(model) + add_tag_from_tags!(labeling,"ux0",[5,7,13,15,17,19,25]) + add_tag_from_tags!(labeling,"ux1",[2,4,6,8,14,16,18,20,26]) + add_tag_from_tags!(labeling,"uxyz0",[1]) + add_tag_from_tags!(labeling,"uxz0",[3]) + + order = 1 + + V = TestFESpace( + reffe=:Lagrangian, valuetype=VectorValue{3,Float64}, order=order, + model=model, + labels = labeling, + dirichlet_tags=["ux0","ux1","uxyz0","uxz0"], + dirichlet_masks=[(true,false,false),(true,false,false),(true,true,true),(true,false,true)]) + + trian = Triangulation(model) + degree = 2*order + quad = CellQuadrature(trian,degree) + + r = QPointCellField(r_0,trian,quad) + d = QPointCellField(0.0,trian,quad) + + nls = NLSolver(show_trace=true, method=:newton) + solver = FESolver(nls) + + function step(uh_in,factor) + + udx = factor*udx_max + u0 = VectorValue(0.0,0.0,0.0) + ud = VectorValue(udx,0.0,0.0) + U = TrialFESpace(V,[u0,ud,u0,u0]) + + res(u,v) = inner( ε(v), σ(ε(u),r,d) ) + jac(u,du,v) = inner( ε(v), dσ(ε(du),ε(u),r,d) ) + t_Ω = FETerm(res,jac,trian,quad) + op = FEOperator(U,V,t_Ω) + + free_values = get_free_values(uh_in) + uh0 = FEFunction(U,free_values) + + uh_out, = solve!(uh0,solver,op) + + update_state_variables!(quad,update,ε(uh_out),r,d) + + uh_out + end + + factors = collect(1:nsteps)*(1/nsteps) + uh = zero(V) + + for (istep,factor) in enumerate(factors) + + uh = step(uh,factor) + + end + +end + +main(n=5,nsteps=20) + + +end diff --git a/test/GridapTests/runtests.jl b/test/GridapTests/runtests.jl index 034e3fca0..4f6421deb 100644 --- a/test/GridapTests/runtests.jl +++ b/test/GridapTests/runtests.jl @@ -18,4 +18,6 @@ using Test @testset "SurfaceCoupling" begin include("SurfaceCouplingTests.jl") end +@testset "IsotropicDamage" begin include("IsotropicDamageTests.jl") end + end # module From 7f51f819f00edabcfb9013e4c9e16da2e71d5073 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Tue, 10 Mar 2020 08:42:14 +0100 Subject: [PATCH 2/4] Testing error norm in IsotropicDamageTest --- test/GridapTests/IsotropicDamageTests.jl | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/test/GridapTests/IsotropicDamageTests.jl b/test/GridapTests/IsotropicDamageTests.jl index 2f562efb3..a7f080230 100644 --- a/test/GridapTests/IsotropicDamageTests.jl +++ b/test/GridapTests/IsotropicDamageTests.jl @@ -3,10 +3,13 @@ module IsotropicDamageTests using Gridap using Gridap.FESpaces using LinearAlgebra +using Test # max imposed displacement const udx_max = 0.05 +const L = 1 + # Elastic model const E = 2.1e4 # Pa const ν = 0.3 # dim-less @@ -15,6 +18,7 @@ const μ = E/(2*(1+ν)) σe(ε) = λ*tr(ε)*one(ε) + 2*μ*ε # Pa τ(ε) = sqrt(inner(ε,σe(ε))) # Pa^(1/2) + # Damage model const σ_u = 5.5 # Pa const r_0 = σ_u / sqrt(E) # Pa^(1/2) @@ -64,9 +68,11 @@ end end end +u(x) = VectorValue( udx_max*x[1]/L, -ν*udx_max*x[2]/L, -ν*udx_max*x[3]/L ) + function main(;n,nsteps) - domain = (0,1,0,1,0,1) + domain = (0,L,0,L,0,L) partition = (n,n,n) model = CartesianDiscreteModel(domain,partition) @@ -126,9 +132,14 @@ function main(;n,nsteps) end + e = uh - u + + e_l2 = sqrt(sum(integrate(e*e,trian,quad))) + @test e_l2 < 1.0e-9 + end -main(n=5,nsteps=20) +main(n=5,nsteps=10) -end +end # module From f0244958e11c71eea6b5436f9d65e61dd75102e3 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Tue, 10 Mar 2020 09:20:07 +0100 Subject: [PATCH 3/4] Minor improvements in IsotropicDamageTest --- src/Arrays/Kernels.jl | 5 ++- src/Exports.jl | 1 + src/Geometry/QPointCellFields.jl | 6 ++++ src/Inference/Inference.jl | 41 ++++++++++++++++++++++++ test/GridapTests/IsotropicDamageTests.jl | 18 ++++++----- 5 files changed, 62 insertions(+), 9 deletions(-) diff --git a/src/Arrays/Kernels.jl b/src/Arrays/Kernels.jl index 6890e0678..e55a15a3f 100644 --- a/src/Arrays/Kernels.jl +++ b/src/Arrays/Kernels.jl @@ -220,7 +220,10 @@ function kernel_cache(f::BCasted,x::NumberOrArray...) bs = Base.Broadcast.broadcast_shape(s...) Te = map(numbertype,x) T = return_type(f.f,Te...) - r = zeros(T,bs...) + N = length(bs) + r = Array{T,N}(undef,bs) + ri = testvalue(T) + fill!(r,ri) cache = CachedArray(r) _prepare_cache(cache,x...) end diff --git a/src/Exports.jl b/src/Exports.jl index ada90f936..3367e51c1 100644 --- a/src/Exports.jl +++ b/src/Exports.jl @@ -126,6 +126,7 @@ end @publish FESpaces FEOperator @publish FESpaces FESolver @publish FESpaces apply_statelaw +@publish FESpaces update_state_variables! using Gridap.FESpaces: @law; export @law using Gridap.FESpaces: @statelaw; export @statelaw diff --git a/src/Geometry/QPointCellFields.jl b/src/Geometry/QPointCellFields.jl index aa98a9d2e..a6d22b7b1 100644 --- a/src/Geometry/QPointCellFields.jl +++ b/src/Geometry/QPointCellFields.jl @@ -9,6 +9,12 @@ function QPointCellField(value::Number,trian::Triangulation,quad::CellQuadrature GenericCellField(array, cell_map) end +""" +""" +function CellField(value::Number,trian::Triangulation,quad::CellQuadrature) + QPointCellField(value,trian,quad) +end + struct EvaluatedField{A<:AbstractArray,P<:AbstractArray} <:Field array::A points::P diff --git a/src/Inference/Inference.jl b/src/Inference/Inference.jl index d15766312..53fb44196 100644 --- a/src/Inference/Inference.jl +++ b/src/Inference/Inference.jl @@ -13,6 +13,7 @@ module Inference using DocStringExtensions using FillArrays +using Gridap.Helpers export testvalue export testvalues @@ -128,6 +129,46 @@ function testvalue(::Type{T}) where T<:Fill{E,N} where {E,N} Fill(zero(E),fill(0,N)...) end +function testvalue(::Type{<:Tuple}) + @notimplemented "testvalue on Tuple type only implemented up to 8 tuple elements" +end + +function testvalue(::Type{Tuple{T1,T2,T3,T4,T5,T6,T7,T8}}) where {T1,T2,T3,T4,T5,T6,T7,T8} + (testvalue(T1),testvalue(T2),testvalue(T3),testvalue(T4),testvalue(T5),testvalue(T6),testvalue(T7),testvalue(T8)) +end + +function testvalue(::Type{Tuple{T1,T2,T3,T4,T5,T6,T7}}) where {T1,T2,T3,T4,T5,T6,T7} + (testvalue(T1),testvalue(T2),testvalue(T3),testvalue(T4),testvalue(T5),testvalue(T6),testvalue(T7)) +end + +function testvalue(::Type{Tuple{T1,T2,T3,T4,T5,T6}}) where {T1,T2,T3,T4,T5,T6} + (testvalue(T1),testvalue(T2),testvalue(T3),testvalue(T4),testvalue(T5),testvalue(T6)) +end + +function testvalue(::Type{Tuple{T1,T2,T3,T4,T5}}) where {T1,T2,T3,T4,T5} + (testvalue(T1),testvalue(T2),testvalue(T3),testvalue(T4),testvalue(T5)) +end + +function testvalue(::Type{Tuple{T1,T2,T3,T4}}) where {T1,T2,T3,T4} + (testvalue(T1),testvalue(T2),testvalue(T3),testvalue(T4)) +end + +function testvalue(::Type{Tuple{T1,T2,T3}}) where {T1,T2,T3} + (testvalue(T1),testvalue(T2),testvalue(T3)) +end + +function testvalue(::Type{Tuple{T1,T2}}) where {T1,T2} + (testvalue(T1),testvalue(T2)) +end + +function testvalue(::Type{Tuple{T1}}) where {T1} + (testvalue(T1),) +end + +function testvalue(::Type{Tuple{}}) + () +end + """ testvalues(Ts::DataType...) -> Tuple diff --git a/test/GridapTests/IsotropicDamageTests.jl b/test/GridapTests/IsotropicDamageTests.jl index a7f080230..75d9bf919 100644 --- a/test/GridapTests/IsotropicDamageTests.jl +++ b/test/GridapTests/IsotropicDamageTests.jl @@ -1,7 +1,6 @@ module IsotropicDamageTests using Gridap -using Gridap.FESpaces using LinearAlgebra using Test @@ -32,7 +31,7 @@ function q(r) r_0 + H*(r-r_0) end -function update(ε_in,r_in,d_in) +@law function update(ε_in,r_in,d_in) τ_in = τ(ε_in) if τ_in <= r_in r_out = r_in @@ -54,9 +53,9 @@ end end -@law function dσ(dε_in,ε_in,r_in,d_in) +@law function dσ(dε_in,ε_in,state) - damaged, r_out, d_out = update(ε_in,r_in,d_in) + damaged, r_out, d_out = state if ! damaged return (1-d_out)*σe(dε_in) @@ -95,10 +94,10 @@ function main(;n,nsteps) degree = 2*order quad = CellQuadrature(trian,degree) - r = QPointCellField(r_0,trian,quad) - d = QPointCellField(0.0,trian,quad) + r = CellField(r_0,trian,quad) + d = CellField(0.0,trian,quad) - nls = NLSolver(show_trace=true, method=:newton) + nls = NLSolver(show_trace=false, method=:newton) solver = FESolver(nls) function step(uh_in,factor) @@ -109,7 +108,10 @@ function main(;n,nsteps) U = TrialFESpace(V,[u0,ud,u0,u0]) res(u,v) = inner( ε(v), σ(ε(u),r,d) ) - jac(u,du,v) = inner( ε(v), dσ(ε(du),ε(u),r,d) ) + function jac(u,du,v) + state = update(ε(u),r,d) + inner( ε(v), dσ(ε(du),ε(u),state) ) + end t_Ω = FETerm(res,jac,trian,quad) op = FEOperator(U,V,t_Ω) From 9676cd0b077effafde1ca1baa43a20e1d236bf05 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Tue, 10 Mar 2020 10:00:48 +0100 Subject: [PATCH 4/4] Minor in tests --- test/FESpacesTests/runtests.jl | 2 +- test/GeometryTests/QPointCellFieldsTests.jl | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/test/FESpacesTests/runtests.jl b/test/FESpacesTests/runtests.jl index a18bdf3ec..ac4b8f236 100644 --- a/test/FESpacesTests/runtests.jl +++ b/test/FESpacesTests/runtests.jl @@ -28,7 +28,7 @@ using Test @testset "DivConformingFESpaces" begin include("DivConformingFESpacesTests.jl") end -@testset "CurlConformingFESpaces" begin include("DivConformingFESpacesTests.jl") end +@testset "CurlConformingFESpaces" begin include("CurlConformingFESpacesTests.jl") end @testset "DiscontinuousFESpaces" begin include("DiscontinuousFESpacesTests.jl") end diff --git a/test/GeometryTests/QPointCellFieldsTests.jl b/test/GeometryTests/QPointCellFieldsTests.jl index 79645967d..85fd4e065 100644 --- a/test/GeometryTests/QPointCellFieldsTests.jl +++ b/test/GeometryTests/QPointCellFieldsTests.jl @@ -21,7 +21,6 @@ test_array_of_fields(a,q,s_q) r = QPointCellField(0.0,trian,quad) r_q = evaluate(r,q) -@show r_q === r.array.array test_cell_field(r,q,r_q) end # module