p4est  1.1
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
Functions | Variables
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

static int refine_fn (p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
 Callback function to decide on refinement. More...
 
static double func_rhs (const double vxyz[3])
 Right hand side function for the 2D Poisson problem. More...
 
static double func_uexact (const double vxyz[3])
 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]lnodesThis 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]vxyzx, 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[3])
static

Exact solution for the 2D Poisson problem.

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

Parameters
[in]vxyzx, 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]p4estThe forest is not changed.
[in]lnodesThe node numbering is not changed.
[out]rhs_evalIs allocated and filled with function values.
[out]uexact_evalIs allocated and filled with function values.
[out]pbcBoolean 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_codeThis number encodes the child id of the quadrant and the hanging status of faces and edges.
[in,out]inplaceOn 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]p4estCan be used to access the connectivity.
[in]ttThe tree number (always zero for the unit square).
[in]nodeThe 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.

See Also
p4est_lnodes.h for an in-depth discussion of the encoding.
Parameters
[in]face_codeBit code as defined in p{4,8}est_lnodes.h.
[out]hanging_cornerUndefined 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]p4estThe forest is not changed.
[in]lnodesThe node numbering is not changed.
[in]bcBoolean flags for Dirichlet boundary nodes. If NULL, no special action is taken.
[in]stiffnessIf false use scaling for the mass matrix, if true use the scaling for stiffness matrix.
[in]matrixA 4x4 matrix computed on the reference element.
[in]inInput vector x.
[out]outOutput 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]lnodesThe node numbering is not changed.
[in]bcBoolean flags for Dirichlet boundary nodes. If NULL, this function does nothing.
[in,out]vDirichlet 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]p4estThe mesh is not changed.
[in]lnodesThe node numbering is not changed.
[in,out]vOn 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]p4estThe forest is not changed.
[in]lnodesThe node numbering is not changed.
[in]bcBoolean flags for Dirichlet boundary nodes.
[in]stiffnessIf false use scaling for the mass matrix, if true use the scaling for stiffness matrix.
[in]matrixThe mass matrix should be passed in here.
[in]bThe right hand side vector.
[out]xThe 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]p4estSolve 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]p4estThe forest is not changed.
[in]lnodesThe node numbering is not changed.
[in]matrixThe mass matrix should be passed in here.
[in]tmpMust be allocated, entries are undefined.
[in,out]lumpMust 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]p4estThe forest is not changed.
[in]lnodesThe node numbering is not changed.
[in]aThe scalar.
[in]xFirst node vector.
[in,out]ySecond node vector.
static void vector_copy ( p4est_t p4est,
p4est_lnodes_t lnodes,
const double *  x,
double *  y 
)
static

copy a vector.

Parameters
[in]p4estThe forest is not changed.
[in]lnodesThe node numbering is not changed.
[in]xInput node vector.
[out]youtput 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]p4estThe forest is not changed.
[in]lnodesThe node numbering is not changed.
[in]v1First node vector.
[in]v2Second 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]p4estThe forest is not changed.
[in]lnodesThe node numbering is not changed.
[in]xFirst node vector.
[in]bThe scalar.
[in,out]ySecond node vector.
static void vector_zero ( p4est_t p4est,
p4est_lnodes_t lnodes,
double *  x 
)
static

Zero all entries of a vector.

Parameters
[in]p4estThe forest is not changed.
[in]lnodesThe node numbering is not changed.
[out]xThe 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.