p4est  1.1
p4est_step4.c File Reference

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>`

## Functions

Callback function to decide on refinement. More...

static double func_rhs (const double vxyz)
Right hand side function for the 2D Poisson problem. More...

static double func_uexact (const double vxyz)
Exact solution for the 2D Poisson problem. More...

static int is_boundary_unitsquare (p4est_t *p4est, p4est_topidx_t tt, p4est_quadrant_t *node)
Determine the boundary status on the unit square/cube. More...

static int lnodes_decode2 (p4est_lnodes_code_t face_code, int hanging_corner[P4EST_CHILDREN])
Decode the information from p{4,8}est_lnodes_t for a given element. More...

static void interpolate_hanging_nodes (p4est_lnodes_code_t face_code, double inplace[P4EST_CHILDREN])
Compute values at hanging nodes by interpolation. More...

static void share_sum (p4est_t *p4est, p4est_lnodes_t *lnodes, double *v)
Parallel sum of values in node vector across all sharers. More...

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. More...

static void vector_axpy (p4est_t *p4est, p4est_lnodes_t *lnodes, double a, const double *x, double *y)
Compute y := y + a * x. More...

static void vector_xpby (p4est_t *p4est, p4est_lnodes_t *lnodes, const double *x, double b, double *y)
Compute y := x + b * y. More...

static void vector_zero (p4est_t *p4est, p4est_lnodes_t *lnodes, double *x)
Zero all entries of a vector. More...

static void vector_copy (p4est_t *p4est, p4est_lnodes_t *lnodes, const double *x, double *y)
copy a vector. More...

static double * allocate_vector (p4est_lnodes_t *lnodes)
Allocate storage for processor-relevant nodal degrees of freedom. More...

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. More...

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. More...

static void set_dirichlet (p4est_lnodes_t *lnodes, const int8_t *bc, double *v)
Set Dirichlet boundary values of a node vector to zero. More...

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. More...

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. More...

static void solve_poisson (p4est_t *p4est)
Execute the numerical part of the example: Solve Poisson's equation. More...

int main (int argc, char **argv)
The main function of the step4 example program. More...

## Variables

static const int corner_num_hanging [P4EST_CHILDREN]
List number of possible independent nodes for each hanging node. More...

static const int zero = 0
Constant zero. More...

static const int ones = P4EST_CHILDREN - 1
One bit per dimension. More...

static const int * corner_to_hanging [P4EST_CHILDREN]
For each node i of the reference quadrant, corner_num_hanging[i] many. More...

## 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 )
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. Even for the 2D case z is well defined, but ignored here.
Returns
Scalar function value at vxyz.
 static double func_uexact ( const double vxyz )
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. Even for the 2D case z is well defined, but ignored here.
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 void interpolate_hanging_nodes ( p4est_lnodes_code_t face_code, double inplace[P4EST_CHILDREN] )
static

Compute values at hanging nodes by interpolation.

A face hanging node in 3D depends on the four corner nodes of the face, edge hanging nodes or face hanging nodes in 2D depend on two nodes. This function works in place, we have to be careful about the ordering. Face hanging node values are not reused, so they are overwritten first.

Parameters
 [in] face_code This number encodes the child id of the quadrant and the hanging status of faces and edges. [in,out] inplace On input, the values at the independent nodes. On output, interpolated to hanging node locations.
 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/cube.

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 p{4,8}est_lnodes_t for a given element.

p4est_lnodes.h for an in-depth discussion of the encoding.
Parameters
 [in] face_code Bit code as defined in p{4,8}est_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/edge otherwise. For faces in 3D, it is diagonally opposite.
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

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.

## Variable Documentation

 const int corner_num_hanging[P4EST_CHILDREN]
static
Initial value:
=
{ 1, 2, 2, 1 }

List number of possible independent nodes for each hanging node.

 const int* corner_to_hanging[P4EST_CHILDREN]
static

For each node i of the reference quadrant, corner_num_hanging[i] many.

 const int ones = P4EST_CHILDREN - 1
static

One bit per dimension.

 const int zero = 0
static

Constant zero.