Final retouches for historical variables

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

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

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

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
 

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
 

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

test/GridapTests/IsotropicDamageTests.jl

@@ -0,0 +1,147 @@
+module IsotropicDamageTests
+
+using Gridap
+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
+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
+
+@law 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,state)
+
+  damaged, r_out, d_out = state
+
+  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
+
+u(x) = VectorValue( udx_max*x[1]/L, -ν*udx_max*x[2]/L, -ν*udx_max*x[3]/L )
+
+function main(;n,nsteps)
+
+  domain = (0,L,0,L,0,L)
+  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 = CellField(r_0,trian,quad)
+  d = CellField(0.0,trian,quad)
+
+  nls = NLSolver(show_trace=false, 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) )
+    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_Ω)
+
+    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
+
+  e = uh - u
+
+  e_l2 = sqrt(sum(integrate(e*e,trian,quad)))
+  @test e_l2 < 1.0e-9
+
+end
+
+main(n=5,nsteps=10)
+
+
+end # module

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