Get the first and second derivative of a function

[WaremF]

Hallo,
I am trying to compute the first an the second derivative of w and gd see the code below. I am not sure about to do that. The picture show you the K and f that I would like to implement.
w(ϕ) = ϕ and
g(ϕ) = (1 - ϕ)^2.
How to compute the 1st and 2nd derivative of the function w and g i.e. w'(ϕ), w''(ϕ)and g'(ϕ), g''(ϕ). I need the derivative for computing the contribution to fe and to Ke

code:

#Ke, fe, cellvalues_ϕ, CellScalarValues, ϕe are already defined
#Element stiffness matrix
nbase_ϕ = getnbasefunctions(cellvalues_ϕ)

# Loop over integration points

for q_point in 1:getnquadpoints(cellvalues_ϕ)
    # Get the quadrature weight
    dΩᵩ =getdetJdV(cellvalues_ϕ,q_point)
    #get pf value
    ϕ = function_value(cellvalues_ϕ,q_point,ϕe)
    ∇ϕ = function_gradient(cellvalues_ϕ,q_point,ϕe)
    # define w function
    w = ϕ
    gd = (1.0-ϕ)^2
    # derivative of the function gd
    gd′ =  -2.0(1.0-ϕ)
    #Loop over test shape functions
    for i in 1:nbase_ϕ
        δϕ = shape_value(cellvalues_ϕ,q_point,i)
        δ∇ϕ = shape_gradient(cellvalues_ϕ,q_point,i)
        # Add contribution to fe
        fe[i] += (gd′δϕ + δϕ * w′ + δ∇ϕ ⋅ ∇ϕ)*dΩᵩ ##### Question: how to compute the derivative of w i.e w′ = ϕ′
        # Loop over trial shape functions
        for j in 1:i
            # I am not sure if the following is correct
            ϕ′ = shape_value(cellvalues_ϕ,q_point,j)
            ∇ϕ′= shape_gradient(cellvalues_ϕ,q_point,j)
            gd′′ = 2.0ϕ′ # 2nd derivative of gd 
            # 
            Ke[i,j] += (gd′′*δϕ+ δϕ*w′′ + δ∇ϕ⋅∇ϕ′)*dΩᵩ # Question how to get the 2nd derivative of w i.e w′′ = ϕ′′
        end
    end
end
symmetrize_lower!(Ke)
return Ke, fe

[KnutAM]

Hi 👋

In your equations, you have indices i (for f) and (i,j) (for K) on the left side, but not on the right side.

Hopefully the following notes can help, that you want to have,

$$\begin{align} f_i^\phi &= \int_\Omega N^\phi_i [g'(\phi) + w'(\phi)] + \nabla N^\phi_i \cdot \nabla \phi\ \mathrm{d}\Omega \\\ \phi &= N^\phi_j a^\phi_j\quad \left(\text{using Einstein summation: } \phi =\sum_i N^\phi_j a^\phi_j\right) \\\ K_{ij}^\phi &= \frac{\partial f_i^\phi}{\partial a_j^\phi} \end{align}$$

Hence you already have already calculated most derivatives correctly, $w'(\phi) = 1$, $w''(\phi) = 0$, $g'(\phi) = -2(1 - \phi)$, and $g''(\phi) = 2$. However, to get $K_{ij}$, you now need to differentiate wrt. the unknown dof values, $a^\phi_j$. (Denoted ϕe in your code).

Update:
The derivatives of $\phi$ and $\nabla\phi$ wrt $a_j$ are $N_j$ and $\nabla N_j$, respectively, which you get from shape_value(cv, q_point, j) and shape_gradient(cv, q_point, j).

In the current v1 of Ferrite, CellScalarValues are replaced by CellValues. If you work on an old version, make sure that you look at the documentation for that version.

[WaremF]

Hi
Thank you for the fast reply. You are right, I forgot the indices i for f and (i,j) for K on the right hand side of the equations. I have corrected the equation and the code. Would you please tell me if I implemented it correctly?
I'm currently working with an older version of ferrite, but will consider upgrading.

#Ke, fe, cellvalues_ϕ, CellScalarValues, ϕe are already defined
#Element stiffness matrix
nbase_ϕ = getnbasefunctions(cellvalues_ϕ)
# Loop over integration points
for q_point in 1:getnquadpoints(cellvalues_ϕ)
    # Get the quadrature weight
    dΩᵩ =getdetJdV(cellvalues_ϕ,q_point)
    #get pf value
    ϕ = function_value(cellvalues_ϕ,q_point,ϕe)
    ∇ϕ = function_gradient(cellvalues_ϕ,q_point,ϕe)
    # define w function
    w = ϕ
    dw = 1.0
    ddw = 0.0
    g = (1.0-ϕ)^2
    dg =  -2.0(1.0-ϕ)
    ddg = 2.0
    #Loop over test shape functions
    for i in 1:nbase_ϕ
        δϕ = shape_value(cellvalues_ϕ,q_point,i)
        δ∇ϕ = shape_gradient(cellvalues_ϕ,q_point,i)
        # Add contribution to fe
        fe[i] += (δϕ*(dg + dw) + δ∇ϕ ⋅ ∇ϕ)*dΩᵩ
        # Loop over trial shape functions
        for j in 1:i
            ϕ′ = shape_value(cellvalues_ϕ,q_point,j)
            ∇ϕ′= shape_gradient(cellvalues_ϕ,q_point,j)
            Ke[i,j] += (δϕ⋅ϕ′*(ddg + ddw) + δ∇ϕ⋅∇ϕ′)*dΩᵩ
        end
    end
end
symmetrize_lower!(Ke)
return Ke, fe