|
ViennaGrid - The Vienna Grid Library
2.1.0
|
00001 #ifndef VIENNAGRID_ALGORITHM_DETAIL_REFINE_TET_HPP 00002 #define VIENNAGRID_ALGORITHM_DETAIL_REFINE_TET_HPP 00003 00004 /* ======================================================================= 00005 Copyright (c) 2011-2014, Institute for Microelectronics, 00006 Institute for Analysis and Scientific Computing, 00007 TU Wien. 00008 00009 ----------------- 00010 ViennaGrid - The Vienna Grid Library 00011 ----------------- 00012 00013 License: MIT (X11), see file LICENSE in the base directory 00014 ======================================================================= */ 00015 00016 #include "viennagrid/forwards.hpp" 00017 #include "viennagrid/mesh/segmentation.hpp" 00018 #include "viennagrid/topology/vertex.hpp" 00019 #include "viennagrid/topology/line.hpp" 00020 #include "viennagrid/topology/simplex.hpp" 00021 #include "viennagrid/algorithm/norm.hpp" 00022 00027 namespace viennagrid 00028 { 00029 namespace detail 00030 { 00031 00037 template<typename MeshT, typename VertexHandleT> 00038 bool stable_line_is_longer(MeshT const & mesh, 00039 VertexHandleT vh1_1, VertexHandleT vh1_2, 00040 VertexHandleT vh2_1, VertexHandleT vh2_2) 00041 { 00042 typedef typename viennagrid::detail::result_of::value_type< VertexHandleT >::type VertexType; 00043 typedef typename viennagrid::result_of::point< MeshT >::type PointType; 00044 typedef typename viennagrid::result_of::coord< PointType >::type ScalarType; 00045 00046 const VertexType & v1_1 = viennagrid::dereference_handle( mesh, vh1_1 ); 00047 const VertexType & v1_2 = viennagrid::dereference_handle( mesh, vh1_2 ); 00048 const VertexType & v2_1 = viennagrid::dereference_handle( mesh, vh2_1 ); 00049 const VertexType & v2_2 = viennagrid::dereference_handle( mesh, vh2_2 ); 00050 00051 const VertexType & v1_1_ptr = (v1_1.id() < v1_2.id()) ? v1_1 : v1_2; //v1_1 carries smaller ID 00052 const VertexType & v1_2_ptr = (v1_1.id() < v1_2.id()) ? v1_2 : v1_1; //v1_2 carries larger ID 00053 00054 const VertexType & v2_1_ptr = (v2_1.id() < v2_2.id()) ? v2_1 : v2_2; //v2_1 carries smaller ID 00055 const VertexType & v2_2_ptr = (v2_1.id() < v2_2.id()) ? v2_2 : v2_1; //v2_2 carries larger ID 00056 00057 ScalarType line1 = viennagrid::norm( viennagrid::point(mesh, v1_1) - viennagrid::point(mesh, v1_2) ); 00058 ScalarType line2 = viennagrid::norm( viennagrid::point(mesh, v2_1) - viennagrid::point(mesh, v2_2) ); 00059 00060 00061 if (line1 > line2) 00062 { 00063 return true; 00064 } 00065 else if (line1 >= line2) // use of == leads to floating-point comparison warning in Clang 00066 { 00067 //compare IDs: 00068 if (v1_1_ptr.id() > v2_1_ptr.id()) 00069 { 00070 return true; 00071 } 00072 else if (v1_1_ptr.id() == v2_1_ptr.id()) 00073 { 00074 if (v1_2_ptr.id() > v2_2_ptr.id()) 00075 return true; 00076 else if (v1_2_ptr.id() == v2_2_ptr.id()) 00077 return false; //identical lines are compared! 00078 } 00079 } 00080 00081 return false; 00082 } 00083 00084 template<typename MeshT, typename VertexHandleContainer> 00085 bool stable_line_is_longer(MeshT const & mesh, VertexHandleContainer vertices, 00086 unsigned int i0, unsigned int i1, unsigned int i2, unsigned int i3) 00087 { 00088 return stable_line_is_longer(mesh, 00089 *viennagrid::advance(vertices.begin(), i0), *viennagrid::advance(vertices.begin(), i1), 00090 *viennagrid::advance(vertices.begin(), i2), *viennagrid::advance(vertices.begin(), i3)); 00091 } 00092 00093 00094 template<typename ElementsVerticesHandleContainerT, typename VertexHandleContainerT> 00095 void add_refinement_element(ElementsVerticesHandleContainerT & elements_vertices, 00096 VertexHandleContainerT vertex_handle_container, 00097 unsigned int i0, unsigned int i1, unsigned int i2, unsigned int i3) 00098 { 00099 elements_vertices.resize( elements_vertices.size()+1 ); 00100 00101 elements_vertices.back().resize(4); 00102 elements_vertices.back()[0] = *viennagrid::advance(vertex_handle_container.begin(), i0); 00103 elements_vertices.back()[1] = *viennagrid::advance(vertex_handle_container.begin(), i1); 00104 elements_vertices.back()[2] = *viennagrid::advance(vertex_handle_container.begin(), i2); 00105 elements_vertices.back()[3] = *viennagrid::advance(vertex_handle_container.begin(), i3); 00106 } 00107 00108 00109 00110 00116 template<> 00117 struct element_refinement<tetrahedron_tag> 00118 { 00119 00121 template<typename ElementT, typename MeshT, 00122 typename ElementsVerticesHandleContainerT, typename VertexCopyMapT, 00123 typename EdgeRefinementFlagAccessor, typename EdgeToVertexHandleAccessor> 00124 static void apply0(ElementT const & element_in, MeshT const &, 00125 ElementsVerticesHandleContainerT & elements_vertices, 00126 VertexCopyMapT & vertex_copy_map_, 00127 EdgeRefinementFlagAccessor const &, EdgeToVertexHandleAccessor const &) 00128 { 00129 typedef typename viennagrid::result_of::const_element_range<ElementT, vertex_tag>::type VertexOnCellRange; 00130 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator; 00131 00132 typedef typename viennagrid::result_of::handle<MeshT, vertex_tag>::type VertexHandleType; 00133 00134 static_array< VertexHandleType, boundary_elements<tetrahedron_tag, vertex_tag>::num > vertex_handles; 00135 00136 // Step 1: grab existing vertices: 00137 VertexOnCellRange vertices_on_cell = viennagrid::elements<vertex_tag>(element_in); 00138 VertexOnCellIterator vocit = vertices_on_cell.begin(); 00139 00140 vertex_handles[0] = vertex_copy_map_(*vocit); ++vocit; 00141 vertex_handles[1] = vertex_copy_map_(*vocit); ++vocit; 00142 vertex_handles[2] = vertex_copy_map_(*vocit); ++vocit; 00143 vertex_handles[3] = vertex_copy_map_(*vocit); 00144 00145 // // Step 2: Add new cells to new mesh: 00146 add_refinement_element(elements_vertices, vertex_handles, 0, 1, 2, 3); 00147 } 00148 00150 template<typename ElementT, typename MeshT, 00151 typename ElementsVerticesHandleContainerT, typename VertexCopyMapT, 00152 typename EdgeRefinementFlagAccessor, typename EdgeToVertexHandleAccessor> 00153 static void apply1(ElementT const & element_in, MeshT const &, 00154 ElementsVerticesHandleContainerT & elements_vertices, 00155 VertexCopyMapT & vertex_copy_map_, 00156 EdgeRefinementFlagAccessor const & edge_refinement_flag_accessor, EdgeToVertexHandleAccessor const & edge_to_vertex_handle_accessor) 00157 { 00158 00159 typedef typename viennagrid::result_of::const_element_range<ElementT, vertex_tag>::type VertexOnCellRange; 00160 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator; 00161 typedef typename viennagrid::result_of::const_element_range<ElementT, line_tag>::type EdgeOnCellRange; 00162 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator; 00163 00164 typedef typename viennagrid::result_of::element<ElementT, line_tag>::type EdgeType; 00165 00166 typedef typename viennagrid::result_of::handle<MeshT, vertex_tag>::type VertexHandleType; 00167 00168 static_array< VertexHandleType, boundary_elements<tetrahedron_tag, vertex_tag>::num > vertices; 00169 00170 00171 // 00172 // Step 1: Get vertices from input cell 00173 // 00174 VertexOnCellRange vertices_on_cell(element_in); 00175 VertexOnCellIterator vocit = vertices_on_cell.begin(); 00176 00177 00178 vertices[0] = vertex_copy_map_(*vocit); ++vocit; 00179 vertices[1] = vertex_copy_map_(*vocit); ++vocit; 00180 vertices[2] = vertex_copy_map_(*vocit); ++vocit; 00181 vertices[3] = vertex_copy_map_(*vocit); 00182 00183 // 00184 // Step 2: Bring vertices in correct order, such that refined edge is on {0,1}-edge 00185 // 00186 00187 static_array< VertexHandleType, boundary_elements<tetrahedron_tag, vertex_tag>::num + 1 > ordered_vertices; 00188 00189 EdgeOnCellRange edges_on_cell(element_in); 00190 EdgeOnCellIterator eocit = edges_on_cell.begin(); 00191 EdgeType const & e0 = *eocit; ++eocit; 00192 EdgeType const & e1 = *eocit; ++eocit; 00193 EdgeType const & e2 = *eocit; ++eocit; 00194 EdgeType const & e3 = *eocit; ++eocit; 00195 EdgeType const & e4 = *eocit; ++eocit; 00196 EdgeType const & e5 = *eocit; 00197 00198 if (edge_refinement_flag_accessor(e0) == true) 00199 { 00200 ordered_vertices[0] = vertices[0]; 00201 ordered_vertices[1] = vertices[1]; 00202 ordered_vertices[2] = vertices[2]; 00203 ordered_vertices[3] = vertices[3]; 00204 ordered_vertices[4] = edge_to_vertex_handle_accessor(e0); 00205 } 00206 else if (edge_refinement_flag_accessor(e1) == true) 00207 { 00208 ordered_vertices[0] = vertices[2]; 00209 ordered_vertices[1] = vertices[0]; 00210 ordered_vertices[2] = vertices[1]; 00211 ordered_vertices[3] = vertices[3]; 00212 ordered_vertices[4] = edge_to_vertex_handle_accessor(e1); 00213 } 00214 else if (edge_refinement_flag_accessor(e2) == true) 00215 { 00216 ordered_vertices[0] = vertices[0]; 00217 ordered_vertices[1] = vertices[3]; 00218 ordered_vertices[2] = vertices[1]; 00219 ordered_vertices[3] = vertices[2]; 00220 ordered_vertices[4] = edge_to_vertex_handle_accessor(e2); 00221 } 00222 else if (edge_refinement_flag_accessor(e3) == true) 00223 { 00224 ordered_vertices[0] = vertices[1]; 00225 ordered_vertices[1] = vertices[2]; 00226 ordered_vertices[2] = vertices[0]; 00227 ordered_vertices[3] = vertices[3]; 00228 ordered_vertices[4] = edge_to_vertex_handle_accessor(e3); 00229 } 00230 else if (edge_refinement_flag_accessor(e4) == true) 00231 { 00232 ordered_vertices[0] = vertices[3]; 00233 ordered_vertices[1] = vertices[1]; 00234 ordered_vertices[2] = vertices[0]; 00235 ordered_vertices[3] = vertices[2]; 00236 ordered_vertices[4] = edge_to_vertex_handle_accessor(e4); 00237 } 00238 else if (edge_refinement_flag_accessor(e5) == true) 00239 { 00240 ordered_vertices[0] = vertices[3]; 00241 ordered_vertices[1] = vertices[2]; 00242 ordered_vertices[2] = vertices[1]; 00243 ordered_vertices[3] = vertices[0]; 00244 ordered_vertices[4] = edge_to_vertex_handle_accessor(e5); 00245 } 00246 else 00247 { 00248 assert(false && "Logic error: No edge for refinement found!"); 00249 } 00250 00251 // 00252 // Step 3: Write new cells to mesh_out 00253 // 00254 00255 //cell containing vertex 0: 00256 add_refinement_element( elements_vertices, ordered_vertices, 0, 4, 2, 3); 00257 00258 //cell without vertex 0: 00259 add_refinement_element( elements_vertices, ordered_vertices, 4, 1, 2, 3); 00260 } 00261 00262 00264 00265 00266 // 00267 // Refinement of a tetrahedron, bisecting two edges 00268 // 00269 00283 template<typename MeshT, typename ElementsVerticesHandleContainerT, typename VertexHandleIteratorType> 00284 static void apply2_1(MeshT const & mesh, 00285 ElementsVerticesHandleContainerT & element_vertices, 00286 VertexHandleIteratorType vertices) 00287 { 00288 add_refinement_element( element_vertices, vertices, 4, 1, 5, 3); 00289 00290 if (stable_line_is_longer(mesh, vertices, 0, 1, 2, 1)) 00291 { 00292 add_refinement_element( element_vertices, vertices, 0, 4, 2, 3); 00293 add_refinement_element( element_vertices, vertices, 4, 5, 2, 3); 00294 } 00295 else //split edge 12, introduce line 05 00296 { 00297 add_refinement_element( element_vertices, vertices, 0, 4, 5, 3); 00298 add_refinement_element( element_vertices, vertices, 0, 5, 2, 3); 00299 } 00300 } 00301 00315 template<typename MeshT, typename ElementsVerticesHandleContainerT, typename VertexHandleIteratorType> 00316 static void apply2_2(MeshT const &, 00317 ElementsVerticesHandleContainerT & element_vertices, 00318 VertexHandleIteratorType vertices) 00319 { 00320 add_refinement_element( element_vertices, vertices, 4, 1, 2, 5); 00321 add_refinement_element( element_vertices, vertices, 4, 1, 5, 3); 00322 add_refinement_element( element_vertices, vertices, 0, 4, 2, 5); 00323 add_refinement_element( element_vertices, vertices, 0, 4, 5, 3); 00324 } 00325 00329 template<typename ElementType, typename MeshT, 00330 typename ElementsVerticesHandleContainerT, typename VertexCopyMapT, 00331 typename EdgeRefinementFlagAccessor, typename EdgeToVertexHandleAccessor> 00332 static void apply2(ElementType const & element_in, MeshT const & mesh, 00333 ElementsVerticesHandleContainerT & element_vertices, 00334 VertexCopyMapT & vertex_copy_map_, 00335 EdgeRefinementFlagAccessor const & edge_refinement_flag_accessor, EdgeToVertexHandleAccessor const & edge_to_vertex_handle_accessor) 00336 { 00337 typedef typename viennagrid::result_of::const_element_range<ElementType, vertex_tag>::type VertexOnCellRange; 00338 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator; 00339 typedef typename viennagrid::result_of::const_element_range<ElementType, line_tag>::type EdgeOnCellRange; 00340 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator; 00341 00342 typedef typename viennagrid::result_of::handle<ElementType, vertex_tag>::type VertexHandleType; 00343 typedef typename viennagrid::result_of::element<ElementType, line_tag>::type EdgeType; 00344 00345 static_array< VertexHandleType, boundary_elements<tetrahedron_tag, vertex_tag>::num > vertices; 00346 00347 // 00348 // Step 1: Get vertices from input cell 00349 // 00350 VertexOnCellRange vertices_on_cell = viennagrid::elements<vertex_tag>(element_in); 00351 VertexOnCellIterator vocit = vertices_on_cell.begin(); 00352 vertices[0] = vertex_copy_map_(*vocit); ++vocit; 00353 vertices[1] = vertex_copy_map_(*vocit); ++vocit; 00354 vertices[2] = vertex_copy_map_(*vocit); ++vocit; 00355 vertices[3] = vertex_copy_map_(*vocit); 00356 00357 // 00358 // Step 2: Bring vertices in correct order, such that refined edge is on {0,1}-edge 00359 // 00360 00361 static_array< VertexHandleType, boundary_elements<tetrahedron_tag, vertex_tag>::num + 2 > ordered_vertices; 00362 EdgeOnCellRange edges_on_cell = viennagrid::elements<line_tag>(element_in); 00363 EdgeOnCellIterator eocit = edges_on_cell.begin(); 00364 EdgeType const & e0 = *eocit; ++eocit; 00365 EdgeType const & e1 = *eocit; ++eocit; 00366 EdgeType const & e2 = *eocit; ++eocit; 00367 EdgeType const & e3 = *eocit; ++eocit; 00368 EdgeType const & e4 = *eocit; ++eocit; 00369 EdgeType const & e5 = *eocit; 00370 00371 //with e0 00372 if (edge_refinement_flag_accessor(e0) == true) 00373 { 00374 if (edge_refinement_flag_accessor(e1) == true) 00375 { 00376 ordered_vertices[0] = vertices[2]; 00377 ordered_vertices[1] = vertices[0]; 00378 ordered_vertices[2] = vertices[1]; 00379 ordered_vertices[3] = vertices[3]; 00380 ordered_vertices[4] = edge_to_vertex_handle_accessor(e1); 00381 ordered_vertices[5] = edge_to_vertex_handle_accessor(e0); 00382 00383 apply2_1(mesh, element_vertices, ordered_vertices); 00384 } 00385 else if (edge_refinement_flag_accessor(e2) == true) 00386 { 00387 ordered_vertices[0] = vertices[1]; 00388 ordered_vertices[1] = vertices[0]; 00389 ordered_vertices[2] = vertices[3]; 00390 ordered_vertices[3] = vertices[2]; 00391 ordered_vertices[4] = edge_to_vertex_handle_accessor(e0); 00392 ordered_vertices[5] = edge_to_vertex_handle_accessor(e2); 00393 00394 apply2_1(mesh, element_vertices, ordered_vertices); 00395 } 00396 else if (edge_refinement_flag_accessor(e3) == true) 00397 { 00398 ordered_vertices[0] = vertices[0]; 00399 ordered_vertices[1] = vertices[1]; 00400 ordered_vertices[2] = vertices[2]; 00401 ordered_vertices[3] = vertices[3]; 00402 ordered_vertices[4] = edge_to_vertex_handle_accessor(e0); 00403 ordered_vertices[5] = edge_to_vertex_handle_accessor(e3); 00404 00405 apply2_1(mesh, element_vertices, ordered_vertices); 00406 } 00407 else if (edge_refinement_flag_accessor(e4) == true) 00408 { 00409 ordered_vertices[0] = vertices[3]; 00410 ordered_vertices[1] = vertices[1]; 00411 ordered_vertices[2] = vertices[0]; 00412 ordered_vertices[3] = vertices[2]; 00413 ordered_vertices[4] = edge_to_vertex_handle_accessor(e4); 00414 ordered_vertices[5] = edge_to_vertex_handle_accessor(e0); 00415 00416 apply2_1(mesh, element_vertices, ordered_vertices); 00417 } 00418 else if (edge_refinement_flag_accessor(e5) == true) 00419 { 00420 ordered_vertices[0] = vertices[0]; 00421 ordered_vertices[1] = vertices[1]; 00422 ordered_vertices[2] = vertices[2]; 00423 ordered_vertices[3] = vertices[3]; 00424 ordered_vertices[4] = edge_to_vertex_handle_accessor(e0); 00425 ordered_vertices[5] = edge_to_vertex_handle_accessor(e5); 00426 00427 apply2_2(mesh, element_vertices, ordered_vertices); 00428 } 00429 else 00430 { 00431 assert(false && "Logic error: No edge for refinement found!"); 00432 } 00433 } 00434 else if (edge_refinement_flag_accessor(e1) == true) 00435 { 00436 if (edge_refinement_flag_accessor(e2) == true) 00437 { 00438 ordered_vertices[0] = vertices[3]; 00439 ordered_vertices[1] = vertices[0]; 00440 ordered_vertices[2] = vertices[2]; 00441 ordered_vertices[3] = vertices[1]; 00442 ordered_vertices[4] = edge_to_vertex_handle_accessor(e2); 00443 ordered_vertices[5] = edge_to_vertex_handle_accessor(e1); 00444 00445 apply2_1(mesh, element_vertices, ordered_vertices); 00446 } 00447 else if (edge_refinement_flag_accessor(e3) == true) 00448 { 00449 ordered_vertices[0] = vertices[1]; 00450 ordered_vertices[1] = vertices[2]; 00451 ordered_vertices[2] = vertices[0]; 00452 ordered_vertices[3] = vertices[3]; 00453 ordered_vertices[4] = edge_to_vertex_handle_accessor(e3); 00454 ordered_vertices[5] = edge_to_vertex_handle_accessor(e1); 00455 00456 00457 apply2_1(mesh, element_vertices, ordered_vertices); 00458 } 00459 else if (edge_refinement_flag_accessor(e4) == true) 00460 { 00461 ordered_vertices[0] = vertices[2]; 00462 ordered_vertices[1] = vertices[0]; 00463 ordered_vertices[2] = vertices[1]; 00464 ordered_vertices[3] = vertices[3]; 00465 ordered_vertices[4] = edge_to_vertex_handle_accessor(e1); 00466 ordered_vertices[5] = edge_to_vertex_handle_accessor(e4); 00467 00468 apply2_2(mesh, element_vertices, ordered_vertices); 00469 } 00470 else if (edge_refinement_flag_accessor(e5) == true) 00471 { 00472 ordered_vertices[0] = vertices[0]; 00473 ordered_vertices[1] = vertices[2]; 00474 ordered_vertices[2] = vertices[3]; 00475 ordered_vertices[3] = vertices[1]; 00476 ordered_vertices[4] = edge_to_vertex_handle_accessor(e1); 00477 ordered_vertices[5] = edge_to_vertex_handle_accessor(e5); 00478 00479 apply2_1(mesh, element_vertices, ordered_vertices); 00480 } 00481 else 00482 { 00483 assert(false && "Logic error: No edge for refinement found!"); 00484 } 00485 } 00486 else if (edge_refinement_flag_accessor(e2) == true) 00487 { 00488 if (edge_refinement_flag_accessor(e3) == true) 00489 { 00490 ordered_vertices[0] = vertices[3]; 00491 ordered_vertices[1] = vertices[0]; 00492 ordered_vertices[2] = vertices[2]; 00493 ordered_vertices[3] = vertices[1]; 00494 ordered_vertices[4] = edge_to_vertex_handle_accessor(e2); 00495 ordered_vertices[5] = edge_to_vertex_handle_accessor(e3); 00496 00497 apply2_2(mesh, element_vertices, ordered_vertices); 00498 } 00499 else if (edge_refinement_flag_accessor(e4) == true) 00500 { 00501 ordered_vertices[0] = vertices[0]; 00502 ordered_vertices[1] = vertices[3]; 00503 ordered_vertices[2] = vertices[1]; 00504 ordered_vertices[3] = vertices[2]; 00505 ordered_vertices[4] = edge_to_vertex_handle_accessor(e2); 00506 ordered_vertices[5] = edge_to_vertex_handle_accessor(e4); 00507 00508 apply2_1(mesh, element_vertices, ordered_vertices); 00509 } 00510 else if (edge_refinement_flag_accessor(e5) == true) 00511 { 00512 ordered_vertices[0] = vertices[2]; 00513 ordered_vertices[1] = vertices[3]; 00514 ordered_vertices[2] = vertices[0]; 00515 ordered_vertices[3] = vertices[1]; 00516 ordered_vertices[4] = edge_to_vertex_handle_accessor(e5); 00517 ordered_vertices[5] = edge_to_vertex_handle_accessor(e2); 00518 00519 apply2_1(mesh, element_vertices, ordered_vertices); 00520 } 00521 else 00522 { 00523 assert(false && "Logic error: No edge for refinement found!"); 00524 } 00525 } 00526 else if (edge_refinement_flag_accessor(e3) == true) 00527 { 00528 if (edge_refinement_flag_accessor(e4) == true) 00529 { 00530 ordered_vertices[0] = vertices[2]; 00531 ordered_vertices[1] = vertices[1]; 00532 ordered_vertices[2] = vertices[3]; 00533 ordered_vertices[3] = vertices[0]; 00534 ordered_vertices[4] = edge_to_vertex_handle_accessor(e3); 00535 ordered_vertices[5] = edge_to_vertex_handle_accessor(e4); 00536 00537 apply2_1(mesh, element_vertices, ordered_vertices); 00538 } 00539 else if (edge_refinement_flag_accessor(e5) == true) 00540 { 00541 ordered_vertices[0] = vertices[3]; 00542 ordered_vertices[1] = vertices[2]; 00543 ordered_vertices[2] = vertices[1]; 00544 ordered_vertices[3] = vertices[0]; 00545 ordered_vertices[4] = edge_to_vertex_handle_accessor(e5); 00546 ordered_vertices[5] = edge_to_vertex_handle_accessor(e3); 00547 00548 apply2_1(mesh, element_vertices, ordered_vertices); 00549 } 00550 else 00551 { 00552 assert(false && "Logic error: No edge for refinement found!"); 00553 } 00554 } 00555 else if (edge_refinement_flag_accessor(e4) == true) 00556 { 00557 if (edge_refinement_flag_accessor(e5) == true) 00558 { 00559 ordered_vertices[0] = vertices[1]; 00560 ordered_vertices[1] = vertices[3]; 00561 ordered_vertices[2] = vertices[2]; 00562 ordered_vertices[3] = vertices[0]; 00563 ordered_vertices[4] = edge_to_vertex_handle_accessor(e4); 00564 ordered_vertices[5] = edge_to_vertex_handle_accessor(e5); 00565 00566 apply2_1(mesh, element_vertices, ordered_vertices); 00567 } 00568 else 00569 { 00570 assert(false && "Logic error: No edge for refinement found!"); 00571 } 00572 } 00573 else 00574 { 00575 assert(false && "Logic error: No edge for refinement found!"); 00576 } 00577 } 00578 00579 00580 00582 00583 00597 template<typename MeshT, typename ElementsVerticesHandleContainerT, typename VertexHandleIteratorType> 00598 static void apply3_1(MeshT const & mesh, 00599 ElementsVerticesHandleContainerT & element_vertices, 00600 VertexHandleIteratorType vertices) 00601 { 00602 add_refinement_element( element_vertices, vertices, 4, 1, 5, 6); 00603 00604 // Strategy: The two longest edges of the common vertex are split 'twice', 00605 // i.e. two new edges start from the center of the two longest edges 00606 00607 if (stable_line_is_longer(mesh, vertices, 0, 1, 1, 2)) 00608 { 00609 //if (diag13_len > diag12_len) //split edge 13 again, introduce line 62 00610 if (stable_line_is_longer(mesh, vertices, 1, 3, 1, 2)) 00611 { 00612 if (stable_line_is_longer(mesh, vertices, 1, 3, 0, 1)) 00613 { 00614 add_refinement_element( element_vertices, vertices, 0, 6, 2, 3); 00615 add_refinement_element( element_vertices, vertices, 0, 4, 2, 6); 00616 add_refinement_element( element_vertices, vertices, 4, 5, 2, 6); 00617 } 00618 else //split edge 01 again, introduce line 43 00619 { 00620 add_refinement_element( element_vertices, vertices, 0, 4, 2, 3); 00621 add_refinement_element( element_vertices, vertices, 3, 4, 2, 6); 00622 add_refinement_element( element_vertices, vertices, 4, 5, 2, 6); 00623 } 00624 } 00625 else //split edge 12 again, introduce lines 43 and 53 00626 { 00627 if (stable_line_is_longer(mesh, vertices, 1, 3, 0, 1)) 00628 { 00629 assert(false && "diag13_len > diag01_len impossible!"); 00630 } 00631 else //split edge 01 again, introduce line 43 00632 { 00633 add_refinement_element( element_vertices, vertices, 0, 4, 2, 3); 00634 add_refinement_element( element_vertices, vertices, 4, 5, 2, 3); 00635 add_refinement_element( element_vertices, vertices, 4, 6, 5, 3); 00636 } 00637 } 00638 } 00639 else //split edge 12, introduce line 05 00640 { 00641 if (stable_line_is_longer(mesh, vertices, 1, 3, 1, 2)) //split edge 13 again, introduce line 62 00642 { 00643 if (stable_line_is_longer(mesh, vertices, 1, 3, 0, 1)) //split edge 13 again, introduce line 60 00644 { 00645 add_refinement_element( element_vertices, vertices, 0, 4, 5, 6); 00646 add_refinement_element( element_vertices, vertices, 0, 6, 5, 2); 00647 add_refinement_element( element_vertices, vertices, 0, 6, 2, 3); 00648 } 00649 else //split edge 01 again, introduce line 43 00650 { 00651 assert(false && "diag13_len > diag01_len impossible!"); 00652 } 00653 } 00654 else //split edge 12 again, introduce line 53 00655 { 00656 if (stable_line_is_longer(mesh, vertices, 1, 3, 0, 1)) //split edge 13 again, introduce line 60 00657 { 00658 add_refinement_element( element_vertices, vertices, 0, 4, 5, 6); 00659 add_refinement_element( element_vertices, vertices, 0, 6, 5, 3); 00660 add_refinement_element( element_vertices, vertices, 0, 5, 2, 3); 00661 } 00662 else //split edge 01 again, introduce line 43 00663 { 00664 //std::cout << "apply_3_1_4" << std::endl; 00665 add_refinement_element( element_vertices, vertices, 0, 4, 5, 3); 00666 add_refinement_element( element_vertices, vertices, 4, 5, 3, 6); 00667 add_refinement_element( element_vertices, vertices, 0, 5, 2, 3); 00668 } 00669 } 00670 } 00671 00672 } 00673 00674 00689 template<typename MeshT, typename ElementsVerticesHandleContainerT, typename VertexHandleIteratorType> 00690 static void apply3_2(MeshT const &, 00691 ElementsVerticesHandleContainerT & element_vertices, 00692 VertexHandleIteratorType vertices) 00693 { 00694 add_refinement_element( element_vertices, vertices, 0, 4, 6, 3); 00695 add_refinement_element( element_vertices, vertices, 4, 5, 6, 3); 00696 add_refinement_element( element_vertices, vertices, 4, 1, 5, 3); 00697 add_refinement_element( element_vertices, vertices, 6, 5, 2, 3); 00698 } 00699 00700 00714 template<typename MeshT, typename ElementsVerticesHandleContainerT, typename VertexHandleIteratorType> 00715 static void apply3_3(MeshT const & mesh, 00716 ElementsVerticesHandleContainerT & element_vertices, 00717 VertexHandleIteratorType vertices) 00718 { 00719 // Strategy: The two longest edges of the common vertex are split 'twice', 00720 // i.e. two new edges start from the center of the two longest edges 00721 //if (diag01_len > diag03_len) //split edge 01 again, introduce line 43 00722 if (stable_line_is_longer(mesh, vertices, 0, 1, 0, 3)) //split edge 01 again, introduce line 43 00723 { 00724 add_refinement_element( element_vertices, vertices, 4, 1, 5, 3); 00725 00726 //if (diag01_len > diag12_len) //split edge 01 again, introduce line 42 00727 if (stable_line_is_longer(mesh, vertices, 0, 1, 1, 2)) //split edge 01 again, introduce line 42 00728 { 00729 add_refinement_element( element_vertices, vertices, 0, 4, 2, 6); 00730 add_refinement_element( element_vertices, vertices, 6, 4, 2, 3); 00731 add_refinement_element( element_vertices, vertices, 4, 5, 2, 3); 00732 } 00733 else //split edge 12 again, introduce line 50 00734 { 00735 add_refinement_element( element_vertices, vertices, 0, 4, 5, 6); 00736 add_refinement_element( element_vertices, vertices, 0, 5, 2, 6); 00737 add_refinement_element( element_vertices, vertices, 6, 4, 5, 3); 00738 add_refinement_element( element_vertices, vertices, 6, 5, 2, 3); 00739 } 00740 } 00741 else //split edge 03 again, introduce line 61 00742 { 00743 add_refinement_element( element_vertices, vertices, 4, 1, 5, 6); 00744 add_refinement_element( element_vertices, vertices, 6, 1, 5, 3); 00745 add_refinement_element( element_vertices, vertices, 6, 5, 2, 3); 00746 00747 if (stable_line_is_longer(mesh, vertices, 0, 1, 1, 2)) //split edge 01 again, introduce line 42 00748 { 00749 add_refinement_element( element_vertices, vertices, 0, 4, 2, 6); 00750 add_refinement_element( element_vertices, vertices, 6, 4, 5, 2); 00751 } 00752 else //split edge 12 again, introduce line 50 00753 { 00754 add_refinement_element( element_vertices, vertices, 0, 4, 5, 6); 00755 add_refinement_element( element_vertices, vertices, 0, 5, 2, 6); 00756 } 00757 } 00758 00759 } 00760 00774 template<typename MeshT, typename ElementsVerticesHandleContainerT, typename VertexHandleIteratorType> 00775 static void apply3_4(MeshT const & mesh, 00776 ElementsVerticesHandleContainerT & element_vertices, 00777 VertexHandleIteratorType vertices) 00778 { 00779 00780 // Strategy: The two longest edges of the common vertex are split 'twice', 00781 // i.e. two new edges start from the center of the two longest edges 00782 00783 if (stable_line_is_longer(mesh, vertices, 0, 1, 1, 2)) //split edge 01 again, introduce line 42 00784 { 00785 add_refinement_element( element_vertices, vertices, 0, 4, 2, 6); 00786 add_refinement_element( element_vertices, vertices, 0, 4, 6, 3); 00787 add_refinement_element( element_vertices, vertices, 4, 5, 2, 6); 00788 00789 //if (diag12_len > diag23_len) //split edge 12 again, introduce line 53 00790 if (stable_line_is_longer(mesh, vertices, 1, 2, 2, 3)) //split edge 12 again, introduce line 53 00791 { 00792 add_refinement_element( element_vertices, vertices, 4, 1, 5, 3); 00793 add_refinement_element( element_vertices, vertices, 4, 5, 6, 3); 00794 } 00795 else //split edge 23 again, introduce line 61 00796 { 00797 add_refinement_element( element_vertices, vertices, 4, 1, 6, 3); 00798 add_refinement_element( element_vertices, vertices, 4, 1, 5, 6); 00799 } 00800 } 00801 else //split edge 12, introduce line 50 00802 { 00803 //if (diag12_len > diag23_len) //split edge 12 again, introduce line 53 00804 if (stable_line_is_longer(mesh, vertices, 1, 2, 2, 3)) //split edge 12 again, introduce line 53 00805 { 00806 add_refinement_element( element_vertices, vertices, 0, 4, 5, 3); 00807 add_refinement_element( element_vertices, vertices, 0, 5, 6, 3); 00808 add_refinement_element( element_vertices, vertices, 0, 5, 2, 6); 00809 add_refinement_element( element_vertices, vertices, 4, 1, 5, 3); 00810 } 00811 else //split edge 23 again, introduce line 61 00812 { 00813 add_refinement_element( element_vertices, vertices, 0, 4, 5, 6); 00814 add_refinement_element( element_vertices, vertices, 0, 4, 6, 3); 00815 add_refinement_element( element_vertices, vertices, 0, 5, 2, 6); 00816 add_refinement_element( element_vertices, vertices, 4, 1, 5, 6); 00817 add_refinement_element( element_vertices, vertices, 4, 1, 6, 3); 00818 } 00819 } 00820 00821 } 00822 00823 00824 00826 template<typename ElementType, typename MeshT, 00827 typename ElementsVerticesHandleContainerT, typename VertexCopyMapT, 00828 typename EdgeRefinementFlagAccessor, typename EdgeToVertexHandleAccessor> 00829 static void apply3(ElementType const & element_in, MeshT const & mesh, 00830 ElementsVerticesHandleContainerT & element_vertices, 00831 VertexCopyMapT & vertex_copy_map_, 00832 EdgeRefinementFlagAccessor const & edge_refinement_flag_accessor, EdgeToVertexHandleAccessor const & edge_to_vertex_handle_accessor) 00833 { 00834 typedef typename viennagrid::result_of::const_element_range<ElementType, vertex_tag>::type VertexOnCellRange; 00835 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator; 00836 typedef typename viennagrid::result_of::const_element_range<ElementType, line_tag>::type EdgeOnCellRange; 00837 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator; 00838 00839 typedef typename viennagrid::result_of::handle<ElementType, vertex_tag>::type VertexHandleType; 00840 typedef typename viennagrid::result_of::element<ElementType, line_tag>::type EdgeType; 00841 00842 static_array< VertexHandleType, boundary_elements<tetrahedron_tag, vertex_tag>::num > vertices; 00843 00844 // 00845 // Step 1: Get vertices from input cell 00846 // 00847 VertexOnCellRange vertices_on_cell = viennagrid::elements<vertex_tag>(element_in); 00848 VertexOnCellIterator vocit = vertices_on_cell.begin(); 00849 vertices[0] = vertex_copy_map_(*vocit); ++vocit; 00850 vertices[1] = vertex_copy_map_(*vocit); ++vocit; 00851 vertices[2] = vertex_copy_map_(*vocit); ++vocit; 00852 vertices[3] = vertex_copy_map_(*vocit); 00853 00854 00855 // 00856 // Step 2: Bring vertices in correct order 00857 // 00858 static_array< VertexHandleType, boundary_elements<tetrahedron_tag, vertex_tag>::num + 3 > ordered_vertices; 00859 EdgeOnCellRange edges_on_cell = viennagrid::elements<line_tag>(element_in); 00860 EdgeOnCellIterator eocit = edges_on_cell.begin(); 00861 EdgeType const & e0 = *eocit; ++eocit; 00862 EdgeType const & e1 = *eocit; ++eocit; 00863 EdgeType const & e2 = *eocit; ++eocit; 00864 EdgeType const & e3 = *eocit; ++eocit; 00865 EdgeType const & e4 = *eocit; ++eocit; 00866 EdgeType const & e5 = *eocit; 00867 00868 //with e0 00869 if (edge_refinement_flag_accessor(e0) == true) 00870 { 00871 if (edge_refinement_flag_accessor(e1) == true) 00872 { 00873 ordered_vertices[0] = vertices[2]; 00874 ordered_vertices[1] = vertices[0]; 00875 ordered_vertices[2] = vertices[1]; 00876 ordered_vertices[3] = vertices[3]; 00877 ordered_vertices[4] = edge_to_vertex_handle_accessor(e1); 00878 ordered_vertices[5] = edge_to_vertex_handle_accessor(e0); 00879 00880 if (edge_refinement_flag_accessor(e2) == true) 00881 { 00882 ordered_vertices[6] = edge_to_vertex_handle_accessor(e2); 00883 apply3_1(mesh, element_vertices, ordered_vertices); 00884 } 00885 else if (edge_refinement_flag_accessor(e3) == true) 00886 { 00887 ordered_vertices[6] = edge_to_vertex_handle_accessor(e3); 00888 apply3_2(mesh, element_vertices, ordered_vertices); 00889 } 00890 else if (edge_refinement_flag_accessor(e4) == true) 00891 { 00892 ordered_vertices[6] = edge_to_vertex_handle_accessor(e4); 00893 apply3_4(mesh, element_vertices, ordered_vertices); 00894 } 00895 else if (edge_refinement_flag_accessor(e5) == true) 00896 { 00897 ordered_vertices[6] = edge_to_vertex_handle_accessor(e5); 00898 apply3_3(mesh, element_vertices, ordered_vertices); 00899 } 00900 else 00901 { 00902 assert(false && "Logic error: No edge for refinement found!"); 00903 } 00904 } 00905 else if (edge_refinement_flag_accessor(e2) == true) 00906 { 00907 ordered_vertices[0] = vertices[1]; 00908 ordered_vertices[1] = vertices[0]; 00909 ordered_vertices[2] = vertices[3]; 00910 ordered_vertices[3] = vertices[2]; 00911 ordered_vertices[4] = edge_to_vertex_handle_accessor(e0); 00912 ordered_vertices[5] = edge_to_vertex_handle_accessor(e2); 00913 00914 if (edge_refinement_flag_accessor(e3) == true) 00915 { 00916 ordered_vertices[6] = edge_to_vertex_handle_accessor(e3); 00917 apply3_3(mesh, element_vertices, ordered_vertices); 00918 } 00919 else if (edge_refinement_flag_accessor(e4) == true) 00920 { 00921 ordered_vertices[6] = edge_to_vertex_handle_accessor(e4); 00922 apply3_2(mesh, element_vertices, ordered_vertices); 00923 } 00924 else if (edge_refinement_flag_accessor(e5) == true) 00925 { 00926 ordered_vertices[6] = edge_to_vertex_handle_accessor(e5); 00927 apply3_4(mesh, element_vertices, ordered_vertices); 00928 } 00929 else 00930 { 00931 assert(false && "Logic error: No edge for refinement found!"); 00932 } 00933 } 00934 else if (edge_refinement_flag_accessor(e3) == true) 00935 { 00936 ordered_vertices[0] = vertices[0]; 00937 ordered_vertices[1] = vertices[1]; 00938 ordered_vertices[2] = vertices[2]; 00939 ordered_vertices[3] = vertices[3]; 00940 ordered_vertices[4] = edge_to_vertex_handle_accessor(e0); 00941 ordered_vertices[5] = edge_to_vertex_handle_accessor(e3); 00942 00943 if (edge_refinement_flag_accessor(e4) == true) 00944 { 00945 ordered_vertices[6] = edge_to_vertex_handle_accessor(e4); 00946 apply3_1(mesh, element_vertices, ordered_vertices); 00947 } 00948 else if (edge_refinement_flag_accessor(e5) == true) 00949 { 00950 ordered_vertices[6] = edge_to_vertex_handle_accessor(e5); 00951 apply3_4(mesh, element_vertices, ordered_vertices); 00952 } 00953 else 00954 { 00955 assert(false && "Logic error: No edge for refinement found!"); 00956 } 00957 } 00958 else if (edge_refinement_flag_accessor(e4) == true) 00959 { 00960 ordered_vertices[0] = vertices[3]; 00961 ordered_vertices[1] = vertices[1]; 00962 ordered_vertices[2] = vertices[0]; 00963 ordered_vertices[3] = vertices[2]; 00964 ordered_vertices[4] = edge_to_vertex_handle_accessor(e4); 00965 ordered_vertices[5] = edge_to_vertex_handle_accessor(e0); 00966 00967 if (edge_refinement_flag_accessor(e5) == true) 00968 { 00969 ordered_vertices[6] = edge_to_vertex_handle_accessor(e5); 00970 apply3_3(mesh, element_vertices, ordered_vertices); 00971 } 00972 else 00973 { 00974 assert(false && "Logic error: No edge for refinement found!"); 00975 } 00976 } 00977 else //no second edge 00978 { 00979 assert(false && "Logic error: No edge for refinement found!"); 00980 } 00981 } 00982 else if (edge_refinement_flag_accessor(e1) == true) 00983 { 00984 if (edge_refinement_flag_accessor(e2) == true) 00985 { 00986 ordered_vertices[0] = vertices[3]; 00987 ordered_vertices[1] = vertices[0]; 00988 ordered_vertices[2] = vertices[2]; 00989 ordered_vertices[3] = vertices[1]; 00990 ordered_vertices[4] = edge_to_vertex_handle_accessor(e2); 00991 ordered_vertices[5] = edge_to_vertex_handle_accessor(e1); 00992 00993 if (edge_refinement_flag_accessor(e3) == true) 00994 { 00995 ordered_vertices[6] = edge_to_vertex_handle_accessor(e3); 00996 apply3_4(mesh, element_vertices, ordered_vertices); 00997 } 00998 else if (edge_refinement_flag_accessor(e4) == true) 00999 { 01000 ordered_vertices[6] = edge_to_vertex_handle_accessor(e4); 01001 apply3_3(mesh, element_vertices, ordered_vertices); 01002 } 01003 else if (edge_refinement_flag_accessor(e5) == true) 01004 { 01005 ordered_vertices[6] = edge_to_vertex_handle_accessor(e5); 01006 apply3_2(mesh, element_vertices, ordered_vertices); 01007 } 01008 else 01009 { 01010 assert(false && "Logic error: No edge for refinement found!"); 01011 } 01012 } 01013 else if (edge_refinement_flag_accessor(e3) == true) 01014 { 01015 ordered_vertices[0] = vertices[1]; 01016 ordered_vertices[1] = vertices[2]; 01017 ordered_vertices[2] = vertices[0]; 01018 ordered_vertices[3] = vertices[3]; 01019 ordered_vertices[4] = edge_to_vertex_handle_accessor(e3); 01020 ordered_vertices[5] = edge_to_vertex_handle_accessor(e1); 01021 01022 if (edge_refinement_flag_accessor(e4) == true) 01023 { 01024 ordered_vertices[6] = edge_to_vertex_handle_accessor(e4); 01025 apply3_3(mesh, element_vertices, ordered_vertices); 01026 } 01027 else if (edge_refinement_flag_accessor(e5) == true) 01028 { 01029 ordered_vertices[6] = edge_to_vertex_handle_accessor(e5); 01030 apply3_1(mesh, element_vertices, ordered_vertices); 01031 } 01032 else 01033 { 01034 assert(false && "Logic error: No edge for refinement found!"); 01035 } 01036 } 01037 else if (edge_refinement_flag_accessor(e4) == true) 01038 { 01039 if (edge_refinement_flag_accessor(e5) == true) 01040 { 01041 //make edges 4 and 5 the references 01042 ordered_vertices[0] = vertices[1]; 01043 ordered_vertices[1] = vertices[3]; 01044 ordered_vertices[2] = vertices[2]; 01045 ordered_vertices[3] = vertices[0]; 01046 ordered_vertices[4] = edge_to_vertex_handle_accessor(e4); 01047 ordered_vertices[5] = edge_to_vertex_handle_accessor(e5); 01048 ordered_vertices[6] = edge_to_vertex_handle_accessor(e1); 01049 01050 apply3_4(mesh, element_vertices, ordered_vertices); 01051 } 01052 else 01053 { 01054 assert(false && "Logic error: No edge for refinement found!"); 01055 } 01056 } 01057 else //no second edge 01058 { 01059 assert(false && "Logic error: No edge for refinement found!"); 01060 } 01061 } 01062 else if (edge_refinement_flag_accessor(e2) == true) 01063 { 01064 if (edge_refinement_flag_accessor(e3) == true) 01065 { 01066 //NOTE: edges 2 and 3 don't have a common vertex, therefore 'base facet' is chosen depending on the third edge 01067 01068 if (edge_refinement_flag_accessor(e4) == true) 01069 { 01070 // take edges 2 and 4 as reference 01071 ordered_vertices[0] = vertices[0]; 01072 ordered_vertices[1] = vertices[3]; 01073 ordered_vertices[2] = vertices[1]; 01074 ordered_vertices[3] = vertices[2]; 01075 ordered_vertices[4] = edge_to_vertex_handle_accessor(e2); 01076 ordered_vertices[5] = edge_to_vertex_handle_accessor(e4); 01077 ordered_vertices[6] = edge_to_vertex_handle_accessor(e3); 01078 01079 apply3_4(mesh, element_vertices, ordered_vertices); 01080 } 01081 else if (edge_refinement_flag_accessor(e5) == true) 01082 { 01083 // take edges 5 and 3 as reference 01084 ordered_vertices[0] = vertices[3]; 01085 ordered_vertices[1] = vertices[2]; 01086 ordered_vertices[2] = vertices[1]; 01087 ordered_vertices[3] = vertices[0]; 01088 ordered_vertices[4] = edge_to_vertex_handle_accessor(e5); 01089 ordered_vertices[5] = edge_to_vertex_handle_accessor(e3); 01090 ordered_vertices[6] = edge_to_vertex_handle_accessor(e2); 01091 01092 apply3_3(mesh, element_vertices, ordered_vertices); 01093 } 01094 else 01095 { 01096 assert(false && "Logic error: No edge for refinement found!"); 01097 } 01098 } 01099 else if (edge_refinement_flag_accessor(e4) == true) 01100 { 01101 ordered_vertices[0] = vertices[0]; 01102 ordered_vertices[1] = vertices[3]; 01103 ordered_vertices[2] = vertices[1]; 01104 ordered_vertices[3] = vertices[2]; 01105 ordered_vertices[4] = edge_to_vertex_handle_accessor(e2); 01106 ordered_vertices[5] = edge_to_vertex_handle_accessor(e4); 01107 01108 if (edge_refinement_flag_accessor(e5) == true) 01109 { 01110 ordered_vertices[6] = edge_to_vertex_handle_accessor(e5); 01111 apply3_1(mesh, element_vertices, ordered_vertices); 01112 } 01113 else 01114 { 01115 assert(false && "Logic error: No edge for refinement found!"); 01116 } 01117 } 01118 else //no second edge 01119 { 01120 assert(false && "Logic error: No edge for refinement found!"); 01121 } 01122 } 01123 else if (edge_refinement_flag_accessor(e3) == true) 01124 { 01125 if (edge_refinement_flag_accessor(e4) == true) 01126 { 01127 ordered_vertices[0] = vertices[2]; 01128 ordered_vertices[1] = vertices[1]; 01129 ordered_vertices[2] = vertices[3]; 01130 ordered_vertices[3] = vertices[0]; 01131 ordered_vertices[4] = edge_to_vertex_handle_accessor(e3); 01132 ordered_vertices[5] = edge_to_vertex_handle_accessor(e4); 01133 01134 if (edge_refinement_flag_accessor(e5) == true) 01135 { 01136 ordered_vertices[6] = edge_to_vertex_handle_accessor(e5); 01137 apply3_2(mesh, element_vertices, ordered_vertices); 01138 } 01139 else 01140 { 01141 assert(false && "Logic error: No edge for refinement found!"); 01142 } 01143 } 01144 else //no second edge 01145 { 01146 assert(false && "Logic error: No edge for refinement found!"); 01147 } 01148 } 01149 else //no first edge found 01150 { 01151 assert(false && "Logic error: No edge for refinement found!"); 01152 } 01153 01154 } //apply3() 01155 01156 01157 01158 01160 01174 template<typename MeshT, typename ElementsVerticesHandleContainerT, typename VertexHandleIteratorType> 01175 static void apply4_1(MeshT const & mesh, 01176 ElementsVerticesHandleContainerT & element_vertices, 01177 VertexHandleIteratorType vertices) 01178 { 01179 01180 //if (diag03_len > diag13_len) //split edge 03, introduce line 71 01181 if (stable_line_is_longer(mesh, vertices, 0, 3, 1, 3)) //split edge 03, introduce line 71 01182 { 01183 add_refinement_element( element_vertices, vertices, 0, 1, 4, 7); 01184 add_refinement_element( element_vertices, vertices, 7, 5, 6, 3); 01185 01186 //if (diag13_len > diag23_len) //split edge 13, introduce line 52 01187 if (stable_line_is_longer(mesh, vertices, 1, 3, 2, 3)) //split edge 13, introduce line 52 01188 { 01189 //std::cout << "apply_4_1_1" << std::endl; 01190 add_refinement_element( element_vertices, vertices, 7, 1, 4, 5); 01191 add_refinement_element( element_vertices, vertices, 7, 5, 4, 6); 01192 add_refinement_element( element_vertices, vertices, 1, 2, 4, 5); 01193 add_refinement_element( element_vertices, vertices, 4, 5, 2, 6); 01194 } 01195 else //split edge 23, introduce line 61 01196 { 01197 //std::cout << "apply_4_1_2" << std::endl; 01198 add_refinement_element( element_vertices, vertices, 7, 1, 6, 5); 01199 add_refinement_element( element_vertices, vertices, 7, 1, 4, 6); 01200 add_refinement_element( element_vertices, vertices, 1, 2, 4, 6); 01201 } 01202 } 01203 else //split edge 13, introduce line 50 01204 { 01205 add_refinement_element( element_vertices, vertices, 0, 1, 4, 5); 01206 add_refinement_element( element_vertices, vertices, 0, 5, 4, 7); 01207 add_refinement_element( element_vertices, vertices, 7, 5, 6, 3); 01208 add_refinement_element( element_vertices, vertices, 7, 5, 4, 6); 01209 01210 //if (diag13_len > diag23_len) //split edge 13 again, introduce line 52 01211 if (stable_line_is_longer(mesh, vertices, 1, 3, 2, 3)) //split edge 13 again, introduce line 52 01212 { 01213 //std::cout << "apply_4_1_3" << std::endl; 01214 add_refinement_element( element_vertices, vertices, 1, 2, 4, 5); 01215 add_refinement_element( element_vertices, vertices, 4, 5, 2, 6); 01216 } 01217 else //split edge 23, introduce line 61 01218 { 01219 //std::cout << "apply_4_1_4" << std::endl; 01220 add_refinement_element( element_vertices, vertices, 5, 1, 4, 6); 01221 add_refinement_element( element_vertices, vertices, 4, 1, 2, 6); 01222 } 01223 } 01224 01225 } 01226 01241 template<typename MeshT, typename ElementsVerticesHandleContainerT, typename VertexHandleIteratorType> 01242 static void apply4_2(MeshT const & mesh, 01243 ElementsVerticesHandleContainerT & element_vertices, 01244 VertexHandleIteratorType vertices) 01245 { 01246 //if (diag03_len > diag13_len) //split edge 03, introduce line 61 01247 if (stable_line_is_longer(mesh, vertices, 0, 3, 1, 3)) //split edge 03, introduce line 61 01248 { 01249 //if (diag13_len > diag12_len) //split edge 13, introduce line 72 01250 if (stable_line_is_longer(mesh, vertices, 1, 3, 1, 2)) //split edge 13, introduce line 72 01251 { 01252 //if (diag02_len > diag03_len) //split edge 02, introduce line 53 01253 if (stable_line_is_longer(mesh, vertices, 0, 2, 0, 3)) //split edge 02, introduce line 53 01254 { 01255 //if (diag02_len > diag12_len) //split edge 02, introduce line 51 01256 if (stable_line_is_longer(mesh, vertices, 0, 2, 1, 2)) //split edge 02, introduce line 51 01257 { 01258 //std::cout << "apply_4_2_1_" << std::endl; 01259 add_refinement_element( element_vertices, vertices, 0, 1, 5, 6); 01260 add_refinement_element( element_vertices, vertices, 6, 1, 5, 7); 01261 add_refinement_element( element_vertices, vertices, 1, 4, 5, 7); 01262 add_refinement_element( element_vertices, vertices, 4, 2, 5, 7); 01263 add_refinement_element( element_vertices, vertices, 6, 7, 5, 3); 01264 add_refinement_element( element_vertices, vertices, 7, 2, 5, 3); 01265 } 01266 else //split edge 12, introduce line 40 01267 { 01268 assert( false && "Logic error: diag02 < diag12 not possible here!"); 01269 } 01270 } 01271 else //split edge 03, introduce line 62 01272 { 01273 //std::cout << "split!" << std::endl; 01274 //if (diag02_len > diag12_len) //split edge 02, introduce line 51 01275 if (stable_line_is_longer(mesh, vertices, 0, 2, 1, 2)) //split edge 02, introduce line 51 01276 { 01277 //std::cout << "split!" << std::endl; 01278 //std::cout << "apply_4_2_2" << std::endl; 01279 add_refinement_element( element_vertices, vertices, 0, 1, 5, 6); 01280 add_refinement_element( element_vertices, vertices, 1, 4, 5, 6); 01281 add_refinement_element( element_vertices, vertices, 1, 4, 6, 7); 01282 add_refinement_element( element_vertices, vertices, 7, 4, 6, 2); 01283 add_refinement_element( element_vertices, vertices, 4, 2, 5, 6); 01284 add_refinement_element( element_vertices, vertices, 6, 7, 2, 3); 01285 //std::cout << "done!" << std::endl; 01286 } 01287 else //split edge 12, introduce line 40 01288 { 01289 //std::cout << "apply_4_2_3" << std::endl; 01290 add_refinement_element( element_vertices, vertices, 0, 1, 4, 6); 01291 add_refinement_element( element_vertices, vertices, 0, 4, 5, 6); 01292 add_refinement_element( element_vertices, vertices, 6, 1, 4, 7); 01293 add_refinement_element( element_vertices, vertices, 6, 7, 4, 2); 01294 add_refinement_element( element_vertices, vertices, 6, 4, 5, 2); 01295 add_refinement_element( element_vertices, vertices, 6, 7, 2, 3); 01296 } 01297 } 01298 } 01299 else //split edge 12, introduce line 43 01300 { 01301 //if (diag02_len > diag03_len) //split edge 02, introduce line 53 01302 if (stable_line_is_longer(mesh, vertices, 0, 2, 0, 3)) //split edge 02, introduce line 53 01303 { 01304 //if (diag02_len > diag12_len) //split edge 02, introduce line 51 01305 if (stable_line_is_longer(mesh, vertices, 0, 2, 1, 2)) //split edge 02, introduce line 51 01306 { 01307 //std::cout << "apply_4_2_4" << std::endl; 01308 add_refinement_element( element_vertices, vertices, 0, 1, 5, 6); 01309 add_refinement_element( element_vertices, vertices, 6, 1, 5, 7); 01310 add_refinement_element( element_vertices, vertices, 1, 4, 5, 7); 01311 add_refinement_element( element_vertices, vertices, 5, 4, 2, 3); 01312 add_refinement_element( element_vertices, vertices, 5, 7, 4, 3); 01313 add_refinement_element( element_vertices, vertices, 6, 7, 5, 3); 01314 } 01315 else //split edge 12, introduce line 40 01316 { 01317 //std::cout << "apply_4_2_5" << std::endl; 01318 add_refinement_element( element_vertices, vertices, 0, 1, 4, 6); 01319 add_refinement_element( element_vertices, vertices, 6, 1, 4, 7); 01320 add_refinement_element( element_vertices, vertices, 6, 7, 4, 3); 01321 add_refinement_element( element_vertices, vertices, 6, 4, 5, 3); 01322 add_refinement_element( element_vertices, vertices, 0, 4, 5, 6); 01323 add_refinement_element( element_vertices, vertices, 5, 4, 2, 3); 01324 } 01325 } 01326 else //split edge 03, introduce line 62 01327 { 01328 //if (diag02_len > diag12_len) //split edge 02, introduce line 51 01329 if (stable_line_is_longer(mesh, vertices, 0, 2, 1, 2)) //split edge 02, introduce line 51 01330 { 01331 //std::cout << "apply_4_2_6" << std::endl; 01332 add_refinement_element( element_vertices, vertices, 0, 1, 5, 6); 01333 add_refinement_element( element_vertices, vertices, 1, 4, 5, 6); 01334 add_refinement_element( element_vertices, vertices, 1, 4, 6, 7); 01335 add_refinement_element( element_vertices, vertices, 7, 4, 6, 3); 01336 add_refinement_element( element_vertices, vertices, 4, 2, 5, 6); 01337 add_refinement_element( element_vertices, vertices, 6, 4, 2, 3); 01338 } 01339 else //split edge 12, introduce line 40 01340 { 01341 //std::cout << "apply_4_2_7" << std::endl; 01342 add_refinement_element( element_vertices, vertices, 0, 1, 4, 6); 01343 add_refinement_element( element_vertices, vertices, 6, 1, 4, 7); 01344 add_refinement_element( element_vertices, vertices, 6, 7, 4, 3); 01345 add_refinement_element( element_vertices, vertices, 0, 4, 5, 6); 01346 add_refinement_element( element_vertices, vertices, 5, 4, 2, 6); 01347 add_refinement_element( element_vertices, vertices, 6, 4, 2, 3); 01348 } 01349 } 01350 } 01351 } 01352 else //split edge 13, introduce line 70 01353 { 01354 //if (diag13_len > diag12_len) //split edge 13, introduce line 72 01355 if (stable_line_is_longer(mesh, vertices, 1, 3, 1, 2)) //split edge 13, introduce line 72 01356 { 01357 //if (diag02_len > diag03_len) //split edge 02, introduce line 53 01358 if (stable_line_is_longer(mesh, vertices, 0, 2, 0, 3)) //split edge 02, introduce line 53 01359 { 01360 //if (diag02_len > diag12_len) //split edge 02, introduce line 51 01361 if (stable_line_is_longer(mesh, vertices, 0, 2, 1, 2)) //split edge 02, introduce line 51 01362 { 01363 //std::cout << "apply_4_2_8" << std::endl; 01364 add_refinement_element( element_vertices, vertices, 0, 1, 5, 7); 01365 add_refinement_element( element_vertices, vertices, 0, 7, 5, 6); 01366 add_refinement_element( element_vertices, vertices, 1, 4, 5, 7); 01367 add_refinement_element( element_vertices, vertices, 5, 4, 2, 7); 01368 add_refinement_element( element_vertices, vertices, 5, 7, 2, 3); 01369 add_refinement_element( element_vertices, vertices, 6, 7, 5, 3); 01370 } 01371 else //split edge 12, introduce line 40 01372 { 01373 //std::cout << "apply_4_2_9" << std::endl; 01374 add_refinement_element( element_vertices, vertices, 0, 1, 4, 7); 01375 add_refinement_element( element_vertices, vertices, 0, 4, 5, 7); 01376 add_refinement_element( element_vertices, vertices, 0, 7, 5, 6); 01377 add_refinement_element( element_vertices, vertices, 6, 7, 5, 3); 01378 add_refinement_element( element_vertices, vertices, 5, 4, 2, 7); 01379 add_refinement_element( element_vertices, vertices, 5, 7, 2, 3); 01380 } 01381 } 01382 else //split edge 03, introduce line 62 01383 { 01384 //if (diag02_len > diag12_len) //split edge 02, introduce line 51 01385 if (stable_line_is_longer(mesh, vertices, 0, 2, 1, 2)) //split edge 02, introduce line 51 01386 { 01387 //std::cout << "apply_4_2_10" << std::endl; 01388 add_refinement_element( element_vertices, vertices, 0, 1, 5, 7); 01389 add_refinement_element( element_vertices, vertices, 1, 4, 5, 7); 01390 add_refinement_element( element_vertices, vertices, 0, 7, 5, 6); 01391 add_refinement_element( element_vertices, vertices, 4, 2, 5, 7); 01392 add_refinement_element( element_vertices, vertices, 7, 2, 5, 6); 01393 add_refinement_element( element_vertices, vertices, 7, 2, 6, 3); 01394 } 01395 else //split edge 12, introduce line 40 01396 { 01397 //std::cout << "apply_4_2_11" << std::endl; 01398 add_refinement_element( element_vertices, vertices, 0, 1, 4, 7); 01399 add_refinement_element( element_vertices, vertices, 0, 7, 4, 6); 01400 add_refinement_element( element_vertices, vertices, 0, 4, 5, 6); 01401 add_refinement_element( element_vertices, vertices, 6, 7, 4, 2); 01402 add_refinement_element( element_vertices, vertices, 6, 4, 5, 2); 01403 add_refinement_element( element_vertices, vertices, 6, 7, 2, 3); 01404 } 01405 } 01406 } 01407 else //split edge 12, introduce line 43 01408 { 01409 //if (diag02_len > diag03_len) //split edge 02, introduce line 53 01410 if (stable_line_is_longer(mesh, vertices, 0, 2, 0, 3)) //split edge 02, introduce line 53 01411 { 01412 //if (diag02_len > diag12_len) //split edge 02, introduce line 51 01413 if (stable_line_is_longer(mesh, vertices, 0, 2, 1, 2)) //split edge 02, introduce line 51 01414 { 01415 //std::cout << "apply_4_2_12" << std::endl; 01416 add_refinement_element( element_vertices, vertices, 0, 1, 5, 7); 01417 add_refinement_element( element_vertices, vertices, 0, 7, 5, 6); 01418 add_refinement_element( element_vertices, vertices, 6, 7, 5, 3); 01419 add_refinement_element( element_vertices, vertices, 1, 4, 5, 7); 01420 add_refinement_element( element_vertices, vertices, 7, 4, 5, 3); 01421 add_refinement_element( element_vertices, vertices, 4, 2, 5, 3); 01422 } 01423 else //split edge 12, introduce line 40 01424 { 01425 //std::cout << "apply_4_2_13" << std::endl; 01426 add_refinement_element( element_vertices, vertices, 0, 1, 4, 7); 01427 add_refinement_element( element_vertices, vertices, 0, 4, 5, 7); 01428 add_refinement_element( element_vertices, vertices, 0, 7, 5, 6); 01429 add_refinement_element( element_vertices, vertices, 6, 7, 5, 3); 01430 add_refinement_element( element_vertices, vertices, 7, 4, 5, 3); 01431 add_refinement_element( element_vertices, vertices, 4, 2, 5, 3); 01432 } 01433 } 01434 else //split edge 03, introduce line 62 01435 { 01436 //if (diag02_len > diag12_len) //split edge 02, introduce line 51 01437 if (stable_line_is_longer(mesh, vertices, 0, 2, 1, 2)) //split edge 02, introduce line 51 01438 { 01439 //we have diag12_len > diag13_len > diag03_len > diag02_len alreday, hence this case is bogus! 01440 assert( false && "Logic error: diag02 > diag12 not possible here!"); 01441 } 01442 else //split edge 12, introduce line 40 01443 { 01444 //std::cout << "apply_4_2_14" << std::endl; 01445 add_refinement_element( element_vertices, vertices, 0, 1, 4, 7); 01446 add_refinement_element( element_vertices, vertices, 0, 7, 4, 6); 01447 add_refinement_element( element_vertices, vertices, 0, 4, 5, 6); 01448 add_refinement_element( element_vertices, vertices, 6, 7, 4, 3); 01449 add_refinement_element( element_vertices, vertices, 6, 4, 2, 3); 01450 add_refinement_element( element_vertices, vertices, 5, 4, 2, 6); 01451 } 01452 } 01453 } 01454 } 01455 } 01456 01457 01459 template<typename ElementType, typename MeshT, 01460 typename ElementsVerticesHandleContainerT, typename VertexCopyMapT, 01461 typename EdgeRefinementFlagAccessor, typename EdgeToVertexHandleAccessor> 01462 static void apply4(ElementType const & element_in, MeshT const & mesh, 01463 ElementsVerticesHandleContainerT & element_vertices, 01464 VertexCopyMapT & vertex_copy_map_, 01465 EdgeRefinementFlagAccessor const & edge_refinement_flag_accessor, EdgeToVertexHandleAccessor const & edge_to_vertex_handle_accessor) 01466 { 01467 typedef typename viennagrid::result_of::const_element_range<ElementType, vertex_tag>::type VertexOnCellRange; 01468 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator; 01469 typedef typename viennagrid::result_of::const_element_range<ElementType, line_tag>::type EdgeOnCellRange; 01470 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator; 01471 01472 typedef typename viennagrid::result_of::handle<ElementType, vertex_tag>::type VertexHandleType; 01473 typedef typename viennagrid::result_of::element<ElementType, line_tag>::type EdgeType; 01474 01475 static_array< VertexHandleType, boundary_elements<tetrahedron_tag, vertex_tag>::num > vertices; 01476 01477 // 01478 // Step 1: Get vertices from input cell 01479 // 01480 VertexOnCellRange vertices_on_cell = viennagrid::elements<vertex_tag>(element_in); 01481 VertexOnCellIterator vocit = vertices_on_cell.begin(); 01482 vertices[0] = vertex_copy_map_(*vocit); ++vocit; 01483 vertices[1] = vertex_copy_map_(*vocit); ++vocit; 01484 vertices[2] = vertex_copy_map_(*vocit); ++vocit; 01485 vertices[3] = vertex_copy_map_(*vocit); 01486 01487 01488 // 01489 // Step 2: Bring vertices in correct order, such that refined edge is on {0,1}-edge 01490 // 01491 static_array< VertexHandleType, boundary_elements<tetrahedron_tag, vertex_tag>::num + 4 > ordered_vertices; 01492 EdgeOnCellRange edges_on_cell = viennagrid::elements<line_tag>(element_in); 01493 EdgeOnCellIterator eocit = edges_on_cell.begin(); 01494 EdgeType const & e0 = *eocit; ++eocit; 01495 EdgeType const & e1 = *eocit; ++eocit; 01496 EdgeType const & e2 = *eocit; ++eocit; 01497 EdgeType const & e3 = *eocit; ++eocit; 01498 EdgeType const & e4 = *eocit; ++eocit; 01499 EdgeType const & e5 = *eocit; 01500 01501 //with e0 01502 if (edge_refinement_flag_accessor(e0) == false) 01503 { 01504 if (edge_refinement_flag_accessor(e1) == false) 01505 { 01506 ordered_vertices[0] = vertices[2]; 01507 ordered_vertices[1] = vertices[0]; 01508 ordered_vertices[2] = vertices[1]; 01509 ordered_vertices[3] = vertices[3]; 01510 ordered_vertices[4] = edge_to_vertex_handle_accessor(e3); 01511 ordered_vertices[5] = edge_to_vertex_handle_accessor(e2); 01512 ordered_vertices[6] = edge_to_vertex_handle_accessor(e4); 01513 ordered_vertices[7] = edge_to_vertex_handle_accessor(e5); 01514 01515 apply4_1(mesh, element_vertices, ordered_vertices); 01516 } 01517 else if (edge_refinement_flag_accessor(e2) == false) 01518 { 01519 ordered_vertices[0] = vertices[1]; 01520 ordered_vertices[1] = vertices[0]; 01521 ordered_vertices[2] = vertices[3]; 01522 ordered_vertices[3] = vertices[2]; 01523 ordered_vertices[4] = edge_to_vertex_handle_accessor(e4); 01524 ordered_vertices[5] = edge_to_vertex_handle_accessor(e1); 01525 ordered_vertices[6] = edge_to_vertex_handle_accessor(e5); 01526 ordered_vertices[7] = edge_to_vertex_handle_accessor(e3); 01527 01528 apply4_1(mesh, element_vertices, ordered_vertices); 01529 } 01530 else if (edge_refinement_flag_accessor(e3) == false) 01531 { 01532 ordered_vertices[0] = vertices[0]; 01533 ordered_vertices[1] = vertices[1]; 01534 ordered_vertices[2] = vertices[2]; 01535 ordered_vertices[3] = vertices[3]; 01536 ordered_vertices[4] = edge_to_vertex_handle_accessor(e1); 01537 ordered_vertices[5] = edge_to_vertex_handle_accessor(e4); 01538 ordered_vertices[6] = edge_to_vertex_handle_accessor(e5); 01539 ordered_vertices[7] = edge_to_vertex_handle_accessor(e2); 01540 01541 apply4_1(mesh, element_vertices, ordered_vertices); 01542 } 01543 else if (edge_refinement_flag_accessor(e4) == false) 01544 { 01545 ordered_vertices[0] = vertices[3]; 01546 ordered_vertices[1] = vertices[1]; 01547 ordered_vertices[2] = vertices[0]; 01548 ordered_vertices[3] = vertices[2]; 01549 ordered_vertices[4] = edge_to_vertex_handle_accessor(e2); 01550 ordered_vertices[5] = edge_to_vertex_handle_accessor(e3); 01551 ordered_vertices[6] = edge_to_vertex_handle_accessor(e1); 01552 ordered_vertices[7] = edge_to_vertex_handle_accessor(e5); 01553 01554 apply4_1(mesh, element_vertices, ordered_vertices); 01555 } 01556 else if (edge_refinement_flag_accessor(e5) == false) 01557 { 01558 ordered_vertices[0] = vertices[0]; 01559 ordered_vertices[1] = vertices[1]; 01560 ordered_vertices[2] = vertices[2]; 01561 ordered_vertices[3] = vertices[3]; 01562 ordered_vertices[4] = edge_to_vertex_handle_accessor(e3); 01563 ordered_vertices[5] = edge_to_vertex_handle_accessor(e1); 01564 ordered_vertices[6] = edge_to_vertex_handle_accessor(e2); 01565 ordered_vertices[7] = edge_to_vertex_handle_accessor(e4); 01566 01567 apply4_2(mesh, element_vertices, ordered_vertices); 01568 } 01569 else 01570 { 01571 assert(false && "Logic error: No edge for refinement found!"); 01572 } 01573 } 01574 else if (edge_refinement_flag_accessor(e1) == false) 01575 { 01576 if (edge_refinement_flag_accessor(e2) == false) 01577 { 01578 ordered_vertices[0] = vertices[3]; 01579 ordered_vertices[1] = vertices[0]; 01580 ordered_vertices[2] = vertices[2]; 01581 ordered_vertices[3] = vertices[1]; 01582 ordered_vertices[4] = edge_to_vertex_handle_accessor(e5); 01583 ordered_vertices[5] = edge_to_vertex_handle_accessor(e0); 01584 ordered_vertices[6] = edge_to_vertex_handle_accessor(e3); 01585 ordered_vertices[7] = edge_to_vertex_handle_accessor(e4); 01586 01587 apply4_1(mesh, element_vertices, ordered_vertices); 01588 } 01589 else if (edge_refinement_flag_accessor(e3) == false) 01590 { 01591 ordered_vertices[0] = vertices[1]; 01592 ordered_vertices[1] = vertices[2]; 01593 ordered_vertices[2] = vertices[0]; 01594 ordered_vertices[3] = vertices[3]; 01595 ordered_vertices[4] = edge_to_vertex_handle_accessor(e0); 01596 ordered_vertices[5] = edge_to_vertex_handle_accessor(e5); 01597 ordered_vertices[6] = edge_to_vertex_handle_accessor(e2); 01598 ordered_vertices[7] = edge_to_vertex_handle_accessor(e4); 01599 01600 apply4_1(mesh, element_vertices, ordered_vertices); 01601 } 01602 else if (edge_refinement_flag_accessor(e4) == false) 01603 { 01604 ordered_vertices[0] = vertices[2]; 01605 ordered_vertices[1] = vertices[0]; 01606 ordered_vertices[2] = vertices[1]; 01607 ordered_vertices[3] = vertices[3]; 01608 ordered_vertices[4] = edge_to_vertex_handle_accessor(e0); 01609 ordered_vertices[5] = edge_to_vertex_handle_accessor(e3); 01610 ordered_vertices[6] = edge_to_vertex_handle_accessor(e5); 01611 ordered_vertices[7] = edge_to_vertex_handle_accessor(e2); 01612 01613 apply4_2(mesh, element_vertices, ordered_vertices); 01614 } 01615 else if (edge_refinement_flag_accessor(e5) == false) 01616 { 01617 ordered_vertices[0] = vertices[0]; 01618 ordered_vertices[1] = vertices[2]; 01619 ordered_vertices[2] = vertices[3]; 01620 ordered_vertices[3] = vertices[1]; 01621 ordered_vertices[4] = edge_to_vertex_handle_accessor(e2); 01622 ordered_vertices[5] = edge_to_vertex_handle_accessor(e3); 01623 ordered_vertices[6] = edge_to_vertex_handle_accessor(e4); 01624 ordered_vertices[7] = edge_to_vertex_handle_accessor(e0); 01625 01626 apply4_1(mesh, element_vertices, ordered_vertices); 01627 } 01628 else 01629 { 01630 assert(false && "Logic error: No edge for refinement found!"); 01631 } 01632 } 01633 else if (edge_refinement_flag_accessor(e2) == false) 01634 { 01635 if (edge_refinement_flag_accessor(e3) == false) 01636 { 01637 ordered_vertices[0] = vertices[3]; 01638 ordered_vertices[1] = vertices[0]; 01639 ordered_vertices[2] = vertices[2]; 01640 ordered_vertices[3] = vertices[1]; 01641 ordered_vertices[4] = edge_to_vertex_handle_accessor(e1); 01642 ordered_vertices[5] = edge_to_vertex_handle_accessor(e5); 01643 ordered_vertices[6] = edge_to_vertex_handle_accessor(e4); 01644 ordered_vertices[7] = edge_to_vertex_handle_accessor(e0); 01645 01646 apply4_2(mesh, element_vertices, ordered_vertices); 01647 } 01648 else if (edge_refinement_flag_accessor(e4) == false) 01649 { 01650 ordered_vertices[0] = vertices[0]; 01651 ordered_vertices[1] = vertices[3]; 01652 ordered_vertices[2] = vertices[1]; 01653 ordered_vertices[3] = vertices[2]; 01654 ordered_vertices[4] = edge_to_vertex_handle_accessor(e0); 01655 ordered_vertices[5] = edge_to_vertex_handle_accessor(e5); 01656 ordered_vertices[6] = edge_to_vertex_handle_accessor(e3); 01657 ordered_vertices[7] = edge_to_vertex_handle_accessor(e1); 01658 01659 apply4_1(mesh, element_vertices, ordered_vertices); 01660 } 01661 else if (edge_refinement_flag_accessor(e5) == false) 01662 { 01663 ordered_vertices[0] = vertices[2]; 01664 ordered_vertices[1] = vertices[3]; 01665 ordered_vertices[2] = vertices[0]; 01666 ordered_vertices[3] = vertices[1]; 01667 ordered_vertices[4] = edge_to_vertex_handle_accessor(e1); 01668 ordered_vertices[5] = edge_to_vertex_handle_accessor(e4); 01669 ordered_vertices[6] = edge_to_vertex_handle_accessor(e0); 01670 ordered_vertices[7] = edge_to_vertex_handle_accessor(e3); 01671 01672 apply4_1(mesh, element_vertices, ordered_vertices); 01673 } 01674 else 01675 { 01676 assert(false && "Logic error: No edge for refinement found!"); 01677 } 01678 } 01679 else if (edge_refinement_flag_accessor(e3) == false) 01680 { 01681 if (edge_refinement_flag_accessor(e4) == false) 01682 { 01683 ordered_vertices[0] = vertices[2]; 01684 ordered_vertices[1] = vertices[1]; 01685 ordered_vertices[2] = vertices[3]; 01686 ordered_vertices[3] = vertices[0]; 01687 ordered_vertices[4] = edge_to_vertex_handle_accessor(e5); 01688 ordered_vertices[5] = edge_to_vertex_handle_accessor(e0); 01689 ordered_vertices[6] = edge_to_vertex_handle_accessor(e2); 01690 ordered_vertices[7] = edge_to_vertex_handle_accessor(e1); 01691 01692 apply4_1(mesh, element_vertices, ordered_vertices); 01693 } 01694 else if (edge_refinement_flag_accessor(e5) == false) 01695 { 01696 ordered_vertices[0] = vertices[3]; 01697 ordered_vertices[1] = vertices[2]; 01698 ordered_vertices[2] = vertices[1]; 01699 ordered_vertices[3] = vertices[0]; 01700 ordered_vertices[4] = edge_to_vertex_handle_accessor(e4); 01701 ordered_vertices[5] = edge_to_vertex_handle_accessor(e1); 01702 ordered_vertices[6] = edge_to_vertex_handle_accessor(e0); 01703 ordered_vertices[7] = edge_to_vertex_handle_accessor(e2); 01704 01705 apply4_1(mesh, element_vertices, ordered_vertices); 01706 } 01707 else 01708 { 01709 assert(false && "Logic error: No edge for refinement found!"); 01710 } 01711 } 01712 else if (edge_refinement_flag_accessor(e4) == false) 01713 { 01714 if (edge_refinement_flag_accessor(e5) == false) 01715 { 01716 ordered_vertices[0] = vertices[1]; 01717 ordered_vertices[1] = vertices[3]; 01718 ordered_vertices[2] = vertices[2]; 01719 ordered_vertices[3] = vertices[0]; 01720 ordered_vertices[4] = edge_to_vertex_handle_accessor(e3); 01721 ordered_vertices[5] = edge_to_vertex_handle_accessor(e2); 01722 ordered_vertices[6] = edge_to_vertex_handle_accessor(e1); 01723 ordered_vertices[7] = edge_to_vertex_handle_accessor(e0); 01724 01725 apply4_1(mesh, element_vertices, ordered_vertices); 01726 } 01727 else 01728 { 01729 assert(false && "Logic error: No edge for refinement found!"); 01730 } 01731 } 01732 else 01733 { 01734 assert(false && "Logic error: No edge for refinement found!"); 01735 } 01736 } 01737 01738 01739 01741 01742 01756 template<typename MeshT, typename ElementsVerticesHandleContainerT, typename VertexHandleIteratorType> 01757 static void apply5_1(MeshT const & mesh, 01758 ElementsVerticesHandleContainerT & element_vertices, 01759 VertexHandleIteratorType vertices) 01760 { 01761 add_refinement_element( element_vertices, vertices, 6, 7, 8, 3); 01762 add_refinement_element( element_vertices, vertices, 5, 4, 2, 8); 01763 01764 01765 //if (diag03_len > diag13_len) //split edge 03, introduce line 61 01766 if (stable_line_is_longer(mesh, vertices, 0, 3, 1, 3)) //split edge 03, introduce line 61 01767 { 01768 add_refinement_element( element_vertices, vertices, 6, 4, 5, 8); 01769 add_refinement_element( element_vertices, vertices, 6, 7, 4, 8); 01770 add_refinement_element( element_vertices, vertices, 1, 4, 6, 7); 01771 01772 //if (diag02_len > diag12_len) //split edge 02, introduce line 51 01773 if (stable_line_is_longer(mesh, vertices, 0, 2, 1, 2)) //split edge 02, introduce line 51 01774 { 01775 add_refinement_element( element_vertices, vertices, 0, 1, 5, 6); 01776 add_refinement_element( element_vertices, vertices, 1, 4, 5, 6); 01777 } 01778 else //split edge 12, introduce line 40 01779 { 01780 add_refinement_element( element_vertices, vertices, 0, 1, 4, 6); 01781 add_refinement_element( element_vertices, vertices, 0, 4, 5, 6); 01782 } 01783 } 01784 else //split edge 13, introduce line 70 01785 { 01786 //if (diag02_len > diag12_len) //split edge 02, introduce line 51 01787 if (stable_line_is_longer(mesh, vertices, 0, 2, 1, 2)) //split edge 02, introduce line 51 01788 { 01789 add_refinement_element( element_vertices, vertices, 0, 1, 5, 7); 01790 add_refinement_element( element_vertices, vertices, 0, 7, 5, 6); 01791 add_refinement_element( element_vertices, vertices, 1, 4, 5, 7); 01792 add_refinement_element( element_vertices, vertices, 7, 4, 5, 8); 01793 add_refinement_element( element_vertices, vertices, 6, 7, 5, 8); 01794 } 01795 else //split edge 12, introduce line 40 01796 { 01797 add_refinement_element( element_vertices, vertices, 0, 1, 4, 7); 01798 add_refinement_element( element_vertices, vertices, 0, 7, 4, 6); 01799 add_refinement_element( element_vertices, vertices, 0, 4, 5, 6); 01800 add_refinement_element( element_vertices, vertices, 6, 4, 5, 8); 01801 add_refinement_element( element_vertices, vertices, 6, 7, 4, 8); 01802 } 01803 } 01804 } 01805 01807 template<typename ElementType, typename MeshT, 01808 typename ElementsVerticesHandleContainerT, typename VertexCopyMapT, 01809 typename EdgeRefinementFlagAccessor, typename EdgeToVertexHandleAccessor> 01810 static void apply5(ElementType const & element_in, MeshT const & mesh, 01811 ElementsVerticesHandleContainerT & element_vertices, 01812 VertexCopyMapT & vertex_copy_map_, 01813 EdgeRefinementFlagAccessor const & edge_refinement_flag_accessor, EdgeToVertexHandleAccessor const & edge_to_vertex_handle_accessor) 01814 { 01815 typedef typename viennagrid::result_of::const_element_range<ElementType, vertex_tag>::type VertexOnCellRange; 01816 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator; 01817 typedef typename viennagrid::result_of::const_element_range<ElementType, line_tag>::type EdgeOnCellRange; 01818 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator; 01819 01820 typedef typename viennagrid::result_of::handle<ElementType, vertex_tag>::type VertexHandleType; 01821 typedef typename viennagrid::result_of::element<ElementType, line_tag>::type EdgeType; 01822 01823 static_array< VertexHandleType, boundary_elements<tetrahedron_tag, vertex_tag>::num > vertices; 01824 01825 // 01826 // Step 1: Get vertices from input cell 01827 // 01828 VertexOnCellRange vertices_on_cell = viennagrid::elements<vertex_tag>(element_in); 01829 VertexOnCellIterator vocit = vertices_on_cell.begin(); 01830 vertices[0] = vertex_copy_map_(*vocit); ++vocit; 01831 vertices[1] = vertex_copy_map_(*vocit); ++vocit; 01832 vertices[2] = vertex_copy_map_(*vocit); ++vocit; 01833 vertices[3] = vertex_copy_map_(*vocit); 01834 01835 01836 // 01837 // Step 2: Bring vertices in correct order, such that refined edge is on {0,1}-edge 01838 // 01839 static_array< VertexHandleType, boundary_elements<tetrahedron_tag, vertex_tag>::num + 5 > ordered_vertices; 01840 EdgeOnCellRange edges_on_cell = viennagrid::elements<line_tag>(element_in); 01841 EdgeOnCellIterator eocit = edges_on_cell.begin(); 01842 EdgeType const & e0 = *eocit; ++eocit; 01843 EdgeType const & e1 = *eocit; ++eocit; 01844 EdgeType const & e2 = *eocit; ++eocit; 01845 EdgeType const & e3 = *eocit; ++eocit; 01846 EdgeType const & e4 = *eocit; ++eocit; 01847 EdgeType const & e5 = *eocit; 01848 01849 if (edge_refinement_flag_accessor(e0) == false) 01850 { 01851 ordered_vertices[0] = vertices[0]; 01852 ordered_vertices[1] = vertices[1]; 01853 ordered_vertices[2] = vertices[2]; 01854 ordered_vertices[3] = vertices[3]; 01855 ordered_vertices[4] = edge_to_vertex_handle_accessor(e3); 01856 ordered_vertices[5] = edge_to_vertex_handle_accessor(e1); 01857 ordered_vertices[6] = edge_to_vertex_handle_accessor(e2); 01858 ordered_vertices[7] = edge_to_vertex_handle_accessor(e4); 01859 ordered_vertices[8] = edge_to_vertex_handle_accessor(e5); 01860 01861 apply5_1(mesh, element_vertices, ordered_vertices); 01862 } 01863 else if (edge_refinement_flag_accessor(e1) == false) 01864 { 01865 ordered_vertices[0] = vertices[2]; 01866 ordered_vertices[1] = vertices[0]; 01867 ordered_vertices[2] = vertices[1]; 01868 ordered_vertices[3] = vertices[3]; 01869 ordered_vertices[4] = edge_to_vertex_handle_accessor(e0); 01870 ordered_vertices[5] = edge_to_vertex_handle_accessor(e3); 01871 ordered_vertices[6] = edge_to_vertex_handle_accessor(e5); 01872 ordered_vertices[7] = edge_to_vertex_handle_accessor(e2); 01873 ordered_vertices[8] = edge_to_vertex_handle_accessor(e4); 01874 01875 apply5_1(mesh, element_vertices, ordered_vertices); 01876 } 01877 else if (edge_refinement_flag_accessor(e2) == false) 01878 { 01879 ordered_vertices[0] = vertices[0]; 01880 ordered_vertices[1] = vertices[3]; 01881 ordered_vertices[2] = vertices[1]; 01882 ordered_vertices[3] = vertices[2]; 01883 ordered_vertices[4] = edge_to_vertex_handle_accessor(e4); 01884 ordered_vertices[5] = edge_to_vertex_handle_accessor(e0); 01885 ordered_vertices[6] = edge_to_vertex_handle_accessor(e1); 01886 ordered_vertices[7] = edge_to_vertex_handle_accessor(e5); 01887 ordered_vertices[8] = edge_to_vertex_handle_accessor(e3); 01888 01889 apply5_1(mesh, element_vertices, ordered_vertices); 01890 } 01891 else if (edge_refinement_flag_accessor(e3) == false) 01892 { 01893 ordered_vertices[0] = vertices[1]; 01894 ordered_vertices[1] = vertices[2]; 01895 ordered_vertices[2] = vertices[0]; 01896 ordered_vertices[3] = vertices[3]; 01897 ordered_vertices[4] = edge_to_vertex_handle_accessor(e1); 01898 ordered_vertices[5] = edge_to_vertex_handle_accessor(e0); 01899 ordered_vertices[6] = edge_to_vertex_handle_accessor(e4); 01900 ordered_vertices[7] = edge_to_vertex_handle_accessor(e5); 01901 ordered_vertices[8] = edge_to_vertex_handle_accessor(e2); 01902 01903 apply5_1(mesh, element_vertices, ordered_vertices); 01904 } 01905 else if (edge_refinement_flag_accessor(e4) == false) 01906 { 01907 ordered_vertices[0] = vertices[1]; 01908 ordered_vertices[1] = vertices[3]; 01909 ordered_vertices[2] = vertices[2]; 01910 ordered_vertices[3] = vertices[0]; 01911 ordered_vertices[4] = edge_to_vertex_handle_accessor(e5); 01912 ordered_vertices[5] = edge_to_vertex_handle_accessor(e3); 01913 ordered_vertices[6] = edge_to_vertex_handle_accessor(e0); 01914 ordered_vertices[7] = edge_to_vertex_handle_accessor(e2); 01915 ordered_vertices[8] = edge_to_vertex_handle_accessor(e1); 01916 01917 apply5_1(mesh, element_vertices, ordered_vertices); 01918 } 01919 else if (edge_refinement_flag_accessor(e5) == false) 01920 { 01921 ordered_vertices[0] = vertices[3]; 01922 ordered_vertices[1] = vertices[2]; 01923 ordered_vertices[2] = vertices[1]; 01924 ordered_vertices[3] = vertices[0]; 01925 ordered_vertices[4] = edge_to_vertex_handle_accessor(e3); 01926 ordered_vertices[5] = edge_to_vertex_handle_accessor(e4); 01927 ordered_vertices[6] = edge_to_vertex_handle_accessor(e2); 01928 ordered_vertices[7] = edge_to_vertex_handle_accessor(e1); 01929 ordered_vertices[8] = edge_to_vertex_handle_accessor(e0); 01930 01931 apply5_1(mesh, element_vertices, ordered_vertices); 01932 } 01933 else 01934 { 01935 assert(false && "Logic error: No edge for refinement found!"); 01936 } 01937 } 01938 01939 01940 01941 01943 01944 01945 01947 template<typename ElementType, typename MeshT, 01948 typename ElementsVerticesHandleContainerT, typename VertexCopyMapT, 01949 typename EdgeRefinementFlagAccessor, typename EdgeToVertexHandleAccessor> 01950 static void apply6(ElementType const & element_in, MeshT const & mesh, 01951 ElementsVerticesHandleContainerT & element_vertices, 01952 VertexCopyMapT & vertex_copy_map_, 01953 EdgeRefinementFlagAccessor const &, EdgeToVertexHandleAccessor const & edge_to_vertex_handle_accessor) 01954 { 01955 typedef typename viennagrid::result_of::const_element_range<ElementType, vertex_tag>::type VertexOnCellRange; 01956 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator; 01957 typedef typename viennagrid::result_of::const_element_range<ElementType, line_tag>::type EdgeOnCellRange; 01958 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator; 01959 01960 typedef typename viennagrid::result_of::handle<ElementType, vertex_tag>::type VertexHandleType; 01961 01962 static_array< VertexHandleType, boundary_elements<tetrahedron_tag, vertex_tag>::num + 01963 boundary_elements<tetrahedron_tag, line_tag>::num> vertices; 01964 01965 // 01966 // Step 1: Get vertices on the new mesh 01967 // 01968 01969 //grab existing vertices: 01970 VertexOnCellRange vertices_on_cell = viennagrid::elements<vertex_tag>(element_in); 01971 VertexOnCellIterator vocit = vertices_on_cell.begin(); 01972 vertices[0] = vertex_copy_map_(*vocit); ++vocit; 01973 vertices[1] = vertex_copy_map_(*vocit); ++vocit; 01974 vertices[2] = vertex_copy_map_(*vocit); ++vocit; 01975 vertices[3] = vertex_copy_map_(*vocit); 01976 01977 //add vertices from edge 01978 EdgeOnCellRange edges_on_cell = viennagrid::elements<line_tag>(element_in); 01979 EdgeOnCellIterator eocit = edges_on_cell.begin(); 01980 01981 01982 01983 vertices[4] = edge_to_vertex_handle_accessor(*eocit); ++eocit; 01984 vertices[5] = edge_to_vertex_handle_accessor(*eocit); ++eocit; 01985 vertices[6] = edge_to_vertex_handle_accessor(*eocit); ++eocit; 01986 vertices[7] = edge_to_vertex_handle_accessor(*eocit); ++eocit; 01987 vertices[8] = edge_to_vertex_handle_accessor(*eocit); ++eocit; 01988 vertices[9] = edge_to_vertex_handle_accessor(*eocit); 01989 01990 // 01991 // Step 2: Add new cells to new mesh: 01992 // 01993 01994 //0-4-5-6: 01995 add_refinement_element( element_vertices, vertices, 0, 4, 5, 6); 01996 01997 //1-7-4-8: 01998 add_refinement_element( element_vertices, vertices, 1, 7, 4, 8); 01999 02000 //2-5-7-9: 02001 add_refinement_element( element_vertices, vertices, 2, 5, 7, 9); 02002 02003 //3-8-6-9: 02004 add_refinement_element( element_vertices, vertices, 3, 8, 6, 9); 02005 02006 double diag58 = viennagrid::norm( viennagrid::point(mesh, vertices[5]) - viennagrid::point(mesh, vertices[8]) ); 02007 double diag67 = viennagrid::norm( viennagrid::point(mesh, vertices[6]) - viennagrid::point(mesh, vertices[7]) ); 02008 double diag49 = viennagrid::norm( viennagrid::point(mesh, vertices[4]) - viennagrid::point(mesh, vertices[9]) ); 02009 02010 if ( (diag58 <= diag67) && (diag58 <= diag49) ) //diag58 is shortest: keep it, split others 02011 { 02012 //4-8-5-6: 02013 add_refinement_element( element_vertices, vertices, 4, 8, 5, 6); 02014 02015 //4-8-7-5: 02016 add_refinement_element( element_vertices, vertices, 4, 8, 7, 5); 02017 02018 //7-8-9-5: 02019 add_refinement_element( element_vertices, vertices, 7, 8, 9, 5); 02020 02021 //8-6-9-5: 02022 add_refinement_element( element_vertices, vertices, 8, 6, 9, 5); 02023 } 02024 else if ( (diag67 <= diag58) && (diag67 <= diag49) ) //diag67 is shortest: keep it, split others 02025 { 02026 //4-7-6-8: 02027 add_refinement_element( element_vertices, vertices, 4, 7, 6, 8); 02028 02029 //4-7-5-6: 02030 add_refinement_element( element_vertices, vertices, 4, 7, 5, 6); 02031 02032 //7-9-6-8: 02033 add_refinement_element( element_vertices, vertices, 7, 9, 6, 8); 02034 02035 //7-9-5-6: 02036 add_refinement_element( element_vertices, vertices, 7, 9, 5, 6); 02037 } 02038 else //keep shortest diagonal diag49 02039 { 02040 //4-9-6-8: 02041 add_refinement_element( element_vertices, vertices, 4, 9, 6, 8); 02042 02043 //4-9-5-6: 02044 add_refinement_element( element_vertices, vertices, 4, 9, 5, 6); 02045 02046 //4-7-9-8: 02047 add_refinement_element( element_vertices, vertices, 4, 7, 9, 8); 02048 02049 //4-7-5-9: 02050 add_refinement_element( element_vertices, vertices, 4, 7, 5, 9); 02051 } 02052 02053 } //apply6() 02054 02055 02056 // 02066 template<typename ElementType, 02067 typename MeshT, typename ElementsVerticesHandleContainerT, 02068 typename VertexCopyMapT, 02069 typename EdgeRefinementFlagAccessor, typename EdgeToVertexHandleAccessor> 02070 static void apply(ElementType const & element_in, MeshT const & mesh, 02071 ElementsVerticesHandleContainerT & element_vertices, 02072 VertexCopyMapT & vertex_copy_map_, 02073 EdgeRefinementFlagAccessor const & edge_refinement_flag_accessor, 02074 EdgeToVertexHandleAccessor const & edge_to_vertex_handle_accessor) 02075 { 02076 typedef typename viennagrid::result_of::const_element_range<ElementType, line_tag>::type EdgeOnCellRange; 02077 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator; 02078 02079 std::size_t edges_to_refine = 0; 02080 EdgeOnCellRange edges_on_cell(element_in); 02081 for (EdgeOnCellIterator eocit = edges_on_cell.begin(); 02082 eocit != edges_on_cell.end(); 02083 ++eocit) 02084 { 02085 if (edge_refinement_flag_accessor(*eocit) == true) 02086 ++edges_to_refine; 02087 } 02088 02089 switch (edges_to_refine) 02090 { 02091 case 0: apply0(element_in, mesh, element_vertices, vertex_copy_map_, edge_refinement_flag_accessor, edge_to_vertex_handle_accessor); break; 02092 case 1: apply1(element_in, mesh, element_vertices, vertex_copy_map_, edge_refinement_flag_accessor, edge_to_vertex_handle_accessor); break; 02093 case 2: apply2(element_in, mesh, element_vertices, vertex_copy_map_, edge_refinement_flag_accessor, edge_to_vertex_handle_accessor); break; 02094 case 3: apply3(element_in, mesh, element_vertices, vertex_copy_map_, edge_refinement_flag_accessor, edge_to_vertex_handle_accessor); break; 02095 case 4: apply4(element_in, mesh, element_vertices, vertex_copy_map_, edge_refinement_flag_accessor, edge_to_vertex_handle_accessor); break; 02096 case 5: apply5(element_in, mesh, element_vertices, vertex_copy_map_, edge_refinement_flag_accessor, edge_to_vertex_handle_accessor); break; 02097 case 6: apply6(element_in, mesh, element_vertices, vertex_copy_map_, edge_refinement_flag_accessor, edge_to_vertex_handle_accessor); break; 02098 default: //nothing to do... 02099 break; 02100 } 02101 } //apply() 02102 02103 }; 02104 02105 } // namespace detail 02106 02107 } 02108 02109 #endif 02110
1.7.6.1