1 /*
2 * proj1.cpp
3 *
4 * project on a sphere of radius R, delete sets if needed, and delete edges between parts
5 * (created by resolve shared ents)
6 */
7
8 #include "moab/Core.hpp"
9 #include "moab/Interface.hpp"
10 #include <iostream>
11 #include <math.h>
12
13 #include "moab/IntxMesh/IntxUtils.hpp"
14 #include <assert.h>
15 using namespace moab;
16
17 double radius = 1.;// in m: 6371220.
18
main(int argc,char ** argv)19 int main(int argc, char **argv)
20 {
21
22 bool delete_partition_sets = false;
23
24 if (argc < 3)
25 return 1;
26
27 int index = 1;
28 char * input_mesh1 = argv[1];
29 char * output = argv[2];
30 while (index < argc)
31 {
32 if (!strcmp(argv[index], "-R")) // this is for radius to project
33 {
34 radius = atof(argv[++index]);
35 }
36 if (!strcmp(argv[index], "-DS")) // delete partition sets
37 {
38 delete_partition_sets = true;
39 }
40
41 if (!strcmp(argv[index], "-h"))
42 {
43 std::cout << " usage: proj1 <input> <output> -R <value> -DS (delete partition sets)\n";
44 return 1;
45 }
46 index++;
47 }
48
49 Core moab;
50 Interface & mb = moab;
51
52 ErrorCode rval = mb.load_mesh(input_mesh1);
53
54 std::cout << " -R " << radius << " input: " << input_mesh1 <<
55 " output: " << output << "\n";
56
57 Range verts;
58 rval = mb.get_entities_by_dimension(0, 0, verts);
59 if (MB_SUCCESS != rval)
60 return 1;
61
62 double *x_ptr, *y_ptr, *z_ptr;
63 int count;
64 rval = mb.coords_iterate(verts.begin(), verts.end(), x_ptr, y_ptr, z_ptr, count);
65 if (MB_SUCCESS != rval)
66 return 1;
67 assert(count == (int) verts.size()); // should end up with just one contiguous chunk of vertices
68
69 for (int v = 0; v < count; v++) {
70 //EntityHandle v = verts[v];
71 CartVect pos( x_ptr[v], y_ptr[v] , z_ptr[v]);
72 pos = pos/pos.length();
73 pos = radius*pos;
74 x_ptr[v] = pos[0];
75 y_ptr[v] = pos[1];
76 z_ptr[v] = pos[2];
77 }
78
79 Range edges;
80 rval = mb.get_entities_by_dimension(0, 1, edges);
81 if (MB_SUCCESS != rval)
82 return 1;
83 // write edges to a new set, and after that, write the set, delete the edges and the set
84 EntityHandle sf1;
85 rval = mb.create_meshset(MESHSET_SET, sf1);
86 if (MB_SUCCESS != rval)
87 return 1;
88 rval = mb.add_entities(sf1, edges);
89 if (MB_SUCCESS != rval)
90 return 1;
91 rval = mb.write_mesh("edgesOnly.h5m", &sf1, 1);
92 if (MB_SUCCESS != rval)
93 return 1;
94 rval = mb.delete_entities(&sf1, 1);
95 if (MB_SUCCESS != rval)
96 return 1;
97 mb.delete_entities(edges);
98 mb.write_file(output);
99
100 if (delete_partition_sets)
101 {
102 Tag par_tag;
103 rval = mb.tag_get_handle("PARALLEL_PARTITION", par_tag);
104 if (MB_SUCCESS == rval)
105
106 {
107 Range par_sets;
108 rval = mb.get_entities_by_type_and_tag(0, MBENTITYSET, &par_tag, NULL, 1, par_sets,
109 moab::Interface::UNION);
110 if (!par_sets.empty())
111 mb.delete_entities(par_sets);
112 mb.tag_delete(par_tag);
113 }
114 }
115
116
117 return 0;
118 }
119