• Main Page
  • Namespaces
  • Data Structures
  • Files
  • File List
  • Globals

/export/development/ViennaGrid/release/ViennaGrid-1.0.0/viennagrid/topology/tetrahedron.hpp

Go to the documentation of this file.
00001 #ifndef VIENNAGRID_TOPOLOGY_TETRAHEDRON_HPP
00002 #define VIENNAGRID_TOPOLOGY_TETRAHEDRON_HPP
00003 
00004 /* =======================================================================
00005    Copyright (c) 2011, Institute for Microelectronics,
00006                        Institute for Analysis and Scientific Computing,
00007                        TU Wien.
00008 
00009                             -----------------
00010                      ViennaGrid - The Vienna Grid Library
00011                             -----------------
00012 
00013    Authors:      Karl Rupp                           rupp@iue.tuwien.ac.at
00014                  Josef Weinbub                    weinbub@iue.tuwien.ac.at
00015                
00016    (A list of additional contributors can be found in the PDF manual)
00017 
00018    License:      MIT (X11), see file LICENSE in the base directory
00019 ======================================================================= */
00020 
00021 #include "viennagrid/forwards.h"
00022 #include "viennagrid/topology/point.hpp"
00023 #include "viennagrid/topology/line.hpp"
00024 #include "viennagrid/topology/triangle.hpp"
00025 #include "viennadata/api.hpp"
00026 #include "viennagrid/algorithm/norm.hpp"
00027 
00032 namespace viennagrid
00033 {
00034 
00036   template <>
00037   struct simplex_tag<3>
00038   {
00039     enum{ dim = 3 };
00040     static std::string name() { return "Tetrahedron"; }
00041   };
00042 
00043 
00044   namespace topology
00045   {
00046 
00047     //Tetrahedron:
00049     template <>
00050     struct bndcells<tetrahedron_tag, 0>
00051     {
00052       typedef point_tag             tag;
00053 
00054       enum{ num = 4 };     //4 vertices
00055     };
00056 
00058     template <>
00059     struct bndcells<tetrahedron_tag, 1>
00060     {
00061       typedef simplex_tag<1>              tag;
00062 
00063       enum{ num = 6 };     //6 edges
00064     };
00065 
00067     template <>
00068     struct bndcells<tetrahedron_tag, 2>
00069     {
00070       typedef triangle_tag          tag;
00071 
00072       enum{ num = 4 };     //4 facets
00073     };
00074     
00076 
00077     template <>
00078     struct bndcell_filler<tetrahedron_tag, 1>
00079     {
00080       //fill edges:
00081       template <typename ElementType, typename Vertices, typename Orientations, typename Segment>
00082       static void fill(ElementType ** elements, Vertices ** vertices, Orientations * orientations, Segment & seg)
00083       {
00084         Vertices * edgevertices[2];
00085         ElementType edge;
00086 
00087         //fill edges according to orientation and ordering induced by k-tuple-metafunction:
00088         edgevertices[0] = vertices[0];
00089         edgevertices[1] = vertices[1];
00090         edge.vertices(edgevertices);
00091         elements[0] = seg.push_back(edge, (orientations == NULL) ? NULL : orientations);
00092 
00093         edgevertices[0] = vertices[0];
00094         edgevertices[1] = vertices[2];
00095         edge.vertices(edgevertices);
00096         elements[1] = seg.push_back(edge, (orientations == NULL) ? NULL : orientations + 1);
00097 
00098         edgevertices[0] = vertices[0];
00099         edgevertices[1] = vertices[3];
00100         edge.vertices(edgevertices);
00101         elements[2] = seg.push_back(edge, (orientations == NULL) ? NULL : orientations + 2);
00102 
00103         edgevertices[0] = vertices[1];
00104         edgevertices[1] = vertices[2];
00105         edge.vertices(edgevertices);
00106         elements[3] = seg.push_back(edge, (orientations == NULL) ? NULL : orientations + 3);
00107 
00108         edgevertices[0] = vertices[1];
00109         edgevertices[1] = vertices[3];
00110         edge.vertices(edgevertices);
00111         elements[4] = seg.push_back(edge, (orientations == NULL) ? NULL : orientations + 4);
00112 
00113         edgevertices[0] = vertices[2];
00114         edgevertices[1] = vertices[3];
00115         edge.vertices(edgevertices);
00116         elements[5] = seg.push_back(edge, (orientations == NULL) ? NULL : orientations + 5);
00117       }
00118     };
00119     
00121     template <>
00122     struct bndcell_filler<tetrahedron_tag, 2>
00123     {
00124       //fill facets:
00125       template <typename ElementType, typename Vertices, typename Orientations, typename Segment>
00126       static void fill(ElementType ** elements, Vertices ** vertices, Orientations * orientations, Segment & seg)
00127       {
00128         Vertices * facetvertices[3];
00129         ElementType facet;
00130 
00131         //fill edges according to orientation and ordering induced by k-tuple-metafunction:
00132 
00133         facetvertices[0] = vertices[0];
00134         facetvertices[1] = vertices[1];
00135         facetvertices[2] = vertices[2];
00136         facet.vertices(facetvertices);
00137         elements[0] = seg.push_back(facet, orientations );
00138         //this new facet must be initialized too:
00139         elements[0]->fill(seg);
00140 
00141         facetvertices[0] = vertices[0];
00142         facetvertices[1] = vertices[1];
00143         facetvertices[2] = vertices[3];
00144         facet.vertices(facetvertices);
00145         elements[1] = seg.push_back(facet, (orientations == NULL) ? NULL : orientations + 1 );
00146         elements[1]->fill(seg);
00147 
00148         facetvertices[0] = vertices[0];
00149         facetvertices[1] = vertices[2];
00150         facetvertices[2] = vertices[3];
00151         facet.vertices(facetvertices);
00152         elements[2] = seg.push_back(facet, (orientations == NULL) ? NULL : orientations + 2 );
00153         elements[2]->fill(seg);
00154 
00155         facetvertices[0] = vertices[1];
00156         facetvertices[1] = vertices[2];
00157         facetvertices[2] = vertices[3];
00158         facet.vertices(facetvertices);
00159         elements[3] = seg.push_back(facet, (orientations == NULL) ? NULL : orientations + 3 );
00160         elements[3]->fill(seg);
00161 
00162       }
00163     };
00164 
00165   } //topology
00166   
00167   
00168   
00174   template <typename VertexType>
00175   bool stable_line_is_longer(VertexType const * v1_1, VertexType const * v1_2,  
00176                              VertexType const * v2_1, VertexType const * v2_2)
00177   {
00178     typedef typename VertexType::config_type::numeric_type       ScalarType;
00179     
00180     VertexType const * v1_1_ptr = (v1_1->id() < v1_2->id()) ? v1_1 : v1_2; //v1_1 carries smaller ID
00181     VertexType const * v1_2_ptr = (v1_1->id() < v1_2->id()) ? v1_2 : v1_1; //v1_2 carries larger ID
00182     
00183     VertexType const * v2_1_ptr = (v2_1->id() < v2_2->id()) ? v2_1 : v2_2; //v2_1 carries smaller ID
00184     VertexType const * v2_2_ptr = (v2_1->id() < v2_2->id()) ? v2_2 : v2_1; //v2_2 carries larger ID
00185     
00186     ScalarType line1 = viennagrid::norm(v1_1->point() - v1_2->point());
00187     ScalarType line2 = viennagrid::norm(v2_1->point() - v2_2->point());
00188     
00189     if (line1 > line2)
00190     {
00191       return true;
00192     }
00193     else if (line1 == line2)
00194     {
00195       //compare IDs:
00196       if (v1_1_ptr->id() > v2_1_ptr->id())
00197       {
00198         return true;
00199       }
00200       else if (v1_1_ptr->id() == v2_1_ptr->id())
00201       {
00202         if (v1_2_ptr->id() > v2_2_ptr->id())
00203           return true;
00204         else if (v1_2_ptr->id() == v2_2_ptr->id())
00205           return false; //identical lines are compared!
00206       }
00207     }
00208 
00209     return false;
00210   }
00211   
00212   
00213   
00219   template <>
00220   struct element_refinement<tetrahedron_tag>
00221   {
00222 
00224     template <typename CellType, typename DomainTypeOut>
00225     static void apply0(CellType const & cell_in, DomainTypeOut & segment_out)
00226     {
00227       //std::cout << "tetrahedron::apply0()" << std::endl;
00228       typedef typename CellType::config_type        ConfigTypeIn;
00229       typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type            VertexOnCellRange;
00230       typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type         VertexOnCellIterator;            
00231       typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type            EdgeOnCellRange;
00232       typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type           EdgeOnCellIterator;            
00233       
00234       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type             VertexType;
00235 
00236       VertexType * vertices[topology::bndcells<tetrahedron_tag, 0>::num];
00237       
00238       // Step 1: grab existing vertices:
00239       VertexOnCellRange vertices_on_cell = viennagrid::ncells<0>(cell_in);
00240       VertexOnCellIterator vocit = vertices_on_cell.begin();
00241       vertices[0] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00242       vertices[1] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00243       vertices[2] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00244       vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]);
00245 
00246       // Step 2: Add new cells to new domain:
00247       CellType new_cell;
00248       new_cell.vertices(vertices);
00249       segment_out.push_back(new_cell);
00250 
00251     }    
00252 
00254     template <typename CellType, typename DomainTypeOut>
00255     static void apply1(CellType const & cell_in, DomainTypeOut & segment_out)
00256     {
00257       typedef typename CellType::config_type        ConfigTypeIn;
00258       typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type            VertexOnCellRange;
00259       typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type         VertexOnCellIterator;            
00260       typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type            EdgeOnCellRange;
00261       typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type           EdgeOnCellIterator;            
00262       
00263       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type             VertexType;
00264       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 1>::type             EdgeType;
00265 
00266       VertexType * vertices[topology::bndcells<tetrahedron_tag, 0>::num];
00267       //std::cout << "apply1()" << std::endl;
00268       
00269       //
00270       // Step 1: Get vertices from input cell
00271       //
00272       VertexOnCellRange vertices_on_cell = viennagrid::ncells<0>(cell_in);
00273       VertexOnCellIterator vocit = vertices_on_cell.begin();
00274       vertices[0] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00275       vertices[1] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00276       vertices[2] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00277       vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]);
00278       
00279 
00280       //
00281       // Step 2: Bring vertices in correct order, such that refined edge is on {0,1}-edge
00282       //
00283       VertexType * ordered_vertices[topology::bndcells<tetrahedron_tag, 0>::num + 1];
00284       EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(cell_in);
00285       EdgeOnCellIterator eocit = edges_on_cell.begin();
00286       EdgeType const & e0 = *eocit; ++eocit;
00287       EdgeType const & e1 = *eocit; ++eocit;
00288       EdgeType const & e2 = *eocit; ++eocit;
00289       EdgeType const & e3 = *eocit; ++eocit;
00290       EdgeType const & e4 = *eocit; ++eocit;
00291       EdgeType const & e5 = *eocit;
00292       
00293       if (viennadata::access<refinement_key, bool>(refinement_key())(e0) == true)
00294       {
00295         ordered_vertices[0] = vertices[0];
00296         ordered_vertices[1] = vertices[1];
00297         ordered_vertices[2] = vertices[2];
00298         ordered_vertices[3] = vertices[3];
00299         ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
00300       }
00301       else if (viennadata::access<refinement_key, bool>(refinement_key())(e1) == true)
00302       {
00303         ordered_vertices[0] = vertices[2];
00304         ordered_vertices[1] = vertices[0];
00305         ordered_vertices[2] = vertices[1];
00306         ordered_vertices[3] = vertices[3];
00307         ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
00308       }
00309       else if (viennadata::access<refinement_key, bool>(refinement_key())(e2) == true)
00310       {
00311         ordered_vertices[0] = vertices[0];
00312         ordered_vertices[1] = vertices[3];
00313         ordered_vertices[2] = vertices[1];
00314         ordered_vertices[3] = vertices[2];
00315         ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
00316       }
00317       else if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == true)
00318       {
00319         ordered_vertices[0] = vertices[1];
00320         ordered_vertices[1] = vertices[2];
00321         ordered_vertices[2] = vertices[0];
00322         ordered_vertices[3] = vertices[3];
00323         ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
00324       }
00325       else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true)
00326       {
00327         ordered_vertices[0] = vertices[3];
00328         ordered_vertices[1] = vertices[1];
00329         ordered_vertices[2] = vertices[0];
00330         ordered_vertices[3] = vertices[2];
00331         ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
00332       }
00333       else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true)
00334       {
00335         ordered_vertices[0] = vertices[3];
00336         ordered_vertices[1] = vertices[2];
00337         ordered_vertices[2] = vertices[1];
00338         ordered_vertices[3] = vertices[0];
00339         ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
00340       }
00341       else
00342       {
00343         assert(false && "Logic error: No edge for refinement found!"); 
00344       }
00345       
00346       
00347       //
00348       // Step 3: Write new cells to domain_out
00349       //
00350       CellType new_cell;
00351       VertexType * cellvertices[topology::bndcells<tetrahedron_tag, 0>::num];
00352       
00353       //cell containing vertex 0:
00354       cellvertices[0] = ordered_vertices[0];
00355       cellvertices[1] = ordered_vertices[4];
00356       cellvertices[2] = ordered_vertices[2];
00357       cellvertices[3] = ordered_vertices[3];
00358       new_cell.vertices(cellvertices);
00359       segment_out.push_back(new_cell);
00360 
00361       //cell without vertex 0:
00362       cellvertices[0] = ordered_vertices[4];
00363       cellvertices[1] = ordered_vertices[1];
00364       cellvertices[2] = ordered_vertices[2];
00365       cellvertices[3] = ordered_vertices[3];
00366       new_cell.vertices(cellvertices);
00367       segment_out.push_back(new_cell);
00368       
00369     }    
00370 
00371 
00373 
00374 
00375     //
00376     // Refinement of a tetrahedron, bisecting two edges
00377     //
00378     
00392     template <typename CellType, typename DomainTypeOut, typename VertexType>
00393     static void apply2_1(CellType const & cell_in, 
00394                          DomainTypeOut & segment_out,
00395                          VertexType ** vertices
00396                         )
00397     {
00398       CellType new_cell;
00399       VertexType * cellvertices[topology::bndcells<tetrahedron_tag, 0>::num];
00400 
00401       cellvertices[0] = vertices[4];
00402       cellvertices[1] = vertices[1];
00403       cellvertices[2] = vertices[5];
00404       cellvertices[3] = vertices[3];
00405       new_cell.vertices(cellvertices);
00406       segment_out.push_back(new_cell);
00407       
00408       //double diag01_len = viennagrid::norm(vertices[0]->point() - vertices[1]->point());
00409       //double diag12_len = viennagrid::norm(vertices[2]->point() - vertices[1]->point());
00410 
00411       //if (diag01_len > diag12_len) //split edge 01, introduce line 42
00412       if (stable_line_is_longer(vertices[0], vertices[1],
00413                                 vertices[2], vertices[1])) //split edge 01, introduce line 42
00414       {
00415         /*std::cout << "Norm " << vertices[0]->id() << vertices[1]->id() << ": " 
00416                              << viennagrid::norm(vertices[0]->point() - vertices[1]->point()) << std::endl;
00417         std::cout << "Norm " << vertices[2]->id() << vertices[1]->id() << ": " 
00418                              << viennagrid::norm(vertices[2]->point() - vertices[1]->point()) << std::endl;
00419         std::cout << "Splitting " << vertices[0]->id() << vertices[1]->id() << std::endl;*/
00420         cellvertices[0] = vertices[0];
00421         cellvertices[1] = vertices[4];
00422         cellvertices[2] = vertices[2];
00423         cellvertices[3] = vertices[3];
00424         new_cell.vertices(cellvertices);
00425         segment_out.push_back(new_cell);
00426         
00427         cellvertices[0] = vertices[4];
00428         cellvertices[1] = vertices[5];
00429         cellvertices[2] = vertices[2];
00430         cellvertices[3] = vertices[3];
00431         new_cell.vertices(cellvertices);
00432         segment_out.push_back(new_cell);
00433       }
00434       else //split edge 12, introduce line 05
00435       {
00436         /*std::cout << "Norm " << vertices[0]->id() << vertices[1]->id() << ": " 
00437                              << viennagrid::norm(vertices[0]->point() - vertices[1]->point()) << std::endl;
00438         std::cout << "Norm " << vertices[2]->id() << vertices[1]->id() << ": " 
00439                              << viennagrid::norm(vertices[2]->point() - vertices[1]->point()) << std::endl;
00440         std::cout << "Splitting " << vertices[2]->id() << vertices[1]->id() << std::endl;*/
00441         cellvertices[0] = vertices[0];
00442         cellvertices[1] = vertices[4];
00443         cellvertices[2] = vertices[5];
00444         cellvertices[3] = vertices[3];
00445         new_cell.vertices(cellvertices);
00446         segment_out.push_back(new_cell);
00447         
00448         cellvertices[0] = vertices[0];
00449         cellvertices[1] = vertices[5];
00450         cellvertices[2] = vertices[2];
00451         cellvertices[3] = vertices[3];
00452         new_cell.vertices(cellvertices);
00453         segment_out.push_back(new_cell);
00454       }
00455       
00456     }
00457 
00471     template <typename CellType, typename DomainTypeOut, typename VertexType>
00472     static void apply2_2(CellType const & cell_in, 
00473                          DomainTypeOut & segment_out,
00474                          VertexType ** vertices
00475                         )
00476     {
00477       CellType new_cell;
00478       VertexType * cellvertices[topology::bndcells<tetrahedron_tag, 0>::num];
00479 
00480       cellvertices[0] = vertices[4];
00481       cellvertices[1] = vertices[1];
00482       cellvertices[2] = vertices[2];
00483       cellvertices[3] = vertices[5];
00484       new_cell.vertices(cellvertices);
00485       segment_out.push_back(new_cell);
00486       
00487       cellvertices[0] = vertices[4];
00488       cellvertices[1] = vertices[1];
00489       cellvertices[2] = vertices[5];
00490       cellvertices[3] = vertices[3];
00491       new_cell.vertices(cellvertices);
00492       segment_out.push_back(new_cell);
00493       
00494       cellvertices[0] = vertices[0];
00495       cellvertices[1] = vertices[4];
00496       cellvertices[2] = vertices[2];
00497       cellvertices[3] = vertices[5];
00498       new_cell.vertices(cellvertices);
00499       segment_out.push_back(new_cell);
00500       
00501       cellvertices[0] = vertices[0];
00502       cellvertices[1] = vertices[4];
00503       cellvertices[2] = vertices[5];
00504       cellvertices[3] = vertices[3];
00505       new_cell.vertices(cellvertices);
00506       segment_out.push_back(new_cell);
00507     }
00508 
00512     template <typename CellType, typename DomainTypeOut>
00513     static void apply2(CellType const & cell_in, DomainTypeOut & segment_out)
00514     {
00515       typedef typename CellType::config_type        ConfigTypeIn;
00516       typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type            VertexOnCellRange;
00517       typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type         VertexOnCellIterator;            
00518       typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type            EdgeOnCellRange;
00519       typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type           EdgeOnCellIterator;            
00520       
00521       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type             VertexType;
00522       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 1>::type             EdgeType;
00523 
00524       VertexType * vertices[topology::bndcells<tetrahedron_tag, 0>::num];
00525       //std::cout << "apply2()" << std::endl;
00526       
00527       //
00528       // Step 1: Get vertices from input cell
00529       //
00530       VertexOnCellRange vertices_on_cell = viennagrid::ncells<0>(cell_in);
00531       VertexOnCellIterator vocit = vertices_on_cell.begin();
00532       vertices[0] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00533       vertices[1] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00534       vertices[2] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00535       vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]);
00536       
00537 
00538       //
00539       // Step 2: Bring vertices in correct order, such that refined edge is on {0,1}-edge
00540       //
00541       VertexType * ordered_vertices[topology::bndcells<tetrahedron_tag, 0>::num + 2];
00542       EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(cell_in);
00543       EdgeOnCellIterator eocit = edges_on_cell.begin();
00544       EdgeType const & e0 = *eocit; ++eocit;
00545       EdgeType const & e1 = *eocit; ++eocit;
00546       EdgeType const & e2 = *eocit; ++eocit;
00547       EdgeType const & e3 = *eocit; ++eocit;
00548       EdgeType const & e4 = *eocit; ++eocit;
00549       EdgeType const & e5 = *eocit;
00550       
00551       //with e0
00552       if (viennadata::access<refinement_key, bool>(refinement_key())(e0) == true)
00553       {
00554         if (viennadata::access<refinement_key, bool>(refinement_key())(e1) == true)
00555         {
00556           ordered_vertices[0] = vertices[2];
00557           ordered_vertices[1] = vertices[0];
00558           ordered_vertices[2] = vertices[1];
00559           ordered_vertices[3] = vertices[3];
00560           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
00561           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
00562           
00563           apply2_1(cell_in, segment_out, ordered_vertices);
00564         }
00565         else if (viennadata::access<refinement_key, bool>(refinement_key())(e2) == true)
00566         {
00567           ordered_vertices[0] = vertices[1];
00568           ordered_vertices[1] = vertices[0];
00569           ordered_vertices[2] = vertices[3];
00570           ordered_vertices[3] = vertices[2];
00571           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
00572           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
00573           
00574           apply2_1(cell_in, segment_out, ordered_vertices);
00575         }
00576         else if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == true) 
00577         {
00578           ordered_vertices[0] = vertices[0];
00579           ordered_vertices[1] = vertices[1];
00580           ordered_vertices[2] = vertices[2];
00581           ordered_vertices[3] = vertices[3];
00582           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
00583           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
00584           
00585           apply2_1(cell_in, segment_out, ordered_vertices);
00586         }
00587         else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true) 
00588         {
00589           ordered_vertices[0] = vertices[3];
00590           ordered_vertices[1] = vertices[1];
00591           ordered_vertices[2] = vertices[0];
00592           ordered_vertices[3] = vertices[2];
00593           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
00594           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
00595           
00596           apply2_1(cell_in, segment_out, ordered_vertices);
00597         }
00598         else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true) 
00599         {
00600           ordered_vertices[0] = vertices[0];
00601           ordered_vertices[1] = vertices[1];
00602           ordered_vertices[2] = vertices[2];
00603           ordered_vertices[3] = vertices[3];
00604           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
00605           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
00606           
00607           apply2_2(cell_in, segment_out, ordered_vertices);
00608         }
00609         else
00610         {
00611           assert(false && "Logic error: No edge for refinement found!"); 
00612         }
00613       }
00614       else if (viennadata::access<refinement_key, bool>(refinement_key())(e1) == true)
00615       {
00616         if (viennadata::access<refinement_key, bool>(refinement_key())(e2) == true)
00617         {
00618           ordered_vertices[0] = vertices[3];
00619           ordered_vertices[1] = vertices[0];
00620           ordered_vertices[2] = vertices[2];
00621           ordered_vertices[3] = vertices[1];
00622           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
00623           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
00624           
00625           apply2_1(cell_in, segment_out, ordered_vertices);
00626         }
00627         else if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == true) 
00628         {
00629           ordered_vertices[0] = vertices[1];
00630           ordered_vertices[1] = vertices[2];
00631           ordered_vertices[2] = vertices[0];
00632           ordered_vertices[3] = vertices[3];
00633           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
00634           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
00635           
00636           apply2_1(cell_in, segment_out, ordered_vertices);
00637         }
00638         else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true) 
00639         {
00640           ordered_vertices[0] = vertices[2];
00641           ordered_vertices[1] = vertices[0];
00642           ordered_vertices[2] = vertices[1];
00643           ordered_vertices[3] = vertices[3];
00644           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
00645           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
00646           
00647           apply2_2(cell_in, segment_out, ordered_vertices);
00648         }
00649         else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true) 
00650         {
00651           ordered_vertices[0] = vertices[0];
00652           ordered_vertices[1] = vertices[2];
00653           ordered_vertices[2] = vertices[3];
00654           ordered_vertices[3] = vertices[1];
00655           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
00656           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
00657           
00658           apply2_1(cell_in, segment_out, ordered_vertices);
00659         }
00660         else
00661         {
00662           assert(false && "Logic error: No edge for refinement found!"); 
00663         }
00664       }
00665       else if (viennadata::access<refinement_key, bool>(refinement_key())(e2) == true)
00666       {
00667         if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == true) 
00668         {
00669           ordered_vertices[0] = vertices[3];
00670           ordered_vertices[1] = vertices[0];
00671           ordered_vertices[2] = vertices[2];
00672           ordered_vertices[3] = vertices[1];
00673           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
00674           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
00675           
00676           apply2_2(cell_in, segment_out, ordered_vertices);
00677         }
00678         else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true) 
00679         {
00680           ordered_vertices[0] = vertices[0];
00681           ordered_vertices[1] = vertices[3];
00682           ordered_vertices[2] = vertices[1];
00683           ordered_vertices[3] = vertices[2];
00684           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
00685           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
00686           
00687           apply2_1(cell_in, segment_out, ordered_vertices);
00688         }
00689         else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true) 
00690         {
00691           ordered_vertices[0] = vertices[2];
00692           ordered_vertices[1] = vertices[3];
00693           ordered_vertices[2] = vertices[0];
00694           ordered_vertices[3] = vertices[1];
00695           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
00696           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
00697           
00698           apply2_1(cell_in, segment_out, ordered_vertices);
00699         }
00700         else
00701         {
00702           assert(false && "Logic error: No edge for refinement found!"); 
00703         }
00704       }
00705       else if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == true)
00706       {
00707         if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true) 
00708         {
00709           ordered_vertices[0] = vertices[2];
00710           ordered_vertices[1] = vertices[1];
00711           ordered_vertices[2] = vertices[3];
00712           ordered_vertices[3] = vertices[0];
00713           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
00714           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
00715           
00716           apply2_1(cell_in, segment_out, ordered_vertices);
00717         }
00718         else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true) 
00719         {
00720           ordered_vertices[0] = vertices[3];
00721           ordered_vertices[1] = vertices[2];
00722           ordered_vertices[2] = vertices[1];
00723           ordered_vertices[3] = vertices[0];
00724           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
00725           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
00726           
00727           apply2_1(cell_in, segment_out, ordered_vertices);
00728         }
00729         else
00730         {
00731           assert(false && "Logic error: No edge for refinement found!"); 
00732         }
00733       }
00734       else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true)
00735       {
00736         if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true) 
00737         {
00738           ordered_vertices[0] = vertices[1];
00739           ordered_vertices[1] = vertices[3];
00740           ordered_vertices[2] = vertices[2];
00741           ordered_vertices[3] = vertices[0];
00742           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
00743           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
00744           
00745           apply2_1(cell_in, segment_out, ordered_vertices);
00746         }
00747         else
00748         {
00749           assert(false && "Logic error: No edge for refinement found!"); 
00750         }
00751       }
00752       else
00753       {
00754         assert(false && "Logic error: No edge for refinement found!"); 
00755       }
00756     }    
00757 
00758 
00759 
00761 
00762 
00776     template <typename CellType, typename DomainTypeOut, typename VertexType>
00777     static void apply3_1(CellType const & cell_in, 
00778                          DomainTypeOut & segment_out,
00779                          VertexType ** vertices
00780                         )
00781     {
00782       CellType new_cell;
00783       VertexType * cellvertices[topology::bndcells<tetrahedron_tag, 0>::num];
00784 
00785       cellvertices[0] = vertices[4];
00786       cellvertices[1] = vertices[1];
00787       cellvertices[2] = vertices[5];
00788       cellvertices[3] = vertices[6];
00789       new_cell.vertices(cellvertices);
00790       segment_out.push_back(new_cell);
00791       
00792       //double diag01_len = viennagrid::norm(vertices[0]->point() - vertices[1]->point());
00793       //double diag12_len = viennagrid::norm(vertices[2]->point() - vertices[1]->point());
00794       //double diag13_len = viennagrid::norm(vertices[3]->point() - vertices[1]->point());
00795 
00796       // Strategy: The two longest edges of the common vertex are split 'twice',
00797       //           i.e. two new edges start from the center of the two longest edges
00798       //if (diag01_len > diag12_len) //split edge 01 again, introduce line 42
00799       if (stable_line_is_longer(vertices[0], vertices[1],
00800                                 vertices[1], vertices[2])) //split edge 01 again, introduce line 42
00801       {
00802         /*std::cout << "Norm " << vertices[0]->id() << vertices[1]->id() << ": " 
00803                              << viennagrid::norm(vertices[0]->point() - vertices[1]->point()) << std::endl;
00804         std::cout << "Norm " << vertices[2]->id() << vertices[1]->id() << ": " 
00805                              << viennagrid::norm(vertices[2]->point() - vertices[1]->point()) << std::endl;
00806         std::cout << "Splitting " << vertices[0]->id() << vertices[1]->id() << std::endl;*/
00807         //if (diag13_len > diag12_len) //split edge 13 again, introduce line 62
00808         if (stable_line_is_longer(vertices[1], vertices[3],
00809                                   vertices[1], vertices[2])) //split edge 13 again, introduce line 62
00810         {
00811           /*std::cout << "Norm " << vertices[1]->id() << vertices[3]->id() << ": " 
00812                                 << viennagrid::norm(vertices[1]->point() - vertices[3]->point()) << std::endl;
00813           std::cout << "Norm " << vertices[2]->id() << vertices[1]->id() << ": " 
00814                                 << viennagrid::norm(vertices[2]->point() - vertices[1]->point()) << std::endl;
00815           std::cout << "Splitting " << vertices[1]->id() << vertices[3]->id() << std::endl; */
00816           
00817           if (stable_line_is_longer(vertices[1], vertices[3],
00818                                     vertices[0], vertices[1])) //split edge 13 again, introduce line 60
00819           {
00820             cellvertices[0] = vertices[0];
00821             cellvertices[1] = vertices[6];
00822             cellvertices[2] = vertices[2];
00823             cellvertices[3] = vertices[3];
00824             new_cell.vertices(cellvertices);
00825             segment_out.push_back(new_cell);
00826             
00827             cellvertices[0] = vertices[0];
00828             cellvertices[1] = vertices[4];
00829             cellvertices[2] = vertices[2];
00830             cellvertices[3] = vertices[6];
00831             new_cell.vertices(cellvertices);
00832             segment_out.push_back(new_cell);
00833             
00834             cellvertices[0] = vertices[4];
00835             cellvertices[1] = vertices[5];
00836             cellvertices[2] = vertices[2];
00837             cellvertices[3] = vertices[6];
00838             new_cell.vertices(cellvertices);
00839             segment_out.push_back(new_cell);
00840           }
00841           else  //split edge 01 again, introduce line 43
00842           {
00843             //std::cout << "apply_3_1_1" << std::endl;
00844             cellvertices[0] = vertices[0];
00845             cellvertices[1] = vertices[4];
00846             cellvertices[2] = vertices[2];
00847             cellvertices[3] = vertices[3];
00848             new_cell.vertices(cellvertices);
00849             segment_out.push_back(new_cell);
00850             
00851             cellvertices[0] = vertices[3];
00852             cellvertices[1] = vertices[4];
00853             cellvertices[2] = vertices[2];
00854             cellvertices[3] = vertices[6];
00855             new_cell.vertices(cellvertices);
00856             segment_out.push_back(new_cell);
00857             
00858             cellvertices[0] = vertices[4];
00859             cellvertices[1] = vertices[5];
00860             cellvertices[2] = vertices[2];
00861             cellvertices[3] = vertices[6];
00862             new_cell.vertices(cellvertices);
00863             segment_out.push_back(new_cell);
00864           }
00865         }
00866         else //split edge 12 again, introduce lines 43 and 53
00867         {
00868           if (stable_line_is_longer(vertices[1], vertices[3],
00869                                     vertices[0], vertices[1])) //split edge 13 again, introduce line 60
00870           {
00871             assert(false && "diag13_len > diag01_len impossible!");
00872           }
00873           else  //split edge 01 again, introduce line 43
00874           {
00875             //std::cout << "apply_3_1_2" << std::endl;
00876             cellvertices[0] = vertices[0];
00877             cellvertices[1] = vertices[4];
00878             cellvertices[2] = vertices[2];
00879             cellvertices[3] = vertices[3];
00880             new_cell.vertices(cellvertices);
00881             segment_out.push_back(new_cell);
00882             
00883             cellvertices[0] = vertices[4];
00884             cellvertices[1] = vertices[5];
00885             cellvertices[2] = vertices[2];
00886             cellvertices[3] = vertices[3];
00887             new_cell.vertices(cellvertices);
00888             segment_out.push_back(new_cell);
00889             
00890             cellvertices[0] = vertices[4];
00891             cellvertices[1] = vertices[6];
00892             cellvertices[2] = vertices[5];
00893             cellvertices[3] = vertices[3];
00894             new_cell.vertices(cellvertices);
00895             segment_out.push_back(new_cell);
00896           }
00897         }
00898       }
00899       else //split edge 12, introduce line 05
00900       {
00901         if (stable_line_is_longer(vertices[1], vertices[3],
00902                                   vertices[1], vertices[2])) //split edge 13 again, introduce line 62
00903         {
00904           if (stable_line_is_longer(vertices[1], vertices[3],
00905                                     vertices[0], vertices[1])) //split edge 13 again, introduce line 60
00906           {
00907             //std::cout << "apply_3_1_3" << std::endl;
00908             cellvertices[0] = vertices[0];
00909             cellvertices[1] = vertices[4];
00910             cellvertices[2] = vertices[5];
00911             cellvertices[3] = vertices[6];
00912             new_cell.vertices(cellvertices);
00913             segment_out.push_back(new_cell);
00914             
00915             cellvertices[0] = vertices[0];
00916             cellvertices[1] = vertices[6];
00917             cellvertices[2] = vertices[5];
00918             cellvertices[3] = vertices[2];
00919             new_cell.vertices(cellvertices);
00920             segment_out.push_back(new_cell);
00921             
00922             cellvertices[0] = vertices[0];
00923             cellvertices[1] = vertices[6];
00924             cellvertices[2] = vertices[2];
00925             cellvertices[3] = vertices[3];
00926             new_cell.vertices(cellvertices);
00927             segment_out.push_back(new_cell);
00928           }
00929           else  //split edge 01 again, introduce line 43
00930           {
00931             assert(false && "diag13_len > diag01_len impossible!");
00932           }
00933         }
00934         else //split edge 12 again, introduce line 53
00935         {
00936           if (stable_line_is_longer(vertices[1], vertices[3],
00937                                     vertices[0], vertices[1])) //split edge 13 again, introduce line 60
00938           {
00939             cellvertices[0] = vertices[0];
00940             cellvertices[1] = vertices[4];
00941             cellvertices[2] = vertices[5];
00942             cellvertices[3] = vertices[6];
00943             new_cell.vertices(cellvertices);
00944             segment_out.push_back(new_cell);
00945             
00946             cellvertices[0] = vertices[0];
00947             cellvertices[1] = vertices[6];
00948             cellvertices[2] = vertices[5];
00949             cellvertices[3] = vertices[3];
00950             new_cell.vertices(cellvertices);
00951             segment_out.push_back(new_cell);
00952             
00953             cellvertices[0] = vertices[0];
00954             cellvertices[1] = vertices[5];
00955             cellvertices[2] = vertices[2];
00956             cellvertices[3] = vertices[3];
00957             new_cell.vertices(cellvertices);
00958             segment_out.push_back(new_cell);
00959           }
00960           else  //split edge 01 again, introduce line 43
00961           {
00962             //std::cout << "apply_3_1_4" << std::endl;
00963             cellvertices[0] = vertices[0];
00964             cellvertices[1] = vertices[4];
00965             cellvertices[2] = vertices[5];
00966             cellvertices[3] = vertices[3];
00967             new_cell.vertices(cellvertices);
00968             segment_out.push_back(new_cell);
00969             
00970             cellvertices[0] = vertices[4];
00971             cellvertices[1] = vertices[5];
00972             cellvertices[2] = vertices[3];
00973             cellvertices[3] = vertices[6];
00974             new_cell.vertices(cellvertices);
00975             segment_out.push_back(new_cell);
00976             
00977             cellvertices[0] = vertices[0];
00978             cellvertices[1] = vertices[5];
00979             cellvertices[2] = vertices[2];
00980             cellvertices[3] = vertices[3];
00981             new_cell.vertices(cellvertices);
00982             segment_out.push_back(new_cell);
00983           }
00984         }
00985       }
00986       
00987     }
00988 
00989 
01004     template <typename CellType, typename DomainTypeOut, typename VertexType>
01005     static void apply3_2(CellType const & cell_in, 
01006                          DomainTypeOut & segment_out,
01007                          VertexType ** vertices
01008                         )
01009     {
01010       CellType new_cell;
01011       VertexType * cellvertices[topology::bndcells<tetrahedron_tag, 0>::num];
01012 
01013       //std::cout << "apply_3_2" << std::endl;
01014       
01015       cellvertices[0] = vertices[0];
01016       cellvertices[1] = vertices[4];
01017       cellvertices[2] = vertices[6];
01018       cellvertices[3] = vertices[3];
01019       new_cell.vertices(cellvertices);
01020       segment_out.push_back(new_cell);
01021       
01022       cellvertices[0] = vertices[4];
01023       cellvertices[1] = vertices[5];
01024       cellvertices[2] = vertices[6];
01025       cellvertices[3] = vertices[3];
01026       new_cell.vertices(cellvertices);
01027       segment_out.push_back(new_cell);
01028       
01029       cellvertices[0] = vertices[4];
01030       cellvertices[1] = vertices[1];
01031       cellvertices[2] = vertices[5];
01032       cellvertices[3] = vertices[3];
01033       new_cell.vertices(cellvertices);
01034       segment_out.push_back(new_cell);
01035       
01036       cellvertices[0] = vertices[6];
01037       cellvertices[1] = vertices[5];
01038       cellvertices[2] = vertices[2];
01039       cellvertices[3] = vertices[3];
01040       new_cell.vertices(cellvertices);
01041       segment_out.push_back(new_cell);
01042       
01043     }
01044 
01045 
01059     template <typename CellType, typename DomainTypeOut, typename VertexType>
01060     static void apply3_3(CellType const & cell_in, 
01061                          DomainTypeOut & segment_out,
01062                          VertexType ** vertices
01063                         )
01064     {
01065       //std::cout << "Found!" << std::endl;
01066       CellType new_cell;
01067       VertexType * cellvertices[topology::bndcells<tetrahedron_tag, 0>::num];
01068 
01069       //double diag01_len = viennagrid::norm(vertices[0]->point() - vertices[1]->point());
01070       //double diag12_len = viennagrid::norm(vertices[1]->point() - vertices[2]->point());
01071       //double diag03_len = viennagrid::norm(vertices[0]->point() - vertices[3]->point());
01072 
01073       // Strategy: The two longest edges of the common vertex are split 'twice',
01074       //           i.e. two new edges start from the center of the two longest edges
01075       //if (diag01_len > diag03_len) //split edge 01 again, introduce line 43
01076       if (stable_line_is_longer(vertices[0], vertices[1],
01077                                 vertices[0], vertices[3])) //split edge 01 again, introduce line 43
01078       {
01079         cellvertices[0] = vertices[4];
01080         cellvertices[1] = vertices[1];
01081         cellvertices[2] = vertices[5];
01082         cellvertices[3] = vertices[3];
01083         new_cell.vertices(cellvertices);
01084         segment_out.push_back(new_cell);
01085           
01086         //if (diag01_len > diag12_len) //split edge 01 again, introduce line 42
01087         if (stable_line_is_longer(vertices[0], vertices[1],
01088                                   vertices[1], vertices[2])) //split edge 01 again, introduce line 42
01089         {
01090           //std::cout << "apply_3_3_1" << std::endl;
01091           cellvertices[0] = vertices[0];
01092           cellvertices[1] = vertices[4];
01093           cellvertices[2] = vertices[2];
01094           cellvertices[3] = vertices[6];
01095           new_cell.vertices(cellvertices);
01096           segment_out.push_back(new_cell);
01097           
01098           cellvertices[0] = vertices[6];
01099           cellvertices[1] = vertices[4];
01100           cellvertices[2] = vertices[2];
01101           cellvertices[3] = vertices[3];
01102           new_cell.vertices(cellvertices);
01103           segment_out.push_back(new_cell);
01104           
01105           cellvertices[0] = vertices[4];
01106           cellvertices[1] = vertices[5];
01107           cellvertices[2] = vertices[2];
01108           cellvertices[3] = vertices[3];
01109           new_cell.vertices(cellvertices);
01110           segment_out.push_back(new_cell);
01111         }
01112         else //split edge 12 again, introduce line 50
01113         {
01114           //std::cout << "apply_3_3_2" << std::endl;
01115           cellvertices[0] = vertices[0];
01116           cellvertices[1] = vertices[4];
01117           cellvertices[2] = vertices[5];
01118           cellvertices[3] = vertices[6];
01119           new_cell.vertices(cellvertices);
01120           segment_out.push_back(new_cell);
01121           
01122           cellvertices[0] = vertices[0];
01123           cellvertices[1] = vertices[5];
01124           cellvertices[2] = vertices[2];
01125           cellvertices[3] = vertices[6];
01126           new_cell.vertices(cellvertices);
01127           segment_out.push_back(new_cell);
01128           
01129           cellvertices[0] = vertices[6];
01130           cellvertices[1] = vertices[4];
01131           cellvertices[2] = vertices[5];
01132           cellvertices[3] = vertices[3];
01133           new_cell.vertices(cellvertices);
01134           segment_out.push_back(new_cell);
01135           
01136           cellvertices[0] = vertices[6];
01137           cellvertices[1] = vertices[5];
01138           cellvertices[2] = vertices[2];
01139           cellvertices[3] = vertices[3];
01140           new_cell.vertices(cellvertices);
01141           segment_out.push_back(new_cell);
01142         }
01143       }
01144       else  //split edge 03 again, introduce line 61
01145       {
01146         cellvertices[0] = vertices[4];
01147         cellvertices[1] = vertices[1];
01148         cellvertices[2] = vertices[5];
01149         cellvertices[3] = vertices[6];
01150         new_cell.vertices(cellvertices);
01151         segment_out.push_back(new_cell);
01152           
01153         cellvertices[0] = vertices[6];
01154         cellvertices[1] = vertices[1];
01155         cellvertices[2] = vertices[5];
01156         cellvertices[3] = vertices[3];
01157         new_cell.vertices(cellvertices);
01158         segment_out.push_back(new_cell);
01159           
01160         cellvertices[0] = vertices[6];
01161         cellvertices[1] = vertices[5];
01162         cellvertices[2] = vertices[2];
01163         cellvertices[3] = vertices[3];
01164         new_cell.vertices(cellvertices);
01165         segment_out.push_back(new_cell);
01166           
01167         //if (diag01_len > diag12_len) //split edge 01 again, introduce line 42
01168         if (stable_line_is_longer(vertices[0], vertices[1],
01169                                   vertices[1], vertices[2])) //split edge 01 again, introduce line 42
01170         {
01171           //std::cout << "apply_3_3_3" << std::endl;
01172           cellvertices[0] = vertices[0];
01173           cellvertices[1] = vertices[4];
01174           cellvertices[2] = vertices[2];
01175           cellvertices[3] = vertices[6];
01176           new_cell.vertices(cellvertices);
01177           segment_out.push_back(new_cell);
01178           
01179           cellvertices[0] = vertices[6];
01180           cellvertices[1] = vertices[4];
01181           cellvertices[2] = vertices[5];
01182           cellvertices[3] = vertices[2];
01183           new_cell.vertices(cellvertices);
01184           segment_out.push_back(new_cell);
01185         }
01186         else //split edge 12 again, introduce line 50
01187         {
01188           //std::cout << "apply_3_3_4" << std::endl;
01189           cellvertices[0] = vertices[0];
01190           cellvertices[1] = vertices[4];
01191           cellvertices[2] = vertices[5];
01192           cellvertices[3] = vertices[6];
01193           new_cell.vertices(cellvertices);
01194           segment_out.push_back(new_cell);
01195           
01196           cellvertices[0] = vertices[0];
01197           cellvertices[1] = vertices[5];
01198           cellvertices[2] = vertices[2];
01199           cellvertices[3] = vertices[6];
01200           new_cell.vertices(cellvertices);
01201           segment_out.push_back(new_cell);
01202         }
01203       }
01204       
01205     }
01206 
01220     template <typename CellType, typename DomainTypeOut, typename VertexType>
01221     static void apply3_4(CellType const & cell_in, 
01222                          DomainTypeOut & segment_out,
01223                          VertexType ** vertices
01224                         )
01225     {
01226       CellType new_cell;
01227       VertexType * cellvertices[topology::bndcells<tetrahedron_tag, 0>::num];
01228 
01229       //double diag01_len = viennagrid::norm(vertices[0]->point() - vertices[1]->point());
01230       //double diag12_len = viennagrid::norm(vertices[1]->point() - vertices[2]->point());
01231       //double diag23_len = viennagrid::norm(vertices[2]->point() - vertices[3]->point());
01232 
01233       // Strategy: The two longest edges of the common vertex are split 'twice',
01234       //           i.e. two new edges start from the center of the two longest edges
01235       //if (diag01_len > diag12_len) //split edge 01 again, introduce line 42
01236       if (stable_line_is_longer(vertices[0], vertices[1],
01237                                 vertices[1], vertices[2])) //split edge 01 again, introduce line 42
01238       {
01239         cellvertices[0] = vertices[0];
01240         cellvertices[1] = vertices[4];
01241         cellvertices[2] = vertices[2];
01242         cellvertices[3] = vertices[6];
01243         new_cell.vertices(cellvertices);
01244         segment_out.push_back(new_cell);
01245       
01246         cellvertices[0] = vertices[0];
01247         cellvertices[1] = vertices[4];
01248         cellvertices[2] = vertices[6];
01249         cellvertices[3] = vertices[3];
01250         new_cell.vertices(cellvertices);
01251         segment_out.push_back(new_cell);
01252         
01253         cellvertices[0] = vertices[4];
01254         cellvertices[1] = vertices[5];
01255         cellvertices[2] = vertices[2];
01256         cellvertices[3] = vertices[6];
01257         new_cell.vertices(cellvertices);
01258         segment_out.push_back(new_cell);
01259           
01260         //if (diag12_len > diag23_len) //split edge 12 again, introduce line 53
01261         if (stable_line_is_longer(vertices[1], vertices[2],
01262                                   vertices[2], vertices[3])) //split edge 12 again, introduce line 53
01263         {
01264           //std::cout << "apply_3_4_1" << std::endl;
01265           cellvertices[0] = vertices[4];
01266           cellvertices[1] = vertices[1];
01267           cellvertices[2] = vertices[5];
01268           cellvertices[3] = vertices[3];
01269           new_cell.vertices(cellvertices);
01270           segment_out.push_back(new_cell);
01271           
01272           cellvertices[0] = vertices[4];
01273           cellvertices[1] = vertices[5];
01274           cellvertices[2] = vertices[6];
01275           cellvertices[3] = vertices[3];
01276           new_cell.vertices(cellvertices);
01277           segment_out.push_back(new_cell);
01278         }
01279         else //split edge 23 again, introduce line 61
01280         {
01281           //std::cout << "apply_3_4_2" << std::endl;
01282           cellvertices[0] = vertices[4];
01283           cellvertices[1] = vertices[1];
01284           cellvertices[2] = vertices[6];
01285           cellvertices[3] = vertices[3];
01286           new_cell.vertices(cellvertices);
01287           segment_out.push_back(new_cell);
01288           
01289           cellvertices[0] = vertices[4];
01290           cellvertices[1] = vertices[1];
01291           cellvertices[2] = vertices[5];
01292           cellvertices[3] = vertices[6];
01293           new_cell.vertices(cellvertices);
01294           segment_out.push_back(new_cell);
01295         }
01296       }
01297       else //split edge 12, introduce line 50
01298       {
01299         //if (diag12_len > diag23_len) //split edge 12 again, introduce line 53
01300         if (stable_line_is_longer(vertices[1], vertices[2],
01301                                   vertices[2], vertices[3])) //split edge 12 again, introduce line 53
01302         {
01303           //std::cout << "apply_3_4_3" << std::endl;
01304           cellvertices[0] = vertices[0];
01305           cellvertices[1] = vertices[4];
01306           cellvertices[2] = vertices[5];
01307           cellvertices[3] = vertices[3];
01308           new_cell.vertices(cellvertices);
01309           segment_out.push_back(new_cell);
01310           
01311           cellvertices[0] = vertices[0];
01312           cellvertices[1] = vertices[5];
01313           cellvertices[2] = vertices[6];
01314           cellvertices[3] = vertices[3];
01315           new_cell.vertices(cellvertices);
01316           segment_out.push_back(new_cell);
01317           
01318           cellvertices[0] = vertices[0];
01319           cellvertices[1] = vertices[5];
01320           cellvertices[2] = vertices[2];
01321           cellvertices[3] = vertices[6];
01322           new_cell.vertices(cellvertices);
01323           segment_out.push_back(new_cell);
01324             
01325           cellvertices[0] = vertices[4];
01326           cellvertices[1] = vertices[1];
01327           cellvertices[2] = vertices[5];
01328           cellvertices[3] = vertices[3];
01329           new_cell.vertices(cellvertices);
01330           segment_out.push_back(new_cell);
01331         }
01332         else //split edge 23 again, introduce line 61
01333         {
01334           //std::cout << "apply_3_4_4" << std::endl;
01335           cellvertices[0] = vertices[0];
01336           cellvertices[1] = vertices[4];
01337           cellvertices[2] = vertices[5];
01338           cellvertices[3] = vertices[6];
01339           new_cell.vertices(cellvertices);
01340           segment_out.push_back(new_cell);
01341           
01342           cellvertices[0] = vertices[0];
01343           cellvertices[1] = vertices[4];
01344           cellvertices[2] = vertices[6];
01345           cellvertices[3] = vertices[3];
01346           new_cell.vertices(cellvertices);
01347           segment_out.push_back(new_cell);
01348           
01349           cellvertices[0] = vertices[0];
01350           cellvertices[1] = vertices[5];
01351           cellvertices[2] = vertices[2];
01352           cellvertices[3] = vertices[6];
01353           new_cell.vertices(cellvertices);
01354           segment_out.push_back(new_cell);
01355           
01356           cellvertices[0] = vertices[4];
01357           cellvertices[1] = vertices[1];
01358           cellvertices[2] = vertices[5];
01359           cellvertices[3] = vertices[6];
01360           new_cell.vertices(cellvertices);
01361           segment_out.push_back(new_cell);
01362           
01363           cellvertices[0] = vertices[4];
01364           cellvertices[1] = vertices[1];
01365           cellvertices[2] = vertices[6];
01366           cellvertices[3] = vertices[3];
01367           new_cell.vertices(cellvertices);
01368           segment_out.push_back(new_cell);          
01369         }
01370       }
01371       
01372     }
01373 
01374 
01375 
01377     template <typename CellType, typename DomainTypeOut>
01378     static void apply3(CellType const & cell_in, DomainTypeOut & segment_out)
01379     {
01380       typedef typename CellType::config_type        ConfigTypeIn;
01381       typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type            VertexOnCellRange;
01382       typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type         VertexOnCellIterator;            
01383       typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type            EdgeOnCellRange;
01384       typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type           EdgeOnCellIterator;            
01385       
01386       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type             VertexType;
01387       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 1>::type             EdgeType;
01388 
01389       VertexType * vertices[topology::bndcells<tetrahedron_tag, 0>::num];
01390       //std::cout << "apply3()" << std::endl;
01391       
01392       //
01393       // Step 1: Get vertices from input cell
01394       //
01395       VertexOnCellRange vertices_on_cell = viennagrid::ncells<0>(cell_in);
01396       VertexOnCellIterator vocit = vertices_on_cell.begin();
01397       vertices[0] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
01398       vertices[1] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
01399       vertices[2] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
01400       vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]);
01401       
01402 
01403       //
01404       // Step 2: Bring vertices in correct order
01405       //
01406       VertexType * ordered_vertices[topology::bndcells<tetrahedron_tag, 0>::num + 3];
01407       EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(cell_in);
01408       EdgeOnCellIterator eocit = edges_on_cell.begin();
01409       EdgeType const & e0 = *eocit; ++eocit;
01410       EdgeType const & e1 = *eocit; ++eocit;
01411       EdgeType const & e2 = *eocit; ++eocit;
01412       EdgeType const & e3 = *eocit; ++eocit;
01413       EdgeType const & e4 = *eocit; ++eocit;
01414       EdgeType const & e5 = *eocit;
01415       
01416       //with e0
01417       if (viennadata::access<refinement_key, bool>(refinement_key())(e0) == true)
01418       {
01419         if (viennadata::access<refinement_key, bool>(refinement_key())(e1) == true)
01420         {
01421           ordered_vertices[0] = vertices[2];
01422           ordered_vertices[1] = vertices[0];
01423           ordered_vertices[2] = vertices[1];
01424           ordered_vertices[3] = vertices[3];
01425           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
01426           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
01427           
01428           if (viennadata::access<refinement_key, bool>(refinement_key())(e2) == true)
01429           {
01430             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
01431             apply3_1(cell_in, segment_out, ordered_vertices);
01432           }
01433           else if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == true)
01434           {
01435             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
01436             apply3_2(cell_in, segment_out, ordered_vertices);
01437           }
01438           else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true)
01439           {
01440             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
01441             apply3_4(cell_in, segment_out, ordered_vertices);
01442           }
01443           else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true)
01444           {
01445             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
01446             apply3_3(cell_in, segment_out, ordered_vertices);
01447           }
01448           else
01449           {
01450             assert(false && "Logic error: No edge for refinement found!"); 
01451           }
01452         }
01453         else if (viennadata::access<refinement_key, bool>(refinement_key())(e2) == true)
01454         {
01455           ordered_vertices[0] = vertices[1];
01456           ordered_vertices[1] = vertices[0];
01457           ordered_vertices[2] = vertices[3];
01458           ordered_vertices[3] = vertices[2];
01459           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
01460           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
01461           
01462           if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == true)
01463           {
01464             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
01465             apply3_3(cell_in, segment_out, ordered_vertices);
01466           }
01467           else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true)
01468           {
01469             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
01470             apply3_2(cell_in, segment_out, ordered_vertices);
01471           }
01472           else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true)
01473           {
01474             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
01475             apply3_4(cell_in, segment_out, ordered_vertices);
01476           }
01477           else
01478           {
01479             assert(false && "Logic error: No edge for refinement found!"); 
01480           }
01481         }
01482         else if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == true)
01483         {
01484           ordered_vertices[0] = vertices[0];
01485           ordered_vertices[1] = vertices[1];
01486           ordered_vertices[2] = vertices[2];
01487           ordered_vertices[3] = vertices[3];
01488           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
01489           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
01490           
01491           if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true)
01492           {
01493             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
01494             apply3_1(cell_in, segment_out, ordered_vertices);
01495           }
01496           else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true)
01497           {
01498             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
01499             apply3_4(cell_in, segment_out, ordered_vertices);
01500           }
01501           else
01502           {
01503             assert(false && "Logic error: No edge for refinement found!"); 
01504           }
01505         }
01506         else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true)
01507         {
01508           ordered_vertices[0] = vertices[3];
01509           ordered_vertices[1] = vertices[1];
01510           ordered_vertices[2] = vertices[0];
01511           ordered_vertices[3] = vertices[2];
01512           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
01513           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
01514           
01515           if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true)
01516           {
01517             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
01518             apply3_3(cell_in, segment_out, ordered_vertices);
01519           }
01520           else
01521           {
01522             assert(false && "Logic error: No edge for refinement found!"); 
01523           }
01524         }
01525         else //no second edge
01526         {
01527           assert(false && "Logic error: No edge for refinement found!"); 
01528         }
01529       } 
01530       else if (viennadata::access<refinement_key, bool>(refinement_key())(e1) == true)
01531       {
01532         if (viennadata::access<refinement_key, bool>(refinement_key())(e2) == true)
01533         {
01534           ordered_vertices[0] = vertices[3];
01535           ordered_vertices[1] = vertices[0];
01536           ordered_vertices[2] = vertices[2];
01537           ordered_vertices[3] = vertices[1];
01538           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
01539           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
01540           
01541           if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == true)
01542           {
01543             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
01544             apply3_4(cell_in, segment_out, ordered_vertices);
01545           }
01546           else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true)
01547           {
01548             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
01549             apply3_3(cell_in, segment_out, ordered_vertices);
01550           }
01551           else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true)
01552           {
01553             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
01554             apply3_2(cell_in, segment_out, ordered_vertices);
01555           }
01556           else
01557           {
01558             assert(false && "Logic error: No edge for refinement found!"); 
01559           }
01560         }
01561         else if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == true)
01562         {
01563           ordered_vertices[0] = vertices[1];
01564           ordered_vertices[1] = vertices[2];
01565           ordered_vertices[2] = vertices[0];
01566           ordered_vertices[3] = vertices[3];
01567           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
01568           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
01569           
01570           if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true)
01571           {
01572             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
01573             apply3_3(cell_in, segment_out, ordered_vertices);
01574           }
01575           else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true)
01576           {
01577             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
01578             apply3_1(cell_in, segment_out, ordered_vertices);
01579           }
01580           else
01581           {
01582             assert(false && "Logic error: No edge for refinement found!"); 
01583           }
01584         }
01585         else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true)
01586         {
01587           if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true)
01588           {
01589             //make edges 4 and 5 the references
01590             ordered_vertices[0] = vertices[1];
01591             ordered_vertices[1] = vertices[3];
01592             ordered_vertices[2] = vertices[2];
01593             ordered_vertices[3] = vertices[0];
01594             ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
01595             ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
01596             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
01597             
01598             apply3_4(cell_in, segment_out, ordered_vertices);
01599           }
01600           else
01601           {
01602             assert(false && "Logic error: No edge for refinement found!"); 
01603           }
01604         }
01605         else //no second edge
01606         {
01607           assert(false && "Logic error: No edge for refinement found!"); 
01608         }
01609       }
01610       else if (viennadata::access<refinement_key, bool>(refinement_key())(e2) == true)
01611       {
01612         if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == true)
01613         {
01614           //NOTE: edges 2 and 3 don't have a common vertex, therefore 'base facet' is chosen depending on the third edge
01615           
01616           if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true)
01617           {
01618             // take edges 2 and 4 as reference
01619             ordered_vertices[0] = vertices[0];
01620             ordered_vertices[1] = vertices[3];
01621             ordered_vertices[2] = vertices[1];
01622             ordered_vertices[3] = vertices[2];
01623             ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
01624             ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
01625             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
01626             
01627             apply3_4(cell_in, segment_out, ordered_vertices);
01628           }
01629           else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true)
01630           {
01631             // take edges 5 and 3 as reference
01632             ordered_vertices[0] = vertices[3];
01633             ordered_vertices[1] = vertices[2];
01634             ordered_vertices[2] = vertices[1];
01635             ordered_vertices[3] = vertices[0];
01636             ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
01637             ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
01638             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
01639             
01640             apply3_3(cell_in, segment_out, ordered_vertices);
01641           }
01642           else
01643           {
01644             assert(false && "Logic error: No edge for refinement found!"); 
01645           }
01646         }
01647         else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true)
01648         {
01649           ordered_vertices[0] = vertices[0];
01650           ordered_vertices[1] = vertices[3];
01651           ordered_vertices[2] = vertices[1];
01652           ordered_vertices[3] = vertices[2];
01653           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
01654           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
01655           
01656           if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true)
01657           {
01658             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
01659             apply3_1(cell_in, segment_out, ordered_vertices);
01660           }
01661           else
01662           {
01663             assert(false && "Logic error: No edge for refinement found!"); 
01664           }
01665         }
01666         else //no second edge
01667         {
01668           assert(false && "Logic error: No edge for refinement found!"); 
01669         }
01670       }
01671       else if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == true)
01672       {
01673         if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == true)
01674         {
01675           ordered_vertices[0] = vertices[2];
01676           ordered_vertices[1] = vertices[1];
01677           ordered_vertices[2] = vertices[3];
01678           ordered_vertices[3] = vertices[0];
01679           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
01680           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
01681           
01682           if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == true)
01683           {
01684             ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
01685             apply3_2(cell_in, segment_out, ordered_vertices);
01686           }
01687           else
01688           {
01689             assert(false && "Logic error: No edge for refinement found!"); 
01690           }
01691         }
01692         else //no second edge
01693         {
01694           assert(false && "Logic error: No edge for refinement found!"); 
01695         }
01696       }
01697       else //no first edge found
01698       {
01699         assert(false && "Logic error: No edge for refinement found!"); 
01700       }
01701       
01702     } //apply3()
01703 
01704 
01705 
01706 
01708 
01722     template <typename CellType, typename DomainTypeOut, typename VertexType>
01723     static void apply4_1(CellType const & cell_in, 
01724                          DomainTypeOut & segment_out,
01725                          VertexType ** vertices
01726                         )
01727     {
01728       CellType new_cell;
01729       VertexType * cellvertices[topology::bndcells<tetrahedron_tag, 0>::num];
01730 
01731       //double diag03_len = viennagrid::norm(vertices[0]->point() - vertices[3]->point());
01732       //double diag13_len = viennagrid::norm(vertices[1]->point() - vertices[3]->point());
01733       //double diag23_len = viennagrid::norm(vertices[2]->point() - vertices[3]->point());
01734 
01735       //if (diag03_len > diag13_len) //split edge 03, introduce line 71
01736       if (stable_line_is_longer(vertices[0], vertices[3],
01737                                 vertices[1], vertices[3])) //split edge 03, introduce line 71
01738       {
01739         cellvertices[0] = vertices[0];
01740         cellvertices[1] = vertices[1];
01741         cellvertices[2] = vertices[4];
01742         cellvertices[3] = vertices[7];
01743         new_cell.vertices(cellvertices);
01744         segment_out.push_back(new_cell);
01745           
01746         cellvertices[0] = vertices[7];
01747         cellvertices[1] = vertices[5];
01748         cellvertices[2] = vertices[6];
01749         cellvertices[3] = vertices[3];
01750         new_cell.vertices(cellvertices);
01751         segment_out.push_back(new_cell);
01752         
01753         //if (diag13_len > diag23_len) //split edge 13, introduce line 52
01754         if (stable_line_is_longer(vertices[1], vertices[3],
01755                                   vertices[2], vertices[3])) //split edge 13, introduce line 52
01756         {
01757           //std::cout << "apply_4_1_1" << std::endl;
01758           cellvertices[0] = vertices[7];
01759           cellvertices[1] = vertices[1];
01760           cellvertices[2] = vertices[4];
01761           cellvertices[3] = vertices[5];
01762           new_cell.vertices(cellvertices);
01763           segment_out.push_back(new_cell);
01764 
01765           cellvertices[0] = vertices[7];
01766           cellvertices[1] = vertices[5];
01767           cellvertices[2] = vertices[4];
01768           cellvertices[3] = vertices[6];
01769           new_cell.vertices(cellvertices);
01770           segment_out.push_back(new_cell);
01771           
01772           cellvertices[0] = vertices[1];
01773           cellvertices[1] = vertices[2];
01774           cellvertices[2] = vertices[4];
01775           cellvertices[3] = vertices[5];
01776           new_cell.vertices(cellvertices);
01777           segment_out.push_back(new_cell);
01778           
01779           cellvertices[0] = vertices[4];
01780           cellvertices[1] = vertices[5];
01781           cellvertices[2] = vertices[2];
01782           cellvertices[3] = vertices[6];
01783           new_cell.vertices(cellvertices);
01784           segment_out.push_back(new_cell);
01785         }
01786         else //split edge 23, introduce line 61
01787         {
01788           //std::cout << "apply_4_1_2" << std::endl;
01789           cellvertices[0] = vertices[7];
01790           cellvertices[1] = vertices[1];
01791           cellvertices[2] = vertices[6];
01792           cellvertices[3] = vertices[5];
01793           new_cell.vertices(cellvertices);
01794           segment_out.push_back(new_cell);
01795           
01796           cellvertices[0] = vertices[7];
01797           cellvertices[1] = vertices[1];
01798           cellvertices[2] = vertices[4];
01799           cellvertices[3] = vertices[6];
01800           new_cell.vertices(cellvertices);
01801           segment_out.push_back(new_cell);
01802 
01803           cellvertices[0] = vertices[1];
01804           cellvertices[1] = vertices[2];
01805           cellvertices[2] = vertices[4];
01806           cellvertices[3] = vertices[6];
01807           new_cell.vertices(cellvertices);
01808           segment_out.push_back(new_cell);
01809         }
01810       }
01811       else //split edge 13, introduce line 50
01812       {
01813         cellvertices[0] = vertices[0];
01814         cellvertices[1] = vertices[1];
01815         cellvertices[2] = vertices[4];
01816         cellvertices[3] = vertices[5];
01817         new_cell.vertices(cellvertices);
01818         segment_out.push_back(new_cell);
01819           
01820         cellvertices[0] = vertices[0];
01821         cellvertices[1] = vertices[5];
01822         cellvertices[2] = vertices[4];
01823         cellvertices[3] = vertices[7];
01824         new_cell.vertices(cellvertices);
01825         segment_out.push_back(new_cell);
01826 
01827         cellvertices[0] = vertices[7];
01828         cellvertices[1] = vertices[5];
01829         cellvertices[2] = vertices[6];
01830         cellvertices[3] = vertices[3];
01831         new_cell.vertices(cellvertices);
01832         segment_out.push_back(new_cell);
01833         
01834         cellvertices[0] = vertices[7];
01835         cellvertices[1] = vertices[5];
01836         cellvertices[2] = vertices[4];
01837         cellvertices[3] = vertices[6];
01838         new_cell.vertices(cellvertices);
01839         segment_out.push_back(new_cell);
01840 
01841         //if (diag13_len > diag23_len) //split edge 13 again, introduce line 52
01842         if (stable_line_is_longer(vertices[1], vertices[3],
01843                                   vertices[2], vertices[3])) //split edge 13 again, introduce line 52
01844         {
01845           //std::cout << "apply_4_1_3" << std::endl;
01846           cellvertices[0] = vertices[1];
01847           cellvertices[1] = vertices[2];
01848           cellvertices[2] = vertices[4];
01849           cellvertices[3] = vertices[5];
01850           new_cell.vertices(cellvertices);
01851           segment_out.push_back(new_cell);
01852           
01853           cellvertices[0] = vertices[4];
01854           cellvertices[1] = vertices[5];
01855           cellvertices[2] = vertices[2];
01856           cellvertices[3] = vertices[6];
01857           new_cell.vertices(cellvertices);
01858           segment_out.push_back(new_cell);
01859         }
01860         else //split edge 23, introduce line 61
01861         {
01862           //std::cout << "apply_4_1_4" << std::endl;
01863           cellvertices[0] = vertices[5];
01864           cellvertices[1] = vertices[1];
01865           cellvertices[2] = vertices[4];
01866           cellvertices[3] = vertices[6];
01867           new_cell.vertices(cellvertices);
01868           segment_out.push_back(new_cell);
01869 
01870           cellvertices[0] = vertices[4];
01871           cellvertices[1] = vertices[1];
01872           cellvertices[2] = vertices[2];
01873           cellvertices[3] = vertices[6];
01874           new_cell.vertices(cellvertices);
01875           segment_out.push_back(new_cell);
01876         }
01877       }
01878       
01879     }
01880     
01895     template <typename CellType, typename DomainTypeOut, typename VertexType>
01896     static void apply4_2(CellType const & cell_in, 
01897                          DomainTypeOut & segment_out,
01898                          VertexType ** vertices
01899                         )
01900     {
01901       CellType new_cell;
01902       VertexType * cellvertices[topology::bndcells<tetrahedron_tag, 0>::num];
01903 
01904       //double diag02_len = viennagrid::norm(vertices[0]->point() - vertices[2]->point());
01905       //double diag03_len = viennagrid::norm(vertices[0]->point() - vertices[3]->point());
01906       //double diag12_len = viennagrid::norm(vertices[1]->point() - vertices[2]->point());
01907       //double diag13_len = viennagrid::norm(vertices[1]->point() - vertices[3]->point());
01908 
01909       //if (diag03_len > diag13_len) //split edge 03, introduce line 61
01910       if (stable_line_is_longer(vertices[0], vertices[3],
01911                                 vertices[1], vertices[3])) //split edge 03, introduce line 61
01912       {
01913         //std::cout << "split!" << std::endl;
01914         //if (diag13_len > diag12_len) //split edge 13, introduce line 72
01915         if (stable_line_is_longer(vertices[1], vertices[3],
01916                                   vertices[1], vertices[2])) //split edge 13, introduce line 72
01917         {
01918           //if (diag02_len > diag03_len) //split edge 02, introduce line 53
01919           if (stable_line_is_longer(vertices[0], vertices[2],
01920                                     vertices[0], vertices[3])) //split edge 02, introduce line 53
01921           {
01922             //if (diag02_len > diag12_len) //split edge 02, introduce line 51
01923             if (stable_line_is_longer(vertices[0], vertices[2],
01924                                       vertices[1], vertices[2])) //split edge 02, introduce line 51
01925             {
01926               //std::cout << "apply_4_2_1_" << std::endl;
01927               cellvertices[0] = vertices[0];
01928               cellvertices[1] = vertices[1];
01929               cellvertices[2] = vertices[5];
01930               cellvertices[3] = vertices[6];
01931               new_cell.vertices(cellvertices);
01932               segment_out.push_back(new_cell);
01933 
01934               cellvertices[0] = vertices[6];
01935               cellvertices[1] = vertices[1];
01936               cellvertices[2] = vertices[5];
01937               cellvertices[3] = vertices[7];
01938               new_cell.vertices(cellvertices);
01939               segment_out.push_back(new_cell);
01940 
01941               cellvertices[0] = vertices[1];
01942               cellvertices[1] = vertices[4];
01943               cellvertices[2] = vertices[5];
01944               cellvertices[3] = vertices[7];
01945               new_cell.vertices(cellvertices);
01946               segment_out.push_back(new_cell);
01947 
01948               cellvertices[0] = vertices[4];
01949               cellvertices[1] = vertices[2];
01950               cellvertices[2] = vertices[5];
01951               cellvertices[3] = vertices[7];
01952               new_cell.vertices(cellvertices);
01953               segment_out.push_back(new_cell);
01954 
01955               cellvertices[0] = vertices[6];
01956               cellvertices[1] = vertices[7];
01957               cellvertices[2] = vertices[5];
01958               cellvertices[3] = vertices[3];
01959               new_cell.vertices(cellvertices);
01960               segment_out.push_back(new_cell);
01961 
01962               cellvertices[0] = vertices[7];
01963               cellvertices[1] = vertices[2];
01964               cellvertices[2] = vertices[5];
01965               cellvertices[3] = vertices[3];
01966               new_cell.vertices(cellvertices);
01967               segment_out.push_back(new_cell);
01968             }
01969             else //split edge 12, introduce line 40
01970             {
01971               assert( false && "Logic error: diag02 < diag12 not possible here!");
01972             }
01973           }
01974           else //split edge 03, introduce line 62
01975           {
01976             //std::cout << "split!" << std::endl;
01977             //if (diag02_len > diag12_len) //split edge 02, introduce line 51
01978             if (stable_line_is_longer(vertices[0], vertices[2],
01979                                       vertices[1], vertices[2])) //split edge 02, introduce line 51
01980             {
01981               //std::cout << "split!" << std::endl;
01982               //std::cout << "apply_4_2_2" << std::endl;
01983               cellvertices[0] = vertices[0];
01984               cellvertices[1] = vertices[1];
01985               cellvertices[2] = vertices[5];
01986               cellvertices[3] = vertices[6];
01987               new_cell.vertices(cellvertices);
01988               segment_out.push_back(new_cell);
01989 
01990               cellvertices[0] = vertices[1];
01991               cellvertices[1] = vertices[4];
01992               cellvertices[2] = vertices[5];
01993               cellvertices[3] = vertices[6];
01994               new_cell.vertices(cellvertices);
01995               segment_out.push_back(new_cell);
01996 
01997               cellvertices[0] = vertices[1];
01998               cellvertices[1] = vertices[4];
01999               cellvertices[2] = vertices[6];
02000               cellvertices[3] = vertices[7];
02001               new_cell.vertices(cellvertices);
02002               segment_out.push_back(new_cell);
02003 
02004               cellvertices[0] = vertices[7];
02005               cellvertices[1] = vertices[4];
02006               cellvertices[2] = vertices[6];
02007               cellvertices[3] = vertices[2];
02008               new_cell.vertices(cellvertices);
02009               segment_out.push_back(new_cell);
02010 
02011               cellvertices[0] = vertices[4];
02012               cellvertices[1] = vertices[2];
02013               cellvertices[2] = vertices[5];
02014               cellvertices[3] = vertices[6];
02015               new_cell.vertices(cellvertices);
02016               segment_out.push_back(new_cell);
02017 
02018               cellvertices[0] = vertices[6];
02019               cellvertices[1] = vertices[7];
02020               cellvertices[2] = vertices[2];
02021               cellvertices[3] = vertices[3];
02022               new_cell.vertices(cellvertices);
02023               segment_out.push_back(new_cell);
02024               //std::cout << "done!" << std::endl;
02025             }
02026             else //split edge 12, introduce line 40
02027             {
02028               //std::cout << "apply_4_2_3" << std::endl;
02029               cellvertices[0] = vertices[0];
02030               cellvertices[1] = vertices[1];
02031               cellvertices[2] = vertices[4];
02032               cellvertices[3] = vertices[6];
02033               new_cell.vertices(cellvertices);
02034               segment_out.push_back(new_cell);
02035 
02036               cellvertices[0] = vertices[0];
02037               cellvertices[1] = vertices[4];
02038               cellvertices[2] = vertices[5];
02039               cellvertices[3] = vertices[6];
02040               new_cell.vertices(cellvertices);
02041               segment_out.push_back(new_cell);
02042 
02043               cellvertices[0] = vertices[6];
02044               cellvertices[1] = vertices[1];
02045               cellvertices[2] = vertices[4];
02046               cellvertices[3] = vertices[7];
02047               new_cell.vertices(cellvertices);
02048               segment_out.push_back(new_cell);
02049 
02050               cellvertices[0] = vertices[6];
02051               cellvertices[1] = vertices[7];
02052               cellvertices[2] = vertices[4];
02053               cellvertices[3] = vertices[2];
02054               new_cell.vertices(cellvertices);
02055               segment_out.push_back(new_cell);
02056 
02057               cellvertices[0] = vertices[6];
02058               cellvertices[1] = vertices[4];
02059               cellvertices[2] = vertices[5];
02060               cellvertices[3] = vertices[2];
02061               new_cell.vertices(cellvertices);
02062               segment_out.push_back(new_cell);
02063 
02064               cellvertices[0] = vertices[6];
02065               cellvertices[1] = vertices[7];
02066               cellvertices[2] = vertices[2];
02067               cellvertices[3] = vertices[3];
02068               new_cell.vertices(cellvertices);
02069               segment_out.push_back(new_cell);
02070             }
02071           }
02072         }
02073         else //split edge 12, introduce line 43
02074         {
02075           //if (diag02_len > diag03_len) //split edge 02, introduce line 53
02076           if (stable_line_is_longer(vertices[0], vertices[2],
02077                                     vertices[0], vertices[3])) //split edge 02, introduce line 53
02078           {
02079             //if (diag02_len > diag12_len) //split edge 02, introduce line 51
02080             if (stable_line_is_longer(vertices[0], vertices[2],
02081                                       vertices[1], vertices[2])) //split edge 02, introduce line 51
02082             {
02083               //std::cout << "apply_4_2_4" << std::endl;
02084               cellvertices[0] = vertices[0];
02085               cellvertices[1] = vertices[1];
02086               cellvertices[2] = vertices[5];
02087               cellvertices[3] = vertices[6];
02088               new_cell.vertices(cellvertices);
02089               segment_out.push_back(new_cell);
02090 
02091               cellvertices[0] = vertices[6];
02092               cellvertices[1] = vertices[1];
02093               cellvertices[2] = vertices[5];
02094               cellvertices[3] = vertices[7];
02095               new_cell.vertices(cellvertices);
02096               segment_out.push_back(new_cell);
02097 
02098               cellvertices[0] = vertices[1];
02099               cellvertices[1] = vertices[4];
02100               cellvertices[2] = vertices[5];
02101               cellvertices[3] = vertices[7];
02102               new_cell.vertices(cellvertices);
02103               segment_out.push_back(new_cell);
02104 
02105               cellvertices[0] = vertices[5];
02106               cellvertices[1] = vertices[4];
02107               cellvertices[2] = vertices[2];
02108               cellvertices[3] = vertices[3];
02109               new_cell.vertices(cellvertices);
02110               segment_out.push_back(new_cell);
02111 
02112               cellvertices[0] = vertices[5];
02113               cellvertices[1] = vertices[7];
02114               cellvertices[2] = vertices[4];
02115               cellvertices[3] = vertices[3];
02116               new_cell.vertices(cellvertices);
02117               segment_out.push_back(new_cell);
02118 
02119               cellvertices[0] = vertices[6];
02120               cellvertices[1] = vertices[7];
02121               cellvertices[2] = vertices[5];
02122               cellvertices[3] = vertices[3];
02123               new_cell.vertices(cellvertices);
02124               segment_out.push_back(new_cell);
02125             }
02126             else //split edge 12, introduce line 40
02127             {
02128               //std::cout << "apply_4_2_5" << std::endl;
02129               cellvertices[0] = vertices[0];
02130               cellvertices[1] = vertices[1];
02131               cellvertices[2] = vertices[4];
02132               cellvertices[3] = vertices[6];
02133               new_cell.vertices(cellvertices);
02134               segment_out.push_back(new_cell);
02135 
02136               cellvertices[0] = vertices[6];
02137               cellvertices[1] = vertices[1];
02138               cellvertices[2] = vertices[4];
02139               cellvertices[3] = vertices[7];
02140               new_cell.vertices(cellvertices);
02141               segment_out.push_back(new_cell);
02142 
02143               cellvertices[0] = vertices[6];
02144               cellvertices[1] = vertices[7];
02145               cellvertices[2] = vertices[4];
02146               cellvertices[3] = vertices[3];
02147               new_cell.vertices(cellvertices);
02148               segment_out.push_back(new_cell);
02149 
02150               cellvertices[0] = vertices[6];
02151               cellvertices[1] = vertices[4];
02152               cellvertices[2] = vertices[5];
02153               cellvertices[3] = vertices[3];
02154               new_cell.vertices(cellvertices);
02155               segment_out.push_back(new_cell);
02156 
02157               cellvertices[0] = vertices[0];
02158               cellvertices[1] = vertices[4];
02159               cellvertices[2] = vertices[5];
02160               cellvertices[3] = vertices[6];
02161               new_cell.vertices(cellvertices);
02162               segment_out.push_back(new_cell);
02163 
02164               cellvertices[0] = vertices[5];
02165               cellvertices[1] = vertices[4];
02166               cellvertices[2] = vertices[2];
02167               cellvertices[3] = vertices[3];
02168               new_cell.vertices(cellvertices);
02169               segment_out.push_back(new_cell);
02170             }
02171           }
02172           else //split edge 03, introduce line 62
02173           {
02174             //if (diag02_len > diag12_len) //split edge 02, introduce line 51
02175             if (stable_line_is_longer(vertices[0], vertices[2],
02176                                       vertices[1], vertices[2])) //split edge 02, introduce line 51
02177             {
02178               //std::cout << "apply_4_2_6" << std::endl;
02179               cellvertices[0] = vertices[0];
02180               cellvertices[1] = vertices[1];
02181               cellvertices[2] = vertices[5];
02182               cellvertices[3] = vertices[6];
02183               new_cell.vertices(cellvertices);
02184               segment_out.push_back(new_cell);
02185 
02186               cellvertices[0] = vertices[1];
02187               cellvertices[1] = vertices[4];
02188               cellvertices[2] = vertices[5];
02189               cellvertices[3] = vertices[6];
02190               new_cell.vertices(cellvertices);
02191               segment_out.push_back(new_cell);
02192 
02193               cellvertices[0] = vertices[1];
02194               cellvertices[1] = vertices[4];
02195               cellvertices[2] = vertices[6];
02196               cellvertices[3] = vertices[7];
02197               new_cell.vertices(cellvertices);
02198               segment_out.push_back(new_cell);
02199 
02200               cellvertices[0] = vertices[7];
02201               cellvertices[1] = vertices[4];
02202               cellvertices[2] = vertices[6];
02203               cellvertices[3] = vertices[3];
02204               new_cell.vertices(cellvertices);
02205               segment_out.push_back(new_cell);
02206 
02207               cellvertices[0] = vertices[4];
02208               cellvertices[1] = vertices[2];
02209               cellvertices[2] = vertices[5];
02210               cellvertices[3] = vertices[6];
02211               new_cell.vertices(cellvertices);
02212               segment_out.push_back(new_cell);
02213 
02214               cellvertices[0] = vertices[6];
02215               cellvertices[1] = vertices[4];
02216               cellvertices[2] = vertices[2];
02217               cellvertices[3] = vertices[3];
02218               new_cell.vertices(cellvertices);
02219               segment_out.push_back(new_cell);
02220             }
02221             else //split edge 12, introduce line 40
02222             {
02223               //std::cout << "apply_4_2_7" << std::endl;
02224               cellvertices[0] = vertices[0];
02225               cellvertices[1] = vertices[1];
02226               cellvertices[2] = vertices[4];
02227               cellvertices[3] = vertices[6];
02228               new_cell.vertices(cellvertices);
02229               segment_out.push_back(new_cell);
02230 
02231               cellvertices[0] = vertices[6];
02232               cellvertices[1] = vertices[1];
02233               cellvertices[2] = vertices[4];
02234               cellvertices[3] = vertices[7];
02235               new_cell.vertices(cellvertices);
02236               segment_out.push_back(new_cell);
02237 
02238               cellvertices[0] = vertices[6];
02239               cellvertices[1] = vertices[7];
02240               cellvertices[2] = vertices[4];
02241               cellvertices[3] = vertices[3];
02242               new_cell.vertices(cellvertices);
02243               segment_out.push_back(new_cell);
02244 
02245               cellvertices[0] = vertices[0];
02246               cellvertices[1] = vertices[4];
02247               cellvertices[2] = vertices[5];
02248               cellvertices[3] = vertices[6];
02249               new_cell.vertices(cellvertices);
02250               segment_out.push_back(new_cell);
02251 
02252               cellvertices[0] = vertices[5];
02253               cellvertices[1] = vertices[4];
02254               cellvertices[2] = vertices[2];
02255               cellvertices[3] = vertices[6];
02256               new_cell.vertices(cellvertices);
02257               segment_out.push_back(new_cell);
02258 
02259               cellvertices[0] = vertices[6];
02260               cellvertices[1] = vertices[4];
02261               cellvertices[2] = vertices[2];
02262               cellvertices[3] = vertices[3];
02263               new_cell.vertices(cellvertices);
02264               segment_out.push_back(new_cell);
02265             }
02266           }
02267         }
02268       }
02269       else //split edge 13, introduce line 70
02270       {
02271         //if (diag13_len > diag12_len) //split edge 13, introduce line 72
02272         if (stable_line_is_longer(vertices[1], vertices[3],
02273                                   vertices[1], vertices[2])) //split edge 13, introduce line 72
02274         {
02275           //if (diag02_len > diag03_len) //split edge 02, introduce line 53
02276           if (stable_line_is_longer(vertices[0], vertices[2],
02277                                     vertices[0], vertices[3])) //split edge 02, introduce line 53
02278           {
02279             //if (diag02_len > diag12_len) //split edge 02, introduce line 51
02280             if (stable_line_is_longer(vertices[0], vertices[2],
02281                                       vertices[1], vertices[2])) //split edge 02, introduce line 51
02282             {
02283               //std::cout << "apply_4_2_8" << std::endl;
02284               cellvertices[0] = vertices[0];
02285               cellvertices[1] = vertices[1];
02286               cellvertices[2] = vertices[5];
02287               cellvertices[3] = vertices[7];
02288               new_cell.vertices(cellvertices);
02289               segment_out.push_back(new_cell);
02290 
02291               cellvertices[0] = vertices[0];
02292               cellvertices[1] = vertices[7];
02293               cellvertices[2] = vertices[5];
02294               cellvertices[3] = vertices[6];
02295               new_cell.vertices(cellvertices);
02296               segment_out.push_back(new_cell);
02297 
02298               cellvertices[0] = vertices[1];
02299               cellvertices[1] = vertices[4];
02300               cellvertices[2] = vertices[5];
02301               cellvertices[3] = vertices[7];
02302               new_cell.vertices(cellvertices);
02303               segment_out.push_back(new_cell);
02304 
02305               cellvertices[0] = vertices[5];
02306               cellvertices[1] = vertices[4];
02307               cellvertices[2] = vertices[2];
02308               cellvertices[3] = vertices[7];
02309               new_cell.vertices(cellvertices);
02310               segment_out.push_back(new_cell);
02311 
02312               cellvertices[0] = vertices[5];
02313               cellvertices[1] = vertices[7];
02314               cellvertices[2] = vertices[2];
02315               cellvertices[3] = vertices[3];
02316               new_cell.vertices(cellvertices);
02317               segment_out.push_back(new_cell);
02318 
02319               cellvertices[0] = vertices[6];
02320               cellvertices[1] = vertices[7];
02321               cellvertices[2] = vertices[5];
02322               cellvertices[3] = vertices[3];
02323               new_cell.vertices(cellvertices);
02324               segment_out.push_back(new_cell);
02325             }
02326             else //split edge 12, introduce line 40
02327             {
02328               //std::cout << "apply_4_2_9" << std::endl;
02329               cellvertices[0] = vertices[0];
02330               cellvertices[1] = vertices[1];
02331               cellvertices[2] = vertices[4];
02332               cellvertices[3] = vertices[7];
02333               new_cell.vertices(cellvertices);
02334               segment_out.push_back(new_cell);
02335 
02336               cellvertices[0] = vertices[0];
02337               cellvertices[1] = vertices[4];
02338               cellvertices[2] = vertices[5];
02339               cellvertices[3] = vertices[7];
02340               new_cell.vertices(cellvertices);
02341               segment_out.push_back(new_cell);
02342 
02343               cellvertices[0] = vertices[0];
02344               cellvertices[1] = vertices[7];
02345               cellvertices[2] = vertices[5];
02346               cellvertices[3] = vertices[6];
02347               new_cell.vertices(cellvertices);
02348               segment_out.push_back(new_cell);
02349 
02350               cellvertices[0] = vertices[6];
02351               cellvertices[1] = vertices[7];
02352               cellvertices[2] = vertices[5];
02353               cellvertices[3] = vertices[3];
02354               new_cell.vertices(cellvertices);
02355               segment_out.push_back(new_cell);
02356 
02357               cellvertices[0] = vertices[5];
02358               cellvertices[1] = vertices[4];
02359               cellvertices[2] = vertices[2];
02360               cellvertices[3] = vertices[7];
02361               new_cell.vertices(cellvertices);
02362               segment_out.push_back(new_cell);
02363 
02364               cellvertices[0] = vertices[5];
02365               cellvertices[1] = vertices[7];
02366               cellvertices[2] = vertices[2];
02367               cellvertices[3] = vertices[3];
02368               new_cell.vertices(cellvertices);
02369               segment_out.push_back(new_cell);
02370             }
02371           }
02372           else //split edge 03, introduce line 62
02373           {
02374             //if (diag02_len > diag12_len) //split edge 02, introduce line 51
02375             if (stable_line_is_longer(vertices[0], vertices[2],
02376                                       vertices[1], vertices[2])) //split edge 02, introduce line 51
02377             {
02378               //std::cout << "apply_4_2_10" << std::endl;
02379               cellvertices[0] = vertices[0];
02380               cellvertices[1] = vertices[1];
02381               cellvertices[2] = vertices[5];
02382               cellvertices[3] = vertices[7];
02383               new_cell.vertices(cellvertices);
02384               segment_out.push_back(new_cell);
02385 
02386               cellvertices[0] = vertices[1];
02387               cellvertices[1] = vertices[4];
02388               cellvertices[2] = vertices[5];
02389               cellvertices[3] = vertices[7];
02390               new_cell.vertices(cellvertices);
02391               segment_out.push_back(new_cell);
02392 
02393               cellvertices[0] = vertices[0];
02394               cellvertices[1] = vertices[7];
02395               cellvertices[2] = vertices[5];
02396               cellvertices[3] = vertices[6];
02397               new_cell.vertices(cellvertices);
02398               segment_out.push_back(new_cell);
02399 
02400               cellvertices[0] = vertices[4];
02401               cellvertices[1] = vertices[2];
02402               cellvertices[2] = vertices[5];
02403               cellvertices[3] = vertices[7];
02404               new_cell.vertices(cellvertices);
02405               segment_out.push_back(new_cell);
02406 
02407               cellvertices[0] = vertices[7];
02408               cellvertices[1] = vertices[2];
02409               cellvertices[2] = vertices[5];
02410               cellvertices[3] = vertices[6];
02411               new_cell.vertices(cellvertices);
02412               segment_out.push_back(new_cell);
02413 
02414               cellvertices[0] = vertices[7];
02415               cellvertices[1] = vertices[2];
02416               cellvertices[2] = vertices[6];
02417               cellvertices[3] = vertices[3];
02418               new_cell.vertices(cellvertices);
02419               segment_out.push_back(new_cell);
02420             }
02421             else //split edge 12, introduce line 40
02422             {
02423               //std::cout << "apply_4_2_11" << std::endl;
02424               cellvertices[0] = vertices[0];
02425               cellvertices[1] = vertices[1];
02426               cellvertices[2] = vertices[4];
02427               cellvertices[3] = vertices[7];
02428               new_cell.vertices(cellvertices);
02429               segment_out.push_back(new_cell);
02430 
02431               cellvertices[0] = vertices[0];
02432               cellvertices[1] = vertices[7];
02433               cellvertices[2] = vertices[4];
02434               cellvertices[3] = vertices[6];
02435               new_cell.vertices(cellvertices);
02436               segment_out.push_back(new_cell);
02437 
02438               cellvertices[0] = vertices[0];
02439               cellvertices[1] = vertices[4];
02440               cellvertices[2] = vertices[5];
02441               cellvertices[3] = vertices[6];
02442               new_cell.vertices(cellvertices);
02443               segment_out.push_back(new_cell);
02444 
02445               cellvertices[0] = vertices[6];
02446               cellvertices[1] = vertices[7];
02447               cellvertices[2] = vertices[4];
02448               cellvertices[3] = vertices[2];
02449               new_cell.vertices(cellvertices);
02450               segment_out.push_back(new_cell);
02451 
02452               cellvertices[0] = vertices[6];
02453               cellvertices[1] = vertices[4];
02454               cellvertices[2] = vertices[5];
02455               cellvertices[3] = vertices[2];
02456               new_cell.vertices(cellvertices);
02457               segment_out.push_back(new_cell);
02458 
02459               cellvertices[0] = vertices[6];
02460               cellvertices[1] = vertices[7];
02461               cellvertices[2] = vertices[2];
02462               cellvertices[3] = vertices[3];
02463               new_cell.vertices(cellvertices);
02464               segment_out.push_back(new_cell);
02465             }
02466           }
02467         }
02468         else //split edge 12, introduce line 43
02469         {
02470           //if (diag02_len > diag03_len) //split edge 02, introduce line 53
02471           if (stable_line_is_longer(vertices[0], vertices[2],
02472                                     vertices[0], vertices[3])) //split edge 02, introduce line 53
02473           {
02474             //if (diag02_len > diag12_len) //split edge 02, introduce line 51
02475             if (stable_line_is_longer(vertices[0], vertices[2],
02476                                       vertices[1], vertices[2])) //split edge 02, introduce line 51
02477             {
02478               //std::cout << "apply_4_2_12" << std::endl;
02479               cellvertices[0] = vertices[0];
02480               cellvertices[1] = vertices[1];
02481               cellvertices[2] = vertices[5];
02482               cellvertices[3] = vertices[7];
02483               new_cell.vertices(cellvertices);
02484               segment_out.push_back(new_cell);
02485 
02486               cellvertices[0] = vertices[0];
02487               cellvertices[1] = vertices[7];
02488               cellvertices[2] = vertices[5];
02489               cellvertices[3] = vertices[6];
02490               new_cell.vertices(cellvertices);
02491               segment_out.push_back(new_cell);
02492 
02493               cellvertices[0] = vertices[6];
02494               cellvertices[1] = vertices[7];
02495               cellvertices[2] = vertices[5];
02496               cellvertices[3] = vertices[3];
02497               new_cell.vertices(cellvertices);
02498               segment_out.push_back(new_cell);
02499 
02500               cellvertices[0] = vertices[1];
02501               cellvertices[1] = vertices[4];
02502               cellvertices[2] = vertices[5];
02503               cellvertices[3] = vertices[7];
02504               new_cell.vertices(cellvertices);
02505               segment_out.push_back(new_cell);
02506 
02507               cellvertices[0] = vertices[7];
02508               cellvertices[1] = vertices[4];
02509               cellvertices[2] = vertices[5];
02510               cellvertices[3] = vertices[3];
02511               new_cell.vertices(cellvertices);
02512               segment_out.push_back(new_cell);
02513 
02514               cellvertices[0] = vertices[4];
02515               cellvertices[1] = vertices[2];
02516               cellvertices[2] = vertices[5];
02517               cellvertices[3] = vertices[3];
02518               new_cell.vertices(cellvertices);
02519               segment_out.push_back(new_cell);
02520             }
02521             else //split edge 12, introduce line 40
02522             {
02523               //std::cout << "apply_4_2_13" << std::endl;
02524               cellvertices[0] = vertices[0];
02525               cellvertices[1] = vertices[1];
02526               cellvertices[2] = vertices[4];
02527               cellvertices[3] = vertices[7];
02528               new_cell.vertices(cellvertices);
02529               segment_out.push_back(new_cell);
02530 
02531               cellvertices[0] = vertices[0];
02532               cellvertices[1] = vertices[4];
02533               cellvertices[2] = vertices[5];
02534               cellvertices[3] = vertices[7];
02535               new_cell.vertices(cellvertices);
02536               segment_out.push_back(new_cell);
02537 
02538               cellvertices[0] = vertices[0];
02539               cellvertices[1] = vertices[7];
02540               cellvertices[2] = vertices[5];
02541               cellvertices[3] = vertices[6];
02542               new_cell.vertices(cellvertices);
02543               segment_out.push_back(new_cell);
02544 
02545               cellvertices[0] = vertices[6];
02546               cellvertices[1] = vertices[7];
02547               cellvertices[2] = vertices[5];
02548               cellvertices[3] = vertices[3];
02549               new_cell.vertices(cellvertices);
02550               segment_out.push_back(new_cell);
02551 
02552               cellvertices[0] = vertices[7];
02553               cellvertices[1] = vertices[4];
02554               cellvertices[2] = vertices[5];
02555               cellvertices[3] = vertices[3];
02556               new_cell.vertices(cellvertices);
02557               segment_out.push_back(new_cell);
02558 
02559               cellvertices[0] = vertices[4];
02560               cellvertices[1] = vertices[2];
02561               cellvertices[2] = vertices[5];
02562               cellvertices[3] = vertices[3];
02563               new_cell.vertices(cellvertices);
02564               segment_out.push_back(new_cell);
02565             }
02566           }
02567           else //split edge 03, introduce line 62
02568           {
02569             //if (diag02_len > diag12_len) //split edge 02, introduce line 51
02570             if (stable_line_is_longer(vertices[0], vertices[2],
02571                                       vertices[1], vertices[2])) //split edge 02, introduce line 51
02572             {
02573               //we have diag12_len > diag13_len > diag03_len > diag02_len alreday, hence this case is bogus!
02574               assert( false && "Logic error: diag02 > diag12 not possible here!");
02575             }
02576             else //split edge 12, introduce line 40
02577             {
02578               //std::cout << "apply_4_2_14" << std::endl;
02579               cellvertices[0] = vertices[0];
02580               cellvertices[1] = vertices[1];
02581               cellvertices[2] = vertices[4];
02582               cellvertices[3] = vertices[7];
02583               new_cell.vertices(cellvertices);
02584               segment_out.push_back(new_cell);
02585 
02586               cellvertices[0] = vertices[0];
02587               cellvertices[1] = vertices[7];
02588               cellvertices[2] = vertices[4];
02589               cellvertices[3] = vertices[6];
02590               new_cell.vertices(cellvertices);
02591               segment_out.push_back(new_cell);
02592 
02593               cellvertices[0] = vertices[0];
02594               cellvertices[1] = vertices[4];
02595               cellvertices[2] = vertices[5];
02596               cellvertices[3] = vertices[6];
02597               new_cell.vertices(cellvertices);
02598               segment_out.push_back(new_cell);
02599 
02600               cellvertices[0] = vertices[6];
02601               cellvertices[1] = vertices[7];
02602               cellvertices[2] = vertices[4];
02603               cellvertices[3] = vertices[3];
02604               new_cell.vertices(cellvertices);
02605               segment_out.push_back(new_cell);
02606 
02607               cellvertices[0] = vertices[6];
02608               cellvertices[1] = vertices[4];
02609               cellvertices[2] = vertices[2];
02610               cellvertices[3] = vertices[3];
02611               new_cell.vertices(cellvertices);
02612               segment_out.push_back(new_cell);
02613 
02614               cellvertices[0] = vertices[5];
02615               cellvertices[1] = vertices[4];
02616               cellvertices[2] = vertices[2];
02617               cellvertices[3] = vertices[6];
02618               new_cell.vertices(cellvertices);
02619               segment_out.push_back(new_cell);
02620             }
02621           }
02622         }
02623       }
02624     }
02625     
02626     
02628     template <typename CellType, typename DomainTypeOut>
02629     static void apply4(CellType const & cell_in, DomainTypeOut & segment_out)
02630     {
02631       typedef typename CellType::config_type        ConfigTypeIn;
02632       typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type            VertexOnCellRange;
02633       typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type         VertexOnCellIterator;            
02634       typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type            EdgeOnCellRange;
02635       typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type           EdgeOnCellIterator;            
02636       
02637       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type             VertexType;
02638       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 1>::type             EdgeType;
02639 
02640       VertexType * vertices[topology::bndcells<tetrahedron_tag, 0>::num];
02641       //std::cout << "apply4()" << std::endl;
02642       
02643       //
02644       // Step 1: Get vertices from input cell
02645       //
02646       VertexOnCellRange vertices_on_cell = viennagrid::ncells<0>(cell_in);
02647       VertexOnCellIterator vocit = vertices_on_cell.begin();
02648       vertices[0] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
02649       vertices[1] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
02650       vertices[2] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
02651       vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]);
02652       
02653 
02654       //
02655       // Step 2: Bring vertices in correct order, such that refined edge is on {0,1}-edge
02656       //
02657       VertexType * ordered_vertices[topology::bndcells<tetrahedron_tag, 0>::num + 4];
02658       EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(cell_in);
02659       EdgeOnCellIterator eocit = edges_on_cell.begin();
02660       EdgeType const & e0 = *eocit; ++eocit;
02661       EdgeType const & e1 = *eocit; ++eocit;
02662       EdgeType const & e2 = *eocit; ++eocit;
02663       EdgeType const & e3 = *eocit; ++eocit;
02664       EdgeType const & e4 = *eocit; ++eocit;
02665       EdgeType const & e5 = *eocit;
02666       
02667       //with e0
02668       if (viennadata::access<refinement_key, bool>(refinement_key())(e0) == false)
02669       {
02670         if (viennadata::access<refinement_key, bool>(refinement_key())(e1) == false)
02671         {
02672           ordered_vertices[0] = vertices[2];
02673           ordered_vertices[1] = vertices[0];
02674           ordered_vertices[2] = vertices[1];
02675           ordered_vertices[3] = vertices[3];
02676           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
02677           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
02678           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
02679           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
02680           
02681           apply4_1(cell_in, segment_out, ordered_vertices);
02682         }
02683         else if (viennadata::access<refinement_key, bool>(refinement_key())(e2) == false)
02684         {
02685           ordered_vertices[0] = vertices[1];
02686           ordered_vertices[1] = vertices[0];
02687           ordered_vertices[2] = vertices[3];
02688           ordered_vertices[3] = vertices[2];
02689           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
02690           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
02691           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
02692           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
02693           
02694           apply4_1(cell_in, segment_out, ordered_vertices);
02695         }
02696         else if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == false) 
02697         {
02698           ordered_vertices[0] = vertices[0];
02699           ordered_vertices[1] = vertices[1];
02700           ordered_vertices[2] = vertices[2];
02701           ordered_vertices[3] = vertices[3];
02702           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
02703           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
02704           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
02705           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
02706           
02707           apply4_1(cell_in, segment_out, ordered_vertices);
02708         }
02709         else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == false) 
02710         {
02711           ordered_vertices[0] = vertices[3];
02712           ordered_vertices[1] = vertices[1];
02713           ordered_vertices[2] = vertices[0];
02714           ordered_vertices[3] = vertices[2];
02715           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
02716           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
02717           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
02718           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
02719           
02720           apply4_1(cell_in, segment_out, ordered_vertices);
02721         }
02722         else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == false) 
02723         {
02724           ordered_vertices[0] = vertices[0];
02725           ordered_vertices[1] = vertices[1];
02726           ordered_vertices[2] = vertices[2];
02727           ordered_vertices[3] = vertices[3];
02728           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
02729           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
02730           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
02731           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
02732           
02733           apply4_2(cell_in, segment_out, ordered_vertices);
02734         }
02735         else
02736         {
02737           assert(false && "Logic error: No edge for refinement found!"); 
02738         }
02739       }
02740       else if (viennadata::access<refinement_key, bool>(refinement_key())(e1) == false)
02741       {
02742         if (viennadata::access<refinement_key, bool>(refinement_key())(e2) == false)
02743         {
02744           ordered_vertices[0] = vertices[3];
02745           ordered_vertices[1] = vertices[0];
02746           ordered_vertices[2] = vertices[2];
02747           ordered_vertices[3] = vertices[1];
02748           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
02749           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
02750           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
02751           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
02752           
02753           apply4_1(cell_in, segment_out, ordered_vertices);
02754         }
02755         else if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == false) 
02756         {
02757           ordered_vertices[0] = vertices[1];
02758           ordered_vertices[1] = vertices[2];
02759           ordered_vertices[2] = vertices[0];
02760           ordered_vertices[3] = vertices[3];
02761           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
02762           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
02763           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
02764           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
02765           
02766           apply4_1(cell_in, segment_out, ordered_vertices);
02767         }
02768         else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == false) 
02769         {
02770           ordered_vertices[0] = vertices[2];
02771           ordered_vertices[1] = vertices[0];
02772           ordered_vertices[2] = vertices[1];
02773           ordered_vertices[3] = vertices[3];
02774           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
02775           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
02776           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
02777           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
02778           
02779           apply4_2(cell_in, segment_out, ordered_vertices);
02780         }
02781         else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == false) 
02782         {
02783           ordered_vertices[0] = vertices[0];
02784           ordered_vertices[1] = vertices[2];
02785           ordered_vertices[2] = vertices[3];
02786           ordered_vertices[3] = vertices[1];
02787           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
02788           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
02789           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
02790           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
02791           
02792           apply4_1(cell_in, segment_out, ordered_vertices);
02793         }
02794         else
02795         {
02796           assert(false && "Logic error: No edge for refinement found!"); 
02797         }
02798       }
02799       else if (viennadata::access<refinement_key, bool>(refinement_key())(e2) == false)
02800       {
02801         if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == false) 
02802         {
02803           ordered_vertices[0] = vertices[3];
02804           ordered_vertices[1] = vertices[0];
02805           ordered_vertices[2] = vertices[2];
02806           ordered_vertices[3] = vertices[1];
02807           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
02808           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
02809           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
02810           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
02811           
02812           apply4_2(cell_in, segment_out, ordered_vertices);
02813         }
02814         else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == false) 
02815         {
02816           ordered_vertices[0] = vertices[0];
02817           ordered_vertices[1] = vertices[3];
02818           ordered_vertices[2] = vertices[1];
02819           ordered_vertices[3] = vertices[2];
02820           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
02821           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
02822           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
02823           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
02824           
02825           apply4_1(cell_in, segment_out, ordered_vertices);
02826         }
02827         else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == false) 
02828         {
02829           ordered_vertices[0] = vertices[2];
02830           ordered_vertices[1] = vertices[3];
02831           ordered_vertices[2] = vertices[0];
02832           ordered_vertices[3] = vertices[1];
02833           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
02834           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
02835           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
02836           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
02837           
02838           apply4_1(cell_in, segment_out, ordered_vertices);
02839         }
02840         else
02841         {
02842           assert(false && "Logic error: No edge for refinement found!"); 
02843         }
02844       }
02845       else if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == false)
02846       {
02847         if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == false) 
02848         {
02849           ordered_vertices[0] = vertices[2];
02850           ordered_vertices[1] = vertices[1];
02851           ordered_vertices[2] = vertices[3];
02852           ordered_vertices[3] = vertices[0];
02853           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
02854           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
02855           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
02856           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
02857           
02858           apply4_1(cell_in, segment_out, ordered_vertices);
02859         }
02860         else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == false) 
02861         {
02862           ordered_vertices[0] = vertices[3];
02863           ordered_vertices[1] = vertices[2];
02864           ordered_vertices[2] = vertices[1];
02865           ordered_vertices[3] = vertices[0];
02866           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
02867           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
02868           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
02869           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
02870           
02871           apply4_1(cell_in, segment_out, ordered_vertices);
02872         }
02873         else
02874         {
02875           assert(false && "Logic error: No edge for refinement found!"); 
02876         }
02877       }
02878       else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == false)
02879       {
02880         if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == false) 
02881         {
02882           ordered_vertices[0] = vertices[1];
02883           ordered_vertices[1] = vertices[3];
02884           ordered_vertices[2] = vertices[2];
02885           ordered_vertices[3] = vertices[0];
02886           ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
02887           ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
02888           ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
02889           ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
02890           
02891           apply4_1(cell_in, segment_out, ordered_vertices);
02892         }
02893         else
02894         {
02895           assert(false && "Logic error: No edge for refinement found!"); 
02896         }
02897       }
02898       else
02899       {
02900         assert(false && "Logic error: No edge for refinement found!"); 
02901       }
02902     }    
02903 
02904 
02905 
02907 
02908     
02922     template <typename CellType, typename DomainTypeOut, typename VertexType>
02923     static void apply5_1(CellType const & cell_in, 
02924                          DomainTypeOut & segment_out,
02925                          VertexType ** vertices
02926                         )
02927     {
02928       CellType new_cell;
02929       VertexType * cellvertices[topology::bndcells<tetrahedron_tag, 0>::num];
02930 
02931       cellvertices[0] = vertices[6];
02932       cellvertices[1] = vertices[7];
02933       cellvertices[2] = vertices[8];
02934       cellvertices[3] = vertices[3];
02935       new_cell.vertices(cellvertices);
02936       segment_out.push_back(new_cell);
02937 
02938       cellvertices[0] = vertices[5];
02939       cellvertices[1] = vertices[4];
02940       cellvertices[2] = vertices[2];
02941       cellvertices[3] = vertices[8];
02942       new_cell.vertices(cellvertices);
02943       segment_out.push_back(new_cell);
02944       
02945 
02946       //if (diag03_len > diag13_len) //split edge 03, introduce line 61
02947       if (stable_line_is_longer(vertices[0], vertices[3],
02948                                 vertices[1], vertices[3])) //split edge 03, introduce line 61
02949       {
02950         cellvertices[0] = vertices[6];
02951         cellvertices[1] = vertices[4];
02952         cellvertices[2] = vertices[5];
02953         cellvertices[3] = vertices[8];
02954         new_cell.vertices(cellvertices);
02955         segment_out.push_back(new_cell);
02956         
02957         cellvertices[0] = vertices[6];
02958         cellvertices[1] = vertices[7];
02959         cellvertices[2] = vertices[4];
02960         cellvertices[3] = vertices[8];
02961         new_cell.vertices(cellvertices);
02962         segment_out.push_back(new_cell);
02963         
02964         cellvertices[0] = vertices[1];
02965         cellvertices[1] = vertices[4];
02966         cellvertices[2] = vertices[6];
02967         cellvertices[3] = vertices[7];
02968         new_cell.vertices(cellvertices);
02969         segment_out.push_back(new_cell);
02970         
02971         //if (diag02_len > diag12_len) //split edge 02, introduce line 51
02972         if (stable_line_is_longer(vertices[0], vertices[2],
02973                                   vertices[1], vertices[2])) //split edge 02, introduce line 51
02974         {
02975           cellvertices[0] = vertices[0];
02976           cellvertices[1] = vertices[1];
02977           cellvertices[2] = vertices[5];
02978           cellvertices[3] = vertices[6];
02979           new_cell.vertices(cellvertices);
02980           segment_out.push_back(new_cell);
02981           
02982           cellvertices[0] = vertices[1];
02983           cellvertices[1] = vertices[4];
02984           cellvertices[2] = vertices[5];
02985           cellvertices[3] = vertices[6];
02986           new_cell.vertices(cellvertices);
02987           segment_out.push_back(new_cell);
02988         }
02989         else //split edge 12, introduce line 40
02990         {
02991           cellvertices[0] = vertices[0];
02992           cellvertices[1] = vertices[1];
02993           cellvertices[2] = vertices[4];
02994           cellvertices[3] = vertices[6];
02995           new_cell.vertices(cellvertices);
02996           segment_out.push_back(new_cell);
02997           
02998           cellvertices[0] = vertices[0];
02999           cellvertices[1] = vertices[4];
03000           cellvertices[2] = vertices[5];
03001           cellvertices[3] = vertices[6];
03002           new_cell.vertices(cellvertices);
03003           segment_out.push_back(new_cell);
03004         }
03005       }
03006       else  //split edge 13, introduce line 70
03007       {
03008         //if (diag02_len > diag12_len) //split edge 02, introduce line 51
03009         if (stable_line_is_longer(vertices[0], vertices[2],
03010                                   vertices[1], vertices[2])) //split edge 02, introduce line 51
03011         {
03012           cellvertices[0] = vertices[0];
03013           cellvertices[1] = vertices[1];
03014           cellvertices[2] = vertices[5];
03015           cellvertices[3] = vertices[7];
03016           new_cell.vertices(cellvertices);
03017           segment_out.push_back(new_cell);
03018           
03019           cellvertices[0] = vertices[0];
03020           cellvertices[1] = vertices[7];
03021           cellvertices[2] = vertices[5];
03022           cellvertices[3] = vertices[6];
03023           new_cell.vertices(cellvertices);
03024           segment_out.push_back(new_cell);
03025           
03026           cellvertices[0] = vertices[1];
03027           cellvertices[1] = vertices[4];
03028           cellvertices[2] = vertices[5];
03029           cellvertices[3] = vertices[7];
03030           new_cell.vertices(cellvertices);
03031           segment_out.push_back(new_cell);
03032           
03033           cellvertices[0] = vertices[7];
03034           cellvertices[1] = vertices[4];
03035           cellvertices[2] = vertices[5];
03036           cellvertices[3] = vertices[8];
03037           new_cell.vertices(cellvertices);
03038           segment_out.push_back(new_cell);
03039           
03040           cellvertices[0] = vertices[6];
03041           cellvertices[1] = vertices[7];
03042           cellvertices[2] = vertices[5];
03043           cellvertices[3] = vertices[8];
03044           new_cell.vertices(cellvertices);
03045           segment_out.push_back(new_cell);
03046         }
03047         else //split edge 12, introduce line 40
03048         {
03049           cellvertices[0] = vertices[0];
03050           cellvertices[1] = vertices[1];
03051           cellvertices[2] = vertices[4];
03052           cellvertices[3] = vertices[7];
03053           new_cell.vertices(cellvertices);
03054           segment_out.push_back(new_cell);
03055           
03056           cellvertices[0] = vertices[0];
03057           cellvertices[1] = vertices[7];
03058           cellvertices[2] = vertices[4];
03059           cellvertices[3] = vertices[6];
03060           new_cell.vertices(cellvertices);
03061           segment_out.push_back(new_cell);
03062           
03063           cellvertices[0] = vertices[0];
03064           cellvertices[1] = vertices[4];
03065           cellvertices[2] = vertices[5];
03066           cellvertices[3] = vertices[6];
03067           new_cell.vertices(cellvertices);
03068           segment_out.push_back(new_cell);
03069           
03070           cellvertices[0] = vertices[6];
03071           cellvertices[1] = vertices[4];
03072           cellvertices[2] = vertices[5];
03073           cellvertices[3] = vertices[8];
03074           new_cell.vertices(cellvertices);
03075           segment_out.push_back(new_cell);
03076           
03077           cellvertices[0] = vertices[6];
03078           cellvertices[1] = vertices[7];
03079           cellvertices[2] = vertices[4];
03080           cellvertices[3] = vertices[8];
03081           new_cell.vertices(cellvertices);
03082           segment_out.push_back(new_cell);
03083         }
03084       }
03085     }
03086     
03088     template <typename CellType, typename DomainTypeOut>
03089     static void apply5(CellType const & cell_in, DomainTypeOut & segment_out)
03090     {
03091       typedef typename CellType::config_type        ConfigTypeIn;
03092       typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type            VertexOnCellRange;
03093       typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type         VertexOnCellIterator;            
03094       typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type            EdgeOnCellRange;
03095       typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type           EdgeOnCellIterator;            
03096       
03097       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type             VertexType;
03098       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 1>::type             EdgeType;
03099 
03100       VertexType * vertices[topology::bndcells<tetrahedron_tag, 0>::num];
03101       //std::cout << "apply5()" << std::endl;
03102       
03103       //
03104       // Step 1: Get vertices from input cell
03105       //
03106       VertexOnCellRange vertices_on_cell = viennagrid::ncells<0>(cell_in);
03107       VertexOnCellIterator vocit = vertices_on_cell.begin();
03108       vertices[0] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
03109       vertices[1] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
03110       vertices[2] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
03111       vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]);
03112       
03113 
03114       //
03115       // Step 2: Bring vertices in correct order, such that refined edge is on {0,1}-edge
03116       //
03117       VertexType * ordered_vertices[topology::bndcells<tetrahedron_tag, 0>::num + 5];
03118       EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(cell_in);
03119       EdgeOnCellIterator eocit = edges_on_cell.begin();
03120       EdgeType const & e0 = *eocit; ++eocit;
03121       EdgeType const & e1 = *eocit; ++eocit;
03122       EdgeType const & e2 = *eocit; ++eocit;
03123       EdgeType const & e3 = *eocit; ++eocit;
03124       EdgeType const & e4 = *eocit; ++eocit;
03125       EdgeType const & e5 = *eocit;
03126       
03127       if (viennadata::access<refinement_key, bool>(refinement_key())(e0) == false)
03128       {
03129         ordered_vertices[0] = vertices[0];
03130         ordered_vertices[1] = vertices[1];
03131         ordered_vertices[2] = vertices[2];
03132         ordered_vertices[3] = vertices[3];
03133         ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
03134         ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
03135         ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
03136         ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
03137         ordered_vertices[8] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
03138         
03139         apply5_1(cell_in, segment_out, ordered_vertices);
03140       }
03141       else if (viennadata::access<refinement_key, bool>(refinement_key())(e1) == false)
03142       {
03143         ordered_vertices[0] = vertices[2];
03144         ordered_vertices[1] = vertices[0];
03145         ordered_vertices[2] = vertices[1];
03146         ordered_vertices[3] = vertices[3];
03147         ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
03148         ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
03149         ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
03150         ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
03151         ordered_vertices[8] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
03152         
03153         apply5_1(cell_in, segment_out, ordered_vertices);
03154       }
03155       else if (viennadata::access<refinement_key, bool>(refinement_key())(e2) == false)
03156       {
03157         ordered_vertices[0] = vertices[0];
03158         ordered_vertices[1] = vertices[3];
03159         ordered_vertices[2] = vertices[1];
03160         ordered_vertices[3] = vertices[2];
03161         ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
03162         ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
03163         ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
03164         ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
03165         ordered_vertices[8] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
03166         
03167         apply5_1(cell_in, segment_out, ordered_vertices);
03168       }
03169       else if (viennadata::access<refinement_key, bool>(refinement_key())(e3) == false)
03170       {
03171         ordered_vertices[0] = vertices[1];
03172         ordered_vertices[1] = vertices[2];
03173         ordered_vertices[2] = vertices[0];
03174         ordered_vertices[3] = vertices[3];
03175         ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
03176         ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
03177         ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
03178         ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
03179         ordered_vertices[8] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
03180         
03181         apply5_1(cell_in, segment_out, ordered_vertices);
03182       }
03183       else if (viennadata::access<refinement_key, bool>(refinement_key())(e4) == false)
03184       {
03185         ordered_vertices[0] = vertices[1];
03186         ordered_vertices[1] = vertices[3];
03187         ordered_vertices[2] = vertices[2];
03188         ordered_vertices[3] = vertices[0];
03189         ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e5)]);
03190         ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
03191         ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
03192         ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
03193         ordered_vertices[8] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
03194         
03195         apply5_1(cell_in, segment_out, ordered_vertices);
03196       }
03197       else if (viennadata::access<refinement_key, bool>(refinement_key())(e5) == false)
03198       {
03199         ordered_vertices[0] = vertices[3];
03200         ordered_vertices[1] = vertices[2];
03201         ordered_vertices[2] = vertices[1];
03202         ordered_vertices[3] = vertices[0];
03203         ordered_vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e3)]);
03204         ordered_vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e4)]);
03205         ordered_vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
03206         ordered_vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
03207         ordered_vertices[8] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
03208         
03209         apply5_1(cell_in, segment_out, ordered_vertices);
03210       }
03211       else
03212       {
03213         assert(false && "Logic error: No edge for refinement found!"); 
03214       }
03215     }    
03216 
03217 
03218 
03219 
03221 
03222 
03223 
03225     template <typename CellType, typename DomainTypeOut>
03226     static void apply6(CellType const & cell_in, DomainTypeOut & segment_out)
03227     {
03228       typedef typename CellType::config_type        ConfigTypeIn;
03229       typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type            VertexOnCellRange;
03230       typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type         VertexOnCellIterator;            
03231       typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type            EdgeOnCellRange;
03232       typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type           EdgeOnCellIterator;            
03233       
03234       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type             VertexType;
03235 
03236       VertexType * vertices[topology::bndcells<tetrahedron_tag, 0>::num
03237                             + topology::bndcells<tetrahedron_tag, 1>::num];
03238       //std::cout << "apply6()" << std::endl;
03239       
03240       //
03241       // Step 1: Get vertices on the new domain
03242       //
03243       
03244       //grab existing vertices:
03245       VertexOnCellRange vertices_on_cell = viennagrid::ncells<0>(cell_in);
03246       VertexOnCellIterator vocit = vertices_on_cell.begin();
03247       vertices[0] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
03248       vertices[1] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
03249       vertices[2] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
03250       vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]);
03251 
03252       //add vertices from edge
03253       EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(cell_in);
03254       EdgeOnCellIterator eocit = edges_on_cell.begin();
03255       vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(*eocit)]); ++eocit;
03256       vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(*eocit)]); ++eocit;
03257       vertices[6] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(*eocit)]); ++eocit;
03258       vertices[7] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(*eocit)]); ++eocit;
03259       vertices[8] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(*eocit)]); ++eocit;
03260       vertices[9] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(*eocit)]);
03261       
03262       //
03263       // Step 2: Add new cells to new domain:
03264       //
03265       CellType new_cell;
03266       VertexType * cellvertices[topology::bndcells<tetrahedron_tag, 0>::num];
03267       
03268       //0-4-5-6:
03269       cellvertices[0] = vertices[0];
03270       cellvertices[1] = vertices[4];
03271       cellvertices[2] = vertices[5];
03272       cellvertices[3] = vertices[6];
03273       new_cell.vertices(cellvertices);
03274       segment_out.push_back(new_cell);
03275 
03276       //1-7-4-8:
03277       cellvertices[0] = vertices[1];
03278       cellvertices[1] = vertices[7];
03279       cellvertices[2] = vertices[4];
03280       cellvertices[3] = vertices[8];
03281       new_cell.vertices(cellvertices);
03282       segment_out.push_back(new_cell);
03283 
03284       //2-5-7-9:
03285       cellvertices[0] = vertices[2];
03286       cellvertices[1] = vertices[5];
03287       cellvertices[2] = vertices[7];
03288       cellvertices[3] = vertices[9];
03289       new_cell.vertices(cellvertices);
03290       segment_out.push_back(new_cell);
03291 
03292       //3-8-6-9:
03293       cellvertices[0] = vertices[3];
03294       cellvertices[1] = vertices[8];
03295       cellvertices[2] = vertices[6];
03296       cellvertices[3] = vertices[9];
03297       new_cell.vertices(cellvertices);
03298       segment_out.push_back(new_cell);
03299       
03300       double diag58 = viennagrid::norm(vertices[5]->point() - vertices[8]->point());
03301       double diag67 = viennagrid::norm(vertices[6]->point() - vertices[7]->point());
03302       double diag49 = viennagrid::norm(vertices[4]->point() - vertices[9]->point());
03303       
03304       if ( (diag58 <= diag67) && (diag58 <= diag49) )  //diag58 is shortest: keep it, split others
03305       {
03306         //4-8-5-6:
03307         cellvertices[0] = vertices[4];
03308         cellvertices[1] = vertices[8];
03309         cellvertices[2] = vertices[5];
03310         cellvertices[3] = vertices[6];
03311         new_cell.vertices(cellvertices);
03312         segment_out.push_back(new_cell);
03313         
03314         //4-8-7-5:
03315         cellvertices[0] = vertices[4];
03316         cellvertices[1] = vertices[8];
03317         cellvertices[2] = vertices[7];
03318         cellvertices[3] = vertices[5];
03319         new_cell.vertices(cellvertices);
03320         segment_out.push_back(new_cell);
03321 
03322         //7-8-9-5:
03323         cellvertices[0] = vertices[7];
03324         cellvertices[1] = vertices[8];
03325         cellvertices[2] = vertices[9];
03326         cellvertices[3] = vertices[5];
03327         new_cell.vertices(cellvertices);
03328         segment_out.push_back(new_cell);
03329 
03330         //8-6-9-5:
03331         cellvertices[0] = vertices[8];
03332         cellvertices[1] = vertices[6];
03333         cellvertices[2] = vertices[9];
03334         cellvertices[3] = vertices[5];
03335         new_cell.vertices(cellvertices);
03336         segment_out.push_back(new_cell);
03337       }
03338       else if ( (diag67 <= diag58) && (diag67 <= diag49) ) //diag67 is shortest: keep it, split others
03339       {
03340         //4-7-6-8:
03341         cellvertices[0] = vertices[4];
03342         cellvertices[1] = vertices[7];
03343         cellvertices[2] = vertices[6];
03344         cellvertices[3] = vertices[8];
03345         new_cell.vertices(cellvertices);
03346         segment_out.push_back(new_cell);
03347         
03348         //4-7-5-6:
03349         cellvertices[0] = vertices[4];
03350         cellvertices[1] = vertices[7];
03351         cellvertices[2] = vertices[5];
03352         cellvertices[3] = vertices[6];
03353         new_cell.vertices(cellvertices);
03354         segment_out.push_back(new_cell);
03355 
03356         //7-9-6-8:
03357         cellvertices[0] = vertices[7];
03358         cellvertices[1] = vertices[9];
03359         cellvertices[2] = vertices[6];
03360         cellvertices[3] = vertices[8];
03361         new_cell.vertices(cellvertices);
03362         segment_out.push_back(new_cell);
03363 
03364         //7-9-5-6:
03365         cellvertices[0] = vertices[7];
03366         cellvertices[1] = vertices[9];
03367         cellvertices[2] = vertices[5];
03368         cellvertices[3] = vertices[6];
03369         new_cell.vertices(cellvertices);
03370         segment_out.push_back(new_cell);
03371       }
03372       else //keep shortest diagonal diag49
03373       {
03374         //4-9-6-8:
03375         cellvertices[0] = vertices[4];
03376         cellvertices[1] = vertices[9];
03377         cellvertices[2] = vertices[6];
03378         cellvertices[3] = vertices[8];
03379         new_cell.vertices(cellvertices);
03380         segment_out.push_back(new_cell);
03381         
03382         //4-9-5-6:
03383         cellvertices[0] = vertices[4];
03384         cellvertices[1] = vertices[9];
03385         cellvertices[2] = vertices[5];
03386         cellvertices[3] = vertices[6];
03387         new_cell.vertices(cellvertices);
03388         segment_out.push_back(new_cell);
03389 
03390         //4-7-9-8:
03391         cellvertices[0] = vertices[4];
03392         cellvertices[1] = vertices[7];
03393         cellvertices[2] = vertices[9];
03394         cellvertices[3] = vertices[8];
03395         new_cell.vertices(cellvertices);
03396         segment_out.push_back(new_cell);
03397 
03398         //4-7-5-9:
03399         cellvertices[0] = vertices[4];
03400         cellvertices[1] = vertices[7];
03401         cellvertices[2] = vertices[5];
03402         cellvertices[3] = vertices[9];
03403         new_cell.vertices(cellvertices);
03404         segment_out.push_back(new_cell);
03405       }
03406       
03407     } //apply6()
03408     
03409     
03410     //
03416     template <typename CellType, typename DomainTypeOut>
03417     static void apply(CellType const & cell_in, DomainTypeOut & segment_out)
03418     {
03419       typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type            EdgeOnCellRange;
03420       typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type                 EdgeOnCellIterator;            
03421       
03422       //std::cout << "tetrahedron_tag::apply()" << std::endl;
03423       std::size_t edges_to_refine = 0;
03424       EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(cell_in);
03425       for (EdgeOnCellIterator eocit = edges_on_cell.begin();
03426                               eocit != edges_on_cell.end();
03427                             ++eocit)
03428       {
03429         if (viennadata::access<refinement_key, bool>(refinement_key())(*eocit) == true)
03430           ++edges_to_refine;
03431       }
03432       
03433       switch (edges_to_refine)
03434       {
03435         case 0: apply0(cell_in, segment_out); break;
03436         case 1: apply1(cell_in, segment_out); break;
03437         case 2: apply2(cell_in, segment_out); break;
03438         case 3: apply3(cell_in, segment_out); break;
03439         case 4: apply4(cell_in, segment_out); break;
03440         case 5: apply5(cell_in, segment_out); break;
03441         case 6: apply6(cell_in, segment_out); break;
03442         default: //nothing to do...
03443                 break;
03444       }
03445     } //apply()
03446     
03447   };
03448   
03449     
03450   
03451 }
03452 
03453 #endif
03454 

Generated on Wed Sep 14 2011 19:21:30 for ViennaGrid by  doxygen 1.7.1