Skip to content

Commit 2efebff

Browse files
sebastiangrimbergsebastiangrimberg
authored
Merge pull request #1156 from sebastiangrimberg/sjg/hcurl-basis-dev
H(curl) basis and `CEED_EVAL_CURL`
2 parents f724464 + b2e3f8e commit 2efebff

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

48 files changed

+1897
-613
lines changed

‎backends/blocked/ceed-blocked-operator.c

+5-39
Original file line numberDiff line numberDiff line change
@@ -77,6 +77,7 @@ static int CeedOperatorSetupFields_Blocked(CeedQFunction qf, CeedOperator op, bo
7777
case CEED_EVAL_INTERP:
7878
case CEED_EVAL_GRAD:
7979
case CEED_EVAL_DIV:
80+
case CEED_EVAL_CURL:
8081
CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis));
8182
CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size));
8283
CeedCallBackend(CeedBasisGetNumNodes(basis, &P));
@@ -92,8 +93,6 @@ static int CeedOperatorSetupFields_Blocked(CeedQFunction qf, CeedOperator op, bo
9293
CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
9394
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i]));
9495
break;
95-
case CEED_EVAL_CURL:
96-
break; // Not implemented
9796
}
9897
}
9998
return CEED_ERROR_SUCCESS;
@@ -227,33 +226,16 @@ static inline int CeedOperatorInputBasis_Blocked(CeedInt e, CeedInt Q, CeedQFunc
227226
CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * Q * size]));
228227
break;
229228
case CEED_EVAL_INTERP:
230-
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
231-
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
232-
CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp]));
233-
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, impl->e_vecs_in[i], impl->q_vecs_in[i]));
234-
break;
235229
case CEED_EVAL_GRAD:
236-
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
237-
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
238-
CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp]));
239-
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, impl->e_vecs_in[i], impl->q_vecs_in[i]));
240-
break;
241230
case CEED_EVAL_DIV:
231+
case CEED_EVAL_CURL:
242232
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
243233
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
244234
CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp]));
245-
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_DIV, impl->e_vecs_in[i], impl->q_vecs_in[i]));
235+
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs_in[i], impl->q_vecs_in[i]));
246236
break;
247237
case CEED_EVAL_WEIGHT:
248238
break; // No action
249-
// LCOV_EXCL_START
250-
case CEED_EVAL_CURL: {
251-
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
252-
Ceed ceed;
253-
CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
254-
return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented");
255-
// LCOV_EXCL_STOP
256-
}
257239
}
258240
}
259241
return CEED_ERROR_SUCCESS;
@@ -280,36 +262,20 @@ static inline int CeedOperatorOutputBasis_Blocked(CeedInt e, CeedInt Q, CeedQFun
280262
case CEED_EVAL_NONE:
281263
break; // No action
282264
case CEED_EVAL_INTERP:
283-
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
284-
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
285-
CeedCallBackend(
286-
CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp]));
287-
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_INTERP, impl->q_vecs_out[i], impl->e_vecs_out[i]));
288-
break;
289265
case CEED_EVAL_GRAD:
290-
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
291-
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
292-
CeedCallBackend(
293-
CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp]));
294-
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_GRAD, impl->q_vecs_out[i], impl->e_vecs_out[i]));
295-
break;
296266
case CEED_EVAL_DIV:
267+
case CEED_EVAL_CURL:
297268
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
298269
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
299270
CeedCallBackend(
300271
CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp]));
301-
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_DIV, impl->q_vecs_out[i], impl->e_vecs_out[i]));
272+
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i]));
302273
break;
303274
// LCOV_EXCL_START
304275
case CEED_EVAL_WEIGHT: {
305276
Ceed ceed;
306277
CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
307278
return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode");
308-
}
309-
case CEED_EVAL_CURL: {
310-
Ceed ceed;
311-
CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
312-
return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented");
313279
// LCOV_EXCL_STOP
314280
}
315281
}

‎backends/opt/ceed-opt-operator.c

