text:
#
# This file is part of module
#   ModiaMedia (ModiaMedia/src/ModiaMedia.jl)
#
using LinearAlgebra

const MoistAir_excludeEnthalpyOfFormation = true                            # If true, enthalpy of formation Hf is not included in specific enthalpy hc
const MoistAir_referenceEnthalpy          = ReferenceEnthalpy_UserDefined   # Choice of reference enthalpy
const MoistAir_h_offset                   = 0.0                             # User defined offset for reference enthalpy, if ReferenceEnthalpy = ReferenceEnthalpy_UserDefined

### Data structures

"""
    data = MoistAirData(steam::SingleGasNasaData, dryair::SingleGasNasaData)

Generate an `MoistAirData` object containing the data
for a moist air model.
"""
mutable struct MoistAirData
    Water::Int                  # "Index of water (in substanceNames, massFractions X, etc.)"
    Air::Int                    # "Index of air (in substanceNames, massFractions X, etc.)";
    k_mair::Float64             # "Ratio of molar weights";
    dryair::SingleGasNasaData   # Data of dry air
    steam::SingleGasNasaData    # Data of steam
    MMX::SVector{2,Float64}     # "Molar masses of components";

    MoistAirData(steam::SingleGasNasaData, dryair::SingleGasNasaData) =
        new(1,2,steam.MM/dryair.MM, dryair, steam, SVector{2,Float64}(steam.MM,dryair.MM))
end

"""
    medium = MoistAir()

Generate an `MoistAir <: CondensingGases` medium object.
Valid for T = 190 K ... 647 K.
"""
struct MoistAir <: CondensingGases
    infos::FluidInfos
    fluidConstants::SVector{2,IdealGasFluidConstants}
    fluidLimits::FluidLimits
    data::MoistAirData
end

"""
    state = MoistAirState(Medium, p, T, X)

Generate a `MoistAirState <: MixtureThermodynamicState` object containing
pressure `p` [Pa], temperature `T` [K], and a vector of mass fractions as states,
where `X[1]` is the mass fraction of Steam and `X[2]` is the mass fraction of dry air.
If argument `X` has only one element, `X[2]` is computed from `X[1]` and stored in the state.
"""
mutable struct MoistAirState <: MixtureThermodynamicState
    Medium::MoistAir
    p::Float64
    T::Float64
    X::MVector{2,Float64}

    MoistAirState(Medium::MoistAir, p::Float64, T::Float64, X::AbstractVector) =
        new(Medium, p, T, MVector{2,Float64}(X[1], length(X)>1 ? X[2] : max(1.0 - X[1], 0.0)))
end

function MoistAir(fluidConstants_H2O::IdealGasFluidConstants,
                  fluidConstants_Air::IdealGasFluidConstants,
                  steam::SingleGasNasaData,
                  dryair::SingleGasNasaData)
    infos = FluidInfos(mediumName           = "MoistAir",
                       substanceNames       = ["water", "air"],
                       extraPropertiesNames = fill("",0),
                       ThermoStates         = IndependentVariables_pTX,
                       singleState          = false,
                       reducedX             = true,
                       fixedX               = false,
                       reference_p          = 101325.0,
                       reference_T          = 298.15,
                       reference_X          = [0.01, 0.99],
                       p_default            = 101325.0,
                       T_default            = 293.15,
                       X_default            = [0.01, 0.99])

    # these limits should be fixed after the specific enthalpy can be calculated below 273.15K; the high end also has problems as well.
    fluidLimits = FluidLimits(TMIN = 190,
                              TMAX = 647)
    data = MoistAirData(steam, dryair)

    infos.h_default = h_pTX(data, infos.p_default, infos.T_default, infos.X_default)

    # HMIN = h_pTX(MoistAirData(data),p_default,280.0,X_default),
    # HMAX = h_pTX(MoistAirData(data),p_default,600.0,X_default))

    Medium = MoistAir(infos,
                      SVector{2,IdealGasFluidConstants}([fluidConstants_H2O,
                                                         fluidConstants_Air]),
                      fluidLimits, data)

    return Medium
end

"Spline interpolation of two functions"
function spliceFunction(
    pos::Float64,                    # Returned value for x-deltax >= 0
    neg::Float64,                    # Returned value for x+deltax <= 0
    x::Float64,                      # Function argument
    deltax::Float64                  # Region around x with spline interpolation
    )::Float64

    scaledX1 = x/deltax
    scaledX = scaledX1*asin(1)
    if scaledX1<=-0.999999999
        y = 0.0
    elseif scaledX1>=0.999999999
        y = 1.0
    else
        y = (tanh(tan(scaledX))+1.0)/2
    end
    return pos*y+(1.0-y)*neg
end

"""
    saturationPressureOfLiquidWater(Tsat)

Return saturation pressure of liquid water as a function of saturation temperature Tsat
in the range of 273.16 to 647.096 K
"""
function saturationPressureOfLiquidWater(Tsat::Float64)::Float64
    Tcritical=647.096
    pcritical=22.064e6
    r1=(1 - Tsat/Tcritical)

    a = SVector{6,Float64}(-7.85951783,1.84408259,-11.7866497,22.6807411,-15.9618719,1.80122502)
    n = SVector{6,Float64}(1.0,1.5,3.0,3.5,4.0,7.5)
    psat = exp(((a[1]*r1^n[1] + a[2]*r1^n[2] + a[3]*r1^n[3] + a[4]*r1^n[4]
     + a[5]*r1^n[5] + a[6]*r1^n[6])*Tcritical)/Tsat)*pcritical

     return psat
end

"""
    sublimationPressureIce(Tsat)

Return sublimation pressure of water as a function of saturation temperature `Tsat`
between 190 and 273.16 K
"""
function sublimationPressureIce(Tsat::Float64)::Float64
    Ttriple=273.16                # Triple point temperature
    ptriple=611.657                # Triple point pressure
    r1=Tsat/Ttriple                     # Common subexpression
    a1=-13.9281690
    a2=34.7078238
    n1=-1.5
    n2=-1.25
    psat = exp(a1-a1*r1^n1+a2-a2*r1^n2)*ptriple

    return psat
end

"""
    saturationPressureOfWater(Tsat)

Return saturation pressure of water as a function of
saturation temperature `Tsat` between 190 K and 647.096 K
"""
saturationPressureOfWater(Tsat::Float64) = spliceFunction(saturationPressureOfLiquidWater(Tsat), sublimationPressureIce(Tsat), Tsat-273.16, 1.0)

"Return specific enthalpy of water (solid/liquid) near atmospheric pressure from temperature T"
enthalpyOfWater(T::Float64) = spliceFunction(4200.0*(T-273.15), 2050.0*(T-273.15)-333000.0, T-273.16, 0.1)

function h_pTX(data::MoistAirData,p::Float64,T::Float64,X::AbstractVector)
    # steam.MM/dryAir.MM
    k_mair = data.k_mair

    p_steam_sat = saturationPressureOfWater(T)
    X_sat = min(p_steam_sat*k_mair/max(100*eps(), p-p_steam_sat)*(1-X[1]), 1.0)
    X_liquid = max(X[1]-X_sat, 0.0);
    X_steam = X[1]-X_liquid;
    X_air = 1-X[1];

    h = h_Tlow(data.steam , T, refChoice=ReferenceEnthalpy_UserDefined, h_off=46479.819+2501014.5)*X_steam +
        h_Tlow(data.dryair, T, refChoice=ReferenceEnthalpy_UserDefined, h_off=25104.684)*X_air +
        enthalpyOfWater(T)*X_liquid

    return h
end

function u_pTX(data::MoistAirData,p::Float64,T::Float64,X::AbstractVector)
    # steam.MM/dryAir.MM
    k_mair = data.k_mair

    p_steam_sat = saturationPressureOfWater(T)
    X_sat = min(p_steam_sat*k_mair/max(100*eps(), p-p_steam_sat)*(1-X[1]), 1.0)
    X_liquid = max(X[1]-X_sat, 0.0);
    X_steam = X[1]-X_liquid;
    X_air = 1-X[1];
    R_gas = data.dryair.R*X_air/(1-X_liquid) + data.steam.R*X_steam/(1-X_liquid);

    u = h_Tlow(data.steam , T, refChoice=ReferenceEnthalpy_UserDefined, h_off=46479.819+2501014.5)*X_steam +
        h_Tlow(data.dryair, T, refChoice=ReferenceEnthalpy_UserDefined, h_off=25104.684)*X_air +
        enthalpyOfWater(T)*X_liquid - R_gas*T
    return u
end

T_phX(data::MoistAirData, p::Float64, h::Float64, X::AbstractVector) = ModiaMath.solveOneNonlinearEquation(T->h-h_pTX(data,p,T,X), 190.0, 647; u_nominal=300.0)
p_dTX(data::MoistAirData,d,T,X::AbstractVector) = d*(data.steam.R*X[1]+data.dryair.R*(1.0-X[1]))*T
T_pTX(data::MoistAirData, p_a::Float64, dp::Float64, T_a::Float64, X_a::AbstractVector) = T_phX(data, p_a+dp, h_pTX(data,p_a,T_a,X_a), X_a)

setState_pTX(m::MoistAir,p,T,X::AbstractVector) = MoistAirState(m,p,T,X)
setState_phX(m::MoistAir,p,h,X::AbstractVector) = MoistAirState(m,p,T_phX(m.data,p,h,X),X)
setState_dTX(m::MoistAir,d,T,X::AbstractVector) = MoistAirState(m,p_dTX(m.data,d,T,X),T,X)
isenthalpicState(m::MoistAir, state::MoistAirState, dp::Float64) = MoistAirState(m, state.p+dp, T_pTX(m.data,state.p,dp,state.T,state.X), state.X)

setState_pTX!(state::MoistAirState,p,T,X::AbstractVector) = begin state.p=p; state.T=T; state.X[1]=X[1]; nothing end
setState_phX!(state::MoistAirState,p,h,X::AbstractVector) = begin state.p=p; state.T=T_phX(state.Medium.data,p,h,X); state.X[1]=X[1]; nothing end
setState_psX!(state::MoistAirState,p,s,X::AbstractVector) = begin state.p=p; state.T=T_ps(state.Medium,p,s); nothing end
setState_dTX!(state::MoistAirState,d,T,X::AbstractVector) = begin state.p=p_dTX(state.Medium.data,d,T,X); state.T=T; state.X[1]=X[1]; nothing end
isenthalpicState!(state_b::MoistAirState, state_a::MoistAirState, dp::Float64) = begin state_b.p    = state_a.p+dp;
                                                                                state_b.T    = T_pTX(state.Medium.data,state_a.p,dp,state_a.T,state_a.X);
                                                                                state_b.X[1] = X[1]; nothing end

pressure(               m::MoistAir, state::MoistAirState)::Float64 = state.p
temperature(            m::MoistAir, state::MoistAirState)::Float64 = state.T
gasConstant(            m::MoistAir, state::MoistAirState)::Float64 = m.data.dryair.R*(1-state.X[1]) + m.data.steam.R*state.X[1]
density(                m::MoistAir, state::MoistAirState)::Float64 = state.p/(gasConstant(m,state)*state.T)
specificEnthalpy(       m::MoistAir, state::MoistAirState)::Float64 = h_pTX(m.data,state.p,state.T,state.X)
specificInternalEnergy( m::MoistAir, state::MoistAirState)::Float64 = u_pTX(m.data,state.p,state.T,state.X)

#=
            function dynamicViscosity(fluid::Any,
                )
                    eta = 1e-6*Polynomials_Temp.evaluateWithRange([9.7391102886305869E-15,-3.1353724870333906E-11,4.3004876595642225E-08,-3.8228016291758240E-05,5.0427874367180762E-02,1.7239260139242528E+01], Cv.to_degC(123.15), Cv.to_degC(1273.15), Cv.to_degC(state.T))

                return
                #= @extends dynamicViscosity
                 =#
            end

=#

function standardCharacteristics(m::MoistAir)::Dict{AbstractString,Any}
    p     = m.infos.reference_p
    T     = collect( range(m.fluidLimits.TMIN, m.fluidLimits.TMAX, length=101) )
    X1    = [0.01, 0.02, 0.03]
    nT    = length(T)
    nX1   = length(X1)
    h     = zeros(nT,nX1)
    u     = zeros(nT,nX1)
    #cp    = zeros(nT,nX1)
    d     = zeros(nT,nX1)

    for i in 1:nT
        for j in 1:nX1
            state   = setState_pTX(m,p,T[i],[X1[j]])
            h[i,j]  = specificEnthalpy(      state)
            u[i,j]  = specificInternalEnergy(state)
            #cp[i,j] = specificHeatCapacityCp(state)
            d[i,j]  = to_DensityDisplayUnit( density(state) )
        end
    end

    mediumDict = Dict{AbstractString,Any}()
    mediumDict["p"]  = p
    mediumDict["T"]  = uconvert.(u"°C", T*1u"K")
    mediumDict["h"]  = h*1u"J/kg"
    mediumDict["u"]  = u*1u"J/kg"
    #mediumDict["cp"] = cp*1u"J/(kg*K)"
    mediumDict["d"]  = d*1u"g/cm^3"

    return mediumDict
