Draft: Possible anisotropy concatenation

[miguelmaso]

FYI @jmartfrut @RogelioOrtigosa , just a reminder to keep the code at some place. This is the cleanest implementation I've reached. To be discussed:

abstract type Elasto end

abstract type IsoElastic <: Elasto end

abstract type AnisoElastic <: Elasto end

struct Foo <: IsoElastic end

struct Bar <: AnisoElastic end

struct MultiAniso <: AnisoElastic  # Not exported
    aniso::NTuple{N,AnisoElastic} where N
end

struct IsoMultiAniso <: AnisoElastic  # Not exported
    iso::IsoElastic
    aniso::MultiAniso
end

Base.hcat(a::AnisoElastic...) = MultiAniso(a)

+(a::IsoElastic, b::AnisoElastic) = IsoMultiAniso(a,MultiAniso((b,)))

+(a::IsoElastic, b::MultiAniso) = IsoMultiAniso(a,b)

transpose(x::NTuple{N, <:Tuple{<:Function, <:Function, <:Function}}) where N = map(i -> getindex.(x, i), 1:3)

function (obj::Foo)(args...)
    Ψ(F) = 0.0
    ∂Ψ∂F(F) = [0.0, 0.0, 0.0]
    ∂Ψ∂FF(F) = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    (Ψ, ∂Ψ∂F, ∂Ψ∂FF)
end

function (obj::Bar)(args...)
    Ψ(F,N) = 0.0
    ∂Ψ∂F(F,N) = [0.0, 0.0, 0.0]
    ∂Ψ∂FF(F,N) = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    (Ψ, ∂Ψ∂F, ∂Ψ∂FF)
end

function (obj::MultiAniso)(args...)
    DΨ = map(a -> a(args...), obj.aniso)
    Ψα, ∂Ψα∂F, ∂Ψα∂FF = transpose(DΨ)
    Ψ(F,N) = mapreduce(Ψi -> Ψi(F,N), Ψα)
    ∂Ψ∂F(F,N) = mapreduce(∂Ψi∂F -> ∂Ψi∂F(F,N), ∂Ψα∂F)
    ∂Ψ∂FF(F,N) = mapreduce(∂Ψi∂FF -> ∂Ψi∂FF(F,N), ∂Ψα∂FF)
    (Ψ, ∂Ψ∂F, ∂Ψ∂FF)
end

function (obj::IsoMultiAniso)(args...)
    Ψi, ∂Ψi∂F, ∂Ψi∂FF = obj.iso(args...)
    Ψa, ∂Ψa∂F, ∂Ψa∂FF = obj.aniso(args...)
    Ψ(F,N) = Ψi(F) + Ψa(F,N)
    ∂Ψ∂F(F,N) = ∂Ψi∂F(F) + ∂Ψa∂F(F,N)
    ∂Ψ∂FF(F,N) = ∂Ψi∂FF(F) + ∂Ψa∂FF(F,N)
    (Ψ, ∂Ψ∂F, ∂Ψ∂FF)
end

f = Foo()
b1 = Bar()
b2 = Bar()

fb = f + b1
bb = [b1 b2]
fbb = f + [b1 b2]

∂Ψ, ∂Ψ∂F, ∂Ψ∂FF = fb()
∂Ψ, ∂Ψ∂F, ∂Ψ∂FF = bb()
∂Ψ, ∂Ψ∂F, ∂Ψ∂FF = fbb()

[jmartfrut]

What about the following proposal

struct MultiIsoElastic <: IsoElastic
  Model1::IsoElastic
  Model2::IsoElastic
  Kinematic
  function MultiIsoElastic(model1::IsoElastic,model2::IsoElastic)
    @assert model1.Kinematic == model2.Kinematic
    new(model1,model2,model1.Kinematic)
  end
  function (obj::MultiIsoElastic)(Λ::Float64=1.0)
    DΨ1 = obj.Model1(Λ)
    DΨ2 = obj.Model2(Λ)
    Ψ, ∂Ψ, ∂∂Ψ = map((ψ1,ψ2) -> (x) -> ψ1(x) + ψ2(x), DΨ1, DΨ2)
    return (Ψ, ∂Ψ, ∂∂Ψ)
  end
