MeshKit
1.0
|
00001 #include "meshkit/PostBL.hpp" 00002 namespace MeshKit 00003 { 00004 int PostBL:: compute_normals() 00005 // --------------------------------------------------------------------------- 00009 // --------------------------------------------------------------------------- 00010 { 00011 // COMPUTE NORMALS 00012 // get tag data and print 00013 all_bl.resize(nodes.size(), 0); 00014 mb->tag_get_data(MNTag, nodes, &all_bl[0]); 00015 int count = -1; 00016 for (Range::iterator kter = nodes.begin(); kter != nodes.end(); ++kter){ 00017 ++count; 00018 // only one normal in case of single material 00019 double xdisp = 0.0, ydisp = 0.0, zdisp = 0.0; 00020 00021 // check if node belongs to one or more materials 00022 if(all_bl[count] == 0){ 00023 moab::Range adj_for_normal; 00024 MBERRCHK(mb->get_adjacencies(&(*kter), 1, m_BLDim, false, adj_for_normal, Interface::UNION), mb); 00025 //create the coordinates of the innermost node corresponding to this node 00026 moab::Range adj_for_norm = intersect(quads, adj_for_normal); 00027 // now compute the average normal direction for this vertex 00028 moab::CartVect rt(0.0, 0.0, 0.0), v(0.0, 0.0, 0.0); 00029 00030 for(Range::iterator qter = adj_for_norm.begin(); qter != adj_for_norm.end(); ++qter){ 00031 MBERRCHK(mb->get_connectivity(&(*qter), 1, adj_qconn),mb); 00032 00033 int side_number = 0, sense = 1, offset = 0; 00034 mb->side_number(old_hex[0], (*qter), side_number, sense, offset); 00035 00036 if(m_GD==3){ 00037 get_normal_quad (adj_qconn, v); 00038 if(sense == 1){ 00039 // do nothing 00040 } 00041 else{ 00042 v=-v; 00043 } 00044 } 00045 else if(m_GD==2){ 00046 if(sense == 1){ 00047 get_normal_edge(adj_qconn, surf_normal, v); 00048 } 00049 else{ 00050 get_normal_edge(adj_qconn, -surf_normal, v); 00051 } 00052 } 00053 rt = rt + v; 00054 } 00055 if(rt.length() !=0){ 00056 xdisp=rt[0]/rt.length(); 00057 ydisp=rt[1]/rt.length(); 00058 zdisp=rt[2]/rt.length(); 00059 } 00060 else{ 00061 xdisp=0.0; 00062 ydisp=0.0; 00063 zdisp=0.0; 00064 } 00065 } 00066 else if(all_bl[count] > 0 && fixmat != -1){ // node belongs to more than one material and fixmat specified 00067 // NODE B/W MATERIALS 00068 moab::Range adj_for_normal; 00069 MBERRCHK(mb->get_adjacencies(&(*kter), 1, m_BLDim, false, adj_for_normal, Interface::UNION), mb); 00070 //create the coordinates of the innermost node corresponding to this node 00071 moab::Range adj_for_norm = intersect(quads, adj_for_normal); 00072 // now compute the average normal direction for this vertex 00073 moab::CartVect rt(0.0, 0.0, 0.0), v(0.0, 0.0, 0.0); 00074 00075 for(Range::iterator qter = adj_for_norm.begin(); qter != adj_for_norm.end(); ++qter){ 00076 MBERRCHK(mb->get_connectivity(&(*qter), 1, adj_qconn),mb); 00077 00078 moab::Range this_quad_hex; 00079 MBERRCHK(mb->get_adjacencies(&(*qter), 1, m_GD, false, this_quad_hex, moab::Interface::UNION),mb); 00080 moab::Range quad_hex = intersect(fixmat_ents, this_quad_hex); 00081 00082 int side_number = 0, sense = 1, offset = 0; 00083 Range::iterator hexter = quad_hex.begin(); 00084 mb->side_number(*hexter, (*qter), side_number, sense, offset); 00085 00086 if(m_GD==3){ 00087 get_normal_quad (adj_qconn, v); 00088 if(sense == 1 ){ 00089 v=-v; 00090 } 00091 else{ 00092 // do nothing 00093 } 00094 } 00095 else if(m_GD==2){ 00096 if(sense == 1 ){ 00097 get_normal_edge(adj_qconn, -surf_normal, v); 00098 } 00099 else{ 00100 get_normal_edge(adj_qconn, surf_normal, v); 00101 } 00102 } 00103 rt = rt + v; 00104 } 00105 if(rt.length() !=0){ 00106 xdisp=rt[0]/rt.length(); 00107 ydisp=rt[1]/rt.length(); 00108 zdisp=rt[2]/rt.length(); 00109 } 00110 else{ 00111 xdisp=0.0; 00112 ydisp=0.0; 00113 zdisp=0.0; 00114 } 00115 } 00116 else if(all_bl[count] > 0 && fixmat == -1){ // node belongs to more than one material and fixmat not specified 00117 // NODE ON BOUNDARY 00118 // get the edge that is not on the boundary and count how many such edges we have 00119 moab::Range adj_for_normal; 00120 int nEdgeDim = 1; 00121 MBERRCHK(mb->get_adjacencies(&(*kter), 1, nEdgeDim, true, adj_for_normal, Interface::UNION), mb); 00122 moab::Range edge_normal; 00123 if(m_GD == 2) 00124 edge_normal = subtract(adj_for_normal, quads); 00125 else if(m_GD == 3) 00126 edge_normal = subtract(adj_for_normal, edges); 00127 00128 00129 if(edge_normal.size() > 1){ 00130 double ncoord[3]; 00131 MBERRCHK(mb->get_coords(&(*kter), 1, ncoord),mb); 00132 00133 m_LogFile << "Multiple normals are needed at: " << ncoord[0] 00134 << ", " << ncoord[1] << ", " << ncoord[2] << " #normals " << edge_normal.size() << std::endl; 00135 exit(0); 00136 } 00137 else{ 00138 m_LogFile << "We've one edge seperating materials 1, edge normal is needed" << edge_normal.size() << std::endl; 00139 moab::Range edge_conn; 00140 MBERRCHK(mb->get_connectivity(&(*edge_normal.begin()), 1, edge_conn),mb); 00141 // now get the normal direction for this edge 00142 moab::CartVect coords[2], bl_coord[1]; 00143 moab::CartVect AB; 00144 MBERRCHK(mb->get_coords(&(*kter), 1, (double*) &bl_coord[0]), mb); 00145 MBERRCHK(mb->get_coords(edge_conn, (double*) &coords[0]), mb); 00146 for(int d = 0; d<2; d++){ 00147 if(bl_coord[0][0] == coords[d][0] && 00148 bl_coord[0][1] == coords[d][1] && 00149 bl_coord[0][2] == coords[d][2]){ 00150 // do nothing 00151 } 00152 else{ 00153 AB = (bl_coord[0] - coords[d]); 00154 } 00155 } 00156 00157 xdisp=AB[0]/AB.length(); 00158 ydisp=AB[1]/AB.length(); 00159 zdisp=AB[2]/AB.length(); 00160 } 00161 } 00162 else if(all_bl[count] < 0 ){ 00163 m_LogFile << "Material must have an associated BLNode: Error! - shouldn't have gotten here: " << count << std::endl; 00164 } 00165 00166 // after the normal computation is done create new BL nodes 00167 double coords_bl_quad[3], move = 0.0; 00168 MBERRCHK(mb->get_coords(&(*kter), 1, coords_bl_quad),mb); 00169 00170 double temp = 0; 00171 double num = m_Thickness*(m_Bias-1)*(pow(m_Bias, m_Intervals -1)); 00172 double deno = pow(m_Bias, m_Intervals) - 1; 00173 if (deno !=0) 00174 temp = num/deno; 00175 else 00176 temp = m_Thickness/m_Intervals; 00177 00178 // loop thru intervals to create BL nodes 00179 for(int j=0; j< m_Intervals; j++){ 00180 00181 move+= temp/pow(m_Bias,j); 00182 if (debug){ 00183 m_LogFile << " move:" << move; 00184 } 00185 // now compute the coords of the new vertex 00186 coords_new_quad[0] = coords_bl_quad[0]-move*xdisp; 00187 coords_new_quad[1] = coords_bl_quad[1]-move*ydisp; 00188 coords_new_quad[2] = coords_bl_quad[2]-move*zdisp; 00189 00190 int nid = count*m_Intervals+j; 00191 // Possible TODO's 00192 //TODO: Check if this vertex is possible (detect possible collision with geometry) 00193 // TODO: See posibility of using ray tracing 00194 // TODO: Parallize: Profile T-junction model and try to device an algorithm 00195 // TODO: Modularize node creation part and use doxygen for all code and design of code, python design and test cases - current functions in code: 00196 // Setup this, Execute this -- break info sub functions and classes, 00197 // prepareIO --make this optional when using python, 00198 // get normal (2d and 3d) -- can be combined to one function 00199 // get det jacobian (hex elements) --needs check for other elements 00200 // 00201 MBERRCHK(mb->create_vertex(coords_new_quad, new_vert[nid]), mb); 00202 if (debug){ 00203 m_LogFile << std::setprecision (3) << std::scientific << " : created node:" << (nid + 1) 00204 << " of " << new_vert.size() << " new nodes:- coords: " << coords_new_quad[0] 00205 << ", " << coords_new_quad[1] << ", " << coords_new_quad[2] << std::endl; 00206 } 00207 } 00208 } 00209 00210 return 0; 00211 } 00212 00213 void PostBL::get_normal_quad (std::vector<EntityHandle>conn, moab::CartVect &v) 00214 // --------------------------------------------------------------------------- 00218 // --------------------------------------------------------------------------- 00219 { 00220 moab::CartVect coords[3]; 00221 MBERRCHK(mb->get_coords(&conn[0], 3, (double*) &coords[0]), mb); 00222 moab::CartVect AB(coords[1] - coords[0]); 00223 moab::CartVect BC(coords[2] - coords[1]); 00224 moab::CartVect normal = AB*BC; 00225 normal = normal/normal.length(); 00226 v = normal; 00227 } 00228 00229 void PostBL::get_normal_edge (std::vector<EntityHandle>conn, moab::CartVect BC, moab::CartVect &v) 00230 // --------------------------------------------------------------------------- 00234 // --------------------------------------------------------------------------- 00235 { 00236 moab::CartVect coords[2]; 00237 MBERRCHK(mb->get_coords(&conn[0], 2, (double*) &coords[0]), mb); 00238 moab::CartVect AB(coords[1] - coords[0]); 00239 moab::CartVect normal = AB*BC; 00240 normal = normal/normal.length(); 00241 v = normal; 00242 } 00243 00244 void PostBL::get_det_jacobian(std::vector<moab::EntityHandle> conn, int offset, double &AvgJ) 00245 // --------------------------------------------------------------------------- 00249 // --------------------------------------------------------------------------- 00250 { 00251 //TODO: Add quality check for tri/quad and pyramids 00252 if(m_Conn ==8){ 00253 ++m_JacCalls; 00254 moab::CartVect vertex[8], xi; 00255 mstream m_LogFile; 00256 MBERRCHK(mb->get_coords(&conn[offset], 8, (double*) &vertex[0]), mb); 00257 00258 double corner[8][3] = { { -1, -1, -1 }, 00259 { 1, -1, -1 }, 00260 { 1, 1, -1 }, 00261 { -1, 1, -1 }, 00262 { -1, -1, 1 }, 00263 { 1, -1, 1 }, 00264 { 1, 1, 1 }, 00265 { -1, 1, 1 } }; 00266 00267 for (unsigned j = 0; j < 8; ++j) { 00268 xi[0] = corner[j][0]; 00269 xi[1] = corner[j][1]; 00270 xi[2] = corner[j][2]; 00271 Matrix3 J(0.0); 00272 double detJ = 0; 00273 for (unsigned i = 0; i < 8; ++i) { 00274 const double xi_p = 1 + xi[0]*corner[i][0]; 00275 const double eta_p = 1 + xi[1]*corner[i][1]; 00276 const double zeta_p = 1 + xi[2]*corner[i][2]; 00277 const double dNi_dxi = corner[i][0] * eta_p * zeta_p; 00278 const double dNi_deta = corner[i][1] * xi_p * zeta_p; 00279 const double dNi_dzeta = corner[i][2] * xi_p * eta_p; 00280 J(0,0) += dNi_dxi * vertex[i][0]; 00281 J(1,0) += dNi_dxi * vertex[i][1]; 00282 J(2,0) += dNi_dxi * vertex[i][2]; 00283 J(0,1) += dNi_deta * vertex[i][0]; 00284 J(1,1) += dNi_deta * vertex[i][1]; 00285 J(2,1) += dNi_deta * vertex[i][2]; 00286 J(0,2) += dNi_dzeta * vertex[i][0]; 00287 J(1,2) += dNi_dzeta * vertex[i][1]; 00288 J(2,2) += dNi_dzeta * vertex[i][2]; 00289 } 00290 J *= 0.125; 00291 detJ = J.determinant(); 00292 if(detJ <= 0.0){ 00293 m_LogFile << "We've negative jacobian at the hex corner: "<< j+1 << std::endl; 00294 exit(0); 00295 } 00296 AvgJ+=detJ; 00297 } 00298 AvgJ/=8; 00299 if(m_JacCalls == 1){ 00300 m_JLo = AvgJ; 00301 m_JHi = AvgJ; 00302 } 00303 else if(AvgJ < m_JLo){ 00304 m_JLo = AvgJ; 00305 } 00306 else if(AvgJ > m_JHi){ 00307 m_JHi = AvgJ; 00308 } 00309 } 00310 } 00311 void PostBL::find_min_edge_length(moab::Range adj_qn, moab::EntityHandle node , moab::Range bl_nodes, double &e_len) 00312 // --------------------------------------------------------------------------- 00316 // --------------------------------------------------------------------------- 00317 { 00318 // get nodes adj(a) to BL node 00319 moab::CartVect coords[1]; 00320 double len = 0; e_len = 0; 00321 MBERRCHK(mb->get_coords(&node, 1, (double*) &coords[0]), mb); 00322 moab::Range non_bl; 00323 non_bl = subtract(adj_qn, bl_nodes); 00324 int n_non_bl = (int) non_bl.size(); 00325 00326 // if there are no bl nodes - this case has already been dealt with 00327 if(non_bl.size() > 0){ 00328 moab::CartVect non_bl_coords[4]; 00329 for(int i=0; i< n_non_bl; i++){ 00330 MBERRCHK(mb->get_coords(non_bl,(double*) &non_bl_coords[0]), mb); 00331 moab::CartVect edge(coords[0] - non_bl_coords[0]); 00332 len = edge.length(); 00333 if(i==0) 00334 e_len = len; 00335 if(len < e_len) 00336 e_len = len; 00337 } 00338 } 00339 // m_LogFile << " node minimum edge length" << e_len << std::endl; 00340 } 00341 }