text:
# Tok-Tok Tokenizer
# Replace non-breaking spaces with normal spaces.
const NON_BREAKING = ("\u00A0",) => " "

# Pad some funky punctuation.
const FUNKY_PUNCT_1 = string.(Tuple("،;؛¿!])}»›”؟¡%٪°±©®।॥…"))
# Pad more funky punctuation.
const FUNKY_PUNCT_2 = string.(Tuple("[“‘„‚«‹「『"))
# Pad En dash and em dash
const EN_EM_DASHES = ("–", "—")

# Replace problematic character with numeric character reference.
const AMPERCENT = ("&",) => "&"
const TAB = ("\t",) => "	"
const PIPE = ("|",) => "|"

# Just pad problematic (often neurotic) hyphen/single quote, etc.
const PROB_SINGLE_QUOTES = ("'", "’", "`")
# Group ` ` stupid quotes ' ' into a single token.
const STUPID_QUOTES_1 = ("` `",) => "``"
const STUPID_QUOTES_2 = ("' '",) => "''"

# This is the \p{Open_Punctuation} from Perl's perluniprops
# see http://perldoc.perl.org/perluniprops.html
const OPEN_PUNCT = (
        "\u0f3a", "\u0f3c", "\u169b", "\u201a", "\u201e", "\u2045", "\u207d", "\uff62",
        "\u208d", "\u2329", "\u2768", "\u276a", "\u276c", "\u276e", "\u2770", "\u2772",
        "\u2774", "\u27c5", "\u27e6", "\u27e8", "\u27ea", "\u27ec", "\u27ee", "\u2983",
        "\u2985", "\u2987", "\u2989", "\u298b", "\u298d", "\u298f", "\u2991", "\u2993",
        "\u2995", "\u2997", "\u29d8", "\u29da", "\u29fc", "\u2e22", "\u2e24", "\u2e26",
        "\u2e28", "\u3008", "\u300a", "\u300c", "\u300e", "\u3010", "\u3014", "\u3016",
        "\u3018", "\u301a", "\u301d", "\ufd3e", "\ufe17", "\ufe35", "\ufe37", "\ufe39",
        "\ufe3b", "\ufe3d", "\ufe3f", "\ufe41", "\ufe43", "\ufe47", "\ufe59", "\ufe5b",
        "\ufe5d", "\uff08", "\uff3b", "\uff5b", "\uff5f")

# This is the \p{Close_Punctuation} from Perl's perluniprops
const CLOSE_PUNCT = (
        "\u0f3b", "\u0f3d", "\u169c", "\u2046", "\u207e", "\u208e", "\u232a", "\uff63",
        "\u2769", "\u276b", "\u276d", "\u276f", "\u2771", "\u2773", "\u2775", "\u27c6",
        "\u27e7", "\u27e9", "\u27eb", "\u27ed", "\u27ef", "\u2984", "\u2986", "\u2988",
        "\u298a", "\u298c", "\u298e", "\u2990", "\u2992", "\u2994", "\u2996", "\u2998",
        "\u29d9", "\u29db", "\u29fd", "\u2e23", "\u2e25", "\u2e27", "\u2e29", "\u3009",
        "\u300b", "\u300d", "\u300f", "\u3011", "\u3015", "\u3017", "\u3019", "\u301b",
        "\u301e", "\u301f", "\ufd3f", "\ufe18", "\ufe36", "\ufe38", "\ufe3a", "\ufe3c",
        "\ufe3e", "\ufe40", "\ufe42", "\ufe44", "\ufe48", "\ufe5a", "\ufe5c", "\ufe5e",
        "\uff09", "\uff3d", "\uff5d", "\uff60")

# This is the \p{Close_Punctuation} from Perl's perluniprops
const CURRENCY_SYM = (
        "\u00a2", "\u00a3", "\u00a4", "\u00a5", "\u058f", "\u060b", "\u09f2", "\u09f3",
        "\u0af1", "\u0bf9", "\u0e3f", "\u17db", "\u20a0", "\u20a1", "\u20a2", "\u20a3",
        "\u20a4", "\u20a5", "\u20a6", "\u20a7", "\u20a8", "\u20a9", "\u20aa", "\u20ab",
        "\u20ac", "\u20ad", "\u20ae", "\u20af", "\u20b0", "\u20b1", "\u20b2", "\u20b3",
        "\u20b4", "\u20b5", "\u20b6", "\u20b7", "\u20b8", "\u20b9", "\u20ba", "\ua838",
        "\ufdfc", "\ufe69", "\uff04", "\uffe0", "\uffe1", "\uffe5", "\uffe6", "\u09fb")

