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 #include "libmesh/petsc_solver_exception.h" 55 mooseError(
" PETSc external solver TS does not exist or this is not a petsc external solver");
58 std::unique_ptr<MooseMesh>
74 return i +
j * (nx + 1);
87 BoundaryInfo & boundary_info =
mesh.get_boundary_info();
93 PetscInt Mx, My, xp, yp;
94 LibmeshPetscCallA(
mesh.comm().get(),
110 const PetscInt *lx, *ly;
113 #if PETSC_VERSION_LESS_THAN(3, 9, 0) 114 LibmeshPetscCallA(
mesh.comm().get(), DMGetWorkArray(da, xp + yp + 2, PETSC_INT, &lxo));
117 LibmeshPetscCallA(
mesh.comm().get(), DMGetWorkArray(da, xp + yp + 2, MPIU_INT, &lxo));
123 LibmeshPetscCallA(
mesh.comm().get(), DMDAGetOwnershipRanges(da, &lx, &ly, NULL));
125 for (PetscInt i = 0; i < xp; i++)
126 lxo[i + 1] = lxo[i] + lx[i];
130 for (PetscInt i = 0; i < yp; i++)
131 lyo[i + 1] = lyo[i] + ly[i];
134 PetscInt xpid, ypid, xpidplus, ypidplus;
138 LibmeshPetscCallA(
mesh.comm().get(), PetscFindInt(i, xp + 1, lxo, &xpid));
140 xpid = xpid < 0 ? -xpid - 1 - 1 : xpid;
142 LibmeshPetscCallA(
mesh.comm().get(), PetscFindInt(i + 1, xp + 1, lxo, &xpidplus));
144 xpidplus = xpidplus < 0 ? -xpidplus - 1 - 1 : xpidplus;
146 LibmeshPetscCallA(
mesh.comm().get(), PetscFindInt(
j, yp + 1, lyo, &ypid));
148 ypid = ypid < 0 ? -ypid - 1 - 1 : ypid;
150 LibmeshPetscCallA(
mesh.comm().get(), PetscFindInt(
j + 1, yp + 1, lyo, &ypidplus));
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));
156 LibmeshPetscCallA(
mesh.comm().get(), DMRestoreWorkArray(da, xp + yp + 2, MPIU_INT, &lxo));
160 auto node0_ptr =
mesh.add_point(Point(static_cast<Real>(i) / nx, static_cast<Real>(
j) / ny, 0),
163 node0_ptr->set_id() = node0_ptr->unique_id();
165 node0_ptr->processor_id() = xpid + ypid * xp;
169 mesh.add_point(Point(static_cast<Real>(i + 1) / nx, static_cast<Real>(
j) / ny, 0),
172 node1_ptr->set_id() = node1_ptr->unique_id();
173 node1_ptr->processor_id() = xpidplus + ypid * xp;
177 mesh.add_point(Point(static_cast<Real>(i + 1) / nx, static_cast<Real>(
j + 1) / ny, 0),
180 node2_ptr->set_id() = node2_ptr->unique_id();
181 node2_ptr->processor_id() = xpidplus + ypidplus * xp;
185 mesh.add_point(Point(static_cast<Real>(i) / nx, static_cast<Real>(
j + 1) / ny, 0),
188 node3_ptr->set_id() = node3_ptr->unique_id();
189 node3_ptr->processor_id() = xpid + ypidplus * xp;
192 Elem * elem =
new Quad4;
193 elem->set_id(elem_id);
194 elem->processor_id() = pid;
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);
205 boundary_info.add_side(elem, 0, 0);
209 boundary_info.add_side(elem, 1, 1);
213 boundary_info.add_side(elem, 2, 2);
217 boundary_info.add_side(elem, 3, 3);
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";
237 j = (elem_id - i) / nx;
255 std::vector<dof_id_type> & neighbors)
257 std::fill(neighbors.begin(), neighbors.end(), Elem::invalid_id);
279 const MeshBase & mesh,
280 std::set<dof_id_type> & ghost_elems)
282 auto & boundary_info =
mesh.get_boundary_info();
286 std::vector<dof_id_type> neighbors(4);
288 for (
auto elem_ptr :
mesh.element_ptr_range())
290 for (
unsigned int s = 0; s < elem_ptr->n_sides(); s++)
293 if (!elem_ptr->neighbor_ptr(s))
296 if (!boundary_info.n_boundary_ids(elem_ptr, s))
298 auto elem_id = elem_ptr->id();
304 ghost_elems.insert(neighbors[s]);
320 auto node0_ptr =
mesh.add_point(Point(static_cast<Real>(i) / nx, static_cast<Real>(
j) / ny, 0),
323 node0_ptr->processor_id() = pid;
329 const auto pid =
mesh.comm().rank();
331 BoundaryInfo & boundary_info =
mesh.get_boundary_info();
340 PetscInt xs, ys, xm, ym, Mx, My, xp, yp;
343 LibmeshPetscCallA(
mesh.comm().get(),
344 DMDAGetCorners(da, &xs, &ys, PETSC_IGNORE, &xm, &ym, PETSC_IGNORE));
345 LibmeshPetscCallA(
mesh.comm().get(),
361 for (PetscInt
j = ys;
j < ys + ym;
j++)
362 for (PetscInt i = xs; i < xs + xm; i++)
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++)
385 mesh.find_neighbors();
387 mesh.find_neighbors(
true);
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));
397 mesh.skip_partitioning(
true);
401 mesh.set_next_unique_id(Mx * My + (Mx - 1) * (My - 1));
404 mesh.allow_renumbering(
false);
405 mesh.allow_find_neighbors(
false);
406 mesh.prepare_for_use();
407 mesh.allow_find_neighbors(
true);
static InputParameters validParams()
static InputParameters validParams()
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)
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
std::unique_ptr< T > copyConstruct(const T &object)
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)
uint8_t processor_id_type
Generate a parallel (distributed) mesh from PETSc DMDA.
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)
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 mooseError(Args &&... args) const
virtual void buildMesh() override
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
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_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 build_cube_Quad4(UnstructuredMesh &mesh, DM da)
PETScDMDAMesh(const InputParameters ¶meters)
const RemoteElem * remote_elem
registerMooseObject("ExternalPetscSolverApp", PETScDMDAMesh)