MFront material

.travis.yml

@@ -22,7 +22,7 @@ jobs:
     - stage: "Documentation"
       os: linux
       before_script:
-        - julia --project=docs/ -e 'using Pkg; Pkg.pin(PackageSpec(name="CxxWrap",version="0.8.1"))'
+        - julia --project=docs/ -e 'using Pkg; Pkg.pin(PackageSpec(name="CxxWrap",version="0.8.1")); Pkg.add(PackageSpec(url="https://github.com/JuliaFEM/FEMMaterials.jl", rev="master"))'
       script:
         - julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate(); Pkg.build()'
         - julia --project=docs/ docs/make.jl
@@ -33,4 +33,4 @@ after_success:
   - julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; Coveralls.submit(process_folder())'
 
 before_script:
-  - julia -e 'using Pkg; Pkg.pin(PackageSpec(name="CxxWrap",version="0.8.1"))'
+  - julia -e 'using Pkg; Pkg.pin(PackageSpec(name="CxxWrap", version="0.8.1")); Pkg.add(PackageSpec(url="https://github.com/JuliaFEM/FEMMaterials.jl", rev="master"))'

Project.toml

@@ -7,14 +7,18 @@ version = "0.1.0"
 BinaryProvider = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
 CxxWrap = "1f15a43c-97ca-5a2a-ae31-89f07a497df4"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+FEMBase = "fbcbbc08-f1bf-5204-9233-b69f5d396135"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+Materials = "35fa313e-aa48-52ea-bc82-ddb91c7db87b"
+Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
 
 [compat]
 julia = "1"
 
 [extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Test", "Suppressor"]

src/MFrontInterface.jl

@@ -91,4 +91,10 @@ export load, BehaviourData, get_variable_offset, get_internal_state_variables
 export get_hypothesis, set_time_increment!, set_external_state_variable!
 export get_final_state, update, get_gradients, get_initial_state, integrate
 export get_initial_state, get_parameters, get_external_state_variables
+
+include("mfront_material.jl")
+export MFrontMaterial, MFrontDriverState, MFrontVariableState, MFrontExternalVariableState
+export integrate_material!, update_material!, reset_material!
+export material_preprocess_analysis!, material_preprocess_increment!
+
 end # module MFront`

src/mfront_material.jl

@@ -0,0 +1,148 @@
+using Parameters
+using Tensors
+using Materials
+using FEMMaterials
+
+mgis_bv = MFrontInterface.behaviour
+
+"""
+Variables updated by MFront.
+"""
+@with_kw struct MFrontVariableState <: AbstractMaterialState
+    stress :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
+    jacobian :: SymmetricTensor{4,3} = zero(SymmetricTensor{4,3,Float64})
+end
+
+"""
+Variables passed in for information.
+These drive evolution of the material state.
+"""
+@with_kw mutable struct MFrontDriverState <: AbstractMaterialState
+    time :: Float64 = zero(Float64)
+    strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
+end
+
+"""
+Material external variables in order that is specific to chosen MFront behaviour.
+"""
+@with_kw struct MFrontExternalVariableState <: AbstractMaterialState
+    names :: Array{String,1} = [""]
+    values :: Array{Float64,1} = zeros(Float64, 1)
+end
+
+"""
+MFront material structure.
+
+`lib_path` is the path to the compiled shared library.
+`behaviour` is the loaded material behaviour from the shared library.
+"""
+@with_kw mutable struct MFrontMaterial <: AbstractMaterial
+    drivers :: MFrontDriverState = MFrontDriverState()
+    ddrivers :: MFrontDriverState = MFrontDriverState()
+    variables :: MFrontVariableState = MFrontVariableState()
+    variables_new :: MFrontVariableState = MFrontVariableState()
+
+    external_variables :: MFrontExternalVariableState = MFrontExternalVariableState()
+
+    behaviour :: MFrontInterface.behaviour.BehaviourAllocated
+    behaviour_data :: MFrontInterface.behaviour.BehaviourDataAllocated
+end
+
+function Materials.integrate_material!(material::MFrontMaterial)
+    behaviour = material.behaviour
+    behaviour_data = material.behaviour_data
+
+    mgis_bv.set_time_increment!(behaviour_data, material.ddrivers.time)
+    mgis_bv.revert(material.behaviour_data)
+
+    # setting the external variables (like temperature)
+    for j in 1:length(material.external_variables.names)
+        if material.external_variables.names[j] != ""
+            mgis_bv.set_external_state_variable!(mgis_bv.get_s1(behaviour_data), material.external_variables.names[j], material.external_variables.values[j])
+        end
+    end
+
+    # passing strain from material struct to the mfront interface
+    dstrain = tovoigt(material.ddrivers.strain; offdiagscale=2.0)
+    # now reorder from voigt 11, 22, 33, 23, 13, 12 -> 11, 22, 33, 12, 13, 23
+    # and use Mandel notation (scaling with sqrt(2.0))
+    # https://en.wikipedia.org/wiki/Voigt_notation#Mandel_notation
+    mfront_dstrain = tovoigt(frommandel(SymmetricTensor{2, 3}, [dstrain[1], dstrain[2], dstrain[3], dstrain[6], dstrain[5], dstrain[4]]))
+    gradients = mgis_bv.get_gradients(mgis_bv.get_final_state(behaviour_data))
+    for j in 1:6
+        gradients[j] += mfront_dstrain[j]
+    end
+
+    # tell mfront interface to calculate the tangent
+    # if K[0] is greater than 3.5, the consistent tangent operator must be computed.
+    dummy = zeros(36)
+    dummy[1] = 4.0
+    mgis_bv.set_tangent_operator!(behaviour_data, dummy)
+
+    mgis_bv.integrate(behaviour_data, behaviour)
+
+    stress = [mgis_bv.get_thermodynamic_forces(mgis_bv.get_final_state(behaviour_data))[k] for k in 1:6]
+    jacobian = reshape([mgis_bv.get_tangent_operator(behaviour_data)[k] for k in 1:36], 6, 6)
+
+    # now reorder to voigt 11, 22, 33, 12, 13, 23 -> 11, 22, 33, 23, 13, 12
+    stress = [stress[1], stress[2], stress[3], stress[6], stress[5], stress[4]]
+    stress = frommandel(SymmetricTensor{2, 3}, stress)
+
+    r1 = jacobian[4, :]
+    r3 = jacobian[6, :]
+    jacobian[4, :] .= r3
+    jacobian[6, :] .= r1
+    c1 = jacobian[:, 4]
+    c3 = jacobian[:, 6]
+    jacobian[:, 4] .= c3
+    jacobian[:, 6] .= c1
+    jacobian = frommandel(SymmetricTensor{4, 3}, jacobian)
+
+    variables_new = MFrontVariableState(stress=stress, jacobian=jacobian)
+    material.variables_new = variables_new
+end
+
+function Materials.update_material!(material::MFrontMaterial)
+    mgis_bv.update(material.behaviour_data)
+
+    material.drivers += material.ddrivers
+    material.variables = material.variables_new
+
+    reset_material!(material)
+end
+
+function Materials.reset_material!(material::MFrontMaterial)
+    material.ddrivers = typeof(material.ddrivers)()
+    material.variables_new = typeof(material.variables_new)()
+end
+
+# other material_* functions are used from FEMMaterials
+
+import FEMMaterials: update_ip!
+import FEMBase: Element, Hex8
+
+"""
+Initializes integration point `ip` for data storage of both `variables` and `drivers` at simulation start `time`.
+"""
+function FEMMaterials.material_preprocess_analysis!(material::MFrontMaterial, element::Element{Hex8}, ip, time)
+    update_ip!(material, ip, time)
+    # Read parameter values
+    values = element("external_variables", ip, time)
+    material.external_variables = MFrontExternalVariableState(material.external_variables.names, values)
+end
+
+"""
+Initializes integration point `ip` for data storage of both `variables` and `drivers` at simulation start `time`.
+Updates external variables, e.g. temperature, stored in `ip` to material
+"""
+function FEMMaterials.material_preprocess_increment!(material::MFrontMaterial, element::Element{Hex8}, ip, time)
+
+    values = element("external_variables", ip, time)
+    material.external_variables = MFrontExternalVariableState(material.external_variables.names, values)
+
+    # Update time increment
+    dtime = time - material.drivers.time
+    material.ddrivers.time = dtime
+
+    return nothing
+end

test/runtests.jl

@@ -15,5 +15,17 @@ eps = 1.e-12
     if Sys.islinux()
         @testset "Norton model" begin include("test_norton_model.jl") end
         include("test_show_methods.jl")
+
+        @testset "test MFront ideal plastic material model" begin
+        include("test_isotropic_linear_hardening_plasticity.jl")
+        end
+
+        @testset "test MFront ideal plastic material model with shear strain" begin
+        include("test_isotropic_linear_hardening_plasticity_shear.jl")
+        end
+
+        @testset "test MFront together with FEMMaterials" begin
+        include("test_mfront_mecamatso.jl")
+        end
     end
 end

test/test_isotropic_linear_hardening_plasticity.jl

@@ -0,0 +1,68 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/MFrontInterface.jl/blob/master/LICENSE
+
+using MFrontInterface
+using Test
+using Tensors
+
+mgis_bv = MFrontInterface.behaviour
+
+function MaterialTest()
+    lib_path = "test_plasticity/libBehaviour.so"
+    behaviour_name = "IsotropicLinearHardeningPlasticity"
+    hypothesis = mgis_bv.Tridimensional
+
+    behaviour = load(lib_path, behaviour_name, hypothesis)
+    behaviour_data = BehaviourData(behaviour)
+
+    ext_variable_names = [mgis_bv.get_name(mgis_bv.get_external_state_variables(behaviour)[i]) for i in 1:mgis_bv.length(mgis_bv.get_external_state_variables(behaviour))]
+    ext_variable_values = zeros(length(ext_variable_names))
+    ext_vatiable_state = MFrontExternalVariableState(names=ext_variable_names, values=ext_variable_values)
+
+    return MFrontMaterial(behaviour=behaviour, behaviour_data=behaviour_data, external_variables=ext_vatiable_state)
+end
+
+mat = MaterialTest()
+
+d = mat.behaviour_data
+o = get_variable_offset(get_internal_state_variables(mat.behaviour),
+                        "EquivalentPlasticStrain",
+                        get_hypothesis(mat.behaviour))
+
+dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0)
+ddrivers = MFrontDriverState(time = 0.25, strain = 0.25*dstrain_dtime)
+mat.ddrivers = ddrivers
+
+integrate_material!(mat)
+update_material!(mat)
+@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, 1.0e6*[50.0, 0.0, 0.0, 0.0, 0.0, 0.0]))
+
+mat.ddrivers = ddrivers
+integrate_material!(mat)
+update_material!(mat)
+@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, 1.0e6*[100.0, 0.0, 0.0, 0.0, 0.0, 0.0]))
+@test isapprox(get_internal_state_variables(get_initial_state(d))[o], 0.0; atol=1.0e-6)
+
+dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.5, -0.5, 0.0, 0.0, 0.0]; offdiagscale=2.0)
+ddrivers = MFrontDriverState(time = 0.25, strain = 0.25*dstrain_dtime)
+mat.ddrivers = ddrivers
+integrate_material!(mat)
+update_material!(mat)
+@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, 1.0e6*[100.0, 0.0, 0.0, 0.0, 0.0, 0.0]); atol=1.0e-6)
+@test isapprox(get_internal_state_variables(get_initial_state(d))[o], 0.25*1.0e-3)
+
+dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64}, -1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0)
+ddrivers = MFrontDriverState(time = 0.25, strain = 0.25*dstrain_dtime)
+mat.ddrivers = ddrivers
+integrate_material!(mat)
+update_material!(mat)
+@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, 1.0e6*[50.0, 0.0, 0.0, 0.0, 0.0, 0.0]); atol=1.0e-6)
+
+dstrain_dtime = (-0.75*fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0)
+                -0.25*fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.5, -0.5, 0.0, 0.0, 0.0]; offdiagscale=2.0))
+ddrivers = MFrontDriverState(time = 1.0, strain = dstrain_dtime)
+mat.ddrivers = ddrivers
+integrate_material!(mat)
+integrate_material!(mat)
+update_material!(mat)
+@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, 1.0e6*[-100.0, 0.0, 0.0, 0.0, 0.0, 0.0]))