const rules_atoms = collect.(Tuple(Iterators.flatten([
        CURRENCY_SYM,
        FUNKY_PUNCT_1,
        FUNKY_PUNCT_2,
        EN_EM_DASHES,
        PROB_SINGLE_QUOTES,
        OPEN_PUNCT,
        CLOSE_PUNCT])))

const rules_replaces = Tuple(Iterators.flatten([
        [AMPERCENT],
        [TAB],
        [PIPE],
        [STUPID_QUOTES_1],
        [STUPID_QUOTES_2]]))

"""
    toktok_tokenize(instring::AstractString)

This tokenizer is a simple, general tokenizer, where the input has one sentence per line; thus only final period is tokenized.
This is an enhanced version of the [original toktok Tokenizer](https://github.com/jonsafari/tok-tok)
It has been tested on and gives reasonably good results for English, Persian, Russian, Czech, French, German, Vietnamese,
Tajik, and a few others.
"""
function toktok_tokenize(instring::AbstractString)
    ts = TokenBuffer(instring)
    isempty(ts.input) && return ts.tokens

    effective_end, flush_later1, flush_later2 = handle_final_periods(ts)

    while !isdone(ts) && ts.idx <= effective_end
        if string(ts.input[ts.idx]) == NON_BREAKING[1]
            ts.input[ts.idx] = NON_BREAKING[2]
        end

        url_handler4(ts) ||   # these url handlers have priority over others
        url_handler3(ts) ||
        url_handler2(ts) ||
        url_handler1(ts) ||
        repeated_character_seq(ts, ',', 2) ||
        repeated_character_seq(ts, '-', 2) ||
        repeated_character_seq(ts, '.', 2) ||
        number(ts) ||
        spaces(ts) ||      # Handles ONE_SPACE rules from original toktok perl script
        replaces(ts, rules_replaces) ||    # most expensive steps, keep low priority
        atoms(ts, rules_atoms) ||
        character(ts)
    end
    if flush_later1 == nothing
        flush!(ts)
    elseif flush_later2 == nothing
        flush!(ts, flush_later1)
    else
        flush!(ts, flush_later1, flush_later2)
    end

    return ts.tokens
end

"""
    handle_final_periods(::TokenBuffer)

Handles the following rules from original toktok perl script:
Don't tokenize period unless it ends the line and that it isn't preceded by another period (FINAL_PERIOD_1)
Don't tokenize period unless it ends the line (FINAL_PERIOD_2)
"""
function handle_final_periods(ts::TokenBuffer)
    effective_end = length(ts.input)

    # handles spaces
    while effective_end >=1 && isspace(ts.input[effective_end])
        effective_end -= 1
    end

    # handles FINAL_PERIOD_1 = r"(?<!\.)\.$"
    if effective_end > 1 && length(ts.input) >= 2 && ts.input[effective_end] == '.' && ts.input[effective_end-1] != '.'
        effective_end -= 1
        return effective_end, ".", nothing
    end

    # handles FINAL_PERIOD_2 = r"(?<!\.)\.\s*(["'’»›”]) *$"
    if effective_end > 1 && ts.input[effective_end] in ('\"', '“', '”', '‘', '’', '›')
        token_position = effective_end
        effective_end -= 1

        while effective_end >=1 && isspace(ts.input[effective_end])
            effective_end -= 1
        end

        if effective_end > 1 && ts.input[effective_end] == '.'
            if effective_end >= 2 && ts.input[effective_end - 1] == '.'
                return token_position + 1, nothing, nothing # No use iterating over spaces again.
            else
                effective_end -= 1
                return effective_end, ".", string(ts.input[token_position])
            end
        end
        return effective_end, string(ts.input[token_position]), nothing
    end

    return effective_end, nothing, nothing
end

