1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #include "solvers-atpmg.hpp"
13 #include "util.hpp"
14
15 #ifdef MFEM_USE_CEED
16 #include <ceed/backend.h>
17
18 #include <math.h>
19 // todo: should probably use Ceed memory wrappers instead of calloc/free?
20 #include <stdlib.h>
21
22 namespace mfem
23 {
24
25 namespace ceed
26 {
27
28 // In one dimension, return corresponding coarse edof index for
29 // given fine index, with -1 meaning the edof disappears on the
30 // coarse grid
coarse_1d_edof(int i,int P1d,int coarse_P1d)31 int coarse_1d_edof(int i, int P1d, int coarse_P1d)
32 {
33 int coarse_i = (i < coarse_P1d - 1) ? i : -1;
34 if (i == P1d - 1)
35 {
36 coarse_i = coarse_P1d - 1;
37 }
38 return coarse_i;
39 }
40
reverse_coarse_1d_edof(int i,int P1d,int coarse_P1d)41 int reverse_coarse_1d_edof(int i, int P1d, int coarse_P1d)
42 {
43 int coarse_i;
44 if (i > P1d - coarse_P1d)
45 {
46 coarse_i = i - (P1d - coarse_P1d);
47 }
48 else
49 {
50 coarse_i = -1;
51 }
52 if (i == 0)
53 {
54 coarse_i = 0;
55 }
56 return coarse_i;
57 }
58
min4(int a,int b,int c,int d)59 int min4(int a, int b, int c, int d)
60 {
61 if (a <= b && a <= c && a <= d)
62 {
63 return a;
64 }
65 else if (b <= a && b <= c && b <= d)
66 {
67 return b;
68 }
69 else if (c <= a && c <= b && c <= d)
70 {
71 return c;
72 }
73 else
74 {
75 return d;
76 }
77 }
78
CeedATPMGElemRestriction(int order,int order_reduction,CeedElemRestriction er_in,CeedElemRestriction * er_out,CeedInt * & dof_map)79 int CeedATPMGElemRestriction(int order,
80 int order_reduction,
81 CeedElemRestriction er_in,
82 CeedElemRestriction* er_out,
83 CeedInt *&dof_map)
84 {
85 int ierr;
86 Ceed ceed;
87 ierr = CeedElemRestrictionGetCeed(er_in, &ceed); CeedChk(ierr);
88
89 CeedInt numelem, numnodes, numcomp, elemsize;
90 ierr = CeedElemRestrictionGetNumElements(er_in, &numelem); CeedChk(ierr);
91 ierr = CeedElemRestrictionGetLVectorSize(er_in, &numnodes); CeedChk(ierr);
92 ierr = CeedElemRestrictionGetElementSize(er_in, &elemsize); CeedChk(ierr);
93 ierr = CeedElemRestrictionGetNumComponents(er_in, &numcomp); CeedChk(ierr);
94 if (numcomp != 1)
95 {
96 // todo: multi-component will require more thought
97 return CeedError(ceed, 1, "Algebraic element restriction not "
98 "implemented for multiple components.");
99 }
100
101 int P1d = order + 1;
102 int coarse_P1d = P1d - order_reduction;
103 int dim = (log((double) elemsize) / log((double) P1d)) + 1.e-3;
104
105 CeedVector in_lvec, in_evec;
106 ierr = CeedElemRestrictionCreateVector(er_in, &in_lvec, &in_evec);
107 CeedChk(ierr);
108
109 // Create the elem_dof array from the given high-order ElemRestriction
110 // by using it to map the L-vector indices to an E-vector
111 CeedScalar * lvec_data;
112 ierr = CeedVectorGetArray(in_lvec, CEED_MEM_HOST, &lvec_data); CeedChk(ierr);
113 for (int i = 0; i < numnodes; ++i)
114 {
115 lvec_data[i] = (CeedScalar) i;
116 }
117 ierr = CeedVectorRestoreArray(in_lvec, &lvec_data); CeedChk(ierr);
118 CeedInt in_layout[3];
119 ierr = CeedElemRestrictionGetELayout(er_in, &in_layout); CeedChk(ierr);
120 if (in_layout[0] == 0 && in_layout[1] == 0 && in_layout[2] == 0)
121 {
122 return CeedError(ceed, 1, "Cannot interpret e-vector ordering of given"
123 "CeedElemRestriction!");
124 }
125 ierr = CeedElemRestrictionApply(er_in, CEED_NOTRANSPOSE, in_lvec, in_evec,
126 CEED_REQUEST_IMMEDIATE); CeedChk(ierr);
127 ierr = CeedVectorDestroy(&in_lvec); CeedChk(ierr);
128 const CeedScalar * in_elem_dof;
129 ierr = CeedVectorGetArrayRead(in_evec, CEED_MEM_HOST, &in_elem_dof);
130 CeedChk(ierr);
131
132 // Create a map (dof_map) that maps high-order ldof indices to
133 // low-order ldof indices, with -1 indicating no correspondence
134 // (NOTE: it is the caller's responsibility to free dof_map)
135 dof_map = new CeedInt[numnodes];
136 for (int i = 0; i < numnodes; ++i)
137 {
138 dof_map[i] = -1;
139 }
140 CeedInt coarse_elemsize = pow(coarse_P1d, dim);
141 CeedInt * out_elem_dof = new CeedInt[coarse_elemsize * numelem];
142 const double rounding_guard = 1.e-10;
143 int running_out_ldof_count = 0;
144 if (dim == 2)
145 {
146 for (int e = 0; e < numelem; ++e)
147 {
148 // Loop over edofs in element
149 for (int i = 0; i < P1d; ++i)
150 {
151 for (int j = 0; j < P1d; ++j)
152 {
153 // Determine topology; is this edof on the outside of the element
154 // in the i or j direction?
155 int in_edof = i*P1d + j;
156 const int edof_index = in_edof*in_layout[0] + e*in_layout[2];
157 int in_ldof = in_elem_dof[edof_index] + rounding_guard;
158 bool i_edge = (i == 0 || i == P1d - 1);
159 bool j_edge = (j == 0 || j == P1d - 1);
160
161 // Determine corresponding coarse 1D edof indices
162 // We do this systematically, orienting edges and faces based on ldof
163 // orientation, so that the choices are consistent when we visit a
164 // shared dof multiple times
165 int coarse_i, coarse_j;
166 if (i_edge == j_edge) // edof is a vertex or interior
167 {
168 // note that interiors could be done with elements in parallel
169 // (you'd have to rethink numbering but it could be done in advance)
170 coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
171 coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
172 }
173 else // edof is on an edge but not a vertex
174 {
175 // Orient coarse_i, coarse_j based on numbering of ldofs on vertices
176 int left_in_edof, left_in_ldof, right_in_edof, right_in_ldof;
177 if (i_edge)
178 {
179 left_in_edof = i*P1d + 0;
180 right_in_edof = i*P1d + (P1d - 1);
181 left_in_ldof = in_elem_dof[left_in_edof*in_layout[0]
182 + e*in_layout[2]] + rounding_guard;
183 right_in_ldof = in_elem_dof[right_in_edof*in_layout[0]
184 + e*in_layout[2]] + rounding_guard;
185 coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
186 coarse_j = (left_in_ldof < right_in_ldof) ?
187 coarse_1d_edof(j, P1d, coarse_P1d) : reverse_coarse_1d_edof(j, P1d, coarse_P1d);
188 }
189 else
190 {
191 left_in_edof = 0*P1d + j;
192 right_in_edof = (P1d - 1)*P1d + j;
193 left_in_ldof = in_elem_dof[left_in_edof*in_layout[0]
194 + e*in_layout[2]] + rounding_guard;
195 right_in_ldof = in_elem_dof[right_in_edof*in_layout[0]
196 + e*in_layout[2]] + rounding_guard;
197 coarse_i = (left_in_ldof < right_in_ldof) ?
198 coarse_1d_edof(i, P1d, coarse_P1d) : reverse_coarse_1d_edof(i, P1d, coarse_P1d);
199 coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
200 }
201 }
202
203 // Select edof to be on coarse grid and assign numbering and maps
204 if (coarse_i >= 0 && coarse_j >= 0)
205 {
206 int out_edof = coarse_i*coarse_P1d + coarse_j;
207 if (dof_map[in_ldof] >= 0)
208 {
209 out_elem_dof[e*coarse_elemsize + out_edof] = dof_map[in_ldof];
210 }
211 else
212 {
213 out_elem_dof[e*coarse_elemsize + out_edof] = running_out_ldof_count;
214 dof_map[in_ldof] = running_out_ldof_count;
215 running_out_ldof_count++;
216 }
217 }
218 }
219 }
220 }
221 }
222 else if (dim == 3)
223 {
224 // The 3D code is perhaps overly complicated and could be optimized
225 for (int e = 0; e < numelem; ++e)
226 {
227 // Loop over edofs in element
228 for (int i = 0; i < P1d; ++i)
229 {
230 for (int j = 0; j < P1d; ++j)
231 {
232 for (int k = 0; k < P1d; ++k)
233 {
234 // Determine topology; is this edof on the outside of the element
235 // in the i, j, or k direction?
236 int in_edof = i*P1d*P1d + j*P1d + k;
237 int in_ldof = in_elem_dof[in_edof*in_layout[0] +
238 e*in_layout[2]] + rounding_guard;
239 bool i_edge = (i == 0 || i == P1d - 1);
240 bool j_edge = (j == 0 || j == P1d - 1);
241 bool k_edge = (k == 0 || k == P1d - 1);
242 int topo = 0;
243 if (i_edge) { topo++; }
244 if (j_edge) { topo++; }
245 if (k_edge) { topo++; }
246
247 // Determine corresponding coarse 1D edof indices
248 // We do this systematically, orienting edges and faces based on ldof
249 // orientation, so that the choices are consistent when we visit a
250 // shared dof multiple times
251 int coarse_i, coarse_j, coarse_k;
252 if (topo == 0 || topo == 3)
253 {
254 // edof is a vertex or interior
255 coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
256 coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
257 coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
258 }
259 else if (topo == 2)
260 {
261 // edof is on an edge, not a vertex
262 // Orient based on ldof numbering of vertices that define edge
263 int left_in_edof, left_in_ldof, right_in_edof, right_in_ldof;
264 if (!i_edge)
265 {
266 left_in_edof = 0*P1d*P1d + j*P1d + k;
267 right_in_edof = (P1d - 1)*P1d*P1d + j*P1d + k;
268 left_in_ldof = in_elem_dof[left_in_edof*in_layout[0] +
269 e*in_layout[2]] + rounding_guard;
270 right_in_ldof = in_elem_dof[right_in_edof*in_layout[0] +
271 e*in_layout[2]] + rounding_guard;
272 coarse_i = (left_in_ldof < right_in_ldof) ?
273 coarse_1d_edof(i, P1d, coarse_P1d) : reverse_coarse_1d_edof(i, P1d, coarse_P1d);
274 coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
275 coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
276 }
277 else if (!j_edge)
278 {
279 left_in_edof = i*P1d*P1d + 0*P1d + k;
280 right_in_edof = i*P1d*P1d + (P1d - 1)*P1d + k;
281 left_in_ldof = in_elem_dof[left_in_edof*in_layout[0] +
282 e*in_layout[2]] + rounding_guard;
283 right_in_ldof = in_elem_dof[right_in_edof*in_layout[0] +
284 e*in_layout[2]] + rounding_guard;
285 coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
286 coarse_j = (left_in_ldof < right_in_ldof) ?
287 coarse_1d_edof(j, P1d, coarse_P1d) : reverse_coarse_1d_edof(j, P1d, coarse_P1d);
288 coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
289 }
290 else
291 {
292 if (k_edge)
293 {
294 return CeedError(ceed, 1,
295 "Element connectivity does not make sense!");
296 }
297 left_in_edof = i*P1d*P1d + j*P1d + 0;
298 right_in_edof = i*P1d*P1d + j*P1d + (P1d - 1);
299 left_in_ldof = in_elem_dof[left_in_edof*in_layout[0] +
300 e*in_layout[2]] + rounding_guard;
301 right_in_ldof = in_elem_dof[right_in_edof*in_layout[0] +
302 e*in_layout[2]] + rounding_guard;
303 coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
304 coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
305 coarse_k = (left_in_ldof < right_in_ldof) ?
306 coarse_1d_edof(k, P1d, coarse_P1d) : reverse_coarse_1d_edof(k, P1d, coarse_P1d);
307 }
308 }
309 else
310 {
311 // edof is on a face, not an edge
312 // Orient based on four vertices that define the face
313 if (topo != 1)
314 {
315 return CeedError(ceed, 1,
316 "Element connectivity does not match topology!");
317 }
318 int bottom_left_edof, bottom_right_edof, top_left_edof, top_right_edof;
319 int bottom_left_ldof, bottom_right_ldof, top_left_ldof, top_right_ldof;
320 if (i_edge)
321 {
322 bottom_left_edof = i*P1d*P1d + 0*P1d + 0;
323 bottom_right_edof = i*P1d*P1d + 0*P1d + (P1d - 1);
324 top_right_edof = i*P1d*P1d + (P1d - 1)*P1d + (P1d - 1);
325 top_left_edof = i*P1d*P1d + (P1d - 1)*P1d + 0;
326 bottom_left_ldof = in_elem_dof[bottom_left_edof*in_layout[0] +
327 e*in_layout[2]] + rounding_guard;
328 bottom_right_ldof = in_elem_dof[bottom_right_edof*in_layout[0] +
329 e*in_layout[2]] + rounding_guard;
330 top_right_ldof = in_elem_dof[top_right_edof*in_layout[0] +
331 e*in_layout[2]] + rounding_guard;
332 top_left_ldof = in_elem_dof[top_left_edof*in_layout[0] +
333 e*in_layout[2]] + rounding_guard;
334 int m = min4(bottom_left_ldof, bottom_right_ldof, top_right_ldof,
335 top_left_ldof);
336 coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
337 if (m == bottom_left_ldof)
338 {
339 coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
340 coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
341 }
342 else if (m == bottom_right_ldof) // j=0, k=P1d-1
343 {
344 coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
345 coarse_k = reverse_coarse_1d_edof(k, P1d, coarse_P1d);
346 }
347 else if (m == top_right_ldof)
348 {
349 coarse_j = reverse_coarse_1d_edof(j, P1d, coarse_P1d);
350 coarse_k = reverse_coarse_1d_edof(k, P1d, coarse_P1d);
351 }
352 else // j=P1d-1, k=0
353 {
354 coarse_j = reverse_coarse_1d_edof(j, P1d, coarse_P1d);
355 coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
356 }
357 }
358 else if (j_edge)
359 {
360 bottom_left_edof = 0*P1d*P1d + j*P1d + 0;
361 bottom_right_edof = 0*P1d*P1d + j*P1d + (P1d - 1);
362 top_right_edof = (P1d - 1)*P1d*P1d + j*P1d + (P1d - 1);
363 top_left_edof = (P1d - 1)*P1d*P1d + j*P1d + 0;
364 bottom_left_ldof = in_elem_dof[bottom_left_edof*in_layout[0] +
365 e*in_layout[2]] + rounding_guard;
366 bottom_right_ldof = in_elem_dof[bottom_right_edof*in_layout[0] +
367 e*in_layout[2]] + rounding_guard;
368 top_right_ldof = in_elem_dof[top_right_edof*in_layout[0] +
369 e*in_layout[2]] + rounding_guard;
370 top_left_ldof = in_elem_dof[top_left_edof*in_layout[0] +
371 e*in_layout[2]] + rounding_guard;
372 int m = min4(bottom_left_ldof, bottom_right_ldof, top_right_ldof,
373 top_left_ldof);
374 coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
375 if (m == bottom_left_ldof)
376 {
377 coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
378 coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
379 }
380 else if (m == bottom_right_ldof) // i=0, k=P1d-1
381 {
382 coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
383 coarse_k = reverse_coarse_1d_edof(k, P1d, coarse_P1d);
384 }
385 else if (m == top_right_ldof)
386 {
387 coarse_i = reverse_coarse_1d_edof(i, P1d, coarse_P1d);
388 coarse_k = reverse_coarse_1d_edof(k, P1d, coarse_P1d);
389 }
390 else // i=P1d-1, k=0
391 {
392 coarse_i = reverse_coarse_1d_edof(i, P1d, coarse_P1d);
393 coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
394 }
395 }
396 else
397 {
398 if (!k_edge)
399 {
400 return CeedError(ceed, 1,
401 "Element connectivity does not make sense!");
402 }
403 bottom_left_edof = 0*P1d*P1d + 0*P1d + k;
404 bottom_right_edof = 0*P1d*P1d + (P1d - 1)*P1d + k;
405 top_right_edof = (P1d - 1)*P1d*P1d + (P1d - 1)*P1d + k;
406 top_left_edof = (P1d - 1)*P1d*P1d + 0*P1d + k;
407 bottom_left_ldof = in_elem_dof[bottom_left_edof*in_layout[0] +
408 e*in_layout[2]] + rounding_guard;
409 bottom_right_ldof = in_elem_dof[bottom_right_edof*in_layout[0] +
410 e*in_layout[2]] + rounding_guard;
411 top_right_ldof = in_elem_dof[top_right_edof*in_layout[0] +
412 e*in_layout[2]] + rounding_guard;
413 top_left_ldof = in_elem_dof[top_left_edof*in_layout[0] +
414 e*in_layout[2]] + rounding_guard;
415 int m = min4(bottom_left_ldof, bottom_right_ldof,
416 top_right_ldof, top_left_ldof);
417 coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
418 if (m == bottom_left_ldof)
419 {
420 coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
421 coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
422 }
423 else if (m == bottom_right_ldof) // i=0, j=P1d-1
424 {
425 coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
426 coarse_j = reverse_coarse_1d_edof(j, P1d, coarse_P1d);
427 }
428 else if (m == top_right_ldof)
429 {
430 coarse_i = reverse_coarse_1d_edof(i, P1d, coarse_P1d);
431 coarse_j = reverse_coarse_1d_edof(j, P1d, coarse_P1d);
432 }
433 else // i=P1d-1, j=0
434 {
435 coarse_i = reverse_coarse_1d_edof(i, P1d, coarse_P1d);
436 coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
437 }
438 }
439 }
440
441 // Select edof to be on coarse grid and assign numbering and maps
442 if (coarse_i >= 0 && coarse_j >= 0 && coarse_k >= 0)
443 {
444 int out_edof = coarse_i*coarse_P1d*coarse_P1d + coarse_j*coarse_P1d + coarse_k;
445 if (dof_map[in_ldof] >= 0)
446 {
447 out_elem_dof[e*coarse_elemsize + out_edof] = dof_map[in_ldof];
448 }
449 else
450 {
451 out_elem_dof[e*coarse_elemsize + out_edof] = running_out_ldof_count;
452 dof_map[in_ldof] = running_out_ldof_count;
453 running_out_ldof_count++;
454 }
455 }
456 }
457 }
458 }
459 }
460
461 }
462 else
463 {
464 return CeedError(ceed, 1,
465 "CeedATPMGElemRestriction does not yet support this dimension.");
466 }
467
468 ierr = CeedVectorRestoreArrayRead(in_evec, &in_elem_dof); CeedChk(ierr);
469 ierr = CeedVectorDestroy(&in_evec); CeedChk(ierr);
470
471 ierr = CeedElemRestrictionCreate(ceed, numelem, coarse_elemsize, numcomp,
472 0, running_out_ldof_count,
473 CEED_MEM_HOST, CEED_COPY_VALUES, out_elem_dof,
474 er_out); CeedChk(ierr);
475
476 delete [] out_elem_dof;
477
478 return 0;
479 }
480
481
CeedBasisATPMGCoarseToFine(Ceed ceed,int P1d,int dim,int order_reduction,CeedBasis * basisc2f)482 int CeedBasisATPMGCoarseToFine(Ceed ceed, int P1d, int dim, int order_reduction,
483 CeedBasis *basisc2f)
484 {
485 // this assumes Lobatto nodes on fine and coarse again
486 // (not so hard to generalize, but we would have to write it ourselves instead of
487 // calling the following Ceed function)
488 int ierr;
489 ierr = CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P1d - order_reduction, P1d,
490 CEED_GAUSS_LOBATTO, basisc2f); CeedChk(ierr);
491 return 0;
492 }
493
CeedBasisATPMGCoarseToFine(CeedBasis basisin,CeedBasis * basisc2f,int order_reduction)494 int CeedBasisATPMGCoarseToFine(CeedBasis basisin,
495 CeedBasis *basisc2f,
496 int order_reduction)
497 {
498 int ierr;
499 Ceed ceed;
500 ierr = CeedBasisGetCeed(basisin, &ceed); CeedChk(ierr);
501
502 CeedInt dim, P1d;
503 ierr = CeedBasisGetDimension(basisin, &dim); CeedChk(ierr);
504 ierr = CeedBasisGetNumNodes1D(basisin, &P1d); CeedChk(ierr);
505 ierr = CeedBasisATPMGCoarseToFine(ceed, P1d, dim, order_reduction,
506 basisc2f); CeedChk(ierr);
507 return 0;
508 }
509
CeedBasisATPMGCoarsen(CeedBasis basisin,CeedBasis basisc2f,CeedBasis * basisout,int order_reduction)510 int CeedBasisATPMGCoarsen(CeedBasis basisin,
511 CeedBasis basisc2f,
512 CeedBasis* basisout,
513 int order_reduction)
514 {
515 int ierr;
516 Ceed ceed;
517 ierr = CeedBasisGetCeed(basisin, &ceed); CeedChk(ierr);
518
519 CeedInt dim, ncomp, P1d, Q1d;
520 ierr = CeedBasisGetDimension(basisin, &dim); CeedChk(ierr);
521 ierr = CeedBasisGetNumComponents(basisin, &ncomp); CeedChk(ierr);
522 ierr = CeedBasisGetNumNodes1D(basisin, &P1d); CeedChk(ierr);
523 ierr = CeedBasisGetNumQuadraturePoints1D(basisin, &Q1d); CeedChk(ierr);
524
525 CeedInt coarse_P1d = P1d - order_reduction;
526
527 const CeedScalar *interp1d;
528 ierr = CeedBasisGetInterp1D(basisin, &interp1d); CeedChk(ierr);
529 const CeedScalar * grad1d;
530 ierr = CeedBasisGetGrad1D(basisin, &grad1d); CeedChk(ierr);
531
532 CeedScalar * coarse_interp1d = new CeedScalar[coarse_P1d * Q1d];
533 CeedScalar * coarse_grad1d = new CeedScalar[coarse_P1d * Q1d];
534 CeedScalar * fine_nodal_points = new CeedScalar[P1d];
535
536 // these things are in [-1, 1], not [0, 1], which matters
537 // (todo: how can we determine this or something related, algebraically?)
538 /* one way you might be able to tell is to just run this algorithm
539 with coarse_P1d = 2 (i.e., linear) and look for symmetry in the coarse
540 basis matrix? */
541 ierr = CeedLobattoQuadrature(P1d, fine_nodal_points, NULL); CeedChk(ierr);
542 for (int i = 0; i < P1d; ++i)
543 {
544 fine_nodal_points[i] = 0.5 * fine_nodal_points[i] + 0.5; // cheating
545 }
546
547 const CeedScalar *interp_ctof;
548 ierr = CeedBasisGetInterp1D(basisc2f, &interp_ctof); CeedChk(ierr);
549
550 for (int i = 0; i < Q1d; ++i)
551 {
552 for (int j = 0; j < coarse_P1d; ++j)
553 {
554 coarse_interp1d[i * coarse_P1d + j] = 0.0;
555 coarse_grad1d[i * coarse_P1d + j] = 0.0;
556 for (int k = 0; k < P1d; ++k)
557 {
558 coarse_interp1d[i * coarse_P1d + j] += interp_ctof[k * coarse_P1d + j] *
559 interp1d[i * P1d + k];
560 coarse_grad1d[i * coarse_P1d + j] += interp_ctof[k * coarse_P1d + j] *
561 grad1d[i * P1d + k];
562 }
563 }
564 }
565
566 const CeedScalar * qref1d;
567 ierr = CeedBasisGetQRef(basisin, &qref1d); CeedChk(ierr);
568 const CeedScalar * qweight1d;
569 ierr = CeedBasisGetQWeights(basisin, &qweight1d); CeedChk(ierr);
570 ierr = CeedBasisCreateTensorH1(ceed, dim, ncomp,
571 coarse_P1d, Q1d, coarse_interp1d, coarse_grad1d,
572 qref1d, qweight1d, basisout); CeedChk(ierr);
573
574 delete [] fine_nodal_points;
575 delete [] coarse_interp1d;
576 delete [] coarse_grad1d;
577
578 return 0;
579 }
580
CeedATPMGOperator(CeedOperator oper,int order_reduction,CeedElemRestriction coarse_er,CeedBasis coarse_basis_in,CeedBasis basis_ctof_in,CeedOperator * out)581 int CeedATPMGOperator(CeedOperator oper, int order_reduction,
582 CeedElemRestriction coarse_er,
583 CeedBasis coarse_basis_in,
584 CeedBasis basis_ctof_in,
585 CeedOperator* out)
586 {
587 int ierr;
588 Ceed ceed;
589 ierr = CeedOperatorGetCeed(oper, &ceed); CeedChk(ierr);
590
591 CeedQFunction qf;
592 ierr = CeedOperatorGetQFunction(oper, &qf); CeedChk(ierr);
593 CeedInt numinputfields, numoutputfields;
594 ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
595 CeedQFunctionField *inputqfields, *outputqfields;
596 ierr = CeedQFunctionGetFields(qf, &inputqfields, &outputqfields); CeedChk(ierr);
597 CeedOperatorField *inputfields, *outputfields;
598 ierr = CeedOperatorGetFields(oper, &inputfields, &outputfields); CeedChk(ierr);
599
600 CeedElemRestriction * er_input = new CeedElemRestriction[numinputfields];
601 CeedElemRestriction * er_output = new CeedElemRestriction[numoutputfields];
602 CeedVector * if_vector = new CeedVector[numinputfields];
603 CeedVector * of_vector = new CeedVector[numoutputfields];
604 CeedBasis * basis_input = new CeedBasis[numinputfields];
605 CeedBasis * basis_output = new CeedBasis[numoutputfields];
606 CeedBasis cbasis = coarse_basis_in;
607
608 int active_input_basis = -1;
609 for (int i = 0; i < numinputfields; ++i)
610 {
611 ierr = CeedOperatorFieldGetElemRestriction(inputfields[i],
612 &er_input[i]); CeedChk(ierr);
613 ierr = CeedOperatorFieldGetVector(inputfields[i], &if_vector[i]); CeedChk(ierr);
614 ierr = CeedOperatorFieldGetBasis(inputfields[i], &basis_input[i]);
615 CeedChk(ierr);
616 if (if_vector[i] == CEED_VECTOR_ACTIVE)
617 {
618 if (active_input_basis < 0)
619 {
620 active_input_basis = i;
621 }
622 else if (basis_input[i] != basis_input[active_input_basis])
623 {
624 return CeedError(ceed, 1, "Two different active input basis!");
625 }
626 }
627 }
628 for (int i = 0; i < numoutputfields; ++i)
629 {
630 ierr = CeedOperatorFieldGetElemRestriction(outputfields[i],
631 &er_output[i]); CeedChk(ierr);
632 ierr = CeedOperatorFieldGetVector(outputfields[i], &of_vector[i]);
633 CeedChk(ierr);
634 ierr = CeedOperatorFieldGetBasis(outputfields[i], &basis_output[i]);
635 CeedChk(ierr);
636 if (of_vector[i] == CEED_VECTOR_ACTIVE)
637 {
638 // should already be coarsened
639 if (basis_output[i] != basis_input[active_input_basis])
640 {
641 return CeedError(ceed, 1, "Input and output basis do not match!");
642 }
643 if (er_output[i] != er_input[active_input_basis])
644 {
645 return CeedError(ceed, 1, "Input and output elem-restriction do not match!");
646 }
647 }
648 }
649
650 CeedOperator coper;
651 ierr = CeedOperatorCreate(ceed, qf, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
652 &coper); CeedChk(ierr);
653
654 for (int i = 0; i < numinputfields; ++i)
655 {
656 char * fieldname;
657 ierr = CeedQFunctionFieldGetName(inputqfields[i], &fieldname); CeedChk(ierr);
658 if (if_vector[i] == CEED_VECTOR_ACTIVE)
659 {
660 ierr = CeedOperatorSetField(coper, fieldname, coarse_er, cbasis,
661 if_vector[i]); CeedChk(ierr);
662 }
663 else
664 {
665 ierr = CeedOperatorSetField(coper, fieldname, er_input[i], basis_input[i],
666 if_vector[i]); CeedChk(ierr);
667 }
668 }
669 for (int i = 0; i < numoutputfields; ++i)
670 {
671 char * fieldname;
672 ierr = CeedQFunctionFieldGetName(outputqfields[i], &fieldname); CeedChk(ierr);
673 if (of_vector[i] == CEED_VECTOR_ACTIVE)
674 {
675 ierr = CeedOperatorSetField(coper, fieldname, coarse_er, cbasis,
676 of_vector[i]); CeedChk(ierr);
677 }
678 else
679 {
680 ierr = CeedOperatorSetField(coper, fieldname, er_output[i], basis_output[i],
681 of_vector[i]); CeedChk(ierr);
682 }
683 }
684 delete [] er_input;
685 delete [] er_output;
686 delete [] if_vector;
687 delete [] of_vector;
688 delete [] basis_input;
689 delete [] basis_output;
690
691 *out = coper;
692 return 0;
693 }
694
CeedOperatorGetActiveBasis(CeedOperator oper,CeedBasis * basis)695 int CeedOperatorGetActiveBasis(CeedOperator oper, CeedBasis *basis)
696 {
697 int ierr;
698 Ceed ceed;
699 ierr = CeedOperatorGetCeed(oper, &ceed); CeedChk(ierr);
700 CeedQFunction qf;
701 ierr = CeedOperatorGetQFunction(oper, &qf); CeedChk(ierr);
702 CeedInt numinputfields, numoutputfields;
703 ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
704 CeedOperatorField *inputfields, *outputfields;
705 ierr = CeedOperatorGetFields(oper, &inputfields, &outputfields); CeedChk(ierr);
706
707 *basis = NULL;
708 for (int i = 0; i < numinputfields; ++i)
709 {
710 CeedVector if_vector;
711 CeedBasis basis_in;
712 ierr = CeedOperatorFieldGetVector(inputfields[i], &if_vector); CeedChk(ierr);
713 ierr = CeedOperatorFieldGetBasis(inputfields[i], &basis_in); CeedChk(ierr);
714 if (if_vector == CEED_VECTOR_ACTIVE)
715 {
716 if (*basis == NULL)
717 {
718 *basis = basis_in;
719 }
720 else if (*basis != basis_in)
721 {
722 return CeedError(ceed, 1, "Two different active input basis!");
723 }
724 }
725 }
726 for (int i = 0; i < numoutputfields; ++i)
727 {
728 CeedVector of_vector;
729 CeedBasis basis_out;
730 ierr = CeedOperatorFieldGetVector(outputfields[i], &of_vector); CeedChk(ierr);
731 ierr = CeedOperatorFieldGetBasis(outputfields[i], &basis_out); CeedChk(ierr);
732 if (of_vector == CEED_VECTOR_ACTIVE)
733 {
734 if (*basis != basis_out)
735 {
736 return CeedError(ceed, 1, "Input and output basis do not match!");
737 }
738 }
739 }
740 return 0;
741 }
742
CeedATPMGOperator(CeedOperator oper,int order_reduction,CeedElemRestriction coarse_er,CeedBasis * coarse_basis_out,CeedBasis * basis_ctof_out,CeedOperator * out)743 int CeedATPMGOperator(CeedOperator oper, int order_reduction,
744 CeedElemRestriction coarse_er,
745 CeedBasis *coarse_basis_out,
746 CeedBasis *basis_ctof_out,
747 CeedOperator *out)
748 {
749 int ierr;
750
751 CeedQFunction qf;
752 ierr = CeedOperatorGetQFunction(oper, &qf); CeedChk(ierr);
753 CeedInt numinputfields, numoutputfields;
754 ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
755 CeedOperatorField *inputfields;
756 ierr = CeedOperatorGetFields(oper, &inputfields, NULL); CeedChk(ierr);
757
758 CeedBasis basis;
759 ierr = CeedOperatorGetActiveBasis(oper, &basis); CeedChk(ierr);
760 ierr = CeedBasisATPMGCoarseToFine(basis, basis_ctof_out, order_reduction);
761 CeedChk(ierr);
762 ierr = CeedBasisATPMGCoarsen(basis, *basis_ctof_out, coarse_basis_out,
763 order_reduction); CeedChk(ierr);
764 ierr = CeedATPMGOperator(oper, order_reduction, coarse_er, *coarse_basis_out,
765 *basis_ctof_out, out); CeedChk(ierr);
766 return 0;
767 }
768
CeedOperatorGetOrder(CeedOperator oper,CeedInt * order)769 int CeedOperatorGetOrder(CeedOperator oper, CeedInt * order)
770 {
771 int ierr;
772
773 CeedOperatorField active_field;
774 ierr = CeedOperatorGetActiveField(oper, &active_field); CeedChk(ierr);
775 CeedBasis basis;
776 ierr = CeedOperatorFieldGetBasis(active_field, &basis); CeedChk(ierr);
777 int P1d;
778 ierr = CeedBasisGetNumNodes1D(basis, &P1d); CeedChk(ierr);
779 *order = P1d - 1;
780
781 return 0;
782 }
783
CeedATPMGBundle(CeedOperator oper,int order_reduction,CeedBasis * coarse_basis_out,CeedBasis * basis_ctof_out,CeedElemRestriction * er_out,CeedOperator * coarse_oper,CeedInt * & dof_map)784 int CeedATPMGBundle(CeedOperator oper, int order_reduction,
785 CeedBasis* coarse_basis_out,
786 CeedBasis* basis_ctof_out,
787 CeedElemRestriction* er_out,
788 CeedOperator* coarse_oper,
789 CeedInt *&dof_map)
790 {
791 int ierr;
792 CeedInt order;
793 ierr = CeedOperatorGetOrder(oper, &order); CeedChk(ierr);
794 CeedElemRestriction ho_er;
795 ierr = CeedOperatorGetActiveElemRestriction(oper, &ho_er); CeedChk(ierr);
796 ierr = CeedATPMGElemRestriction(order, order_reduction, ho_er, er_out, dof_map);
797 CeedChk(ierr);
798 ierr = CeedATPMGOperator(oper, order_reduction, *er_out, coarse_basis_out,
799 basis_ctof_out, coarse_oper); CeedChk(ierr);
800 return 0;
801 }
802
803 } // namespace ceed
804
805 } // namespace mfem
806
807 #endif // MFEM_USE_CEED
808