This 2D example program solves the Poisson equation using finite elements. More...

#include <p4est_bits.h>
#include <p4est_ghost.h>
#include <p4est_lnodes.h>
#include <p4est_vtk.h>

Include dependency graph for p4est_step4.c:

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

Functions
static int	refine_fn (p4est_t p4est, p4est_topidx_t which_tree, p4est_quadrant_t quadrant)
	Callback function to decide on refinement.
static double	func_rhs (const double vxyz[3])
	Right hand side function for the 2D Poisson problem.
static double	func_uexact (const double vxyz[3])
	Exact solution for the 2D Poisson problem.
static int	is_boundary_unitsquare (p4est_t p4est, p4est_topidx_t tt, p4est_quadrant_t node)
	Determine the boundary status on the unit square.
static int	lnodes_decode2 (p4est_lnodes_code_t face_code, int hanging_corner[P4EST_CHILDREN])
	Decode the information from p4est_lnodes_t for a given element.
static void	share_sum (p4est_t p4est, p4est_lnodes_t lnodes, double *v)
	Parallel sum of values in node vector across all sharers.
static double	vector_dot (p4est_t p4est, p4est_lnodes_t lnodes, const double v1, const double v2)
	Compute the inner product of two node vectors in parallel.
static void	vector_axpy (p4est_t p4est, p4est_lnodes_t lnodes, double a, const double x, double y)
	Compute y := y + a * x.
static void	vector_xpby (p4est_t p4est, p4est_lnodes_t lnodes, const double x, double b, double y)
	Compute y := x + b * y.
static void	vector_zero (p4est_t p4est, p4est_lnodes_t lnodes, double *x)
	Zero all entries of a vector.
static void	vector_copy (p4est_t p4est, p4est_lnodes_t lnodes, const double x, double y)
	copy a vector.
static double *	allocate_vector (p4est_lnodes_t *lnodes)
	Allocate storage for processor-relevant nodal degrees of freedom.
static void	interpolate_functions (p4est_t p4est, p4est_lnodes_t lnodes, double rhs_eval, double uexact_eval, int8_t **pbc)
	Interpolate right hand side and exact solution onto mesh nodes.
static void	multiply_matrix (p4est_t p4est, p4est_lnodes_t lnodes, const int8_t bc, int stiffness, double(matrix)[P4EST_CHILDREN][P4EST_CHILDREN], const double in, double out)
	Apply a finite element matrix to a node vector, y = Mx.
static void	set_dirichlet (p4est_lnodes_t lnodes, const int8_t bc, double *v)
	Set Dirichlet boundary values of a node vector to zero.
static void	test_area (p4est_t p4est, p4est_lnodes_t lnodes, double(matrix)[P4EST_CHILDREN][P4EST_CHILDREN], double tmp, double *lump)
	Multiply the mass matrix with a vector of ones to compute the volume.
static void	solve_by_cg (p4est_t p4est, p4est_lnodes_t lnodes, const int8_t bc, int stiffness, double(matrix)[P4EST_CHILDREN][P4EST_CHILDREN], const double b, double x)
	Execute the conjugate gradient method.
static void	solve_poisson (p4est_t *p4est)
	Execute the numerical part of the example: Solve Poisson's equation.
int	main (int argc, char **argv)
	The main function of the step4 example program.

Detailed Description

This 2D example program solves the Poisson equation using finite elements.

Currently, it works on the unit square. For more general domains, a coordinate transformation to physical space would need to be implemented. This will usually entail using a quadrature instead of exact integration. The check for boundary nodes would need to be adapted to the new geometry.

Function Documentation

static double* allocate_vector ( p4est_lnodes_t * lnodes ) [static]

Allocate storage for processor-relevant nodal degrees of freedom.

Parameters:

[in] lnodes This structure is queried for the node count.

Returns:: Allocated double array; must be freed with P4EST_FREE.

static double func_rhs ( const double vxyz[3] ) [static]

Right hand side function for the 2D Poisson problem.

