13d8e8822SJeremy L Thompson // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 23d8e8822SJeremy L Thompson // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 350c301a5SRezgar Shakeri // 43d8e8822SJeremy L Thompson // SPDX-License-Identifier: BSD-2-Clause 550c301a5SRezgar Shakeri // 63d8e8822SJeremy L Thompson // This file is part of CEED: http://github.com/ceed 750c301a5SRezgar Shakeri 82b730f8bSJeremy L Thompson #include <ceed.h> 92b730f8bSJeremy L Thompson 10*b5404d5dSSebastian Grimberg // H(div) basis for quadrilateral linear BDM element in 2D 11*b5404d5dSSebastian Grimberg // Local numbering is as follows (each edge has 2 vector dof) 1250c301a5SRezgar Shakeri // b4 b5 1350c301a5SRezgar Shakeri // 2---------3 1450c301a5SRezgar Shakeri // b7| |b3 1550c301a5SRezgar Shakeri // | | 1650c301a5SRezgar Shakeri // b6| |b2 1750c301a5SRezgar Shakeri // 0---------1 1850c301a5SRezgar Shakeri // b0 b1 1950c301a5SRezgar Shakeri // Bx[0-->7] = b0_x-->b7_x, By[0-->7] = b0_y-->b7_y 2050c301a5SRezgar Shakeri // To see how the nodal basis is constructed visit: 2150c301a5SRezgar Shakeri // https://github.com/rezgarshakeri/H-div-Tests 224fee36f0SJeremy L Thompson int BuildNodalHdivQuadrilateral(CeedScalar *x, CeedScalar *Bx, CeedScalar *By) { 234fee36f0SJeremy L Thompson CeedScalar x_hat = x[0]; 244fee36f0SJeremy L Thompson CeedScalar y_hat = x[1]; 2550c301a5SRezgar Shakeri Bx[0] = -0.125 + 0.125 * x_hat * x_hat; 2650c301a5SRezgar Shakeri By[0] = -0.25 + 0.25 * x_hat + 0.25 * y_hat + -0.25 * x_hat * y_hat; 2750c301a5SRezgar Shakeri Bx[1] = 0.125 + -0.125 * x_hat * x_hat; 2850c301a5SRezgar Shakeri By[1] = -0.25 + -0.25 * x_hat + 0.25 * y_hat + 0.25 * x_hat * y_hat; 2950c301a5SRezgar Shakeri Bx[2] = 0.25 + 0.25 * x_hat + -0.25 * y_hat + -0.25 * x_hat * y_hat; 3050c301a5SRezgar Shakeri By[2] = -0.125 + 0.125 * y_hat * y_hat; 3150c301a5SRezgar Shakeri Bx[3] = 0.25 + 0.25 * x_hat + 0.25 * y_hat + 0.25 * x_hat * y_hat; 3250c301a5SRezgar Shakeri By[3] = 0.125 + -0.125 * y_hat * y_hat; 3350c301a5SRezgar Shakeri Bx[4] = -0.125 + 0.125 * x_hat * x_hat; 3450c301a5SRezgar Shakeri By[4] = 0.25 + -0.25 * x_hat + 0.25 * y_hat + -0.25 * x_hat * y_hat; 3550c301a5SRezgar Shakeri Bx[5] = 0.125 + -0.125 * x_hat * x_hat; 3650c301a5SRezgar Shakeri By[5] = 0.25 + 0.25 * x_hat + 0.25 * y_hat + 0.25 * x_hat * y_hat; 3750c301a5SRezgar Shakeri Bx[6] = -0.25 + 0.25 * x_hat + 0.25 * y_hat + -0.25 * x_hat * y_hat; 3850c301a5SRezgar Shakeri By[6] = -0.125 + 0.125 * y_hat * y_hat; 3950c301a5SRezgar Shakeri Bx[7] = -0.25 + 0.25 * x_hat + -0.25 * y_hat + 0.25 * x_hat * y_hat; 4050c301a5SRezgar Shakeri By[7] = 0.125 + -0.125 * y_hat * y_hat; 4150c301a5SRezgar Shakeri return 0; 4250c301a5SRezgar Shakeri } 43*b5404d5dSSebastian Grimberg 444fee36f0SJeremy L Thompson static void BuildHdivQuadrilateral(CeedInt q, CeedScalar *q_ref, CeedScalar *q_weights, CeedScalar *interp, CeedScalar *div, CeedQuadMode quad_mode) { 4550c301a5SRezgar Shakeri // Get 1D quadrature on [-1,1] 464fee36f0SJeremy L Thompson CeedScalar q_ref_1d[q], q_weight_1d[q]; 4750c301a5SRezgar Shakeri switch (quad_mode) { 4850c301a5SRezgar Shakeri case CEED_GAUSS: 494fee36f0SJeremy L Thompson CeedGaussQuadrature(q, q_ref_1d, q_weight_1d); 5050c301a5SRezgar Shakeri break; 5159058f14SJeremy L Thompson // LCOV_EXCL_START 5250c301a5SRezgar Shakeri case CEED_GAUSS_LOBATTO: 534fee36f0SJeremy L Thompson CeedLobattoQuadrature(q, q_ref_1d, q_weight_1d); 5450c301a5SRezgar Shakeri break; 5550c301a5SRezgar Shakeri } 5659058f14SJeremy L Thompson // LCOV_EXCL_STOP 5750c301a5SRezgar Shakeri 5850c301a5SRezgar Shakeri // Divergence operator; Divergence of nodal basis for ref element 5950c301a5SRezgar Shakeri CeedScalar D[8] = {0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}; 6050c301a5SRezgar Shakeri CeedScalar Bx[8], By[8]; 6150c301a5SRezgar Shakeri CeedScalar X[2]; 6250c301a5SRezgar Shakeri 63*b5404d5dSSebastian Grimberg // Loop over quadrature points 644fee36f0SJeremy L Thompson for (CeedInt i = 0; i < q; i++) { 654fee36f0SJeremy L Thompson for (CeedInt j = 0; j < q; j++) { 664fee36f0SJeremy L Thompson CeedInt k1 = q * i + j; 6750c301a5SRezgar Shakeri q_ref[k1] = q_ref_1d[j]; 684fee36f0SJeremy L Thompson q_ref[k1 + q * q] = q_ref_1d[i]; 6950c301a5SRezgar Shakeri q_weights[k1] = q_weight_1d[j] * q_weight_1d[i]; 7050c301a5SRezgar Shakeri X[0] = q_ref_1d[j]; 7150c301a5SRezgar Shakeri X[1] = q_ref_1d[i]; 724fee36f0SJeremy L Thompson BuildNodalHdivQuadrilateral(X, Bx, By); 7350c301a5SRezgar Shakeri for (CeedInt k = 0; k < 8; k++) { 7450c301a5SRezgar Shakeri interp[k1 * 8 + k] = Bx[k]; 754fee36f0SJeremy L Thompson interp[k1 * 8 + k + 8 * q * q] = By[k]; 7650c301a5SRezgar Shakeri div[k1 * 8 + k] = D[k]; 7750c301a5SRezgar Shakeri } 7850c301a5SRezgar Shakeri } 7950c301a5SRezgar Shakeri } 8050c301a5SRezgar Shakeri } 81