diff --git a/backends/blocked/ceed-blocked-operator.c b/backends/blocked/ceed-blocked-operator.c index 0d0f30cd68..f0d01b5508 100644 --- a/backends/blocked/ceed-blocked-operator.c +++ b/backends/blocked/ceed-blocked-operator.c @@ -77,6 +77,7 @@ static int CeedOperatorSetupFields_Blocked(CeedQFunction qf, CeedOperator op, bo case CEED_EVAL_INTERP: case CEED_EVAL_GRAD: case CEED_EVAL_DIV: + case CEED_EVAL_CURL: CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); @@ -92,8 +93,6 @@ static int CeedOperatorSetupFields_Blocked(CeedQFunction qf, CeedOperator op, bo CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i])); break; - case CEED_EVAL_CURL: - break; // Not implemented } } return CEED_ERROR_SUCCESS; @@ -227,33 +226,16 @@ static inline int CeedOperatorInputBasis_Blocked(CeedInt e, CeedInt Q, CeedQFunc CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * Q * size])); break; case CEED_EVAL_INTERP: - CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); - CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); - CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp])); - CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, impl->e_vecs_in[i], impl->q_vecs_in[i])); - break; case CEED_EVAL_GRAD: - CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); - CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); - CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp])); - CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, impl->e_vecs_in[i], impl->q_vecs_in[i])); - break; case CEED_EVAL_DIV: + case CEED_EVAL_CURL: CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp])); - CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_DIV, impl->e_vecs_in[i], impl->q_vecs_in[i])); + CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs_in[i], impl->q_vecs_in[i])); break; case CEED_EVAL_WEIGHT: break; // No action - // LCOV_EXCL_START - case CEED_EVAL_CURL: { - CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); - Ceed ceed; - CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); - return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented"); - // LCOV_EXCL_STOP - } } } return CEED_ERROR_SUCCESS; @@ -280,36 +262,20 @@ static inline int CeedOperatorOutputBasis_Blocked(CeedInt e, CeedInt Q, CeedQFun case CEED_EVAL_NONE: break; // No action case CEED_EVAL_INTERP: - CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); - CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); - CeedCallBackend( - CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp])); - CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_INTERP, impl->q_vecs_out[i], impl->e_vecs_out[i])); - break; case CEED_EVAL_GRAD: - CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); - CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); - CeedCallBackend( - CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp])); - CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_GRAD, impl->q_vecs_out[i], impl->e_vecs_out[i])); - break; case CEED_EVAL_DIV: + case CEED_EVAL_CURL: CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); CeedCallBackend( CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp])); - CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_DIV, impl->q_vecs_out[i], impl->e_vecs_out[i])); + CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i])); break; // LCOV_EXCL_START case CEED_EVAL_WEIGHT: { Ceed ceed; CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); - } - case CEED_EVAL_CURL: { - Ceed ceed; - CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); - return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented"); // LCOV_EXCL_STOP } } diff --git a/backends/opt/ceed-opt-operator.c b/backends/opt/ceed-opt-operator.c index 23051b4295..b49f4f2b39 100644 --- a/backends/opt/ceed-opt-operator.c +++ b/backends/opt/ceed-opt-operator.c @@ -80,6 +80,7 @@ static int CeedOperatorSetupFields_Opt(CeedQFunction qf, CeedOperator op, bool i case CEED_EVAL_INTERP: case CEED_EVAL_GRAD: case CEED_EVAL_DIV: + case CEED_EVAL_CURL: CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); @@ -95,8 +96,6 @@ static int CeedOperatorSetupFields_Opt(CeedQFunction qf, CeedOperator op, bool i CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i])); break; - case CEED_EVAL_CURL: - break; // Not implemented } if (is_input && !!e_vecs[i]) { 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 } break; case CEED_EVAL_INTERP: - CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); - if (!active_in) { - CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); - CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][e * elem_size * num_comp])); - } - CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, impl->e_vecs_in[i], impl->q_vecs_in[i])); - break; case CEED_EVAL_GRAD: - CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); - if (!active_in) { - CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); - CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][e * elem_size * num_comp])); - } - CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, impl->e_vecs_in[i], impl->q_vecs_in[i])); - break; case CEED_EVAL_DIV: + case CEED_EVAL_CURL: CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); if (!active_in) { CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][e * elem_size * num_comp])); } - CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_DIV, impl->e_vecs_in[i], impl->q_vecs_in[i])); + CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs_in[i], impl->q_vecs_in[i])); break; case CEED_EVAL_WEIGHT: break; // No action - // LCOV_EXCL_START - case CEED_EVAL_CURL: { - CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); - Ceed ceed; - CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); - return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented"); - // LCOV_EXCL_STOP - } } } return CEED_ERROR_SUCCESS; @@ -306,16 +284,11 @@ static inline int CeedOperatorOutputBasis_Opt(CeedInt e, CeedInt Q, CeedQFunctio case CEED_EVAL_NONE: break; // No action case CEED_EVAL_INTERP: - CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); - CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_INTERP, impl->q_vecs_out[i], impl->e_vecs_out[i])); - break; case CEED_EVAL_GRAD: - CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); - CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_GRAD, impl->q_vecs_out[i], impl->e_vecs_out[i])); - break; case CEED_EVAL_DIV: + case CEED_EVAL_CURL: CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); - CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_DIV, impl->q_vecs_out[i], impl->e_vecs_out[i])); + CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i])); break; // LCOV_EXCL_START case CEED_EVAL_WEIGHT: { @@ -324,11 +297,6 @@ static inline int CeedOperatorOutputBasis_Opt(CeedInt e, CeedInt Q, CeedQFunctio return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output " "evaluation mode"); - } - case CEED_EVAL_CURL: { - Ceed ceed; - CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); - return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented"); // LCOV_EXCL_STOP } } diff --git a/backends/ref/ceed-ref-basis.c b/backends/ref/ceed-ref-basis.c index 18e9dd007d..775dd82c67 100644 --- a/backends/ref/ceed-ref-basis.c +++ b/backends/ref/ceed-ref-basis.c @@ -19,12 +19,12 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U, CeedVector V) { Ceed ceed; CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); - CeedInt dim, num_comp, num_nodes, num_qpts, Q_comp; + CeedInt dim, num_comp, q_comp, num_nodes, num_qpts; CeedCallBackend(CeedBasisGetDimension(basis, &dim)); CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); + CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); CeedCallBackend(CeedBasisGetNumNodes(basis, &num_nodes)); CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); - CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, &Q_comp)); CeedTensorContract contract; CeedCallBackend(CeedBasisGetTensorContract(basis, &contract)); const CeedInt add = (t_mode == CEED_TRANSPOSE); @@ -46,8 +46,8 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo } bool tensor_basis; CeedCallBackend(CeedBasisIsTensor(basis, &tensor_basis)); - // Tensor basis if (tensor_basis) { + // Tensor basis CeedInt P_1d, Q_1d; CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); @@ -191,36 +191,31 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo } } else { // Non-tensor basis + CeedInt P = num_nodes, Q = num_qpts; switch (eval_mode) { // Interpolate to/from quadrature points case CEED_EVAL_INTERP: { - CeedInt P = num_nodes, Q = Q_comp * num_qpts; const CeedScalar *interp; CeedCallBackend(CeedBasisGetInterp(basis, &interp)); - if (t_mode == CEED_TRANSPOSE) { - P = Q_comp * num_qpts; - Q = num_nodes; - } - CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, interp, t_mode, add, u, v)); + CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, interp, t_mode, add, u, v)); } break; // Evaluate the gradient to/from quadrature points case CEED_EVAL_GRAD: { - CeedInt P = num_nodes, Q = num_qpts; - CeedInt dim_stride = num_qpts * num_comp * num_elem; - CeedInt grad_stride = num_qpts * num_nodes; const CeedScalar *grad; CeedCallBackend(CeedBasisGetGrad(basis, &grad)); - if (t_mode == CEED_TRANSPOSE) { - P = num_qpts; - Q = num_nodes; - for (CeedInt d = 0; d < dim; d++) { - CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, grad + d * grad_stride, t_mode, add, u + d * dim_stride, v)); - } - } else { - for (CeedInt d = 0; d < dim; d++) { - CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, grad + d * grad_stride, t_mode, add, u, v + d * dim_stride)); - } - } + CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, grad, t_mode, add, u, v)); + } break; + // Evaluate the divergence to/from the quadrature points + case CEED_EVAL_DIV: { + const CeedScalar *div; + CeedCallBackend(CeedBasisGetDiv(basis, &div)); + CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, div, t_mode, add, u, v)); + } break; + // Evaluate the curl to/from the quadrature points + case CEED_EVAL_CURL: { + const CeedScalar *curl; + CeedCallBackend(CeedBasisGetCurl(basis, &curl)); + CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, curl, t_mode, add, u, v)); } break; // Retrieve interpolation weights case CEED_EVAL_WEIGHT: { @@ -235,21 +230,7 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo for (CeedInt e = 0; e < num_elem; e++) v[i * num_elem + e] = q_weight[i]; } } break; - // Evaluate the divergence to/from the quadrature points - case CEED_EVAL_DIV: { - CeedInt P = num_nodes, Q = num_qpts; - const CeedScalar *div; - CeedCallBackend(CeedBasisGetDiv(basis, &div)); - if (t_mode == CEED_TRANSPOSE) { - P = num_qpts; - Q = num_nodes; - } - CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, div, t_mode, add, u, v)); - } break; // LCOV_EXCL_START - // Evaluate the curl to/from the quadrature points - case CEED_EVAL_CURL: - return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_CURL not supported"); // Take no action, BasisApply should not have been called case CEED_EVAL_NONE: 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 return CEED_ERROR_SUCCESS; } +//------------------------------------------------------------------------------ +// Basis Create Non-Tensor H(curl) +//------------------------------------------------------------------------------ +int CeedBasisCreateHcurl_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, + const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) { + Ceed ceed; + CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); + + Ceed parent; + CeedCallBackend(CeedGetParent(ceed, &parent)); + CeedTensorContract contract; + CeedCallBackend(CeedTensorContractCreate(parent, basis, &contract)); + CeedCallBackend(CeedBasisSetTensorContract(basis, contract)); + + CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref)); + + return CEED_ERROR_SUCCESS; +} + //------------------------------------------------------------------------------ // Basis Destroy Tensor //------------------------------------------------------------------------------ diff --git a/backends/ref/ceed-ref-operator.c b/backends/ref/ceed-ref-operator.c index 2d7e1e20bd..6b0a061979 100644 --- a/backends/ref/ceed-ref-operator.c +++ b/backends/ref/ceed-ref-operator.c @@ -53,6 +53,7 @@ static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op, bool i case CEED_EVAL_INTERP: case CEED_EVAL_GRAD: case CEED_EVAL_DIV: + case CEED_EVAL_CURL: CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); @@ -68,8 +69,6 @@ static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op, bool i CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i])); break; - case CEED_EVAL_CURL: - break; // Not implemented } } return CEED_ERROR_SUCCESS; @@ -204,33 +203,16 @@ static inline int CeedOperatorInputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunction CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * Q * size])); break; case CEED_EVAL_INTERP: - CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); - CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); - CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp])); - CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, impl->e_vecs_in[i], impl->q_vecs_in[i])); - break; case CEED_EVAL_GRAD: - CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); - CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); - CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp])); - CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, impl->e_vecs_in[i], impl->q_vecs_in[i])); - break; case CEED_EVAL_DIV: + case CEED_EVAL_CURL: CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp])); - CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_DIV, impl->e_vecs_in[i], impl->q_vecs_in[i])); + CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs_in[i], impl->q_vecs_in[i])); break; case CEED_EVAL_WEIGHT: break; // No action - // LCOV_EXCL_START - case CEED_EVAL_CURL: { - CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); - Ceed ceed; - CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); - return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented"); - // LCOV_EXCL_STOP - } } } return CEED_ERROR_SUCCESS; @@ -257,25 +239,14 @@ static inline int CeedOperatorOutputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctio case CEED_EVAL_NONE: break; // No action case CEED_EVAL_INTERP: - CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); - CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); - CeedCallBackend( - CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp])); - CeedCallBackend(CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, impl->q_vecs_out[i], impl->e_vecs_out[i])); - break; case CEED_EVAL_GRAD: - CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); - CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); - CeedCallBackend( - CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp])); - CeedCallBackend(CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, impl->q_vecs_out[i], impl->e_vecs_out[i])); - break; case CEED_EVAL_DIV: + case CEED_EVAL_CURL: CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); CeedCallBackend( CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp])); - CeedCallBackend(CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_DIV, impl->q_vecs_out[i], impl->e_vecs_out[i])); + CeedCallBackend(CeedBasisApply(basis, 1, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i])); break; // LCOV_EXCL_START case CEED_EVAL_WEIGHT: { @@ -284,11 +255,6 @@ static inline int CeedOperatorOutputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctio return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output " "evaluation mode"); - } - case CEED_EVAL_CURL: { - Ceed ceed; - CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); - return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented"); // LCOV_EXCL_STOP } } diff --git a/backends/ref/ceed-ref.c b/backends/ref/ceed-ref.c index 3377a2213c..e052e1643d 100644 --- a/backends/ref/ceed-ref.c +++ b/backends/ref/ceed-ref.c @@ -26,6 +26,7 @@ static int CeedInit_Ref(const char *resource, Ceed ceed) { CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "BasisCreateTensorH1", CeedBasisCreateTensorH1_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "BasisCreateH1", CeedBasisCreateH1_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "BasisCreateHdiv", CeedBasisCreateHdiv_Ref)); + CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "BasisCreateHcurl", CeedBasisCreateHcurl_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "TensorContractCreate", CeedTensorContractCreate_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "ElemRestrictionCreate", CeedElemRestrictionCreate_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "ElemRestrictionCreateOriented", CeedElemRestrictionCreateOriented_Ref)); diff --git a/backends/ref/ceed-ref.h b/backends/ref/ceed-ref.h index 017110c468..325accfb13 100644 --- a/backends/ref/ceed-ref.h +++ b/backends/ref/ceed-ref.h @@ -70,6 +70,8 @@ CEED_INTERN int CeedBasisCreateH1_Ref(CeedElemTopology topo, CeedInt dim, CeedIn const CeedScalar *grad, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis); CEED_INTERN int CeedBasisCreateHdiv_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *div, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis); +CEED_INTERN int CeedBasisCreateHcurl_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, + const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis); CEED_INTERN int CeedTensorContractCreate_Ref(CeedBasis basis, CeedTensorContract contract); diff --git a/include/ceed-impl.h b/include/ceed-impl.h index 7ad39ab463..3b9a728cee 100644 --- a/include/ceed-impl.h +++ b/include/ceed-impl.h @@ -105,6 +105,8 @@ struct Ceed_private { CeedBasis); int (*BasisCreateHdiv)(CeedElemTopology, CeedInt, CeedInt, CeedInt, const CeedScalar *, const CeedScalar *, const CeedScalar *, const CeedScalar *, CeedBasis); + int (*BasisCreateHcurl)(CeedElemTopology, CeedInt, CeedInt, CeedInt, const CeedScalar *, const CeedScalar *, const CeedScalar *, const CeedScalar *, + CeedBasis); int (*TensorContractCreate)(CeedBasis, CeedTensorContract); int (*QFunctionCreate)(CeedQFunction); int (*QFunctionContextCreate)(CeedQFunctionContext); @@ -171,27 +173,28 @@ struct CeedBasis_private { Ceed ceed; int (*Apply)(CeedBasis, CeedInt, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector); int (*Destroy)(CeedBasis); - int ref_count; - bool tensor_basis; /* flag for tensor basis */ - CeedInt dim; /* topological dimension */ - CeedElemTopology topo; /* element topology */ - CeedInt num_comp; /* number of field components (1 for scalar fields) */ - CeedInt Q_comp; /* number of Q-vector components (1 for H^1, dim for H(div)) */ - CeedInt P_1d; /* number of nodes in one dimension */ - CeedInt Q_1d; /* number of quadrature points in one dimension */ - CeedInt P; /* total number of nodes */ - CeedInt Q; /* total number of quadrature points */ - CeedScalar *q_ref_1d; /* Array of length Q1d holding the locations of quadrature points on the 1D reference element [-1, 1] */ - CeedScalar *q_weight_1d; /* array of length Q1d holding the quadrature weights on the reference element */ - CeedScalar *interp; /* row-major matrix of shape [Q_comp*Q, P] expressing the values of nodal basis functions at quadrature points */ - CeedScalar *interp_1d; /* row-major matrix of shape [Q1d, P1d] expressing the values of nodal basis functions at quadrature points */ - CeedScalar *grad; /* row-major matrix of shape [dim*Q_comp*Q, P] matrix expressing derivatives of nodal basis functions at quadrature points */ - CeedScalar *grad_1d; /* row-major matrix of shape [Q1d, P1d] matrix expressing derivatives of nodal basis functions at quadrature points */ - CeedTensorContract contract; /* tensor contraction object */ - CeedFESpace basis_space; /* Initialize in basis constructor with 1,2 for H^1, H(div) FE space */ - CeedScalar - *div; /* row-major matrix of shape [Q, P] expressing the divergence of nodal basis functions at quadrature points for H(div) discretizations */ - void *data; /* place for the backend to store any data */ + int ref_count; + bool tensor_basis; /* flag for tensor basis */ + CeedInt dim; /* topological dimension */ + CeedElemTopology topo; /* element topology */ + CeedInt num_comp; /* number of field components (1 for scalar fields) */ + CeedInt P_1d; /* number of nodes in one dimension */ + CeedInt Q_1d; /* number of quadrature points in one dimension */ + CeedInt P; /* total number of nodes */ + CeedInt Q; /* total number of quadrature points */ + CeedFESpace fe_space; /* initialized in basis constructor with 1, 2, 3 for H^1, H(div), and H(curl) FE space */ + CeedTensorContract contract; /* tensor contraction object */ + CeedScalar *q_ref_1d; /* array of length Q1d holding the locations of quadrature points on the 1D reference element [-1, 1] */ + CeedScalar *q_weight_1d; /* array of length Q1d holding the quadrature weights on the reference element */ + CeedScalar *interp; /* row-major matrix of shape [Q, P] or [dim * Q, P] expressing the values of nodal basis functions or vector basis functions at + quadrature points */ + CeedScalar *interp_1d; /* row-major matrix of shape [Q1d, P1d] expressing the values of nodal basis functions at quadrature points */ + CeedScalar *grad; /* row-major matrix of shape [dim * Q, P] matrix expressing derivatives of nodal basis functions at quadrature points */ + CeedScalar *grad_1d; /* row-major matrix of shape [Q1d, P1d] matrix expressing derivatives of nodal basis functions at quadrature points */ + CeedScalar *div; /* row-major matrix of shape [Q, P] expressing the divergence of basis functions at quadrature points for H(div) discretizations */ + CeedScalar *curl; /* row-major matrix of shape [curl_dim * Q, P], curl_dim = 1 if dim < 3 else dim, expressing the curl of basis functions at + quadrature points for H(curl) discretizations */ + void *data; /* place for the backend to store any data */ }; struct CeedTensorContract_private { diff --git a/include/ceed/backend.h b/include/ceed/backend.h index c915d52d98..e2ef15a751 100644 --- a/include/ceed/backend.h +++ b/include/ceed/backend.h @@ -173,10 +173,12 @@ CEED_EXTERN int CeedElemRestrictionGetFlopsEstimate(CeedElemRestriction rstr, Ce /// Type of FE space; /// @ingroup CeedBasis typedef enum { - /// H1 FE space + /// H^1 FE space CEED_FE_SPACE_H1 = 1, /// H(div) FE space CEED_FE_SPACE_HDIV = 2, + /// H(curl) FE space + CEED_FE_SPACE_HCURL = 3, } CeedFESpace; CEED_EXTERN const char *const CeedFESpaces[]; @@ -185,15 +187,19 @@ CEED_EXTERN int CeedBasisIsTensor(CeedBasis basis, bool *is_tensor); CEED_EXTERN int CeedBasisGetData(CeedBasis basis, void *data); CEED_EXTERN int CeedBasisSetData(CeedBasis basis, void *data); CEED_EXTERN int CeedBasisReference(CeedBasis basis); +CEED_EXTERN int CeedBasisGetNumQuadratureComponents(CeedBasis basis, CeedEvalMode eval_mode, CeedInt *q_comp); CEED_EXTERN int CeedBasisGetFlopsEstimate(CeedBasis basis, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedSize *flops); - +CEED_EXTERN int CeedBasisGetFESpace(CeedBasis basis, CeedFESpace *fe_space); CEED_EXTERN int CeedBasisGetTopologyDimension(CeedElemTopology topo, CeedInt *dim); - CEED_EXTERN int CeedBasisGetTensorContract(CeedBasis basis, CeedTensorContract *contract); CEED_EXTERN int CeedBasisSetTensorContract(CeedBasis basis, CeedTensorContract contract); + CEED_EXTERN int CeedTensorContractCreate(Ceed ceed, CeedBasis basis, CeedTensorContract *contract); CEED_EXTERN int CeedTensorContractApply(CeedTensorContract contract, CeedInt A, CeedInt B, CeedInt C, CeedInt J, const CeedScalar *__restrict__ t, CeedTransposeMode t_mode, const CeedInt Add, const CeedScalar *__restrict__ u, CeedScalar *__restrict__ v); +CEED_EXTERN int CeedTensorContractStridedApply(CeedTensorContract contract, CeedInt A, CeedInt B, CeedInt C, CeedInt D, CeedInt J, + const CeedScalar *__restrict__ t, CeedTransposeMode t_mode, const CeedInt add, + const CeedScalar *__restrict__ u, CeedScalar *__restrict__ v); CEED_EXTERN int CeedTensorContractGetCeed(CeedTensorContract contract, Ceed *ceed); CEED_EXTERN int CeedTensorContractGetData(CeedTensorContract contract, void *data); CEED_EXTERN int CeedTensorContractSetData(CeedTensorContract contract, void *data); diff --git a/include/ceed/ceed.h b/include/ceed/ceed.h index ee064cdddf..e104f01417 100644 --- a/include/ceed/ceed.h +++ b/include/ceed/ceed.h @@ -274,6 +274,8 @@ CEED_EXTERN int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt num_ const CeedScalar *grad, const CeedScalar *q_ref, const CeedScalar *q_weights, CeedBasis *basis); CEED_EXTERN int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt nqpts, const CeedScalar *interp, const CeedScalar *div, const CeedScalar *q_ref, const CeedScalar *q_weights, CeedBasis *basis); +CEED_EXTERN int CeedBasisCreateHcurl(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt nqpts, const CeedScalar *interp, + const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weights, CeedBasis *basis); CEED_EXTERN int CeedBasisCreateProjection(CeedBasis basis_from, CeedBasis basis_to, CeedBasis *basis_project); CEED_EXTERN int CeedBasisReferenceCopy(CeedBasis basis, CeedBasis *basis_copy); CEED_EXTERN int CeedBasisView(CeedBasis basis, FILE *stream); @@ -281,7 +283,6 @@ CEED_EXTERN int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeM CEED_EXTERN int CeedBasisGetCeed(CeedBasis basis, Ceed *ceed); CEED_EXTERN int CeedBasisGetDimension(CeedBasis basis, CeedInt *dim); CEED_EXTERN int CeedBasisGetTopology(CeedBasis basis, CeedElemTopology *topo); -CEED_EXTERN int CeedBasisGetNumQuadratureComponents(CeedBasis basis, CeedInt *Q_comp); CEED_EXTERN int CeedBasisGetNumComponents(CeedBasis basis, CeedInt *num_comp); CEED_EXTERN int CeedBasisGetNumNodes(CeedBasis basis, CeedInt *P); CEED_EXTERN int CeedBasisGetNumNodes1D(CeedBasis basis, CeedInt *P_1d); @@ -294,6 +295,7 @@ CEED_EXTERN int CeedBasisGetInterp1D(CeedBasis basis, const CeedScalar **interp_ CEED_EXTERN int CeedBasisGetGrad(CeedBasis basis, const CeedScalar **grad); CEED_EXTERN int CeedBasisGetGrad1D(CeedBasis basis, const CeedScalar **grad_1d); CEED_EXTERN int CeedBasisGetDiv(CeedBasis basis, const CeedScalar **div); +CEED_EXTERN int CeedBasisGetCurl(CeedBasis basis, const CeedScalar **curl); CEED_EXTERN int CeedBasisDestroy(CeedBasis *basis); CEED_EXTERN int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d, CeedScalar *q_weight_1d); diff --git a/include/ceed/types.h b/include/ceed/types.h index 00e53fbd4d..cf9045d012 100644 --- a/include/ceed/types.h +++ b/include/ceed/types.h @@ -194,11 +194,11 @@ typedef enum { CEED_EVAL_NONE = 0, /// Interpolate from nodes to quadrature points CEED_EVAL_INTERP = 1, - /// Evaluate gradients at quadrature points from input in a nodal basis + /// Evaluate gradients at quadrature points from input in the basis CEED_EVAL_GRAD = 2, - /// Evaluate divergence at quadrature points from input in a nodal basis + /// Evaluate divergence at quadrature points from input in the basis CEED_EVAL_DIV = 4, - /// Evaluate curl at quadrature points from input in a nodal basis + /// Evaluate curl at quadrature points from input in the basis CEED_EVAL_CURL = 8, /// Using no input, evaluate quadrature weights on the reference element CEED_EVAL_WEIGHT = 16, @@ -213,7 +213,7 @@ typedef enum { CEED_GAUSS_LOBATTO = 1, } CeedQuadMode; -/// Type of basis shape to create non-tensor H1 element basis. +/// Type of basis shape to create non-tensor element basis. /// Dimension can be extracted with bitwise AND (CeedElemTopology & 2**(dim + 2)) == TRUE /// @ingroup CeedBasis typedef enum { diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c index 1338d84c42..ecf1310485 100644 --- a/interface/ceed-basis.c +++ b/interface/ceed-basis.c @@ -124,8 +124,9 @@ static int CeedScalarView(const char *name, const char *fp_fmt, CeedInt m, CeedI /** @brief Create the interpolation and gradient matrices for projection from the nodes of `basis_from` to the nodes of `basis_to`. - The interpolation is given by `interp_project = interp_to^+ * interp_from`, where the pesudoinverse `interp_to^+` is given by QR factorization. - The gradient is given by `grad_project = interp_to^+ * grad_from`. + The interpolation is given by `interp_project = interp_to^+ * interp_from`, where the pseudoinverse `interp_to^+` is given by QR factorization. + The gradient is given by `grad_project = interp_to^+ * grad_from`, and is only computed for H^1 spaces otherwise it should not be used. + Note: `basis_from` and `basis_to` must have compatible quadrature spaces. @param[in] basis_from CeedBasis to project from @@ -169,28 +170,62 @@ static int CeedBasisCreateProjectionMatrices(CeedBasis basis_from, CeedBasis bas // LCOV_EXCL_STOP } + // Check for matching FE space + CeedFESpace fe_space_to, fe_space_from; + CeedCall(CeedBasisGetFESpace(basis_to, &fe_space_to)); + CeedCall(CeedBasisGetFESpace(basis_from, &fe_space_from)); + if (fe_space_to != fe_space_from) { + // LCOV_EXCL_START + return CeedError(ceed, CEED_ERROR_MINOR, "Bases must both be the same FE space type"); + // LCOV_EXCL_STOP + } + // Get source matrices - CeedInt dim; - CeedScalar *interp_to, *interp_from, *tau; + CeedInt dim, q_comp = 1; + const CeedScalar *interp_to_source = NULL, *interp_from_source = NULL; + CeedScalar *interp_to, *interp_from, *tau; CeedCall(CeedBasisGetDimension(basis_to, &dim)); - CeedCall(CeedMalloc(Q * P_from, &interp_from)); - CeedCall(CeedMalloc(Q * P_to, &interp_to)); - CeedCall(CeedCalloc(P_to * P_from, interp_project)); - CeedCall(CeedCalloc(P_to * P_from * (is_tensor_to ? 1 : dim), grad_project)); - CeedCall(CeedMalloc(Q, &tau)); - const CeedScalar *interp_to_source = NULL, *interp_from_source = NULL, *grad_from_source; if (is_tensor_to) { CeedCall(CeedBasisGetInterp1D(basis_to, &interp_to_source)); CeedCall(CeedBasisGetInterp1D(basis_from, &interp_from_source)); - CeedCall(CeedBasisGetGrad1D(basis_from, &grad_from_source)); } else { + CeedCall(CeedBasisGetNumQuadratureComponents(basis_from, CEED_EVAL_INTERP, &q_comp)); CeedCall(CeedBasisGetInterp(basis_to, &interp_to_source)); CeedCall(CeedBasisGetInterp(basis_from, &interp_from_source)); - CeedCall(CeedBasisGetGrad(basis_from, &grad_from_source)); + } + CeedCall(CeedMalloc(Q * P_from * q_comp, &interp_from)); + CeedCall(CeedMalloc(Q * P_to * q_comp, &interp_to)); + CeedCall(CeedCalloc(P_to * P_from, interp_project)); + CeedCall(CeedMalloc(Q * q_comp, &tau)); + + // `grad_project = interp_to^+ * grad_from` is computed for the H^1 space case so the + // projection basis will have a gradient operation + const CeedScalar *grad_from_source = NULL; + if (fe_space_to == CEED_FE_SPACE_H1) { + if (is_tensor_to) { + CeedCall(CeedBasisGetGrad1D(basis_from, &grad_from_source)); + } else { + CeedCall(CeedBasisGetGrad(basis_from, &grad_from_source)); + } + CeedCall(CeedCalloc(P_to * P_from * (is_tensor_to ? 1 : dim), grad_project)); + } else { + // Projected derivate operator is only calculated for H^1 but we allocate it for the + // basis construction later + CeedInt q_comp_deriv = 1; + if (fe_space_to == CEED_FE_SPACE_HDIV) { + CeedCall(CeedBasisGetNumQuadratureComponents(basis_from, CEED_EVAL_DIV, &q_comp_deriv)); + } else if (fe_space_to == CEED_FE_SPACE_HCURL) { + CeedCall(CeedBasisGetNumQuadratureComponents(basis_from, CEED_EVAL_CURL, &q_comp_deriv)); + } + CeedCall(CeedCalloc(P_to * P_from * q_comp_deriv, grad_project)); } + // QR Factorization, interp_to = Q R + memcpy(interp_to, interp_to_source, Q * P_to * q_comp * sizeof(interp_to_source[0])); + CeedCall(CeedQRFactorization(ceed, interp_to, tau, Q * q_comp, P_to)); + // Build matrices - CeedInt num_matrices = 1 + (is_tensor_to ? 1 : dim); + CeedInt num_matrices = 1 + (fe_space_to == CEED_FE_SPACE_H1) * (is_tensor_to ? 1 : dim); CeedScalar *input_from[num_matrices], *output_project[num_matrices]; input_from[0] = (CeedScalar *)interp_from_source; output_project[0] = *interp_project; @@ -199,15 +234,11 @@ static int CeedBasisCreateProjectionMatrices(CeedBasis basis_from, CeedBasis bas output_project[m] = &((*grad_project)[(m - 1) * P_to * P_from]); } for (CeedInt m = 0; m < num_matrices; m++) { - // -- QR Factorization, interp_to = Q R - memcpy(interp_to, interp_to_source, Q * P_to * sizeof(interp_to_source[0])); - CeedCall(CeedQRFactorization(ceed, interp_to, tau, Q, P_to)); - - // -- Apply Qtranspose, interp_to = Qtranspose interp_from - memcpy(interp_from, input_from[m], Q * P_from * sizeof(input_from[m][0])); - CeedCall(CeedHouseholderApplyQ(interp_from, interp_to, tau, CEED_TRANSPOSE, Q, P_from, P_to, P_from, 1)); + // Apply Q^T, interp_from = Q^T interp_from + memcpy(interp_from, input_from[m], Q * P_from * q_comp * sizeof(input_from[m][0])); + CeedCall(CeedHouseholderApplyQ(interp_from, interp_to, tau, CEED_TRANSPOSE, Q * q_comp, P_from, P_to, P_from, 1)); - // -- Apply Rinv, interp_project = Rinv interp_c + // Apply Rinv, output_project = Rinv interp_from for (CeedInt j = 0; j < P_from; j++) { // Column j output_project[m][j + P_from * (P_to - 1)] = interp_from[j + P_from * (P_to - 1)] / interp_to[P_to * P_to - 1]; for (CeedInt i = P_to - 2; i >= 0; i--) { // Row i @@ -274,7 +305,7 @@ int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collo_grad_1d) { for (j = P_1d; j < Q_1d; j++) collo_grad_1d[j + Q_1d * i] = 0; } - // Apply Qtranspose, collo_grad = collo_grad Q_transpose + // Apply Q^T, collo_grad_1d = collo_grad_1d Q^T CeedCall(CeedHouseholderApplyQ(collo_grad_1d, interp_1d, tau, CEED_NOTRANSPOSE, Q_1d, Q_1d, P_1d, 1, Q_1d)); CeedCall(CeedFree(&interp_1d)); @@ -342,6 +373,42 @@ int CeedBasisReference(CeedBasis basis) { return CEED_ERROR_SUCCESS; } +/** + @brief Get number of Q-vector components for given CeedBasis + + @param[in] basis CeedBasis + @param[in] eval_mode \ref CEED_EVAL_INTERP to use interpolated values, + \ref CEED_EVAL_GRAD to use gradients, + \ref CEED_EVAL_DIV to use divergence, + \ref CEED_EVAL_CURL to use curl. + @param[out] q_comp Variable to store number of Q-vector components of basis + + @return An error code: 0 - success, otherwise - failure + + @ref Backend +**/ +int CeedBasisGetNumQuadratureComponents(CeedBasis basis, CeedEvalMode eval_mode, CeedInt *q_comp) { + switch (eval_mode) { + case CEED_EVAL_INTERP: + *q_comp = (basis->fe_space == CEED_FE_SPACE_H1) ? 1 : basis->dim; + break; + case CEED_EVAL_GRAD: + *q_comp = basis->dim; + break; + case CEED_EVAL_DIV: + *q_comp = 1; + break; + case CEED_EVAL_CURL: + *q_comp = (basis->dim < 3) ? 1 : basis->dim; + break; + case CEED_EVAL_NONE: + case CEED_EVAL_WEIGHT: + *q_comp = 1; + break; + } + return CEED_ERROR_SUCCESS; +} + /** @brief Estimate number of FLOPs required to apply CeedBasis in t_mode and eval_mode @@ -395,27 +462,21 @@ int CeedBasisGetFlopsEstimate(CeedBasis basis, CeedTransposeMode t_mode, CeedEva break; } } else { - CeedInt dim, num_comp, num_nodes, num_qpts, Q_comp; + CeedInt dim, num_comp, q_comp, num_nodes, num_qpts; CeedCall(CeedBasisGetDimension(basis, &dim)); CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); + CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); - CeedCall(CeedBasisGetNumQuadratureComponents(basis, &Q_comp)); switch (eval_mode) { case CEED_EVAL_NONE: *flops = 0; break; case CEED_EVAL_INTERP: - *flops = num_nodes * num_qpts * num_comp; - break; case CEED_EVAL_GRAD: - *flops = num_nodes * num_qpts * num_comp * dim; - break; case CEED_EVAL_DIV: - *flops = num_nodes * num_qpts; - break; case CEED_EVAL_CURL: - *flops = num_nodes * num_qpts * dim; + *flops = num_nodes * num_qpts * num_comp * q_comp; break; case CEED_EVAL_WEIGHT: *flops = 0; @@ -426,6 +487,21 @@ int CeedBasisGetFlopsEstimate(CeedBasis basis, CeedTransposeMode t_mode, CeedEva return CEED_ERROR_SUCCESS; } +/** + @brief Get CeedFESpace for a CeedBasis + + @param[in] basis CeedBasis + @param[out] fe_space Variable to store CeedFESpace + + @return An error code: 0 - success, otherwise - failure + + @ref Backend +**/ +int CeedBasisGetFESpace(CeedBasis basis, CeedFESpace *fe_space) { + *fe_space = basis->fe_space; + return CEED_ERROR_SUCCESS; +} + /** @brief Get dimension for given CeedElemTopology @@ -570,7 +646,7 @@ int CeedQRFactorization(Ceed ceed, CeedScalar *mat, CeedScalar *tau, CeedInt m, @return An error code: 0 - success, otherwise - failure - @ref Developer + @ref Utility **/ int CeedHouseholderApplyQ(CeedScalar *mat_A, const CeedScalar *mat_Q, const CeedScalar *tau, CeedTransposeMode t_mode, CeedInt m, CeedInt n, CeedInt k, CeedInt row, CeedInt col) { @@ -901,8 +977,7 @@ int CeedBasisCreateTensorH1(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P_ (*basis)->Q_1d = Q_1d; (*basis)->P = CeedIntPow(P_1d, dim); (*basis)->Q = CeedIntPow(Q_1d, dim); - (*basis)->Q_comp = 1; - (*basis)->basis_space = 1; // 1 for H^1 space + (*basis)->fe_space = CEED_FE_SPACE_H1; CeedCall(CeedCalloc(Q_1d, &(*basis)->q_ref_1d)); CeedCall(CeedCalloc(Q_1d, &(*basis)->q_weight_1d)); if (q_ref_1d) memcpy((*basis)->q_ref_1d, q_ref_1d, Q_1d * sizeof(q_ref_1d[0])); @@ -1020,7 +1095,7 @@ int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt num_comp, Ce @param[in] num_nodes Total number of nodes @param[in] num_qpts Total number of quadrature points @param[in] interp Row-major (num_qpts * num_nodes) matrix expressing the values of nodal basis functions at quadrature points - @param[in] grad Row-major (num_qpts * dim * num_nodes) matrix expressing derivatives of nodal basis functions at quadrature points + @param[in] grad Row-major (dim * num_qpts * num_nodes) matrix expressing derivatives of nodal basis functions at quadrature points @param[in] q_ref Array of length num_qpts * dim holding the locations of quadrature points on the reference element @param[in] q_weight Array of length num_qpts holding the quadrature weights on the reference element @param[out] basis Address of the variable where the newly created CeedBasis will be stored. @@ -1078,8 +1153,7 @@ int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedIn (*basis)->num_comp = num_comp; (*basis)->P = P; (*basis)->Q = Q; - (*basis)->Q_comp = 1; - (*basis)->basis_space = 1; // 1 for H^1 space + (*basis)->fe_space = CEED_FE_SPACE_H1; CeedCall(CeedCalloc(Q * dim, &(*basis)->q_ref_1d)); CeedCall(CeedCalloc(Q, &(*basis)->q_weight_1d)); if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); @@ -1100,8 +1174,8 @@ int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedIn @param[in] num_comp Number of components (usually 1 for vectors in H(div) bases) @param[in] num_nodes Total number of nodes (dofs per element) @param[in] num_qpts Total number of quadrature points - @param[in] interp Row-major (dim*num_qpts * num_nodes) matrix expressing the values of nodal basis functions at quadrature points - @param[in] div Row-major (num_qpts * num_nodes) matrix expressing divergence of nodal basis functions at quadrature points + @param[in] interp Row-major (dim * num_qpts * num_nodes) matrix expressing the values of basis functions at quadrature points + @param[in] div Row-major (num_qpts * num_nodes) matrix expressing divergence of basis functions at quadrature points @param[in] q_ref Array of length num_qpts * dim holding the locations of quadrature points on the reference element @param[in] q_weight Array of length num_qpts holding the quadrature weights on the reference element @param[out] basis Address of the variable where the newly created CeedBasis will be stored. @@ -1113,7 +1187,7 @@ int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedIn int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *div, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis *basis) { CeedInt Q = num_qpts, P = num_nodes, dim = 0; - CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); + if (!ceed->BasisCreateHdiv) { Ceed delegate; CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis")); @@ -1148,6 +1222,8 @@ int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, Ceed CeedCall(CeedCalloc(1, basis)); + CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); + (*basis)->ceed = ceed; CeedCall(CeedReference(ceed)); (*basis)->ref_count = 1; @@ -1157,8 +1233,7 @@ int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, Ceed (*basis)->num_comp = num_comp; (*basis)->P = P; (*basis)->Q = Q; - (*basis)->Q_comp = dim; - (*basis)->basis_space = 2; // 2 for H(div) space + (*basis)->fe_space = CEED_FE_SPACE_HDIV; CeedCall(CeedMalloc(Q * dim, &(*basis)->q_ref_1d)); CeedCall(CeedMalloc(Q, &(*basis)->q_weight_1d)); if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); @@ -1171,13 +1246,97 @@ int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, Ceed return CEED_ERROR_SUCCESS; } +/** + @brief Create a non tensor-product basis for H(curl) discretizations + + @param[in] ceed Ceed object where the CeedBasis will be created + @param[in] topo Topology of element (`CEED_TOPOLOGY_QUAD`, `CEED_TOPOLOGY_PRISM`, etc.), dimension of which is used in some array sizes below + @param[in] num_comp Number of components (usually 1 for vectors in H(curl) bases) + @param[in] num_nodes Total number of nodes (dofs per element) + @param[in] num_qpts Total number of quadrature points + @param[in] interp Row-major (dim * num_qpts * num_nodes) matrix expressing the values of basis functions at quadrature points + @param[in] curl Row-major (curl_comp * num_qpts * num_nodes, curl_comp = 1 if dim < 3 else dim) matrix expressing curl of basis functions at +quadrature points + @param[in] q_ref Array of length num_qpts * dim holding the locations of quadrature points on the reference element + @param[in] q_weight Array of length num_qpts holding the quadrature weights on the reference element + @param[out] basis Address of the variable where the newly created CeedBasis will be stored. + + @return An error code: 0 - success, otherwise - failure + + @ref User +**/ +int CeedBasisCreateHcurl(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, + const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis *basis) { + CeedInt Q = num_qpts, P = num_nodes, dim = 0, curl_comp = 0; + + if (!ceed->BasisCreateHdiv) { + Ceed delegate; + CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis")); + + if (!delegate) { + // LCOV_EXCL_START + return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateHcurl"); + // LCOV_EXCL_STOP + } + + CeedCall(CeedBasisCreateHcurl(delegate, topo, num_comp, num_nodes, num_qpts, interp, curl, q_ref, q_weight, basis)); + return CEED_ERROR_SUCCESS; + } + + if (num_comp < 1) { + // LCOV_EXCL_START + return CeedError(ceed, CEED_ERROR_DIMENSION, "Basis must have at least 1 component"); + // LCOV_EXCL_STOP + } + + if (num_nodes < 1) { + // LCOV_EXCL_START + return CeedError(ceed, CEED_ERROR_DIMENSION, "Basis must have at least 1 node"); + // LCOV_EXCL_STOP + } + + if (num_qpts < 1) { + // LCOV_EXCL_START + return CeedError(ceed, CEED_ERROR_DIMENSION, "Basis must have at least 1 quadrature point"); + // LCOV_EXCL_STOP + } + + CeedCall(CeedCalloc(1, basis)); + + CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); + curl_comp = (dim < 3) ? 1 : dim; + + (*basis)->ceed = ceed; + CeedCall(CeedReference(ceed)); + (*basis)->ref_count = 1; + (*basis)->tensor_basis = 0; + (*basis)->dim = dim; + (*basis)->topo = topo; + (*basis)->num_comp = num_comp; + (*basis)->P = P; + (*basis)->Q = Q; + (*basis)->fe_space = CEED_FE_SPACE_HCURL; + CeedCall(CeedMalloc(Q * dim, &(*basis)->q_ref_1d)); + CeedCall(CeedMalloc(Q, &(*basis)->q_weight_1d)); + if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); + if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0])); + CeedCall(CeedMalloc(dim * Q * P, &(*basis)->interp)); + CeedCall(CeedMalloc(curl_comp * Q * P, &(*basis)->curl)); + if (interp) memcpy((*basis)->interp, interp, dim * Q * P * sizeof(interp[0])); + if (curl) memcpy((*basis)->curl, curl, curl_comp * Q * P * sizeof(curl[0])); + CeedCall(ceed->BasisCreateHcurl(topo, dim, P, Q, interp, curl, q_ref, q_weight, *basis)); + return CEED_ERROR_SUCCESS; +} + /** @brief Create a CeedBasis for projection from the nodes of `basis_from` to the nodes of `basis_to`. - Only `CEED_EVAL_INTERP` and `CEED_EVAL_GRAD` will be valid for the new basis, `basis_project`. + Only `CEED_EVAL_INTERP` will be valid for the new basis, `basis_project`. For H^1 spaces, `CEED_EVAL_GRAD` will also be valid. The interpolation is given by `interp_project = interp_to^+ * interp_from`, where the pesudoinverse `interp_to^+` is given by QR -factorization. The gradient is given by `grad_project = interp_to^+ * grad_from`. Note: `basis_from` and `basis_to` must have compatible quadrature -spaces. +factorization. The gradient (for the H^1 case) is given by `grad_project = interp_to^+ * grad_from`. + + Note: `basis_from` and `basis_to` must have compatible quadrature spaces. + Note: `basis_project` will have the same number of components as `basis_from`, regardless of the number of components that `basis_to` has. If `basis_from` has 3 components and `basis_to` has 5 components, then `basis_project` will have 3 components. @@ -1199,9 +1358,11 @@ int CeedBasisCreateProjection(CeedBasis basis_from, CeedBasis basis_to, CeedBasi // Build basis bool is_tensor; + CeedFESpace fe_space; CeedInt dim, num_comp; CeedScalar *q_ref, *q_weight; CeedCall(CeedBasisIsTensor(basis_to, &is_tensor)); + CeedCall(CeedBasisGetFESpace(basis_to, &fe_space)); CeedCall(CeedBasisGetDimension(basis_to, &dim)); CeedCall(CeedBasisGetNumComponents(basis_from, &num_comp)); if (is_tensor) { @@ -1219,7 +1380,19 @@ int CeedBasisCreateProjection(CeedBasis basis_from, CeedBasis basis_to, CeedBasi CeedCall(CeedBasisGetNumNodes(basis_to, &num_nodes_to)); CeedCall(CeedCalloc(num_nodes_to * dim, &q_ref)); CeedCall(CeedCalloc(num_nodes_to, &q_weight)); - CeedCall(CeedBasisCreateH1(ceed, topo, num_comp, num_nodes_from, num_nodes_to, interp_project, grad_project, q_ref, q_weight, basis_project)); + switch (fe_space) { + case CEED_FE_SPACE_H1: + CeedCall(CeedBasisCreateH1(ceed, topo, num_comp, num_nodes_from, num_nodes_to, interp_project, grad_project, q_ref, q_weight, basis_project)); + break; + case CEED_FE_SPACE_HDIV: + CeedCall( + CeedBasisCreateHdiv(ceed, topo, num_comp, num_nodes_from, num_nodes_to, interp_project, grad_project, q_ref, q_weight, basis_project)); + break; + case CEED_FE_SPACE_HCURL: + CeedCall( + CeedBasisCreateHcurl(ceed, topo, num_comp, num_nodes_from, num_nodes_to, interp_project, grad_project, q_ref, q_weight, basis_project)); + break; + } } // Cleanup @@ -1262,15 +1435,16 @@ int CeedBasisReferenceCopy(CeedBasis basis, CeedBasis *basis_copy) { @ref User **/ int CeedBasisView(CeedBasis basis, FILE *stream) { - CeedFESpace FE_space = basis->basis_space; CeedElemTopology topo = basis->topo; + CeedFESpace fe_space = basis->fe_space; + CeedInt q_comp = 0; // Print FE space and element topology of the basis if (basis->tensor_basis) { - fprintf(stream, "CeedBasis (%s on a %s element): dim=%" CeedInt_FMT " P=%" CeedInt_FMT " Q=%" CeedInt_FMT "\n", CeedFESpaces[FE_space], + fprintf(stream, "CeedBasis (%s on a %s element): dim=%" CeedInt_FMT " P=%" CeedInt_FMT " Q=%" CeedInt_FMT "\n", CeedFESpaces[fe_space], CeedElemTopologies[topo], basis->dim, basis->P_1d, basis->Q_1d); } else { - fprintf(stream, "CeedBasis (%s on a %s element): dim=%" CeedInt_FMT " P=%" CeedInt_FMT " Q=%" CeedInt_FMT "\n", CeedFESpaces[FE_space], + fprintf(stream, "CeedBasis (%s on a %s element): dim=%" CeedInt_FMT " P=%" CeedInt_FMT " Q=%" CeedInt_FMT "\n", CeedFESpaces[fe_space], CeedElemTopologies[topo], basis->dim, basis->P, basis->Q); } // Print quadrature data, interpolation/gradient/divergence/curl of the basis @@ -1282,12 +1456,19 @@ int CeedBasisView(CeedBasis basis, FILE *stream) { } else { // non-tensor basis CeedCall(CeedScalarView("qref", "\t% 12.8f", 1, basis->Q * basis->dim, basis->q_ref_1d, stream)); CeedCall(CeedScalarView("qweight", "\t% 12.8f", 1, basis->Q, basis->q_weight_1d, stream)); - CeedCall(CeedScalarView("interp", "\t% 12.8f", basis->Q_comp * basis->Q, basis->P, basis->interp, stream)); + CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp)); + CeedCall(CeedScalarView("interp", "\t% 12.8f", q_comp * basis->Q, basis->P, basis->interp, stream)); if (basis->grad) { - CeedCall(CeedScalarView("grad", "\t% 12.8f", basis->dim * basis->Q, basis->P, basis->grad, stream)); + CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp)); + CeedCall(CeedScalarView("grad", "\t% 12.8f", q_comp * basis->Q, basis->P, basis->grad, stream)); } if (basis->div) { - CeedCall(CeedScalarView("div", "\t% 12.8f", basis->Q, basis->P, basis->div, stream)); + CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp)); + CeedCall(CeedScalarView("div", "\t% 12.8f", q_comp * basis->Q, basis->P, basis->div, stream)); + } + if (basis->curl) { + CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp)); + CeedCall(CeedScalarView("curl", "\t% 12.8f", q_comp * basis->Q, basis->P, basis->curl, stream)); } } return CEED_ERROR_SUCCESS; @@ -1304,6 +1485,8 @@ int CeedBasisView(CeedBasis basis, FILE *stream) { @param[in] eval_mode \ref CEED_EVAL_NONE to use values directly, \ref CEED_EVAL_INTERP to use interpolated values, \ref CEED_EVAL_GRAD to use gradients, + \ref CEED_EVAL_DIV to use divergence, + \ref CEED_EVAL_CURL to use curl, \ref CEED_EVAL_WEIGHT to use quadrature weights. @param[in] u Input CeedVector @param[out] v Output CeedVector @@ -1314,9 +1497,10 @@ int CeedBasisView(CeedBasis basis, FILE *stream) { **/ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) { CeedSize u_length = 0, v_length; - CeedInt dim, num_comp, num_nodes, num_qpts; + CeedInt dim, num_comp, q_comp, num_nodes, num_qpts; CeedCall(CeedBasisGetDimension(basis, &dim)); CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); + CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); CeedCall(CeedVectorGetLength(v, &v_length)); @@ -1343,26 +1527,15 @@ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, switch (eval_mode) { case CEED_EVAL_NONE: case CEED_EVAL_INTERP: - bad_dims = ((t_mode == CEED_TRANSPOSE && (u_length < num_elem * num_comp * num_qpts || v_length < num_elem * num_comp * num_nodes)) || - (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem * num_qpts * num_comp || u_length < num_elem * num_comp * num_nodes))); - break; case CEED_EVAL_GRAD: - bad_dims = ((t_mode == CEED_TRANSPOSE && (u_length < num_elem * num_comp * num_qpts * dim || v_length < num_elem * num_comp * num_nodes)) || - (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem * num_qpts * num_comp * dim || u_length < num_elem * num_comp * num_nodes))); + case CEED_EVAL_DIV: + case CEED_EVAL_CURL: + bad_dims = ((t_mode == CEED_TRANSPOSE && (u_length < num_elem * num_comp * num_qpts * q_comp || v_length < num_elem * num_comp * num_nodes)) || + (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem * num_qpts * num_comp * q_comp || u_length < num_elem * num_comp * num_nodes))); break; case CEED_EVAL_WEIGHT: bad_dims = v_length < num_elem * num_qpts; break; - // LCOV_EXCL_START - case CEED_EVAL_DIV: - bad_dims = ((t_mode == CEED_TRANSPOSE && (u_length < num_elem * num_comp * num_qpts || v_length < num_elem * num_comp * num_nodes)) || - (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem * num_qpts * num_comp || u_length < num_elem * num_comp * num_nodes))); - break; - case CEED_EVAL_CURL: - bad_dims = ((t_mode == CEED_TRANSPOSE && (u_length < num_elem * num_comp * num_qpts || v_length < num_elem * num_comp * num_nodes)) || - (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem * num_qpts * num_comp || u_length < num_elem * num_comp * num_nodes))); - break; - // LCOV_EXCL_STOP } if (bad_dims) { // LCOV_EXCL_START @@ -1419,21 +1592,6 @@ int CeedBasisGetTopology(CeedBasis basis, CeedElemTopology *topo) { return CEED_ERROR_SUCCESS; } -/** - @brief Get number of Q-vector components for given CeedBasis - - @param[in] basis CeedBasis - @param[out] Q_comp Variable to store number of Q-vector components of basis - - @return An error code: 0 - success, otherwise - failure - - @ref Advanced -**/ -int CeedBasisGetNumQuadratureComponents(CeedBasis basis, CeedInt *Q_comp) { - *Q_comp = basis->Q_comp; - return CEED_ERROR_SUCCESS; -} - /** @brief Get number of components for given CeedBasis @@ -1683,6 +1841,27 @@ int CeedBasisGetDiv(CeedBasis basis, const CeedScalar **div) { return CEED_ERROR_SUCCESS; } +/** + @brief Get curl matrix of a CeedBasis + + @param[in] basis CeedBasis + @param[out] curl Variable to store curl matrix + + @return An error code: 0 - success, otherwise - failure + + @ref Advanced +**/ +int CeedBasisGetCurl(CeedBasis basis, const CeedScalar **curl) { + if (!basis->curl) { + // LCOV_EXCL_START + return CeedError(basis->ceed, CEED_ERROR_MINOR, "CeedBasis does not have curl matrix."); + // LCOV_EXCL_STOP + } + + *curl = basis->curl; + return CEED_ERROR_SUCCESS; +} + /** @brief Destroy a CeedBasis @@ -1699,13 +1878,14 @@ int CeedBasisDestroy(CeedBasis *basis) { } if ((*basis)->Destroy) CeedCall((*basis)->Destroy(*basis)); if ((*basis)->contract) CeedCall(CeedTensorContractDestroy(&(*basis)->contract)); + CeedCall(CeedFree(&(*basis)->q_ref_1d)); + CeedCall(CeedFree(&(*basis)->q_weight_1d)); CeedCall(CeedFree(&(*basis)->interp)); CeedCall(CeedFree(&(*basis)->interp_1d)); CeedCall(CeedFree(&(*basis)->grad)); - CeedCall(CeedFree(&(*basis)->div)); CeedCall(CeedFree(&(*basis)->grad_1d)); - CeedCall(CeedFree(&(*basis)->q_ref_1d)); - CeedCall(CeedFree(&(*basis)->q_weight_1d)); + CeedCall(CeedFree(&(*basis)->div)); + CeedCall(CeedFree(&(*basis)->curl)); CeedCall(CeedDestroy(&(*basis)->ceed)); CeedCall(CeedFree(basis)); return CEED_ERROR_SUCCESS; diff --git a/interface/ceed-fortran.c b/interface/ceed-fortran.c index f1062ffe89..e97020b733 100644 --- a/interface/ceed-fortran.c +++ b/interface/ceed-fortran.c @@ -449,6 +449,40 @@ CEED_EXTERN void fCeedBasisCreateH1(int *ceed, int *topo, int *num_comp, int *nn } } +#define fCeedBasisCreateHdiv FORTRAN_NAME(ceedbasiscreatehdiv, CEEDBASISCREATEHDIV) +CEED_EXTERN void fCeedBasisCreateHdiv(int *ceed, int *topo, int *num_comp, int *nnodes, int *nqpts, const CeedScalar *interp, const CeedScalar *div, + const CeedScalar *qref, const CeedScalar *qweight, int *basis, int *err) { + if (CeedBasis_count == CeedBasis_count_max) { + CeedBasis_count_max += CeedBasis_count_max / 2 + 1; + CeedRealloc(CeedBasis_count_max, &CeedBasis_dict); + } + + *err = CeedBasisCreateHdiv(Ceed_dict[*ceed], (CeedElemTopology)*topo, *num_comp, *nnodes, *nqpts, interp, div, qref, qweight, + &CeedBasis_dict[CeedBasis_count]); + + if (*err == 0) { + *basis = CeedBasis_count++; + CeedBasis_n++; + } +} + +#define fCeedBasisCreateHcurl FORTRAN_NAME(ceedbasiscreatehcurl, CEEDBASISCREATEHCURL) +CEED_EXTERN void fCeedBasisCreateHcurl(int *ceed, int *topo, int *num_comp, int *nnodes, int *nqpts, const CeedScalar *interp, const CeedScalar *curl, + const CeedScalar *qref, const CeedScalar *qweight, int *basis, int *err) { + if (CeedBasis_count == CeedBasis_count_max) { + CeedBasis_count_max += CeedBasis_count_max / 2 + 1; + CeedRealloc(CeedBasis_count_max, &CeedBasis_dict); + } + + *err = CeedBasisCreateHcurl(Ceed_dict[*ceed], (CeedElemTopology)*topo, *num_comp, *nnodes, *nqpts, interp, curl, qref, qweight, + &CeedBasis_dict[CeedBasis_count]); + + if (*err == 0) { + *basis = CeedBasis_count++; + CeedBasis_n++; + } +} + #define fCeedBasisView FORTRAN_NAME(ceedbasisview, CEEDBASISVIEW) CEED_EXTERN void fCeedBasisView(int *basis, int *err) { *err = CeedBasisView(CeedBasis_dict[*basis], stdout); } diff --git a/interface/ceed-operator.c b/interface/ceed-operator.c index 8f2bf988eb..595ada819e 100644 --- a/interface/ceed-operator.c +++ b/interface/ceed-operator.c @@ -35,7 +35,7 @@ **/ static int CeedOperatorCheckField(Ceed ceed, CeedQFunctionField qf_field, CeedElemRestriction r, CeedBasis b) { CeedEvalMode eval_mode = qf_field->eval_mode; - CeedInt dim = 1, num_comp = 1, Q_comp = 1, restr_num_comp = 1, size = qf_field->size; + CeedInt dim = 1, num_comp = 1, q_comp = 1, restr_num_comp = 1, size = qf_field->size; // Restriction if (r != CEED_ELEMRESTRICTION_NONE) { @@ -61,7 +61,7 @@ static int CeedOperatorCheckField(Ceed ceed, CeedQFunctionField qf_field, CeedEl } CeedCall(CeedBasisGetDimension(b, &dim)); CeedCall(CeedBasisGetNumComponents(b, &num_comp)); - CeedCall(CeedBasisGetNumQuadratureComponents(b, &Q_comp)); + CeedCall(CeedBasisGetNumQuadratureComponents(b, eval_mode, &q_comp)); if (r != CEED_ELEMRESTRICTION_NONE && restr_num_comp != num_comp) { // LCOV_EXCL_START return CeedError(ceed, CEED_ERROR_DIMENSION, @@ -88,38 +88,19 @@ static int CeedOperatorCheckField(Ceed ceed, CeedQFunctionField qf_field, CeedEl } break; case CEED_EVAL_INTERP: - if (size != num_comp * Q_comp) { - // LCOV_EXCL_START - return CeedError(ceed, CEED_ERROR_DIMENSION, - "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: ElemRestriction/Basis has %" CeedInt_FMT " components", - qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], num_comp * Q_comp); - // LCOV_EXCL_STOP - } - break; case CEED_EVAL_GRAD: - if (size != num_comp * dim) { - // LCOV_EXCL_START - return CeedError(ceed, CEED_ERROR_DIMENSION, - "Field '%s' of size %" CeedInt_FMT " and EvalMode %s in %" CeedInt_FMT " dimensions: ElemRestriction/Basis has %" CeedInt_FMT - " components", - qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], dim, num_comp); - // LCOV_EXCL_STOP - } - break; - case CEED_EVAL_WEIGHT: - // No additional checks required - break; case CEED_EVAL_DIV: - if (size != num_comp) { + case CEED_EVAL_CURL: + if (size != num_comp * q_comp) { // LCOV_EXCL_START return CeedError(ceed, CEED_ERROR_DIMENSION, "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: ElemRestriction/Basis has %" CeedInt_FMT " components", - qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], num_comp); + qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], num_comp * q_comp); // LCOV_EXCL_STOP } break; - case CEED_EVAL_CURL: - // Not implemented + case CEED_EVAL_WEIGHT: + // No additional checks required break; } return CEED_ERROR_SUCCESS; diff --git a/interface/ceed-preconditioning.c b/interface/ceed-preconditioning.c index 6b4cbe4d1c..07ec08a7f7 100644 --- a/interface/ceed-preconditioning.c +++ b/interface/ceed-preconditioning.c @@ -180,32 +180,36 @@ int CeedOperatorGetFallback(CeedOperator op, CeedOperator *op_fallback) { /** @brief Select correct basis matrix pointer based on CeedEvalMode + @param[in] basis CeedBasis from which to get the basis matrix @param[in] eval_mode Current basis evaluation mode @param[in] identity Pointer to identity matrix - @param[in] interp Pointer to interpolation matrix - @param[in] grad Pointer to gradient matrix @param[out] basis_ptr Basis pointer to set @ref Developer **/ -static inline void CeedOperatorGetBasisPointer(CeedEvalMode eval_mode, const CeedScalar *identity, const CeedScalar *interp, const CeedScalar *grad, - const CeedScalar **basis_ptr) { +static inline int CeedOperatorGetBasisPointer(CeedBasis basis, CeedEvalMode eval_mode, const CeedScalar *identity, const CeedScalar **basis_ptr) { switch (eval_mode) { case CEED_EVAL_NONE: *basis_ptr = identity; break; case CEED_EVAL_INTERP: - *basis_ptr = interp; + CeedCall(CeedBasisGetInterp(basis, basis_ptr)); break; case CEED_EVAL_GRAD: - *basis_ptr = grad; + CeedCall(CeedBasisGetGrad(basis, basis_ptr)); break; - case CEED_EVAL_WEIGHT: case CEED_EVAL_DIV: + CeedCall(CeedBasisGetDiv(basis, basis_ptr)); + break; case CEED_EVAL_CURL: + CeedCall(CeedBasisGetCurl(basis, basis_ptr)); + break; + case CEED_EVAL_WEIGHT: break; // Caught by QF Assembly } assert(*basis_ptr != NULL); + + return CEED_ERROR_SUCCESS; } /** @@ -287,7 +291,6 @@ static inline int CeedSingleOperatorAssembleAddDiagonal_Core(CeedOperator op, Ce CeedSize **eval_mode_offsets_in, **eval_mode_offsets_out, num_output_components; CeedBasis *active_bases; CeedElemRestriction *active_elem_rstrs; - CeedCall(CeedOperatorGetOperatorAssemblyData(op, &data)); CeedCall(CeedOperatorAssemblyDataGetEvalModes(data, &num_active_bases, &num_eval_modes_in, &eval_modes_in, &eval_mode_offsets_in, &num_eval_modes_out, &eval_modes_out, &eval_mode_offsets_out, &num_output_components)); @@ -321,10 +324,9 @@ static inline int CeedSingleOperatorAssembleAddDiagonal_Core(CeedOperator op, Ce CeedCall(CeedBasisGetNumComponents(active_bases[b], &num_components)); CeedCall(CeedBasisGetNumQuadraturePoints(active_bases[b], &num_qpts)); - // Basis matrices - const CeedScalar *interp, *grad; - CeedScalar *identity = NULL; - bool has_eval_none = false; + // Construct identity matrix for basis if required + bool has_eval_none = false; + CeedScalar *identity = NULL; for (CeedInt i = 0; i < num_eval_modes_in[b]; i++) { has_eval_none = has_eval_none || (eval_modes_in[b][i] == CEED_EVAL_NONE); } @@ -335,22 +337,35 @@ static inline int CeedSingleOperatorAssembleAddDiagonal_Core(CeedOperator op, Ce CeedCall(CeedCalloc(num_qpts * num_nodes, &identity)); for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) identity[i * num_nodes + i] = 1.0; } - CeedCall(CeedBasisGetInterp(active_bases[b], &interp)); - CeedCall(CeedBasisGetGrad(active_bases[b], &grad)); + // Compute the diagonal of B^T D B // Each element for (CeedInt e = 0; e < num_elem; e++) { - CeedInt d_out = -1; // Each basis eval mode pair + CeedInt d_out = 0, q_comp_out; + CeedEvalMode eval_mode_out_prev = CEED_EVAL_NONE; for (CeedInt e_out = 0; e_out < num_eval_modes_out[b]; e_out++) { const CeedScalar *B_t = NULL; - if (eval_modes_out[b][e_out] == CEED_EVAL_GRAD) d_out += 1; - CeedOperatorGetBasisPointer(eval_modes_out[b][e_out], identity, interp, &grad[d_out * num_qpts * num_nodes], &B_t); - CeedInt d_in = -1; + CeedOperatorGetBasisPointer(active_bases[b], eval_modes_out[b][e_out], identity, &B_t); + CeedCall(CeedBasisGetNumQuadratureComponents(active_bases[b], eval_modes_out[b][e_out], &q_comp_out)); + if (q_comp_out > 1) { + if (e_out == 0 || eval_modes_out[b][e_out] != eval_mode_out_prev) d_out = 0; + else B_t = &B_t[(++d_out) * num_qpts * num_nodes]; + } + eval_mode_out_prev = eval_modes_out[b][e_out]; + + CeedInt d_in = 0, q_comp_in; + CeedEvalMode eval_mode_in_prev = CEED_EVAL_NONE; for (CeedInt e_in = 0; e_in < num_eval_modes_in[b]; e_in++) { const CeedScalar *B = NULL; - if (eval_modes_in[b][e_in] == CEED_EVAL_GRAD) d_in += 1; - CeedOperatorGetBasisPointer(eval_modes_in[b][e_in], identity, interp, &grad[d_in * num_qpts * num_nodes], &B); + CeedOperatorGetBasisPointer(active_bases[b], eval_modes_in[b][e_in], identity, &B); + CeedCall(CeedBasisGetNumQuadratureComponents(active_bases[b], eval_modes_in[b][e_in], &q_comp_in)); + if (q_comp_in > 1) { + if (e_in == 0 || eval_modes_in[b][e_in] != eval_mode_in_prev) d_in = 0; + else B = &B[(++d_in) * num_qpts * num_nodes]; + } + eval_mode_in_prev = eval_modes_in[b][e_in]; + // Each component for (CeedInt c_out = 0; c_out < num_components; c_out++) { // Each qpt/node pair @@ -603,10 +618,6 @@ static int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset, CeedVecto CeedInt local_num_entries = elem_size * num_comp * elem_size * num_comp * num_elem; // loop over elements and put in data structure - const CeedScalar *interp_in, *grad_in; - CeedCall(CeedBasisGetInterp(basis_in, &interp_in)); - CeedCall(CeedBasisGetGrad(basis_in, &grad_in)); - const CeedScalar *assembled_qf_array; CeedCall(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_HOST, &assembled_qf_array)); @@ -1151,26 +1162,23 @@ int CeedOperatorAssemblyDataCreate(Ceed ceed, CeedOperator op, CeedOperatorAssem CeedInt *num_eval_modes_in = NULL, *num_eval_modes_out = NULL, offset = 0; CeedEvalMode **eval_modes_in = NULL, **eval_modes_out = NULL; CeedSize **eval_mode_offsets_in = NULL, **eval_mode_offsets_out = NULL; - for (CeedInt i = 0; i < num_input_fields; i++) { CeedVector vec; - CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); if (vec == CEED_VECTOR_ACTIVE) { - CeedInt index = -1, dim = 1, num_components; CeedBasis basis_in = NULL; CeedEvalMode eval_mode; - + CeedInt index = -1, dim, num_comp, q_comp; CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis_in)); - CeedCall(CeedBasisGetDimension(basis_in, &dim)); - CeedCall(CeedBasisGetNumComponents(basis_in, &num_components)); CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); + CeedCall(CeedBasisGetDimension(basis_in, &dim)); + CeedCall(CeedBasisGetNumComponents(basis_in, &num_comp)); + CeedCall(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); for (CeedInt i = 0; i < num_active_bases; i++) { if ((*data)->active_bases[i] == basis_in) index = i; } if (index == -1) { CeedElemRestriction elem_rstr_in; - index = num_active_bases; CeedCall(CeedRealloc(num_active_bases + 1, &(*data)->active_bases)); (*data)->active_bases[num_active_bases] = NULL; @@ -1197,30 +1205,16 @@ int CeedOperatorAssemblyDataCreate(Ceed ceed, CeedOperator op, CeedOperatorAssem (*data)->assembled_bases_out[index] = NULL; num_active_bases++; } - switch (eval_mode) { - case CEED_EVAL_NONE: - case CEED_EVAL_INTERP: - CeedCall(CeedRealloc(num_eval_modes_in[index] + 1, &eval_modes_in[index])); - CeedCall(CeedRealloc(num_eval_modes_in[index] + 1, &eval_mode_offsets_in[index])); - eval_modes_in[index][num_eval_modes_in[index]] = eval_mode; - eval_mode_offsets_in[index][num_eval_modes_in[index]] = offset; - offset += num_components; - num_eval_modes_in[index] += 1; - break; - case CEED_EVAL_GRAD: - CeedCall(CeedRealloc(num_eval_modes_in[index] + dim, &eval_modes_in[index])); - CeedCall(CeedRealloc(num_eval_modes_in[index] + dim, &eval_mode_offsets_in[index])); - for (CeedInt d = 0; d < dim; d++) { - eval_modes_in[index][num_eval_modes_in[index] + d] = eval_mode; - eval_mode_offsets_in[index][num_eval_modes_in[index] + d] = offset; - offset += num_components; - } - num_eval_modes_in[index] += dim; - break; - case CEED_EVAL_WEIGHT: - case CEED_EVAL_DIV: - case CEED_EVAL_CURL: - break; // Caught by QF Assembly + if (eval_mode != CEED_EVAL_WEIGHT) { + // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly + CeedCall(CeedRealloc(num_eval_modes_in[index] + q_comp, &eval_modes_in[index])); + CeedCall(CeedRealloc(num_eval_modes_in[index] + q_comp, &eval_mode_offsets_in[index])); + for (CeedInt d = 0; d < q_comp; d++) { + eval_modes_in[index][num_eval_modes_in[index] + d] = eval_mode; + eval_mode_offsets_in[index][num_eval_modes_in[index] + d] = offset; + offset += num_comp; + } + num_eval_modes_in[index] += q_comp; } } } @@ -1230,23 +1224,21 @@ int CeedOperatorAssemblyDataCreate(Ceed ceed, CeedOperator op, CeedOperatorAssem // Determine active output basis CeedInt num_output_fields; - CeedCall(CeedQFunctionGetFields(qf, NULL, NULL, &num_output_fields, &qf_fields)); CeedCall(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); offset = 0; for (CeedInt i = 0; i < num_output_fields; i++) { CeedVector vec; - CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); if (vec == CEED_VECTOR_ACTIVE) { - CeedInt index = -1, dim = 1, num_components; CeedBasis basis_out = NULL; CeedEvalMode eval_mode; - + CeedInt index = -1, dim, num_comp, q_comp; CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis_out)); - CeedCall(CeedBasisGetDimension(basis_out, &dim)); - CeedCall(CeedBasisGetNumComponents(basis_out, &num_components)); CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); + CeedCall(CeedBasisGetDimension(basis_out, &dim)); + CeedCall(CeedBasisGetNumComponents(basis_out, &num_comp)); + CeedCall(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); for (CeedInt i = 0; i < num_active_bases; i++) { if ((*data)->active_bases[i] == basis_out) index = i; } @@ -1279,30 +1271,16 @@ int CeedOperatorAssemblyDataCreate(Ceed ceed, CeedOperator op, CeedOperatorAssem (*data)->assembled_bases_out[index] = NULL; num_active_bases++; } - switch (eval_mode) { - case CEED_EVAL_NONE: - case CEED_EVAL_INTERP: - CeedCall(CeedRealloc(num_eval_modes_out[index] + 1, &eval_modes_out[index])); - CeedCall(CeedRealloc(num_eval_modes_out[index] + 1, &eval_mode_offsets_out[index])); - eval_modes_out[index][num_eval_modes_out[index]] = eval_mode; - eval_mode_offsets_out[index][num_eval_modes_out[index]] = offset; - offset += num_components; - num_eval_modes_out[index] += 1; - break; - case CEED_EVAL_GRAD: - CeedCall(CeedRealloc(num_eval_modes_out[index] + dim, &eval_modes_out[index])); - CeedCall(CeedRealloc(num_eval_modes_out[index] + dim, &eval_mode_offsets_out[index])); - for (CeedInt d = 0; d < dim; d++) { - eval_modes_out[index][num_eval_modes_out[index] + d] = eval_mode; - eval_mode_offsets_out[index][num_eval_modes_out[index] + d] = offset; - offset += num_components; - } - num_eval_modes_out[index] += dim; - break; - case CEED_EVAL_WEIGHT: - case CEED_EVAL_DIV: - case CEED_EVAL_CURL: - break; // Caught by QF Assembly + if (eval_mode != CEED_EVAL_WEIGHT) { + // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly + CeedCall(CeedRealloc(num_eval_modes_out[index] + q_comp, &eval_modes_out[index])); + CeedCall(CeedRealloc(num_eval_modes_out[index] + q_comp, &eval_mode_offsets_out[index])); + for (CeedInt d = 0; d < q_comp; d++) { + eval_modes_out[index][num_eval_modes_out[index] + d] = eval_mode; + eval_mode_offsets_out[index][num_eval_modes_out[index] + d] = offset; + offset += num_comp; + } + num_eval_modes_out[index] += q_comp; } } } @@ -1375,36 +1353,38 @@ int CeedOperatorAssemblyDataGetBases(CeedOperatorAssemblyData data, CeedInt *num CeedCall(CeedBasisGetNumQuadraturePoints(data->active_bases[0], &num_qpts)); for (CeedInt b = 0; b < data->num_active_bases; b++) { - CeedInt elem_size; - CeedScalar *B_in = NULL, *identity = NULL; - const CeedScalar *interp_in, *grad_in; - bool has_eval_none = false; + CeedInt num_nodes; + CeedScalar *B_in = NULL, *identity = NULL; + bool has_eval_none = false; - CeedCall(CeedBasisGetNumNodes(data->active_bases[b], &elem_size)); - CeedCall(CeedCalloc(num_qpts * elem_size * data->num_eval_modes_in[b], &B_in)); + CeedCall(CeedBasisGetNumNodes(data->active_bases[b], &num_nodes)); + CeedCall(CeedCalloc(num_qpts * num_nodes * data->num_eval_modes_in[b], &B_in)); for (CeedInt i = 0; i < data->num_eval_modes_in[b]; i++) { has_eval_none = has_eval_none || (data->eval_modes_in[b][i] == CEED_EVAL_NONE); } if (has_eval_none) { - CeedCall(CeedCalloc(num_qpts * elem_size, &identity)); - for (CeedInt i = 0; i < (elem_size < num_qpts ? elem_size : num_qpts); i++) { - identity[i * elem_size + i] = 1.0; + CeedCall(CeedCalloc(num_qpts * num_nodes, &identity)); + for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) { + identity[i * num_nodes + i] = 1.0; } } - CeedCall(CeedBasisGetInterp(data->active_bases[b], &interp_in)); - CeedCall(CeedBasisGetGrad(data->active_bases[b], &grad_in)); for (CeedInt q = 0; q < num_qpts; q++) { - for (CeedInt n = 0; n < elem_size; n++) { - CeedInt d_in = -1; + for (CeedInt n = 0; n < num_nodes; n++) { + CeedInt d_in = 0, q_comp_in; + CeedEvalMode eval_mode_in_prev = CEED_EVAL_NONE; for (CeedInt e_in = 0; e_in < data->num_eval_modes_in[b]; e_in++) { const CeedInt qq = data->num_eval_modes_in[b] * q; const CeedScalar *B = NULL; - - if (data->eval_modes_in[b][e_in] == CEED_EVAL_GRAD) d_in++; - CeedOperatorGetBasisPointer(data->eval_modes_in[b][e_in], identity, interp_in, &grad_in[d_in * num_qpts * elem_size], &B); - B_in[(qq + e_in) * elem_size + n] = B[q * elem_size + n]; + CeedOperatorGetBasisPointer(data->active_bases[b], data->eval_modes_in[b][e_in], identity, &B); + CeedCall(CeedBasisGetNumQuadratureComponents(data->active_bases[b], data->eval_modes_in[b][e_in], &q_comp_in)); + if (q_comp_in > 1) { + if (e_in == 0 || data->eval_modes_in[b][e_in] != eval_mode_in_prev) d_in = 0; + else B = &B[(++d_in) * num_qpts * num_nodes]; + } + eval_mode_in_prev = data->eval_modes_in[b][e_in]; + B_in[(qq + e_in) * num_nodes + n] = B[q * num_nodes + n]; } } } @@ -1418,36 +1398,38 @@ int CeedOperatorAssemblyDataGetBases(CeedOperatorAssemblyData data, CeedInt *num CeedCall(CeedBasisGetNumQuadraturePoints(data->active_bases[0], &num_qpts)); for (CeedInt b = 0; b < data->num_active_bases; b++) { - CeedInt elem_size; - const CeedScalar *interp_out, *grad_out; - bool has_eval_none = false; - CeedScalar *B_out = NULL, *identity = NULL; + CeedInt num_nodes; + bool has_eval_none = false; + CeedScalar *B_out = NULL, *identity = NULL; - CeedCall(CeedBasisGetNumNodes(data->active_bases[b], &elem_size)); - CeedCall(CeedCalloc(num_qpts * elem_size * data->num_eval_modes_out[b], &B_out)); + CeedCall(CeedBasisGetNumNodes(data->active_bases[b], &num_nodes)); + CeedCall(CeedCalloc(num_qpts * num_nodes * data->num_eval_modes_out[b], &B_out)); for (CeedInt i = 0; i < data->num_eval_modes_out[b]; i++) { has_eval_none = has_eval_none || (data->eval_modes_out[b][i] == CEED_EVAL_NONE); } if (has_eval_none) { - CeedCall(CeedCalloc(num_qpts * elem_size, &identity)); - for (CeedInt i = 0; i < (elem_size < num_qpts ? elem_size : num_qpts); i++) { - identity[i * elem_size + i] = 1.0; + CeedCall(CeedCalloc(num_qpts * num_nodes, &identity)); + for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) { + identity[i * num_nodes + i] = 1.0; } } - CeedCall(CeedBasisGetInterp(data->active_bases[b], &interp_out)); - CeedCall(CeedBasisGetGrad(data->active_bases[b], &grad_out)); for (CeedInt q = 0; q < num_qpts; q++) { - for (CeedInt n = 0; n < elem_size; n++) { - CeedInt d_out = -1; + for (CeedInt n = 0; n < num_nodes; n++) { + CeedInt d_out = 0, q_comp_out; + CeedEvalMode eval_mode_out_prev = CEED_EVAL_NONE; for (CeedInt e_out = 0; e_out < data->num_eval_modes_out[b]; e_out++) { const CeedInt qq = data->num_eval_modes_out[b] * q; const CeedScalar *B = NULL; - - if (data->eval_modes_out[b][e_out] == CEED_EVAL_GRAD) d_out++; - CeedOperatorGetBasisPointer(data->eval_modes_out[b][e_out], identity, interp_out, &grad_out[d_out * num_qpts * elem_size], &B); - B_out[(qq + e_out) * elem_size + n] = B[q * elem_size + n]; + CeedOperatorGetBasisPointer(data->active_bases[b], data->eval_modes_out[b][e_out], identity, &B); + CeedCall(CeedBasisGetNumQuadratureComponents(data->active_bases[b], data->eval_modes_out[b][e_out], &q_comp_out)); + if (q_comp_out > 1) { + if (e_out == 0 || data->eval_modes_out[b][e_out] != eval_mode_out_prev) d_out = 0; + else B = &B[(++d_out) * num_qpts * num_nodes]; + } + eval_mode_out_prev = data->eval_modes_out[b][e_out]; + B_out[(qq + e_out) * num_nodes + n] = B[q * num_nodes + n]; } } } @@ -2354,9 +2336,9 @@ int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdm_inv, // LCOV_EXCL_STOP } CeedSize l_size = 1; - CeedInt P_1d, Q_1d, elem_size, num_qpts, dim, num_comp = 1, num_elem = 1; + CeedInt P_1d, Q_1d, num_nodes, num_qpts, dim, num_comp = 1, num_elem = 1; CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); - CeedCall(CeedBasisGetNumNodes(basis, &elem_size)); + CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); CeedCall(CeedBasisGetDimension(basis, &dim)); @@ -2439,26 +2421,26 @@ int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdm_inv, // Build FDM diagonal CeedVector q_data; CeedScalar *q_data_array, *fdm_diagonal; - CeedCall(CeedCalloc(num_comp * elem_size, &fdm_diagonal)); - const CeedScalar fdm_diagonal_bound = elem_size * CEED_EPSILON; + CeedCall(CeedCalloc(num_comp * num_nodes, &fdm_diagonal)); + const CeedScalar fdm_diagonal_bound = num_nodes * CEED_EPSILON; for (CeedInt c = 0; c < num_comp; c++) { - for (CeedInt n = 0; n < elem_size; n++) { - if (interp) fdm_diagonal[c * elem_size + n] = 1.0; + for (CeedInt n = 0; n < num_nodes; n++) { + if (interp) fdm_diagonal[c * num_nodes + n] = 1.0; if (grad) { for (CeedInt d = 0; d < dim; d++) { CeedInt i = (n / CeedIntPow(P_1d, d)) % P_1d; - fdm_diagonal[c * elem_size + n] += lambda[i]; + fdm_diagonal[c * num_nodes + n] += lambda[i]; } } - if (fabs(fdm_diagonal[c * elem_size + n]) < fdm_diagonal_bound) fdm_diagonal[c * elem_size + n] = fdm_diagonal_bound; + if (fabs(fdm_diagonal[c * num_nodes + n]) < fdm_diagonal_bound) fdm_diagonal[c * num_nodes + n] = fdm_diagonal_bound; } } - CeedCall(CeedVectorCreate(ceed_parent, num_elem * num_comp * elem_size, &q_data)); + CeedCall(CeedVectorCreate(ceed_parent, num_elem * num_comp * num_nodes, &q_data)); CeedCall(CeedVectorSetValue(q_data, 0.0)); CeedCall(CeedVectorGetArrayWrite(q_data, CEED_MEM_HOST, &q_data_array)); for (CeedInt e = 0; e < num_elem; e++) { for (CeedInt c = 0; c < num_comp; c++) { - for (CeedInt n = 0; n < elem_size; n++) q_data_array[(e * num_comp + c) * elem_size + n] = 1. / (elem_avg[e] * fdm_diagonal[c * elem_size + n]); + for (CeedInt n = 0; n < num_nodes; n++) q_data_array[(e * num_comp + c) * num_nodes + n] = 1. / (elem_avg[e] * fdm_diagonal[c * num_nodes + n]); } } CeedCall(CeedFree(&elem_avg)); @@ -2481,8 +2463,8 @@ int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdm_inv, // -- Restriction CeedElemRestriction rstr_qd_i; - CeedInt strides[3] = {1, elem_size, elem_size * num_comp}; - CeedCall(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, elem_size, num_comp, num_elem * num_comp * elem_size, strides, &rstr_qd_i)); + CeedInt strides[3] = {1, num_nodes, num_nodes * num_comp}; + CeedCall(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, num_nodes, num_comp, num_elem * num_comp * num_nodes, strides, &rstr_qd_i)); // -- QFunction CeedQFunction qf_fdm; CeedCall(CeedQFunctionCreateInteriorByName(ceed_parent, "Scale", &qf_fdm)); diff --git a/interface/ceed-qfunction.c b/interface/ceed-qfunction.c index d2e52ebfcc..b4e664a7d6 100644 --- a/interface/ceed-qfunction.c +++ b/interface/ceed-qfunction.c @@ -89,12 +89,16 @@ int CeedQFunctionRegister(const char *name, const char *source, CeedInt vec_leng @param[out] f CeedQFunctionField @param[in] field_name Name of QFunction field - @param[in] size Size of QFunction field, (num_comp * dim) for @ref CEED_EVAL_GRAD or (num_comp * 1) for @ref CEED_EVAL_NONE, @ref -CEED_EVAL_INTERP, and @ref CEED_EVAL_WEIGHT + @param[in] size Size of QFunction field, (num_comp * 1) for @ref CEED_EVAL_NONE and @ref CEED_EVAL_WEIGHT, +(num_comp * 1) for @ref CEED_EVAL_INTERP for an H^1 space or (num_comp * dim) for an H(div) or H(curl) space, +(num_comp * dim) for @ref CEED_EVAL_GRAD, or (num_comp * 1) for @ref CEED_EVAL_DIV, and +(num_comp * curl_dim) with curl_dim = 1 if dim < 3 else dim for @ref CEED_EVAL_CURL. @param[in] eval_mode \ref CEED_EVAL_NONE to use values directly, + \ref CEED_EVAL_WEIGHT to use quadrature weights, \ref CEED_EVAL_INTERP to use interpolated values, \ref CEED_EVAL_GRAD to use gradients, - \ref CEED_EVAL_WEIGHT to use quadrature weights. + \ref CEED_EVAL_DIV to use divergence, + \ref CEED_EVAL_CURL to use curl. @return An error code: 0 - success, otherwise - failure @@ -722,13 +726,17 @@ int CeedQFunctionReferenceCopy(CeedQFunction qf, CeedQFunction *qf_copy) { /** @brief Add a CeedQFunction input - @param[in,out] qf CeedQFunction - @param[in] field_name Name of QFunction field - @param[in] size Size of QFunction field, (num_comp * dim) for @ref CEED_EVAL_GRAD or (num_comp * 1) for @ref CEED_EVAL_NONE and @ref -CEED_EVAL_INTERP - @param[in] eval_mode \ref CEED_EVAL_NONE to use values directly, - \ref CEED_EVAL_INTERP to use interpolated values, - \ref CEED_EVAL_GRAD to use gradients. + @param[in,out] qf CeedQFunction + @param[in] field_name Name of QFunction field + @param[in] size Size of QFunction field, (num_comp * 1) for @ref CEED_EVAL_NONE, +(num_comp * 1) for @ref CEED_EVAL_INTERP for an H^1 space or (num_comp * dim) for an H(div) or H(curl) space, +(num_comp * dim) for @ref CEED_EVAL_GRAD, or (num_comp * 1) for @ref CEED_EVAL_DIV, and +(num_comp * curl_dim) with curl_dim = 1 if dim < 3 else dim for @ref CEED_EVAL_CURL. + @param[in] eval_mode \ref CEED_EVAL_NONE to use values directly, + \ref CEED_EVAL_INTERP to use interpolated values, + \ref CEED_EVAL_GRAD to use gradients, + \ref CEED_EVAL_DIV to use divergence, + \ref CEED_EVAL_CURL to use curl. @return An error code: 0 - success, otherwise - failure @@ -769,11 +777,15 @@ int CeedQFunctionAddInput(CeedQFunction qf, const char *field_name, CeedInt size @param[in,out] qf CeedQFunction @param[in] field_name Name of QFunction field - @param[in] size Size of QFunction field, (num_comp * dim) for @ref CEED_EVAL_GRAD or (num_comp * 1) for @ref CEED_EVAL_NONE and @ref -CEED_EVAL_INTERP + @param[in] size Size of QFunction field, (num_comp * 1) for @ref CEED_EVAL_NONE, +(num_comp * 1) for @ref CEED_EVAL_INTERP for an H^1 space or (num_comp * dim) for an H(div) or H(curl) space, +(num_comp * dim) for @ref CEED_EVAL_GRAD, or (num_comp * 1) for @ref CEED_EVAL_DIV, and +(num_comp * curl_dim) with curl_dim = 1 if dim < 3 else dim for @ref CEED_EVAL_CURL. @param[in] eval_mode \ref CEED_EVAL_NONE to use values directly, \ref CEED_EVAL_INTERP to use interpolated values, - \ref CEED_EVAL_GRAD to use gradients. + \ref CEED_EVAL_GRAD to use gradients, + \ref CEED_EVAL_DIV to use divergence, + \ref CEED_EVAL_CURL to use curl. @return An error code: 0 - success, otherwise - failure diff --git a/interface/ceed-tensor.c b/interface/ceed-tensor.c index 6ed6cea00c..a8e1b21bf0 100644 --- a/interface/ceed-tensor.c +++ b/interface/ceed-tensor.c @@ -82,6 +82,44 @@ int CeedTensorContractApply(CeedTensorContract contract, CeedInt A, CeedInt B, C return CEED_ERROR_SUCCESS; } +/** + @brief Apply tensor contraction + + Contracts on the middle index + NOTRANSPOSE: v_dajc = t_djb u_abc + TRANSPOSE: v_ajc = t_dbj u_dabc + If add != 0, "=" is replaced by "+=" + + @param[in] contract CeedTensorContract to use + @param[in] A First index of u, second index of v + @param[in] B Middle index of u, one of last two indices of t + @param[in] C Last index of u, v + @param[in] D First index of v, first index of t + @param[in] J Third index of v, one of last two indices of t + @param[in] t Tensor array to contract against + @param[in] t_mode Transpose mode for t, \ref CEED_NOTRANSPOSE for t_jb \ref CEED_TRANSPOSE for t_bj + @param[in] add Add mode + @param[in] u Input array + @param[out] v Output array + + @return An error code: 0 - success, otherwise - failure + + @ref Backend +**/ +int CeedTensorContractStridedApply(CeedTensorContract contract, CeedInt A, CeedInt B, CeedInt C, CeedInt D, CeedInt J, const CeedScalar *restrict t, + CeedTransposeMode t_mode, const CeedInt add, const CeedScalar *restrict u, CeedScalar *restrict v) { + if (t_mode == CEED_TRANSPOSE) { + for (CeedInt d = 0; d < D; d++) { + CeedCall(contract->Apply(contract, A, J, C, B, t + d * B * J, t_mode, add, u + d * A * J * C, v)); + } + } else { + for (CeedInt d = 0; d < D; d++) { + CeedCall(contract->Apply(contract, A, B, C, J, t + d * B * J, t_mode, add, u, v + d * A * J * C)); + } + } + return CEED_ERROR_SUCCESS; +} + /** @brief Get Ceed associated with a CeedTensorContract diff --git a/interface/ceed-types.c b/interface/ceed-types.c index 1b362084ce..0e603bead8 100644 --- a/interface/ceed-types.c +++ b/interface/ceed-types.c @@ -59,6 +59,7 @@ const char *const CeedContextFieldTypes[] = { }; const char *const CeedFESpaces[] = { - [CEED_FE_SPACE_H1] = "H^1 space", - [CEED_FE_SPACE_HDIV] = "H(div) space", + [CEED_FE_SPACE_H1] = "H^1 space", + [CEED_FE_SPACE_HDIV] = "H(div) space", + [CEED_FE_SPACE_HCURL] = "H(curl) space", }; diff --git a/interface/ceed.c b/interface/ceed.c index 1d5eb382f9..af49b3ff03 100644 --- a/interface/ceed.c +++ b/interface/ceed.c @@ -820,6 +820,7 @@ int CeedInit(const char *resource, Ceed *ceed) { CEED_FTABLE_ENTRY(Ceed, BasisCreateTensorH1), CEED_FTABLE_ENTRY(Ceed, BasisCreateH1), CEED_FTABLE_ENTRY(Ceed, BasisCreateHdiv), + CEED_FTABLE_ENTRY(Ceed, BasisCreateHcurl), CEED_FTABLE_ENTRY(Ceed, TensorContractCreate), CEED_FTABLE_ENTRY(Ceed, QFunctionCreate), CEED_FTABLE_ENTRY(Ceed, QFunctionContextCreate), diff --git a/julia/LibCEED.jl/docs/src/Basis.md b/julia/LibCEED.jl/docs/src/Basis.md index 1a482991a6..cbf690b68a 100644 --- a/julia/LibCEED.jl/docs/src/Basis.md +++ b/julia/LibCEED.jl/docs/src/Basis.md @@ -14,11 +14,13 @@ BasisCollocated create_tensor_h1_lagrange_basis create_tensor_h1_basis create_h1_basis +create_hdiv_basis +create_hcurl_basis apply!(b::Basis, nelem, tmode::TransposeMode, emode::EvalMode, u::LibCEED.AbstractCeedVector, v::LibCEED.AbstractCeedVector) apply(b::Basis, u::AbstractVector; nelem=1, tmode=NOTRANSPOSE, emode=EVAL_INTERP) getdimension gettopology -getnumcomponents(b::Basis) +getnumcomponents getnumnodes getnumnodes1d getnumqpts @@ -29,4 +31,6 @@ getinterp getinterp1d getgrad getgrad1d +getdiv +getcurl ``` diff --git a/julia/LibCEED.jl/src/Basis.jl b/julia/LibCEED.jl/src/Basis.jl index 05b25c7c35..6769b72902 100644 --- a/julia/LibCEED.jl/src/Basis.jl +++ b/julia/LibCEED.jl/src/Basis.jl @@ -17,6 +17,8 @@ created using one of: - [`create_tensor_h1_lagrange_basis`](@ref) - [`create_tensor_h1_basis`](@ref) - [`create_h1_basis`](@ref) +- [`create_hdiv_basis`](@ref) +- [`create_hcurl_basis`](@ref) """ mutable struct Basis <: AbstractBasis ref::RefValue{C.CeedBasis} @@ -112,7 +114,7 @@ end @doc raw""" create_h1_basis(c::Ceed, topo::Topology, ncomp, nnodes, nqpts, interp, grad, qref, qweight) -Create a non tensor-product basis for H^1 discretizations +Create a non tensor-product basis for $H^1$ discretizations # Arguments: - `ceed`: A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created. @@ -124,7 +126,7 @@ Create a non tensor-product basis for H^1 discretizations at quadrature points. - `grad`: Array of size `(dim, nqpts, nnodes)` expressing derivatives of nodal basis functions at quadrature points. -- `qref`: Array of length `nqpts` holding the locations of quadrature points on the +- `qref`: Matrix of size `(dim, nqpts)` holding the locations of quadrature points on the reference element $[-1, 1]$. - `qweight`: Array of length `nqpts` holding the quadrature weights on the reference element. @@ -143,12 +145,13 @@ function create_h1_basis( dim = getdimension(topo) @assert size(interp) == (nqpts, nnodes) @assert size(grad) == (dim, nqpts, nnodes) - @assert length(qref) == nqpts + @assert size(qref) == (dim, nqpts) @assert length(qweight) == nqpts # Convert from Julia matrices and tensors (column-major) to row-major format interp_rowmajor = collect(interp') grad_rowmajor = permutedims(grad, [3, 2, 1]) + qref_rowmajor = collect(qref') ref = Ref{C.CeedBasis}() C.CeedBasisCreateH1( @@ -159,7 +162,124 @@ function create_h1_basis( nqpts, interp_rowmajor, grad_rowmajor, - qref, + qref_rowmajor, + qweight, + ref, + ) + Basis(ref) +end + +@doc raw""" + create_hdiv_basis(c::Ceed, topo::Topology, ncomp, nnodes, nqpts, interp, div, qref, qweight) + +Create a non tensor-product basis for H(div) discretizations + +# Arguments: +- `ceed`: A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created. +- `topo`: [`Topology`](@ref) of element, e.g. hypercube, simplex, etc. +- `ncomp`: Number of field components (1 for scalar fields). +- `nnodes`: Total number of nodes. +- `nqpts`: Total number of quadrature points. +- `interp`: Matrix of size `(dim, nqpts, nnodes)` expressing the values of basis functions + at quadrature points. +- `div`: Array of size `(nqpts, nnodes)` expressing divergence of basis functions at + quadrature points. +- `qref`: Matrix of size `(dim, nqpts)` holding the locations of quadrature points on the + reference element $[-1, 1]$. +- `qweight`: Array of length `nqpts` holding the quadrature weights on the reference + element. +""" +function create_hdiv_basis( + c::Ceed, + topo::Topology, + ncomp, + nnodes, + nqpts, + interp::AbstractArray{CeedScalar}, + div::AbstractArray{CeedScalar}, + qref::AbstractArray{CeedScalar}, + qweight::AbstractArray{CeedScalar}, +) + dim = getdimension(topo) + @assert size(interp) == (dim, nqpts, nnodes) + @assert size(div) == (nqpts, nnodes) + @assert size(qref) == (dim, nqpts) + @assert length(qweight) == nqpts + + # Convert from Julia matrices and tensors (column-major) to row-major format + interp_rowmajor = permutedims(interp, [3, 2, 1]) + div_rowmajor = collect(div') + qref_rowmajor = collect(qref') + + ref = Ref{C.CeedBasis}() + C.CeedBasisCreateHdiv( + c[], + topo, + ncomp, + nnodes, + nqpts, + interp_rowmajor, + div_rowmajor, + qref_rowmajor, + qweight, + ref, + ) + Basis(ref) +end + +@doc raw""" + create_hcurl_basis(c::Ceed, topo::Topology, ncomp, nnodes, nqpts, interp, curl, qref, qweight) + +Create a non tensor-product basis for H(curl) discretizations + +# Arguments: +- `ceed`: A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created. +- `topo`: [`Topology`](@ref) of element, e.g. hypercube, simplex, etc. +- `ncomp`: Number of field components (1 for scalar fields). +- `nnodes`: Total number of nodes. +- `nqpts`: Total number of quadrature points. +- `interp`: Matrix of size `(dim, nqpts, nnodes)` expressing the values of basis functions + at quadrature points. +- `curl`: Matrix of size `(curlcomp, nqpts, nnodes)`, `curlcomp = 1 if dim < 3 else dim`) + matrix expressing curl of basis functions at quadrature points. +- `qref`: Matrix of size `(dim, nqpts)` holding the locations of quadrature points on the + reference element $[-1, 1]$. +- `qweight`: Array of length `nqpts` holding the quadrature weights on the reference + element. +""" +function create_hcurl_basis( + c::Ceed, + topo::Topology, + ncomp, + nnodes, + nqpts, + interp::AbstractArray{CeedScalar}, + curl::AbstractArray{CeedScalar}, + qref::AbstractArray{CeedScalar}, + qweight::AbstractArray{CeedScalar}, +) + dim = getdimension(topo) + curlcomp = dim < 3 ? 1 : dim + @assert size(interp) == (dim, nqpts, nnodes) + @assert size(curl) == (curlcomp, nqpts, nnodes) + @assert size(qref) == (dim, nqpts) + @assert length(qweight) == nqpts + + # Convert from Julia matrices and tensors (column-major) to row-major format + interp_rowmajor = permutedims(interp, [3, 2, 1]) + curl_rowmajor = permutedims(curl, [3, 2, 1]) + qref_rowmajor = collect(qref') + + ref = Ref{C.CeedBasis}() + C.CeedBasisCreateHcurl( + c[], + topo, + ncomp, + nnodes, + nqpts, + interp_rowmajor, + curl_rowmajor, + qref_rowmajor, qweight, ref, ) @@ -211,10 +331,9 @@ function apply(b::Basis, u::AbstractVector; nelem=1, tmode=NOTRANSPOSE, emode=EV u_vec = CeedVector(c, u) - len_v = (tmode == TRANSPOSE) ? getnumnodes(b) : getnumqpts(b) - if emode == EVAL_GRAD - len_v *= getdimension(b) - end + qcomp = Ref{CeedInt}() + C.CeedBasisGetNumQuadratureComponents(b[], emode, qcomp) + len_v = (tmode == TRANSPOSE) ? getnumnodes(b) : qcomp[]*getnumqpts(b) v_vec = CeedVector(c, len_v) @@ -320,11 +439,8 @@ function getqref(b::Basis) ref = Ref{Ptr{CeedScalar}}() C.CeedBasisGetQRef(b[], ref) copy( - unsafe_wrap( - Array, - ref[], - istensor[] ? getnumqpts1d(b) : (getnumqpts(b)*getdimension(b)), - ), + istensor[] ? unsafe_wrap(Array, ref[], getnumqpts1d(b)) : + unsafe_wrap(Array, ref[], (getnumqpts(b), getdimension(b)))', ) end @@ -342,18 +458,25 @@ function getqweights(b::Basis) copy(unsafe_wrap(Array, ref[], istensor[] ? getnumqpts1d(b) : getnumqpts(b))) end -""" +@doc raw""" getinterp(b::Basis) Get the interpolation matrix of the given [`Basis`](@ref). Returns a matrix of size -`(getnumqpts(b), getnumnodes(b))`. +`(getnumqpts(b), getnumnodes(b))` for a given $H^1$ basis or `(getdimension(b), +getnumqpts(b), getnumnodes(b))` for a given vector $H(div)$ or $H(curl)$ basis. """ function getinterp(b::Basis) ref = Ref{Ptr{CeedScalar}}() C.CeedBasisGetInterp(b[], ref) q = getnumqpts(b) p = getnumnodes(b) - collect(unsafe_wrap(Array, ref[], (p, q))') + qcomp = Ref{CeedInt}() + C.CeedBasisGetNumQuadratureComponents(b[], C.CEED_EVAL_INTERP, qcomp) + if qcomp[] == 1 + collect(unsafe_wrap(Array, ref[], (p, q))') + else + permutedims(unsafe_wrap(Array, ref[], (p, q, qcomp[])), [3, 2, 1]) + end end """ @@ -399,3 +522,34 @@ function getgrad1d(b::Basis) p = getnumnodes1d(b) collect(unsafe_wrap(Array, ref[], (p, q))') end + +""" + getdiv(b::Basis) + +Get the divergence matrix of the given [`Basis`](@ref). Returns a tensor of size +`(getnumqpts(b), getnumnodes(b))`. +""" +function getdiv(b::Basis) + ref = Ref{Ptr{CeedScalar}}() + C.CeedBasisGetDiv(b[], ref) + q = getnumqpts(b) + p = getnumnodes(b) + collect(unsafe_wrap(Array, ref[], (p, q))') +end + +""" + getcurl(b::Basis) + +Get the curl matrix of the given [`Basis`](@ref). Returns a tensor of size +`(curlcomp, getnumqpts(b), getnumnodes(b))`, `curlcomp = 1 if getdimension(b) < 3 else +getdimension(b)`. +""" +function getcurl(b::Basis) + ref = Ref{Ptr{CeedScalar}}() + C.CeedBasisGetCurl(b[], ref) + q = getnumqpts(b) + p = getnumnodes(b) + qcomp = Ref{CeedInt}() + C.CeedBasisGetNumQuadratureComponents(b[], C.CEED_EVAL_CURL, qcomp) + permutedims(unsafe_wrap(Array, ref[], (p, q, qcomp[])), [3, 2, 1]) +end diff --git a/julia/LibCEED.jl/src/LibCEED.jl b/julia/LibCEED.jl/src/LibCEED.jl index aaa8d5d88d..fd127aa051 100644 --- a/julia/LibCEED.jl/src/LibCEED.jl +++ b/julia/LibCEED.jl/src/LibCEED.jl @@ -79,6 +79,8 @@ export @interior_qf, create_elem_restriction_strided, create_evector, create_h1_basis, + create_hdiv_basis, + create_hcurl_basis, create_identity_qfunction, create_interior_qfunction, create_lvector, @@ -150,7 +152,7 @@ include("Request.jl") include("Operator.jl") include("Misc.jl") -const minimum_libceed_version = v"0.10.0" +const minimum_libceed_version = v"0.12.0" function __init__() if !ceedversion_ge(minimum_libceed_version) diff --git a/julia/LibCEED.jl/src/generated/libceed_bindings.jl b/julia/LibCEED.jl/src/generated/libceed_bindings.jl index 7e954c2d40..ee71bd8974 100644 --- a/julia/LibCEED.jl/src/generated/libceed_bindings.jl +++ b/julia/LibCEED.jl/src/generated/libceed_bindings.jl @@ -24,6 +24,57 @@ end CEED_ERROR_UNSUPPORTED = -3 end +@cenum CeedMemType::UInt32 begin + CEED_MEM_HOST = 0 + CEED_MEM_DEVICE = 1 +end + +@cenum CeedCopyMode::UInt32 begin + CEED_COPY_VALUES = 0 + CEED_USE_POINTER = 1 + CEED_OWN_POINTER = 2 +end + +@cenum CeedNormType::UInt32 begin + CEED_NORM_1 = 0 + CEED_NORM_2 = 1 + CEED_NORM_MAX = 2 +end + +@cenum CeedTransposeMode::UInt32 begin + CEED_NOTRANSPOSE = 0 + CEED_TRANSPOSE = 1 +end + +@cenum CeedEvalMode::UInt32 begin + CEED_EVAL_NONE = 0 + CEED_EVAL_INTERP = 1 + CEED_EVAL_GRAD = 2 + CEED_EVAL_DIV = 4 + CEED_EVAL_CURL = 8 + CEED_EVAL_WEIGHT = 16 +end + +@cenum CeedQuadMode::UInt32 begin + CEED_GAUSS = 0 + CEED_GAUSS_LOBATTO = 1 +end + +@cenum CeedElemTopology::UInt32 begin + CEED_TOPOLOGY_LINE = 65536 + CEED_TOPOLOGY_TRIANGLE = 131073 + CEED_TOPOLOGY_QUAD = 131074 + CEED_TOPOLOGY_TET = 196611 + CEED_TOPOLOGY_PYRAMID = 196612 + CEED_TOPOLOGY_PRISM = 196613 + CEED_TOPOLOGY_HEX = 196614 +end + +@cenum CeedContextFieldType::UInt32 begin + CEED_CONTEXT_FIELD_DOUBLE = 1 + CEED_CONTEXT_FIELD_INT32 = 2 +end + mutable struct Ceed_private end const Ceed = Ptr{Ceed_private} @@ -123,27 +174,10 @@ function CeedGetScalarType(scalar_type) ccall((:CeedGetScalarType, libceed), Cint, (Ptr{CeedScalarType},), scalar_type) end -@cenum CeedMemType::UInt32 begin - CEED_MEM_HOST = 0 - CEED_MEM_DEVICE = 1 -end - function CeedGetPreferredMemType(ceed, type) ccall((:CeedGetPreferredMemType, libceed), Cint, (Ceed, Ptr{CeedMemType}), ceed, type) end -@cenum CeedCopyMode::UInt32 begin - CEED_COPY_VALUES = 0 - CEED_USE_POINTER = 1 - CEED_OWN_POINTER = 2 -end - -@cenum CeedNormType::UInt32 begin - CEED_NORM_1 = 0 - CEED_NORM_2 = 1 - CEED_NORM_MAX = 2 -end - function CeedVectorCreate(ceed, len, vec) ccall((:CeedVectorCreate, libceed), Cint, (Ceed, CeedSize, Ptr{CeedVector}), ceed, len, vec) end @@ -152,6 +186,10 @@ function CeedVectorReferenceCopy(vec, vec_copy) ccall((:CeedVectorReferenceCopy, libceed), Cint, (CeedVector, Ptr{CeedVector}), vec, vec_copy) end +function CeedVectorCopy(vec, vec_copy) + ccall((:CeedVectorCopy, libceed), Cint, (CeedVector, CeedVector), vec, vec_copy) +end + function CeedVectorSetArray(vec, mem_type, copy_mode, array) ccall((:CeedVectorSetArray, libceed), Cint, (CeedVector, CeedMemType, CeedCopyMode, Ptr{CeedScalar}), vec, mem_type, copy_mode, array) end @@ -200,6 +238,10 @@ function CeedVectorAXPY(y, alpha, x) ccall((:CeedVectorAXPY, libceed), Cint, (CeedVector, CeedScalar, CeedVector), y, alpha, x) end +function CeedVectorAXPBY(y, alpha, beta, x) + ccall((:CeedVectorAXPBY, libceed), Cint, (CeedVector, CeedScalar, CeedScalar, CeedVector), y, alpha, beta, x) +end + function CeedVectorPointwiseMult(w, x, y) ccall((:CeedVectorPointwiseMult, libceed), Cint, (CeedVector, CeedVector, CeedVector), w, x, y) end @@ -208,6 +250,10 @@ function CeedVectorReciprocal(vec) ccall((:CeedVectorReciprocal, libceed), Cint, (CeedVector,), vec) end +function CeedVectorViewRange(vec, start, stop, step, fp_fmt, stream) + ccall((:CeedVectorViewRange, libceed), Cint, (CeedVector, CeedSize, CeedSize, CeedInt, Ptr{Cchar}, Ptr{Libc.FILE}), vec, start, stop, step, fp_fmt, stream) +end + function CeedVectorView(vec, fp_fmt, stream) ccall((:CeedVectorView, libceed), Cint, (CeedVector, Ptr{Cchar}, Ptr{Libc.FILE}), vec, fp_fmt, stream) end @@ -228,11 +274,6 @@ function CeedRequestWait(req) ccall((:CeedRequestWait, libceed), Cint, (Ptr{CeedRequest},), req) end -@cenum CeedTransposeMode::UInt32 begin - CEED_NOTRANSPOSE = 0 - CEED_TRANSPOSE = 1 -end - function CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, rstr) ccall((:CeedElemRestrictionCreate, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedSize, CeedMemType, CeedCopyMode, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, rstr) end @@ -313,30 +354,6 @@ function CeedElemRestrictionDestroy(rstr) ccall((:CeedElemRestrictionDestroy, libceed), Cint, (Ptr{CeedElemRestriction},), rstr) end -@cenum CeedEvalMode::UInt32 begin - CEED_EVAL_NONE = 0 - CEED_EVAL_INTERP = 1 - CEED_EVAL_GRAD = 2 - CEED_EVAL_DIV = 4 - CEED_EVAL_CURL = 8 - CEED_EVAL_WEIGHT = 16 -end - -@cenum CeedQuadMode::UInt32 begin - CEED_GAUSS = 0 - CEED_GAUSS_LOBATTO = 1 -end - -@cenum CeedElemTopology::UInt32 begin - CEED_TOPOLOGY_LINE = 65536 - CEED_TOPOLOGY_TRIANGLE = 131073 - CEED_TOPOLOGY_QUAD = 131074 - CEED_TOPOLOGY_TET = 196611 - CEED_TOPOLOGY_PYRAMID = 196612 - CEED_TOPOLOGY_PRISM = 196613 - CEED_TOPOLOGY_HEX = 196614 -end - function CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, P, Q, quad_mode, basis) ccall((:CeedBasisCreateTensorH1Lagrange, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedQuadMode, Ptr{CeedBasis}), ceed, dim, num_comp, P, Q, quad_mode, basis) end @@ -353,6 +370,10 @@ function CeedBasisCreateHdiv(ceed, topo, num_comp, num_nodes, nqpts, interp, div ccall((:CeedBasisCreateHdiv, libceed), Cint, (Ceed, CeedElemTopology, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedBasis}), ceed, topo, num_comp, num_nodes, nqpts, interp, div, q_ref, q_weights, basis) end +function CeedBasisCreateHcurl(ceed, topo, num_comp, num_nodes, nqpts, interp, curl, q_ref, q_weights, basis) + ccall((:CeedBasisCreateHcurl, libceed), Cint, (Ceed, CeedElemTopology, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedBasis}), ceed, topo, num_comp, num_nodes, nqpts, interp, curl, q_ref, q_weights, basis) +end + function CeedBasisCreateProjection(basis_from, basis_to, basis_project) ccall((:CeedBasisCreateProjection, libceed), Cint, (CeedBasis, CeedBasis, Ptr{CeedBasis}), basis_from, basis_to, basis_project) end @@ -381,10 +402,6 @@ function CeedBasisGetTopology(basis, topo) ccall((:CeedBasisGetTopology, libceed), Cint, (CeedBasis, Ptr{CeedElemTopology}), basis, topo) end -function CeedBasisGetNumQuadratureComponents(basis, Q_comp) - ccall((:CeedBasisGetNumQuadratureComponents, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, Q_comp) -end - function CeedBasisGetNumComponents(basis, num_comp) ccall((:CeedBasisGetNumComponents, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, num_comp) end @@ -433,6 +450,10 @@ function CeedBasisGetDiv(basis, div) ccall((:CeedBasisGetDiv, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, div) end +function CeedBasisGetCurl(basis, curl) + ccall((:CeedBasisGetCurl, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, curl) +end + function CeedBasisDestroy(basis) ccall((:CeedBasisDestroy, libceed), Cint, (Ptr{CeedBasis},), basis) end @@ -516,11 +537,6 @@ function CeedQFunctionFieldGetEvalMode(qf_field, eval_mode) ccall((:CeedQFunctionFieldGetEvalMode, libceed), Cint, (CeedQFunctionField, Ptr{CeedEvalMode}), qf_field, eval_mode) end -@cenum CeedContextFieldType::UInt32 begin - CEED_CONTEXT_FIELD_DOUBLE = 1 - CEED_CONTEXT_FIELD_INT32 = 2 -end - # typedef int ( * CeedQFunctionContextDataDestroyUser ) ( void * data ) const CeedQFunctionContextDataDestroyUser = Ptr{Cvoid} @@ -716,32 +732,36 @@ function CeedOperatorGetFlopsEstimate(op, flops) ccall((:CeedOperatorGetFlopsEstimate, libceed), Cint, (CeedOperator, Ptr{CeedSize}), op, flops) end -function CeedOperatorContextGetFieldLabel(op, field_name, field_label) - ccall((:CeedOperatorContextGetFieldLabel, libceed), Cint, (CeedOperator, Ptr{Cchar}, Ptr{CeedContextFieldLabel}), op, field_name, field_label) +function CeedOperatorGetContext(op, ctx) + ccall((:CeedOperatorGetContext, libceed), Cint, (CeedOperator, Ptr{CeedQFunctionContext}), op, ctx) +end + +function CeedOperatorGetContextFieldLabel(op, field_name, field_label) + ccall((:CeedOperatorGetContextFieldLabel, libceed), Cint, (CeedOperator, Ptr{Cchar}, Ptr{CeedContextFieldLabel}), op, field_name, field_label) end -function CeedOperatorContextSetDouble(op, field_label, values) - ccall((:CeedOperatorContextSetDouble, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Cdouble}), op, field_label, values) +function CeedOperatorSetContextDouble(op, field_label, values) + ccall((:CeedOperatorSetContextDouble, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Cdouble}), op, field_label, values) end -function CeedOperatorContextGetDoubleRead(op, field_label, num_values, values) - ccall((:CeedOperatorContextGetDoubleRead, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Csize_t}, Ptr{Ptr{Cdouble}}), op, field_label, num_values, values) +function CeedOperatorGetContextDoubleRead(op, field_label, num_values, values) + ccall((:CeedOperatorGetContextDoubleRead, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Csize_t}, Ptr{Ptr{Cdouble}}), op, field_label, num_values, values) end -function CeedOperatorContextRestoreDoubleRead(op, field_label, values) - ccall((:CeedOperatorContextRestoreDoubleRead, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Ptr{Cdouble}}), op, field_label, values) +function CeedOperatorRestoreContextDoubleRead(op, field_label, values) + ccall((:CeedOperatorRestoreContextDoubleRead, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Ptr{Cdouble}}), op, field_label, values) end -function CeedOperatorContextSetInt32(op, field_label, values) - ccall((:CeedOperatorContextSetInt32, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Cint}), op, field_label, values) +function CeedOperatorSetContextInt32(op, field_label, values) + ccall((:CeedOperatorSetContextInt32, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Cint}), op, field_label, values) end -function CeedOperatorContextGetInt32Read(op, field_label, num_values, values) - ccall((:CeedOperatorContextGetInt32Read, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Csize_t}, Ptr{Ptr{Cint}}), op, field_label, num_values, values) +function CeedOperatorGetContextInt32Read(op, field_label, num_values, values) + ccall((:CeedOperatorGetContextInt32Read, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Csize_t}, Ptr{Ptr{Cint}}), op, field_label, num_values, values) end -function CeedOperatorContextRestoreInt32Read(op, field_label, values) - ccall((:CeedOperatorContextRestoreInt32Read, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Ptr{Cint}}), op, field_label, values) +function CeedOperatorRestoreContextInt32Read(op, field_label, values) + ccall((:CeedOperatorRestoreContextInt32Read, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Ptr{Cint}}), op, field_label, values) end function CeedOperatorApply(op, in, out, request) @@ -987,6 +1007,7 @@ end @cenum CeedFESpace::UInt32 begin CEED_FE_SPACE_H1 = 1 CEED_FE_SPACE_HDIV = 2 + CEED_FE_SPACE_HCURL = 3 end function CeedBasisGetCollocatedGrad(basis, colo_grad_1d) @@ -1009,10 +1030,18 @@ function CeedBasisReference(basis) ccall((:CeedBasisReference, libceed), Cint, (CeedBasis,), basis) end +function CeedBasisGetNumQuadratureComponents(basis, eval_mode, q_comp) + ccall((:CeedBasisGetNumQuadratureComponents, libceed), Cint, (CeedBasis, CeedEvalMode, Ptr{CeedInt}), basis, eval_mode, q_comp) +end + function CeedBasisGetFlopsEstimate(basis, t_mode, eval_mode, flops) ccall((:CeedBasisGetFlopsEstimate, libceed), Cint, (CeedBasis, CeedTransposeMode, CeedEvalMode, Ptr{CeedSize}), basis, t_mode, eval_mode, flops) end +function CeedBasisGetFESpace(basis, fe_space) + ccall((:CeedBasisGetFESpace, libceed), Cint, (CeedBasis, Ptr{CeedFESpace}), basis, fe_space) +end + function CeedBasisGetTopologyDimension(topo, dim) ccall((:CeedBasisGetTopologyDimension, libceed), Cint, (CeedElemTopology, Ptr{CeedInt}), topo, dim) end @@ -1033,6 +1062,10 @@ function CeedTensorContractApply(contract, A, B, C, J, t, t_mode, Add, u, v) ccall((:CeedTensorContractApply, libceed), Cint, (CeedTensorContract, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, CeedTransposeMode, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}), contract, A, B, C, J, t, t_mode, Add, u, v) end +function CeedTensorContractStridedApply(contract, A, B, C, D, J, t, t_mode, add, u, v) + ccall((:CeedTensorContractStridedApply, libceed), Cint, (CeedTensorContract, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, CeedTransposeMode, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}), contract, A, B, C, D, J, t, t_mode, add, u, v) +end + function CeedTensorContractGetCeed(contract, ceed) ccall((:CeedTensorContractGetCeed, libceed), Cint, (CeedTensorContract, Ptr{Ceed}), contract, ceed) end @@ -1331,7 +1364,9 @@ end # Skipping MacroDefinition: CEED_EXTERN extern CEED_VISIBILITY ( default ) -# Skipping MacroDefinition: CEED_QFUNCTION_HELPER CEED_QFUNCTION_ATTR static inline +# Skipping MacroDefinition: CEED_QFUNCTION_HELPER_ATTR CEED_QFUNCTION_ATTR __attribute__ ( ( always_inline ) ) + +# Skipping MacroDefinition: CEED_QFUNCTION_HELPER CEED_QFUNCTION_HELPER_ATTR static inline const CeedInt_FMT = "d" @@ -1341,7 +1376,7 @@ const CEED_VERSION_MINOR = 11 const CEED_VERSION_PATCH = 0 -const CEED_VERSION_RELEASE = true +const CEED_VERSION_RELEASE = false # Skipping MacroDefinition: CEED_INTERN extern CEED_VISIBILITY ( hidden ) diff --git a/julia/LibCEED.jl/test/buildmats.jl b/julia/LibCEED.jl/test/buildmats.jl new file mode 100644 index 0000000000..0e5d5e7758 --- /dev/null +++ b/julia/LibCEED.jl/test/buildmats.jl @@ -0,0 +1,78 @@ +function build_mats_hdiv(qref, qweight, ::Type{T}=Float64) where {T} + P, Q, dim = 4, 4, 2 + interp = Array{T}(undef, dim, Q, P) + div = Array{T}(undef, Q, P) + + qref[1, 1] = -1.0/sqrt(3.0) + qref[1, 2] = qref[1, 1] + qref[1, 3] = qref[1, 1] + qref[1, 4] = -qref[1, 1] + qref[2, 1] = -qref[1, 1] + qref[2, 2] = -qref[1, 1] + qref[2, 3] = qref[1, 1] + qref[2, 4] = qref[1, 1] + qweight[1] = 1.0 + qweight[2] = 1.0 + qweight[3] = 1.0 + qweight[4] = 1.0 + + # Loop over quadrature points + for i = 1:Q + x1 = qref[1, i] + x2 = qref[2, i] + # Interp + interp[1, i, 1] = 0.0 + interp[2, i, 1] = 1.0 - x2 + interp[1, i, 2] = x1 - 1.0 + interp[2, i, 2] = 0.0 + interp[1, i, 3] = -x1 + interp[2, i, 3] = 0.0 + interp[1, i, 4] = 0.0 + interp[2, i, 4] = x2 + # Div + div[i, 1] = -1.0 + div[i, 2] = 1.0 + div[i, 3] = -1.0 + div[i, 4] = 1.0 + end + + return interp, div +end + +function build_mats_hcurl(qref, qweight, ::Type{T}=Float64) where {T} + P, Q, dim = 3, 4, 2 + interp = Array{T}(undef, dim, Q, P) + curl = Array{T}(undef, 1, Q, P) + + qref[1, 1] = 0.2 + qref[1, 2] = 0.6 + qref[1, 3] = 1.0/3.0 + qref[1, 4] = 0.2 + qref[2, 1] = 0.2 + qref[2, 2] = 0.2 + qref[2, 3] = 1.0/3.0 + qref[2, 4] = 0.6 + qweight[1] = 25.0/96.0 + qweight[2] = 25.0/96.0 + qweight[3] = -27.0/96.0 + qweight[4] = 25.0/96.0 + + # Loop over quadrature points + for i = 1:Q + x1 = qref[1, i] + x2 = qref[2, i] + # Interp + interp[1, i, 1] = -x2 + interp[2, i, 1] = x1 + interp[1, i, 2] = x2 + interp[2, i, 2] = 1.0 - x1 + interp[1, i, 3] = 1.0 - x2 + interp[2, i, 3] = x1 + # Curl + curl[1, i, 1] = 2.0 + curl[1, i, 2] = -2.0 + curl[1, i, 3] = -2.0 + end + + return interp, curl +end diff --git a/julia/LibCEED.jl/test/rundevtests.jl b/julia/LibCEED.jl/test/rundevtests.jl index 75d0f78410..bc2a232e33 100644 --- a/julia/LibCEED.jl/test/rundevtests.jl +++ b/julia/LibCEED.jl/test/rundevtests.jl @@ -1,3 +1,48 @@ using Test, LibCEED, LinearAlgebra, StaticArrays -@testset "LibCEED Development Tests" begin end +include("buildmats.jl") + +@testset "LibCEED Development Tests" begin + @testset "Basis" begin + c = Ceed() + dim = 2 + ncomp = 1 + p1 = 4 + q1 = 4 + qref1 = Array{Float64}(undef, dim, q1) + qweight1 = Array{Float64}(undef, q1) + interp1, div1 = build_mats_hdiv(qref1, qweight1) + b1 = create_hdiv_basis(c, QUAD, ncomp, p1, q1, interp1, div1, qref1, qweight1) + + u1 = ones(Float64, p1) + v1 = apply(b1, u1) + + for i = 1:q1 + @test v1[i] ≈ -1.0 + @test v1[q1+i] ≈ 1.0 + end + + p2 = 3 + q2 = 4 + qref2 = Array{Float64}(undef, dim, q2) + qweight2 = Array{Float64}(undef, q2) + interp2, curl2 = build_mats_hcurl(qref2, qweight2) + b2 = create_hcurl_basis(c, TRIANGLE, ncomp, p2, q2, interp2, curl2, qref2, qweight2) + + u2 = [1.0, 2.0, 1.0] + v2 = apply(b2, u2) + + for i = 1:q2 + @test v2[i] ≈ 1.0 + end + + u2[1] = -1.0 + u2[2] = 1.0 + u2[3] = 2.0 + v2 = apply(b2, u2) + + for i = 1:q2 + @test v2[q2+i] ≈ 1.0 + end + end +end diff --git a/julia/LibCEED.jl/test/runtests.jl b/julia/LibCEED.jl/test/runtests.jl index e5e3f58637..8dbeed6b30 100644 --- a/julia/LibCEED.jl/test/runtests.jl +++ b/julia/LibCEED.jl/test/runtests.jl @@ -174,8 +174,18 @@ else @test getgrad1d(b2) == d1d @test checkoutput(showstr(b2), "b2.out") - b3 = create_h1_basis(c, LINE, 1, p, q, b1d, reshape(d1d, 1, q, p), q1d, w1d) - @test getqref(b3) == q1d + b3 = create_h1_basis( + c, + LINE, + 1, + p, + q, + b1d, + reshape(d1d, 1, q, p), + reshape(q1d, 1, q), + w1d, + ) + @test getqref(b3) == reshape(q1d, 1, q) @test getqweights(b3) == w1d @test checkoutput(showstr(b3), "b3.out") diff --git a/python/__init__.py b/python/__init__.py index 311b1af3b7..40cc969118 100644 --- a/python/__init__.py +++ b/python/__init__.py @@ -7,7 +7,7 @@ from .ceed import Ceed from .ceed_vector import Vector -from .ceed_basis import Basis, BasisTensorH1, BasisTensorH1Lagrange, BasisH1 +from .ceed_basis import Basis, BasisTensorH1, BasisTensorH1Lagrange, BasisH1, BasisHdiv, BasisHcurl from .ceed_elemrestriction import ElemRestriction, StridedElemRestriction, BlockedElemRestriction, BlockedStridedElemRestriction from .ceed_qfunction import QFunction, QFunctionByName, IdentityQFunction from .ceed_operator import Operator, CompositeOperator @@ -18,7 +18,7 @@ # ------------------------------------------------------------------------------ __all__ = ["Ceed", "Vector", - "Basis", "BasisTensorH1", "BasisTensorH1Lagrange", "BasisH1", + "Basis", "BasisTensorH1", "BasisTensorH1Lagrange", "BasisH1", "BasisHdiv", "BasisHcurl", "ElemRestriction", "StridedElemRestriction", "BlockedElemRestriction", "BlockedStridedelemRestriction", "QFunction", "QFunctionByName", "IdentityQFunction", "Operator", "CompositeOperator", diff --git a/python/ceed.py b/python/ceed.py index 6274c66f8a..9af5c34346 100644 --- a/python/ceed.py +++ b/python/ceed.py @@ -13,7 +13,7 @@ import tempfile from abc import ABC from .ceed_vector import Vector -from .ceed_basis import BasisTensorH1, BasisTensorH1Lagrange, BasisH1 +from .ceed_basis import BasisTensorH1, BasisTensorH1Lagrange, BasisH1, BasisHdiv, BasisHcurl from .ceed_elemrestriction import ElemRestriction, StridedElemRestriction, BlockedElemRestriction, BlockedStridedElemRestriction from .ceed_qfunction import QFunction, QFunctionByName, IdentityQFunction from .ceed_qfunctioncontext import QFunctionContext @@ -356,7 +356,7 @@ def BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight): *interp: Numpy array holding the row-major (nqpts * nnodes) matrix expressing the values of nodal basis functions at quadrature points - *grad: Numpy array holding the row-major (nqpts * dim * nnodes) + *grad: Numpy array holding the row-major (dim * nqpts * nnodes) matrix expressing the derivatives of nodal basis functions at quadrature points *qref: Numpy array of length (nqpts * dim) holding the locations of @@ -370,6 +370,59 @@ def BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight): return BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight) + def BasisHdiv(self, topo, ncomp, nnodes, nqpts, interp, div, qref, qweight): + """Ceed Hdiv Basis: finite element non tensor-product basis for H(div) + discretizations. + + Args: + topo: topology of the element, e.g. hypercube, simplex, etc + ncomp: number of field components (1 for scalar fields) + nnodes: total number of nodes + nqpts: total number of quadrature points + *interp: Numpy array holding the row-major (dim * nqpts * nnodes) + matrix expressing the values of basis functions at + quadrature points + *div: Numpy array holding the row-major (nqpts * nnodes) matrix + expressing the divergence of basis functions at + quadrature points + *qref: Numpy array of length (nqpts * dim) holding the locations of + quadrature points on the reference element [-1, 1] + *qweight: Numpy array of length nnodes holding the quadrature + weights on the reference element + + Returns: + basis: Ceed Basis""" + + return BasisHdiv(self, topo, ncomp, nnodes, nqpts, + interp, div, qref, qweight) + + def BasisHcurl(self, topo, ncomp, nnodes, nqpts, + interp, curl, qref, qweight): + """Ceed Hcurl Basis: finite element non tensor-product basis for H(curl) + discretizations. + + Args: + topo: topology of the element, e.g. hypercube, simplex, etc + ncomp: number of field components (1 for scalar fields) + nnodes: total number of nodes + nqpts: total number of quadrature points + *interp: Numpy array holding the row-major (dim * nqpts * nnodes) + matrix expressing the values of basis functions at + quadrature points + *curl: Numpy array holding the row-major (curlcomp * nqpts * nnodes), + curlcomp = 1 if dim < 3 else dim, matrix expressing the curl + of basis functions at quadrature points + *qref: Numpy array of length (nqpts * dim) holding the locations of + quadrature points on the reference element [-1, 1] + *qweight: Numpy array of length nnodes holding the quadrature + weights on the reference element + + Returns: + basis: Ceed Basis""" + + return BasisHcurl(self, topo, ncomp, nnodes, nqpts, + interp, curl, qref, qweight) + # CeedQFunction def QFunction(self, vlength, f, source): """Ceed QFunction: point-wise operation at quadrature points for diff --git a/python/ceed_basis.py b/python/ceed_basis.py index 6be8765b2d..02a9eb34e7 100644 --- a/python/ceed_basis.py +++ b/python/ceed_basis.py @@ -339,6 +339,88 @@ def __init__(self, ceed, topo, ncomp, nnodes, # ------------------------------------------------------------------------------ +class BasisHdiv(Basis): + """Ceed Hdiv Basis: finite element non tensor-product basis for H(div) discretizations.""" + + # Constructor + def __init__(self, ceed, topo, ncomp, nnodes, + nqpts, interp, grad, qref, qweight): + + # Setup arguments + self._pointer = ffi.new("CeedBasis *") + + self._ceed = ceed + + interp_pointer = ffi.new("CeedScalar *") + interp_pointer = ffi.cast( + "CeedScalar *", + interp.__array_interface__['data'][0]) + + div_pointer = ffi.new("CeedScalar *") + div_pointer = ffi.cast( + "CeedScalar *", + grad.__array_interface__['data'][0]) + + qref_pointer = ffi.new("CeedScalar *") + qref_pointer = ffi.cast( + "CeedScalar *", + qref.__array_interface__['data'][0]) + + qweight_pointer = ffi.new("CeedScalar *") + qweight_pointer = ffi.cast( + "CeedScalar *", + qweight.__array_interface__['data'][0]) + + # libCEED call + err_code = lib.CeedBasisCreateHdiv(self._ceed._pointer[0], topo, ncomp, + nnodes, nqpts, interp_pointer, + div_pointer, qref_pointer, + qweight_pointer, self._pointer) + +# ------------------------------------------------------------------------------ + + +class BasisHcurl(Basis): + """Ceed Hcurl Basis: finite element non tensor-product basis for H(curl) discretizations.""" + + # Constructor + def __init__(self, ceed, topo, ncomp, nnodes, + nqpts, interp, grad, qref, qweight): + + # Setup arguments + self._pointer = ffi.new("CeedBasis *") + + self._ceed = ceed + + interp_pointer = ffi.new("CeedScalar *") + interp_pointer = ffi.cast( + "CeedScalar *", + interp.__array_interface__['data'][0]) + + curl_pointer = ffi.new("CeedScalar *") + curl_pointer = ffi.cast( + "CeedScalar *", + grad.__array_interface__['data'][0]) + + qref_pointer = ffi.new("CeedScalar *") + qref_pointer = ffi.cast( + "CeedScalar *", + qref.__array_interface__['data'][0]) + + qweight_pointer = ffi.new("CeedScalar *") + qweight_pointer = ffi.cast( + "CeedScalar *", + qweight.__array_interface__['data'][0]) + + # libCEED call + err_code = lib.CeedBasisCreateHcurl(self._ceed._pointer[0], topo, ncomp, + nnodes, nqpts, interp_pointer, + curl_pointer, qref_pointer, + qweight_pointer, self._pointer) + +# ------------------------------------------------------------------------------ + + class TransposeBasis(): """Transpose Ceed Basis: transpose of finite element tensor-product basis objects.""" diff --git a/python/tests/buildmats.py b/python/tests/buildmats.py index e444c7d6e1..a361c1def9 100644 --- a/python/tests/buildmats.py +++ b/python/tests/buildmats.py @@ -47,3 +47,80 @@ def buildmats(qref, qweight, mat_dtype="float64"): grad[(i + Q) * P + 5] = 2. * (1. * (x2 - 1. / 2.) + x2 * 1.) return interp, grad + + +def buildmatshdiv(qref, qweight, mat_dtype="float64"): + P, Q, dim = 4, 4, 2 + interp = np.empty(dim * P * Q, dtype=mat_dtype) + div = np.empty(P * Q, dtype=mat_dtype) + + qref[0] = -1. / np.sqrt(3.) + qref[1] = qref[0] + qref[2] = qref[0] + qref[3] = -qref[0] + qref[4] = -qref[0] + qref[5] = -qref[0] + qref[6] = qref[0] + qref[7] = qref[0] + qweight[0] = 1. + qweight[1] = 1. + qweight[2] = 1. + qweight[3] = 1. + + # Loop over quadrature points + for i in range(Q): + x1 = qref[0 * Q + i] + x2 = qref[1 * Q + i] + # Interp + interp[(i + 0) * P + 0] = 0. + interp[(i + Q) * P + 0] = 1. - x2 + interp[(i + 0) * P + 1] = x1 - 1. + interp[(i + Q) * P + 1] = 0. + interp[(i + 0) * P + 2] = -x1 + interp[(i + Q) * P + 2] = 0. + interp[(i + 0) * P + 3] = 0. + interp[(i + Q) * P + 3] = x2 + # Div + div[i * P + 0] = -1. + div[i * P + 1] = 1. + div[i * P + 2] = -1. + div[i * P + 3] = 1. + + return interp, div + + +def buildmatshcurl(qref, qweight, mat_dtype="float64"): + P, Q, dim = 3, 4, 2 + interp = np.empty(dim * P * Q, dtype=mat_dtype) + curl = np.empty(P * Q, dtype=mat_dtype) + + qref[0] = 0.2 + qref[1] = 0.6 + qref[2] = 1. / 3. + qref[3] = 0.2 + qref[4] = 0.2 + qref[5] = 0.2 + qref[6] = 1. / 3. + qref[7] = 0.6 + qweight[0] = 25. / 96. + qweight[1] = 25. / 96. + qweight[2] = -27. / 96. + qweight[3] = 25. / 96. + + # Loop over quadrature points + for i in range(Q): + x1 = qref[0 * Q + i] + x2 = qref[1 * Q + i] + # Interp + interp[(i + 0) * P + 0] = -x2 + interp[(i + Q) * P + 0] = x1 + interp[(i + 0) * P + 1] = x2 + interp[(i + Q) * P + 1] = 1. - x1 + interp[(i + 0) * P + 2] = 1. - x2 + interp[(i + Q) * P + 2] = x1 + # Curl + curl[i * P + 0] = 2. + curl[i * P + 1] = -2. + curl[i * P + 2] = -2. + + return interp, curl diff --git a/python/tests/test-3-basis.py b/python/tests/test-3-basis.py index 2a08c23ed3..fe77428dde 100644 --- a/python/tests/test-3-basis.py +++ b/python/tests/test-3-basis.py @@ -303,3 +303,86 @@ def test_323(ceed_resource): assert math.fabs(out_array[1 * Q + i] - value) < 10. * TOL # ------------------------------------------------------------------------------- +# Test interpolation with a 2D Quad non-tensor H(div) basis +# ------------------------------------------------------------------------------- + + +def test_330(ceed_resource): + ceed = libceed.Ceed(ceed_resource) + + P, Q, dim = 4, 4, 2 + + in_array = np.ones(P, dtype=ceed.scalar_type()) + qref = np.empty(dim * Q, dtype=ceed.scalar_type()) + qweight = np.empty(Q, dtype=ceed.scalar_type()) + + interp, div = bm.buildmatshdiv(qref, qweight, libceed.scalar_types[ + libceed.lib.CEED_SCALAR_TYPE]) + + b = ceed.BasisHdiv(libceed.QUAD, 1, P, Q, interp, div, qref, qweight) + + # Interpolate function to quadrature points + in_vec = ceed.Vector(P) + in_vec.set_array(in_array, cmode=libceed.USE_POINTER) + out_vec = ceed.Vector(Q * dim) + out_vec.set_value(0) + + b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec) + + # Check values at quadrature points + with out_vec.array_read() as out_array: + for i in range(Q): + assert math.fabs(out_array[0 * Q + i] + 1.) < 10. * TOL + assert math.fabs(out_array[1 * Q + i] - 1.) < 10. * TOL + +# ------------------------------------------------------------------------------- +# Test interpolation with a 2D Simplex non-tensor H(curl) basis +# ------------------------------------------------------------------------------- + + +def test_340(ceed_resource): + ceed = libceed.Ceed(ceed_resource) + + P, Q, dim = 3, 4, 2 + + in_array = np.empty(P, dtype=ceed.scalar_type()) + qref = np.empty(dim * Q, dtype=ceed.scalar_type()) + qweight = np.empty(Q, dtype=ceed.scalar_type()) + + interp, curl = bm.buildmatshcurl(qref, qweight, libceed.scalar_types[ + libceed.lib.CEED_SCALAR_TYPE]) + + b = ceed.BasisHcurl(libceed.TRIANGLE, 1, P, Q, interp, curl, qref, qweight) + + # Interpolate function to quadrature points + in_array[0] = 1. + in_array[1] = 2. + in_array[2] = 1. + + in_vec = ceed.Vector(P) + in_vec.set_array(in_array, cmode=libceed.USE_POINTER) + out_vec = ceed.Vector(Q * dim) + out_vec.set_value(0) + + b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec) + + # Check values at quadrature points + with out_vec.array_read() as out_array: + for i in range(Q): + assert math.fabs(out_array[0 * Q + i] - 1.) < 10. * TOL + + # Interpolate function to quadrature points + in_array[0] = -1. + in_array[1] = 1. + in_array[2] = 2. + + out_vec.set_value(0) + + b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec) + + # Check values at quadrature points + with out_vec.array_read() as out_array: + for i in range(Q): + assert math.fabs(out_array[1 * Q + i] - 1.) < 10. * TOL + +# ------------------------------------------------------------------------------- diff --git a/rust/libceed/src/basis.rs b/rust/libceed/src/basis.rs index d3fe83edaa..ea3b2b24c0 100644 --- a/rust/libceed/src/basis.rs +++ b/rust/libceed/src/basis.rs @@ -238,6 +238,84 @@ impl<'a> Basis<'a> { }) } + pub fn create_Hdiv( + ceed: &crate::Ceed, + topo: crate::ElemTopology, + ncomp: usize, + nnodes: usize, + nqpts: usize, + interp: &[crate::Scalar], + div: &[crate::Scalar], + qref: &[crate::Scalar], + qweight: &[crate::Scalar], + ) -> crate::Result { + let mut ptr = std::ptr::null_mut(); + let (topo, ncomp, nnodes, nqpts) = ( + topo as bind_ceed::CeedElemTopology, + i32::try_from(ncomp).unwrap(), + i32::try_from(nnodes).unwrap(), + i32::try_from(nqpts).unwrap(), + ); + let ierr = unsafe { + bind_ceed::CeedBasisCreateHdiv( + ceed.ptr, + topo, + ncomp, + nnodes, + nqpts, + interp.as_ptr(), + div.as_ptr(), + qref.as_ptr(), + qweight.as_ptr(), + &mut ptr, + ) + }; + ceed.check_error(ierr)?; + Ok(Self { + ptr, + _lifeline: PhantomData, + }) + } + + pub fn create_Hcurl( + ceed: &crate::Ceed, + topo: crate::ElemTopology, + ncomp: usize, + nnodes: usize, + nqpts: usize, + interp: &[crate::Scalar], + curl: &[crate::Scalar], + qref: &[crate::Scalar], + qweight: &[crate::Scalar], + ) -> crate::Result { + let mut ptr = std::ptr::null_mut(); + let (topo, ncomp, nnodes, nqpts) = ( + topo as bind_ceed::CeedElemTopology, + i32::try_from(ncomp).unwrap(), + i32::try_from(nnodes).unwrap(), + i32::try_from(nqpts).unwrap(), + ); + let ierr = unsafe { + bind_ceed::CeedBasisCreateHcurl( + ceed.ptr, + topo, + ncomp, + nnodes, + nqpts, + interp.as_ptr(), + curl.as_ptr(), + qref.as_ptr(), + qweight.as_ptr(), + &mut ptr, + ) + }; + ceed.check_error(ierr)?; + Ok(Self { + ptr, + _lifeline: PhantomData, + }) + } + // Error handling #[doc(hidden)] fn check_error(&self, ierr: i32) -> crate::Result { diff --git a/rust/libceed/src/lib.rs b/rust/libceed/src/lib.rs index c11bc63da0..a884585169 100755 --- a/rust/libceed/src/lib.rs +++ b/rust/libceed/src/lib.rs @@ -586,7 +586,7 @@ impl Ceed { /// * `nqpts` - Total number of quadrature points /// * `interp` - Row-major `(nqpts * nnodes)` matrix expressing the values of /// nodal basis functions at quadrature points - /// * `grad` - Row-major `(nqpts * dim * nnodes)` matrix expressing + /// * `grad` - Row-major `(dim * nqpts * nnodes)` matrix expressing /// derivatives of nodal basis functions at quadrature points /// * `qref` - Array of length `nqpts` holding the locations of quadrature /// points on the reference element `[-1, 1]` @@ -707,6 +707,198 @@ impl Ceed { ) } + /// Returns an $H(div)$ Basis + /// + /// # arguments + /// + /// * `topo` - Topology of element, e.g. hypercube, simplex, ect + /// * `ncomp` - Number of field components (1 for scalar fields) + /// * `nnodes` - Total number of nodes + /// * `nqpts` - Total number of quadrature points + /// * `interp` - Row-major `(dim * nqpts * nnodes)` matrix expressing the + /// values of basis functions at quadrature points + /// * `div` - Row-major `(nqpts * nnodes)` matrix expressing the + /// divergence of basis functions at quadrature points + /// * `qref` - Array of length `nqpts` holding the locations of quadrature + /// points on the reference element `[-1, 1]` + /// * `qweight` - Array of length `nqpts` holding the quadrature weights on + /// the reference element + /// + /// ``` + /// # use libceed::prelude::*; + /// # fn main() -> libceed::Result<()> { + /// # let ceed = libceed::Ceed::default_init(); + /// let interp = [ + /// 0.00000000, + /// -1.57735027, + /// 0.57735027, + /// 0.00000000, + /// 0.00000000, + /// -1.57735027, + /// 0.57735027, + /// 0.00000000, + /// 0.00000000, + /// -1.57735027, + /// 0.57735027, + /// 0.00000000, + /// 0.00000000, + /// -0.42264973, + /// -0.57735027, + /// 0.00000000, + /// 0.42264973, + /// 0.00000000, + /// 0.00000000, + /// 0.57735027, + /// 0.42264973, + /// 0.00000000, + /// 0.00000000, + /// 0.57735027, + /// 1.57735027, + /// 0.00000000, + /// 0.00000000, + /// -0.57735027, + /// 1.57735027, + /// 0.00000000, + /// 0.00000000, + /// -0.57735027, + /// ]; + /// let div = [ + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// ]; + /// let qref = [ + /// 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, + /// 0.57735026, + /// ]; + /// let qweight = [1.00000000, 1.00000000, 1.00000000, 1.00000000]; + /// let b = ceed.basis_Hdiv(ElemTopology::Quad, 1, 4, 4, &interp, &div, &qref, &qweight)?; + /// # Ok(()) + /// # } + /// ``` + pub fn basis_Hdiv<'a>( + &self, + topo: ElemTopology, + ncomp: usize, + nnodes: usize, + nqpts: usize, + interp: &[crate::Scalar], + div: &[crate::Scalar], + qref: &[crate::Scalar], + qweight: &[crate::Scalar], + ) -> Result> { + Basis::create_Hdiv(self, topo, ncomp, nnodes, nqpts, interp, div, qref, qweight) + } + + /// Returns an $H(curl)$ Basis + /// + /// # arguments + /// + /// * `topo` - Topology of element, e.g. hypercube, simplex, ect + /// * `ncomp` - Number of field components (1 for scalar fields) + /// * `nnodes` - Total number of nodes + /// * `nqpts` - Total number of quadrature points + /// * `interp` - Row-major `(dim * nqpts * nnodes)` matrix expressing the + /// values of basis functions at quadrature points + /// * `curl` - Row-major `(curl_comp * nqpts * nnodes)`, `curl_comp = 1 if + /// dim < 3 else dim` matrix expressing the curl of basis + /// functions at quadrature points + /// * `qref` - Array of length `nqpts` holding the locations of quadrature + /// points on the reference element `[-1, 1]` + /// * `qweight` - Array of length `nqpts` holding the quadrature weights on + /// the reference element + /// + /// ``` + /// # use libceed::prelude::*; + /// # fn main() -> libceed::Result<()> { + /// # let ceed = libceed::Ceed::default_init(); + /// let interp = [ + /// -0.20000000, + /// 0.20000000, + /// 0.80000000, + /// -0.20000000, + /// 0.20000000, + /// 0.80000000, + /// -0.33333333, + /// 0.33333333, + /// 0.66666667, + /// -0.60000000, + /// 0.60000000, + /// 0.40000000, + /// 0.20000000, + /// 0.80000000, + /// 0.20000000, + /// 0.60000000, + /// 0.40000000, + /// 0.60000000, + /// 0.33333333, + /// 0.66666667, + /// 0.33333333, + /// 0.20000000, + /// 0.80000000, + /// 0.20000000, + /// ]; + /// let curl = [ + /// 2.00000000, + /// -2.00000000, + /// -2.00000000, + /// 2.00000000, + /// -2.00000000, + /// -2.00000000, + /// 2.00000000, + /// -2.00000000, + /// -2.00000000, + /// 2.00000000, + /// -2.00000000, + /// -2.00000000, + /// ]; + /// let qref = [ + /// 0.20000000, 0.60000000, 0.33333333, 0.20000000, 0.20000000, 0.20000000, 0.33333333, + /// 0.60000000, + /// ]; + /// let qweight = [0.26041667, 0.26041667, -0.28125000, 0.26041667]; + /// let b = ceed.basis_Hcurl( + /// ElemTopology::Triangle, + /// 1, + /// 3, + /// 4, + /// &interp, + /// &curl, + /// &qref, + /// &qweight, + /// )?; + /// # Ok(()) + /// # } + /// ``` + pub fn basis_Hcurl<'a>( + &self, + topo: ElemTopology, + ncomp: usize, + nnodes: usize, + nqpts: usize, + interp: &[crate::Scalar], + curl: &[crate::Scalar], + qref: &[crate::Scalar], + qweight: &[crate::Scalar], + ) -> Result> { + Basis::create_Hcurl( + self, topo, ncomp, nnodes, nqpts, interp, curl, qref, qweight, + ) + } + /// Returns a QFunction for evaluating interior (volumetric) terms /// of a weak form corresponding to the $L^2$ inner product /// diff --git a/tests/README.md b/tests/README.md index c000c4f60d..e461403046 100644 --- a/tests/README.md +++ b/tests/README.md @@ -4,22 +4,25 @@ This page provides a brief description of the tests for the libCEED library. The tests are organized by API object, and some tests are further organized, as required. -0. Ceed Tests +0. Ceed Tests 1. CeedVector Tests 1.0. CeedVector general tests - 1.1. CeedVector error tests -2. CeedElemRestriction Tests + 1.1. CeedVector error tests +2. CeedElemRestriction Tests 3. CeedBasis Tests 3.0. CeedBasis utility tests 3.1. CeedBasis tensor basis tests - 3.2. CeedBasis simplex basis tests + 3.2. CeedBasis simplex basis tests + 3.3. CeedBasis non-tensor H(div) basis tests + 3.4. CeedBasis non-tensor H(curl) basis tests 4. CeedQFunction Tests 4.0. CeedQFunction user code tests - 4.1. CeedQFunction gallery code tests + 4.1. CeedQFunction gallery code tests 5. CeedOperator Tests 5.0. CeedOperator with tensor bases tests 5.1. CeedOperator with simplex bases tests 5.2. CeedOperator with operator composition tests - 5.3. CeedOperator and CeedQFunction assembly tests + 5.3. CeedOperator diagonal and CeedQFunction assembly tests 5.4. CeedOperator element inverse tests - 5.5. CeedOperator multigrid level tests + 5.5. CeedOperator multigrid level tests + 5.6. CeedOperator full assembly tests diff --git a/tests/output/t340-basis.out b/tests/output/t340-basis.out new file mode 100644 index 0000000000..33dca069a9 --- /dev/null +++ b/tests/output/t340-basis.out @@ -0,0 +1,15 @@ +CeedBasis (H(curl) space on a triangle element): dim=2 P=8 Q=4 + qref: 0.33333333 0.20000000 0.20000000 0.60000000 0.33333333 0.20000000 0.60000000 0.20000000 + qweight: -0.26041667 0.26041667 0.28125000 0.26041667 + interp[0]: -0.22222222 0.44444444 0.22222222 -0.44444444 -0.22222222 -0.22222222 2.66666667 0.00000000 + interp[1]: 0.08000000 0.48000000 0.56000000 -0.48000000 1.04000000 -0.72000000 2.24000000 -0.64000000 + interp[2]: 0.24000000 -0.48000000 -0.24000000 0.48000000 -0.56000000 -0.56000000 2.88000000 0.00000000 + interp[3]: -0.56000000 0.48000000 -0.08000000 -0.48000000 -0.72000000 1.04000000 1.60000000 0.64000000 + interp[4]: -0.44444444 0.22222222 -0.22222222 -0.22222222 0.22222222 -0.44444444 0.00000000 2.66666667 + interp[5]: -0.48000000 -0.08000000 1.04000000 -0.72000000 0.56000000 -0.48000000 -0.64000000 2.24000000 + interp[6]: -0.48000000 0.56000000 -0.72000000 1.04000000 -0.08000000 -0.48000000 0.64000000 1.60000000 + interp[7]: 0.48000000 -0.24000000 -0.56000000 -0.56000000 -0.24000000 0.48000000 0.00000000 2.88000000 + curl[0]: 2.00000000 2.00000000 -2.00000000 -2.00000000 2.00000000 2.00000000 0.00000000 0.00000000 + curl[1]: -1.20000000 -1.20000000 -8.40000000 1.20000000 8.40000000 -1.20000000 -9.60000000 9.60000000 + curl[2]: -1.20000000 8.40000000 1.20000000 -8.40000000 -1.20000000 -1.20000000 9.60000000 9.60000000 + curl[3]: 8.40000000 -1.20000000 1.20000000 1.20000000 -1.20000000 8.40000000 0.00000000 -19.20000000 diff --git a/tests/t320-basis.c b/tests/t320-basis.c index 86534ef4d7..c028fcd0a5 100644 --- a/tests/t320-basis.c +++ b/tests/t320-basis.c @@ -1,6 +1,6 @@ /// @file -/// Test creation and destruction of a 2D Simplex non-tensor H1 basis -/// \test Test creation and destruction of a 2D Simplex non-tensor H1 basis +/// Test creation and destruction of a 2D Simplex non-tensor H^1 basis +/// \test Test creation and destruction of a 2D Simplex non-tensor H^1 basis #include "t320-basis.h" #include diff --git a/tests/t321-basis.c b/tests/t321-basis.c index a490f883b0..1274e4ef89 100644 --- a/tests/t321-basis.c +++ b/tests/t321-basis.c @@ -1,6 +1,6 @@ /// @file -/// Test interpolation with a 2D Simplex non-tensor H1 basis -/// \test Test interpolation with a 2D Simplex non-tensor H1 basis +/// Test interpolation with a 2D Simplex non-tensor H^1 basis +/// \test Test interpolation with a 2D Simplex non-tensor H^1 basis #include #include #include diff --git a/tests/t322-basis.c b/tests/t322-basis.c index fd739934c0..3b636ca9de 100644 --- a/tests/t322-basis.c +++ b/tests/t322-basis.c @@ -1,6 +1,6 @@ /// @file -/// Test integration with a 2D Simplex non-tensor H1 basis -/// \test Test integration with a 2D Simplex non-tensor H1 basis +/// Test integration with a 2D Simplex non-tensor H^1 basis +/// \test Test integration with a 2D Simplex non-tensor H^1 basis #include #include #include diff --git a/tests/t323-basis.c b/tests/t323-basis.c index e3408a378d..3a6adae84c 100644 --- a/tests/t323-basis.c +++ b/tests/t323-basis.c @@ -1,6 +1,6 @@ /// @file -/// Test grad with a 2D Simplex non-tensor H1 basis -/// \test Test grad with a 2D Simplex non-tensor H1 basis +/// Test grad with a 2D Simplex non-tensor H^1 basis +/// \test Test grad with a 2D Simplex non-tensor H^1 basis #include #include #include diff --git a/tests/t324-basis.c b/tests/t324-basis.c index 1d7c6f7862..6a4673f83b 100644 --- a/tests/t324-basis.c +++ b/tests/t324-basis.c @@ -1,6 +1,6 @@ /// @file -/// Test grad transpose with a 2D Simplex non-tensor H1 basis -/// \test Test grad transpose with a 2D Simplex non-tensor H1 basis +/// Test grad transpose with a 2D Simplex non-tensor H^1 basis +/// \test Test grad transpose with a 2D Simplex non-tensor H^1 basis #include #include #include diff --git a/tests/t325-basis.c b/tests/t325-basis.c index 844cb0e7d6..92479810ae 100644 --- a/tests/t325-basis.c +++ b/tests/t325-basis.c @@ -1,6 +1,6 @@ /// @file -/// Test grad transpose with a 2D Simplex non-tensor H1 basis -/// \test Test grad transpose with a 2D Simplex non-tensor H1 basis +/// Test grad transpose with a 2D Simplex non-tensor H^1 basis +/// \test Test grad transpose with a 2D Simplex non-tensor H^1 basis #include #include #include diff --git a/tests/t330-basis.c b/tests/t330-basis.c index 3133ecbe51..dfbf3373ff 100644 --- a/tests/t330-basis.c +++ b/tests/t330-basis.c @@ -1,6 +1,6 @@ /// @file -/// Test creation and destruction of a 2D Quad non-tensor Hdiv basis -/// \test Test creation and destruction of a 2D Quad non-tensor Hdiv basis +/// Test creation and destruction of a 2D Quad non-tensor H(div) basis +/// \test Test creation and destruction of a 2D Quad non-tensor H(div) basis #include "t330-basis.h" #include @@ -8,9 +8,7 @@ int main(int argc, char **argv) { Ceed ceed; - const CeedInt q = 3, dim = 2, num_qpts = q * q, elem_nodes = 4; - CeedInt num_comp = 1; - CeedInt p = dim * elem_nodes; // DoF per element, DoF are vector in H(div) + const CeedInt p = 8, q = 3, dim = 2, num_qpts = q * q; CeedBasis basis; CeedScalar q_ref[dim * num_qpts], q_weights[num_qpts]; CeedScalar interp[dim * p * num_qpts], div[p * num_qpts]; @@ -18,15 +16,10 @@ int main(int argc, char **argv) { CeedInit(argv[1], &ceed); // Test skipped if using single precision - if (CEED_SCALAR_TYPE == CEED_SCALAR_FP32) - // LCOV_EXCL_START - return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Test not implemented in single precision"); - // LCOV_EXCL_STOP + if (CEED_SCALAR_TYPE == CEED_SCALAR_FP32) return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Test not implemented in single precision"); BuildHdivQuadrilateral(q, q_ref, q_weights, interp, div, CEED_GAUSS); - CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, num_comp, p, num_qpts, interp, div, q_ref, q_weights, &basis); - // interp[0]--.interp[num_qpts-1] ==> basis in x-direction - // interp[num_qpts]--.interp[dim*num_qpts-1] ==> basis in y-direction + CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, 1, p, num_qpts, interp, div, q_ref, q_weights, &basis); CeedBasisView(basis, stdout); CeedBasisDestroy(&basis); diff --git a/tests/t330-basis.h b/tests/t330-basis.h index fd7069cb61..af973ff14f 100644 --- a/tests/t330-basis.h +++ b/tests/t330-basis.h @@ -7,8 +7,8 @@ #include -// Hdiv basis for quadrilateral linear BDM element in 2D -// Local numbering is as follow (each edge has 2 vector DoF) +// H(div) basis for quadrilateral linear BDM element in 2D +// Local numbering is as follows (each edge has 2 vector dof) // b4 b5 // 2---------3 // b7| |b3 @@ -40,6 +40,7 @@ int BuildNodalHdivQuadrilateral(CeedScalar *x, CeedScalar *Bx, CeedScalar *By) { By[7] = 0.125 + -0.125 * y_hat * y_hat; return 0; } + static void BuildHdivQuadrilateral(CeedInt q, CeedScalar *q_ref, CeedScalar *q_weights, CeedScalar *interp, CeedScalar *div, CeedQuadMode quad_mode) { // Get 1D quadrature on [-1,1] CeedScalar q_ref_1d[q], q_weight_1d[q]; @@ -56,10 +57,10 @@ static void BuildHdivQuadrilateral(CeedInt q, CeedScalar *q_ref, CeedScalar *q_w // Divergence operator; Divergence of nodal basis for ref element CeedScalar D[8] = {0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}; - // Loop over quadrature points CeedScalar Bx[8], By[8]; CeedScalar X[2]; + // Loop over quadrature points for (CeedInt i = 0; i < q; i++) { for (CeedInt j = 0; j < q; j++) { CeedInt k1 = q * i + j; @@ -76,4 +77,4 @@ static void BuildHdivQuadrilateral(CeedInt q, CeedScalar *q_ref, CeedScalar *q_w } } } -} \ No newline at end of file +} diff --git a/tests/t331-basis.c b/tests/t331-basis.c index ff423c30f2..75722b2b03 100644 --- a/tests/t331-basis.c +++ b/tests/t331-basis.c @@ -1,6 +1,6 @@ /// @file -/// Test GetInterp and BasisApply for a 2D Quad non-tensor H(div) basis -/// \test Test GetInterp and BasisApply for a 2D Quad non-tensor H(div) basis +/// Test interpolation with a 2D Quad non-tensor H(div) basis +/// \test Test interpolation with a 2D Quad non-tensor H(div) basis #include #include #include @@ -9,38 +9,32 @@ int main(int argc, char **argv) { Ceed ceed; - const CeedInt num_nodes = 4, q = 3, dim = 2, num_qpts = q * q; - CeedInt num_comp = 1; // one vector component - CeedInt p = dim * num_nodes; // DoF per element + CeedVector u, v; + const CeedInt p = 8, q = 3, dim = 2, num_qpts = q * q; CeedBasis basis; CeedScalar q_ref[dim * num_qpts], q_weights[num_qpts]; - CeedScalar div[p * num_qpts], interp[p * dim * num_qpts]; - CeedVector u, v; + CeedScalar interp[dim * p * num_qpts], div[p * num_qpts]; CeedInit(argv[1], &ceed); BuildHdivQuadrilateral(q, q_ref, q_weights, interp, div, CEED_GAUSS); - CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, num_comp, p, num_qpts, interp, div, q_ref, q_weights, &basis); + CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, 1, p, num_qpts, interp, div, q_ref, q_weights, &basis); - // Test GetInterp for H(div) + // Test interpolation for H(div) { - const CeedScalar *interp_2; - CeedBasisGetInterp(basis, &interp_2); - for (CeedInt i = 0; i < p * dim * num_qpts; i++) { - if (fabs(interp[i] - interp_2[i]) > 100. * CEED_EPSILON) { - // LCOV_EXCL_START - printf("%f != %f\n", interp[i], interp_2[i]); - // LCOV_EXCL_STOP - } + const CeedScalar *interp_in_basis; + + CeedBasisGetInterp(basis, &interp_in_basis); + for (CeedInt i = 0; i < dim * p * num_qpts; i++) { + if (fabs(interp[i] - interp_in_basis[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", interp[i], interp_in_basis[i]); } } CeedVectorCreate(ceed, p, &u); CeedVectorSetValue(u, 1.0); - CeedVectorCreate(ceed, num_qpts * dim, &v); - CeedVectorSetValue(v, 0.); + CeedVectorCreate(ceed, dim * num_qpts, &v); + CeedVectorSetValue(v, 0.0); - // BasisApply for H(div): CEED_EVAL_INTERP, NOTRANSPOSE case CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, v); { @@ -48,11 +42,7 @@ int main(int argc, char **argv) { CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); for (CeedInt i = 0; i < dim * num_qpts; i++) { - if (fabs(q_ref[i] - v_array[i]) > 100. * CEED_EPSILON) { - // LCOV_EXCL_START - printf("%f != %f\n", q_ref[i], v_array[i]); - // LCOV_EXCL_STOP - } + if (fabs(q_ref[i] - v_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", q_ref[i], v_array[i]); } CeedVectorRestoreArrayRead(v, &v_array); } @@ -60,7 +50,6 @@ int main(int argc, char **argv) { CeedVectorSetValue(v, 1.0); CeedVectorSetValue(u, 0.0); - // BasisApply for Hdiv: CEED_EVAL_INTERP, TRANSPOSE case CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, v, u); { diff --git a/tests/t332-basis.c b/tests/t332-basis.c index 49bbebb80d..3af8bdf49d 100644 --- a/tests/t332-basis.c +++ b/tests/t332-basis.c @@ -1,6 +1,6 @@ /// @file -/// Test GetDiv and BasisApply for a 2D Quad non-tensor H(div) basis -/// \test Test GetDiv and BasisApply for a 2D Quad non-tensor H(div) basis +/// Test div with a 2D Quad non-tensor H(div) basis +/// \test Test div with a 2D Quad non-tensor H(div) basis #include #include #include @@ -9,35 +9,32 @@ int main(int argc, char **argv) { Ceed ceed; - const CeedInt num_nodes = 4, q = 3, dim = 2, num_qpts = q * q; - CeedInt num_comp = 1; // one vector component - CeedInt p = dim * num_nodes; // DoF per element + CeedVector u, v; + const CeedInt p = 8, q = 3, dim = 2, num_qpts = q * q; CeedBasis basis; CeedScalar q_ref[dim * num_qpts], q_weights[num_qpts]; - CeedScalar div[p * num_qpts], interp[p * dim * num_qpts]; - CeedVector u, v; + CeedScalar interp[dim * p * num_qpts], div[p * num_qpts]; CeedInit(argv[1], &ceed); BuildHdivQuadrilateral(q, q_ref, q_weights, interp, div, CEED_GAUSS); - CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, num_comp, p, num_qpts, interp, div, q_ref, q_weights, &basis); + CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, 1, p, num_qpts, interp, div, q_ref, q_weights, &basis); - // Test GetDiv + // Test div for H(div) { - const CeedScalar *div_2; + const CeedScalar *div_in_basis; - CeedBasisGetDiv(basis, &div_2); + CeedBasisGetDiv(basis, &div_in_basis); for (CeedInt i = 0; i < p * num_qpts; i++) { - if (fabs(div[i] - div_2[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", div[i], div_2[i]); + if (fabs(div[i] - div_in_basis[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", div[i], div_in_basis[i]); } } CeedVectorCreate(ceed, p, &u); - CeedVectorSetValue(u, 1); + CeedVectorSetValue(u, 1.0); CeedVectorCreate(ceed, num_qpts, &v); - CeedVectorSetValue(v, 0); + CeedVectorSetValue(v, 0.0); - // BasisApply for H(div): CEED_EVAL_DIV, NOTRANSPOSE case CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_DIV, u, v); { @@ -53,7 +50,6 @@ int main(int argc, char **argv) { CeedVectorSetValue(v, 1.0); CeedVectorSetValue(u, 0.0); - // BasisApply for H(div): CEED_EVAL_DIV, TRANSPOSE case CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_DIV, v, u); { diff --git a/tests/t340-basis.c b/tests/t340-basis.c new file mode 100644 index 0000000000..35bd6eefeb --- /dev/null +++ b/tests/t340-basis.c @@ -0,0 +1,27 @@ +/// @file +/// Test creation and destruction of a 2D Simplex non-tensor H(curl) basis +/// \test Test creation and distruction of a 2D Simplex non-tensor H(curl) basis +#include "t340-basis.h" + +#include + +int main(int argc, char **argv) { + Ceed ceed; + const CeedInt p = 8, q = 4, dim = 2; + CeedBasis basis; + CeedScalar q_ref[dim * q], q_weight[q]; + CeedScalar interp[dim * p * q], curl[p * q]; + + CeedInit(argv[1], &ceed); + + // Test skipped if using single precision + if (CEED_SCALAR_TYPE == CEED_SCALAR_FP32) return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Test not implemented in single precision"); + + BuildHcurl2DSimplex(q_ref, q_weight, interp, curl); + CeedBasisCreateHcurl(ceed, CEED_TOPOLOGY_TRIANGLE, 1, p, q, interp, curl, q_ref, q_weight, &basis); + CeedBasisView(basis, stdout); + + CeedBasisDestroy(&basis); + CeedDestroy(&ceed); + return 0; +} diff --git a/tests/t340-basis.h b/tests/t340-basis.h new file mode 100644 index 0000000000..b98ca6100c --- /dev/null +++ b/tests/t340-basis.h @@ -0,0 +1,58 @@ +// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. +// All Rights Reserved. See the top-level LICENSE and NOTICE files for details. +// +// SPDX-License-Identifier: BSD-2-Clause +// +// This file is part of CEED: http://github.com/ceed + +#include + +// H(curl) basis for order 2 Nédélec (first kind) triangle element in 2D +// See: https://defelement.com/elements/examples/triangle-N1curl-2.html +static void BuildHcurl2DSimplex(CeedScalar *q_ref, CeedScalar *q_weight, CeedScalar *interp, CeedScalar *curl) { + CeedInt P = 8, Q = 4; + + q_ref[0] = 1. / 3.; + q_ref[4] = 1. / 3.; + q_ref[1] = 0.2; + q_ref[5] = 0.2; + q_ref[2] = 0.2; + q_ref[6] = 0.6; + q_ref[3] = 0.6; + q_ref[7] = 0.2; + q_weight[0] = -25. / 96.; + q_weight[1] = 25. / 96.; + q_weight[2] = 27. / 96.; + q_weight[3] = 25. / 96.; + + // Loop over quadrature points + for (int i = 0; i < Q; i++) { + CeedScalar x1 = q_ref[0 * Q + i], x2 = q_ref[1 * Q + i]; + // Interp + interp[(i + 0) * P + 0] = 2. * x2 * (1. - 4. * x1); + interp[(i + Q) * P + 0] = 4. * x1 * (2. * x1 - 1.); + interp[(i + 0) * P + 1] = 4. * x2 * (1. - 2. * x2); + interp[(i + Q) * P + 1] = 2. * x1 * (4. * x2 - 1.); + interp[(i + 0) * P + 2] = 2. * x2 * (-4. * x1 - 4. * x2 + 3.); + interp[(i + Q) * P + 2] = 8. * x1 * x1 + 8. * x1 * x2 - 12. * x1 - 6. * x2 + 4.; + interp[(i + 0) * P + 3] = 4. * x2 * (2. * x2 - 1.); + interp[(i + Q) * P + 3] = -8. * x1 * x2 + 2. * x1 + 6. * x2 - 2.; + interp[(i + 0) * P + 4] = 8. * x1 * x2 - 6. * x1 + 8. * x2 * x2 - 12. * x2 + 4.; + interp[(i + Q) * P + 4] = 2. * x1 * (-4. * x1 - 4. * x2 + 3.); + interp[(i + 0) * P + 5] = -8. * x1 * x2 + 6. * x1 + 2. * x2 - 2.; + interp[(i + Q) * P + 5] = 4. * x1 * (2. * x1 - 1.); + interp[(i + 0) * P + 6] = 8. * x2 * (-x1 - 2. * x2 + 2.); + interp[(i + Q) * P + 6] = 8. * x1 * (x1 + 2. * x2 - 1.); + interp[(i + 0) * P + 7] = 8. * x2 * (2. * x1 + x2 - 1.); + interp[(i + Q) * P + 7] = 8. * x1 * (-2. * x1 - x2 + 2.); + // Curl + curl[i * P + 0] = 24. * x1 - 6.; + curl[i * P + 1] = 24. * x2 - 6.; + curl[i * P + 2] = 24. * x1 + 24. * x2 - 18.; + curl[i * P + 3] = -24. * x2 + 6.; + curl[i * P + 4] = -24. * x1 - 24. * x2 + 18.; + curl[i * P + 5] = 24. * x1 - 6.; + curl[i * P + 6] = 24. * x1 + 48. * x2 - 24.; + curl[i * P + 7] = -48. * x1 - 24. * x1 + 24.; + } +} diff --git a/tests/t341-basis.c b/tests/t341-basis.c new file mode 100644 index 0000000000..dcf7ac0765 --- /dev/null +++ b/tests/t341-basis.c @@ -0,0 +1,83 @@ +/// @file +/// Test interpolation with a 2D Simplex non-tensor H(curl) basis +/// \test Test interpolation with a 2D Simplex non-tensor H(curl) basis +#include +#include + +#include "t340-basis.h" + +int main(int argc, char **argv) { + Ceed ceed; + CeedVector u, v; + const CeedInt p = 8, q = 4, dim = 2; + CeedBasis basis; + CeedScalar q_ref[dim * q], q_weight[q]; + CeedScalar interp[dim * p * q], curl[p * q]; + CeedScalar row_sum[dim * q], column_sum[p]; + + CeedInit(argv[1], &ceed); + + BuildHcurl2DSimplex(q_ref, q_weight, interp, curl); + CeedBasisCreateHcurl(ceed, CEED_TOPOLOGY_TRIANGLE, 1, p, q, interp, curl, q_ref, q_weight, &basis); + + // Test interpolation for H(curl) + { + const CeedScalar *interp_in_basis; + + CeedBasisGetInterp(basis, &interp_in_basis); + for (CeedInt i = 0; i < dim * p * q; i++) { + if (fabs(interp[i] - interp_in_basis[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", interp[i], interp_in_basis[i]); + } + } + + for (int i = 0; i < dim * q; i++) { + row_sum[i] = 0.0; + for (int j = 0; j < p; j++) { + row_sum[i] += interp[j + i * p]; + } + } + for (int i = 0; i < p; i++) { + column_sum[i] = 0.0; + for (int j = 0; j < dim * q; j++) { + column_sum[i] += interp[i + j * p]; + } + } + + CeedVectorCreate(ceed, p, &u); + CeedVectorSetValue(u, 1.0); + CeedVectorCreate(ceed, dim * q, &v); + CeedVectorSetValue(v, 0.0); + + CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, v); + + { + const CeedScalar *v_array; + + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + for (CeedInt i = 0; i < dim * q; i++) { + if (fabs(row_sum[i] - v_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", row_sum[i], v_array[i]); + } + CeedVectorRestoreArrayRead(v, &v_array); + } + + CeedVectorSetValue(v, 1.0); + CeedVectorSetValue(u, 0.0); + + CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, v, u); + + { + const CeedScalar *u_array; + + CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); + for (CeedInt i = 0; i < p; i++) { + if (fabs(column_sum[i] - u_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", column_sum[i], u_array[i]); + } + CeedVectorRestoreArrayRead(u, &u_array); + } + + CeedBasisDestroy(&basis); + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedDestroy(&ceed); + return 0; +} diff --git a/tests/t342-basis.c b/tests/t342-basis.c new file mode 100644 index 0000000000..8e30e5c344 --- /dev/null +++ b/tests/t342-basis.c @@ -0,0 +1,83 @@ +/// @file +/// Test curl with a 2D Simplex non-tensor H(curl) basis +/// \test Test curl with a 2D Simplex non-tensor H(curl) basis +#include +#include + +#include "t340-basis.h" + +int main(int argc, char **argv) { + Ceed ceed; + CeedVector u, v; + const CeedInt p = 8, q = 4, dim = 2; + CeedBasis basis; + CeedScalar q_ref[dim * q], q_weight[q]; + CeedScalar interp[dim * p * q], curl[p * q]; + CeedScalar row_sum[dim * q], column_sum[p]; + + CeedInit(argv[1], &ceed); + + BuildHcurl2DSimplex(q_ref, q_weight, interp, curl); + CeedBasisCreateHcurl(ceed, CEED_TOPOLOGY_TRIANGLE, 1, p, q, interp, curl, q_ref, q_weight, &basis); + + // Test curl for H(curl) + { + const CeedScalar *curl_in_basis; + + CeedBasisGetCurl(basis, &curl_in_basis); + for (CeedInt i = 0; i < p * q; i++) { + if (fabs(curl[i] - curl_in_basis[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", curl[i], curl_in_basis[i]); + } + } + + for (int i = 0; i < q; i++) { + row_sum[i] = 0; + for (int j = 0; j < p; j++) { + row_sum[i] += curl[j + i * p]; + } + } + for (int i = 0; i < p; i++) { + column_sum[i] = 0; + for (int j = 0; j < q; j++) { + column_sum[i] += curl[i + j * p]; + } + } + + CeedVectorCreate(ceed, p, &u); + CeedVectorSetValue(u, 1.0); + CeedVectorCreate(ceed, q, &v); + CeedVectorSetValue(v, 0.0); + + CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_CURL, u, v); + + { + const CeedScalar *v_array; + + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + for (CeedInt i = 0; i < q; i++) { + if (fabs(row_sum[i] - v_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", row_sum[i], v_array[i]); + } + CeedVectorRestoreArrayRead(v, &v_array); + } + + CeedVectorSetValue(v, 1.0); + CeedVectorSetValue(u, 0.0); + + CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_CURL, v, u); + + { + const CeedScalar *u_array; + + CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); + for (CeedInt i = 0; i < p; i++) { + if (fabs(column_sum[i] - u_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", column_sum[i], u_array[i]); + } + CeedVectorRestoreArrayRead(u, &u_array); + } + + CeedBasisDestroy(&basis); + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedDestroy(&ceed); + return 0; +}