1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2023 by the deal.II authors 4 // 5 // This file is part of the deal.II library. 6 // 7 // The deal.II library is free software; you can use it, redistribute 8 // it, and/or modify it under the terms of the GNU Lesser General 9 // Public License as published by the Free Software Foundation; either 10 // version 2.1 of the License, or (at your option) any later version. 11 // The full text of the license can be found in the file LICENSE.md at 12 // the top level directory of deal.II. 13 // 14 // Authors: Peter Munch, Martin Kronbichler 15 // 16 // --------------------------------------------------------------------- 17 18 #include <ceed.h> 19 20 21 22 /** 23 * Context passed to libCEED Q-function. 24 */ 25 struct BuildContext 26 { 27 CeedInt dim, space_dim; 28 }; 29 30 31 32 /** 33 * libCEED Q-function for building quadrature data for a mass operator 34 */ 35 CEED_QFUNCTION(f_build_mass) 36 (void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) 37 { 38 BuildContext *bc = (BuildContext *)ctx; 39 const CeedScalar *J = in[0], *w = in[1]; 40 CeedScalar *qdata = out[0]; 41 42 switch (bc->dim + 10 * bc->space_dim) 43 { 44 case 11: 45 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i) 46 { 47 qdata[i] = J[i] * w[i]; 48 } 49 break; 50 case 22: 51 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i) 52 { 53 qdata[i] = (J[i + Q * 0] * J[i + Q * 3] - J[i + Q * 1] * J[i + Q * 2]) * w[i]; 54 } 55 break; 56 case 33: 57 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i) 58 { 59 qdata[i] = (J[i + Q * 0] * (J[i + Q * 4] * J[i + Q * 8] - J[i + Q * 5] * J[i + Q * 7]) - 60 J[i + Q * 1] * (J[i + Q * 3] * J[i + Q * 8] - J[i + Q * 5] * J[i + Q * 6]) + 61 J[i + Q * 2] * (J[i + Q * 3] * J[i + Q * 7] - J[i + Q * 4] * J[i + Q * 6])) * 62 w[i]; 63 } 64 break; 65 } 66 return 0; 67 } 68 69 70 71 /** 72 * libCEED Q-function for applying a mass operator 73 */ 74 CEED_QFUNCTION(f_apply_mass) 75 (void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) 76 { 77 (void)ctx; 78 79 const CeedScalar *u = in[0]; 80 const CeedScalar *JxW = in[1]; 81 CeedScalar *v = out[0]; 82 83 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i) 84 { 85 v[i] = JxW[i] * u[i]; 86 } 87 return 0; 88 } 89 90 91 92 /** 93 * libCEED Q-function for applying a vector mass operator 94 */ 95 CEED_QFUNCTION(f_apply_mass_vec) 96 (void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) 97 { 98 BuildContext *bc = (BuildContext *)ctx; 99 100 const CeedScalar *u = in[0]; 101 const CeedScalar *JxW = in[1]; 102 CeedScalar *v = out[0]; 103 104 switch (bc->dim + 10 * bc->space_dim) 105 { 106 case 11: 107 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i) 108 { 109 v[i + Q * 0] = JxW[i] * u[i + Q * 0]; 110 } 111 break; 112 case 22: 113 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i) 114 { 115 v[i + Q * 0] = JxW[i] * u[i + Q * 0]; 116 v[i + Q * 1] = JxW[i] * u[i + Q * 1]; 117 } 118 break; 119 case 33: 120 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i) 121 { 122 v[i + Q * 0] = JxW[i] * u[i + Q * 0]; 123 v[i + Q * 1] = JxW[i] * u[i + Q * 1]; 124 v[i + Q * 2] = JxW[i] * u[i + Q * 2]; 125 } 126 break; 127 } 128 return 0; 129 } 130 131 132 133 /** 134 * libCEED Q-function for building quadrature data for a Poisson operator 135 */ 136 CEED_QFUNCTION(f_build_poisson) 137 (void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) 138 { 139 BuildContext *bc = (BuildContext *)ctx; 140 const CeedScalar *J = in[0], *w = in[1]; 141 CeedScalar *qdata = out[0]; 142 143 switch (bc->dim + 10 * bc->space_dim) 144 { 145 case 11: 146 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) 147 { 148 qdata[i] = w[i] / J[i]; 149 } 150 break; 151 case 22: 152 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) 153 { 154 const CeedScalar J11 = J[i + Q * 0]; 155 const CeedScalar J21 = J[i + Q * 1]; 156 const CeedScalar J12 = J[i + Q * 2]; 157 const CeedScalar J22 = J[i + Q * 3]; 158 const CeedScalar qw = w[i] / (J11 * J22 - J21 * J12); 159 qdata[i + Q * 0] = qw * (J12 * J12 + J22 * J22); 160 qdata[i + Q * 1] = qw * (J11 * J11 + J21 * J21); 161 qdata[i + Q * 2] = -qw * (J11 * J12 + J21 * J22); 162 } 163 break; 164 case 33: 165 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) 166 { 167 const CeedScalar J11 = J[i + Q * 0]; 168 const CeedScalar J21 = J[i + Q * 1]; 169 const CeedScalar J31 = J[i + Q * 2]; 170 const CeedScalar J12 = J[i + Q * 3]; 171 const CeedScalar J22 = J[i + Q * 4]; 172 const CeedScalar J32 = J[i + Q * 5]; 173 const CeedScalar J13 = J[i + Q * 6]; 174 const CeedScalar J23 = J[i + Q * 7]; 175 const CeedScalar J33 = J[i + Q * 8]; 176 const CeedScalar A11 = J22 * J33 - J23 * J32; 177 const CeedScalar A12 = J13 * J32 - J12 * J33; 178 const CeedScalar A13 = J12 * J23 - J13 * J22; 179 const CeedScalar A21 = J23 * J31 - J21 * J33; 180 const CeedScalar A22 = J11 * J33 - J13 * J31; 181 const CeedScalar A23 = J13 * J21 - J11 * J23; 182 const CeedScalar A31 = J21 * J32 - J22 * J31; 183 const CeedScalar A32 = J12 * J31 - J11 * J32; 184 const CeedScalar A33 = J11 * J22 - J12 * J21; 185 const CeedScalar qw = w[i] / (J11 * A11 + J21 * A12 + J31 * A13); 186 qdata[i + Q * 0] = qw * (A11 * A11 + A12 * A12 + A13 * A13); 187 qdata[i + Q * 1] = qw * (A21 * A21 + A22 * A22 + A23 * A23); 188 qdata[i + Q * 2] = qw * (A31 * A31 + A32 * A32 + A33 * A33); 189 qdata[i + Q * 3] = qw * (A21 * A31 + A22 * A32 + A23 * A33); 190 qdata[i + Q * 4] = qw * (A11 * A31 + A12 * A32 + A13 * A33); 191 qdata[i + Q * 5] = qw * (A11 * A21 + A12 * A22 + A13 * A23); 192 } 193 break; 194 } 195 return 0; 196 } 197 198 199 200 /** 201 * libCEED Q-function for applying a Poisson operator 202 */ 203 CEED_QFUNCTION(f_apply_poisson) 204 (void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) 205 { 206 BuildContext *bc = (BuildContext *)ctx; 207 const CeedScalar *ug = in[0], *qdata = in[1]; 208 CeedScalar *vg = out[0]; 209 210 switch (bc->dim) 211 { 212 case 1: 213 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) 214 { 215 vg[i] = ug[i] * qdata[i]; 216 } 217 break; 218 case 2: 219 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) 220 { 221 const CeedScalar ug0 = ug[i + Q * 0]; 222 const CeedScalar ug1 = ug[i + Q * 1]; 223 vg[i + Q * 0] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 2] * ug1; 224 vg[i + Q * 1] = qdata[i + Q * 2] * ug0 + qdata[i + Q * 1] * ug1; 225 } 226 break; 227 case 3: 228 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) 229 { 230 const CeedScalar ug0 = ug[i + Q * 0]; 231 const CeedScalar ug1 = ug[i + Q * 1]; 232 const CeedScalar ug2 = ug[i + Q * 2]; 233 vg[i + Q * 0] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 5] * ug1 + qdata[i + Q * 4] * ug2; 234 vg[i + Q * 1] = qdata[i + Q * 5] * ug0 + qdata[i + Q * 1] * ug1 + qdata[i + Q * 3] * ug2; 235 vg[i + Q * 2] = qdata[i + Q * 4] * ug0 + qdata[i + Q * 3] * ug1 + qdata[i + Q * 2] * ug2; 236 } 237 break; 238 } 239 return 0; 240 } 241 242 243 244 /** 245 *libCEED Q-function for applying a vector Poisson operator 246 */ 247 CEED_QFUNCTION(f_apply_poisson_vec) 248 (void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) 249 { 250 BuildContext *bc = (BuildContext *)ctx; 251 const CeedScalar *ug = in[0], *qdata = in[1]; 252 CeedScalar *vg = out[0]; 253 254 switch (bc->dim) 255 { 256 case 1: 257 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) 258 { 259 vg[i] = ug[i] * qdata[i]; 260 } 261 break; 262 case 2: 263 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) 264 { 265 { 266 const CeedScalar ug0 = ug[i + Q * 0 + Q * 2 * 0]; 267 const CeedScalar ug1 = ug[i + Q * 1 + Q * 2 * 0]; 268 vg[i + Q * 0 + Q * 2 * 0] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 2] * ug1; 269 vg[i + Q * 1 + Q * 2 * 0] = qdata[i + Q * 2] * ug0 + qdata[i + Q * 1] * ug1; 270 } 271 { 272 const CeedScalar ug0 = ug[i + Q * 0 + Q * 2 * 1]; 273 const CeedScalar ug1 = ug[i + Q * 1 + Q * 2 * 1]; 274 vg[i + Q * 0 + Q * 2 * 1] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 2] * ug1; 275 vg[i + Q * 1 + Q * 2 * 1] = qdata[i + Q * 2] * ug0 + qdata[i + Q * 1] * ug1; 276 } 277 } 278 break; 279 case 3: 280 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) 281 { 282 { 283 const CeedScalar ug0 = ug[i + Q * 0 + Q * 3 * 0]; 284 const CeedScalar ug1 = ug[i + Q * 1 + Q * 3 * 0]; 285 const CeedScalar ug2 = ug[i + Q * 2 + Q * 3 * 0]; 286 vg[i + Q * 0 + Q * 3 * 0] = 287 qdata[i + Q * 0] * ug0 + qdata[i + Q * 5] * ug1 + qdata[i + Q * 4] * ug2; 288 vg[i + Q * 1 + Q * 3 * 0] = 289 qdata[i + Q * 5] * ug0 + qdata[i + Q * 1] * ug1 + qdata[i + Q * 3] * ug2; 290 vg[i + Q * 2 + Q * 3 * 0] = 291 qdata[i + Q * 4] * ug0 + qdata[i + Q * 3] * ug1 + qdata[i + Q * 2] * ug2; 292 } 293 { 294 const CeedScalar ug0 = ug[i + Q * 0 + Q * 3 * 1]; 295 const CeedScalar ug1 = ug[i + Q * 1 + Q * 3 * 1]; 296 const CeedScalar ug2 = ug[i + Q * 2 + Q * 3 * 1]; 297 vg[i + Q * 0 + Q * 3 * 1] = 298 qdata[i + Q * 0] * ug0 + qdata[i + Q * 5] * ug1 + qdata[i + Q * 4] * ug2; 299 vg[i + Q * 1 + Q * 3 * 1] = 300 qdata[i + Q * 5] * ug0 + qdata[i + Q * 1] * ug1 + qdata[i + Q * 3] * ug2; 301 vg[i + Q * 2 + Q * 3 * 1] = 302 qdata[i + Q * 4] * ug0 + qdata[i + Q * 3] * ug1 + qdata[i + Q * 2] * ug2; 303 } 304 { 305 const CeedScalar ug0 = ug[i + Q * 0 + Q * 3 * 2]; 306 const CeedScalar ug1 = ug[i + Q * 1 + Q * 3 * 2]; 307 const CeedScalar ug2 = ug[i + Q * 2 + Q * 3 * 2]; 308 vg[i + Q * 0 + Q * 3 * 2] = 309 qdata[i + Q * 0] * ug0 + qdata[i + Q * 5] * ug1 + qdata[i + Q * 4] * ug2; 310 vg[i + Q * 1 + Q * 3 * 2] = 311 qdata[i + Q * 5] * ug0 + qdata[i + Q * 1] * ug1 + qdata[i + Q * 3] * ug2; 312 vg[i + Q * 2 + Q * 3 * 2] = 313 qdata[i + Q * 4] * ug0 + qdata[i + Q * 3] * ug1 + qdata[i + Q * 2] * ug2; 314 } 315 } 316 break; 317 } 318 return 0; 319 } 320