|
MeshKit
1.0
|
|
MeshKit
1.0
|
00001 /********************************************* 00002 July 31,2013 00003 PostBL: Post-Mesh BOundary Layer Improved Algorithms 2 00004 Argonne National Laboratory 00005 *********************************************/ 00006 00007 #include "meshkit/PostBL.hpp" 00008 using namespace MeshKit; 00009 #include <string.h> 00010 00011 void PostBL:: Algo2(){ 00012 00013 m_LogFile << "\nIn execute this : creating boundary layer elements.." << std::endl; 00014 // start the timer 00015 CClock Timer; 00016 clock_t sTime = clock(); 00017 std::string szDateTime; 00018 Timer.GetDateTime (szDateTime); 00019 00020 m_LogFile << "\nStarting out at : " << szDateTime << std::endl; 00021 m_LogFile << "\n Loading meshfile: " << m_MeshFile << ".." << std::endl; 00022 00023 // load specified mesh file 00024 IBERRCHK(imesh->load(0, m_MeshFile.c_str(),0), *imesh); 00025 moab::Range all_elems, all_verts; 00026 MBERRCHK(mb->get_entities_by_dimension(0, 3, all_elems,true),mb); 00027 if (all_elems.size() == 0) 00028 m_GD = 2; 00029 else if (all_elems.size() > 0) 00030 m_GD = 3; 00031 else 00032 exit(0); 00033 all_elems.clear(); 00034 m_LogFile << "Geometric dimension of meshfile = "<< m_GD <<std::endl; 00035 00036 // obtain existing tag handles 00037 moab::Tag GDTag, GIDTag, NTag, MTag, STag, FTag, MNTag, MatIDTag, BLNodeIDTag; 00038 MBERRCHK(mb->tag_get_handle("GEOM_DIMENSION", 1, moab::MB_TYPE_INTEGER, GDTag),mb); 00039 MBERRCHK(mb->tag_get_handle("NEUMANN_SET", 1, moab::MB_TYPE_INTEGER, NTag),mb); 00040 MBERRCHK(mb->tag_get_handle("MATERIAL_SET", 1, moab::MB_TYPE_INTEGER, MTag),mb); 00041 MBERRCHK(mb->tag_get_handle("GLOBAL_ID", 1, moab::MB_TYPE_INTEGER, GIDTag),mb); 00042 // create smoothset and fixed tag for mesquite 00043 MBERRCHK(mb->tag_get_handle("SMOOTHSET", 1, moab::MB_TYPE_INTEGER, STag, 00044 moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT),mb); 00045 MBERRCHK(mb->tag_get_handle("fixed", 1, moab::MB_TYPE_INTEGER, FTag, 00046 moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT),mb); 00047 MBERRCHK(mb->tag_get_handle("mnode", 1, moab::MB_TYPE_INTEGER, MNTag, 00048 moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT),mb); 00049 MBERRCHK(mb->tag_get_handle("matid", 1, moab::MB_TYPE_INTEGER, MatIDTag, 00050 moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT),mb); 00051 MBERRCHK(mb->tag_get_handle("BLNODEID", 1, moab::MB_TYPE_INTEGER, BLNodeIDTag, 00052 moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT),mb); 00053 00054 // get all the entity sets with boundary layer geom dimension, neumann sets and material sets 00055 moab::Range sets, n_sets, m_sets; 00056 m_BLDim = m_GD - 1; 00057 00058 const void* gdim[] = {&m_BLDim}; 00059 00060 MBERRCHK(mb->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &GDTag, 00061 gdim, 1 , sets, moab::Interface::INTERSECT, false), mb); 00062 MBERRCHK(mb->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &NTag, 0, 1 , n_sets),mb); 00063 MBERRCHK(mb->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &MTag, 0, 1 , m_sets),mb); 00064 00065 // Handling NeumannSets (if BL surf in input via NS) 00066 moab::Range::iterator set_it; 00067 moab::EntityHandle this_set = 0; 00068 for (set_it = n_sets.begin(); set_it != n_sets.end(); set_it++) { 00069 00070 this_set = *set_it; 00071 00072 // get entity handle of NS specified in the input file 00073 int set_id; 00074 MBERRCHK(mb->tag_get_data(NTag, &this_set, 1, &set_id), mb); 00075 if(set_id == m_NeumannSet) 00076 break; 00077 this_set = 0; 00078 } 00079 if (debug && m_NeumannSet != -1 && this_set != 0){ 00080 m_LogFile << "Looking for NS with id " << m_NeumannSet << 00081 ". Total NS found are: "<< n_sets.size() << std::endl; 00082 } 00083 00084 // For specified surface: get the all the quads and nodes in a range 00085 moab::EntityHandle s1; 00086 moab::Range quads, nodes,edges, fixmat_ents; 00087 int dims; // variable to store global id of boundary layer specified in the input file 00088 00089 // Method 1: INPUT by NeumannSet 00090 if(m_NeumannSet != -1 && this_set != 0){ 00091 MBERRCHK(mb->get_entities_by_dimension(this_set, m_BLDim, quads,true),mb); 00092 if (quads.size() <=0){ 00093 m_LogFile << " No quads found, aborting.. " << std::endl; 00094 exit(0); 00095 } 00096 00097 MBERRCHK(mb->get_adjacencies(quads, 0, false, nodes, moab::Interface::UNION),mb); 00098 if(m_GD == 3) 00099 MBERRCHK(mb->get_adjacencies(quads, 1, true, edges, moab::Interface::UNION),mb); 00100 00101 if (debug) { 00102 m_LogFile << "Found NeumannSet with id : " << m_NeumannSet << std::endl; 00103 m_LogFile << "#Quads in this surface: " << quads.size() << std::endl; 00104 m_LogFile << "#Nodes in this surface: " << nodes.size() << std::endl; 00105 m_LogFile << "#New nodes to be created:" << m_Intervals*nodes.size() << std::endl; 00106 } 00107 } 00108 // Method 2: INPUT by surface id (geom dimension) 00109 else if (m_SurfId !=-1){ 00110 for(moab::Range::iterator rit=sets.begin(); rit != sets.end(); ++rit){ 00111 s1 = *rit; 00112 MBERRCHK(mb->tag_get_data(GIDTag, &s1, 1, &dims),mb); 00113 00114 if(dims == m_SurfId && m_SurfId != -1){ 00115 MBERRCHK(mb->get_entities_by_dimension(s1, m_BLDim, quads,true),mb); 00116 if (quads.size() <=0){ 00117 m_LogFile << " No quads found, aborting.. " << std::endl; 00118 exit(0); 00119 } 00120 00121 MBERRCHK(mb->get_adjacencies(quads, 0, false, nodes, moab::Interface::UNION),mb); 00122 if(m_GD == 3) 00123 MBERRCHK(mb->get_adjacencies(edges, 0, false, nodes, moab::Interface::UNION),mb); 00124 if (debug) { 00125 m_LogFile << "Found surface with id : " << m_SurfId << std::endl; 00126 m_LogFile << "#Quads in this surface: " << quads.size() << std::endl; 00127 m_LogFile << "#Nodes in this surface: " << nodes.size() << std::endl; 00128 m_LogFile << "#New nodes to be created:" << m_Intervals*nodes.size() << std::endl; 00129 } 00130 } 00131 } 00132 } 00133 00134 if (quads.size() == 0 || nodes.size() == 0) { 00135 m_LogFile << "Invalid boundary layer specification, aborting.." << std::endl; 00136 exit(0); 00137 } 00138 00139 // placeholder for storing smoothing entities 00140 moab::EntityHandle smooth_set; 00141 int s_id = 100; 00142 MBERRCHK(mb->create_meshset(moab::MESHSET_SET, smooth_set, 1), mb); 00143 MBERRCHK(mb->tag_set_data(STag, &smooth_set, 1, &s_id), mb); 00144 00145 // placeholder for storing gd on new entities 00146 moab::EntityHandle geom_set; 00147 MBERRCHK(mb->create_meshset(moab::MESHSET_SET, geom_set, 1), mb); 00148 MBERRCHK(mb->tag_set_data(GDTag, &geom_set, 1, &m_GD), mb); 00149 00150 // declare variables before starting BL creation 00151 std::vector <bool> node_status(false); // size of verts of bl surface 00152 node_status.resize(nodes.size()); 00153 moab::Range hexes, hex_edge, quad_verts; 00154 double coords_new_quad[3]; 00155 moab::EntityHandle hex, hex1; 00156 int qcount = 0; 00157 00158 //size of the following is based on element type 00159 std::vector<moab::EntityHandle> conn, qconn, adj_qconn, tri_conn, 00160 new_vert(m_Intervals*nodes.size()), old_hex_conn, adj_hexes, adj_quads, adj_hex_nodes1; 00161 moab::CartVect surf_normal; 00162 00163 // element on the boundary 00164 Range::iterator kter = quads.begin(); 00165 00166 std::vector<moab::EntityHandle> old_hex; 00167 MBERRCHK(mb->get_adjacencies(&(*kter), 1, m_GD, false, old_hex),mb); 00168 if((int) old_hex.size() == 0){ 00169 m_LogFile << "unable to find adjacent hex for BL quad, aborting..."; 00170 exit(0); 00171 } 00172 00173 // allocate space for connectivity/adjacency during the first pass of this loop 00174 if(mb->type_from_handle(old_hex[0]) == MBHEX){ 00175 m_Conn = 8; 00176 m_BElemNodes = 4; 00177 m_HConn = 8; 00178 //allocating based on element type 00179 conn.resize(m_Intervals*m_Conn), qconn.resize(m_BElemNodes), adj_qconn.resize(m_BElemNodes), 00180 old_hex_conn.resize(m_Conn), adj_hex_nodes1.resize(m_Conn); 00181 } 00182 else if(mb->type_from_handle(old_hex[0]) == MBTET){ 00183 m_Conn = 4; 00184 m_HConn = 6; 00185 m_BElemNodes = 3; 00186 //allocating based on element type - thrice the number of elements 00187 if(hybrid) 00188 conn.resize(m_Intervals*6); 00189 else 00190 conn.resize(2*m_Intervals*m_Conn); 00191 qconn.resize(m_BElemNodes), adj_qconn.resize(m_BElemNodes), 00192 old_hex_conn.resize(m_Conn), adj_hex_nodes1.resize(m_Conn); 00193 } 00194 else if(mb->type_from_handle(old_hex[0]) == MBQUAD){ 00195 m_Conn = 4; 00196 m_HConn = 4; 00197 m_BElemNodes = 2; 00198 //allocating based on element type 00199 conn.resize(m_Intervals*m_Conn), qconn.resize(m_BElemNodes), adj_qconn.resize(m_BElemNodes), 00200 old_hex_conn.resize(m_Conn), adj_hex_nodes1.resize(m_Conn); 00201 } 00202 else if(mb->type_from_handle(old_hex[0]) == MBTRI){ 00203 m_Conn = 3; 00204 m_HConn = 4; 00205 m_BElemNodes = 2; 00206 //allocating based on element type - twice the number of elements 00207 if(hybrid){ 00208 m_HConn = 4; 00209 conn.resize(m_Intervals*m_HConn); 00210 } 00211 else{ 00212 tri_conn.resize(2*m_Intervals*m_Conn); 00213 conn.resize(2*m_Intervals*m_Conn); 00214 } 00215 qconn.resize(m_BElemNodes), adj_qconn.resize(m_BElemNodes), 00216 old_hex_conn.resize(m_Conn), adj_hex_nodes1.resize(m_Conn); 00217 } 00218 else if(m_Conn == 0 || m_BElemNodes == 0){ 00219 m_LogFile << "This mesh type is not supported by this tool" << std::endl; 00220 exit(0); 00221 } 00222 00223 // Tag all nodes on outer boundary with a unique number 00224 int node_id = 0; 00225 std::vector<int> NId(nodes.size()); 00226 for(moab::Range::iterator nodes_iter = nodes.begin(); nodes_iter != nodes.end(); nodes_iter++){ 00227 NId[node_id] = node_id; 00228 MBERRCHK(mb->tag_set_data(BLNodeIDTag, &(*nodes_iter),1, &NId[node_id]), mb); 00229 ++node_id; 00230 } 00231 00232 // Handling MaterialSet 00233 moab::Range::iterator mset_it; 00234 moab::EntityHandle mthis_set; 00235 int mset_id = 0, found = 0; 00236 for (mset_it = m_sets.begin(); mset_it != m_sets.end(); mset_it++) { 00237 00238 mthis_set = *mset_it; 00239 00240 // get entity handle of MS specified in the input file 00241 MBERRCHK(mb->tag_get_data(MTag, &mthis_set, 1, &mset_id), mb); 00242 00243 // if no material set is specified, we'll have to resolve 00244 00245 // else just set all the MNTag to 0 00246 00247 if(mset_id == m_Material){ 00248 found = 1; 00249 break; 00250 } 00251 else if(mset_id == fixmat){ 00252 MBERRCHK(mb->get_entities_by_dimension(mthis_set, m_GD, fixmat_ents ,true),mb); 00253 } 00254 00255 // get all the nodes in the material and tag bl nodes 00256 moab::Range mat_nodes, mat_hexes; 00257 if(m_GD == 3){ 00258 if(m_Conn == 8) 00259 MBERRCHK(mb->get_entities_by_type(mthis_set, moab::MBHEX, mat_hexes, true), mb); 00260 else if(m_Conn ==4) 00261 MBERRCHK(mb->get_entities_by_type(mthis_set, moab::MBTET, mat_hexes, true), mb); 00262 } 00263 else if(m_GD == 2){ 00264 if(m_Conn == 4) 00265 MBERRCHK(mb->get_entities_by_type(mthis_set, moab::MBQUAD, mat_hexes, true), mb); 00266 else if(m_Conn == 3) 00267 MBERRCHK(mb->get_entities_by_type(mthis_set, moab::MBTRI, mat_hexes, true), mb); 00268 } 00269 // tag all the mat_hexes with matid 00270 std::vector<int> matID(mat_hexes.size(), mset_id); 00271 MBERRCHK(mb->tag_set_data(MatIDTag, mat_hexes, &matID[0]), mb); 00272 // 00273 MBERRCHK(mb->get_adjacencies(mat_hexes, 0, false, mat_nodes, Interface::UNION), mb); 00274 moab::Range mat_b_nodes = intersect(nodes, mat_nodes); 00275 00276 // 00277 std::vector<int> bl_node_data(mat_b_nodes.size(), 0); 00278 std::vector<int> node_tag_data(mat_b_nodes.size(),-1); 00279 // don't error check, as it is supposed to give error when multiple material case is encountered 00280 mb->tag_get_data(MNTag, mat_b_nodes, &node_tag_data[0]); 00281 for(int i=0; i< (int)mat_b_nodes.size(); i++){ 00282 // already a part of some material 00283 if(node_tag_data[i] != -1){ 00284 bl_node_data[i] = node_tag_data[i]+1; 00285 } 00286 } 00287 MBERRCHK(mb->tag_set_data(MNTag, mat_b_nodes, &bl_node_data[0]), mb); 00288 mat_hexes.clear(); 00289 mat_b_nodes.clear(); 00290 mat_nodes.clear(); 00291 mthis_set = 0; 00292 } // end handling material set 00293 00294 // if fixmat specified, filter old hex, we don't have to correct both sides of the boundary 00295 if (fixmat !=0 && (int) old_hex.size() > 1){ 00296 moab::EntityHandle old_hex_set; 00297 MBERRCHK(mb->create_meshset(moab::MESHSET_SET, old_hex_set, 1), mb); 00298 MBERRCHK(mb->add_entities(old_hex_set,&old_hex[0], (int) old_hex.size()), mb); 00299 // TODO: Find a faster way of doing this 00300 MBERRCHK(mb->remove_entities(old_hex_set, fixmat_ents), mb); 00301 old_hex.clear(); 00302 old_hex.empty(); 00303 // the the old hex to be modified 00304 MBERRCHK(mb->get_entities_by_dimension(old_hex_set, m_GD, old_hex), mb); 00305 } 00306 else if(fixmat ==0 && old_hex.size()>1){ 00307 m_LogFile << "FIXMAT not defined, elements found on either side of specified BL surface, aborting..."; 00308 exit(0); 00309 } 00310 00311 MBERRCHK(mb->get_connectivity(&(*kter), 1, qconn),mb); 00312 MBERRCHK(mb->get_connectivity(&old_hex[0], 1, old_hex_conn),mb); 00313 // get the normal to the plane of the 2D mesh 00314 if(m_GD==2) 00315 get_normal_quad (old_hex_conn, surf_normal); 00316 if(found == 1 && m_Material !=999){ 00317 m_LogFile << "Found material set with id " << m_Material << std::endl; 00318 } 00319 else if (m_Material !=999 && found == 0){ 00320 // material set found, but non-existant. Now create this material set 00321 m_LogFile << "Creating material set with id " << m_Material << std::endl; 00322 MBERRCHK(mb->create_meshset(moab::MESHSET_SET, mthis_set, 1), mb); 00323 MBERRCHK(mb->tag_set_data(MTag, &mthis_set, 1, &m_Material), mb); 00324 } 00325 else if (m_Material == 999 && found == 0){ 00326 // 00327 m_LogFile << "Use old materials for new boundary layer elements. No material specified." << std::endl; 00328 } 00329 else{ 00330 m_LogFile << "Unhandled case" << std::endl; 00331 exit(0); 00332 } 00333 00334 // COMPUTE NORMALS 00335 // get tag data and print 00336 std::vector<int> all_bl(nodes.size(), 0); 00337 mb->tag_get_data(MNTag, nodes, &all_bl[0]); 00338 int count = -1; 00339 for (Range::iterator kter = nodes.begin(); kter != nodes.end(); ++kter){ 00340 ++count; 00341 // only one normal in case of single material 00342 double xdisp = 0.0, ydisp = 0.0, zdisp = 0.0; 00343 // check if node belongs to one or more materials 00344 if(all_bl[count] == 0){ 00345 moab::Range adj_for_normal; 00346 MBERRCHK(mb->get_adjacencies(&(*kter), 1, m_BLDim, false, adj_for_normal, Interface::UNION), mb); 00347 //create the coordinates of the innermost node corresponding to this node 00348 moab::Range adj_for_norm = intersect(quads, adj_for_normal); 00349 // now compute the average normal direction for this vertex 00350 moab::CartVect rt(0.0, 0.0, 0.0), v(0.0, 0.0, 0.0); 00351 00352 for(Range::iterator qter = adj_for_norm.begin(); qter != adj_for_norm.end(); ++qter){ 00353 MBERRCHK(mb->get_connectivity(&(*qter), 1, adj_qconn),mb); 00354 00355 int side_number = 0, sense = 1, offset = 0; 00356 MBERRCHK(mb->side_number(old_hex[0], (*qter), side_number, sense, offset), mb); 00357 00358 if(m_GD==3){ 00359 get_normal_quad (adj_qconn, v); 00360 if(sense == 1){ 00361 // do nothing 00362 } 00363 else{ 00364 v=-v; 00365 } 00366 } 00367 else if(m_GD==2){ 00368 if(sense == 1){ 00369 get_normal_edge(adj_qconn, surf_normal, v); 00370 } 00371 else{ 00372 get_normal_edge(adj_qconn, -surf_normal, v); 00373 } 00374 } 00375 rt = rt + v; 00376 } 00377 if(rt.length() !=0){ 00378 xdisp=rt[0]/rt.length(); 00379 ydisp=rt[1]/rt.length(); 00380 zdisp=rt[2]/rt.length(); 00381 } 00382 else{ 00383 xdisp=0.0; 00384 ydisp=0.0; 00385 zdisp=0.0; 00386 } 00387 } 00388 else if(all_bl[count] > 0 && fixmat != -1){ // node belongs to more than one material and fixmat specified 00389 // NODE B/W MATERIALS 00390 moab::Range adj_for_normal; 00391 MBERRCHK(mb->get_adjacencies(&(*kter), 1, m_BLDim, false, adj_for_normal, Interface::UNION), mb); 00392 //create the coordinates of the innermost node corresponding to this node 00393 moab::Range adj_for_norm = intersect(quads, adj_for_normal); 00394 // now compute the average normal direction for this vertex 00395 moab::CartVect rt(0.0, 0.0, 0.0), v(0.0, 0.0, 0.0); 00396 00397 for(Range::iterator qter = adj_for_norm.begin(); qter != adj_for_norm.end(); ++qter){ 00398 MBERRCHK(mb->get_connectivity(&(*qter), 1, adj_qconn),mb); 00399 00400 moab::Range this_quad_hex; 00401 MBERRCHK(mb->get_adjacencies(&(*qter), 1, m_GD, false, this_quad_hex, moab::Interface::UNION),mb); 00402 moab::Range quad_hex = intersect(fixmat_ents, this_quad_hex); 00403 00404 int side_number = 0, sense = 1, offset = 0; 00405 Range::iterator hexter = quad_hex.begin(); 00406 MBERRCHK(mb->side_number(*hexter, (*qter), side_number, sense, offset), mb); 00407 00408 if(m_GD==3){ 00409 get_normal_quad (adj_qconn, v); 00410 if(sense == 1 ){ 00411 v=-v; 00412 } 00413 else{ 00414 // do nothing 00415 } 00416 } 00417 else if(m_GD==2){ 00418 if(sense == 1 ){ 00419 get_normal_edge(adj_qconn, -surf_normal, v); 00420 } 00421 else{ 00422 get_normal_edge(adj_qconn, surf_normal, v); 00423 } 00424 } 00425 rt = rt + v; 00426 } 00427 if(rt.length() !=0){ 00428 xdisp=rt[0]/rt.length(); 00429 ydisp=rt[1]/rt.length(); 00430 zdisp=rt[2]/rt.length(); 00431 } 00432 else{ 00433 xdisp=0.0; 00434 ydisp=0.0; 00435 zdisp=0.0; 00436 } 00437 } 00438 else if(all_bl[count] > 0 && fixmat == -1){ // node belongs to more than one material and fixmat not specified 00439 // NODE ON BOUNDARY 00440 // get the edge that is not on the boundary and count how many such edges we have 00441 moab::Range adj_for_normal; 00442 int nEdgeDim = 1; 00443 MBERRCHK(mb->get_adjacencies(&(*kter), 1, nEdgeDim, true, adj_for_normal, Interface::UNION), mb); 00444 moab::Range edge_normal; 00445 if(m_GD == 2) 00446 edge_normal = subtract(adj_for_normal, quads); 00447 else if(m_GD == 3) 00448 edge_normal = subtract(adj_for_normal, edges); 00449 00450 00451 if(edge_normal.size() > 1){ 00452 double ncoord[3]; 00453 MBERRCHK(mb->get_coords(&(*kter), 1, ncoord),mb); 00454 00455 m_LogFile << "Work in progress !! :- MULTIPLE NORMALS ARE NEEDED AT" << ncoord[0] 00456 << ", " << ncoord[1] << ", " << ncoord[2] << " #normals " << edge_normal.size() << std::endl; 00457 exit(0); 00458 } 00459 else{ 00460 m_LogFile << "We've one edge seperating materials 1 NORMAL IS NEEDED" << edge_normal.size() << std::endl; 00461 moab::Range edge_conn; 00462 MBERRCHK(mb->get_connectivity(&(*edge_normal.begin()), 1, edge_conn),mb); 00463 // now get the normal direction for this edge 00464 moab::CartVect coords[2], bl_coord[1]; 00465 moab::CartVect AB; 00466 MBERRCHK(mb->get_coords(&(*kter), 1, (double*) &bl_coord[0]), mb); 00467 MBERRCHK(mb->get_coords(edge_conn, (double*) &coords[0]), mb); 00468 for(int d = 0; d<2; d++){ 00469 if(bl_coord[0][0] == coords[d][0] && 00470 bl_coord[0][1] == coords[d][1] && 00471 bl_coord[0][2] == coords[d][2]){ 00472 // do nothing 00473 } 00474 else{ 00475 AB = (bl_coord[0] - coords[d]); 00476 } 00477 } 00478 00479 xdisp=AB[0]/AB.length(); 00480 ydisp=AB[1]/AB.length(); 00481 zdisp=AB[2]/AB.length(); 00482 } 00483 } 00484 else if(all_bl[count] < 0 ){ 00485 m_LogFile << "Material must have associated with BLNode: Error, shouldn't have gotten here: " << count << std::endl; 00486 } 00487 00488 00489 // after the normal computation is done create new BL nodes 00490 double coords_bl_quad[3], move = 0.0; 00491 MBERRCHK(mb->get_coords(&(*kter), 1, coords_bl_quad),mb); 00492 00493 double temp = 0; 00494 double num = m_Thickness*(m_Bias-1)*(pow(m_Bias, m_Intervals -1)); 00495 double deno = pow(m_Bias, m_Intervals) - 1; 00496 if (deno !=0) 00497 temp = num/deno; 00498 else 00499 temp = m_Thickness/m_Intervals; 00500 00501 // loop thru intervals to create BL nodes 00502 for(int j=0; j< m_Intervals; j++){ 00503 00504 move+= temp/pow(m_Bias,j); 00505 if (debug){ 00506 m_LogFile << " move:" << move; 00507 } 00508 // now compute the coords of the new vertex 00509 coords_new_quad[0] = coords_bl_quad[0]-move*xdisp; 00510 coords_new_quad[1] = coords_bl_quad[1]-move*ydisp; 00511 coords_new_quad[2] = coords_bl_quad[2]-move*zdisp; 00512 00513 int nid = count*m_Intervals+j; 00514 // Possible TODO's 00515 //TODO: Check if this vertex is possible (detect possible collision with geometry) 00516 // TODO: See posibility of using ray tracing 00517 // TODO: Parallize: Profile T-junction model and try to device an algorithm 00518 // TODO: Modularize node creation part and use doxygen for all code and design of code, python design and test cases - current functions in code: 00519 // Setup this, Execute this -- break info sub functions and classes, 00520 // prepareIO --make this optional when using python, 00521 // get normal (2d and 3d) -- can be combined to one function 00522 // get det jacobian (hex elements) --needs check for other elements 00523 // 00524 MBERRCHK(mb->create_vertex(coords_new_quad, new_vert[nid]), mb); 00525 if (debug){ 00526 m_LogFile << std::setprecision (3) << std::scientific << " : created node:" << (nid + 1) 00527 << " of " << new_vert.size() << " new nodes:- coords: " << coords_new_quad[0] 00528 << ", " << coords_new_quad[1] << ", " << coords_new_quad[2] << std::endl; 00529 } 00530 } 00531 } 00532 00533 // swap nodal coordinates of input nodes with innermost BL nodes to push the bulk mesh 00534 count = -1; 00535 for (Range::iterator kter = nodes.begin(); kter != nodes.end(); ++kter){ 00536 ++count; 00537 if(all_bl[count] == 0 ){ 00538 double coords_new_quad[3]; 00539 double coords_old_quad[3]; 00540 00541 int nid = (count+1)*m_Intervals - 1; 00542 MBERRCHK(mb->get_coords(&new_vert[nid], 1, coords_new_quad),mb); 00543 00544 MBERRCHK(mb->get_coords(&(*kter), 1, coords_old_quad),mb); 00545 00546 MBERRCHK(mb->set_coords(&(*kter), 1, coords_new_quad),mb); 00547 00548 MBERRCHK(mb->set_coords(&new_vert[nid], 1, coords_old_quad),mb); 00549 if (debug) { 00550 m_LogFile << std::setprecision (3) << std::scientific << " : NID:" << (nid) 00551 << coords_old_quad[0] 00552 << ", " << coords_old_quad[1] << ", " << coords_old_quad[2] << " OLD:- coords: NEW" << coords_new_quad[0] 00553 << ", " << coords_new_quad[1] << ", " << coords_new_quad[2] << std::endl; 00554 } 00555 } 00556 else if(all_bl[count] > 0 && fixmat != -1){ // node belongs to more than one material and fixmat specified 00557 // NODE B/W MATERIALS 00558 double coords_new_quad[3]; 00559 double coords_old_quad[3]; 00560 00561 int nid = (count+1)*m_Intervals - 1; 00562 MBERRCHK(mb->get_coords(&new_vert[nid], 1, coords_new_quad),mb); 00563 00564 MBERRCHK(mb->get_coords(&(*kter), 1, coords_old_quad),mb); 00565 00566 MBERRCHK(mb->set_coords(&(*kter), 1, coords_new_quad),mb); 00567 00568 MBERRCHK(mb->set_coords(&new_vert[nid], 1, coords_old_quad),mb); 00569 00570 m_LogFile << std::setprecision (3) << std::scientific << "FM : NID:" << (nid) 00571 << coords_old_quad[0] 00572 << ", " << coords_old_quad[1] << ", " << coords_old_quad[2] << " OLD:- coords: NEW" << coords_new_quad[0] 00573 << ", " << coords_new_quad[1] << ", " << coords_new_quad[2] << std::endl; 00574 00575 // now find the hex in fixmat_ents that was affected by this swap and UNswap the coords 00576 moab::Range fhex; 00577 MBERRCHK(mb->get_adjacencies(&(*kter), 1, m_GD, false, fhex, Interface::UNION), mb); 00578 moab::Range fmhex= intersect(fhex,fixmat_ents); 00579 std::vector<moab::EntityHandle> fmconn; 00580 for(Range::iterator fmter = fmhex.begin(); fmter != fmhex.end(); ++fmter){ 00581 MBERRCHK(mb->get_connectivity(&(*fmter), 1, fmconn),mb); 00582 for(int i=0; i < (int)fmconn.size(); i++){ 00583 if((*kter) == fmconn[i]){ 00584 //we have a node to unswap or reset connectivity for fixmat 00585 fmconn[i] = new_vert[nid]; 00586 } 00587 } 00588 MBERRCHK(mb->set_connectivity(*fmter, &fmconn[0], fmconn.size()), mb); 00589 // MBERRCHK(mb->add_adjacencies((*fmter), hexes, true), mb); 00590 // MBERRCHK(mb->add_adjacencies((*fmter), fixmat_ents, true), mb); 00591 00592 } 00593 00594 m_LogFile << " We're here in MM case, now go along the edge --- have fun in the process !!" << std::endl; 00595 // material boundary - get edge direction and length 00596 00597 // find_min_edge_length(adj_qconn_r, qconn[i], nodes, m_MinEdgeLength); 00598 00599 // check to see if this is a fixmat case 00600 } 00601 else if(all_bl[count] > 0 && fixmat == -1){ // node belongs to more than one material and fixmat not specified 00602 // NODE ON BOUNDARY 00603 double coords_new_quad[3]; 00604 double coords_old_quad[3]; 00605 00606 int nid = (count+1)*m_Intervals - 1; 00607 MBERRCHK(mb->get_coords(&new_vert[nid], 1, coords_new_quad),mb); 00608 00609 MBERRCHK(mb->get_coords(&(*kter), 1, coords_old_quad),mb); 00610 00611 MBERRCHK(mb->set_coords(&(*kter), 1, coords_new_quad),mb); 00612 00613 MBERRCHK(mb->set_coords(&new_vert[nid], 1, coords_old_quad),mb); 00614 00615 m_LogFile << std::setprecision (3) << std::scientific << " : NID:" << (nid) 00616 << coords_old_quad[0] 00617 << ", " << coords_old_quad[1] << ", " << coords_old_quad[2] << " OLD:- coords: NEW" << coords_new_quad[0] 00618 << ", " << coords_new_quad[1] << ", " << coords_new_quad[2] << std::endl; 00619 } 00620 } 00621 00622 // check for the volume of penultimate elements if -ve volume encountered. Report. 00623 count = -1; 00624 // Try to move another layer attached to it. 00625 for (Range::iterator kter = nodes.begin(); kter != nodes.end(); ++kter){ 00626 ++count; 00627 for(int j=0; j< m_Intervals; j++){ 00628 int nid = count*m_Intervals+j; 00629 double coords_new_quad[3]; 00630 MBERRCHK(mb->get_coords(&new_vert[nid], 1, coords_new_quad),mb); 00631 m_LogFile << std::setprecision (3) << std::scientific << " : NID:" << (nid) 00632 << " of " << new_vert.size() << " new nodes:- coords: " << coords_new_quad[0] 00633 << ", " << coords_new_quad[1] << ", " << coords_new_quad[2] << std::endl; 00634 } 00635 } 00636 // shoot multiple normals for multiple materials case only. 00637 00638 // This can be used of regular case also how to invoke it. mention in algorithm 00639 00640 // mention local and global smoothing. 00641 00642 // qcount = -1; 00643 // int flag[quads.size()]; 00644 00645 // Now start creating New elements 00646 for (Range::iterator kter = quads.begin(); kter != quads.end(); ++kter){ 00647 qcount++; 00648 // if(flag[qcount] != 1){ 00649 00650 MBERRCHK(mb->get_connectivity(&(*kter), 1, qconn),mb); 00651 double one_node_in_quad[3]; 00652 for (int i=0; i<m_BElemNodes; i++){ 00653 00654 int node_tag_id = 0; 00655 MBERRCHK(mb->tag_get_data(BLNodeIDTag, &qconn[i], 1, &node_tag_id) ,mb); 00656 MBERRCHK(mb->get_coords(&qconn[i], 1, one_node_in_quad),mb); 00657 m_LogFile << std::setprecision (3) << std::scientific << " new nodes:- coords: " << one_node_in_quad[0] 00658 << ", " << one_node_in_quad[1] << ", " << one_node_in_quad[2] << std::endl; 00659 00660 //populate the connectivity after creating nodes for this BL node 00661 for(int j=0; j< m_Intervals; j++){ 00662 if(m_Conn == 8 && m_BElemNodes == 4){ // hex 00663 int nid = node_tag_id*m_Intervals + j; 00664 if(m_Intervals == 1){ 00665 conn[m_Conn*j +i] = qconn[i]; 00666 conn[m_Conn*j + i+m_BElemNodes] = new_vert[nid]; 00667 } 00668 else if(j==0){ 00669 conn[m_Conn*j +i] = qconn[i]; 00670 conn[m_Conn*j + i+m_BElemNodes] = new_vert[nid + m_Intervals - 2]; 00671 } 00672 else if(j==(m_Intervals-1)){ 00673 conn[m_Conn*j +i] = new_vert[nid - m_Intervals + 1]; 00674 conn[m_Conn*j + i+m_BElemNodes] = new_vert[nid]; 00675 } 00676 else { 00677 conn[m_Conn*j +i] = new_vert[nid + m_Intervals - 2*j -1]; 00678 conn[m_Conn*j + i+m_BElemNodes] = new_vert[nid + m_Intervals - 2*j -2]; 00679 } 00680 } 00681 else if(m_Conn == 4 && m_BElemNodes == 2){ // Quads 00682 int nid = node_tag_id*m_Intervals+j; 00683 if(m_Intervals == 1){ 00684 conn[m_Conn*j +i] = new_vert[nid]; 00685 conn[m_Conn*j + i+m_BElemNodes] = qconn[m_BElemNodes-i-1]; 00686 } 00687 else if(j==0){ 00688 conn[m_Conn*j +i] = qconn[m_BElemNodes-i-1]; 00689 conn[m_Conn*j + i+m_BElemNodes] = new_vert[nid + m_Intervals - 2]; 00690 } 00691 else if(j==(m_Intervals-1)){ 00692 conn[m_Conn*j +i] = new_vert[nid]; 00693 conn[m_Conn*j + m_BElemNodes + 1 -i] = new_vert[nid - m_Intervals + 1]; 00694 } 00695 else { 00696 conn[m_Conn*j +i] = new_vert[nid + m_Intervals - 2*j -2]; 00697 conn[m_Conn*j + m_BElemNodes + 1 -i] = new_vert[nid + m_Intervals - 2*j -1]; 00698 } 00699 } 00700 else if(m_Conn == 4 && m_BElemNodes == 3 && hybrid == true){ // make wedges aka prisms for tet mesh 00701 int nid = node_tag_id*m_Intervals+j; 00702 if(m_Intervals == 1){ 00703 conn[m_Conn*j +i] = qconn[i]; 00704 conn[m_Conn*j + i+m_BElemNodes] = new_vert[nid]; 00705 } 00706 else if(j==0){ 00707 conn[m_HConn*j +i] = qconn[i]; 00708 conn[m_HConn*j + i+m_BElemNodes] = new_vert[nid + m_Intervals - 2]; 00709 } 00710 else if(j==(m_Intervals-1)){ 00711 conn[m_HConn*j +i] = new_vert[nid - m_Intervals + 1]; 00712 conn[m_HConn*j + i+m_BElemNodes] = new_vert[nid]; 00713 } 00714 else { 00715 conn[m_HConn*j +i] = new_vert[nid + m_Intervals - 2*j -1]; 00716 conn[m_HConn*j + i+m_BElemNodes] = new_vert[nid + m_Intervals - 2*j -2]; 00717 } 00718 } 00719 else if(m_Conn == 3 && m_BElemNodes == 2){ // make quads for tri mesh 00720 int nid = node_tag_id*m_Intervals+j; 00721 if(m_Intervals == 1){ 00722 conn[m_Conn*j +i] = new_vert[nid]; 00723 conn[m_Conn*j + i+m_BElemNodes] = qconn[m_BElemNodes-i-1]; 00724 } 00725 else if(j==0){ 00726 conn[m_HConn*j +i] = qconn[m_BElemNodes-i-1]; 00727 conn[m_HConn*j + i+m_BElemNodes] = new_vert[nid + m_Intervals - 2]; 00728 } 00729 else if(j==(m_Intervals-1)){ 00730 conn[m_HConn*j +i] = new_vert[nid]; 00731 conn[m_HConn*j + m_BElemNodes + 1 - i] = new_vert[nid - m_Intervals + 1]; 00732 } 00733 else { 00734 conn[m_HConn*j +i] = new_vert[nid + m_Intervals - 2*j -2]; 00735 conn[m_HConn*j + m_BElemNodes + 1 - i] = new_vert[nid + m_Intervals - 2*j -1]; 00736 } 00737 } 00738 } 00739 } 00740 00741 //TODO: Set Connectivity of tet's, break prisms into 3 tets, Another loop is required. 00742 if(m_Conn == 3 && m_BElemNodes == 2 && hybrid == false){ 00743 for(int c=0; c<m_Intervals; c++){ 00744 if(tri_sch == 1){ 00745 // lower triangle 00746 tri_conn[m_Conn*c*2] = conn[c*m_HConn + 0]; 00747 tri_conn[m_Conn*c*2 + 1] = conn[c*m_HConn + 1]; 00748 tri_conn[m_Conn*c*2 + 2] = conn[c*m_HConn + 2]; 00749 // upper triangl 00750 tri_conn[m_Conn*c*2 + 3] = conn[c*m_HConn + 2]; 00751 tri_conn[m_Conn*c*2 + 4] = conn[c*m_HConn + 3]; 00752 tri_conn[m_Conn*c*2 + 5] = conn[c*m_HConn + 0]; 00753 } 00754 else if(tri_sch == 2){ 00755 // lower triangle 00756 tri_conn[m_Conn*c*2] = conn[c*m_HConn + 0]; 00757 tri_conn[m_Conn*c*2 + 1] = conn[c*m_HConn + 1]; 00758 tri_conn[m_Conn*c*2 + 2] = conn[c*m_HConn + 2]; 00759 // upper triangle 00760 tri_conn[m_Conn*c*2 + 3] = conn[c*m_HConn + 2]; 00761 tri_conn[m_Conn*c*2 + 4] = conn[c*m_HConn + 3]; 00762 tri_conn[m_Conn*c*2 + 5] = conn[c*m_HConn + 0]; 00763 } 00764 } 00765 } 00766 00767 // create boundary layer hexes 00768 for(int j=0; j< m_Intervals; j++){ 00769 //double j_hex = 0.0; 00770 // get_det_jacobian(conn, j*m_Conn, j_hex); 00771 if(m_Conn == 8){ 00772 MBERRCHK(mb->create_element(MBHEX, &conn[j*m_Conn], m_Conn, hex),mb); 00773 } 00774 else if(m_Conn==4 && m_GD ==3 && hybrid == true){ 00775 MBERRCHK(mb->create_element(MBPRISM, &conn[j*6], 6, hex),mb); 00776 } 00777 else if(m_Conn==4 && m_GD ==2){ 00778 MBERRCHK(mb->create_element(MBQUAD, &conn[j*m_Conn], m_Conn, hex),mb); 00779 } 00780 else if(m_Conn==3 && m_GD ==2 && hybrid == true){ 00781 MBERRCHK(mb->create_element(MBQUAD, &conn[j*m_HConn], m_HConn, hex),mb); 00782 } 00783 else if(m_Conn==3 && m_GD ==2 && hybrid == false){ 00784 MBERRCHK(mb->create_element(MBTRI, &tri_conn[j*m_Conn*2], m_Conn, hex),mb); 00785 MBERRCHK(mb->create_element(MBTRI, &tri_conn[j*m_Conn*2+3], m_Conn, hex1),mb); 00786 // MBERRCHK(mb->add_entities(mthis_set, &hex1, 1), mb); 00787 // MBERRCHK(mb->add_entities(smooth_set, &hex1, 1), mb); 00788 } 00789 // add this hex to a block 00790 moab::Range adj_hex_for_mat; 00791 // int hmat_id = 0; 00792 00793 if(mthis_set == 0){ 00794 // std::vector<int> hmat_id(qconn.size(), 0); 00795 // MBERRCHK(mb->tag_get_data(MNTag, &qconn[0], 1, &hmat_id) ,mb); 00796 // MBERRCHK(mb->get_adjacencies(&qconn[0], 1, m_GD, false, adj_hex_for_mat, moab::Interface::INTERSECT), mb); 00797 MBERRCHK(mb->get_adjacencies(&(*kter), 1, m_GD, false, adj_hex_for_mat, moab::Interface::INTERSECT), mb); 00798 MBERRCHK(mb->add_adjacencies(hex, adj_hex_for_mat, true), mb); 00799 00800 std::vector<int> hmat_id(adj_hex_for_mat.size(), 0); 00801 00802 // this will lead to an error, so no error checking, new adj hexes don't have matidtag 00803 mb->tag_get_data(MatIDTag, adj_hex_for_mat, &hmat_id[0]);//, mb); 00804 for(int p=0; p< (int)hmat_id.size(); p++){ 00805 if(hmat_id[p] !=0){ 00806 // this is our mat id for this hex 00807 moab::EntityHandle mat_set = 0; 00808 for (set_it = m_sets.begin(); set_it != m_sets.end(); set_it++) { 00809 mat_set = *set_it; 00810 int set_id; 00811 MBERRCHK(mb->tag_get_data(MTag, &mat_set, 1, &set_id), mb); 00812 if(set_id == hmat_id[p]) 00813 break; 00814 } 00815 MBERRCHK(mb->add_entities(mat_set, &hex, 1), mb); 00816 if(m_Conn==3 && m_GD ==2 && hybrid == false) 00817 MBERRCHK(mb->add_entities(mthis_set, &hex1, 1), mb); 00818 } 00819 } 00820 00821 00822 // here find the material set for this hex 00823 00824 // m_LogFile << "Find out the material that this hex corresponds to?? bailing out" << std::endl; 00825 // exit(0); 00826 } 00827 else { 00828 MBERRCHK(mb->add_entities(mthis_set, &hex, 1), mb); 00829 if(m_Conn==3 && m_GD ==2 && hybrid == false) 00830 MBERRCHK(mb->add_entities(mthis_set, &hex1, 1), mb); 00831 00832 } 00833 // mark entities for smoothing 00834 // MBERRCHK(mb->add_entities(smooth_set, &hex, 1), mb); 00835 // add geom dim tag 00836 MBERRCHK(mb->add_entities(geom_set, &hex, 1), mb); 00837 // // TODO: Add Local Smooting 00838 } 00839 00840 } 00841 00842 // all_elems.clear(); 00843 // moab::Range skin_verts; 00844 // MBERRCHK(mb->get_entities_by_dimension(0, 3, all_elems,true),mb); 00845 00846 // moab::Skinner skinner(mb); 00847 // skinner.find_skin(0, all_elems, 0, skin_verts); 00848 00849 // m_LogFile << "setting 'fixed'' tag = 1 on verts in the skin = " << skin_verts.size() << std::endl; 00850 00851 // // set fixed tag = 1 on all skin verts 00852 // std::vector<int> all_skin_data(skin_verts.size(), 1); 00853 // MBERRCHK(mb->tag_set_data(FTag, skin_verts, &all_skin_data[0]), mb); 00854 00855 m_LogFile << "\nTotal Jacobian calls/Min/Max: " << m_JacCalls << ", " << m_JLo << ", " << m_JHi << std::endl; 00856 00857 // save the final boundary layer mesh 00858 MBERRCHK(mb->write_mesh(m_OutFile.c_str()),mb); 00859 m_LogFile << "\n\nWrote Mesh File: " << m_OutFile << std::endl; 00860 // get the current date and time 00861 Timer.GetDateTime (szDateTime); 00862 m_LogFile << "Ending at : " << szDateTime; 00863 // report/compute the elapsed time 00864 m_LogFile << "Elapsed wall clock time: " << Timer.DiffTime () 00865 << " seconds or " << (Timer.DiffTime ())/60.0 << " mins\n"; 00866 m_LogFile << "AL2 Total CPU time used: " << (double) (clock() - sTime)/CLOCKS_PER_SEC \ 00867 << " seconds" << std::endl; 00868 00869 }
1.7.6.1 
