CutFEM Method

class cutfem_method.CutFemMethod(level_set, level_set_space, space_displacement, ds, bc, bc_velocity, parameters, shift)

This is the CutFEM class.

Some details about the initialization of linear elasticity problem with CutFEM method.

Definition of Primal problem

Linear elasticity problem is given by: Find \(u:\Omega \rightarrow \mathbb{R}^{d}\)

\[\begin{split}\begin{align} \begin{cases} -\text{div}( \sigma(u)) & \!\!\!\!=0 \text{ in }\Omega\\ u& \!\!\!\!=0\text{ on }\Gamma_{D}\\ \sigma(u)\cdot n & \!\!\!\!=g\text{ on }\Gamma_{N} \end{cases} \end{align}\end{split}\]

Where \(d\) the dimension of the problem. We assume small deformations and zero volumetric forces.

This yields to the following weak formulation: Find \(u \in V\), such that for all \(v \in V\) we have

\[a\left(u,v\right)=l\left(v\right)\]

with:

\[\begin{split}\begin{align} a\left(u,v\right) &= 2\mu\left(\varepsilon(u),\varepsilon(v)\right)_{L^{2}(\Omega)} + \lambda\left(\nabla\cdot u,\nabla\cdot v\right)_{L^{2}(\Omega)} \\ l\left(v\right) &= \left(g,v\right)_{L^{2}\left(\Gamma_{N}\right)}, \end{align}\end{split}\]

Bilinear form (primal):

import ufl

u =ufl.TrialFunction(self.space_displacement)
v =ufl.TestFunction(self.space_displacement)

self.gamma = 1e-5*(self.lame_mu + self.lame_lambda)

self.h = CellDiameter(self.mesh)

self.bc = bc

self.a_primal =  2.0*self.lame_mu  * ufl.inner(mecanics_tool.strain(u), mecanics_tool.strain(v)) * self.dxq \
    + self.lame_lambda *  ufl.inner(ufl.nabla_div(u), ufl.nabla_div(v)) * self.dxq
# Stabilization:
self.a_primal += avg(self.gamma) * avg(self.h)**3*ufl.inner(ufl.jump(ufl.grad(u),self.n),\
    ufl.jump(ufl.grad(v),self.n))*self.dS(0)

Linear form (primal):

self.L_primal = ufl.dot(self.shift,v) * self.ds(2)

Definition of Dual problem

Some details about the initialization of adjoint problem with CutFEM.

Bilinear form (dual):

import ufl

u =ufl.TrialFunction(self.space_displacement)
v =ufl.TestFunction(self.space_displacement)

self.gamma = 1e-5*(self.lame_mu + self.lame_lambda)

self.h = CellDiameter(self.mesh)

self.bc = bc

self.a_adjoint =  2.0*self.lame_mu  * ufl.inner(mecanics_tool.strain(u), mecanics_tool.strain(v)) * self.dxq \
    + self.lame_lambda *  ufl.inner(ufl.nabla_div(u), ufl.nabla_div(v)) * self.dxq
# Stabilization:
self.a_adjoint += avg(self.gamma) * avg(self.h)**3*ufl.inner(ufl.jump(ufl.grad(u),self.n),\
    ufl.jump(ufl.grad(v),self.n))*self.dS(0)

Linear form (dual):

The dual operator is compute using the automatic differentiation :

## Exemple for Lp nom of VonMises constraint minimization:

self.J = ((mechanics_tool.von_mises(self.uh,self.lame_mu,self.lame_lambda,self.dim)/parameters.elasticity_limit)**self.p_const)*self.dxq

self.L_adj = ufl.derivative(self.J,self.uh,v_adj)

adjoint_problem(u, parameters, level_set, adjoint=0)

Resolution of the dual problem with the CutFEM method.

Parameters:

• u (fem.Function) – The displacement field function, \(u_{h}\).

• parameters (Parameters) – The object parameters.

• adjoint (ufl.Expression) – The adjoint operator if needed.

