Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include <vector>
34 : #include <algorithm>
35 : #include "grid_objects.h"
36 : #include "common/common.h"
37 : #include "tetrahedron_rules.h"
38 : #include "octahedron_rules.h"
39 : #include "pyramid_rules.h"
40 : #include "prism_rules.h"
41 : #include "hexahedron_rules.h"
42 : #include "grid_object_ids.h"
43 :
44 : //#include "../algorithms/geom_obj_util/geom_obj_util.h"
45 :
46 : using namespace std;
47 :
48 : namespace ug
49 : {
50 :
51 : class TetrahedronClass{
52 : public:
53 : enum{
54 : NUM_VERTICES = tet_rules::NUM_VERTICES,
55 : NUM_EDGES = tet_rules::NUM_EDGES,
56 : NUM_FACES = tet_rules::NUM_FACES,
57 : MAX_NUM_INDS_OUT = tet_rules::MAX_NUM_INDS_OUT
58 : };
59 : };
60 :
61 : class OctahedronClass{
62 : public:
63 : enum{
64 : NUM_VERTICES = oct_rules::NUM_VERTICES,
65 : NUM_EDGES = oct_rules::NUM_EDGES,
66 : NUM_FACES = oct_rules::NUM_FACES,
67 : MAX_NUM_INDS_OUT = oct_rules::MAX_NUM_INDS_OUT
68 : };
69 : };
70 :
71 : class PyramidClass{
72 : public:
73 : enum{
74 : NUM_VERTICES = pyra_rules::NUM_VERTICES,
75 : NUM_EDGES = pyra_rules::NUM_EDGES,
76 : NUM_FACES = pyra_rules::NUM_FACES,
77 : MAX_NUM_INDS_OUT = pyra_rules::MAX_NUM_INDS_OUT
78 : };
79 : };
80 :
81 : class PrismClass{
82 : public:
83 : enum{
84 : NUM_VERTICES = prism_rules::NUM_VERTICES,
85 : NUM_EDGES = prism_rules::NUM_EDGES,
86 : NUM_FACES = prism_rules::NUM_FACES,
87 : MAX_NUM_INDS_OUT = prism_rules::MAX_NUM_INDS_OUT
88 : };
89 : };
90 :
91 : class HexahedronClass{
92 : public:
93 : enum{
94 : NUM_VERTICES = hex_rules::NUM_VERTICES,
95 : NUM_EDGES = hex_rules::NUM_EDGES,
96 : NUM_FACES = hex_rules::NUM_FACES,
97 : MAX_NUM_INDS_OUT = hex_rules::MAX_NUM_INDS_OUT
98 : };
99 : };
100 :
101 : struct GridObjectInfo{
102 : public:
103 0 : static size_t num_vertices(int gridObjectID) {return inst().m_numVertices[gridObjectID];}
104 :
105 : private:
106 0 : static GridObjectInfo& inst(){
107 0 : static GridObjectInfo goi;
108 0 : return goi;
109 : }
110 :
111 0 : GridObjectInfo(){
112 0 : for(size_t i = 0; i < GOID_NUM_GRID_OBJECT_IDS; ++i)
113 0 : m_numVertices[i] = 0;
114 :
115 0 : m_numVertices[GOID_TETRAHEDRON] = 4;
116 0 : m_numVertices[GOID_PYRAMID] = 5;
117 0 : m_numVertices[GOID_PRISM] = 6;
118 0 : m_numVertices[GOID_OCTAHEDRON] = 6;
119 0 : m_numVertices[GOID_HEXAHEDRON] = 8;
120 : }
121 :
122 : size_t m_numVertices[GOID_NUM_GRID_OBJECT_IDS];
123 : };
124 :
125 : ////////////////////////////////////////////////////////////////////////
126 : ////////////////////////////////////////////////////////////////////////
127 : // TOOLS
128 : static
129 0 : void CreateVolumesFromElementIndexList (
130 : vector<Volume*>& volsOut,
131 : int* elemIndexList,
132 : int elemIndexListSize,
133 : Vertex** vrts)
134 : {
135 0 : VolumeDescriptor vd;
136 : volsOut.clear();
137 :
138 0 : for(int i = 0; i < elemIndexListSize;){
139 0 : int gridObjectID = elemIndexList[i++];
140 : size_t num = GridObjectInfo::num_vertices(gridObjectID);
141 0 : vd.set_num_vertices(num);
142 0 : for(size_t j = 0; j < num; ++j){
143 : assert(vrts[elemIndexList[i]]);
144 0 : vd.set_vertex(j, vrts[elemIndexList[i++]]);
145 : }
146 :
147 0 : switch(gridObjectID){
148 0 : case GOID_TETRAHEDRON: volsOut.push_back(new Tetrahedron(vd)); break;
149 0 : case GOID_PYRAMID: volsOut.push_back(new Pyramid(vd)); break;
150 0 : case GOID_PRISM: volsOut.push_back(new Prism(vd)); break;
151 0 : case GOID_HEXAHEDRON: volsOut.push_back(new Hexahedron(vd)); break;
152 0 : case GOID_OCTAHEDRON: volsOut.push_back(new Octahedron(vd)); break;
153 : }
154 : }
155 0 : }
156 :
157 : /** This refinement helper is called by the different refine implementations.
158 : * The last parameter is the actual refinement procedure as defined in
159 : * ug::tet_rules, ug::pyra_rules, ug::hex_rules or ug::prism_rules.
160 : */
161 : template <class TElemClass>
162 0 : static bool Refine(std::vector<Volume*>& vNewVolumesOut,
163 : Vertex** ppNewVertexOut,
164 : Vertex** newEdgeVertices,
165 : Vertex** newFaceVertices,
166 : Vertex* newVolumeVertex,
167 : const Vertex& prototypeVertex,
168 : Vertex** vrts,
169 : int (*funcRefine)(int*, int*, bool&, vector3*, bool*),
170 : vector3* corners = NULL,
171 : bool* isSnapPoint = NULL)
172 : {
173 : vNewVolumesOut.clear();
174 0 : *ppNewVertexOut = NULL;
175 :
176 : // allVrts is an array holding both, the vertices and the new edge-vertices.
177 : // we will index it with the results later on to create the new elements.
178 : const int allVrtsSize = TElemClass::NUM_VERTICES + TElemClass::NUM_EDGES
179 : + TElemClass::NUM_FACES + 1;
180 : Vertex* allVrts[allVrtsSize];
181 0 : for(int i = 0; i < TElemClass::NUM_VERTICES; ++i)
182 0 : allVrts[i] = vrts[i];
183 :
184 : // check which edge has to be refined, and which not
185 : int newEdgeVrts[TElemClass::NUM_EDGES];
186 0 : for(int i = 0; i < TElemClass::NUM_EDGES; ++i){
187 0 : allVrts[TElemClass::NUM_VERTICES + i] = newEdgeVertices[i];
188 0 : if(newEdgeVertices[i])
189 0 : newEdgeVrts[i] = 1;
190 : else
191 0 : newEdgeVrts[i] = 0;
192 : }
193 :
194 : // copy new face vertices to the allVrts array
195 0 : if(newFaceVertices){
196 0 : for(int i = 0; i < TElemClass::NUM_FACES; ++i){
197 0 : allVrts[TElemClass::NUM_VERTICES + TElemClass::NUM_EDGES + i] =
198 0 : newFaceVertices[i];
199 : }
200 : }
201 :
202 : // in this array we'll receive the new indices
203 : int newElemInds[TElemClass::MAX_NUM_INDS_OUT];
204 :
205 : // perform refine
206 0 : bool centerVrtRequired = false;
207 0 : int numElemInds = funcRefine(newElemInds, newEdgeVrts, centerVrtRequired, corners, isSnapPoint);
208 :
209 : assert(numElemInds != 0 && "PROBLEM in Refine(...): "
210 : "refine with 1 new edge vertex failed.");
211 :
212 0 : if(numElemInds == 0)
213 : return false;
214 :
215 : // if a new center vertex is required, then we'll create one now.
216 0 : if(centerVrtRequired){
217 0 : if(!newVolumeVertex)
218 : newVolumeVertex = static_cast<Vertex*>(
219 0 : prototypeVertex.create_empty_instance());
220 0 : *ppNewVertexOut = newVolumeVertex;
221 0 : allVrts[allVrtsSize - 1] = *ppNewVertexOut;
222 : }
223 :
224 : // debug log of inds
225 : /*
226 : UG_LOG("newElemInds:");
227 : for(int i = 0; i < numElemInds; ++i){
228 : UG_LOG(" " << newElemInds[i]);
229 : }
230 : UG_LOG(endl);
231 :
232 : UG_LOG("allVrts, (allVrtsSize: " << allVrtsSize << ") ; ");
233 : for(int i = 0; i < allVrtsSize; ++i){
234 : UG_LOG(" " << allVrts[i]);
235 : }
236 : UG_LOG(endl);
237 : */
238 :
239 0 : CreateVolumesFromElementIndexList(vNewVolumesOut, newElemInds, numElemInds, allVrts);
240 : // for(int i = 0; i < numElemInds;){
241 : // int gridObjectID = newElemInds[i++];
242 : // size_t num = GridObjectInfo::num_vertices(gridObjectID);
243 : // vd.set_num_vertices(num);
244 : // for(size_t j = 0; j < num; ++j){
245 : // assert(allVrts[newElemInds[i]]);
246 : // vd.set_vertex(j, allVrts[newElemInds[i++]]);
247 : // }
248 :
249 : // switch(gridObjectID){
250 : // case GOID_TETRAHEDRON: vNewVolumesOut.push_back(new Tetrahedron(vd)); break;
251 : // case GOID_PYRAMID: vNewVolumesOut.push_back(new Pyramid(vd)); break;
252 : // case GOID_PRISM: vNewVolumesOut.push_back(new Prism(vd)); break;
253 : // case GOID_HEXAHEDRON: vNewVolumesOut.push_back(new Hexahedron(vd)); break;
254 : // case GOID_OCTAHEDRON: vNewVolumesOut.push_back(new Octahedron(vd)); break;
255 : // }
256 : // }
257 :
258 : return true;
259 : }
260 :
261 : ////////////////////////////////////////////////////////////////////////
262 : ////////////////////////////////////////////////////////////////////////
263 : // VOLUMES
264 :
265 : ////////////////////////////////////////////////////////////////////////
266 : // TetrahedronDescriptor
267 0 : TetrahedronDescriptor::TetrahedronDescriptor(const TetrahedronDescriptor& td)
268 : {
269 0 : m_vertex[0] = td.vertex(0);
270 0 : m_vertex[1] = td.vertex(1);
271 0 : m_vertex[2] = td.vertex(2);
272 0 : m_vertex[3] = td.vertex(3);
273 0 : }
274 :
275 0 : TetrahedronDescriptor::TetrahedronDescriptor(const VolumeVertices& vv)
276 : {
277 : assert((vv.num_vertices() == 4) && "Bad number of vertices in volume-descriptor. Should be 4.");
278 0 : m_vertex[0] = vv.vertex(0);
279 0 : m_vertex[1] = vv.vertex(1);
280 0 : m_vertex[2] = vv.vertex(2);
281 0 : m_vertex[3] = vv.vertex(3);
282 0 : }
283 :
284 0 : TetrahedronDescriptor::TetrahedronDescriptor(Vertex* v1, Vertex* v2, Vertex* v3, Vertex* v4)
285 : {
286 0 : m_vertex[0] = v1;
287 0 : m_vertex[1] = v2;
288 0 : m_vertex[2] = v3;
289 0 : m_vertex[3] = v4;
290 0 : }
291 :
292 : ////////////////////////////////////////////////////////////////////////
293 : // Tetrahedron
294 0 : Tetrahedron::Tetrahedron(const TetrahedronDescriptor& td)
295 : {
296 0 : m_vertices[0] = td.vertex(0);
297 0 : m_vertices[1] = td.vertex(1);
298 0 : m_vertices[2] = td.vertex(2);
299 0 : m_vertices[3] = td.vertex(3);
300 0 : }
301 :
302 0 : Tetrahedron::Tetrahedron(Vertex* v1, Vertex* v2, Vertex* v3, Vertex* v4)
303 : {
304 0 : m_vertices[0] = v1;
305 0 : m_vertices[1] = v2;
306 0 : m_vertices[2] = v3;
307 0 : m_vertices[3] = v4;
308 0 : }
309 :
310 0 : EdgeDescriptor Tetrahedron::edge_desc(int index) const
311 : {
312 0 : EdgeDescriptor ed;
313 0 : edge_desc(index, ed);
314 0 : return ed;
315 : }
316 :
317 0 : void Tetrahedron::edge_desc(int index, EdgeDescriptor& edOut) const
318 : {
319 : using namespace tet_rules;
320 : assert(index >= 0 && index <= NUM_EDGES);
321 0 : edOut.set_vertices(m_vertices[EDGE_VRT_INDS[index][0]],
322 0 : m_vertices[EDGE_VRT_INDS[index][1]]);
323 0 : }
324 :
325 0 : uint Tetrahedron::num_edges() const
326 : {
327 0 : return 6;
328 : }
329 :
330 0 : FaceDescriptor Tetrahedron::face_desc(int index) const
331 : {
332 0 : FaceDescriptor fd;
333 0 : face_desc(index, fd);
334 0 : return fd;
335 : }
336 :
337 0 : void Tetrahedron::face_desc(int index, FaceDescriptor& fdOut) const
338 : {
339 : using namespace tet_rules;
340 : assert(index >= 0 && index < NUM_FACES);
341 :
342 : fdOut.set_num_vertices(3);
343 0 : fdOut.set_vertex(0, m_vertices[FACE_VRT_INDS[index][0]]);
344 0 : fdOut.set_vertex(1, m_vertices[FACE_VRT_INDS[index][2]]);
345 0 : fdOut.set_vertex(2, m_vertices[FACE_VRT_INDS[index][1]]);
346 0 : }
347 :
348 0 : uint Tetrahedron::num_faces() const
349 : {
350 0 : return 4;
351 : }
352 :
353 0 : Edge* Tetrahedron::create_edge(int index)
354 : {
355 : using namespace tet_rules;
356 : assert(index >= 0 && index < NUM_EDGES);
357 0 : const int* e = EDGE_VRT_INDS[index];
358 0 : return new RegularEdge(m_vertices[e[0]], m_vertices[e[1]]);
359 : }
360 :
361 0 : Face* Tetrahedron::create_face(int index)
362 : {
363 : using namespace tet_rules;
364 : assert(index >= 0 && index < NUM_FACES);
365 :
366 0 : const int* f = FACE_VRT_INDS[index];
367 0 : return new Triangle(m_vertices[f[0]], m_vertices[f[2]], m_vertices[f[1]]);
368 : }
369 :
370 0 : void Tetrahedron::
371 : get_vertex_indices_of_edge (size_t& ind1Out,
372 : size_t& ind2Out,
373 : size_t edgeInd) const
374 : {
375 : assert(edgeInd >= 0 && edgeInd < 6);
376 0 : ind1Out = tet_rules::EDGE_VRT_INDS[edgeInd][0];
377 0 : ind2Out = tet_rules::EDGE_VRT_INDS[edgeInd][1];
378 0 : }
379 :
380 0 : void Tetrahedron::
381 : get_vertex_indices_of_face (std::vector<size_t>& indsOut,
382 : size_t side) const
383 : {
384 : assert(side >= 0 && side < 4);
385 0 : indsOut.resize(3);
386 0 : indsOut[0] = tet_rules::FACE_VRT_INDS[side][0];
387 0 : indsOut[1] = tet_rules::FACE_VRT_INDS[side][2];
388 0 : indsOut[2] = tet_rules::FACE_VRT_INDS[side][1];
389 0 : }
390 :
391 0 : int Tetrahedron::
392 : get_edge_index_from_vertices( const size_t vi0,
393 : const size_t vi1) const
394 : {
395 0 : return tet_rules::EDGE_FROM_VRTS[vi0][vi1];
396 : }
397 :
398 0 : int Tetrahedron::
399 : get_face_edge_index(const size_t faceInd,
400 : const size_t faceEdgeInd) const
401 : {
402 0 : return tet_rules::FACE_EDGE_INDS[faceInd][2 - faceEdgeInd];
403 : }
404 :
405 0 : bool Tetrahedron::collapse_edge(std::vector<Volume*>& vNewVolumesOut,
406 : int edgeIndex, Vertex* newVertex,
407 : std::vector<Vertex*>* pvSubstituteVertices)
408 : {
409 : // if an edge of a tetrahedron is collapsed, nothing remains.
410 : vNewVolumesOut.clear();
411 0 : return true;
412 : }
413 :
414 0 : std::pair<GridBaseObjectId, int> Tetrahedron::
415 : get_opposing_object(Vertex* vrt) const
416 : {
417 : using namespace tet_rules;
418 0 : for(int i = 0; i < tet_rules::NUM_VERTICES; ++i){
419 0 : if(vrt == m_vertices[i]){
420 0 : return make_pair(static_cast<GridBaseObjectId>(OPPOSED_OBJECT[i][0]),
421 0 : OPPOSED_OBJECT[i][1]);
422 : }
423 : }
424 :
425 0 : UG_THROW("Specified vertex is not part of this element.");
426 : }
427 :
428 0 : bool Tetrahedron::refine(std::vector<Volume*>& vNewVolumesOut,
429 : Vertex** ppNewVertexOut,
430 : Vertex** newEdgeVertices,
431 : Vertex** newFaceVertices,
432 : Vertex* newVolumeVertex,
433 : const Vertex& prototypeVertex,
434 : Vertex** pSubstituteVertices,
435 : vector3* corners,
436 : bool* isSnapPoint)
437 : {
438 : // handle substitute vertices.
439 : Vertex** vrts;
440 0 : if(pSubstituteVertices)
441 : vrts = pSubstituteVertices;
442 : else
443 0 : vrts = m_vertices;
444 :
445 0 : return Refine<TetrahedronClass>(vNewVolumesOut, ppNewVertexOut,
446 : newEdgeVertices, newFaceVertices,
447 : newVolumeVertex, prototypeVertex,
448 : vrts, tet_rules::Refine, corners,
449 0 : isSnapPoint);
450 : }
451 :
452 :
453 0 : bool Tetrahedron::is_regular_ref_rule(int edgeMarks) const
454 : {
455 0 : return tet_rules::IsRegularRefRule(edgeMarks);
456 : }
457 :
458 :
459 0 : void Tetrahedron::get_flipped_orientation(VolumeDescriptor& vdOut) const
460 : {
461 : // in order to flip a tetrahedrons orientation, we have to invert the order
462 : // of the base-vertices
463 : vdOut.set_num_vertices(4);
464 0 : vdOut.set_vertex(0, vertex(2));
465 0 : vdOut.set_vertex(1, vertex(1));
466 0 : vdOut.set_vertex(2, vertex(0));
467 0 : vdOut.set_vertex(3, vertex(3));
468 0 : }
469 :
470 :
471 : ////////////////////////////////////////////////////////////////////////
472 : // HexahedronDescriptor
473 0 : HexahedronDescriptor::HexahedronDescriptor(const HexahedronDescriptor& td)
474 : {
475 0 : m_vertex[0] = td.vertex(0);
476 0 : m_vertex[1] = td.vertex(1);
477 0 : m_vertex[2] = td.vertex(2);
478 0 : m_vertex[3] = td.vertex(3);
479 0 : m_vertex[4] = td.vertex(4);
480 0 : m_vertex[5] = td.vertex(5);
481 0 : m_vertex[6] = td.vertex(6);
482 0 : m_vertex[7] = td.vertex(7);
483 0 : }
484 :
485 0 : HexahedronDescriptor::HexahedronDescriptor(const VolumeVertices& vv)
486 : {
487 : assert((vv.num_vertices() == 8) && "Bad number of vertices in volume-descriptor. Should be 8.");
488 0 : m_vertex[0] = vv.vertex(0);
489 0 : m_vertex[1] = vv.vertex(1);
490 0 : m_vertex[2] = vv.vertex(2);
491 0 : m_vertex[3] = vv.vertex(3);
492 0 : m_vertex[4] = vv.vertex(4);
493 0 : m_vertex[5] = vv.vertex(5);
494 0 : m_vertex[6] = vv.vertex(6);
495 0 : m_vertex[7] = vv.vertex(7);
496 0 : }
497 :
498 0 : HexahedronDescriptor::HexahedronDescriptor(Vertex* v1, Vertex* v2, Vertex* v3, Vertex* v4,
499 0 : Vertex* v5, Vertex* v6, Vertex* v7, Vertex* v8)
500 : {
501 0 : m_vertex[0] = v1;
502 0 : m_vertex[1] = v2;
503 0 : m_vertex[2] = v3;
504 0 : m_vertex[3] = v4;
505 0 : m_vertex[4] = v5;
506 0 : m_vertex[5] = v6;
507 0 : m_vertex[6] = v7;
508 0 : m_vertex[7] = v8;
509 0 : }
510 :
511 : ////////////////////////////////////////////////////////////////////////
512 : // Hexahedron
513 0 : Hexahedron::Hexahedron(const HexahedronDescriptor& td)
514 : {
515 0 : m_vertices[0] = td.vertex(0);
516 0 : m_vertices[1] = td.vertex(1);
517 0 : m_vertices[2] = td.vertex(2);
518 0 : m_vertices[3] = td.vertex(3);
519 0 : m_vertices[4] = td.vertex(4);
520 0 : m_vertices[5] = td.vertex(5);
521 0 : m_vertices[6] = td.vertex(6);
522 0 : m_vertices[7] = td.vertex(7);
523 0 : }
524 :
525 0 : Hexahedron::Hexahedron(Vertex* v1, Vertex* v2, Vertex* v3, Vertex* v4,
526 0 : Vertex* v5, Vertex* v6, Vertex* v7, Vertex* v8)
527 : {
528 0 : m_vertices[0] = v1;
529 0 : m_vertices[1] = v2;
530 0 : m_vertices[2] = v3;
531 0 : m_vertices[3] = v4;
532 0 : m_vertices[4] = v5;
533 0 : m_vertices[5] = v6;
534 0 : m_vertices[6] = v7;
535 0 : m_vertices[7] = v8;
536 0 : }
537 :
538 0 : EdgeDescriptor Hexahedron::edge_desc(int index) const
539 : {
540 0 : EdgeDescriptor ed;
541 0 : edge_desc(index, ed);
542 0 : return ed;
543 : }
544 :
545 0 : void Hexahedron::edge_desc(int index, EdgeDescriptor& edOut) const
546 : {
547 : using namespace hex_rules;
548 : assert(index >= 0 && index < NUM_EDGES);
549 0 : const int* e = EDGE_VRT_INDS[index];
550 0 : edOut.set_vertices(m_vertices[e[0]],
551 0 : m_vertices[e[1]]);
552 0 : }
553 :
554 0 : uint Hexahedron::num_edges() const
555 : {
556 0 : return 12;
557 : }
558 :
559 0 : FaceDescriptor Hexahedron::face_desc(int index) const
560 : {
561 0 : FaceDescriptor fd;
562 0 : face_desc(index, fd);
563 0 : return fd;
564 : }
565 :
566 0 : void Hexahedron::face_desc(int index, FaceDescriptor& fdOut) const
567 : {
568 : using namespace hex_rules;
569 : assert(index >= 0 && index < NUM_FACES);
570 :
571 0 : const int* f = FACE_VRT_INDS[index];
572 :
573 : fdOut.set_num_vertices(4);
574 0 : fdOut.set_vertex(0, m_vertices[f[0]]);
575 0 : fdOut.set_vertex(1, m_vertices[f[3]]);
576 0 : fdOut.set_vertex(2, m_vertices[f[2]]);
577 0 : fdOut.set_vertex(3, m_vertices[f[1]]);
578 0 : }
579 :
580 0 : uint Hexahedron::num_faces() const
581 : {
582 0 : return 6;
583 : }
584 :
585 0 : Edge* Hexahedron::create_edge(int index)
586 : {
587 : using namespace hex_rules;
588 : assert(index >= 0 && index < NUM_EDGES);
589 0 : const int* e = EDGE_VRT_INDS[index];
590 0 : return new RegularEdge(m_vertices[e[0]], m_vertices[e[1]]);
591 : }
592 :
593 0 : Face* Hexahedron::create_face(int index)
594 : {
595 : using namespace hex_rules;
596 : assert(index >= 0 && index < NUM_FACES);
597 :
598 0 : const int* f = FACE_VRT_INDS[index];
599 0 : return new Quadrilateral(m_vertices[f[0]], m_vertices[f[3]],
600 0 : m_vertices[f[2]], m_vertices[f[1]]);
601 : }
602 :
603 0 : void Hexahedron::
604 : get_vertex_indices_of_edge (size_t& ind1Out,
605 : size_t& ind2Out,
606 : size_t edgeInd) const
607 : {
608 : assert(edgeInd >= 0 && edgeInd < hex_rules::NUM_EDGES);
609 0 : ind1Out = hex_rules::EDGE_VRT_INDS[edgeInd][0];
610 0 : ind2Out = hex_rules::EDGE_VRT_INDS[edgeInd][1];
611 0 : }
612 :
613 0 : void Hexahedron::
614 : get_vertex_indices_of_face (std::vector<size_t>& indsOut,
615 : size_t side) const
616 : {
617 : assert(side >= 0 && side < hex_rules::NUM_FACES);
618 0 : indsOut.resize(4);
619 0 : indsOut[0] = hex_rules::FACE_VRT_INDS[side][0];
620 0 : indsOut[1] = hex_rules::FACE_VRT_INDS[side][3];
621 0 : indsOut[2] = hex_rules::FACE_VRT_INDS[side][2];
622 0 : indsOut[3] = hex_rules::FACE_VRT_INDS[side][1];
623 0 : }
624 :
625 0 : int Hexahedron::
626 : get_edge_index_from_vertices( const size_t vi0,
627 : const size_t vi1) const
628 : {
629 0 : return hex_rules::EDGE_FROM_VRTS[vi0][vi1];
630 : }
631 :
632 0 : int Hexahedron::
633 : get_face_edge_index(const size_t faceInd,
634 : const size_t faceEdgeInd) const
635 : {
636 0 : return hex_rules::FACE_EDGE_INDS[faceInd][3 - faceEdgeInd];
637 : }
638 :
639 0 : bool Hexahedron::get_opposing_side(FaceVertices* f, FaceDescriptor& fdOut) const
640 : {
641 : using namespace hex_rules;
642 0 : int localInd = get_local_side_index(f);
643 0 : if(localInd == -1)
644 : return false;
645 :
646 0 : face_desc(OPPOSED_FACE[localInd], fdOut);
647 0 : return true;
648 : }
649 :
650 0 : std::pair<GridBaseObjectId, int> Hexahedron::
651 : get_opposing_object(Vertex* vrt) const
652 : {
653 : using namespace hex_rules;
654 0 : for(int i = 0; i < hex_rules::NUM_VERTICES; ++i){
655 0 : if(vrt == m_vertices[i]){
656 0 : return make_pair(static_cast<GridBaseObjectId>(OPPOSED_OBJECT[i][0]),
657 0 : OPPOSED_OBJECT[i][1]);
658 : }
659 : }
660 :
661 0 : UG_THROW("Specified vertex is not part of this element.");
662 : }
663 :
664 0 : bool Hexahedron::collapse_edge(std::vector<Volume*>& vNewVolumesOut,
665 : int edgeIndex, Vertex* newVertex,
666 : std::vector<Vertex*>* pvSubstituteVertices)
667 : {
668 : // NOT YET SUPPORTED!
669 : //TODO: implement Hexahedron::collapse_edge
670 : vNewVolumesOut.clear();
671 : UG_LOG("edge-collapse for hexahedrons not yet implemented... sorry\n");
672 0 : return false;
673 : }
674 :
675 0 : bool Hexahedron::refine(std::vector<Volume*>& vNewVolumesOut,
676 : Vertex** ppNewVertexOut,
677 : Vertex** newEdgeVertices,
678 : Vertex** newFaceVertices,
679 : Vertex* newVolumeVertex,
680 : const Vertex& prototypeVertex,
681 : Vertex** pSubstituteVertices,
682 : vector3* corners,
683 : bool* isSnapPoint)
684 : {
685 : // handle substitute vertices.
686 : Vertex** vrts;
687 0 : if(pSubstituteVertices)
688 : vrts = pSubstituteVertices;
689 : else
690 0 : vrts = m_vertices;
691 :
692 0 : return Refine<HexahedronClass>(vNewVolumesOut, ppNewVertexOut,
693 : newEdgeVertices, newFaceVertices,
694 : newVolumeVertex, prototypeVertex,
695 : vrts, hex_rules::Refine, corners,
696 0 : isSnapPoint);
697 : }
698 :
699 0 : bool Hexahedron::is_regular_ref_rule(int edgeMarks) const
700 : {
701 0 : return hex_rules::IsRegularRefRule(edgeMarks);
702 : }
703 :
704 0 : void Hexahedron::get_flipped_orientation(VolumeDescriptor& vdOut) const
705 : {
706 : // in order to flip a hexahedrons orientation, we have to swap
707 : // the bottom and top vertices
708 : vdOut.set_num_vertices(8);
709 0 : vdOut.set_vertex(0, vertex(4));
710 0 : vdOut.set_vertex(1, vertex(5));
711 0 : vdOut.set_vertex(2, vertex(6));
712 0 : vdOut.set_vertex(3, vertex(7));
713 0 : vdOut.set_vertex(4, vertex(0));
714 0 : vdOut.set_vertex(5, vertex(1));
715 0 : vdOut.set_vertex(6, vertex(2));
716 0 : vdOut.set_vertex(7, vertex(3));
717 0 : }
718 :
719 :
720 :
721 : ////////////////////////////////////////////////////////////////////////
722 : ////////////////////////////////////////////////////////////////////////
723 : // PrismDescriptor
724 0 : PrismDescriptor::PrismDescriptor(const PrismDescriptor& td)
725 : {
726 0 : m_vertex[0] = td.vertex(0);
727 0 : m_vertex[1] = td.vertex(1);
728 0 : m_vertex[2] = td.vertex(2);
729 0 : m_vertex[3] = td.vertex(3);
730 0 : m_vertex[4] = td.vertex(4);
731 0 : m_vertex[5] = td.vertex(5);
732 0 : }
733 :
734 0 : PrismDescriptor::PrismDescriptor(const VolumeVertices& vv)
735 : {
736 : assert((vv.num_vertices() == 6) && "Bad number of vertices in volume-descriptor. Should be 6.");
737 0 : m_vertex[0] = vv.vertex(0);
738 0 : m_vertex[1] = vv.vertex(1);
739 0 : m_vertex[2] = vv.vertex(2);
740 0 : m_vertex[3] = vv.vertex(3);
741 0 : m_vertex[4] = vv.vertex(4);
742 0 : m_vertex[5] = vv.vertex(5);
743 0 : }
744 :
745 0 : PrismDescriptor::PrismDescriptor(Vertex* v1, Vertex* v2, Vertex* v3,
746 0 : Vertex* v4, Vertex* v5, Vertex* v6)
747 : {
748 0 : m_vertex[0] = v1;
749 0 : m_vertex[1] = v2;
750 0 : m_vertex[2] = v3;
751 0 : m_vertex[3] = v4;
752 0 : m_vertex[4] = v5;
753 0 : m_vertex[5] = v6;
754 0 : }
755 :
756 : ////////////////////////////////////////////////////////////////////////
757 : // Prism
758 0 : Prism::Prism(const PrismDescriptor& td)
759 : {
760 0 : m_vertices[0] = td.vertex(0);
761 0 : m_vertices[1] = td.vertex(1);
762 0 : m_vertices[2] = td.vertex(2);
763 0 : m_vertices[3] = td.vertex(3);
764 0 : m_vertices[4] = td.vertex(4);
765 0 : m_vertices[5] = td.vertex(5);
766 0 : }
767 :
768 0 : Prism::Prism(Vertex* v1, Vertex* v2, Vertex* v3,
769 0 : Vertex* v4, Vertex* v5, Vertex* v6)
770 : {
771 0 : m_vertices[0] = v1;
772 0 : m_vertices[1] = v2;
773 0 : m_vertices[2] = v3;
774 0 : m_vertices[3] = v4;
775 0 : m_vertices[4] = v5;
776 0 : m_vertices[5] = v6;
777 0 : }
778 :
779 0 : EdgeDescriptor Prism::edge_desc(int index) const
780 : {
781 0 : EdgeDescriptor ed;
782 0 : edge_desc(index, ed);
783 0 : return ed;
784 : }
785 :
786 0 : void Prism::edge_desc(int index, EdgeDescriptor& edOut) const
787 : {
788 : using namespace prism_rules;
789 : assert(index >= 0 && index < NUM_EDGES);
790 0 : const int* e = EDGE_VRT_INDS[index];
791 0 : edOut.set_vertices(m_vertices[e[0]],
792 0 : m_vertices[e[1]]);
793 0 : }
794 :
795 0 : uint Prism::num_edges() const
796 : {
797 0 : return 9;
798 : }
799 :
800 0 : FaceDescriptor Prism::face_desc(int index) const
801 : {
802 0 : FaceDescriptor fd;
803 0 : face_desc(index, fd);
804 0 : return fd;
805 : }
806 :
807 0 : void Prism::face_desc(int index, FaceDescriptor& fdOut) const
808 : {
809 : using namespace prism_rules;
810 : assert(index >= 0 && index < NUM_FACES);
811 :
812 0 : const int* f = FACE_VRT_INDS[index];
813 0 : if(f[3] == -1){
814 : fdOut.set_num_vertices(3);
815 0 : fdOut.set_vertex(0, m_vertices[f[0]]);
816 0 : fdOut.set_vertex(1, m_vertices[f[2]]);
817 0 : fdOut.set_vertex(2, m_vertices[f[1]]);
818 : }
819 : else{
820 : fdOut.set_num_vertices(4);
821 0 : fdOut.set_vertex(0, m_vertices[f[0]]);
822 0 : fdOut.set_vertex(1, m_vertices[f[3]]);
823 0 : fdOut.set_vertex(2, m_vertices[f[2]]);
824 0 : fdOut.set_vertex(3, m_vertices[f[1]]);
825 : }
826 0 : }
827 :
828 0 : uint Prism::num_faces() const
829 : {
830 0 : return 5;
831 : }
832 :
833 0 : Edge* Prism::create_edge(int index)
834 : {
835 : using namespace prism_rules;
836 : assert(index >= 0 && index < NUM_EDGES);
837 0 : const int* e = EDGE_VRT_INDS[index];
838 0 : return new RegularEdge(m_vertices[e[0]], m_vertices[e[1]]);
839 : }
840 :
841 0 : Face* Prism::create_face(int index)
842 : {
843 : using namespace prism_rules;
844 : assert(index >= 0 && index < NUM_FACES);
845 :
846 0 : const int* f = FACE_VRT_INDS[index];
847 0 : if(f[3] == -1){
848 0 : return new Triangle(m_vertices[f[0]], m_vertices[f[2]],
849 0 : m_vertices[f[1]]);
850 : }
851 : else{
852 0 : return new Quadrilateral(m_vertices[f[0]], m_vertices[f[3]],
853 0 : m_vertices[f[2]], m_vertices[f[1]]);
854 : }
855 : }
856 :
857 0 : void Prism::
858 : get_vertex_indices_of_edge (size_t& ind1Out,
859 : size_t& ind2Out,
860 : size_t edgeInd) const
861 : {
862 : assert(edgeInd >= 0 && edgeInd < prism_rules::NUM_EDGES);
863 0 : ind1Out = prism_rules::EDGE_VRT_INDS[edgeInd][0];
864 0 : ind2Out = prism_rules::EDGE_VRT_INDS[edgeInd][1];
865 0 : }
866 :
867 0 : void Prism::
868 : get_vertex_indices_of_face (std::vector<size_t>& indsOut,
869 : size_t side) const
870 : {
871 : assert(side >= 0 && side < prism_rules::NUM_FACES);
872 :
873 0 : if(prism_rules::FACE_VRT_INDS[side][3] == -1){
874 0 : indsOut.resize(3);
875 0 : indsOut[0] = prism_rules::FACE_VRT_INDS[side][0];
876 0 : indsOut[1] = prism_rules::FACE_VRT_INDS[side][2];
877 0 : indsOut[2] = prism_rules::FACE_VRT_INDS[side][1];
878 : }
879 : else{
880 0 : indsOut.resize(4);
881 0 : indsOut[0] = prism_rules::FACE_VRT_INDS[side][0];
882 0 : indsOut[1] = prism_rules::FACE_VRT_INDS[side][3];
883 0 : indsOut[2] = prism_rules::FACE_VRT_INDS[side][2];
884 0 : indsOut[3] = prism_rules::FACE_VRT_INDS[side][1];
885 : }
886 0 : }
887 :
888 0 : int Prism::
889 : get_edge_index_from_vertices( const size_t vi0,
890 : const size_t vi1) const
891 : {
892 0 : return prism_rules::EDGE_FROM_VRTS[vi0][vi1];
893 : }
894 :
895 0 : int Prism::
896 : get_face_edge_index(const size_t faceInd,
897 : const size_t faceEdgeInd) const
898 : {
899 0 : if(prism_rules::FACE_EDGE_INDS[faceInd][3] == -1)
900 0 : return prism_rules::FACE_EDGE_INDS[faceInd][2 - faceEdgeInd];
901 : else
902 0 : return prism_rules::FACE_EDGE_INDS[faceInd][3 - faceEdgeInd];
903 : }
904 :
905 0 : bool Prism::get_opposing_side(FaceVertices* f, FaceDescriptor& fdOut) const
906 : {
907 : using namespace prism_rules;
908 0 : int localInd = get_local_side_index(f);
909 0 : if(localInd == -1)
910 : return false;
911 :
912 0 : int opposedInd = OPPOSED_FACE[localInd];
913 0 : if(opposedInd == -1)
914 : return false;
915 :
916 0 : face_desc(opposedInd, fdOut);
917 0 : return true;
918 : }
919 :
920 0 : std::pair<GridBaseObjectId, int> Prism::
921 : get_opposing_object(Vertex* vrt) const
922 : {
923 : using namespace prism_rules;
924 0 : for(int i = 0; i < prism_rules::NUM_VERTICES; ++i){
925 0 : if(vrt == m_vertices[i]){
926 0 : return make_pair(static_cast<GridBaseObjectId>(OPPOSED_OBJECT[i][0]),
927 0 : OPPOSED_OBJECT[i][1]);
928 : }
929 : }
930 0 : UG_THROW("Specified vertex is not part of this element.");
931 : }
932 :
933 0 : bool Prism::collapse_edge(std::vector<Volume*>& vNewVolumesOut,
934 : int edgeIndex, Vertex* newVertex,
935 : std::vector<Vertex*>* pvSubstituteVertices)
936 : {
937 : using namespace prism_rules;
938 :
939 : int elemInds[MAX_NUM_COLLAPSE_INDS_OUT];
940 0 : int elemIndsSize = CollapseEdge(
941 : elemInds,
942 0 : EDGE_VRT_INDS[edgeIndex][0],
943 0 : EDGE_VRT_INDS[edgeIndex][1]);
944 :
945 0 : if(elemIndsSize > 0){
946 : Vertex** vrts;
947 0 : if(pvSubstituteVertices)
948 : vrts = &pvSubstituteVertices->front();
949 : else
950 0 : vrts = m_vertices;
951 :
952 0 : CreateVolumesFromElementIndexList(
953 : vNewVolumesOut,
954 : elemInds,
955 : elemIndsSize,
956 : vrts);
957 0 : return !vNewVolumesOut.empty();
958 : }
959 : else{
960 : vNewVolumesOut.clear();
961 0 : return false;
962 : }
963 : }
964 :
965 0 : bool Prism::refine(std::vector<Volume*>& vNewVolumesOut,
966 : Vertex** ppNewVertexOut,
967 : Vertex** newEdgeVertices,
968 : Vertex** newFaceVertices,
969 : Vertex* newVolumeVertex,
970 : const Vertex& prototypeVertex,
971 : Vertex** pSubstituteVertices,
972 : vector3* corners,
973 : bool* isSnapPoint)
974 : {
975 : // handle substitute vertices.
976 : Vertex** vrts;
977 0 : if(pSubstituteVertices)
978 : vrts = pSubstituteVertices;
979 : else
980 0 : vrts = m_vertices;
981 :
982 0 : return Refine<PrismClass>(vNewVolumesOut, ppNewVertexOut,
983 : newEdgeVertices, newFaceVertices,
984 : newVolumeVertex, prototypeVertex,
985 : vrts, prism_rules::Refine, corners,
986 0 : isSnapPoint);
987 : }
988 :
989 0 : bool Prism::is_regular_ref_rule(int edgeMarks) const
990 : {
991 0 : return prism_rules::IsRegularRefRule(edgeMarks);
992 : }
993 :
994 0 : void Prism::get_flipped_orientation(VolumeDescriptor& vdOut) const
995 : {
996 : // in order to flip a prisms orientation, we have to swap
997 : // the bottom and top vertices
998 : vdOut.set_num_vertices(6);
999 0 : vdOut.set_vertex(0, vertex(3));
1000 0 : vdOut.set_vertex(1, vertex(4));
1001 0 : vdOut.set_vertex(2, vertex(5));
1002 0 : vdOut.set_vertex(3, vertex(0));
1003 0 : vdOut.set_vertex(4, vertex(1));
1004 0 : vdOut.set_vertex(5, vertex(2));
1005 0 : }
1006 :
1007 :
1008 : ////////////////////////////////////////////////////////////////////////
1009 : ////////////////////////////////////////////////////////////////////////
1010 : // PyramidDescriptor
1011 0 : PyramidDescriptor::PyramidDescriptor(const PyramidDescriptor& td)
1012 : {
1013 0 : m_vertex[0] = td.vertex(0);
1014 0 : m_vertex[1] = td.vertex(1);
1015 0 : m_vertex[2] = td.vertex(2);
1016 0 : m_vertex[3] = td.vertex(3);
1017 0 : m_vertex[4] = td.vertex(4);
1018 0 : }
1019 :
1020 0 : PyramidDescriptor::PyramidDescriptor(const VolumeVertices& vv)
1021 : {
1022 : assert((vv.num_vertices() == 5) && "Bad number of vertices in volume-descriptor. Should be 5.");
1023 0 : m_vertex[0] = vv.vertex(0);
1024 0 : m_vertex[1] = vv.vertex(1);
1025 0 : m_vertex[2] = vv.vertex(2);
1026 0 : m_vertex[3] = vv.vertex(3);
1027 0 : m_vertex[4] = vv.vertex(4);
1028 0 : }
1029 :
1030 0 : PyramidDescriptor::PyramidDescriptor(Vertex* v1, Vertex* v2, Vertex* v3,
1031 0 : Vertex* v4, Vertex* v5)
1032 : {
1033 0 : m_vertex[0] = v1;
1034 0 : m_vertex[1] = v2;
1035 0 : m_vertex[2] = v3;
1036 0 : m_vertex[3] = v4;
1037 0 : m_vertex[4] = v5;
1038 0 : }
1039 :
1040 : ////////////////////////////////////////////////////////////////////////
1041 : // Pyramid
1042 0 : Pyramid::Pyramid(const PyramidDescriptor& td)
1043 : {
1044 0 : m_vertices[0] = td.vertex(0);
1045 0 : m_vertices[1] = td.vertex(1);
1046 0 : m_vertices[2] = td.vertex(2);
1047 0 : m_vertices[3] = td.vertex(3);
1048 0 : m_vertices[4] = td.vertex(4);
1049 0 : }
1050 :
1051 0 : Pyramid::Pyramid(Vertex* v1, Vertex* v2, Vertex* v3,
1052 0 : Vertex* v4, Vertex* v5)
1053 : {
1054 0 : m_vertices[0] = v1;
1055 0 : m_vertices[1] = v2;
1056 0 : m_vertices[2] = v3;
1057 0 : m_vertices[3] = v4;
1058 0 : m_vertices[4] = v5;
1059 0 : }
1060 :
1061 0 : EdgeDescriptor Pyramid::edge_desc(int index) const
1062 : {
1063 0 : EdgeDescriptor ed;
1064 0 : edge_desc(index, ed);
1065 0 : return ed;
1066 : }
1067 :
1068 0 : void Pyramid::edge_desc(int index, EdgeDescriptor& edOut) const
1069 : {
1070 : using namespace pyra_rules;
1071 : assert(index >= 0 && index < NUM_EDGES);
1072 0 : const int* e = EDGE_VRT_INDS[index];
1073 0 : edOut.set_vertices(m_vertices[e[0]],
1074 0 : m_vertices[e[1]]);
1075 0 : }
1076 :
1077 0 : uint Pyramid::num_edges() const
1078 : {
1079 0 : return 8;
1080 : }
1081 :
1082 0 : FaceDescriptor Pyramid::face_desc(int index) const
1083 : {
1084 0 : FaceDescriptor fd;
1085 0 : face_desc(index, fd);
1086 0 : return fd;
1087 : }
1088 :
1089 0 : void Pyramid::face_desc(int index, FaceDescriptor& fdOut) const
1090 : {
1091 : using namespace pyra_rules;
1092 : assert(index >= 0 && index < NUM_FACES);
1093 :
1094 0 : const int* f = FACE_VRT_INDS[index];
1095 0 : if(f[3] == -1){
1096 : fdOut.set_num_vertices(3);
1097 0 : fdOut.set_vertex(0, m_vertices[f[0]]);
1098 0 : fdOut.set_vertex(1, m_vertices[f[2]]);
1099 0 : fdOut.set_vertex(2, m_vertices[f[1]]);
1100 : }
1101 : else{
1102 : fdOut.set_num_vertices(4);
1103 0 : fdOut.set_vertex(0, m_vertices[f[0]]);
1104 0 : fdOut.set_vertex(1, m_vertices[f[3]]);
1105 0 : fdOut.set_vertex(2, m_vertices[f[2]]);
1106 0 : fdOut.set_vertex(3, m_vertices[f[1]]);
1107 : }
1108 0 : }
1109 :
1110 0 : uint Pyramid::num_faces() const
1111 : {
1112 0 : return 5;
1113 : }
1114 :
1115 0 : Edge* Pyramid::create_edge(int index)
1116 : {
1117 : using namespace pyra_rules;
1118 : assert(index >= 0 && index < NUM_EDGES);
1119 0 : const int* e = EDGE_VRT_INDS[index];
1120 0 : return new RegularEdge(m_vertices[e[0]], m_vertices[e[1]]);
1121 : }
1122 :
1123 0 : Face* Pyramid::create_face(int index)
1124 : {
1125 : using namespace pyra_rules;
1126 : assert(index >= 0 && index < NUM_FACES);
1127 :
1128 0 : const int* f = FACE_VRT_INDS[index];
1129 0 : if(f[3] == -1){
1130 0 : return new Triangle(m_vertices[f[0]], m_vertices[f[2]],
1131 0 : m_vertices[f[1]]);
1132 : }
1133 : else{
1134 0 : return new Quadrilateral(m_vertices[f[0]], m_vertices[f[3]],
1135 0 : m_vertices[f[2]], m_vertices[f[1]]);
1136 : }
1137 : }
1138 :
1139 0 : void Pyramid::
1140 : get_vertex_indices_of_edge (size_t& ind1Out,
1141 : size_t& ind2Out,
1142 : size_t edgeInd) const
1143 : {
1144 : assert(edgeInd >= 0 && edgeInd < pyra_rules::NUM_EDGES);
1145 0 : ind1Out = pyra_rules::EDGE_VRT_INDS[edgeInd][0];
1146 0 : ind2Out = pyra_rules::EDGE_VRT_INDS[edgeInd][1];
1147 0 : }
1148 :
1149 0 : void Pyramid::
1150 : get_vertex_indices_of_face (std::vector<size_t>& indsOut,
1151 : size_t side) const
1152 : {
1153 : assert(side >= 0 && side < pyra_rules::NUM_FACES);
1154 :
1155 0 : if(pyra_rules::FACE_VRT_INDS[side][3] == -1){
1156 0 : indsOut.resize(3);
1157 0 : indsOut[0] = pyra_rules::FACE_VRT_INDS[side][0];
1158 0 : indsOut[1] = pyra_rules::FACE_VRT_INDS[side][2];
1159 0 : indsOut[2] = pyra_rules::FACE_VRT_INDS[side][1];
1160 : }
1161 : else{
1162 0 : indsOut.resize(4);
1163 0 : indsOut[0] = pyra_rules::FACE_VRT_INDS[side][0];
1164 0 : indsOut[1] = pyra_rules::FACE_VRT_INDS[side][3];
1165 0 : indsOut[2] = pyra_rules::FACE_VRT_INDS[side][2];
1166 0 : indsOut[3] = pyra_rules::FACE_VRT_INDS[side][1];
1167 : }
1168 0 : }
1169 :
1170 0 : int Pyramid::
1171 : get_edge_index_from_vertices( const size_t vi0,
1172 : const size_t vi1) const
1173 : {
1174 0 : return pyra_rules::EDGE_FROM_VRTS[vi0][vi1];
1175 : }
1176 :
1177 0 : int Pyramid::
1178 : get_face_edge_index(const size_t faceInd,
1179 : const size_t faceEdgeInd) const
1180 : {
1181 0 : if(pyra_rules::FACE_EDGE_INDS[faceInd][3] == -1)
1182 0 : return pyra_rules::FACE_EDGE_INDS[faceInd][2 - faceEdgeInd];
1183 : else
1184 0 : return pyra_rules::FACE_EDGE_INDS[faceInd][3 - faceEdgeInd];
1185 : }
1186 :
1187 0 : std::pair<GridBaseObjectId, int> Pyramid::
1188 : get_opposing_object(Vertex* vrt) const
1189 : {
1190 : using namespace pyra_rules;
1191 0 : for(int i = 0; i < pyra_rules::NUM_VERTICES; ++i){
1192 0 : if(vrt == m_vertices[i]){
1193 0 : return make_pair(static_cast<GridBaseObjectId>(OPPOSED_OBJECT[i][0]),
1194 0 : OPPOSED_OBJECT[i][1]);
1195 : }
1196 : }
1197 0 : UG_THROW("Specified vertex is not part of this element.");
1198 : }
1199 :
1200 0 : bool Pyramid::collapse_edge(std::vector<Volume*>& vNewVolumesOut,
1201 : int edgeIndex, Vertex* newVertex,
1202 : std::vector<Vertex*>* pvSubstituteVertices)
1203 : {
1204 : using namespace pyra_rules;
1205 :
1206 : int elemInds[MAX_NUM_COLLAPSE_INDS_OUT];
1207 0 : int elemIndsSize = CollapseEdge(
1208 : elemInds,
1209 0 : EDGE_VRT_INDS[edgeIndex][0],
1210 0 : EDGE_VRT_INDS[edgeIndex][1]);
1211 :
1212 0 : if(elemIndsSize > 0){
1213 : Vertex** vrts;
1214 0 : if(pvSubstituteVertices)
1215 : vrts = &pvSubstituteVertices->front();
1216 : else
1217 0 : vrts = m_vertices;
1218 :
1219 0 : CreateVolumesFromElementIndexList(
1220 : vNewVolumesOut,
1221 : elemInds,
1222 : elemIndsSize,
1223 : vrts);
1224 0 : return !vNewVolumesOut.empty();
1225 : }
1226 : else{
1227 : vNewVolumesOut.clear();
1228 0 : return false;
1229 : }
1230 : }
1231 :
1232 0 : bool Pyramid::refine(std::vector<Volume*>& vNewVolumesOut,
1233 : Vertex** ppNewVertexOut,
1234 : Vertex** newEdgeVertices,
1235 : Vertex** newFaceVertices,
1236 : Vertex* newVolumeVertex,
1237 : const Vertex& prototypeVertex,
1238 : Vertex** pSubstituteVertices,
1239 : vector3* corners,
1240 : bool* isSnapPoint)
1241 : {
1242 : // handle substitute vertices.
1243 : Vertex** vrts;
1244 0 : if(pSubstituteVertices)
1245 : vrts = pSubstituteVertices;
1246 : else
1247 0 : vrts = m_vertices;
1248 :
1249 0 : return Refine<PyramidClass>(vNewVolumesOut, ppNewVertexOut,
1250 : newEdgeVertices, newFaceVertices,
1251 : newVolumeVertex, prototypeVertex,
1252 : vrts, pyra_rules::Refine, corners,
1253 0 : isSnapPoint);
1254 : }
1255 :
1256 0 : bool Pyramid::is_regular_ref_rule(int edgeMarks) const
1257 : {
1258 0 : return pyra_rules::IsRegularRefRule(edgeMarks);
1259 : }
1260 :
1261 0 : void Pyramid::get_flipped_orientation(VolumeDescriptor& vdOut) const
1262 : {
1263 : // in order to flip a pyramids orientation, we have to invert the order
1264 : // of the base-vertices
1265 : vdOut.set_num_vertices(5);
1266 0 : vdOut.set_vertex(0, vertex(3));
1267 0 : vdOut.set_vertex(1, vertex(2));
1268 0 : vdOut.set_vertex(2, vertex(1));
1269 0 : vdOut.set_vertex(3, vertex(0));
1270 0 : vdOut.set_vertex(4, vertex(4));
1271 0 : }
1272 :
1273 :
1274 : ////////////////////////////////////////////////////////////////////////
1275 : // OctahedronDescriptor
1276 0 : OctahedronDescriptor::OctahedronDescriptor(const OctahedronDescriptor& td)
1277 : {
1278 0 : m_vertex[0] = td.vertex(0);
1279 0 : m_vertex[1] = td.vertex(1);
1280 0 : m_vertex[2] = td.vertex(2);
1281 0 : m_vertex[3] = td.vertex(3);
1282 0 : m_vertex[4] = td.vertex(4);
1283 0 : m_vertex[5] = td.vertex(5);
1284 0 : }
1285 :
1286 0 : OctahedronDescriptor::OctahedronDescriptor(const VolumeVertices& vv)
1287 : {
1288 : assert((vv.num_vertices() == 6) && "Bad number of vertices in volume-descriptor. Should be 6.");
1289 0 : m_vertex[0] = vv.vertex(0);
1290 0 : m_vertex[1] = vv.vertex(1);
1291 0 : m_vertex[2] = vv.vertex(2);
1292 0 : m_vertex[3] = vv.vertex(3);
1293 0 : m_vertex[4] = vv.vertex(4);
1294 0 : m_vertex[5] = vv.vertex(5);
1295 0 : }
1296 :
1297 0 : OctahedronDescriptor::OctahedronDescriptor(Vertex* v1, Vertex* v2, Vertex* v3, Vertex* v4, Vertex* v5, Vertex* v6)
1298 : {
1299 0 : m_vertex[0] = v1;
1300 0 : m_vertex[1] = v2;
1301 0 : m_vertex[2] = v3;
1302 0 : m_vertex[3] = v4;
1303 0 : m_vertex[4] = v5;
1304 0 : m_vertex[5] = v6;
1305 0 : }
1306 :
1307 : ////////////////////////////////////////////////////////////////////////
1308 : // Octahedron
1309 0 : Octahedron::Octahedron(const OctahedronDescriptor& td)
1310 : {
1311 0 : m_vertices[0] = td.vertex(0);
1312 0 : m_vertices[1] = td.vertex(1);
1313 0 : m_vertices[2] = td.vertex(2);
1314 0 : m_vertices[3] = td.vertex(3);
1315 0 : m_vertices[4] = td.vertex(4);
1316 0 : m_vertices[5] = td.vertex(5);
1317 0 : }
1318 :
1319 0 : Octahedron::Octahedron(Vertex* v1, Vertex* v2, Vertex* v3, Vertex* v4, Vertex* v5, Vertex* v6)
1320 : {
1321 0 : m_vertices[0] = v1;
1322 0 : m_vertices[1] = v2;
1323 0 : m_vertices[2] = v3;
1324 0 : m_vertices[3] = v4;
1325 0 : m_vertices[4] = v5;
1326 0 : m_vertices[5] = v6;
1327 0 : }
1328 :
1329 0 : EdgeDescriptor Octahedron::edge_desc(int index) const
1330 : {
1331 0 : EdgeDescriptor ed;
1332 0 : edge_desc(index, ed);
1333 0 : return ed;
1334 : }
1335 :
1336 0 : void Octahedron::edge_desc(int index, EdgeDescriptor& edOut) const
1337 : {
1338 : using namespace oct_rules;
1339 : assert(index >= 0 && index <= NUM_EDGES);
1340 0 : edOut.set_vertices(m_vertices[EDGE_VRT_INDS[index][0]],
1341 0 : m_vertices[EDGE_VRT_INDS[index][1]]);
1342 0 : }
1343 :
1344 0 : uint Octahedron::num_edges() const
1345 : {
1346 0 : return 12;
1347 : }
1348 :
1349 0 : FaceDescriptor Octahedron::face_desc(int index) const
1350 : {
1351 0 : FaceDescriptor fd;
1352 0 : face_desc(index, fd);
1353 0 : return fd;
1354 : }
1355 :
1356 0 : void Octahedron::face_desc(int index, FaceDescriptor& fdOut) const
1357 : {
1358 : using namespace oct_rules;
1359 : assert(index >= 0 && index < NUM_FACES);
1360 :
1361 : fdOut.set_num_vertices(3);
1362 0 : fdOut.set_vertex(0, m_vertices[FACE_VRT_INDS[index][0]]);
1363 0 : fdOut.set_vertex(1, m_vertices[FACE_VRT_INDS[index][2]]);
1364 0 : fdOut.set_vertex(2, m_vertices[FACE_VRT_INDS[index][1]]);
1365 0 : }
1366 :
1367 0 : uint Octahedron::num_faces() const
1368 : {
1369 0 : return 8;
1370 : }
1371 :
1372 0 : Edge* Octahedron::create_edge(int index)
1373 : {
1374 : using namespace oct_rules;
1375 : assert(index >= 0 && index < NUM_EDGES);
1376 0 : const int* e = EDGE_VRT_INDS[index];
1377 0 : return new RegularEdge(m_vertices[e[0]], m_vertices[e[1]]);
1378 : }
1379 :
1380 0 : Face* Octahedron::create_face(int index)
1381 : {
1382 : using namespace oct_rules;
1383 : assert(index >= 0 && index < NUM_FACES);
1384 :
1385 0 : const int* f = FACE_VRT_INDS[index];
1386 0 : return new Triangle(m_vertices[f[0]], m_vertices[f[2]], m_vertices[f[1]]);
1387 : }
1388 :
1389 0 : void Octahedron::
1390 : get_vertex_indices_of_edge (size_t& ind1Out,
1391 : size_t& ind2Out,
1392 : size_t edgeInd) const
1393 : {
1394 : assert(edgeInd >= 0 && edgeInd < oct_rules::NUM_EDGES);
1395 0 : ind1Out = oct_rules::EDGE_VRT_INDS[edgeInd][0];
1396 0 : ind2Out = oct_rules::EDGE_VRT_INDS[edgeInd][1];
1397 0 : }
1398 :
1399 0 : void Octahedron::
1400 : get_vertex_indices_of_face (std::vector<size_t>& indsOut,
1401 : size_t side) const
1402 : {
1403 : assert(side >= 0 && side < oct_rules::NUM_FACES);
1404 :
1405 0 : indsOut.resize(3);
1406 0 : indsOut[0] = oct_rules::FACE_VRT_INDS[side][0];
1407 0 : indsOut[1] = oct_rules::FACE_VRT_INDS[side][2];
1408 0 : indsOut[2] = oct_rules::FACE_VRT_INDS[side][1];
1409 0 : }
1410 :
1411 0 : int Octahedron::
1412 : get_edge_index_from_vertices( const size_t vi0,
1413 : const size_t vi1) const
1414 : {
1415 0 : return oct_rules::EDGE_FROM_VRTS[vi0][vi1];
1416 : }
1417 :
1418 0 : int Octahedron::
1419 : get_face_edge_index(const size_t faceInd,
1420 : const size_t faceEdgeInd) const
1421 : {
1422 0 : return oct_rules::FACE_EDGE_INDS[faceInd][2 - faceEdgeInd];
1423 : }
1424 :
1425 0 : std::pair<GridBaseObjectId, int> Octahedron::
1426 : get_opposing_object(Vertex* vrt) const
1427 : {
1428 : using namespace oct_rules;
1429 0 : for(int i = 0; i < oct_rules::NUM_VERTICES; ++i){
1430 0 : if(vrt == m_vertices[i]){
1431 0 : return make_pair(static_cast<GridBaseObjectId>(OPPOSED_OBJECT[i][0]),
1432 0 : OPPOSED_OBJECT[i][1]);
1433 : }
1434 : }
1435 :
1436 0 : UG_THROW("Specified vertex is not part of this element.");
1437 : }
1438 :
1439 0 : bool Octahedron::collapse_edge(std::vector<Volume*>& vNewVolumesOut,
1440 : int edgeIndex, Vertex* newVertex,
1441 : std::vector<Vertex*>* pvSubstituteVertices)
1442 : {
1443 : // NOT YET SUPPORTED!
1444 : //TODO: implement octahedron::collapse_edge
1445 : vNewVolumesOut.clear();
1446 : UG_LOG("edge-collapse for octahedrons not yet implemented... sorry\n");
1447 0 : return false;
1448 : }
1449 :
1450 :
1451 0 : bool Octahedron::refine(std::vector<Volume*>& vNewVolumesOut,
1452 : Vertex** ppNewVertexOut,
1453 : Vertex** newEdgeVertices,
1454 : Vertex** newFaceVertices,
1455 : Vertex* newVolumeVertex,
1456 : const Vertex& prototypeVertex,
1457 : Vertex** pSubstituteVertices,
1458 : vector3* corners,
1459 : bool* isSnapPoint)
1460 : {
1461 : // handle substitute vertices.
1462 : Vertex** vrts;
1463 0 : if(pSubstituteVertices)
1464 : vrts = pSubstituteVertices;
1465 : else
1466 0 : vrts = m_vertices;
1467 :
1468 0 : return Refine<OctahedronClass>(vNewVolumesOut, ppNewVertexOut,
1469 : newEdgeVertices, newFaceVertices,
1470 : newVolumeVertex, prototypeVertex,
1471 : vrts, oct_rules::Refine, corners,
1472 0 : isSnapPoint);
1473 : }
1474 :
1475 0 : bool Octahedron::is_regular_ref_rule(int edgeMarks) const
1476 : {
1477 0 : return oct_rules::IsRegularRefRule(edgeMarks);
1478 : }
1479 :
1480 0 : void Octahedron::get_flipped_orientation(VolumeDescriptor& vdOut) const
1481 : {
1482 : // in order to flip a pyramids orientation, we have to invert the order
1483 : // of the base-vertices
1484 : vdOut.set_num_vertices(6);
1485 0 : vdOut.set_vertex(0, vertex(0));
1486 0 : vdOut.set_vertex(1, vertex(4));
1487 0 : vdOut.set_vertex(2, vertex(3));
1488 0 : vdOut.set_vertex(3, vertex(2));
1489 0 : vdOut.set_vertex(4, vertex(1));
1490 0 : vdOut.set_vertex(5, vertex(5));
1491 0 : }
1492 :
1493 : }// end of namespace
|