end


struct MultiAnisoElastic{A,B} <: AnisoElastic
  Model1::A
  Model2::B
  Kinematic
  function MultiAnisoElastic(model1::A,model2::B) where{A,B}
    @assert model1.Kinematic == model2.Kinematic
    new{A,B}(model1,model2,model1.Kinematic)
  end
  function (obj::MultiAnisoElastic{<:AnisoElastic,<:AnisoElastic})(Λ::Float64=1.0)
    DΨ1 = obj.Model1(Λ)
    DΨ2 = obj.Model2(Λ)
    Ψ, ∂Ψ, ∂∂Ψ = map((ψ1,ψ2) -> (x,N) -> ψ1(x,N) + ψ2(x,N), DΨ1, DΨ2)
    return (Ψ, ∂Ψ, ∂∂Ψ)
  end
  function (obj::MultiAnisoElastic{<:IsoElastic,<:AnisoElastic})(Λ::Float64=1.0)
    DΨ1 = obj.Model1(Λ)
    DΨ2 = obj.Model2(Λ)
    Ψ, ∂Ψ, ∂∂Ψ = map((ψ1,ψ2) -> (x,N) -> ψ1(x) + ψ2(x,N), DΨ1, DΨ2)
    return (Ψ, ∂Ψ, ∂∂Ψ)
  end
  function (obj::MultiAnisoElastic{<:AnisoElastic,<:IsoElastic})(Λ::Float64=1.0)
    DΨ1 = obj.Model1(Λ)
    DΨ2 = obj.Model2(Λ)
    Ψ, ∂Ψ, ∂∂Ψ = map((ψ1,ψ2) -> (x,N) -> ψ1(x,N) + ψ2(x), DΨ1, DΨ2)
    return (Ψ, ∂Ψ, ∂∂Ψ)
  end
end
 
function (+)(Model1::IsoElastic, Model2::AnisoElastic)
  MultiAnisoElastic(Model1,Model2)
end
function (+)(Model1::AnisoElastic, Model2::IsoElastic)
  MultiAnisoElastic(Model1,Model2)
end
function (+)(Model1::AnisoElastic, Model2::AnisoElastic)
  MultiAnisoElastic(Model1,Model2)
end
function (+)(Model1::IsoElastic, Model2::IsoElastic)
  MultiIsoElastic(Model1,Model2)
end

∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
N = VectorValue(1.0, 2.0, 3.0)
model1 = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=2.0)
model2 = TransverseIsotropy3D(μ=μParams[5], α=μParams[6], β=μParams[7])
model=model1+model2+model1+model2

[miguelmaso]

I find the second option more complicated

function (obj::MultiAnisoElastic{<:AnisoElastic,<:IsoElastic})(Λ::Float64=1.0)

You need to deep into the internals of the model to do multiple dispatch. In the first verion, multiple dispatch is easier

Ψ, ∂Ψ, ∂∂Ψ = map((ψ1,ψ2) -> (x) -> ψ1(x) + ψ2(x), DΨ1, DΨ2)

When there are only two models, I find better to use the sum, easier to undertand.

From the point of view of usage, which is more likely to occur: multiple fiber directions or multiple isotropic models? The first approach is reflecting the multiple fiber directions (analog to multiple viscous branches).

Finally, when state variables are involved, you'll need to definitely go for a Tuple. Recursion is difficult to unpack...

[miguelmaso]

Latest iteration: updated function (obj::MultiAniso)(args...) with a tuple of anisotropy directions N:

abstract type Elasto end

abstract type IsoElastic <: Elasto end

abstract type AnisoElastic <: Elasto end

struct Foo <: IsoElastic end

struct Bar <: AnisoElastic end

struct MultiAniso <: AnisoElastic
    aniso::NTuple{N,AnisoElastic} where N
end

