1 // Copyright (c) 2017-2026, 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-impl.h>
9 #include <ceed.h>
10 #include <ceed/backend.h>
11 #include <math.h>
12 #include <stdbool.h>
13 #include <stdio.h>
14 #include <string.h>
15
16 /// @file
17 /// Implementation of CeedBasis interfaces
18
19 /// @cond DOXYGEN_SKIP
20 static struct CeedBasis_private ceed_basis_none;
21 /// @endcond
22
23 /// @addtogroup CeedBasisUser
24 /// @{
25
26 /// Argument for @ref CeedOperatorSetField() indicating that the field does not require a `CeedBasis`
27 const CeedBasis CEED_BASIS_NONE = &ceed_basis_none;
28
29 /// @}
30
31 /// ----------------------------------------------------------------------------
32 /// CeedBasis Library Internal Functions
33 /// ----------------------------------------------------------------------------
34 /// @addtogroup CeedBasisDeveloper
35 /// @{
36
37 /**
38 @brief Compute Chebyshev polynomial values at a point
39
40 @param[in] x Coordinate to evaluate Chebyshev polynomials at
41 @param[in] n Number of Chebyshev polynomials to evaluate, `n >= 2`
42 @param[out] chebyshev_x Array of Chebyshev polynomial values
43
44 @return An error code: 0 - success, otherwise - failure
45
46 @ref Developer
47 **/
CeedChebyshevPolynomialsAtPoint(CeedScalar x,CeedInt n,CeedScalar * chebyshev_x)48 static int CeedChebyshevPolynomialsAtPoint(CeedScalar x, CeedInt n, CeedScalar *chebyshev_x) {
49 chebyshev_x[0] = 1.0;
50 chebyshev_x[1] = 2 * x;
51 for (CeedInt i = 2; i < n; i++) chebyshev_x[i] = 2 * x * chebyshev_x[i - 1] - chebyshev_x[i - 2];
52 return CEED_ERROR_SUCCESS;
53 }
54
55 /**
56 @brief Compute values of the derivative of Chebyshev polynomials at a point
57
58 @param[in] x Coordinate to evaluate derivative of Chebyshev polynomials at
59 @param[in] n Number of Chebyshev polynomials to evaluate, `n >= 2`
60 @param[out] chebyshev_dx Array of Chebyshev polynomial derivative values
61
62 @return An error code: 0 - success, otherwise - failure
63
64 @ref Developer
65 **/
CeedChebyshevDerivativeAtPoint(CeedScalar x,CeedInt n,CeedScalar * chebyshev_dx)66 static int CeedChebyshevDerivativeAtPoint(CeedScalar x, CeedInt n, CeedScalar *chebyshev_dx) {
67 CeedScalar chebyshev_x[3];
68
69 chebyshev_x[1] = 1.0;
70 chebyshev_x[2] = 2 * x;
71 chebyshev_dx[0] = 0.0;
72 chebyshev_dx[1] = 2.0;
73 for (CeedInt i = 2; i < n; i++) {
74 chebyshev_x[0] = chebyshev_x[1];
75 chebyshev_x[1] = chebyshev_x[2];
76 chebyshev_x[2] = 2 * x * chebyshev_x[1] - chebyshev_x[0];
77 chebyshev_dx[i] = 2 * x * chebyshev_dx[i - 1] + 2 * chebyshev_x[1] - chebyshev_dx[i - 2];
78 }
79 return CEED_ERROR_SUCCESS;
80 }
81
82 /**
83 @brief Compute Householder reflection.
84
85 Computes \f$A = (I - b v v^T) A\f$, where \f$A\f$ is an \f$m \times n\f$ matrix indexed as `A[i*row + j*col]`.
86
87 @param[in,out] A Matrix to apply Householder reflection to, in place
88 @param[in] v Householder vector
89 @param[in] b Scaling factor
90 @param[in] m Number of rows in `A`
91 @param[in] n Number of columns in `A`
92 @param[in] row Row stride
93 @param[in] col Col stride
94
95 @return An error code: 0 - success, otherwise - failure
96
97 @ref Developer
98 **/
CeedHouseholderReflect(CeedScalar * A,const CeedScalar * v,CeedScalar b,CeedInt m,CeedInt n,CeedInt row,CeedInt col)99 static int CeedHouseholderReflect(CeedScalar *A, const CeedScalar *v, CeedScalar b, CeedInt m, CeedInt n, CeedInt row, CeedInt col) {
100 for (CeedInt j = 0; j < n; j++) {
101 CeedScalar w = A[0 * row + j * col];
102
103 for (CeedInt i = 1; i < m; i++) w += v[i] * A[i * row + j * col];
104 A[0 * row + j * col] -= b * w;
105 for (CeedInt i = 1; i < m; i++) A[i * row + j * col] -= b * w * v[i];
106 }
107 return CEED_ERROR_SUCCESS;
108 }
109
110 /**
111 @brief Compute Givens rotation
112
113 Computes \f$A = G A\f$ (or \f$G^T A\f$ in transpose mode), where \f$A\f$ is an \f$m \times n\f$ matrix indexed as `A[i*n + j*m]`.
114
115 @param[in,out] A Row major matrix to apply Givens rotation to, in place
116 @param[in] c Cosine factor
117 @param[in] s Sine factor
118 @param[in] t_mode @ref CEED_NOTRANSPOSE to rotate the basis counter-clockwise, which has the effect of rotating columns of `A` clockwise;
119 @ref CEED_TRANSPOSE for the opposite rotation
120 @param[in] i First row/column to apply rotation
121 @param[in] k Second row/column to apply rotation
122 @param[in] m Number of rows in `A`
123 @param[in] n Number of columns in `A`
124
125 @return An error code: 0 - success, otherwise - failure
126
127 @ref Developer
128 **/
CeedGivensRotation(CeedScalar * A,CeedScalar c,CeedScalar s,CeedTransposeMode t_mode,CeedInt i,CeedInt k,CeedInt m,CeedInt n)129 static int CeedGivensRotation(CeedScalar *A, CeedScalar c, CeedScalar s, CeedTransposeMode t_mode, CeedInt i, CeedInt k, CeedInt m, CeedInt n) {
130 CeedInt stride_j = 1, stride_ik = m, num_its = n;
131
132 if (t_mode == CEED_NOTRANSPOSE) {
133 stride_j = n;
134 stride_ik = 1;
135 num_its = m;
136 }
137
138 // Apply rotation
139 for (CeedInt j = 0; j < num_its; j++) {
140 CeedScalar tau1 = A[i * stride_ik + j * stride_j], tau2 = A[k * stride_ik + j * stride_j];
141
142 A[i * stride_ik + j * stride_j] = c * tau1 - s * tau2;
143 A[k * stride_ik + j * stride_j] = s * tau1 + c * tau2;
144 }
145 return CEED_ERROR_SUCCESS;
146 }
147
148 /**
149 @brief View an array stored in a `CeedBasis`
150
151 @param[in] name Name of array
152 @param[in] fp_fmt Printing format
153 @param[in] m Number of rows in array
154 @param[in] n Number of columns in array
155 @param[in] a Array to be viewed
156 @param[in] tabs Tabs to append before each new line
157 @param[in] stream Stream to view to, e.g., `stdout`
158
159 @return An error code: 0 - success, otherwise - failure
160
161 @ref Developer
162 **/
CeedScalarView(const char * name,const char * fp_fmt,CeedInt m,CeedInt n,const CeedScalar * a,const char * tabs,FILE * stream)163 static int CeedScalarView(const char *name, const char *fp_fmt, CeedInt m, CeedInt n, const CeedScalar *a, const char *tabs, FILE *stream) {
164 if (m > 1) {
165 fprintf(stream, "%s %s:\n", tabs, name);
166 } else {
167 char padded_name[12];
168
169 snprintf(padded_name, 11, "%s:", name);
170 fprintf(stream, "%s %-10s", tabs, padded_name);
171 }
172 for (CeedInt i = 0; i < m; i++) {
173 if (m > 1) fprintf(stream, "%s [%" CeedInt_FMT "]", tabs, i);
174 for (CeedInt j = 0; j < n; j++) fprintf(stream, fp_fmt, fabs(a[i * n + j]) > 1E-14 ? a[i * n + j] : 0);
175 fputs("\n", stream);
176 }
177 return CEED_ERROR_SUCCESS;
178 }
179
180 /**
181 @brief View a `CeedBasis` passed as a `CeedObject`
182
183 @param[in] basis `CeedBasis` to view
184 @param[in] stream Filestream to write to
185
186 @return An error code: 0 - success, otherwise - failure
187
188 @ref Developer
189 **/
CeedBasisView_Object(CeedObject basis,FILE * stream)190 static int CeedBasisView_Object(CeedObject basis, FILE *stream) {
191 CeedCall(CeedBasisView((CeedBasis)basis, stream));
192 return CEED_ERROR_SUCCESS;
193 }
194
195 /**
196 @brief Destroy a `CeedBasis` passed as a `CeedObject`
197
198 @param[in,out] basis Address of `CeedBasis` to destroy
199
200 @return An error code: 0 - success, otherwise - failure
201
202 @ref Developer
203 **/
CeedBasisDestroy_Object(CeedObject * basis)204 static int CeedBasisDestroy_Object(CeedObject *basis) {
205 CeedCall(CeedBasisDestroy((CeedBasis *)basis));
206 return CEED_ERROR_SUCCESS;
207 }
208
209 /**
210 @brief Create the interpolation and gradient matrices for projection from the nodes of `basis_from` to the nodes of `basis_to`.
211
212 The interpolation is given by `interp_project = interp_to^+ * interp_from`, where the pseudoinverse `interp_to^+` is given by QR factorization.
213 The gradient is given by `grad_project = interp_to^+ * grad_from`, and is only computed for \f$H^1\f$ spaces otherwise it should not be used.
214
215 Note: `basis_from` and `basis_to` must have compatible quadrature spaces.
216
217 @param[in] basis_from `CeedBasis` to project from
218 @param[in] basis_to `CeedBasis` to project to
219 @param[out] interp_project Address of the variable where the newly created interpolation matrix will be stored
220 @param[out] grad_project Address of the variable where the newly created gradient matrix will be stored
221
222 @return An error code: 0 - success, otherwise - failure
223
224 @ref Developer
225 **/
CeedBasisCreateProjectionMatrices(CeedBasis basis_from,CeedBasis basis_to,CeedScalar ** interp_project,CeedScalar ** grad_project)226 static int CeedBasisCreateProjectionMatrices(CeedBasis basis_from, CeedBasis basis_to, CeedScalar **interp_project, CeedScalar **grad_project) {
227 bool are_both_tensor;
228 CeedInt Q, Q_to, Q_from, P_to, P_from;
229
230 // Check for compatible quadrature spaces
231 CeedCall(CeedBasisGetNumQuadraturePoints(basis_to, &Q_to));
232 CeedCall(CeedBasisGetNumQuadraturePoints(basis_from, &Q_from));
233 CeedCheck(Q_to == Q_from, CeedBasisReturnCeed(basis_to), CEED_ERROR_DIMENSION,
234 "Bases must have compatible quadrature spaces."
235 " 'basis_from' has %" CeedInt_FMT " points and 'basis_to' has %" CeedInt_FMT,
236 Q_from, Q_to);
237 Q = Q_to;
238
239 // Check for matching tensor or non-tensor
240 {
241 bool is_tensor_to, is_tensor_from;
242
243 CeedCall(CeedBasisIsTensor(basis_to, &is_tensor_to));
244 CeedCall(CeedBasisIsTensor(basis_from, &is_tensor_from));
245 are_both_tensor = is_tensor_to && is_tensor_from;
246 }
247 if (are_both_tensor) {
248 CeedCall(CeedBasisGetNumNodes1D(basis_to, &P_to));
249 CeedCall(CeedBasisGetNumNodes1D(basis_from, &P_from));
250 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis_from, &Q));
251 } else {
252 CeedCall(CeedBasisGetNumNodes(basis_to, &P_to));
253 CeedCall(CeedBasisGetNumNodes(basis_from, &P_from));
254 }
255
256 // Check for matching FE space
257 CeedFESpace fe_space_to, fe_space_from;
258
259 CeedCall(CeedBasisGetFESpace(basis_to, &fe_space_to));
260 CeedCall(CeedBasisGetFESpace(basis_from, &fe_space_from));
261 CeedCheck(fe_space_to == fe_space_from, CeedBasisReturnCeed(basis_to), CEED_ERROR_MINOR,
262 "Bases must both be the same FE space type."
263 " 'basis_from' is a %s and 'basis_to' is a %s",
264 CeedFESpaces[fe_space_from], CeedFESpaces[fe_space_to]);
265
266 // Get source matrices
267 CeedInt dim, q_comp = 1;
268 CeedScalar *interp_to_inv, *interp_from;
269 const CeedScalar *interp_to_source = NULL, *interp_from_source = NULL, *grad_from_source = NULL;
270
271 CeedCall(CeedBasisGetDimension(basis_from, &dim));
272 if (are_both_tensor) {
273 CeedCall(CeedBasisGetInterp1D(basis_to, &interp_to_source));
274 CeedCall(CeedBasisGetInterp1D(basis_from, &interp_from_source));
275 } else {
276 CeedCall(CeedBasisGetNumQuadratureComponents(basis_from, CEED_EVAL_INTERP, &q_comp));
277 CeedCall(CeedBasisGetInterp(basis_to, &interp_to_source));
278 CeedCall(CeedBasisGetInterp(basis_from, &interp_from_source));
279 }
280 CeedCall(CeedMalloc(Q * P_from * q_comp, &interp_from));
281 CeedCall(CeedCalloc(P_to * P_from, interp_project));
282
283 // `grad_project = interp_to^+ * grad_from` is computed for the H^1 space case so the
284 // projection basis will have a gradient operation (allocated even if not H^1 for the
285 // basis construction later on)
286 if (fe_space_to == CEED_FE_SPACE_H1) {
287 if (are_both_tensor) {
288 CeedCall(CeedBasisGetGrad1D(basis_from, &grad_from_source));
289 } else {
290 CeedCall(CeedBasisGetGrad(basis_from, &grad_from_source));
291 }
292 }
293 CeedCall(CeedCalloc(P_to * P_from * (are_both_tensor ? 1 : dim), grad_project));
294
295 // Compute interp_to^+, pseudoinverse of interp_to
296 CeedCall(CeedCalloc(Q * q_comp * P_to, &interp_to_inv));
297 CeedCall(CeedMatrixPseudoinverse(CeedBasisReturnCeed(basis_to), interp_to_source, Q * q_comp, P_to, interp_to_inv));
298 // Build matrices
299 CeedInt num_matrices = 1 + (fe_space_to == CEED_FE_SPACE_H1) * (are_both_tensor ? 1 : dim);
300 CeedScalar *input_from[num_matrices], *output_project[num_matrices];
301
302 input_from[0] = (CeedScalar *)interp_from_source;
303 output_project[0] = *interp_project;
304 for (CeedInt m = 1; m < num_matrices; m++) {
305 input_from[m] = (CeedScalar *)&grad_from_source[(m - 1) * Q * P_from];
306 output_project[m] = &((*grad_project)[(m - 1) * P_to * P_from]);
307 }
308 for (CeedInt m = 0; m < num_matrices; m++) {
309 // output_project = interp_to^+ * interp_from
310 memcpy(interp_from, input_from[m], Q * P_from * q_comp * sizeof(input_from[m][0]));
311 CeedCall(CeedMatrixMatrixMultiply(CeedBasisReturnCeed(basis_to), interp_to_inv, input_from[m], output_project[m], P_to, P_from, Q * q_comp));
312 // Round zero to machine precision
313 for (CeedInt i = 0; i < P_to * P_from; i++) {
314 if (fabs(output_project[m][i]) < 10 * CEED_EPSILON) output_project[m][i] = 0.0;
315 }
316 }
317
318 // Cleanup
319 CeedCall(CeedFree(&interp_to_inv));
320 CeedCall(CeedFree(&interp_from));
321 return CEED_ERROR_SUCCESS;
322 }
323
324 /**
325 @brief Check input vector dimensions for CeedBasisApply[Add]
326
327 @param[in] basis `CeedBasis` to evaluate
328 @param[in] num_elem The number of elements to apply the basis evaluation to;
329 the backend will specify the ordering in @ref CeedElemRestrictionCreate()
330 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to quadrature points;
331 @ref CEED_TRANSPOSE to apply the transpose, mapping from quadrature points to nodes
332 @param[in] eval_mode @ref CEED_EVAL_NONE to use values directly,
333 @ref CEED_EVAL_INTERP to use interpolated values,
334 @ref CEED_EVAL_GRAD to use gradients,
335 @ref CEED_EVAL_DIV to use divergence,
336 @ref CEED_EVAL_CURL to use curl,
337 @ref CEED_EVAL_WEIGHT to use quadrature weights
338 @param[in] u Input `CeedVector`
339 @param[out] v Output `CeedVector`
340
341 @return An error code: 0 - success, otherwise - failure
342
343 @ref Developer
344 **/
CeedBasisApplyCheckDims(CeedBasis basis,CeedInt num_elem,CeedTransposeMode t_mode,CeedEvalMode eval_mode,CeedVector u,CeedVector v)345 static int CeedBasisApplyCheckDims(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) {
346 CeedInt dim, num_comp, q_comp, num_nodes, num_qpts;
347 CeedSize u_length = 0, v_length;
348
349 CeedCall(CeedBasisGetDimension(basis, &dim));
350 CeedCall(CeedBasisGetNumComponents(basis, &num_comp));
351 CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp));
352 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes));
353 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
354 CeedCall(CeedVectorGetLength(v, &v_length));
355 if (u) CeedCall(CeedVectorGetLength(u, &u_length));
356
357 // Check vector lengths to prevent out of bounds issues
358 bool has_good_dims = true;
359 switch (eval_mode) {
360 case CEED_EVAL_NONE:
361 case CEED_EVAL_INTERP:
362 case CEED_EVAL_GRAD:
363 case CEED_EVAL_DIV:
364 case CEED_EVAL_CURL:
365 has_good_dims = ((t_mode == CEED_TRANSPOSE && u_length >= (CeedSize)num_elem * (CeedSize)num_comp * (CeedSize)num_qpts * (CeedSize)q_comp &&
366 v_length >= (CeedSize)num_elem * (CeedSize)num_comp * (CeedSize)num_nodes) ||
367 (t_mode == CEED_NOTRANSPOSE && v_length >= (CeedSize)num_elem * (CeedSize)num_qpts * (CeedSize)num_comp * (CeedSize)q_comp &&
368 u_length >= (CeedSize)num_elem * (CeedSize)num_comp * (CeedSize)num_nodes));
369 break;
370 case CEED_EVAL_WEIGHT:
371 has_good_dims = v_length >= (CeedSize)num_elem * (CeedSize)num_qpts;
372 break;
373 }
374 CeedCheck(has_good_dims, CeedBasisReturnCeed(basis), CEED_ERROR_DIMENSION, "Input/output vectors too short for basis and evaluation mode");
375 return CEED_ERROR_SUCCESS;
376 }
377
378 /**
379 @brief Check input vector dimensions for CeedBasisApply[Add]AtPoints
380
381 @param[in] basis `CeedBasis` to evaluate
382 @param[in] num_elem The number of elements to apply the basis evaluation to;
383 the backend will specify the ordering in @ref CeedElemRestrictionCreate()
384 @param[in] num_points Array of the number of points to apply the basis evaluation to in each element, size `num_elem`
385 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points;
386 @ref CEED_TRANSPOSE to apply the transpose, mapping from points to nodes
387 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values,
388 @ref CEED_EVAL_GRAD to use gradients,
389 @ref CEED_EVAL_WEIGHT to use quadrature weights
390 @param[in] x_ref `CeedVector` holding reference coordinates of each point
391 @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE
392 @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP
393
394 @return An error code: 0 - success, otherwise - failure
395
396 @ref Developer
397 **/
CeedBasisApplyAtPointsCheckDims(CeedBasis basis,CeedInt num_elem,const CeedInt * num_points,CeedTransposeMode t_mode,CeedEvalMode eval_mode,CeedVector x_ref,CeedVector u,CeedVector v)398 static int CeedBasisApplyAtPointsCheckDims(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode,
399 CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) {
400 CeedInt dim, num_comp, num_q_comp, num_nodes, P_1d = 1, Q_1d = 1, total_num_points = 0;
401 CeedSize x_length = 0, u_length = 0, v_length;
402
403 CeedCall(CeedBasisGetDimension(basis, &dim));
404 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d));
405 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
406 CeedCall(CeedBasisGetNumComponents(basis, &num_comp));
407 CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &num_q_comp));
408 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes));
409 CeedCall(CeedVectorGetLength(v, &v_length));
410 if (x_ref != CEED_VECTOR_NONE) CeedCall(CeedVectorGetLength(x_ref, &x_length));
411 if (u != CEED_VECTOR_NONE) CeedCall(CeedVectorGetLength(u, &u_length));
412
413 // Check compatibility coordinates vector
414 for (CeedInt i = 0; i < num_elem; i++) total_num_points += num_points[i];
415 CeedCheck((x_length >= (CeedSize)total_num_points * (CeedSize)dim) || (eval_mode == CEED_EVAL_WEIGHT), CeedBasisReturnCeed(basis),
416 CEED_ERROR_DIMENSION,
417 "Length of reference coordinate vector incompatible with basis dimension and number of points."
418 " Found reference coordinate vector of length %" CeedSize_FMT ", not of length %" CeedSize_FMT ".",
419 x_length, (CeedSize)total_num_points * (CeedSize)dim);
420
421 // Check CEED_EVAL_WEIGHT only on CEED_NOTRANSPOSE
422 CeedCheck(eval_mode != CEED_EVAL_WEIGHT || t_mode == CEED_NOTRANSPOSE, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED,
423 "CEED_EVAL_WEIGHT only supported with CEED_NOTRANSPOSE");
424
425 // Check vector lengths to prevent out of bounds issues
426 bool has_good_dims = true;
427 switch (eval_mode) {
428 case CEED_EVAL_INTERP:
429 has_good_dims = ((t_mode == CEED_TRANSPOSE && (u_length >= (CeedSize)total_num_points * (CeedSize)num_q_comp ||
430 v_length >= (CeedSize)num_elem * (CeedSize)num_nodes * (CeedSize)num_comp)) ||
431 (t_mode == CEED_NOTRANSPOSE && (v_length >= (CeedSize)total_num_points * (CeedSize)num_q_comp ||
432 u_length >= (CeedSize)num_elem * (CeedSize)num_nodes * (CeedSize)num_comp)));
433 break;
434 case CEED_EVAL_GRAD:
435 has_good_dims = ((t_mode == CEED_TRANSPOSE && (u_length >= (CeedSize)total_num_points * (CeedSize)num_q_comp * (CeedSize)dim ||
436 v_length >= (CeedSize)num_elem * (CeedSize)num_nodes * (CeedSize)num_comp)) ||
437 (t_mode == CEED_NOTRANSPOSE && (v_length >= (CeedSize)total_num_points * (CeedSize)num_q_comp * (CeedSize)dim ||
438 u_length >= (CeedSize)num_elem * (CeedSize)num_nodes * (CeedSize)num_comp)));
439 break;
440 case CEED_EVAL_WEIGHT:
441 has_good_dims = t_mode == CEED_NOTRANSPOSE && (v_length >= total_num_points);
442 break;
443 // LCOV_EXCL_START
444 case CEED_EVAL_NONE:
445 case CEED_EVAL_DIV:
446 case CEED_EVAL_CURL:
447 return CeedError(CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Evaluation at arbitrary points not supported for %s",
448 CeedEvalModes[eval_mode]);
449 // LCOV_EXCL_STOP
450 }
451 CeedCheck(has_good_dims, CeedBasisReturnCeed(basis), CEED_ERROR_DIMENSION, "Input/output vectors too short for basis and evaluation mode");
452 return CEED_ERROR_SUCCESS;
453 }
454
455 /**
456 @brief Default implimentation to apply basis evaluation from nodes to arbitrary points
457
458 @param[in] basis `CeedBasis` to evaluate
459 @param[in] apply_add Sum result into target vector or overwrite
460 @param[in] num_elem The number of elements to apply the basis evaluation to;
461 the backend will specify the ordering in @ref CeedElemRestrictionCreate()
462 @param[in] num_points Array of the number of points to apply the basis evaluation to in each element, size `num_elem`
463 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points;
464 @ref CEED_TRANSPOSE to apply the transpose, mapping from points to nodes
465 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values,
466 @ref CEED_EVAL_GRAD to use gradients,
467 @ref CEED_EVAL_WEIGHT to use quadrature weights
468 @param[in] x_ref `CeedVector` holding reference coordinates of each point
469 @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE
470 @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP
471
472 @return An error code: 0 - success, otherwise - failure
473
474 @ref Developer
475 **/
CeedBasisApplyAtPoints_Core(CeedBasis basis,bool apply_add,CeedInt num_elem,const CeedInt * num_points,CeedTransposeMode t_mode,CeedEvalMode eval_mode,CeedVector x_ref,CeedVector u,CeedVector v)476 static int CeedBasisApplyAtPoints_Core(CeedBasis basis, bool apply_add, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode,
477 CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) {
478 CeedInt dim, num_comp, P_1d = 1, Q_1d = 1, total_num_points = num_points[0];
479
480 CeedCall(CeedBasisGetDimension(basis, &dim));
481 // Inserting check because clang-tidy doesn't understand this cannot occur
482 CeedCheck(dim > 0, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Malformed CeedBasis, dim > 0 is required");
483 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d));
484 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
485 CeedCall(CeedBasisGetNumComponents(basis, &num_comp));
486
487 // Default implementation
488 {
489 bool is_tensor_basis;
490
491 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis));
492 CeedCheck(is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED,
493 "Evaluation at arbitrary points only supported for tensor product bases");
494 }
495 CeedCheck(num_elem == 1, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED,
496 "Evaluation at arbitrary points only supported for a single element at a time");
497 if (eval_mode == CEED_EVAL_WEIGHT) {
498 CeedCall(CeedVectorSetValue(v, 1.0));
499 return CEED_ERROR_SUCCESS;
500 }
501 if (!basis->basis_chebyshev) {
502 // Build basis mapping from nodes to Chebyshev coefficients
503 CeedScalar *chebyshev_interp_1d, *chebyshev_grad_1d, *chebyshev_q_weight_1d;
504 const CeedScalar *q_ref_1d;
505 Ceed ceed;
506
507 CeedCall(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d));
508 CeedCall(CeedCalloc(P_1d * Q_1d, &chebyshev_grad_1d));
509 CeedCall(CeedCalloc(Q_1d, &chebyshev_q_weight_1d));
510 CeedCall(CeedBasisGetQRef(basis, &q_ref_1d));
511 CeedCall(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d));
512
513 CeedCall(CeedBasisGetCeed(basis, &ceed));
514 CeedCall(CeedVectorCreate(ceed, num_comp * CeedIntPow(Q_1d, dim), &basis->vec_chebyshev));
515 CeedCall(CeedBasisCreateTensorH1(ceed, dim, num_comp, P_1d, Q_1d, chebyshev_interp_1d, chebyshev_grad_1d, q_ref_1d, chebyshev_q_weight_1d,
516 &basis->basis_chebyshev));
517
518 // Cleanup
519 CeedCall(CeedFree(&chebyshev_interp_1d));
520 CeedCall(CeedFree(&chebyshev_grad_1d));
521 CeedCall(CeedFree(&chebyshev_q_weight_1d));
522 CeedCall(CeedDestroy(&ceed));
523 }
524
525 // Create TensorContract object if needed, such as a basis from the GPU backends
526 if (!basis->contract) {
527 Ceed ceed_ref;
528 CeedBasis basis_ref = NULL;
529
530 CeedCall(CeedInit("/cpu/self", &ceed_ref));
531 // Only need matching tensor contraction dimensions, any type of basis will work
532 CeedCall(CeedBasisCreateTensorH1Lagrange(ceed_ref, dim, num_comp, P_1d, Q_1d, CEED_GAUSS, &basis_ref));
533 // Note - clang-tidy doesn't know basis_ref->contract must be valid here
534 CeedCheck(basis_ref && basis_ref->contract, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED,
535 "Reference CPU ceed failed to create a tensor contraction object");
536 CeedCall(CeedTensorContractReferenceCopy(basis_ref->contract, &basis->contract));
537 CeedCall(CeedBasisDestroy(&basis_ref));
538 CeedCall(CeedDestroy(&ceed_ref));
539 }
540
541 // Basis evaluation
542 switch (t_mode) {
543 case CEED_NOTRANSPOSE: {
544 // Nodes to arbitrary points
545 CeedScalar *v_array;
546 const CeedScalar *chebyshev_coeffs, *x_array_read;
547
548 // -- Interpolate to Chebyshev coefficients
549 CeedCall(CeedBasisApply(basis->basis_chebyshev, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, basis->vec_chebyshev));
550
551 // -- Evaluate Chebyshev polynomials at arbitrary points
552 CeedCall(CeedVectorGetArrayRead(basis->vec_chebyshev, CEED_MEM_HOST, &chebyshev_coeffs));
553 CeedCall(CeedVectorGetArrayRead(x_ref, CEED_MEM_HOST, &x_array_read));
554 CeedCall(CeedVectorGetArrayWrite(v, CEED_MEM_HOST, &v_array));
555 switch (eval_mode) {
556 case CEED_EVAL_INTERP: {
557 CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d];
558
559 // ---- Values at point
560 for (CeedInt p = 0; p < total_num_points; p++) {
561 CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1;
562
563 for (CeedInt d = 0; d < dim; d++) {
564 // ------ Tensor contract with current Chebyshev polynomial values
565 CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x));
566 CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false,
567 d == 0 ? chebyshev_coeffs : tmp[d % 2], tmp[(d + 1) % 2]));
568 pre /= Q_1d;
569 post *= 1;
570 }
571 for (CeedInt c = 0; c < num_comp; c++) v_array[c * total_num_points + p] = tmp[dim % 2][c];
572 }
573 break;
574 }
575 case CEED_EVAL_GRAD: {
576 CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d];
577
578 // ---- Values at point
579 for (CeedInt p = 0; p < total_num_points; p++) {
580 // Dim**2 contractions, apply grad when pass == dim
581 for (CeedInt pass = 0; pass < dim; pass++) {
582 CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1;
583
584 for (CeedInt d = 0; d < dim; d++) {
585 // ------ Tensor contract with current Chebyshev polynomial values
586 if (pass == d) CeedCall(CeedChebyshevDerivativeAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x));
587 else CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x));
588 CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false,
589 d == 0 ? chebyshev_coeffs : tmp[d % 2], tmp[(d + 1) % 2]));
590 pre /= Q_1d;
591 post *= 1;
592 }
593 for (CeedInt c = 0; c < num_comp; c++) v_array[(pass * num_comp + c) * total_num_points + p] = tmp[dim % 2][c];
594 }
595 }
596 break;
597 }
598 default:
599 // Nothing to do, excluded above
600 break;
601 }
602 CeedCall(CeedVectorRestoreArrayRead(basis->vec_chebyshev, &chebyshev_coeffs));
603 CeedCall(CeedVectorRestoreArrayRead(x_ref, &x_array_read));
604 CeedCall(CeedVectorRestoreArray(v, &v_array));
605 break;
606 }
607 case CEED_TRANSPOSE: {
608 // Note: No switch on e_mode here because only CEED_EVAL_INTERP is supported at this time
609 // Arbitrary points to nodes
610 CeedScalar *chebyshev_coeffs;
611 const CeedScalar *u_array, *x_array_read;
612
613 // -- Transpose of evaluation of Chebyshev polynomials at arbitrary points
614 CeedCall(CeedVectorGetArrayWrite(basis->vec_chebyshev, CEED_MEM_HOST, &chebyshev_coeffs));
615 CeedCall(CeedVectorGetArrayRead(x_ref, CEED_MEM_HOST, &x_array_read));
616 CeedCall(CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array));
617
618 switch (eval_mode) {
619 case CEED_EVAL_INTERP: {
620 CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d];
621
622 // ---- Values at point
623 for (CeedInt p = 0; p < total_num_points; p++) {
624 CeedInt pre = num_comp * 1, post = 1;
625
626 for (CeedInt c = 0; c < num_comp; c++) tmp[0][c] = u_array[c * total_num_points + p];
627 for (CeedInt d = 0; d < dim; d++) {
628 // ------ Tensor contract with current Chebyshev polynomial values
629 CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x));
630 CeedCall(CeedTensorContractApply(basis->contract, pre, 1, post, Q_1d, chebyshev_x, t_mode, p > 0 && d == (dim - 1), tmp[d % 2],
631 d == (dim - 1) ? chebyshev_coeffs : tmp[(d + 1) % 2]));
632 pre /= 1;
633 post *= Q_1d;
634 }
635 }
636 break;
637 }
638 case CEED_EVAL_GRAD: {
639 CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d];
640
641 // ---- Values at point
642 for (CeedInt p = 0; p < total_num_points; p++) {
643 // Dim**2 contractions, apply grad when pass == dim
644 for (CeedInt pass = 0; pass < dim; pass++) {
645 CeedInt pre = num_comp * 1, post = 1;
646
647 for (CeedInt c = 0; c < num_comp; c++) tmp[0][c] = u_array[(pass * num_comp + c) * total_num_points + p];
648 for (CeedInt d = 0; d < dim; d++) {
649 // ------ Tensor contract with current Chebyshev polynomial values
650 if (pass == d) CeedCall(CeedChebyshevDerivativeAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x));
651 else CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x));
652 CeedCall(CeedTensorContractApply(basis->contract, pre, 1, post, Q_1d, chebyshev_x, t_mode,
653 (p > 0 || (p == 0 && pass > 0)) && d == (dim - 1), tmp[d % 2],
654 d == (dim - 1) ? chebyshev_coeffs : tmp[(d + 1) % 2]));
655 pre /= 1;
656 post *= Q_1d;
657 }
658 }
659 }
660 break;
661 }
662 default:
663 // Nothing to do, excluded above
664 break;
665 }
666 CeedCall(CeedVectorRestoreArray(basis->vec_chebyshev, &chebyshev_coeffs));
667 CeedCall(CeedVectorRestoreArrayRead(x_ref, &x_array_read));
668 CeedCall(CeedVectorRestoreArrayRead(u, &u_array));
669
670 // -- Interpolate transpose from Chebyshev coefficients
671 if (apply_add) CeedCall(CeedBasisApplyAdd(basis->basis_chebyshev, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, basis->vec_chebyshev, v));
672 else CeedCall(CeedBasisApply(basis->basis_chebyshev, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, basis->vec_chebyshev, v));
673 break;
674 }
675 }
676 return CEED_ERROR_SUCCESS;
677 }
678
679 /// @}
680
681 /// ----------------------------------------------------------------------------
682 /// Ceed Backend API
683 /// ----------------------------------------------------------------------------
684 /// @addtogroup CeedBasisBackend
685 /// @{
686
687 /**
688 @brief Fallback to a reference implementation for a non tensor-product basis for \f$H^1\f$ discretizations.
689 This function may only be called inside of a backend `BasisCreateH1` function.
690 This is used by a backend when the specific parameters for a `CeedBasis` exceed the backend's support, such as
691 when a `interp` and `grad` matrices require too many bytes to fit into shared memory on a GPU.
692
693 @param[in] ceed `Ceed` object used to create the `CeedBasis`
694 @param[in] topo Topology of element, e.g. hypercube, simplex, etc
695 @param[in] num_comp Number of field components (1 for scalar fields)
696 @param[in] num_nodes Total number of nodes
697 @param[in] num_qpts Total number of quadrature points
698 @param[in] interp Row-major (`num_qpts * num_nodes`) matrix expressing the values of nodal basis functions at quadrature points
699 @param[in] grad Row-major (`dim * num_qpts * num_nodes`) matrix expressing derivatives of nodal basis functions at quadrature points
700 @param[in] q_ref Array of length `num_qpts * dim` holding the locations of quadrature points on the reference element
701 @param[in] q_weight Array of length `num_qpts` holding the quadrature weights on the reference element
702 @param[out] basis Newly created `CeedBasis`
703
704 @return An error code: 0 - success, otherwise - failure
705
706 @ref User
707 **/
CeedBasisCreateH1Fallback(Ceed ceed,CeedElemTopology topo,CeedInt num_comp,CeedInt num_nodes,CeedInt num_qpts,const CeedScalar * interp,const CeedScalar * grad,const CeedScalar * q_ref,const CeedScalar * q_weight,CeedBasis basis)708 int CeedBasisCreateH1Fallback(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp,
709 const CeedScalar *grad, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
710 CeedInt P = num_nodes, Q = num_qpts, dim = 0;
711 Ceed delegate;
712
713 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis"));
714 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateH1");
715
716 CeedCall(CeedReferenceCopy(delegate, &(basis)->obj.ceed));
717 CeedCall(CeedBasisGetTopologyDimension(topo, &dim));
718 CeedCall(delegate->BasisCreateH1(topo, dim, P, Q, interp, grad, q_ref, q_weight, basis));
719 CeedCall(CeedDestroy(&delegate));
720 return CEED_ERROR_SUCCESS;
721 }
722
723 /**
724 @brief Return collocated gradient matrix
725
726 @param[in] basis `CeedBasis`
727 @param[out] collo_grad_1d Row-major (`Q_1d * Q_1d`) matrix expressing derivatives of basis functions at quadrature points
728
729 @return An error code: 0 - success, otherwise - failure
730
731 @ref Backend
732 **/
CeedBasisGetCollocatedGrad(CeedBasis basis,CeedScalar * collo_grad_1d)733 int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collo_grad_1d) {
734 Ceed ceed;
735 CeedInt P_1d, Q_1d;
736 CeedScalar *interp_1d_pinv;
737 const CeedScalar *grad_1d, *interp_1d;
738
739 // Note: This function is for backend use, so all errors are terminal and we do not need to clean up memory on failure.
740 CeedCall(CeedBasisGetCeed(basis, &ceed));
741 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d));
742 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
743
744 // Compute interp_1d^+, pseudoinverse of interp_1d
745 CeedCall(CeedCalloc(P_1d * Q_1d, &interp_1d_pinv));
746 CeedCall(CeedBasisGetInterp1D(basis, &interp_1d));
747 CeedCall(CeedMatrixPseudoinverse(ceed, interp_1d, Q_1d, P_1d, interp_1d_pinv));
748 CeedCall(CeedBasisGetGrad1D(basis, &grad_1d));
749 CeedCall(CeedMatrixMatrixMultiply(ceed, grad_1d, (const CeedScalar *)interp_1d_pinv, collo_grad_1d, Q_1d, Q_1d, P_1d));
750
751 CeedCall(CeedFree(&interp_1d_pinv));
752 CeedCall(CeedDestroy(&ceed));
753 return CEED_ERROR_SUCCESS;
754 }
755
756 /**
757 @brief Return 1D interpolation matrix to Chebyshev polynomial coefficients on quadrature space
758
759 @param[in] basis `CeedBasis`
760 @param[out] chebyshev_interp_1d Row-major (`P_1d * Q_1d`) matrix interpolating from basis nodes to Chebyshev polynomial coefficients
761
762 @return An error code: 0 - success, otherwise - failure
763
764 @ref Backend
765 **/
CeedBasisGetChebyshevInterp1D(CeedBasis basis,CeedScalar * chebyshev_interp_1d)766 int CeedBasisGetChebyshevInterp1D(CeedBasis basis, CeedScalar *chebyshev_interp_1d) {
767 CeedInt P_1d, Q_1d;
768 CeedScalar *C, *chebyshev_coeffs_1d_inv;
769 const CeedScalar *interp_1d, *q_ref_1d;
770 Ceed ceed;
771
772 CeedCall(CeedBasisGetCeed(basis, &ceed));
773 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d));
774 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
775
776 // Build coefficient matrix
777 // -- Note: Clang-tidy needs this check
778 CeedCheck((P_1d > 0) && (Q_1d > 0), ceed, CEED_ERROR_INCOMPATIBLE, "CeedBasis dimensions are malformed");
779 CeedCall(CeedCalloc(Q_1d * Q_1d, &C));
780 CeedCall(CeedBasisGetQRef(basis, &q_ref_1d));
781 for (CeedInt i = 0; i < Q_1d; i++) CeedCall(CeedChebyshevPolynomialsAtPoint(q_ref_1d[i], Q_1d, &C[i * Q_1d]));
782
783 // Compute C^+, pseudoinverse of coefficient matrix
784 CeedCall(CeedCalloc(Q_1d * Q_1d, &chebyshev_coeffs_1d_inv));
785 CeedCall(CeedMatrixPseudoinverse(ceed, C, Q_1d, Q_1d, chebyshev_coeffs_1d_inv));
786
787 // Build mapping from nodes to Chebyshev coefficients
788 CeedCall(CeedBasisGetInterp1D(basis, &interp_1d));
789 CeedCall(CeedMatrixMatrixMultiply(ceed, chebyshev_coeffs_1d_inv, interp_1d, chebyshev_interp_1d, Q_1d, P_1d, Q_1d));
790
791 // Cleanup
792 CeedCall(CeedFree(&C));
793 CeedCall(CeedFree(&chebyshev_coeffs_1d_inv));
794 CeedCall(CeedDestroy(&ceed));
795 return CEED_ERROR_SUCCESS;
796 }
797
798 /**
799 @brief Get tensor status for given `CeedBasis`
800
801 @param[in] basis `CeedBasis`
802 @param[out] is_tensor Variable to store tensor status
803
804 @return An error code: 0 - success, otherwise - failure
805
806 @ref Backend
807 **/
CeedBasisIsTensor(CeedBasis basis,bool * is_tensor)808 int CeedBasisIsTensor(CeedBasis basis, bool *is_tensor) {
809 *is_tensor = basis->is_tensor_basis;
810 return CEED_ERROR_SUCCESS;
811 }
812
813 /**
814 @brief Determine if given `CeedBasis` has nodes collocated with quadrature points
815
816 @param[in] basis `CeedBasis`
817 @param[out] is_collocated Variable to store collocated status
818
819 @return An error code: 0 - success, otherwise - failure
820
821 @ref Backend
822 **/
CeedBasisIsCollocated(CeedBasis basis,bool * is_collocated)823 int CeedBasisIsCollocated(CeedBasis basis, bool *is_collocated) {
824 if (basis->is_tensor_basis && (basis->Q_1d == basis->P_1d)) {
825 *is_collocated = true;
826
827 for (CeedInt i = 0; i < basis->P_1d; i++) {
828 *is_collocated = *is_collocated && (fabs(basis->interp_1d[i + basis->P_1d * i] - 1.0) < 10 * CEED_EPSILON);
829 for (CeedInt j = 0; j < basis->Q_1d; j++) {
830 if (j != i) *is_collocated = *is_collocated && (fabs(basis->interp_1d[j + basis->P_1d * i]) < 10 * CEED_EPSILON);
831 }
832 }
833 } else {
834 *is_collocated = false;
835 }
836 return CEED_ERROR_SUCCESS;
837 }
838
839 /**
840 @brief Get backend data of a `CeedBasis`
841
842 @param[in] basis `CeedBasis`
843 @param[out] data Variable to store data
844
845 @return An error code: 0 - success, otherwise - failure
846
847 @ref Backend
848 **/
CeedBasisGetData(CeedBasis basis,void * data)849 int CeedBasisGetData(CeedBasis basis, void *data) {
850 *(void **)data = basis->data;
851 return CEED_ERROR_SUCCESS;
852 }
853
854 /**
855 @brief Set backend data of a `CeedBasis`
856
857 @param[in,out] basis `CeedBasis`
858 @param[in] data Data to set
859
860 @return An error code: 0 - success, otherwise - failure
861
862 @ref Backend
863 **/
CeedBasisSetData(CeedBasis basis,void * data)864 int CeedBasisSetData(CeedBasis basis, void *data) {
865 basis->data = data;
866 return CEED_ERROR_SUCCESS;
867 }
868
869 /**
870 @brief Increment the reference counter for a `CeedBasis`
871
872 @param[in,out] basis `CeedBasis` to increment the reference counter
873
874 @return An error code: 0 - success, otherwise - failure
875
876 @ref Backend
877 **/
CeedBasisReference(CeedBasis basis)878 int CeedBasisReference(CeedBasis basis) {
879 CeedCall(CeedObjectReference((CeedObject)basis));
880 return CEED_ERROR_SUCCESS;
881 }
882
883 /**
884 @brief Get number of Q-vector components for given `CeedBasis`
885
886 @param[in] basis `CeedBasis`
887 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values,
888 @ref CEED_EVAL_GRAD to use gradients,
889 @ref CEED_EVAL_DIV to use divergence,
890 @ref CEED_EVAL_CURL to use curl
891 @param[out] q_comp Variable to store number of Q-vector components of basis
892
893 @return An error code: 0 - success, otherwise - failure
894
895 @ref Backend
896 **/
CeedBasisGetNumQuadratureComponents(CeedBasis basis,CeedEvalMode eval_mode,CeedInt * q_comp)897 int CeedBasisGetNumQuadratureComponents(CeedBasis basis, CeedEvalMode eval_mode, CeedInt *q_comp) {
898 CeedInt dim;
899
900 CeedCall(CeedBasisGetDimension(basis, &dim));
901 switch (eval_mode) {
902 case CEED_EVAL_INTERP: {
903 CeedFESpace fe_space;
904
905 CeedCall(CeedBasisGetFESpace(basis, &fe_space));
906 *q_comp = (fe_space == CEED_FE_SPACE_H1) ? 1 : dim;
907 } break;
908 case CEED_EVAL_GRAD:
909 *q_comp = dim;
910 break;
911 case CEED_EVAL_DIV:
912 *q_comp = 1;
913 break;
914 case CEED_EVAL_CURL:
915 *q_comp = (dim < 3) ? 1 : dim;
916 break;
917 case CEED_EVAL_NONE:
918 case CEED_EVAL_WEIGHT:
919 *q_comp = 1;
920 break;
921 }
922 return CEED_ERROR_SUCCESS;
923 }
924
925 /**
926 @brief Estimate number of FLOPs required to apply `CeedBasis` in `t_mode` and `eval_mode`
927
928 @param[in] basis `CeedBasis` to estimate FLOPs for
929 @param[in] t_mode Apply basis or transpose
930 @param[in] eval_mode @ref CeedEvalMode
931 @param[in] is_at_points Evaluate the basis at points or quadrature points
932 @param[in] num_points Number of points basis is evaluated at
933 @param[out] flops Address of variable to hold FLOPs estimate
934
935 @ref Backend
936 **/
CeedBasisGetFlopsEstimate(CeedBasis basis,CeedTransposeMode t_mode,CeedEvalMode eval_mode,bool is_at_points,CeedInt num_points,CeedSize * flops)937 int CeedBasisGetFlopsEstimate(CeedBasis basis, CeedTransposeMode t_mode, CeedEvalMode eval_mode, bool is_at_points, CeedInt num_points,
938 CeedSize *flops) {
939 bool is_tensor;
940
941 CeedCall(CeedBasisIsTensor(basis, &is_tensor));
942 CeedCheck(!is_at_points || is_tensor, CeedBasisReturnCeed(basis), CEED_ERROR_INCOMPATIBLE, "Can only evaluate tensor-product bases at points");
943 if (is_tensor) {
944 CeedInt dim, num_comp, P_1d, Q_1d;
945
946 CeedCall(CeedBasisGetDimension(basis, &dim));
947 CeedCall(CeedBasisGetNumComponents(basis, &num_comp));
948 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d));
949 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
950 if (t_mode == CEED_TRANSPOSE) {
951 P_1d = Q_1d;
952 Q_1d = P_1d;
953 }
954 CeedInt tensor_flops = 0, pre = num_comp * CeedIntPow(P_1d, dim - 1), post = 1;
955
956 for (CeedInt d = 0; d < dim; d++) {
957 tensor_flops += 2 * pre * P_1d * post * Q_1d;
958 pre /= P_1d;
959 post *= Q_1d;
960 }
961 if (is_at_points) {
962 bool is_gpu = false;
963
964 {
965 CeedMemType mem_type;
966
967 CeedCall(CeedGetPreferredMemType(CeedBasisReturnCeed(basis), &mem_type));
968 is_gpu = mem_type == CEED_MEM_DEVICE;
969 }
970
971 CeedInt chebyshev_flops = (Q_1d - 2) * 3 + 1, d_chebyshev_flops = (Q_1d - 2) * 8 + 1;
972 CeedInt point_tensor_flops = 0, pre = CeedIntPow(Q_1d, dim - 1), post = 1;
973
974 for (CeedInt d = 0; d < dim; d++) {
975 point_tensor_flops += 2 * pre * Q_1d * post * 1;
976 pre /= P_1d;
977 post *= Q_1d;
978 }
979
980 switch (eval_mode) {
981 case CEED_EVAL_NONE:
982 *flops = 0;
983 break;
984 case CEED_EVAL_INTERP: {
985 *flops = tensor_flops + num_points * num_comp * (point_tensor_flops + (t_mode == CEED_TRANSPOSE ? CeedIntPow(Q_1d, dim) : 0));
986 if (dim == 3 && is_gpu) {
987 *flops += num_points * Q_1d *
988 (chebyshev_flops + num_comp * (2 * chebyshev_flops + 2 * Q_1d * Q_1d + (t_mode == CEED_TRANSPOSE ? 2 * Q_1d + 1 : 3 * Q_1d)));
989 } else {
990 *flops += num_points * (is_gpu ? num_comp : 1) * dim * chebyshev_flops;
991 }
992 break;
993 }
994 case CEED_EVAL_GRAD: {
995 *flops = tensor_flops + num_points * num_comp * (point_tensor_flops + (t_mode == CEED_TRANSPOSE ? CeedIntPow(Q_1d, dim) : 0));
996 if (dim == 3 && is_gpu) {
997 CeedInt inner_flops =
998 dim * (2 * Q_1d * Q_1d + (t_mode == CEED_TRANSPOSE ? 2 : 3) * Q_1d) + (dim - 1) * (2 * chebyshev_flops + d_chebyshev_flops);
999
1000 *flops += num_points * Q_1d * (chebyshev_flops + d_chebyshev_flops + num_comp * (inner_flops + (t_mode == CEED_TRANSPOSE ? 1 : 0)));
1001 } else {
1002 *flops += num_points * (is_gpu ? num_comp : 1) * dim * (d_chebyshev_flops + (dim - 1) * chebyshev_flops);
1003 }
1004 break;
1005 }
1006 case CEED_EVAL_DIV:
1007 case CEED_EVAL_CURL: {
1008 // LCOV_EXCL_START
1009 return CeedError(CeedBasisReturnCeed(basis), CEED_ERROR_INCOMPATIBLE, "Tensor basis evaluation for %s not supported at points",
1010 CeedEvalModes[eval_mode]);
1011 break;
1012 // LCOV_EXCL_STOP
1013 }
1014 case CEED_EVAL_WEIGHT:
1015 *flops = num_points;
1016 break;
1017 }
1018 } else {
1019 switch (eval_mode) {
1020 case CEED_EVAL_NONE:
1021 *flops = 0;
1022 break;
1023 case CEED_EVAL_INTERP:
1024 *flops = tensor_flops;
1025 break;
1026 case CEED_EVAL_GRAD:
1027 *flops = tensor_flops * 2;
1028 break;
1029 case CEED_EVAL_DIV:
1030 case CEED_EVAL_CURL: {
1031 // LCOV_EXCL_START
1032 return CeedError(CeedBasisReturnCeed(basis), CEED_ERROR_INCOMPATIBLE, "Tensor basis evaluation for %s not supported",
1033 CeedEvalModes[eval_mode]);
1034 break;
1035 // LCOV_EXCL_STOP
1036 }
1037 case CEED_EVAL_WEIGHT:
1038 *flops = dim * CeedIntPow(Q_1d, dim);
1039 break;
1040 }
1041 }
1042 } else {
1043 CeedInt dim, num_comp, q_comp, num_nodes, num_qpts;
1044
1045 CeedCall(CeedBasisGetDimension(basis, &dim));
1046 CeedCall(CeedBasisGetNumComponents(basis, &num_comp));
1047 CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp));
1048 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes));
1049 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
1050 switch (eval_mode) {
1051 case CEED_EVAL_NONE:
1052 *flops = 0;
1053 break;
1054 case CEED_EVAL_INTERP:
1055 case CEED_EVAL_GRAD:
1056 case CEED_EVAL_DIV:
1057 case CEED_EVAL_CURL:
1058 *flops = num_nodes * num_qpts * num_comp * q_comp;
1059 break;
1060 case CEED_EVAL_WEIGHT:
1061 *flops = 0;
1062 break;
1063 }
1064 }
1065 return CEED_ERROR_SUCCESS;
1066 }
1067
1068 /**
1069 @brief Get `CeedFESpace` for a `CeedBasis`
1070
1071 @param[in] basis `CeedBasis`
1072 @param[out] fe_space Variable to store `CeedFESpace`
1073
1074 @return An error code: 0 - success, otherwise - failure
1075
1076 @ref Backend
1077 **/
CeedBasisGetFESpace(CeedBasis basis,CeedFESpace * fe_space)1078 int CeedBasisGetFESpace(CeedBasis basis, CeedFESpace *fe_space) {
1079 *fe_space = basis->fe_space;
1080 return CEED_ERROR_SUCCESS;
1081 }
1082
1083 /**
1084 @brief Get dimension for given `CeedElemTopology`
1085
1086 @param[in] topo `CeedElemTopology`
1087 @param[out] dim Variable to store dimension of topology
1088
1089 @return An error code: 0 - success, otherwise - failure
1090
1091 @ref Backend
1092 **/
CeedBasisGetTopologyDimension(CeedElemTopology topo,CeedInt * dim)1093 int CeedBasisGetTopologyDimension(CeedElemTopology topo, CeedInt *dim) {
1094 *dim = (CeedInt)topo >> 16;
1095 return CEED_ERROR_SUCCESS;
1096 }
1097
1098 /**
1099 @brief Get `CeedTensorContract` of a `CeedBasis`
1100
1101 @param[in] basis `CeedBasis`
1102 @param[out] contract Variable to store `CeedTensorContract`
1103
1104 @return An error code: 0 - success, otherwise - failure
1105
1106 @ref Backend
1107 **/
CeedBasisGetTensorContract(CeedBasis basis,CeedTensorContract * contract)1108 int CeedBasisGetTensorContract(CeedBasis basis, CeedTensorContract *contract) {
1109 *contract = basis->contract;
1110 return CEED_ERROR_SUCCESS;
1111 }
1112
1113 /**
1114 @brief Set `CeedTensorContract` of a `CeedBasis`
1115
1116 @param[in,out] basis `CeedBasis`
1117 @param[in] contract `CeedTensorContract` to set
1118
1119 @return An error code: 0 - success, otherwise - failure
1120
1121 @ref Backend
1122 **/
CeedBasisSetTensorContract(CeedBasis basis,CeedTensorContract contract)1123 int CeedBasisSetTensorContract(CeedBasis basis, CeedTensorContract contract) {
1124 basis->contract = contract;
1125 CeedCall(CeedTensorContractReference(contract));
1126 return CEED_ERROR_SUCCESS;
1127 }
1128
1129 /**
1130 @brief Return a reference implementation of matrix multiplication \f$C = A B\f$.
1131
1132 Note: This is a reference implementation for CPU `CeedScalar` pointers that is not intended for high performance.
1133
1134 @param[in] ceed `Ceed` context for error handling
1135 @param[in] mat_A Row-major matrix `A`
1136 @param[in] mat_B Row-major matrix `B`
1137 @param[out] mat_C Row-major output matrix `C`
1138 @param[in] m Number of rows of `C`
1139 @param[in] n Number of columns of `C`
1140 @param[in] kk Number of columns of `A`/rows of `B`
1141
1142 @return An error code: 0 - success, otherwise - failure
1143
1144 @ref Utility
1145 **/
CeedMatrixMatrixMultiply(Ceed ceed,const CeedScalar * mat_A,const CeedScalar * mat_B,CeedScalar * mat_C,CeedInt m,CeedInt n,CeedInt kk)1146 int CeedMatrixMatrixMultiply(Ceed ceed, const CeedScalar *mat_A, const CeedScalar *mat_B, CeedScalar *mat_C, CeedInt m, CeedInt n, CeedInt kk) {
1147 for (CeedInt i = 0; i < m; i++) {
1148 for (CeedInt j = 0; j < n; j++) {
1149 CeedScalar sum = 0;
1150
1151 for (CeedInt k = 0; k < kk; k++) sum += mat_A[k + i * kk] * mat_B[j + k * n];
1152 mat_C[j + i * n] = sum;
1153 }
1154 }
1155 return CEED_ERROR_SUCCESS;
1156 }
1157
1158 /**
1159 @brief Return QR Factorization of a matrix
1160
1161 @param[in] ceed `Ceed` context for error handling
1162 @param[in,out] mat Row-major matrix to be factorized in place
1163 @param[in,out] tau Vector of length `m` of scaling factors
1164 @param[in] m Number of rows
1165 @param[in] n Number of columns
1166
1167 @return An error code: 0 - success, otherwise - failure
1168
1169 @ref Utility
1170 **/
CeedQRFactorization(Ceed ceed,CeedScalar * mat,CeedScalar * tau,CeedInt m,CeedInt n)1171 int CeedQRFactorization(Ceed ceed, CeedScalar *mat, CeedScalar *tau, CeedInt m, CeedInt n) {
1172 CeedScalar v[m];
1173
1174 // Check matrix shape
1175 CeedCheck(n <= m, ceed, CEED_ERROR_UNSUPPORTED, "Cannot compute QR factorization with n > m");
1176
1177 for (CeedInt i = 0; i < n; i++) {
1178 CeedScalar sigma = 0.0;
1179
1180 if (i >= m - 1) { // last row of matrix, no reflection needed
1181 tau[i] = 0.;
1182 break;
1183 }
1184 // Calculate Householder vector, magnitude
1185 v[i] = mat[i + n * i];
1186 for (CeedInt j = i + 1; j < m; j++) {
1187 v[j] = mat[i + n * j];
1188 sigma += v[j] * v[j];
1189 }
1190 const CeedScalar norm = sqrt(v[i] * v[i] + sigma); // norm of v[i:m]
1191 const CeedScalar R_ii = -copysign(norm, v[i]);
1192
1193 v[i] -= R_ii;
1194 // norm of v[i:m] after modification above and scaling below
1195 // norm = sqrt(v[i]*v[i] + sigma) / v[i];
1196 // tau = 2 / (norm*norm)
1197 tau[i] = 2 * v[i] * v[i] / (v[i] * v[i] + sigma);
1198 for (CeedInt j = i + 1; j < m; j++) v[j] /= v[i];
1199
1200 // Apply Householder reflector to lower right panel
1201 CeedHouseholderReflect(&mat[i * n + i + 1], &v[i], tau[i], m - i, n - i - 1, n, 1);
1202 // Save v
1203 mat[i + n * i] = R_ii;
1204 for (CeedInt j = i + 1; j < m; j++) mat[i + n * j] = v[j];
1205 }
1206 return CEED_ERROR_SUCCESS;
1207 }
1208
1209 /**
1210 @brief Apply Householder Q matrix
1211
1212 Compute `mat_A = mat_Q mat_A`, where `mat_Q` is \f$m \times m\f$ and `mat_A` is \f$m \times n\f$.
1213
1214 @param[in,out] mat_A Matrix to apply Householder Q to, in place
1215 @param[in] mat_Q Householder Q matrix
1216 @param[in] tau Householder scaling factors
1217 @param[in] t_mode Transpose mode for application
1218 @param[in] m Number of rows in `A`
1219 @param[in] n Number of columns in `A`
1220 @param[in] k Number of elementary reflectors in Q, `k < m`
1221 @param[in] row Row stride in `A`
1222 @param[in] col Col stride in `A`
1223
1224 @return An error code: 0 - success, otherwise - failure
1225
1226 @ref Utility
1227 **/
CeedHouseholderApplyQ(CeedScalar * mat_A,const CeedScalar * mat_Q,const CeedScalar * tau,CeedTransposeMode t_mode,CeedInt m,CeedInt n,CeedInt k,CeedInt row,CeedInt col)1228 int CeedHouseholderApplyQ(CeedScalar *mat_A, const CeedScalar *mat_Q, const CeedScalar *tau, CeedTransposeMode t_mode, CeedInt m, CeedInt n,
1229 CeedInt k, CeedInt row, CeedInt col) {
1230 CeedScalar *v;
1231
1232 CeedCall(CeedMalloc(m, &v));
1233 for (CeedInt ii = 0; ii < k; ii++) {
1234 CeedInt i = t_mode == CEED_TRANSPOSE ? ii : k - 1 - ii;
1235 for (CeedInt j = i + 1; j < m; j++) v[j] = mat_Q[j * k + i];
1236 // Apply Householder reflector (I - tau v v^T) collo_grad_1d^T
1237 CeedCall(CeedHouseholderReflect(&mat_A[i * row], &v[i], tau[i], m - i, n, row, col));
1238 }
1239 CeedCall(CeedFree(&v));
1240 return CEED_ERROR_SUCCESS;
1241 }
1242
1243 /**
1244 @brief Return pseudoinverse of a matrix
1245
1246 @param[in] ceed Ceed context for error handling
1247 @param[in] mat Row-major matrix to compute pseudoinverse of
1248 @param[in] m Number of rows
1249 @param[in] n Number of columns
1250 @param[out] mat_pinv Row-major pseudoinverse matrix
1251
1252 @return An error code: 0 - success, otherwise - failure
1253
1254 @ref Utility
1255 **/
CeedMatrixPseudoinverse(Ceed ceed,const CeedScalar * mat,CeedInt m,CeedInt n,CeedScalar * mat_pinv)1256 int CeedMatrixPseudoinverse(Ceed ceed, const CeedScalar *mat, CeedInt m, CeedInt n, CeedScalar *mat_pinv) {
1257 CeedScalar *tau, *I, *mat_copy;
1258
1259 CeedCall(CeedCalloc(m, &tau));
1260 CeedCall(CeedCalloc(m * m, &I));
1261 CeedCall(CeedCalloc(m * n, &mat_copy));
1262 memcpy(mat_copy, mat, m * n * sizeof mat[0]);
1263
1264 // QR Factorization, mat = Q R
1265 CeedCall(CeedQRFactorization(ceed, mat_copy, tau, m, n));
1266
1267 // -- Apply Q^T, I = Q^T * I
1268 for (CeedInt i = 0; i < m; i++) I[i * m + i] = 1.0;
1269 CeedCall(CeedHouseholderApplyQ(I, mat_copy, tau, CEED_TRANSPOSE, m, m, n, m, 1));
1270 // -- Apply R_inv, mat_pinv = R_inv * Q^T
1271 for (CeedInt j = 0; j < m; j++) { // Column j
1272 mat_pinv[j + m * (n - 1)] = I[j + m * (n - 1)] / mat_copy[n * n - 1];
1273 for (CeedInt i = n - 2; i >= 0; i--) { // Row i
1274 mat_pinv[j + m * i] = I[j + m * i];
1275 for (CeedInt k = i + 1; k < n; k++) mat_pinv[j + m * i] -= mat_copy[k + n * i] * mat_pinv[j + m * k];
1276 mat_pinv[j + m * i] /= mat_copy[i + n * i];
1277 }
1278 }
1279
1280 // Cleanup
1281 CeedCall(CeedFree(&I));
1282 CeedCall(CeedFree(&tau));
1283 CeedCall(CeedFree(&mat_copy));
1284 return CEED_ERROR_SUCCESS;
1285 }
1286
1287 /**
1288 @brief Return symmetric Schur decomposition of the symmetric matrix mat via symmetric QR factorization
1289
1290 @param[in] ceed `Ceed` context for error handling
1291 @param[in,out] mat Row-major matrix to be factorized in place
1292 @param[out] lambda Vector of length n of eigenvalues
1293 @param[in] n Number of rows/columns
1294
1295 @return An error code: 0 - success, otherwise - failure
1296
1297 @ref Utility
1298 **/
1299 CeedPragmaOptimizeOff
CeedSymmetricSchurDecomposition(Ceed ceed,CeedScalar * mat,CeedScalar * lambda,CeedInt n)1300 int CeedSymmetricSchurDecomposition(Ceed ceed, CeedScalar *mat, CeedScalar *lambda, CeedInt n) {
1301 // Check bounds for clang-tidy
1302 CeedCheck(n > 1, ceed, CEED_ERROR_UNSUPPORTED, "Cannot compute symmetric Schur decomposition of scalars");
1303
1304 CeedScalar v[n - 1], tau[n - 1], mat_T[n * n];
1305
1306 // Copy mat to mat_T and set mat to I
1307 memcpy(mat_T, mat, n * n * sizeof(mat[0]));
1308 for (CeedInt i = 0; i < n; i++) {
1309 for (CeedInt j = 0; j < n; j++) mat[j + n * i] = (i == j) ? 1 : 0;
1310 }
1311
1312 // Reduce to tridiagonal
1313 for (CeedInt i = 0; i < n - 1; i++) {
1314 // Calculate Householder vector, magnitude
1315 CeedScalar sigma = 0.0;
1316
1317 v[i] = mat_T[i + n * (i + 1)];
1318 for (CeedInt j = i + 1; j < n - 1; j++) {
1319 v[j] = mat_T[i + n * (j + 1)];
1320 sigma += v[j] * v[j];
1321 }
1322 const CeedScalar norm = sqrt(v[i] * v[i] + sigma); // norm of v[i:n-1]
1323 const CeedScalar R_ii = -copysign(norm, v[i]);
1324
1325 v[i] -= R_ii;
1326 // norm of v[i:m] after modification above and scaling below
1327 // norm = sqrt(v[i]*v[i] + sigma) / v[i];
1328 // tau = 2 / (norm*norm)
1329 tau[i] = i == n - 2 ? 2 : 2 * v[i] * v[i] / (v[i] * v[i] + sigma);
1330 for (CeedInt j = i + 1; j < n - 1; j++) v[j] /= v[i];
1331
1332 // Update sub and super diagonal
1333 for (CeedInt j = i + 2; j < n; j++) {
1334 mat_T[i + n * j] = 0;
1335 mat_T[j + n * i] = 0;
1336 }
1337 // Apply symmetric Householder reflector to lower right panel
1338 CeedHouseholderReflect(&mat_T[(i + 1) + n * (i + 1)], &v[i], tau[i], n - (i + 1), n - (i + 1), n, 1);
1339 CeedHouseholderReflect(&mat_T[(i + 1) + n * (i + 1)], &v[i], tau[i], n - (i + 1), n - (i + 1), 1, n);
1340
1341 // Save v
1342 mat_T[i + n * (i + 1)] = R_ii;
1343 mat_T[(i + 1) + n * i] = R_ii;
1344 for (CeedInt j = i + 1; j < n - 1; j++) {
1345 mat_T[i + n * (j + 1)] = v[j];
1346 }
1347 }
1348 // Backwards accumulation of Q
1349 for (CeedInt i = n - 2; i >= 0; i--) {
1350 if (tau[i] > 0.0) {
1351 v[i] = 1;
1352 for (CeedInt j = i + 1; j < n - 1; j++) {
1353 v[j] = mat_T[i + n * (j + 1)];
1354 mat_T[i + n * (j + 1)] = 0;
1355 }
1356 CeedHouseholderReflect(&mat[(i + 1) + n * (i + 1)], &v[i], tau[i], n - (i + 1), n - (i + 1), n, 1);
1357 }
1358 }
1359
1360 // Reduce sub and super diagonal
1361 CeedInt p = 0, q = 0, itr = 0, max_itr = n * n * n * n;
1362 CeedScalar tol = CEED_EPSILON;
1363
1364 while (itr < max_itr) {
1365 // Update p, q, size of reduced portions of diagonal
1366 p = 0;
1367 q = 0;
1368 for (CeedInt i = n - 2; i >= 0; i--) {
1369 if (fabs(mat_T[i + n * (i + 1)]) < tol) q += 1;
1370 else break;
1371 }
1372 for (CeedInt i = 0; i < n - q - 1; i++) {
1373 if (fabs(mat_T[i + n * (i + 1)]) < tol) p += 1;
1374 else break;
1375 }
1376 if (q == n - 1) break; // Finished reducing
1377
1378 // Reduce tridiagonal portion
1379 CeedScalar t_nn = mat_T[(n - 1 - q) + n * (n - 1 - q)], t_nnm1 = mat_T[(n - 2 - q) + n * (n - 1 - q)];
1380 CeedScalar d = (mat_T[(n - 2 - q) + n * (n - 2 - q)] - t_nn) / 2;
1381 CeedScalar mu = t_nn - t_nnm1 * t_nnm1 / (d + copysign(sqrt(d * d + t_nnm1 * t_nnm1), d));
1382 CeedScalar x = mat_T[p + n * p] - mu;
1383 CeedScalar z = mat_T[p + n * (p + 1)];
1384
1385 for (CeedInt k = p; k < n - q - 1; k++) {
1386 // Compute Givens rotation
1387 CeedScalar c = 1, s = 0;
1388
1389 if (fabs(z) > tol) {
1390 if (fabs(z) > fabs(x)) {
1391 const CeedScalar tau = -x / z;
1392
1393 s = 1 / sqrt(1 + tau * tau);
1394 c = s * tau;
1395 } else {
1396 const CeedScalar tau = -z / x;
1397
1398 c = 1 / sqrt(1 + tau * tau);
1399 s = c * tau;
1400 }
1401 }
1402
1403 // Apply Givens rotation to T
1404 CeedGivensRotation(mat_T, c, s, CEED_NOTRANSPOSE, k, k + 1, n, n);
1405 CeedGivensRotation(mat_T, c, s, CEED_TRANSPOSE, k, k + 1, n, n);
1406
1407 // Apply Givens rotation to Q
1408 CeedGivensRotation(mat, c, s, CEED_NOTRANSPOSE, k, k + 1, n, n);
1409
1410 // Update x, z
1411 if (k < n - q - 2) {
1412 x = mat_T[k + n * (k + 1)];
1413 z = mat_T[k + n * (k + 2)];
1414 }
1415 }
1416 itr++;
1417 }
1418
1419 // Save eigenvalues
1420 for (CeedInt i = 0; i < n; i++) lambda[i] = mat_T[i + n * i];
1421
1422 // Check convergence
1423 CeedCheck(itr < max_itr || q > n, ceed, CEED_ERROR_MINOR, "Symmetric QR failed to converge");
1424 return CEED_ERROR_SUCCESS;
1425 }
1426 CeedPragmaOptimizeOn
1427
1428 /**
1429 @brief Return Simultaneous Diagonalization of two matrices.
1430
1431 This solves the generalized eigenvalue problem `A x = lambda B x`, where `A` and `B` are symmetric and `B` is positive definite.
1432 We generate the matrix `X` and vector `Lambda` such that `X^T A X = Lambda` and `X^T B X = I`.
1433 This is equivalent to the LAPACK routine 'sygv' with `TYPE = 1`.
1434
1435 @param[in] ceed `Ceed` context for error handling
1436 @param[in] mat_A Row-major matrix to be factorized with eigenvalues
1437 @param[in] mat_B Row-major matrix to be factorized to identity
1438 @param[out] mat_X Row-major orthogonal matrix
1439 @param[out] lambda Vector of length `n` of generalized eigenvalues
1440 @param[in] n Number of rows/columns
1441
1442 @return An error code: 0 - success, otherwise - failure
1443
1444 @ref Utility
1445 **/
1446 CeedPragmaOptimizeOff
CeedSimultaneousDiagonalization(Ceed ceed,CeedScalar * mat_A,CeedScalar * mat_B,CeedScalar * mat_X,CeedScalar * lambda,CeedInt n)1447 int CeedSimultaneousDiagonalization(Ceed ceed, CeedScalar *mat_A, CeedScalar *mat_B, CeedScalar *mat_X, CeedScalar *lambda, CeedInt n) {
1448 CeedScalar *mat_C, *mat_G, *vec_D;
1449
1450 CeedCall(CeedCalloc(n * n, &mat_C));
1451 CeedCall(CeedCalloc(n * n, &mat_G));
1452 CeedCall(CeedCalloc(n, &vec_D));
1453
1454 // Compute B = G D G^T
1455 memcpy(mat_G, mat_B, n * n * sizeof(mat_B[0]));
1456 CeedCall(CeedSymmetricSchurDecomposition(ceed, mat_G, vec_D, n));
1457
1458 // Sort eigenvalues
1459 for (CeedInt i = n - 1; i >= 0; i--) {
1460 for (CeedInt j = 0; j < i; j++) {
1461 if (fabs(vec_D[j]) > fabs(vec_D[j + 1])) {
1462 CeedScalarSwap(vec_D[j], vec_D[j + 1]);
1463 for (CeedInt k = 0; k < n; k++) CeedScalarSwap(mat_G[k * n + j], mat_G[k * n + j + 1]);
1464 }
1465 }
1466 }
1467
1468 // Compute C = (G D^1/2)^-1 A (G D^1/2)^-T
1469 // = D^-1/2 G^T A G D^-1/2
1470 // -- D = D^-1/2
1471 for (CeedInt i = 0; i < n; i++) vec_D[i] = 1. / sqrt(vec_D[i]);
1472 // -- G = G D^-1/2
1473 // -- C = D^-1/2 G^T
1474 for (CeedInt i = 0; i < n; i++) {
1475 for (CeedInt j = 0; j < n; j++) {
1476 mat_G[i * n + j] *= vec_D[j];
1477 mat_C[j * n + i] = mat_G[i * n + j];
1478 }
1479 }
1480 // -- X = (D^-1/2 G^T) A
1481 CeedCall(CeedMatrixMatrixMultiply(ceed, (const CeedScalar *)mat_C, (const CeedScalar *)mat_A, mat_X, n, n, n));
1482 // -- C = (D^-1/2 G^T A) (G D^-1/2)
1483 CeedCall(CeedMatrixMatrixMultiply(ceed, (const CeedScalar *)mat_X, (const CeedScalar *)mat_G, mat_C, n, n, n));
1484
1485 // Compute Q^T C Q = lambda
1486 CeedCall(CeedSymmetricSchurDecomposition(ceed, mat_C, lambda, n));
1487
1488 // Sort eigenvalues
1489 for (CeedInt i = n - 1; i >= 0; i--) {
1490 for (CeedInt j = 0; j < i; j++) {
1491 if (fabs(lambda[j]) > fabs(lambda[j + 1])) {
1492 CeedScalarSwap(lambda[j], lambda[j + 1]);
1493 for (CeedInt k = 0; k < n; k++) CeedScalarSwap(mat_C[k * n + j], mat_C[k * n + j + 1]);
1494 }
1495 }
1496 }
1497
1498 // Set X = (G D^1/2)^-T Q
1499 // = G D^-1/2 Q
1500 CeedCall(CeedMatrixMatrixMultiply(ceed, (const CeedScalar *)mat_G, (const CeedScalar *)mat_C, mat_X, n, n, n));
1501
1502 // Cleanup
1503 CeedCall(CeedFree(&mat_C));
1504 CeedCall(CeedFree(&mat_G));
1505 CeedCall(CeedFree(&vec_D));
1506 return CEED_ERROR_SUCCESS;
1507 }
1508 CeedPragmaOptimizeOn
1509
1510 /// @}
1511
1512 /// ----------------------------------------------------------------------------
1513 /// CeedBasis Public API
1514 /// ----------------------------------------------------------------------------
1515 /// @addtogroup CeedBasisUser
1516 /// @{
1517
1518 /**
1519 @brief Create a tensor-product basis for \f$H^1\f$ discretizations
1520
1521 @param[in] ceed `Ceed` object used to create the `CeedBasis`
1522 @param[in] dim Topological dimension
1523 @param[in] num_comp Number of field components (1 for scalar fields)
1524 @param[in] P_1d Number of nodes in one dimension
1525 @param[in] Q_1d Number of quadrature points in one dimension
1526 @param[in] interp_1d Row-major (`Q_1d * P_1d`) matrix expressing the values of nodal basis functions at quadrature points
1527 @param[in] grad_1d Row-major (`Q_1d * P_1d`) matrix expressing derivatives of nodal basis functions at quadrature points
1528 @param[in] q_ref_1d Array of length `Q_1d` holding the locations of quadrature points on the 1D reference element `[-1, 1]`
1529 @param[in] q_weight_1d Array of length `Q_1d` holding the quadrature weights on the reference element
1530 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored
1531
1532 @return An error code: 0 - success, otherwise - failure
1533
1534 @ref User
1535 **/
CeedBasisCreateTensorH1(Ceed ceed,CeedInt dim,CeedInt num_comp,CeedInt P_1d,CeedInt Q_1d,const CeedScalar * interp_1d,const CeedScalar * grad_1d,const CeedScalar * q_ref_1d,const CeedScalar * q_weight_1d,CeedBasis * basis)1536 int CeedBasisCreateTensorH1(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d,
1537 const CeedScalar *grad_1d, const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis *basis) {
1538 if (!ceed->BasisCreateTensorH1) {
1539 Ceed delegate;
1540
1541 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis"));
1542 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateTensorH1");
1543 CeedCall(CeedBasisCreateTensorH1(delegate, dim, num_comp, P_1d, Q_1d, interp_1d, grad_1d, q_ref_1d, q_weight_1d, basis));
1544 CeedCall(CeedDestroy(&delegate));
1545 return CEED_ERROR_SUCCESS;
1546 }
1547
1548 CeedCheck(dim > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis dimension must be a positive value");
1549 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component");
1550 CeedCheck(P_1d > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node");
1551 CeedCheck(Q_1d > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point");
1552
1553 CeedElemTopology topo = dim == 1 ? CEED_TOPOLOGY_LINE : dim == 2 ? CEED_TOPOLOGY_QUAD : CEED_TOPOLOGY_HEX;
1554
1555 CeedCall(CeedCalloc(1, basis));
1556 CeedCall(CeedObjectCreate(ceed, CeedBasisView_Object, CeedBasisDestroy_Object, &(*basis)->obj));
1557 (*basis)->is_tensor_basis = true;
1558 (*basis)->dim = dim;
1559 (*basis)->topo = topo;
1560 (*basis)->num_comp = num_comp;
1561 (*basis)->P_1d = P_1d;
1562 (*basis)->Q_1d = Q_1d;
1563 (*basis)->P = CeedIntPow(P_1d, dim);
1564 (*basis)->Q = CeedIntPow(Q_1d, dim);
1565 (*basis)->fe_space = CEED_FE_SPACE_H1;
1566 CeedCall(CeedCalloc(Q_1d, &(*basis)->q_ref_1d));
1567 CeedCall(CeedCalloc(Q_1d, &(*basis)->q_weight_1d));
1568 if (q_ref_1d) memcpy((*basis)->q_ref_1d, q_ref_1d, Q_1d * sizeof(q_ref_1d[0]));
1569 if (q_weight_1d) memcpy((*basis)->q_weight_1d, q_weight_1d, Q_1d * sizeof(q_weight_1d[0]));
1570 CeedCall(CeedCalloc(Q_1d * P_1d, &(*basis)->interp_1d));
1571 CeedCall(CeedCalloc(Q_1d * P_1d, &(*basis)->grad_1d));
1572 if (interp_1d) memcpy((*basis)->interp_1d, interp_1d, Q_1d * P_1d * sizeof(interp_1d[0]));
1573 if (grad_1d) memcpy((*basis)->grad_1d, grad_1d, Q_1d * P_1d * sizeof(grad_1d[0]));
1574 CeedCall(ceed->BasisCreateTensorH1(dim, P_1d, Q_1d, interp_1d, grad_1d, q_ref_1d, q_weight_1d, *basis));
1575 return CEED_ERROR_SUCCESS;
1576 }
1577
1578 /**
1579 @brief Create a tensor-product \f$H^1\f$ Lagrange basis
1580
1581 @param[in] ceed `Ceed` object used to create the `CeedBasis`
1582 @param[in] dim Topological dimension of element
1583 @param[in] num_comp Number of field components (1 for scalar fields)
1584 @param[in] P Number of Gauss-Lobatto nodes in one dimension.
1585 The polynomial degree of the resulting `Q_k` element is `k = P - 1`.
1586 @param[in] Q Number of quadrature points in one dimension.
1587 @param[in] quad_mode Distribution of the `Q` quadrature points (affects order of accuracy for the quadrature)
1588 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored
1589
1590 @return An error code: 0 - success, otherwise - failure
1591
1592 @ref User
1593 **/
CeedBasisCreateTensorH1Lagrange(Ceed ceed,CeedInt dim,CeedInt num_comp,CeedInt P,CeedInt Q,CeedQuadMode quad_mode,CeedBasis * basis)1594 int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P, CeedInt Q, CeedQuadMode quad_mode, CeedBasis *basis) {
1595 // Allocate
1596 int ierr = CEED_ERROR_SUCCESS;
1597 CeedScalar c1, c2, c3, c4, dx, *nodes, *interp_1d, *grad_1d, *q_ref_1d, *q_weight_1d;
1598
1599 CeedCheck(dim > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis dimension must be a positive value");
1600 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component");
1601 CeedCheck(P > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node");
1602 CeedCheck(Q > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point");
1603
1604 // Get Nodes and Weights
1605 CeedCall(CeedCalloc(P * Q, &interp_1d));
1606 CeedCall(CeedCalloc(P * Q, &grad_1d));
1607 CeedCall(CeedCalloc(P, &nodes));
1608 CeedCall(CeedCalloc(Q, &q_ref_1d));
1609 CeedCall(CeedCalloc(Q, &q_weight_1d));
1610 if (CeedLobattoQuadrature(P, nodes, NULL) != CEED_ERROR_SUCCESS) goto cleanup;
1611 switch (quad_mode) {
1612 case CEED_GAUSS:
1613 ierr = CeedGaussQuadrature(Q, q_ref_1d, q_weight_1d);
1614 break;
1615 case CEED_GAUSS_LOBATTO:
1616 ierr = CeedLobattoQuadrature(Q, q_ref_1d, q_weight_1d);
1617 break;
1618 }
1619 if (ierr != CEED_ERROR_SUCCESS) goto cleanup;
1620
1621 // Build B, D matrix
1622 // Fornberg, 1998
1623 for (CeedInt i = 0; i < Q; i++) {
1624 c1 = 1.0;
1625 c3 = nodes[0] - q_ref_1d[i];
1626 interp_1d[i * P + 0] = 1.0;
1627 for (CeedInt j = 1; j < P; j++) {
1628 c2 = 1.0;
1629 c4 = c3;
1630 c3 = nodes[j] - q_ref_1d[i];
1631 for (CeedInt k = 0; k < j; k++) {
1632 dx = nodes[j] - nodes[k];
1633 c2 *= dx;
1634 if (k == j - 1) {
1635 grad_1d[i * P + j] = c1 * (interp_1d[i * P + k] - c4 * grad_1d[i * P + k]) / c2;
1636 interp_1d[i * P + j] = -c1 * c4 * interp_1d[i * P + k] / c2;
1637 }
1638 grad_1d[i * P + k] = (c3 * grad_1d[i * P + k] - interp_1d[i * P + k]) / dx;
1639 interp_1d[i * P + k] = c3 * interp_1d[i * P + k] / dx;
1640 }
1641 c1 = c2;
1642 }
1643 }
1644 // Pass to CeedBasisCreateTensorH1
1645 CeedCall(CeedBasisCreateTensorH1(ceed, dim, num_comp, P, Q, interp_1d, grad_1d, q_ref_1d, q_weight_1d, basis));
1646 cleanup:
1647 CeedCall(CeedFree(&interp_1d));
1648 CeedCall(CeedFree(&grad_1d));
1649 CeedCall(CeedFree(&nodes));
1650 CeedCall(CeedFree(&q_ref_1d));
1651 CeedCall(CeedFree(&q_weight_1d));
1652 return CEED_ERROR_SUCCESS;
1653 }
1654
1655 /**
1656 @brief Create a non tensor-product basis for \f$H^1\f$ discretizations
1657
1658 @param[in] ceed `Ceed` object used to create the `CeedBasis`
1659 @param[in] topo Topology of element, e.g. hypercube, simplex, etc
1660 @param[in] num_comp Number of field components (1 for scalar fields)
1661 @param[in] num_nodes Total number of nodes
1662 @param[in] num_qpts Total number of quadrature points
1663 @param[in] interp Row-major (`num_qpts * num_nodes`) matrix expressing the values of nodal basis functions at quadrature points
1664 @param[in] grad Row-major (`dim * num_qpts * num_nodes`) matrix expressing derivatives of nodal basis functions at quadrature points
1665 @param[in] q_ref Array of length `num_qpts * dim` holding the locations of quadrature points on the reference element
1666 @param[in] q_weight Array of length `num_qpts` holding the quadrature weights on the reference element
1667 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored
1668
1669 @return An error code: 0 - success, otherwise - failure
1670
1671 @ref User
1672 **/
CeedBasisCreateH1(Ceed ceed,CeedElemTopology topo,CeedInt num_comp,CeedInt num_nodes,CeedInt num_qpts,const CeedScalar * interp,const CeedScalar * grad,const CeedScalar * q_ref,const CeedScalar * q_weight,CeedBasis * basis)1673 int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp,
1674 const CeedScalar *grad, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis *basis) {
1675 CeedInt P = num_nodes, Q = num_qpts, dim = 0;
1676
1677 if (!ceed->BasisCreateH1) {
1678 Ceed delegate;
1679
1680 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis"));
1681 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateH1");
1682 CeedCall(CeedBasisCreateH1(delegate, topo, num_comp, num_nodes, num_qpts, interp, grad, q_ref, q_weight, basis));
1683 CeedCall(CeedDestroy(&delegate));
1684 return CEED_ERROR_SUCCESS;
1685 }
1686
1687 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component");
1688 CeedCheck(num_nodes > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node");
1689 CeedCheck(num_qpts > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point");
1690
1691 CeedCall(CeedBasisGetTopologyDimension(topo, &dim));
1692
1693 CeedCall(CeedCalloc(1, basis));
1694 CeedCall(CeedObjectCreate(ceed, CeedBasisView_Object, CeedBasisDestroy_Object, &(*basis)->obj));
1695 (*basis)->is_tensor_basis = false;
1696 (*basis)->dim = dim;
1697 (*basis)->topo = topo;
1698 (*basis)->num_comp = num_comp;
1699 (*basis)->P = P;
1700 (*basis)->Q = Q;
1701 (*basis)->fe_space = CEED_FE_SPACE_H1;
1702 CeedCall(CeedCalloc(Q * dim, &(*basis)->q_ref_1d));
1703 CeedCall(CeedCalloc(Q, &(*basis)->q_weight_1d));
1704 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0]));
1705 if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0]));
1706 CeedCall(CeedCalloc(Q * P, &(*basis)->interp));
1707 CeedCall(CeedCalloc(dim * Q * P, &(*basis)->grad));
1708 if (interp) memcpy((*basis)->interp, interp, Q * P * sizeof(interp[0]));
1709 if (grad) memcpy((*basis)->grad, grad, dim * Q * P * sizeof(grad[0]));
1710 CeedCall(ceed->BasisCreateH1(topo, dim, P, Q, interp, grad, q_ref, q_weight, *basis));
1711 return CEED_ERROR_SUCCESS;
1712 }
1713
1714 /**
1715 @brief Create a non tensor-product basis for \f$H(\mathrm{div})\f$ discretizations
1716
1717 @param[in] ceed `Ceed` object used to create the `CeedBasis`
1718 @param[in] topo Topology of element (`CEED_TOPOLOGY_QUAD`, `CEED_TOPOLOGY_PRISM`, etc.), dimension of which is used in some array sizes below
1719 @param[in] num_comp Number of components (usually 1 for vectors in H(div) bases)
1720 @param[in] num_nodes Total number of nodes (DoFs per element)
1721 @param[in] num_qpts Total number of quadrature points
1722 @param[in] interp Row-major (`dim * num_qpts * num_nodes`) matrix expressing the values of basis functions at quadrature points
1723 @param[in] div Row-major (`num_qpts * num_nodes`) matrix expressing divergence of basis functions at quadrature points
1724 @param[in] q_ref Array of length `num_qpts` * dim holding the locations of quadrature points on the reference element
1725 @param[in] q_weight Array of length `num_qpts` holding the quadrature weights on the reference element
1726 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored
1727
1728 @return An error code: 0 - success, otherwise - failure
1729
1730 @ref User
1731 **/
CeedBasisCreateHdiv(Ceed ceed,CeedElemTopology topo,CeedInt num_comp,CeedInt num_nodes,CeedInt num_qpts,const CeedScalar * interp,const CeedScalar * div,const CeedScalar * q_ref,const CeedScalar * q_weight,CeedBasis * basis)1732 int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp,
1733 const CeedScalar *div, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis *basis) {
1734 CeedInt Q = num_qpts, P = num_nodes, dim = 0;
1735
1736 if (!ceed->BasisCreateHdiv) {
1737 Ceed delegate;
1738
1739 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis"));
1740 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateHdiv");
1741 CeedCall(CeedBasisCreateHdiv(delegate, topo, num_comp, num_nodes, num_qpts, interp, div, q_ref, q_weight, basis));
1742 CeedCall(CeedDestroy(&delegate));
1743 return CEED_ERROR_SUCCESS;
1744 }
1745
1746 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component");
1747 CeedCheck(num_nodes > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node");
1748 CeedCheck(num_qpts > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point");
1749
1750 CeedCall(CeedBasisGetTopologyDimension(topo, &dim));
1751
1752 CeedCall(CeedCalloc(1, basis));
1753 CeedCall(CeedObjectCreate(ceed, CeedBasisView_Object, CeedBasisDestroy_Object, &(*basis)->obj));
1754 (*basis)->is_tensor_basis = false;
1755 (*basis)->dim = dim;
1756 (*basis)->topo = topo;
1757 (*basis)->num_comp = num_comp;
1758 (*basis)->P = P;
1759 (*basis)->Q = Q;
1760 (*basis)->fe_space = CEED_FE_SPACE_HDIV;
1761 CeedCall(CeedMalloc(Q * dim, &(*basis)->q_ref_1d));
1762 CeedCall(CeedMalloc(Q, &(*basis)->q_weight_1d));
1763 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0]));
1764 if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0]));
1765 CeedCall(CeedMalloc(dim * Q * P, &(*basis)->interp));
1766 CeedCall(CeedMalloc(Q * P, &(*basis)->div));
1767 if (interp) memcpy((*basis)->interp, interp, dim * Q * P * sizeof(interp[0]));
1768 if (div) memcpy((*basis)->div, div, Q * P * sizeof(div[0]));
1769 CeedCall(ceed->BasisCreateHdiv(topo, dim, P, Q, interp, div, q_ref, q_weight, *basis));
1770 return CEED_ERROR_SUCCESS;
1771 }
1772
1773 /**
1774 @brief Create a non tensor-product basis for \f$H(\mathrm{curl})\f$ discretizations
1775
1776 @param[in] ceed `Ceed` object used to create the `CeedBasis`
1777 @param[in] topo Topology of element (`CEED_TOPOLOGY_QUAD`, `CEED_TOPOLOGY_PRISM`, etc.), dimension of which is used in some array sizes below
1778 @param[in] num_comp Number of components (usually 1 for vectors in \f$H(\mathrm{curl})\f$ bases)
1779 @param[in] num_nodes Total number of nodes (DoFs per element)
1780 @param[in] num_qpts Total number of quadrature points
1781 @param[in] interp Row-major (`dim * num_qpts * num_nodes`) matrix expressing the values of basis functions at quadrature points
1782 @param[in] curl Row-major (`curl_comp * num_qpts * num_nodes`, `curl_comp = 1` if `dim < 3` otherwise `curl_comp = dim`) matrix expressing curl of basis functions at quadrature points
1783 @param[in] q_ref Array of length `num_qpts * dim` holding the locations of quadrature points on the reference element
1784 @param[in] q_weight Array of length `num_qpts` holding the quadrature weights on the reference element
1785 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored
1786
1787 @return An error code: 0 - success, otherwise - failure
1788
1789 @ref User
1790 **/
CeedBasisCreateHcurl(Ceed ceed,CeedElemTopology topo,CeedInt num_comp,CeedInt num_nodes,CeedInt num_qpts,const CeedScalar * interp,const CeedScalar * curl,const CeedScalar * q_ref,const CeedScalar * q_weight,CeedBasis * basis)1791 int CeedBasisCreateHcurl(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp,
1792 const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis *basis) {
1793 CeedInt Q = num_qpts, P = num_nodes, dim = 0, curl_comp = 0;
1794
1795 if (!ceed->BasisCreateHcurl) {
1796 Ceed delegate;
1797
1798 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis"));
1799 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateHcurl");
1800 CeedCall(CeedBasisCreateHcurl(delegate, topo, num_comp, num_nodes, num_qpts, interp, curl, q_ref, q_weight, basis));
1801 CeedCall(CeedDestroy(&delegate));
1802 return CEED_ERROR_SUCCESS;
1803 }
1804
1805 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component");
1806 CeedCheck(num_nodes > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node");
1807 CeedCheck(num_qpts > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point");
1808
1809 CeedCall(CeedBasisGetTopologyDimension(topo, &dim));
1810 curl_comp = (dim < 3) ? 1 : dim;
1811
1812 CeedCall(CeedCalloc(1, basis));
1813 CeedCall(CeedObjectCreate(ceed, CeedBasisView_Object, CeedBasisDestroy_Object, &(*basis)->obj));
1814 (*basis)->is_tensor_basis = false;
1815 (*basis)->dim = dim;
1816 (*basis)->topo = topo;
1817 (*basis)->num_comp = num_comp;
1818 (*basis)->P = P;
1819 (*basis)->Q = Q;
1820 (*basis)->fe_space = CEED_FE_SPACE_HCURL;
1821 CeedCall(CeedMalloc(Q * dim, &(*basis)->q_ref_1d));
1822 CeedCall(CeedMalloc(Q, &(*basis)->q_weight_1d));
1823 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0]));
1824 if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0]));
1825 CeedCall(CeedMalloc(dim * Q * P, &(*basis)->interp));
1826 CeedCall(CeedMalloc(curl_comp * Q * P, &(*basis)->curl));
1827 if (interp) memcpy((*basis)->interp, interp, dim * Q * P * sizeof(interp[0]));
1828 if (curl) memcpy((*basis)->curl, curl, curl_comp * Q * P * sizeof(curl[0]));
1829 CeedCall(ceed->BasisCreateHcurl(topo, dim, P, Q, interp, curl, q_ref, q_weight, *basis));
1830 return CEED_ERROR_SUCCESS;
1831 }
1832
1833 /**
1834 @brief Create a `CeedBasis` for projection from the nodes of `basis_from` to the nodes of `basis_to`.
1835
1836 Only @ref CEED_EVAL_INTERP will be valid for the new basis, `basis_project`.
1837 For \f$H^1\f$ spaces, @ref CEED_EVAL_GRAD will also be valid.
1838 The interpolation is given by `interp_project = interp_to^+ * interp_from`, where the pseudoinverse `interp_to^+` is given by QR factorization.
1839 The gradient (for the \f$H^1\f$ case) is given by `grad_project = interp_to^+ * grad_from`.
1840
1841 Note: `basis_from` and `basis_to` must have compatible quadrature spaces.
1842
1843 Note: `basis_project` will have the same number of components as `basis_from`, regardless of the number of components that `basis_to` has.
1844 If `basis_from` has 3 components and `basis_to` has 5 components, then `basis_project` will have 3 components.
1845
1846 Note: If either `basis_from` or `basis_to` are non-tensor, then `basis_project` will also be non-tensor
1847
1848 @param[in] basis_from `CeedBasis` to prolong from
1849 @param[in] basis_to `CeedBasis` to prolong to
1850 @param[out] basis_project Address of the variable where the newly created `CeedBasis` will be stored
1851
1852 @return An error code: 0 - success, otherwise - failure
1853
1854 @ref User
1855 **/
CeedBasisCreateProjection(CeedBasis basis_from,CeedBasis basis_to,CeedBasis * basis_project)1856 int CeedBasisCreateProjection(CeedBasis basis_from, CeedBasis basis_to, CeedBasis *basis_project) {
1857 Ceed ceed;
1858 bool create_tensor;
1859 CeedInt dim, num_comp;
1860 CeedScalar *interp_project, *grad_project;
1861
1862 CeedCall(CeedBasisGetCeed(basis_to, &ceed));
1863
1864 // Create projection matrix
1865 CeedCall(CeedBasisCreateProjectionMatrices(basis_from, basis_to, &interp_project, &grad_project));
1866
1867 // Build basis
1868 {
1869 bool is_tensor_to, is_tensor_from;
1870
1871 CeedCall(CeedBasisIsTensor(basis_to, &is_tensor_to));
1872 CeedCall(CeedBasisIsTensor(basis_from, &is_tensor_from));
1873 create_tensor = is_tensor_from && is_tensor_to;
1874 }
1875 CeedCall(CeedBasisGetDimension(basis_to, &dim));
1876 CeedCall(CeedBasisGetNumComponents(basis_from, &num_comp));
1877 if (create_tensor) {
1878 CeedInt P_1d_to, P_1d_from;
1879
1880 CeedCall(CeedBasisGetNumNodes1D(basis_from, &P_1d_from));
1881 CeedCall(CeedBasisGetNumNodes1D(basis_to, &P_1d_to));
1882 CeedCall(CeedBasisCreateTensorH1(ceed, dim, num_comp, P_1d_from, P_1d_to, interp_project, grad_project, NULL, NULL, basis_project));
1883 } else {
1884 // Even if basis_to and basis_from are not H1, the resulting basis is H1 for interpolation to work
1885 CeedInt num_nodes_to, num_nodes_from;
1886 CeedElemTopology topo;
1887
1888 CeedCall(CeedBasisGetTopology(basis_from, &topo));
1889 CeedCall(CeedBasisGetNumNodes(basis_from, &num_nodes_from));
1890 CeedCall(CeedBasisGetNumNodes(basis_to, &num_nodes_to));
1891 CeedCall(CeedBasisCreateH1(ceed, topo, num_comp, num_nodes_from, num_nodes_to, interp_project, grad_project, NULL, NULL, basis_project));
1892 }
1893
1894 // Cleanup
1895 CeedCall(CeedFree(&interp_project));
1896 CeedCall(CeedFree(&grad_project));
1897 CeedCall(CeedDestroy(&ceed));
1898 return CEED_ERROR_SUCCESS;
1899 }
1900
1901 /**
1902 @brief Copy the pointer to a `CeedBasis`.
1903
1904 Note: If the value of `*basis_copy` passed into this function is non-`NULL`, then it is assumed that `*basis_copy` is a pointer to a `CeedBasis`.
1905 This `CeedBasis` will be destroyed if `*basis_copy` is the only reference to this `CeedBasis`.
1906
1907 @param[in] basis `CeedBasis` to copy reference to
1908 @param[in,out] basis_copy Variable to store copied reference
1909
1910 @return An error code: 0 - success, otherwise - failure
1911
1912 @ref User
1913 **/
CeedBasisReferenceCopy(CeedBasis basis,CeedBasis * basis_copy)1914 int CeedBasisReferenceCopy(CeedBasis basis, CeedBasis *basis_copy) {
1915 if (basis != CEED_BASIS_NONE) CeedCall(CeedBasisReference(basis));
1916 CeedCall(CeedBasisDestroy(basis_copy));
1917 *basis_copy = basis;
1918 return CEED_ERROR_SUCCESS;
1919 }
1920
1921 /**
1922 @brief Set the number of tabs to indent for @ref CeedBasisView() output
1923
1924 @param[in] basis `CeedBasis` to set the number of view tabs
1925 @param[in] num_tabs Number of view tabs to set
1926
1927 @return Error code: 0 - success, otherwise - failure
1928
1929 @ref User
1930 **/
CeedBasisSetNumViewTabs(CeedBasis basis,CeedInt num_tabs)1931 int CeedBasisSetNumViewTabs(CeedBasis basis, CeedInt num_tabs) {
1932 CeedCall(CeedObjectSetNumViewTabs((CeedObject)basis, num_tabs));
1933 return CEED_ERROR_SUCCESS;
1934 }
1935
1936 /**
1937 @brief Get the number of tabs to indent for @ref CeedBasisView() output
1938
1939 @param[in] basis `CeedBasis` to get the number of view tabs
1940 @param[out] num_tabs Number of view tabs
1941
1942 @return Error code: 0 - success, otherwise - failure
1943
1944 @ref User
1945 **/
CeedBasisGetNumViewTabs(CeedBasis basis,CeedInt * num_tabs)1946 int CeedBasisGetNumViewTabs(CeedBasis basis, CeedInt *num_tabs) {
1947 CeedCall(CeedObjectGetNumViewTabs((CeedObject)basis, num_tabs));
1948 return CEED_ERROR_SUCCESS;
1949 }
1950
1951 /**
1952 @brief View a `CeedBasis`
1953
1954 @param[in] basis `CeedBasis` to view
1955 @param[in] stream Stream to view to, e.g., `stdout`
1956
1957 @return An error code: 0 - success, otherwise - failure
1958
1959 @ref User
1960 **/
CeedBasisView(CeedBasis basis,FILE * stream)1961 int CeedBasisView(CeedBasis basis, FILE *stream) {
1962 bool is_tensor_basis;
1963 char *tabs = NULL;
1964 CeedElemTopology topo;
1965 CeedFESpace fe_space;
1966
1967 // Basis data
1968 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis));
1969 CeedCall(CeedBasisGetTopology(basis, &topo));
1970 CeedCall(CeedBasisGetFESpace(basis, &fe_space));
1971
1972 {
1973 CeedInt num_tabs = 0;
1974
1975 CeedCall(CeedBasisGetNumViewTabs(basis, &num_tabs));
1976 CeedCall(CeedCalloc(CEED_TAB_WIDTH * num_tabs + 1, &tabs));
1977 for (CeedInt i = 0; i < CEED_TAB_WIDTH * num_tabs; i++) tabs[i] = ' ';
1978 }
1979
1980 // Print FE space and element topology of the basis
1981 fprintf(stream, "%sCeedBasis in a %s on a %s element\n", tabs, CeedFESpaces[fe_space], CeedElemTopologies[topo]);
1982 if (is_tensor_basis) {
1983 fprintf(stream, "%s P: %" CeedInt_FMT "\n%s Q: %" CeedInt_FMT "\n", tabs, basis->P_1d, tabs, basis->Q_1d);
1984 } else {
1985 fprintf(stream, "%s P: %" CeedInt_FMT "\n%s Q: %" CeedInt_FMT "\n", tabs, basis->P, tabs, basis->Q);
1986 }
1987 fprintf(stream, "%s dimension: %" CeedInt_FMT "\n%s field components: %" CeedInt_FMT "\n", tabs, basis->dim, tabs, basis->num_comp);
1988 // Print quadrature data, interpolation/gradient/divergence/curl of the basis
1989 if (is_tensor_basis) { // tensor basis
1990 CeedInt P_1d, Q_1d;
1991 const CeedScalar *q_ref_1d, *q_weight_1d, *interp_1d, *grad_1d;
1992
1993 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d));
1994 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
1995 CeedCall(CeedBasisGetQRef(basis, &q_ref_1d));
1996 CeedCall(CeedBasisGetQWeights(basis, &q_weight_1d));
1997 CeedCall(CeedBasisGetInterp1D(basis, &interp_1d));
1998 CeedCall(CeedBasisGetGrad1D(basis, &grad_1d));
1999
2000 CeedCall(CeedScalarView("qref1d", "\t% 12.8f", 1, Q_1d, q_ref_1d, tabs, stream));
2001 CeedCall(CeedScalarView("qweight1d", "\t% 12.8f", 1, Q_1d, q_weight_1d, tabs, stream));
2002 CeedCall(CeedScalarView("interp1d", "\t% 12.8f", Q_1d, P_1d, interp_1d, tabs, stream));
2003 CeedCall(CeedScalarView("grad1d", "\t% 12.8f", Q_1d, P_1d, grad_1d, tabs, stream));
2004 } else { // non-tensor basis
2005 CeedInt P, Q, dim, q_comp;
2006 const CeedScalar *q_ref, *q_weight, *interp, *grad, *div, *curl;
2007
2008 CeedCall(CeedBasisGetNumNodes(basis, &P));
2009 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &Q));
2010 CeedCall(CeedBasisGetDimension(basis, &dim));
2011 CeedCall(CeedBasisGetQRef(basis, &q_ref));
2012 CeedCall(CeedBasisGetQWeights(basis, &q_weight));
2013 CeedCall(CeedBasisGetInterp(basis, &interp));
2014 CeedCall(CeedBasisGetGrad(basis, &grad));
2015 CeedCall(CeedBasisGetDiv(basis, &div));
2016 CeedCall(CeedBasisGetCurl(basis, &curl));
2017
2018 CeedCall(CeedScalarView("qref", "\t% 12.8f", 1, Q * dim, q_ref, tabs, stream));
2019 CeedCall(CeedScalarView("qweight", "\t% 12.8f", 1, Q, q_weight, tabs, stream));
2020 CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp));
2021 CeedCall(CeedScalarView("interp", "\t% 12.8f", q_comp * Q, P, interp, tabs, stream));
2022 if (grad) {
2023 CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp));
2024 CeedCall(CeedScalarView("grad", "\t% 12.8f", q_comp * Q, P, grad, tabs, stream));
2025 }
2026 if (div) {
2027 CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp));
2028 CeedCall(CeedScalarView("div", "\t% 12.8f", q_comp * Q, P, div, tabs, stream));
2029 }
2030 if (curl) {
2031 CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp));
2032 CeedCall(CeedScalarView("curl", "\t% 12.8f", q_comp * Q, P, curl, tabs, stream));
2033 }
2034 }
2035 CeedCall(CeedFree(&tabs));
2036 return CEED_ERROR_SUCCESS;
2037 }
2038
2039 /**
2040 @brief Apply basis evaluation from nodes to quadrature points or vice versa
2041
2042 @param[in] basis `CeedBasis` to evaluate
2043 @param[in] num_elem The number of elements to apply the basis evaluation to;
2044 the backend will specify the ordering in @ref CeedElemRestrictionCreate()
2045 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to quadrature points;
2046 @ref CEED_TRANSPOSE to apply the transpose, mapping from quadrature points to nodes
2047 @param[in] eval_mode @ref CEED_EVAL_NONE to use values directly,
2048 @ref CEED_EVAL_INTERP to use interpolated values,
2049 @ref CEED_EVAL_GRAD to use gradients,
2050 @ref CEED_EVAL_DIV to use divergence,
2051 @ref CEED_EVAL_CURL to use curl,
2052 @ref CEED_EVAL_WEIGHT to use quadrature weights
2053 @param[in] u Input `CeedVector`
2054 @param[out] v Output `CeedVector`
2055
2056 @return An error code: 0 - success, otherwise - failure
2057
2058 @ref User
2059 **/
CeedBasisApply(CeedBasis basis,CeedInt num_elem,CeedTransposeMode t_mode,CeedEvalMode eval_mode,CeedVector u,CeedVector v)2060 int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) {
2061 CeedCall(CeedBasisApplyCheckDims(basis, num_elem, t_mode, eval_mode, u, v));
2062 CeedCheck(basis->Apply, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Backend does not support CeedBasisApply");
2063 CeedCall(basis->Apply(basis, num_elem, t_mode, eval_mode, u, v));
2064 return CEED_ERROR_SUCCESS;
2065 }
2066
2067 /**
2068 @brief Apply basis evaluation from quadrature points to nodes and sum into target vector
2069
2070 @param[in] basis `CeedBasis` to evaluate
2071 @param[in] num_elem The number of elements to apply the basis evaluation to;
2072 the backend will specify the ordering in @ref CeedElemRestrictionCreate()
2073 @param[in] t_mode @ref CEED_TRANSPOSE to apply the transpose, mapping from quadrature points to nodes;
2074 @ref CEED_NOTRANSPOSE is not valid for `CeedBasisApplyAdd()`
2075 @param[in] eval_mode @ref CEED_EVAL_NONE to use values directly,
2076 @ref CEED_EVAL_INTERP to use interpolated values,
2077 @ref CEED_EVAL_GRAD to use gradients,
2078 @ref CEED_EVAL_DIV to use divergence,
2079 @ref CEED_EVAL_CURL to use curl,
2080 @ref CEED_EVAL_WEIGHT to use quadrature weights
2081 @param[in] u Input `CeedVector`
2082 @param[out] v Output `CeedVector` to sum into
2083
2084 @return An error code: 0 - success, otherwise - failure
2085
2086 @ref User
2087 **/
CeedBasisApplyAdd(CeedBasis basis,CeedInt num_elem,CeedTransposeMode t_mode,CeedEvalMode eval_mode,CeedVector u,CeedVector v)2088 int CeedBasisApplyAdd(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) {
2089 CeedCheck(t_mode == CEED_TRANSPOSE, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "CeedBasisApplyAdd only supports CEED_TRANSPOSE");
2090 CeedCall(CeedBasisApplyCheckDims(basis, num_elem, t_mode, eval_mode, u, v));
2091 CeedCheck(basis->ApplyAdd, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedBasisApplyAdd");
2092 CeedCall(basis->ApplyAdd(basis, num_elem, t_mode, eval_mode, u, v));
2093 return CEED_ERROR_SUCCESS;
2094 }
2095
2096 /**
2097 @brief Apply basis evaluation from nodes to arbitrary points
2098
2099 @param[in] basis `CeedBasis` to evaluate
2100 @param[in] num_elem The number of elements to apply the basis evaluation to;
2101 the backend will specify the ordering in @ref CeedElemRestrictionCreate()
2102 @param[in] num_points Array of the number of points to apply the basis evaluation to in each element, size `num_elem`
2103 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points;
2104 @ref CEED_TRANSPOSE to apply the transpose, mapping from points to nodes
2105 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values,
2106 @ref CEED_EVAL_GRAD to use gradients,
2107 @ref CEED_EVAL_WEIGHT to use quadrature weights
2108 @param[in] x_ref `CeedVector` holding reference coordinates of each point
2109 @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE
2110 @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP
2111
2112 @return An error code: 0 - success, otherwise - failure
2113
2114 @ref User
2115 **/
CeedBasisApplyAtPoints(CeedBasis basis,CeedInt num_elem,const CeedInt * num_points,CeedTransposeMode t_mode,CeedEvalMode eval_mode,CeedVector x_ref,CeedVector u,CeedVector v)2116 int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode,
2117 CeedVector x_ref, CeedVector u, CeedVector v) {
2118 CeedCall(CeedBasisApplyAtPointsCheckDims(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v));
2119 if (basis->ApplyAtPoints) {
2120 CeedCall(basis->ApplyAtPoints(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v));
2121 } else {
2122 CeedCall(CeedBasisApplyAtPoints_Core(basis, false, num_elem, num_points, t_mode, eval_mode, x_ref, u, v));
2123 }
2124 return CEED_ERROR_SUCCESS;
2125 }
2126
2127 /**
2128 @brief Apply basis evaluation from nodes to arbitrary points and sum into target vector
2129
2130 @param[in] basis `CeedBasis` to evaluate
2131 @param[in] num_elem The number of elements to apply the basis evaluation to;
2132 the backend will specify the ordering in @ref CeedElemRestrictionCreate()
2133 @param[in] num_points Array of the number of points to apply the basis evaluation to in each element, size `num_elem`
2134 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points;
2135 @ref CEED_NOTRANSPOSE is not valid for `CeedBasisApplyAddAtPoints()`
2136 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values,
2137 @ref CEED_EVAL_GRAD to use gradients,
2138 @ref CEED_EVAL_WEIGHT to use quadrature weights
2139 @param[in] x_ref `CeedVector` holding reference coordinates of each point
2140 @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE
2141 @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP
2142
2143 @return An error code: 0 - success, otherwise - failure
2144
2145 @ref User
2146 **/
CeedBasisApplyAddAtPoints(CeedBasis basis,CeedInt num_elem,const CeedInt * num_points,CeedTransposeMode t_mode,CeedEvalMode eval_mode,CeedVector x_ref,CeedVector u,CeedVector v)2147 int CeedBasisApplyAddAtPoints(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode,
2148 CeedVector x_ref, CeedVector u, CeedVector v) {
2149 CeedCheck(t_mode == CEED_TRANSPOSE, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "CeedBasisApplyAddAtPoints only supports CEED_TRANSPOSE");
2150 CeedCall(CeedBasisApplyAtPointsCheckDims(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v));
2151 if (basis->ApplyAddAtPoints) {
2152 CeedCall(basis->ApplyAddAtPoints(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v));
2153 } else {
2154 CeedCall(CeedBasisApplyAtPoints_Core(basis, true, num_elem, num_points, t_mode, eval_mode, x_ref, u, v));
2155 }
2156 return CEED_ERROR_SUCCESS;
2157 }
2158
2159 /**
2160 @brief Get the `Ceed` associated with a `CeedBasis`
2161
2162 @param[in] basis `CeedBasis`
2163 @param[out] ceed Variable to store `Ceed`
2164
2165 @return An error code: 0 - success, otherwise - failure
2166
2167 @ref Advanced
2168 **/
CeedBasisGetCeed(CeedBasis basis,Ceed * ceed)2169 int CeedBasisGetCeed(CeedBasis basis, Ceed *ceed) {
2170 CeedCall(CeedObjectGetCeed((CeedObject)basis, ceed));
2171 return CEED_ERROR_SUCCESS;
2172 }
2173
2174 /**
2175 @brief Return the `Ceed` associated with a `CeedBasis`
2176
2177 @param[in] basis `CeedBasis`
2178
2179 @return `Ceed` associated with the `basis`
2180
2181 @ref Advanced
2182 **/
CeedBasisReturnCeed(CeedBasis basis)2183 Ceed CeedBasisReturnCeed(CeedBasis basis) { return CeedObjectReturnCeed((CeedObject)basis); }
2184
2185 /**
2186 @brief Get dimension for given `CeedBasis`
2187
2188 @param[in] basis `CeedBasis`
2189 @param[out] dim Variable to store dimension of basis
2190
2191 @return An error code: 0 - success, otherwise - failure
2192
2193 @ref Advanced
2194 **/
CeedBasisGetDimension(CeedBasis basis,CeedInt * dim)2195 int CeedBasisGetDimension(CeedBasis basis, CeedInt *dim) {
2196 *dim = basis->dim;
2197 return CEED_ERROR_SUCCESS;
2198 }
2199
2200 /**
2201 @brief Get topology for given `CeedBasis`
2202
2203 @param[in] basis `CeedBasis`
2204 @param[out] topo Variable to store topology of basis
2205
2206 @return An error code: 0 - success, otherwise - failure
2207
2208 @ref Advanced
2209 **/
CeedBasisGetTopology(CeedBasis basis,CeedElemTopology * topo)2210 int CeedBasisGetTopology(CeedBasis basis, CeedElemTopology *topo) {
2211 *topo = basis->topo;
2212 return CEED_ERROR_SUCCESS;
2213 }
2214
2215 /**
2216 @brief Get number of components for given `CeedBasis`
2217
2218 @param[in] basis `CeedBasis`
2219 @param[out] num_comp Variable to store number of components
2220
2221 @return An error code: 0 - success, otherwise - failure
2222
2223 @ref Advanced
2224 **/
CeedBasisGetNumComponents(CeedBasis basis,CeedInt * num_comp)2225 int CeedBasisGetNumComponents(CeedBasis basis, CeedInt *num_comp) {
2226 *num_comp = basis->num_comp;
2227 return CEED_ERROR_SUCCESS;
2228 }
2229
2230 /**
2231 @brief Get total number of nodes (in `dim` dimensions) of a `CeedBasis`
2232
2233 @param[in] basis `CeedBasis`
2234 @param[out] P Variable to store number of nodes
2235
2236 @return An error code: 0 - success, otherwise - failure
2237
2238 @ref Utility
2239 **/
CeedBasisGetNumNodes(CeedBasis basis,CeedInt * P)2240 int CeedBasisGetNumNodes(CeedBasis basis, CeedInt *P) {
2241 *P = basis->P;
2242 return CEED_ERROR_SUCCESS;
2243 }
2244
2245 /**
2246 @brief Get total number of nodes (in 1 dimension) of a `CeedBasis`
2247
2248 @param[in] basis `CeedBasis`
2249 @param[out] P_1d Variable to store number of nodes
2250
2251 @return An error code: 0 - success, otherwise - failure
2252
2253 @ref Advanced
2254 **/
CeedBasisGetNumNodes1D(CeedBasis basis,CeedInt * P_1d)2255 int CeedBasisGetNumNodes1D(CeedBasis basis, CeedInt *P_1d) {
2256 CeedCheck(basis->is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "Cannot supply P_1d for non-tensor CeedBasis");
2257 *P_1d = basis->P_1d;
2258 return CEED_ERROR_SUCCESS;
2259 }
2260
2261 /**
2262 @brief Get total number of quadrature points (in `dim` dimensions) of a `CeedBasis`
2263
2264 @param[in] basis `CeedBasis`
2265 @param[out] Q Variable to store number of quadrature points
2266
2267 @return An error code: 0 - success, otherwise - failure
2268
2269 @ref Utility
2270 **/
CeedBasisGetNumQuadraturePoints(CeedBasis basis,CeedInt * Q)2271 int CeedBasisGetNumQuadraturePoints(CeedBasis basis, CeedInt *Q) {
2272 *Q = basis->Q;
2273 return CEED_ERROR_SUCCESS;
2274 }
2275
2276 /**
2277 @brief Get total number of quadrature points (in 1 dimension) of a `CeedBasis`
2278
2279 @param[in] basis `CeedBasis`
2280 @param[out] Q_1d Variable to store number of quadrature points
2281
2282 @return An error code: 0 - success, otherwise - failure
2283
2284 @ref Advanced
2285 **/
CeedBasisGetNumQuadraturePoints1D(CeedBasis basis,CeedInt * Q_1d)2286 int CeedBasisGetNumQuadraturePoints1D(CeedBasis basis, CeedInt *Q_1d) {
2287 CeedCheck(basis->is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "Cannot supply Q_1d for non-tensor CeedBasis");
2288 *Q_1d = basis->Q_1d;
2289 return CEED_ERROR_SUCCESS;
2290 }
2291
2292 /**
2293 @brief Get reference coordinates of quadrature points (in `dim` dimensions) of a `CeedBasis`
2294
2295 @param[in] basis `CeedBasis`
2296 @param[out] q_ref Variable to store reference coordinates of quadrature points
2297
2298 @return An error code: 0 - success, otherwise - failure
2299
2300 @ref Advanced
2301 **/
CeedBasisGetQRef(CeedBasis basis,const CeedScalar ** q_ref)2302 int CeedBasisGetQRef(CeedBasis basis, const CeedScalar **q_ref) {
2303 *q_ref = basis->q_ref_1d;
2304 return CEED_ERROR_SUCCESS;
2305 }
2306
2307 /**
2308 @brief Get quadrature weights of quadrature points (in `dim` dimensions) of a `CeedBasis`
2309
2310 @param[in] basis `CeedBasis`
2311 @param[out] q_weight Variable to store quadrature weights
2312
2313 @return An error code: 0 - success, otherwise - failure
2314
2315 @ref Advanced
2316 **/
CeedBasisGetQWeights(CeedBasis basis,const CeedScalar ** q_weight)2317 int CeedBasisGetQWeights(CeedBasis basis, const CeedScalar **q_weight) {
2318 *q_weight = basis->q_weight_1d;
2319 return CEED_ERROR_SUCCESS;
2320 }
2321
2322 /**
2323 @brief Get interpolation matrix of a `CeedBasis`
2324
2325 @param[in] basis `CeedBasis`
2326 @param[out] interp Variable to store interpolation matrix
2327
2328 @return An error code: 0 - success, otherwise - failure
2329
2330 @ref Advanced
2331 **/
CeedBasisGetInterp(CeedBasis basis,const CeedScalar ** interp)2332 int CeedBasisGetInterp(CeedBasis basis, const CeedScalar **interp) {
2333 if (!basis->interp && basis->is_tensor_basis) {
2334 // Allocate
2335 CeedCall(CeedMalloc(basis->Q * basis->P, &basis->interp));
2336
2337 // Initialize
2338 for (CeedInt i = 0; i < basis->Q * basis->P; i++) basis->interp[i] = 1.0;
2339
2340 // Calculate
2341 for (CeedInt d = 0; d < basis->dim; d++) {
2342 for (CeedInt qpt = 0; qpt < basis->Q; qpt++) {
2343 for (CeedInt node = 0; node < basis->P; node++) {
2344 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d;
2345 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d;
2346
2347 basis->interp[qpt * (basis->P) + node] *= basis->interp_1d[q * basis->P_1d + p];
2348 }
2349 }
2350 }
2351 }
2352 *interp = basis->interp;
2353 return CEED_ERROR_SUCCESS;
2354 }
2355
2356 /**
2357 @brief Get 1D interpolation matrix of a tensor product `CeedBasis`
2358
2359 @param[in] basis `CeedBasis`
2360 @param[out] interp_1d Variable to store interpolation matrix
2361
2362 @return An error code: 0 - success, otherwise - failure
2363
2364 @ref Backend
2365 **/
CeedBasisGetInterp1D(CeedBasis basis,const CeedScalar ** interp_1d)2366 int CeedBasisGetInterp1D(CeedBasis basis, const CeedScalar **interp_1d) {
2367 bool is_tensor_basis;
2368
2369 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis));
2370 CeedCheck(is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "CeedBasis is not a tensor product CeedBasis");
2371 *interp_1d = basis->interp_1d;
2372 return CEED_ERROR_SUCCESS;
2373 }
2374
2375 /**
2376 @brief Get gradient matrix of a `CeedBasis`
2377
2378 @param[in] basis `CeedBasis`
2379 @param[out] grad Variable to store gradient matrix
2380
2381 @return An error code: 0 - success, otherwise - failure
2382
2383 @ref Advanced
2384 **/
CeedBasisGetGrad(CeedBasis basis,const CeedScalar ** grad)2385 int CeedBasisGetGrad(CeedBasis basis, const CeedScalar **grad) {
2386 if (!basis->grad && basis->is_tensor_basis) {
2387 // Allocate
2388 CeedCall(CeedMalloc(basis->dim * basis->Q * basis->P, &basis->grad));
2389
2390 // Initialize
2391 for (CeedInt i = 0; i < basis->dim * basis->Q * basis->P; i++) basis->grad[i] = 1.0;
2392
2393 // Calculate
2394 for (CeedInt d = 0; d < basis->dim; d++) {
2395 for (CeedInt i = 0; i < basis->dim; i++) {
2396 for (CeedInt qpt = 0; qpt < basis->Q; qpt++) {
2397 for (CeedInt node = 0; node < basis->P; node++) {
2398 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d;
2399 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d;
2400
2401 if (i == d) basis->grad[(i * basis->Q + qpt) * (basis->P) + node] *= basis->grad_1d[q * basis->P_1d + p];
2402 else basis->grad[(i * basis->Q + qpt) * (basis->P) + node] *= basis->interp_1d[q * basis->P_1d + p];
2403 }
2404 }
2405 }
2406 }
2407 }
2408 *grad = basis->grad;
2409 return CEED_ERROR_SUCCESS;
2410 }
2411
2412 /**
2413 @brief Get 1D gradient matrix of a tensor product `CeedBasis`
2414
2415 @param[in] basis `CeedBasis`
2416 @param[out] grad_1d Variable to store gradient matrix
2417
2418 @return An error code: 0 - success, otherwise - failure
2419
2420 @ref Advanced
2421 **/
CeedBasisGetGrad1D(CeedBasis basis,const CeedScalar ** grad_1d)2422 int CeedBasisGetGrad1D(CeedBasis basis, const CeedScalar **grad_1d) {
2423 bool is_tensor_basis;
2424
2425 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis));
2426 CeedCheck(is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "CeedBasis is not a tensor product CeedBasis");
2427 *grad_1d = basis->grad_1d;
2428 return CEED_ERROR_SUCCESS;
2429 }
2430
2431 /**
2432 @brief Get divergence matrix of a `CeedBasis`
2433
2434 @param[in] basis `CeedBasis`
2435 @param[out] div Variable to store divergence matrix
2436
2437 @return An error code: 0 - success, otherwise - failure
2438
2439 @ref Advanced
2440 **/
CeedBasisGetDiv(CeedBasis basis,const CeedScalar ** div)2441 int CeedBasisGetDiv(CeedBasis basis, const CeedScalar **div) {
2442 *div = basis->div;
2443 return CEED_ERROR_SUCCESS;
2444 }
2445
2446 /**
2447 @brief Get curl matrix of a `CeedBasis`
2448
2449 @param[in] basis `CeedBasis`
2450 @param[out] curl Variable to store curl matrix
2451
2452 @return An error code: 0 - success, otherwise - failure
2453
2454 @ref Advanced
2455 **/
CeedBasisGetCurl(CeedBasis basis,const CeedScalar ** curl)2456 int CeedBasisGetCurl(CeedBasis basis, const CeedScalar **curl) {
2457 *curl = basis->curl;
2458 return CEED_ERROR_SUCCESS;
2459 }
2460
2461 /**
2462 @brief Destroy a @ref CeedBasis
2463
2464 @param[in,out] basis `CeedBasis` to destroy
2465
2466 @return An error code: 0 - success, otherwise - failure
2467
2468 @ref User
2469 **/
CeedBasisDestroy(CeedBasis * basis)2470 int CeedBasisDestroy(CeedBasis *basis) {
2471 if (!*basis || *basis == CEED_BASIS_NONE || CeedObjectDereference((CeedObject)*basis) > 0) {
2472 *basis = NULL;
2473 return CEED_ERROR_SUCCESS;
2474 }
2475 if ((*basis)->Destroy) CeedCall((*basis)->Destroy(*basis));
2476 CeedCall(CeedTensorContractDestroy(&(*basis)->contract));
2477 CeedCall(CeedFree(&(*basis)->q_ref_1d));
2478 CeedCall(CeedFree(&(*basis)->q_weight_1d));
2479 CeedCall(CeedFree(&(*basis)->interp));
2480 CeedCall(CeedFree(&(*basis)->interp_1d));
2481 CeedCall(CeedFree(&(*basis)->grad));
2482 CeedCall(CeedFree(&(*basis)->grad_1d));
2483 CeedCall(CeedFree(&(*basis)->div));
2484 CeedCall(CeedFree(&(*basis)->curl));
2485 CeedCall(CeedVectorDestroy(&(*basis)->vec_chebyshev));
2486 CeedCall(CeedBasisDestroy(&(*basis)->basis_chebyshev));
2487 CeedCall(CeedObjectDestroy_Private(&(*basis)->obj));
2488 CeedCall(CeedFree(basis));
2489 return CEED_ERROR_SUCCESS;
2490 }
2491
2492 /**
2493 @brief Construct a Gauss-Legendre quadrature
2494
2495 @param[in] Q Number of quadrature points (integrates polynomials of degree `2*Q-1` exactly)
2496 @param[out] q_ref_1d Array of length `Q` to hold the abscissa on `[-1, 1]`
2497 @param[out] q_weight_1d Array of length `Q` to hold the weights
2498
2499 @return An error code: 0 - success, otherwise - failure
2500
2501 @ref Utility
2502 **/
CeedGaussQuadrature(CeedInt Q,CeedScalar * q_ref_1d,CeedScalar * q_weight_1d)2503 int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d, CeedScalar *q_weight_1d) {
2504 CeedScalar P0, P1, P2, dP2, xi, wi, PI = 4.0 * atan(1.0);
2505
2506 // Build q_ref_1d, q_weight_1d
2507 for (CeedInt i = 0; i <= Q / 2; i++) {
2508 // Guess
2509 xi = cos(PI * (CeedScalar)(2 * i + 1) / ((CeedScalar)(2 * Q)));
2510 // Pn(xi)
2511 P0 = 1.0;
2512 P1 = xi;
2513 P2 = 0.0;
2514 for (CeedInt j = 2; j <= Q; j++) {
2515 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j));
2516 P0 = P1;
2517 P1 = P2;
2518 }
2519 // First Newton Step
2520 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0);
2521 xi = xi - P2 / dP2;
2522 // Newton to convergence
2523 for (CeedInt k = 0; k < 100 && fabs(P2) > 10 * CEED_EPSILON; k++) {
2524 P0 = 1.0;
2525 P1 = xi;
2526 for (CeedInt j = 2; j <= Q; j++) {
2527 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j));
2528 P0 = P1;
2529 P1 = P2;
2530 }
2531 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0);
2532 xi = xi - P2 / dP2;
2533 }
2534 // Save xi, wi
2535 wi = 2.0 / ((1.0 - xi * xi) * dP2 * dP2);
2536 q_weight_1d[i] = wi;
2537 q_weight_1d[Q - 1 - i] = wi;
2538 q_ref_1d[i] = -xi;
2539 q_ref_1d[Q - 1 - i] = xi;
2540 }
2541 return CEED_ERROR_SUCCESS;
2542 }
2543
2544 /**
2545 @brief Construct a Gauss-Legendre-Lobatto quadrature
2546
2547 @param[in] Q Number of quadrature points (integrates polynomials of degree `2*Q-3` exactly)
2548 @param[out] q_ref_1d Array of length `Q` to hold the abscissa on `[-1, 1]`
2549 @param[out] q_weight_1d Array of length `Q` to hold the weights
2550
2551 @return An error code: 0 - success, otherwise - failure
2552
2553 @ref Utility
2554 **/
CeedLobattoQuadrature(CeedInt Q,CeedScalar * q_ref_1d,CeedScalar * q_weight_1d)2555 int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d, CeedScalar *q_weight_1d) {
2556 CeedScalar P0, P1, P2, dP2, d2P2, xi, wi, PI = 4.0 * atan(1.0);
2557
2558 // Build q_ref_1d, q_weight_1d
2559 // Set endpoints
2560 CeedCheck(Q > 1, NULL, CEED_ERROR_DIMENSION, "Cannot create Lobatto quadrature with Q=%" CeedInt_FMT " < 2 points", Q);
2561 wi = 2.0 / ((CeedScalar)(Q * (Q - 1)));
2562 if (q_weight_1d) {
2563 q_weight_1d[0] = wi;
2564 q_weight_1d[Q - 1] = wi;
2565 }
2566 q_ref_1d[0] = -1.0;
2567 q_ref_1d[Q - 1] = 1.0;
2568 // Interior
2569 for (CeedInt i = 1; i <= (Q - 1) / 2; i++) {
2570 // Guess
2571 xi = cos(PI * (CeedScalar)(i) / (CeedScalar)(Q - 1));
2572 // Pn(xi)
2573 P0 = 1.0;
2574 P1 = xi;
2575 P2 = 0.0;
2576 for (CeedInt j = 2; j < Q; j++) {
2577 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j));
2578 P0 = P1;
2579 P1 = P2;
2580 }
2581 // First Newton step
2582 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0);
2583 d2P2 = (2 * xi * dP2 - (CeedScalar)(Q * (Q - 1)) * P2) / (1.0 - xi * xi);
2584 xi = xi - dP2 / d2P2;
2585 // Newton to convergence
2586 for (CeedInt k = 0; k < 100 && fabs(dP2) > 10 * CEED_EPSILON; k++) {
2587 P0 = 1.0;
2588 P1 = xi;
2589 for (CeedInt j = 2; j < Q; j++) {
2590 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j));
2591 P0 = P1;
2592 P1 = P2;
2593 }
2594 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0);
2595 d2P2 = (2 * xi * dP2 - (CeedScalar)(Q * (Q - 1)) * P2) / (1.0 - xi * xi);
2596 xi = xi - dP2 / d2P2;
2597 }
2598 // Save xi, wi
2599 wi = 2.0 / (((CeedScalar)(Q * (Q - 1))) * P2 * P2);
2600 if (q_weight_1d) {
2601 q_weight_1d[i] = wi;
2602 q_weight_1d[Q - 1 - i] = wi;
2603 }
2604 q_ref_1d[i] = -xi;
2605 q_ref_1d[Q - 1 - i] = xi;
2606 }
2607 return CEED_ERROR_SUCCESS;
2608 }
2609
2610 /// @}
2611