|
MeshKit
1.0
|
|
MeshKit
1.0
|
00001 #include "meshkit/PostBL.hpp" 00002 00003 namespace MeshKit 00004 { 00005 00006 // static registration of this mesh scheme 00007 moab::EntityType PostBL_tps[] = {moab::MBHEX, 00008 moab::MBMAXTYPE}; 00009 const moab::EntityType* PostBL::output_types() 00010 { return PostBL_tps; } 00011 00012 PostBL::PostBL( MKCore *mk, const MEntVector &me_vec) 00013 : MeshScheme( mk, me_vec), 00014 igeom(mk->igeom_instance()), imesh(mk->imesh_instance()), 00015 mb (mk->moab_instance()) 00016 // --------------------------------------------------------------------------- 00020 // --------------------------------------------------------------------------- 00021 { 00022 tri_sch = 2; 00023 m_Conn = 0; 00024 m_BElemNodes = 0; 00025 m_SurfId = -1; 00026 check_bl_edge_length = false; 00027 debug = false; 00028 hybrid = false; 00029 m_NeumannSet = -1; 00030 m_Material = 999; 00031 m_nLineNumber = 0; 00032 szComment = "!"; 00033 MAXCHARS = 300; 00034 m_JacCalls = 0; 00035 m_JLo = 0.0; 00036 m_JHi = 0.0; 00037 err = 0; 00038 fixmat = -1; 00039 } 00040 00041 PostBL::~PostBL() 00042 // --------------------------------------------------------------------------- 00046 // --------------------------------------------------------------------------- 00047 {} 00048 00049 bool PostBL::add_modelent(ModelEnt *model_ent) 00050 // --------------------------------------------------------------------------- 00054 // --------------------------------------------------------------------------- 00055 { 00056 return MeshOp::add_modelent(model_ent); 00057 } 00058 00059 void PostBL::setup_this() 00060 // --------------------------------------------------------------------------- 00064 // --------------------------------------------------------------------------- 00065 { 00066 if (debug) { 00067 m_LogFile << "\nIn setup this : " << std::endl; 00068 } 00069 00070 } 00071 00072 void PostBL::execute_this() 00073 // --------------------------------------------------------------------------- 00077 // --------------------------------------------------------------------------- 00078 { 00079 00082 int algo_no = 2; 00083 if(algo_no == 1){ 00084 m_LogFile << "\nIn execute this : creating boundary layer elements.." << std::endl; 00085 // start the timer 00086 CClock Timer; 00087 clock_t sTime = clock(); 00088 std::string szDateTime; 00089 Timer.GetDateTime (szDateTime); 00090 00091 m_LogFile << "\nStarting out at : " << szDateTime << std::endl; 00092 m_LogFile << "\n Loading meshfile: " << m_MeshFile << ".." << std::endl; 00093 00094 // load specified mesh file 00095 IBERRCHK(imesh->load(0, m_MeshFile.c_str(),0), *imesh); 00096 // m_GD = imesh->getGeometricDimension(); Doesn't work ! 00097 moab::Range all_elems, all_verts; 00098 MBERRCHK(mb->get_entities_by_dimension(0, 3, all_elems,true),mb); 00099 if (all_elems.size() == 0) 00100 m_GD = 2; 00101 else if (all_elems.size() > 0) 00102 m_GD = 3; 00103 else 00104 exit(0); 00105 all_elems.clear(); 00106 m_LogFile << "Geometric dimension of meshfile = "<< m_GD <<std::endl; 00107 00108 // obtain existing tag handles 00109 moab::Tag GDTag, GIDTag, NTag, MTag; //, STag, FTag, MNTag 00110 MBERRCHK(mb->tag_get_handle("GEOM_DIMENSION", 1, moab::MB_TYPE_INTEGER, GDTag),mb); 00111 MBERRCHK(mb->tag_get_handle("NEUMANN_SET", 1, moab::MB_TYPE_INTEGER, NTag),mb); 00112 MBERRCHK(mb->tag_get_handle("MATERIAL_SET", 1, moab::MB_TYPE_INTEGER, MTag),mb); 00113 MBERRCHK(mb->tag_get_handle("GLOBAL_ID", 1, moab::MB_TYPE_INTEGER, GIDTag),mb); 00114 // create smoothset and fixed tag for mesquite 00115 // MBERRCHK(mb->tag_get_handle("SMOOTHSET", 1, moab::MB_TYPE_INTEGER, STag, 00116 // moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT),mb); 00117 // MBERRCHK(mb->tag_get_handle("fixed", 1, moab::MB_TYPE_INTEGER, FTag, 00118 // moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT),mb); 00119 // MBERRCHK(mb->tag_get_handle("mnode", 1, moab::MB_TYPE_INTEGER, MNTag, 00120 // moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT),mb); 00121 00122 // get all the entity sets with boundary layer geom dimension, neumann sets and material sets 00123 moab::Range sets, n_sets, m_sets; 00124 m_BLDim = m_GD - 1; 00125 00126 const void* gdim[] = {&m_BLDim}; 00127 00128 MBERRCHK(mb->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &GDTag, 00129 gdim, 1 , sets, moab::Interface::INTERSECT, false), mb); 00130 MBERRCHK(mb->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &NTag, 0, 1 , n_sets),mb); 00131 MBERRCHK(mb->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &MTag, 0, 1 , m_sets),mb); 00132 00133 // Handling NeumannSets (if BL surf in input via NS) 00134 moab::Range::iterator set_it; 00135 moab::EntityHandle this_set = 0; 00136 for (set_it = n_sets.begin(); set_it != n_sets.end(); set_it++) { 00137 00138 this_set = *set_it; 00139 00140 // get entity handle of NS specified in the input file 00141 int set_id; 00142 MBERRCHK(mb->tag_get_data(NTag, &this_set, 1, &set_id), mb); 00143 if(set_id == m_NeumannSet) 00144 break; 00145 this_set = 0; 00146 } 00147 if (debug && m_NeumannSet != -1 && this_set != 0){ 00148 m_LogFile << "Looking for NS with id " << m_NeumannSet << 00149 ". Total NS found are: "<< n_sets.size() << std::endl; 00150 } 00151 00152 // For specified surface: get the all the quads and nodes in a range 00153 moab::EntityHandle s1; 00154 moab::Range quads, nodes, fixmat_ents; 00155 int dims; // variable to store global id of boundary layer specified in the input file 00156 00157 // Method 1: INPUT by NeumannSet 00158 if(m_NeumannSet != -1 && this_set != 0){ 00159 MBERRCHK(mb->get_entities_by_dimension(this_set, m_BLDim, quads,true),mb); 00160 if (quads.size() <=0){ 00161 m_LogFile << " No quads found, aborting.. " << std::endl; 00162 exit(0); 00163 } 00164 00165 if(mb->type_from_handle(*quads.rbegin()) != MBQUAD){ 00166 m_Conn = 4; 00167 m_BElemNodes = 3; 00168 } 00169 00170 MBERRCHK(mb->get_adjacencies(quads, 0, false, nodes, moab::Interface::UNION),mb); 00171 if (debug) { 00172 m_LogFile << "Found NeumannSet with id : " << m_NeumannSet << std::endl; 00173 m_LogFile << "#Quads in this surface: " << quads.size() << std::endl; 00174 m_LogFile << "#Nodes in this surface: " << nodes.size() << std::endl; 00175 m_LogFile << "#New nodes to be created:" << m_Intervals*nodes.size() << std::endl; 00176 } 00177 } 00178 // Method 2: INPUT by surface id (geom dimension) 00179 else if (m_SurfId !=-1){ 00180 for(moab::Range::iterator rit=sets.begin(); rit != sets.end(); ++rit){ 00181 s1 = *rit; 00182 MBERRCHK(mb->tag_get_data(GIDTag, &s1, 1, &dims),mb); 00183 00184 if(dims == m_SurfId && m_SurfId != -1){ 00185 MBERRCHK(mb->get_entities_by_dimension(s1, m_BLDim, quads,true),mb); 00186 if (quads.size() <=0){ 00187 m_LogFile << " No quads found, aborting.. " << std::endl; 00188 exit(0); 00189 } 00190 00191 MBERRCHK(mb->get_adjacencies(quads, 0, false, nodes, moab::Interface::UNION),mb); 00192 if (debug) { 00193 m_LogFile << "Found surface with id : " << m_SurfId << std::endl; 00194 m_LogFile << "#Quads in this surface: " << quads.size() << std::endl; 00195 m_LogFile << "#Nodes in this surface: " << nodes.size() << std::endl; 00196 m_LogFile << "#New nodes to be created:" << m_Intervals*nodes.size() << std::endl; 00197 } 00198 } 00199 } 00200 } 00201 00202 if (quads.size() == 0 || nodes.size() == 0) { 00203 m_LogFile << "Invalid boundary layer specification, aborting.." << std::endl; 00204 exit(0); 00205 } 00206 00207 // set fixed tag on all the BL nodes 00208 // MBERRCHK(mb->get_entities_by_dimension(0, 0, all_verts, true),mb); 00209 // std::vector<int> all_node_data(all_verts.size(), 0); 00210 // MBERRCHK(mb->tag_set_data(FTag, all_verts, &all_node_data[0]), mb); 00211 00212 // Handling MaterialSet 00213 moab::Range::iterator mset_it; 00214 moab::EntityHandle mthis_set; 00215 int mset_id = 0, found = 0; 00216 for (mset_it = m_sets.begin(); mset_it != m_sets.end(); mset_it++) { 00217 00218 mthis_set = *mset_it; 00219 00220 // get entity handle of MS specified in the input file 00221 MBERRCHK(mb->tag_get_data(MTag, &mthis_set, 1, &mset_id), mb); 00222 if(mset_id == m_Material){ 00223 found = 1; 00224 break; 00225 } 00226 else if(mset_id == fixmat){ 00227 MBERRCHK(mb->get_entities_by_dimension(mthis_set, m_GD, fixmat_ents ,true),mb); 00228 } 00229 00230 00231 // // get all the nodes in the material and tag bl nodes 00232 // moab::Range mat_nodes, mat_hexes; 00233 // if(m_GD == 3) 00234 // MBERRCHK(mb->get_entities_by_type(mthis_set, moab::MBHEX, mat_hexes, true), mb); 00235 // MBERRCHK(mb->get_adjacencies(mat_hexes, 0, false, mat_nodes, Interface::UNION), mb); 00236 00237 // moab::Range mat_b_nodes = intersect(nodes, mat_nodes); 00238 // std::cout << "MAT NO:" << mset_id << " has " << mat_hexes.size() << "HEXES " << mat_nodes.size() << "NODES AND " 00239 // << mat_b_nodes.size() << " NODES ON BOUNDARY" << std::endl; 00240 00241 // std::vector<int> bl_node_data(mat_b_nodes.size(), 0); 00242 // std::vector<int> node_tag_data(mat_b_nodes.size(),-1); 00243 // // don't error check, as it is supposed to give error when multiple material case is encountered 00244 // mb->tag_get_data(MNTag, mat_b_nodes, &node_tag_data[0]); 00245 // for(int i=0; i<mat_b_nodes.size(); i++){ 00246 // std::cout << node_tag_data[i] << std::endl; 00247 // // already a part of some material 00248 // if(node_tag_data[i] != -1){ 00249 // bl_node_data[i] = node_tag_data[i]+1; 00250 // } 00251 // } 00252 // MBERRCHK(mb->tag_set_data(MNTag, mat_b_nodes, &bl_node_data[0]), mb); 00253 // mat_hexes.clear(); 00254 // mat_b_nodes.clear(); 00255 // mat_nodes.clear(); 00256 mthis_set = 0; 00257 } 00258 // get tag data and print 00259 // std::vector<int> all_bl(nodes.size(),-100); 00260 // mb->tag_get_data(MNTag, nodes, &all_bl[0]); 00261 // for(int i=0; i<nodes.size(); i++){ 00262 // std::cout << " Node " << i << " tag value " << all_bl[i] << std::endl; 00263 // } 00264 if(found == 1 && m_Material !=999){ 00265 m_LogFile << "Found material set with id " << m_Material << std::endl; 00266 } 00267 else{ 00268 // No material set found, creating material set 999 for BL elements 00269 m_LogFile << "Creating material set with id " << m_Material << std::endl; 00270 MBERRCHK(mb->create_meshset(moab::MESHSET_SET, mthis_set, 1), mb); 00271 MBERRCHK(mb->tag_set_data(MTag, &mthis_set, 1, &m_Material), mb); 00272 } 00273 00274 // // placeholder for storing smoothing entities 00275 // moab::EntityHandle smooth_set; 00276 // int s_id = 100; 00277 // MBERRCHK(mb->create_meshset(moab::MESHSET_SET, smooth_set, 1), mb); 00278 // MBERRCHK(mb->tag_set_data(STag, &smooth_set, 1, &s_id), mb); 00279 00280 // placeholder for storing gd on new entities 00281 moab::EntityHandle geom_set; 00282 MBERRCHK(mb->create_meshset(moab::MESHSET_SET, geom_set, 1), mb); 00283 MBERRCHK(mb->tag_set_data(GDTag, &geom_set, 1, &m_GD), mb); 00284 00285 // declare variables before starting BL creation 00286 std::vector <bool> node_status(false); // size of verts of bl surface 00287 node_status.resize(nodes.size()); 00288 moab::Range edges, hexes, hex_edge, quad_verts; 00289 double coords_bl_quad[3], coords_new_quad[3], xdisp = 0.0, ydisp = 0.0, zdisp = 0.0; 00290 moab::EntityHandle hex, hex1; 00291 int qcount = 0; 00292 00293 //size of the following is based on element type 00294 std::vector<moab::EntityHandle> conn, qconn, adj_qconn, tri_conn, 00295 new_vert(m_Intervals*nodes.size()), old_hex_conn, adj_hexes, adj_quads, adj_hex_nodes1; 00296 moab::CartVect surf_normal(3); 00297 00298 // Now start creating New elements 00299 for (Range::iterator kter = quads.begin(); kter != quads.end(); ++kter){ 00300 qcount++; 00301 00302 if (debug){ 00303 m_LogFile << "\n\n*** QUAD: " << qcount << std::endl; 00304 } 00305 00306 std::vector<moab::EntityHandle> old_hex; 00307 MBERRCHK(mb->get_adjacencies(&(*kter), 1, m_GD, false, old_hex),mb); 00308 if((int) old_hex.size() == 0){ 00309 m_LogFile << "unable to find adjacent hex for BL quad, aborting..."; 00310 exit(0); 00311 } 00312 00313 // if fixmat specified, filter old hex, we don't have to correct both sides of the boundary 00314 if (fixmat !=0 && (int) old_hex.size() > 1){ 00315 moab::EntityHandle old_hex_set; 00316 MBERRCHK(mb->create_meshset(moab::MESHSET_SET, old_hex_set, 1), mb); 00317 MBERRCHK(mb->add_entities(old_hex_set,&old_hex[0], (int) old_hex.size()), mb); 00318 // TODO: Find a faster way of doing this 00319 MBERRCHK(mb->remove_entities(old_hex_set, fixmat_ents), mb); 00320 old_hex.clear(); 00321 old_hex.empty(); 00322 MBERRCHK(mb->get_entities_by_dimension(old_hex_set, m_GD, old_hex), mb); 00323 } 00324 // allocate space for connectivity/adjacency during the first pass of this loop 00325 if(qcount ==1){ 00326 if(mb->type_from_handle(old_hex[0]) == MBHEX){ 00327 m_Conn = 8; 00328 m_BElemNodes = 4; 00329 m_HConn = 8; 00330 //allocating based on element type 00331 conn.resize(m_Intervals*m_Conn), qconn.resize(m_BElemNodes), adj_qconn.resize(m_BElemNodes), 00332 old_hex_conn.resize(m_Conn), adj_hex_nodes1.resize(m_Conn); 00333 } 00334 else if(mb->type_from_handle(old_hex[0]) == MBTET){ 00335 m_Conn = 4; 00336 m_HConn = 6; 00337 m_BElemNodes = 3; 00338 //allocating based on element type - thrice the number of elements 00339 if(hybrid) 00340 conn.resize(m_Intervals*6); 00341 else 00342 conn.resize(2*m_Intervals*m_Conn); 00343 qconn.resize(m_BElemNodes), adj_qconn.resize(m_BElemNodes), 00344 old_hex_conn.resize(m_Conn), adj_hex_nodes1.resize(m_Conn); 00345 } 00346 else if(mb->type_from_handle(old_hex[0]) == MBQUAD){ 00347 m_Conn = 4; 00348 m_HConn = 4; 00349 m_BElemNodes = 2; 00350 //allocating based on element type 00351 conn.resize(m_Intervals*m_Conn), qconn.resize(m_BElemNodes), adj_qconn.resize(m_BElemNodes), 00352 old_hex_conn.resize(m_Conn), adj_hex_nodes1.resize(m_Conn); 00353 } 00354 else if(mb->type_from_handle(old_hex[0]) == MBTRI){ 00355 m_Conn = 3; 00356 m_HConn = 4; 00357 m_BElemNodes = 2; 00358 //allocating based on element type - twice the number of elements 00359 if(hybrid){ 00360 m_HConn = 4; 00361 conn.resize(m_Intervals*m_HConn); 00362 } 00363 else{ 00364 tri_conn.resize(2*m_Intervals*m_Conn); 00365 conn.resize(2*m_Intervals*m_Conn); 00366 } 00367 qconn.resize(m_BElemNodes), adj_qconn.resize(m_BElemNodes), 00368 old_hex_conn.resize(m_Conn), adj_hex_nodes1.resize(m_Conn); 00369 } 00370 else if(m_Conn == 0 || m_BElemNodes == 0){ 00371 m_LogFile << "MeshType is not supported by this tool" << std::endl; 00372 exit(0); 00373 } 00374 } 00375 qconn.clear(); 00376 old_hex_conn.clear(); 00377 MBERRCHK(mb->get_connectivity(&(*kter), 1, qconn),mb); 00378 MBERRCHK(mb->get_connectivity(&old_hex[0], 1, old_hex_conn),mb); 00379 00380 // get the normal to the plane of the 2D mesh 00381 if(m_GD==2 && qcount == 1) 00382 get_normal_quad (old_hex_conn, surf_normal); 00383 for (int i=0; i<m_BElemNodes; i++){ 00384 MBERRCHK(mb->get_coords(&qconn[i], 1, coords_bl_quad),mb); 00385 00386 // check to see if this node is already dealt with 00387 int blNodeId = 0; 00388 for(int n=0; n< (int) nodes.size(); n++){ 00389 if(nodes[n] == qconn[i]){ 00390 blNodeId = n; 00391 break; 00392 } 00393 else if(n == (int) nodes.size()){ 00394 m_LogFile << "QUAD doesn't have a node in BL NODES, aborting.." << std::endl; 00395 exit(0); 00396 } 00397 } 00398 if (debug){ 00399 m_LogFile << "\n*Working on Node: " << blNodeId << " coords: " << coords_bl_quad[0] 00400 << ", " << coords_bl_quad[1] << ", " << coords_bl_quad[2] << std::endl; 00401 } 00402 double temp; 00403 //create new nodes 00404 if(node_status[blNodeId] == false){ 00405 adj_hexes.clear(); 00406 MBERRCHK(mb->get_adjacencies(&qconn[i], 1, m_GD, false, adj_hexes, moab::Interface::UNION), mb); 00407 MBERRCHK(mb->get_adjacencies(&qconn[i], 1, m_BLDim, false, adj_quads, moab::Interface::UNION), mb); 00408 00409 // if fixmat specified, filter old hex, we don't have to correct both sides of the boundary 00410 if (fixmat !=0 && (int) adj_hexes.size() > 1){ 00411 moab::EntityHandle adj_hex_set; 00412 MBERRCHK(mb->create_meshset(moab::MESHSET_SET, adj_hex_set, 1), mb); 00413 MBERRCHK(mb->add_entities(adj_hex_set,&adj_hexes[0], (int) adj_hexes.size()), mb); 00414 MBERRCHK(mb->remove_entities(adj_hex_set, fixmat_ents), mb); 00415 adj_hexes.clear(); 00416 adj_hexes.empty(); 00417 MBERRCHK(mb->get_entities_by_dimension(adj_hex_set, m_GD, adj_hexes), mb); 00418 } 00419 00420 int side_number = 0, sense = 1, offset = 0; 00421 MBERRCHK(mb->side_number(old_hex[0], (*kter), side_number, sense, offset), mb); 00422 00423 moab::CartVect rt(0.0, 0.0, 0.0), v(3); 00424 moab::Range adj_qconn_r; 00425 // TODO: Add feature to add element on both sides of the boundary layer 00426 // find the normal direction where new nodes are to be created - xdisp, ydisp and zdisp 00427 for (int q=0; q < (int) quads.size(); q++){ 00428 for(int r=0; r < (int) adj_quads.size(); r++){ 00429 if (adj_quads[r] == quads[q]){ 00430 // it's a BL quad, get the normal and prepare to compute average 00431 MBERRCHK(mb->get_connectivity(&adj_quads[r], 1, adj_qconn),mb); 00432 if(m_GD==3){ 00433 get_normal_quad (adj_qconn, v); 00434 if (sense == -1) 00435 v=-v; 00436 } 00437 else if(m_GD==2){ 00438 if(sense == 1) 00439 get_normal_edge(adj_qconn, surf_normal, v); 00440 else 00441 get_normal_edge(adj_qconn, -surf_normal, v); 00442 } 00443 //TODO: Check to make sure that normal is inwards 00444 rt = rt + v; 00445 xdisp=rt[0]/rt.length(); 00446 ydisp=rt[1]/rt.length(); 00447 zdisp=rt[2]/rt.length(); 00448 } 00449 else{ 00450 // it's not a BL quad, instead of shooting the normal find the max. BL distance allowable 00451 if(check_bl_edge_length){ 00452 MBERRCHK(mb->get_connectivity(&adj_quads[r],1,adj_qconn_r), mb); 00453 find_min_edge_length(adj_qconn_r, qconn[i], nodes, m_MinEdgeLength); 00454 } 00455 } 00456 adj_qconn_r.clear(); 00457 } 00458 } 00459 // TODO: shoot normal in direction xdisp, yd.. from point coords_bl_quad 00460 // and find the point of intersection and distance with existing mesh 00461 // range nodes, coords_bl_quad and the nodes of the hex of interest here is old_hex_conn 00462 00463 // Half check of thickness validity 00464 if(m_Thickness > m_MinEdgeLength && m_MinEdgeLength > 0){ 00465 // This BL creation might fail 00466 m_LogFile << "Specified thickess is " << m_Thickness 00467 << " and BL elements edge length is " << m_MinEdgeLength 00468 << " BL creation might fail, " 00469 << "\n Set check_bl_edge_length to false to avoid this check, aborting.. " << std::endl; 00470 exit(0); 00471 } 00472 00473 00474 adj_quads.clear(); 00475 double num = m_Thickness*(m_Bias-1)*(pow(m_Bias, m_Intervals -1)); 00476 double deno = pow(m_Bias, m_Intervals) - 1; 00477 if (deno !=0) 00478 temp = num/deno; 00479 else 00480 temp = m_Thickness/m_Intervals; 00481 00482 // created BL nodes 00483 double move = 0.0; 00484 for(int j=0; j< m_Intervals; j++){ 00485 00486 move+= temp/pow(m_Bias,j); 00487 if (debug){ 00488 m_LogFile << " move:" << move; 00489 } 00490 // now compute the coords of the new vertex 00491 coords_new_quad[0] = coords_bl_quad[0]-move*xdisp; 00492 coords_new_quad[1] = coords_bl_quad[1]-move*ydisp; 00493 coords_new_quad[2] = coords_bl_quad[2]-move*zdisp; 00494 00495 int nid = blNodeId*m_Intervals+j; 00496 // Possible TODO's 00497 //TODO: Check if this vertex is possible (detect possible collision with geometry) 00498 // TODO: See posibility of using ray tracing 00499 // TODO: Parallize: Profile T-junction model and try to device an algorithm 00500 // TODO: Modularize node creation part and use doxygen for all code and design of code, python design and test cases - current functions in code: 00501 // Setup this, Execute this -- break info sub functions and classes, 00502 // prepareIO --make this optional when using python, 00503 // get normal (2d and 3d) -- can be combined to one function 00504 // get det jacobian (hex elements) --needs check for other elements 00505 // 00506 MBERRCHK(mb->create_vertex(coords_new_quad, new_vert[nid]), mb); 00507 if (debug){ 00508 m_LogFile << std::setprecision (3) << std::scientific << " : created node:" << (nid + 1) 00509 << " of " << new_vert.size() << " new nodes:- coords: " << coords_new_quad[0] 00510 << ", " << coords_new_quad[1] << ", " << coords_new_quad[2] << std::endl; 00511 } 00512 } 00513 node_status[blNodeId] = true; 00514 } 00515 00516 //populate the connectivity after creating nodes for this BL node 00517 for(int j=0; j< m_Intervals; j++){ 00518 if(m_Conn == 8 && m_BElemNodes == 4){ 00519 int nid = blNodeId*m_Intervals+j; 00520 if(j==0) // set connectivity of boundary layer hex 00521 conn[m_Conn*j + i+m_BElemNodes] = qconn[i]; 00522 else 00523 conn[m_Conn*j + i+m_BElemNodes] = new_vert[nid-1]; 00524 conn[m_Conn*j +i] = new_vert[nid]; 00525 } 00526 else if(m_Conn == 4 && m_BElemNodes == 2){ 00527 int nid = blNodeId*m_Intervals+j; 00528 if(j==0) // set connectivity of boundary layer hex 00529 conn[m_Conn*j + i+m_BElemNodes] = qconn[m_BElemNodes-i-1]; 00530 else 00531 conn[m_Conn*j + m_BElemNodes + 1 - i] = new_vert[nid-1]; 00532 conn[m_Conn*j +i] = new_vert[nid]; 00533 } 00534 else if(m_Conn == 4 && m_BElemNodes == 3 && hybrid == true){ // make wedges aka prisms for tet mesh 00535 int nid = blNodeId*m_Intervals+j; 00536 if(j==0) // set connectivity of boundary layer hex 00537 conn[m_HConn*j + i+m_BElemNodes] = qconn[i]; 00538 else 00539 conn[m_HConn*j + i+m_BElemNodes] = new_vert[nid-1]; 00540 conn[m_HConn*j +i] = new_vert[nid]; 00541 } 00542 else if(m_Conn == 3 && m_BElemNodes == 2){ // make quads for tri mesh 00543 int nid = blNodeId*m_Intervals+j; 00544 if(j==0) // set connectivity of boundary layer hex 00545 conn[m_HConn*j + i+m_BElemNodes] = qconn[m_BElemNodes-i-1]; 00546 else 00547 conn[m_HConn*j + m_BElemNodes + 1 - i] = new_vert[nid-1]; 00548 conn[m_HConn*j +i] = new_vert[nid]; 00549 } 00550 00551 // Another block for setting up fixed tag for smoothing using Mesquite 00552 // int fix0 = 0; 00553 // int nid = blNodeId*m_Intervals+j; 00554 // MBERRCHK(mb->tag_set_data(FTag, &new_vert[nid], 1, &fix0), mb); 00555 } 00556 00557 // if a hex does have a quad on BL, but, has a node or edge on the boundary layer 00558 for(int k=0; k < (int) adj_hexes.size(); k++){ 00559 adj_hex_nodes1.clear(); 00560 MBERRCHK(mb->get_connectivity(&adj_hexes[k], 1, adj_hex_nodes1), mb); 00561 std::vector<EntityHandle> inodes; 00562 int var = 0; 00563 for(int n = 0; n < (int) nodes.size(); n++){ 00564 for(int m = 0; m < m_Conn; m++){ 00565 if(nodes[n] == adj_hex_nodes1[m]){ 00566 inodes.push_back(adj_hex_nodes1[m]); 00567 var++; 00568 } 00569 } 00570 } 00571 // doesn't have a quad or when inodes is 0 it mean this is a newly created hex 00572 if((int) inodes.size() != m_BElemNodes && (int) inodes.size() != 0){ 00573 if(debug){ 00574 m_LogFile << "Hex on BL surface is connected by a node or edge" << 00575 "\n nodes on BL - " << inodes.size() << std::endl; 00576 } 00577 // mark this hex and set it's connectivity later 00578 for(int p = 0; p < m_Conn; p++){ 00579 if(qconn[i] == adj_hex_nodes1[p]){ 00580 adj_hex_nodes1[p] = conn[m_HConn*(m_Intervals-1)+i]; 00581 } 00582 } 00583 // check to see if this node position change also changes lower dimension entities 00584 std::vector<EntityHandle> adj1, adj2; 00585 EntityHandle one = qconn[i]; 00586 EntityHandle two = adj_hexes[k]; 00587 std::vector<EntityHandle> from_entities; 00588 from_entities.push_back(one); 00589 from_entities.push_back(two); 00590 if(m_BLDim == 2){ 00591 MBERRCHK(mb->get_adjacencies(&from_entities[0], 2, 1, false, adj1), mb); 00592 MBERRCHK(mb->delete_entities(&adj1[0], (int) adj1.size()), mb); 00593 } 00594 MBERRCHK(mb->get_adjacencies(&from_entities[0], 2, m_BLDim, false, adj2), mb); 00595 MBERRCHK(mb->delete_entities(&adj2[0], (int) adj2.size()), mb); 00596 MBERRCHK(mb->set_connectivity(adj_hexes[k], &adj_hex_nodes1[0], m_Conn), mb); 00597 } 00598 double j_ahex = 0.0; 00599 get_det_jacobian(adj_hex_nodes1, 0, j_ahex); 00600 // mark entities for smoothing 00601 // MBERRCHK(mb->add_entities(smooth_set, &adj_hexes[k], 1), mb); 00602 inodes.clear(); 00603 } 00604 } // Loop thru BL element nodes ends 00605 00606 //TODO: Set Connectivity of tet's, break prisms into 3 tets, Another loop is required. 00607 if(m_Conn == 3 && m_BElemNodes == 2 && hybrid == false){ 00608 for(int c=0; c<m_Intervals; c++){ 00609 if(tri_sch == 1){ 00610 // lower triangle 00611 tri_conn[m_Conn*c*2] = conn[c*m_HConn]; 00612 tri_conn[m_Conn*c*2 + 1] = conn[c*m_HConn + 1]; 00613 tri_conn[m_Conn*c*2 + 2] = conn[c*m_HConn + 3]; 00614 // upper triangle 00615 tri_conn[m_Conn*c*2 + 3] = conn[c*m_HConn + 1]; 00616 tri_conn[m_Conn*c*2 + 4] = conn[c*m_HConn + 2]; 00617 tri_conn[m_Conn*c*2 + 5] = conn[c*m_HConn + 3]; 00618 } 00619 else if(tri_sch == 2){ 00620 // lower triangle 00621 tri_conn[m_Conn*c*2] = conn[c*m_HConn]; 00622 tri_conn[m_Conn*c*2 + 1] = conn[c*m_HConn + 1]; 00623 tri_conn[m_Conn*c*2 + 2] = conn[c*m_HConn + 2]; 00624 // upper triangle 00625 tri_conn[m_Conn*c*2 + 3] = conn[c*m_HConn + 2]; 00626 tri_conn[m_Conn*c*2 + 4] = conn[c*m_HConn + 3]; 00627 tri_conn[m_Conn*c*2 + 5] = conn[c*m_HConn]; 00628 } 00629 } 00630 } 00631 00632 if (debug) 00633 { 00634 m_LogFile << "\nsetting connectivity of the old BL hex " << std::endl; 00635 } 00636 // First replace BL nodes (part of old hex) with newly created nodes, then set connectivity of the old_hex 00637 // deleting old connectivities 00638 //TODO: Check to see if this step can be avoided for some cases 00639 for(int p=0; p<m_Conn; p++){ 00640 for(int q=0; q<m_BElemNodes; q++){ 00641 if (old_hex_conn[p] == qconn[q]){ 00642 old_hex_conn[p] = conn[m_HConn*(m_Intervals-1) + q]; 00643 // check to see if this node position change also changes lower dimension entities 00644 std::vector<EntityHandle> adj1, adj2; 00645 EntityHandle one = qconn[q]; 00646 EntityHandle two = old_hex[0]; 00647 std::vector<EntityHandle> from_entities; 00648 from_entities.push_back(one); 00649 from_entities.push_back(two); 00650 if(m_BLDim == 2){ 00651 MBERRCHK(mb->get_adjacencies(&from_entities[0], 2, 1, false, adj1), mb); 00652 MBERRCHK(mb->delete_entities(&adj1[0], (int) adj1.size()), mb); 00653 } 00654 MBERRCHK(mb->get_adjacencies(&from_entities[0], 2, m_BLDim, false, adj2), mb); 00655 MBERRCHK(mb->delete_entities(&adj2[0], (int) adj2.size()), mb); 00656 } 00657 } 00658 } 00659 double j_old_hex = 0.0; 00660 get_det_jacobian(old_hex_conn, 0, j_old_hex); 00661 MBERRCHK(mb->set_connectivity(old_hex[0], &old_hex_conn[0], m_Conn), mb); 00662 old_hex.clear(); 00663 00664 // mark entities for smoothing 00665 // MBERRCHK(mb->add_entities(smooth_set, &old_hex[0], 1), mb); 00666 00667 if (debug){ 00668 m_LogFile << "creating new boundary layer hexes" << std::endl; 00669 } 00670 // create boundary layer hexes 00671 for(int j=0; j< m_Intervals; j++){ 00672 double j_hex = 0.0; 00673 get_det_jacobian(conn, j*m_Conn, j_hex); 00674 if(m_Conn == 8){ 00675 MBERRCHK(mb->create_element(MBHEX, &conn[j*m_Conn], m_Conn, hex),mb); 00676 } 00677 else if(m_Conn==4 && m_GD ==3 && hybrid == true){ 00678 MBERRCHK(mb->create_element(MBPRISM, &conn[j*6], 6, hex),mb); 00679 } 00680 else if(m_Conn==4 && m_GD ==2){ 00681 MBERRCHK(mb->create_element(MBQUAD, &conn[j*m_Conn], m_Conn, hex),mb); 00682 } 00683 else if(m_Conn==3 && m_GD ==2 && hybrid == true){ 00684 MBERRCHK(mb->create_element(MBQUAD, &conn[j*m_HConn], m_HConn, hex),mb); 00685 } 00686 else if(m_Conn==3 && m_GD ==2 && hybrid == false){ 00687 MBERRCHK(mb->create_element(MBTRI, &tri_conn[j*m_Conn*2], m_Conn, hex),mb); 00688 MBERRCHK(mb->create_element(MBTRI, &tri_conn[j*m_Conn*2+3], m_Conn, hex1),mb); 00689 MBERRCHK(mb->add_entities(mthis_set, &hex1, 1), mb); 00690 // MBERRCHK(mb->add_entities(smooth_set, &hex1, 1), mb); 00691 } 00692 // add this hex to a block 00693 MBERRCHK(mb->add_entities(mthis_set, &hex, 1), mb); 00694 // mark entities for smoothing 00695 // MBERRCHK(mb->add_entities(smooth_set, &hex, 1), mb); 00696 // add geom dim tag 00697 MBERRCHK(mb->add_entities(geom_set, &hex, 1), mb); 00698 // TODO: Add Local Smooting 00699 } 00700 } // Loop thru quads ends 00701 00702 // TODO: Global smoothing works after PostBL step; make global smoothing available as an option inside 00703 // get the skin of the entities 00704 // moab::Range skin_verts; 00705 // all_elems.clear(); 00706 // MBERRCHK(mb->get_entities_by_dimension(0, 3, all_elems,true),mb); 00707 00708 // moab::Skinner skinner(mb); 00709 // skinner.find_skin(all_elems, 0, skin_verts); 00710 00711 // m_LogFile << "setting 'fixed'' tag = 1 on verts in the skin = " << skin_verts.size() << std::endl; 00712 00713 // // set fixed tag = 1 on all skin verts 00714 // std::vector<int> all_skin_data(skin_verts.size(), 1); 00715 // MBERRCHK(mb->tag_set_data(FTag, skin_verts, &all_skin_data[0]), mb); 00716 00717 m_LogFile << "\nTotal Jacobian calls/Min/Max: " << m_JacCalls << ", " << m_JLo << ", " << m_JHi << std::endl; 00718 00719 // save the final boundary layer mesh 00720 MBERRCHK(mb->write_mesh(m_OutFile.c_str()),mb); 00721 m_LogFile << "\n\nWrote Mesh File: " << m_OutFile << std::endl; 00722 // get the current date and time 00723 Timer.GetDateTime (szDateTime); 00724 m_LogFile << "Ending at : " << szDateTime; 00725 // report/compute the elapsed time 00726 m_LogFile << "Elapsed wall clock time: " << Timer.DiffTime () 00727 << " seconds or " << (Timer.DiffTime ())/60.0 << " mins\n"; 00728 m_LogFile << "Total CPU time used: " << (double) (clock() - sTime)/CLOCKS_PER_SEC \ 00729 << " seconds" << std::endl; 00730 } 00731 else{ 00732 Algo2(); 00733 } 00734 } 00735 00736 void PostBL::PrepareIO (int argc, char *argv[], std::string TestDir) 00737 // --------------------------------------------------------------------------- 00741 // --------------------------------------------------------------------------- 00742 { 00743 // set and open input output files 00744 bool bDone = false; 00745 do{ 00746 if (2 == argc) { 00747 m_InputFile = argv[1]; 00748 m_LogName = m_InputFile + ".log"; 00749 } 00750 else if (1 == argc){ 00751 m_LogFile << "\nRunning default case:\n" << std::endl; 00752 00753 m_InputFile = TestDir + "/" + (char *)DEFAULT_TEST_POSTBL; 00754 m_LogName = m_InputFile + ".log"; 00755 } 00756 00757 // open input file for reading 00758 m_FileInput.open (m_InputFile.c_str(), std::ios::in); 00759 if (!m_FileInput){ 00760 m_LogFile << "Unable to open file: " << m_InputFile << std::endl; 00761 m_FileInput.clear (); 00762 exit(1); 00763 } 00764 else 00765 bDone = true; // file opened successfully 00766 00767 // open the log file for dumping debug/output statements 00768 m_LogFile.coss.open (m_LogName.c_str(), std::ios::out); 00769 if (!m_LogFile.coss){ 00770 m_LogFile << "Unable to open file: " << m_LogName << std::endl; 00771 m_LogFile.coss.clear (); 00772 exit(1); 00773 } 00774 else 00775 bDone = true; // file opened successfully 00776 m_LogFile << '\n'; 00777 m_LogFile << "\t\t---------------------------------------------------------" << '\n'; 00778 m_LogFile << "\t\t Tool to generate Post-mesh Boundary Layers " << '\n'; 00779 m_LogFile << "\t\t\t\tArgonne National Laboratory" << '\n'; 00780 m_LogFile << "\t\t\t\t 2012 " << '\n'; 00781 m_LogFile << "\t\t---------------------------------------------------------" << '\n'; 00782 m_LogFile << "\nsee README file for using the program and details on various cards.\n"<< std::endl; 00783 00784 }while (!bDone); 00785 00786 // Get the meshfile name, surface(s), thickness, intervals and bias 00787 CParser Parse; 00788 00789 // count the total number of cylinder commands in each pincellh 00790 for(;;){ 00791 if (!Parse.ReadNextLine (m_FileInput, m_nLineNumber, szInputString, 00792 MAXCHARS, szComment)) 00793 IOErrorHandler (INVALIDINPUT); 00794 00795 // Get tri scheme 00796 if (szInputString.substr(0,9) == "trischeme"){ 00797 std::istringstream szFormatString (szInputString); 00798 szFormatString >> m_Card >> tri_sch; 00799 if(szFormatString.fail()) 00800 IOErrorHandler(INVALIDINPUT); 00801 m_LogFile << m_Card << " name read: "<< tri_sch << std::endl; 00802 } 00803 // Get hybrid 00804 if (szInputString.substr(0,6) == "hybrid"){ 00805 std::istringstream szFormatString (szInputString); 00806 szFormatString >> m_Card >> hybrid; 00807 if(szFormatString.fail()) 00808 IOErrorHandler(INVALIDINPUT); 00809 m_LogFile << m_Card << " name read: "<< hybrid << std::endl; 00810 } 00811 // Get hybrid 00812 if (szInputString.substr(0,6) == "fixmat"){ 00813 std::istringstream szFormatString (szInputString); 00814 szFormatString >> m_Card >> fixmat; 00815 if(szFormatString.fail()) 00816 IOErrorHandler(INVALIDINPUT); 00817 m_LogFile << m_Card << " name read: "<< fixmat << std::endl; 00818 } 00819 // Get MeshFile name 00820 if (szInputString.substr(0,8) == "meshfile"){ 00821 std::istringstream szFormatString (szInputString); 00822 szFormatString >> m_Card >> m_MeshFile; 00823 if(szFormatString.fail()) 00824 IOErrorHandler(INVALIDINPUT); 00825 m_LogFile << m_Card << " name read: "<< m_MeshFile << std::endl; 00826 if (argc == 1){ 00827 m_MeshFile = TestDir + "/" + m_MeshFile; 00828 } 00829 } 00830 // Get BL surface 00831 if (szInputString.substr(0,8) == "surfaces"){ 00832 std::istringstream szFormatString (szInputString); 00833 szFormatString >> m_Card >> m_SurfId; 00834 if(m_SurfId < 0 || szFormatString.fail()) 00835 IOErrorHandler(INVALIDINPUT); 00836 m_LogFile << m_Card << " read: "<< m_SurfId <<std::endl; 00837 } 00838 // Get BL surface via neumann set or sideset 00839 if (szInputString.substr(0,10) == "neumannset"){ 00840 std::istringstream szFormatString (szInputString); 00841 szFormatString >> m_Card >> m_NeumannSet; 00842 if(m_NeumannSet < 0 || szFormatString.fail()) 00843 IOErrorHandler(INVALIDINPUT); 00844 m_LogFile << m_Card << " read: "<< m_NeumannSet <<std::endl; 00845 } 00846 // Get BL material (block) number 00847 if (szInputString.substr(0,8) == "material"){ 00848 std::istringstream szFormatString (szInputString); 00849 szFormatString >> m_Card >> m_Material; 00850 if(m_Material < 0 || szFormatString.fail()) 00851 IOErrorHandler(INVALIDINPUT); 00852 m_LogFile << m_Card << " read: "<< m_Material <<std::endl; 00853 } 00854 00855 // Get thickness 00856 if (szInputString.substr(0,9) == "thickness"){ 00857 std::istringstream szFormatString (szInputString); 00858 szFormatString >> m_Card >> m_Thickness; 00859 if(m_Thickness < 0 || szFormatString.fail()) 00860 IOErrorHandler(INVALIDINPUT); 00861 m_LogFile << m_Card << " read: "<< m_Thickness <<std::endl; 00862 } 00863 // Get intervals 00864 if (szInputString.substr(0,9) == "intervals"){ 00865 std::istringstream szFormatString (szInputString); 00866 szFormatString >> m_Card >> m_Intervals; 00867 if(m_Intervals < 0 || szFormatString.fail()) 00868 IOErrorHandler(INVALIDINPUT); 00869 m_LogFile << m_Card << " read: "<< m_Intervals <<std::endl; 00870 } 00871 // Get bias 00872 if (szInputString.substr(0,4) == "bias"){ 00873 std::istringstream szFormatString (szInputString); 00874 szFormatString >> m_Card >> m_Bias; 00875 if(m_Bias < 0 || szFormatString.fail()) 00876 IOErrorHandler(INVALIDINPUT); 00877 m_LogFile << m_Card << " read: "<< m_Bias <<std::endl; 00878 } 00879 // Output file name 00880 if (szInputString.substr(0,7) == "outfile"){ 00881 std::istringstream szFormatString (szInputString); 00882 szFormatString >> m_Card >> m_OutFile; 00883 if(szFormatString.fail()) 00884 IOErrorHandler(INVALIDINPUT); 00885 m_LogFile << m_Card << " name read: "<< m_OutFile <<std::endl; 00886 } 00887 // Debug flag 00888 if (szInputString.substr(0,5) == "debug"){ 00889 std::istringstream szFormatString (szInputString); 00890 szFormatString >> m_Card >> debug; 00891 if(szFormatString.fail()) 00892 IOErrorHandler(INVALIDINPUT); 00893 m_LogFile << m_Card << " name read: "<< debug <<std::endl; 00894 } 00895 if (szInputString.substr(0,3) == "end"){ 00896 break; 00897 } 00898 } 00899 } 00900 00901 void PostBL::IOErrorHandler (ErrorStates ECode) const 00902 // --------------------------------------------------------------------------- 00906 // --------------------------------------------------------------------------- 00907 { 00908 std::cerr << '\n'; 00909 if (ECode == INVALIDINPUT) // invalid input 00910 std::cerr << "Invalid input."; 00911 else 00912 std::cerr << "Unknown error ...?"; 00913 00914 std::cerr << '\n' << "Error in input file line : " << m_nLineNumber; 00915 std::cerr << std::endl; 00916 exit (1); 00917 } 00918 00919 void PostBL::get_normal_quad (std::vector<EntityHandle>conn, moab::CartVect &v) 00920 // --------------------------------------------------------------------------- 00924 // --------------------------------------------------------------------------- 00925 { 00926 moab::CartVect coords[3]; 00927 MBERRCHK(mb->get_coords(&conn[0], 3, (double*) &coords[0]), mb); 00928 moab::CartVect AB(coords[1] - coords[0]); 00929 moab::CartVect BC(coords[2] - coords[1]); 00930 moab::CartVect normal = AB*BC; 00931 normal = normal/normal.length(); 00932 v = normal; 00933 } 00934 00935 void PostBL::get_normal_edge (std::vector<EntityHandle>conn, moab::CartVect BC, moab::CartVect &v) 00936 // --------------------------------------------------------------------------- 00940 // --------------------------------------------------------------------------- 00941 { 00942 moab::CartVect coords[2]; 00943 MBERRCHK(mb->get_coords(&conn[0], 2, (double*) &coords[0]), mb); 00944 moab::CartVect AB(coords[1] - coords[0]); 00945 moab::CartVect normal = AB*BC; 00946 normal = normal/normal.length(); 00947 v = normal; 00948 } 00949 00950 void PostBL::find_min_edge_length(moab::Range adj_qn, moab::EntityHandle node , moab::Range bl_nodes, double &e_len) 00951 // --------------------------------------------------------------------------- 00955 // --------------------------------------------------------------------------- 00956 { 00957 // get nodes adj(a) to BL node 00958 moab::CartVect coords[1]; 00959 double len = 0; e_len = 0; 00960 MBERRCHK(mb->get_coords(&node, 1, (double*) &coords[0]), mb); 00961 moab::Range non_bl; 00962 non_bl = subtract(adj_qn, bl_nodes); 00963 int n_non_bl = (int) non_bl.size(); 00964 00965 // if there are no bl nodes - this case has already been dealt with 00966 if(non_bl.size() > 0){ 00967 moab::CartVect non_bl_coords[4]; 00968 for(int i=0; i< n_non_bl; i++){ 00969 MBERRCHK(mb->get_coords(non_bl,(double*) &non_bl_coords[0]), mb); 00970 moab::CartVect edge(coords[0] - non_bl_coords[0]); 00971 len = edge.length(); 00972 if(i==0) 00973 e_len = len; 00974 if(len < e_len) 00975 e_len = len; 00976 } 00977 } 00978 // m_LogFile << " node minimum edge length" << e_len << std::endl; 00979 } 00980 00981 void PostBL::get_det_jacobian(std::vector<moab::EntityHandle> conn, int offset, double &AvgJ) 00982 // --------------------------------------------------------------------------- 00986 // --------------------------------------------------------------------------- 00987 { 00988 //TODO: Add quality check for tri/quad and pyramids 00989 if(m_Conn ==8){ 00990 ++m_JacCalls; 00991 moab::CartVect vertex[8], xi; 00992 mstream m_LogFile; 00993 MBERRCHK(mb->get_coords(&conn[offset], 8, (double*) &vertex[0]), mb); 00994 00995 double corner[8][3] = { { -1, -1, -1 }, 00996 { 1, -1, -1 }, 00997 { 1, 1, -1 }, 00998 { -1, 1, -1 }, 00999 { -1, -1, 1 }, 01000 { 1, -1, 1 }, 01001 { 1, 1, 1 }, 01002 { -1, 1, 1 } }; 01003 01004 for (unsigned j = 0; j < 8; ++j) { 01005 xi[0] = corner[j][0]; 01006 xi[1] = corner[j][1]; 01007 xi[2] = corner[j][2]; 01008 Matrix3 J(0.0); 01009 double detJ = 0; 01010 for (unsigned i = 0; i < 8; ++i) { 01011 const double xi_p = 1 + xi[0]*corner[i][0]; 01012 const double eta_p = 1 + xi[1]*corner[i][1]; 01013 const double zeta_p = 1 + xi[2]*corner[i][2]; 01014 const double dNi_dxi = corner[i][0] * eta_p * zeta_p; 01015 const double dNi_deta = corner[i][1] * xi_p * zeta_p; 01016 const double dNi_dzeta = corner[i][2] * xi_p * eta_p; 01017 J(0,0) += dNi_dxi * vertex[i][0]; 01018 J(1,0) += dNi_dxi * vertex[i][1]; 01019 J(2,0) += dNi_dxi * vertex[i][2]; 01020 J(0,1) += dNi_deta * vertex[i][0]; 01021 J(1,1) += dNi_deta * vertex[i][1]; 01022 J(2,1) += dNi_deta * vertex[i][2]; 01023 J(0,2) += dNi_dzeta * vertex[i][0]; 01024 J(1,2) += dNi_dzeta * vertex[i][1]; 01025 J(2,2) += dNi_dzeta * vertex[i][2]; 01026 } 01027 J *= 0.125; 01028 detJ = J.determinant(); 01029 if(detJ <= 0.0){ 01030 m_LogFile << "We've negative jacobian at the hex corner: "<< j+1 << std::endl; 01031 exit(0); 01032 } 01033 AvgJ+=detJ; 01034 } 01035 AvgJ/=8; 01036 if(m_JacCalls == 1){ 01037 m_JLo = AvgJ; 01038 m_JHi = AvgJ; 01039 } 01040 else if(AvgJ < m_JLo){ 01041 m_JLo = AvgJ; 01042 } 01043 else if(AvgJ > m_JHi){ 01044 m_JHi = AvgJ; 01045 } 01046 } 01047 } 01048 } // namespace MeshKit 01049 01050
1.7.6.1 