end

function standardPlot(m::MoistAir; figure=1)
    mediumDict = standardCharacteristics(m)
    p = mediumDict["p"]
    ModiaMath.plot(mediumDict, [("h", "u"), "d"], xAxis="T", heading=m.infos.mediumName*" (p=$p Pa, X[Water] = 0.01, 0.02, 0.03)", figure=figure)
end


text:
function write_data(	file_path::String,
					ax::Float64,
					ay::Float64,
					number_of_frames::Array{Int64, 1},
					deltat::Float64,
					K::Array{Int64, 1},
					DE::Array{Float64, 1})

	file_stream::IOStream = open(file_path, "w")

	@printf(file_stream, "%s", "<data>\n")

	write_key(file_stream, "ax", ax)
	write_key(file_stream, "ay", ay)
	write_key(file_stream, "number_of_frames", number_of_frames)
	write_key(file_stream, "deltat", deltat)
	write_key(file_stream, "K", K)
	write_key(file_stream, "DE", DE)

	@printf(file_stream, "%s", "</data>")

	close(file_stream)

	nothing
end

ext: jl


text:
"""
    ITensor([::Type{ElT} = Float64, ][flux::QN = QN(), ]inds)
    ITensor([::Type{ElT} = Float64, ][flux::QN = QN(), ]inds::Index...)

Construct an ITensor with BlockSparse storage filled with `zero(ElT)` where the nonzero blocks are determined by `flux`.

If `ElT` is not specified it defaults to `Float64`.
"""
function ITensor(::Type{ElT},
                 flux::QN,
                 inds::Indices) where {ElT <: Number}
  blocks = nzblocks(flux, IndexSet(inds))
  T = BlockSparseTensor(ElT, blocks, inds)
  return itensor(T)
end

function ITensor(::Type{ElT},
                 flux::QN,
                 inds::Index...) where {ElT <: Number}
  return ITensor(ElT, flux, IndexSet(inds...))
end

ITensor(flux::QN, inds::Indices) = ITensor(Float64, flux, inds)

ITensor(flux::QN,
        inds::Index...) = ITensor(Float64, flux, IndexSet(inds...))

function ITensor(::Type{ElT},
                 inds::QNIndices) where {ElT <: Number}
  return ITensor(ElT, QN(), inds)
end

ITensor(inds::QNIndices) = ITensor(Float64, QN(), inds)

function ITensor(::Type{ElT}, inds::QNIndex...) where {ElT<:Number}
  return ITensor(ElT, QN(), IndexSet(inds...))
end

ITensor(inds::QNIndex...) = ITensor(Float64, QN(), IndexSet(inds...))

"""
    emptyITensor([::Type{ElT} = Float64, ]inds)
    emptyITensor([::Type{ElT} = Float64, ]inds::QNIndex...)

Construct an ITensor with `NDTensors.BlockSparse` storage of element type `ElT` with the no blocks.

If `ElT` is not specified it defaults to `Float64`.

In the future, this will use the storage `NDTensors.EmptyBlockSparse`.
"""
function emptyITensor(::Type{ElT},
                      inds::QNIndices) where {ElT <: Number}
  return itensor(EmptyBlockSparseTensor(ElT, inds))
end

emptyITensor(inds::QNIndices) = emptyITensor(Float64, inds)

"""
    randomITensor([::Type{ElT} = Float64, ][flux::QN = QN(), ]inds)
    randomITensor([::Type{ElT} = Float64, ][flux::QN = QN(), ]inds::Index...)

Construct an ITensor with `NDTensors.BlockSparse` storage filled with random elements of type `ElT` where the nonzero blocks are determined by `flux`.

If `ElT` is not specified it defaults to `Float64`. If the flux is not specified it defaults to `QN()`.
"""
function randomITensor(::Type{ElT},
                       flux::QN,
                       inds::Indices) where {ElT <: Number}
  T = ITensor(ElT, flux, inds)
  randn!(T)
  return T
end

function randomITensor(::Type{ElT},
                       flux::QN,
                       inds::Index...) where {ElT<:Number}
  return randomITensor(ElT, flux, IndexSet(inds...))
end

function randomITensor(::Type{ElT},
                       inds::QNIndices) where {ElT<:Number}
  return randomITensor(ElT, QN(), inds)
end

randomITensor(flux::QN,
              inds::Indices) = randomITensor(Float64, flux, inds)

randomITensor(flux::QN,
              inds::Index...) = randomITensor(Float64,
                                              flux,
                                              IndexSet(inds...))

function randomITensor(::Type{ElT},
                       inds::QNIndex...) where {ElT<:Number}
  return randomITensor(ElT, QN(), IndexSet(inds...))
end

randomITensor(inds::QNIndices) = randomITensor(Float64, QN(), inds)

randomITensor(inds::QNIndex...) = randomITensor(Float64, QN(), IndexSet(inds...))

function combiner(inds::QNIndices; kwargs...)
  # TODO: support combining multiple set of indices
  tags = get(kwargs, :tags, "CMB,Link")
  new_ind = ⊗(inds...)
  if all(i->dir(i)!=Out,inds) && all(i->dir(i)!=In,inds)
    new_ind = dag(new_ind)
    new_ind = replaceqns(new_ind,-qnblocks(new_ind))
  end
  new_ind = settags(new_ind,tags)
  comb_ind,perm,comb = combineblocks(new_ind)
  return itensor(Combiner(perm,comb),
                 IndexSet(comb_ind, dag.(inds)...))
end

#
# DiagBlockSparse ITensor constructors
#

"""
    diagITensor([::Type{ElT} = Float64, ][flux::QN = QN(), ]is)
    diagITensor([::Type{ElT} = Float64, ][flux::QN = QN(), ]is::Index...)

Make an ITensor with storage type `NDTensors.DiagBlockSparse` with elements `zero(ElT)`. The ITensor only has diagonal blocks consistent with the specified `flux`.

If the element type is not specified, it defaults to `Float64`. If theflux is not specified, it defaults to `QN()`.
"""
function diagITensor(::Type{ElT},
                     flux::QN,
                     is::Indices) where {ElT <: Number}
  blocks = nzdiagblocks(flux, IndexSet(is))
  T = DiagBlockSparseTensor(ElT, blocks, is)
  return itensor(T)
end

function diagITensor(::Type{ElT},
                     flux::QN,
                     inds::Index...) where {ElT <: Number}
  return diagITensor(ElT, flux, IndexSet(inds...))
end

diagITensor(flux::QN,
            is::Indices) = diagITensor(Float64, flux, is)

diagITensor(flux::QN,
            inds::Index...) = diagITensor(Float64,
                                          flux,
                                          IndexSet(inds...))

function diagITensor(::Type{ElT},
                     inds::QNIndices) where {ElT <: Number}
  return diagITensor(ElT, QN(), inds)
end

function diagITensor(inds::QNIndices)
  return diagITensor(Float64, QN(), inds)
end

"""
    delta([::Type{ElT} = Float64, ][flux::QN = QN(), ]is)
    delta([::Type{ElT} = Float64, ][flux::QN = QN(), ]is::Index...)

Make an ITensor with storage type `NDTensors.DiagBlockSparse` with uniform elements `one(ElT)`. The ITensor only has diagonal blocks consistent with the specified `flux`.

If the element type is not specified, it defaults to `Float64`. If theflux is not specified, it defaults to `QN()`.
"""
function delta(::Type{ElT},
               flux::QN,
               inds::Indices) where {ElT <: Number}
  blocks = nzdiagblocks(flux, IndexSet(inds))
  T = DiagBlockSparseTensor(one(ElT), blocks, inds)
  return itensor(T)
end

function delta(::Type{ElT},
               flux::QN,
               inds::Index...) where {ElT <: Number}
  return delta(ElT, flux, IndexSet(inds...))
end

delta(flux::QN,
      inds::Indices) = delta(Float64, flux, is)

delta(flux::QN,
      inds::Index...) = delta(Float64, flux, IndexSet(inds...))

function delta(::Type{ElT},
               inds::QNIndices) where {ElT <: Number}
  return delta(ElT, QN(), inds)
end

delta(inds::QNIndices) = delta(Float64, QN(), inds)


text:
# Template for implementing an invertible neural network layer and its logdet.
# We will take an affine layer as an example to explain all necessary steps.
# The affine layer is defined as:
#
# Y = f(X) = S .* X .+ B
#
# Here, X is the 4D input tensor of dimensions (nx, ny, num_channel, batchsize).
# S and B are a scaling and bias term, both with dimensions (nx, ny, num_channel).
#
# The trainable parameters of this layer are S, B. The input is X and the output is Y.
#

# You need to expoert whichever functions you want to be able to access
export AffineLayer

# This immutable structure defines our network layer. The structure contains the
# parameters of our layer (in this case S and B) or any other building blocks that
# you want to use in your layer (for example 1x1 convolutions). However, in this
# case we only have parameters S and B.
struct AffineLayer <: NeuralNetLayer
    S::Parameter    # trainable parameters are defined as Parameters.
    B::Parameter    # both S and B have two fields: S.data and S.grad
    logdet::Bool    # bool to indicate whether you want to compute the logdet
end

# Functor the layer for gpu/cpu offloading
@Flux.functor AffineLayer

# The constructor builds and returns a new network layer for given input dimensions.
function AffineLayer(nx, ny, nc; logdet=false)

    # Create S and B
    S = Parameter(glorot_uniform(nx, ny, nc))   # initiliaze S with random values and make it a parameter
    B = Parameter(zeros(Float32, nx, ny, nc))   # initilize B with zeros and make it a parameter

    # Build an affine layer
    return AffineLayer(S, B, logdet)
end

# Foward pass: Input X, Output Y
# The forward pass for the affine layer is:
# Y = X .* S .+ B
function forward(X::AbstractArray{T, N}, AL::AffineLayer) where {T, N}

    Y = X .* AL.S.data .+ AL.B.data   # S and B are Parameters, so access their values via the .data field

    # If logdet is true, also compute the logdet and return it as a second output argument.
    # Otherwise only return Y.
    if AL.logdet == true
        return Y, logdet_forward(S)
    else
        return Y
    end
end

# Inverse pass: Input Y, Output X
# The inverse pass for our affine layer is:
# X = (Y .- B) ./ S
# To avoid division by zero, we add numerical noise to S in the division.
function inverse(Y::AbstractArray{T, N}, AL::AffineLayer) where {T, N}
    X = (Y .- AL.B.data) ./ (AL.S.data .+ eps(T))   # avoid division by 0
    return X
end

# Backward pass: Input (ΔY, Y), Output (ΔY, Y)
# Assuming that the layer is invertible, the backward function takes
# 2 input arguments: the data residual ΔY and the original output Y.
# The first step is to recompute the original input X by calling the
# inverse function on Y. If the layer is not invertible, this layer
# needs X as an input instead of Y.
# Second of all, we compute the partial derivatives of our layer
# with respect to X, S, and B:
# ΔX = S .* ΔY (corresponds to df/dX * dY)
# ΔS = X .* ΔY (corresponds to df/dS * dY)
# ΔB = 1 .* ΔY (corresponds to df/dB * dY)
function backward(ΔY::AbstractArray{T, N}, Y::AbstractArray{T, N}, AL::AffineLayer) where {T, N}
    nx, ny, n_in, batchsize = size(Y)

    # Recompute X from Y
    X = inverse(Y, AL)

    # Gradient w.r.t. X
    ΔX = ΔY .* AL.S.data

    # Gradient w.r.t. S (sum over the batchsize, as S is only a 3D tensor)
    ΔS = sum(ΔY .* X, dims=4)[:,:,:,1]

    # If the logdet is computed, in the forward pass, also compute the gradient of the
    # logdet term and subtract from S (assuming that the logdet term is subtracted in the
    # objective function)
    AL.logdet == true && (ΔS -= logdet_backward(S))

    # Gradient w.r.t. B (sum over the batchsize)
    ΔB = sum(ΔY, dims=4)[:,:,:,1]

    # Set the gradient fields of the parameters S and B
    AL.S.grad = ΔS
    AL.B.grad = ΔB

    # Return the backpropagated data residual and the re-computed X
    return ΔX, X
end

# For optimization, we need a function that clears all the gradients.
# I.e. we set the .grad fields of all parameters to nothing.
function clear_grad!(AL::AffineLayer)
    AL.S.grad = nothing
    AL.B.grad = nothing
