fractions.h File Reference

#include "geometry.h"
#include "myc2d.h"

Include dependency graph for fractions.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Macros
#define	intersection(a, b)   min(a,b)
#define	union(a, b)   max(a,b)
#define	difference(a, b)   min(a,-(b))

Functions
auto macro coord	interface_normal (Point p, scalar c)
	By default the interface normal is computed using the MYC approximation.
void	fraction_refine (Point point, scalar c)
static void	alpha_refine (Point point, scalar alpha)
trace void	fractions (scalar Phi, scalar c, vector s={0}, double val=0.)
macro	fraction (scalar f, double func)
	The convenience macros below can be used to define volume and surface fraction fields directly from a function.
macro	solid (scalar cs, vector fs, double func)
coord	youngs_normal (Point point, scalar c)
coord	facet_normal (Point point, scalar c, vector s)
trace void	reconstruction (const scalar c, vector n, scalar alpha)
trace void	output_facets (scalar c, FILE *fp=stdout, vector s={{-1}})
trace double	interface_area (scalar c)
	for (int _d=0;_d< 2;_d++) static double interface_fraction_x(coord m
	if (m.y< 0.)
return	clamp ((alpha - m.x *xo)/m.y, 0., 1.)
trace void	face_fraction (scalar f, vector s, double tol=1.e-6)
	This function calculates a single surface fraction value in a given face.

Variables
	attribute
	Finally, we also need to prolongate the reconstructed value of \(\alpha\).
double	alpha = 0.5*(m.x + m.y)
double bool	right
double	xo = right ? 1. : 0.

Macro Definition Documentation

difference

#define difference	(		a,
			b
	)		min(a,-(b))

Definition at line 389 of file fractions.h.

intersection

#define intersection	(		a,
			b
	)		min(a,b)

Boolean operations

Implicit surface representations have the advantage of allowing simple constructive solid geometry operations.

Definition at line 387 of file fractions.h.

union

#define union	(		a,
			b
	)		max(a,b)

Definition at line 388 of file fractions.h.

Function Documentation

alpha_refine()

static void alpha_refine ( Point  point, scalar  alpha  )

static

Definition at line 86 of file fractions.h.

References alpha, dimension, m(), n, x, and coord::x.

Referenced by reconstruction().

Here is the call graph for this function:

Here is the caller graph for this function:

clamp()

return clamp	(	(alpha - m.x *xo)/m.	y,
		0.	,
		1.
	)

face_fraction()

trace void face_fraction	(	scalar	f,
		vector	s,
		double	tol = 1.e-6
	)

This function calculates a single surface fraction value in a given face.

The PLIC interface representation cannot guarantee the continuity of the planar segments across the faces. Therefore, we obtain a single value using a geometric mean between the left and right sides. The tolerance defines the interfacial cells and it can be modified.

Note that adaptation of the face fraction s is currently not performed through a dedicated refinement procedure. Consequently, s must be recalculated after each grid adaptation and update of the volume fraction.

We calculate the interface normal field.

Case 1: the face is shared between a full and an empty cell. The surface fraction is null.

Case 2: if both cells are full, the surface fraction is unitary.

Case 3: if both cells are interfacial, the face contains the interface, and we proceed with the calculation of the surface fraction.

Definition at line 628 of file fractions.h.

References _i, alpha, dimension, f, m(), mycs(), n, plane_alpha, point, s, x, and coord::x.

Here is the call graph for this function:

facet_normal()

coord facet_normal	(	Point	point,
		scalar	c,
		vector	s
	)

Normal approximation using MYC or face fractions

Definition at line 426 of file fractions.h.

References c, dimension, scalar::i, interface_normal(), n, nn, point, s, x, and coord::x.

Referenced by embed_flux(), embed_fraction_refine(), embed_geometry(), and refine_face_x_axi().

Here is the call graph for this function:

Here is the caller graph for this function:

for()

for

(

)

Face fraction

We wish to calculate the fraction of face surface occupied by a phase defined by a volume fraction field. This operation can be useful in different contexts, for example the solution of the diffusion equation in a specific phase.

The surface fraction is computed as the intersection between the face of the cell under investigation and the PLIC interface fragment. The problem boils down to the calculation of the intersection between the interfacial plane and the coordinate of the cell face ( \(x = x_o\)). In 2D:

\[ y = \frac{\alpha - m_x x_o}{m_y} \]

