ceed-basis.c (38d0029a08e84b16bbdde77c27446234f5f09cf2) ceed-basis.c (7a982d89c751e293e39d23a7c44a161cef1fcd06)
1// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3// reserved. See files LICENSE and NOTICE for details.
4//
5// This file is part of CEED, a collection of benchmarks, miniapps, software
6// libraries and APIs for efficient high-order finite element and spectral
7// element discretizations for exascale applications. For more information and
8// source code availability see http://github.com/ceed.

--- 7 unchanged lines hidden (view full) ---

16
17#include <ceed-impl.h>
18#include <ceed-backend.h>
19#include <math.h>
20#include <stdio.h>
21#include <stdlib.h>
22#include <string.h>
23
1// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3// reserved. See files LICENSE and NOTICE for details.
4//
5// This file is part of CEED, a collection of benchmarks, miniapps, software
6// libraries and APIs for efficient high-order finite element and spectral
7// element discretizations for exascale applications. For more information and
8// source code availability see http://github.com/ceed.

--- 7 unchanged lines hidden (view full) ---

16
17#include <ceed-impl.h>
18#include <ceed-backend.h>
19#include <math.h>
20#include <stdio.h>
21#include <stdlib.h>
22#include <string.h>
23
24/// @file
25/// Implementation of CeedBasis interfaces
26
24/// @cond DOXYGEN_SKIP
25static struct CeedBasis_private ceed_basis_collocated;
26/// @endcond
27
27/// @cond DOXYGEN_SKIP
28static struct CeedBasis_private ceed_basis_collocated;
29/// @endcond
30
28/// @file
29/// Implementation of public CeedBasis interfaces
30///
31/// @addtogroup CeedBasis
31/// @addtogroup CeedBasisUser
32/// @{
33
32/// @{
33
34/// Indicate that the quadrature points are collocated with the nodes
35const CeedBasis CEED_BASIS_COLLOCATED = &ceed_basis_collocated;
36
37/// @}
38
39/// ----------------------------------------------------------------------------
40/// CeedBasis Library Internal Functions
41/// ----------------------------------------------------------------------------
42/// @addtogroup CeedBasisDeveloper
43/// @{
44
34/**
45/**
46 @brief Compute Householder reflection
47
48 Computes A = (I - b v v^T) A
49 where A is an mxn matrix indexed as A[i*row + j*col]
50
51 @param[in,out] A Matrix to apply Householder reflection to, in place
52 @param v Householder vector
53 @param b Scaling factor
54 @param m Number of rows in A
55 @param n Number of columns in A
56 @param row Row stride
57 @param col Col stride
58
59 @return An error code: 0 - success, otherwise - failure
60
61 @ref Developer
62**/
63static int CeedHouseholderReflect(CeedScalar *A, const CeedScalar *v,
64 CeedScalar b, CeedInt m, CeedInt n,
65 CeedInt row, CeedInt col) {
66 for (CeedInt j=0; j<n; j++) {
67 CeedScalar w = A[0*row + j*col];
68 for (CeedInt i=1; i<m; i++)
69 w += v[i] * A[i*row + j*col];
70 A[0*row + j*col] -= b * w;
71 for (CeedInt i=1; i<m; i++)
72 A[i*row + j*col] -= b * w * v[i];
73 }
74 return 0;
75}
76
77/**
78 @brief Apply Householder Q matrix
79
80 Compute A = Q A where Q is mxm and A is mxn.
81
82 @param[in,out] A Matrix to apply Householder Q to, in place
83 @param Q Householder Q matrix
84 @param tau Householder scaling factors
85 @param tmode Transpose mode for application
86 @param m Number of rows in A
87 @param n Number of columns in A
88 @param k Number of elementary reflectors in Q, k<m
89 @param row Row stride in A
90 @param col Col stride in A
91
92 @return An error code: 0 - success, otherwise - failure
93
94 @ref Developer
95**/
96static int CeedHouseholderApplyQ(CeedScalar *A, const CeedScalar *Q,
97 const CeedScalar *tau, CeedTransposeMode tmode,
98 CeedInt m, CeedInt n, CeedInt k,
99 CeedInt row, CeedInt col) {
100 CeedScalar v[m];
101 for (CeedInt ii=0; ii<k; ii++) {
102 CeedInt i = tmode == CEED_TRANSPOSE ? ii : k-1-ii;
103 for (CeedInt j=i+1; j<m; j++)
104 v[j] = Q[j*k+i];
105 // Apply Householder reflector (I - tau v v^T) collograd1d^T
106 CeedHouseholderReflect(&A[i*row], &v[i], tau[i], m-i, n, row, col);
107 }
108 return 0;
109}
110
111/**
112 @brief Compute Givens rotation
113
114 Computes A = G A (or G^T A in transpose mode)
115 where A is an mxn matrix indexed as A[i*n + j*m]
116
117 @param[in,out] A Row major matrix to apply Givens rotation to, in place
118 @param c Cosine factor
119 @param s Sine factor
120 @param i First row/column to apply rotation
121 @param k Second row/column to apply rotation
122 @param m Number of rows in A
123 @param n Number of columns in A
124
125 @return An error code: 0 - success, otherwise - failure
126
127 @ref Developer
128**/
129static int CeedGivensRotation(CeedScalar *A, CeedScalar c, CeedScalar s,
130 CeedTransposeMode tmode, CeedInt i, CeedInt k,
131 CeedInt m, CeedInt n) {
132 CeedInt stridej = 1, strideik = m, numits = n;
133 if (tmode == CEED_NOTRANSPOSE) {
134 stridej = n; strideik = 1; numits = m;
135 }
136
137 // Apply rotation
138 for (CeedInt j=0; j<numits; j++) {
139 CeedScalar tau1 = A[i*strideik+j*stridej], tau2 = A[k*strideik+j*stridej];
140 A[i*strideik+j*stridej] = c*tau1 - s*tau2;
141 A[k*strideik+j*stridej] = s*tau1 + c*tau2;
142 }
143
144 return 0;
145}
146
147/**
148 @brief View an array stored in a CeedBasis
149
150 @param name Name of array
151 @param fpformat Printing format
152 @param m Number of rows in array
153 @param n Number of columns in array
154 @param a Array to be viewed
155 @param stream Stream to view to, e.g., stdout
156
157 @return An error code: 0 - success, otherwise - failure
158
159 @ref Developer
160**/
161static int CeedScalarView(const char *name, const char *fpformat, CeedInt m,
162 CeedInt n, const CeedScalar *a, FILE *stream) {
163 for (int i=0; i<m; i++) {
164 if (m > 1)
165 fprintf(stream, "%12s[%d]:", name, i);
166 else
167 fprintf(stream, "%12s:", name);
168 for (int j=0; j<n; j++)
169 fprintf(stream, fpformat, fabs(a[i*n+j]) > 1E-14 ? a[i*n+j] : 0);
170 fputs("\n", stream);
171 }
172 return 0;
173}
174
175/// @}
176
177/// ----------------------------------------------------------------------------
178/// Ceed Backend API
179/// ----------------------------------------------------------------------------
180/// @addtogroup CeedBasisBackend
181/// @{
182
183/**
184 @brief Return collocated grad matrix
185
186 @param basis CeedBasis
187 @param[out] collograd1d Row-major (Q1d * Q1d) matrix expressing derivatives of
188 basis functions at quadrature points
189
190 @return An error code: 0 - success, otherwise - failure
191
192 @ref Backend
193**/
194int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collograd1d) {
195 int i, j, k;
196 Ceed ceed;
197 CeedInt ierr, P1d=(basis)->P1d, Q1d=(basis)->Q1d;
198 CeedScalar *interp1d, *grad1d, tau[Q1d];
199
200 ierr = CeedMalloc(Q1d*P1d, &interp1d); CeedChk(ierr);
201 ierr = CeedMalloc(Q1d*P1d, &grad1d); CeedChk(ierr);
202 memcpy(interp1d, (basis)->interp1d, Q1d*P1d*sizeof(basis)->interp1d[0]);
203 memcpy(grad1d, (basis)->grad1d, Q1d*P1d*sizeof(basis)->interp1d[0]);
204
205 // QR Factorization, interp1d = Q R
206 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr);
207 ierr = CeedQRFactorization(ceed, interp1d, tau, Q1d, P1d); CeedChk(ierr);
208
209 // Apply Rinv, collograd1d = grad1d Rinv
210 for (i=0; i<Q1d; i++) { // Row i
211 collograd1d[Q1d*i] = grad1d[P1d*i]/interp1d[0];
212 for (j=1; j<P1d; j++) { // Column j
213 collograd1d[j+Q1d*i] = grad1d[j+P1d*i];
214 for (k=0; k<j; k++)
215 collograd1d[j+Q1d*i] -= interp1d[j+P1d*k]*collograd1d[k+Q1d*i];
216 collograd1d[j+Q1d*i] /= interp1d[j+P1d*j];
217 }
218 for (j=P1d; j<Q1d; j++)
219 collograd1d[j+Q1d*i] = 0;
220 }
221
222 // Apply Qtranspose, collograd = collograd Qtranspose
223 CeedHouseholderApplyQ(collograd1d, interp1d, tau, CEED_NOTRANSPOSE,
224 Q1d, Q1d, P1d, 1, Q1d);
225
226 ierr = CeedFree(&interp1d); CeedChk(ierr);
227 ierr = CeedFree(&grad1d); CeedChk(ierr);
228
229 return 0;
230}
231
232/**
233 @brief Get Ceed associated with a CeedBasis
234
235 @param basis CeedBasis
236 @param[out] ceed Variable to store Ceed
237
238 @return An error code: 0 - success, otherwise - failure
239
240 @ref Backend
241**/
242int CeedBasisGetCeed(CeedBasis basis, Ceed *ceed) {
243 *ceed = basis->ceed;
244 return 0;
245}
246
247/**
248 @brief Get tensor status for given CeedBasis
249
250 @param basis CeedBasis
251 @param[out] tensor Variable to store tensor status
252
253 @return An error code: 0 - success, otherwise - failure
254
255 @ref Backend
256**/
257int CeedBasisGetTensorStatus(CeedBasis basis, bool *tensor) {
258 *tensor = basis->tensorbasis;
259 return 0;
260}
261
262/**
263 @brief Get dimension for given CeedBasis
264
265 @param basis CeedBasis
266 @param[out] dim Variable to store dimension of basis
267
268 @return An error code: 0 - success, otherwise - failure
269
270 @ref Backend
271**/
272int CeedBasisGetDimension(CeedBasis basis, CeedInt *dim) {
273 *dim = basis->dim;
274 return 0;
275}
276
277/**
278 @brief Get number of components for given CeedBasis
279
280 @param basis CeedBasis
281 @param[out] numcomp Variable to store number of components of basis
282
283 @return An error code: 0 - success, otherwise - failure
284
285 @ref Backend
286**/
287int CeedBasisGetNumComponents(CeedBasis basis, CeedInt *numcomp) {
288 *numcomp = basis->ncomp;
289 return 0;
290}
291
292/**
293 @brief Get total number of nodes (in 1 dimension) of a CeedBasis
294
295 @param basis CeedBasis
296 @param[out] P1d Variable to store number of nodes
297
298 @return An error code: 0 - success, otherwise - failure
299
300 @ref Backend
301**/
302int CeedBasisGetNumNodes1D(CeedBasis basis, CeedInt *P1d) {
303 if (!basis->tensorbasis)
304 // LCOV_EXCL_START
305 return CeedError(basis->ceed, 1, "Cannot supply P1d for non-tensor basis");
306 // LCOV_EXCL_STOP
307
308 *P1d = basis->P1d;
309 return 0;
310}
311
312/**
313 @brief Get total number of quadrature points (in 1 dimension) of a CeedBasis
314
315 @param basis CeedBasis
316 @param[out] Q1d Variable to store number of quadrature points
317
318 @return An error code: 0 - success, otherwise - failure
319
320 @ref Backend
321**/
322int CeedBasisGetNumQuadraturePoints1D(CeedBasis basis, CeedInt *Q1d) {
323 if (!basis->tensorbasis)
324 // LCOV_EXCL_START
325 return CeedError(basis->ceed, 1, "Cannot supply Q1d for non-tensor basis");
326 // LCOV_EXCL_STOP
327
328 *Q1d = basis->Q1d;
329 return 0;
330}
331
332/**
333 @brief Get reference coordinates of quadrature points (in dim dimensions)
334 of a CeedBasis
335
336 @param basis CeedBasis
337 @param[out] qref Variable to store reference coordinates of quadrature points
338
339 @return An error code: 0 - success, otherwise - failure
340
341 @ref Backend
342**/
343int CeedBasisGetQRef(CeedBasis basis, CeedScalar **qref) {
344 *qref = basis->qref1d;
345 return 0;
346}
347
348/**
349 @brief Get quadrature weights of quadrature points (in dim dimensions)
350 of a CeedBasis
351
352 @param basis CeedBasis
353 @param[out] qweight Variable to store quadrature weights
354
355 @return An error code: 0 - success, otherwise - failure
356
357 @ref Backend
358**/
359int CeedBasisGetQWeights(CeedBasis basis, CeedScalar **qweight) {
360 *qweight = basis->qweight1d;
361 return 0;
362}
363
364/**
365 @brief Get interpolation matrix of a CeedBasis
366
367 @param basis CeedBasis
368 @param[out] interp Variable to store interpolation matrix
369
370 @return An error code: 0 - success, otherwise - failure
371
372 @ref Backend
373**/
374int CeedBasisGetInterp(CeedBasis basis, CeedScalar **interp) {
375 if (!basis->interp && basis->tensorbasis) {
376 // Allocate
377 int ierr;
378 ierr = CeedMalloc(basis->Q*basis->P, &basis->interp); CeedChk(ierr);
379
380 // Initialize
381 for (CeedInt i=0; i<basis->Q*basis->P; i++)
382 basis->interp[i] = 1.0;
383
384 // Calculate
385 for (CeedInt d=0; d<basis->dim; d++)
386 for (CeedInt qpt=0; qpt<basis->Q; qpt++)
387 for (CeedInt node=0; node<basis->P; node++) {
388 CeedInt p = (node / CeedIntPow(basis->P1d, d)) % basis->P1d;
389 CeedInt q = (qpt / CeedIntPow(basis->Q1d, d)) % basis->Q1d;
390 basis->interp[qpt*(basis->P)+node] *= basis->interp1d[q*basis->P1d+p];
391 }
392 }
393
394 *interp = basis->interp;
395
396 return 0;
397}
398
399/**
400 @brief Get 1D interpolation matrix of a tensor product CeedBasis
401
402 @param basis CeedBasis
403 @param[out] interp1d Variable to store interpolation matrix
404
405 @return An error code: 0 - success, otherwise - failure
406
407 @ref Backend
408**/
409int CeedBasisGetInterp1D(CeedBasis basis, CeedScalar **interp1d) {
410 if (!basis->tensorbasis)
411 // LCOV_EXCL_START
412 return CeedError(basis->ceed, 1, "CeedBasis is not a tensor product basis.");
413 // LCOV_EXCL_STOP
414
415 *interp1d = basis->interp1d;
416
417 return 0;
418}
419
420/**
421 @brief Get gradient matrix of a CeedBasis
422
423 @param basis CeedBasis
424 @param[out] grad Variable to store gradient matrix
425
426 @return An error code: 0 - success, otherwise - failure
427
428 @ref Backend
429**/
430int CeedBasisGetGrad(CeedBasis basis, CeedScalar **grad) {
431 if (!basis->grad && basis->tensorbasis) {
432 // Allocate
433 int ierr;
434 ierr = CeedMalloc(basis->dim*basis->Q*basis->P, &basis->grad);
435 CeedChk(ierr);
436
437 // Initialize
438 for (CeedInt i=0; i<basis->dim*basis->Q*basis->P; i++)
439 basis->grad[i] = 1.0;
440
441 // Calculate
442 for (CeedInt d=0; d<basis->dim; d++)
443 for (CeedInt i=0; i<basis->dim; i++)
444 for (CeedInt qpt=0; qpt<basis->Q; qpt++)
445 for (CeedInt node=0; node<basis->P; node++) {
446 CeedInt p = (node / CeedIntPow(basis->P1d, d)) % basis->P1d;
447 CeedInt q = (qpt / CeedIntPow(basis->Q1d, d)) % basis->Q1d;
448 if (i == d)
449 basis->grad[(i*basis->Q+qpt)*(basis->P)+node] *=
450 basis->grad1d[q*basis->P1d+p];
451 else
452 basis->grad[(i*basis->Q+qpt)*(basis->P)+node] *=
453 basis->interp1d[q*basis->P1d+p];
454 }
455 }
456
457 *grad = basis->grad;
458
459 return 0;
460}
461
462/**
463 @brief Get 1D gradient matrix of a tensor product CeedBasis
464
465 @param basis CeedBasis
466 @param[out] grad1d Variable to store gradient matrix
467
468 @return An error code: 0 - success, otherwise - failure
469
470 @ref Backend
471**/
472int CeedBasisGetGrad1D(CeedBasis basis, CeedScalar **grad1d) {
473 if (!basis->tensorbasis)
474 // LCOV_EXCL_START
475 return CeedError(basis->ceed, 1, "CeedBasis is not a tensor product basis.");
476 // LCOV_EXCL_STOP
477
478 *grad1d = basis->grad1d;
479
480 return 0;
481}
482
483/**
484 @brief Get backend data of a CeedBasis
485
486 @param basis CeedBasis
487 @param[out] data Variable to store data
488
489 @return An error code: 0 - success, otherwise - failure
490
491 @ref Backend
492**/
493int CeedBasisGetData(CeedBasis basis, void **data) {
494 *data = basis->data;
495 return 0;
496}
497
498/**
499 @brief Set backend data of a CeedBasis
500
501 @param[out] basis CeedBasis
502 @param data Data to set
503
504 @return An error code: 0 - success, otherwise - failure
505
506 @ref Backend
507**/
508int CeedBasisSetData(CeedBasis basis, void **data) {
509 basis->data = *data;
510 return 0;
511}
512
513/**
514 @brief Get dimension for given CeedElemTopology
515
516 @param topo CeedElemTopology
517 @param[out] dim Variable to store dimension of topology
518
519 @return An error code: 0 - success, otherwise - failure
520
521 @ref Backend
522**/
523int CeedBasisGetTopologyDimension(CeedElemTopology topo, CeedInt *dim) {
524 *dim = (CeedInt) topo >> 16;
525 return 0;
526}
527
528/**
529 @brief Get CeedTensorContract of a CeedBasis
530
531 @param basis CeedBasis
532 @param[out] contract Variable to store CeedTensorContract
533
534 @return An error code: 0 - success, otherwise - failure
535
536 @ref Backend
537**/
538int CeedBasisGetTensorContract(CeedBasis basis, CeedTensorContract *contract) {
539 *contract = basis->contract;
540 return 0;
541}
542
543/**
544 @brief Set CeedTensorContract of a CeedBasis
545
546 @param[out] basis CeedBasis
547 @param contract CeedTensorContract to set
548
549 @return An error code: 0 - success, otherwise - failure
550
551 @ref Backend
552**/
553int CeedBasisSetTensorContract(CeedBasis basis, CeedTensorContract *contract) {
554 basis->contract = *contract;
555 return 0;
556}
557
558/**
559 @brief Return a reference implementation of matrix multiplication C = A B.
560 Note, this is a reference implementation for CPU CeedScalar pointers
561 that is not intended for high performance.
562
563 @param ceed A Ceed context for error handling
564 @param[in] matA Row-major matrix A
565 @param[in] matB Row-major matrix B
566 @param[out] matC Row-major output matrix C
567 @param m Number of rows of C
568 @param n Number of columns of C
569 @param kk Number of columns of A/rows of B
570
571 @return An error code: 0 - success, otherwise - failure
572
573 @ref Utility
574**/
575int CeedMatrixMultiply(Ceed ceed, const CeedScalar *matA,
576 const CeedScalar *matB, CeedScalar *matC, CeedInt m,
577 CeedInt n, CeedInt kk) {
578 for (CeedInt i=0; i<m; i++)
579 for (CeedInt j=0; j<n; j++) {
580 CeedScalar sum = 0;
581 for (CeedInt k=0; k<kk; k++)
582 sum += matA[k+i*kk]*matB[j+k*n];
583 matC[j+i*n] = sum;
584 }
585 return 0;
586}
587
588/// @}
589
590/// ----------------------------------------------------------------------------
591/// CeedBasis Public API
592/// ----------------------------------------------------------------------------
593/// @addtogroup CeedBasisUser
594/// @{
595
596/**
35 @brief Create a tensor-product basis for H^1 discretizations
36
37 @param ceed A Ceed object where the CeedBasis will be created
38 @param dim Topological dimension
39 @param ncomp Number of field components (1 for scalar fields)
40 @param P1d Number of nodes in one dimension
41 @param Q1d Number of quadrature points in one dimension
42 @param interp1d Row-major (Q1d * P1d) matrix expressing the values of nodal

--- 4 unchanged lines hidden (view full) ---

47 on the 1D reference element [-1, 1]
48 @param qweight1d Array of length Q1d holding the quadrature weights on the
49 reference element
50 @param[out] basis Address of the variable where the newly created
51 CeedBasis will be stored.
52
53 @return An error code: 0 - success, otherwise - failure
54
597 @brief Create a tensor-product basis for H^1 discretizations
598
599 @param ceed A Ceed object where the CeedBasis will be created
600 @param dim Topological dimension
601 @param ncomp Number of field components (1 for scalar fields)
602 @param P1d Number of nodes in one dimension
603 @param Q1d Number of quadrature points in one dimension
604 @param interp1d Row-major (Q1d * P1d) matrix expressing the values of nodal

--- 4 unchanged lines hidden (view full) ---

609 on the 1D reference element [-1, 1]
610 @param qweight1d Array of length Q1d holding the quadrature weights on the
611 reference element
612 @param[out] basis Address of the variable where the newly created
613 CeedBasis will be stored.
614
615 @return An error code: 0 - success, otherwise - failure
616
55 @ref Basic
617 @ref User
56**/
57int CeedBasisCreateTensorH1(Ceed ceed, CeedInt dim, CeedInt ncomp, CeedInt P1d,
58 CeedInt Q1d, const CeedScalar *interp1d,
59 const CeedScalar *grad1d, const CeedScalar *qref1d,
60 const CeedScalar *qweight1d, CeedBasis *basis) {
61 int ierr;
62
63 if (dim<1)

--- 50 unchanged lines hidden (view full) ---

114 @param Q Number of quadrature points in one dimension.
115 @param qmode Distribution of the Q quadrature points (affects order of
116 accuracy for the quadrature)
117 @param[out] basis Address of the variable where the newly created
118 CeedBasis will be stored.
119
120 @return An error code: 0 - success, otherwise - failure
121
618**/
619int CeedBasisCreateTensorH1(Ceed ceed, CeedInt dim, CeedInt ncomp, CeedInt P1d,
620 CeedInt Q1d, const CeedScalar *interp1d,
621 const CeedScalar *grad1d, const CeedScalar *qref1d,
622 const CeedScalar *qweight1d, CeedBasis *basis) {
623 int ierr;
624
625 if (dim<1)

--- 50 unchanged lines hidden (view full) ---

676 @param Q Number of quadrature points in one dimension.
677 @param qmode Distribution of the Q quadrature points (affects order of
678 accuracy for the quadrature)
679 @param[out] basis Address of the variable where the newly created
680 CeedBasis will be stored.
681
682 @return An error code: 0 - success, otherwise - failure
683
122 @ref Basic
684 @ref User
123**/
124int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt ncomp,
125 CeedInt P, CeedInt Q, CeedQuadMode qmode,
126 CeedBasis *basis) {
127 // Allocate
128 int ierr, i, j, k;
129 CeedScalar c1, c2, c3, c4, dx, *nodes, *interp1d, *grad1d, *qref1d, *qweight1d;
130

--- 67 unchanged lines hidden (view full) ---

198 points on the reference element [-1, 1]
199 @param qweight Array of length nqpts holding the quadrature weights on the
200 reference element
201 @param[out] basis Address of the variable where the newly created
202 CeedBasis will be stored.
203
204 @return An error code: 0 - success, otherwise - failure
205
685**/
686int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt ncomp,
687 CeedInt P, CeedInt Q, CeedQuadMode qmode,
688 CeedBasis *basis) {
689 // Allocate
690 int ierr, i, j, k;
691 CeedScalar c1, c2, c3, c4, dx, *nodes, *interp1d, *grad1d, *qref1d, *qweight1d;
692

--- 67 unchanged lines hidden (view full) ---

760 points on the reference element [-1, 1]
761 @param qweight Array of length nqpts holding the quadrature weights on the
762 reference element
763 @param[out] basis Address of the variable where the newly created
764 CeedBasis will be stored.
765
766 @return An error code: 0 - success, otherwise - failure
767
206 @ref Basic
768 @ref User
207**/
208int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt ncomp,
209 CeedInt nnodes, CeedInt nqpts, const CeedScalar *interp,
210 const CeedScalar *grad, const CeedScalar *qref,
211 const CeedScalar *qweight, CeedBasis *basis) {
212 int ierr;
213 CeedInt P = nnodes, Q = nqpts, dim = 0;
214

--- 33 unchanged lines hidden (view full) ---

248 memcpy((*basis)->interp, interp, Q*P*sizeof(interp[0]));
249 memcpy((*basis)->grad, grad, dim*Q*P*sizeof(grad[0]));
250 ierr = ceed->BasisCreateH1(topo, dim, P, Q, interp, grad, qref,
251 qweight, *basis); CeedChk(ierr);
252 return 0;
253}
254
255/**
769**/
770int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt ncomp,
771 CeedInt nnodes, CeedInt nqpts, const CeedScalar *interp,
772 const CeedScalar *grad, const CeedScalar *qref,
773 const CeedScalar *qweight, CeedBasis *basis) {
774 int ierr;
775 CeedInt P = nnodes, Q = nqpts, dim = 0;
776

--- 33 unchanged lines hidden (view full) ---

810 memcpy((*basis)->interp, interp, Q*P*sizeof(interp[0]));
811 memcpy((*basis)->grad, grad, dim*Q*P*sizeof(grad[0]));
812 ierr = ceed->BasisCreateH1(topo, dim, P, Q, interp, grad, qref,
813 qweight, *basis); CeedChk(ierr);
814 return 0;
815}
816
817/**
818 @brief View a CeedBasis
819
820 @param basis CeedBasis to view
821 @param stream Stream to view to, e.g., stdout
822
823 @return An error code: 0 - success, otherwise - failure
824
825 @ref User
826**/
827int CeedBasisView(CeedBasis basis, FILE *stream) {
828 int ierr;
829
830 if (basis->tensorbasis) {
831 fprintf(stream, "CeedBasis: dim=%d P=%d Q=%d\n", basis->dim, basis->P1d,
832 basis->Q1d);
833 ierr = CeedScalarView("qref1d", "\t% 12.8f", 1, basis->Q1d, basis->qref1d,
834 stream); CeedChk(ierr);
835 ierr = CeedScalarView("qweight1d", "\t% 12.8f", 1, basis->Q1d,
836 basis->qweight1d, stream); CeedChk(ierr);
837 ierr = CeedScalarView("interp1d", "\t% 12.8f", basis->Q1d, basis->P1d,
838 basis->interp1d, stream); CeedChk(ierr);
839 ierr = CeedScalarView("grad1d", "\t% 12.8f", basis->Q1d, basis->P1d,
840 basis->grad1d, stream); CeedChk(ierr);
841 } else {
842 fprintf(stream, "CeedBasis: dim=%d P=%d Q=%d\n", basis->dim, basis->P,
843 basis->Q);
844 ierr = CeedScalarView("qref", "\t% 12.8f", 1, basis->Q*basis->dim,
845 basis->qref1d,
846 stream); CeedChk(ierr);
847 ierr = CeedScalarView("qweight", "\t% 12.8f", 1, basis->Q, basis->qweight1d,
848 stream); CeedChk(ierr);
849 ierr = CeedScalarView("interp", "\t% 12.8f", basis->Q, basis->P,
850 basis->interp, stream); CeedChk(ierr);
851 ierr = CeedScalarView("grad", "\t% 12.8f", basis->dim*basis->Q, basis->P,
852 basis->grad, stream); CeedChk(ierr);
853 }
854 return 0;
855}
856
857/**
858 @brief Get total number of nodes (in dim dimensions) of a CeedBasis
859
860 @param basis CeedBasis
861 @param[out] P Variable to store number of nodes
862
863 @return An error code: 0 - success, otherwise - failure
864
865 @ref Utility
866**/
867int CeedBasisGetNumNodes(CeedBasis basis, CeedInt *P) {
868 *P = basis->P;
869 return 0;
870}
871
872/**
873 @brief Get total number of quadrature points (in dim dimensions) of a CeedBasis
874
875 @param basis CeedBasis
876 @param[out] Q Variable to store number of quadrature points
877
878 @return An error code: 0 - success, otherwise - failure
879
880 @ref Utility
881**/
882int CeedBasisGetNumQuadraturePoints(CeedBasis basis, CeedInt *Q) {
883 *Q = basis->Q;
884 return 0;
885}
886
887/**
888 @brief Apply basis evaluation from nodes to quadrature points or vice versa
889
890 @param basis CeedBasis to evaluate
891 @param nelem The number of elements to apply the basis evaluation to;
892 the backend will specify the ordering in
893 ElemRestrictionCreateBlocked
894 @param tmode \ref CEED_NOTRANSPOSE to evaluate from nodes to quadrature
895 points, \ref CEED_TRANSPOSE to apply the transpose, mapping
896 from quadrature points to nodes
897 @param emode \ref CEED_EVAL_NONE to use values directly,
898 \ref CEED_EVAL_INTERP to use interpolated values,
899 \ref CEED_EVAL_GRAD to use gradients,
900 \ref CEED_EVAL_WEIGHT to use quadrature weights.
901 @param[in] u Input CeedVector
902 @param[out] v Output CeedVector
903
904 @return An error code: 0 - success, otherwise - failure
905
906 @ref User
907**/
908int CeedBasisApply(CeedBasis basis, CeedInt nelem, CeedTransposeMode tmode,
909 CeedEvalMode emode, CeedVector u, CeedVector v) {
910 int ierr;
911 CeedInt ulength = 0, vlength, nnodes, nqpt;
912 if (!basis->Apply)
913 // LCOV_EXCL_START
914 return CeedError(basis->ceed, 1, "Backend does not support BasisApply");
915 // LCOV_EXCL_STOP
916
917 // Check compatibility of topological and geometrical dimensions
918 ierr = CeedBasisGetNumNodes(basis, &nnodes); CeedChk(ierr);
919 ierr = CeedBasisGetNumQuadraturePoints(basis, &nqpt); CeedChk(ierr);
920 ierr = CeedVectorGetLength(v, &vlength); CeedChk(ierr);
921
922 if (u) {
923 ierr = CeedVectorGetLength(u, &ulength); CeedChk(ierr);
924 }
925
926 if ((tmode == CEED_TRANSPOSE && (vlength%nnodes != 0 || ulength%nqpt != 0)) ||
927 (tmode == CEED_NOTRANSPOSE && (ulength%nnodes != 0 || vlength%nqpt != 0)))
928 return CeedError(basis->ceed, 1, "Length of input/output vectors "
929 "incompatible with basis dimensions");
930
931 ierr = basis->Apply(basis, nelem, tmode, emode, u, v); CeedChk(ierr);
932 return 0;
933}
934
935/**
936 @brief Destroy a CeedBasis
937
938 @param basis CeedBasis to destroy
939
940 @return An error code: 0 - success, otherwise - failure
941
942 @ref User
943**/
944int CeedBasisDestroy(CeedBasis *basis) {
945 int ierr;
946
947 if (!*basis || --(*basis)->refcount > 0)
948 return 0;
949 if ((*basis)->Destroy) {
950 ierr = (*basis)->Destroy(*basis); CeedChk(ierr);
951 }
952 ierr = CeedFree(&(*basis)->interp); CeedChk(ierr);
953 ierr = CeedFree(&(*basis)->interp1d); CeedChk(ierr);
954 ierr = CeedFree(&(*basis)->grad); CeedChk(ierr);
955 ierr = CeedFree(&(*basis)->grad1d); CeedChk(ierr);
956 ierr = CeedFree(&(*basis)->qref1d); CeedChk(ierr);
957 ierr = CeedFree(&(*basis)->qweight1d); CeedChk(ierr);
958 ierr = CeedDestroy(&(*basis)->ceed); CeedChk(ierr);
959 ierr = CeedFree(basis); CeedChk(ierr);
960 return 0;
961}
962
963/**
256 @brief Construct a Gauss-Legendre quadrature
257
258 @param Q Number of quadrature points (integrates polynomials of
259 degree 2*Q-1 exactly)
260 @param[out] qref1d Array of length Q to hold the abscissa on [-1, 1]
261 @param[out] qweight1d Array of length Q to hold the weights
262
263 @return An error code: 0 - success, otherwise - failure

--- 107 unchanged lines hidden (view full) ---

371 }
372 qref1d[i] = -xi;
373 qref1d[Q-1-i]= xi;
374 }
375 return 0;
376}
377
378/**
964 @brief Construct a Gauss-Legendre quadrature
965
966 @param Q Number of quadrature points (integrates polynomials of
967 degree 2*Q-1 exactly)
968 @param[out] qref1d Array of length Q to hold the abscissa on [-1, 1]
969 @param[out] qweight1d Array of length Q to hold the weights
970
971 @return An error code: 0 - success, otherwise - failure

--- 107 unchanged lines hidden (view full) ---

1079 }
1080 qref1d[i] = -xi;
1081 qref1d[Q-1-i]= xi;
1082 }
1083 return 0;
1084}
1085
1086/**
379 @brief View an array stored in a CeedBasis
380
381 @param name Name of array
382 @param fpformat Printing format
383 @param m Number of rows in array
384 @param n Number of columns in array
385 @param a Array to be viewed
386 @param stream Stream to view to, e.g., stdout
387
388 @return An error code: 0 - success, otherwise - failure
389
390 @ref Utility
391**/
392static int CeedScalarView(const char *name, const char *fpformat, CeedInt m,
393 CeedInt n, const CeedScalar *a, FILE *stream) {
394 for (int i=0; i<m; i++) {
395 if (m > 1)
396 fprintf(stream, "%12s[%d]:", name, i);
397 else
398 fprintf(stream, "%12s:", name);
399 for (int j=0; j<n; j++)
400 fprintf(stream, fpformat, fabs(a[i*n+j]) > 1E-14 ? a[i*n+j] : 0);
401 fputs("\n", stream);
402 }
403 return 0;
404}
405
406/**
407 @brief View a CeedBasis
408
409 @param basis CeedBasis to view
410 @param stream Stream to view to, e.g., stdout
411
412 @return An error code: 0 - success, otherwise - failure
413
414 @ref Utility
415**/
416int CeedBasisView(CeedBasis basis, FILE *stream) {
417 int ierr;
418
419 if (basis->tensorbasis) {
420 fprintf(stream, "CeedBasis: dim=%d P=%d Q=%d\n", basis->dim, basis->P1d,
421 basis->Q1d);
422 ierr = CeedScalarView("qref1d", "\t% 12.8f", 1, basis->Q1d, basis->qref1d,
423 stream); CeedChk(ierr);
424 ierr = CeedScalarView("qweight1d", "\t% 12.8f", 1, basis->Q1d,
425 basis->qweight1d, stream); CeedChk(ierr);
426 ierr = CeedScalarView("interp1d", "\t% 12.8f", basis->Q1d, basis->P1d,
427 basis->interp1d, stream); CeedChk(ierr);
428 ierr = CeedScalarView("grad1d", "\t% 12.8f", basis->Q1d, basis->P1d,
429 basis->grad1d, stream); CeedChk(ierr);
430 } else {
431 fprintf(stream, "CeedBasis: dim=%d P=%d Q=%d\n", basis->dim, basis->P,
432 basis->Q);
433 ierr = CeedScalarView("qref", "\t% 12.8f", 1, basis->Q*basis->dim,
434 basis->qref1d,
435 stream); CeedChk(ierr);
436 ierr = CeedScalarView("qweight", "\t% 12.8f", 1, basis->Q, basis->qweight1d,
437 stream); CeedChk(ierr);
438 ierr = CeedScalarView("interp", "\t% 12.8f", basis->Q, basis->P,
439 basis->interp, stream); CeedChk(ierr);
440 ierr = CeedScalarView("grad", "\t% 12.8f", basis->dim*basis->Q, basis->P,
441 basis->grad, stream); CeedChk(ierr);
442 }
443 return 0;
444}
445
446/**
447 @brief Compute Householder reflection
448
449 Computes A = (I - b v v^T) A
450 where A is an mxn matrix indexed as A[i*row + j*col]
451
452 @param[in,out] A Matrix to apply Householder reflection to, in place
453 @param v Householder vector
454 @param b Scaling factor
455 @param m Number of rows in A
456 @param n Number of columns in A
457 @param row Row stride
458 @param col Col stride
459
460 @return An error code: 0 - success, otherwise - failure
461
462 @ref Developer
463**/
464static int CeedHouseholderReflect(CeedScalar *A, const CeedScalar *v,
465 CeedScalar b, CeedInt m, CeedInt n,
466 CeedInt row, CeedInt col) {
467 for (CeedInt j=0; j<n; j++) {
468 CeedScalar w = A[0*row + j*col];
469 for (CeedInt i=1; i<m; i++)
470 w += v[i] * A[i*row + j*col];
471 A[0*row + j*col] -= b * w;
472 for (CeedInt i=1; i<m; i++)
473 A[i*row + j*col] -= b * w * v[i];
474 }
475 return 0;
476}
477
478/**
479 @brief Apply Householder Q matrix
480
481 Compute A = Q A where Q is mxm and A is mxn.
482
483 @param[in,out] A Matrix to apply Householder Q to, in place
484 @param Q Householder Q matrix
485 @param tau Householder scaling factors
486 @param tmode Transpose mode for application
487 @param m Number of rows in A
488 @param n Number of columns in A
489 @param k Number of elementary reflectors in Q, k<m
490 @param row Row stride in A
491 @param col Col stride in A
492
493 @return An error code: 0 - success, otherwise - failure
494
495 @ref Developer
496**/
497static int CeedHouseholderApplyQ(CeedScalar *A, const CeedScalar *Q,
498 const CeedScalar *tau, CeedTransposeMode tmode,
499 CeedInt m, CeedInt n, CeedInt k,
500 CeedInt row, CeedInt col) {
501 CeedScalar v[m];
502 for (CeedInt ii=0; ii<k; ii++) {
503 CeedInt i = tmode == CEED_TRANSPOSE ? ii : k-1-ii;
504 for (CeedInt j=i+1; j<m; j++)
505 v[j] = Q[j*k+i];
506 // Apply Householder reflector (I - tau v v^T) collograd1d^T
507 CeedHouseholderReflect(&A[i*row], &v[i], tau[i], m-i, n, row, col);
508 }
509 return 0;
510}
511
512/**
513 @brief Compute Givens rotation
514
515 Computes A = G A (or G^T A in transpose mode)
516 where A is an mxn matrix indexed as A[i*n + j*m]
517
518 @param[in,out] A Row major matrix to apply Givens rotation to, in place
519 @param c Cosine factor
520 @param s Sine factor
521 @param i First row/column to apply rotation
522 @param k Second row/column to apply rotation
523 @param m Number of rows in A
524 @param n Number of columns in A
525
526 @return An error code: 0 - success, otherwise - failure
527
528 @ref Developer
529**/
530static int CeedGivensRotation(CeedScalar *A, CeedScalar c, CeedScalar s,
531 CeedTransposeMode tmode, CeedInt i, CeedInt k,
532 CeedInt m, CeedInt n) {
533 CeedInt stridej = 1, strideik = m, numits = n;
534 if (tmode == CEED_NOTRANSPOSE) {
535 stridej = n; strideik = 1; numits = m;
536 }
537
538 // Apply rotation
539 for (CeedInt j=0; j<numits; j++) {
540 CeedScalar tau1 = A[i*strideik+j*stridej], tau2 = A[k*strideik+j*stridej];
541 A[i*strideik+j*stridej] = c*tau1 - s*tau2;
542 A[k*strideik+j*stridej] = s*tau1 + c*tau2;
543 }
544
545 return 0;
546}
547
548/**
549 @brief Return QR Factorization of a matrix
550
551 @param ceed A Ceed context for error handling
552 @param[in,out] mat Row-major matrix to be factorized in place
553 @param[in,out] tau Vector of length m of scaling factors
554 @param m Number of rows
555 @param n Number of columns
556

--- 187 unchanged lines hidden (view full) ---

744 // LCOV_EXCL_START
745 return CeedError(ceed, 1, "Symmetric QR failed to converge");
746 // LCOV_EXCL_STOP
747
748 return 0;
749}
750
751/**
1087 @brief Return QR Factorization of a matrix
1088
1089 @param ceed A Ceed context for error handling
1090 @param[in,out] mat Row-major matrix to be factorized in place
1091 @param[in,out] tau Vector of length m of scaling factors
1092 @param m Number of rows
1093 @param n Number of columns
1094

--- 187 unchanged lines hidden (view full) ---

1282 // LCOV_EXCL_START
1283 return CeedError(ceed, 1, "Symmetric QR failed to converge");
1284 // LCOV_EXCL_STOP
1285
1286 return 0;
1287}
1288
1289/**
752 @brief Return a reference implementation of matrix multiplication C = A B.
753 Note, this is a reference implementation for CPU CeedScalar pointers
754 that is not intended for high performance.
755
756 @param ceed A Ceed context for error handling
757 @param[in] matA Row-major matrix A
758 @param[in] matB Row-major matrix B
759 @param[out] matC Row-major output matrix C
760 @param m Number of rows of C
761 @param n Number of columns of C
762 @param kk Number of columns of A/rows of B
763
764 @return An error code: 0 - success, otherwise - failure
765
766 @ref Utility
767**/
768int CeedMatrixMultiply(Ceed ceed, const CeedScalar *matA,
769 const CeedScalar *matB, CeedScalar *matC, CeedInt m,
770 CeedInt n, CeedInt kk) {
771 for (CeedInt i=0; i<m; i++)
772 for (CeedInt j=0; j<n; j++) {
773 CeedScalar sum = 0;
774 for (CeedInt k=0; k<kk; k++)
775 sum += matA[k+i*kk]*matB[j+k*n];
776 matC[j+i*n] = sum;
777 }
778 return 0;
779}
780
781/**
782 @brief Return Simultaneous Diagonalization of two matrices. This solves the
783 generalized eigenvalue problem A x = lambda B x, where A and B
784 are symmetric and B is positive definite. We generate the matrix X
785 and vector Lambda such that X^T A X = Lambda and X^T B X = I. This
786 is equivalent to the LAPACK routine 'sygv' with TYPE = 1.
787
788 @param ceed A Ceed context for error handling
789 @param[in] matA Row-major matrix to be factorized with eigenvalues

--- 40 unchanged lines hidden (view full) ---

830 // = G D^-1/2 Q
831 ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)matG,
832 (const CeedScalar *)matC, x, n, n, n);
833 CeedChk(ierr);
834
835 return 0;
836}
837
1290 @brief Return Simultaneous Diagonalization of two matrices. This solves the
1291 generalized eigenvalue problem A x = lambda B x, where A and B
1292 are symmetric and B is positive definite. We generate the matrix X
1293 and vector Lambda such that X^T A X = Lambda and X^T B X = I. This
1294 is equivalent to the LAPACK routine 'sygv' with TYPE = 1.
1295
1296 @param ceed A Ceed context for error handling
1297 @param[in] matA Row-major matrix to be factorized with eigenvalues

--- 40 unchanged lines hidden (view full) ---

1338 // = G D^-1/2 Q
1339 ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)matG,
1340 (const CeedScalar *)matC, x, n, n, n);
1341 CeedChk(ierr);
1342
1343 return 0;
1344}
1345
838/**
839 @brief Return collocated grad matrix
840
841 @param basis CeedBasis
842 @param[out] collograd1d Row-major (Q1d * Q1d) matrix expressing derivatives of
843 basis functions at quadrature points
844
845 @return An error code: 0 - success, otherwise - failure
846
847 @ref Advanced
848**/
849int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collograd1d) {
850 int i, j, k;
851 Ceed ceed;
852 CeedInt ierr, P1d=(basis)->P1d, Q1d=(basis)->Q1d;
853 CeedScalar *interp1d, *grad1d, tau[Q1d];
854
855 ierr = CeedMalloc(Q1d*P1d, &interp1d); CeedChk(ierr);
856 ierr = CeedMalloc(Q1d*P1d, &grad1d); CeedChk(ierr);
857 memcpy(interp1d, (basis)->interp1d, Q1d*P1d*sizeof(basis)->interp1d[0]);
858 memcpy(grad1d, (basis)->grad1d, Q1d*P1d*sizeof(basis)->interp1d[0]);
859
860 // QR Factorization, interp1d = Q R
861 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr);
862 ierr = CeedQRFactorization(ceed, interp1d, tau, Q1d, P1d); CeedChk(ierr);
863
864 // Apply Rinv, collograd1d = grad1d Rinv
865 for (i=0; i<Q1d; i++) { // Row i
866 collograd1d[Q1d*i] = grad1d[P1d*i]/interp1d[0];
867 for (j=1; j<P1d; j++) { // Column j
868 collograd1d[j+Q1d*i] = grad1d[j+P1d*i];
869 for (k=0; k<j; k++)
870 collograd1d[j+Q1d*i] -= interp1d[j+P1d*k]*collograd1d[k+Q1d*i];
871 collograd1d[j+Q1d*i] /= interp1d[j+P1d*j];
872 }
873 for (j=P1d; j<Q1d; j++)
874 collograd1d[j+Q1d*i] = 0;
875 }
876
877 // Apply Qtranspose, collograd = collograd Qtranspose
878 CeedHouseholderApplyQ(collograd1d, interp1d, tau, CEED_NOTRANSPOSE,
879 Q1d, Q1d, P1d, 1, Q1d);
880
881 ierr = CeedFree(&interp1d); CeedChk(ierr);
882 ierr = CeedFree(&grad1d); CeedChk(ierr);
883
884 return 0;
885}
886
887/**
888 @brief Apply basis evaluation from nodes to quadrature points or vice versa
889
890 @param basis CeedBasis to evaluate
891 @param nelem The number of elements to apply the basis evaluation to;
892 the backend will specify the ordering in
893 ElemRestrictionCreateBlocked
894 @param tmode \ref CEED_NOTRANSPOSE to evaluate from nodes to quadrature
895 points, \ref CEED_TRANSPOSE to apply the transpose, mapping
896 from quadrature points to nodes
897 @param emode \ref CEED_EVAL_NONE to use values directly,
898 \ref CEED_EVAL_INTERP to use interpolated values,
899 \ref CEED_EVAL_GRAD to use gradients,
900 \ref CEED_EVAL_WEIGHT to use quadrature weights.
901 @param[in] u Input CeedVector
902 @param[out] v Output CeedVector
903
904 @return An error code: 0 - success, otherwise - failure
905
906 @ref Advanced
907**/
908int CeedBasisApply(CeedBasis basis, CeedInt nelem, CeedTransposeMode tmode,
909 CeedEvalMode emode, CeedVector u, CeedVector v) {
910 int ierr;
911 CeedInt ulength = 0, vlength, nnodes, nqpt;
912 if (!basis->Apply)
913 // LCOV_EXCL_START
914 return CeedError(basis->ceed, 1, "Backend does not support BasisApply");
915 // LCOV_EXCL_STOP
916
917 // Check compatibility of topological and geometrical dimensions
918 ierr = CeedBasisGetNumNodes(basis, &nnodes); CeedChk(ierr);
919 ierr = CeedBasisGetNumQuadraturePoints(basis, &nqpt); CeedChk(ierr);
920 ierr = CeedVectorGetLength(v, &vlength); CeedChk(ierr);
921
922 if (u) {
923 ierr = CeedVectorGetLength(u, &ulength); CeedChk(ierr);
924 }
925
926 if ((tmode == CEED_TRANSPOSE && (vlength%nnodes != 0 || ulength%nqpt != 0)) ||
927 (tmode == CEED_NOTRANSPOSE && (ulength%nnodes != 0 || vlength%nqpt != 0)))
928 return CeedError(basis->ceed, 1, "Length of input/output vectors "
929 "incompatible with basis dimensions");
930
931 ierr = basis->Apply(basis, nelem, tmode, emode, u, v); CeedChk(ierr);
932 return 0;
933}
934
935/**
936 @brief Get Ceed associated with a CeedBasis
937
938 @param basis CeedBasis
939 @param[out] ceed Variable to store Ceed
940
941 @return An error code: 0 - success, otherwise - failure
942
943 @ref Advanced
944**/
945int CeedBasisGetCeed(CeedBasis basis, Ceed *ceed) {
946 *ceed = basis->ceed;
947 return 0;
948};
949
950/**
951 @brief Get dimension for given CeedBasis
952
953 @param basis CeedBasis
954 @param[out] dim Variable to store dimension of basis
955
956 @return An error code: 0 - success, otherwise - failure
957
958 @ref Advanced
959**/
960int CeedBasisGetDimension(CeedBasis basis, CeedInt *dim) {
961 *dim = basis->dim;
962 return 0;
963};
964
965/**
966 @brief Get tensor status for given CeedBasis
967
968 @param basis CeedBasis
969 @param[out] tensor Variable to store tensor status
970
971 @return An error code: 0 - success, otherwise - failure
972
973 @ref Advanced
974**/
975int CeedBasisGetTensorStatus(CeedBasis basis, bool *tensor) {
976 *tensor = basis->tensorbasis;
977 return 0;
978};
979
980/**
981 @brief Get number of components for given CeedBasis
982
983 @param basis CeedBasis
984 @param[out] numcomp Variable to store number of components of basis
985
986 @return An error code: 0 - success, otherwise - failure
987
988 @ref Advanced
989**/
990int CeedBasisGetNumComponents(CeedBasis basis, CeedInt *numcomp) {
991 *numcomp = basis->ncomp;
992 return 0;
993};
994
995/**
996 @brief Get total number of nodes (in 1 dimension) of a CeedBasis
997
998 @param basis CeedBasis
999 @param[out] P1d Variable to store number of nodes
1000
1001 @return An error code: 0 - success, otherwise - failure
1002
1003 @ref Advanced
1004**/
1005int CeedBasisGetNumNodes1D(CeedBasis basis, CeedInt *P1d) {
1006 if (!basis->tensorbasis)
1007 // LCOV_EXCL_START
1008 return CeedError(basis->ceed, 1, "Cannot supply P1d for non-tensor basis");
1009 // LCOV_EXCL_STOP
1010
1011 *P1d = basis->P1d;
1012 return 0;
1013}
1014
1015/**
1016 @brief Get total number of quadrature points (in 1 dimension) of a CeedBasis
1017
1018 @param basis CeedBasis
1019 @param[out] Q1d Variable to store number of quadrature points
1020
1021 @return An error code: 0 - success, otherwise - failure
1022
1023 @ref Advanced
1024**/
1025int CeedBasisGetNumQuadraturePoints1D(CeedBasis basis, CeedInt *Q1d) {
1026 if (!basis->tensorbasis)
1027 // LCOV_EXCL_START
1028 return CeedError(basis->ceed, 1, "Cannot supply Q1d for non-tensor basis");
1029 // LCOV_EXCL_STOP
1030
1031 *Q1d = basis->Q1d;
1032 return 0;
1033}
1034
1035/**
1036 @brief Get total number of nodes (in dim dimensions) of a CeedBasis
1037
1038 @param basis CeedBasis
1039 @param[out] P Variable to store number of nodes
1040
1041 @return An error code: 0 - success, otherwise - failure
1042
1043 @ref Utility
1044**/
1045int CeedBasisGetNumNodes(CeedBasis basis, CeedInt *P) {
1046 *P = basis->P;
1047 return 0;
1048}
1049
1050/**
1051 @brief Get total number of quadrature points (in dim dimensions) of a CeedBasis
1052
1053 @param basis CeedBasis
1054 @param[out] Q Variable to store number of quadrature points
1055
1056 @return An error code: 0 - success, otherwise - failure
1057
1058 @ref Utility
1059**/
1060int CeedBasisGetNumQuadraturePoints(CeedBasis basis, CeedInt *Q) {
1061 *Q = basis->Q;
1062 return 0;
1063}
1064
1065/**
1066 @brief Get reference coordinates of quadrature points (in dim dimensions)
1067 of a CeedBasis
1068
1069 @param basis CeedBasis
1070 @param[out] qref Variable to store reference coordinates of quadrature points
1071
1072 @return An error code: 0 - success, otherwise - failure
1073
1074 @ref Advanced
1075**/
1076int CeedBasisGetQRef(CeedBasis basis, CeedScalar **qref) {
1077 *qref = basis->qref1d;
1078 return 0;
1079}
1080
1081/**
1082 @brief Get quadrature weights of quadrature points (in dim dimensions)
1083 of a CeedBasis
1084
1085 @param basis CeedBasis
1086 @param[out] qweight Variable to store quadrature weights
1087
1088 @return An error code: 0 - success, otherwise - failure
1089
1090 @ref Advanced
1091**/
1092int CeedBasisGetQWeights(CeedBasis basis, CeedScalar **qweight) {
1093 *qweight = basis->qweight1d;
1094 return 0;
1095}
1096
1097/**
1098 @brief Get interpolation matrix of a CeedBasis
1099
1100 @param basis CeedBasis
1101 @param[out] interp Variable to store interpolation matrix
1102
1103 @return An error code: 0 - success, otherwise - failure
1104
1105 @ref Advanced
1106**/
1107int CeedBasisGetInterp(CeedBasis basis, CeedScalar **interp) {
1108 if (!basis->interp && basis->tensorbasis) {
1109 // Allocate
1110 int ierr;
1111 ierr = CeedMalloc(basis->Q*basis->P, &basis->interp); CeedChk(ierr);
1112
1113 // Initialize
1114 for (CeedInt i=0; i<basis->Q*basis->P; i++)
1115 basis->interp[i] = 1.0;
1116
1117 // Calculate
1118 for (CeedInt d=0; d<basis->dim; d++)
1119 for (CeedInt qpt=0; qpt<basis->Q; qpt++)
1120 for (CeedInt node=0; node<basis->P; node++) {
1121 CeedInt p = (node / CeedIntPow(basis->P1d, d)) % basis->P1d;
1122 CeedInt q = (qpt / CeedIntPow(basis->Q1d, d)) % basis->Q1d;
1123 basis->interp[qpt*(basis->P)+node] *= basis->interp1d[q*basis->P1d+p];
1124 }
1125 }
1126
1127 *interp = basis->interp;
1128
1129 return 0;
1130}
1131
1132/**
1133 @brief Get 1D interpolation matrix of a tensor product CeedBasis
1134
1135 @param basis CeedBasis
1136 @param[out] interp1d Variable to store interpolation matrix
1137
1138 @return An error code: 0 - success, otherwise - failure
1139
1140 @ref Advanced
1141**/
1142int CeedBasisGetInterp1D(CeedBasis basis, CeedScalar **interp1d) {
1143 if (!basis->tensorbasis)
1144 // LCOV_EXCL_START
1145 return CeedError(basis->ceed, 1, "CeedBasis is not a tensor product basis.");
1146 // LCOV_EXCL_STOP
1147
1148 *interp1d = basis->interp1d;
1149
1150 return 0;
1151}
1152
1153/**
1154 @brief Get gradient matrix of a CeedBasis
1155
1156 @param basis CeedBasis
1157 @param[out] grad Variable to store gradient matrix
1158
1159 @return An error code: 0 - success, otherwise - failure
1160
1161 @ref Advanced
1162**/
1163int CeedBasisGetGrad(CeedBasis basis, CeedScalar **grad) {
1164 if (!basis->grad && basis->tensorbasis) {
1165 // Allocate
1166 int ierr;
1167 ierr = CeedMalloc(basis->dim*basis->Q*basis->P, &basis->grad);
1168 CeedChk(ierr);
1169
1170 // Initialize
1171 for (CeedInt i=0; i<basis->dim*basis->Q*basis->P; i++)
1172 basis->grad[i] = 1.0;
1173
1174 // Calculate
1175 for (CeedInt d=0; d<basis->dim; d++)
1176 for (CeedInt i=0; i<basis->dim; i++)
1177 for (CeedInt qpt=0; qpt<basis->Q; qpt++)
1178 for (CeedInt node=0; node<basis->P; node++) {
1179 CeedInt p = (node / CeedIntPow(basis->P1d, d)) % basis->P1d;
1180 CeedInt q = (qpt / CeedIntPow(basis->Q1d, d)) % basis->Q1d;
1181 if (i == d)
1182 basis->grad[(i*basis->Q+qpt)*(basis->P)+node] *=
1183 basis->grad1d[q*basis->P1d+p];
1184 else
1185 basis->grad[(i*basis->Q+qpt)*(basis->P)+node] *=
1186 basis->interp1d[q*basis->P1d+p];
1187 }
1188 }
1189
1190 *grad = basis->grad;
1191
1192 return 0;
1193}
1194
1195/**
1196 @brief Get 1D gradient matrix of a tensor product CeedBasis
1197
1198 @param basis CeedBasis
1199 @param[out] grad1d Variable to store gradient matrix
1200
1201 @return An error code: 0 - success, otherwise - failure
1202
1203 @ref Advanced
1204**/
1205int CeedBasisGetGrad1D(CeedBasis basis, CeedScalar **grad1d) {
1206 if (!basis->tensorbasis)
1207 // LCOV_EXCL_START
1208 return CeedError(basis->ceed, 1, "CeedBasis is not a tensor product basis.");
1209 // LCOV_EXCL_STOP
1210
1211 *grad1d = basis->grad1d;
1212
1213 return 0;
1214}
1215
1216/**
1217 @brief Get backend data of a CeedBasis
1218
1219 @param basis CeedBasis
1220 @param[out] data Variable to store data
1221
1222 @return An error code: 0 - success, otherwise - failure
1223
1224 @ref Advanced
1225**/
1226int CeedBasisGetData(CeedBasis basis, void **data) {
1227 *data = basis->data;
1228 return 0;
1229}
1230
1231/**
1232 @brief Set backend data of a CeedBasis
1233
1234 @param[out] basis CeedBasis
1235 @param data Data to set
1236
1237 @return An error code: 0 - success, otherwise - failure
1238
1239 @ref Advanced
1240**/
1241int CeedBasisSetData(CeedBasis basis, void **data) {
1242 basis->data = *data;
1243 return 0;
1244}
1245
1246/**
1247 @brief Get CeedTensorContract of a CeedBasis
1248
1249 @param basis CeedBasis
1250 @param[out] contract Variable to store CeedTensorContract
1251
1252 @return An error code: 0 - success, otherwise - failure
1253
1254 @ref Advanced
1255**/
1256int CeedBasisGetTensorContract(CeedBasis basis, CeedTensorContract *contract) {
1257 *contract = basis->contract;
1258 return 0;
1259}
1260
1261/**
1262 @brief Set CeedTensorContract of a CeedBasis
1263
1264 @param[out] basis CeedBasis
1265 @param contract CeedTensorContract to set
1266
1267 @return An error code: 0 - success, otherwise - failure
1268
1269 @ref Advanced
1270**/
1271int CeedBasisSetTensorContract(CeedBasis basis, CeedTensorContract *contract) {
1272 basis->contract = *contract;
1273 return 0;
1274}
1275
1276/**
1277 @brief Get dimension for given CeedElemTopology
1278
1279 @param topo CeedElemTopology
1280 @param[out] dim Variable to store dimension of topology
1281
1282 @return An error code: 0 - success, otherwise - failure
1283
1284 @ref Advanced
1285**/
1286int CeedBasisGetTopologyDimension(CeedElemTopology topo, CeedInt *dim) {
1287 *dim = (CeedInt) topo >> 16;
1288 return 0;
1289};
1290
1291/**
1292 @brief Destroy a CeedBasis
1293
1294 @param basis CeedBasis to destroy
1295
1296 @return An error code: 0 - success, otherwise - failure
1297
1298 @ref Basic
1299**/
1300int CeedBasisDestroy(CeedBasis *basis) {
1301 int ierr;
1302
1303 if (!*basis || --(*basis)->refcount > 0)
1304 return 0;
1305 if ((*basis)->Destroy) {
1306 ierr = (*basis)->Destroy(*basis); CeedChk(ierr);
1307 }
1308 ierr = CeedFree(&(*basis)->interp); CeedChk(ierr);
1309 ierr = CeedFree(&(*basis)->interp1d); CeedChk(ierr);
1310 ierr = CeedFree(&(*basis)->grad); CeedChk(ierr);
1311 ierr = CeedFree(&(*basis)->grad1d); CeedChk(ierr);
1312 ierr = CeedFree(&(*basis)->qref1d); CeedChk(ierr);
1313 ierr = CeedFree(&(*basis)->qweight1d); CeedChk(ierr);
1314 ierr = CeedDestroy(&(*basis)->ceed); CeedChk(ierr);
1315 ierr = CeedFree(basis); CeedChk(ierr);
1316 return 0;
1317}
1318
1319/// @cond DOXYGEN_SKIP
1320// Indicate that the quadrature points are collocated with the nodes
1321const CeedBasis CEED_BASIS_COLLOCATED = &ceed_basis_collocated;
1322/// @endcond
1323/// @}
1346/// @}