Skip to content

Commit 1fefc7a

Browse files
jedbrownjedbrown
authored
Merge pull request #1105 from CEED/jeremy/composite-multiplicity
Add CeedCompositeOperatorGetMultiplicity
2 parents e3fc958 + a00f0c5 commit 1fefc7a

File tree

6 files changed

+346
-62
lines changed

6 files changed

+346
-62
lines changed

‎doc/sphinx/source/releasenotes.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ On this page we provide a summary of the main API changes, new features and exam
1010

1111
- Added {c:func}`CeedOperatorSetName` for more readable {c:func}`CeedOperatorView` output.
1212
- Added {c:func}`CeedBasisCreateProjection` to facilitate interpolation between nodes for separate `CeedBases`.
13+
- Rename and move {c:func}`CeedCompositeOperatorGetNumSub` and {c:func}`CeedCompositeOperatorGetSubList` to public interface.
1314

1415
### New features
1516

‎include/ceed/backend.h

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -262,8 +262,6 @@ CEED_EXTERN int CeedOperatorGetNumArgs(CeedOperator op, CeedInt *num_args);
262262
CEED_EXTERN int CeedOperatorIsSetupDone(CeedOperator op, bool *is_setup_done);
263263
CEED_EXTERN int CeedOperatorGetQFunction(CeedOperator op, CeedQFunction *qf);
264264
CEED_EXTERN int CeedOperatorIsComposite(CeedOperator op, bool *is_composite);
265-
CEED_EXTERN int CeedOperatorGetNumSub(CeedOperator op, CeedInt *num_suboperators);
266-
CEED_EXTERN int CeedOperatorGetSubList(CeedOperator op, CeedOperator **sub_operators);
267265
CEED_EXTERN int CeedOperatorGetData(CeedOperator op, void *data);
268266
CEED_EXTERN int CeedOperatorSetData(CeedOperator op, void *data);
269267
CEED_EXTERN int CeedOperatorReference(CeedOperator op);

‎include/ceed/ceed.h

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -474,6 +474,8 @@ CEED_EXTERN int CeedOperatorSetField(CeedOperator op, const char *field_name, Ce
474474
CEED_EXTERN int CeedOperatorGetFields(CeedOperator op, CeedInt *num_input_fields, CeedOperatorField **input_fields, CeedInt *num_output_fields,
475475
CeedOperatorField **output_fields);
476476
CEED_EXTERN int CeedCompositeOperatorAddSub(CeedOperator composite_op, CeedOperator sub_op);
477+
CEED_EXTERN int CeedCompositeOperatorGetNumSub(CeedOperator op, CeedInt *num_suboperators);
478+
CEED_EXTERN int CeedCompositeOperatorGetSubList(CeedOperator op, CeedOperator **sub_operators);
477479
CEED_EXTERN int CeedOperatorCheckReady(CeedOperator op);
478480
CEED_EXTERN int CeedOperatorGetActiveVectorLengths(CeedOperator op, CeedSize *input_size, CeedSize *output_size);
479481
CEED_EXTERN int CeedOperatorSetQFunctionAssemblyReuse(CeedOperator op, bool reuse_assembly_data);
@@ -487,6 +489,7 @@ CEED_EXTERN int CeedOperatorLinearAssemblePointBlockDiagonal(CeedOperator op, Ce
487489
CEED_EXTERN int CeedOperatorLinearAssembleAddPointBlockDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request);
488490
CEED_EXTERN int CeedOperatorLinearAssembleSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols);
489491
CEED_EXTERN int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values);
492+
CEED_EXTERN int CeedCompositeOperatorGetMultiplicity(CeedOperator op, CeedInt num_skip_indices, CeedInt *skip_indices, CeedVector mult);
490493
CEED_EXTERN int CeedOperatorMultigridLevelCreate(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse,
491494
CeedBasis basis_coarse, CeedOperator *op_coarse, CeedOperator *op_prolong,
492495
CeedOperator *op_restrict);

‎interface/ceed-operator.c

