ceed-ref-basis.c (dc1dbf0744997d47eb98656f475bfd131b09dcf8) ceed-ref-basis.c (84a01de5ce080ac9cdd243d9d64da2df0ae9cb77)
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.

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

16
17#include <string.h>
18#include "ceed-ref.h"
19
20// Contracts on the middle index
21// NOTRANSPOSE: V_ajc = T_jb U_abc
22// TRANSPOSE: V_ajc = T_bj U_abc
23// If Add != 0, "=" is replaced by "+="
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.

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

16
17#include <string.h>
18#include "ceed-ref.h"
19
20// Contracts on the middle index
21// NOTRANSPOSE: V_ajc = T_jb U_abc
22// TRANSPOSE: V_ajc = T_bj U_abc
23// If Add != 0, "=" is replaced by "+="
24static int CeedTensorContract_Ref(Ceed ceed, CeedInt A, CeedInt B, CeedInt C,
25 CeedInt J,
26 const CeedScalar *restrict t, CeedTransposeMode tmode,
27 const CeedInt Add,
28 const CeedScalar *restrict u, CeedScalar *restrict v) {
24static int CeedTensorContract_Ref(Ceed ceed, CeedInt A, CeedInt B,
25 CeedInt C, CeedInt J,
26 const CeedScalar *restrict t,
27 CeedTransposeMode tmode, const CeedInt Add,
28 const CeedScalar *restrict u,
29 CeedScalar *restrict v) {
29 CeedInt tstride0 = B, tstride1 = 1;
30 if (tmode == CEED_TRANSPOSE) {
31 tstride0 = 1; tstride1 = J;
32 }
33
34 if (!Add)
35 for (CeedInt q=0; q<A*J*C; q++)
36 v[q] = (CeedScalar) 0.0;

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

62 if (U) {
63 ierr = CeedVectorGetArrayRead(U, CEED_MEM_HOST, &u); CeedChk(ierr);
64 } else if (emode != CEED_EVAL_WEIGHT) {
65 return CeedError(ceed, 1,
66 "An input vector is required for this CeedEvalMode");
67 }
68 ierr = CeedVectorGetArray(V, CEED_MEM_HOST, &v); CeedChk(ierr);
69
30 CeedInt tstride0 = B, tstride1 = 1;
31 if (tmode == CEED_TRANSPOSE) {
32 tstride0 = 1; tstride1 = J;
33 }
34
35 if (!Add)
36 for (CeedInt q=0; q<A*J*C; q++)
37 v[q] = (CeedScalar) 0.0;

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

63 if (U) {
64 ierr = CeedVectorGetArrayRead(U, CEED_MEM_HOST, &u); CeedChk(ierr);
65 } else if (emode != CEED_EVAL_WEIGHT) {
66 return CeedError(ceed, 1,
67 "An input vector is required for this CeedEvalMode");
68 }
69 ierr = CeedVectorGetArray(V, CEED_MEM_HOST, &v); CeedChk(ierr);
70
70 if (nelem != 1)
71 return CeedError(ceed, 1,
72 "This backend does not support BasisApply for multiple elements");
73
74 // Clear v if operating in transpose
75 if (tmode == CEED_TRANSPOSE) {
71 // Clear v if operating in transpose
72 if (tmode == CEED_TRANSPOSE) {
76 const CeedInt vsize = ncomp*ndof;
73 const CeedInt vsize = nelem*ncomp*ndof;
77 for (CeedInt i = 0; i < vsize; i++)
74 for (CeedInt i = 0; i < vsize; i++)
78 v[i] = (CeedScalar) 0;
75 v[i] = (CeedScalar) 0.0;
79 }
76 }
80 // Tensor basis
81 bool tensorbasis;
82 ierr = CeedBasisGetTensorStatus(basis, &tensorbasis); CeedChk(ierr);
77 bool tensorbasis;
78 ierr = CeedBasisGetTensorStatus(basis, &tensorbasis); CeedChk(ierr);
79 // Tensor basis
83 if (tensorbasis) {
84 CeedInt P1d, Q1d;
85 ierr = CeedBasisGetNumNodes1D(basis, &P1d); CeedChk(ierr);
86 ierr = CeedBasisGetNumQuadraturePoints1D(basis, &Q1d); CeedChk(ierr);
87 switch (emode) {
88 // Interpolate to/from quadrature points
89 case CEED_EVAL_INTERP: {
90 CeedInt P = P1d, Q = Q1d;
91 if (tmode == CEED_TRANSPOSE) {
92 P = Q1d; Q = P1d;
93 }
80 if (tensorbasis) {
81 CeedInt P1d, Q1d;
82 ierr = CeedBasisGetNumNodes1D(basis, &P1d); CeedChk(ierr);
83 ierr = CeedBasisGetNumQuadraturePoints1D(basis, &Q1d); CeedChk(ierr);
84 switch (emode) {
85 // Interpolate to/from quadrature points
86 case CEED_EVAL_INTERP: {
87 CeedInt P = P1d, Q = Q1d;
88 if (tmode == CEED_TRANSPOSE) {
89 P = Q1d; Q = P1d;
90 }
94 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = 1;
95 CeedScalar tmp[2][ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)];
91 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = nelem;
92 CeedScalar tmp[2][nelem*ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)];
96 CeedScalar *interp1d;
97 ierr = CeedBasisGetInterp(basis, &interp1d); CeedChk(ierr);
98 for (CeedInt d=0; d<dim; d++) {
93 CeedScalar *interp1d;
94 ierr = CeedBasisGetInterp(basis, &interp1d); CeedChk(ierr);
95 for (CeedInt d=0; d<dim; d++) {
99 ierr = CeedTensorContract_Ref(ceed, pre, P, post, Q, interp1d,
100 tmode, add&&(d==dim-1),
96 ierr = CeedTensorContract_Ref(ceed, pre, P, post, Q,
97 interp1d, tmode, add&&(d==dim-1),
101 d==0?u:tmp[d%2], d==dim-1?v:tmp[(d+1)%2]);
102 CeedChk(ierr);
103 pre /= P;
104 post *= Q;
105 }
106 } break;
107 // Evaluate the gradient to/from quadrature points
108 case CEED_EVAL_GRAD: {
98 d==0?u:tmp[d%2], d==dim-1?v:tmp[(d+1)%2]);
99 CeedChk(ierr);
100 pre /= P;
101 post *= Q;
102 }
103 } break;
104 // Evaluate the gradient to/from quadrature points
105 case CEED_EVAL_GRAD: {
109 CeedInt P = P1d, Q = Q1d;
110 // In CEED_NOTRANSPOSE mode:
111 // u has shape [dim, ncomp, P^dim, nelem], row-major layout
112 // v has shape [dim, ncomp, Q^dim, nelem], row-major layout
113 // In CEED_TRANSPOSE mode, the sizes of u and v are switched.
106 // In CEED_NOTRANSPOSE mode:
107 // u has shape [dim, ncomp, P^dim, nelem], row-major layout
108 // v has shape [dim, ncomp, Q^dim, nelem], row-major layout
109 // In CEED_TRANSPOSE mode, the sizes of u and v are switched.
110 CeedInt P = P1d, Q = Q1d;
114 if (tmode == CEED_TRANSPOSE) {
111 if (tmode == CEED_TRANSPOSE) {
115 P = Q1d, Q = P1d;
112 P = Q1d, Q = Q1d;
116 }
113 }
117 CeedScalar tmp[2][ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)];
118 CeedScalar *interp1d, *grad1d;
114 CeedBasis_Ref *impl;
115 ierr = CeedBasisGetData(basis, (void *)&impl); CeedChk(ierr);
116 CeedScalar interp[nelem*ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)];
117 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = nelem;
118 CeedScalar tmp[2][nelem*ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)];
119 CeedScalar *interp1d;
119 ierr = CeedBasisGetInterp(basis, &interp1d); CeedChk(ierr);
120 ierr = CeedBasisGetInterp(basis, &interp1d); CeedChk(ierr);
120 ierr = CeedBasisGetGrad(basis, &grad1d); CeedChk(ierr);
121 for (CeedInt p = 0; p < dim; p++) {
122 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = 1;
123 for (CeedInt d=0; d<dim; d++) {
124 ierr = CeedTensorContract_Ref(ceed, pre, P, post, Q,
125 (p==d)?grad1d:interp1d,
126 tmode, add&&(d==dim-1),
127 (d == 0
128 ? (tmode==CEED_NOTRANSPOSE?u:u+p*ncomp*nqpt)
129 : tmp[d%2]),
130 (d == dim-1
131 ? (tmode==CEED_TRANSPOSE?v:v+p*ncomp*nqpt)
132 : tmp[(d+1)%2]));
133 CeedChk(ierr);
134 pre /= P;
135 post *= Q;
136 }
121 // Interpolate to quadrature points (NoTranspose)
122 // or Grad to quadrature points (Transpose)
123 for (CeedInt d=0; d<dim; d++) {
124 ierr = CeedTensorContract_Ref(ceed, pre, P, post, Q,
125 (tmode == CEED_NOTRANSPOSE
126 ? interp1d
127 : impl->colograd1d),
128 tmode, add&&(d>0),
129 (tmode == CEED_NOTRANSPOSE
130 ? (d==0?u:tmp[d%2])
131 : u + d*nqpt*ncomp*nelem),
132 (tmode == CEED_NOTRANSPOSE
133 ? (d==dim-1?interp:tmp[(d+1)%2])
134 : interp));
135 CeedChk(ierr);
136 pre /= P;
137 post *= Q;
137 }
138 }
139 // Grad to quadrature points (NoTranspose)
140 // or Interpolate to dofs (Transpose)
141 P = Q1d, Q = Q1d;
142 if (tmode == CEED_TRANSPOSE) {
143 P = Q1d, Q = P1d;
144 }
145 pre = ncomp*CeedIntPow(P, dim-1), post = nelem;
146 for (CeedInt d=0; d<dim; d++) {
147 ierr = CeedTensorContract_Ref(ceed, pre, P, post, Q,
148 (tmode == CEED_NOTRANSPOSE
149 ? impl->colograd1d
150 : interp1d),
151 tmode, add&&(d==dim-1),
152 (tmode == CEED_NOTRANSPOSE
153 ? interp
154 : (d==0?interp:tmp[d%2])),
155 (tmode == CEED_NOTRANSPOSE
156 ? v + d*nqpt*ncomp*nelem
157 : (d==dim-1?v:tmp[(d+1)%2])));
158 CeedChk(ierr);
159 pre /= P;
160 post *= Q;
161 }
138 } break;
139 // Retrieve interpolation weights
140 case CEED_EVAL_WEIGHT: {
141 if (tmode == CEED_TRANSPOSE)
142 return CeedError(ceed, 1,
143 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE");
144 CeedInt Q = Q1d;
145 CeedScalar *qweight1d;
146 ierr = CeedBasisGetQWeights(basis, &qweight1d); CeedChk(ierr);
147 for (CeedInt d=0; d<dim; d++) {
148 CeedInt pre = CeedIntPow(Q, dim-d-1), post = CeedIntPow(Q, d);
149 for (CeedInt i=0; i<pre; i++)
150 for (CeedInt j=0; j<Q; j++)
162 } break;
163 // Retrieve interpolation weights
164 case CEED_EVAL_WEIGHT: {
165 if (tmode == CEED_TRANSPOSE)
166 return CeedError(ceed, 1,
167 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE");
168 CeedInt Q = Q1d;
169 CeedScalar *qweight1d;
170 ierr = CeedBasisGetQWeights(basis, &qweight1d); CeedChk(ierr);
171 for (CeedInt d=0; d<dim; d++) {
172 CeedInt pre = CeedIntPow(Q, dim-d-1), post = CeedIntPow(Q, d);
173 for (CeedInt i=0; i<pre; i++)
174 for (CeedInt j=0; j<Q; j++)
151 for (CeedInt k=0; k<post; k++)
152 v[(i*Q + j)*post + k] = qweight1d[j]
153 * (d == 0 ? 1 : v[(i*Q + j)*post + k]);
175 for (CeedInt k=0; k<post; k++) {
176 CeedScalar w = qweight1d[j]
177 * (d == 0 ? 1 : v[((i*Q + j)*post + k)*nelem]);
178 for (CeedInt e=0; e<nelem; e++)
179 v[((i*Q + j)*post + k)*nelem + e] = w;
180 }
154 }
155 } break;
156 // Evaluate the divergence to/from the quadrature points
157 case CEED_EVAL_DIV:
158 return CeedError(ceed, 1, "CEED_EVAL_DIV not supported");
159 // Evaluate the curl to/from the quadrature points
160 case CEED_EVAL_CURL:
161 return CeedError(ceed, 1, "CEED_EVAL_CURL not supported");
162 // Take no action, BasisApply should not have been called
163 case CEED_EVAL_NONE:
164 return CeedError(ceed, 1,
165 "CEED_EVAL_NONE does not make sense in this context");
166 }
167 } else {
168 // Non-tensor basis
169 switch (emode) {
181 }
182 } break;
183 // Evaluate the divergence to/from the quadrature points
184 case CEED_EVAL_DIV:
185 return CeedError(ceed, 1, "CEED_EVAL_DIV not supported");
186 // Evaluate the curl to/from the quadrature points
187 case CEED_EVAL_CURL:
188 return CeedError(ceed, 1, "CEED_EVAL_CURL not supported");
189 // Take no action, BasisApply should not have been called
190 case CEED_EVAL_NONE:
191 return CeedError(ceed, 1,
192 "CEED_EVAL_NONE does not make sense in this context");
193 }
194 } else {
195 // Non-tensor basis
196 switch (emode) {
197 // Interpolate to/from quadrature points
170 case CEED_EVAL_INTERP: {
171 CeedInt P = ndof, Q = nqpt;
172 CeedScalar *interp;
173 ierr = CeedBasisGetInterp(basis, &interp); CeedChk(ierr);
174 if (tmode == CEED_TRANSPOSE) {
175 P = nqpt; Q = ndof;
176 }
198 case CEED_EVAL_INTERP: {
199 CeedInt P = ndof, Q = nqpt;
200 CeedScalar *interp;
201 ierr = CeedBasisGetInterp(basis, &interp); CeedChk(ierr);
202 if (tmode == CEED_TRANSPOSE) {
203 P = nqpt; Q = ndof;
204 }
177 ierr = CeedTensorContract_Ref(ceed, ncomp, P, 1, Q, interp,
178 tmode, add, u, v);
205 ierr = CeedTensorContract_Ref(ceed, ncomp, P, nelem, Q,
206 interp, tmode, add, u, v);
179 CeedChk(ierr);
180 }
181 break;
207 CeedChk(ierr);
208 }
209 break;
210 // Evaluate the gradient to/from quadrature points
182 case CEED_EVAL_GRAD: {
183 CeedInt P = ndof, Q = dim*nqpt;
184 CeedScalar *grad;
185 ierr = CeedBasisGetGrad(basis, &grad); CeedChk(ierr);
186 if (tmode == CEED_TRANSPOSE) {
187 P = dim*nqpt; Q = ndof;
188 }
211 case CEED_EVAL_GRAD: {
212 CeedInt P = ndof, Q = dim*nqpt;
213 CeedScalar *grad;
214 ierr = CeedBasisGetGrad(basis, &grad); CeedChk(ierr);
215 if (tmode == CEED_TRANSPOSE) {
216 P = dim*nqpt; Q = ndof;
217 }
189 ierr = CeedTensorContract_Ref(ceed, ncomp, P, 1, Q, grad,
190 tmode, add, u, v);
218 ierr = CeedTensorContract_Ref(ceed, ncomp, P, nelem, Q,
219 grad, tmode, add, u, v);
191 CeedChk(ierr);
192 }
193 break;
220 CeedChk(ierr);
221 }
222 break;
223 // Retrieve interpolation weights
194 case CEED_EVAL_WEIGHT: {
195 if (tmode == CEED_TRANSPOSE)
196 return CeedError(ceed, 1,
197 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE");
198 CeedScalar *qweight;
199 ierr = CeedBasisGetQWeights(basis, &qweight); CeedChk(ierr);
200 for (CeedInt i=0; i<nqpt; i++)
224 case CEED_EVAL_WEIGHT: {
225 if (tmode == CEED_TRANSPOSE)
226 return CeedError(ceed, 1,
227 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE");
228 CeedScalar *qweight;
229 ierr = CeedBasisGetQWeights(basis, &qweight); CeedChk(ierr);
230 for (CeedInt i=0; i<nqpt; i++)
201 v[i] = qweight[i];
231 for (CeedInt e=0; e<nelem; e++)
232 v[i*nelem + e] = qweight[i];
202 } break;
233 } break;
234 // Evaluate the divergence to/from the quadrature points
203 case CEED_EVAL_DIV:
204 return CeedError(ceed, 1, "CEED_EVAL_DIV not supported");
235 case CEED_EVAL_DIV:
236 return CeedError(ceed, 1, "CEED_EVAL_DIV not supported");
237 // Evaluate the curl to/from the quadrature points
205 case CEED_EVAL_CURL:
206 return CeedError(ceed, 1, "CEED_EVAL_CURL not supported");
238 case CEED_EVAL_CURL:
239 return CeedError(ceed, 1, "CEED_EVAL_CURL not supported");
240 // Take no action, BasisApply should not have been called
207 case CEED_EVAL_NONE:
208 return CeedError(ceed, 1,
209 "CEED_EVAL_NONE does not make sense in this context");
210 }
211 }
212 if (U) {
213 ierr = CeedVectorRestoreArrayRead(U, &u); CeedChk(ierr);
214 }
215 ierr = CeedVectorRestoreArray(V, &v); CeedChk(ierr);
216 return 0;
217}
218
241 case CEED_EVAL_NONE:
242 return CeedError(ceed, 1,
243 "CEED_EVAL_NONE does not make sense in this context");
244 }
245 }
246 if (U) {
247 ierr = CeedVectorRestoreArrayRead(U, &u); CeedChk(ierr);
248 }
249 ierr = CeedVectorRestoreArray(V, &v); CeedChk(ierr);
250 return 0;
251}
252
219static int CeedBasisDestroy_Ref(CeedBasis basis) {
253static int CeedBasisDestroyNonTensor_Ref(CeedBasis basis) {
220 return 0;
221}
222
254 return 0;
255}
256
257static int CeedBasisDestroyTensor_Ref(CeedBasis basis) {
258 int ierr;
259 CeedBasis_Ref *impl;
260 ierr = CeedBasisGetData(basis, (void *)&impl); CeedChk(ierr);
261
262 ierr = CeedFree(&impl->colograd1d); CeedChk(ierr);
263 ierr = CeedFree(&impl); CeedChk(ierr);
264
265 return 0;
266}
267
223int CeedBasisCreateTensorH1_Ref(CeedInt dim, CeedInt P1d,
224 CeedInt Q1d, const CeedScalar *interp1d,
225 const CeedScalar *grad1d,
226 const CeedScalar *qref1d,
227 const CeedScalar *qweight1d,
228 CeedBasis basis) {
229 int ierr;
230 Ceed ceed;
231 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr);
268int CeedBasisCreateTensorH1_Ref(CeedInt dim, CeedInt P1d,
269 CeedInt Q1d, const CeedScalar *interp1d,
270 const CeedScalar *grad1d,
271 const CeedScalar *qref1d,
272 const CeedScalar *qweight1d,
273 CeedBasis basis) {
274 int ierr;
275 Ceed ceed;
276 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr);
277 CeedBasis_Ref *impl;
278 ierr = CeedCalloc(1, &impl); CeedChk(ierr);
279 ierr = CeedMalloc(Q1d*Q1d, &impl->colograd1d); CeedChk(ierr);
280 ierr = CeedBasisGetCollocatedGrad(basis, impl->colograd1d); CeedChk(ierr);
281 ierr = CeedBasisSetData(basis, (void *)&impl); CeedChk(ierr);
232
233 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Apply",
234 CeedBasisApply_Ref); CeedChk(ierr);
235 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Destroy",
282
283 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Apply",
284 CeedBasisApply_Ref); CeedChk(ierr);
285 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Destroy",
236 CeedBasisDestroy_Ref); CeedChk(ierr);
286 CeedBasisDestroyTensor_Ref); CeedChk(ierr);
237 return 0;
238}
239
287 return 0;
288}
289
290
291
240int CeedBasisCreateH1_Ref(CeedElemTopology topo, CeedInt dim,
241 CeedInt ndof, CeedInt nqpts,
242 const CeedScalar *interp,
243 const CeedScalar *grad,
244 const CeedScalar *qref,
245 const CeedScalar *qweight,
246 CeedBasis basis) {
247 int ierr;
248 Ceed ceed;
249 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr);
250
251 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Apply",
252 CeedBasisApply_Ref); CeedChk(ierr);
253 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Destroy",
292int CeedBasisCreateH1_Ref(CeedElemTopology topo, CeedInt dim,
293 CeedInt ndof, CeedInt nqpts,
294 const CeedScalar *interp,
295 const CeedScalar *grad,
296 const CeedScalar *qref,
297 const CeedScalar *qweight,
298 CeedBasis basis) {
299 int ierr;
300 Ceed ceed;
301 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr);
302
303 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Apply",
304 CeedBasisApply_Ref); CeedChk(ierr);
305 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Destroy",
254 CeedBasisDestroy_Ref); CeedChk(ierr);
306 CeedBasisDestroyNonTensor_Ref); CeedChk(ierr);
307
255 return 0;
256}
308 return 0;
309}