test/test_isotropic_linear_hardening_plasticity_shear.jl

@@ -0,0 +1,68 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/MFrontInterface.jl/blob/master/LICENSE
+
+using MFrontInterface
+using Test
+using Tensors
+
+mgis_bv = MFrontInterface.behaviour
+
+function MaterialTest()
+    lib_path = "test_plasticity/libBehaviour.so"
+    behaviour_name = "IsotropicLinearHardeningPlasticity"
+    hypothesis = mgis_bv.Tridimensional
+
+    behaviour = load(lib_path, behaviour_name, hypothesis)
+    behaviour_data = BehaviourData(behaviour)
+
+    ext_variable_names = [mgis_bv.get_name(mgis_bv.get_external_state_variables(behaviour)[i]) for i in 1:mgis_bv.length(mgis_bv.get_external_state_variables(behaviour))]
+    ext_variable_values = zeros(length(ext_variable_names))
+    ext_vatiable_state = MFrontExternalVariableState(names=ext_variable_names, values=ext_variable_values)
+
+    return MFrontMaterial(behaviour=behaviour, behaviour_data=behaviour_data, external_variables=ext_vatiable_state)
+end
+
+mat = MaterialTest()
+
+E = 200.0e9
+nu = 0.3
+syield = 100.0e6
+
+times = [0.0]
+loads = [0.0]
+dt = 0.5
+G = 0.5*E/(1+nu)
+
+ea = 2*syield/(sqrt(3)*G)
+# Go to elastic border
+push!(times, times[end]+dt)
+push!(loads, loads[end] + ea*dt)
+ # Proceed to plastic flow
+push!(times, times[end]+dt)
+push!(loads, loads[end] + ea*dt)
+ # Reverse direction
+push!(times, times[end]+dt)
+push!(loads, loads[end] - ea*dt)
+ # Continue and pass yield criterion
+push!(times, times[end]+dt)
+push!(loads, loads[end] - 2*ea*dt)
+stresses = [copy(tovoigt(mat.variables.stress))]
+for i=2:length(times)
+    dtime = times[i]-times[i-1]
+    dstrain12 = loads[i]-loads[i-1]
+    dstrain = [0.0, 0.0, 0.0, 0.0, 0.0, dstrain12]
+    dstrain_ = fromvoigt(SymmetricTensor{2,3,Float64}, dstrain; offdiagscale=2.0)
+    ddrivers = MFrontDriverState(time = dtime, strain = dstrain_)
+    mat.ddrivers = ddrivers
+    integrate_material!(mat)
+    update_material!(mat)
+    push!(stresses, copy(tovoigt(mat.variables.stress)))
+end
+
+for i in 1:length(times)
+    @test isapprox(stresses[i][1:5], zeros(5); atol=1e-6)
+end
+s12 = [s[6] for s in stresses]
+
+s12_expected = [0.0, syield/sqrt(3.0), syield/sqrt(3.0), 0.0, -syield/sqrt(3.0)]
+@test isapprox(s12, s12_expected; rtol=1.0e-2)