# In below functions flush!() is used when some given string needs to be a seperate token
# and push!() is used when it needs to appended to some processsing token in buffer i.e not as seperate token
"""
    url_handler1(::TokenBuffer)

Handles this rule from python using TokenBuffer API:
URL_FOE_1 = re.compile(r':(?!//)'), r' : '
"""
function url_handler1(ts::TokenBuffer)
    str = ':'
    pattern = ['/', '/']
    if ts.idx + length(pattern) <= length(ts.input) && ts.input[ts.idx] == str
        if ts.input[ts.idx + 1 : ts.idx + length(pattern)] != pattern
            flush!(ts, string(str))
            ts.idx += 1
            return true
        else
            for i in 1 : length(pattern) + 1
                push!(ts.buffer, ts.input[ts.idx])
                ts.idx += 1
            end
            return true
        end
    end
    return false
end

"""
    url_handler2(::TokenBuffer)

Handles this rule from python using TokenBuffer API:
URL_FOE_2 = re.compile(r'\\?(?!\\S)'), r' ? '
"""
function url_handler2(ts::TokenBuffer)
    str = '?'
    if ts.idx + 1 <= length(ts.input) && ts.input[ts.idx] == str
        if isspace(ts.input[ts.idx + 1])
            flush!(ts, string(str))
            ts.idx += 1
            return true
        else
            for i in 1 : 2
                push!(ts.buffer, ts.input[ts.idx])
                ts.idx += 1
            end
            return true
        end
    end
    return false
end

"""
    url_handler3(::TokenBuffer)

Handles this rule from python using TokenBuffer API:
URL_FOE_3 = re.compile(r'(://)[\\S+\\.\\S+/\\S+][/]'), ' / '
"""
function url_handler3(ts::TokenBuffer)
    str = "://"
    pattern = ['+', '/', '.']
    if lookahead(ts, str) && ts.idx + 4 <= length(ts.input)
        if (ts.input[ts.idx + length(str)] in pattern || isspace(ts.input[ts.idx + length(str)])) &&
            ts.input[ts.idx + length(str) + 1] == '/'
            flush!(ts, "/")
            ts.idx += 5
            return true
        end
    end
    return false
end

"""
    url_handler4(::TokenBuffer)

Handles this rule from python using TokenBuffer API:
URL_FOE_4 = re.compile(r' /'), r' / '
"""
function url_handler4(ts::TokenBuffer)
    if lookahead(ts, " /")
        flush!(ts, "/")
        ts.idx += 2
        return true
    elseif lookahead(ts, "/")
        push!(ts.buffer, '/')
        ts.idx += 1
        return true
    end
    return false
end

"""
    repeated_character_seq(::TokenBuffer, char, min_repeats=2)

Matches sequences of characters that are repreated at least `min_repeats` times.
Treat them as fake characters and ignores them.
"""
function repeated_character_seq(ts, char, min_repeats = 2)
    if ts.input[ts.idx] != char
        return false
    end
    ind = ts.idx
    while ind <= length(ts.input) && ts.input[ind] == char
        ind += 1
    end

    if ind - ts.idx >= min_repeats
        ts.idx = ind
        return true
    end

    return false
end


text:
using UnitlessFlatten, Unitful, Test
import UnitlessFlatten: flattenable

@flattenable struct Partial{T}
    a::T | true
    b::T | true
    c::T | false
end

@flattenable struct NestedPartial{P,T}
    np::P | true
    nb::T | true
    nc::T | false
end

# With units
partialunits = Partial(1.0u"s", 2.0u"s", 3.0u"s")
nestedunits = NestedPartial(Partial(1.0u"km", 2.0u"km", 3.0u"km"), 4.0u"g", 5.0u"g")
@test flatten(Vector, partialunits) == [1.0, 2.0]
@test flatten(Vector, reconstruct(partialunits, flatten(Vector, partialunits))) == flatten(Vector, partialunits)
@test flatten(Tuple, reconstruct(partialunits, flatten(Tuple, partialunits))) == flatten(Tuple, partialunits)
@test flatten(Vector, reconstruct(nestedunits, flatten(Vector, nestedunits))) == flatten(Vector, nestedunits)
@test flatten(Tuple, reconstruct(nestedunits, flatten(Tuple, nestedunits))) == flatten(Tuple, nestedunits)
@inferred flatten(Tuple, reconstruct(nestedunits, flatten(Tuple, nestedunits))) == flatten(Tuple, nestedunits)