end

# Also we define a get_params function that returns an array of all
# the parameters. In this case, our parameters are S and B
get_params(AL::AffineLayer) = [AL.S, AL.B]

# Function for the logdet and for computing the gradient of the logdet.
# For our affine layer consisting of an element-wise multiplication of S
# and X, the Jacobian is given by S, and the logdet is the sum of the logarithm
# of the (absolute) values.
logdet_forward(S) = sum(log.(abs.(S.data)))

# The gradient of the forward logdet function is given by 1/S
logdet_backward(S) = 1f0 ./ S.data


text:
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("color codes print in REPL") do

    include(Pkg.dir("Watchmaker/src/colors.jl")) # pulls in palette

    context("palette list is intact") do
        @fact typeof(palette)                                  --> Array{ASCIIString,1}
        @fact length(palette)                                  --> 661
        @fact sum([contains(p, "goldenrod") for p in palette]) --> 17
    end
end


text:
import Tables
import TableOperations
import SentinelArrays

import Base: fetch, show, length

export DTable, tabletype, tabletype!, trim, trim!, leftjoin, innerjoin

const VTYPE = Vector{Union{Dagger.Chunk,Dagger.EagerThunk}}

"""
    DTable

Structure representing the distributed table based on Dagger.

The table is stored as a vector of `Chunk` structures which hold partitions of the table.
That vector can also store `EagerThunk` structures when an operation that modifies
the underlying partitions was applied to it (currently only `filter`).
"""
mutable struct DTable
    chunks::VTYPE
    tabletype
    schema::Union{Nothing,Tables.Schema}
    DTable(chunks::VTYPE, tabletype) = new(chunks, tabletype, nothing)
end

DTable(chunks::Vector{Dagger.EagerThunk}, args...) = DTable(VTYPE(chunks), args...)
DTable(chunks::Vector{Dagger.Chunk}, args...) = DTable(VTYPE(chunks), args...)

"""
    DTable(table; tabletype=nothing) -> DTable

Constructs a `DTable` using a `Tables.jl`-compatible input `table`.
Calls `Tables.partitions` on `table` and assumes the provided partitioning.
"""
function DTable(table; tabletype=nothing)
    chunks = Vector{Dagger.Chunk}()
    type = nothing
    sink = nothing
    for partition in Tables.partitions(table)
        if sink === nothing
            sink = Tables.materializer(tabletype !== nothing ? tabletype() : partition)
        end

        tpart = sink(partition)
        push!(chunks, Dagger.tochunk(tpart))

        if type === nothing
            type = typeof(tpart).name.wrapper
        end
    end
    return DTable(chunks, type)
end

"""
    DTable(table, chunksize; tabletype=nothing) -> DTable

Constructs a `DTable` using a `Tables.jl` compatible `table` input.
It assumes no initial partitioning of the table and uses the `chunksize`
argument to partition the table (based on row count).

Providing `tabletype` kwarg overrides the internal table partition type.
"""
function DTable(table, chunksize::Integer; tabletype=nothing)
    chunks = Vector{Dagger.Chunk}()
    type = nothing
    sink = Tables.materializer(tabletype !== nothing ? tabletype() : table)
    for outer_partition in Tables.partitions(table)
        for inner_partition in Tables.partitions(TableOperations.makepartitions(outer_partition, chunksize))
            tpart = sink(inner_partition)
            push!(chunks, Dagger.tochunk(tpart))
            if type === nothing
                type = typeof(tpart).name.wrapper
            end
        end
    end
    return DTable(chunks, type)
end

"""
    DTable(loader_function, files::Vector{String}; tabletype=nothing)

Constructs a `DTable` using a list of filenames and a `loader_function`.
Partitioning is based on the contents of the files provided, which means that
one file is used to create one partition.

Providing `tabletype` kwarg overrides the internal table partition type.
"""
function DTable(loader_function, files::Vector{String}; tabletype=nothing)
    chunks = Vector{Dagger.EagerThunk}()
    sizehint!(chunks, length(files))

    append!(chunks, map(file -> Dagger.spawn(_file_load, file, loader_function, tabletype), files))

    return DTable(chunks, tabletype)
end

function _file_load(filename, loader_function, tabletype)
    part = loader_function(filename)
    sink = Tables.materializer(tabletype === nothing ? part : tabletype())
    tpart = sink(part)
    return tpart
end

"""
    fetch(d::DTable)

Collects all the chunks in the `DTable` into a single, non-distributed
instance of the underlying table type.

Fetching an empty DTable results in returning an empty `NamedTuple` regardless of the underlying `tabletype`.
"""
function fetch(d::DTable)
    sink = Tables.materializer(tabletype(d)())
    sink(_retrieve_partitions(d))
end

"""
    fetch(d::DTable, sink)

Collects all the chunks in the `DTable` into a single, non-distributed
instance of table type created using the provided `sink` function.
"""
fetch(d::DTable, sink) = sink(_retrieve_partitions(d))

function _retrieve_partitions(d::DTable)
    d2 = trim(d)
    return nchunks(d2) > 0 ?
        TableOperations.joinpartitions(Tables.partitioner(_retrieve, d2.chunks)) : NamedTuple()
end

_retrieve(x::Dagger.EagerThunk) = fetch(x)
_retrieve(x::Dagger.Chunk) = collect(x)

"""
    tabletype!(d::DTable)

Provides the type of the underlying table partition and caches it in `d`.

In case the tabletype cannot be obtained the default return value is `NamedTuple`.
"""
tabletype!(d::DTable) = d.tabletype = resolve_tabletype(d)

"""
    tabletype(d::DTable)

Provides the type of the underlying table partition.
Uses the cached tabletype if available.

In case the tabletype cannot be obtained the default return value is `NamedTuple`.
"""
tabletype(d::DTable) = d.tabletype === nothing ? resolve_tabletype(d) : d.tabletype

function resolve_tabletype(d::DTable)
    _type = c -> isnonempty(c) ? typeof(c).name.wrapper : nothing
    t = nothing

    if nchunks(d) > 0
        for chunk in d.chunks
            t = fetch(Dagger.@spawn _type(chunk))
            t !== nothing && break
        end
    end
    t !== nothing ? t : NamedTuple
end

function isnonempty(chunk)
    length(Tables.rows(chunk)) > 0 && length(Tables.columnnames(chunk)) > 0
end

"""
    trim!(d::DTable) -> DTable

Removes empty chunks from `d`.
"""
function trim!(d::DTable)
    check_result = [Dagger.@spawn isnonempty(c) for c in d.chunks]
    d.chunks = getindex.(filter(x -> fetch(check_result[x[1]]), collect(enumerate(d.chunks))), 2)
    d
end

"""
    trim(d::DTable) -> DTable

Returns `d` with empty chunks removed.
"""
trim(d::DTable) = trim!(DTable(d.chunks, d.tabletype))

show(io::IO, d::DTable) = show(io, MIME"text/plain"(), d)

function show(io::IO, ::MIME"text/plain", d::DTable)
    tabletype = d.tabletype === nothing ? "unknown (use `tabletype!(::DTable)`)" : d.tabletype
    println(io, "DTable with $(nchunks(d)) partitions")
    print(io, "Tabletype: $tabletype")
    nothing
end

function length(table::DTable)
    f = x -> length(Tables.rows(x))
    sum(fetch.([Dagger.@spawn f(c) for c in table.chunks]))
end

function _columnnames_svector(d::DTable)
    colnames_tuple = determine_columnnames(d)
    colnames_tuple !== nothing ? [sym for sym in colnames_tuple] : nothing
end

@inline nchunks(d::DTable) = length(d.chunks)

merge_chunks(sink, chunks) = sink(TableOperations.joinpartitions(Tables.partitioner(_retrieve, chunks)))


text:
# This page was generated on DATEOFTODAY

using GeoPhyInv
using Test
using Random
#jl @init_parallel_stencil(3, false, Float32, 4)

# # Intro

#md # ```@docs
#md # Records
#md # ```

# # Examples
# Define an acquisition geometry.
mgrid = [range(0.0, stop=10.,step=0.1), range(0.0, stop=30.,step=0.2)];
ageom=AGeom(mgrid, SSrcs(10), Srcs(10), Recs(10));

# Need a time grid.
tgrid=range(0, stop=1.0, step=0.1);

# Lets initialize records for `:P` and `:vz` fields.
data=Records(tgrid, ageom, [:p,:vz]);

# Fill the `:P` field of 3rd supersource with random numbers.
Random.randn!(data[3][:p]);

# # Methods
# Methods listed for `SrcWav` can used on instances of `Records` too.


text:
abstract type AbstractPrimRec end

struct Succ <: AbstractPrimRec end

(::Succ)(xs::Int...)::Int = xs[1] + 1

const succ = Succ()

struct Zero <: AbstractPrimRec end

(::Zero)(xs::Int...)::Int = 0

const zro = Zero()

struct Proj{T} <: AbstractPrimRec end

function (::Proj{T})(xs::Int...)::Int where {T}; xs[T] end

struct Comb <: AbstractPrimRec
    f::AbstractPrimRec
    gs::Tuple{Vararg{T where T<:AbstractPrimRec}}
    Comb(f, gs...) = new(f, gs)
end

(c::Comb)(xs::Int...)::Int = c.f((g(xs...) for g in c.gs)...)

struct PrimRec <: AbstractPrimRec
    basic::AbstractPrimRec
    step::AbstractPrimRec
end

(pr::PrimRec)(xs::Int...)::Int = begin
    t = xs[1]
    y = xs[2:end]
    cur = pr.basic(y...)
    for i in 0:(t - 1)
        cur = pr.step(cur, i, y...)
    end
    cur
end

(f::AbstractPrimRec)(gs::AbstractPrimRec...) = Comb(f, gs...)


text:
import JSON
import UUIDs: UUID
import HTTP
import Sockets

connectedclients = Dict{Symbol,Client}()
notebooks = Dict{UUID,Notebook}()

"Send `messages` to all clients connected to the `notebook`."
function putnotebookupdates!(notebook::Notebook, messages::UpdateMessage...)
    listeners = filter(collect(values(connectedclients))) do c
        c.connected_notebook !== nothing &&
        c.connected_notebook.uuid == notebook.uuid
    end
    if isempty(listeners)
        @info "no clients connected to this notebook!"
    else
        for next_to_send in messages, client in listeners
            put!(client.pendingupdates, next_to_send)
        end
    end
    flushallclients(listeners)
    listeners
end

"Send `messages` to all connected clients."
function putplutoupdates!(messages::UpdateMessage...)
    listeners = collect(values(connectedclients))
    if isempty(listeners)
        @info "no clients connected to pluto!"
    else
        for next_to_send in messages, client in listeners
            put!(client.pendingupdates, next_to_send)
        end
    end
    flushallclients(listeners)
    listeners
end

"Send `messages` to a `client`."
function putclientupdates!(client::Client, messages::UpdateMessage...)
    for next_to_send in messages
        put!(client.pendingupdates, next_to_send)
    end
    flushclient(client)
end

"Send `messages` to the `Client` who initiated."
function putclientupdates!(initiator::Initiator, messages::UpdateMessage...)
    putclientupdates!(connectedclients[initiator.clientID], messages...)
end

# flushtoken = Channel{Nothing}(1)
# put!(flushtoken, nothing)

function flushclient(client::Client)
    # take!(flushtoken)
    didsomething = false
    while isready(client.pendingupdates)
        next_to_send = take!(client.pendingupdates)
        didsomething = true

        token = take!(client.stream_accesstoken)
        try
            if client.stream !== nothing
                if isopen(client.stream)
                    write(client.stream, serialize_message(next_to_send))
                else
                    @info "Client $(client.id) stream closed."
                    return false
                end
            end
        catch ex
            bt = stacktrace(catch_backtrace())
            @warn "Failed to write to WebSocket of $(client.id) " exception=(ex,bt)
            return false
        end
        put!(client.stream_accesstoken, token)
    end
    # put!(flushtoken, nothing)
    true
end

function flushallclients(subset::Union{Set{Client}, AbstractArray{Client}})
    disconnected = Set{Symbol}()
    for client in subset
        stillconnected = flushclient(client)
        if !stillconnected
            push!(disconnected, client.id)
        end
    end
    for to_deleteID in disconnected
        delete!(connectedclients, to_deleteID)
    end
end

function flushallclients()
    flushallclients(values(connectedclients))
end

"Will hold all 'response handlers': functions that respond to a WebSocket request from the client. These are defined in `src/webserver/Dynamic.jl`."
responses = Dict{Symbol,Function}()

