Skip to content

Commit 81dbaa2

Browse files
jeremyltjeremylt
committed
wip
1 parent 7efa755 commit 81dbaa2

5 files changed

Lines changed: 460 additions & 376 deletions

File tree

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

Lines changed: 29 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -1527,7 +1527,7 @@ static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref(CeedOperator op, Ce
15271527
//------------------------------------------------------------------------------
15281528
// Assemble Operator AtPoints
15291529
//------------------------------------------------------------------------------
1530-
static int CeedSingleOperatorAssemble_Ref(CeedOperator op, CeedInt offset, CeedVector values) {
1530+
static int CeedSingleOperatorAssembleAtPoints_Ref(CeedOperator op, CeedInt offset, CeedVector values) {
15311531
CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem, num_comp_active = 1;
15321532
CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}, *assembled;
15331533
Ceed ceed;
@@ -1573,16 +1573,16 @@ static int CeedSingleOperatorAssemble_Ref(CeedOperator op, CeedInt offset, CeedV
15731573
// Assembled array
15741574
CeedCallBackend(CeedVectorGetArray(values, CEED_MEM_HOST, &assembled));
15751575

1576-
// Clear input Qvecs
1576+
// Clear input Evecs
15771577
for (CeedInt i = 0; i < num_input_fields; i++) {
15781578
bool is_active;
15791579
CeedVector vec;
15801580

15811581
CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
15821582
is_active = vec == CEED_VECTOR_ACTIVE;
15831583
CeedCallBackend(CeedVectorDestroy(&vec));
1584-
if (!is_active) continue;
1585-
CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0));
1584+
if (!is_active || impl->skip_rstr_in[i]) continue;
1585+
CeedCallBackend(CeedVectorSetValue(impl->e_vecs_in[i], 0.0));
15861586
}
15871587

15881588
// Input Evecs and Restriction
@@ -1599,7 +1599,7 @@ static int CeedSingleOperatorAssemble_Ref(CeedOperator op, CeedInt offset, CeedV
15991599

16001600
// Input basis apply for non-active bases
16011601
CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec,
1602-
impl->point_coords_elem, true, e_data, impl, CEED_REQUEST_IMMEDIATE));
1602+
impl->point_coords_elem, true, false, e_data, impl, CEED_REQUEST_IMMEDIATE));
16031603

16041604
// Loop over points on element
16051605
for (CeedInt i = 0; i < num_input_fields; i++) {
@@ -1613,8 +1613,7 @@ static int CeedSingleOperatorAssemble_Ref(CeedOperator op, CeedInt offset, CeedV
16131613
CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
16141614
is_active = vec == CEED_VECTOR_ACTIVE;
16151615
CeedCallBackend(CeedVectorDestroy(&vec));
1616-
if (!is_active) continue;
1617-
if (impl->skip_rstr_in[i]) continue;
1616+
if (!is_active || impl->skip_rstr_in[i]) continue;
16181617

16191618
// -- Get active restriction type
16201619
CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
@@ -1627,8 +1626,8 @@ static int CeedSingleOperatorAssemble_Ref(CeedOperator op, CeedInt offset, CeedV
16271626

16281627
e_vec_size = elem_size_active * num_comp_active;
16291628
for (CeedInt s = 0; s < e_vec_size; s++) {
1630-
CeedEvalMode eval_mode;
1631-
CeedBasis basis;
1629+
const CeedInt comp_in = s / elem_size_active;
1630+
const CeedInt node_in = s % elem_size_active;
16321631

16331632
// -- Update unit vector
16341633
{
@@ -1640,50 +1639,33 @@ static int CeedSingleOperatorAssemble_Ref(CeedOperator op, CeedInt offset, CeedV
16401639
if (s > 0) array[s - 1] = 0.0;
16411640
CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_in[i], &array));
16421641
}
1643-
// -- Basis action
1644-
CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
1645-
switch (eval_mode) {
1646-
case CEED_EVAL_NONE:
1647-
break;
1648-
// Note - these basis eval modes require FEM fields
1649-
case CEED_EVAL_INTERP:
1650-
case CEED_EVAL_GRAD:
1651-
case CEED_EVAL_DIV:
1652-
case CEED_EVAL_CURL:
1653-
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
1654-
CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_NOTRANSPOSE, eval_mode, impl->point_coords_elem, impl->e_vecs_in[i],
1655-
impl->q_vecs_in[i]));
1656-
CeedCallBackend(CeedBasisDestroy(&basis));
1657-
break;
1658-
case CEED_EVAL_WEIGHT:
1659-
break; // No action
1660-
}
1642+
// Input basis apply for active bases
1643+
CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields,
1644+
in_vec, impl->point_coords_elem, false, true, e_data, impl, CEED_REQUEST_IMMEDIATE));
16611645