test/test_mfront_mecamatso.jl

@@ -0,0 +1,39 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/MFrontInterface.jl/blob/master/LICENSE
+
+using MFrontInterface
+using Test
+using Materials, Tensors, FEMMaterials, FEMBase, Test
+
+mgis_bv = MFrontInterface.behaviour
+
+function MaterialTest()
+    lib_path = "test_plasticity/libBehaviour.so"
+    behaviour_name = "IsotropicLinearHardeningPlasticity"
+    hypothesis = mgis_bv.Tridimensional
+
+    behaviour = load(lib_path, behaviour_name, hypothesis)
+    behaviour_data = BehaviourData(behaviour)
+
+    ext_variable_names = [mgis_bv.get_name(mgis_bv.get_external_state_variables(behaviour)[i]) for i in 1:mgis_bv.length(mgis_bv.get_external_state_variables(behaviour))]
+    ext_variable_values = zeros(length(ext_variable_names))
+    ext_vatiable_state = MFrontExternalVariableState(names=ext_variable_names, values=ext_variable_values)
+
+    return MFrontMaterial(behaviour=behaviour, behaviour_data=behaviour_data, external_variables=ext_vatiable_state)
+end
+
+analysis, problem, element, bc_elements, ip = get_one_element_material_analysis(:(Main.MaterialTest))
+
+temperature = 293.15
+update!(element, "external_variables", [temperature])
+
+times = [0.0, 1.0, 2.0, 3.0]
+loads = [0.0, 1.0e-3, -1.0e-3, 1.0e-3]
+loading = AxialStrainLoading(times, loads)
+update_bc_elements!(bc_elements, loading)
+analysis.properties.t1 = maximum(times)
+
+run!(analysis)
+s33 = [tovoigt(ip("stress", t))[3] for t in times]
+s33_expected = 1e6*[0.0, 100.0, -100.0, 100.0]
+@test isapprox(s33, s33_expected; rtol=1.0e-2)