Defining a new PDE

ATOmiCS wraps FEniCS for the solution of the PDEs. The input for FEniCS is the variational (weak) form. Thus, to implement a physics problem that is not included in ATOmiCS, the users need to implement the variational form using FEniCS unified form language (UFL) and add it to the ATOmiCS PDE folder (atomics/atomics/pdes). We demonstrate this process using an example that adds the st. Venant–Kirchhoff model to the ATOmiCS the PDE folder.

1.1 Derive the variational form

St. Venant–Kirchhoff model is a type of simple hyperelastic problem. The variational form of this kind of hyperelastic problem can be derived by minimizing the potential energy on domain \(\Omega\)

\[\Psi = \int_{\Omega} \psi(u) \, {\rm d} x - \int_{\Omega} b \cdot u \, {\rm d} x - \int_{\partial\Omega} t \cdot u \, {\rm d} s,\]

where \(\psi\) is the strain energy density for st. Venant–Kirchhoff model; \(b\) and \(t\) are the body and the traction force, respectively; \(u\) is the displacements function.

We minimize this potential energy by making its directional derivative with respect to the change in \(u\).

\[D_{v} \Psi = \left. \frac{d \Psi(u + \epsilon v)}{d\epsilon} \right|_{\epsilon = 0}=0.\]

We compute the \(\psi(u)\) according to the st. Venant–Kirchhoff model. We consider the deformation gradient

\[F = I + \nabla u;\]

the right Cauchy–Green tensor

\[C = F^{T} F;\]

the Green–Lagrange strain tensor

\[E = \frac{1}{2} (C - I)\]

the stress tensor

\[S = 2 \mu E + \lambda {\rm trace}(E)I;\]

and finally, the strain energy density

\[\psi(u) = \frac{1}{2} S:E\]

this final weak form is also equivalent to

\[D_{v} \int_{\Omega} \psi(u) \, {\rm d} x - \int_{\Omega} b \cdot v \ - \int_{\partial\Omega} t \cdot v \, {\rm d} s =0\]

1.2 Write the variational form in atomics/atomics/pdes/<filename>.py using UFL

We write the bilinear term as below:

import dolfin as df
import numpy as np

def get_residual_form(u, v, rho_e, method='RAMP'):
    df.dx = df.dx(metadata={"quadrature_degree": 4})
    if method == 'SIMP':
        stiffness = rho_e**3
    else:  #RAMP
        stiffness = rho_e / (1 + 8. * (1. - rho_e))

    # Kinematics
    k = 3e1
    E = k * stiffness
    nu = 0.3
    mu, lmbda = (E / (2 * (1 + nu))), (E * nu / ((1 + nu) * (1 - 2 * nu)))

    d = len(u)
    I = df.Identity(d)  # Identity tensor
    F = I + df.grad(u)  # Deformation gradient
    C = F.T * F  # Right Cauchy-Green tensor

    E_ = 0.5 * (C - I)  # Green--Lagrange strain
    S = 2.0 * mu * E_ + lmbda * df.tr(E_) * df.Identity(
        d)  # stress tensor (C:eps)
    psi = 0.5 * df.inner(S, E_)  # 0.5*eps:C:eps

    # Total potential energy
    Pi = psi * df.dx
    # Solve weak problem obtained by differentiating Pi:
    res = df.derivative(Pi, u, v)
    return res

Since the number of integration terms in the linear term

\[\int_{\Omega} b \cdot u \, {\rm d} x + \int_{\partial\Omega} t \cdot u \, {\rm d} s,\]

and the subdomain where they are applied on subject to are problem specific. We define this term in the run_file. The user can refer to Steps for solving a topology optimization problem in ATOMiCS on how to define a subdomain.

residual_form -= df.dot(f, v) * dss(6)