xref: /libCEED/backends/ref/ceed-ref-basis.c (revision 91216f769229b3793e221c647f948ce6d0249ef6) 
1 // Copyright (c) 2017-2025, 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 CeedBasisApplyCore_Ref(CeedBasis basis, bool apply_add, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U,
20                                   CeedVector V) {
21   bool               is_tensor_basis, add = apply_add || (t_mode == CEED_TRANSPOSE);
22   CeedInt            dim, num_comp, q_comp, num_nodes, num_qpts;
23   const CeedScalar  *u;
24   CeedScalar        *v;
25   CeedTensorContract contract;
26   CeedBasis_Ref     *impl;
27 
28   CeedCallBackend(CeedBasisGetData(basis, &impl));
29   CeedCallBackend(CeedBasisGetDimension(basis, &dim));
30   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
31   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp));
32   CeedCallBackend(CeedBasisGetNumNodes(basis, &num_nodes));
33   CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
34   CeedCallBackend(CeedBasisGetTensorContract(basis, &contract));
35   if (U != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(U, CEED_MEM_HOST, &u));
36   else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, CeedBasisReturnCeed(basis), CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode");
37   // Clear v if operating in transpose
38   if (apply_add) CeedCallBackend(CeedVectorGetArray(V, CEED_MEM_HOST, &v));
39   else CeedCallBackend(CeedVectorGetArrayWrite(V, CEED_MEM_HOST, &v));
40 
41   if (t_mode == CEED_TRANSPOSE && !apply_add) {
42     CeedSize len;
43 
44     CeedCallBackend(CeedVectorGetLength(V, &len));
45     for (CeedInt i = 0; i < len; i++) v[i] = 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                                                     (t_mode == CEED_TRANSPOSE) && (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,
121                 (t_mode == CEED_NOTRANSPOSE && apply_add) || (t_mode == CEED_TRANSPOSE && (d == dim - 1)),
122                 (t_mode == CEED_NOTRANSPOSE ? interp : (d == 0 ? interp : tmp[d % 2])),
123                 (t_mode == CEED_NOTRANSPOSE ? &v[d * num_qpts * num_comp * num_elem] : (d == dim - 1 ? v : tmp[(d + 1) % 2]))));
124             pre /= P;
125             post *= Q;
126           }
127         } else if (impl->has_collo_interp) {  // Qpts collocated with nodes
128           const CeedScalar *grad_1d;
129 
130           CeedCallBackend(CeedBasisGetGrad1D(basis, &grad_1d));
131 
132           // Dim contractions, identity in other directions
133           CeedInt pre = num_comp * CeedIntPow(P, dim - 1), post = num_elem;
134 
135           for (CeedInt d = 0; d < dim; d++) {
136             CeedCallBackend(CeedTensorContractApply(contract, pre, P, post, Q, grad_1d, t_mode, add && (d > 0),
137                                                     t_mode == CEED_NOTRANSPOSE ? u : &u[d * num_comp * num_qpts * num_elem],
138                                                     t_mode == CEED_TRANSPOSE ? v : &v[d * num_comp * num_qpts * num_elem]));
139             pre /= P;
140             post *= Q;
141           }
142         } else {  // Underintegration, P > Q
143           const CeedScalar *grad_1d;
144 
145           CeedCallBackend(CeedBasisGetGrad1D(basis, &grad_1d));
146 
147           if (t_mode == CEED_TRANSPOSE) {
148             P = Q_1d;
149             Q = P_1d;
150           }
151           CeedScalar tmp[2][num_elem * num_comp * Q * CeedIntPow(P > Q ? P : Q, dim - 1)];
152 
153           // Dim**2 contractions, apply grad when pass == dim
154           for (CeedInt p = 0; p < dim; p++) {
155             CeedInt pre = num_comp * CeedIntPow(P, dim - 1), post = num_elem;
156 
157             for (CeedInt d = 0; d < dim; d++) {
158               CeedCallBackend(CeedTensorContractApply(
159                   contract, pre, P, post, Q, (p == d) ? grad_1d : interp_1d, t_mode, add && (d == dim - 1),
160                   (d == 0 ? (t_mode == CEED_NOTRANSPOSE ? u : &u[p * num_comp * num_qpts * num_elem]) : tmp[d % 2]),
161                   (d == dim - 1 ? (t_mode == CEED_TRANSPOSE ? v : &v[p * num_comp * num_qpts * num_elem]) : tmp[(d + 1) % 2])));
162               pre /= P;
163               post *= Q;
164             }
165           }
166         }
167       } break;
168       // Retrieve interpolation weights
169       case CEED_EVAL_WEIGHT: {
170         CeedInt           Q = Q_1d;
171         const CeedScalar *q_weight_1d;
172 
173         CeedCheck(t_mode == CEED_NOTRANSPOSE, CeedBasisReturnCeed(basis), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE");
174         CeedCallBackend(CeedBasisGetQWeights(basis, &q_weight_1d));
175         for (CeedInt d = 0; d < dim; d++) {
176           CeedInt pre = CeedIntPow(Q, dim - d - 1), post = CeedIntPow(Q, d);
177 
178           for (CeedInt i = 0; i < pre; i++) {
179             for (CeedInt j = 0; j < Q; j++) {
180               for (CeedInt k = 0; k < post; k++) {
181                 const CeedScalar w = q_weight_1d[j] * (d == 0 ? 1 : v[((i * Q + j) * post + k) * num_elem]);
182 
183                 for (CeedInt e = 0; e < num_elem; e++) v[((i * Q + j) * post + k) * num_elem + e] = w;
184               }
185             }
186           }
187         }
188       } break;
189       // LCOV_EXCL_START
190       case CEED_EVAL_DIV:
191       case CEED_EVAL_CURL:
192         return CeedError(CeedBasisReturnCeed(basis), CEED_ERROR_BACKEND, "%s not supported", CeedEvalModes[eval_mode]);
193       case CEED_EVAL_NONE:
194         return CeedError(CeedBasisReturnCeed(basis), CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context");
195         // LCOV_EXCL_STOP
196     }
197   } else {
198     // Non-tensor basis
199     CeedInt P = num_nodes, Q = num_qpts;
200 
201     switch (eval_mode) {
202       // Interpolate to/from quadrature points
203       case CEED_EVAL_INTERP: {
204         const CeedScalar *interp;
205 
206         CeedCallBackend(CeedBasisGetInterp(basis, &interp));
207         CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, interp, t_mode, add, u, v));
208       } break;
209       // Evaluate the gradient to/from quadrature points
210       case CEED_EVAL_GRAD: {
211         const CeedScalar *grad;
212 
213         CeedCallBackend(CeedBasisGetGrad(basis, &grad));
214         CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, grad, t_mode, add, u, v));
215       } break;
216       // Evaluate the divergence to/from the quadrature points
217       case CEED_EVAL_DIV: {
218         const CeedScalar *div;
219 
220         CeedCallBackend(CeedBasisGetDiv(basis, &div));
221         CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, div, t_mode, add, u, v));
222       } break;
223       // Evaluate the curl to/from the quadrature points
224       case CEED_EVAL_CURL: {
225         const CeedScalar *curl;
226 
227         CeedCallBackend(CeedBasisGetCurl(basis, &curl));
228         CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, curl, t_mode, add, u, v));
229       } break;
230       // Retrieve interpolation weights
231       case CEED_EVAL_WEIGHT: {
232         const CeedScalar *q_weight;
233 
234         CeedCheck(t_mode == CEED_NOTRANSPOSE, CeedBasisReturnCeed(basis), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE");
235         CeedCallBackend(CeedBasisGetQWeights(basis, &q_weight));
236         for (CeedInt i = 0; i < num_qpts; i++) {
237           for (CeedInt e = 0; e < num_elem; e++) v[i * num_elem + e] = q_weight[i];
238         }
239       } break;
240       // LCOV_EXCL_START
241       case CEED_EVAL_NONE:
242         return CeedError(CeedBasisReturnCeed(basis), CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context");
243         // LCOV_EXCL_STOP
244     }
245   }
246   if (U != CEED_VECTOR_NONE) {
247     CeedCallBackend(CeedVectorRestoreArrayRead(U, &u));
248   }
249   CeedCallBackend(CeedVectorRestoreArray(V, &v));
250   return CEED_ERROR_SUCCESS;
251 }
252 
253 static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U, CeedVector V) {
254   CeedCallBackend(CeedBasisApplyCore_Ref(basis, false, num_elem, t_mode, eval_mode, U, V));
255   return CEED_ERROR_SUCCESS;
256 }
257 
258 static int CeedBasisApplyAdd_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U, CeedVector V) {
259   CeedCallBackend(CeedBasisApplyCore_Ref(basis, true, num_elem, t_mode, eval_mode, U, V));
260   return CEED_ERROR_SUCCESS;
261 }
262 
263 //------------------------------------------------------------------------------
264 // Basis Destroy Tensor
265 //------------------------------------------------------------------------------
266 static int CeedBasisDestroyTensor_Ref(CeedBasis basis) {
267   CeedBasis_Ref *impl;
268 
269   CeedCallBackend(CeedBasisGetData(basis, &impl));
270   CeedCallBackend(CeedFree(&impl->collo_grad_1d));
271   CeedCallBackend(CeedFree(&impl));
272   return CEED_ERROR_SUCCESS;
273 }
274 
275 //------------------------------------------------------------------------------
276 // Basis Create Tensor
277 //------------------------------------------------------------------------------
278 int CeedBasisCreateTensorH1_Ref(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d,
279                                 const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis) {
280   Ceed               ceed, ceed_parent;
281   CeedBasis_Ref     *impl;
282   CeedTensorContract contract;
283 
284   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
285   CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
286 
287   CeedCallBackend(CeedCalloc(1, &impl));
288   // Check for collocated interp
289   if (Q_1d == P_1d) {
290     bool has_collocated = true;
291 
292     for (CeedInt i = 0; i < P_1d; i++) {
293       has_collocated = has_collocated && (fabs(interp_1d[i + P_1d * i] - 1.0) < 1e-14);
294       for (CeedInt j = 0; j < P_1d; j++) {
295         if (j != i) has_collocated = has_collocated && (fabs(interp_1d[j + P_1d * i]) < 1e-14);
296       }
297     }
298     impl->has_collo_interp = has_collocated;
299   }
300   // Calculate collocated grad
301   if (Q_1d >= P_1d && !impl->has_collo_interp) {
302     CeedCallBackend(CeedMalloc(Q_1d * Q_1d, &impl->collo_grad_1d));
303     CeedCallBackend(CeedBasisGetCollocatedGrad(basis, impl->collo_grad_1d));
304   }
305   CeedCallBackend(CeedBasisSetData(basis, impl));
306 
307   CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract));
308   CeedCallBackend(CeedBasisSetTensorContract(basis, contract));
309 
310   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref));
311   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAdd_Ref));
312   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroyTensor_Ref));
313   CeedCallBackend(CeedDestroy(&ceed));
314   CeedCallBackend(CeedDestroy(&ceed_parent));
315   return CEED_ERROR_SUCCESS;
316 }
317 
318 //------------------------------------------------------------------------------
319 // Basis Create Non-Tensor H^1
320 //------------------------------------------------------------------------------
321 int CeedBasisCreateH1_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *grad,
322                           const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
323   Ceed               ceed, ceed_parent;
324   CeedTensorContract contract;
325 
326   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
327   CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
328 
329   CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract));
330   CeedCallBackend(CeedBasisSetTensorContract(basis, contract));
331 
332   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref));
333   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAdd_Ref));
334   CeedCallBackend(CeedDestroy(&ceed));
335   CeedCallBackend(CeedDestroy(&ceed_parent));
336   return CEED_ERROR_SUCCESS;
337 }
338 
339 //------------------------------------------------------------------------------
340 // Basis Create Non-Tensor H(div)
341 //------------------------------------------------------------------------------
342 int CeedBasisCreateHdiv_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *div,
343                             const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
344   Ceed               ceed, ceed_parent;
345   CeedTensorContract contract;
346 
347   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
348   CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
349 
350   CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract));
351   CeedCallBackend(CeedBasisSetTensorContract(basis, contract));
352 
353   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref));
354   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAdd_Ref));
355   CeedCallBackend(CeedDestroy(&ceed));
356   CeedCallBackend(CeedDestroy(&ceed_parent));
357   return CEED_ERROR_SUCCESS;
358 }
359 
360 //------------------------------------------------------------------------------
361 // Basis Create Non-Tensor H(curl)
362 //------------------------------------------------------------------------------
363 int CeedBasisCreateHcurl_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp,
364                              const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
365   Ceed               ceed, ceed_parent;
366   CeedTensorContract contract;
367 
368   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
369   CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
370 
371   CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract));
372   CeedCallBackend(CeedBasisSetTensorContract(basis, contract));
373 
374   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref));
375   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAdd_Ref));
376   CeedCallBackend(CeedDestroy(&ceed));
377   CeedCallBackend(CeedDestroy(&ceed_parent));
378   return CEED_ERROR_SUCCESS;
379 }
380 
381 //------------------------------------------------------------------------------
382