1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // Authors: Peter Munch, Martin Kronbichler
15 //
16 // ---------------------------------------------------------------------
17
18 #include <ceed/types.h>
19
20
21
22 /**
23 * Context passed to libCEED Q-function.
24 */
25 struct BuildContext
26 {
27 CeedInt dim, space_dim;
28 };
29
30
31
32 /**
33 * libCEED Q-function for building quadrature data for a mass operator
34 */
f_build_mass(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)35 CEED_QFUNCTION(f_build_mass)
36 (void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out)
37 {
38 BuildContext *bc = (BuildContext *)ctx;
39 const CeedScalar *J = in[0], *w = in[1];
40 CeedScalar *qdata = out[0];
41
42 switch (bc->dim + 10 * bc->space_dim)
43 {
44 case 11:
45 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
46 {
47 qdata[i] = J[i] * w[i];
48 }
49 break;
50 case 22:
51 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
52 {
53 qdata[i] = (J[i + Q * 0] * J[i + Q * 3] - J[i + Q * 1] * J[i + Q * 2]) * w[i];
54 }
55 break;
56 case 33:
57 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
58 {
59 qdata[i] = (J[i + Q * 0] * (J[i + Q * 4] * J[i + Q * 8] - J[i + Q * 5] * J[i + Q * 7]) -
60 J[i + Q * 1] * (J[i + Q * 3] * J[i + Q * 8] - J[i + Q * 5] * J[i + Q * 6]) +
61 J[i + Q * 2] * (J[i + Q * 3] * J[i + Q * 7] - J[i + Q * 4] * J[i + Q * 6])) *
62 w[i];
63 }
64 break;
65 }
66 return 0;
67 }
68
69
70
71 /**
72 * libCEED Q-function for applying a mass operator
73 */
f_apply_mass(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)74 CEED_QFUNCTION(f_apply_mass)
75 (void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out)
76 {
77 (void)ctx;
78
79 const CeedScalar *u = in[0];
80 const CeedScalar *JxW = in[1];
81 CeedScalar *v = out[0];
82
83 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
84 {
85 v[i] = JxW[i] * u[i];
86 }
87 return 0;
88 }
89
90
91
92 /**
93 * libCEED Q-function for applying a vector mass operator
94 */
f_apply_mass_vec(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)95 CEED_QFUNCTION(f_apply_mass_vec)
96 (void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out)
97 {
98 BuildContext *bc = (BuildContext *)ctx;
99
100 const CeedScalar *u = in[0];
101 const CeedScalar *JxW = in[1];
102 CeedScalar *v = out[0];
103
104 switch (bc->dim + 10 * bc->space_dim)
105 {
106 case 11:
107 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
108 {
109 v[i + Q * 0] = JxW[i] * u[i + Q * 0];
110 }
111 break;
112 case 22:
113 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
114 {
115 v[i + Q * 0] = JxW[i] * u[i + Q * 0];
116 v[i + Q * 1] = JxW[i] * u[i + Q * 1];
117 }
118 break;
119 case 33:
120 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
121 {
122 v[i + Q * 0] = JxW[i] * u[i + Q * 0];
123 v[i + Q * 1] = JxW[i] * u[i + Q * 1];
124 v[i + Q * 2] = JxW[i] * u[i + Q * 2];
125 }
126 break;
127 }
128 return 0;
129 }
130
131
132
133 /**
134 * libCEED Q-function for building quadrature data for a Poisson operator
135 */
f_build_poisson(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)136 CEED_QFUNCTION(f_build_poisson)
137 (void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out)
138 {
139 BuildContext *bc = (BuildContext *)ctx;
140 const CeedScalar *J = in[0], *w = in[1];
141 CeedScalar *qdata = out[0];
142
143 switch (bc->dim + 10 * bc->space_dim)
144 {
145 case 11:
146 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++)
147 {
148 qdata[i] = w[i] / J[i];
149 }
150 break;
151 case 22:
152 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++)
153 {
154 const CeedScalar J11 = J[i + Q * 0];
155 const CeedScalar J21 = J[i + Q * 1];
156 const CeedScalar J12 = J[i + Q * 2];
157 const CeedScalar J22 = J[i + Q * 3];
158 const CeedScalar qw = w[i] / (J11 * J22 - J21 * J12);
159 qdata[i + Q * 0] = qw * (J12 * J12 + J22 * J22);
160 qdata[i + Q * 1] = qw * (J11 * J11 + J21 * J21);
161 qdata[i + Q * 2] = -qw * (J11 * J12 + J21 * J22);
162 }
163 break;
164 case 33:
165 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++)
166 {
167 const CeedScalar J11 = J[i + Q * 0];
168 const CeedScalar J21 = J[i + Q * 1];
169 const CeedScalar J31 = J[i + Q * 2];
170 const CeedScalar J12 = J[i + Q * 3];
171 const CeedScalar J22 = J[i + Q * 4];
172 const CeedScalar J32 = J[i + Q * 5];
173 const CeedScalar J13 = J[i + Q * 6];
174 const CeedScalar J23 = J[i + Q * 7];
175 const CeedScalar J33 = J[i + Q * 8];
176 const CeedScalar A11 = J22 * J33 - J23 * J32;
177 const CeedScalar A12 = J13 * J32 - J12 * J33;
178 const CeedScalar A13 = J12 * J23 - J13 * J22;
179 const CeedScalar A21 = J23 * J31 - J21 * J33;
180 const CeedScalar A22 = J11 * J33 - J13 * J31;
181 const CeedScalar A23 = J13 * J21 - J11 * J23;
182 const CeedScalar A31 = J21 * J32 - J22 * J31;
183 const CeedScalar A32 = J12 * J31 - J11 * J32;
184 const CeedScalar A33 = J11 * J22 - J12 * J21;
185 const CeedScalar qw = w[i] / (J11 * A11 + J21 * A12 + J31 * A13);
186 qdata[i + Q * 0] = qw * (A11 * A11 + A12 * A12 + A13 * A13);
187 qdata[i + Q * 1] = qw * (A21 * A21 + A22 * A22 + A23 * A23);
188 qdata[i + Q * 2] = qw * (A31 * A31 + A32 * A32 + A33 * A33);
189 qdata[i + Q * 3] = qw * (A21 * A31 + A22 * A32 + A23 * A33);
190 qdata[i + Q * 4] = qw * (A11 * A31 + A12 * A32 + A13 * A33);
191 qdata[i + Q * 5] = qw * (A11 * A21 + A12 * A22 + A13 * A23);
192 }
193 break;
194 }
195 return 0;
196 }
197
198
199
200 /**
201 * libCEED Q-function for applying a Poisson operator
202 */
f_apply_poisson(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)203 CEED_QFUNCTION(f_apply_poisson)
204 (void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out)
205 {
206 BuildContext *bc = (BuildContext *)ctx;
207 const CeedScalar *ug = in[0], *qdata = in[1];
208 CeedScalar *vg = out[0];
209
210 switch (bc->dim)
211 {
212 case 1:
213 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++)
214 {
215 vg[i] = ug[i] * qdata[i];
216 }
217 break;
218 case 2:
219 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++)
220 {
221 const CeedScalar ug0 = ug[i + Q * 0];
222 const CeedScalar ug1 = ug[i + Q * 1];
223 vg[i + Q * 0] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 2] * ug1;
224 vg[i + Q * 1] = qdata[i + Q * 2] * ug0 + qdata[i + Q * 1] * ug1;
225 }
226 break;
227 case 3:
228 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++)
229 {
230 const CeedScalar ug0 = ug[i + Q * 0];
231 const CeedScalar ug1 = ug[i + Q * 1];
232 const CeedScalar ug2 = ug[i + Q * 2];
233 vg[i + Q * 0] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 5] * ug1 + qdata[i + Q * 4] * ug2;
234 vg[i + Q * 1] = qdata[i + Q * 5] * ug0 + qdata[i + Q * 1] * ug1 + qdata[i + Q * 3] * ug2;
235 vg[i + Q * 2] = qdata[i + Q * 4] * ug0 + qdata[i + Q * 3] * ug1 + qdata[i + Q * 2] * ug2;
236 }
237 break;
238 }
239 return 0;
240 }
241
242
243
244 /**
245 *libCEED Q-function for applying a vector Poisson operator
246 */
f_apply_poisson_vec(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)247 CEED_QFUNCTION(f_apply_poisson_vec)
248 (void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out)
249 {
250 BuildContext *bc = (BuildContext *)ctx;
251 const CeedScalar *ug = in[0], *qdata = in[1];
252 CeedScalar *vg = out[0];
253
254 switch (bc->dim)
255 {
256 case 1:
257 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++)
258 {
259 vg[i] = ug[i] * qdata[i];
260 }
261 break;
262 case 2:
263 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++)
264 {
265 {
266 const CeedScalar ug0 = ug[i + Q * 0 + Q * 2 * 0];
267 const CeedScalar ug1 = ug[i + Q * 1 + Q * 2 * 0];
268 vg[i + Q * 0 + Q * 2 * 0] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 2] * ug1;
269 vg[i + Q * 1 + Q * 2 * 0] = qdata[i + Q * 2] * ug0 + qdata[i + Q * 1] * ug1;
270 }
271 {
272 const CeedScalar ug0 = ug[i + Q * 0 + Q * 2 * 1];
273 const CeedScalar ug1 = ug[i + Q * 1 + Q * 2 * 1];
274 vg[i + Q * 0 + Q * 2 * 1] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 2] * ug1;
275 vg[i + Q * 1 + Q * 2 * 1] = qdata[i + Q * 2] * ug0 + qdata[i + Q * 1] * ug1;
276 }
277 }
278 break;
279 case 3:
280 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++)
281 {
282 {
283 const CeedScalar ug0 = ug[i + Q * 0 + Q * 3 * 0];
284 const CeedScalar ug1 = ug[i + Q * 1 + Q * 3 * 0];
285 const CeedScalar ug2 = ug[i + Q * 2 + Q * 3 * 0];
286 vg[i + Q * 0 + Q * 3 * 0] =
287 qdata[i + Q * 0] * ug0 + qdata[i + Q * 5] * ug1 + qdata[i + Q * 4] * ug2;
288 vg[i + Q * 1 + Q * 3 * 0] =
289 qdata[i + Q * 5] * ug0 + qdata[i + Q * 1] * ug1 + qdata[i + Q * 3] * ug2;
290 vg[i + Q * 2 + Q * 3 * 0] =
291 qdata[i + Q * 4] * ug0 + qdata[i + Q * 3] * ug1 + qdata[i + Q * 2] * ug2;
292 }
293 {
294 const CeedScalar ug0 = ug[i + Q * 0 + Q * 3 * 1];
295 const CeedScalar ug1 = ug[i + Q * 1 + Q * 3 * 1];
296 const CeedScalar ug2 = ug[i + Q * 2 + Q * 3 * 1];
297 vg[i + Q * 0 + Q * 3 * 1] =
298 qdata[i + Q * 0] * ug0 + qdata[i + Q * 5] * ug1 + qdata[i + Q * 4] * ug2;
299 vg[i + Q * 1 + Q * 3 * 1] =
300 qdata[i + Q * 5] * ug0 + qdata[i + Q * 1] * ug1 + qdata[i + Q * 3] * ug2;
301 vg[i + Q * 2 + Q * 3 * 1] =
302 qdata[i + Q * 4] * ug0 + qdata[i + Q * 3] * ug1 + qdata[i + Q * 2] * ug2;
303 }
304 {
305 const CeedScalar ug0 = ug[i + Q * 0 + Q * 3 * 2];
306 const CeedScalar ug1 = ug[i + Q * 1 + Q * 3 * 2];
307 const CeedScalar ug2 = ug[i + Q * 2 + Q * 3 * 2];
308 vg[i + Q * 0 + Q * 3 * 2] =
309 qdata[i + Q * 0] * ug0 + qdata[i + Q * 5] * ug1 + qdata[i + Q * 4] * ug2;
310 vg[i + Q * 1 + Q * 3 * 2] =
311 qdata[i + Q * 5] * ug0 + qdata[i + Q * 1] * ug1 + qdata[i + Q * 3] * ug2;
312 vg[i + Q * 2 + Q * 3 * 2] =
313 qdata[i + Q * 4] * ug0 + qdata[i + Q * 3] * ug1 + qdata[i + Q * 2] * ug2;
314 }
315 }
316 break;
317 }
318 return 0;
319 }
320