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