after shifting the reference frame and handling degenerate cases ( \(m_y = 0\)). In 3D, this concept is extended by computing the area occupied by the phase. The plane is considered either on the left or on the right side of the face, as controlled using the boolean right.

fraction()

macro fraction	(	scalar	f,
		double	func
	)

The convenience macros below can be used to define volume and surface fraction fields directly from a function.

Definition at line 359 of file fractions.h.

References _i, f, fractions(), func, phi, and x.

Here is the call graph for this function:

fraction_refine()

void fraction_refine	(	Point	point,
		scalar	c
	)

Coarsening and refinement of a volume fraction field

On trees, we need to define how to coarsen (i.e. "restrict") or refine (i.e. "prolongate") interface definitions (see [geometry.h]() for a basic explanation of how interfaces are defined).

If the parent cell is empty or full, we just use the same value for the fine cell.

Otherwise, we reconstruct the interface in the parent cell.

And compute the volume fraction in the quadrant of the coarse cell matching the fine cells. We use symmetries to simplify the combinations.

Definition at line 42 of file fractions.h.

References a, alpha, b, c, cc, dimension, mycs(), n, nc, plane_alpha, point, rectangle_fraction(), x, and coord::x.

Referenced by event_acceleration(), event_defaults(), event_metric(), and event_properties().

Here is the call graph for this function:

Here is the caller graph for this function:

fractions()

trace void fractions	(	scalar	Phi,
		scalar	c,
		vector	s = {0},
		double	val = 0.
	)

Computing volume fractions from a "levelset" function

Initialising a volume fraction field representing an interface is not trivial since it involves the numerical evaluation of surface integrals.

Here we define a function which allows the approximation of these surface integrals in the case of an interface defined by a "levelset" function \(\Phi\) sampled on the vertices of the grid.

By convention the "inside" of the interface corresponds to \(\Phi > 0\).

The function takes the scalar field \(\Phi\) as input and fills c with the volume fraction and, optionally if it is given, s with the surface fractions i.e. the fractions of the faces of the cell which are inside the interface.

We store the positions of the intersections of the surface with the edges of the cell in vector field p. In two dimensions, this field is just the transpose of the line fractions s, in 3D we need to allocate a new field.

Line fraction computation

We start by computing the line fractions i.e. the (normalised) lengths of the edges of the cell within the surface.

If the values of \(\Phi\) on the vertices of the edge have opposite signs, we know that the edge is cut by the interface.

In that case we can find an approximation of the interface position by simple linear interpolation. We also check the sign of one of the vertices to orient the interface properly.

If the values of \(\Phi\) on the vertices of the edge have the same sign (or are zero), then the edge is either entirely outside or entirely inside the interface. We check the sign of both vertices to treat limit cases properly (when the interface intersects the edge exactly on one of the vertices).

Surface fraction computation

We can now compute the surface fractions. In 3D they will be computed for each face (in the z, x and y directions) and stored in the face field s. In 2D the surface fraction in the z-direction is the volume fraction c.

We first compute the normal to the interface. This can be done easily using the line fractions. The idea is to compute the circulation of the normal along the boundary \(\partial\Omega\) of the fraction of the cell \(\Omega\) inside the interface. Since this is a closed curve, we have

\[ \oint_{\partial\Omega}\mathbf{n}\;dl = 0 \]

We can further decompose the integral into its parts along the edges of the square and the part along the interface. For the case pictured above, we get for one component (and similarly for the other)

\[ - s_x[] + \oint_{\Phi=0}n_x\;dl = 0 \]

If we now define the average normal to the interface as

\[ \overline{\mathbf{n}} = \oint_{\Phi=0}\mathbf{n}\;dl \]

We have in the general case

\[ \overline{\mathbf{n}}_x = s_x[] - s_x[1,0] \]

and

\[ |\overline{\mathbf{n}}| = \oint_{\Phi=0}\;dl \]

Note also that this average normal is exact in the case of a linear interface.

If the norm is zero, the cell is full or empty and the surface fraction is identical to one of the line fractions.

Otherwise we are in a cell containing the interface. We first normalise the normal.

To find the intercept \(\alpha\), we look for edges which are cut by the interface, find the coordinate \(a\) of the intersection and use it to derive \(\alpha\). We take the average of \(\alpha\) for all intersections.

Once we have \(\mathbf{n}\) and \(\alpha\), the (linear) interface is fully defined and we can compute the surface fraction using our pre-defined function. For marginal cases, the cell is full or empty (ni == 0) and we look at the line fractions to decide.

Volume fraction computation

To compute the volume fraction in 3D, we use the same approach.

