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