MeshKit  1.0
computenormals.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines