Skip to content

Commit

Permalink
SSAFEM CFBC test passes with the new fem::BoundaryQuadrature
Browse files Browse the repository at this point in the history 
Next: implement fem::p1::BoundaryQuadrature2 (2-point Gaussian
quadrature on sides of P1 elements).
  • Loading branch information
@ckhroulev
ckhroulev committed Oct 14, 2020
1 parent 0dd8413 commit e451324
Show file tree
Hide file tree
Showing 3 changed files with 153 additions and 86 deletions.
25 changes: 13 additions & 12 deletions src/stressbalance/ssa/SSAFEM.cc
View file
Expand Up @@ -547,12 +547,15 @@ void SSAFEM::cache_residual_cfbc(const Inputs &inputs) {

const double dx = m_grid->dx(), dy = m_grid->dy();

// Q1 boundary quadrature
fem::q1::BoundaryQuadrature2 bq(dx, dy, 1.0);

// Q1 element geometry information
fem::q1::Q1ElementGeometry q1;

// Q1 boundary quadratures
fem::q1::BoundaryQuadrature2 q1_bq(dx, dy, 1.0);

// P1 element geometry information (one per P1 element type)
std::vector<fem::p1::P1ElementGeometry> p1;

for (unsigned int k = 0; k < 4; ++k) {
p1.push_back(fem::p1::P1ElementGeometry(k, dx, dy));
}
Expand Down Expand Up @@ -622,9 +625,11 @@ void SSAFEM::cache_residual_cfbc(const Inputs &inputs) {
const bool p1_interior_node = n_exterior_nodes == 1;

fem::ElementGeometry *E = NULL;
fem::BoundaryQuadrature *Q = NULL;

if (q1_interior_element) {
E = &q1;
Q = &q1_bq;
} else if (p1_interior_node) {
int type = -1;

Expand All @@ -637,8 +642,6 @@ void SSAFEM::cache_residual_cfbc(const Inputs &inputs) {

E = &p1[type];
continue;

// quadrature = bq_p1[p1_number];
} else {
// an exterior node
continue;
Expand All @@ -660,7 +663,7 @@ void SSAFEM::cache_residual_cfbc(const Inputs &inputs) {
// second
double psi[2] = {0.0, 0.0};

const unsigned int Nq = bq.n();
const unsigned int Nq = Q->n();
const unsigned int n_sides = E->n_sides();
// loop over element sides
for (unsigned int side = 0; side < n_sides; ++side) {
Expand All @@ -676,16 +679,14 @@ void SSAFEM::cache_residual_cfbc(const Inputs &inputs) {
continue;
}

// in our case (i.e. uniform spacing in x and y directions) W is the same at all
// quadrature points along a side.
const double W = bq.weights(side);

for (unsigned int q = 0; q < Nq; ++q) {

const double W = Q->weight(side, q);

// test functions at nodes incident to the current side, evaluated at the quadrature point
// q
psi[0] = bq.germ(side, q, n0).val;
psi[1] = bq.germ(side, q, n1).val;
psi[0] = Q->germ(side, q, n0).val;
psi[1] = Q->germ(side, q, n1).val;

// Compute ice thickness and bed elevation at a quadrature point. This uses a 1D basis
// expansion on the current side.
Expand Down
155 changes: 113 additions & 42 deletions src/util/FETools.cc
View file
Expand Up @@ -66,6 +66,23 @@ static double determinant(const double J[2][2]) {
return J[0][0] * J[1][1] - J[1][0] * J[0][1];
}

// Multiply a matrix by a vector.
static Vector2 multiply(const double A[2][2], const Vector2 &v) {
Vector2 result;
result.u = v.u * A[0][0] + v.v * A[0][1];
result.v = v.u * A[1][0] + v.v * A[1][1];
return result;
}

//! Compute derivatives with respect to x,y using J^{-1} and derivatives with respect to xi, eta.
static Germ multiply(const double A[2][2], const Germ &v) {
Germ result;
result.val = v.val;
result.dx = v.dx * A[0][0] + v.dy * A[0][1];
result.dy = v.dx * A[1][0] + v.dy * A[1][1];
return result;
}

//! Compute the inverse of a two by two matrix.
static void invert(const double A[2][2], double A_inv[2][2]) {
const double det_A = determinant(A);
Expand All @@ -78,15 +95,6 @@ static void invert(const double A[2][2], double A_inv[2][2]) {
A_inv[1][1] = A[0][0] / det_A;
}

//! Compute derivatives with respect to x,y using J^{-1} and derivatives with respect to xi, eta.
static Germ apply_jacobian_inverse(const double J_inv[2][2], const Germ &f) {
Germ result;
result.val = f.val;
result.dx = f.dx * J_inv[0][0] + f.dy * J_inv[0][1];
result.dy = f.dx * J_inv[1][0] + f.dy * J_inv[1][1];
return result;
}

ElementGeometry::ElementGeometry(unsigned int n)
: m_n_sides(n) {
// empty
Expand All @@ -106,6 +114,37 @@ Vector2 ElementGeometry::normal(unsigned int side) const {
return m_normals[side];
}

BoundaryQuadrature::BoundaryQuadrature(unsigned int size)
: m_Nq(size) {
// empty
}

BoundaryQuadrature::~BoundaryQuadrature() {
// empty
}

unsigned int BoundaryQuadrature::n() const {
return m_Nq;
}

double BoundaryQuadrature::weight(unsigned int side,
unsigned int q) const {
assert(side < q1::n_sides);
assert(q < m_Nq);

return this->weight_impl(side, q);
}

const Germ& BoundaryQuadrature::germ(unsigned int side,
unsigned int q,
unsigned int test_function) const {
assert(side < q1::n_sides);
assert(q < m_Nq);
assert(test_function < q1::n_chi);

return this->germ_impl(side, q, test_function);
}

namespace q1 {

Q1ElementGeometry::Q1ElementGeometry()
Expand All @@ -124,13 +163,14 @@ unsigned int Q1ElementGeometry::incident_node_impl(unsigned int side, unsigned i
return nodes[side][k];
}

// coordinates of reference element nodes
static const double xi[n_chi] = {-1.0, 1.0, 1.0, -1.0};
static const double eta[n_chi] = {-1.0, -1.0, 1.0, 1.0};

//! Q1 basis functions on the reference element with nodes (-1,-1), (1,-1), (1,1), (-1,1).
Germ chi(unsigned int k, const QuadPoint &pt) {
assert(k < q1::n_chi);

const double xi[] = {-1.0, 1.0, 1.0, -1.0};
const double eta[] = {-1.0, -1.0, 1.0, 1.0};

Germ result;

result.val = 0.25 * (1.0 + xi[k] * pt.xi) * (1.0 + eta[k] * pt.eta);
Expand All @@ -140,45 +180,74 @@ Germ chi(unsigned int k, const QuadPoint &pt) {
return result;
}

BoundaryQuadrature2::BoundaryQuadrature2(double dx, double dy, double L) {
// Parameterization of sides of the Q1 reference element (t \in [-1, 1]).
static QuadPoint r_star(unsigned int side, double t) {

// Map t (in [-1, 1]) to [0, 1] to simplify interpolation
const double L = 0.5 * (t + 1.0);

const unsigned int
k = n_chi - 1, // largest allowed index
j0 = side,
j1 = side % k + 1;

QuadPoint result;
result.xi = (1.0 - L) * xi[j0] + L * xi[j1];
result.eta = (1.0 - L) * eta[j0] + L * eta[j1];

return result;
}

BoundaryQuadrature2::BoundaryQuadrature2(double dx, double dy, double L)
: BoundaryQuadrature(m_size) {

// The Jacobian of the map from the reference element to a physical element.
const double J[2][2] = {{0.5 * dx / L, 0.0},
{0.0, 0.5 * dy / L}};

// The inverse of the Jacobian.
// derivative of r_star(t) = (xi(t), eta(t)) (the parameterization of the selected side of the
// reference element) with respect to t. See fem_q1_boundary.mac for a derivation.
Vector2 dr_star[n_sides] = {{1.0, 0.0}, {0.0, 1.0}, {-1.0, 0.0}, {0.0, -1.0}};

// 2-point Gaussian quadrature on [-1, 1].
const double points[m_size] = {-1.0 / sqrt(3.0), 1.0 / sqrt(3.0)};
const double weights[m_size] = {1.0, 1.0};

// The inverse of the Jacobian
double J_inv[2][2];
invert(J, J_inv);

// Note that all quadrature weights are 1.0 (and so they are implicitly included below).
//
// bottom
m_W[0] = J[0][0];
// right
m_W[1] = J[1][1];
// top
m_W[2] = J[0][0];
// left
m_W[3] = J[1][1];

const double C = 1.0 / sqrt(3);
const QuadPoint points[q1::n_sides][m_Nq] = {
{{ -C, -1.0}, { C, -1.0}}, // South
{{ 1.0, -C}, { 1.0, C}}, // East
{{ -C, 1.0}, { C, 1.0}}, // North
{{-1.0, -C}, {-1.0, C}} // West
};
for (unsigned int side = 0; side < n_sides; ++side) {
// Magnitude of the derivative r(t) = (x(t), y(t)) (the parameterization of the current side of
// a physical element) with respect to t, computed using the chain rule.
const Vector2 dr = multiply(J, dr_star[side]);

memset(m_germs, 0, q1::n_sides*m_Nq*q1::n_chi*sizeof(Germ));
for (unsigned int q = 0; q < m_size; ++q) {
QuadPoint pt = r_star(side, points[q]);

for (unsigned int side = 0; side < q1::n_sides; ++side) {
for (unsigned int q = 0; q < m_Nq; ++q) {
for (unsigned int k = 0; k < q1::n_chi; ++k) {
Germ phi = q1::chi(k, points[side][q]);
m_W[side][q] = weights[q] * dr.magnitude();

m_germs[side][q][k] = apply_jacobian_inverse(J_inv, phi);
for (unsigned int k = 0; k < n_chi; ++k) {
// Compute the value of the current shape function and convert derivatives with respect to
// xi and eta into derivatives with respect to x and y:
m_germs[side][q][k] = multiply(J_inv, q1::chi(k, pt));
}
}
}
} // end of loop over quadrature points
} // end of loop over sides
}

double BoundaryQuadrature2::weight_impl(unsigned int side, unsigned int q) const {
return m_W[side][q];
}


//! @brief Return the "germ" (value and partial derivatives) of a basis function @f$ \chi_k @f$
//! evaluated at the point `pt` on the side `side` of an element.
const Germ& BoundaryQuadrature2::germ_impl(unsigned int side,
unsigned int q,
unsigned int k) const {

return m_germs[side][q][k];
}

} // end of namespace q1
Expand Down Expand Up @@ -226,6 +295,7 @@ P1ElementGeometry::P1ElementGeometry(unsigned int type, double dx, double dy)
m_normals[0] = n23;
m_normals[1] = n30;
m_normals[2] = -1.0 * n20;
break;
}
}

Expand Down Expand Up @@ -267,7 +337,8 @@ Germ chi(unsigned int k, const QuadPoint &pt) {
result.val = 0.0;
result.dx = 0.0;
result.dy = 0.0;
}
break;
}
return result;
}

Expand Down Expand Up @@ -594,7 +665,7 @@ void Quadrature::initialize(ShapeFunction2 f,

for (unsigned int q = 0; q < m_Nq; q++) {
for (unsigned int k = 0; k < n_chi; k++) {
m_germs[q][k] = apply_jacobian_inverse(J_inv, f(k, points[q]));
m_germs[q][k] = multiply(J_inv, f(k, points[q]));
}
}

Expand Down
59 changes: 27 additions & 32 deletions src/util/FETools.hh
View file
Expand Up @@ -207,6 +207,26 @@ private:
unsigned int m_n_sides;
};

class BoundaryQuadrature {
public:
BoundaryQuadrature(unsigned int size);
virtual ~BoundaryQuadrature();

unsigned int n() const;

double weight(unsigned int side, unsigned int q) const;

const Germ& germ(unsigned int side, unsigned int q, unsigned int test_function) const;
protected:
virtual double weight_impl(unsigned int side,
unsigned int q) const = 0;
virtual const Germ& germ_impl(unsigned int side,
unsigned int q,
unsigned int test_function) const = 0;
private:
unsigned int m_Nq;
};

//! Q1 element information.
namespace q1 {

Expand All @@ -228,45 +248,20 @@ private:


//! 2-point quadrature for sides of Q1 quadrilateral elements.
class BoundaryQuadrature2 {
class BoundaryQuadrature2 : public BoundaryQuadrature {
public:
BoundaryQuadrature2(double dx, double dy, double L);

inline unsigned int n() const;

inline double weights(unsigned int side) const;

inline const Germ& germ(unsigned int side,
unsigned int func,
unsigned int pt) const;
protected:
double weight_impl(unsigned int side, unsigned int q) const;
const Germ& germ_impl(unsigned int side, unsigned int q, unsigned int test_function) const;
private:
//! Number of quadrature points per side.
static const unsigned int m_Nq = 2;
Germ m_germs[q1::n_sides][m_Nq][q1::n_chi];
double m_W[q1::n_sides];
static const unsigned int m_size = 2;
double m_W[n_sides][m_size];
Germ m_germs[n_sides][m_size][q1::n_chi];
};

inline unsigned int BoundaryQuadrature2::n() const {
return m_Nq;
}

inline double BoundaryQuadrature2::weights(unsigned int side) const {
assert(side < q1::n_sides);
return m_W[side];
}

//! @brief Return the "germ" (value and partial derivatives) of a basis function @f$ \chi_k @f$
//! evaluated at the point `pt` on the side `side` of an element.
inline const Germ& BoundaryQuadrature2::germ(unsigned int side,
unsigned int q,
unsigned int k) const {
assert(side < q1::n_sides);
assert(k < q1::n_chi);
assert(q < 2);

return m_germs[side][q][k];
}

} // end of namespace q1

//! @brief P1 element embedded in a Q1 element.
Expand Down

0 comments on commit e451324

Please sign in to comment.