This is the negative Laplace operator acting on the function uexact.

Parameters:

[in] vxyz x, y, and z coordinates in physical space.

Returns:: Scalar function value at vxyz.

static double func_uexact ( const double vxyz[3] ) [static]

Exact solution for the 2D Poisson problem.

We pick a function with zero Dirichlet boundary conditions on the unit square.

Parameters:

[in] vxyz x, y, and z coordinates in physical space.

Returns:: Scalar function value at vxyz.

static void interpolate_functions	(	p4est_t *	p4est,
		p4est_lnodes_t *	lnodes,
		double **	rhs_eval,
		double **	uexact_eval,
		int8_t **	pbc
	)		`[static]`

Interpolate right hand side and exact solution onto mesh nodes.

Parameters:

[in]	p4est	The forest is not changed.
[in]	lnodes	The node numbering is not changed.
[out]	rhs_eval	Is allocated and filled with function values.
[out]	uexact_eval	Is allocated and filled with function values.
[out]	pbc	Boolean flags for Dirichlet boundary nodes.

static int is_boundary_unitsquare	(	p4est_t *	p4est,
		p4est_topidx_t	tt,
		p4est_quadrant_t *	node
	)		`[static]`

Determine the boundary status on the unit square.

Parameters:

[in]	p4est	Can be used to access the connectivity.
[in]	tt	The tree number (always zero for the unit square).
[in]	node	The corner node of an element to be examined.

Returns:: True for Dirichlet boundary, false otherwise.

static int lnodes_decode2	(	p4est_lnodes_code_t	face_code,
		int	hanging_corner[P4EST_CHILDREN]
	)		`[static]`

Decode the information from p4est_lnodes_t for a given element.

See also:: p4est_lnodes.h for an in-depth discussion of the encoding.

Parameters:

[in]	face_code	Bit code as defined in p4est_lnodes.h.
[out]	hanging_corner	Undefined if no node is hanging. If any node is hanging, this contains one integer per corner, which is -1 for corners that are not hanging, and the number of the non-hanging corner on the hanging face otherwise.

Returns:: true if any node is hanging, false otherwise.

int main	(	int	argc,
		char **	argv
	)

The main function of the step4 example program.

It creates a connectivity and forest, refines it, and solves the Poisson equation with piecewise linear finite elements.

static void multiply_matrix	(	p4est_t *	p4est,
		p4est_lnodes_t *	lnodes,
		const int8_t *	bc,
		int	stiffness,
		double(*)	matrix[P4EST_CHILDREN][P4EST_CHILDREN],
		const double *	in,
		double *	out
	)		`[static]`

Apply a finite element matrix to a node vector, y = Mx.

Parameters:

[in]	p4est	The forest is not changed.
[in]	lnodes	The node numbering is not changed.
[in]	bc	Boolean flags for Dirichlet boundary nodes. If NULL, no special action is taken.
[in]	stiffness	If false use scaling for the mass matrix, if true use the scaling for stiffness matrix.
[in]	matrix	A 4x4 matrix computed on the reference element.
[in]	in	Input vector x.
[out]	out	Output vector y = Mx.

static int refine_fn	(	p4est_t *	p4est,
		p4est_topidx_t	which_tree,
		p4est_quadrant_t *	quadrant
	)		`[static]`

Callback function to decide on refinement.

This function is called for every processor-local quadrant in order; its return value is understood as a boolean refinement flag. We refine around a h = 1/8 block with left front lower corner (5/8, 2/8, 6/8).

static void set_dirichlet	(	p4est_lnodes_t *	lnodes,
		const int8_t *	bc,
		double *	v
	)		`[static]`

Set Dirichlet boundary values of a node vector to zero.

Parameters:

[in]	lnodes	The node numbering is not changed.
[in]	bc	Boolean flags for Dirichlet boundary nodes. If NULL, this function does nothing.
[in,out]	v	Dirichlet nodes are overwritten with zero.

