MFEM  v4.5.2
Finite element discretization library
solvers-atpmg.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 #include <ceed/backend.h>
19 
20 #include <math.h>
21 // todo: should probably use Ceed memory wrappers instead of calloc/free?
22 #include <stdlib.h>
23 
24 namespace mfem
25 {
26 
27 namespace 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
33 int 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 
43 int 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 
61 int 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); CeedChk(ierr);
90 
91  CeedInt numelem, numcomp, elemsize;
92  CeedSize numnodes;
93  ierr = CeedElemRestrictionGetNumElements(er_in, &numelem); CeedChk(ierr);
94  ierr = CeedElemRestrictionGetLVectorSize(er_in, &numnodes); CeedChk(ierr);
95  ierr = CeedElemRestrictionGetElementSize(er_in, &elemsize); CeedChk(ierr);
96  ierr = CeedElemRestrictionGetNumComponents(er_in, &numcomp); CeedChk(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  CeedChk(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  CeedChk(ierr);
117  for (CeedSize i = 0; i < numnodes; ++i)
118  {
119  lvec_data[i] = (CeedScalar) i;
120  }
121  ierr = CeedVectorRestoreArray(in_lvec, &lvec_data); CeedChk(ierr);
122  CeedInt in_layout[3];
123  ierr = CeedElemRestrictionGetELayout(er_in, &in_layout); CeedChk(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); CeedChk(ierr);
131  ierr = CeedVectorDestroy(&in_lvec); CeedChk(ierr);
132  const CeedScalar * in_elem_dof;
133  ierr = CeedVectorGetArrayRead(in_evec, CEED_MEM_HOST, &in_elem_dof);
134  CeedChk(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); CeedChk(ierr);
473  ierr = CeedVectorDestroy(&in_evec); CeedChk(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); CeedChk(ierr);
479 
480  delete [] out_elem_dof;
481 
482  return 0;
483 }
484 
485 
486 int 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); CeedChk(ierr);
495  return 0;
496 }
497 
498 int CeedBasisATPMGCoarseToFine(CeedBasis basisin,
499  CeedBasis *basisc2f,
500  int order_reduction)
501 {
502  int ierr;
503  Ceed ceed;
504  ierr = CeedBasisGetCeed(basisin, &ceed); CeedChk(ierr);
505 
506  CeedInt dim, P1d;
507  ierr = CeedBasisGetDimension(basisin, &dim); CeedChk(ierr);
508  ierr = CeedBasisGetNumNodes1D(basisin, &P1d); CeedChk(ierr);
509  ierr = CeedBasisATPMGCoarseToFine(ceed, P1d, dim, order_reduction,
510  basisc2f); CeedChk(ierr);
511  return 0;
512 }
513 
514 int 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); CeedChk(ierr);
522 
523  CeedInt dim, ncomp, P1d, Q1d;
524  ierr = CeedBasisGetDimension(basisin, &dim); CeedChk(ierr);
525  ierr = CeedBasisGetNumComponents(basisin, &ncomp); CeedChk(ierr);
526  ierr = CeedBasisGetNumNodes1D(basisin, &P1d); CeedChk(ierr);
527  ierr = CeedBasisGetNumQuadraturePoints1D(basisin, &Q1d); CeedChk(ierr);
528 
529  CeedInt coarse_P1d = P1d - order_reduction;
530 
531  const CeedScalar *interp1d;
532  ierr = CeedBasisGetInterp1D(basisin, &interp1d); CeedChk(ierr);
533  const CeedScalar * grad1d;
534  ierr = CeedBasisGetGrad1D(basisin, &grad1d); CeedChk(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); CeedChk(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); CeedChk(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); CeedChk(ierr);
572  const CeedScalar * qweight1d;
573  ierr = CeedBasisGetQWeights(basisin, &qweight1d); CeedChk(ierr);
574  ierr = CeedBasisCreateTensorH1(ceed, dim, ncomp,
575  coarse_P1d, Q1d, coarse_interp1d, coarse_grad1d,
576  qref1d, qweight1d, basisout); CeedChk(ierr);
577 
578  delete [] fine_nodal_points;
579  delete [] coarse_interp1d;
580  delete [] coarse_grad1d;
581 
582  return 0;
583 }
584 
585 int 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); CeedChk(ierr);
597 
598  CeedQFunction qf;
599  ierr = CeedOperatorGetQFunction(oper, &qf); CeedChk(ierr);
600  CeedInt numinputfields, numoutputfields;
601  CeedQFunctionField *inputqfields, *outputqfields;
602  ierr = CeedQFunctionGetFields(qf, &numinputfields, &inputqfields,
603  &numoutputfields, &outputqfields);
604  CeedChk(ierr);
605  CeedOperatorField *inputfields, *outputfields;
606  ierr = CeedOperatorGetFields(oper, &numinputfields, &inputfields,
607  &numoutputfields, &outputfields);
608  CeedChk(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]); CeedChk(ierr);
623  ierr = CeedOperatorFieldGetVector(inputfields[i], &if_vector[i]); CeedChk(ierr);
624  ierr = CeedOperatorFieldGetBasis(inputfields[i], &basis_input[i]);
625  CeedChk(ierr);
626  if (if_vector[i] == CEED_VECTOR_ACTIVE)
627  {
628  if (active_input_basis < 0)
629  {
630  active_input_basis = i;
631  }
632  else if (basis_input[i] != basis_input[active_input_basis])
633  {
634  return CeedError(ceed, 1, "Two different active input basis!");
635  }
636  }
637  }
638  for (int i = 0; i < numoutputfields; ++i)
639  {
640  ierr = CeedOperatorFieldGetElemRestriction(outputfields[i],
641  &er_output[i]); CeedChk(ierr);
642  ierr = CeedOperatorFieldGetVector(outputfields[i], &of_vector[i]);
643  CeedChk(ierr);
644  ierr = CeedOperatorFieldGetBasis(outputfields[i], &basis_output[i]);
645  CeedChk(ierr);
646  if (of_vector[i] == CEED_VECTOR_ACTIVE)
647  {
648  // should already be coarsened
649  if (basis_output[i] != basis_input[active_input_basis])
650  {
651  return CeedError(ceed, 1, "Input and output basis do not match!");
652  }
653  if (er_output[i] != er_input[active_input_basis])
654  {
655  return CeedError(ceed, 1, "Input and output elem-restriction do not match!");
656  }
657  }
658  }
659 
660  CeedOperator coper;
661  ierr = CeedOperatorCreate(ceed, qf, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
662  &coper); CeedChk(ierr);
663 
664  for (int i = 0; i < numinputfields; ++i)
665  {
666  char * fieldname;
667  ierr = CeedQFunctionFieldGetName(inputqfields[i], &fieldname); CeedChk(ierr);
668  if (if_vector[i] == CEED_VECTOR_ACTIVE)
669  {
670  ierr = CeedOperatorSetField(coper, fieldname, coarse_er, cbasis,
671  if_vector[i]); CeedChk(ierr);
672  }
673  else
674  {
675  ierr = CeedOperatorSetField(coper, fieldname, er_input[i], basis_input[i],
676  if_vector[i]); CeedChk(ierr);
677  }
678  }
679  for (int i = 0; i < numoutputfields; ++i)
680  {
681  char * fieldname;
682  ierr = CeedQFunctionFieldGetName(outputqfields[i], &fieldname); CeedChk(ierr);
683  if (of_vector[i] == CEED_VECTOR_ACTIVE)
684  {
685  ierr = CeedOperatorSetField(coper, fieldname, coarse_er, cbasis,
686  of_vector[i]); CeedChk(ierr);
687  }
688  else
689  {
690  ierr = CeedOperatorSetField(coper, fieldname, er_output[i], basis_output[i],
691  of_vector[i]); CeedChk(ierr);
692  }
693  }
694  delete [] er_input;
695  delete [] er_output;
696  delete [] if_vector;
697  delete [] of_vector;
698  delete [] basis_input;
699  delete [] basis_output;
700 
701  *out = coper;
702  return 0;
703 }
704 
705 int CeedATPMGOperator(CeedOperator oper, int order_reduction,
706  CeedElemRestriction coarse_er,
707  CeedBasis *coarse_basis_out,
708  CeedBasis *basis_ctof_out,
709  CeedOperator *out)
710 {
711  int ierr;
712 
713  CeedQFunction qf;
714  ierr = CeedOperatorGetQFunction(oper, &qf); CeedChk(ierr);
715  CeedInt numinputfields, numoutputfields;
716  CeedOperatorField *inputfields;
717  ierr = CeedOperatorGetFields(oper, &numinputfields, &inputfields,
718  &numoutputfields, NULL);
719  CeedChk(ierr);
720 
721  CeedBasis basis;
722  ierr = CeedOperatorGetActiveBasis(oper, &basis); CeedChk(ierr);
723  ierr = CeedBasisATPMGCoarseToFine(basis, basis_ctof_out, order_reduction);
724  CeedChk(ierr);
725  ierr = CeedBasisATPMGCoarsen(basis, *basis_ctof_out, coarse_basis_out,
726  order_reduction); CeedChk(ierr);
727  ierr = CeedATPMGOperator(oper, order_reduction, coarse_er, *coarse_basis_out,
728  *basis_ctof_out, out); CeedChk(ierr);
729  return 0;
730 }
731 
732 int CeedOperatorGetOrder(CeedOperator oper, CeedInt * order)
733 {
734  int ierr;
735 
736  CeedOperatorField active_field;
737  ierr = CeedOperatorGetActiveField(oper, &active_field); CeedChk(ierr);
738  CeedBasis basis;
739  ierr = CeedOperatorFieldGetBasis(active_field, &basis); CeedChk(ierr);
740  int P1d;
741  ierr = CeedBasisGetNumNodes1D(basis, &P1d); CeedChk(ierr);
742  *order = P1d - 1;
743 
744  return 0;
745 }
746 
747 int CeedATPMGBundle(CeedOperator oper, int order_reduction,
748  CeedBasis* coarse_basis_out,
749  CeedBasis* basis_ctof_out,
750  CeedElemRestriction* er_out,
751  CeedOperator* coarse_oper,
752  CeedInt *&dof_map)
753 {
754  int ierr;
755  CeedInt order;
756  ierr = CeedOperatorGetOrder(oper, &order); CeedChk(ierr);
757  CeedElemRestriction ho_er;
758  ierr = CeedOperatorGetActiveElemRestriction(oper, &ho_er); CeedChk(ierr);
759  ierr = CeedATPMGElemRestriction(order, order_reduction, ho_er, er_out, dof_map);
760  CeedChk(ierr);
761  ierr = CeedATPMGOperator(oper, order_reduction, *er_out, coarse_basis_out,
762  basis_ctof_out, coarse_oper); CeedChk(ierr);
763  return 0;
764 }
765 
766 } // namespace ceed
767 
768 } // namespace mfem
769 
770 #endif // MFEM_USE_CEED
int CeedATPMGOperator(CeedOperator oper, int order_reduction, CeedElemRestriction coarse_er, CeedBasis coarse_basis_in, CeedBasis basis_ctof_in, CeedOperator *out)
int min4(int a, int b, int c, int d)
int reverse_coarse_1d_edof(int i, int P1d, int coarse_P1d)
int CeedBasisATPMGCoarsen(CeedBasis basisin, CeedBasis basisc2f, CeedBasis *basisout, int order_reduction)
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 CeedOperatorGetActiveField(CeedOperator oper, CeedOperatorField *field)
Definition: util.cpp:131
double b
Definition: lissajous.cpp:42
int CeedOperatorGetOrder(CeedOperator oper, CeedInt *order)
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
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.
double a
Definition: lissajous.cpp:41
int dim
Definition: ex24.cpp:53
int coarse_1d_edof(int i, int P1d, int coarse_P1d)