Definition at line 123 of file fractions.h.

Referenced by event_init(), fraction(), levelset_to_vof(), solid(), and zarea().

Here is the caller graph for this function:

if()

if

(

)

Definition at line 596 of file fractions.h.

References alpha, and m().

Here is the call graph for this function:

interface_area()

trace double interface_area

(

scalar

c

)

Interfacial area

This function returns the surface area of the interface as estimated using its VOF reconstruction.

Definition at line 553 of file fractions.h.

References alpha, area(), c, dimension, interface_normal(), n, p, plane_alpha, plane_area_center, point, pow(), and x.

Here is the call graph for this function:

interface_normal()

auto macro coord interface_normal	(	Point	p,
		scalar	c
	)

By default the interface normal is computed using the MYC approximation.

Volume fractions

These functions are used to maintain or define volume and surface fractions either from an initial geometric definition or from an existing volume fraction field.

We will use basic geometric functions for square cut cells and the "Mixed-Youngs-Centered" (MYC) normal approximation of Ruben Scardovelli. This can be overloaded by redefining this macro.

Definition at line 29 of file fractions.h.

References c, mycs(), and p.

Referenced by facet_normal(), interface_area(), and reconstruction().

Here is the call graph for this function:

Here is the caller graph for this function:

output_facets()

trace void output_facets	(	scalar	c,
		FILE *	fp = stdout,
		vector	s = {{-1}}
	)

Interface output

This function "draws" interface facets in a file. The segment endpoints are defined by pairs of coordinates. Each pair of endpoints is separated from the next pair by a newline, so that the resulting file is directly visualisable with gnuplot.

The input parameters are a volume fraction field c, an optional file pointer fp (which defaults to stdout) and an optional face vector field s containing the surface fractions.

If s is specified, the surface fractions are used to compute the interface normals which leads to a continuous interface representation in most cases. Otherwise the interface normals are approximated from the volume fraction field, which results in a piecewise continuous (i.e. geometric VOF) interface representation.

Definition at line 518 of file fractions.h.

reconstruction()

trace void reconstruction	(	const scalar	c,
		vector	n,
		scalar	alpha
	)

Interface reconstruction

The reconstruction function takes a volume fraction field c and returns the corresponding normal vector field n and intercept field \(\alpha\).

If the cell is empty or full, we set \(\mathbf{n}\) and \(\alpha\) only to avoid using uninitialised values in alpha_refine().

Otherwise, we compute the interface normal using the Mixed-Youngs-Centered scheme, copy the result into the normal field and compute the intercept \(\alpha\) using our predefined function.

On a tree grid, for the normal to the interface, we don't use any interpolation from coarse to fine i.e. we use straight "injection".

We set our refinement function for alpha.

Definition at line 454 of file fractions.h.

References _i, alpha, alpha_refine(), c, dimension, interface_normal(), m(), n, plane_alpha, point, refine_injection(), and x.

Here is the call graph for this function:

solid()

macro solid	(	scalar	cs,
		vector	fs,
		double	func
	)

Definition at line 369 of file fractions.h.

References _i, cs, fractions(), fs, func, phi, and x.

Here is the call graph for this function:

youngs_normal()

coord youngs_normal	(	Point	point,
		scalar	c
	)

Interface reconstruction from volume fractions

The reconstruction of the interface geometry from the volume fraction field requires computing an approximation to the interface normal.

Youngs normal approximation

This a simple, but relatively inaccurate way of approximating the normal. It is simply a weighted average of centered volume fraction gradients. We include it as an example but it is not used.

Definition at line 403 of file fractions.h.

References assert, c, dimension, n, nn, and x.

Variable Documentation

alpha

alpha = 0.5*(m.x + m.y)

Definition at line 587 of file fractions.h.

Referenced by alpha_refine(), face_fraction(), fraction_refine(), if(), interface_area(), and reconstruction().

attribute

attribute

Initial value:

{
  vector n

Finally, we also need to prolongate the reconstructed value of \(\alpha\).

This is done with the simple formula below. We add an attribute so that we can access the normal from the refinement function.

Definition at line 82 of file fractions.h.

right

double bool right

Initial value:

{
 
 
 
  if (fabs (m.y) < 1e-4) 
    return right ? (m.x < 0. ? 1. : 0.)
                 : (m.x > 0. ? 1. : 0.)

Definition at line 587 of file fractions.h.

xo

double xo = right ? 1. : 0.

Definition at line 600 of file fractions.h.