www.mooseframework.org
PETScDMDAMesh.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "PETScDMDAMesh.h"
11 
12 #include "libmesh/mesh_generation.h"
13 #include "libmesh/string_to_enum.h"
14 #include "libmesh/periodic_boundaries.h"
15 #include "libmesh/periodic_boundary_base.h"
16 #include "libmesh/unstructured_mesh.h"
17 #include "libmesh/partitioner.h"
18 #include "libmesh/metis_csr_graph.h"
19 #include "libmesh/edge_edge2.h"
20 #include "libmesh/edge_edge3.h"
21 #include "libmesh/edge_edge4.h"
22 #include "libmesh/mesh_communication.h"
23 #include "libmesh/remote_elem.h"
24 #include "libmesh/face_quad4.h"
25 #include "libmesh/cell_hex8.h"
26 
27 #include "ExternalPetscSolverApp.h"
28 
29 // C++ includes
30 #include <cmath> // provides round, not std::round (see http://www.cplusplus.com/reference/cmath/round/)
31 
32 registerMooseObject("ExternalPetscSolverApp", PETScDMDAMesh);
33 
34 template <>
35 InputParameters
37 {
38  InputParameters params = validParams<MooseMesh>();
39 
40  MooseEnum elem_types("EDGE2 QUAD4 HEX8"); // no default
41 
42  MooseEnum dims("1=1 2 3", "2");
43  params.addRequiredParam<MooseEnum>(
44  "dim", dims, "The dimension of the mesh to be generated"); // Make this parameter required
45 
46  params.addParam<dof_id_type>("nx", 11, "Number of elements in the X direction");
47  params.addParam<dof_id_type>("ny", 11, "Number of elements in the Y direction");
48  params.addParam<dof_id_type>("nz", 11, "Number of elements in the Z direction");
49  params.addParam<Real>("xmin", 0.0, "Lower X Coordinate of the generated mesh");
50  params.addParam<Real>("ymin", 0.0, "Lower Y Coordinate of the generated mesh");
51  params.addParam<Real>("zmin", 0.0, "Lower Z Coordinate of the generated mesh");
52  params.addParam<Real>("xmax", 1.0, "Upper X Coordinate of the generated mesh");
53  params.addParam<Real>("ymax", 1.0, "Upper Y Coordinate of the generated mesh");
54  params.addParam<Real>("zmax", 1.0, "Upper Z Coordinate of the generated mesh");
55  params.addParam<MooseEnum>("elem_type",
56  elem_types,
57  "The type of element from libMesh to "
58  "generate (default: linear element for "
59  "requested dimension)");
60 
61  params.addParamNamesToGroup("dim", "Main");
62 
63  params.addClassDescription(
64  "Create a line, square, or cube mesh with uniformly spaced memsh using PETSc DMDA.");
65 
66  // This mesh is always distributed
67  params.set<MooseEnum>("parallel_type") = "DISTRIBUTED";
68 
69  return params;
70 }
71 
72 PETScDMDAMesh::PETScDMDAMesh(const InputParameters & parameters)
73  : MooseMesh(parameters),
74  _dim(getParam<MooseEnum>("dim")),
75  _nx(getParam<dof_id_type>("nx")),
76  _ny(getParam<dof_id_type>("ny")),
77  _nz(getParam<dof_id_type>("nz")),
78  _xmin(getParam<Real>("xmin")),
79  _xmax(getParam<Real>("xmax")),
80  _ymin(getParam<Real>("ymin")),
81  _ymax(getParam<Real>("ymax")),
82  _zmin(getParam<Real>("zmin")),
83  _zmax(getParam<Real>("zmax"))
84 {
85  // All generated meshes are regular orthogonal meshes
86  _regular_orthogonal_mesh = true;
87 
88  // We support 2D mesh only at this point. If you need 3D or 1D
89  // Please contact MOOSE developers
90  if (_dim != 2)
91  mooseError("Support 2 dimensional mesh only");
92 
93 #if LIBMESH_HAVE_PETSC
94  // If we are going to couple PETSc external solver,
95  // take a DA from PETSc
96  ExternalPetscSolverApp * petsc_app = dynamic_cast<ExternalPetscSolverApp *>(&_app);
97  if (petsc_app)
98  {
99  TS & ts = petsc_app->getExternalPETScTS();
100  // Retrieve mesh from TS
101  TSGetDM(ts, &_dmda);
102  _need_to_destroy_dmda = false;
103  }
104  // This mesh object can be used independently for any MOOSE apps.
105  // We create one from scratch
106  else
107  {
108  DMDACreate2d(_communicator.get(),
109  DM_BOUNDARY_NONE,
110  DM_BOUNDARY_NONE,
111  DMDA_STENCIL_BOX,
112  _nx,
113  _ny,
114  PETSC_DECIDE,
115  PETSC_DECIDE,
116  1,
117  1,
118  NULL,
119  NULL,
120  &_dmda);
121  // Let DMDA take PETSc options
122  DMSetFromOptions(_dmda);
123  // Finalize mesh setup
124  DMSetUp(_dmda);
125  _need_to_destroy_dmda = true;
126  }
127 #else
128  mooseError("You need PETSc installed to use PETScDMDAMesh");
129 #endif
130 }
131 
132 Real
134 {
135  switch (component)
136  {
137  case 0:
138  return _xmin;
139  case 1:
140  return _dim > 1 ? _ymin : 0;
141  case 2:
142  return _dim > 2 ? _zmin : 0;
143  default:
144  mooseError("Invalid component");
145  }
146 }
147 
148 Real
150 {
151  switch (component)
152  {
153  case 0:
154  return _xmax;
155  case 1:
156  return _dim > 1 ? _ymax : 0;
157  case 2:
158  return _dim > 2 ? _zmax : 0;
159  default:
160  mooseError("Invalid component");
161  }
162 }
163 
164 std::unique_ptr<MooseMesh>
166 {
167  return libmesh_make_unique<PETScDMDAMesh>(*this);
168 }
169 
170 inline dof_id_type
171 node_id_Quad4(const ElemType /*type*/,
172  const dof_id_type nx,
173  const dof_id_type /*ny*/,
174  const dof_id_type i,
175  const dof_id_type j,
176  const dof_id_type /*k*/)
177 
178 {
179  // Transform a grid coordinate (i, j) to its global node ID
180  // This match what PETSc does
181  return i + j * (nx + 1);
182 }
183 
184 void
186  const dof_id_type nx,
187  const dof_id_type ny,
188  const dof_id_type i,
189  const dof_id_type j,
190  const dof_id_type elem_id,
191  const processor_id_type pid,
192  const ElemType type,
193  MeshBase & mesh)
194 {
195  BoundaryInfo & boundary_info = mesh.get_boundary_info();
196 
197 #if LIBMESH_HAVE_PETSC
198  // Mx: number of grid points in x direction for all processors
199  // My: number of grid points in y direction for all processors
200  // xp: number of processors in x direction
201  // yp: number of processors in y direction
202  PetscInt Mx, My, xp, yp;
203  DMDAGetInfo(da,
204  PETSC_IGNORE,
205  &Mx,
206  &My,
207  PETSC_IGNORE,
208  &xp,
209  &yp,
210  PETSC_IGNORE,
211  PETSC_IGNORE,
212  PETSC_IGNORE,
213  PETSC_IGNORE,
214  PETSC_IGNORE,
215  PETSC_IGNORE,
216  PETSC_IGNORE);
217 
218  const PetscInt *lx, *ly;
219  PetscInt *lxo, *lyo;
220  // PETSc-3.8.x or older use PetscDataType
221 #if PETSC_VERSION_LESS_THAN(3, 9, 0)
222  DMGetWorkArray(da, xp + yp + 2, PETSC_INT, &lxo);
223 #else
224  // PETSc-3.9.x or newer use MPI_DataType
225  DMGetWorkArray(da, xp + yp + 2, MPIU_INT, &lxo);
226 #endif
227  // Gets the ranges of indices in the x, y and z direction that are owned by each process
228  // Ranges here are different from what we have in Mat and Vec.
229  // It means how many points each processor holds
230  DMDAGetOwnershipRanges(da, &lx, &ly, NULL);
231  lxo[0] = 0;
232  for (PetscInt i = 0; i < xp; i++)
233  lxo[i + 1] = lxo[i] + lx[i];
234 
235  lyo = lxo + xp + 1;
236  lyo[0] = 0;
237  for (PetscInt i = 0; i < yp; i++)
238  lyo[i + 1] = lyo[i] + ly[i];
239 
240  // Try to calculate processor-grid coordinate (xpid, ypid)
241  PetscInt xpid, ypid, xpidplus, ypidplus;
242  // Finds integer in a sorted array of integers
243  // Loc: the location if found, otherwise -(slot+1)
244  // where slot is the place the value would go
245  PetscFindInt(i, xp + 1, lxo, &xpid);
246 
247  xpid = xpid < 0 ? -xpid - 1 - 1 : xpid;
248 
249  PetscFindInt(i + 1, xp + 1, lxo, &xpidplus);
250 
251  xpidplus = xpidplus < 0 ? -xpidplus - 1 - 1 : xpidplus;
252 
253  PetscFindInt(j, yp + 1, lyo, &ypid);
254 
255  ypid = ypid < 0 ? -ypid - 1 - 1 : ypid;
256 
257  PetscFindInt(j + 1, yp + 1, lyo, &ypidplus);
258 
259  ypidplus = ypidplus < 0 ? -ypidplus - 1 - 1 : ypidplus;
260 #if PETSC_VERSION_LESS_THAN(3, 9, 0)
261  DMRestoreWorkArray(da, xp + yp + 2, PETSC_INT, &lxo);
262 #else
263  DMRestoreWorkArray(da, xp + yp + 2, MPIU_INT, &lxo);
264 #endif
265 #endif
266  // Bottom Left
267  auto node0_ptr = mesh.add_point(Point(static_cast<Real>(i) / nx, static_cast<Real>(j) / ny, 0),
268  node_id_Quad4(type, nx, 0, i, j, 0));
269  node0_ptr->set_unique_id() = node_id_Quad4(type, nx, 0, i, j, 0);
270  node0_ptr->set_id() = node0_ptr->unique_id();
271  // xpid + ypid * xp is the global processor ID
272  node0_ptr->processor_id() = xpid + ypid * xp;
273 
274  // Bottom Right
275  auto node1_ptr =
276  mesh.add_point(Point(static_cast<Real>(i + 1) / nx, static_cast<Real>(j) / ny, 0),
277  node_id_Quad4(type, nx, 0, i + 1, j, 0));
278  node1_ptr->set_unique_id() = node_id_Quad4(type, nx, 0, i + 1, j, 0);
279  node1_ptr->set_id() = node1_ptr->unique_id();
280  node1_ptr->processor_id() = xpidplus + ypid * xp;
281 
282  // Top Right
283  auto node2_ptr =
284  mesh.add_point(Point(static_cast<Real>(i + 1) / nx, static_cast<Real>(j + 1) / ny, 0),
285  node_id_Quad4(type, nx, 0, i + 1, j + 1, 0));
286  node2_ptr->set_unique_id() = node_id_Quad4(type, nx, 0, i + 1, j + 1, 0);
287  node2_ptr->set_id() = node2_ptr->unique_id();
288  node2_ptr->processor_id() = xpidplus + ypidplus * xp;
289 
290  // Top Left
291  auto node3_ptr =
292  mesh.add_point(Point(static_cast<Real>(i) / nx, static_cast<Real>(j + 1) / ny, 0),
293  node_id_Quad4(type, nx, 0, i, j + 1, 0));
294  node3_ptr->set_unique_id() = node_id_Quad4(type, nx, 0, i, j + 1, 0);
295  node3_ptr->set_id() = node3_ptr->unique_id();
296  node3_ptr->processor_id() = xpid + ypidplus * xp;
297 
298  // New an element and attach four nodes to it
299  Elem * elem = new Quad4;
300  elem->set_id(elem_id);
301  elem->processor_id() = pid;
302  elem->set_unique_id() = elem_id;
303  elem = mesh.add_elem(elem);
304  elem->set_node(0) = node0_ptr;
305  elem->set_node(1) = node1_ptr;
306  elem->set_node(2) = node2_ptr;
307  elem->set_node(3) = node3_ptr;
308 
309  // Bottom
310  if (j == 0)
311  boundary_info.add_side(elem, 0, 0);
312 
313  // Right
314  if (i == nx - 1)
315  boundary_info.add_side(elem, 1, 1);
316 
317  // Top
318  if (j == ny - 1)
319  boundary_info.add_side(elem, 2, 2);
320 
321  // Left
322  if (i == 0)
323  boundary_info.add_side(elem, 3, 3);
324 }
325 
326 void
327 set_boundary_names_Quad4(BoundaryInfo & boundary_info)
328 {
329  boundary_info.sideset_name(0) = "bottom";
330  boundary_info.sideset_name(1) = "right";
331  boundary_info.sideset_name(2) = "top";
332  boundary_info.sideset_name(3) = "left";
333 }
334 
335 void
336 get_indices_Quad4(const dof_id_type nx,
337  const dof_id_type /*ny*/,
338  const dof_id_type elem_id,
339  dof_id_type & i,
340  dof_id_type & j)
341 {
342  i = elem_id % nx;
343  j = (elem_id - i) / nx;
344 }
345 
346 dof_id_type
347 elem_id_Quad4(const dof_id_type nx,
348  const dof_id_type /*nx*/,
349  const dof_id_type i,
350  const dof_id_type j,
351  const dof_id_type /*k*/)
352 {
353  return (j * nx) + i;
354 }
355 
356 inline void
357 get_neighbors_Quad4(const dof_id_type nx,
358  const dof_id_type ny,
359  const dof_id_type i,
360  const dof_id_type j,
361  std::vector<dof_id_type> & neighbors)
362 {
363  std::fill(neighbors.begin(), neighbors.end(), Elem::invalid_id);
364 
365  // Bottom
366  if (j != 0)
367  neighbors[0] = elem_id_Quad4(nx, 0, i, j - 1, 0);
368 
369  // Right
370  if (i != nx - 1)
371  neighbors[1] = elem_id_Quad4(nx, 0, i + 1, j, 0);
372 
373  // Top
374  if (j != ny - 1)
375  neighbors[2] = elem_id_Quad4(nx, 0, i, j + 1, 0);
376 
377  // Left
378  if (i != 0)
379  neighbors[3] = elem_id_Quad4(nx, 0, i - 1, j, 0);
380 }
381 
382 void
383 get_ghost_neighbors_Quad4(const dof_id_type nx,
384  const dof_id_type ny,
385  const MeshBase & mesh,
386  std::set<dof_id_type> & ghost_elems)
387 {
388  auto & boundary_info = mesh.get_boundary_info();
389 
390  dof_id_type i, j;
391 
392  std::vector<dof_id_type> neighbors(4);
393 
394  for (auto elem_ptr : mesh.element_ptr_range())
395  {
396  for (unsigned int s = 0; s < elem_ptr->n_sides(); s++)
397  {
398  // No current neighbor
399  if (!elem_ptr->neighbor_ptr(s))
400  {
401  // Not on a boundary
402  if (!boundary_info.n_boundary_ids(elem_ptr, s))
403  {
404  auto elem_id = elem_ptr->id();
405 
406  get_indices_Quad4(nx, 0, elem_id, i, j);
407 
408  get_neighbors_Quad4(nx, ny, i, j, neighbors);
409 
410  ghost_elems.insert(neighbors[s]);
411  }
412  }
413  }
414  }
415 }
416 
417 void
418 add_node_Qua4(dof_id_type nx,
419  dof_id_type ny,
420  dof_id_type i,
421  dof_id_type j,
422  processor_id_type pid,
423  ElemType type,
424  MeshBase & mesh)
425 {
426  // Bottom Left
427  auto node0_ptr = mesh.add_point(Point(static_cast<Real>(i) / nx, static_cast<Real>(j) / ny, 0),
428  node_id_Quad4(type, nx, 0, i, j, 0));
429  node0_ptr->set_unique_id() = node_id_Quad4(type, nx, 0, i, j, 0);
430  node0_ptr->processor_id() = pid;
431 }
432 
433 #if LIBMESH_HAVE_PETSC
434 void
435 build_cube_Quad4(UnstructuredMesh & mesh, DM da, const ElemType type)
436 {
437  const auto pid = mesh.comm().rank();
438 
439  BoundaryInfo & boundary_info = mesh.get_boundary_info();
440  // xs: start grid point (not element) index on local in x direction
441  // ys: start grid point index on local in y direction
442  // xm: number of grid points owned by the local processor in x direction
443  // ym: number of grid points owned by the local processor in y direction
444  // Mx: number of grid points on all processors in x direction
445  // My: number of grid points on all processors in y direction
446  // xp: number of processor cores in x direction
447  // yp: number of processor cores in y direction
448  PetscInt xs, ys, xm, ym, Mx, My, xp, yp;
449 
450  /* Get local grid boundaries */
451  DMDAGetCorners(da, &xs, &ys, PETSC_IGNORE, &xm, &ym, PETSC_IGNORE);
452  DMDAGetInfo(da,
453  PETSC_IGNORE,
454  &Mx,
455  &My,
456  PETSC_IGNORE,
457  &xp,
458  &yp,
459  PETSC_IGNORE,
460  PETSC_IGNORE,
461  PETSC_IGNORE,
462  PETSC_IGNORE,
463  PETSC_IGNORE,
464  PETSC_IGNORE,
465  PETSC_IGNORE);
466 
467  for (PetscInt j = ys; j < ys + ym; j++)
468  for (PetscInt i = xs; i < xs + xm; i++)
469  {
470  // We loop over grid points, but we are
471  // here building elements. So that we just
472  // simply skip the first x and y points since the
473  // number of grid ponts is one more than
474  // the number of grid elements
475  if (!i || !j)
476  continue;
477 
478  dof_id_type ele_id = (i - 1) + (j - 1) * (Mx - 1);
479 
480  add_element_Quad4(da, Mx - 1, My - 1, i - 1, j - 1, ele_id, pid, type, mesh);
481  }
482 
483  // If there is no element at the given processor
484  // We need to manually add all mesh nodes
485  if ((ys == 0 && ym == 1) || (xs == 0 && xm == 1))
486  for (PetscInt j = ys; j < ys + ym; j++)
487  for (PetscInt i = xs; i < xs + xm; i++)
488  add_node_Qua4(Mx, My, i, j, pid, type, mesh);
489 
490  // Need to link up the local elements before we can know what's missing
491  mesh.find_neighbors();
492 
493  mesh.find_neighbors(true);
494 
495  // Set RemoteElem neighbors
496  for (auto & elem_ptr : mesh.element_ptr_range())
497  for (unsigned int s = 0; s < elem_ptr->n_sides(); s++)
498  if (!elem_ptr->neighbor_ptr(s) && !boundary_info.n_boundary_ids(elem_ptr, s))
499  elem_ptr->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
500 
501  set_boundary_names_Quad4(boundary_info);
502  // Already partitioned!
503  mesh.skip_partitioning(true);
504 
505  // No need to renumber or find neighbors - done did it.
506  // Avoid deprecation message/error by _also_ setting
507  // allow_renumbering(false). This is a bit silly, but we want to
508  // catch cases where people are purely using the old "skip"
509  // interface and not the new flag setting one.
510  mesh.allow_renumbering(false);
511  mesh.prepare_for_use(/*skip_renumber (ignored!) = */ false,
512  /*skip_find_neighbors = */ true);
513 }
514 #endif
515 
516 void
518 {
519  // Reference to the libmesh mesh
520  MeshBase & mesh = getMesh();
521 
522  MooseEnum elem_type_enum = getParam<MooseEnum>("elem_type");
523 
524  if (!isParamValid("elem_type"))
525  {
526  // Switching on MooseEnum
527  switch (_dim)
528  {
529  case 2:
530  elem_type_enum = "QUAD4";
531  break;
532 
533  default:
534  mooseError("Does not support dimension ", _dim, "yet");
535  }
536  }
537 
538  _elem_type = Utility::string_to_enum<ElemType>(elem_type_enum);
539 
540  mesh.set_mesh_dimension(_dim);
541  mesh.set_spatial_dimension(_dim);
542 
543  // Switching on MooseEnum
544  switch (_dim)
545  {
546 #if LIBMESH_HAVE_PETSC
547  case 2:
548  build_cube_Quad4(dynamic_cast<UnstructuredMesh &>(getMesh()), _dmda, _elem_type);
549  break;
550 #endif
551  default:
552  mooseError("Does not support dimension ", _dim, "yet");
553  }
554 }
dof_id_type node_id_Quad4(const ElemType, const dof_id_type nx, const dof_id_type, const dof_id_type i, const dof_id_type j, const dof_id_type)
dof_id_type _nx
Number of elements in x, y, z direction.
Definition: PETScDMDAMesh.h:54
void add_element_Quad4(DM da, const dof_id_type nx, const dof_id_type ny, const dof_id_type i, const dof_id_type j, const dof_id_type elem_id, const processor_id_type pid, const ElemType type, MeshBase &mesh)
dof_id_type _ny
Definition: PETScDMDAMesh.h:54
Real component(const SymmTensor &symm_tensor, unsigned int index)
InputParameters validParams< PETScDMDAMesh >()
Definition: PETScDMDAMesh.C:36
ElemType _elem_type
The type of element to build.
Definition: PETScDMDAMesh.h:60
void add_node_Qua4(dof_id_type nx, dof_id_type ny, dof_id_type i, dof_id_type j, processor_id_type pid, ElemType type, MeshBase &mesh)
This is a demo used to demonstrate how to couple an external app through the MOOSE wrapper APP...
virtual std::unique_ptr< MooseMesh > safeClone() const override
void get_indices_Quad4(const dof_id_type nx, const dof_id_type, const dof_id_type elem_id, dof_id_type &i, dof_id_type &j)
virtual Real getMinInDimension(unsigned int component) const override
void get_neighbors_Quad4(const dof_id_type nx, const dof_id_type ny, const dof_id_type i, const dof_id_type j, std::vector< dof_id_type > &neighbors)
void build_cube_Quad4(UnstructuredMesh &mesh, DM da, const ElemType type)
TS & getExternalPETScTS()
Return a time-stepping (TS) component that holds all the ingredients of applicaiton.
Generate a parallel (distributed) mesh from PETSc DMDA.
Definition: PETScDMDAMesh.h:27
void set_boundary_names_Quad4(BoundaryInfo &boundary_info)
dof_id_type elem_id_Quad4(const dof_id_type nx, const dof_id_type, const dof_id_type i, const dof_id_type j, const dof_id_type)
MooseEnum _dim
The dimension of the mesh.
Definition: PETScDMDAMesh.h:51
DM _dmda
Mesh object.
Definition: PETScDMDAMesh.h:66
bool _need_to_destroy_dmda
If DMDA is created on the fly, we should destroy it.
Definition: PETScDMDAMesh.h:63
Real _xmin
The min/max values for x,y,z component.
Definition: PETScDMDAMesh.h:57
virtual void buildMesh() override
void get_ghost_neighbors_Quad4(const dof_id_type nx, const dof_id_type ny, const MeshBase &mesh, std::set< dof_id_type > &ghost_elems)
virtual Real getMaxInDimension(unsigned int component) const override
PETScDMDAMesh(const InputParameters &parameters)
Definition: PETScDMDAMesh.C:72
registerMooseObject("ExternalPetscSolverApp", PETScDMDAMesh)