1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 #include <ceed.h> 9 #include <ceed/backend.h> 10 #include <math.h> 11 #include <stdbool.h> 12 #include <string.h> 13 14 #include "ceed-ref.h" 15 16 //------------------------------------------------------------------------------ 17 // Basis Apply 18 //------------------------------------------------------------------------------ 19 static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U, CeedVector V) { 20 Ceed ceed; 21 bool is_tensor_basis; 22 CeedInt dim, num_comp, q_comp, num_nodes, num_qpts; 23 const CeedInt add = (t_mode == CEED_TRANSPOSE); 24 const CeedScalar *u; 25 CeedScalar *v; 26 CeedTensorContract contract; 27 CeedBasis_Ref *impl; 28 29 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 30 CeedCallBackend(CeedBasisGetData(basis, &impl)); 31 CeedCallBackend(CeedBasisGetDimension(basis, &dim)); 32 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 33 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); 34 CeedCallBackend(CeedBasisGetNumNodes(basis, &num_nodes)); 35 CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); 36 CeedCallBackend(CeedBasisGetTensorContract(basis, &contract)); 37 if (U != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(U, CEED_MEM_HOST, &u)); 38 else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode"); 39 CeedCallBackend(CeedVectorGetArrayWrite(V, CEED_MEM_HOST, &v)); 40 41 // Clear v if operating in transpose 42 if (t_mode == CEED_TRANSPOSE) { 43 const CeedInt v_size = num_elem * num_comp * num_nodes; 44 45 for (CeedInt i = 0; i < v_size; i++) v[i] = (CeedScalar)0.0; 46 } 47 48 CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor_basis)); 49 if (is_tensor_basis) { 50 // Tensor basis 51 CeedInt P_1d, Q_1d; 52 53 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 54 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 55 switch (eval_mode) { 56 // Interpolate to/from quadrature points 57 case CEED_EVAL_INTERP: { 58 if (impl->has_collo_interp) { 59 memcpy(v, u, num_elem * num_comp * num_nodes * sizeof(u[0])); 60 } else { 61 CeedInt P = P_1d, Q = Q_1d; 62 63 if (t_mode == CEED_TRANSPOSE) { 64 P = Q_1d; 65 Q = P_1d; 66 } 67 CeedInt pre = num_comp * CeedIntPow(P, dim - 1), post = num_elem; 68 CeedScalar tmp[2][num_elem * num_comp * Q * CeedIntPow(P > Q ? P : Q, dim - 1)]; 69 const CeedScalar *interp_1d; 70 71 CeedCallBackend(CeedBasisGetInterp1D(basis, &interp_1d)); 72 for (CeedInt d = 0; d < dim; d++) { 73 CeedCallBackend(CeedTensorContractApply(contract, pre, P, post, Q, interp_1d, t_mode, add && (d == dim - 1), d == 0 ? u : tmp[d % 2], 74 d == dim - 1 ? v : tmp[(d + 1) % 2])); 75 pre /= P; 76 post *= Q; 77 } 78 } 79 } break; 80 // Evaluate the gradient to/from quadrature points 81 case CEED_EVAL_GRAD: { 82 // In CEED_NOTRANSPOSE mode: 83 // u has shape [dim, num_comp, P^dim, num_elem], row-major layout 84 // v has shape [dim, num_comp, Q^dim, num_elem], row-major layout 85 // In CEED_TRANSPOSE mode, the sizes of u and v are switched. 86 CeedInt P = P_1d, Q = Q_1d; 87 88 if (t_mode == CEED_TRANSPOSE) { 89 P = Q_1d; 90 Q = Q_1d; 91 } 92 CeedInt pre = num_comp * CeedIntPow(P, dim - 1), post = num_elem; 93 const CeedScalar *interp_1d; 94 95 CeedCallBackend(CeedBasisGetInterp1D(basis, &interp_1d)); 96 if (impl->collo_grad_1d) { 97 CeedScalar tmp[2][num_elem * num_comp * Q * CeedIntPow(P > Q ? P : Q, dim - 1)]; 98 CeedScalar interp[num_elem * num_comp * Q * CeedIntPow(P > Q ? P : Q, dim - 1)]; 99 100 // Interpolate to quadrature points (NoTranspose) 101 // or Grad to quadrature points (Transpose) 102 for (CeedInt d = 0; d < dim; d++) { 103 CeedCallBackend(CeedTensorContractApply(contract, pre, P, post, Q, (t_mode == CEED_NOTRANSPOSE ? interp_1d : impl->collo_grad_1d), t_mode, 104 add && (d > 0), 105 (t_mode == CEED_NOTRANSPOSE ? (d == 0 ? u : tmp[d % 2]) : u + d * num_qpts * num_comp * num_elem), 106 (t_mode == CEED_NOTRANSPOSE ? (d == dim - 1 ? interp : tmp[(d + 1) % 2]) : interp))); 107 pre /= P; 108 post *= Q; 109 } 110 // Grad to quadrature points (NoTranspose) 111 // or Interpolate to nodes (Transpose) 112 P = Q_1d, Q = Q_1d; 113 if (t_mode == CEED_TRANSPOSE) { 114 P = Q_1d; 115 Q = P_1d; 116 } 117 pre = num_comp * CeedIntPow(P, dim - 1), post = num_elem; 118 for (CeedInt d = 0; d < dim; d++) { 119 CeedCallBackend(CeedTensorContractApply( 120 contract, pre, P, post, Q, (t_mode == CEED_NOTRANSPOSE ? impl->collo_grad_1d : interp_1d), t_mode, add && (d == dim - 1), 121 (t_mode == CEED_NOTRANSPOSE ? interp : (d == 0 ? interp : tmp[d % 2])), 122 (t_mode == CEED_NOTRANSPOSE ? v + d * num_qpts * num_comp * num_elem : (d == dim - 1 ? v : tmp[(d + 1) % 2])))); 123 pre /= P; 124 post *= Q; 125 } 126 } else if (impl->has_collo_interp) { // Qpts collocated with nodes 127 const CeedScalar *grad_1d; 128 129 CeedCallBackend(CeedBasisGetGrad1D(basis, &grad_1d)); 130 131 // Dim contractions, identity in other directions 132 CeedInt pre = num_comp * CeedIntPow(P, dim - 1), post = num_elem; 133 134 for (CeedInt d = 0; d < dim; d++) { 135 CeedCallBackend(CeedTensorContractApply(contract, pre, P, post, Q, grad_1d, t_mode, add && (d > 0), 136 t_mode == CEED_NOTRANSPOSE ? u : u + d * num_comp * num_qpts * num_elem, 137 t_mode == CEED_TRANSPOSE ? v : v + d * num_comp * num_qpts * num_elem)); 138 pre /= P; 139 post *= Q; 140 } 141 } else { // Underintegration, P > Q 142 const CeedScalar *grad_1d; 143 144 CeedCallBackend(CeedBasisGetGrad1D(basis, &grad_1d)); 145 146 if (t_mode == CEED_TRANSPOSE) { 147 P = Q_1d; 148 Q = P_1d; 149 } 150 CeedScalar tmp[2][num_elem * num_comp * Q * CeedIntPow(P > Q ? P : Q, dim - 1)]; 151 152 // Dim**2 contractions, apply grad when pass == dim 153 for (CeedInt p = 0; p < dim; p++) { 154 CeedInt pre = num_comp * CeedIntPow(P, dim - 1), post = num_elem; 155 156 for (CeedInt d = 0; d < dim; d++) { 157 CeedCallBackend(CeedTensorContractApply( 158 contract, pre, P, post, Q, (p == d) ? grad_1d : interp_1d, t_mode, add && (d == dim - 1), 159 (d == 0 ? (t_mode == CEED_NOTRANSPOSE ? u : u + p * num_comp * num_qpts * num_elem) : tmp[d % 2]), 160 (d == dim - 1 ? (t_mode == CEED_TRANSPOSE ? v : v + p * num_comp * num_qpts * num_elem) : tmp[(d + 1) % 2]))); 161 pre /= P; 162 post *= Q; 163 } 164 } 165 } 166 } break; 167 // Retrieve interpolation weights 168 case CEED_EVAL_WEIGHT: { 169 CeedInt Q = Q_1d; 170 const CeedScalar *q_weight_1d; 171 172 CeedCheck(t_mode == CEED_NOTRANSPOSE, ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 173 CeedCallBackend(CeedBasisGetQWeights(basis, &q_weight_1d)); 174 for (CeedInt d = 0; d < dim; d++) { 175 CeedInt pre = CeedIntPow(Q, dim - d - 1), post = CeedIntPow(Q, d); 176 177 for (CeedInt i = 0; i < pre; i++) { 178 for (CeedInt j = 0; j < Q; j++) { 179 for (CeedInt k = 0; k < post; k++) { 180 const CeedScalar w = q_weight_1d[j] * (d == 0 ? 1 : v[((i * Q + j) * post + k) * num_elem]); 181 182 for (CeedInt e = 0; e < num_elem; e++) v[((i * Q + j) * post + k) * num_elem + e] = w; 183 } 184 } 185 } 186 } 187 } break; 188 // LCOV_EXCL_START 189 case CEED_EVAL_DIV: 190 case CEED_EVAL_CURL: 191 return CeedError(ceed, CEED_ERROR_BACKEND, "%s not supported", CeedEvalModes[eval_mode]); 192 case CEED_EVAL_NONE: 193 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context"); 194 // LCOV_EXCL_STOP 195 } 196 } else { 197 // Non-tensor basis 198 CeedInt P = num_nodes, Q = num_qpts; 199 200 switch (eval_mode) { 201 // Interpolate to/from quadrature points 202 case CEED_EVAL_INTERP: { 203 const CeedScalar *interp; 204 205 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 206 CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, interp, t_mode, add, u, v)); 207 } break; 208 // Evaluate the gradient to/from quadrature points 209 case CEED_EVAL_GRAD: { 210 const CeedScalar *grad; 211 212 CeedCallBackend(CeedBasisGetGrad(basis, &grad)); 213 CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, grad, t_mode, add, u, v)); 214 } break; 215 // Evaluate the divergence to/from the quadrature points 216 case CEED_EVAL_DIV: { 217 const CeedScalar *div; 218 219 CeedCallBackend(CeedBasisGetDiv(basis, &div)); 220 CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, div, t_mode, add, u, v)); 221 } break; 222 // Evaluate the curl to/from the quadrature points 223 case CEED_EVAL_CURL: { 224 const CeedScalar *curl; 225 226 CeedCallBackend(CeedBasisGetCurl(basis, &curl)); 227 CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, curl, t_mode, add, u, v)); 228 } break; 229 // Retrieve interpolation weights 230 case CEED_EVAL_WEIGHT: { 231 const CeedScalar *q_weight; 232 233 CeedCheck(t_mode == CEED_NOTRANSPOSE, ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 234 CeedCallBackend(CeedBasisGetQWeights(basis, &q_weight)); 235 for (CeedInt i = 0; i < num_qpts; i++) { 236 for (CeedInt e = 0; e < num_elem; e++) v[i * num_elem + e] = q_weight[i]; 237 } 238 } break; 239 // LCOV_EXCL_START 240 case CEED_EVAL_NONE: 241 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context"); 242 // LCOV_EXCL_STOP 243 } 244 } 245 if (U != CEED_VECTOR_NONE) { 246 CeedCallBackend(CeedVectorRestoreArrayRead(U, &u)); 247 } 248 CeedCallBackend(CeedVectorRestoreArray(V, &v)); 249 return CEED_ERROR_SUCCESS; 250 } 251 252 //------------------------------------------------------------------------------ 253 // Basis Destroy Tensor 254 //------------------------------------------------------------------------------ 255 static int CeedBasisDestroyTensor_Ref(CeedBasis basis) { 256 CeedBasis_Ref *impl; 257 258 CeedCallBackend(CeedBasisGetData(basis, &impl)); 259 CeedCallBackend(CeedFree(&impl->collo_grad_1d)); 260 CeedCallBackend(CeedFree(&impl)); 261 return CEED_ERROR_SUCCESS; 262 } 263 264 //------------------------------------------------------------------------------ 265 // Basis Create Tensor 266 //------------------------------------------------------------------------------ 267 int CeedBasisCreateTensorH1_Ref(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d, 268 const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis) { 269 Ceed ceed, ceed_parent; 270 CeedBasis_Ref *impl; 271 CeedTensorContract contract; 272 273 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 274 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 275 276 CeedCallBackend(CeedCalloc(1, &impl)); 277 // Check for collocated interp 278 if (Q_1d == P_1d) { 279 bool has_collocated = true; 280 281 for (CeedInt i = 0; i < P_1d; i++) { 282 has_collocated = has_collocated && (fabs(interp_1d[i + P_1d * i] - 1.0) < 1e-14); 283 for (CeedInt j = 0; j < P_1d; j++) { 284 if (j != i) has_collocated = has_collocated && (fabs(interp_1d[j + P_1d * i]) < 1e-14); 285 } 286 } 287 impl->has_collo_interp = has_collocated; 288 } 289 // Calculate collocated grad 290 if (Q_1d >= P_1d && !impl->has_collo_interp) { 291 CeedCallBackend(CeedMalloc(Q_1d * Q_1d, &impl->collo_grad_1d)); 292 CeedCallBackend(CeedBasisGetCollocatedGrad(basis, impl->collo_grad_1d)); 293 } 294 CeedCallBackend(CeedBasisSetData(basis, impl)); 295 296 CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract)); 297 CeedCallBackend(CeedBasisSetTensorContract(basis, contract)); 298 299 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref)); 300 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroyTensor_Ref)); 301 return CEED_ERROR_SUCCESS; 302 } 303 304 //------------------------------------------------------------------------------ 305 // Basis Create Non-Tensor H^1 306 //------------------------------------------------------------------------------ 307 int CeedBasisCreateH1_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *grad, 308 const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) { 309 Ceed ceed, ceed_parent; 310 CeedTensorContract contract; 311 312 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 313 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 314 315 CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract)); 316 CeedCallBackend(CeedBasisSetTensorContract(basis, contract)); 317 318 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref)); 319 return CEED_ERROR_SUCCESS; 320 } 321 322 //------------------------------------------------------------------------------ 323 // Basis Create Non-Tensor H(div) 324 //------------------------------------------------------------------------------ 325 int CeedBasisCreateHdiv_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *div, 326 const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) { 327 Ceed ceed, ceed_parent; 328 CeedTensorContract contract; 329 330 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 331 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 332 333 CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract)); 334 CeedCallBackend(CeedBasisSetTensorContract(basis, contract)); 335 336 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref)); 337 return CEED_ERROR_SUCCESS; 338 } 339 340 //------------------------------------------------------------------------------ 341 // Basis Create Non-Tensor H(curl) 342 //------------------------------------------------------------------------------ 343 int CeedBasisCreateHcurl_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, 344 const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) { 345 Ceed ceed, ceed_parent; 346 CeedTensorContract contract; 347 348 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 349 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 350 351 CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract)); 352 CeedCallBackend(CeedBasisSetTensorContract(basis, contract)); 353 354 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref)); 355 return CEED_ERROR_SUCCESS; 356 } 357 358 //------------------------------------------------------------------------------ 359