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_2d.h"
36 : //#include "../algorithms/geom_obj_util/geom_obj_util.h"
37 :
38 : using namespace std;
39 :
40 : namespace ug
41 : {
42 : ////////////////////////////////////////////////////////////////////////
43 : ////////////////////////////////////////////////////////////////////////
44 : // TOOLS
45 : /// helpful if a local vertex-order is required
46 : /**
47 : * cornersOut and cornersIn both have to be of size numCorners.
48 : * After termination cornersOut will contain the vertices of
49 : * cornersIn, starting from firstCorner, taking vertices modulo numCorners.
50 : * If cornersOut == cornersIn, the method will fail! This is ok since
51 : * the method is used locally and has been created for a special case.
52 : */
53 : static inline
54 : bool ReorderCornersCCW(Vertex** cornersOut, Vertex** const cornersIn,
55 : int numCorners, int firstCorner)
56 : {
57 0 : cornersOut[0] = cornersIn[firstCorner];
58 0 : for(int i = 1; i < numCorners; ++i)
59 0 : cornersOut[i] = cornersIn[(firstCorner + i) % numCorners];
60 : return true;
61 : }
62 :
63 :
64 : ////////////////////////////////////////////////////////////////////////
65 : ////////////////////////////////////////////////////////////////////////
66 : // FACES
67 :
68 : ////////////////////////////////////////////////////////////////////////
69 : // TriangleDescriptor
70 0 : TriangleDescriptor::TriangleDescriptor(const TriangleDescriptor& td)
71 : {
72 0 : m_vertex[0] = td.vertex(0);
73 0 : m_vertex[1] = td.vertex(1);
74 0 : m_vertex[2] = td.vertex(2);
75 0 : }
76 :
77 1 : TriangleDescriptor::TriangleDescriptor(Vertex* v1, Vertex* v2, Vertex* v3)
78 : {
79 1 : m_vertex[0] = v1;
80 1 : m_vertex[1] = v2;
81 1 : m_vertex[2] = v3;
82 1 : }
83 :
84 : ////////////////////////////////////////////////////////////////////////
85 : // CustomTriangle
86 : template <class ConcreteTriangleType, class BaseClass, class RefTriType, class RefQuadType>
87 1 : CustomTriangle<ConcreteTriangleType, BaseClass, RefTriType, RefQuadType>::
88 1 : CustomTriangle(const TriangleDescriptor& td)
89 : {
90 1 : m_vertices[0] = td.vertex(0);
91 1 : m_vertices[1] = td.vertex(1);
92 1 : m_vertices[2] = td.vertex(2);
93 1 : }
94 :
95 : template <class ConcreteTriangleType, class BaseClass, class RefTriType, class RefQuadType>
96 0 : CustomTriangle<ConcreteTriangleType, BaseClass, RefTriType, RefQuadType>::
97 0 : CustomTriangle(Vertex* v1, Vertex* v2, Vertex* v3)
98 : {
99 0 : m_vertices[0] = v1;
100 0 : m_vertices[1] = v2;
101 0 : m_vertices[2] = v3;
102 0 : }
103 :
104 : template <class ConcreteTriangleType, class BaseClass, class RefTriType, class RefQuadType>
105 : std::pair<GridBaseObjectId, int>
106 0 : CustomTriangle<ConcreteTriangleType, BaseClass, RefTriType, RefQuadType>::
107 : get_opposing_object(Vertex* vrt) const
108 : {
109 0 : for(int i = 0; i < 3; ++i){
110 0 : if(vrt == m_vertices[i]){
111 0 : return make_pair(EDGE, (i + 1) % 3);
112 : }
113 : }
114 :
115 0 : UG_THROW("The given vertex is not contained in the given face.");
116 : }
117 :
118 : template <class ConcreteTriangleType, class BaseClass, class RefTriType, class RefQuadType>
119 : bool
120 0 : CustomTriangle<ConcreteTriangleType, BaseClass, RefTriType, RefQuadType>::
121 : refine(std::vector<Face*>& vNewFacesOut,
122 : Vertex** newFaceVertexOut,
123 : Vertex** newEdgeVertices,
124 : Vertex* newFaceVertex,
125 : Vertex** pSubstituteVertices,
126 : int snapPointIndex)
127 : {
128 : //TODO: complete triangle refine
129 :
130 0 : *newFaceVertexOut = newFaceVertex;
131 : vNewFacesOut.clear();
132 :
133 : // handle substitute vertices.
134 : Vertex** vrts;
135 0 : if(pSubstituteVertices)
136 : vrts = pSubstituteVertices;
137 : else
138 0 : vrts = m_vertices;
139 :
140 : // get the number of new vertices.
141 : uint numNewVrts = 0;
142 0 : for(uint i = 0; i < 3; ++i)
143 : {
144 0 : if(newEdgeVertices[i] != NULL)
145 0 : ++numNewVrts;
146 : }
147 :
148 : // if newFaceVertex is specified, then create three sub-triangles and
149 : // refine each. If not then refine the triangle depending
150 0 : if(newFaceVertex)
151 : {
152 0 : if(numNewVrts > 0){
153 : assert(!"Problem in CustomTriangle::refine: refine with newFaceVertex and newEdgeVertices is not yet supported.");
154 : return false;
155 : }
156 :
157 : // create three new triangles
158 0 : vNewFacesOut.push_back(new RefTriType(vrts[0], vrts[1], newFaceVertex));
159 0 : vNewFacesOut.push_back(new RefTriType(vrts[1], vrts[2], newFaceVertex));
160 0 : vNewFacesOut.push_back(new RefTriType(vrts[2], vrts[0], newFaceVertex));
161 0 : return true;
162 : }
163 : else
164 : {
165 0 : switch(numNewVrts)
166 : {
167 0 : case 0: // this may happen when the triangle belongs to a prism being anisotropically refined
168 : // and the volume on the other side is not being refined
169 0 : vNewFacesOut.push_back(new RefTriType(vrts[0], vrts[1], vrts[2]));
170 0 : return true;
171 :
172 : case 1:
173 : {
174 : // get the index of the edge that will be refined
175 : int iNew = -1;;
176 0 : for(int i = 0; i < 3; ++i){
177 0 : if(newEdgeVertices[i]){
178 : iNew = i;
179 : break;
180 : }
181 : }
182 :
183 : // the corners. The first corner is the corner on the opposite side of iNew.
184 : // Other follow in ccw order
185 : int iCorner[3];
186 0 : iCorner[0] = (iNew + 2) % 3;
187 0 : iCorner[1] = (iCorner[0] + 1) % 3;
188 0 : iCorner[2] = (iCorner[1] + 1) % 3;
189 :
190 : // create the new triangles.
191 0 : vNewFacesOut.push_back(new RefTriType(vrts[iCorner[0]], vrts[iCorner[1]],
192 0 : newEdgeVertices[iNew]));
193 0 : vNewFacesOut.push_back(new RefTriType(vrts[iCorner[0]], newEdgeVertices[iNew],
194 0 : vrts[iCorner[2]]));
195 :
196 : return true;
197 : }
198 :
199 : case 2:
200 : {
201 : // get the index of the edge that won't be refined
202 : int iFree = -1;
203 0 : for(int i = 0; i < 3; ++i){
204 0 : if(!newEdgeVertices[i]){
205 : iFree = i;
206 : break;
207 : }
208 : }
209 :
210 : // the refined edges
211 : int iNew[2];
212 0 : iNew[0] = (iFree + 1) % 3;
213 0 : iNew[1] = (iFree + 2) % 3;
214 :
215 : // the corners
216 : int iCorner[3];
217 : iCorner[0] = iFree;
218 : iCorner[1] = (iFree + 1) % 3;
219 : iCorner[2] = (iFree + 2) % 3;
220 :
221 : // create the faces
222 0 : vNewFacesOut.push_back(new RefTriType(newEdgeVertices[iNew[0]],
223 0 : vrts[iCorner[2]],
224 0 : newEdgeVertices[iNew[1]]));
225 0 : vNewFacesOut.push_back(new RefQuadType(vrts[iCorner[0]], vrts[iCorner[1]],
226 : newEdgeVertices[iNew[0]], newEdgeVertices[iNew[1]]));
227 : return true;
228 : }
229 :
230 0 : case 3:
231 : {
232 : // perform regular refine.
233 0 : vNewFacesOut.push_back(new RefTriType(vrts[0], newEdgeVertices[0], newEdgeVertices[2]));
234 0 : vNewFacesOut.push_back(new RefTriType(vrts[1], newEdgeVertices[1], newEdgeVertices[0]));
235 0 : vNewFacesOut.push_back(new RefTriType(vrts[2], newEdgeVertices[2], newEdgeVertices[1]));
236 0 : vNewFacesOut.push_back(new RefTriType(newEdgeVertices[0], newEdgeVertices[1], newEdgeVertices[2]));
237 0 : return true;
238 : }
239 :
240 : }
241 : }
242 :
243 : return false;
244 : }
245 :
246 : template <class ConcreteTriangleType, class BaseClass, class RefTriType, class RefQuadType>
247 : bool
248 0 : CustomTriangle<ConcreteTriangleType, BaseClass, RefTriType, RefQuadType>::
249 : is_regular_ref_rule(int edgeMarks) const
250 : {
251 0 : return edgeMarks == 7;
252 : }
253 :
254 : template <class ConcreteTriangleType, class BaseClass, class RefTriType, class RefQuadType>
255 : bool
256 0 : CustomTriangle<ConcreteTriangleType, BaseClass, RefTriType, RefQuadType>::
257 : collapse_edge(std::vector<Face*>& vNewFacesOut,
258 : int edgeIndex, Vertex* newVertex,
259 : Vertex** pSubstituteVertices)
260 : {
261 : // if an edge of the triangle is collapsed, nothing remains
262 : vNewFacesOut.clear();
263 0 : return true;
264 : }
265 :
266 : template <class ConcreteTriangleType, class BaseClass, class RefTriType, class RefQuadType>
267 : bool
268 0 : CustomTriangle<ConcreteTriangleType, BaseClass, RefTriType, RefQuadType>::
269 : collapse_edges(std::vector<Face*>& vNewFacesOut,
270 : std::vector<Vertex*>& vNewEdgeVertices,
271 : Vertex** pSubstituteVertices)
272 : {
273 0 : if(vNewEdgeVertices.size() > BaseClass::num_edges())
274 : {
275 : assert(!"WARNING in Triangle::collapse_edges(...): bad number of newEdgeVertices.");
276 : LOG("WARNING in Triangle::collapse_edges(...): bad number of newEdgeVertices.");
277 0 : return false;
278 : }
279 :
280 : // check if there is a valid entry in vNewEdgeVertices
281 : bool bGotOne = false;
282 0 : for(uint i = 0; i < vNewEdgeVertices.size(); ++i)
283 : {
284 0 : if(vNewEdgeVertices[i] != NULL)
285 : {
286 : bGotOne = true;
287 : break;
288 : }
289 : }
290 :
291 0 : if(!bGotOne)
292 : {
293 : assert(!"WARNING in Triangle::collapse:edges(...): no new vertex was specified.");
294 : LOG("WARNING in Triangle::collapse:edges(...): no new vertex was specified.");
295 0 : return false;
296 : }
297 :
298 : // if an edge of the triangle is collapsed, nothing remains
299 : vNewFacesOut.clear();
300 : return true;
301 : }
302 :
303 : // BEGIN Depreciated
304 : template <class ConcreteTriangleType, class BaseClass, class RefTriType, class RefQuadType>
305 : void
306 0 : CustomTriangle<ConcreteTriangleType, BaseClass, RefTriType, RefQuadType>::
307 : create_faces_by_edge_split(int splitEdgeIndex,
308 : Vertex* newVertex,
309 : std::vector<Face*>& vNewFacesOut,
310 : Vertex** pSubstituteVertices)
311 : {
312 : assert(((splitEdgeIndex >= 0) && (splitEdgeIndex < 3)) && "ERROR in Triangle::create_faces_by_edge_split(...): bad edge index!");
313 :
314 : // if pvSubstituteVertices is supplied, we will use the vertices in
315 : // pvSubstituteVertices instead the ones of 'this'.
316 : // If not, we will redirect the pointer to a local local vector,
317 : // that holds the vertices of 'this'
318 : Vertex** vrts;
319 0 : if(pSubstituteVertices)
320 : vrts = pSubstituteVertices;
321 : else
322 0 : vrts = m_vertices;
323 :
324 : // we have to find the indices ind0, ind1, ind2, where
325 : // ind0 is the index of the vertex on e before newVertex,
326 : // ind1 is the index of the vertex on e after newVertex
327 : // and ind2 is the index of the vertex not located on e.
328 :
329 : int ind0 = splitEdgeIndex;
330 0 : int ind1 = (ind0 + 1) % 3;
331 0 : int ind2 = (ind1 + 1) % 3;
332 :
333 0 : vNewFacesOut.push_back(new Triangle(vrts[ind0], newVertex, vrts[ind2]));
334 0 : vNewFacesOut.push_back(new Triangle(newVertex, vrts[ind1], vrts[ind2]));
335 0 : }
336 :
337 :
338 : ////////////////////////////////////////////////////////////////////////
339 : ////////////////////////////////////////////////////////////////////////
340 : // QUADRILATERALS
341 :
342 : ////////////////////////////////////////////////////////////////////////
343 : // QuadrilateralDescriptor
344 0 : QuadrilateralDescriptor::QuadrilateralDescriptor(const QuadrilateralDescriptor& qd)
345 : {
346 0 : m_vertex[0] = qd.vertex(0);
347 0 : m_vertex[1] = qd.vertex(1);
348 0 : m_vertex[2] = qd.vertex(2);
349 0 : m_vertex[3] = qd.vertex(3);
350 0 : }
351 :
352 0 : QuadrilateralDescriptor::QuadrilateralDescriptor(Vertex* v1, Vertex* v2, Vertex* v3, Vertex* v4)
353 : {
354 0 : m_vertex[0] = v1;
355 0 : m_vertex[1] = v2;
356 0 : m_vertex[2] = v3;
357 0 : m_vertex[3] = v4;
358 0 : }
359 :
360 : ////////////////////////////////////////////////////////////////////////
361 : // Quad
362 :
363 : template <class ConcreteQuadrilateralType, class BaseClass, class RefTriType, class RefQuadType>
364 0 : CustomQuadrilateral<ConcreteQuadrilateralType, BaseClass, RefTriType, RefQuadType>::
365 0 : CustomQuadrilateral(const QuadrilateralDescriptor& qd)
366 : {
367 0 : m_vertices[0] = qd.vertex(0);
368 0 : m_vertices[1] = qd.vertex(1);
369 0 : m_vertices[2] = qd.vertex(2);
370 0 : m_vertices[3] = qd.vertex(3);
371 0 : }
372 :
373 : template <class ConcreteQuadrilateralType, class BaseClass, class RefTriType, class RefQuadType>
374 0 : CustomQuadrilateral<ConcreteQuadrilateralType, BaseClass, RefTriType, RefQuadType>::
375 0 : CustomQuadrilateral(Vertex* v1, Vertex* v2, Vertex* v3, Vertex* v4)
376 : {
377 0 : m_vertices[0] = v1;
378 0 : m_vertices[1] = v2;
379 0 : m_vertices[2] = v3;
380 0 : m_vertices[3] = v4;
381 0 : }
382 :
383 : template <class ConcreteQuadrilateralType, class BaseClass, class RefTriType, class RefQuadType>
384 : bool
385 0 : CustomQuadrilateral<ConcreteQuadrilateralType, BaseClass, RefTriType, RefQuadType>::
386 : get_opposing_side(EdgeVertices* e, EdgeDescriptor& edOut) const
387 : {
388 0 : int localInd = Face::get_local_side_index(e);
389 0 : if(localInd == -1){
390 : return false;
391 : }
392 :
393 0 : edge_desc((localInd + 2) % 4, edOut);
394 0 : return true;
395 : }
396 :
397 : template <class ConcreteQuadrilateralType, class BaseClass, class RefTriType, class RefQuadType>
398 : std::pair<GridBaseObjectId, int>
399 0 : CustomQuadrilateral<ConcreteQuadrilateralType, BaseClass, RefTriType, RefQuadType>::
400 : get_opposing_object(Vertex* vrt) const
401 : {
402 0 : for(int i = 0; i < 4; ++i){
403 0 : if(vrt == m_vertices[i]){
404 0 : return make_pair(VERTEX, (i + 2) % 4);
405 : }
406 : }
407 :
408 0 : UG_THROW("The given vertex is not contained in the given face.");
409 : }
410 :
411 : template <class ConcreteQuadrilateralType, class BaseClass, class RefTriType, class RefQuadType>
412 : void
413 0 : CustomQuadrilateral<ConcreteQuadrilateralType, BaseClass, RefTriType, RefQuadType>::
414 : create_faces_by_edge_split(int splitEdgeIndex,
415 : Vertex* newVertex,
416 : std::vector<Face*>& vNewFacesOut,
417 : Vertex** pSubstituteVertices)
418 : {
419 : assert(((splitEdgeIndex >= 0) && (splitEdgeIndex < 4)) && "ERROR in Quadrilateral::create_faces_by_edge_split(...): bad edge index!");
420 :
421 : // if pvSubstituteVertices is supplied, we will use the vertices in
422 : // pvSubstituteVertices instead the ones of 'this'.
423 : // If not, we will redirect the pointer to a local local vector,
424 : // that holds the vertices of 'this'
425 : Vertex** vrts;
426 0 : if(pSubstituteVertices)
427 : vrts = pSubstituteVertices;
428 : else
429 0 : vrts = m_vertices;
430 :
431 : // we have to find the indices ind0, ind1, ind2, where
432 : // ind0 is the index of the vertex on e before newVertex,
433 : // ind1 is the index of the vertex on e after newVertex
434 : // and ind2 and ind3 are the indices of the vertices not located on e.
435 : int ind0 = splitEdgeIndex; //edge-index equals the index of its first vertex.
436 0 : int ind1 = (ind0 + 1) % 4;
437 0 : int ind2 = (ind0 + 2) % 4;
438 0 : int ind3 = (ind0 + 3) % 4;
439 :
440 : TriangleDescriptor td;
441 :
442 : // edge-split generates 3 triangles
443 0 : vNewFacesOut.push_back(new RefTriType(vrts[ind0], newVertex, vrts[ind3]));
444 0 : vNewFacesOut.push_back(new RefTriType(newVertex, vrts[ind1], vrts[ind2]));
445 0 : vNewFacesOut.push_back(new RefTriType(vrts[ind3], newVertex, vrts[ind2]));
446 :
447 : // we're done.
448 0 : }
449 :
450 : ////////////////////////////////////////////////////////////////////////
451 : // Quadrilateral::refine
452 : template <class ConcreteQuadrilateralType, class BaseClass, class RefTriType, class RefQuadType>
453 : bool
454 0 : CustomQuadrilateral<ConcreteQuadrilateralType, BaseClass, RefTriType, RefQuadType>::
455 : refine(std::vector<Face*>& vNewFacesOut,
456 : Vertex** newFaceVertexOut,
457 : Vertex** edgeVrts,
458 : Vertex* newFaceVertex,
459 : Vertex** pSubstituteVertices,
460 : int snapPointIndex)
461 : {
462 : //TODO: complete quad refine
463 0 : *newFaceVertexOut = newFaceVertex;
464 : vNewFacesOut.clear();
465 :
466 : // handle substitute vertices.
467 : Vertex** vrts;
468 0 : if(pSubstituteVertices)
469 : vrts = pSubstituteVertices;
470 : else
471 0 : vrts = m_vertices;
472 :
473 : // check which edges have to be refined and perform the required operations.
474 : // get the number of new vertices.
475 : uint numNewVrts = 0;
476 0 : for(uint i = 0; i < 4; ++i)
477 : {
478 0 : if(edgeVrts[i] != NULL)
479 0 : ++numNewVrts;
480 : }
481 :
482 0 : switch(numNewVrts)
483 : {
484 0 : case 0:
485 : // refine with mid point if it exists
486 0 : if (newFaceVertex)
487 : {
488 : // create four new triangles
489 0 : vNewFacesOut.push_back(new RefTriType(vrts[0], vrts[1], newFaceVertex));
490 0 : vNewFacesOut.push_back(new RefTriType(vrts[1], vrts[2], newFaceVertex));
491 0 : vNewFacesOut.push_back(new RefTriType(vrts[2], vrts[3], newFaceVertex));
492 0 : vNewFacesOut.push_back(new RefTriType(vrts[3], vrts[0], newFaceVertex));
493 0 : return true;
494 : }
495 :
496 : // in case the mid point does not exists, we need a simple copy
497 : // This may happen when the quad belongs to a hexahedron being anisotropically refined
498 : // and the volume on the other side is not being refined.
499 0 : vNewFacesOut.push_back(new RefQuadType(vrts[0], vrts[1], vrts[2], vrts[3]));
500 0 : return true;
501 :
502 0 : case 1:
503 : {
504 0 : if(newFaceVertex){
505 : assert(!"Problem in CustomQuadrilateral::refine: refine with newFaceVertex and one newEdgeVertex is not yet supported.");
506 : return false;
507 : }
508 :
509 : int iNew = -1;
510 0 : for(int i = 0; i < 4; ++i){
511 0 : if(edgeVrts[i]){
512 : iNew = i;
513 : break;
514 : }
515 : }
516 :
517 0 : if(snapPointIndex != -1 && (snapPointIndex == iNew || snapPointIndex == (iNew + 1) % 4)){
518 : UG_LOG("WARNING: Invalid snap-point distribution detected. Ignoring snap-points for this element.\n");
519 : snapPointIndex = -1;
520 : }
521 :
522 : // the corners in a local ordering relative to iNew. Keeping ccw order.
523 : Vertex* corner[4];
524 0 : const int rot = (iNew + 3) % 4;
525 : ReorderCornersCCW(corner, vrts, 4, rot);
526 :
527 : // create the new elements
528 0 : if(snapPointIndex == -1){
529 0 : vNewFacesOut.push_back(new RefTriType(corner[0], corner[1], edgeVrts[iNew]));
530 0 : vNewFacesOut.push_back(new RefTriType(corner[0], edgeVrts[iNew], corner[3]));
531 0 : vNewFacesOut.push_back(new RefTriType(corner[3], edgeVrts[iNew], corner[2]));
532 : }
533 : else{
534 0 : snapPointIndex = (snapPointIndex + 4 - rot) % 4;
535 0 : if(snapPointIndex == 0){
536 0 : vNewFacesOut.push_back(new RefTriType(corner[0], corner[1], edgeVrts[iNew]));
537 0 : vNewFacesOut.push_back(new RefQuadType(corner[0], edgeVrts[iNew], corner[2], corner[3]));
538 : }
539 0 : else if(snapPointIndex == 3){
540 0 : vNewFacesOut.push_back(new RefQuadType(corner[0], corner[1], edgeVrts[iNew], corner[3]));
541 0 : vNewFacesOut.push_back(new RefTriType(corner[3], edgeVrts[iNew], corner[2]));
542 : }
543 : else{
544 0 : UG_THROW("Unexpected snap-point index: " << snapPointIndex << ". This is an implementation error!");
545 : }
546 : }
547 :
548 : return true;
549 : }
550 :
551 0 : case 2:
552 : {
553 0 : if(newFaceVertex){
554 : assert(!"Problem in CustomQuadrilateral::refine: refine with newFaceVertex and two newEdgeVertices is not yet supported.");
555 : return false;
556 : }
557 :
558 : // there are two cases (refined edges are adjacent or opposite to each other)
559 : int iNew[2];
560 : int counter = 0;
561 0 : for(int i = 0; i < 4; ++i){
562 0 : if(edgeVrts[i]){
563 0 : iNew[counter] = i;
564 0 : ++counter;
565 : }
566 : }
567 :
568 : // corners will be filled later on
569 : Vertex* corner[4];
570 :
571 : // check which case applies
572 0 : if(iNew[1] - iNew[0] == 2){
573 : // edges are opposite to each other
574 : // the corners in a local ordering relative to iNew[0]. Keeping ccw order.
575 0 : ReorderCornersCCW(corner, vrts, 4, (iNew[0] + 3) % 4);
576 :
577 : // create new faces
578 0 : vNewFacesOut.push_back(new RefQuadType(corner[0], corner[1],
579 0 : edgeVrts[iNew[0]], edgeVrts[iNew[1]]));
580 :
581 0 : vNewFacesOut.push_back(new RefQuadType(edgeVrts[iNew[1]], edgeVrts[iNew[0]],
582 : corner[2], corner[3]));
583 : }
584 : else{
585 : // edges are adjacent
586 : // we want that (iNew[0] + 1) % 4 == iNew[1]
587 0 : if((iNew[0] + 1) % 4 != iNew[1])
588 : swap(iNew[0], iNew[1]);
589 :
590 : // the corners in a local ordering relative to iNew[0]. Keeping ccw order.
591 0 : const int rot = (iNew[0] + 3) % 4;
592 : ReorderCornersCCW(corner, vrts, 4, rot);
593 :
594 0 : if(snapPointIndex != -1 && ((snapPointIndex + 4 - rot) % 4) != 0){
595 : snapPointIndex = -1;
596 : UG_LOG("WARNING: Invalid snap-point distribution detected. Ignoring snap-points for this element.\n");
597 : }
598 :
599 :
600 : // create new faces
601 : if(snapPointIndex == -1){
602 0 : vNewFacesOut.push_back(new RefTriType(corner[0], corner[1], edgeVrts[iNew[0]]));
603 0 : vNewFacesOut.push_back(new RefTriType(edgeVrts[iNew[0]], corner[2], edgeVrts[iNew[1]]));
604 0 : vNewFacesOut.push_back(new RefTriType(corner[3], corner[0], edgeVrts[iNew[1]]));
605 0 : vNewFacesOut.push_back(new RefTriType(corner[0], edgeVrts[iNew[0]], edgeVrts[iNew[1]]));
606 : }
607 : else{
608 0 : vNewFacesOut.push_back(new RefTriType(corner[0], corner[1], edgeVrts[iNew[0]]));
609 0 : vNewFacesOut.push_back(new RefTriType(corner[3], corner[0], edgeVrts[iNew[1]]));
610 0 : vNewFacesOut.push_back(new RefQuadType(corner[0], edgeVrts[iNew[0]], corner[2], edgeVrts[iNew[1]]));
611 : }
612 : }
613 :
614 : return true;
615 : }
616 :
617 0 : case 3:
618 : {
619 0 : if(newFaceVertex){
620 : assert(!"Problem in CustomQuadrilateral::refine: refine with newFaceVertex and three newEdgeVertices is not yet supported.");
621 : return false;
622 : }
623 :
624 : int iFree = -1;
625 0 : for(int i = 0; i < 4; ++i){
626 0 : if(!edgeVrts[i]){
627 : iFree = i;
628 : break;
629 : }
630 : }
631 :
632 : // the vertices on the edges:
633 : Vertex* nvrts[3];
634 0 : nvrts[0] = edgeVrts[(iFree + 1) % 4];
635 0 : nvrts[1] = edgeVrts[(iFree + 2) % 4];
636 0 : nvrts[2] = edgeVrts[(iFree + 3) % 4];
637 :
638 : // the corners in a local ordering relative to iNew. Keeping ccw order.
639 : Vertex* corner[4];
640 : ReorderCornersCCW(corner, vrts, 4, (iFree + 1) % 4);
641 :
642 : // create the faces
643 0 : vNewFacesOut.push_back(new RefQuadType(corner[0], nvrts[0], nvrts[2], corner[3]));
644 0 : vNewFacesOut.push_back(new RefTriType(corner[1], nvrts[1], nvrts[0]));
645 0 : vNewFacesOut.push_back(new RefTriType(corner[2], nvrts[2], nvrts[1]));
646 0 : vNewFacesOut.push_back(new RefTriType(nvrts[0], nvrts[1], nvrts[2]));
647 :
648 0 : return true;
649 : }
650 :
651 0 : case 4:
652 : {
653 : // we'll create 4 new quads. create a new center if required.
654 0 : if(!newFaceVertex)
655 0 : newFaceVertex = new RegularVertex;
656 :
657 0 : *newFaceVertexOut = newFaceVertex;
658 :
659 0 : vNewFacesOut.push_back(new RefQuadType(vrts[0], edgeVrts[0], newFaceVertex, edgeVrts[3]));
660 0 : vNewFacesOut.push_back(new RefQuadType(vrts[1], edgeVrts[1], newFaceVertex, edgeVrts[0]));
661 0 : vNewFacesOut.push_back(new RefQuadType(vrts[2], edgeVrts[2], newFaceVertex, edgeVrts[1]));
662 0 : vNewFacesOut.push_back(new RefQuadType(vrts[3], edgeVrts[3], newFaceVertex, edgeVrts[2]));
663 0 : return true;
664 : }
665 : }
666 :
667 : return false;
668 : }
669 :
670 : template <class ConcreteQuadrilateralType, class BaseClass, class RefTriType, class RefQuadType>
671 : bool
672 0 : CustomQuadrilateral<ConcreteQuadrilateralType, BaseClass, RefTriType, RefQuadType>::
673 : is_regular_ref_rule(int edgeMarks) const
674 : {
675 : static const int allEdges = 15; // 1111
676 : static const int hEdges = 5; // 0101
677 : static const int vEdges = 10; // 1010
678 :
679 0 : return (edgeMarks == allEdges)
680 0 : || (edgeMarks == hEdges)
681 0 : || (edgeMarks == vEdges);
682 : }
683 :
684 :
685 : template <class ConcreteQuadrilateralType, class BaseClass, class RefTriType, class RefQuadType>
686 : bool
687 0 : CustomQuadrilateral<ConcreteQuadrilateralType, BaseClass, RefTriType, RefQuadType>::
688 : collapse_edge(std::vector<Face*>& vNewFacesOut,
689 : int edgeIndex, Vertex* newVertex,
690 : Vertex** pSubstituteVertices)
691 : {
692 : // if an edge of the triangle is collapsed, nothing remains
693 : vNewFacesOut.clear();
694 :
695 : // handle substitute vertices.
696 : Vertex** vrts;
697 0 : if(pSubstituteVertices)
698 : vrts = pSubstituteVertices;
699 : else
700 0 : vrts = m_vertices;
701 :
702 : // the collapsed edge connects vertices at ind0 and ind1.
703 : int ind0 = edgeIndex; //edge-index equals the index of its first vertex.
704 0 : int ind1 = (ind0 + 1) % 4;
705 0 : int ind2 = (ind1 + 1) % 4;
706 0 : int ind3 = (ind2 + 1) % 4;
707 :
708 : // ind0 and ind1 will be replaced by newVertex.
709 0 : vNewFacesOut.push_back(new RefTriType(newVertex, vrts[ind2], vrts[ind3]));
710 0 : return true;
711 : }
712 :
713 : template <class ConcreteQuadrilateralType, class BaseClass, class RefTriType, class RefQuadType>
714 : bool
715 0 : CustomQuadrilateral<ConcreteQuadrilateralType, BaseClass, RefTriType, RefQuadType>::
716 : collapse_edges(std::vector<Face*>& vNewFacesOut,
717 : std::vector<Vertex*>& vNewEdgeVertices,
718 : Vertex** pSubstituteVertices)
719 : {
720 0 : if(vNewEdgeVertices.size() > BaseClass::num_edges())
721 : {
722 : assert(!"WARNING in Quadrilateral::collapse_edges(...): bad number of newEdgeVertices.");
723 : LOG("WARNING in Quadrilateral::collapse_edges(...): bad number of newEdgeVertices.");
724 0 : return false;
725 : }
726 :
727 : vNewFacesOut.clear();
728 :
729 : // check if there is a valid entry in vNewEdgeVertices
730 : int collapseIndex = -1;
731 : int numCollapses = 0;
732 0 : for(uint i = 0; i < 4; ++i)
733 : {
734 0 : if(i < vNewEdgeVertices.size())
735 : {
736 0 : if(vNewEdgeVertices[i] != NULL)
737 : {
738 0 : ++numCollapses;
739 0 : collapseIndex = i;
740 : }
741 : }
742 : }
743 :
744 : // if more than 1 edge is collapsed, nothing remains.
745 0 : if(numCollapses == 0)
746 : {
747 : assert(!"WARNING in Quadrilateral::collapse:edges(...): no new vertex was specified.");
748 : LOG("WARNING in Quadrilateral::collapse:edges(...): no new vertex was specified.");
749 0 : return false;
750 : }
751 0 : else if(numCollapses == 1)
752 : {
753 : // call collapse_edge with the edges index.
754 0 : collapse_edge(vNewFacesOut, collapseIndex, vNewEdgeVertices[collapseIndex], pSubstituteVertices);
755 : }
756 :
757 : return true;
758 : }
759 :
760 : // explicit instantiation
761 : template class CustomTriangle<Triangle, Face, Triangle, Quadrilateral>;
762 : template class CustomTriangle<ConstrainedTriangle, ConstrainedFace,
763 : ConstrainedTriangle, ConstrainedQuadrilateral>;
764 : template class CustomTriangle<ConstrainingTriangle, ConstrainingFace,
765 : ConstrainingTriangle, ConstrainingQuadrilateral>;
766 :
767 : template class CustomQuadrilateral<Quadrilateral, Face, Triangle, Quadrilateral>;
768 : template class CustomQuadrilateral<ConstrainedQuadrilateral, ConstrainedFace,
769 : ConstrainedTriangle, ConstrainedQuadrilateral>;
770 : template class CustomQuadrilateral<ConstrainingQuadrilateral, ConstrainingFace,
771 : ConstrainingTriangle, ConstrainingQuadrilateral>;
772 :
773 :
774 : ////////////////////////////////////////////////////////////////////////////////
775 : // CONSTRAINING FACE
776 : template <> size_t
777 0 : ConstrainingFace::
778 : num_constrained<Vertex>() const
779 : {
780 0 : return num_constrained_vertices();
781 : }
782 :
783 : template <> size_t
784 0 : ConstrainingFace::
785 : num_constrained<Edge>() const
786 : {
787 0 : return num_constrained_edges();
788 : }
789 :
790 : template <> size_t
791 0 : ConstrainingFace::
792 : num_constrained<Face>() const
793 : {
794 0 : return num_constrained_faces();
795 : }
796 :
797 :
798 : template <> Vertex*
799 0 : ConstrainingFace::
800 : constrained<Vertex>(size_t ind) const
801 : {
802 0 : return constrained_vertex(ind);
803 : }
804 :
805 : template <> Edge*
806 0 : ConstrainingFace::
807 : constrained<Edge>(size_t ind) const
808 : {
809 0 : return constrained_edge(ind);
810 : }
811 :
812 : template <> Face*
813 0 : ConstrainingFace::
814 : constrained<Face>(size_t ind) const
815 : {
816 0 : return constrained_face(ind);
817 : }
818 :
819 : }// end of namespace
|
