AMDiS - Adaptive Multi-Dimensional Simulations

class: center, middle

# Session 6
## Complex geometries and Surface PDEs

---

# Overview

- Local refinement
- Implicit geometry description
- Surface FEM
- Branched meshes

---

# Motivation

Granular structur in cross-sections of paper coating

---

# Motivation

Granular structur in cross-sections of paper coating

---

# Local refinement

- mesh adaption based on a-priori information
- use an indicator-function to select elements for refinement/coarsening

```
// using namespace AMDiS::extensions;
class RefinementExpression
{
public:
  RefinementExpression(Mesh *mesh);

template <class EXPRESSION>
  bool refine(EXPRESSION e);

template <class EXPRESSION>
  bool refine(int numRefinements, EXPRESSION e);

// ...
};
```
with `EXPRESSION -> int` gives the local refinement level.
- Level `0` means: macro mesh level
- Level `10` means: 10 bisections of the element
- Refine level by level until `numRefinements` is reached. Default is `numRefinements=15`.

---

# Examples

- Let `C` be a phase-field functions with values in `\([0,1]\)`.
- Refine the mesh up to level `10` whenever an element is in the region `\(C\in[0.1, 0.9]\)`.
- Everything else should be as coarse as possible.
--

```
RefinementExpression adaption(prob.getMesh());
adaption.refine( func([](double c) -> int
  {
    return  0.1 <= c && c <= 0.9 ? 10 : 0;
  },
  valueOf(C)) );
```
--

- Refinement around a corner of the domain:

```
adaption.refine( func([](WorldVector<double> const& x) -> int
  {
    return two_norm(x) < 0.1 ? 10 : 0;
  },
  X()) );
```

---

# Combination of refinement indicators

- You want to refine at the interface of a phase-field and in a corner?
- Take the maximum of two refinement functions:

```
auto f1 = func([](double c) -> int
               { return  0.1 <= c && c <= 0.9 ? 10 : 0; }, valueOf(C));
auto f2 = func([](WorldVector<double> const& x) -> int
               { return two_norm(x) < 0.1 ? 10 : 0; }, X());

adaption.refine( max(f1, f2) );
```

---

# Diffuse Domain approach

Describe your domain implicitly, using a smeared out phase-field function:

\\[
\Omega := \\{ \mathbf{x} \in \bar{\Omega}\;:\; d(\mathbf{x}) < 0 \\}
\\]
with `\(d(\mathbf{x}) = \pm \operatorname{dist}(\mathbf{x},\partial\Omega)\)` a
signed distance function (negative inside of `\(\Omega\)` and positive outside).

Then a phase-field can be defined as
\\[
\phi(\mathbf{x}) := \frac{1}{2}(1 - \tanh(3 d(\mathbf{x})/\epsilon) )
\\]
with `\(\epsilon\ll 1\)` a smoothing parameter.

---

# Example: Circular domain

The signed distance function is defined by
\\[
d(\mathbf{x}) = ||\mathbf{x}|| - R
\\]
with `\(R\)` the radius of the circle around the coordinate center `\((0,0)\)`.

This can be implemented in AMDiS using predefined functors:
```
#include "SignedDistFunctors.h"
auto* D1 = new Circle(radius); // or
auto* D2 = new Circle(radius, center);
```
The result is an `AbstractFunction<double, WorldVector<double> >`.

---

# Example: polygonal domain

2D Domain is bounded by polygonal chain, given as a vector of coordinates:
```
std::vector< WorldVector<double> > points;
double angle = 2*M_PI/n;
for (size_t i = 0; i < n; ++i) {
  WorldVector<double> x;
  x[0] = std::cos(i*angle); x[1] = std::sin(i*angle);
  points.push_back(x);
}

auto* D3 = new Polygon(points);
```

This polygon can be deformed using a modifier method:
```
D3->move( velocity );
```
where velocity is given as `DOFVector<WorldVector<double> >*` or as `AbstractFunction`.

---

# Create a phase-field from the distance function

- Create a new `AbstractFunction` that represents a phase-field by wrapping the distance function
```
#include "PhaseFieldConvert.h"
auto* P1 = new SignedDistFctToPhaseField(epsilon, D1);
```
- Use EXPRESSIONS as wrappers:
```
auto p1 = 0.5*( 1 - tanh(3*eval(D1) / epsilon) );
```
- Or interpolate distance function and phase-field directly to `DOFVectors`:

```
DOFVector<double> dist1(prob.getFeSpace(), "dist1");
DOFVector<double> phase1(prob.getFeSpace(), "phase1");

dist1 << eval(D1);
phase1 << func( wrap(new SignedDistToPhaseField), valueOf(dist1) ); // or
phase1 << func( wrap(new SignedDistToPhaseField), eval(D1) );
```

---

# Create a distance function from phase-field

The opposite way is more complicated:
- Function `tanh` is rounded to -1, 1.
- Only in the interface region good approximation

Need **redistancing**, i.e. calc a distance function from a levelset-function
```
auto p = max(1.e-12, min(1.0 - 1.e-12, valueOf(phase1)));
dist1 << epsilon/3 * atanh(1 - 2*p) );
```
--

Redistancing only available for `P1` FeSpace:
```
#include "HL_SignedDistTraverse.h"
#include "Recovery.h"

feSpaceP1 = ...;
DOFVector<double> tmp(feSpaceP1, "tmp");
tmp.interpol(*dist1);

HL_SignedDistTraverse reinit("reinit", dim);
reinit.calcSignedDistFct(adaptInfo, &tmp);

Recovery rec(L2_NORM, 1);
rec.recoveryUh(&tmp, *dist1);
// or dist1->interpol(tmp);
```

---

# Calc integrals