"""Start a Pluto server _synchronously_ (i.e. blocking call) on `http://localhost:[port]/`.

This will start the static HTTP server and a WebSocket server. Pluto Notebooks will be started dynamically (by user input)."""
function run(port = 1234, launchbrowser = false)
    serversocket = Sockets.listen(UInt16(port))
    @async HTTP.serve(Sockets.localhost, UInt16(port), stream = true, server = serversocket) do http::HTTP.Stream
        # messy messy code so that we can use the websocket on the same port as the HTTP server

        if HTTP.WebSockets.is_upgrade(http.message)
            try
                HTTP.WebSockets.upgrade(http) do clientstream
                    if !isopen(clientstream)
                        return
                    end
                    while !eof(clientstream)
                        # This stream contains data received over the WebSocket.
                        # It is formatted and JSON-encoded by send(...) in editor.html
                        try
                            parentbody = let
                                data = String(readavailable(clientstream))
                                JSON.parse(data)
                            end

                            client = let
                                clientID = Symbol(parentbody["clientID"])
                                Client(clientID, clientstream)
                            end
                            # add to our list of connected clients
                            connectedclients[client.id] = client

                            messagetype = Symbol(parentbody["type"])
                            requestID = Symbol(parentbody["requestID"])

                            if messagetype == :disconnect
                                delete!(connectedclients, client.id)
                                close(clientstream)
                            else
                                body = parentbody["body"]

                                args = []
                                if haskey(parentbody, "notebookID")
                                    notebook = let
                                        notebookID = UUID(parentbody["notebookID"])
                                        get(notebooks, notebookID, nothing)
                                    end
                                    if notebook === nothing
                                        messagetype === :connect || @warn "Remote notebook not found locally!"
                                    else
                                        client.connected_notebook = notebook
                                        push!(args, notebook)
                                    end
                                end

                                if haskey(parentbody, "cellID")
                                    cell = let
                                        cellID = UUID(parentbody["cellID"])
                                        selectcell_byuuid(notebook, cellID)
                                    end
                                    if cell === nothing
                                        @warn "Remote cell not found locally!"
                                    else
                                        push!(args, cell)
                                    end
                                end

                                if haskey(responses, messagetype)
                                    responsefunc = responses[messagetype]
                                    try
                                        responsefunc(body, args..., initiator=Initiator(client.id, requestID))
                                    catch ex
                                        @warn "Response function to message of type $(messagetype) failed"
                                        rethrow(ex)
                                    end
                                else
                                    @warn "Message of type $(messagetype) not recognised"
                                end
                            end
                        catch ex
                            if ex isa InterruptException
                                rethrow(ex)
                            elseif ex isa HTTP.WebSockets.WebSocketError
                                # that's fine!
                            elseif ex isa InexactError
                                # that's fine! this is a (fixed) HTTP.jl bug: https://github.com/JuliaWeb/HTTP.jl/issues/471
                                # TODO: remove this switch
                            else
                                bt = stacktrace(catch_backtrace())
                                @warn "Reading WebSocket client stream failed for unknown reason:" exception=(ex,bt)
                            end
                        end
                    end
                end
            catch ex
                if ex isa InterruptException
                    rethrow(ex)
                elseif ex isa Base.IOError
                    # that's fine!
                elseif ex isa ArgumentError && occursin("stream is closed", ex.msg)
                    # that's fine!
                else
                    bt = stacktrace(catch_backtrace())
                    @warn "HTTP upgrade failed for unknown reason" exception=(ex,bt)
                end
            end
        else
            request::HTTP.Request = http.message
            request.body = read(http)
            HTTP.closeread(http)

            request_body = IOBuffer(HTTP.payload(request))
            if eof(request_body)
                # no request body
                response_body = HTTP.handle(PLUTOROUTER, request)
            else
                # there's a body, so pass it on to the handler we dispatch to
                response_body = HTTP.handle(PLUTOROUTER, request, JSON.parse(request_body))
            end

            request.response::HTTP.Response = response_body
            request.response.request = request
            try
                HTTP.startwrite(http)
                write(http, request.response.body)
            catch e
                if isa(e, HTTP.IOError) || isa(e, ArgumentError)
                    @warn "Attempted to write to a closed stream at $(request.target)"
                else
                    rethrow(e)
                end
            end
        end
    end

    println("Go to http://localhost:$(port)/ to start writing! ⚙")
    println()
    controlkey = Sys.isapple() ? "Command" : "Ctrl"
    println("Press $controlkey+C to stop Pluto")
    println()

    launchbrowser && @warn "Not implemented yet"

    # create blocking call:
    try
        while true
            sleep(typemax(UInt64))
        end
    catch e
        if isa(e, InterruptException)
            println("\nClosing Pluto... Restart Julia for a fresh session. \n\nHave a nice day! 🎈")
            close(serversocket)
            for (uuid, ws) in WorkspaceManager.workspaces
                WorkspaceManager.unmake_workspace(ws)
            end
        else
            rethrow(e)
        end
    end
end


text:
# Author: Sylvain Barbot (sbarbot@ntu.edu.sg) - May 19, 2018, Los Angeles.
# Translated by Pengcheng Shi (shipengcheng1230@gmail.com) 08/2019

# greens function ϵ vs u for tet4 element
export _disp_vol_tet4, _disp_vol_tet4!
export _strain_vol_tet4, _strain_vol_tet4!
export _stress_vol_tet4, _stress_vol_tet4!

export disp_vol_tet4, disp_vol_tet4!
export stress_vol_tet4, stress_vol_tet4!

function _swap_coord!(x::AbstractVector)
    x[1], x[2], x[3] = x[2], x[1], -x[3]
end

function disp_vol_tet4!(u::AbstractVector{R}, quadrature::Q, x::R, y::R, z::R, A::U, B::U, C::U, D::U, ϵxx::R, ϵxy::R, ϵxz::R, ϵyy::R, ϵyz::R, ϵzz::R, ν::R) where {R, Q, U}
    foreach(x -> _swap_coord!(x), [A, B, C, D])
    _disp_vol_tet4!(u, quadrature, y, x, -z, A, B, C, D, ϵyy, ϵxy, -ϵyz, ϵxx, -ϵxz, ϵzz, ν)
    foreach(x -> _swap_coord!(x), [A, B, C, D])
    u[1], u[2], u[3] = u[2], u[1], -u[3]
end

function disp_vol_tet4(quadrature::Q, x::R, y::R, z::R, A::U, B::U, C::U, D::U, ϵxx::R, ϵxy::R, ϵxz::R, ϵyy::R, ϵyz::R, ϵzz::R, ν::R) where {R, Q, U}
    u = Vector{R}(undef, 3)
    disp_vol_tet4!(u, quadrature, x, y, z, A, B, C, D, ϵxx, ϵxy, ϵxz, ϵyy, ϵyz, ϵzz, ν)
    return u
end

function stress_vol_tet4!(σ::AbstractVector{R}, quadrature::Q, x::R, y::R, z::R, A::U, B::U, C::U, D::U, ϵxx::R, ϵxy::R, ϵxz::R, ϵyy::R, ϵyz::R, ϵzz::R, G::R, ν::R) where {R, Q, U}
    foreach(x -> _swap_coord!(x), [A, B, C, D])
    _stress_vol_tet4!(σ, quadrature, y, x, -z, A, B, C, D, ϵyy, ϵxy, -ϵyz, ϵxx, -ϵxz, ϵzz, G, ν)
    foreach(x -> _swap_coord!(x), [A, B, C, D])
    σ[1], σ[4] = σ[4], σ[1]
    σ[3], σ[5] = -σ[5], -σ[3]
end

function stress_vol_tet4(quadrature::Q, x::R, y::R, z::R, A::U, B::U, C::U, D::U, ϵxx::R, ϵxy::R, ϵxz::R, ϵyy::R, ϵyz::R, ϵzz::R, G::R, ν::R) where {R, Q, U}
    σ = Vector{R}(undef, 6)
    stress_vol_tet4!(σ, quadrature, x, y, z, A, B, C, D, ϵxx, ϵxy, ϵxz, ϵyy, ϵyz, ϵzz, G, ν)
    return σ
end

@doc raw"""
    _disp_vol_tet4(quadrature::Q,
        x1::R, x2::R, x3::R, A::U, B::U, C::U, D::U,
        e11::R, e12::R, e13::R, e22::R, e23::R, e33::R, nu::R
        ) where {R, U, Q}

Compute displacement arisen from inelastic strain in Tet4 elements.
    Please see [original version](https://bitbucket.org/sbarbot/bssa-2018058/src/default/)
    for complete details, especially the **coordinate system** used here.

## Arguments
- `quadrature`: quadrature rule for integration, see
    [FastGaussQuadrature.jl](https://github.com/JuliaApproximation/FastGaussQuadrature.jl)
- `x1`, `x2`, `x3`: observational position, where ``x_{3} ≥ 0``
- `A`, `B`, `C`, `D`: a list of 3 numbers for each, each of which represents
    coordinates of the vertex. All depth coordinates must be greater or
    equal to 0 (no checking is performed here)
- `e**`: strain components, each is ``ϵ_{11}``, ``ϵ_{12}``, ``ϵ_{13}``,
    ``ϵ_{22}``, ``ϵ_{23}``, ``ϵ_{33}``
- `nu`: poisson ratio

## Output
A vector of 3 numbers, each represents ``u_{1}``, ``u_{2}``, ``u_{3}``

## Notice
- Inplace version: `_disp_vol_tet4!(u, args...)` where u is a vector of 3 numbers.
"""
function _disp_vol_tet4(quadrature::Q,
    x1::R, x2::R, x3::R, A::U, B::U, C::U, D::U,
    e11::R, e12::R, e13::R, e22::R, e23::R, e33::R, nu::R
    ) where {R, U, Q}

    u = Vector{R}(undef, 3)
    _disp_vol_tet4!(u, quadrature, x1, x2, x3, A, B, C, D, e11, e12, e13, e22, e23, e33, nu)
    return u
end