+5-37
Original file line numberDiff line numberDiff line change
@@ -80,6 +80,7 @@ static int CeedOperatorSetupFields_Opt(CeedQFunction qf, CeedOperator op, bool i
8080
case CEED_EVAL_INTERP:
8181
case CEED_EVAL_GRAD:
8282
case CEED_EVAL_DIV:
83+
case CEED_EVAL_CURL:
8384
CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis));
8485
CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size));
8586
CeedCallBackend(CeedBasisGetNumNodes(basis, &P));
@@ -95,8 +96,6 @@ static int CeedOperatorSetupFields_Opt(CeedQFunction qf, CeedOperator op, bool i
9596
CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
9697
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i]));
9798
break;
98-
case CEED_EVAL_CURL:
99-
break; // Not implemented
10099
}
101100
if (is_input && !!e_vecs[i]) {
102101
CeedCallBackend(CeedVectorSetArray(e_vecs[i], CEED_MEM_HOST, CEED_COPY_VALUES, NULL));
@@ -248,39 +247,18 @@ static inline int CeedOperatorInputBasis_Opt(CeedInt e, CeedInt Q, CeedQFunction
248247
}
249248
break;
250249
case CEED_EVAL_INTERP:
251-
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
252-
if (!active_in) {
253-
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
254-
CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][e * elem_size * num_comp]));
255-
}
256-
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, impl->e_vecs_in[i], impl->q_vecs_in[i]));
257-
break;
258250
case CEED_EVAL_GRAD:
259-
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
260-
if (!active_in) {
261-
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
262-
CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][e * elem_size * num_comp]));
263-
}
264-
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, impl->e_vecs_in[i], impl->q_vecs_in[i]));
265-
break;
266251
case CEED_EVAL_DIV:
252+
case CEED_EVAL_CURL:
267253
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
268254
if (!active_in) {
269255
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
270256
CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][e * elem_size * num_comp]));
271257
}
272-
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_DIV, impl->e_vecs_in[i], impl->q_vecs_in[i]));
258+
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs_in[i], impl->q_vecs_in[i]));
273259
break;
274260
case CEED_EVAL_WEIGHT:
275261
break; // No action
276-
// LCOV_EXCL_START
277-
case CEED_EVAL_CURL: {
278-
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
279-
Ceed ceed;
280-
CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
281-
return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented");
282-
// LCOV_EXCL_STOP
283-
}
284262
}
285263
}
286264
return CEED_ERROR_SUCCESS;
@@ -306,16 +284,11 @@ static inline int CeedOperatorOutputBasis_Opt(CeedInt e, CeedInt Q, CeedQFunctio
306284
case CEED_EVAL_NONE:
307285
break; // No action
308286
case CEED_EVAL_INTERP:
309-
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
310-
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_INTERP, impl->q_vecs_out[i], impl->e_vecs_out[i]));
311-
break;
312287
case CEED_EVAL_GRAD:
313-
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
314-
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_GRAD, impl->q_vecs_out[i], impl->e_vecs_out[i]));
315-
break;
316288
case CEED_EVAL_DIV:
289+
case CEED_EVAL_CURL:
317290
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
318-
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_DIV, impl->q_vecs_out[i], impl->e_vecs_out[i]));
291+
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i]));
319292
break;
320293
// LCOV_EXCL_START
321294
case CEED_EVAL_WEIGHT: {
@@ -324,11 +297,6 @@ static inline int CeedOperatorOutputBasis_Opt(CeedInt e, CeedInt Q, CeedQFunctio
324297
return CeedError(ceed, CEED_ERROR_BACKEND,
325298
"CEED_EVAL_WEIGHT cannot be an output "
326299
"evaluation mode");
327-
}
328-
case CEED_EVAL_CURL: {
329-
Ceed ceed;
330-
CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
331-
return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented");
332300
// LCOV_EXCL_STOP
333301
}
334302
}

‎backends/ref/ceed-ref-basis.c

