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 1050c301a5SRezgar Shakeri // Hdiv basis for quadrilateral linear BDM element in 2D 11*4fee36f0SJeremy L Thompson // Local numbering is as follow (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 22*4fee36f0SJeremy L Thompson int BuildNodalHdivQuadrilateral(CeedScalar *x, CeedScalar *Bx, CeedScalar *By) { 23*4fee36f0SJeremy L Thompson CeedScalar x_hat = x[0]; 24*4fee36f0SJeremy 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*4fee36f0SJeremy L Thompson static void BuildHdivQuadrilateral(CeedInt q, CeedScalar *q_ref, CeedScalar *q_weights, CeedScalar *interp, CeedScalar *div, CeedQuadMode quad_mode) { 4450c301a5SRezgar Shakeri // Get 1D quadrature on [-1,1] 45*4fee36f0SJeremy L Thompson CeedScalar q_ref_1d[q], q_weight_1d[q]; 4650c301a5SRezgar Shakeri switch (quad_mode) { 4750c301a5SRezgar Shakeri case CEED_GAUSS: 48*4fee36f0SJeremy L Thompson CeedGaussQuadrature(q, q_ref_1d, q_weight_1d); 4950c301a5SRezgar Shakeri break; 5059058f14SJeremy L Thompson // LCOV_EXCL_START 5150c301a5SRezgar Shakeri case CEED_GAUSS_LOBATTO: 52*4fee36f0SJeremy L Thompson CeedLobattoQuadrature(q, q_ref_1d, q_weight_1d); 5350c301a5SRezgar Shakeri break; 5450c301a5SRezgar Shakeri } 5559058f14SJeremy L Thompson // LCOV_EXCL_STOP 5650c301a5SRezgar Shakeri 5750c301a5SRezgar Shakeri // Divergence operator; Divergence of nodal basis for ref element 5850c301a5SRezgar Shakeri CeedScalar D[8] = {0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}; 5950c301a5SRezgar Shakeri // Loop over quadrature points 6050c301a5SRezgar Shakeri CeedScalar Bx[8], By[8]; 6150c301a5SRezgar Shakeri CeedScalar X[2]; 6250c301a5SRezgar Shakeri 63*4fee36f0SJeremy L Thompson for (CeedInt i = 0; i < q; i++) { 64*4fee36f0SJeremy L Thompson for (CeedInt j = 0; j < q; j++) { 65*4fee36f0SJeremy L Thompson CeedInt k1 = q * i + j; 6650c301a5SRezgar Shakeri q_ref[k1] = q_ref_1d[j]; 67*4fee36f0SJeremy L Thompson q_ref[k1 + q * q] = q_ref_1d[i]; 6850c301a5SRezgar Shakeri q_weights[k1] = q_weight_1d[j] * q_weight_1d[i]; 6950c301a5SRezgar Shakeri X[0] = q_ref_1d[j]; 7050c301a5SRezgar Shakeri X[1] = q_ref_1d[i]; 71*4fee36f0SJeremy L Thompson BuildNodalHdivQuadrilateral(X, Bx, By); 7250c301a5SRezgar Shakeri for (CeedInt k = 0; k < 8; k++) { 7350c301a5SRezgar Shakeri interp[k1 * 8 + k] = Bx[k]; 74*4fee36f0SJeremy L Thompson interp[k1 * 8 + k + 8 * q * q] = By[k]; 7550c301a5SRezgar Shakeri div[k1 * 8 + k] = D[k]; 7650c301a5SRezgar Shakeri } 7750c301a5SRezgar Shakeri } 7850c301a5SRezgar Shakeri } 7950c301a5SRezgar Shakeri }