ext: jl


text:
# abstract type and required methods for a shock sensor

"""
  Abstract type for all shock sensors.  The purpose of shock sensors it to tell
  whether or not there is a shock, and if so, what the viscoscity coefficient
  should be.  The type of artificial viscoscity to add is determined by the
  [`AbstractShockCapturing`](@ref).  Each type must implement:

  ```
  getShockSensor
  getShockSensor_diff
  setAlpha
  ```

  Each type may optionally implement

  ```
  isShockElement
  updateMetrics
  ```
"""
abstract type _AbstractShockSensor end  # this is a dummy type, the real one
                                        # is declared in defs.jl, this one
                                        # is here for documentation purposes
                                        # only.

"""
  Returns true if the element has non-zero diffusion in it, false otherwise.

  This function has a fallback defined in terms of `getShockSensor`, but types
  may be able to implement this more efficiently.

  **Inputs**

   * params: a `AbstractParamType`
   * sbp: SBP operator
   * sensor: the `AbstractShockSensor` object
   * q: solution on a particular element, `numDofPerNode` x `numNodesPerElement`
   * elnum: element number
   * coords: `dim` x `numNodesPerElement` array containing the coordinates
             of each node of the element
   * dxidx: `dim` x `dim` x `numNodesPerElement` array containing the (scaled)
            mapping jacobian dxi/dx
   * jac: mapping jacobian determinant for each node of the element, length
          `numNodesPerElement`

  **Outputs**

   * is_shock: bool
"""
function isShockElement(params::AbstractParamType{Tdim}, sbp::AbstractOperator,
                        sensor::AbstractShockSensor,
                        q::AbstractMatrix{Tsol}, elnum::Integer,
                        coords::AbstractMatrix,
                        dxidx::Abstract3DArray, jac::AbstractVector{Tmsh},
                        ) where {Tsol, Tmsh, Tdim}

  Tres = promote_type(Tsol, Tmsh)
  numDofPerNode, numNodesPerElement = size(q)
  dim = size(dxidx, 1)
  Se = EmptyArray{Tres}(dim, numNodesPerElement)
  ee = EmptyArray{Tres}(dim, numNodesPerElement)

  return getShockSensor(params, sbp, sensor, q, elnum, coords, dxidx, jac, Se,
                        ee)
end

"""
  Computes the shock sensor value and the associated viscoscity.  The
  viscoscity can be different in each Cartesian and at each node of the
  element,  A full diffusion tensor is not supported.

  **Inputs**

   * params: a `AbstractParamType`
   * sbp: SBP operator
   * sensor: the `AbstractShockSensor` object
   * q: solution on a particular element, `numDofPerNode` x `numNodesPerElement`
   * elnum: element number
   * coords: `dim` x `numNodesPerElement` array containing the coordinates
             of each node of the element
   * dxidx: `dim` x `dim` x `numNodesPerElement` array containing the (scaled)
            mapping jacobian dxi/dx
   * jac: mapping jacobian determinant for each node of the element, length
          `numNodesPerElement`

  **Inputs/Outputs**

   * Se: `dim` x `numNodesPerElement` containing the shock sensor in each
         direction at each node.
   * ee: the viscoscity coefficient, same size as `Se`

  **Outputs**

   * is_nonzero: Bool, true if ee is zero at all nodes, false otherwise

"""
function getShockSensor(params::AbstractParamType{Tdim}, sbp::AbstractOperator,
                          sensor::AbstractShockSensor,
                          q::AbstractMatrix{Tsol}, elnum::Integer,
                          coords::AbstractMatrix,
                          dxidx::Abstract3DArray, jac::AbstractVector{Tmsh},
                          Se::AbstractMatrix, ee::AbstractMatrix
                         ) where {Tsol, Tmsh, Tdim}

  error("abstract method for getShockSensor() called: did you forget to extend it with a new method for shock sensor $(typeof(sensor))")

end

