Base functions

Calculation of base functions at integration points. More...

Collaboration diagram for Base functions:

Functions
PetscErrorCode	Legendre_polynomials (int p, double s, double *diff_s, double *L, double *diffL, const int dim)
	Calculate Legendre approximation basis. More...
PetscErrorCode	Lobatto_polynomials (int p, double s, double *diff_s, double *L, double *diffL, const int dim)
	Calculate Lobatto base functions [29]. More...
PetscErrorCode	LobattoKernel_polynomials (int p, double s, double *diff_s, double *L, double *diffL, const int dim)
	Calculate Kernel Lobatto base functions. More...

Detailed Description

Calculation of base functions at integration points.

Function Documentation

Legendre_polynomials()

PetscErrorCode Legendre_polynomials	(	int	p,
		double	s,
		double *	diff_s,
		double *	L,
		double *	diffL,
		const int	dim
	)

Calculate Legendre approximation basis.

Lagrange polynomial is given by

\[ L_0(s)=1;\quad L_1(s) = s \]

and following terms are generated inductively

\[ L_{l+1}=\frac{2l+1}{l+1}sL_l(s)-\frac{l}{l+1}L_{l-1}(s) \]

Note that:

\[ s\in[-1,1] \quad \textrm{and}\; s=s(\xi_0,\xi_1,\xi_2) \]

where \(\xi_i\) are barycentric coordinates of element.

Parameters

p	is approximation order
s	is position \(s\in[-1,1]\)
diff_s	derivatives of shape functions, i.e. \(\frac{\partial s}{\partial \xi_i}\)

Return values

L	approximation functions
diffL	derivatives, i.e. \(\frac{\partial L}{\partial \xi_i}\)

Parameters

dim

dimension

Returns

error code

Examples

edge_and_bubble_shape_functions_on_quad.cpp.

Definition at line 15 of file base_functions.c.

                                                                  {
  MoFEMFunctionBeginHot;
#ifndef NDEBUG
  if (dim < 1)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
  if (dim > 3)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
  if (p < 0)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 0");
#endif // NDEBUG
 
  L[0] = 1;
  if (diffL != NULL) {
    for (int d = 0; d != dim; ++d) {
      diffL[d * (p + 1) + 0] = 0;
    }
  }
  if (p == 0)
    MoFEMFunctionReturnHot(0);
 
  L[1] = s;
  if (diffL != NULL) {
#ifndef NDEBUG
    if (diff_s == NULL) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
    }
#endif // NDEBUG
    for (int d = 0; d != dim; ++d) {
      diffL[d * (p + 1) + 1] = diff_s[d];
    }
  }
  if (p == 1)
    MoFEMFunctionReturnHot(0);
 
  int l = 1;
  for (; l < p; l++) {
    double A = ((2 * (double)l + 1) / ((double)l + 1));
    double B = ((double)l / ((double)l + 1));
    L[l + 1] = A * s * L[l] - B * L[l - 1];
    if (diffL != NULL) {
      for (int d = 0; d != dim; ++d) {
        diffL[d * (p + 1) + l + 1] =
            A * (diff_s[d] * L[l] + s * diffL[d * (p + 1) + l]) -
            B * diffL[d * (p + 1) + l - 1];
      }      
    }
  }
 
  MoFEMFunctionReturnHot(0);
}

Lobatto_polynomials()

PetscErrorCode Lobatto_polynomials	(	int	p,
		double	s,
		double *	diff_s,
		double *	L,
		double *	diffL,
		const int	dim
	)

Calculate Lobatto base functions [29].

Order of first function is 2 and goes to p.

Parameters

p	is approximation order
s	is a mapping of coordinates of edge to \([-1, 1]\), i.e., \(s(\xi_1,\cdot,\xi_{dim})\in[-1,1]\)
diff_s	jacobian of the transformation, i.e. \(\frac{\partial s}{\partial \xi_i}\) output

Return values

L	values basis functions at s
diffL	derivatives of basis functions at s, i.e. \(\frac{\partial L}{\partial \xi_i}\)

Parameters

dim

dimension

Returns

error code

Definition at line 199 of file base_functions.c.

                                                                 {
 
  MoFEMFunctionBeginHot;
#ifndef NDEBUG
  if (dim < 1)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
  if (dim > 3)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
  if (p < 2)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 2");
#endif // NDEBUG
  double l[p + 1];
  ierr = Legendre_polynomials(p, s, NULL, l, NULL, 1);
  CHKERRQ(ierr);
 
  L[0] = 1;
  if (diffL != NULL) {
    for (int d = 0; d != dim; ++d) {
      diffL[d * (p + 1) + 0] = 0;
    }
  }
  L[1] = s;
  if (diffL != NULL) {
#ifndef NDEBUG
    if (diff_s == NULL) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
    }
#endif // NDEBUG
    for (int d = 0; d != dim; ++d) {
      diffL[d * (p + 1) + 1] = diff_s[d];
    }
  }
 
  // Integrated Legendre
  for (int k = 2; k <= p; k++) {
    const double factor = 2 * (2 * k - 1);
    L[k] = 1.0 / factor * (l[k] - l[k - 2]);
  }
 
  if (diffL != NULL) {
    for (int k = 2; k <= p; k++) {
      double a = l[k - 1] / 2.;
      for (int d = 0; d != dim; ++d) {
        diffL[d * (p + 1) + k] = a * diff_s[d];
      }
    }
  }
  MoFEMFunctionReturnHot(0);
}

LobattoKernel_polynomials()

PetscErrorCode LobattoKernel_polynomials	(	int	p,
		double	s,
		double *	diff_s,
		double *	L,
		double *	diffL,
		const int	dim
	)

Calculate Kernel Lobatto base functions.

This is implemented using definitions from Hermes2d https://github.com/hpfem/hermes following book by Pavel Solin et al [solin2003higher].

Parameters

p	is approximation order
s	is position \(s\in[-1,1]\)
diff_s	derivatives of shape functions, i.e. \(\frac{\partial s}{\partial \xi_i}\)

Return values

L	approximation functions
diffL	derivatives, i.e. \(\frac{\partial L}{\partial \xi_i}\)

Parameters

dim

dimension

Returns

error code

Definition at line 345 of file base_functions.c.

                                                        {
  MoFEMFunctionBeginHot;
#ifndef NDEBUG
  if (dim < 1)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
  if (dim > 3)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
  if (p < 0)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 0");
  if (p > 9)
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Polynomial beyond order 9 is not implemented");
#endif // NDEBUG
  if (L) {
    int l = 0;
    for (; l != p + 1; l++) {
      L[l] = f_phi[l](s);
    }
  }
  if (diffL != NULL) {
#ifndef NDEBUG
    if (diff_s == NULL) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
    }
#endif // NDEBUG
    int l = 0;
    for (; l != p + 1; l++) {
      double a = f_phix[l](s);
      diffL[0 * (p + 1) + l] = diff_s[0] * a;
      if (dim >= 2) {
        diffL[1 * (p + 1) + l] = diff_s[1] * a;
        if (dim == 3) {
          diffL[2 * (p + 1) + l] = diff_s[2] * a;
        }
      }
    }
  }
  MoFEMFunctionReturnHot(0);
}