Select Git revision
-
Christophe Geuzaine authoredChristophe Geuzaine authored
orthogonalBasis.cpp 5.69 KiB
// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
//
// Contributed by Amaury Johnen
#include <cmath>
#include <vector>
#include <algorithm>
#include "orthogonalBasis.h"
#include "FuncSpaceData.h"
#include "Numeric.h"
orthogonalBasis::orthogonalBasis(const FuncSpaceData &_data)
: _type(_data.elementType()), _order(_data.spaceOrder())
{
}
void orthogonalBasis::f(double u, double v, double w, double *sf) const
{
static double sf0[100];
static double sf1[100];
int k = 0;
switch(_type) {
case TYPE_LIN: LegendrePolynomials::f(_order, u, sf); return;
case TYPE_TRI:
if(u >= 1.) {
for(int k = 0; k <= _order; ++k) {
sf[k] = k + 1;
}
for(int k = _order; k < (_order + 1) * (_order + 2) / 2; ++k) {
sf[k] = 0;
}
return;
}
LegendrePolynomials::f(_order, 2 * v / (1 - u) - 1, sf0);
for(int i = 0; i <= _order; ++i) {
JacobiPolynomials::f(_order - i, 2 * i + 1, 0, 2 * u - 1, sf1);
double coeff = pow_int(1 - u, i);
for(int j = 0; j <= _order - i; ++j) {
sf1[j] *= coeff;
}
for(int j = 0; j <= _order - i; ++j) {
sf[k++] = sf0[i] * sf1[j];
}
}
return;
case TYPE_QUA:
LegendrePolynomials::f(_order, u, sf0);
LegendrePolynomials::f(_order, v, sf1);
for(int i = 0; i <= _order; ++i) {
for(int j = 0; j <= _order; ++j) {
sf[k++] = sf0[i] * sf1[j];
}
}
return;
}
}
void orthogonalBasis::integralfSquared(double *val) const
{
int k = 0;
switch(_type) {
case TYPE_LIN:
for(int i = 0; i <= _order; ++i) {
val[i] = 2. / (1 + 2 * i);
}
return;
case TYPE_TRI:
for(int i = 0; i <= _order; ++i) {
for(int j = 0; j <= _order - i; ++j) {
val[k++] = .5 / (1 + i + j) / (1 + 2 * i);
}
}
return;
case TYPE_QUA:
for(int i = 0; i <= _order; ++i) {
for(int j = 0; j <= _order; ++j) {
val[k++] = 4. / (1 + 2 * i) / (1 + 2 * j);
}
}
}
}
namespace LegendrePolynomials {
void f(int n, double u, double *val)
{
val[0] = 1;
for(int i = 0; i < n; i++) {
double a1i = i + 1;
double a3i = 2. * i + 1;
double a4i = i;
val[i + 1] = a3i * u * val[i];
if(i > 0) val[i + 1] -= a4i * val[i - 1];
val[i + 1] /= a1i;
}
}
void fc(int n, double u, double *val)
{
f(n, u, val);
for(int i = 2; i < n + 1; ++i) {
if(i % 2)
val[i] -= u;
else
val[i] -= 1;
}
}
void df(int n, double u, double *val)
{
// Indeterminate form for u == -1 and u == 1
if((u == 1.) || (u == -1.)) {
for(int k = 0; k <= n; k++) val[k] = 0.5 * k * (k + 1);
if((u == -1.) && (n >= 2))
for(int k = 2; k <= n; k += 2) val[k] = -val[k];
return;
}
// Now general case
// Values of the Legendre polynomials
std::vector<double> tmp(n + 1);
f(n, u, &(tmp[0]));
// First value of the derivative
val[0] = 0;
double g2 = (1. - u * u);
// Values of the derivative for orders > 1 computed from the values of the
// polynomials
for(int i = 1; i <= n; i++) {
double g1 = -u * i;
double g0 = (double)i;
val[i] = (g1 * tmp[i] + g0 * tmp[i - 1]) / g2;
}
}
} // namespace LegendrePolynomials
namespace JacobiPolynomials {
void f(int n, double alpha, double beta, double u, double *val)
{
const double alphaPlusBeta = alpha + beta;
const double a2MinusB2 = alpha * alpha - beta * beta;
val[0] = 1.;
if(n >= 1)
val[1] = 0.5 * (2. * (alpha + 1.) + (alphaPlusBeta + 2.) * (u - 1.));
for(int i = 1; i < n; i++) {
double ii = (double)i;
double twoI = 2. * ii;
double a1i =
2. * (ii + 1.) * (ii + alphaPlusBeta + 1.) * (twoI + alphaPlusBeta);
double a2i = (twoI + alphaPlusBeta + 1.) * (a2MinusB2);
double a3i = Pochhammer(twoI + alphaPlusBeta, 3);
double a4i =
2. * (ii + alpha) * (ii + beta) * (twoI + alphaPlusBeta + 2.);
val[i + 1] = ((a2i + a3i * u) * val[i] - a4i * val[i - 1]) / a1i;
}
}
void df(int n, double alpha, double beta, double u, double *val)
{
const double alphaPlusBeta = alpha + beta;
// Indeterminate form for u == -1 and u == 1
// TODO: Extend to non-integer alpha & beta?
if((u == 1.) || (u == -1.)) {
// alpha or beta in formula, depending on u
int coeff = (u == 1.) ? (int)alpha : (int)beta;
// Compute factorial
const int fMax = std::max(n, 1) + coeff;
std::vector<double> fact(fMax + 1);
fact[0] = 1.;
for(int i = 1; i <= fMax; i++) fact[i] = i * fact[i - 1];
// Compute formula (with appropriate sign at even orders for u == -1)
val[0] = 0.;
for(int k = 1; k <= n; k++)
val[k] = 0.5 * (k + alphaPlusBeta + 1) * fact[k + coeff] /
(fact[coeff + 1] * fact[k - 1]);
if((u == -1.) && (n >= 2))
for(int k = 2; k <= n; k += 2) val[k] = -val[k];
return;
}
// Now general case
// Values of the Jacobi polynomials
std::vector<double> tmp(n + 1);
f(n, alpha, beta, u, &(tmp[0]));
// First 2 values of the derivatives
val[0] = 0;
if(n >= 1) val[1] = 0.5 * (alphaPlusBeta + 2.);
// Values of the derivative for orders > 1 computed from the values of the
// polynomials
for(int i = 2; i <= n; i++) {
double ii = (double)i;
double aa = (2. * ii + alphaPlusBeta);
double g2 = aa * (1. - u * u);
double g1 = ii * (alpha - beta - aa * u);
double g0 = 2. * (ii + alpha) * (ii + beta);
val[i] = (g1 * tmp[i] + g0 * tmp[i - 1]) / g2;
}
}
} // namespace JacobiPolynomials