51 #include "Teuchos_XMLParameterListHelpers.hpp" 54 #include "create_inline_mesh.h" 63 typedef Tpetra::DefaultPlatform::DefaultPlatformType
Platform;
64 typedef Tpetra::MultiVector<double, int, int>
tMVector_t;
72 int main(
int narg,
char *arg[]) {
74 Teuchos::GlobalMPISession mpiSession(&narg, &arg,0);
75 Platform &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
76 RCP<const Teuchos::Comm<int> > CommT = platform.getComm();
78 int me = CommT->getRank();
79 int numProcs = CommT->getSize();
86 std::string xmlMeshInFileName(
"Poisson.xml");
89 Teuchos::CommandLineProcessor cmdp (
false,
false);
90 cmdp.setOption(
"xmlfile", &xmlMeshInFileName,
91 "XML file with PamGen specifications");
92 cmdp.parse(narg, arg);
95 ParameterList inputMeshList;
97 if(xmlMeshInFileName.length()) {
99 cout <<
"\nReading parameter list from the XML file \"" 100 <<xmlMeshInFileName<<
"\" ...\n\n";
102 Teuchos::updateParametersFromXmlFile(xmlMeshInFileName,
103 Teuchos::inoutArg(inputMeshList));
105 inputMeshList.print(cout,2,
true,
true);
110 cout <<
"Cannot read input file: " << xmlMeshInFileName <<
"\n";
115 std::string meshInput = Teuchos::getParameter<std::string>(inputMeshList,
129 if (me == 0) cout <<
"Generating mesh ... \n\n";
132 long long maxInt = 9223372036854775807LL;
133 Create_Pamgen_Mesh(meshInput.c_str(), dim, me, numProcs, maxInt);
136 if (me == 0) cout <<
"Creating mesh adapter ... \n\n";
140 inputAdapter_t ia(*CommT,
"region");
141 inputAdapter_t ia2(*CommT,
"vertex");
142 inputAdapter_t::gno_t
const *adjacencyIds=NULL;
143 inputAdapter_t::lno_t
const *offsets=NULL;
148 int dimension, num_nodes, num_elem;
152 int num_elem_blk, num_node_sets, num_side_sets;
153 error += im_ex_get_init(exoid, title, &dimension, &num_nodes, &num_elem,
154 &num_elem_blk, &num_node_sets, &num_side_sets);
156 int *element_num_map =
new int [num_elem];
157 error += im_ex_get_elem_num_map(exoid, element_num_map);
159 inputAdapter_t::gno_t *node_num_map =
new int [num_nodes];
160 error += im_ex_get_node_num_map(exoid, node_num_map);
162 int *elem_blk_ids =
new int [num_elem_blk];
163 error += im_ex_get_elem_blk_ids(exoid, elem_blk_ids);
165 int *num_nodes_per_elem =
new int [num_elem_blk];
166 int *num_attr =
new int [num_elem_blk];
167 int *num_elem_this_blk =
new int [num_elem_blk];
168 char **elem_type =
new char * [num_elem_blk];
169 int **connect =
new int * [num_elem_blk];
171 for(
int i = 0; i < num_elem_blk; i++){
172 elem_type[i] =
new char [MAX_STR_LENGTH + 1];
173 error += im_ex_get_elem_block(exoid, elem_blk_ids[i], elem_type[i],
174 (
int*)&(num_elem_this_blk[i]),
175 (
int*)&(num_nodes_per_elem[i]),
176 (
int*)&(num_attr[i]));
177 delete[] elem_type[i];
185 for(
int b = 0; b < num_elem_blk; b++) {
186 connect[b] =
new int [num_nodes_per_elem[b]*num_elem_this_blk[b]];
187 error += im_ex_get_elem_conn(exoid, elem_blk_ids[b], connect[b]);
190 delete[] elem_blk_ids;
195 if (ia.availAdjs(primaryEType, adjEType)) {
196 ia.getAdjsView(primaryEType, adjEType, offsets, adjacencyIds);
198 if ((
int)ia.getLocalNumOf(primaryEType) != num_elem) {
199 cout <<
"Number of elements do not match\n";
203 for (
int b = 0; b < num_elem_blk; b++) {
204 for (
int i = 0; i < num_elem_this_blk[b]; i++) {
205 if (offsets[telct + 1] - offsets[telct] != num_nodes_per_elem[b]) {
206 std::cout <<
"Number of adjacencies do not match" << std::endl;
210 for (
int j = 0; j < num_nodes_per_elem[b]; j++) {
211 ssize_t in_list = -1;
213 for(inputAdapter_t::lno_t k=offsets[telct];k<offsets[telct+1];k++) {
214 if(adjacencyIds[k] ==
215 node_num_map[connect[b][i*num_nodes_per_elem[b]+j]-1]) {
222 std::cout <<
"Adjacency missing" << std::endl;
231 if (telct != num_elem) {
232 cout <<
"Number of elements do not match\n";
237 std::cout <<
"Adjacencies not available" << std::endl;
241 primaryEType = ia2.getPrimaryEntityType();
242 adjEType = ia2.getAdjacencyEntityType();
244 if (ia2.availAdjs(primaryEType, adjEType)) {
245 ia2.getAdjsView(primaryEType, adjEType, offsets, adjacencyIds);
247 if ((
int)ia2.getLocalNumOf(primaryEType) != num_nodes) {
248 cout <<
"Number of nodes do not match\n";
253 int *num_adj =
new int[num_nodes];
255 for (
int i = 0; i < num_nodes; i++) {
259 for (
int b = 0; b < num_elem_blk; b++) {
260 for (
int i = 0; i < num_elem_this_blk[b]; i++) {
261 for (
int j = 0; j < num_nodes_per_elem[b]; j++) {
262 ssize_t in_list = -1;
263 ++num_adj[connect[b][i * num_nodes_per_elem[b] + j] - 1];
265 for(inputAdapter_t::lno_t k =
266 offsets[connect[b][i * num_nodes_per_elem[b] + j] - 1];
267 k < offsets[connect[b][i * num_nodes_per_elem[b] + j]]; k++) {
268 if(adjacencyIds[k] == element_num_map[telct]) {
275 std::cout <<
"Adjacency missing" << std::endl;
284 for (
int i = 0; i < num_nodes; i++) {
285 if (offsets[i + 1] - offsets[i] != num_adj[i]) {
286 std::cout <<
"Number of adjacencies do not match" << std::endl;
295 std::cout <<
"Adjacencies not available" << std::endl;
300 if (me == 0) cout <<
"Deleting the mesh ... \n\n";
302 Delete_Pamgen_Mesh();
305 std::cout <<
"PASS" << std::endl;
Defines the PamgenMeshAdapter class.
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
This class represents a mesh.