Branch data Line data Source code
1 : : #include "meshkit/PostBL.hpp"
2 : :
3 : : namespace MeshKit
4 : : {
5 : :
6 : : // static registration of this mesh scheme
7 : : moab::EntityType PostBL_tps[] = {moab::MBHEX,
8 : : moab::MBMAXTYPE};
9 : 40 : const moab::EntityType* PostBL::output_types()
10 : 40 : { return PostBL_tps; }
11 : :
12 : 1 : PostBL::PostBL( MKCore *mk, const MEntVector &me_vec)
13 : : : MeshScheme( mk, me_vec),
14 [ + - ][ + - ]: 1 : igeom(mk->igeom_instance()), imesh(mk->imesh_instance()),
15 [ + - ][ + - ]: 2 : mb (mk->moab_instance())
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
16 : : // ---------------------------------------------------------------------------
17 : : //! Function: Constructor \n
18 : : //! Input: Initialize mesh and geometry instances and parameters \n
19 : : //! Output: none
20 : : // ---------------------------------------------------------------------------
21 : : {
22 : 1 : tri_sch = 2;
23 : 1 : m_Conn = 0;
24 : 1 : m_BElemNodes = 0;
25 : 1 : m_SurfId = -1;
26 : 1 : check_bl_edge_length = false;
27 : 1 : debug = false;
28 : 1 : hybrid = false;
29 : 1 : m_NeumannSet = -1;
30 : 1 : m_Material = 999999;
31 : 1 : m_nLineNumber = 0;
32 [ + - ]: 1 : szComment = "!";
33 : 1 : MAXCHARS = 300;
34 : 1 : m_JacCalls = 0;
35 : 1 : m_JLo = 0.0;
36 : 1 : m_JHi = 0.0;
37 : 1 : err = 0;
38 : 1 : fixmat = -1;
39 : 1 : hex27 = 0;
40 : : // hex = NULL;
41 : : // hex1 = NULL;
42 : : // hex2 = NULL;
43 : 1 : }
44 : :
45 : 3 : PostBL::~PostBL()
46 : : // ---------------------------------------------------------------------------
47 : : //! Function: Destructor, does nothing..\n
48 : : //! Input: none \n
49 : : //! Output: none \n
50 : : // ---------------------------------------------------------------------------
51 [ - + ]: 2 : {}
52 : :
53 : 0 : bool PostBL::add_modelent(ModelEnt *model_ent)
54 : : // ---------------------------------------------------------------------------
55 : : //! Function: Adds entities for PosBL graph node.\n
56 : : //! Input: ModelEnt \n
57 : : //! Output: none \n
58 : : // ---------------------------------------------------------------------------
59 : : {
60 : 0 : return MeshOp::add_modelent(model_ent);
61 : : }
62 : :
63 : 1 : void PostBL::setup_this()
64 : : // ---------------------------------------------------------------------------
65 : : //! Function: Setup the graph node for PostBL \n
66 : : //! Input: none \n
67 : : //! Output: none \n
68 : : // ---------------------------------------------------------------------------
69 : : {
70 [ - + ]: 1 : if (debug) {
71 : 0 : m_LogFile << "\nIn setup this : " << std::endl;
72 : : }
73 : :
74 : 1 : }
75 : :
76 : 1 : void PostBL::execute_this()
77 : : // ---------------------------------------------------------------------------
78 : : //! Function: Read user input from file and run the PostBL algorithm \n
79 : : //! Input: Uses the file name (.inp) with keywords predefined by PosBL algorithm. \n
80 : : //! Output: Resulting mesh file is saved. \n
81 : : // ---------------------------------------------------------------------------
82 : : {
83 [ + - ][ + - ]: 1 : m_LogFile << "\nIn execute this : creating boundary layer elements.." << std::endl;
84 : : // start the timer
85 [ + - ]: 1 : CClock Timer;
86 : 1 : clock_t sTime = clock();
87 [ + - ]: 2 : std::string szDateTime;
88 [ + - ]: 1 : Timer.GetDateTime (szDateTime);
89 [ + - ]: 2 : VerdictWrapper vw(mb);
90 : :
91 [ + - ][ + - ]: 1 : m_LogFile << "\nStarting out at : " << szDateTime << std::endl;
[ + - ][ + - ]
92 [ + - ][ + - ]: 1 : m_LogFile << "\n Loading meshfile: " << m_MeshFile << std::endl;
[ + - ][ + - ]
93 : :
94 : : // load specified mesh file
95 [ + - ][ + - ]: 1 : IBERRCHK(imesh->load(0, m_MeshFile.c_str(),0), *imesh);
96 : :
97 : : //check if intervals is read from input file
98 [ - + ]: 1 : if(m_Intervals == 0){
99 [ # # ][ # # ]: 0 : m_LogFile << "Please specify desired number of intervals using 'Intervals' keyword" << std::endl;
100 : 0 : exit(0);
101 : : }
102 : : // if no NS specified, take all d-1 elements as BL input
103 [ - + ][ # # ]: 1 : if (m_NeumannSet == -1 && m_SurfId == -1) {
104 : : // create a neumann set with id 99999 in the model and set it for BL generation
105 [ # # ]: 0 : moab::Range neuEnts;
106 : : moab::EntityHandle neuEntSet;
107 : : moab::Tag neuTag;
108 [ # # ]: 0 : mb->create_meshset(moab::MESHSET_SET, neuEntSet);
109 [ # # ][ # # ]: 0 : mb->get_entities_by_dimension(mb->get_root_set(), 2, neuEnts, true);
110 [ # # ]: 0 : mb->add_entities(neuEntSet, neuEnts);
111 [ # # ]: 0 : mb->tag_get_handle("NEUMANN_SET", neuTag);
112 : 0 : m_NeumannSet = 999999;
113 [ # # ]: 0 : mb->tag_set_data(neuTag, &neuEntSet, 1, (void*) &m_NeumannSet);
114 : : }
115 : :
116 [ + - ][ + - ]: 2 : moab::Range all_elems, all_verts;
117 [ + - ][ - + ]: 1 : MBERRCHK(mb->get_entities_by_dimension(0, 3, all_elems,true),mb);
[ # # ][ # # ]
[ # # ][ # # ]
118 [ + - ][ - + ]: 1 : if (all_elems.size() == 0)
119 : 0 : m_GD = 2;
120 [ + - ][ + - ]: 1 : else if (all_elems.size() > 0)
121 : 1 : m_GD = 3;
122 : : else
123 : 0 : exit(0);
124 [ + - ]: 1 : all_elems.clear();
125 [ + - ][ + - ]: 1 : m_LogFile << "Geometric dimension of meshfile = "<< m_GD <<std::endl;
[ + - ]
126 : :
127 : : // obtain existing tag handles
128 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_get_handle("GEOM_DIMENSION", 1, moab::MB_TYPE_INTEGER, GDTag),mb);
[ # # ][ # # ]
[ # # ][ # # ]
129 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_get_handle("NEUMANN_SET", 1, moab::MB_TYPE_INTEGER, NTag),mb);
[ # # ][ # # ]
[ # # ][ # # ]
130 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_get_handle("MATERIAL_SET", 1, moab::MB_TYPE_INTEGER, MTag),mb);
[ # # ][ # # ]
[ # # ][ # # ]
131 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_get_handle("GLOBAL_ID", 1, moab::MB_TYPE_INTEGER, GIDTag),mb);
[ # # ][ # # ]
[ # # ][ # # ]
132 : : // create smoothset and fixed tag for mesquite
133 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_get_handle("SMOOTHSET", 1, moab::MB_TYPE_INTEGER, STag,
[ # # ][ # # ]
[ # # ][ # # ]
134 : : moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT),mb);
135 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_get_handle("fixed", 1, moab::MB_TYPE_INTEGER, FTag,
[ # # ][ # # ]
[ # # ][ # # ]
136 : : moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT),mb);
137 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_get_handle("mnode", 1, moab::MB_TYPE_INTEGER, MNTag,
[ # # ][ # # ]
[ # # ][ # # ]
138 : : moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT),mb);
139 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_get_handle("matid", 1, moab::MB_TYPE_INTEGER, MatIDTag,
[ # # ][ # # ]
[ # # ][ # # ]
140 : : moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT),mb);
141 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_get_handle("BLNODEID", 1, moab::MB_TYPE_INTEGER, BLNodeIDTag,
[ # # ][ # # ]
[ # # ][ # # ]
142 : : moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT),mb);
143 : :
144 : : // get all the entity sets with boundary layer geom dimension, neumann sets and material sets
145 : 1 : m_BLDim = m_GD - 1;
146 : :
147 : 1 : const void* gdim[] = {&m_BLDim};
148 : :
149 [ + - ][ - + ]: 1 : MBERRCHK(mb->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &GDTag,
[ # # ][ # # ]
[ # # ][ # # ]
150 : : gdim, 1 , sets, moab::Interface::INTERSECT, false), mb);
151 [ + - ][ - + ]: 1 : MBERRCHK(mb->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &NTag, 0, 1 , n_sets),mb);
[ # # ][ # # ]
[ # # ][ # # ]
152 [ + - ][ - + ]: 1 : MBERRCHK(mb->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &MTag, 0, 1 , m_sets),mb);
[ # # ][ # # ]
[ # # ][ # # ]
153 : :
154 : : // Handling NeumannSets (if BL surf in input via NS)
155 : 1 : moab::EntityHandle this_set = 0;
156 [ + - ][ # # ]: 1 : for (set_it = n_sets.begin(); set_it != n_sets.end(); set_it++) {
[ + - ][ + - ]
[ + - ]
157 : :
158 [ + - ]: 1 : this_set = *set_it;
159 : :
160 : : // get entity handle of NS specified in the input file
161 : : int set_id;
162 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_get_data(NTag, &this_set, 1, &set_id), mb);
[ # # ][ # # ]
[ # # ][ # # ]
163 [ + - ]: 1 : if(set_id == m_NeumannSet)
164 : 1 : break;
165 : 0 : this_set = 0;
166 : : }
167 [ - + ][ # # ]: 1 : if (debug && m_NeumannSet != -1 && this_set != 0){
[ # # ]
168 [ # # ][ # # ]: 0 : m_LogFile << "Looking for NS with id " << m_NeumannSet <<
169 [ # # ][ # # ]: 0 : ". Total NS found are: "<< n_sets.size() << std::endl;
[ # # ][ # # ]
170 : : }
171 : :
172 : : // For specified surface: get the all the quads and nodes in a range
173 : : moab::EntityHandle s1;
174 : : // variable to store global id of boundary layer specified in the input file
175 : : int dims;
176 : :
177 : : // Method 1: INPUT by NeumannSet
178 [ + - ][ + - ]: 1 : if(m_NeumannSet != -1 && this_set != 0){
179 [ + - ][ - + ]: 1 : MBERRCHK(mb->get_entities_by_dimension(this_set, m_BLDim, quads,true),mb);
[ # # ][ # # ]
[ # # ][ # # ]
180 [ + - ][ - + ]: 1 : if (quads.size() <=0){
181 [ # # ][ # # ]: 0 : m_LogFile << " No quads found, aborting.. " << std::endl;
182 : 0 : exit(0);
183 : : }
184 : :
185 [ + - ][ - + ]: 1 : MBERRCHK(mb->get_adjacencies(quads, 0, false, nodes, moab::Interface::UNION),mb);
[ # # ][ # # ]
[ # # ][ # # ]
186 [ + - ]: 1 : if(m_GD == 3)
187 [ + - ][ - + ]: 1 : MBERRCHK(mb->get_adjacencies(quads, 1, true, edges, moab::Interface::UNION),mb);
[ # # ][ # # ]
[ # # ][ # # ]
188 : :
189 [ - + ]: 1 : if (debug) {
190 [ # # ][ # # ]: 0 : m_LogFile << "Found NeumannSet with id : " << m_NeumannSet << std::endl;
[ # # ]
191 [ # # ][ # # ]: 0 : m_LogFile << "#Quads in this surface: " << quads.size() << std::endl;
[ # # ][ # # ]
192 [ # # ][ # # ]: 0 : m_LogFile << "#Nodes in this surface: " << nodes.size() << std::endl;
[ # # ][ # # ]
193 [ # # ][ # # ]: 0 : m_LogFile << "#New nodes to be created:" << m_Intervals*nodes.size() << std::endl;
[ # # ][ # # ]
194 : 1 : }
195 : : }
196 : : // Method 2: INPUT by surface id (geom dimension)
197 [ # # ]: 0 : else if (m_SurfId !=-1){
198 [ # # ][ # # ]: 0 : for(moab::Range::iterator rit=sets.begin(); rit != sets.end(); ++rit){
[ # # ][ # # ]
[ # # ]
199 [ # # ]: 0 : s1 = *rit;
200 [ # # ][ # # ]: 0 : MBERRCHK(mb->tag_get_data(GIDTag, &s1, 1, &dims),mb);
[ # # ][ # # ]
[ # # ][ # # ]
201 : :
202 [ # # ][ # # ]: 0 : if(dims == m_SurfId && m_SurfId != -1){
203 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_entities_by_dimension(s1, m_BLDim, quads,true),mb);
[ # # ][ # # ]
[ # # ][ # # ]
204 [ # # ][ # # ]: 0 : if (quads.size() <=0){
205 [ # # ][ # # ]: 0 : m_LogFile << " No quads found, aborting.. " << std::endl;
206 : 0 : exit(0);
207 : : }
208 : :
209 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_adjacencies(quads, 0, false, nodes, moab::Interface::UNION),mb);
[ # # ][ # # ]
[ # # ][ # # ]
210 [ # # ]: 0 : if(m_GD == 3)
211 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_adjacencies(edges, 0, false, nodes, moab::Interface::UNION),mb);
[ # # ][ # # ]
[ # # ][ # # ]
212 [ # # ]: 0 : if (debug) {
213 [ # # ][ # # ]: 0 : m_LogFile << "Found surface with id : " << m_SurfId << std::endl;
[ # # ]
214 [ # # ][ # # ]: 0 : m_LogFile << "#Quads in this surface: " << quads.size() << std::endl;
[ # # ][ # # ]
215 [ # # ][ # # ]: 0 : m_LogFile << "#Nodes in this surface: " << nodes.size() << std::endl;
[ # # ][ # # ]
216 [ # # ][ # # ]: 0 : m_LogFile << "#New nodes to be created:" << m_Intervals*nodes.size() << std::endl;
[ # # ][ # # ]
217 : : }
218 : : }
219 : : }
220 : : }
221 : :
222 [ + - ][ + - ]: 1 : if (quads.size() == 0 || nodes.size() == 0) {
[ + - ][ - + ]
[ - + ]
223 [ # # ][ # # ]: 0 : m_LogFile << "Invalid boundary layer specification, aborting.." << std::endl;
224 : 0 : exit(0);
225 : : }
226 : :
227 : : // placeholder for storing smoothing entities
228 : : moab::EntityHandle smooth_set;
229 : 1 : int s_id = 100;
230 [ + - ][ - + ]: 1 : MBERRCHK(mb->create_meshset(moab::MESHSET_SET, smooth_set, 1), mb);
[ # # ][ # # ]
[ # # ][ # # ]
231 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_set_data(STag, &smooth_set, 1, &s_id), mb);
[ # # ][ # # ]
[ # # ][ # # ]
232 : :
233 : : // placeholder for storing gd on new entities
234 [ + - ][ - + ]: 1 : MBERRCHK(mb->create_meshset(moab::MESHSET_SET, geom_set, 1), mb);
[ # # ][ # # ]
[ # # ][ # # ]
235 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_set_data(GDTag, &geom_set, 1, &m_GD), mb);
[ # # ][ # # ]
[ # # ][ # # ]
236 : :
237 : : // declare variables before starting BL creation
238 [ + - ]: 2 : std::vector <bool> node_status(false); // size of verts of bl surface
239 [ + - ][ + - ]: 1 : node_status.resize(nodes.size());
240 [ + - ][ + - ]: 1 : blmaterial_id.resize(quads.size());
241 : :
242 : : //size of the following is based on element type
243 [ + - ][ + - ]: 1 : new_vert.resize(m_Intervals*nodes.size());
244 : : // element on the boundary
245 [ + - ]: 1 : Range::iterator kter = quads.begin();
246 : :
247 [ + - ][ + - ]: 1 : MBERRCHK(mb->get_adjacencies(&(*kter), 1, m_GD, false, old_hex),mb);
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
248 [ + - ][ + - ]: 1 : std::cout << "old_hex size is - before filtering fixmat - " << old_hex.size() << std::endl;
[ + - ]
249 [ - + ]: 1 : if((int) old_hex.size() == 0){
250 [ # # ]: 0 : m_LogFile << "unable to find adjacent hex for BL quad, aborting...";
251 : 0 : exit(0);
252 : : }
253 : :
254 : : // allocate space for connectivity/adjacency during the first pass of this loop
255 [ + - ][ + - ]: 1 : if(mb->type_from_handle(old_hex[0]) == moab::MBHEX){
[ + - ]
256 : 1 : m_Conn = 8;
257 : 1 : m_BElemNodes = 4;
258 : 1 : m_HConn = 8;
259 : : //allocating based on element type
260 [ + - ][ + - ]: 1 : conn.resize(m_Intervals*m_Conn), qconn.resize(m_BElemNodes), adj_qconn.resize(m_BElemNodes),
[ + - ]
261 [ + - ][ + - ]: 1 : old_hex_conn.resize(m_Conn), adj_hex_nodes1.resize(m_Conn);
262 : : }
263 [ # # ][ # # ]: 0 : else if(mb->type_from_handle(old_hex[0]) == MBTET){
[ # # ]
264 : 0 : m_Conn = 4;
265 : 0 : m_BElemNodes = 3;
266 : : //allocating based on element type - thrice the number of elements
267 [ # # ]: 0 : if(hybrid){
268 : 0 : m_HConn = 6;
269 [ # # ]: 0 : conn.resize(m_Intervals*6);
270 : : }
271 : : else{
272 : 0 : m_HConn = 6; // we use the prism for making tets
273 [ # # ]: 0 : tet_conn.resize(3*m_Intervals*m_Conn);
274 [ # # ]: 0 : conn.resize(m_Intervals*6);
275 : : }
276 [ # # ][ # # ]: 0 : qconn.resize(m_BElemNodes), adj_qconn.resize(m_BElemNodes),
277 [ # # ][ # # ]: 0 : old_hex_conn.resize(m_Conn), adj_hex_nodes1.resize(m_Conn);
278 : : }
279 [ # # ][ # # ]: 0 : else if(mb->type_from_handle(old_hex[0]) == MBQUAD){
[ # # ]
280 : 0 : m_Conn = 4;
281 : 0 : m_HConn = 4;
282 : 0 : m_BElemNodes = 2;
283 : : //allocating based on element type
284 [ # # ][ # # ]: 0 : conn.resize(m_Intervals*m_Conn), qconn.resize(m_BElemNodes), adj_qconn.resize(m_BElemNodes),
[ # # ]
285 [ # # ][ # # ]: 0 : old_hex_conn.resize(m_Conn), adj_hex_nodes1.resize(m_Conn);
286 : : }
287 [ # # ][ # # ]: 0 : else if(mb->type_from_handle(old_hex[0]) == MBTRI){
[ # # ]
288 : 0 : m_Conn = 3;
289 : 0 : m_HConn = 4;
290 : 0 : m_BElemNodes = 2;
291 : : //allocating based on element type - twice the number of elements
292 [ # # ]: 0 : if(hybrid){
293 : 0 : m_HConn = 4;
294 [ # # ]: 0 : conn.resize(m_Intervals*m_HConn);
295 : : }
296 : : else{
297 [ # # ]: 0 : tri_conn.resize(2*m_Intervals*m_Conn);
298 [ # # ]: 0 : conn.resize(2*m_Intervals*m_Conn);
299 : : }
300 [ # # ][ # # ]: 0 : qconn.resize(m_BElemNodes), adj_qconn.resize(m_BElemNodes),
301 [ # # ][ # # ]: 0 : old_hex_conn.resize(m_Conn), adj_hex_nodes1.resize(m_Conn);
302 : : }
303 [ # # ][ # # ]: 0 : else if(m_Conn == 0 || m_BElemNodes == 0){
304 [ # # ][ # # ]: 0 : m_LogFile << "This mesh type is not supported by this tool" << std::endl;
305 : 0 : exit(0);
306 : : }
307 : :
308 : : // Tag all nodes on outer boundary with a unique number
309 : 1 : int node_id = 0;
310 [ + - ][ + - ]: 2 : std::vector<int> NId(nodes.size());
311 [ + - ][ + - ]: 9 : for(moab::Range::iterator nodes_iter = nodes.begin(); nodes_iter != nodes.end(); nodes_iter++){
[ + - ][ + - ]
[ + + ]
312 [ + - ]: 8 : NId[node_id] = node_id;
313 [ + - ][ + - ]: 8 : MBERRCHK(mb->tag_set_data(BLNodeIDTag, &(*nodes_iter),1, &NId[node_id]), mb);
[ + - ][ - + ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
314 : 8 : ++node_id;
315 : : }
316 : :
317 : : // Handling MaterialSet
318 : 1 : int mset_id = 0, found = 0;
319 [ + - ][ + - ]: 2 : for (mset_it = m_sets.begin(); mset_it != m_sets.end(); mset_it++) {
[ + - ][ + - ]
[ + + ]
320 : :
321 [ + - ]: 1 : mthis_set = *mset_it;
322 : :
323 : : // get entity handle of MS specified in the input file
324 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_get_data(MTag, &mthis_set, 1, &mset_id), mb);
[ # # ][ # # ]
[ # # ][ # # ]
325 : :
326 : : // if no material set is specified, we'll have to resolve else just set all the MNTag to 0
327 [ - + ]: 1 : if(mset_id == m_Material){
328 : 0 : found = 1;
329 : 0 : break;
330 : : }
331 [ - + ]: 1 : else if(mset_id == fixmat){
332 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_entities_by_dimension(mthis_set, m_GD, fixmat_ents ,true),mb);
[ # # ][ # # ]
[ # # ][ # # ]
333 : : }
334 : :
335 : : // get all the nodes in the material and tag bl nodes
336 [ + - ][ + - ]: 2 : moab::Range mat_nodes, mat_hexes;
337 [ + - ]: 1 : if(m_GD == 3){
338 [ + - ]: 1 : if(m_Conn == 8)
339 [ + - ][ - + ]: 1 : MBERRCHK(mb->get_entities_by_type(mthis_set, moab::MBHEX, mat_hexes, true), mb);
[ # # ][ # # ]
[ # # ][ # # ]
340 [ # # ]: 0 : else if(m_Conn ==4)
341 [ # # ][ # # ]: 1 : MBERRCHK(mb->get_entities_by_type(mthis_set, moab::MBTET, mat_hexes, true), mb);
[ # # ][ # # ]
[ # # ][ # # ]
342 : : }
343 [ # # ]: 0 : else if(m_GD == 2){
344 [ # # ]: 0 : if(m_Conn == 4)
345 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_entities_by_type(mthis_set, moab::MBQUAD, mat_hexes, true), mb);
[ # # ][ # # ]
[ # # ][ # # ]
346 [ # # ]: 0 : else if(m_Conn == 3)
347 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_entities_by_type(mthis_set, moab::MBTRI, mat_hexes, true), mb);
[ # # ][ # # ]
[ # # ][ # # ]
348 : : }
349 : :
350 : : // tag all the mat_hexes with matid
351 [ + - ][ + - ]: 2 : std::vector<int> matID(mat_hexes.size(), mset_id);
352 [ + - ][ + - ]: 1 : MBERRCHK(mb->tag_set_data(MatIDTag, mat_hexes, &matID[0]), mb);
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
353 : : //
354 [ + - ][ - + ]: 1 : MBERRCHK(mb->get_adjacencies(mat_hexes, 0, false, mat_nodes, Interface::UNION), mb);
[ # # ][ # # ]
[ # # ][ # # ]
355 [ + - ]: 2 : moab::Range mat_b_nodes = intersect(nodes, mat_nodes);
356 : :
357 [ + - ][ + - ]: 2 : std::vector<int> bl_node_data(mat_b_nodes.size(), 0);
358 [ + - ][ + - ]: 2 : std::vector<int> node_tag_data(mat_b_nodes.size(),-1);
359 : :
360 : : // don't do error check, as it is supposed to give error when multiple material case is encountered
361 [ + - ][ + - ]: 1 : mb->tag_get_data(MNTag, mat_b_nodes, &node_tag_data[0]);
362 [ + - ][ + + ]: 9 : for(int i=0; i< (int)mat_b_nodes.size(); i++){
363 : : // already a part of some material
364 [ + - ][ - + ]: 8 : if(node_tag_data[i] != -1){
365 [ # # ][ # # ]: 0 : bl_node_data[i] = node_tag_data[i]+1;
366 : : }
367 : : }
368 : :
369 [ + - ][ + - ]: 1 : MBERRCHK(mb->tag_set_data(MNTag, mat_b_nodes, &bl_node_data[0]), mb);
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
370 [ + - ]: 1 : mat_hexes.clear();
371 [ + - ]: 1 : mat_b_nodes.clear();
372 [ + - ]: 1 : mat_nodes.clear();
373 : 1 : mthis_set = 0;
374 : 1 : } // end handling material set
375 : :
376 : : // if fixmat specified, filter old hex, we don't have to correct both sides of the boundary
377 [ - + ][ # # ]: 1 : if (fixmat !=-1 && (int) old_hex.size() > 1){
[ - + ]
378 : : moab::EntityHandle old_hex_set;
379 [ # # ][ # # ]: 0 : MBERRCHK(mb->create_meshset(moab::MESHSET_SET, old_hex_set, 1), mb);
[ # # ][ # # ]
[ # # ][ # # ]
380 [ # # ][ # # ]: 0 : MBERRCHK(mb->add_entities(old_hex_set,&old_hex[0], (int) old_hex.size()), mb);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
381 : : // TODO: Find a faster way of doing this
382 [ # # ][ # # ]: 0 : MBERRCHK(mb->remove_entities(old_hex_set, fixmat_ents), mb);
[ # # ][ # # ]
[ # # ][ # # ]
383 : 0 : old_hex.clear();
384 : 0 : old_hex.empty();
385 : : // the the old hex to be modified
386 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_entities_by_dimension(old_hex_set, m_GD, old_hex), mb);
[ # # ][ # # ]
[ # # ][ # # ]
387 : : }
388 [ + - ][ - + ]: 1 : else if(fixmat ==-1 && (int) old_hex.size()>1){
[ - + ]
389 [ # # ]: 0 : m_LogFile << "FIXMAT not defined, elements found on either side of specified BL surface, aborting...";
390 [ # # ][ # # ]: 0 : m_LogFile << "\n\n Define FIXMAT keyword with material id that remains fixed." << std::endl;
391 : 0 : exit(0);
392 : : }
393 : :
394 [ + - ][ + - ]: 1 : MBERRCHK(mb->get_connectivity(&(*kter), 1, qconn),mb);
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
395 [ + - ][ + - ]: 1 : MBERRCHK(mb->get_connectivity(&old_hex[0], 1, old_hex_conn),mb);
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
396 : :
397 : : // get the normal to the plane of the 2D mesh
398 [ - + ]: 1 : if(m_GD==2)
399 [ # # ][ # # ]: 0 : get_normal_quad (old_hex_conn, surf_normal);
400 [ - + ][ # # ]: 1 : if(found == 1 && m_Material !=999999){
401 [ # # ][ # # ]: 0 : m_LogFile << "Found material set with id " << m_Material << std::endl;
[ # # ]
402 : : }
403 [ + - ][ + - ]: 1 : else if (m_Material !=999999 && found == 0){
404 : : // material set found, but non-existant. Now create this material set
405 [ + - ][ + - ]: 1 : m_LogFile << "Creating material set with id " << m_Material << std::endl;
[ + - ]
406 [ + - ][ - + ]: 1 : MBERRCHK(mb->create_meshset(moab::MESHSET_SET, mthis_set, 1), mb);
[ # # ][ # # ]
[ # # ][ # # ]
407 [ + - ][ - + ]: 1 : MBERRCHK(mb->tag_set_data(MTag, &mthis_set, 1, &m_Material), mb);
[ # # ][ # # ]
[ # # ][ # # ]
408 : : }
409 [ # # ][ # # ]: 0 : else if (m_Material == 999999 && found == 0){
410 [ # # ][ # # ]: 0 : m_LogFile << "Use old materials for new boundary layer elements. No material specified." << std::endl;
411 : :
412 : : // If no material tag exists in the model create one now
413 [ # # ][ # # ]: 0 : if((int)m_sets.size() == 0){
414 [ # # ][ # # ]: 0 : m_LogFile << "Creating material set with id " << m_Material << std::endl;
[ # # ]
415 [ # # ][ # # ]: 0 : MBERRCHK(mb->create_meshset(moab::MESHSET_SET, mthis_set, 1), mb);
[ # # ][ # # ]
[ # # ][ # # ]
416 [ # # ][ # # ]: 0 : MBERRCHK(mb->tag_set_data(MTag, &mthis_set, 1, &m_Material), mb);
[ # # ][ # # ]
[ # # ][ # # ]
417 : 0 : }
418 : : }
419 : : else{
420 [ # # ][ # # ]: 0 : m_LogFile << "Unhandled case" << std::endl;
421 : 0 : exit(0);
422 : : }
423 : :
424 [ + - ]: 1 : err = compute_normals();
425 [ - + ]: 1 : if (err!=0){
426 [ # # ][ # # ]: 0 : m_LogFile << "Failed to compute normals" << std::endl;
427 : 0 : exit(0);
428 : : }
429 : :
430 [ + - ]: 1 : err = push_bulk_mesh(vw);
431 [ - + ]: 1 : if (err!=0){
432 [ # # ][ # # ]: 0 : m_LogFile << "Failed to push bulk mesh" << std::endl;
433 : 0 : exit(0);
434 : : }
435 [ + - ]: 1 : err = create_bl_elements(vw);
436 [ - + ]: 1 : if (err!=0){
437 [ # # ][ # # ]: 0 : m_LogFile << "Failed to create boundary layer elements after pushing bulk mesh" << std::endl;
438 : 0 : exit(0);
439 : : }
440 : :
441 [ + - ][ + - ]: 1 : m_LogFile << "\nTotal Jacobian calls/Min/Max of penultimate hex elements:" << m_JacCalls << ", " << m_JLo << ", " << m_JHi << std::endl;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
442 : :
443 : : // convert the final mesh to hex27 for Nek5000
444 [ - + ]: 1 : if(hex27 == 1){
445 [ # # ]: 0 : moab::Range entities;
446 : : moab::EntityHandle meshset;
447 [ # # ]: 0 : mb->get_entities_by_type(0, MBHEX, entities);
448 [ # # ]: 0 : mb->create_meshset(MESHSET_SET, meshset);
449 [ # # ]: 0 : mb->add_entities(meshset, entities);
450 : : // Add nodes along mid- edge, face and region
451 [ # # ]: 0 : mb->convert_entities(meshset, true, true, true);
452 [ # # ][ # # ]: 0 : m_LogFile << "Mesh converted from hex 8 to hex 27 elements" << std::endl;
453 : : }
454 : :
455 : : // save the final mesh with boundary layer elements
456 [ + - ][ - + ]: 1 : MBERRCHK(mb->write_mesh(m_OutFile.c_str()),mb);
[ # # ][ # # ]
[ # # ][ # # ]
457 [ + - ][ + - ]: 1 : m_LogFile << "\n\nWrote Mesh File: " << m_OutFile << std::endl;
[ + - ][ + - ]
458 : : // get the current date and time
459 [ + - ]: 1 : Timer.GetDateTime (szDateTime);
460 [ + - ][ + - ]: 1 : m_LogFile << "Ending at : " << szDateTime;
[ + - ]
461 : : // report/compute the elapsed time
462 [ + - ][ + - ]: 1 : m_LogFile << "Elapsed wall clock time: " << Timer.DiffTime ()
[ + - ]
463 [ + - ][ + - ]: 2 : << " seconds or " << (Timer.DiffTime ())/60.0 << " mins\n";
[ + - ][ + - ]
464 [ + - ][ + - ]: 1 : m_LogFile << "Total CPU time used: " << (double) (clock() - sTime)/CLOCKS_PER_SEC \
465 [ + - ][ + - ]: 2 : << " seconds" << std::endl;
466 : 1 : }
467 : :
468 : 1 : void PostBL::PrepareIO (int argc, char *argv[], std::string TestDir)
469 : : // ---------------------------------------------------------------------------
470 : : //! Function: Parser for reading the PostBL specification (.inp) file. \n
471 : : //! Input: Command line arguments. \n
472 : : //! Output: none \n
473 : : // ---------------------------------------------------------------------------
474 : : {
475 : : // set and open input output files
476 : 1 : bool bDone = false;
477 [ - + ]: 1 : do{
478 [ - + ]: 1 : if (2 == argc) {
479 [ # # ]: 0 : m_InputFile = argv[1];
480 [ # # ][ # # ]: 0 : m_LogName = m_InputFile + ".log";
481 : : }
482 [ + - ]: 1 : else if (1 == argc){
483 [ + - ][ + - ]: 1 : m_InputFile = TestDir + "/" + (char *)DEFAULT_TEST_POSTBL;
[ + - ]
484 [ + - ][ + - ]: 1 : m_LogName = (std::string)DEFAULT_TEST_POSTBL + ".log";
[ + - ]
485 : : }
486 : :
487 : : // open input file for reading
488 [ + - ]: 1 : m_FileInput.open (m_InputFile.c_str(), std::ios::in);
489 [ + - ][ - + ]: 1 : if (!m_FileInput){
490 [ # # ][ # # ]: 0 : m_LogFile << "Usage: postbl <filename.inp> " << std::endl;
491 [ # # ][ # # ]: 0 : m_LogFile << "Default test file can be found here <Meshkit/data>" << std::endl;
492 [ # # ][ # # ]: 0 : m_LogFile << " Examples input and mesh files are located here <MeshKit/test/algs/postbl_examples>" << std::endl;
493 [ # # ]: 0 : m_FileInput.clear ();
494 : 0 : exit(1);
495 : : }
496 : : else
497 : 1 : bDone = true; // file opened successfully
498 : :
499 : : // open the log file for dumping debug/output statements
500 [ + - ]: 1 : m_LogFile.coss.open (m_LogName.c_str(), std::ios::out);
501 [ + - ][ - + ]: 1 : if (!m_LogFile.coss){
502 [ # # ][ # # ]: 0 : m_LogFile << "Unable to open file: " << m_LogName << std::endl;
[ # # ][ # # ]
503 [ # # ]: 0 : m_LogFile.coss.clear ();
504 : 0 : exit(1);
505 : : }
506 : : else
507 : 1 : bDone = true; // file opened successfully
508 [ + - ]: 1 : m_LogFile << '\n';
509 [ + - ][ + - ]: 1 : m_LogFile << "\t\t---------------------------------------------------------" << '\n';
510 [ + - ][ + - ]: 1 : m_LogFile << "\t\t Tool to generate Post-mesh Boundary Layers " << '\n';
511 [ + - ][ + - ]: 1 : m_LogFile << "\t\t\t\tArgonne National Laboratory" << '\n';
512 [ + - ][ + - ]: 1 : m_LogFile << "\t\t\t\t 2012 " << '\n';
513 [ + - ][ + - ]: 1 : m_LogFile << "\t\t---------------------------------------------------------" << '\n';
514 [ + - ][ + - ]: 1 : m_LogFile << "\nsee README file for using the program and details on various cards.\n"<< std::endl;
515 : :
516 : 1 : }while (!bDone);
517 : :
518 : : // Get the meshfile name, surface(s), thickness, intervals and bias
519 [ + - ]: 1 : CParser Parse;
520 : 1 : int maxloop = 1000;
521 : : // count the total number of cylinder commands in each pincellh
522 : : for(;;){
523 [ - + ]: 9 : if(m_nLineNumber > maxloop){
524 [ # # ][ # # ]: 0 : m_LogFile << "Warning: Didn't find End keyword, stopping after reading 1000 lines from input file" << std::endl;
525 : 0 : break;
526 : : }
527 [ - + ]: 9 : if (!Parse.ReadNextLine (m_FileInput, m_nLineNumber, szInputString,
528 [ + - ]: 9 : MAXCHARS, szComment))
529 [ # # ]: 0 : IOErrorHandler (INVALIDINPUT);
530 : :
531 : : // Get tri scheme
532 [ + - ][ + - ]: 9 : if (szInputString.substr(0,9) == "trischeme"){
[ - + ]
533 [ # # ]: 0 : std::istringstream szFormatString (szInputString);
534 [ # # ][ # # ]: 0 : szFormatString >> m_Card >> tri_sch;
535 [ # # ][ # # ]: 0 : if(szFormatString.fail())
536 [ # # ]: 0 : IOErrorHandler(INVALIDINPUT);
537 [ # # ][ # # ]: 0 : m_LogFile << m_Card << " name read: "<< tri_sch << std::endl;
[ # # ][ # # ]
[ # # ]
538 : : }
539 : : // Get hybrid
540 [ + - ][ + - ]: 9 : if (szInputString.substr(0,6) == "hybrid"){
[ - + ]
541 [ # # ]: 0 : std::istringstream szFormatString (szInputString);
542 [ # # ][ # # ]: 0 : szFormatString >> m_Card >> hybrid;
543 [ # # ][ # # ]: 0 : if(szFormatString.fail())
544 [ # # ]: 0 : IOErrorHandler(INVALIDINPUT);
545 [ # # ][ # # ]: 0 : m_LogFile << m_Card << " name read: "<< hybrid << std::endl;
[ # # ][ # # ]
[ # # ]
546 : : }
547 : : // Get hybrid
548 [ + - ][ + - ]: 9 : if (szInputString.substr(0,6) == "fixmat"){
[ - + ]
549 [ # # ]: 0 : std::istringstream szFormatString (szInputString);
550 [ # # ][ # # ]: 0 : szFormatString >> m_Card >> fixmat;
551 [ # # ][ # # ]: 0 : if(szFormatString.fail())
552 [ # # ]: 0 : IOErrorHandler(INVALIDINPUT);
553 [ # # ][ # # ]: 0 : m_LogFile << m_Card << " name read: "<< fixmat << std::endl;
[ # # ][ # # ]
[ # # ]
554 : : }
555 : : // Get MeshFile name
556 [ + - ][ + - ]: 9 : if (szInputString.substr(0,8) == "meshfile"){
[ + + ]
557 [ + - ]: 1 : std::istringstream szFormatString (szInputString);
558 [ + - ][ + - ]: 1 : szFormatString >> m_Card >> m_MeshFile;
559 [ + - ][ - + ]: 1 : if(szFormatString.fail())
560 [ # # ]: 0 : IOErrorHandler(INVALIDINPUT);
561 [ + - ][ + - ]: 1 : m_LogFile << m_Card << " name read: "<< m_MeshFile << std::endl;
[ + - ][ + - ]
[ + - ][ + - ]
562 [ + - ]: 1 : if (argc == 1){
563 [ + - ][ + - ]: 1 : m_MeshFile = TestDir + "/" + m_MeshFile;
[ + - ]
564 : 1 : }
565 : : }
566 : : // Get BL surface
567 [ + - ][ + - ]: 9 : if (szInputString.substr(0,8) == "surfaces"){
[ - + ]
568 [ # # ]: 0 : std::istringstream szFormatString (szInputString);
569 [ # # ][ # # ]: 0 : szFormatString >> m_Card >> m_SurfId;
570 [ # # ][ # # ]: 0 : if(m_SurfId < 0 || szFormatString.fail())
[ # # ][ # # ]
571 [ # # ]: 0 : IOErrorHandler(INVALIDINPUT);
572 [ # # ][ # # ]: 0 : m_LogFile << m_Card << " read: "<< m_SurfId <<std::endl;
[ # # ][ # # ]
[ # # ]
573 : : }
574 : : // Get BL surface via neumann set or sideset
575 [ + - ][ + - ]: 9 : if (szInputString.substr(0,10) == "neumannset"){
[ + + ]
576 [ + - ]: 1 : std::istringstream szFormatString (szInputString);
577 [ + - ][ + - ]: 1 : szFormatString >> m_Card >> m_NeumannSet;
578 [ + - ][ + - ]: 1 : if(m_NeumannSet < 0 || szFormatString.fail())
[ - + ][ - + ]
579 [ # # ]: 0 : IOErrorHandler(INVALIDINPUT);
580 [ + - ][ + - ]: 1 : m_LogFile << m_Card << " read: "<< m_NeumannSet <<std::endl;
[ + - ][ + - ]
[ + - ]
581 : : }
582 : : // Get BL material (block) number
583 [ + - ][ + - ]: 9 : if (szInputString.substr(0,8) == "material"){
[ + + ]
584 [ + - ]: 1 : std::istringstream szFormatString (szInputString);
585 [ + - ][ + - ]: 1 : szFormatString >> m_Card >> m_Material;
586 [ + - ][ + - ]: 1 : if(m_Material < 0 || szFormatString.fail())
[ - + ][ - + ]
587 [ # # ]: 0 : IOErrorHandler(INVALIDINPUT);
588 [ + - ][ + - ]: 1 : m_LogFile << m_Card << " read: "<< m_Material <<std::endl;
[ + - ][ + - ]
[ + - ]
589 : : }
590 : :
591 : : // Get thickness
592 [ + - ][ + - ]: 9 : if (szInputString.substr(0,9) == "thickness"){
[ + + ]
593 [ + - ]: 1 : std::istringstream szFormatString (szInputString);
594 [ + - ][ + - ]: 1 : szFormatString >> m_Card >> m_Thickness;
595 [ + - ][ + - ]: 1 : if(m_Thickness < 0 || szFormatString.fail())
[ - + ][ - + ]
596 [ # # ]: 0 : IOErrorHandler(INVALIDINPUT);
597 [ + - ][ + - ]: 1 : m_LogFile << m_Card << " read: "<< m_Thickness <<std::endl;
[ + - ][ + - ]
[ + - ]
598 : : }
599 : : // Get intervals
600 [ + - ][ + - ]: 9 : if (szInputString.substr(0,9) == "intervals"){
[ + + ]
601 [ + - ]: 1 : std::istringstream szFormatString (szInputString);
602 [ + - ][ + - ]: 1 : szFormatString >> m_Card >> m_Intervals;
603 [ + - ][ + - ]: 1 : if(m_Intervals < 0 || szFormatString.fail())
[ - + ][ - + ]
604 [ # # ]: 0 : IOErrorHandler(INVALIDINPUT);
605 [ + - ][ + - ]: 1 : m_LogFile << m_Card << " read: "<< m_Intervals <<std::endl;
[ + - ][ + - ]
[ + - ]
606 : : }
607 : : // Get bias
608 [ + - ][ + - ]: 9 : if (szInputString.substr(0,4) == "bias"){
[ + + ]
609 [ + - ]: 1 : std::istringstream szFormatString (szInputString);
610 [ + - ][ + - ]: 1 : szFormatString >> m_Card >> m_Bias;
611 [ + - ][ + - ]: 1 : if(m_Bias < 0 || szFormatString.fail())
[ - + ][ - + ]
612 [ # # ]: 0 : IOErrorHandler(INVALIDINPUT);
613 [ + - ][ + - ]: 1 : m_LogFile << m_Card << " read: "<< m_Bias <<std::endl;
[ + - ][ + - ]
[ + - ]
614 : : }
615 : : // Output file name
616 [ + - ][ + - ]: 9 : if (szInputString.substr(0,7) == "outfile"){
[ + + ]
617 [ + - ]: 1 : std::istringstream szFormatString (szInputString);
618 [ + - ][ + - ]: 1 : szFormatString >> m_Card >> m_OutFile;
619 [ + - ][ - + ]: 1 : if(szFormatString.fail())
620 [ # # ]: 0 : IOErrorHandler(INVALIDINPUT);
621 [ + - ][ + - ]: 1 : m_LogFile << m_Card << " name read: "<< m_OutFile <<std::endl;
[ + - ][ + - ]
[ + - ][ + - ]
622 : : }
623 : : // Output file name
624 [ + - ][ + - ]: 9 : if (szInputString.substr(0,7) == "outfile"){
[ + + ]
625 [ + - ]: 1 : std::istringstream szFormatString (szInputString);
626 [ + - ][ + - ]: 1 : szFormatString >> m_Card >> m_OutFile;
627 [ + - ][ - + ]: 1 : if(szFormatString.fail())
628 [ # # ]: 0 : IOErrorHandler(INVALIDINPUT);
629 [ + - ][ + - ]: 1 : m_LogFile << m_Card << " name read: "<< m_OutFile <<std::endl;
[ + - ][ + - ]
[ + - ][ + - ]
630 : : }
631 : : // save hex27 mesh
632 [ + - ][ + - ]: 9 : if (szInputString.substr(0,5) == "hex27"){
[ - + ]
633 [ # # ]: 0 : std::istringstream szFormatString (szInputString);
634 [ # # ][ # # ]: 0 : szFormatString >> m_Card >> hex27;
635 [ # # ][ # # ]: 0 : if(szFormatString.fail())
636 [ # # ]: 0 : IOErrorHandler(INVALIDINPUT);
637 [ # # ][ # # ]: 0 : m_LogFile << m_Card << " hex27 creation: "<< hex27 <<std::endl;
[ # # ][ # # ]
[ # # ]
638 : : }
639 [ + - ][ + - ]: 9 : if (szInputString.substr(0,3) == "end"){
[ + + ]
640 : 1 : break;
641 : : }
642 : 9 : }
643 : 1 : }
644 : :
645 : 0 : void PostBL::IOErrorHandler (ErrorStates ECode) const
646 : : // ---------------------------------------------------------------------------
647 : : //! Function: Displays error messages related to input data \n
648 : : //! Input: Error code \n
649 : : //! Output: none \n
650 : : // ---------------------------------------------------------------------------
651 : : {
652 : 0 : std::cerr << '\n';
653 [ # # ]: 0 : if (ECode == INVALIDINPUT) // invalid input
654 : 0 : std::cerr << "Invalid input.";
655 : : else
656 : 0 : std::cerr << "Unknown error ...?";
657 : :
658 : 0 : std::cerr << '\n' << "Error in input file line : " << m_nLineNumber;
659 : 0 : std::cerr << std::endl;
660 : 0 : exit (1);
661 : : }
662 : :
663 : :
664 [ + - ][ + - ]: 156 : } // namespace MeshKit
665 : :
666 : :
|
