1 // Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC and
2 // other Axom Project Developers. See the top-level LICENSE file for details.
3 //
4 // SPDX-License-Identifier: (BSD-3-Clause)
5
6 /*!
7 * \file
8 *
9 * \brief Illustrates how to construct and use a mint::UniformMesh object.
10 */
11
12 // sphinx_tutorial_basic_example_start
13
14 // sphinx_tutorial_walkthrough_includes_start
15
16 #include "axom/config.hpp" // compile time definitions
17 #include "axom/core/execution/execution_space.hpp" // for execution_space traits
18
19 #include "axom/mint.hpp" // for mint classes and functions
20 #include "axom/core/numerics/Matrix.hpp" // for numerics::Matrix
21
22 // sphinx_tutorial_walkthrough_includes_end
23
24 // namespace aliases
25 namespace mint = axom::mint;
26 namespace numerics = axom::numerics;
27 namespace xargs = mint::xargs;
28
29 using IndexType = axom::IndexType;
30
31 // compile-time switch for execution policy
32 #if defined(AXOM_USE_RAJA) && defined(AXOM_USE_CUDA) && defined(AXOM_USE_UMPIRE)
33 constexpr int NUM_BLOCKS = 512;
34 using ExecPolicy = axom::CUDA_EXEC<NUM_BLOCKS>;
35 #elif defined(AXOM_USE_RAJA) && defined(AXOM_USE_OPENMP)
36 using ExecPolicy = axom::OMP_EXEC;
37 #else
38 using ExecPolicy = axom::SEQ_EXEC;
39 #endif
40
41 constexpr IndexType NUM_NODES_PER_CELL = 4;
42 constexpr double ONE_OVER_4 = 1. / static_cast<double>(NUM_NODES_PER_CELL);
43
44 /*!
45 * \brief Holds command-line arguments
46 */
47 static struct
48 {
49 int res;
50 bool useUnstructured;
51 } Arguments;
52
53 //------------------------------------------------------------------------------
54 // FUNCTION PROTOTYPES
55 //------------------------------------------------------------------------------
56 void parse_args(int argc, char** argv);
57 mint::Mesh* getUniformMesh();
58 mint::Mesh* getUnstructuredMesh();
59
60 //------------------------------------------------------------------------------
61 // PROGRAM MAIN
62 //------------------------------------------------------------------------------
main(int argc,char ** argv)63 int main(int argc, char** argv)
64 {
65 parse_args(argc, argv);
66
67 // sphinx_tutorial_walkthrough_set_memory_start
68 // NOTE: use unified memory if we are using CUDA
69 const int allocID = axom::execution_space<ExecPolicy>::allocatorID();
70 axom::setDefaultAllocator(allocID);
71 // sphinx_tutorial_walkthrough_set_memory_end
72
73 // sphinx_tutorial_walkthrough_construct_mesh_start
74
75 mint::Mesh* mesh =
76 (Arguments.useUnstructured) ? getUnstructuredMesh() : getUniformMesh();
77
78 // sphinx_tutorial_walkthrough_construct_mesh_end
79
80 // sphinx_tutorial_walkthrough_add_fields_start
81
82 // add a cell-centered and a node-centered field
83 double* phi = mesh->createField<double>("phi", mint::NODE_CENTERED);
84 double* hc = mesh->createField<double>("hc", mint::CELL_CENTERED);
85
86 constexpr int NUM_COMPONENTS = 2;
87 double* xc =
88 mesh->createField<double>("xc", mint::CELL_CENTERED, NUM_COMPONENTS);
89
90 // sphinx_tutorial_walkthrough_add_fields_end
91
92 // sphinx_tutorial_walkthrough_compute_hf_start
93
94 // loop over the nodes and evaluate Himmelblaus Function
95 mint::for_all_nodes<ExecPolicy, xargs::xy>(
96 mesh,
97 AXOM_LAMBDA(IndexType nodeIdx, double x, double y) {
98 const double x_2 = x * x;
99 const double y_2 = y * y;
100 const double A = x_2 + y - 11.0;
101 const double B = x + y_2 - 7.0;
102
103 phi[nodeIdx] = A * A + B * B;
104 });
105
106 // sphinx_tutorial_walkthrough_compute_hf_end
107
108 // sphinx_tutorial_walkthrough_cell_centers_start
109
110 // loop over cells and compute cell centers
111 mint::for_all_cells<ExecPolicy, xargs::coords>(
112 mesh,
113 AXOM_LAMBDA(IndexType cellIdx,
114 const numerics::Matrix<double>& coords,
115 const IndexType* nodeIds) {
116 // NOTE: A column vector of the coords matrix corresponds to a nodes coords
117
118 // Sum the cell's nodal coordinates
119 double xsum = 0.0;
120 double ysum = 0.0;
121 double hsum = 0.0;
122
123 const IndexType numNodes = coords.getNumColumns();
124 for(IndexType inode = 0; inode < numNodes; ++inode)
125 {
126 const double* node = coords.getColumn(inode);
127 xsum += node[mint::X_COORDINATE];
128 ysum += node[mint::Y_COORDINATE];
129
130 hsum += phi[nodeIds[inode]];
131 } // END for all cell nodes
132
133 // compute cell centroid by averaging the nodal coordinate sums
134 const IndexType offset = cellIdx * NUM_COMPONENTS;
135 const double invnnodes = 1.f / static_cast<double>(numNodes);
136 xc[offset] = xsum * invnnodes;
137 xc[offset + 1] = ysum * invnnodes;
138
139 hc[cellIdx] = hsum * invnnodes;
140 });
141
142 // sphinx_tutorial_walkthrough_cell_centers_end
143
144 // sphinx_tutorial_walkthrough_vtk_start
145
146 // write the mesh in a VTK file for visualization
147 std::string vtkfile =
148 (Arguments.useUnstructured) ? "unstructured_mesh.vtk" : "uniform_mesh.vtk";
149 mint::write_vtk(mesh, vtkfile);
150
151 // sphinx_tutorial_walkthrough_vtk_end
152
153 delete mesh;
154 mesh = nullptr;
155
156 return 0;
157 }
158
159 //------------------------------------------------------------------------------
160 // FUNCTION PROTOTYPE IMPLEMENTATION
161 //------------------------------------------------------------------------------
parse_args(int argc,char ** argv)162 void parse_args(int argc, char** argv)
163 {
164 Arguments.res = 25;
165 Arguments.useUnstructured = false;
166
167 for(int i = 1; i < argc; ++i)
168 {
169 if(strcmp(argv[i], "--unstructured") == 0)
170 {
171 Arguments.useUnstructured = true;
172 }
173
174 else if(strcmp(argv[i], "--resolution") == 0)
175 {
176 Arguments.res = std::atoi(argv[++i]);
177 }
178
179 } // END for all arguments
180
181 SLIC_ERROR_IF(
182 Arguments.res < 2,
183 "invalid mesh resolution! Please, pick a value greater than 2.");
184 }
185
186 //------------------------------------------------------------------------------
187 // sphinx_tutorial_walkthrough_construct_umesh_start
getUniformMesh()188 mint::Mesh* getUniformMesh()
189 {
190 // construct a N x N grid within a domain defined in [-5.0, 5.0]
191 const double lo[] = {-5.0, -5.0};
192 const double hi[] = {5.0, 5.0};
193 mint::Mesh* m = new mint::UniformMesh(lo, hi, Arguments.res, Arguments.res);
194 return (m);
195 }
196 // sphinx_tutorial_walkthrough_construct_umesh_end
197
198 //------------------------------------------------------------------------------
getUnstructuredMesh()199 mint::Mesh* getUnstructuredMesh()
200 {
201 mint::Mesh* umesh = getUniformMesh();
202 const IndexType umesh_ncells = umesh->getNumberOfCells();
203 const IndexType umesh_nnodes = umesh->getNumberOfNodes();
204
205 const IndexType ncells = umesh_ncells * 4; // split each quad into 4 triangles
206 const IndexType nnodes = umesh_nnodes + umesh_ncells;
207
208 constexpr int DIMENSION = 2;
209 using MeshType = mint::UnstructuredMesh<mint::SINGLE_SHAPE>;
210 MeshType* m = new MeshType(DIMENSION, mint::TRIANGLE, nnodes, ncells);
211 m->resize(nnodes, ncells);
212
213 double* x = m->getCoordinateArray(mint::X_COORDINATE);
214 double* y = m->getCoordinateArray(mint::Y_COORDINATE);
215 IndexType* cells = m->getCellNodesArray();
216
217 // fill coordinates from uniform mesh
218 mint::for_all_nodes<ExecPolicy, xargs::xy>(
219 umesh,
220 AXOM_LAMBDA(IndexType nodeIdx, double nx, double ny) {
221 x[nodeIdx] = nx;
222 y[nodeIdx] = ny;
223 });
224
225 // loop over cells, compute cell centers and fill connectivity
226 mint::for_all_cells<ExecPolicy, xargs::coords>(
227 umesh,
228 AXOM_LAMBDA(IndexType cellIdx,
229 const numerics::Matrix<double>& coords,
230 const IndexType* nodeIds) {
231 // NOTE: A column vector of the coords matrix corresponds to a nodes coords
232
233 // Sum the cell's nodal coordinates
234 double xsum = 0.0;
235 double ysum = 0.0;
236 for(IndexType inode = 0; inode < NUM_NODES_PER_CELL; ++inode)
237 {
238 const double* node = coords.getColumn(inode);
239 xsum += node[mint::X_COORDINATE];
240 ysum += node[mint::Y_COORDINATE];
241 } // END for all cell nodes
242
243 // compute cell centroid by averaging the nodal coordinate sums
244 const IndexType nc = umesh_nnodes + cellIdx; /* centroid index */
245 x[nc] = xsum * ONE_OVER_4;
246 y[nc] = ysum * ONE_OVER_4;
247
248 // triangulate
249 const IndexType& n0 = nodeIds[0];
250 const IndexType& n1 = nodeIds[1];
251 const IndexType& n2 = nodeIds[2];
252 const IndexType& n3 = nodeIds[3];
253
254 const IndexType offset = cellIdx * 12;
255
256 cells[offset] = n0;
257 cells[offset + 1] = nc;
258 cells[offset + 2] = n3;
259
260 cells[offset + 3] = n0;
261 cells[offset + 4] = n1;
262 cells[offset + 5] = nc;
263
264 cells[offset + 6] = n1;
265 cells[offset + 7] = n2;
266 cells[offset + 8] = nc;
267
268 cells[offset + 9] = n2;
269 cells[offset + 10] = n3;
270 cells[offset + 11] = nc;
271 });
272
273 // delete uniform mesh
274 delete umesh;
275 umesh = nullptr;
276
277 return (m);
278 }
279
280 // sphinx_tutorial_basic_example_end
281