face_quad9.C
Go to the documentation of this file.
00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 00003 00004 // This library is free software; you can redistribute it and/or 00005 // modify it under the terms of the GNU Lesser General Public 00006 // License as published by the Free Software Foundation; either 00007 // version 2.1 of the License, or (at your option) any later version. 00008 00009 // This library is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // Lesser General Public License for more details. 00013 00014 // You should have received a copy of the GNU Lesser General Public 00015 // License along with this library; if not, write to the Free Software 00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00017 00018 // C++ includes 00019 00020 // Local includes 00021 #include "libmesh/side.h" 00022 #include "libmesh/edge_edge3.h" 00023 #include "libmesh/face_quad9.h" 00024 00025 namespace libMesh 00026 { 00027 00028 00029 00030 00031 // ------------------------------------------------------------ 00032 // Quad9 class static member initializations 00033 const unsigned int Quad9::side_nodes_map[4][3] = 00034 { 00035 {0, 1, 4}, // Side 0 00036 {1, 2, 5}, // Side 1 00037 {2, 3, 6}, // Side 2 00038 {3, 0, 7} // Side 3 00039 }; 00040 00041 00042 #ifdef LIBMESH_ENABLE_AMR 00043 00044 const float Quad9::_embedding_matrix[4][9][9] = 00045 { 00046 // embedding matrix for child 0 00047 { 00048 // 0 1 2 3 4 5 6 7 8 00049 { 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 0 00050 { 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 1 00051 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 2 00052 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000 }, // 3 00053 { 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 4 00054 { 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.00000, 0.750000 }, // 5 00055 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.750000 }, // 6 00056 { 0.375000, 0.00000, 0.00000, -0.125000, 0.00000, 0.00000, 0.00000, 0.750000, 0.00000 }, // 7 00057 { 0.140625, -0.0468750, 0.0156250, -0.0468750, 0.281250, -0.0937500, -0.0937500, 0.281250, 0.562500 } // 8 00058 }, 00059 00060 // embedding matrix for child 1 00061 { 00062 // 0 1 2 3 4 5 6 7 8 00063 { 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 0 00064 { 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 1 00065 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000 }, // 2 00066 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 3 00067 { -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 4 00068 { 0.00000, 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000 }, // 5 00069 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.750000 }, // 6 00070 { 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.00000, 0.750000 }, // 7 00071 { -0.0468750, 0.140625, -0.0468750, 0.0156250, 0.281250, 0.281250, -0.0937500, -0.0937500, 0.562500 } // 8 00072 }, 00073 00074 // embedding matrix for child 2 00075 { 00076 // 0 1 2 3 4 5 6 7 8 00077 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000 }, // 0 00078 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 1 00079 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000 }, // 2 00080 { 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 3 00081 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.750000 }, // 4 00082 { 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.00000, 0.750000 }, // 5 00083 { 0.00000, 0.00000, -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000 }, // 6 00084 { -0.125000, 0.00000, 0.00000, 0.375000, 0.00000, 0.00000, 0.00000, 0.750000, 0.00000 }, // 7 00085 { -0.0468750, 0.0156250, -0.0468750, 0.140625, -0.0937500, -0.0937500, 0.281250, 0.281250, 0.562500 } // 8 00086 }, 00087 00088 // embedding matrix for child 3 00089 { 00090 // 0 1 2 3 4 5 6 7 8 00091 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 0 00092 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000 }, // 1 00093 { 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 2 00094 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000 }, // 3 00095 { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.750000 }, // 4 00096 { 0.00000, -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000 }, // 5 00097 { 0.00000, 0.00000, 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000 }, // 6 00098 { 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.00000, 0.750000 }, // 7 00099 { 0.0156250, -0.0468750, 0.140625, -0.0468750, -0.0937500, 0.281250, 0.281250, -0.0937500, 0.562500 } // 8 00100 } 00101 }; 00102 00103 #endif 00104 00105 00106 00107 // ------------------------------------------------------------ 00108 // Quad9 class member functions 00109 00110 bool Quad9::is_vertex(const unsigned int i) const 00111 { 00112 if (i < 4) 00113 return true; 00114 return false; 00115 } 00116 00117 bool Quad9::is_edge(const unsigned int i) const 00118 { 00119 if (i < 4) 00120 return false; 00121 if (i > 7) 00122 return false; 00123 return true; 00124 } 00125 00126 bool Quad9::is_face(const unsigned int i) const 00127 { 00128 if (i > 7) 00129 return true; 00130 return false; 00131 } 00132 00133 bool Quad9::is_node_on_side(const unsigned int n, 00134 const unsigned int s) const 00135 { 00136 libmesh_assert_less (s, n_sides()); 00137 for (unsigned int i = 0; i != 3; ++i) 00138 if (side_nodes_map[s][i] == n) 00139 return true; 00140 return false; 00141 } 00142 00143 00144 00145 bool Quad9::has_affine_map() const 00146 { 00147 // make sure corners form a parallelogram 00148 Point v = this->point(1) - this->point(0); 00149 if (!v.relative_fuzzy_equals(this->point(2) - this->point(3))) 00150 return false; 00151 // make sure "horizontal" sides are straight 00152 v /= 2; 00153 if (!v.relative_fuzzy_equals(this->point(4) - this->point(0)) || 00154 !v.relative_fuzzy_equals(this->point(6) - this->point(3))) 00155 return false; 00156 // make sure "vertical" sides are straight 00157 // and the center node is centered 00158 v = (this->point(3) - this->point(0))/2; 00159 if (!v.relative_fuzzy_equals(this->point(7) - this->point(0)) || 00160 !v.relative_fuzzy_equals(this->point(5) - this->point(1)) || 00161 !v.relative_fuzzy_equals(this->point(8) - this->point(4))) 00162 return false; 00163 return true; 00164 } 00165 00166 00167 00168 dof_id_type Quad9::key (const unsigned int s) const 00169 { 00170 libmesh_assert_less (s, this->n_sides()); 00171 00172 switch (s) 00173 { 00174 case 0: 00175 00176 return 00177 this->compute_key (this->node(4)); 00178 00179 case 1: 00180 00181 return 00182 this->compute_key (this->node(5)); 00183 00184 case 2: 00185 00186 return 00187 this->compute_key (this->node(6)); 00188 00189 case 3: 00190 00191 return 00192 this->compute_key (this->node(7)); 00193 00194 default: 00195 libmesh_error_msg("Invalid side s = " << s); 00196 } 00197 00198 libmesh_error_msg("We'll never get here!"); 00199 return 0; 00200 } 00201 00202 00203 00204 UniquePtr<Elem> Quad9::build_side (const unsigned int i, 00205 bool proxy) const 00206 { 00207 libmesh_assert_less (i, this->n_sides()); 00208 00209 if (proxy) 00210 return UniquePtr<Elem>(new Side<Edge3,Quad9>(this,i)); 00211 00212 else 00213 { 00214 Elem* edge = new Edge3; 00215 edge->subdomain_id() = this->subdomain_id(); 00216 00217 switch (i) 00218 { 00219 case 0: 00220 { 00221 edge->set_node(0) = this->get_node(0); 00222 edge->set_node(1) = this->get_node(1); 00223 edge->set_node(2) = this->get_node(4); 00224 break; 00225 } 00226 case 1: 00227 { 00228 edge->set_node(0) = this->get_node(1); 00229 edge->set_node(1) = this->get_node(2); 00230 edge->set_node(2) = this->get_node(5); 00231 break; 00232 } 00233 case 2: 00234 { 00235 edge->set_node(0) = this->get_node(2); 00236 edge->set_node(1) = this->get_node(3); 00237 edge->set_node(2) = this->get_node(6); 00238 break; 00239 } 00240 case 3: 00241 { 00242 edge->set_node(0) = this->get_node(3); 00243 edge->set_node(1) = this->get_node(0); 00244 edge->set_node(2) = this->get_node(7); 00245 break; 00246 } 00247 default: 00248 libmesh_error_msg("Invalid side i = " << i); 00249 } 00250 00251 return UniquePtr<Elem>(edge); 00252 } 00253 00254 libmesh_error_msg("We'll never get here!"); 00255 return UniquePtr<Elem>(); 00256 } 00257 00258 00259 00260 00261 00262 00263 00264 void Quad9::connectivity(const unsigned int sf, 00265 const IOPackage iop, 00266 std::vector<dof_id_type>& conn) const 00267 { 00268 libmesh_assert_less (sf, this->n_sub_elem()); 00269 libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE); 00270 00271 conn.resize(4); 00272 00273 switch (iop) 00274 { 00275 case TECPLOT: 00276 { 00277 switch(sf) 00278 { 00279 case 0: 00280 // linear sub-quad 0 00281 conn[0] = this->node(0)+1; 00282 conn[1] = this->node(4)+1; 00283 conn[2] = this->node(8)+1; 00284 conn[3] = this->node(7)+1; 00285 return; 00286 00287 case 1: 00288 // linear sub-quad 1 00289 conn[0] = this->node(4)+1; 00290 conn[1] = this->node(1)+1; 00291 conn[2] = this->node(5)+1; 00292 conn[3] = this->node(8)+1; 00293 return; 00294 00295 case 2: 00296 // linear sub-quad 2 00297 conn[0] = this->node(7)+1; 00298 conn[1] = this->node(8)+1; 00299 conn[2] = this->node(6)+1; 00300 conn[3] = this->node(3)+1; 00301 return; 00302 00303 case 3: 00304 // linear sub-quad 3 00305 conn[0] = this->node(8)+1; 00306 conn[1] = this->node(5)+1; 00307 conn[2] = this->node(2)+1; 00308 conn[3] = this->node(6)+1; 00309 return; 00310 00311 default: 00312 libmesh_error_msg("Invalid sf = " << sf); 00313 } 00314 } 00315 00316 case VTK: 00317 { 00318 conn.resize(9); 00319 conn[0] = this->node(0); 00320 conn[1] = this->node(1); 00321 conn[2] = this->node(2); 00322 conn[3] = this->node(3); 00323 conn[4] = this->node(4); 00324 conn[5] = this->node(5); 00325 conn[6] = this->node(6); 00326 conn[7] = this->node(7); 00327 conn[8] = this->node(8); 00328 return; 00329 00330 /* 00331 switch(sf) 00332 { 00333 case 0: 00334 // linear sub-quad 0 00335 conn[0] = this->node(0); 00336 conn[1] = this->node(4); 00337 conn[2] = this->node(8); 00338 conn[3] = this->node(7); 00339 00340 return; 00341 00342 case 1: 00343 // linear sub-quad 1 00344 conn[0] = this->node(4); 00345 conn[1] = this->node(1); 00346 conn[2] = this->node(5); 00347 conn[3] = this->node(8); 00348 00349 return; 00350 00351 case 2: 00352 // linear sub-quad 2 00353 conn[0] = this->node(7); 00354 conn[1] = this->node(8); 00355 conn[2] = this->node(6); 00356 conn[3] = this->node(3); 00357 00358 return; 00359 00360 case 3: 00361 // linear sub-quad 3 00362 conn[0] = this->node(8); 00363 conn[1] = this->node(5); 00364 conn[2] = this->node(2); 00365 conn[3] = this->node(6); 00366 00367 return; 00368 00369 default: 00370 libmesh_error_msg("Invalid sf = " << sf); 00371 }*/ 00372 } 00373 00374 default: 00375 libmesh_error_msg("Unsupported IO package " << iop); 00376 } 00377 } 00378 00379 00380 00381 00382 unsigned int Quad9::n_second_order_adjacent_vertices (const unsigned int n) const 00383 { 00384 switch (n) 00385 { 00386 case 4: 00387 case 5: 00388 case 6: 00389 case 7: 00390 return 2; 00391 00392 case 8: 00393 return 4; 00394 00395 default: 00396 libmesh_error_msg("Invalid n = " << n); 00397 } 00398 00399 libmesh_error_msg("We'll never get here!"); 00400 return libMesh::invalid_uint; 00401 } 00402 00403 00404 00405 unsigned short int Quad9::second_order_adjacent_vertex (const unsigned int n, 00406 const unsigned int v) const 00407 { 00408 libmesh_assert_greater_equal (n, this->n_vertices()); 00409 libmesh_assert_less (n, this->n_nodes()); 00410 00411 switch (n) 00412 { 00413 case 8: 00414 { 00415 libmesh_assert_less (v, 4); 00416 return static_cast<unsigned short int>(v); 00417 } 00418 00419 default: 00420 { 00421 libmesh_assert_less (v, 2); 00422 // use the matrix that we inherited from \p Quad 00423 return _second_order_adjacent_vertices[n-this->n_vertices()][v]; 00424 } 00425 } 00426 } 00427 00428 00429 00430 std::pair<unsigned short int, unsigned short int> 00431 Quad9::second_order_child_vertex (const unsigned int n) const 00432 { 00433 libmesh_assert_greater_equal (n, this->n_vertices()); 00434 libmesh_assert_less (n, this->n_nodes()); 00435 /* 00436 * the _second_order_vertex_child_* vectors are 00437 * stored in face_quad.C, since they are identical 00438 * for Quad8 and Quad9 (for the first 4 higher-order nodes) 00439 */ 00440 return std::pair<unsigned short int, unsigned short int> 00441 (_second_order_vertex_child_number[n], 00442 _second_order_vertex_child_index[n]); 00443 } 00444 00445 } // namespace libMesh
generated by
