ceed-ref-basis.c (d7b241e67f6e33d9b297db3da3be4f167f32bbee) ceed-ref-basis.c (b5cf12ee8b5ffe0ca173fbf1c0d42d7d14238b45)
1// Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC.
2// Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707.
3// All Rights 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.

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

47}
48
49static int CeedBasisApply_Ref(CeedBasis basis, CeedInt nelem,
50 CeedTransposeMode tmode, CeedEvalMode emode,
51 const CeedScalar *u, CeedScalar *v) {
52 int ierr;
53 const CeedInt dim = basis->dim;
54 const CeedInt ncomp = basis->ncomp;
1// Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC.
2// Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707.
3// All Rights 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.

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

47}
48
49static int CeedBasisApply_Ref(CeedBasis basis, CeedInt nelem,
50 CeedTransposeMode tmode, CeedEvalMode emode,
51 const CeedScalar *u, CeedScalar *v) {
52 int ierr;
53 const CeedInt dim = basis->dim;
54 const CeedInt ncomp = basis->ncomp;
55 const CeedInt nqpt = CeedPowInt(basis->Q1d, dim);
55 const CeedInt nqpt = CeedIntPow(basis->Q1d, dim);
56 const CeedInt add = (tmode == CEED_TRANSPOSE);
57
58 if (nelem != 1)
59 return CeedError(basis->ceed, 1,
60 "This backend does not support BasisApply for multiple elements");
61
62 // Clear v if operating in transpose
63 if (tmode == CEED_TRANSPOSE) {
56 const CeedInt add = (tmode == CEED_TRANSPOSE);
57
58 if (nelem != 1)
59 return CeedError(basis->ceed, 1,
60 "This backend does not support BasisApply for multiple elements");
61
62 // Clear v if operating in transpose
63 if (tmode == CEED_TRANSPOSE) {
64 const CeedInt vsize = ncomp*CeedPowInt(basis->P1d, dim);
64 const CeedInt vsize = ncomp*CeedIntPow(basis->P1d, dim);
65 for (CeedInt i = 0; i < vsize; i++)
66 v[i] = (CeedScalar) 0;
67 }
68 switch (emode) {
69 // Interpolate to/from quadrature points
70 case CEED_EVAL_INTERP: {
71 CeedInt P = basis->P1d, Q = basis->Q1d;
72 if (tmode == CEED_TRANSPOSE) {
73 P = basis->Q1d; Q = basis->P1d;
74 }
65 for (CeedInt i = 0; i < vsize; i++)
66 v[i] = (CeedScalar) 0;
67 }
68 switch (emode) {
69 // Interpolate to/from quadrature points
70 case CEED_EVAL_INTERP: {
71 CeedInt P = basis->P1d, Q = basis->Q1d;
72 if (tmode == CEED_TRANSPOSE) {
73 P = basis->Q1d; Q = basis->P1d;
74 }
75 CeedInt pre = ncomp*CeedPowInt(P, dim-1), post = 1;
76 CeedScalar tmp[2][ncomp*Q*CeedPowInt(P>Q?P:Q, dim-1)];
75 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = 1;
76 CeedScalar tmp[2][ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)];
77 for (CeedInt d=0; d<dim; d++) {
78 ierr = CeedTensorContract_Ref(basis->ceed, pre, P, post, Q, basis->interp1d,
79 tmode, add&&(d==dim-1),
80 d==0?u:tmp[d%2], d==dim-1?v:tmp[(d+1)%2]);
81 CeedChk(ierr);
82 pre /= P;
83 post *= Q;
84 }
85 } break;
86 // Evaluate the gradient to/from quadrature points
87 case CEED_EVAL_GRAD: {
88 CeedInt P = basis->P1d, Q = basis->Q1d;
89 // In CEED_NOTRANSPOSE mode:
90 // u is [dim, ncomp, P^dim, nelem], row-major layout
91 // v is [dim, ncomp, Q^dim, nelem], row-major layout
92 // In CEED_TRANSPOSE mode, the sizes of u and v are switched.
93 if (tmode == CEED_TRANSPOSE) {
94 P = basis->Q1d, Q = basis->P1d;
95 }
77 for (CeedInt d=0; d<dim; d++) {
78 ierr = CeedTensorContract_Ref(basis->ceed, pre, P, post, Q, basis->interp1d,
79 tmode, add&&(d==dim-1),
80 d==0?u:tmp[d%2], d==dim-1?v:tmp[(d+1)%2]);
81 CeedChk(ierr);
82 pre /= P;
83 post *= Q;
84 }
85 } break;
86 // Evaluate the gradient to/from quadrature points
87 case CEED_EVAL_GRAD: {
88 CeedInt P = basis->P1d, Q = basis->Q1d;
89 // In CEED_NOTRANSPOSE mode:
90 // u is [dim, ncomp, P^dim, nelem], row-major layout
91 // v is [dim, ncomp, Q^dim, nelem], row-major layout
92 // In CEED_TRANSPOSE mode, the sizes of u and v are switched.
93 if (tmode == CEED_TRANSPOSE) {
94 P = basis->Q1d, Q = basis->P1d;
95 }
96 CeedScalar tmp[2][ncomp*Q*CeedPowInt(P>Q?P:Q, dim-1)];
96 CeedScalar tmp[2][ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)];
97 for (CeedInt p = 0; p < dim; p++) {
97 for (CeedInt p = 0; p < dim; p++) {
98 CeedInt pre = ncomp*CeedPowInt(P, dim-1), post = 1;
98 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = 1;
99 for (CeedInt d=0; d<dim; d++) {
100 ierr = CeedTensorContract_Ref(basis->ceed, pre, P, post, Q,
101 (p==d)?basis->grad1d:basis->interp1d,
102 tmode, add&&(d==dim-1),
103 d==0
104 ? (tmode==CEED_NOTRANSPOSE?u:u+p*ncomp*nqpt)
105 : tmp[d%2],
106 d==dim-1

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

114 } break;
115 // Retrieve interpolation weights
116 case CEED_EVAL_WEIGHT: {
117 if (tmode == CEED_TRANSPOSE)
118 return CeedError(basis->ceed, 1,
119 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE");
120 CeedInt Q = basis->Q1d;
121 for (CeedInt d=0; d<dim; d++) {
99 for (CeedInt d=0; d<dim; d++) {
100 ierr = CeedTensorContract_Ref(basis->ceed, pre, P, post, Q,
101 (p==d)?basis->grad1d:basis->interp1d,
102 tmode, add&&(d==dim-1),
103 d==0
104 ? (tmode==CEED_NOTRANSPOSE?u:u+p*ncomp*nqpt)
105 : tmp[d%2],
106 d==dim-1

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

114 } break;
115 // Retrieve interpolation weights
116 case CEED_EVAL_WEIGHT: {
117 if (tmode == CEED_TRANSPOSE)
118 return CeedError(basis->ceed, 1,
119 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE");
120 CeedInt Q = basis->Q1d;
121 for (CeedInt d=0; d<dim; d++) {
122 CeedInt pre = CeedPowInt(Q, dim-d-1), post = CeedPowInt(Q, d);
122 CeedInt pre = CeedIntPow(Q, dim-d-1), post = CeedIntPow(Q, d);
123 for (CeedInt i=0; i<pre; i++)
124 for (CeedInt j=0; j<Q; j++)
125 for (CeedInt k=0; k<post; k++)
126 v[(i*Q + j)*post + k] = basis->qweight1d[j]
127 * (d == 0 ? 1 : v[(i*Q + j)*post + k]);
128 }
129 } break;
130 // Evaluate the divergence to/from the quadrature points

--- 26 unchanged lines hidden --- 		123 for (CeedInt i=0; i<pre; i++)
124 for (CeedInt j=0; j<Q; j++)
125 for (CeedInt k=0; k<post; k++)
126 v[(i*Q + j)*post + k] = basis->qweight1d[j]
127 * (d == 0 ? 1 : v[(i*Q + j)*post + k]);
128 }
129 } break;
130 // Evaluate the divergence to/from the quadrature points

--- 26 unchanged lines hidden ---