Skip to content

Commit a7d38f9

Browse files
jeremyltjeremylt
committed
wip
1 parent 7481867 commit a7d38f9

1 file changed

Lines changed: 256 additions & 0 deletions

File tree

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

Lines changed: 256 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1551,6 +1551,262 @@ static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref(CeedOperator op, Ce
15511551
CeedCallBackend(CeedQFunctionDestroy(&qf));
15521552
return CEED_ERROR_SUCCESS;
15531553
}
1554+
}
1555+
1556+
//------------------------------------------------------------------------------
1557+
// Assemble Operator AtPoints
1558+
//------------------------------------------------------------------------------
1559+
static int CeedSingleOperatorLinearAssembl_Ref(CeedOperator op, CeedInt offset, CeedVector values) {
1560+
CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem, num_comp_active = 1;
1561+
CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}, *assembled;
1562+
Ceed ceed;
1563+
CeedVector point_coords = NULL, in_vec, out_vec;
1564+
CeedElemRestriction rstr_points = NULL;
1565+
CeedQFunctionField *qf_input_fields, *qf_output_fields;
1566+
CeedQFunction qf;
1567+
CeedOperatorField *op_input_fields, *op_output_fields;
1568+
CeedOperator_Ref *impl;
1569+
1570+
CeedCallBackend(CeedOperatorGetData(op, &impl));
1571+
CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem));
1572+
CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
1573+
CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
1574+
CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
1575+
1576+
// Setup
1577+
CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op));
1578+
1579+
// Ceed
1580+
{
1581+
Ceed ceed_parent;
1582+
1583+
CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
1584+
CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
1585+
CeedCallBackend(CeedReferenceCopy(ceed_parent, &ceed));
1586+
CeedCallBackend(CeedDestroy(&ceed_parent));
1587+
}
1588+
1589+
// Point coordinates
1590+
CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords));
1591+
1592+
// Input and output vectors
1593+
{
1594+
CeedSize input_size, output_size;
1595+
1596+
CeedCallBackend(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size));
1597+
CeedCallBackend(CeedVectorCreate(ceed, input_size, &in_vec));
1598+
CeedCallBackend(CeedVectorCreate(ceed, output_size, &out_vec));
1599+
CeedCallBackend(CeedVectorSetValue(out_vec, 0.0));
1600+
}
1601+
1602+
// Assembled array
1603+
CeedCallBackend(CeedVectorGetArray(values, CEED_MEM_HOST, &assembled));
1604+
1605+
// Clear input Qvecs
1606+
for (CeedInt i = 0; i < num_input_fields; i++) {
1607+
bool is_active;
1608+
CeedVector vec;
1609+
1610+
CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
1611+
is_active = vec == CEED_VECTOR_ACTIVE;
1612+
CeedCallBackend(CeedVectorDestroy(&vec));
1613+
if (!is_active) continue;
1614+
CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0));
1615+
}
1616+
1617+
// Input Evecs and Restriction
1618+
CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request));
1619+
1620+
// Loop through elements
1621+
for (CeedInt e = 0; e < num_elem; e++) {
1622+
CeedInt num_points, e_vec_size = 0;
1623+
1624+
// Setup points for element
1625+
CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request));
1626+
CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points));
1627+
1628+
// Input basis apply for non-active bases
1629+
CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec,
1630+
impl->point_coords_elem, true, e_data, impl, request));
1631+
1632+
// ---- Clear output
1633+
for (CeedInt k = 0; k < elem_size * elem_size) {
1634+
assembled[offset + e * elem_size * elem_size + k] = 0.0;
1635+
}
1636+
1637+
// Loop over points on element
1638+
for (CeedInt i = 0; i < num_input_fields; i++) {
1639+
bool is_active_at_points = true, is_active;
1640+
CeedInt elem_size_active = 1;
1641+
CeedRestrictionType rstr_type;
1642+
CeedVector vec;
1643+
CeedElemRestriction elem_rstr;
1644+
1645+
// -- Skip non-active input
1646+
CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
1647+
is_active = vec == CEED_VECTOR_ACTIVE;
1648+
CeedCallBackend(CeedVectorDestroy(&vec));
1649+
if (!is_active) continue;
1650+
1651+
// -- Get active restriction type
1652+
CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
1653+
CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type));
1654+
is_active_at_points = rstr_type == CEED_RESTRICTION_POINTS;
1655+
if (!is_active_at_points) CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size_active));
1656+
else elem_size_active = num_points;
1657+
CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp_active));
1658+
CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1659+
1660+
e_vec_size = elem_size_active * num_comp_active;
1661+
for (CeedInt s = 0; s < e_vec_size; s++) {
1662+
CeedEvalMode eval_mode;
1663+
CeedBasis basis;
1664+
1665+
// -- Update unit vector
1666+
{
1667+
CeedScalar *array;
1668+
1669+
if (s == 0) CeedCallBackend(CeedVectorSetValue(impl->e_vecs_in[i], 0.0));
1670+
CeedCallBackend(CeedVectorGetArray(impl->e_vecs_in[i], CEED_MEM_HOST, &array));
1671+
array[s] = 1.0;
1672+
if (s > 0) array[s - 1] = 0.0;
1673+
CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_in[i], &array));
1674+
}
1675+
// -- Basis action
1676+
CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
1677+
switch (eval_mode) {
1678+
case CEED_EVAL_NONE:
1679+
break;
1680+
// Note - these basis eval modes require FEM fields
1681+
case CEED_EVAL_INTERP:
1682+
case CEED_EVAL_GRAD:
1683+
case CEED_EVAL_DIV:
1684+
case CEED_EVAL_CURL:
1685+
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
1686+
CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_NOTRANSPOSE, eval_mode, impl->point_coords_elem, impl->e_vecs_in[i],
1687+
impl->q_vecs_in[i]));
1688+
CeedCallBackend(CeedBasisDestroy(&basis));
1689+
break;
1690+
case CEED_EVAL_WEIGHT:
1691+
break; // No action
1692+
}
1693+
1694+
// -- Q function
1695+
if (!impl->is_identity_qf) {
1696+
CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out));
1697+
}
1698+
1699+
// -- Output basis apply and restriction
1700+
CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields,
1701+
num_output_fields, impl->apply_add_basis_out, impl->skip_rstr_out, op, out_vec,
1702+
impl->point_coords_elem, impl, request));
1703+
1704+
// -- Build element matrix
1705+
for (CeedInt j = 0; j < num_output_fields; j++) {
1706+
bool is_active;
1707+
CeedInt elem_size = 0;
1708+
CeedRestrictionType rstr_type;
1709+
CeedEvalMode eval_mode;
1710+
CeedVector vec;
1711+
CeedElemRestriction elem_rstr;
1712+
CeedBasis basis;
1713+
1714+
// ---- Skip non-active output
1715+
CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[j], &vec));
1716+
is_active = vec == CEED_VECTOR_ACTIVE;
1717+
CeedCallBackend(CeedVectorDestroy(&vec));
1718+
if (!is_active) continue;
1719+
1720+
// ---- Check if elem size matches
1721+
CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[j], &elem_rstr));
1722+
CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type));
1723+
if (is_active_at_points && rstr_type != CEED_RESTRICTION_POINTS) {
1724+
CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1725+
continue;
1726+
}
1727+
if (rstr_type == CEED_RESTRICTION_POINTS) {
1728+
CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(elem_rstr, e, &elem_size));
1729+
} else {
1730+
CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
1731+
}
1732+
{
1733+
CeedInt num_comp = 0;
1734+
1735+
CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp));
1736+
if (e_vec_size != num_comp * elem_size) {
1737+
CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1738+
continue;
1739+
}
1740+
}
1741+
1742+
// ---- Basis action
1743+
CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[j], &eval_mode));
1744+
switch (eval_mode) {
1745+
case CEED_EVAL_NONE:
1746+
break; // No action
1747+
case CEED_EVAL_INTERP:
1748+
case CEED_EVAL_GRAD:
1749+
case CEED_EVAL_DIV:
1750+
case CEED_EVAL_CURL:
1751+
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[j], &basis));
1752+
CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_TRANSPOSE, eval_mode, impl->point_coords_elem, impl->q_vecs_out[j],
1753+
impl->e_vecs_out[j]));
1754+
CeedCallBackend(CeedBasisDestroy(&basis));
1755+
break;
1756+
// LCOV_EXCL_START
1757+
case CEED_EVAL_WEIGHT: {
1758+
return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode");
1759+
// LCOV_EXCL_STOP
1760+
}
1761+
}
1762+
// ---- Copy output
1763+
{
1764+
const CeedScalar *output;
1765+
1766+
CeedCallBackend(CeedVectorGetArray(impl->e_vecs_out[j], CEED_MEM_HOST, &output));
1767+
for (CeedInt k = 0; k < elem_size) {
1768+
assembled[offset + e * elem_size * elem_size + i * elem_size + k] += output[k];
1769+
}
1770+
CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[j], &output));
1771+
}
1772+
1773+
// COPIED CODE NEED TO UPDATE ----------------------------------------------------------------------------------------------------------
1774+
1775+
/**
1776+
// Put element matrix in coordinate data structure
1777+
for (CeedInt i = 0; i < elem_size_out; i++) {
1778+
for (CeedInt j = 0; j < elem_size_in; j++) {
1779+
assembled[offset + count] = elem_mat[i * elem_size_in + j];
1780+
count++;
1781+
}
1782+
}
1783+
**/
1784+
1785+
// END OF COPIED CODE ------------------------------------------------------------------------------------------------------------------
1786+
1787+
1788+
// -- Reset unit vector
1789+
if (s == e_vec_size - 1) CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0));
1790+
}
1791+
}
1792+
num_points_offset += num_points;
1793+
}
1794+
1795+
// Restore input arrays
1796+
CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl));
1797+
1798+
// Restore assembled values
1799+
CeedCallBackend(CeedVectorRestoreArray(values, &assembled));
1800+
1801+
// Cleanup
1802+
CeedCallBackend(CeedDestroy(&ceed));
1803+
CeedCallBackend(CeedVectorDestroy(&in_vec));
1804+
CeedCallBackend(CeedVectorDestroy(&out_vec));
1805+
CeedCallBackend(CeedVectorDestroy(&point_coords));
1806+
CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
1807+
CeedCallBackend(CeedQFunctionDestroy(&qf));
1808+
return CEED_ERROR_SUCCESS;
1809+
}
15541810

15551811
//------------------------------------------------------------------------------
15561812
// Operator Destroy

0 commit comments

Comments
 (0)