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