To integrate a `DOFVector` over an implicit domain, either multiply with a phase-field
```
integrate( p1 * EXPRESSION ); // or
integrate( valueOf(phase1) * EXPRESSION );
```

or use method for integration over negative distance fucntion:
```
#include "compositeFEM/CFE_Integration.h"

auto* One = new AMDiS::Constant(1.0);
ElementFunctionAnalytic<double> oneFct(One);
ElementFunctionAnalytic<double> elLevelFct(D1);
ElementLevelSet elLevelSet("circle", &elLevelFct, prob.getMesh());

double disk_area = CFE_Integration::integrate_onNegLs(&oneFct, &elLevelSet);
double disk_perimeter = CFE_Integration::integrate_onZeroLs(&oneFct, &elLevelSet);
```

`->` see next session about composite-FEM method

---

# Surface PDEs

.center[
<img src="images/rho_big.png" width="40%" /> <img src="images/Earth_flow.png" width="40%" />
]

---

#  Surface PDEs

- Domain is a surface manifold
- Dimension of mesh `\(\neq\)` dimension of world
- Differential operators are projected to surface `\(\Gamma\)`:

\\[
-\Delta\_\Gamma u = f,\text{ on }\Gamma\subset\mathbb{R}^{n+1}
\\]
or in weak form
\\[
\langle \nabla\_\Gamma u, \nabla\_\Gamma \theta\rangle\_\Gamma = \langle f, \theta\rangle\_\Gamma
\\]
with the tangential grandient `\(\nabla_\Gamma:=\mathbb{P}_\Gamma\nabla\)`

- Using barycentric coordinates, we get
\\[
\Lambda := \begin{pmatrix} - \nabla\lambda\_0 - \\\\ - \nabla\lambda\_1 - \\\\ \cdots \\\\ - \nabla\lambda\_d - \end{pmatrix}\in\mathbb{R}^{dim+1\times dow},\quad
\nabla\_\Gamma\phi(\mathbf{x}) = \Lambda^\top \nabla\_\lambda \phi(\lambda(\mathbf{x}))
\\]

---

#  Surface PDEs

- No change in assembling routine necessary
- Surface mesh can be provided as macro-mesh
```matlab
DIM:          2
DIM_OF_WORLD: 3
...
```

.center[
<img src="images/torus_macro.jpg" width="40%" />
]

---

#  Surface PDEs

- Refinement of surface mesh puts new vertices inside flat triangles!
- **Idea**: surface projection

```
class Projection
{
public:
  Projection(int id, ProjectionType type);

/// Projection method. Must be overriden in sub classes.
  virtual void project(WorldVector<double>& x) = 0;
};
```
with `ProjectionType = [BOUNDARY_PROJECTION|VOLUME_PROJECTION]`.

- Example: `BallProject`:
\\[
\mathbf{x}\mapsto (\mathbf{x}-\mathbf{x}\_0) \frac{1}{||\mathbf{x}-\mathbf{x}\_0||}
\\]

---

#  Surface PDEs

- **Attention!** Projection is stored in boundary faces only. Thus, boundary vertices of elements that do not have a boundary face are not projected:

.center[
<img src="images/local_refinement_0.png" width="40%" />
]

---

# Branched geometries

Outlook: graph-like geometries or branched facets

---

# Branched geometries

Numbering must account for more than one neighbour.

**Idea**: identify left and right neighbours:

<table width="100%"><tr>
<td width="60%"><pre><code>element neighbours:
-1 1
-1 2
-1 0

element neighbours inverse:
-1 2
-1 0
-1 1
</code></pre>
</td><td width="40%" style="text-align:center;">
<img src="images/graph.png" width="100%" alt="A 1d graph mesh with numbers" />
</td>
</tr></table>

**Status in AMDiS**:
- Linear Lagrange elements on fixed mesh works
- Refinement and higher order finite elements do not yet work

---

class: center, middle
# Exercise
## Diffuse-Domain

---

# Exercise: Diffuse-Domain

Solve the diffusion equation in a circular domain:
\\[
  \partial\_t u - \Delta u = f,\quad\text{ in }\Omega,\quad u\big|\_{\partial\Omega} = g
\\]
with `\(\Omega=0.5 \mathcal{B}^2\)` and `\(g\equiv 0, f\equiv 1\)`.

Rewrite this equation, using the Diffuse-Domain approach:
\\[
  \partial\_t (\phi u) - \nabla\cdot(\phi\nabla u) + \frac{1}{\epsilon^3}(1-\phi)(u - g) = \phi f,\quad\text{ in }\bar{\Omega},\quad \partial\_n u\big|\_{\partial\bar\Omega} = 0
\\]
with `\(\bar\Omega=[0,1]\times[0,1]\)`.

1. Define a phase-field that represents `\(\Omega\)`
2. Refine the mesh at the interface of the phase-field
3. solve the diffuse-domain equation in `\(\bar\Omega\)`

---

# Advanced Exercise: Diffuse-Interface

Instead of the phase-field `\(\phi\)`, use the double-well
\\[
B(\phi) := \phi^2(1 - \phi^2)
\\]
and solve the equation in a 3d box:
\\[
  \partial\_t (B(\phi) u) - \nabla\cdot(B(\phi)\nabla u) = B(\phi) f,\quad\text{ in }\bar{\Omega},\quad \partial\_n u\big|\_{\partial\bar\Omega} = 0
\\]
with `\(\bar\Omega=[0,1]\times[0,1]\times[0,1]\)`, where `\(\phi\)` describes a phase-field that represents the interface of a 2d sphere in 3d.

1. Define an expression that represents `\(B(\phi)\)`
2. Add a regularization parameter to B, i.e. `\(B(\phi) +\!\!\!= 10^{-5}\)`
3. Instead of `\(B(\phi)\)` use the gradient norm `\(||\nabla\phi||\)`