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