curvature.h File Reference

#include "heights.h"
#include "parabola.h"

Include dependency graph for curvature.h:

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

Go to the source code of this file.

Data Structures
struct	cstats
	The function below computes the mean curvature kappa of the interface defined by the volume fraction c. More...

Functions
static void	curvature_restriction (Point point, scalar kappa)
static void	curvature_prolongation (Point point, scalar kappa)
	The prolongation function performs a similar averaging, but using the same stencil as that used for bilinear interpolation, so that the symmetries of the volume fraction field and curvature field are preserved.
	for (int _d=0;_d< 2;_d++) static double kappa_y(Point point
	for (int i=-1;i<=1;i++) if(h.y[i]
return hxx	pow (1.+sq(hx), 3/2.)
	if (h.y[]==nodata) return n
	if (h.y[-1] !=nodata &&orientation(h.y[-1])==ori)
else	if (h.y[1] !=nodata &&orientation(h.y[1])==ori) n.x
static double	height_curvature (Point point, scalar c, vector h)
	We now need to choose one of the \(x\), \(y\) or \(z\) height functions to compute the curvature.
coord	height_normal (Point point, scalar c, vector h)
	The function below works in a similar manner to return the normal estimated using height-functions (or a nodata vector if this cannot be done).
static int	independents (coord *p, int n)
	In three dimensions, these functions return the (two) components of the normal projected onto the $(x,y)$ plane (respectively).
static double	height_curvature_fit (Point point, scalar c, vector h)
	Given a volume fraction field c and a height function field h, this function returns the "mixed heights" parabola-fitted curvature (or nodata if the curvature cannot be computed).
static double	centroids_curvature_fit (Point point, scalar c)
static bool	interfacial (Point point, scalar c)
trace cstats	curvature (scalar c, scalar kappa, double sigma=1.[0], bool add=false)
static double	height_position (Point point, scalar f, vector h, coord *G, coord *Z)
	We now need to choose one of the \(x\), \(y\) or \(z\) height functions to compute the position.
void	position (scalar f, scalar pos, coord G={0}, coord Z={0}, bool add=false)
	The position() function fills field pos with.

Variables
vector	h
double	hx = (h.y[1] - h.y[-1])/2.
double	hxx = (h.y[1] + h.y[-1] - 2.*h.y[])/Delta
int	ori = orientation(h.y[])
else return	n = 0
double	nn = (ori ? -1. : 1.)*sqrt(1. + sq(n.x))
n	x = nn
n	y = 1./nn
vector coord *	G
vector coord coord *	Z
coord	o = {x, y, z}
double	pos = 0.

Function Documentation

centroids_curvature_fit()

static double centroids_curvature_fit ( Point  point, scalar  c  )

static

Parabolic fit of centroids

If all else fails, we try a parabolic fit of interface centroids.

We recover the interface normal and the centroid of the interface fragment and initialize the parabolic fit.

We add the interface centroids in a \(3^d\) neighborhood and compute the curvature.

Definition at line 457 of file curvature.h.

References alpha, area(), c, dimension, kappa, m(), mycs(), parabola_fit_add(), parabola_fit_curvature(), parabola_fit_init(), parabola_fit_solve(), plane_alpha, plane_area_center, point, coord::x, x, y, and z.

Referenced by curvature().

Here is the call graph for this function:

Here is the caller graph for this function:

curvature()

trace cstats curvature	(	scalar	c,
		scalar	kappa,
		double	sigma = 1.[0],
		bool	add = false
	)

On trees we set the prolongation and restriction functions for the curvature.

We first compute a temporary curvature k: a "clone" of \(\kappa\).

If we are not in an interfacial cell, we set \(\kappa\) to nodata.

Otherwise we try the standard HF curvature calculation first, and the "mixed heights" HF curvature second.

We then construct the final curvature field using either the computed temporary curvature...

...or the average of the curvatures in the \(3^{d}\) neighborhood of interfacial cells.

Empty neighborhood: we try centroids as a last resort.

We add or set kappa.

Definition at line 543 of file curvature.h.

References _i, a, c, centroids_curvature_fit(), clamp(), curvature_prolongation(), curvature_restriction(), h, height_curvature(), height_curvature_fit(), heights(), interfacial(), k, kappa, nodata, p, point, scalar_clone, sf, sigma, sign(), and x.

Referenced by event_acceleration(), and event_end_timestep().

Here is the call graph for this function:

Here is the caller graph for this function:

curvature_prolongation()

static void curvature_prolongation ( Point  point, scalar  kappa  )

static

The prolongation function performs a similar averaging, but using the same stencil as that used for bilinear interpolation, so that the symmetries of the volume fraction field and curvature field are preserved.

Definition at line 29 of file curvature.h.

References coarse(), dimension, i, j, k, kappa, nodata, s, and x.

Referenced by curvature().

Here is the call graph for this function:

Here is the caller graph for this function:

curvature_restriction()

static void curvature_restriction ( Point  point, scalar  kappa  )

static

Curvature of an interface

The curvature field is defined only in interfacial cells. In all the other cells it takes the value nodata.

On trees, we need to redefine the restriction function to take this into account i.e. the curvature of the parent cell is the average of the curvatures in the interfacial child cells.

Definition at line 14 of file curvature.h.

References k, kappa, nodata, s, and x.

Referenced by curvature().

Here is the caller graph for this function:

for() [1/2]

for

(

)

Height-function curvature and normal

To compute the curvature, we estimate the derivatives of the height functions in a given direction (x, y or z). We first check that all the heights are defined and that their orientations are the same. We then compute the curvature as

\[ \kappa = \frac{h_{xx}}{(1 + h_x^2)^{3/2}} \]

in two dimensions, or

\[ \kappa = \frac{h_{xx}(1 + h_y^2) + h_{yy}(1 + h_x^2) - 2h_{xy}h_xh_y} {(1 + h_x^2 + h_y^2)^{3/2}} \]

in three dimensions.

The normal is computed in a similar way, but also allowing for asymmetric 2-points stencils and taking into account the orientation.

Position of an interface

This is similar to curvature but this time for the position of the interface, defined as

\[ pos = \mathbf{G}\cdot(\mathbf{x} - \mathbf{Z}) \]

with \(\mathbf{G}\) and \(\mathbf{Z}\) two vectors and \(\mathbf{x}\) the coordinates of the interface.

This is defined only in interfacial cells. In all the other cells it takes the value nodata.

We first need a function to compute the position \(\mathbf{x}\) of an interface. For accuracy, we first try to use height functions.

for() [2/2]

for

(

int

i = -1;i<=1;i++

)

height_curvature()

static double height_curvature ( Point  point, scalar  c, vector  h  )

static

We now need to choose one of the \(x\), \(y\) or \(z\) height functions to compute the curvature.

This is done by the function below which returns the HF curvature given a volume fraction field c and a height function field h.

We first define pairs of normal coordinates n (computed by simple differencing of c) and corresponding HF curvature function kappa (defined above).

We sort these pairs in decreasing order of \(|n|\).

We try each curvature function in turn.

We limit the maximum curvature to \(1/\Delta\).

We add the axisymmetric curvature if necessary.

Definition at line 183 of file curvature.h.

References c, dimension, h, height(), hx, kappa, max, n, nodata, nr, orientation(), point, sign(), sq(), swap, vector, vector::x, x, vector::y, y, and z.

Referenced by curvature().

Here is the call graph for this function:

Here is the caller graph for this function:

height_curvature_fit()

static double height_curvature_fit ( Point  point, scalar  c, vector  h  )

static

Given a volume fraction field c and a height function field h, this function returns the "mixed heights" parabola-fitted curvature (or nodata if the curvature cannot be computed).

The coordinates of the interface points and the number of interface points.

We collect the points along all directions.

We don't want to mix heights with different orientations. We first find the "dominant" orientation ori.

We look for height-functions with the dominant orientation and store the corresponding interface coordinates (relative to the center of the cell and normalised by the cell size).

If we don't have enough independent points, we cannot do the parabolic fit.

We recover the interface normal and the centroid of the interface fragment and initialize the parabolic fit.

We add the collected interface positions and compute the curvature.

Definition at line 364 of file curvature.h.

References alpha, area(), c, dimension, h, height(), i, independents(), j, kappa, m(), mycs(), n, nodata, ori, orientation(), parabola_fit_add(), PARABOLA_FIT_CENTER_WEIGHT, parabola_fit_curvature(), parabola_fit_init(), parabola_fit_solve(), plane_alpha, plane_area_center, point, x, vector::y, y, and z.

Referenced by curvature().

Here is the call graph for this function:

Here is the caller graph for this function:

height_normal()

coord height_normal	(	Point	point,
		scalar	c,
		vector	h
	)

The function below works in a similar manner to return the normal estimated using height-functions (or a nodata vector if this cannot be done).

We first define pairs of normal coordinates n (computed by simple differencing of c) and corresponding normal function normal (defined above).

We sort these pairs in decreasing order of \(|n|\).

We try each normal function in turn.

Definition at line 261 of file curvature.h.

References c, dimension, h, n, nodata, normal, point, swap, vector, coord::x, x, y, and z.

Referenced by height_myc_normal().

Here is the caller graph for this function:

height_position()

static double height_position ( Point  point, scalar  f, vector  h, coord *  G, coord *  Z  )

static

We now need to choose one of the \(x\), \(y\) or \(z\) height functions to compute the position.

This is done by the function below which returns the HF position given a volume fraction field f, a height function field h and vectors G and Z.

We first define pairs of normal coordinates n (computed by simple differencing of f) and corresponding HF position function pos (defined above).

We sort these pairs in decreasing order of \(|n|\).

We try each position function in turn.

Definition at line 686 of file curvature.h.

References dimension, f, G, h, n, nodata, point, pos, swap, vector, x, y, z, and Z.

if() [1/3]

if

(

h.y !

[-1] = nodata && orientation(h.y[-1]) == ori

)

Definition at line 92 of file curvature.h.

References h, n, nodata, ori, orientation(), and vector::y.

Here is the call graph for this function:

if() [2/3]

else if

(

h.y !

[1] = nodata &&orientation(h.y[1])==ori

)

if() [3/3]

if

(

h.

y[] = =nodata

)

independents()

static int independents ( coord *  p, int  n  )

static

In three dimensions, these functions return the (two) components of the normal projected onto the $(x,y)$ plane (respectively).

Parabolic fit of "mixed" height-functions

When the standard height function curvature calculation is not possible (for example because not enough heights are available in any given direction), one can try to combine all the available heights (thus using "mixed" directions) to obtain points on the interface. These point locations can then be fitted with a parabola (using least-mean-square optimisation) and the resulting curvature can be computed. The fitting functions are defined in the file included below. Given n (interface) point coordinates, this function returns the number of "independent" points i.e. points which are more than half-a-cell distant from all the other points.

Definition at line 341 of file curvature.h.

References dimension, i, j, n, p, sq(), and x.

Referenced by height_curvature_fit().

Here is the call graph for this function:

Here is the caller graph for this function:

interfacial()

static bool interfacial ( Point  point, scalar  c  )

inlinestatic

General curvature computation

We first need to define "interfacial cells" i.e. cells which contain an interface. A simple test would just be that the volume fraction is neither zero nor one. As usual things are more complicated because of round-off errors. They can cause the interface to be exactly aligned with cell boundaries, so that cells on either side of this interface have fractions exactly equal to zero or one. The function below takes this into account.

Definition at line 506 of file curvature.h.

References c, dimension, i, and x.

Referenced by curvature(), and dphidt().

Here is the caller graph for this function:

position()

void position	(	scalar	f,
		scalar	pos,
		coord	G = {0},
		coord	Z = {0},
		bool	add = false
	)

The position() function fills field pos with.

\[ \mathbf{G}\cdot(\mathbf{x} - \mathbf{Z}) \]

with \(\mathbf{x}\) the position of the interface defined by \(f\).

If add is true, the position is added to pos.

On trees we set the prolongation and restriction functions for the position.

If the height function is not defined, we use the centroid of the reconstructed VOF interface.

Definition at line 737 of file curvature.h.

pow()

return hxx pow	(	1.+	sqhx,
		3/	2.
	)

Referenced by embed_force(), evaluate_expression(), event_face_fields(), geometric_beta(), interface_area(), parabola_fit_add(), parabola_fit_curvature(), and remove_droplets().

Here is the caller graph for this function:

Variable Documentation

G

vector coord* G

Definition at line 668 of file curvature.h.

Referenced by height_position().

h

vector h

Initial value:

{
  int ori = orientation(h.y[])

Definition at line 74 of file curvature.h.

Referenced by curvature(), height_curvature(), height_curvature_fit(), height_normal(), height_position(), and if().

hx

double hx = (h.y[1] - h.y[-1])/2.

Definition at line 80 of file curvature.h.

Referenced by de(), height_curvature(), and parabola_fit_curvature().

hxx

double hxx = (h.y[1] + h.y[-1] - 2.*h.y[])/Delta

Definition at line 81 of file curvature.h.

Referenced by parabola_fit_curvature().

n

u n = 0

Definition at line 101 of file curvature.h.

Referenced by a32_hash_add(), alpha_refine(), area_integral(), args(), assemble_node(), barycenter(), cfilter(), de(), debug_mpi(), distance(), dtnext(), embed_fraction_refine(), eta_restriction(), evaluate_expression(), event_acceleration(), event_face_fields(), event_properties(), event_register(), face_fraction(), facet_normal(), facets(), falsepos(), fine(), for(), for(), fraction_refine(), glvertex_normal3d(), gpu_errors_scan_bytes(), h_relax(), heap_read(), height_curvature(), height_curvature_fit(), height_myc_normal(), height_normal(), height_position(), if(), if(), if(), if(), if(), inc_scan_bytes(), independents(), init_grid(), init_solver(), input0(), input_init_input(), input_netcdf_init_input_netcdf(), input_pgm(), input_read_obs(), interface_area(), interpolate_array(), line_alpha(), line_center(), line_length_center(), log_base2(), main(), matrix_inverse(), matrix_new(), maxruntime(), mpi_boundary_update_buffers(), mtridiagonal_init_tridiagonal(), mtridiagonal_tridiagonal(), mycs(), normalize(), normf(), output_gauges(), output_vtk(), parabola_fit_init(), parse_params(), PointSegmentOrientation(), PointTriangleOrientation(), post_scan_bytes(), qpopen(), quad_neighbor(), quad_neighbor_finest(), quad_x(), quad_y(), query(), query_sum(), reconstruction(), rectangle_fraction(), refine_cm_axi(), refine_face_x_axi(), relax(), relax_hydro(), relax_nh(), remove_droplets(), right_value(), RPE(), sizes(), smatrix_inverse(), split(), tag(), terrain(), time_update_time(), timer_timing(), update_distance(), update_sum(), util_adv_center(), util_diff_center(), util_diff_face(), util_matrix(), youngs_normal(), zarea(), and zvolume().

nn

double nn = (ori ? -1. : 1.)*sqrt(1. + sq(n.x))

Definition at line 102 of file curvature.h.

Referenced by facet_normal(), kpp_bottom_layer(), kpp_do_kpp(), kpp_interior(), kpp_surface_layer(), observations_const_nns(), observations_const_nnt(), turbulence_algebraiclength(), turbulence_alpha_mnb(), turbulence_dissipationeq(), turbulence_do_epsb(), turbulence_do_kb(), turbulence_do_lengthscale(), turbulence_do_tke(), turbulence_do_turbulence(), turbulence_genericeq(), turbulence_internal_wave(), turbulence_ispralength(), turbulence_kbeq(), turbulence_lengthscaleeq(), turbulence_potentialml(), turbulence_production(), turbulence_q2over2eq(), turbulence_tkealgebraic(), turbulence_tkeeq(), update_distance(), and youngs_normal().

o

coord o = {x, y, z}

Definition at line 672 of file curvature.h.

Referenced by advance_generic(), event_viscous_term(), neighbor(), parabola_fit_init(), periodic_function(), relative(), restriction_embed_linear(), stripslash(), and sum_add_point().

ori

int ori = orientation(h.y[])

Definition at line 91 of file curvature.h.

Referenced by height_curvature_fit(), and if().

pos

return pos = 0.

Definition at line 674 of file curvature.h.

Referenced by height_position(), if(), query(), query_sum(), and split().

x

o x = nn

Definition at line 103 of file curvature.h.

Referenced by centroids_curvature_fit(), curvature(), curvature_prolongation(), curvature_restriction(), height_curvature(), height_curvature_fit(), height_normal(), height_position(), independents(), and interfacial().

y

n y = 1./nn

Definition at line 104 of file curvature.h.

Referenced by centroids_curvature_fit(), height_curvature(), height_curvature_fit(), height_normal(), and height_position().

Z

vector coord coord* Z

Initial value:

{
  if (fabs(height(h.x[])) > 1.)
    return nodata

Definition at line 668 of file curvature.h.

Referenced by height_position().