"""
  Differentiated version of [`getShockSensor`](@ref).  Most of the inputs
  and outputs are the same as that function, see its documentation.

  **Inputs**

   * params
   * sbp
   * q
   * elnum
   * coords
   * dxidx
   * jac

  **Inputs/Outputs**

   * Se_jac: `dim` x `numDofPerNode` x `numNodesPerElement` x
             `numNodesPerElement` array to be overwritten with jacobian of
             `Se` wrt `q`
   * ee_jac: array to be overwritten with jacobian of `ee` wrt `q`, same size as
             `Se_jac`

  **Outputs**

   * is_const: true if the derivative of ee wrt q is zero.  This is useful
               for skipping terms in the jacobian calculation

"""
function getShockSensor_diff(params::AbstractParamType{Tdim}, sbp::AbstractOperator,
                      sensor::AbstractShockSensor,
                      q::AbstractMatrix{Tsol}, elnum::Integer,
                      coords::AbstractMatrix, dxidx::Abstract3DArray,
                      jac::AbstractVector{Tmsh},
                      Se_jac::Abstract4DArray{Tres},
                      ee_jac::Abstract4DArray{Tres}
                     ) where {Tsol, Tmsh, Tres, Tdim}

  error("abstract method for getShockSensor_diff() called: did you forget to extend it with a new method for shock sensor $(typeof(sensor))")

end

"""
  This function sets an arbitrary constant in front of the shock viscosity.

  **Inputs**

   * obj: an `AbstractShockSensor`
   * alpha: the constant
"""
function setAlpha(obj::AbstractShockSensor, alpha::Number)

  error("abstract method for setAlpha() called: did you forget to extend it with a new method for shock sensor $(typeof(obj))")
end

"""
  Function to be called after the mesh metrics are updated.  This gives the
 eshock sensor the opportunity to recalculate mesh-wide quantities.

  **Inputs**

   * mesh
   * sbp
   * opts
   * sensor: an [`AbstractShockSensor`](@ref)
"""
function updateMetrics(mesh, sbp, opts, sensor::AbstractShockSensor)

end

"""
  Reverse mode with respect to `q` of the shock sensor

  **Inputs**

   * params
   * sbp
   * sensor: `AbstractShockSensor`
   * q
   * elnum
   * coords
   * dxidx
   * jac
   * ee_bar: adjoint part of `ee`, same size

  **Inputs/Outputs**

   * ee
   * q_bar: to be updated with result, same size as `q`

  **Outputs**

   * is_nonzero: Bool, true if ee is zero at all nodes, false otherwise
"""
function getShockSensor_revq(params::AbstractParamType, sbp::AbstractOperator,
                        sensor::AbstractShockSensor,
                        q::AbstractMatrix, q_bar::AbstractMatrix,
                        elnum::Integer,
                        coords::AbstractMatrix,
                        dxidx::Abstract3DArray, jac::AbstractVector,
                        ee_mat::AbstractMatrix, ee_bar::AbstractMatrix
                        )

  error("abstract method for getShockSensor_revq() called: did you forget to extend it with a new method for shock sensor $(typeof(sensor))")

end

"""
  Reverse mode wrt metrics of the shock sensor

  **Inputs**

   * params
   * sbp
   * sensor
   * q
   * elnum
   * coords
   * dxidx
   * jac
   * ee_bar: seed matrix for back-propigation

  **Inputs/Outputs**

   * coords_bar
   * dxidx_bar
   * jac_bar
"""
function getShockSensor_revm(params::AbstractParamType, sbp::AbstractOperator,
                        sensor::AbstractShockSensor,
                        q::AbstractMatrix,
                        elnum::Integer,
                        coords::AbstractMatrix, coords_bar::AbstractMatrix,
                        dxidx::Abstract3DArray, dxidx_bar::Abstract3DArray,
                        jac::AbstractVector, jac_bar::AbstractVector,
                        ee_bar::AbstractMatrix
                       )

  error("abstract method for getShockSensor_revm() called: did you forget to extend it with a new method for shock sensor $(typeof(sensor))")

end

"""
  Performs any initialization required by the shock sensor at the beginning of
  a reverse mode wrt metrics calculation.  This is mainly used for zeroing out
  arrays inside the shock sensor that will be += during the computation.

  This function has an (empty) default function, so shock sensors only need to
  implement this function if they have to do setup work.

  **Inputs**

   * sensor: an `AbstractShockSensor`
"""
function initForRevm(sensor::AbstractShockSensor)

  return nothing