Lines changed: 52 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -297,8 +297,8 @@ static int CeedOperatorContextSetGeneric(CeedOperator op, CeedContextFieldLabel
297297
CeedInt num_sub;
298298
CeedOperator *sub_operators;
299299

300-
CeedCall(CeedOperatorGetNumSub(op, &num_sub));
301-
CeedCall(CeedOperatorGetSubList(op, &sub_operators));
300+
CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub));
301+
CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
302302
if (num_sub != field_label->num_sub_labels) {
303303
// LCOV_EXCL_START
304304
return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED,
@@ -407,48 +407,6 @@ int CeedOperatorIsComposite(CeedOperator op, bool *is_composite) {
407407
return CEED_ERROR_SUCCESS;
408408
}
409409

410-
/**
411-
@brief Get the number of sub_operators associated with a CeedOperator
412-
413-
@param[in] op CeedOperator
414-
@param[out] num_suboperators Variable to store number of sub_operators
415-
416-
@return An error code: 0 - success, otherwise - failure
417-
418-
@ref Backend
419-
**/
420-
421-
int CeedOperatorGetNumSub(CeedOperator op, CeedInt *num_suboperators) {
422-
if (!op->is_composite) {
423-
// LCOV_EXCL_START
424-
return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator");
425-
// LCOV_EXCL_STOP
426-
}
427-
*num_suboperators = op->num_suboperators;
428-
return CEED_ERROR_SUCCESS;
429-
}
430-
431-
/**
432-
@brief Get the list of sub_operators associated with a CeedOperator
433-
434-
@param[in] op CeedOperator
435-
@param[out] sub_operators Variable to store list of sub_operators
436-
437-
@return An error code: 0 - success, otherwise - failure
438-
439-
@ref Backend
440-
**/
441-
442-
int CeedOperatorGetSubList(CeedOperator op, CeedOperator **sub_operators) {
443-
if (!op->is_composite) {
444-
// LCOV_EXCL_START
445-
return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator");
446-
// LCOV_EXCL_STOP
447-
}
448-
*sub_operators = op->sub_operators;
449-
return CEED_ERROR_SUCCESS;
450-
}
451-
452410
/**
453411
@brief Get the backend data of a CeedOperator
454412
@@ -909,6 +867,48 @@ int CeedCompositeOperatorAddSub(CeedOperator composite_op, CeedOperator sub_op)
909867
return CEED_ERROR_SUCCESS;
910868
}
911869

870+
/**
871+
@brief Get the number of sub_operators associated with a CeedOperator
872+
873+
@param[in] op CeedOperator
874+
@param[out] num_suboperators Variable to store number of sub_operators
875+
876+
@return An error code: 0 - success, otherwise - failure
877+
878+
@ref Backend
879+
**/
880+
881+
int CeedCompositeOperatorGetNumSub(CeedOperator op, CeedInt *num_suboperators) {
882+
if (!op->is_composite) {
883+
// LCOV_EXCL_START
884+
return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator");
885+
// LCOV_EXCL_STOP
886+
}
887+
*num_suboperators = op->num_suboperators;
888+
return CEED_ERROR_SUCCESS;
889+
}
890+
891+
/**
892+
@brief Get the list of sub_operators associated with a CeedOperator
893+
894+
@param op CeedOperator
895+
@param[out] sub_operators Variable to store list of sub_operators
896+
897+
@return An error code: 0 - success, otherwise - failure
898+
899+
@ref Backend
900+
**/
901+
902+
int CeedCompositeOperatorGetSubList(CeedOperator op, CeedOperator **sub_operators) {
903+
if (!op->is_composite) {
904+
// LCOV_EXCL_START
905+
return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator");
906+
// LCOV_EXCL_STOP
907+
}
908+
*sub_operators = op->sub_operators;
909+
return CEED_ERROR_SUCCESS;
910+
}
911+
912912
/**
913913
@brief Check if a CeedOperator is ready to be used.
914914
@@ -1218,9 +1218,9 @@ int CeedOperatorGetFlopsEstimate(CeedOperator op, CeedSize *flops) {
12181218
CeedCall(CeedOperatorIsComposite(op, &is_composite));
12191219
if (is_composite) {
12201220
CeedInt num_suboperators;
1221-
CeedCall(CeedOperatorGetNumSub(op, &num_suboperators));
1221+
CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators));
12221222
CeedOperator *sub_operators;
1223-
CeedCall(CeedOperatorGetSubList(op, &sub_operators));
1223+
CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
12241224

12251225
// FLOPs for each suboperator
12261226
for (CeedInt i = 0; i < num_suboperators; i++) {
@@ -1299,8 +1299,8 @@ int CeedOperatorContextGetFieldLabel(CeedOperator op, const char *field_name, Ce
12991299
CeedContextFieldLabel new_field_label;
13001300

13011301
CeedCall(CeedCalloc(1, &new_field_label));
1302-
CeedCall(CeedOperatorGetNumSub(op, &num_sub));
1303-
CeedCall(CeedOperatorGetSubList(op, &sub_operators));
1302+
CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub));
1303+
CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
13041304
CeedCall(CeedCalloc(num_sub, &new_field_label->sub_labels));
13051305
new_field_label->num_sub_labels = num_sub;
13061306

@@ -1437,9 +1437,9 @@ int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out, CeedReques
14371437
CeedCall(op->ApplyComposite(op, in, out, request));
14381438
} else {
14391439
CeedInt num_suboperators;
1440-
CeedCall(CeedOperatorGetNumSub(op, &num_suboperators));
1440+
CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators));
14411441
CeedOperator *sub_operators;
1442-
CeedCall(CeedOperatorGetSubList(op, &sub_operators));
1442+
CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
14431443

14441444
// Zero all output vectors
14451445
if (out != CEED_VECTOR_NONE) {
@@ -1489,9 +1489,9 @@ int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, CeedReq
14891489
CeedCall(op->ApplyAddComposite(op, in, out, request));
14901490
} else {
14911491
CeedInt num_suboperators;
1492-
CeedCall(CeedOperatorGetNumSub(op, &num_suboperators));
1492+
CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators));
14931493
CeedOperator *sub_operators;
1494-
CeedCall(CeedOperatorGetSubList(op, &sub_operators));
1494+
CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
14951495

14961496
for (CeedInt i = 0; i < num_suboperators; i++) {
14971497
CeedCall(CeedOperatorApplyAdd(sub_operators[i], in, out, request));

‎interface/ceed-preconditioning.c

Lines changed: 79 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -398,8 +398,8 @@ static inline int CeedCompositeOperatorLinearAssembleAddDiagonal(CeedOperator op
398398
CeedVector assembled) {
399399
CeedInt num_sub;
400400
CeedOperator *suboperators;
401-
CeedCall(CeedOperatorGetNumSub(op, &num_sub));
402-
CeedCall(CeedOperatorGetSubList(op, &suboperators));
401+
CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub));
402+
CeedCall(CeedCompositeOperatorGetSubList(op, &suboperators));
403403
for (CeedInt i = 0; i < num_sub; i++) {
404404
if (is_pointblock) {
405405
CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(suboperators[i], assembled, request));
@@ -1670,8 +1670,8 @@ int CeedOperatorLinearAssembleSymbolic(CeedOperator op, CeedSize *num_entries, C
16701670
CeedCall(CeedOperatorIsComposite(op, &is_composite));
16711671
*num_entries = 0;
16721672
if (is_composite) {
1673-
CeedCall(CeedOperatorGetNumSub(op, &num_suboperators));
1674-
CeedCall(CeedOperatorGetSubList(op, &sub_operators));
1673+
CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators));
1674+
CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
16751675
for (CeedInt k = 0; k < num_suboperators; ++k) {
16761676
CeedCall(CeedSingleOperatorAssemblyCountEntries(sub_operators[k], &single_entries));
16771677
*num_entries += single_entries;
@@ -1686,8 +1686,8 @@ int CeedOperatorLinearAssembleSymbolic(CeedOperator op, CeedSize *num_entries, C
16861686
// assemble nonzero locations
16871687
CeedInt offset = 0;
16881688
if (is_composite) {
1689-
CeedCall(CeedOperatorGetNumSub(op, &num_suboperators));
1690-
CeedCall(CeedOperatorGetSubList(op, &sub_operators));
1689+
CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators));
1690+
CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
16911691
for (CeedInt k = 0; k < num_suboperators; ++k) {
16921692
CeedCall(CeedSingleOperatorAssembleSymbolic(sub_operators[k], offset, *rows, *cols));
16931693
CeedCall(CeedSingleOperatorAssemblyCountEntries(sub_operators[k], &single_entries));
@@ -1744,8 +1744,8 @@ int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values) {
17441744

17451745
CeedInt offset = 0;
17461746
if (is_composite) {
1747-
CeedCall(CeedOperatorGetNumSub(op, &num_suboperators));
1748-
CeedCall(CeedOperatorGetSubList(op, &sub_operators));
1747+
CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators));
1748+
CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
17491749
for (CeedInt k = 0; k < num_suboperators; k++) {
17501750
CeedCall(CeedSingleOperatorAssemble(sub_operators[k], offset, values));
17511751
CeedCall(CeedSingleOperatorAssemblyCountEntries(sub_operators[k], &single_entries));
@@ -1758,6 +1758,77 @@ int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values) {
17581758
return CEED_ERROR_SUCCESS;
17591759
}
17601760

1761+
/**
1762+
@brief Get the multiplicity of nodes across suboperators in a composite CeedOperator
1763+
1764+
Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable.
1765+
1766+
@param[in] op Composite CeedOperator
1767+
@param[in] num_skip_indices Number of suboperators to skip
1768+
@param[in] skip_indices Array of indices of suboperators to skip
1769+
@param[out] mult Vector to store multiplicity (of size l_size)
1770+
1771+
@return An error code: 0 - success, otherwise - failure
1772+
1773+
@ref User
1774+
**/
1775+
int CeedCompositeOperatorGetMultiplicity(CeedOperator op, CeedInt num_skip_indices, CeedInt *skip_indices, CeedVector mult) {
1776+
CeedCall(CeedOperatorCheckReady(op));
1777+
1778+
Ceed ceed;
1779+
CeedInt num_sub_ops;
1780+
CeedSize l_vec_len;
1781+
CeedScalar *mult_array;
1782+
CeedVector ones_l_vec;
1783+
CeedElemRestriction elem_restr;
1784+
CeedOperator *sub_ops;
1785+
1786+
CeedCall(CeedOperatorGetCeed(op, &ceed));
1787+
1788+
// Zero mult vector
1789+
CeedCall(CeedVectorSetValue(mult, 0.0));
1790+
1791+
// Get suboperators
1792+
CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub_ops));
1793+
CeedCall(CeedCompositeOperatorGetSubList(op, &sub_ops));
1794+
if (num_sub_ops == 0) return CEED_ERROR_SUCCESS;
1795+
1796+
// Work vector
1797+
CeedCall(CeedVectorGetLength(mult, &l_vec_len));
1798+
CeedCall(CeedVectorCreate(ceed, l_vec_len, &ones_l_vec));
1799+
CeedCall(CeedVectorSetValue(ones_l_vec, 1.0));
1800+
CeedCall(CeedVectorGetArray(mult, CEED_MEM_HOST, &mult_array));
1801+
1802+
// Compute multiplicity across suboperators
1803+
for (CeedInt i = 0; i < num_sub_ops; i++) {
1804+
const CeedScalar *sub_mult_array;
1805+
CeedVector sub_mult_l_vec, ones_e_vec;
1806+
1807+
// -- Check for suboperator to skip
1808+
for (CeedInt j = 0; j < num_skip_indices; j++) {
1809+
if (skip_indices[j] == i) continue;
1810+
}
1811+
1812+
// -- Sub operator multiplicity
1813+
CeedCall(CeedOperatorGetActiveElemRestriction(sub_ops[i], &elem_restr));
1814+
CeedCall(CeedElemRestrictionCreateVector(elem_restr, &sub_mult_l_vec, &ones_e_vec));
1815+
CeedCall(CeedVectorSetValue(sub_mult_l_vec, 0.0));
1816+
CeedCall(CeedElemRestrictionApply(elem_restr, CEED_NOTRANSPOSE, ones_l_vec, ones_e_vec, CEED_REQUEST_IMMEDIATE));
1817+
CeedCall(CeedElemRestrictionApply(elem_restr, CEED_TRANSPOSE, ones_e_vec, sub_mult_l_vec, CEED_REQUEST_IMMEDIATE));
1818+
CeedCall(CeedVectorGetArrayRead(sub_mult_l_vec, CEED_MEM_HOST, &sub_mult_array));
1819+
// ---- Flag every node present in the current suboperator
1820+
for (CeedInt j = 0; j < l_vec_len; j++) {
1821+
if (sub_mult_array[j] > 0.0) mult_array[j] += 1.0;
1822+
}
1823+
CeedCall(CeedVectorRestoreArrayRead(sub_mult_l_vec, &sub_mult_array));
1824+
CeedCall(CeedVectorDestroy(&sub_mult_l_vec));
1825+
CeedCall(CeedVectorDestroy(&ones_e_vec));
1826+
}
1827+
CeedCall(CeedVectorRestoreArray(mult, &mult_array));
1828+
1829+
return CEED_ERROR_SUCCESS;
1830+
}
1831+
17611832
/**
17621833
@brief Create a multigrid coarse operator and level transfer operators for a CeedOperator, creating the prolongation basis from the fine and coarse
17631834
grid interpolation

0 commit comments

Comments
 (0)