00001 #ifndef VIENNAGRID_TOPOLOGY_TRIANGLE_HPP
00002 #define VIENNAGRID_TOPOLOGY_TRIANGLE_HPP
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "viennagrid/forwards.h"
00022 #include "viennagrid/topology/point.hpp"
00023 #include "viennagrid/topology/line.hpp"
00024 #include "viennagrid/detail/element_iterators.hpp"
00025 #include "viennagrid/algorithm/norm.hpp"
00026
00031 namespace viennagrid
00032 {
00033
00035 template <>
00036 struct simplex_tag<2>
00037 {
00038 enum{ dim = 2 };
00039 static std::string name() { return "Triangle"; }
00040 };
00041
00042
00043 namespace topology
00044 {
00046 template <>
00047 struct bndcells<triangle_tag, 0>
00048 {
00049 typedef point_tag tag;
00050
00051 enum{ num = 3 };
00052 };
00053
00055 template <>
00056 struct bndcells<triangle_tag, 1>
00057 {
00058 typedef simplex_tag<1> tag;
00059
00060 enum{ num = 3 };
00061
00062 };
00063
00065
00066
00068 template <>
00069 struct bndcell_filler<triangle_tag, 1>
00070 {
00071
00072 template <typename ElementType, typename Vertices, typename Orientations, typename Segment>
00073 static void fill(ElementType ** elements, Vertices ** vertices, Orientations * orientations, Segment & seg)
00074 {
00075 Vertices * edgevertices[2];
00076 ElementType edge;
00077
00078 edgevertices[0] = vertices[0];
00079 edgevertices[1] = vertices[1];
00080 edge.vertices(edgevertices);
00081 elements[0] = seg.push_back(edge, (orientations == NULL) ? NULL : orientations);
00082
00083 edgevertices[0] = vertices[0];
00084 edgevertices[1] = vertices[2];
00085 edge.vertices(edgevertices);
00086 elements[1] = seg.push_back(edge, (orientations == NULL) ? NULL : orientations + 1 );
00087
00088 edgevertices[0] = vertices[1];
00089 edgevertices[1] = vertices[2];
00090 edge.vertices(edgevertices);
00091 elements[2] = seg.push_back(edge, (orientations == NULL) ? NULL : orientations + 2 );
00092 }
00093 };
00094
00095 }
00096
00097
00098
00100 template <>
00101 struct element_refinement<triangle_tag>
00102 {
00103
00105 template <typename CellType, typename DomainTypeOut>
00106 static void apply0(CellType const & cell_in, DomainTypeOut & segment_out)
00107 {
00108 typedef typename CellType::config_type ConfigTypeIn;
00109 typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type VertexOnCellRange;
00110 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator;
00111 typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type EdgeOnCellRange;
00112 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator;
00113
00114 typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type VertexType;
00115
00116 VertexType * vertices[topology::bndcells<triangle_tag, 0>::num];
00117
00118
00119
00120
00121
00122
00123 VertexOnCellRange vertices_on_cell = viennagrid::ncells<0>(cell_in);
00124 VertexOnCellIterator vocit = vertices_on_cell.begin();
00125 vertices[0] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00126 vertices[1] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00127 vertices[2] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]);
00128
00129
00130
00131
00132 CellType new_cell;
00133 VertexType * cellvertices[topology::bndcells<triangle_tag, 0>::num];
00134
00135
00136 cellvertices[0] = vertices[0];
00137 cellvertices[1] = vertices[1];
00138 cellvertices[2] = vertices[2];
00139 new_cell.vertices(cellvertices);
00140 segment_out.push_back(new_cell);
00141
00142 }
00143
00144
00154 template <typename CellType, typename DomainTypeOut>
00155 static void apply1(CellType const & cell_in, DomainTypeOut & segment_out)
00156 {
00157 typedef typename CellType::config_type ConfigTypeIn;
00158 typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type VertexOnCellRange;
00159 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator;
00160 typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type EdgeOnCellRange;
00161 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator;
00162
00163 typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type VertexType;
00164 typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 1>::type EdgeType;
00165
00166 VertexType * vertices[topology::bndcells<triangle_tag, 0>::num + 1];
00167
00168
00169
00170
00171
00172
00173 VertexOnCellRange vertices_on_cell = viennagrid::ncells<0>(cell_in);
00174 VertexOnCellIterator vocit = vertices_on_cell.begin();
00175 vertices[0] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00176 vertices[1] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00177 vertices[2] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]);
00178
00179
00180 EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(cell_in);
00181 std::size_t offset = 0;
00182 EdgeOnCellIterator eocit = edges_on_cell.begin();
00183 EdgeType const & e0 = *eocit; ++eocit;
00184 EdgeType const & e1 = *eocit; ++eocit;
00185 EdgeType const & e2 = *eocit;
00186
00187 if ( (viennadata::access<refinement_key, bool>(refinement_key())(e0) == true) )
00188 {
00189 vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
00190 offset = 0;
00191 }
00192 else if ( (viennadata::access<refinement_key, bool>(refinement_key())(e1) == true) )
00193 {
00194 vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
00195 offset = 2;
00196 }
00197 else if ( (viennadata::access<refinement_key, bool>(refinement_key())(e2) == true) )
00198 {
00199 vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
00200 offset = 1;
00201 }
00202 else
00203 {
00204 assert( (2 < 2) && "Logic error: Triangle does not have an edges for bisection!");
00205 }
00206
00207
00208
00209
00210 CellType new_cell;
00211 VertexType * cellvertices[topology::bndcells<triangle_tag, 0>::num];
00212
00213
00214 cellvertices[0] = vertices[(offset + 0) % topology::bndcells<triangle_tag, 0>::num];
00215 cellvertices[1] = vertices[3];
00216 cellvertices[2] = vertices[(offset + 2) % topology::bndcells<triangle_tag, 0>::num];
00217 new_cell.vertices(cellvertices);
00218 segment_out.push_back(new_cell);
00219
00220
00221 cellvertices[0] = vertices[3];
00222 cellvertices[1] = vertices[(offset + 1) % topology::bndcells<triangle_tag, 0>::num];
00223 cellvertices[2] = vertices[(offset + 2) % topology::bndcells<triangle_tag, 0>::num];
00224 new_cell.vertices(cellvertices);
00225 segment_out.push_back(new_cell);
00226 }
00227
00228
00239 template <typename CellType, typename DomainTypeOut>
00240 static void apply2(CellType const & cell_in, DomainTypeOut & segment_out)
00241 {
00242 typedef typename CellType::config_type ConfigTypeIn;
00243 typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type VertexOnCellRange;
00244 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator;
00245 typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type EdgeOnCellRange;
00246 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator;
00247
00248 typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type VertexType;
00249 typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 1>::type EdgeType;
00250
00251 VertexType * vertices[topology::bndcells<triangle_tag, 0>::num + 2];
00252
00253
00254
00255
00256
00257
00258 VertexOnCellRange vertices_on_cell = viennagrid::ncells<0>(cell_in);
00259 VertexOnCellIterator vocit = vertices_on_cell.begin();
00260 vertices[0] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00261 vertices[1] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00262 vertices[2] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]);
00263
00264
00265 EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(cell_in);
00266 std::size_t offset = 0;
00267
00268 EdgeOnCellIterator eocit = edges_on_cell.begin();
00269 EdgeType const & e0 = *eocit; ++eocit;
00270 EdgeType const & e1 = *eocit; ++eocit;
00271 EdgeType const & e2 = *eocit;
00272
00273 if ( (viennadata::access<refinement_key, bool>(refinement_key())(e0) == true)
00274 && (viennadata::access<refinement_key, bool>(refinement_key())(e1) == true) )
00275 {
00276 vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
00277 vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
00278 offset = 2;
00279 }
00280 else if ( (viennadata::access<refinement_key, bool>(refinement_key())(e0) == true)
00281 && (viennadata::access<refinement_key, bool>(refinement_key())(e2) == true) )
00282 {
00283 vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e0)]);
00284 vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
00285 offset = 0;
00286 }
00287 else if ( (viennadata::access<refinement_key, bool>(refinement_key())(e1) == true)
00288 && (viennadata::access<refinement_key, bool>(refinement_key())(e2) == true) )
00289 {
00290 vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e2)]);
00291 vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(e1)]);
00292 offset = 1;
00293 }
00294 else
00295 {
00296 assert( (2 < 2) && "Logic error: Triangle does not have two edges for bisection!");
00297 }
00298
00299
00300
00301
00302 CellType new_cell;
00303 VertexType * cellvertices[topology::bndcells<triangle_tag, 0>::num];
00304
00305
00306 cellvertices[0] = vertices[3];
00307 cellvertices[1] = vertices[(offset + 1) % topology::bndcells<triangle_tag, 0>::num];
00308 cellvertices[2] = vertices[4];
00309 new_cell.vertices(cellvertices);
00310 segment_out.push_back(new_cell);
00311
00312
00313 VertexType & v0 = *(vertices[(offset + 0) % topology::bndcells<triangle_tag, 0>::num]);
00314 VertexType & v1 = *(vertices[(offset + 1) % topology::bndcells<triangle_tag, 0>::num]);
00315 VertexType & v2 = *(vertices[(offset + 2) % topology::bndcells<triangle_tag, 0>::num]);
00316 double len_edge1 = viennagrid::norm(v1.point() - v0.point());
00317 double len_edge2 = viennagrid::norm(v2.point() - v1.point());
00318
00319 if (len_edge1 > len_edge2)
00320 {
00321
00322 cellvertices[0] = vertices[(offset + 0) % topology::bndcells<triangle_tag, 0>::num];
00323 cellvertices[1] = vertices[3];
00324 cellvertices[2] = vertices[(offset + 2) % topology::bndcells<triangle_tag, 0>::num];
00325 new_cell.vertices(cellvertices);
00326 segment_out.push_back(new_cell);
00327
00328
00329 cellvertices[0] = vertices[(offset + 2) % topology::bndcells<triangle_tag, 0>::num];
00330 cellvertices[1] = vertices[3];
00331 cellvertices[2] = vertices[4];
00332 new_cell.vertices(cellvertices);
00333 segment_out.push_back(new_cell);
00334 }
00335 else
00336 {
00337
00338 cellvertices[0] = vertices[(offset + 0) % topology::bndcells<triangle_tag, 0>::num];
00339 cellvertices[1] = vertices[3];
00340 cellvertices[2] = vertices[4];
00341 new_cell.vertices(cellvertices);
00342 segment_out.push_back(new_cell);
00343
00344
00345 cellvertices[0] = vertices[(offset + 0) % topology::bndcells<triangle_tag, 0>::num];
00346 cellvertices[1] = vertices[4];
00347 cellvertices[2] = vertices[(offset + 2) % topology::bndcells<triangle_tag, 0>::num];
00348 new_cell.vertices(cellvertices);
00349 segment_out.push_back(new_cell);
00350 }
00351
00352
00353 }
00354
00355
00356
00357
00359 template <typename CellType, typename DomainTypeOut>
00360 static void apply3(CellType const & cell_in, DomainTypeOut & segment_out)
00361 {
00362 typedef typename CellType::config_type ConfigTypeIn;
00363 typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type VertexOnCellRange;
00364 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator;
00365 typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type EdgeOnCellRange;
00366 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator;
00367
00368 typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type VertexType;
00369
00370 VertexType * vertices[topology::bndcells<triangle_tag, 0>::num
00371 + topology::bndcells<triangle_tag, 1>::num];
00372
00373
00374
00375
00376
00377
00378 VertexOnCellRange vertices_on_cell = viennagrid::ncells<0>(cell_in);
00379 VertexOnCellIterator vocit = vertices_on_cell.begin();
00380 vertices[0] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00381 vertices[1] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]); ++vocit;
00382 vertices[2] = &(viennagrid::ncells<0>(segment_out.domain())[vocit->id()]);
00383
00384
00385 EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(cell_in);
00386 EdgeOnCellIterator eocit = edges_on_cell.begin();
00387 vertices[3] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(*eocit)]); ++eocit;
00388 vertices[4] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(*eocit)]); ++eocit;
00389 vertices[5] = &(viennagrid::ncells<0>(segment_out.domain())[viennadata::access<refinement_key, std::size_t>(refinement_key())(*eocit)]);
00390
00391
00392
00393
00394 CellType new_cell;
00395 VertexType * cellvertices[topology::bndcells<triangle_tag, 0>::num];
00396
00397
00398 cellvertices[0] = vertices[0];
00399 cellvertices[1] = vertices[3];
00400 cellvertices[2] = vertices[4];
00401 new_cell.vertices(cellvertices);
00402 segment_out.push_back(new_cell);
00403
00404
00405 cellvertices[0] = vertices[3];
00406 cellvertices[1] = vertices[1];
00407 cellvertices[2] = vertices[5];
00408 new_cell.vertices(cellvertices);
00409 segment_out.push_back(new_cell);
00410
00411
00412 cellvertices[0] = vertices[4];
00413 cellvertices[1] = vertices[5];
00414 cellvertices[2] = vertices[2];
00415 new_cell.vertices(cellvertices);
00416 segment_out.push_back(new_cell);
00417
00418
00419 cellvertices[0] = vertices[4];
00420 cellvertices[1] = vertices[3];
00421 cellvertices[2] = vertices[5];
00422 new_cell.vertices(cellvertices);
00423 segment_out.push_back(new_cell);
00424
00425 }
00426
00427
00433 template <typename CellType, typename DomainTypeOut>
00434 static void apply(CellType const & cell_in, DomainTypeOut & segment_out)
00435 {
00436 typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type EdgeOnCellRange;
00437 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator;
00438
00439 std::size_t edges_to_refine = 0;
00440 EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(cell_in);
00441 for (EdgeOnCellIterator eocit = edges_on_cell.begin();
00442 eocit != edges_on_cell.end();
00443 ++eocit)
00444 {
00445 if (viennadata::access<refinement_key, bool>(refinement_key())(*eocit) == true)
00446 ++edges_to_refine;
00447 }
00448
00449 switch (edges_to_refine)
00450 {
00451 case 0: apply0(cell_in, segment_out); break;
00452 case 1: apply1(cell_in, segment_out); break;
00453 case 2: apply2(cell_in, segment_out); break;
00454 case 3: apply3(cell_in, segment_out); break;
00455 default:
00456 break;
00457 }
00458 }
00459 };
00460
00461
00462
00463
00464
00465
00466
00467 }
00468
00469 #endif
00470