Returns:

The dual solution, \(p_{h}\).

Return type:

fem.Function

cut_fem_adv(level_set, dt, velocity_field)

Resolution of advection equation with CutFEM stabilization.

Parameters:

• level_set (fem.Function) – The level set field wich defined implicitely the domain \(\Omega\).

• dt (float) – The :math: dt time parameters.

• velocity_field (ufl.Expression) – The value of advection velocity_field, in normal direction of the interface \(\partial\Omega\).

Returns:

The values of the advected level set.

Return type:

fem.Function

cutfem_solver(level_set, parameters, problem_topo=0)

Resolution of the primal and dual problem.

Parameters:

• level_set (fem.Function) – The level set field which defined implicitly the domain \(\Omega\).

• parameters (Parameters) – The object parameters.

• adjoint (ufl.Expression) – The adjoint operator if needed.

Returns:

The values of the primal and dual solution.

Return type:

fem.Function, fem.Function

descent_direction(level_set, parameters, rest_constraint, constraint_integrand, cost_integrand, xsi=0)

Determine the descent direction by solving the following equation:

Find \(v'_{\text{reg}}\in H_{\Gamma_{D}}^{1}=\left\{ v\in H^{1}\left(D\right)\text{ such that }v=0\text{ on }\Gamma_{D}\right\}\) such that \(\forall w\in H_{\Gamma_{D}}^{1}\)

\[\alpha\left(\nabla v'_{\text{reg}},\nabla w\right)_{L^{2}\left(D\right)}+\left(v'_{\text{reg}},w\right)_{L^{2}\left(D\right)}=-J'(\Omega)\left(w\right)\]

with \(J\) the cost function and \(\alpha>0\) is a smoothing parameter instantiated in the Parameter class.

Parameters:

• level_set (fem.Function) – The level set field which defined implicitly the domain \(\Omega\).

• parameters (Parameters) – The object parameters.

• rest_constraint (float) – The value of the constraint function \(C(\Omega)\).

• constraint_integrand (fem.Expression) – The integrand of the constraint function.

• cost_integrand (fem.Expression) – The integrand of the cost function.

Returns:

The velocity field, defined in D.

Return type:

fem.Function

euclidean_norm_grad(func)

Calculation of the integrand of the L2-norm of the gradient of the function provided, given by the following equality:

\[\left|\nabla\phi\right| = \sqrt{\nabla\phi\cdot\nabla\phi}\]

Parameters:

func (fem.Function) – Function field.

Returns:

The values of integrand of the L2-norm of the gradient.

Return type:

fem.Expression

primal_problem(level_set, parameters)

Resolution of the primal problem with the CutFEM method.

Parameters:

• level_set (fem.Function) – The level set field which defined implicitly the domain \(\Omega\).

• parameters (Parameters) – The object parameters.

Returns:

The primal solution.

Return type:

fem.Function

set_measure_dxq(level_set)

Set the measure dxq on \(\Omega\).

Parameters:

level_set (fem.Function) – The level_set field function, \(\phi\).

velocity_normalization(v, c)

Normalization of the Velocity field according to the following equation:

\[\overline{v} = \frac{v}{\sqrt{c\left\Vert \nabla\phi\right\Vert _{L^{2}\left(D\right)}^{2}+\left\Vert \phi\right\Vert _{L^{2}\left(D\right)}^{2}}}\]

with \(c>0\) and \(\left\Vert . \right\Vert _{L^{2}\left(D\right)}\) norm defined as:

\[\left\Vert f \right\Vert _{L^{2}\left(D\right)}^{2} = \int_{D} f \cdot f \text{ }dx.\]

Parameters:

• v (fem.Expression or fem.Function) – The scalar velocity field which defined the value of advection in direction of the normal to \(\partial\Omega\).

• c (float) – Value of the smoothing for the velocity normalization. Topically, this value is equal to the smoothing value in the extension equation.

Returns:

The normalized velocity field, defined in D.

Return type:

fem.Expression