MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2007 Sandia National Laboratories. Developed at the 00005 University of Wisconsin--Madison under SNL contract number 00006 624796. The U.S. Government and the University of Wisconsin 00007 retain certain rights to this software. 00008 00009 This library is free software; you can redistribute it and/or 00010 modify it under the terms of the GNU Lesser General Public 00011 License as published by the Free Software Foundation; either 00012 version 2.1 of the License, or (at your option) any later version. 00013 00014 This library is distributed in the hope that it will be useful, 00015 but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 Lesser General Public License for more details. 00018 00019 You should have received a copy of the GNU Lesser General Public License 00020 (lgpl.txt) along with this library; if not, write to the Free Software 00021 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 00023 (2009) kraftche@cae.wisc.edu 00024 00025 ***************************************************************** */ 00026 00027 /** \file SlaveBoundaryVertices.cpp 00028 * \brief 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #include "SlaveBoundaryVertices.hpp" 00033 #include "Settings.hpp" 00034 #include "MsqError.hpp" 00035 #include "MeshInterface.hpp" 00036 #include "MsqVertex.hpp" 00037 00038 #include <vector> 00039 #include <algorithm> 00040 00041 namespace MBMesquite 00042 { 00043 00044 SlaveBoundaryVertices::SlaveBoundaryVertices( unsigned depth, unsigned dim ) : elemDepth( depth ), domainDoF( dim ) {} 00045 00046 std::string SlaveBoundaryVertices::get_name() const 00047 { 00048 return "SlaveBoundaryVertices"; 00049 } 00050 00051 double SlaveBoundaryVertices::loop_over_mesh( MeshDomainAssoc* mesh_and_domain, const Settings* settings, 00052 MsqError& err ) 00053 { 00054 Mesh* mesh = mesh_and_domain->get_mesh(); 00055 MeshDomain* domain = mesh_and_domain->get_domain(); 00056 00057 if( settings->get_slaved_ho_node_mode() != Settings::SLAVE_CALCULATED ) 00058 { 00059 MSQ_SETERR( err ) 00060 ( "Request to calculate higher-order node slaved status " 00061 "when Settings::get_get_slaved_ho_node_mode() " 00062 "!= SLAVE_CALCUALTED", 00063 MsqError::INVALID_STATE ); 00064 return 0.0; 00065 } 00066 00067 // If user said that we should treat fixed vertices as the 00068 // boundary, but specified that fixed vertices are defined 00069 // by the dimension of their geometric domain, then just 00070 // do distance from the geometric domain. 00071 int dim = this->domainDoF; 00072 if( dim >= 4 && settings->get_fixed_vertex_mode() != Settings::FIXED_FLAG ) dim = settings->get_fixed_vertex_mode(); 00073 00074 // Create a map to contain vertex depth. Intiliaze all to 00075 // elemDepth+1. 00076 std::vector< Mesh::VertexHandle > vertices; 00077 mesh->get_all_vertices( vertices, err ); 00078 MSQ_ERRZERO( err ); 00079 if( vertices.empty() ) return 0.0; 00080 std::sort( vertices.begin(), vertices.end() ); 00081 std::vector< unsigned short > depth( vertices.size(), elemDepth + 1 ); 00082 std::vector< bool > fixed; 00083 mesh->vertices_get_fixed_flag( arrptr( vertices ), fixed, vertices.size(), err ); 00084 MSQ_ERRZERO( err ); 00085 00086 // Initialize map with boundary vertices. 00087 if( dim >= 4 ) 00088 { 00089 for( size_t i = 0; i < vertices.size(); ++i ) 00090 if( fixed[i] ) depth[i] = 0; 00091 } 00092 else 00093 { 00094 if( !domain ) 00095 { 00096 MSQ_SETERR( err ) 00097 ( "Request to calculate higher-order node slaved status " 00098 "by distance from bounding domain without a domain.", 00099 MsqError::INVALID_STATE ); 00100 return 0.0; 00101 } 00102 00103 std::vector< unsigned short > dof( vertices.size() ); 00104 domain->domain_DoF( arrptr( vertices ), arrptr( dof ), vertices.size(), err ); 00105 MSQ_ERRZERO( err ); 00106 for( size_t i = 0; i < vertices.size(); ++i ) 00107 if( dof[i] <= dim ) depth[i] = 0; 00108 } 00109 00110 // Now iterate over elements repeatedly until we've found all of the 00111 // elements near the boundary. This could be done much more efficiently 00112 // using vertex-to-element adjacencies, but it is to common for 00113 // applications not to implement that unless doing relaxation smoothing. 00114 // This is O(elemDepth * elements.size() * ln(vertices.size())); 00115 std::vector< Mesh::ElementHandle > elements; 00116 std::vector< Mesh::ElementHandle >::const_iterator j, k; 00117 std::vector< Mesh::VertexHandle > conn; 00118 std::vector< size_t > junk( 2 ); 00119 mesh->get_all_elements( elements, err ); 00120 MSQ_ERRZERO( err ); 00121 if( elements.empty() ) return 0.0; 00122 bool some_changed; 00123 do 00124 { 00125 some_changed = false; 00126 for( j = elements.begin(); j != elements.end(); ++j ) 00127 { 00128 conn.clear(); 00129 junk.clear(); 00130 mesh->elements_get_attached_vertices( &*j, 1, conn, junk, err ); 00131 MSQ_ERRZERO( err ); 00132 unsigned short elem_depth = elemDepth + 1; 00133 for( k = conn.begin(); k != conn.end(); ++k ) 00134 { 00135 size_t i = std::lower_bound( vertices.begin(), vertices.end(), *k ) - vertices.begin(); 00136 if( i == vertices.size() ) 00137 { 00138 MSQ_SETERR( err ) 00139 ( "Invalid vertex handle in element connectivity list.", MsqError::INVALID_MESH ); 00140 return 0.0; 00141 } 00142 if( depth[i] < elem_depth ) elem_depth = depth[i]; 00143 } 00144 if( elem_depth == elemDepth + 1 ) continue; 00145 00146 ++elem_depth; 00147 for( k = conn.begin(); k != conn.end(); ++k ) 00148 { 00149 size_t i = std::lower_bound( vertices.begin(), vertices.end(), *k ) - vertices.begin(); 00150 if( depth[i] > elem_depth ) 00151 { 00152 depth[i] = elem_depth; 00153 some_changed = true; 00154 } 00155 } 00156 } // for(elements) 00157 } while( some_changed ); 00158 00159 // Now remove any corner vertices from the slaved set 00160 std::vector< Mesh::VertexHandle >::iterator p; 00161 std::vector< EntityTopology > types( elements.size() ); 00162 mesh->elements_get_topologies( arrptr( elements ), arrptr( types ), elements.size(), err ); 00163 MSQ_ERRZERO( err ); 00164 for( j = elements.begin(); j != elements.end(); ++j ) 00165 { 00166 const unsigned corners = TopologyInfo::corners( types[j - elements.begin()] ); 00167 conn.clear(); 00168 junk.clear(); 00169 mesh->elements_get_attached_vertices( &*j, 1, conn, junk, err ); 00170 MSQ_ERRZERO( err ); 00171 for( unsigned i = 0; i < corners; ++i ) 00172 { 00173 p = std::lower_bound( vertices.begin(), vertices.end(), conn[i] ); 00174 depth[p - vertices.begin()] = 0; 00175 } 00176 } 00177 00178 // Now mark all vertices *not* within specified depth as slave vertices. 00179 std::vector< unsigned char > bytes( vertices.size() ); 00180 mesh->vertices_get_byte( arrptr( vertices ), arrptr( bytes ), vertices.size(), err ); 00181 MSQ_ERRZERO( err ); 00182 for( size_t i = 0; i < vertices.size(); ++i ) 00183 { 00184 if( depth[i] <= elemDepth || fixed[i] ) 00185 bytes[i] &= ~MsqVertex::MSQ_DEPENDENT; 00186 else 00187 bytes[i] |= MsqVertex::MSQ_DEPENDENT; 00188 } 00189 00190 mesh->vertices_set_byte( arrptr( vertices ), arrptr( bytes ), vertices.size(), err ); 00191 MSQ_ERRZERO( err ); 00192 return 0.0; 00193 } 00194 00195 void SlaveBoundaryVertices::initialize_queue( MeshDomainAssoc*, const Settings*, MsqError& ) {} 00196 00197 } // namespace MBMesquite