static void share_sum	(	p4est_t *	p4est,
		p4est_lnodes_t *	lnodes,
		double *	v
	)		`[static]`

Parallel sum of values in node vector across all sharers.

This function is necessary in the matrix-vector product since elements from multiple processors can contribute to any given node value.

Parameters:

[in]	p4est	The mesh is not changed.
[in]	lnodes	The node numbering is not changed.
[in,out]	v	On input, vector with local contributions. On output, the node-wise sum across all sharers.

static void solve_by_cg	(	p4est_t *	p4est,
		p4est_lnodes_t *	lnodes,
		const int8_t *	bc,
		int	stiffness,
		double(*)	matrix[P4EST_CHILDREN][P4EST_CHILDREN],
		const double *	b,
		double *	x
	)		`[static]`

Execute the conjugate gradient method.

Parameters:

[in]	p4est	The forest is not changed.
[in]	lnodes	The node numbering is not changed.
[in]	bc	Boolean flags for Dirichlet boundary nodes.
[in]	stiffness	If false use scaling for the mass matrix, if true use the scaling for stiffness matrix.
[in]	matrix	The mass matrix should be passed in here.
[in]	b	The right hand side vector.
[out]	x	The result; we use an initial value of zero.

static void solve_poisson ( p4est_t * p4est ) [static]

Execute the numerical part of the example: Solve Poisson's equation.

Parameters:

[in] p4est Solve the PDE with the given mesh refinement.

static void test_area	(	p4est_t *	p4est,
		p4est_lnodes_t *	lnodes,
		double(*)	matrix[P4EST_CHILDREN][P4EST_CHILDREN],
		double *	tmp,
		double *	lump
	)		`[static]`

Multiply the mass matrix with a vector of ones to compute the volume.

Parameters:

[in]	p4est	The forest is not changed.
[in]	lnodes	The node numbering is not changed.
[in]	matrix	The mass matrix should be passed in here.
[in]	tmp	Must be allocated, entries are undefined.
[in,out]	lump	Must be allocated, receives matrix * ones.

static void vector_axpy	(	p4est_t *	p4est,
		p4est_lnodes_t *	lnodes,
		double	a,
		const double *	x,
		double *	y
	)		`[static]`

Compute y := y + a * x.

Parameters:

[in]	p4est	The forest is not changed.
[in]	lnodes	The node numbering is not changed.
[in]	a	The scalar.
[in]	x	First node vector.
[in,out]	y	Second node vector.

static void vector_copy	(	p4est_t *	p4est,
		p4est_lnodes_t *	lnodes,
		const double *	x,
		double *	y
	)		`[static]`

copy a vector.

Parameters:

[in]	p4est	The forest is not changed.
[in]	lnodes	The node numbering is not changed.
[in]	x	Input node vector.
[out]	y	output node vector.

static double vector_dot	(	p4est_t *	p4est,
		p4est_lnodes_t *	lnodes,
		const double *	v1,
		const double *	v2
	)		`[static]`

Compute the inner product of two node vectors in parallel.

Parameters:

[in]	p4est	The forest is not changed.
[in]	lnodes	The node numbering is not changed.
[in]	v1	First node vector.
[in]	v2	Second node vector.

Returns:: Parallel l_2 inner product over the domain.

static void vector_xpby	(	p4est_t *	p4est,
		p4est_lnodes_t *	lnodes,
		const double *	x,
		double	b,
		double *	y
	)		`[static]`

Compute y := x + b * y.

Parameters:

[in]	p4est	The forest is not changed.
[in]	lnodes	The node numbering is not changed.
[in]	x	First node vector.
[in]	b	The scalar.
[in,out]	y	Second node vector.

static void vector_zero	(	p4est_t *	p4est,
		p4est_lnodes_t *	lnodes,
		double *	x
	)		`[static]`

Zero all entries of a vector.

Parameters:

[in]	p4est	The forest is not changed.
[in]	lnodes	The node numbering is not changed.
[out]	x	The vector.

Functions

Detailed Description

Function Documentation