function _disp_vol_tet4!(u::W, quadrature::Q,
    x1::R, x2::R, x3::R, A::U, B::U, C::U, D::U,
    e11::R, e12::R, e13::R, e22::R, e23::R, e33::R, nu::R
    ) where {R, U, Q, W}

    @assert x3 ≥ zero(R) "Depth must be greater than or equal to 0."
    lambda = 2 * nu / (1 - 2 * nu)
    ekk = e11 + e22 + e33

    nA = cross(C - B, D - B)
    nB = cross(D - C, A - C)
    nC = cross(A - D, B - D)
    nD = cross(B - A, C - A)

    nA /= norm(nA)
    nB /= norm(nB)
    nC /= norm(nC)
    nD /= norm(nD)

    if nA' * (A .- (B .+ C .+ D) ./ 3) > zero(R); nA *= -one(R) end
    if nB' * (B .- (A .+ C .+ D) ./ 3) > zero(R); nB *= -one(R) end
    if nC' * (C .- (B .+ A .+ D) ./ 3) > zero(R); nC *= -one(R) end
    if nD' * (D .- (B .+ C .+ A) ./ 3) > zero(R); nD *= -one(R) end

    ABC = norm(cross(C .- A, B .- A)) / 2
    BCD = norm(cross(D .- B, C .- B)) / 2
    CDA = norm(cross(A .- C, D .- C)) / 2
    DAB = norm(cross(B .- D, A .- D)) / 2

    m11 = lambda * ekk + 2 * e11
    m12 = 2 * e12
    m13 = 2 * e13
    m22 = lambda * ekk + 2 * e22
    m23 = 2 * e23
    m33 = lambda * ekk + 2 * e33

    let lambda=lambda, x1=x1, x2=x2, x3=x3, A=A, B=B, C=C, D=D, e11=e11, e12=e12, e13=e13, e22=e22, e23=e23, e33=e33, nu=nu, nA=nA, nB=nB, nC=nC, ABC=ABC, BCD=BCD, CDA=CDA, DAB=DAB, m11=m11, m12=m12, m13=m13, m22=m22, m23=m23, m33=m33

        function r1(y1::R, y2::R, y3::R) where R
            hypot(x1 - y1, x2 - y2, x3 - y3)
        end

        function r2(y1::R, y2::R, y3::R) where R
            hypot(x1 - y1, x2 - y2, x3 + y3)
        end

        function G11(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (
            (3 - 4 * nu) / r1(y1, y2, y3) + 1 / r2(y1, y2, y3) + (x1 - y1) ^ 2 / r1(y1, y2, y3) ^ 3
            + (3 - 4 * nu) * (x1 - y1) ^ 2 / r2(y1, y2, y3) ^ 3 + 2 * x3 * y3 * (r2(y1, y2, y3) ^ 2 - 3 * (x1 - y1) ^ 2) / r2(y1, y2, y3) ^ 5
            + 4 * (1 - 2 * nu) * (1 - nu) * (r2(y1, y2, y3) ^ 2 - (x1 - y1) ^ 2 + r2(y1, y2, y3) * (x3 + y3)) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3) ^ 2)))
        end

        function G12(y1::R, y2::R, y3::R) where R
            ((x1 - y1) * (x2 - y2) / (16 * π * (1 - nu)) * (
            1  / r1(y1, y2, y3) ^ 3 + (3 - 4 * nu) / r2(y1, y2, y3) ^ 3 - 6 * x3 * y3 / r2(y1, y2, y3) ^ 5
            -4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3) ^ 2)))
        end

        function G13(y1::R, y2::R, y3::R) where R
            ((x1 - y1) /(16*π*(1-nu)) *(
              (x3 - y3) / r1(y1, y2, y3) ^ 3 + (3 - 4 * nu) * (x3 - y3) / r2(y1, y2, y3) ^ 3
            -6 * x3 * y3 * (x3 + y3) / r2(y1, y2, y3) ^ 5 + 4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3))))
        end

        function G21(y1::R, y2::R, y3::R) where R
            ((x1 - y1) * (x2 - y2) / (16 * π * (1 - nu)) * (
            1  / r1(y1, y2, y3) ^ 3 + (3 - 4 * nu) / r2(y1, y2, y3) ^ 3 - 6 * x3 * y3 / r2(y1, y2, y3) ^ 5
            -4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3) ^ 2)))
        end

        function G22(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (
            (3 - 4 * nu) / r1(y1, y2, y3) + 1  / r2(y1, y2, y3) + (x2 - y2) ^ 2  / r1(y1, y2, y3) ^ 3
            +(3 - 4 * nu) * (x2 - y2) ^ 2  / r2(y1, y2, y3) ^ 3 + 2 * x3 * y3 * (r2(y1, y2, y3) ^ 2 - 3 * (x2 - y2) ^ 2) / r2(y1, y2, y3) ^ 5
            +4 * (1 - 2 * nu) * (1 - nu) * (r2(y1, y2, y3) ^ 2 - (x2 - y2) ^ 2 + r2(y1, y2, y3) * (x3 + y3)) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3) ^ 2)))
        end

        function G23(y1::R, y2::R, y3::R) where R
            ((x2 - y2) /(16*π*(1-nu)) *(
              (x3 - y3) / r1(y1, y2, y3) ^ 3 + (3 - 4 * nu) * (x3 - y3) / r2(y1, y2, y3) ^ 3
            -6 * x3 * y3 * (x3 + y3) / r2(y1, y2, y3) ^ 5 + 4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3))))
        end

        function G31(y1::R, y2::R, y3::R) where R
            ((x1 - y1) /(16*π*(1-nu)) *(
              (x3 - y3) / r1(y1, y2, y3) ^ 3 + (3 - 4 * nu) * (x3 - y3) / r2(y1, y2, y3) ^ 3
            +6 * x3 * y3 * (x3 + y3) / r2(y1, y2, y3) ^ 5 - 4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3))))
        end

        function G32(y1::R, y2::R, y3::R) where R
            ((x2 - y2) /(16*π*(1-nu)) *(
              (x3 - y3) / r1(y1, y2, y3) ^ 3 + (3 - 4 * nu) * (x3 - y3) / r2(y1, y2, y3) ^ 3
            +6 * x3 * y3 * (x3 + y3) / r2(y1, y2, y3) ^ 5 - 4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3))))
        end

        function G33(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (
            (3 - 4 * nu) / r1(y1, y2, y3) + (5 - 12 * nu + 8 * nu ^ 2) / r2(y1, y2, y3) + (x3 - y3) ^ 2  / r1(y1, y2, y3) ^ 3
            +6 * x3 * y3 * (x3 + y3) ^ 2  / r2(y1, y2, y3) ^ 5 + ((3 - 4 * nu) * (x3 + y3) ^ 2 - 2 * x3 * y3) / r2(y1, y2, y3) ^ 3))
        end

        function y(u::R, v::R, A::R, B::R, C::R) where R
            A * (1 - u) * (1 - v) / 4 + B * (1 + u) * (1 - v) / 4 + C * (1 + v) / 2
        end

        function IU1(u::R, v::R) where R
            (ABC / 4 * (m11 * nD[1] + m12 * nD[2] + m13 * nD[3]) * G11(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
                + ABC / 4 * (m12 * nD[1] + m22 * nD[2] + m23 * nD[3]) * G21(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
                + ABC / 4 * (m13 * nD[1] + m23 * nD[2] + m33 * nD[3]) * G31(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
                + BCD / 4 * (m11 * nA[1] + m12 * nA[2] + m13 * nA[3]) * G11(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
                + BCD / 4 * (m12 * nA[1] + m22 * nA[2] + m23 * nA[3]) * G21(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
                + BCD / 4 * (m13 * nA[1] + m23 * nA[2] + m33 * nA[3]) * G31(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
                + CDA / 4 * (m11 * nB[1] + m12 * nB[2] + m13 * nB[3]) * G11(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
                + CDA / 4 * (m12 * nB[1] + m22 * nB[2] + m23 * nB[3]) * G21(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
                + CDA / 4 * (m13 * nB[1] + m23 * nB[2] + m33 * nB[3]) * G31(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
                + DAB / 4 * (m11 * nC[1] + m12 * nC[2] + m13 * nC[3]) * G11(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
                + DAB / 4 * (m12 * nC[1] + m22 * nC[2] + m23 * nC[3]) * G21(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
                + DAB / 4 * (m13 * nC[1] + m23 * nC[2] + m33 * nC[3]) * G31(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3])))
        end

        function IU2(u::R, v::R) where R
            (ABC / 4 * (m11 * nD[1] + m12 * nD[2] + m13 * nD[3]) * G12(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
                + ABC / 4 * (m12 * nD[1] + m22 * nD[2] + m23 * nD[3]) * G22(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
                + ABC / 4 * (m13 * nD[1] + m23 * nD[2] + m33 * nD[3]) * G32(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
                + BCD / 4 * (m11 * nA[1] + m12 * nA[2] + m13 * nA[3]) * G12(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
                + BCD / 4 * (m12 * nA[1] + m22 * nA[2] + m23 * nA[3]) * G22(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
                + BCD / 4 * (m13 * nA[1] + m23 * nA[2] + m33 * nA[3]) * G32(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
                + CDA / 4 * (m11 * nB[1] + m12 * nB[2] + m13 * nB[3]) * G12(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
                + CDA / 4 * (m12 * nB[1] + m22 * nB[2] + m23 * nB[3]) * G22(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
                + CDA / 4 * (m13 * nB[1] + m23 * nB[2] + m33 * nB[3]) * G32(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
                + DAB / 4 * (m11 * nC[1] + m12 * nC[2] + m13 * nC[3]) * G12(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
                + DAB / 4 * (m12 * nC[1] + m22 * nC[2] + m23 * nC[3]) * G22(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
                + DAB / 4 * (m13 * nC[1] + m23 * nC[2] + m33 * nC[3]) * G32(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3])))
        end

        function IU3(u::R, v::R) where R
            (ABC / 4 * (m11 * nD[1] + m12 * nD[2] + m13 * nD[3]) * G13(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
                + ABC / 4 * (m12 * nD[1] + m22 * nD[2] + m23 * nD[3]) * G23(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
                + ABC / 4 * (m13 * nD[1] + m23 * nD[2] + m33 * nD[3]) * G33(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
                + BCD / 4 * (m11 * nA[1] + m12 * nA[2] + m13 * nA[3]) * G13(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
                + BCD / 4 * (m12 * nA[1] + m22 * nA[2] + m23 * nA[3]) * G23(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
                + BCD / 4 * (m13 * nA[1] + m23 * nA[2] + m33 * nA[3]) * G33(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
                + CDA / 4 * (m11 * nB[1] + m12 * nB[2] + m13 * nB[3]) * G13(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
                + CDA / 4 * (m12 * nB[1] + m22 * nB[2] + m23 * nB[3]) * G23(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
                + CDA / 4 * (m13 * nB[1] + m23 * nB[2] + m33 * nB[3]) * G33(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
                + DAB / 4 * (m11 * nC[1] + m12 * nC[2] + m13 * nC[3]) * G13(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
                + DAB / 4 * (m12 * nC[1] + m22 * nC[2] + m23 * nC[3]) * G23(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
                + DAB / 4 * (m13 * nC[1] + m23 * nC[2] + m33 * nC[3]) * G33(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3])))
        end

        fill!(u, zero(R))
        x, w = quadrature
        @inbounds @fastmath for k = 1: length(x)
            @simd for j = 1: length(x)
                u[1] += w[j] * w[k] * (1 - x[k]) * IU1(x[j], x[k])
                u[2] += w[j] * w[k] * (1 - x[k]) * IU2(x[j], x[k])
                u[3] += w[j] * w[k] * (1 - x[k]) * IU3(x[j], x[k])
            end
        end

        return nothing
    end
end

@doc raw"""
    _stress_vol_tet4(quadrature::Q,
        x1::R, x2::R, x3::R, A::U, B::U, C::U, D::U,
        e11::R, e12::R, e13::R, e22::R, e23::R, e33::R, G::R, nu::R
        ) where {R, U, Q}

Compute stress arisen from inelastic strain in Tet4 elements.
    Please see [original version](https://bitbucket.org/sbarbot/bssa-2018058/src/default/)
    for complete details, especially the **coordinate system** used here.

## Arguments
- `quadrature`: quadrature rule for integration, see
    [FastGaussQuadrature.jl](https://github.com/JuliaApproximation/FastGaussQuadrature.jl)
- `x1`, `x2`, `x3`: observational position
- `A`, `B`, `C`, `D`: a list of 3 numbers for each, each of which represents
    coordinates of the vertex
- `e**`: strain components, each is ``ϵ_{11}``, ``ϵ_{12}``, ``ϵ_{13}``,
    ``ϵ_{22}``, ``ϵ_{23}``, ``ϵ_{33}``
- `G`: shear modulus
- `nu`: poisson ratio

## Output
A vector of 6 numbers, each represents ``σ_{11}``, ``σ_{12}``, ``σ_{13}``,
    ``σ_{22}``, ``σ_{23}``, ``σ_{33}``

## Notice
- Inplace version: `_stress_vol_tet4!(u, args...)` where u is a vector of 6 numbers.
"""
function _stress_vol_tet4(quadrature::Q,
    x1::R, x2::R, x3::R, A::U, B::U, C::U, D::U,
    e11::R, e12::R, e13::R, e22::R, e23::R, e33::R, G::R, nu::R
    ) where {R, U, Q}

    u = Vector{R}(undef, 6)
    _stress_vol_tet4!(u, quadrature, x1, x2, x3, A, B, C, D, e11, e12, e13, e22, e23, e33, G, nu)
    return u
end

function _stress_vol_tet4!(u::W, quadrature::Q,
    x1::R, x2::R, x3::R, A::U, B::U, C::U, D::U,
    e11::R, e12::R, e13::R, e22::R, e23::R, e33::R, G::R, nu::R
    ) where {R, U, Q, W}

    _strain_vol_tet4!(u, quadrature, x1, x2, x3, A, B, C, D, e11, e12, e13, e22, e23, e33, G, nu)

    lambda = 2 * nu / (1 - 2 * nu)
    ekk = u[1] + u[4] + u[6]
    u[1] = 2G * u[1] + G * lambda * ekk
    u[2] *= 2G
    u[3] *= 2G
    u[4] = 2G * u[4] + G * lambda * ekk
    u[5] *= 2G
    u[6] = 2G * u[6] + G * lambda * ekk

    return nothing
end

function _strain_vol_tet4(quadrature::Q,
    x1::R, x2::R, x3::R, A::U, B::U, C::U, D::U,
    e11::R, e12::R, e13::R, e22::R, e23::R, e33::R, G::R, nu::R
    ) where {R, U, Q}

    u = Vector{R}(undef, 6)
    _strain_vol_tet4!(u, quadrature, x1, x2, x3, A, B, C, D, e11, e12, e13, e22, e23, e33, G, nu)
    return u
end

function _strain_vol_tet4!(u::W, quadrature::Q,
    x1::R, x2::R, x3::R, A::U, B::U, C::U, D::U,
    e11::R, e12::R, e13::R, e22::R, e23::R, e33::R, G::R, nu::R
    ) where {R, U, Q, W}

    @assert x3 ≥ zero(R) "Depth must be greater than or equal to 0."
    lambda = 2 * nu / (1 - 2 * nu)
    ekk = e11 + e22 + e33

    nA = cross(C - B, D - B)
    nB = cross(D - C, A - C)
    nC = cross(A - D, B - D)
    nD = cross(B - A, C - A)

    nA /= norm(nA)
    nB /= norm(nB)
    nC /= norm(nC)
    nD /= norm(nD)

    if nA' * (A .- (B .+ C .+ D) ./ 3) > zero(R); nA *= -one(R) end
    if nB' * (B .- (A .+ C .+ D) ./ 3) > zero(R); nB *= -one(R) end
    if nC' * (C .- (B .+ A .+ D) ./ 3) > zero(R); nC *= -one(R) end
    if nD' * (D .- (B .+ C .+ A) ./ 3) > zero(R); nD *= -one(R) end

    ABC = norm(cross(C .- A, B .- A)) / 2
    BCD = norm(cross(D .- B, C .- B)) / 2
    CDA = norm(cross(A .- C, D .- C)) / 2
    DAB = norm(cross(B .- D, A .- D)) / 2

    m11 = lambda * ekk + 2 * e11
    m12 = 2 * e12
    m13 = 2 * e13
    m22 = lambda * ekk + 2 * e22
    m23 = 2 * e23
    m33 = lambda * ekk + 2 * e33

    let lambda=lambda, x1=x1, x2=x2, x3=x3, A=A, B=B, C=C, D=D, e11=e11, e12=e12, e13=e13, e22=e22, e23=e23, e33=e33, G=G, nu=nu, nA=nA, nB=nB, nC=nC, ABC=ABC, BCD=BCD, CDA=CDA, DAB=DAB, m11=m11, m12=m12, m13=m13, m22=m22, m23=m23, m33=m33

        function r1(y1::R, y2::R, y3::R) where R
            hypot(x1 - y1, x2 - y2, x3 - y3)
        end

        function r2(y1::R, y2::R, y3::R) where R
            hypot(x1 - y1, x2 - y2, x3 + y3)
        end

        function y(u::R, v::R, A::R, B::R, C::R) where R
            A * (1 - u) * (1 - v) / 4 + B * (1 + u) * (1 - v) / 4 + C * (1 + v) / 2
        end

        function G11d1(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x1 - y1) * (
            -(3 - 4 * nu) / r1(y1, y2, y3) ^ 3
            -1  / r2(y1, y2, y3) ^ 3
            +(2 * r1(y1, y2, y3) ^ 2 - 3 * (x1 - y1) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (2 * r2(y1, y2, y3) ^ 2 - 3 * (x1 - y1) ^ 2) / r2(y1, y2, y3) ^ 5
            -6 * y3 * x3 * (3 * r2(y1, y2, y3) ^ 2 - 5 * (x1 - y1) ^ 2) / r2(y1, y2, y3) ^ 7
            -4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3) ^ 2)
            -8 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3) ^ 2)
            +4 * (1 - 2 * nu) * (1 - nu) * (x1 - y1) ^ 2  / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 2)
            +8 * (1 - 2 * nu) * (1 - nu) * (x1 - y1) ^ 2  / (r2(y1, y2, y3) ^ 2  * (r2(y1, y2, y3) + x3 + y3) ^ 3) );)
        end

        function G11d2(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x2 - y2) * (
            -(3 - 4 * nu) / r1(y1, y2, y3) ^ 3
            -1  / r2(y1, y2, y3) ^ 3
            -3 * (x1 - y1) ^ 2  / r1(y1, y2, y3) ^ 5
            -3 * (3 - 4 * nu) * (x1 - y1) ^ 2  / r2(y1, y2, y3) ^ 5
            -6 * y3 * x3 * (r2(y1, y2, y3) ^ 2 - 5 * (x1 - y1) ^ 2) / r2(y1, y2, y3) ^ 7
            -4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3) ^ 2)
            +4 * (1 - 2 * nu) * (1 - nu) * (x1 - y1) ^ 2  * (3 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 3) ))
        end

        function G11d3(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (
            -(3 - 4 * nu) * (x3 - y3) / r1(y1, y2, y3) ^ 3
            -(x3 + y3) / r2(y1, y2, y3) ^ 3
            -3 * (x1 - y1) ^ 2  * (x3 - y3) / r1(y1, y2, y3) ^ 5
            -3 * (3 - 4 * nu) * (x1 - y1) ^ 2  * (x3 + y3) / r2(y1, y2, y3) ^ 5
            +2 * y3 * (r2(y1, y2, y3) ^ 2 - 3 * x3 * (x3 + y3)) / r2(y1, y2, y3) ^ 5
            -6 * y3 * (x1 - y1) ^ 2  * (r2(y1, y2, y3) ^ 2 - 5 * x3 * (x3 + y3)) / r2(y1, y2, y3) ^ 7
            -4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3))
            +4 * (1 - 2 * nu) * (1 - nu) * (x1 - y1) ^ 2  * (2 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 2) ))
        end

        function G21d1(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x2 - y2) * (
            +(r1(y1, y2, y3) ^ 2 - 3 * (x1 - y1) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (r2(y1, y2, y3) ^ 2 - 3 * (x1 - y1) ^ 2) / r2(y1, y2, y3) ^ 5
            -6 * y3 * x3 * (r2(y1, y2, y3) ^ 2 - 5 * (x1 - y1) ^ 2) / r2(y1, y2, y3) ^ 7
            -4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3) ^ 2)
            +4 * (1 - 2 * nu) * (1 - nu) * (x1 - y1) ^ 2  * (3 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 3) ))
        end

        function G21d2(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x1 - y1) * (
            +(r1(y1, y2, y3) ^ 2 - 3 * (x2 - y2) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (r2(y1, y2, y3) ^ 2 - 3 * (x2 - y2) ^ 2) / r2(y1, y2, y3) ^ 5
            -6 * y3 * x3 * (r2(y1, y2, y3) ^ 2 - 5 * (x2 - y2) ^ 2) / r2(y1, y2, y3) ^ 7
            -4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3) ^ 2)
            +4 * (1 - 2 * nu) * (1 - nu) * (x2 - y2) ^ 2  * (3 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 3) ))
        end

        function G21d3(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x1 - y1) * (x2 - y2) * (
            -3 * (x3 - y3) / r1(y1, y2, y3) ^ 5
            -3 * (3 - 4 * nu) * (x3 + y3) / r2(y1, y2, y3) ^ 5
            -6 * y3 * (r2(y1, y2, y3) ^ 2 - 5 * x3 * (x3 + y3)) / r2(y1, y2, y3) ^ 7
            +4 * (1 - 2 * nu) * (1 - nu) * (2 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 2) ))
        end

        function G31d1(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (
            +(x3 - y3) * (r1(y1, y2, y3) ^ 2 - 3 * (x1 - y1) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (x3 - y3) * (r2(y1, y2, y3) ^ 2 - 3 * (x1 - y1) ^ 2) / r2(y1, y2, y3) ^ 5
            +6 * y3 * x3 * (x3 + y3) * (r2(y1, y2, y3) ^ 2 - 5 * (x1 - y1) ^ 2) / r2(y1, y2, y3) ^ 7
            -4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3))
            +4 * (1 - 2 * nu) * (1 - nu) * (x1 - y1) ^ 2  * (2 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 2) ))
        end

        function G31d2(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x1 - y1) * (x2 - y2) * (
            -3 * (x3 - y3) / r1(y1, y2, y3) ^ 5
            -3 * (3 - 4 * nu) * (x3 - y3) / r2(y1, y2, y3) ^ 5
            -30 * y3 * x3 * (x3 + y3) / r2(y1, y2, y3) ^ 7
            +4 * (1 - 2 * nu) * (1 - nu) * (2 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 2) ))
        end

        function G31d3(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x1 - y1) * (
            +(r1(y1, y2, y3) ^ 2 - 3 * (x3 - y3) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (r2(y1, y2, y3) ^ 2 - 3 * (x3 ^ 2 - y3 ^ 2)) / r2(y1, y2, y3) ^ 5
            +6 * y3 * (2 * x3 + y3) / r2(y1, y2, y3) ^ 5
            -30 * y3 * x3 * (x3 + y3) ^ 2  / r2(y1, y2, y3) ^ 7
            +4 * (1 - 2 * nu) * (1 - nu) / r2(y1, y2, y3) ^ 3 ))
        end

        function G12d1(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x2 - y2) * (
            +(r1(y1, y2, y3) ^ 2 - 3 * (x1 - y1) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (r2(y1, y2, y3) ^ 2 - 3 * (x1 - y1) ^ 2) / r2(y1, y2, y3) ^ 5
            -6 * y3 * x3 * (r2(y1, y2, y3) ^ 2 - 5 * (x1 - y1) ^ 2) / r2(y1, y2, y3) ^ 7
            -4 * (1 - nu) * (1 - 2 * nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3) ^ 2)
            +4 * (1 - nu) * (1 - 2 * nu) * (x1 - y1) ^ 2  * (3 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 3) ))
        end

        function G12d2(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x1 - y1) * (
            +(r1(y1, y2, y3) ^ 2 - 3 * (x2 - y2) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (r2(y1, y2, y3) ^ 2 - 3 * (x2 - y2) ^ 2) / r2(y1, y2, y3) ^ 5
            -6 * y3 * x3 * (r2(y1, y2, y3) ^ 2 - 5 * (x2 - y2) ^ 2) / r2(y1, y2, y3) ^ 7
            -4 * (1 - nu) * (1 - 2 * nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3) ^ 2)
            +4 * (1 - nu) * (1 - 2 * nu) * (x2 - y2) ^ 2  * (3 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 3) ))
        end

        function G12d3(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x1 - y1) * (x2 - y2) * (
            -3 * (x3 - y3) / r1(y1, y2, y3) ^ 5
            -3 * (3 - 4 * nu) * (x3 + y3) / r2(y1, y2, y3) ^ 5
            -6 * y3 * (r2(y1, y2, y3) ^ 2 - 5 * x3 * (x3 + y3)) / r2(y1, y2, y3) ^ 7
            +4 * (1 - 2 * nu) * (1 - nu) * (2 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 2) ))
        end

        function G22d1(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x1 - y1) * (
            -(3 - 4 * nu) / r1(y1, y2, y3) ^ 3
            -1  / r2(y1, y2, y3) ^ 3
            -3 * (x2 - y2) ^ 2  / r1(y1, y2, y3) ^ 5
            -3 * (3 - 4 * nu) * (x2 - y2) ^ 2  / r2(y1, y2, y3) ^ 5
            -6 * y3 * x3 * (r2(y1, y2, y3) ^ 2 - 5 * (x2 - y2) ^ 2) / r2(y1, y2, y3) ^ 7
            -4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3) ^ 2)
            +4 * (1 - 2 * nu) * (1 - nu) * (x2 - y2) ^ 2  * (3 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 3) ))
        end

        function G22d2(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x2 - y2) * (
            -(3 - 4 * nu) / r1(y1, y2, y3) ^ 3
            -1  / r2(y1, y2, y3) ^ 3
            +(2 * r1(y1, y2, y3) ^ 2 - 3 * (x2 - y2) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (2 * r2(y1, y2, y3) ^ 2 - 3 * (x2 - y2) ^ 2) / r2(y1, y2, y3) ^ 5
            -6 * y3 * x3 * (3 * r2(y1, y2, y3) ^ 2 - 5 * (x2 - y2) ^ 2) / r2(y1, y2, y3) ^ 7
            -12 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3) ^ 2)
            +4 * (1 - 2 * nu) * (1 - nu) * (x2 - y2) ^ 2  * (3 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 3) ))
        end

        function G22d3(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (
            -(3 - 4 * nu) * (x3 - y3) / r1(y1, y2, y3) ^ 3
            -(x3 + y3) / r2(y1, y2, y3) ^ 3
            -3 * (x2 - y2) ^ 2  * (x3 - y3) / r1(y1, y2, y3) ^ 5
            -3 * (3 - 4 * nu) * (x2 - y2) ^ 2  * (x3 + y3) / r2(y1, y2, y3) ^ 5
            +2 * y3 * (r2(y1, y2, y3) ^ 2 - 3 * x3 * (x3 + y3)) / r2(y1, y2, y3) ^ 5
            -6 * y3 * (x2 - y2) ^ 2  * (r2(y1, y2, y3) ^ 2 - 5 * x3 * (x3 + y3)) / r2(y1, y2, y3) ^ 7
            -4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3))
            +4 * (1 - 2 * nu) * (1 - nu) * (x2 - y2) ^ 2  * (2 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 2) ))
        end

        function G32d1(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x1 - y1) * (x2 - y2) * (
            -3 * (x3 - y3) / r1(y1, y2, y3) ^ 5
            -3 * (3 - 4 * nu) * (x3 - y3) / r2(y1, y2, y3) ^ 5
            -30 * y3 * x3 * (x3 + y3) / r2(y1, y2, y3) ^ 7
            +4 * (1 - 2 * nu) * (1 - nu) * (2 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 2) ))
        end

        function G32d2(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (
            +(x3 - y3) * (r1(y1, y2, y3) ^ 2 - 3 * (x2 - y2) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (x3 - y3) * (r2(y1, y2, y3) ^ 2 - 3 * (x2 - y2) ^ 2) / r2(y1, y2, y3) ^ 5
            +6 * y3 * x3 * (x3 + y3) * (r2(y1, y2, y3) ^ 2 - 5 * (x2 - y2) ^ 2) / r2(y1, y2, y3) ^ 7
            -4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3))
            +4 * (1 - 2 * nu) * (1 - nu) * (x2 - y2) ^ 2  * (2 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 2) ))
        end

        function G32d3(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x2 - y2) * (
            +(r1(y1, y2, y3) ^ 2 - 3 * (x3 - y3) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (r2(y1, y2, y3) ^ 2 - 3 * (x3 ^ 2 - y3 ^ 2)) / r2(y1, y2, y3) ^ 5
            +6 * y3 * (2 * x3 + y3) / r2(y1, y2, y3) ^ 5
            -30 * y3 * x3 * (x3 + y3) ^ 2  / r2(y1, y2, y3) ^ 7
            +4 * (1 - 2 * nu) * (1 - nu) / r2(y1, y2, y3) ^ 3 ))
        end

        function G13d1(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (
            +(x3 - y3) * (r1(y1, y2, y3) ^ 2 - 3 * (x1 - y1) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (x3 - y3) * (r2(y1, y2, y3) ^ 2 - 3 * (x1 - y1) ^ 2) / r2(y1, y2, y3) ^ 5
            -6 * y3 * x3 * (x3 + y3) * (r2(y1, y2, y3) ^ 2 - 5 * (x1 - y1) ^ 2) / r2(y1, y2, y3) ^ 7
            +4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3))
            -4 * (1 - 2 * nu) * (1 - nu) * (x1 - y1) ^ 2  * (2 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 2) ))
        end

        function G13d2(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x1 - y1) * (x2 - y2) * (
            -3 * (x3 - y3) / r1(y1, y2, y3) ^ 5
            -3 * (3 - 4 * nu) * (x3 - y3) / r2(y1, y2, y3) ^ 5
            +30 * y3 * x3 * (x3 + y3) / r2(y1, y2, y3) ^ 7
            -4 * (1 - 2 * nu) * (1 - nu) * (2 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 2) ))
        end

        function G13d3(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x1 - y1) * (
            +(r1(y1, y2, y3) ^ 2 - 3 * (x3 - y3) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (r2(y1, y2, y3) ^ 2 - 3 * (x3 ^ 2 - y3 ^ 2)) / r2(y1, y2, y3) ^ 5
            -6 * y3 * (2 * x3 + y3) / r2(y1, y2, y3) ^ 5
            +30 * y3 * x3 * (x3 + y3) ^ 2  / r2(y1, y2, y3) ^ 7
            -4 * (1 - 2 * nu) * (1 - nu) / r2(y1, y2, y3) ^ 3 ))
        end

        function G23d1(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x1 - y1) * (x2 - y2) * (
            -3 * (x3 - y3) / r1(y1, y2, y3) ^ 5
            -3 * (3 - 4 * nu) * (x3 - y3) / r2(y1, y2, y3) ^ 5
            +30 * y3 * x3 * (x3 + y3) / r2(y1, y2, y3) ^ 7
            -4 * (1 - 2 * nu) * (1 - nu) * (2 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 2) ))
        end

        function G23d2(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (
            +(x3 - y3) * (r1(y1, y2, y3) ^ 2 - 3 * (x2 - y2) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (x3 - y3) * (r2(y1, y2, y3) ^ 2 - 3 * (x2 - y2) ^ 2) / r2(y1, y2, y3) ^ 5
            -6 * y3 * x3 * (x3 + y3) * (r2(y1, y2, y3) ^ 2 - 5 * (x2 - y2) ^ 2) / r2(y1, y2, y3) ^ 7
            +4 * (1 - 2 * nu) * (1 - nu) / (r2(y1, y2, y3) * (r2(y1, y2, y3) + x3 + y3))
            -4 * (1 - 2 * nu) * (1 - nu) * (x2 - y2) ^ 2  * (2 * r2(y1, y2, y3) + x3 + y3) / (r2(y1, y2, y3) ^ 3  * (r2(y1, y2, y3) + x3 + y3) ^ 2) ))
        end

        function G23d3(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x2 - y2) * (
            +(r1(y1, y2, y3) ^ 2 - 3 * (x3 - y3) ^ 2) / r1(y1, y2, y3) ^ 5
            +(3 - 4 * nu) * (r2(y1, y2, y3) ^ 2 - 3 * (x3 ^ 2 - y3 ^ 2)) / r2(y1, y2, y3) ^ 5
            -6 * y3 * (2 * x3 + y3) / r2(y1, y2, y3) ^ 5
            +30 * y3 * x3 * (x3 + y3) ^ 2  / r2(y1, y2, y3) ^ 7
            -4 * (1 - 2 * nu) * (1 - nu) / r2(y1, y2, y3) ^ 3 ))
        end

        function G33d1(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x1 - y1) * (
            -(3 - 4 * nu) / r1(y1, y2, y3) ^ 3
            -(5 - 12 * nu + 8 * nu ^ 2) / r2(y1, y2, y3) ^ 3
            -3 * (x3 - y3) ^ 2  / r1(y1, y2, y3) ^ 5
            -30 * y3 * x3 * (x3 + y3) ^ 2  / r2(y1, y2, y3) ^ 7
            -3 * (3 - 4 * nu) * (x3 + y3) ^ 2  / r2(y1, y2, y3) ^ 5
            +6 * y3 * x3 / r2(y1, y2, y3) ^ 5 ))
        end

        function G33d2(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (x2 - y2) * (
            -(3 - 4 * nu) / r1(y1, y2, y3) ^ 3
            -(5 - 12 * nu + 8 * nu ^ 2) / r2(y1, y2, y3) ^ 3
            -3 * (x3 - y3) ^ 2  / r1(y1, y2, y3) ^ 5
            -30 * y3 * x3 * (x3 + y3) ^ 2  / r2(y1, y2, y3) ^ 7
            -3 * (3 - 4 * nu) * (x3 + y3) ^ 2  / r2(y1, y2, y3) ^ 5
            +6 * y3 * x3 / r2(y1, y2, y3) ^ 5 ))
        end

        function G33d3(y1::R, y2::R, y3::R) where R
            (1 / (16 * π * (1 - nu)) * (
            -(3 - 4 * nu) * (x3 - y3) / r1(y1, y2, y3) ^ 3
            -(5 - 12 * nu + 8 * nu ^ 2) * (x3 + y3) / r2(y1, y2, y3) ^ 3
            +(x3 - y3) * (2 * r1(y1, y2, y3) ^ 2 - 3 * (x3 - y3) ^ 2) / r1(y1, y2, y3) ^ 5
            +6 * y3 * (x3 + y3) ^ 2  / r2(y1, y2, y3) ^ 5
            +6 * y3 * x3 * (x3 + y3) * (2 * r2(y1, y2, y3) ^ 2 - 5 * (x3 + y3) ^ 2) / r2(y1, y2, y3) ^ 7
            +(3 - 4 * nu) * (x3 + y3) * (2 * r2(y1, y2, y3) ^ 2 - 3 * (x3 + y3) ^ 2) / r2(y1, y2, y3) ^ 5
            -2 * y3 * (r2(y1, y2, y3) ^ 2 - 3 * x3 * (x3 + y3)) / r2(y1, y2, y3) ^ 5 ))
        end

        function IU11(u::R, v::R) where R
            (+ ABC / 4 * (m11 * nD[1] + m12 * nD[2] + m13 * nD[3]) * G11d1(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m12 * nD[1] + m22 * nD[2] + m23 * nD[3]) * G21d1(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m13 * nD[1] + m23 * nD[2] + m33 * nD[3]) * G31d1(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + BCD / 4 * (m11 * nA[1] + m12 * nA[2] + m13 * nA[3]) * G11d1(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m12 * nA[1] + m22 * nA[2] + m23 * nA[3]) * G21d1(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m13 * nA[1] + m23 * nA[2] + m33 * nA[3]) * G31d1(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + CDA / 4 * (m11 * nB[1] + m12 * nB[2] + m13 * nB[3]) * G11d1(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m12 * nB[1] + m22 * nB[2] + m23 * nB[3]) * G21d1(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m13 * nB[1] + m23 * nB[2] + m33 * nB[3]) * G31d1(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + DAB / 4 * (m11 * nC[1] + m12 * nC[2] + m13 * nC[3]) * G11d1(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m12 * nC[1] + m22 * nC[2] + m23 * nC[3]) * G21d1(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m13 * nC[1] + m23 * nC[2] + m33 * nC[3]) * G31d1(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3])))
        end

        function IU12(u::R, v::R) where R
            (+ ABC / 4 * (m11 * nD[1] + m12 * nD[2] + m13 * nD[3]) * G11d2(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m12 * nD[1] + m22 * nD[2] + m23 * nD[3]) * G21d2(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m13 * nD[1] + m23 * nD[2] + m33 * nD[3]) * G31d2(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + BCD / 4 * (m11 * nA[1] + m12 * nA[2] + m13 * nA[3]) * G11d2(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m12 * nA[1] + m22 * nA[2] + m23 * nA[3]) * G21d2(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m13 * nA[1] + m23 * nA[2] + m33 * nA[3]) * G31d2(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + CDA / 4 * (m11 * nB[1] + m12 * nB[2] + m13 * nB[3]) * G11d2(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m12 * nB[1] + m22 * nB[2] + m23 * nB[3]) * G21d2(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m13 * nB[1] + m23 * nB[2] + m33 * nB[3]) * G31d2(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + DAB / 4 * (m11 * nC[1] + m12 * nC[2] + m13 * nC[3]) * G11d2(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m12 * nC[1] + m22 * nC[2] + m23 * nC[3]) * G21d2(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m13 * nC[1] + m23 * nC[2] + m33 * nC[3]) * G31d2(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3])))
        end

        function IU13(u::R, v::R) where R
            (+ ABC / 4 * (m11 * nD[1] + m12 * nD[2] + m13 * nD[3]) * G11d3(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m12 * nD[1] + m22 * nD[2] + m23 * nD[3]) * G21d3(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m13 * nD[1] + m23 * nD[2] + m33 * nD[3]) * G31d3(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + BCD / 4 * (m11 * nA[1] + m12 * nA[2] + m13 * nA[3]) * G11d3(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m12 * nA[1] + m22 * nA[2] + m23 * nA[3]) * G21d3(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m13 * nA[1] + m23 * nA[2] + m33 * nA[3]) * G31d3(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + CDA / 4 * (m11 * nB[1] + m12 * nB[2] + m13 * nB[3]) * G11d3(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m12 * nB[1] + m22 * nB[2] + m23 * nB[3]) * G21d3(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m13 * nB[1] + m23 * nB[2] + m33 * nB[3]) * G31d3(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + DAB / 4 * (m11 * nC[1] + m12 * nC[2] + m13 * nC[3]) * G11d3(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m12 * nC[1] + m22 * nC[2] + m23 * nC[3]) * G21d3(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m13 * nC[1] + m23 * nC[2] + m33 * nC[3]) * G31d3(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3])))
        end

        function IU21(u::R, v::R) where R
            (+ ABC / 4 * (m11 * nD[1] + m12 * nD[2] + m13 * nD[3]) * G12d1(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m12 * nD[1] + m22 * nD[2] + m23 * nD[3]) * G22d1(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m13 * nD[1] + m23 * nD[2] + m33 * nD[3]) * G32d1(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + BCD / 4 * (m11 * nA[1] + m12 * nA[2] + m13 * nA[3]) * G12d1(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m12 * nA[1] + m22 * nA[2] + m23 * nA[3]) * G22d1(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m13 * nA[1] + m23 * nA[2] + m33 * nA[3]) * G32d1(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + CDA / 4 * (m11 * nB[1] + m12 * nB[2] + m13 * nB[3]) * G12d1(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m12 * nB[1] + m22 * nB[2] + m23 * nB[3]) * G22d1(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m13 * nB[1] + m23 * nB[2] + m33 * nB[3]) * G32d1(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + DAB / 4 * (m11 * nC[1] + m12 * nC[2] + m13 * nC[3]) * G12d1(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m12 * nC[1] + m22 * nC[2] + m23 * nC[3]) * G22d1(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m13 * nC[1] + m23 * nC[2] + m33 * nC[3]) * G32d1(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3])))
        end

        function IU22(u::R, v::R) where R
            (+ ABC / 4 * (m11 * nD[1] + m12 * nD[2] + m13 * nD[3]) * G12d2(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m12 * nD[1] + m22 * nD[2] + m23 * nD[3]) * G22d2(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m13 * nD[1] + m23 * nD[2] + m33 * nD[3]) * G32d2(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + BCD / 4 * (m11 * nA[1] + m12 * nA[2] + m13 * nA[3]) * G12d2(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m12 * nA[1] + m22 * nA[2] + m23 * nA[3]) * G22d2(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m13 * nA[1] + m23 * nA[2] + m33 * nA[3]) * G32d2(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + CDA / 4 * (m11 * nB[1] + m12 * nB[2] + m13 * nB[3]) * G12d2(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m12 * nB[1] + m22 * nB[2] + m23 * nB[3]) * G22d2(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m13 * nB[1] + m23 * nB[2] + m33 * nB[3]) * G32d2(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + DAB / 4 * (m11 * nC[1] + m12 * nC[2] + m13 * nC[3]) * G12d2(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m12 * nC[1] + m22 * nC[2] + m23 * nC[3]) * G22d2(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m13 * nC[1] + m23 * nC[2] + m33 * nC[3]) * G32d2(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3])))
        end

        function IU23(u::R, v::R) where R
            (+ ABC / 4 * (m11 * nD[1] + m12 * nD[2] + m13 * nD[3]) * G12d3(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m12 * nD[1] + m22 * nD[2] + m23 * nD[3]) * G22d3(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m13 * nD[1] + m23 * nD[2] + m33 * nD[3]) * G32d3(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + BCD / 4 * (m11 * nA[1] + m12 * nA[2] + m13 * nA[3]) * G12d3(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m12 * nA[1] + m22 * nA[2] + m23 * nA[3]) * G22d3(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m13 * nA[1] + m23 * nA[2] + m33 * nA[3]) * G32d3(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + CDA / 4 * (m11 * nB[1] + m12 * nB[2] + m13 * nB[3]) * G12d3(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m12 * nB[1] + m22 * nB[2] + m23 * nB[3]) * G22d3(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m13 * nB[1] + m23 * nB[2] + m33 * nB[3]) * G32d3(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + DAB / 4 * (m11 * nC[1] + m12 * nC[2] + m13 * nC[3]) * G12d3(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m12 * nC[1] + m22 * nC[2] + m23 * nC[3]) * G22d3(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m13 * nC[1] + m23 * nC[2] + m33 * nC[3]) * G32d3(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3])))
        end

        function IU31(u::R, v::R) where R
            (+ ABC / 4 * (m11 * nD[1] + m12 * nD[2] + m13 * nD[3]) * G13d1(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m12 * nD[1] + m22 * nD[2] + m23 * nD[3]) * G23d1(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m13 * nD[1] + m23 * nD[2] + m33 * nD[3]) * G33d1(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + BCD / 4 * (m11 * nA[1] + m12 * nA[2] + m13 * nA[3]) * G13d1(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m12 * nA[1] + m22 * nA[2] + m23 * nA[3]) * G23d1(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m13 * nA[1] + m23 * nA[2] + m33 * nA[3]) * G33d1(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + CDA / 4 * (m11 * nB[1] + m12 * nB[2] + m13 * nB[3]) * G13d1(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m12 * nB[1] + m22 * nB[2] + m23 * nB[3]) * G23d1(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m13 * nB[1] + m23 * nB[2] + m33 * nB[3]) * G33d1(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + DAB / 4 * (m11 * nC[1] + m12 * nC[2] + m13 * nC[3]) * G13d1(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m12 * nC[1] + m22 * nC[2] + m23 * nC[3]) * G23d1(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m13 * nC[1] + m23 * nC[2] + m33 * nC[3]) * G33d1(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3])))
        end

        function IU32(u::R, v::R) where R
            (+ ABC / 4 * (m11 * nD[1] + m12 * nD[2] + m13 * nD[3]) * G13d2(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m12 * nD[1] + m22 * nD[2] + m23 * nD[3]) * G23d2(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m13 * nD[1] + m23 * nD[2] + m33 * nD[3]) * G33d2(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + BCD / 4 * (m11 * nA[1] + m12 * nA[2] + m13 * nA[3]) * G13d2(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m12 * nA[1] + m22 * nA[2] + m23 * nA[3]) * G23d2(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m13 * nA[1] + m23 * nA[2] + m33 * nA[3]) * G33d2(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + CDA / 4 * (m11 * nB[1] + m12 * nB[2] + m13 * nB[3]) * G13d2(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m12 * nB[1] + m22 * nB[2] + m23 * nB[3]) * G23d2(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m13 * nB[1] + m23 * nB[2] + m33 * nB[3]) * G33d2(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + DAB / 4 * (m11 * nC[1] + m12 * nC[2] + m13 * nC[3]) * G13d2(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m12 * nC[1] + m22 * nC[2] + m23 * nC[3]) * G23d2(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m13 * nC[1] + m23 * nC[2] + m33 * nC[3]) * G33d2(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3])))
        end

        function IU33(u::R, v::R) where R
            (+ ABC / 4 * (m11 * nD[1] + m12 * nD[2] + m13 * nD[3]) * G13d3(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m12 * nD[1] + m22 * nD[2] + m23 * nD[3]) * G23d3(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + ABC / 4 * (m13 * nD[1] + m23 * nD[2] + m33 * nD[3]) * G33d3(y(u, v, A[1], B[1], C[1]), y(u, v, A[2], B[2], C[2]), y(u, v, A[3], B[3], C[3]))
            + BCD / 4 * (m11 * nA[1] + m12 * nA[2] + m13 * nA[3]) * G13d3(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m12 * nA[1] + m22 * nA[2] + m23 * nA[3]) * G23d3(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + BCD / 4 * (m13 * nA[1] + m23 * nA[2] + m33 * nA[3]) * G33d3(y(u, v, B[1], C[1], D[1]), y(u, v, B[2], C[2], D[2]), y(u, v, B[3], C[3], D[3]))
            + CDA / 4 * (m11 * nB[1] + m12 * nB[2] + m13 * nB[3]) * G13d3(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m12 * nB[1] + m22 * nB[2] + m23 * nB[3]) * G23d3(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + CDA / 4 * (m13 * nB[1] + m23 * nB[2] + m33 * nB[3]) * G33d3(y(u, v, C[1], D[1], A[1]), y(u, v, C[2], D[2], A[2]), y(u, v, C[3], D[3], A[3]))
            + DAB / 4 * (m11 * nC[1] + m12 * nC[2] + m13 * nC[3]) * G13d3(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m12 * nC[1] + m22 * nC[2] + m23 * nC[3]) * G23d3(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3]))
            + DAB / 4 * (m13 * nC[1] + m23 * nC[2] + m33 * nC[3]) * G33d3(y(u, v, D[1], A[1], B[1]), y(u, v, D[2], A[2], B[2]), y(u, v, D[3], A[3], B[3])))
        end

        u11, u12, u13, u21, u22, u23, u31, u32, u33 = zero(R), zero(R), zero(R), zero(R), zero(R), zero(R), zero(R), zero(R), zero(R)

        x, w = quadrature
        @inbounds @fastmath for k = 1: length(x)
            @simd for j = 1: length(x)
                c = w[j] * w[k] * (1 - x[k])
                u11 += c * IU11(x[j], x[k])
                u12 += c * IU12(x[j], x[k])
                u13 += c * IU13(x[j], x[k])
                u21 += c * IU21(x[j], x[k])
                u22 += c * IU22(x[j], x[k])
                u23 += c * IU23(x[j], x[k])
                u31 += c * IU31(x[j], x[k])
                u32 += c * IU32(x[j], x[k])
                u33 += c * IU33(x[j], x[k])
            end
        end

        omega =
            (heaviside(((A[1] + B[1] + C[1]) / 3 - x1) * nD[1] + ((A[2] + B[2] + C[2]) / 3 - x2) * nD[2] + ((A[3] + B[3] + C[3]) / 3 - x3) * nD[3])
            * heaviside(((B[1] + C[1] + D[1]) / 3 - x1) * nA[1] + ((B[2] + C[2] + D[2]) / 3 - x2) * nA[2] + ((B[3] + C[3] + D[3]) / 3 - x3) * nA[3])
            * heaviside(((C[1] + D[1] + A[1]) / 3 - x1) * nB[1] + ((C[2] + D[2] + A[2]) / 3 - x2) * nB[2] + ((C[3] + D[3] + A[3]) / 3 - x3) * nB[3])
            * heaviside(((D[1] + A[1] + B[1]) / 3 - x1) * nC[1] + ((D[2] + A[2] + B[2]) / 3 - x2) * nC[2] + ((D[3] + A[3] + B[3]) / 3 - x3) * nC[3]))

        u[1] = u11 - omega * e11
        u[2] = (u12 + u21) / 2 - omega * e12
        u[3] = (u13 + u31) / 2 - omega * e13
        u[4] = u22 - omega * e22
        u[5] = (u23 + u32) / 2 - omega * e23
        u[6] = u33 - omega * e33

        return nothing
    end
end


text:
# Struct Det 1D functor (4 methods) #

# Return V(x,ti)
function (a::SolveOutDet1D)(ti::AbstractFloat)
    t = findlast(i->i==ti,a.tsaved)
    return a.V[(t-1)*a.N+1:a.N*t]
end
(a::SolveOutDet1D)(i::Int) = a.V[(i-1)*a.N+1:a.N*i]

# Return V(xj,ti)
function (a::SolveOutDet1D)(xj::AbstractFloat,ti::AbstractFloat)
    t  = findlast(i->i==ti,a.tsaved)
    x  = findlast(j->j==xj,a.x)
    Vi = @view a.V[(t-1)*a.N+1:a.N*t]

    return Vi[x]
end
function (a::SolveOutDet1D)(j::Int,i::Int)
    Vi = @view a.V[(i-1)*a.N+1:a.N*i] # solution at ti
    return Vi[j]
end

# Struct Sto 1D functor (8 methods) #

# Return V(x,ti,p)
function (a::SolveOutSto1D)(ti::AbstractFloat,p::Int)
    t = findlast(i->i==ti,a.tsaved)
    return a.V[(t-1)*a.N+1:a.N*t,p]
end
(a::SolveOutSto1D)(i::Int,p::Int) = a.V[(i-1)*a.N+1:a.N*i,p]

# Return Vmean(x,ti)
function (b::SolveOutSto1D)(ti::AbstractFloat)
    t = findlast(i->i==ti,b.tsaved)
    return b.meanV[(t-1)*b.N+1:b.N*t]
end
(b::SolveOutSto1D)(i::Int) = b.meanV[(i-1)*b.N+1:b.N*i]

# Return V(xj,ti,p)
function (a::SolveOutSto1D)(xj::AbstractFloat,ti::AbstractFloat,p::Int)
    t  = findlast(i->i==ti,a.tsaved)
    x  = findlast(j->j==xj,a.x)
    Vi = @view a.V[(t-1)*a.N+1:a.N*t,p] # path at ti

    return Vi[x] # solution at x and ti; i.e. V(x,t)
end
function (a::SolveOutSto1D)(j::Int,i::Int,p::Int)
    Vi = @view a.V[(i-1)*a.N+1:a.N*i,p]
    return Vi[j]
end

# Return Vmean(xj,ti)
function (b::SolveOutSto1D)(xj::AbstractFloat,ti::AbstractFloat)
    t  = findlast(i->i==ti,b.tsaved)
    x  = findlast(j->j==xj,b.x)
    Vi = @view b.meanV[(t-1)*b.N+1:b.N*t]

    return Vi[x]
end


text:
module MaxHelpingHandMultiRoomStrawberry

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/MultiRoomStrawberry" MultiRoomStrawberry(x::Integer, y::Integer,
    name::String="multi_room_strawberry", winged::Bool=false, moon::Bool=false, checkpointID::Integer=-1, order::Integer=-1)

const placements = Ahorn.PlacementDict(
    "Multi-Room Strawberry (max480's Helping Hand)" => Ahorn.EntityPlacement(
        MultiRoomStrawberry
    )
)

# winged, moon
sprites = Dict{Tuple{Bool, Bool}, String}(
    (false, false) => "collectables/strawberry/normal00",
    (true, false) => "collectables/strawberry/wings01",
    (false, true) => "collectables/moonBerry/normal00",
    (true, true) => "collectables/moonBerry/normal00"
)

function Ahorn.selection(entity::MultiRoomStrawberry)
    x, y = Ahorn.position(entity)

    moon = get(entity.data, "moon", false)
    winged = get(entity.data, "winged", false)
    sprite = sprites[(winged, moon)]

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::MultiRoomStrawberry, room::Maple.Room)
    x, y = Ahorn.position(entity)

    moon = get(entity.data, "moon", false)
    winged = get(entity.data, "winged", false)
    sprite = sprites[(winged, moon)]

    Ahorn.drawSprite(ctx, sprite, x, y)
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand