https://mooseframework.inl.gov
Functions
PETScDMDAMesh.C File Reference

Go to the source code of this file.

Functions

 registerMooseObject ("ExternalPetscSolverApp", PETScDMDAMesh)
 
dof_id_type node_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)
 
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, MeshBase &mesh)
 
void set_boundary_names_Quad4 (BoundaryInfo &boundary_info)
 
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)
 
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)
 
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 get_ghost_neighbors_Quad4 (const dof_id_type nx, const dof_id_type ny, const MeshBase &mesh, std::set< dof_id_type > &ghost_elems)
 
void add_node_Qua4 (dof_id_type nx, dof_id_type ny, dof_id_type i, dof_id_type j, processor_id_type pid, MeshBase &mesh)
 
void build_cube_Quad4 (UnstructuredMesh &mesh, DM da)
 

Function Documentation

◆ add_element_Quad4()

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,
MeshBase &  mesh 
)

Definition at line 78 of file PETScDMDAMesh.C.

Referenced by build_cube_Quad4().

86 {
87  BoundaryInfo & boundary_info = mesh.get_boundary_info();
88 
89  // Mx: number of grid points in x direction for all processors
90  // My: number of grid points in y direction for all processors
91  // xp: number of processors in x direction
92  // yp: number of processors in y direction
93  PetscInt Mx, My, xp, yp;
94  LibmeshPetscCallA(mesh.comm().get(),
95  DMDAGetInfo(da,
96  PETSC_IGNORE,
97  &Mx,
98  &My,
99  PETSC_IGNORE,
100  &xp,
101  &yp,
102  PETSC_IGNORE,
103  PETSC_IGNORE,
104  PETSC_IGNORE,
105  PETSC_IGNORE,
106  PETSC_IGNORE,
107  PETSC_IGNORE,
108  PETSC_IGNORE));
109 
110  const PetscInt *lx, *ly;
111  PetscInt *lxo, *lyo;
112  // PETSc-3.8.x or older use PetscDataType
113 #if PETSC_VERSION_LESS_THAN(3, 9, 0)
114  LibmeshPetscCallA(mesh.comm().get(), DMGetWorkArray(da, xp + yp + 2, PETSC_INT, &lxo));
115 #else
116  // PETSc-3.9.x or newer use MPI_DataType
117  LibmeshPetscCallA(mesh.comm().get(), DMGetWorkArray(da, xp + yp + 2, MPIU_INT, &lxo));
118 #endif
119 
120  // Gets the ranges of indices in the x, y and z direction that are owned by each process
121  // Ranges here are different from what we have in Mat and Vec.
122  // It means how many points each processor holds
123  LibmeshPetscCallA(mesh.comm().get(), DMDAGetOwnershipRanges(da, &lx, &ly, NULL));
124  lxo[0] = 0;
125  for (PetscInt i = 0; i < xp; i++)
126  lxo[i + 1] = lxo[i] + lx[i];
127 
128  lyo = lxo + xp + 1;
129  lyo[0] = 0;
130  for (PetscInt i = 0; i < yp; i++)
131  lyo[i + 1] = lyo[i] + ly[i];
132 
133  // Try to calculate processor-grid coordinate (xpid, ypid)
134  PetscInt xpid, ypid, xpidplus, ypidplus;
135  // Finds integer in a sorted array of integers
136  // Loc: the location if found, otherwise -(slot+1)
137  // where slot is the place the value would go
138  LibmeshPetscCallA(mesh.comm().get(), PetscFindInt(i, xp + 1, lxo, &xpid));
139 
140  xpid = xpid < 0 ? -xpid - 1 - 1 : xpid;
141 
142  LibmeshPetscCallA(mesh.comm().get(), PetscFindInt(i + 1, xp + 1, lxo, &xpidplus));
143 
144  xpidplus = xpidplus < 0 ? -xpidplus - 1 - 1 : xpidplus;
145 
146  LibmeshPetscCallA(mesh.comm().get(), PetscFindInt(j, yp + 1, lyo, &ypid));
147 
148  ypid = ypid < 0 ? -ypid - 1 - 1 : ypid;
149 
150  LibmeshPetscCallA(mesh.comm().get(), PetscFindInt(j + 1, yp + 1, lyo, &ypidplus));
151 
152  ypidplus = ypidplus < 0 ? -ypidplus - 1 - 1 : ypidplus;
153 #if PETSC_VERSION_LESS_THAN(3, 9, 0)
154  LibmeshPetscCallA(mesh.comm().get(), DMRestoreWorkArray(da, xp + yp + 2, PETSC_INT, &lxo));
155 #else
156  LibmeshPetscCallA(mesh.comm().get(), DMRestoreWorkArray(da, xp + yp + 2, MPIU_INT, &lxo));
157 #endif
158 
159  // Bottom Left
160  auto node0_ptr = mesh.add_point(Point(static_cast<Real>(i) / nx, static_cast<Real>(j) / ny, 0),
161  node_id_Quad4(nx, 0, i, j, 0));
162  node0_ptr->set_unique_id(node_id_Quad4(nx, 0, i, j, 0));
163  node0_ptr->set_id() = node0_ptr->unique_id();
164  // xpid + ypid * xp is the global processor ID
165  node0_ptr->processor_id() = xpid + ypid * xp;
166 
167  // Bottom Right
168  auto node1_ptr =
169  mesh.add_point(Point(static_cast<Real>(i + 1) / nx, static_cast<Real>(j) / ny, 0),
170  node_id_Quad4(nx, 0, i + 1, j, 0));
171  node1_ptr->set_unique_id(node_id_Quad4(nx, 0, i + 1, j, 0));
172  node1_ptr->set_id() = node1_ptr->unique_id();
173  node1_ptr->processor_id() = xpidplus + ypid * xp;
174 
175  // Top Right
176  auto node2_ptr =
177  mesh.add_point(Point(static_cast<Real>(i + 1) / nx, static_cast<Real>(j + 1) / ny, 0),
178  node_id_Quad4(nx, 0, i + 1, j + 1, 0));
179  node2_ptr->set_unique_id(node_id_Quad4(nx, 0, i + 1, j + 1, 0));
180  node2_ptr->set_id() = node2_ptr->unique_id();
181  node2_ptr->processor_id() = xpidplus + ypidplus * xp;
182 
183  // Top Left
184  auto node3_ptr =
185  mesh.add_point(Point(static_cast<Real>(i) / nx, static_cast<Real>(j + 1) / ny, 0),
186  node_id_Quad4(nx, 0, i, j + 1, 0));
187  node3_ptr->set_unique_id(node_id_Quad4(nx, 0, i, j + 1, 0));
188  node3_ptr->set_id() = node3_ptr->unique_id();
189  node3_ptr->processor_id() = xpid + ypidplus * xp;
190 
191  // New an element and attach four nodes to it
192  Elem * elem = new Quad4;
193  elem->set_id(elem_id);
194  elem->processor_id() = pid;
195  // Make sure our unique_id doesn't overlap any nodes'
196  elem->set_unique_id(elem_id + (nx + 1) * (ny + 1));
197  elem = mesh.add_elem(elem);
198  elem->set_node(0, node0_ptr);
199  elem->set_node(1, node1_ptr);
200  elem->set_node(2, node2_ptr);
201  elem->set_node(3, node3_ptr);
202 
203  // Bottom
204  if (j == 0)
205  boundary_info.add_side(elem, 0, 0);
206 
207  // Right
208  if (i == nx - 1)
209  boundary_info.add_side(elem, 1, 1);
210 
211  // Top
212  if (j == ny - 1)
213  boundary_info.add_side(elem, 2, 2);
214 
215  // Left
216  if (i == 0)
217  boundary_info.add_side(elem, 3, 3);
218 }
MeshBase & mesh
dof_id_type node_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)
Definition: PETScDMDAMesh.C:65
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")

◆ add_node_Qua4()

void add_node_Qua4 ( dof_id_type  nx,
dof_id_type  ny,
dof_id_type  i,
dof_id_type  j,
processor_id_type  pid,
MeshBase &  mesh 
)

Definition at line 312 of file PETScDMDAMesh.C.

Referenced by build_cube_Quad4().

318 {
319  // Bottom Left
320  auto node0_ptr = mesh.add_point(Point(static_cast<Real>(i) / nx, static_cast<Real>(j) / ny, 0),
321  node_id_Quad4(nx, 0, i, j, 0));
322  node0_ptr->set_unique_id(node_id_Quad4(nx, 0, i, j, 0));
323  node0_ptr->processor_id() = pid;
324 }
MeshBase & mesh
dof_id_type node_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)
Definition: PETScDMDAMesh.C:65
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")

◆ build_cube_Quad4()

void build_cube_Quad4 ( UnstructuredMesh &  mesh,
DM  da 
)

Definition at line 327 of file PETScDMDAMesh.C.

Referenced by PETScDMDAMesh::buildMesh().

328 {
329  const auto pid = mesh.comm().rank();
330 
331  BoundaryInfo & boundary_info = mesh.get_boundary_info();
332  // xs: start grid point (not element) index on local in x direction
333  // ys: start grid point index on local in y direction
334  // xm: number of grid points owned by the local processor in x direction
335  // ym: number of grid points owned by the local processor in y direction
336  // Mx: number of grid points on all processors in x direction
337  // My: number of grid points on all processors in y direction
338  // xp: number of processor cores in x direction
339  // yp: number of processor cores in y direction
340  PetscInt xs, ys, xm, ym, Mx, My, xp, yp;
341 
342  /* Get local grid boundaries */
343  LibmeshPetscCallA(mesh.comm().get(),
344  DMDAGetCorners(da, &xs, &ys, PETSC_IGNORE, &xm, &ym, PETSC_IGNORE));
345  LibmeshPetscCallA(mesh.comm().get(),
346  DMDAGetInfo(da,
347  PETSC_IGNORE,
348  &Mx,
349  &My,
350  PETSC_IGNORE,
351  &xp,
352  &yp,
353  PETSC_IGNORE,
354  PETSC_IGNORE,
355  PETSC_IGNORE,
356  PETSC_IGNORE,
357  PETSC_IGNORE,
358  PETSC_IGNORE,
359  PETSC_IGNORE));
360 
361  for (PetscInt j = ys; j < ys + ym; j++)
362  for (PetscInt i = xs; i < xs + xm; i++)
363  {
364  // We loop over grid points, but we are
365  // here building elements. So that we just
366  // simply skip the first x and y points since the
367  // number of grid ponts is one more than
368  // the number of grid elements
369  if (!i || !j)
370  continue;
371 
372  dof_id_type ele_id = (i - 1) + (j - 1) * (Mx - 1);
373 
374  add_element_Quad4(da, Mx - 1, My - 1, i - 1, j - 1, ele_id, pid, mesh);
375  }
376 
377  // If there is no element at the given processor
378  // We need to manually add all mesh nodes
379  if ((ys == 0 && ym == 1) || (xs == 0 && xm == 1))
380  for (PetscInt j = ys; j < ys + ym; j++)
381  for (PetscInt i = xs; i < xs + xm; i++)
382  add_node_Qua4(Mx, My, i, j, pid, mesh);
383 
384  // Need to link up the local elements before we can know what's missing
385  mesh.find_neighbors();
386 
387  mesh.find_neighbors(true);
388 
389  // Set RemoteElem neighbors
390  for (auto & elem_ptr : mesh.element_ptr_range())
391  for (unsigned int s = 0; s < elem_ptr->n_sides(); s++)
392  if (!elem_ptr->neighbor_ptr(s) && !boundary_info.n_boundary_ids(elem_ptr, s))
393  elem_ptr->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
394 
395  set_boundary_names_Quad4(boundary_info);
396  // Already partitioned!
397  mesh.skip_partitioning(true);
398 
399  // We've already set our own unique_id values; now make sure that
400  // future mesh modifications use subsequent values.
401  mesh.set_next_unique_id(Mx * My + (Mx - 1) * (My - 1));
402 
403  // No need to renumber or find neighbors - done did it.
404  mesh.allow_renumbering(false);
405  mesh.allow_find_neighbors(false);
406  mesh.prepare_for_use();
407  mesh.allow_find_neighbors(true);
408 }
MeshBase & mesh
void set_boundary_names_Quad4(BoundaryInfo &boundary_info)
void add_node_Qua4(dof_id_type nx, dof_id_type ny, dof_id_type i, dof_id_type j, processor_id_type pid, MeshBase &mesh)
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
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, MeshBase &mesh)
Definition: PETScDMDAMesh.C:78
uint8_t dof_id_type

◆ elem_id_Quad4()

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   
)

Definition at line 241 of file PETScDMDAMesh.C.

Referenced by get_neighbors_Quad4().

246 {
247  return (j * nx) + i;
248 }
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")

◆ get_ghost_neighbors_Quad4()

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 
)

Definition at line 277 of file PETScDMDAMesh.C.

281 {
282  auto & boundary_info = mesh.get_boundary_info();
283 
284  dof_id_type i, j;
285 
286  std::vector<dof_id_type> neighbors(4);
287 
288  for (auto elem_ptr : mesh.element_ptr_range())
289  {
290  for (unsigned int s = 0; s < elem_ptr->n_sides(); s++)
291  {
292  // No current neighbor
293  if (!elem_ptr->neighbor_ptr(s))
294  {
295  // Not on a boundary
296  if (!boundary_info.n_boundary_ids(elem_ptr, s))
297  {
298  auto elem_id = elem_ptr->id();
299 
300  get_indices_Quad4(nx, 0, elem_id, i, j);
301 
302  get_neighbors_Quad4(nx, ny, i, j, neighbors);
303 
304  ghost_elems.insert(neighbors[s]);
305  }
306  }
307  }
308  }
309 }
MeshBase & mesh
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)
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)
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
uint8_t dof_id_type

◆ get_indices_Quad4()

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 
)

Definition at line 230 of file PETScDMDAMesh.C.

Referenced by get_ghost_neighbors_Quad4().

235 {
236  i = elem_id % nx;
237  j = (elem_id - i) / nx;
238 }
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")

◆ get_neighbors_Quad4()

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 
)
inline

Definition at line 251 of file PETScDMDAMesh.C.

Referenced by get_ghost_neighbors_Quad4().

256 {
257  std::fill(neighbors.begin(), neighbors.end(), Elem::invalid_id);
258 
259  // Bottom
260  if (j != 0)
261  neighbors[0] = elem_id_Quad4(nx, 0, i, j - 1, 0);
262 
263  // Right
264  if (i != nx - 1)
265  neighbors[1] = elem_id_Quad4(nx, 0, i + 1, j, 0);
266 
267  // Top
268  if (j != ny - 1)
269  neighbors[2] = elem_id_Quad4(nx, 0, i, j + 1, 0);
270 
271  // Left
272  if (i != 0)
273  neighbors[3] = elem_id_Quad4(nx, 0, i - 1, j, 0);
274 }
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)
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")

◆ node_id_Quad4()

dof_id_type node_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   
)
inline

Definition at line 65 of file PETScDMDAMesh.C.

Referenced by add_element_Quad4(), and add_node_Qua4().

71 {
72  // Transform a grid coordinate (i, j) to its global node ID
73  // This match what PETSc does
74  return i + j * (nx + 1);
75 }
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")

◆ registerMooseObject()

registerMooseObject ( "ExternalPetscSolverApp"  ,
PETScDMDAMesh   
)

◆ set_boundary_names_Quad4()

void set_boundary_names_Quad4 ( BoundaryInfo &  boundary_info)

Definition at line 221 of file PETScDMDAMesh.C.

Referenced by build_cube_Quad4().

222 {
223  boundary_info.sideset_name(0) = "bottom";
224  boundary_info.sideset_name(1) = "right";
225  boundary_info.sideset_name(2) = "top";
226  boundary_info.sideset_name(3) = "left";
227 }