end

"""
  Finish back propigation with respect to the metrics by back propgating any
  array in the `sensor` to the `mesh`.

  This function has an (empty) default function, so shock sensors only need to
  implement this function if they have work to do.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * sensor: `AbstractShockSensor`
"""
function finishRevm(mesh::AbstractMesh, sbp::AbstractOperator, eqn::EulerData,
                    opts, sensor::AbstractShockSensor)

  return nothing
end

#------------------------------------------------------------------------------
# Required interface

function getShockSensor(eqn::EulerData)

  return getShockSensor(eqn.shock_capturing)
end

function setShockSensor(eqn::EulerData, sensor::AbstractShockSensor)

  setShockSensor(eqn.shock_capturing, sensor)
end

function setShockSensorAlpha(eqn::EulerData, alpha::Number)

  sensor = getShockSensor(eqn)
  setAlpha(sensor, alpha)
end


text:
#
# Combinator-based embedded query language.
#

module DataKnots

export
    @query,
    Collect,
    Count,
    DataKnot,
    Drop,
    Each,
    Exists,
    Filter,
    First,
    Get,
    Given,
    Group,
    Is,
    Is0to1,
    Is0toN,
    Is1to1,
    Is1toN,
    It,
    Join,
    Keep,
    Label,
    Last,
    Let,
    Lift,
    Max,
    Min,
    Mix,
    Nth,
    Record,
    Sort,
    Sum,
    Tag,
    Take,
    Unique,
    unitknot

include("layouts.jl")
include("vectors.jl")
include("shapes.jl")
include("knots.jl")
include("pipelines.jl")
include("rewrites.jl")
include("queries.jl")

end


