Velocity class

The Velocity class sets up and solves the following weak formulation: Find \(\theta \in \mathbb{H}\) such that $$ B(\theta, \xi) = - \int_{\mathcal{D}} S_0\cdot\xi + S_1 : D\xi \quad\forall \; \xi \in \mathbb{H}, $$ where

• \(\mathbb{H} = H^1(\mathcal{D})^d\) or \(\mathbb{H} = H_0^1(\mathcal{D})^d\);
• \(B\) is a bilinear form defined in the model class;
• \(S_0\) and \(S_1\) are shape derivative components of either the cost functional \(J\) or the Lagrangian \(L\) (if there are constraints).

The velocity field \(\theta\) is used to update the level set function.

---

Details

This class build and solve the velocity equation.

Attributes:

Name	Type	Description
biform	Form	Bilinear form.
liform	Form	Linear functional.
bc	List[DirichletBC]	List with the homogeneous Dirichlet condition.
solver	KSP	Krylov Subspace Solver

Methods:

Name	Description
run	Solves the velocity equation.

Source code in code/formopt.py

class Velocity:
    """
    This class build and solve the velocity equation.

    Attributes
    ----------
    biform : Form
        Bilinear form.
    liform : Form
        Linear functional.
    bc : List[DirichletBC]
        List with the homogeneous Dirichlet condition.
    solver : PETSc.KSP
        Krylov Subspace Solver

    Methods
    -------
    run(theta: Function) -> None
        Solves the velocity equation.
    """

    def __init__(
        self,
        dim: int,
        domain: Mesh,
        space: FunctionSpace,
        biform: Callable[[Argument, Argument], Tuple[Expr, bool]],
        S0: Expr,
        S1: Expr,
    ) -> None:
        """
        Sets up the linear system for the velocity.
        Since the bilinear part remains unchanged,
        it is compiled here.

        Parameters
        ----------
        dim : int
            Domain dimension.
        domain : Mesh
            Problem domain.
        space : FunctionSpace
            Space of functions.
        biform : Callable[[Argument, Argument], Tuple[Expr, bool]]
            Bilinear form method of a subclass of Model.
        S0 : Expr
            S0 component of the shape derivative.
        S1 : Expr
            S1 component of the shape derivative.
        """

        th = TrialFunction(space)
        xi = TestFunction(space)
        dx = Measure("dx", domain=domain)

        b, dirbc = biform(th, xi)

        self.biform = form(b)
        self.liform = form(-(dot(S0, xi) + inner(S1, grad(xi))) * dx)
        self.bc = None

        if dirbc == True:
            # Homogeneous Dirichlet boundary condition on the whole boundary
            # Recall that domain dimension = velocity rank
            self.bc = homogeneus_boundary(domain, space, dim, dim)
        elif isinstance(dirbc, tuple):
            # Homogeneous Dirichlet boundary condition on a subset of the boundary.
            # dirbc is a tuple with two components,
            # boundary tags and a list of marks. For instance,
            # dirbc = (tags, [mark0, mark1])
            self.bc = homogeneous_dirichlet(domain, space, dirbc[0], dirbc[1], dim)
        else:
            # By defult,
            # homogeneous Neumann boundary condition on the whole boundary
            self.bc = []

        self.solver = build_solver(domain, self.biform, self.bc)

    def run(self, theta: Function) -> None:
        """
        Solves the velocity equation.
        Only the linear part must be updated.
        """

        L = assemble_vector(self.liform)
        apply_lifting(L, [self.biform], [self.bc])
        L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        set_bc(L, self.bc)

        self.solver.solve(L, theta.x.petsc_vec)
        theta.x.scatter_forward()

__init__(dim, domain, space, biform, S0, S1)

Sets up the linear system for the velocity. Since the bilinear part remains unchanged, it is compiled here.

Parameters:

Name	Type	Description	Default
dim	int	Domain dimension.	required
domain	Mesh	Problem domain.	required
space	FunctionSpace	Space of functions.	required
biform	Callable[[Argument, Argument], Tuple[Expr, bool]]	Bilinear form method of a subclass of Model.	required
S0	Expr	S0 component of the shape derivative.	required
S1	Expr	S1 component of the shape derivative.	required

Source code in code/formopt.py

def __init__(
    self,
    dim: int,
    domain: Mesh,
    space: FunctionSpace,
    biform: Callable[[Argument, Argument], Tuple[Expr, bool]],
    S0: Expr,
    S1: Expr,
) -> None:
    """
    Sets up the linear system for the velocity.
    Since the bilinear part remains unchanged,
    it is compiled here.

    Parameters
    ----------
    dim : int
        Domain dimension.
    domain : Mesh
        Problem domain.
    space : FunctionSpace
        Space of functions.
    biform : Callable[[Argument, Argument], Tuple[Expr, bool]]
        Bilinear form method of a subclass of Model.
    S0 : Expr
        S0 component of the shape derivative.
    S1 : Expr
        S1 component of the shape derivative.
    """

    th = TrialFunction(space)
    xi = TestFunction(space)
    dx = Measure("dx", domain=domain)

    b, dirbc = biform(th, xi)

    self.biform = form(b)
    self.liform = form(-(dot(S0, xi) + inner(S1, grad(xi))) * dx)
    self.bc = None

    if dirbc == True:
        # Homogeneous Dirichlet boundary condition on the whole boundary
        # Recall that domain dimension = velocity rank
        self.bc = homogeneus_boundary(domain, space, dim, dim)
    elif isinstance(dirbc, tuple):
        # Homogeneous Dirichlet boundary condition on a subset of the boundary.
        # dirbc is a tuple with two components,
        # boundary tags and a list of marks. For instance,
        # dirbc = (tags, [mark0, mark1])
        self.bc = homogeneous_dirichlet(domain, space, dirbc[0], dirbc[1], dim)
    else:
        # By defult,
        # homogeneous Neumann boundary condition on the whole boundary
        self.bc = []

    self.solver = build_solver(domain, self.biform, self.bc)

run(theta)

Solves the velocity equation. Only the linear part must be updated.

Source code in code/formopt.py

def run(self, theta: Function) -> None:
    """
    Solves the velocity equation.
    Only the linear part must be updated.
    """

    L = assemble_vector(self.liform)
    apply_lifting(L, [self.biform], [self.bc])
    L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(L, self.bc)

    self.solver.solve(L, theta.x.petsc_vec)
    theta.x.scatter_forward()