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 // Evaluate the divergence to/from the quadrature points 190 case CEED_EVAL_DIV: 191 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_DIV not supported"); 192 // Evaluate the curl to/from the quadrature points 193 case CEED_EVAL_CURL: 194 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_CURL not supported"); 195 // Take no action, BasisApply should not have been called 196 case CEED_EVAL_NONE: 197 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context"); 198 // LCOV_EXCL_STOP 199 } 200 } else { 201 // Non-tensor basis 202 CeedInt P = num_nodes, Q = num_qpts; 203 204 switch (eval_mode) { 205 // Interpolate to/from quadrature points 206 case CEED_EVAL_INTERP: { 207 const CeedScalar *interp; 208 209 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 210 CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, interp, t_mode, add, u, v)); 211 } break; 212 // Evaluate the gradient to/from quadrature points 213 case CEED_EVAL_GRAD: { 214 const CeedScalar *grad; 215 216 CeedCallBackend(CeedBasisGetGrad(basis, &grad)); 217 CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, grad, t_mode, add, u, v)); 218 } break; 219 // Evaluate the divergence to/from the quadrature points 220 case CEED_EVAL_DIV: { 221 const CeedScalar *div; 222 223 CeedCallBackend(CeedBasisGetDiv(basis, &div)); 224 CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, div, t_mode, add, u, v)); 225 } break; 226 // Evaluate the curl to/from the quadrature points 227 case CEED_EVAL_CURL: { 228 const CeedScalar *curl; 229 230 CeedCallBackend(CeedBasisGetCurl(basis, &curl)); 231 CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, curl, t_mode, add, u, v)); 232 } break; 233 // Retrieve interpolation weights 234 case CEED_EVAL_WEIGHT: { 235 const CeedScalar *q_weight; 236 237 CeedCheck(t_mode == CEED_NOTRANSPOSE, ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 238 CeedCallBackend(CeedBasisGetQWeights(basis, &q_weight)); 239 for (CeedInt i = 0; i < num_qpts; i++) { 240 for (CeedInt e = 0; e < num_elem; e++) v[i * num_elem + e] = q_weight[i]; 241 } 242 } break; 243 // LCOV_EXCL_START 244 // Take no action, BasisApply should not have been called 245 case CEED_EVAL_NONE: 246 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context"); 247 // LCOV_EXCL_STOP 248 } 249 } 250 if (U != CEED_VECTOR_NONE) { 251 CeedCallBackend(CeedVectorRestoreArrayRead(U, &u)); 252 } 253 CeedCallBackend(CeedVectorRestoreArray(V, &v)); 254 return CEED_ERROR_SUCCESS; 255 } 256 257 //------------------------------------------------------------------------------ 258 // Basis Destroy Tensor 259 //------------------------------------------------------------------------------ 260 static int CeedBasisDestroyTensor_Ref(CeedBasis basis) { 261 CeedBasis_Ref *impl; 262 263 CeedCallBackend(CeedBasisGetData(basis, &impl)); 264 CeedCallBackend(CeedFree(&impl->collo_grad_1d)); 265 CeedCallBackend(CeedFree(&impl)); 266 return CEED_ERROR_SUCCESS; 267 } 268 269 //------------------------------------------------------------------------------ 270 // Basis Create Tensor 271 //------------------------------------------------------------------------------ 272 int CeedBasisCreateTensorH1_Ref(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d, 273 const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis) { 274 Ceed ceed, ceed_parent; 275 CeedBasis_Ref *impl; 276 CeedTensorContract contract; 277 278 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 279 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 280 281 CeedCallBackend(CeedCalloc(1, &impl)); 282 // Check for collocated interp 283 if (Q_1d == P_1d) { 284 bool has_collocated = true; 285 286 for (CeedInt i = 0; i < P_1d; i++) { 287 has_collocated = has_collocated && (fabs(interp_1d[i + P_1d * i] - 1.0) < 1e-14); 288 for (CeedInt j = 0; j < P_1d; j++) { 289 if (j != i) has_collocated = has_collocated && (fabs(interp_1d[j + P_1d * i]) < 1e-14); 290 } 291 } 292 impl->has_collo_interp = has_collocated; 293 } 294 // Calculate collocated grad 295 if (Q_1d >= P_1d && !impl->has_collo_interp) { 296 CeedCallBackend(CeedMalloc(Q_1d * Q_1d, &impl->collo_grad_1d)); 297 CeedCallBackend(CeedBasisGetCollocatedGrad(basis, impl->collo_grad_1d)); 298 } 299 CeedCallBackend(CeedBasisSetData(basis, impl)); 300 301 CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract)); 302 CeedCallBackend(CeedBasisSetTensorContract(basis, contract)); 303 304 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref)); 305 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroyTensor_Ref)); 306 return CEED_ERROR_SUCCESS; 307 } 308 309 //------------------------------------------------------------------------------ 310 // Basis Create Non-Tensor H^1 311 //------------------------------------------------------------------------------ 312 int CeedBasisCreateH1_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *grad, 313 const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) { 314 Ceed ceed, ceed_parent; 315 CeedTensorContract contract; 316 317 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 318 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 319 320 CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract)); 321 CeedCallBackend(CeedBasisSetTensorContract(basis, contract)); 322 323 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref)); 324 return CEED_ERROR_SUCCESS; 325 } 326 327 //------------------------------------------------------------------------------ 328 // Basis Create Non-Tensor H(div) 329 //------------------------------------------------------------------------------ 330 int CeedBasisCreateHdiv_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *div, 331 const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) { 332 Ceed ceed, ceed_parent; 333 CeedTensorContract contract; 334 335 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 336 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 337 338 CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract)); 339 CeedCallBackend(CeedBasisSetTensorContract(basis, contract)); 340 341 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref)); 342 return CEED_ERROR_SUCCESS; 343 } 344 345 //------------------------------------------------------------------------------ 346 // Basis Create Non-Tensor H(curl) 347 //------------------------------------------------------------------------------ 348 int CeedBasisCreateHcurl_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, 349 const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) { 350 Ceed ceed, ceed_parent; 351 CeedTensorContract contract; 352 353 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 354 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 355 356 CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract)); 357 CeedCallBackend(CeedBasisSetTensorContract(basis, contract)); 358 359 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref)); 360 return CEED_ERROR_SUCCESS; 361 } 362 363 //------------------------------------------------------------------------------ 364