16621646
// -- Q function
16631647
if (!impl->is_identity_qf) {
16641648
CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out));
16651649
}
16661650

1667-
// -- Output basis apply
1651+
// -- Output basis apply and restriction
16681652
CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields,
16691653
num_output_fields, impl->apply_add_basis_out, impl->skip_rstr_out, op, out_vec,
1670-
impl->point_coords_elem, impl, CEED_REQUEST_IMMEDIATE));
1654+
impl->point_coords_elem, true, impl, CEED_REQUEST_IMMEDIATE));
16711655

16721656
// -- Build element matrix
16731657
for (CeedInt j = 0; j < num_output_fields; j++) {
16741658
bool is_active;
16751659
CeedInt elem_size = 0;
16761660
CeedRestrictionType rstr_type;
1677-
CeedEvalMode eval_mode;
16781661
CeedVector vec;
16791662
CeedElemRestriction elem_rstr;
1680-
CeedBasis basis;
16811663

16821664
// ---- Skip non-active output
16831665
CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[j], &vec));
16841666
is_active = vec == CEED_VECTOR_ACTIVE;
16851667
CeedCallBackend(CeedVectorDestroy(&vec));
1686-
if (!is_active) continue;
1668+
if (!is_active || impl->skip_rstr_out[j]) continue;
16871669

16881670
// ---- Check if elem size matches
16891671
CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[j], &elem_rstr));
@@ -1704,39 +1686,28 @@ static int CeedSingleOperatorAssemble_Ref(CeedOperator op, CeedInt offset, CeedV
17041686
CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
17051687
if (e_vec_size != num_comp * elem_size) continue;
17061688
}
1707-
1708-
// ---- Basis action
1709-
CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[j], &eval_mode));
1710-
switch (eval_mode) {
1711-
case CEED_EVAL_NONE:
1712-
break; // No action
1713-
case CEED_EVAL_INTERP:
1714-
case CEED_EVAL_GRAD:
1715-
case CEED_EVAL_DIV:
1716-
case CEED_EVAL_CURL:
1717-
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[j], &basis));
1718-
CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_TRANSPOSE, eval_mode, impl->point_coords_elem, impl->q_vecs_out[j],
1719-
impl->e_vecs_out[j]));
1720-
CeedCallBackend(CeedBasisDestroy(&basis));
1721-
break;
1722-
// LCOV_EXCL_START
1723-
case CEED_EVAL_WEIGHT: {
1724-
return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode");
1725-
// LCOV_EXCL_STOP
1726-
}
1727-
}
17281689
// ---- Copy output
17291690
{
17301691
const CeedScalar *output;
17311692

17321693
CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_out[j], CEED_MEM_HOST, &output));
1733-
for (CeedInt k = 0; k < elem_size_active; k++) {
1734-
assembled[offset + e * e_vec_size * e_vec_size + s * elem_size_active + k] += output[k];
1694+
for (CeedInt k = 0; k < e_vec_size; k++) {
1695+
const CeedInt comp_out = k / elem_size_active;
1696+
const CeedInt node_out = k % elem_size_active;
1697+
1698+
assembled[offset + e * e_vec_size * e_vec_size + (comp_in * num_comp_active + comp_out) * elem_size_active * elem_size_active +
1699+
node_out * elem_size_active + node_in] = output[k];
17351700
}
17361701
CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_out[j], &output));
17371702
}
1738-
// -- Reset unit vector
1739-
if (s == e_vec_size - 1) CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0));
1703+
}
1704+
// -- Reset unit vector
1705+
if (s == e_vec_size - 1) {
1706+
CeedScalar *array;
1707+
1708+
CeedCallBackend(CeedVectorGetArray(impl->e_vecs_in[i], CEED_MEM_HOST, &array));
1709+
array[s] = 0.0;
1710+
CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_in[i], &array));
17401711
}
17411712
}
17421713
}
@@ -1827,7 +1798,7 @@ int CeedOperatorCreateAtPoints_Ref(CeedOperator op) {
18271798
CeedCallBackend(
18281799
CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref));
18291800
CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref));
1830-
CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleSingle", CeedSingleOperatorAssemble_Ref));
1801+
CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleSingle", CeedSingleOperatorAssembleAtPoints_Ref));
18311802
CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Ref));
18321803
CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref));
18331804
CeedCallBackend(CeedDestroy(&ceed));

0 commit comments

Comments
 (0)