1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7 #pragma once
8
9 #include <ceed.h>
10 #include <petscdm.h>
11
12 /**
13 @brief Copy the reference to a `Vec`.
14 Note: If `vec_copy` is non-null, it is assumed to be a valid pointer to a `Vec` and `VecDestroy()` will be called.
15
16 Collective across MPI processes.
17
18 @param[in] vec `Vec` to reference
19 @param[out] vec_copy Copy of reference
20
21 @return An error code: 0 - success, otherwise - failure
22 **/
VecReferenceCopy(Vec vec,Vec * vec_copy)23 static inline PetscErrorCode VecReferenceCopy(Vec vec, Vec *vec_copy) {
24 PetscFunctionBeginUser;
25 PetscCall(PetscObjectReference((PetscObject)vec));
26 PetscCall(VecDestroy(vec_copy));
27 *vec_copy = vec;
28 PetscFunctionReturn(PETSC_SUCCESS);
29 }
30
31 /**
32 @brief Copy the reference to a `DM`.
33 Note: If `dm_copy` is non-null, it is assumed to be a valid pointer to a `DM` and `DMDestroy()` will be called.
34
35 Collective across MPI processes.
36
37 @param[in] dm `DM` to reference
38 @param[out] dm_copy Copy of reference
39
40 @return An error code: 0 - success, otherwise - failure
41 **/
DMReferenceCopy(DM dm,DM * dm_copy)42 static inline PetscErrorCode DMReferenceCopy(DM dm, DM *dm_copy) {
43 PetscFunctionBeginUser;
44 PetscCall(PetscObjectReference((PetscObject)dm));
45 PetscCall(DMDestroy(dm_copy));
46 *dm_copy = dm;
47 PetscFunctionReturn(PETSC_SUCCESS);
48 }
49
50 /**
51 @brief Translate PetscMemType to CeedMemType
52
53 @param[in] mem_type PetscMemType
54
55 @return Equivalent CeedMemType
56 **/
57 /// @ingroup RatelInternal
MemTypePetscToCeed(PetscMemType mem_type)58 static inline CeedMemType MemTypePetscToCeed(PetscMemType mem_type) { return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST; }
59
60 /**
61 @brief Translate array of `PetscInt` to `CeedInt`.
62 If the types differ, `array_petsc` is freed with `PetscFree()` and `array_ceed` is allocated with `PetscMalloc1()`.
63 Caller is responsible for freeing `array_ceed` with `PetscFree()`.
64
65 Not collective across MPI processes.
66
67 @param[in] num_entries Number of array entries
68 @param[in,out] array_petsc Array of `PetscInt`
69 @param[out] array_ceed Array of `CeedInt`
70
71 @return An error code: 0 - success, otherwise - failure
72 **/
IntArrayCeedToPetsc(PetscInt num_entries,CeedInt ** array_ceed,PetscInt ** array_petsc)73 static inline PetscErrorCode IntArrayCeedToPetsc(PetscInt num_entries, CeedInt **array_ceed, PetscInt **array_petsc) {
74 const CeedInt int_c = 0;
75 const PetscInt int_p = 0;
76
77 PetscFunctionBeginUser;
78 if (sizeof(int_c) == sizeof(int_p)) {
79 *array_petsc = (PetscInt *)*array_ceed;
80 } else {
81 *array_petsc = malloc(num_entries * sizeof(PetscInt));
82 for (PetscInt i = 0; i < num_entries; i++) (*array_petsc)[i] = (*array_ceed)[i];
83 free(*array_ceed);
84 }
85 *array_ceed = NULL;
86 PetscFunctionReturn(PETSC_SUCCESS);
87 }
88
89 /**
90 @brief Translate array of `PetscInt` to `CeedInt`.
91 If the types differ, `array_petsc` is freed with `PetscFree()` and `array_ceed` is allocated with `PetscMalloc1()`.
92 Caller is responsible for freeing `array_ceed` with `PetscFree()`.
93
94 Not collective across MPI processes.
95
96 @param[in] num_entries Number of array entries
97 @param[in,out] array_petsc Array of `PetscInt`
98 @param[out] array_ceed Array of `CeedInt`
99
100 @return An error code: 0 - success, otherwise - failure
101 **/
IntArrayPetscToCeed(PetscInt num_entries,PetscInt ** array_petsc,CeedInt ** array_ceed)102 static inline PetscErrorCode IntArrayPetscToCeed(PetscInt num_entries, PetscInt **array_petsc, CeedInt **array_ceed) {
103 const CeedInt int_c = 0;
104 const PetscInt int_p = 0;
105
106 PetscFunctionBeginUser;
107 if (sizeof(int_c) == sizeof(int_p)) {
108 *array_ceed = (CeedInt *)*array_petsc;
109 } else {
110 PetscCall(PetscMalloc1(num_entries, array_ceed));
111 for (PetscInt i = 0; i < num_entries; i++) (*array_ceed)[i] = (*array_petsc)[i];
112 PetscCall(PetscFree(*array_petsc));
113 }
114 *array_petsc = NULL;
115 PetscFunctionReturn(PETSC_SUCCESS);
116 }
117
118 /**
119 @brief Transfer array from PETSc `Vec` to `CeedVector`.
120
121 Collective across MPI processes.
122
123 @param[in] X_petsc PETSc `Vec`
124 @param[out] mem_type PETSc `MemType`
125 @param[out] x_ceed `CeedVector`
126
127 @return An error code: 0 - success, otherwise - failure
128 **/
VecPetscToCeed(Vec X_petsc,PetscMemType * mem_type,CeedVector x_ceed)129 static inline PetscErrorCode VecPetscToCeed(Vec X_petsc, PetscMemType *mem_type, CeedVector x_ceed) {
130 PetscScalar *x;
131
132 PetscFunctionBeginUser;
133 PetscCall(VecGetArrayAndMemType(X_petsc, &x, mem_type));
134 PetscCallCeed(CeedVectorReturnCeed(x_ceed), CeedVectorSetArray(x_ceed, MemTypePetscToCeed(*mem_type), CEED_USE_POINTER, x));
135 PetscFunctionReturn(PETSC_SUCCESS);
136 }
137
138 /**
139 @brief Transfer array from `CeedVector` to PETSc `Vec`.
140
141 Collective across MPI processes.
142
143 @param[in] x_ceed `CeedVector`
144 @param[in] mem_type PETSc `MemType`
145 @param[out] X_petsc PETSc `Vec`
146
147 @return An error code: 0 - success, otherwise - failure
148 **/
VecCeedToPetsc(CeedVector x_ceed,PetscMemType mem_type,Vec X_petsc)149 static inline PetscErrorCode VecCeedToPetsc(CeedVector x_ceed, PetscMemType mem_type, Vec X_petsc) {
150 PetscScalar *x;
151
152 PetscFunctionBeginUser;
153 PetscCallCeed(CeedVectorReturnCeed(x_ceed), CeedVectorTakeArray(x_ceed, MemTypePetscToCeed(mem_type), &x));
154 PetscCall(VecRestoreArrayAndMemType(X_petsc, &x));
155 PetscFunctionReturn(PETSC_SUCCESS);
156 }
157
158 /**
159 @brief Transfer read only array from PETSc `Vec` to `CeedVector`.
160
161 Collective across MPI processes.
162
163 @param[in] X_petsc PETSc `Vec`
164 @param[out] mem_type PETSc `MemType`
165 @param[out] x_ceed `CeedVector`
166
167 @return An error code: 0 - success, otherwise - failure
168 **/
VecReadPetscToCeed(Vec X_petsc,PetscMemType * mem_type,CeedVector x_ceed)169 static inline PetscErrorCode VecReadPetscToCeed(Vec X_petsc, PetscMemType *mem_type, CeedVector x_ceed) {
170 PetscScalar *x;
171
172 PetscFunctionBeginUser;
173 PetscCall(VecGetArrayReadAndMemType(X_petsc, (const PetscScalar **)&x, mem_type));
174 PetscCallCeed(CeedVectorReturnCeed(x_ceed), CeedVectorSetArray(x_ceed, MemTypePetscToCeed(*mem_type), CEED_USE_POINTER, x));
175 PetscFunctionReturn(PETSC_SUCCESS);
176 }
177
178 /**
179 @brief Transfer read only array from `CeedVector` to PETSc `Vec`.
180
181 Collective across MPI processes.
182
183 @param[in] x_ceed `CeedVector`
184 @param[in] mem_type PETSc `MemType`
185 @param[out] X_petsc PETSc `Vec`
186
187 @return An error code: 0 - success, otherwise - failure
188 **/
VecReadCeedToPetsc(CeedVector x_ceed,PetscMemType mem_type,Vec X_petsc)189 static inline PetscErrorCode VecReadCeedToPetsc(CeedVector x_ceed, PetscMemType mem_type, Vec X_petsc) {
190 PetscScalar *x;
191
192 PetscFunctionBeginUser;
193 PetscCallCeed(CeedVectorReturnCeed(x_ceed), CeedVectorTakeArray(x_ceed, MemTypePetscToCeed(mem_type), &x));
194 PetscCall(VecRestoreArrayReadAndMemType(X_petsc, (const PetscScalar **)&x));
195 PetscFunctionReturn(PETSC_SUCCESS);
196 }
197
198 /**
199 @brief Copy PETSc `Vec` data into `CeedVector`
200
201 @param[in] X_petsc PETSc `Vec`
202 @param[out] x_ceed `CeedVector`
203
204 @return An error code: 0 - success, otherwise - failure
205 **/
VecCopyPetscToCeed(Vec X_petsc,CeedVector x_ceed)206 static inline PetscErrorCode VecCopyPetscToCeed(Vec X_petsc, CeedVector x_ceed) {
207 PetscScalar *x;
208 PetscMemType mem_type;
209 PetscInt X_size;
210 CeedSize x_size;
211 Ceed ceed;
212
213 PetscFunctionBeginUser;
214 PetscCall(CeedVectorGetCeed(x_ceed, &ceed));
215 PetscCall(VecGetLocalSize(X_petsc, &X_size));
216 PetscCallCeed(ceed, CeedVectorGetLength(x_ceed, &x_size));
217 PetscCheck(X_size == x_size, PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "X_petsc (%" PetscInt_FMT ") and x_ceed (%" CeedSize_FMT ") must be same size",
218 X_size, x_size);
219
220 PetscCall(VecGetArrayReadAndMemType(X_petsc, (const PetscScalar **)&x, &mem_type));
221 PetscCallCeed(ceed, CeedVectorSetArray(x_ceed, MemTypePetscToCeed(mem_type), CEED_COPY_VALUES, x));
222 PetscCall(VecRestoreArrayReadAndMemType(X_petsc, (const PetscScalar **)&x));
223 PetscCheck(CeedDestroy(&ceed) == CEED_ERROR_SUCCESS, PetscObjectComm((PetscObject)X_petsc), PETSC_ERR_LIB, "Destroying Ceed object failed");
224 PetscFunctionReturn(PETSC_SUCCESS);
225 }
226
227 /**
228 @brief Return the quadrature size from the eval mode, dimension, and number of components
229
230 @param[in] eval_mode The basis evaluation mode
231 @param[in] dim The basis dimension
232 @param[in] num_components The basis number of components
233
234 @return The maximum of the two integers
235
236 **/
237 /// @ingroup RatelInternal
GetCeedQuadratureSize(CeedEvalMode eval_mode,CeedInt dim,CeedInt num_components)238 static inline CeedInt GetCeedQuadratureSize(CeedEvalMode eval_mode, CeedInt dim, CeedInt num_components) {
239 switch (eval_mode) {
240 case CEED_EVAL_INTERP:
241 return num_components;
242 case CEED_EVAL_GRAD:
243 return dim * num_components;
244 default:
245 return -1;
246 }
247 }
248
249 /**
250 @brief Convert from DMPolytopeType to CeedElemTopology
251
252 @param[in] cell_type DMPolytopeType for the cell
253
254 @return CeedElemTopology, or 0 if no equivelent CeedElemTopology was found
255 **/
PolytopeTypePetscToCeed(DMPolytopeType cell_type)256 static inline CeedElemTopology PolytopeTypePetscToCeed(DMPolytopeType cell_type) {
257 switch (cell_type) {
258 case DM_POLYTOPE_TRIANGLE:
259 return CEED_TOPOLOGY_TRIANGLE;
260 case DM_POLYTOPE_QUADRILATERAL:
261 return CEED_TOPOLOGY_QUAD;
262 case DM_POLYTOPE_TETRAHEDRON:
263 return CEED_TOPOLOGY_TET;
264 case DM_POLYTOPE_HEXAHEDRON:
265 return CEED_TOPOLOGY_HEX;
266 default:
267 return 0;
268 }
269 }
270