+37-37
Original file line numberDiff line numberDiff line change
@@ -19,12 +19,12 @@
1919
static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U, CeedVector V) {
2020
Ceed ceed;
2121
CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
22-
CeedInt dim, num_comp, num_nodes, num_qpts, Q_comp;
22+
CeedInt dim, num_comp, q_comp, num_nodes, num_qpts;
2323
CeedCallBackend(CeedBasisGetDimension(basis, &dim));
2424
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
25+
CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp));
2526
CeedCallBackend(CeedBasisGetNumNodes(basis, &num_nodes));
2627
CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
27-
CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, &Q_comp));
2828
CeedTensorContract contract;
2929
CeedCallBackend(CeedBasisGetTensorContract(basis, &contract));
3030
const CeedInt add = (t_mode == CEED_TRANSPOSE);
@@ -46,8 +46,8 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo
4646
}
4747
bool tensor_basis;
4848
CeedCallBackend(CeedBasisIsTensor(basis, &tensor_basis));
49-
// Tensor basis
5049
if (tensor_basis) {
50+
// Tensor basis
5151
CeedInt P_1d, Q_1d;
5252
CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d));
5353
CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
@@ -191,36 +191,31 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo
191191
}
192192
} else {
193193
// Non-tensor basis
194+
CeedInt P = num_nodes, Q = num_qpts;
194195
switch (eval_mode) {
195196
// Interpolate to/from quadrature points
196197
case CEED_EVAL_INTERP: {
197-
CeedInt P = num_nodes, Q = Q_comp * num_qpts;
198198
const CeedScalar *interp;
199199
CeedCallBackend(CeedBasisGetInterp(basis, &interp));
200-
if (t_mode == CEED_TRANSPOSE) {
201-
P = Q_comp * num_qpts;
202-
Q = num_nodes;
203-
}
204-
CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, interp, t_mode, add, u, v));
200+
CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, interp, t_mode, add, u, v));
205201
} break;
206202
// Evaluate the gradient to/from quadrature points
207203
case CEED_EVAL_GRAD: {
208-
CeedInt P = num_nodes, Q = num_qpts;
209-
CeedInt dim_stride = num_qpts * num_comp * num_elem;
210-
CeedInt grad_stride = num_qpts * num_nodes;
211204
const CeedScalar *grad;
212205
CeedCallBackend(CeedBasisGetGrad(basis, &grad));
213-
if (t_mode == CEED_TRANSPOSE) {
214-
P = num_qpts;
215-
Q = num_nodes;
216-
for (CeedInt d = 0; d < dim; d++) {
217-
CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, grad + d * grad_stride, t_mode, add, u + d * dim_stride, v));
218-
}
219-
} else {
220-
for (CeedInt d = 0; d < dim; d++) {
221-
CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, grad + d * grad_stride, t_mode, add, u, v + d * dim_stride));
222-
}
223-
}
206+
CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, grad, t_mode, add, u, v));
207+
} break;
208+
// Evaluate the divergence to/from the quadrature points
209+
case CEED_EVAL_DIV: {
210+
const CeedScalar *div;
211+
CeedCallBackend(CeedBasisGetDiv(basis, &div));
212+
CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, div, t_mode, add, u, v));
213+
} break;
214+
// Evaluate the curl to/from the quadrature points
215+
case CEED_EVAL_CURL: {
216+
const CeedScalar *curl;
217+
CeedCallBackend(CeedBasisGetCurl(basis, &curl));
218+
CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, curl, t_mode, add, u, v));
224219
} break;
225220
// Retrieve interpolation weights
226221
case CEED_EVAL_WEIGHT: {
@@ -235,21 +230,7 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo
235230
for (CeedInt e = 0; e < num_elem; e++) v[i * num_elem + e] = q_weight[i];
236231
}
237232
} break;
238-
// Evaluate the divergence to/from the quadrature points
239-
case CEED_EVAL_DIV: {
240-
CeedInt P = num_nodes, Q = num_qpts;
241-
const CeedScalar *div;
242-
CeedCallBackend(CeedBasisGetDiv(basis, &div));
243-
if (t_mode == CEED_TRANSPOSE) {
244-
P = num_qpts;
245-
Q = num_nodes;
246-
}
247-
CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, div, t_mode, add, u, v));
248-
} break;
249233
// LCOV_EXCL_START
250-
// Evaluate the curl to/from the quadrature points
251-
case CEED_EVAL_CURL:
252-
return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_CURL not supported");
253234
// Take no action, BasisApply should not have been called
254235
case CEED_EVAL_NONE:
255236
return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context");
@@ -301,6 +282,25 @@ int CeedBasisCreateHdiv_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_node
301282
return CEED_ERROR_SUCCESS;
302283
}
303284

285+
//------------------------------------------------------------------------------
286+
// Basis Create Non-Tensor H(curl)
287+
//------------------------------------------------------------------------------
288+
int CeedBasisCreateHcurl_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp,
289+
const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
290+
Ceed ceed;
291+
CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
292+
293+
Ceed parent;
294+
CeedCallBackend(CeedGetParent(ceed, &parent));
295+
CeedTensorContract contract;
296+
CeedCallBackend(CeedTensorContractCreate(parent, basis, &contract));
297+
CeedCallBackend(CeedBasisSetTensorContract(basis, contract));
298+
299+
CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref));
300+
301+
return CEED_ERROR_SUCCESS;
302+
}
303+
304304
//------------------------------------------------------------------------------
305305
// Basis Destroy Tensor
306306
//------------------------------------------------------------------------------

0 commit comments

Comments
 (0)