struct IsoMultiAniso <: AnisoElastic
    iso::IsoElastic
    aniso::MultiAniso
end

Base.hcat(a::AnisoElastic...) = MultiAniso(a)

+(a::IsoElastic, b::AnisoElastic) = IsoMultiAniso(a,MultiAniso((b,)))

+(a::IsoElastic, b::MultiAniso) = IsoMultiAniso(a,b)

transpose(x::NTuple{N, <:Tuple{<:Function, <:Function, <:Function}}) where N = map(i -> getindex.(x, i), 1:3)

function (obj::Foo)(args...)
    Ψ(F) = 0.0
    ∂Ψ∂F(F) = [0.0, 0.0, 0.0]
    ∂Ψ∂FF(F) = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    (Ψ, ∂Ψ∂F, ∂Ψ∂FF)
end

function (obj::Bar)(args...)
    Ψ(F,N) = 0.0
    ∂Ψ∂F(F,N) = [0.0, 0.0, 0.0]
    ∂Ψ∂FF(F,N) = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    (Ψ, ∂Ψ∂F, ∂Ψ∂FF)
end

function (obj::MultiAniso)(args...)
    DΨ = map(a -> a(args...), obj.aniso)
    Ψα, ∂Ψα∂F, ∂Ψα∂FF = transpose(DΨ)
    Ψ(F,N) = mapreduce((Ψi,Ni) -> Ψi(F,Ni), +, Ψα, N)
    ∂Ψ∂F(F,N) = mapreduce((∂Ψi∂F,Ni) -> ∂Ψi∂F(F,Ni), +, ∂Ψα∂F, N)
    ∂Ψ∂FF(F,N) = mapreduce((∂Ψi∂FF, Ni) -> ∂Ψi∂FF(F,Ni), +, ∂Ψα∂FF, N)
    (Ψ, ∂Ψ∂F, ∂Ψ∂FF)
end

function (obj::IsoMultiAniso)(args...)
    Ψi, ∂Ψi∂F, ∂Ψi∂FF = obj.iso(args...)
    Ψa, ∂Ψa∂F, ∂Ψa∂FF = obj.aniso(args...)
    Ψ(F,N) = Ψi(F) + Ψa(F,N)
    ∂Ψ∂F(F,N) = ∂Ψi∂F(F) + ∂Ψa∂F(F,N)
    ∂Ψ∂FF(F,N) = ∂Ψi∂FF(F) + ∂Ψa∂FF(F,N)
    (Ψ, ∂Ψ∂F, ∂Ψ∂FF)
end

f = Foo()
b1 = Bar()
b2 = Bar()

fb = f + b1
bb = [b1 b2]
fbb = f + [b1 b2]

∂Ψ, ∂Ψ∂F, ∂Ψ∂FF = fb()
∂Ψ, ∂Ψ∂F, ∂Ψ∂FF = bb()
∂Ψ, ∂Ψ∂F, ∂Ψ∂FF = fbb()

NOTE: not sure if we can keep Ψ(F,N) = mapreduce((Ψi,Ni) -> Ψi(F,Ni), +, Ψα, N) or we need to unpack the isotropy directions:

- Ψ(F,N) = mapreduce((Ψi,Ni) -> Ψi(F,Ni), +, Ψα, N)
+ Ψ(F,N...) = mapreduce((Ψi,Ni) -> Ψi(F,Ni), +, Ψα, N)

If so, the same change should be applied the isotropic-anisotropic model:

- Ψ(F,N) = Ψi(F) + Ψa(F,N)
+ Ψ(F,N...) = Ψi(F) + Ψa(F,N...)

Finally, if we must unpack the anisotropy vectors N (I'd like to avoid it), unpacking will conflict with the state variables. I guess we'll need to do some magic to manage the variables...

[miguelmaso]

Finally, if we must unpack the anisotropy vectors N (I'd like to avoid it), unpacking will conflict with the state variables. I guess we'll need to do some magic to manage the variables...

Probably we should implement a macro to tackle this issue