text:
# ---
# title: 450. Delete Node in a BST
# id: problem450
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/delete-node-in-a-bst/description/>
# hidden: true
# ---
#
# Given a root node reference of a BST and a key, delete the node with the given
# key in the BST. Return the root node reference (possibly updated) of the BST.
#
# Basically, the deletion can be divided into two stages:
#
#   1. Search for a node to remove.
#   2. If the node is found, delete the node.
#
# **Follow up:**  Can you solve it with time complexity `O(height of tree)`?
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/04/del_node_1.jpg)
#
#
#
#     Input: root = [5,3,6,2,4,null,7], key = 3
#     Output: [5,4,6,2,null,null,7]
#     Explanation: Given key to delete is 3. So we find the node with value 3 and delete it.
#     One valid answer is [5,4,6,2,null,null,7], shown in the above BST.
#     Please notice that another valid answer is [5,2,6,null,4,null,7] and it's also accepted.
#     ![](https://assets.leetcode.com/uploads/2020/09/04/del_node_supp.jpg)
#
#
# **Example 2:**
#
#
#
#     Input: root = [5,3,6,2,4,null,7], key = 0
#     Output: [5,3,6,2,4,null,7]
#     Explanation: The tree does not contain a node with value = 0.
#
#
# **Example 3:**
#
#
#
#     Input: root = [], key = 0
#     Output: []
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[0, 104]`.
#   * `-105 <= Node.val <= 105`
#   * Each node has a **unique** value.
#   * `root` is a valid binary search tree.
#   * `-105 <= key <= 105`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export consistent_tangent_matrix, isotropic_function, strain_voigt_to_tensor, bi_isotropic_function
@doc raw"""
    consistent_tangent_matrix(inputs::Union{Array{Float64, 2}, PyObject},Dc::Union{Array{Float64,2}, PyObject})

Returns the consistent tangent matrices. The size of the return is $N\times 3 \times 3$.

$$D_c^{ep} = D_c - \frac{D_c\frac{\partial g}{\partial \sigma}\left(\frac{\partial f}{\partial \sigma} \right)^T D_c }{E_p + \frac{\partial g}{\partial \sigma} Dc \left(\frac{\partial f}{\partial \sigma} \right)^T}$$

Here `inputs` is a $N\times 7$ matrix, where each row is

$$\left[\frac{\partial g}{\partial \sigma}, \frac{\partial f}{\partial \sigma}, E_p\right]$$

`Dc` is a $3\times 3$ **row-major** matrix; each row is a linear elasticity matrix.

!!! info
    The gradients with respect to both `inputs` and `H` are available.
"""
function consistent_tangent_matrix(inputs::Union{Array{Float64, 2}, PyObject},H::Union{Array{Float64,2}, PyObject})
    @assert size(inputs,2)==7
    @assert size(H,1)==size(H,2)==3
    plasticity_ = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib","plasticity")
    inputs,H = convert_to_tensor([inputs,H], [Float64,Float64])
    out = plasticity_(inputs,H)
    set_shape(out, (size(inputs, 1), 3, 3))
end

@doc raw"""
    isotropic_function(coef::Union{Array{Float64,2}, PyObject},strain::Union{Array{Float64,2}, PyObject})

Isotropic function of a symmetric tensor.

$$T = s_0 I + s_1 A + s_2 A^2$$

Here

$$\texttt{coef(i,:)} = [s_0\ s_1\ s_2]$$

$$A = \left[\begin{matrix}
S_{i1} & S_{i3}/2\\
S_{i3}/2 & S_{i2}
\end{matrix}\right]$$

where $S_{ij}$ is the $i$-th row and $j$-th column of `strain`
"""
function isotropic_function(coef::Union{Array{Float64,2}, PyObject},strain::Union{Array{Float64,2}, PyObject})
    @assert size(strain,2)==size(coef,2)==3
    @assert size(coef,1)==size(strain,1)
    isotropic_ = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib","isotropic")
    coef,strain = convert_to_tensor([coef,strain], [Float64,Float64])
    out = isotropic_(coef,strain)
    set_shape(out, (size(coef,1), 3))
end

@doc raw"""
    strain_voigt_to_tensor(inp::Union{Array{Float64,2}, PyObject})

Converts Voigt strain tensors to matrix form

$$\begin{bmatrix}
\epsilon_{11}\\
\epsilon_{22}\\
2\epsilon_{12}\\
\end{bmatrix} \rightarrow \begin{bmatrix}
\epsilon_{11} &  \epsilon_{12}\\
\epsilon_{12} & \epsilon_{22}
\end{bmatrix}$$

The input is $N\times 3$ and the output is $N\times 2 \times 2$
"""
function strain_voigt_to_tensor(inp::Union{Array{Float64,2}, PyObject})
    @assert size(inp,2)==3
    tensor_rep_ = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib","tensor_rep")
    inp = convert_to_tensor([inp], [Float64]); inp = inp[1]
    set_shape(tensor_rep_(inp), (size(inp,1), 2,2))
end

@doc raw"""
    bi_isotropic_function(coef::Union{Array{Float64,2}, PyObject},strain::Union{Array{Float64,2}, PyObject},
    strainrate::Union{Array{Float64,2}, PyObject})

Isotropic function of two symmetric tensors $A$ and $B$

$$\begin{aligned}
T(A, B) &= \gamma_0 I + \gamma_1 A + \gamma_2 B + \gamma_3 A^2 + \gamma_4 (A*B+B*A)\\
& + \gamma_5 B^2 + \gamma_6 (A^2 B + BA^2) + \gamma_7 (AB^2 + B^2 A) + \gamma_8 (A^2B^2+B^2A^2)
\end{aligned}$$

This is useful for rate-dependent constitutive relation. For example,

$$\sigma^{n+1} = \mathcal{C}(\epsilon^{n+1}, \dot\epsilon^{n+1})$$

Here

- `coef`: $N\times 9$ coefficient array

- `strain`: $N\times 3$ strain tensor in Voigt notation

- `strainrate`: $N\times 3$ strain rate tensor in Voigt notation
"""
function bi_isotropic_function(coef::Union{Array{Float64,2}, PyObject},strain::Union{Array{Float64,2}, PyObject},
    strainrate::Union{Array{Float64,2}, PyObject})
    @assert size(coef, 2)==9
    isotropic_two_ = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib","isotropic_two")
    coef,strain,strainrate = convert_to_tensor([coef,strain,strainrate], [Float64,Float64,Float64])
    T = isotropic_two_(coef,strain,strainrate)
    set_shape(T, (size(strain, 1), 3))
end