diff --git a/tests/t218-elemrestriction.c b/tests/t218-elemrestriction.c index 82ec7fe561..3494292741 100644 --- a/tests/t218-elemrestriction.c +++ b/tests/t218-elemrestriction.c @@ -1,68 +1,100 @@ /// @file -/// Test creation, use, and destruction of a curl-conforming oriented element restriction -/// \test Test creation, use, and destruction of a curl-conforming oriented element restriction +/// Test creation, use, and destruction of a curl-conforming oriented element restriction with and without unsigned application +/// \test Test creation, use, and destruction of a curl-conforming oriented element restriction with and without unsigned application #include #include int main(int argc, char **argv) { Ceed ceed; - CeedVector x, y; - CeedInt num_elem = 6, elem_size = 2; + CeedVector x, y, y_unsigned; + CeedInt num_elem = 6, elem_size = 4; CeedInt ind[elem_size * num_elem]; CeedInt8 curl_orients[3 * elem_size * num_elem]; - CeedScalar x_array[num_elem + 1]; - CeedElemRestriction elem_restriction; + CeedScalar x_array[3 * num_elem + 1]; + CeedElemRestriction elem_restriction, elem_restriction_unsigned; CeedInit(argv[1], &ceed); - CeedVectorCreate(ceed, num_elem + 1, &x); - for (CeedInt i = 0; i < num_elem + 1; i++) x_array[i] = 10 + i; + CeedVectorCreate(ceed, 3 * num_elem + 1, &x); + for (CeedInt i = 0; i < 3 * num_elem + 1; i++) x_array[i] = 10 + i; CeedVectorSetArray(x, CEED_MEM_HOST, CEED_USE_POINTER, x_array); CeedVectorCreate(ceed, num_elem * elem_size, &y); + CeedVectorCreate(ceed, num_elem * elem_size, &y_unsigned); for (CeedInt i = 0; i < num_elem; i++) { - ind[2 * i + 0] = i; - ind[2 * i + 1] = i + 1; - curl_orients[3 * 2 * i] = curl_orients[3 * 2 * (i + 1) - 1] = 0; + ind[4 * i + 0] = 3 * i + 0; + ind[4 * i + 1] = 3 * i + 1; + ind[4 * i + 2] = 3 * i + 2; + ind[4 * i + 3] = 3 * i + 3; if (i % 2 > 0) { - // T = [0 -1] - // [-1 0] - curl_orients[3 * 2 * i + 1] = 0; - curl_orients[3 * 2 * i + 2] = -1; - curl_orients[3 * 2 * i + 3] = -1; - curl_orients[3 * 2 * i + 4] = 0; + // T = [ 1 0 0 0] + // [ 0 1 0 0] + // [ 0 0 0 -1] + // [ 0 0 -1 0] + curl_orients[3 * 4 * i + 0] = 0; + curl_orients[3 * 4 * i + 1] = 1; + curl_orients[3 * 4 * i + 2] = 0; + curl_orients[3 * 3 * i + 3] = 0; + curl_orients[3 * 4 * i + 4] = 1; + curl_orients[3 * 4 * i + 5] = 0; + curl_orients[3 * 4 * i + 6] = 0; + curl_orients[3 * 4 * i + 7] = 0; + curl_orients[3 * 4 * i + 8] = -1; + curl_orients[3 * 4 * i + 9] = -1; + curl_orients[3 * 4 * i + 10] = 0; + curl_orients[3 * 4 * i + 11] = 0; } else { // T = I - curl_orients[3 * 2 * i + 1] = 1; - curl_orients[3 * 2 * i + 2] = 0; - curl_orients[3 * 2 * i + 3] = 0; - curl_orients[3 * 2 * i + 4] = 1; + curl_orients[3 * 4 * i + 0] = 0; + curl_orients[3 * 4 * i + 1] = 1; + curl_orients[3 * 4 * i + 2] = 0; + curl_orients[3 * 4 * i + 3] = 0; + curl_orients[3 * 4 * i + 4] = 1; + curl_orients[3 * 4 * i + 5] = 0; + curl_orients[3 * 4 * i + 6] = 0; + curl_orients[3 * 4 * i + 7] = 1; + curl_orients[3 * 4 * i + 8] = 0; + curl_orients[3 * 4 * i + 9] = 0; + curl_orients[3 * 4 * i + 10] = 1; + curl_orients[3 * 4 * i + 11] = 0; } } - CeedElemRestrictionCreateCurlOriented(ceed, num_elem, elem_size, 1, 1, num_elem + 1, CEED_MEM_HOST, CEED_USE_POINTER, ind, curl_orients, + CeedElemRestrictionCreateCurlOriented(ceed, num_elem, elem_size, 1, 1, 3 * num_elem + 1, CEED_MEM_HOST, CEED_USE_POINTER, ind, curl_orients, &elem_restriction); + CeedElemRestrictionCreateUnsignedCopy(elem_restriction, &elem_restriction_unsigned); // NoTranspose CeedElemRestrictionApply(elem_restriction, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE); + CeedElemRestrictionApply(elem_restriction_unsigned, CEED_NOTRANSPOSE, x, y_unsigned, CEED_REQUEST_IMMEDIATE); { - const CeedScalar *y_array; + const CeedScalar *y_array, *y_unsigned_array; CeedVectorGetArrayRead(y, CEED_MEM_HOST, &y_array); + CeedVectorGetArrayRead(y_unsigned, CEED_MEM_HOST, &y_unsigned_array); for (CeedInt i = 0; i < num_elem; i++) { for (CeedInt j = 0; j < elem_size; j++) { CeedInt k = j + elem_size * i; - if (i % 2 > 0) { - if (j == 0 && 10 + i + 1 != -y_array[k]) { + if (i % 2 > 0 && j >= 2) { + if (j == 2 && 10 + 3 * i + j + 1 != -y_array[k]) { // LCOV_EXCL_START printf("Error in restricted array y[%" CeedInt_FMT "] = %f\n", k, (CeedScalar)y_array[k]); // LCOV_EXCL_STOP - } else if (j == 1 && 10 + i != -y_array[k]) { + } else if (j == 3 && 10 + 3 * i + j - 1 != -y_array[k]) { // LCOV_EXCL_START printf("Error in restricted array y[%" CeedInt_FMT "] = %f\n", k, (CeedScalar)y_array[k]); // LCOV_EXCL_STOP } + if (j == 2 && 10 + 3 * i + j + 1 != y_unsigned_array[k]) { + // LCOV_EXCL_START + printf("Error in unsigned restricted array y[%" CeedInt_FMT "] = %f\n", k, (CeedScalar)y_unsigned_array[k]); + // LCOV_EXCL_STOP + } else if (j == 3 && 10 + 3 * i + j - 1 != y_unsigned_array[k]) { + // LCOV_EXCL_START + printf("Error in unsigned restricted array y[%" CeedInt_FMT "] = %f\n", k, (CeedScalar)y_unsigned_array[k]); + // LCOV_EXCL_STOP + } } else { - if (10 + (k + 1) / 2 != y_array[k]) { + if (10 + 3 * i + j != y_array[k] || 10 + 3 * i + j != y_unsigned_array[k]) { // LCOV_EXCL_START printf("Error in restricted array y[%" CeedInt_FMT "] = %f\n", k, (CeedScalar)y_array[k]); // LCOV_EXCL_STOP @@ -71,6 +103,7 @@ int main(int argc, char **argv) { } } CeedVectorRestoreArrayRead(y, &y_array); + CeedVectorRestoreArrayRead(y_unsigned, &y_unsigned_array); } // Transpose @@ -80,8 +113,22 @@ int main(int argc, char **argv) { const CeedScalar *x_array; CeedVectorGetArrayRead(x, CEED_MEM_HOST, &x_array); - for (CeedInt i = 0; i < num_elem + 1; i++) { - if (x_array[i] != (10 + i) * (i > 0 && i < num_elem ? 2.0 : 1.0)) + for (CeedInt i = 0; i < 3 * num_elem + 1; i++) { + if (x_array[i] != (10 + i) * (i > 0 && i < 3 * num_elem && i % 3 == 0 ? 2.0 : 1.0)) + printf("Error in restricted array x[%" CeedInt_FMT "] = %f\n", i, (CeedScalar)x_array[i]); + } + CeedVectorRestoreArrayRead(x, &x_array); + } + + // Transpose unsigned + CeedVectorSetValue(x, 0); + CeedElemRestrictionApply(elem_restriction_unsigned, CEED_TRANSPOSE, y_unsigned, x, CEED_REQUEST_IMMEDIATE); + { + const CeedScalar *x_array; + + CeedVectorGetArrayRead(x, CEED_MEM_HOST, &x_array); + for (CeedInt i = 0; i < 3 * num_elem + 1; i++) { + if (x_array[i] != (10 + i) * (i > 0 && i < 3 * num_elem && i % 3 == 0 ? 2.0 : 1.0)) printf("Error in restricted array x[%" CeedInt_FMT "] = %f\n", i, (CeedScalar)x_array[i]); } CeedVectorRestoreArrayRead(x, &x_array); @@ -89,7 +136,9 @@ int main(int argc, char **argv) { CeedVectorDestroy(&x); CeedVectorDestroy(&y); + CeedVectorDestroy(&y_unsigned); CeedElemRestrictionDestroy(&elem_restriction); + CeedElemRestrictionDestroy(&elem_restriction_unsigned); CeedDestroy(&ceed); return 0; } diff --git a/tests/t221-elemrestriction.c b/tests/t220-elemrestriction.c similarity index 90% rename from tests/t221-elemrestriction.c rename to tests/t220-elemrestriction.c index 5ab326eb9e..19e1113ccb 100644 --- a/tests/t221-elemrestriction.c +++ b/tests/t220-elemrestriction.c @@ -1,6 +1,6 @@ /// @file -/// Test creation, use, and destruction of an element restriction oriented with unsigned application -/// \test Test creation, use, and destruction of an element restriction oriented with unsigned application +/// Test creation, use, and destruction of an oriented element restriction with unsigned application +/// \test Test creation, use, and destruction of an oriented element restriction with unsigned application #include #include #include @@ -25,10 +25,9 @@ int main(int argc, char **argv) { CeedVectorCreate(ceed, num_elem * 2, &y_unsigned_copy); for (CeedInt i = 0; i < num_elem; i++) { - ind[2 * i + 0] = i; - ind[2 * i + 1] = i + 1; - // flip the dofs on element 1,3,... - orient[2 * i + 0] = (i % (2)) * -1 < 0; + ind[2 * i + 0] = i; + ind[2 * i + 1] = i + 1; + orient[2 * i + 0] = (i % (2)) * -1 < 0; // flip the dofs on element 1, 3, ... orient[2 * i + 1] = (i % (2)) * -1 < 0; } CeedElemRestrictionCreateOriented(ceed, num_elem, p, dim, 1, num_elem + 1, CEED_MEM_HOST, CEED_USE_POINTER, ind, orient, &elem_restriction); @@ -38,7 +37,6 @@ int main(int argc, char **argv) { CeedElemRestrictionApply(elem_restriction, CEED_NOTRANSPOSE, x, y_oriented, CEED_REQUEST_IMMEDIATE); CeedElemRestrictionApply(elem_restriction_unsigned, CEED_NOTRANSPOSE, x, y_unsigned, CEED_REQUEST_IMMEDIATE); CeedElemRestrictionApply(elem_restriction_copy, CEED_NOTRANSPOSE, x, y_unsigned_copy, CEED_REQUEST_IMMEDIATE); - { const CeedScalar *y_oriented_array, *y_unsigned_array, *y_unsigned_copy_array; diff --git a/tests/t571-operator.c b/tests/t571-operator.c new file mode 100644 index 0000000000..dda0c83df2 --- /dev/null +++ b/tests/t571-operator.c @@ -0,0 +1,205 @@ +/// @file +/// Test full assembly of mass matrix operator with oriented and curl-oriented element restrictions (see t560) +/// \test Test full assembly of mass matrix operator with oriented and curl-oriented element restrictions +#include +#include +#include +#include + +#include "t510-operator.h" + +int main(int argc, char **argv) { + Ceed ceed; + CeedElemRestriction elem_restriction_x, elem_restriction_q_data; + CeedElemRestriction elem_restriction_u, oriented_elem_restriction_u, curl_oriented_elem_restriction_u; + CeedBasis basis_x, basis_u; + CeedQFunction qf_setup, qf_mass; + CeedOperator op_setup, op_mass, op_mass_oriented, op_mass_curl_oriented; + CeedVector q_data, x; + CeedInt p = 3, q = 4, dim = 2; + CeedInt n_x = 3, n_y = 2; + CeedInt num_elem = n_x * n_y; + CeedInt num_dofs = (n_x * 2 + 1) * (n_y * 2 + 1), num_qpts = num_elem * q * q; + CeedInt ind_x[num_elem * p * p]; + bool orients_u[num_elem * p * p]; + CeedInt8 curl_orients_u[3 * num_elem * p * p]; + CeedScalar assembled_values[num_dofs * num_dofs]; + CeedScalar assembled_values_oriented[num_dofs * num_dofs]; + CeedScalar assembled_values_curl_oriented[num_dofs * num_dofs]; + + CeedInit(argv[1], &ceed); + + // Vectors + CeedVectorCreate(ceed, dim * num_dofs, &x); + { + CeedScalar x_array[dim * num_dofs]; + + for (CeedInt i = 0; i < n_x * 2 + 1; i++) { + for (CeedInt j = 0; j < n_y * 2 + 1; j++) { + x_array[i + j * (n_x * 2 + 1) + 0 * num_dofs] = (CeedScalar)i / (2 * n_x); + x_array[i + j * (n_x * 2 + 1) + 1 * num_dofs] = (CeedScalar)j / (2 * n_y); + } + } + CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); + } + CeedVectorCreate(ceed, num_qpts, &q_data); + + // Restrictions + for (CeedInt i = 0; i < num_elem; i++) { + CeedInt col, row, offset; + col = i % n_x; + row = i / n_x; + offset = col * (p - 1) + row * (n_x * 2 + 1) * (p - 1); + for (CeedInt j = 0; j < p; j++) { + for (CeedInt k = 0; k < p; k++) { + ind_x[p * (p * i + k) + j] = offset + k * (n_x * 2 + 1) + j; + orients_u[p * (p * i + k) + j] = false; + curl_orients_u[3 * (p * (p * i + k) + j) + 0] = 0; + curl_orients_u[3 * (p * (p * i + k) + j) + 1] = 1; + curl_orients_u[3 * (p * (p * i + k) + j) + 2] = 0; + } + } + } + CeedElemRestrictionCreate(ceed, num_elem, p * p, dim, num_dofs, dim * num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x); + CeedElemRestrictionCreate(ceed, num_elem, p * p, 1, 1, num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_u); + CeedElemRestrictionCreateOriented(ceed, num_elem, p * p, 1, 1, num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, orients_u, + &oriented_elem_restriction_u); + CeedElemRestrictionCreateCurlOriented(ceed, num_elem, p * p, 1, 1, num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, curl_orients_u, + &curl_oriented_elem_restriction_u); + + CeedInt strides_q_data[3] = {1, q * q, q * q}; + CeedElemRestrictionCreateStrided(ceed, num_elem, q * q, 1, num_qpts, strides_q_data, &elem_restriction_q_data); + + // Bases + CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, p, q, CEED_GAUSS, &basis_x); + CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u); + + // QFunctions + CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup); + CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); + CeedQFunctionAddInput(qf_setup, "dx", dim * dim, CEED_EVAL_GRAD); + CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE); + + CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass); + CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE); + CeedQFunctionAddInput(qf_mass, "u", 1, CEED_EVAL_INTERP); + CeedQFunctionAddOutput(qf_mass, "v", 1, CEED_EVAL_INTERP); + + // Operators + CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup); + CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); + CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_setup, "rho", elem_restriction_q_data, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); + + CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass); + CeedOperatorSetField(op_mass, "rho", elem_restriction_q_data, CEED_BASIS_COLLOCATED, q_data); + CeedOperatorSetField(op_mass, "u", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_mass, "v", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + + CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass_oriented); + CeedOperatorSetField(op_mass_oriented, "rho", elem_restriction_q_data, CEED_BASIS_COLLOCATED, q_data); + CeedOperatorSetField(op_mass_oriented, "u", oriented_elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_mass_oriented, "v", oriented_elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + + CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass_curl_oriented); + CeedOperatorSetField(op_mass_curl_oriented, "rho", elem_restriction_q_data, CEED_BASIS_COLLOCATED, q_data); + CeedOperatorSetField(op_mass_curl_oriented, "u", curl_oriented_elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_mass_curl_oriented, "v", curl_oriented_elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + + // Apply Setup Operator + CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE); + + // Fully assemble operator + CeedSize num_entries, num_entries_oriented, num_entries_curl_oriented; + CeedInt *rows, *rows_oriented, *rows_curl_oriented; + CeedInt *cols, *cols_oriented, *cols_curl_oriented; + CeedVector assembled, assembled_oriented, assembled_curl_oriented; + + for (CeedInt k = 0; k < num_dofs * num_dofs; ++k) { + assembled_values[k] = 0.0; + assembled_values_oriented[k] = 0.0; + assembled_values_curl_oriented[k] = 0.0; + } + CeedOperatorLinearAssembleSymbolic(op_mass, &num_entries, &rows, &cols); + CeedOperatorLinearAssembleSymbolic(op_mass_oriented, &num_entries_oriented, &rows_oriented, &cols_oriented); + CeedOperatorLinearAssembleSymbolic(op_mass_curl_oriented, &num_entries_curl_oriented, &rows_curl_oriented, &cols_curl_oriented); + CeedVectorCreate(ceed, num_entries, &assembled); + CeedVectorCreate(ceed, num_entries_oriented, &assembled_oriented); + CeedVectorCreate(ceed, num_entries_curl_oriented, &assembled_curl_oriented); + CeedOperatorLinearAssemble(op_mass, assembled); + CeedOperatorLinearAssemble(op_mass_oriented, assembled_oriented); + CeedOperatorLinearAssemble(op_mass_curl_oriented, assembled_curl_oriented); + { + const CeedScalar *assembled_array; + + CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array); + for (CeedInt k = 0; k < num_entries; ++k) { + assembled_values[rows[k] * num_dofs + cols[k]] += assembled_array[k]; + } + CeedVectorRestoreArrayRead(assembled, &assembled_array); + } + { + const CeedScalar *assembled_array; + + CeedVectorGetArrayRead(assembled_oriented, CEED_MEM_HOST, &assembled_array); + for (CeedInt k = 0; k < num_entries_oriented; ++k) { + assembled_values_oriented[rows_oriented[k] * num_dofs + cols_oriented[k]] += assembled_array[k]; + } + CeedVectorRestoreArrayRead(assembled_oriented, &assembled_array); + } + { + const CeedScalar *assembled_array; + + CeedVectorGetArrayRead(assembled_curl_oriented, CEED_MEM_HOST, &assembled_array); + for (CeedInt k = 0; k < num_entries_curl_oriented; ++k) { + assembled_values_curl_oriented[rows_curl_oriented[k] * num_dofs + cols_curl_oriented[k]] += assembled_array[k]; + } + CeedVectorRestoreArrayRead(assembled_curl_oriented, &assembled_array); + } + + // Check output + for (CeedInt i = 0; i < num_dofs; i++) { + for (CeedInt j = 0; j < num_dofs; j++) { + if (fabs(assembled_values_oriented[j * num_dofs + i] - assembled_values[j * num_dofs + i]) > 100. * CEED_EPSILON) { + // LCOV_EXCL_START + printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] Error in oriented assembly: %f != %f\n", i, j, assembled_values_oriented[j * num_dofs + i], + assembled_values[j * num_dofs + i]); + // LCOV_EXCL_STOP + } + if (fabs(assembled_values_curl_oriented[j * num_dofs + i] - assembled_values[j * num_dofs + i]) > 100. * CEED_EPSILON) { + // LCOV_EXCL_START + printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] Error in curl-oriented assembly: %f != %f\n", i, j, + assembled_values_curl_oriented[j * num_dofs + i], assembled_values[j * num_dofs + i]); + // LCOV_EXCL_STOP + } + } + } + + // Cleanup + free(rows); + free(cols); + free(rows_oriented); + free(cols_oriented); + free(rows_curl_oriented); + free(cols_curl_oriented); + CeedVectorDestroy(&x); + CeedVectorDestroy(&q_data); + CeedVectorDestroy(&assembled); + CeedVectorDestroy(&assembled_oriented); + CeedVectorDestroy(&assembled_curl_oriented); + CeedElemRestrictionDestroy(&elem_restriction_u); + CeedElemRestrictionDestroy(&oriented_elem_restriction_u); + CeedElemRestrictionDestroy(&curl_oriented_elem_restriction_u); + CeedElemRestrictionDestroy(&elem_restriction_x); + CeedElemRestrictionDestroy(&elem_restriction_q_data); + CeedBasisDestroy(&basis_u); + CeedBasisDestroy(&basis_x); + CeedQFunctionDestroy(&qf_setup); + CeedQFunctionDestroy(&qf_mass); + CeedOperatorDestroy(&op_setup); + CeedOperatorDestroy(&op_mass); + CeedOperatorDestroy(&op_mass_oriented); + CeedOperatorDestroy(&op_mass_curl_oriented); + CeedDestroy(&ceed); + return 0; +}