Line data Source code
1 : /*
2 : * DiamondsEstablish3D.cpp
3 : *
4 : * Created on: 08.12.2025
5 : * Author: Markus Knodel
6 : */
7 :
8 : #include <lib_grid/algorithms/extrusion/DiamondsEstablish3D.h>
9 : #include <typeinfo>
10 :
11 : namespace ug
12 : {
13 :
14 : namespace arte
15 : {
16 :
17 : namespace diamonds
18 : {
19 :
20 0 : DiamondsEstablish3D::DiamondsEstablish3D( Grid & grid,
21 : SubsetHandler & sh,
22 : DiamondsEstablish3D::VecVolManifVrtxCombi const & vecVolManifVrtxC
23 0 : )
24 : :
25 0 : m_grid(grid),
26 0 : m_sh(sh),
27 : m_aaPos(Grid::VertexAttachmentAccessor<APosition>()),
28 0 : m_vecVolManifVrtxCombiToShrink4Diams(vecVolManifVrtxC),
29 : m_vecVolElmFac5(VecVolumeElementFaceQuintuplet()),
30 0 : m_diamondsOnlyPreform(false),
31 : m_vecElems2BQuenched(VecElems2BQuenched()),
32 : m_disappearingVols(std::vector<Volume*>()),
33 : m_disappearingFacs(std::vector<Face*>()),
34 : m_vecElemGroupVrtx2BQuenched(VecElemGroupVrtx2BQnchd4D()),
35 : m_attElmGrpVrtx2BQnchd(AttElemGrpVrtx2BQuenchd()),
36 : m_attAccsElmGrpVrtx2BQnchd(Grid::VertexAttachmentAccessor<AttElemGrpVrtx2BQuenchd>()),
37 : m_attMarkTwinFace(ABool()),
38 : m_attAccsFacIsTwinFac(Grid::FaceAttachmentAccessor<ABool>()),
39 : m_attMarkShiftFace(ABool()),
40 : m_attAccsFacIsShiftFac(Grid::FaceAttachmentAccessor<ABool>()),
41 : m_attMarkShiftTriangleFace(ABool()),
42 : m_attAccsFacIsShiftTriangleFac(Grid::FaceAttachmentAccessor<ABool>()),
43 : m_attMarkShiftQuadriliteralFace(ABool()),
44 : m_attAccsFacIsShiftQuadriliteralFac(Grid::FaceAttachmentAccessor<ABool>()),
45 : m_attMarkVolGetsShrinked(ABool()),
46 : m_attAccsVolGetsShrinked(Grid::VolumeAttachmentAccessor<ABool>()),
47 : m_attVrtVec(AttVrtVec()),
48 : m_attAccsVrtVecVol(Grid::VolumeAttachmentAccessor<AttVrtVec>()),
49 : // m_attVrtVecFace(AttVrtVec()),
50 : // m_attAccsVrtVecFace(Grid::VolumeAttachmentAccessor<AttVrtVec>()),
51 : m_attMarkVrtxIsCenterVrtx(ABool()),
52 : m_attAccsVrtxIsCenterVrtx(Grid::VertexAttachmentAccessor<ABool>()),
53 : m_attMarkVrtxIsShiftVrtx(ABool()),
54 : m_attAccsVrtxIsShiftVrtx(Grid::VertexAttachmentAccessor<ABool>()),
55 : m_attMarkEdgeIsShiftEdge(ABool()),
56 : m_attAccsEdgeIsShiftEdge(Grid::EdgeAttachmentAccessor<ABool>()),
57 : m_attInfoVecSudosTouchingVrtx(AttVecInt()),
58 : m_attAccsInfoVecSudosTouchingVrtx(Grid::VertexAttachmentAccessor<AttVecInt>()),
59 : m_attVrtxIsMidPtOfShiftVrtx(ABool()),
60 : m_attAccsVrtxIsMidPtOfShiftVrtcs(Grid::VertexAttachmentAccessor<ABool>()),
61 : m_vecCombiShiftVrtxMidVrtx(VecCombiShiftVrtxMidVrtx()),
62 : m_attMidPtVrtxOfShiftVrtx(AVertex()),
63 : m_attAccsMidPtVrtxOfShiftVrtx(Grid::VertexAttachmentAccessor<AVertex>()),
64 : m_attCenterVrtxOfShiftVrtx(AVertex()),
65 : m_attAccsCenterVrtxOfShiftVrtx(Grid::VertexAttachmentAccessor<AVertex>()),
66 : m_vecCombiNewVolsTwoCross(VecCombiNewVolsProps()),
67 : m_vecCombiNewVolsThreeCross(VecCombiNewVolsProps()),
68 : m_sudosTable(VecIndxVec()),
69 : m_vecCombiCntrVrtxSudo(VecCombiCntrVrtxSudo()),
70 : m_attEdgeCanBeRemoved(ABool()),
71 : m_attAccsEdgeCanBeRemoved(Grid::EdgeAttachmentAccessor<ABool>()),
72 : m_attNewSudoOfVrtx(AInt()),
73 : m_attAccsNewSudoOfVrtx(Grid::VertexAttachmentAccessor<AInt>()),
74 : m_faces2BDeletedAtLastStep(std::vector<Face*>()),
75 : m_edges2BDeletedAtLastStep(std::vector<Edge*>()),
76 : m_attCentersCutPts(AttVec3Vec()),
77 : m_attAccsCentersCutPts(Grid::EdgeAttachmentAccessor<AttVec3Vec>()),
78 : // m_edgesCut4Diams(std::vector<Edge*>()),
79 : m_attDiamCentrVrtx(AVertex()),
80 0 : m_attAccsDiamCentrVrtx(Grid::EdgeAttachmentAccessor<AVertex>())
81 : {
82 : // // Notloesung, nicht in die erste Initialisierung vor geschweifter Klammer, da copy constructor privat
83 0 : m_sel = Selector();
84 0 : }
85 :
86 : //bool DiamondsEstablish3D::createTheDiamonds()
87 : //{
88 : // IndexType sudosVols = m_sh.num_subsets();
89 : //
90 : // IndexType sudosEdges = sudosVols + 1;
91 : //
92 : // IndexType sudosFaces = sudosVols + 2;
93 : //
94 : // for( auto & vmvcd : m_vecVolManifVrtxCombiToShrink4Diams )
95 : // {
96 : // Volume* vol;
97 : // vmvcd.spuckVol(vol);
98 : //
99 : // m_sh.assign_subset(vol, sudosVols);
100 : //
101 : // IndexType numLowdimElmsFnd = vmvcd.computeTheLowdimElm(m_grid);
102 : //
103 : // if( numLowdimElmsFnd != 1 )
104 : // {
105 : // UG_LOG("number of lowdim elems found strange " << numLowdimElmsFnd << std::endl);
106 : // UG_THROW("number of lowdim elems found strange " << numLowdimElmsFnd << std::endl);
107 : // }
108 : //
109 : // Edge* edge;
110 : // vmvcd.spuckLowdimElem( edge );
111 : //
112 : // if( edge == nullptr )
113 : // {
114 : // UG_LOG("Edge nicht gefunden " << std::endl);
115 : // UG_THROW("Edge nicht gefunden " << std::endl);
116 : // }
117 : // m_sh.assign_subset( edge, sudosEdges);
118 : //
119 : // Face * fac;
120 : // vmvcd.spuckManif(fac);
121 : // m_sh.assign_subset( fac, sudosFaces );
122 : // }
123 : //
124 : // UG_LOG("Established diamonds" << std::endl);
125 : //
126 : //
127 : // return true;
128 : //}
129 :
130 0 : DiamondsEstablish3D::~DiamondsEstablish3D()
131 : {
132 : // Auto-generated destructor stub
133 0 : }
134 :
135 0 : bool DiamondsEstablish3D::initialize()
136 : {
137 0 : if(!m_grid.has_vertex_attachment(aPosition) )
138 : {
139 : UG_LOG("Error in ExpandFractures Arte 3D: Missing position attachment");
140 0 : return false;
141 : }
142 :
143 0 : m_aaPos = Grid::VertexAttachmentAccessor<APosition>(m_grid, aPosition);
144 :
145 0 : return true;
146 : }
147 :
148 : /////////////////////////////////////////////////////////////////
149 :
150 :
151 0 : bool DiamondsEstablish3D::createTheDiamonds( bool diamondsOnlyPreform )
152 : {
153 0 : m_diamondsOnlyPreform = diamondsOnlyPreform;
154 :
155 0 : if( ! initialize())
156 : {
157 : UG_LOG("initialization diamonds did not work " << std::endl);
158 0 : return false;
159 : }
160 :
161 0 : if( ! setSelector() )
162 : return false;
163 :
164 : UG_LOG("selektiert" << std::endl);
165 :
166 0 : if( ! figureOutTheEdges() )
167 : {
168 : UG_LOG("Edges not found " << std::endl);
169 0 : return false;
170 : }
171 :
172 : UG_LOG("edges figured out " << std::endl);
173 :
174 0 : if( ! findRegions2BShrinked())
175 : {
176 : UG_LOG("Regions to be shrinked not found " << std::endl);
177 0 : return false;
178 : }
179 :
180 : UG_LOG("regions to be shrinked found " << std::endl);
181 :
182 0 : if( ! establishElems2BeQuenched() )
183 : {
184 : UG_LOG("quenching elements not establishable " << std::endl);
185 0 : return false;
186 : }
187 :
188 : UG_LOG("quenching elements established " << std::endl);
189 :
190 0 : if( ! sortElems2BQuenched())
191 : {
192 : UG_LOG("sorting quench impossible" << std::endl);
193 0 : return false;
194 : }
195 :
196 : UG_LOG("sorted elements finished " << std::endl);
197 :
198 0 : if( ! attachMarkers())
199 : {
200 : UG_LOG("Markers do not want D" << std::endl);
201 0 : return false;
202 : }
203 :
204 : UG_LOG("markers attached" << std::endl);
205 :
206 0 : if( ! assignBasicAtts())
207 : {
208 : UG_LOG("unassignale basic att" << std::endl);
209 0 : return false;
210 : }
211 :
212 : UG_LOG("basic attcs assigned" << std::endl);
213 :
214 0 : if( ! trafoCollectedInfo2Attachments())
215 : {
216 : UG_LOG("untrafoable infos" << std::endl);
217 0 : return false;
218 : }
219 :
220 : UG_LOG("trafo collected infos " << std::endl);
221 :
222 : // disable selection inheritance to avoid infinite recursion.
223 0 : m_sel.enable_selection_inheritance(false);
224 :
225 0 : if( ! createConditionForNewVrtcs())
226 : {
227 : UG_LOG("conditions for new vertices do not work " << std::endl);
228 0 : return false;
229 : }
230 :
231 : UG_LOG("conditions for new vertices created" << std::endl);
232 :
233 0 : if( ! distributeInfosForShrinkingVols())
234 : {
235 : UG_LOG("info distribution did not work " << std::endl);
236 0 : return false;
237 : }
238 :
239 : UG_LOG("info distributed" << std::endl);
240 :
241 0 : if( ! determineShiftFaces())
242 : {
243 : UG_LOG("shift faces not determinable " << std::endl);
244 0 : return false;
245 : }
246 :
247 : UG_LOG("detect removable edges " << std::endl);
248 :
249 0 : if( ! detectRemovableEdges())
250 : {
251 : UG_LOG("removable edges not detected " << std::endl);
252 0 : return false;
253 : }
254 :
255 : UG_LOG("shrink volumes " << std::endl);
256 :
257 0 : if( ! shrinkVolumes())
258 : {
259 : UG_LOG("shrinken schief gegangen " << std::endl);
260 0 : return false;
261 : }
262 :
263 : UG_LOG("volumes shrinked " << std::endl);
264 :
265 : // if( diamondsOnlyPreform )
266 : // {
267 : // UG_LOG("diamonds only preform, stopping diamond postprocessing at this stage " << std::endl);
268 : // }
269 :
270 0 : if( ! assignSudos2DiamsPreform())
271 : {
272 : UG_LOG("postprocessing diam vols not working" << std::endl);
273 0 : return false;
274 : }
275 :
276 : UG_LOG("sudos to preforms assigned " << std::endl);
277 :
278 0 : if( ! m_diamondsOnlyPreform )
279 : {
280 0 : if( ! splitDiamsPreform2Diams())
281 : {
282 : UG_LOG("did not manage to split diams" << std::endl);
283 0 : return false;
284 : }
285 : }
286 : else
287 : {
288 : UG_LOG("No split of preform diamonds, remaining at preform " << std::endl);
289 : }
290 :
291 :
292 : UG_LOG("diam vols postprocessed" << std::endl);
293 :
294 0 : if( ! detachMarkers())
295 : {
296 : UG_LOG("Markers do not get detouched D" << std::endl);
297 0 : return false;
298 : }
299 :
300 : UG_LOG("markers detached" << std::endl);
301 :
302 : UG_LOG("noch nicht soweit: Established diamonds" << std::endl);
303 :
304 :
305 :
306 0 : return true;
307 : }
308 :
309 : /////////////////////////////////////////////////////////////////
310 :
311 :
312 0 : bool DiamondsEstablish3D::setSelector()
313 : {
314 :
315 0 : m_sel.assign_grid(m_grid);
316 :
317 0 : m_sel.enable_autoselection(false);
318 0 : m_sel.enable_selection_inheritance(true); //required for select and mark, disabled later
319 0 : m_sel.enable_strict_inheritance(false);
320 :
321 :
322 0 : return true;
323 : }
324 :
325 : /////////////////////////////////////////////////////////////////
326 :
327 0 : bool DiamondsEstablish3D::figureOutTheEdges()
328 : {
329 :
330 : // compute the edges connecting the old and shift vertex
331 :
332 0 : for( auto & vmvcd : m_vecVolManifVrtxCombiToShrink4Diams )
333 : {
334 : // IndexType numLowdimElmsFnd = vmvcd.checkIntegrity(m_grid);
335 : // if( numLowdimElmsFnd != 1 )
336 0 : if( ! vmvcd.checkIntegrity(m_grid) )
337 : {
338 : UG_LOG("number of lowdim elems found strange " << std::endl);
339 : return false;
340 : }
341 : }
342 :
343 : UG_LOG("Edges found" << std::endl);
344 :
345 0 : return true;
346 : }
347 :
348 : ///////////////////////////////////////////////////////////////////////////////////////
349 :
350 :
351 : ///////////////////////////////////////////////////////////////////////////////
352 :
353 : ///////////////////////////////////////////////////////////////////////////////////////
354 :
355 0 : bool DiamondsEstablish3D::findRegions2BShrinked()
356 : {
357 : UG_LOG("want to find regions to be shrinked" << std::endl);
358 :
359 :
360 0 : VecVolManifVrtxCombi vecVolManifVrtxCopy = m_vecVolManifVrtxCombiToShrink4Diams;
361 :
362 : int d_out = 0;
363 :
364 0 : for( typename VecVolManifVrtxCombi::iterator itVMVOuter = vecVolManifVrtxCopy.begin();
365 0 : itVMVOuter != vecVolManifVrtxCopy.end(); )
366 : {
367 :
368 : // UG_LOG("out round " << d_out << std::endl);
369 :
370 : bool partnerFound = false;
371 :
372 0 : VolManifVrtxCombi outer = *itVMVOuter;
373 :
374 : VrtxPair oldAndShiftVrtxOuter;
375 : outer.spuckOldAndShiftVrtx( oldAndShiftVrtxOuter );
376 :
377 : Vertex * oldVrtxOuter = oldAndShiftVrtxOuter.first;
378 : Vertex * shiftVrtxOuter = oldAndShiftVrtxOuter.second;
379 :
380 : Face * faceOuter;
381 : outer.spuckManif(faceOuter);
382 :
383 : IndexType sudoOuter = outer.spuckSudo();
384 :
385 : int d_in = 0;
386 :
387 0 : for( typename VecVolManifVrtxCombi::iterator itVMVInner = itVMVOuter + 1;
388 0 : itVMVInner != vecVolManifVrtxCopy.end(); )
389 : {
390 : // UG_LOG("start in " << d_in << std::endl);
391 :
392 0 : VolManifVrtxCombi inner = *itVMVInner;
393 :
394 : VrtxPair oldAndShiftVrtxInner;
395 : inner.spuckOldAndShiftVrtx( oldAndShiftVrtxInner );
396 :
397 : Vertex * oldVrtxInner = oldAndShiftVrtxInner.first;
398 :
399 0 : if( oldVrtxInner == oldVrtxOuter )
400 : {
401 : Vertex * shiftVrtxInner = oldAndShiftVrtxInner.second;
402 :
403 : // if connected by a face, twin to outer found, partner may exist or not
404 : // check if the face is identical, else is from the partner twin
405 :
406 : Face * faceInner;
407 : inner.spuckManif(faceInner);
408 :
409 0 : if( faceInner == faceOuter )
410 : {
411 : // twin found
412 :
413 : PairVolFacVrtxCmb prVolFacVrtxC( outer, inner );
414 : VolumeElementFaceQuintuplet vef5;
415 :
416 0 : trafoVolFacVrtxCombiPair2FullLowDimManifQuintuplet( prVolFacVrtxC, vef5 );
417 :
418 0 : if( ! vef5.checkIntegrity() )
419 : {
420 : UG_LOG("strange twin produced" << std::endl);
421 0 : return false;
422 : }
423 :
424 0 : m_vecVolElmFac5.push_back(vef5);
425 :
426 : itVMVInner = vecVolManifVrtxCopy.erase(itVMVInner);
427 :
428 : partnerFound = true;
429 :
430 : }
431 : else
432 : {
433 : itVMVInner++;
434 : }
435 :
436 : if( partnerFound )
437 : {
438 : break;
439 : }
440 : }
441 : else
442 : {
443 : itVMVInner++;
444 : }
445 : if( partnerFound )
446 : {
447 : break;
448 : }
449 :
450 :
451 : // UG_LOG("end in " << d_in << std::endl);
452 : d_in++;
453 : }
454 :
455 : // UG_LOG("first round finished " << std::endl);
456 0 : if( ! partnerFound )
457 : {
458 : UG_LOG("no partner found " << std::endl);
459 : return false;
460 : }
461 : else
462 : {
463 : itVMVOuter = vecVolManifVrtxCopy.erase(itVMVOuter);
464 : }
465 :
466 : // UG_LOG("end out " << d_out << std::endl);
467 : d_out++;
468 : }
469 :
470 : UG_LOG("end of search reached diams " << std::endl);
471 :
472 : // for( auto & v: m_vecVolElmFac5 )
473 : // {
474 : // IndexType sudoNum = m_sh.num_subsets();
475 : //
476 : // IndexType sudoVols = sudoNum;
477 : // IndexType sudoFacs = sudoNum+1;
478 : // IndexType sudoEdgs = sudoNum+2;
479 : // IndexType sudoVrtx = sudoNum+3;
480 : //
481 : // Volume * vol1;
482 : // Volume * vol2;
483 : // Face * fac;
484 : // Edge * edg1;
485 : // Edge * edg2;
486 : // Vertex * vrtC;
487 : // Vertex * vrtE1;
488 : // Vertex * vrtE2;
489 : //
490 : // std::pair<VolumeElementTwin,VolumeElementTwin> pvv;
491 : //
492 : // v.spuckPairFullLowDimTwin(pvv);
493 : //
494 : // pvv.first.spuckFullDimElem(vol1);
495 : // pvv.second.spuckFullDimElem(vol2);
496 : //
497 : // v.spuckManifElem( fac );
498 : //
499 : // pvv.first.spuckLowDimElem(edg1);
500 : // pvv.second.spuckLowDimElem(edg2);
501 : //
502 : // v.spuckCenterVertex(vrtC);
503 : //
504 : // VrtxPair vp;
505 : // v.spuckShiftVrtcs(vp);
506 : //
507 : // vrtE1 = vp.first;
508 : // vrtE2 = vp.second;
509 : //
510 : // m_sh.assign_subset(vol1, sudoVols);
511 : // m_sh.assign_subset(vol2, sudoVols);
512 : // m_sh.assign_subset(fac, sudoFacs);
513 : // m_sh.assign_subset(edg1, sudoEdgs);
514 : // m_sh.assign_subset(edg2, sudoEdgs);
515 : // m_sh.assign_subset(vrtC, sudoVrtx);
516 : // m_sh.assign_subset(vrtE1, sudoVrtx);
517 : // m_sh.assign_subset(vrtE2, sudoVrtx);
518 : //
519 : // }
520 :
521 :
522 : UG_LOG("found regions to be shrinked" << std::endl);
523 :
524 :
525 : return true;
526 0 : }
527 :
528 :
529 : ////////////////////////////////////////////////////////////////////////////
530 :
531 0 : bool DiamondsEstablish3D::trafoVolFacVrtxCombiPair2FullLowDimManifQuintuplet(
532 : PairVolFacVrtxCmb & prVolFacVrtxC, VolumeElementFaceQuintuplet & vef5 )
533 : {
534 :
535 : VolManifVrtxCombi & mvcOne = prVolFacVrtxC.first;
536 : VolManifVrtxCombi & mvcTwo = prVolFacVrtxC.second;
537 :
538 : IndexType sudoOne = mvcOne.spuckSudo();
539 : IndexType sudoTwo = mvcTwo.spuckSudo();
540 :
541 0 : if( sudoOne != sudoTwo )
542 : {
543 : UG_LOG("sudos differ " << std::endl);
544 0 : return false;
545 : }
546 :
547 : IndexType sudo = sudoOne;
548 :
549 : Face * faceOne;
550 : Face * faceTwo;
551 :
552 : mvcOne.spuckManif(faceOne);
553 : mvcTwo.spuckManif(faceTwo);
554 :
555 0 : if( faceOne != faceTwo )
556 : {
557 : UG_LOG("faces differ " << std::endl);
558 0 : return false;
559 : }
560 :
561 : Face * connectingFace = faceOne;
562 :
563 : VrtxPair oldAndShiftVrtxOne;
564 : mvcOne.spuckOldAndShiftVrtx( oldAndShiftVrtxOne );
565 :
566 : VrtxPair oldAndShiftVrtxTwo;
567 : mvcTwo.spuckOldAndShiftVrtx( oldAndShiftVrtxTwo );
568 :
569 : Vertex * oldVrtxOne = oldAndShiftVrtxOne.first;
570 : Vertex * oldVrtxTwo = oldAndShiftVrtxTwo.first;
571 :
572 0 : if( oldVrtxOne != oldVrtxTwo )
573 : {
574 : UG_LOG("center vertices not identical " << std::endl);
575 0 : return false;
576 : }
577 :
578 : Vertex * oldVrtx = oldVrtxOne;
579 :
580 : Vertex * shiftVrtxOne = oldAndShiftVrtxOne.second;
581 : Vertex * shiftVrtxTwo = oldAndShiftVrtxTwo.second;
582 :
583 0 : if( shiftVrtxOne == shiftVrtxTwo )
584 : {
585 : UG_LOG("shift vertices coincide but should not " << std::endl);
586 0 : return false;
587 : }
588 :
589 : Volume * volOne;
590 : Volume * volTwo;
591 :
592 : mvcOne.spuckFulldimElem( volOne );
593 : mvcTwo.spuckFulldimElem( volTwo );
594 :
595 0 : if( volOne == volTwo )
596 : {
597 : UG_LOG("volumes coincide but should not " << std::endl);
598 0 : return false;
599 : }
600 :
601 : Edge * edgeOne;
602 : Edge * edgeTwo;
603 :
604 : mvcOne.spuckLowDimElem( edgeOne );
605 : mvcTwo.spuckLowDimElem( edgeTwo );
606 0 : if( edgeOne == edgeTwo )
607 : {
608 : UG_LOG("edges coincide but should not " << std::endl);
609 0 : return false;
610 : }
611 :
612 : VolumeElementTwin volElTwinOne( volOne, edgeOne, sudo );
613 : VolumeElementTwin volElTwinTwo( volTwo, edgeTwo, sudo );
614 :
615 0 : if( ! volElTwinOne.checkIntegrity() || ! volElTwinTwo.checkIntegrity() )
616 : {
617 : UG_LOG("twins of vol edge not integer " << std::endl);
618 0 : return false;
619 : }
620 :
621 : std::pair<VolumeElementTwin,VolumeElementTwin> volElTwinPair( volElTwinOne, volElTwinTwo );
622 :
623 0 : vef5 = VolumeElementFaceQuintuplet( volElTwinPair, connectingFace );
624 :
625 0 : if( ! vef5.checkIntegrity() )
626 : {
627 : UG_LOG( "quitent not integer " << std::endl );
628 0 : return false;
629 : }
630 :
631 : return true;
632 : }
633 :
634 :
635 : ////////////////////////////////////////////////////////////////////////////
636 :
637 0 : bool DiamondsEstablish3D::establishElems2BeQuenched()
638 : {
639 0 : VecVolumeElementFaceQuintuplet vvef5 = m_vecVolElmFac5;
640 :
641 0 : for( typename VecVolumeElementFaceQuintuplet::iterator itOuter = vvef5.begin();
642 0 : itOuter != vvef5.end();
643 : )
644 : {
645 0 : VolumeElementFaceQuintuplet vef5Outer = *itOuter;
646 :
647 : VecVolumeElementFaceQuintuplet vfld5ThisEdgePr;
648 :
649 0 : vfld5ThisEdgePr.push_back(vef5Outer);
650 :
651 : // now collect all quintuplets belonging to the same edge pair
652 :
653 : EdgePair edgPrOuter;
654 :
655 : vef5Outer.spuckPairLowDimElem(edgPrOuter);
656 :
657 0 : for( typename VecVolumeElementFaceQuintuplet::iterator itInner = vvef5.begin() + 1;
658 0 : itInner != vvef5.end();
659 : )
660 : {
661 0 : VolumeElementFaceQuintuplet vef5Inner = *itInner;
662 :
663 : EdgePair edgPrInner;
664 : vef5Inner.spuckPairLowDimElem(edgPrInner);
665 :
666 : if( edgPrInner == edgPrOuter )
667 : {
668 0 : vfld5ThisEdgePr.push_back(vef5Inner);
669 :
670 : itInner = vvef5.erase(itInner);
671 : }
672 : else // if( edgPrInner != edgPrOuter )
673 : {
674 : EdgePair testEdges = edgPrInner;
675 : std::swap( testEdges.first, testEdges.second );
676 :
677 : if( testEdges == edgPrOuter )
678 : {
679 0 : if( ! vef5Inner.swapEntries() )
680 : {
681 : UG_LOG("swapping not worked " << std::endl);
682 : return false;
683 : }
684 :
685 0 : if( ! vef5Inner.checkIntegrity())
686 : {
687 : UG_LOG("not integer any more after sapping test" << std::endl);
688 : return false;
689 : }
690 :
691 : EdgePair testAgainEdgPr;
692 : vef5Inner.spuckPairLowDimElem(testAgainEdgPr);
693 :
694 : if( testAgainEdgPr == edgPrOuter )
695 : {
696 0 : vfld5ThisEdgePr.push_back(vef5Inner);
697 :
698 : itInner = vvef5.erase(itInner);
699 : }
700 : else
701 : {
702 : UG_LOG("should be correct but is not " << std::endl);
703 : return false;
704 : }
705 : }
706 : else
707 : {
708 : itInner++;
709 : }
710 : }
711 : }
712 :
713 : Elems2BQuenched elem2BQuenched( vfld5ThisEdgePr );
714 :
715 0 : if( ! elem2BQuenched.checkIntegrity())
716 : {
717 : UG_LOG("an elem to be quenched not integer" << std::endl);
718 : return false;
719 : }
720 :
721 0 : if( ! assignMidPointOfShiftVrtcs(elem2BQuenched))
722 : {
723 : UG_LOG("mid point not assignable" << std::endl);
724 : return false;
725 : }
726 :
727 0 : m_vecElems2BQuenched.push_back(elem2BQuenched);
728 :
729 : itOuter = vvef5.erase(itOuter);
730 0 : }
731 :
732 : // int d_q = 10;
733 :
734 : return true;
735 0 : }
736 :
737 : //////////////////////////////////////////////////////////////////////////////7
738 :
739 0 : void DiamondsEstablish3D::debugE2bQ(Elems2BQuenched & e2bq)
740 : {
741 : // for( Elems2BQuenched & e2bq : m_vecElems2BQuenched )
742 : {
743 0 : IndexType sudoNum = m_sh.num_subsets();
744 :
745 : IndexType sudoVols = sudoNum;
746 0 : IndexType sudoFacs = sudoNum+1;
747 0 : IndexType sudoEdgs = sudoNum+2;
748 0 : IndexType sudoVrtx = sudoNum+3;
749 :
750 : VecVolumeElementFaceQuintuplet vve5;
751 :
752 : EdgePair edgP;
753 : e2bq.spuckPairLowDimElem(edgP);
754 :
755 : e2bq.spuckVecFullLowDimManifQuintuplet(vve5);
756 :
757 : Vertex * centerVrtx;
758 : e2bq.spuckCenterVertex(centerVrtx);
759 :
760 : vector3 vertexLocation;
761 :
762 0 : if( centerVrtx != nullptr )
763 : vertexLocation = m_aaPos[centerVrtx];
764 : else
765 : {
766 : UG_LOG("center null " << std::endl);
767 : }
768 :
769 : number x = vertexLocation[0];
770 : number y = vertexLocation[1];
771 : number z = vertexLocation[2];
772 :
773 0 : number distSq = ( x - 0.5 )*( x - 0.5 ) + ( y - 0.5 )*( y - 0.5 ) + ( z - 0.5 )*( z - 0.5 );
774 :
775 0 : if( distSq < 0.02 )
776 : {
777 :
778 0 : for( VolumeElementFaceQuintuplet & v: vve5 )
779 : {
780 : Volume * vol1;
781 : Volume * vol2;
782 : Face * fac;
783 : Edge * edg1;
784 : Edge * edg2;
785 : Vertex * vrtC;
786 : Vertex * vrtE1;
787 : Vertex * vrtE2;
788 :
789 : std::pair<VolumeElementTwin,VolumeElementTwin> pvv;
790 :
791 : v.spuckPairFullLowDimTwin(pvv);
792 :
793 : pvv.first.spuckFullDimElem(vol1);
794 : pvv.second.spuckFullDimElem(vol2);
795 :
796 : v.spuckManifElem( fac );
797 :
798 : pvv.first.spuckLowDimElem(edg1);
799 : pvv.second.spuckLowDimElem(edg2);
800 :
801 : v.spuckCenterVertex(vrtC);
802 :
803 : VrtxPair vp;
804 : v.spuckShiftVrtcs(vp);
805 :
806 : vrtE1 = vp.first;
807 : vrtE2 = vp.second;
808 :
809 :
810 : Volume * newVol1;
811 : Volume * newVol2;
812 :
813 0 : newVol1 = *m_grid.create<Prism>(
814 0 : PrismDescriptor( vol1->vertex(0),
815 0 : vol1->vertex(1),
816 0 : vol1->vertex(2),
817 0 : vol1->vertex(3),
818 0 : vol1->vertex(4),
819 0 : vol1->vertex(5)
820 : )
821 : );
822 :
823 0 : newVol2 = *m_grid.create<Prism>(
824 0 : PrismDescriptor( vol2->vertex(0),
825 0 : vol2->vertex(1),
826 0 : vol2->vertex(2),
827 0 : vol2->vertex(3),
828 0 : vol2->vertex(4),
829 0 : vol2->vertex(5)
830 : )
831 : );
832 :
833 :
834 0 : m_sh.assign_subset(newVol1, sudoVols);
835 0 : m_sh.assign_subset(newVol2, sudoVols);
836 :
837 : Face * newFac;
838 :
839 0 : newFac = *m_grid.create<Triangle>(TriangleDescriptor( fac->vertex(0), fac->vertex(1), fac->vertex(2) ));
840 :
841 0 : m_sh.assign_subset(newFac, sudoFacs);
842 :
843 :
844 : // m_sh.assign_subset(vol1, sudoVols);
845 : // m_sh.assign_subset(vol2, sudoVols);
846 :
847 : // m_sh.assign_subset(fac, sudoFacs);
848 :
849 :
850 : // m_sh.assign_subset(fac, sudoFacs);
851 : // m_sh.assign_subset(edg1, sudoEdgs);
852 : // m_sh.assign_subset(edg2, sudoEdgs);
853 : // m_sh.assign_subset(vrtC, sudoVrtx);
854 : // m_sh.assign_subset(vrtE1, sudoVrtx);
855 : // m_sh.assign_subset(vrtE2, sudoVrtx);
856 :
857 : Edge * edgOne;
858 : Edge * edgTwo;
859 :
860 0 : edgOne = *m_grid.create<RegularEdge>(EdgeDescriptor( edgP.first->vertex(0), edgP.first->vertex(1) ));
861 0 : edgTwo = *m_grid.create<RegularEdge>(EdgeDescriptor( edgP.second->vertex(0), edgP.second->vertex(1) ));
862 :
863 0 : m_sh.assign_subset(edgOne,sudoEdgs);
864 0 : m_sh.assign_subset(edgTwo,sudoEdgs);
865 :
866 : Vertex * centerV;
867 :
868 : e2bq.spuckCenterVertex(centerV);
869 :
870 0 : Vertex * newCenter = *m_grid.create<RegularVertex>();
871 :
872 : m_aaPos[newCenter] = m_aaPos[centerV];
873 :
874 0 : m_sh.assign_subset(newCenter, sudoVrtx);
875 :
876 : }
877 :
878 :
879 :
880 :
881 : // if( d_q == 20 )
882 : // return true;
883 : //
884 : // d_q++;
885 : //
886 : }
887 0 : }
888 0 : }
889 :
890 : ////////////////////////////////////////////////////////////////////////////
891 :
892 0 : bool DiamondsEstablish3D::sortElems2BQuenched()
893 : {
894 0 : VecElems2BQuenched ve2bq = m_vecElems2BQuenched;
895 :
896 0 : for( typename VecElems2BQuenched::iterator itOuter = ve2bq.begin();
897 0 : itOuter != ve2bq.end();
898 : )
899 : {
900 : Elems2BQuenched e2bqOuter = *itOuter;
901 :
902 : Vertex * centerVrtxOuter;
903 : VecElems2BQuenched vecElems2Q4ThisVrtx;
904 : e2bqOuter.spuckOrigCenterVertex( centerVrtxOuter );
905 :
906 0 : vecElems2Q4ThisVrtx.push_back(e2bqOuter);
907 :
908 0 : for( typename VecElems2BQuenched::iterator itInner = ve2bq.begin() + 1;
909 0 : itInner != ve2bq.end();
910 : )
911 : {
912 : Elems2BQuenched e2bqInner = *itInner;
913 :
914 : Vertex * centerVrtxInner;
915 : e2bqInner.spuckOrigCenterVertex( centerVrtxInner );
916 :
917 0 : if( centerVrtxOuter == centerVrtxInner )
918 : {
919 0 : vecElems2Q4ThisVrtx.push_back(e2bqInner);
920 : itInner = ve2bq.erase(itInner);
921 : }
922 : else
923 : {
924 : itInner++;
925 : }
926 : }
927 :
928 : ElemGroupVrtx2BQuenched4Diams egv2b( vecElems2Q4ThisVrtx );
929 :
930 0 : if( ! egv2b.checkIntegrity() )
931 : {
932 : UG_LOG("elem group not integer for vertex" << std::endl);
933 : return false;
934 : }
935 :
936 0 : m_vecElemGroupVrtx2BQuenched.push_back( egv2b );
937 :
938 : itOuter = ve2bq.erase(itOuter);
939 :
940 0 : }
941 :
942 : return true;
943 0 : }
944 :
945 : ////////////////////////////////////////////////////////////////////////////
946 :
947 0 : bool DiamondsEstablish3D::attachMarkers()
948 : {
949 : // m_aAdjMarkerVFP = AttVertFracProp();
950 : //
951 : //// support::VertexFracturePropertiesVol<IndexType> vfp0; // false, 0 );
952 : // VertxFracPropts vfp0; // false, 0 );
953 : // // default value: no boundary fracture, no fractures crossing
954 : //
955 : // m_grid.attach_to_vertices_dv( m_aAdjMarkerVFP, vfp0 );
956 : // m_aaMarkVrtVFP = Grid::VertexAttachmentAccessor<AttVertFracProp> ( m_grid, m_aAdjMarkerVFP );
957 :
958 : m_attElmGrpVrtx2BQnchd = AttElemGrpVrtx2BQuenchd();
959 :
960 : ElemGroupVrtx2BQuenched4Diams egv2bQEmpty;
961 :
962 0 : m_grid.attach_to_vertices_dv( m_attElmGrpVrtx2BQnchd, egv2bQEmpty );
963 :
964 0 : m_attAccsElmGrpVrtx2BQnchd = Grid::VertexAttachmentAccessor<AttElemGrpVrtx2BQuenchd>( m_grid, m_attElmGrpVrtx2BQnchd);
965 :
966 : m_attMarkTwinFace = ABool(); // used to know if an face is twin face
967 :
968 0 : m_grid.attach_to_faces_dv( m_attMarkTwinFace, false );
969 :
970 0 : m_attAccsFacIsTwinFac = Grid::FaceAttachmentAccessor<ABool>( m_grid, m_attMarkTwinFace );
971 :
972 : m_attMarkShiftFace = ABool(); // used to know if an face is shift face
973 :
974 0 : m_grid.attach_to_faces_dv( m_attMarkShiftFace, false );
975 :
976 0 : m_attAccsFacIsShiftFac = Grid::FaceAttachmentAccessor<ABool>( m_grid, m_attMarkShiftFace );
977 :
978 : m_attMarkShiftTriangleFace = ABool(); // used to know if an face is shift face
979 :
980 0 : m_grid.attach_to_faces_dv(m_attMarkShiftTriangleFace, false);
981 :
982 0 : m_attAccsFacIsShiftTriangleFac = Grid::FaceAttachmentAccessor<ABool>(m_grid, m_attMarkShiftTriangleFace);
983 :
984 : m_attMarkShiftQuadriliteralFace = ABool(); // used to know if an face is shift face
985 :
986 0 : m_grid.attach_to_faces_dv(m_attMarkShiftQuadriliteralFace, false);
987 :
988 0 : m_attAccsFacIsShiftQuadriliteralFac = Grid::FaceAttachmentAccessor<ABool>(m_grid, m_attMarkShiftQuadriliteralFace);
989 :
990 : m_attMarkVolGetsShrinked = ABool();
991 :
992 0 : m_grid.attach_to_volumes_dv( m_attMarkVolGetsShrinked, false );
993 :
994 0 : m_attAccsVolGetsShrinked = Grid::VolumeAttachmentAccessor<ABool>( m_grid, m_attMarkVolGetsShrinked );
995 :
996 : m_attVrtVec = AttVrtVec();
997 :
998 0 : m_grid.attach_to_volumes_dv(m_attVrtVec, std::vector<Vertex*>());
999 :
1000 : // m_attVrtVecFace = AttVrtVec();
1001 : //
1002 : // m_grid.attach_to_faces(m_attVrtVecFace, std::vector<Vertex*>());
1003 : //
1004 : // m_attAccsVrtVecFace = Grid::VolumeAttachmentAccessor<AttVrtVec>(m_grid, m_attVrtVecFace);
1005 :
1006 0 : m_attAccsVrtVecVol = Grid::VolumeAttachmentAccessor<AttVrtVec>(m_grid, m_attVrtVec);
1007 :
1008 : m_attMarkVrtxIsCenterVrtx = ABool();
1009 :
1010 0 : m_grid.attach_to_vertices_dv(m_attMarkVrtxIsCenterVrtx,false);
1011 :
1012 0 : m_attAccsVrtxIsCenterVrtx = Grid::VertexAttachmentAccessor<ABool>( m_grid, m_attMarkVrtxIsCenterVrtx);
1013 :
1014 : m_attMarkVrtxIsShiftVrtx = ABool();
1015 :
1016 0 : m_grid.attach_to_vertices_dv(m_attMarkVrtxIsShiftVrtx,false);
1017 :
1018 0 : m_attAccsVrtxIsShiftVrtx = Grid::VertexAttachmentAccessor<ABool>( m_grid, m_attMarkVrtxIsShiftVrtx);
1019 :
1020 : m_attMarkEdgeIsShiftEdge = ABool();
1021 :
1022 0 : m_grid.attach_to_edges_dv(m_attMarkEdgeIsShiftEdge,false);
1023 :
1024 0 : m_attAccsEdgeIsShiftEdge = Grid::EdgeAttachmentAccessor<ABool>( m_grid, m_attMarkEdgeIsShiftEdge );
1025 :
1026 : m_attInfoVecSudosTouchingVrtx = AttVecInt();
1027 :
1028 0 : m_grid.attach_to_vertices_dv(m_attInfoVecSudosTouchingVrtx, std::vector<IndexType>());
1029 :
1030 0 : m_attAccsInfoVecSudosTouchingVrtx = Grid::VertexAttachmentAccessor<AttVecInt>( m_grid, m_attInfoVecSudosTouchingVrtx);
1031 :
1032 : m_attVrtxIsMidPtOfShiftVrtx = ABool();
1033 :
1034 0 : m_grid.attach_to_vertices_dv(m_attVrtxIsMidPtOfShiftVrtx,false);
1035 :
1036 0 : m_attAccsVrtxIsMidPtOfShiftVrtcs = Grid::VertexAttachmentAccessor<ABool>(m_grid, m_attVrtxIsMidPtOfShiftVrtx);
1037 :
1038 : m_attMidPtVrtxOfShiftVrtx = AVertex();
1039 :
1040 0 : m_grid.attach_to_vertices_dv(m_attMidPtVrtxOfShiftVrtx, nullptr);
1041 :
1042 0 : m_attAccsMidPtVrtxOfShiftVrtx = Grid::VertexAttachmentAccessor<AVertex>( m_grid, m_attMidPtVrtxOfShiftVrtx );
1043 :
1044 : m_attCenterVrtxOfShiftVrtx = AVertex();
1045 :
1046 0 : m_grid.attach_to_vertices_dv( m_attCenterVrtxOfShiftVrtx, nullptr);
1047 :
1048 0 : m_attAccsCenterVrtxOfShiftVrtx = Grid::VertexAttachmentAccessor<AVertex>( m_grid, m_attCenterVrtxOfShiftVrtx );
1049 :
1050 : m_attEdgeCanBeRemoved = ABool();
1051 :
1052 0 : m_grid.attach_to_edges_dv(m_attEdgeCanBeRemoved, false);
1053 :
1054 0 : m_attAccsEdgeCanBeRemoved = Grid::EdgeAttachmentAccessor<ABool>( m_grid, m_attEdgeCanBeRemoved);
1055 :
1056 : m_attNewSudoOfVrtx = AInt();
1057 :
1058 0 : m_grid.attach_to_vertices_dv( m_attNewSudoOfVrtx, -1 );
1059 :
1060 0 : m_attAccsNewSudoOfVrtx = Grid::VertexAttachmentAccessor<AInt>( m_grid, m_attNewSudoOfVrtx );
1061 :
1062 : m_attCentersCutPts = AttVec3Vec();
1063 :
1064 0 : m_grid.attach_to_edges_dv( m_attCentersCutPts, std::vector<vector3>() );
1065 :
1066 0 : m_attAccsCentersCutPts = Grid::EdgeAttachmentAccessor<AttVec3Vec>( m_grid, m_attCentersCutPts );
1067 :
1068 : m_attDiamCentrVrtx = AVertex();
1069 :
1070 0 : m_grid.attach_to_edges_dv( m_attDiamCentrVrtx, nullptr);
1071 :
1072 0 : m_attAccsDiamCentrVrtx = Grid::EdgeAttachmentAccessor<AVertex>( m_grid, m_attDiamCentrVrtx);
1073 :
1074 0 : return true;
1075 : }
1076 :
1077 : ////////////////////////////////////////////////////////////////////////////
1078 :
1079 0 : bool DiamondsEstablish3D::detachMarkers()
1080 : {
1081 0 : m_grid.detach_from_vertices(m_attElmGrpVrtx2BQnchd);
1082 :
1083 0 : m_grid.detach_from_faces(m_attMarkTwinFace);
1084 0 : m_grid.detach_from_faces(m_attMarkShiftFace);
1085 :
1086 0 : m_grid.detach_from_faces(m_attMarkShiftTriangleFace);
1087 0 : m_grid.detach_from_faces(m_attMarkShiftQuadriliteralFace);
1088 :
1089 0 : m_grid.detach_from_volumes(m_attMarkVolGetsShrinked);
1090 :
1091 0 : m_grid.detach_from_volumes( m_attVrtVec );
1092 :
1093 : // m_grid.detach_from_volumes(m_attVrtVecFace);
1094 :
1095 0 : m_grid.detach_from_vertices(m_attMarkVrtxIsCenterVrtx);
1096 :
1097 0 : m_grid.detach_from_vertices(m_attMarkVrtxIsShiftVrtx);
1098 :
1099 0 : m_grid.detach_from_edges(m_attMarkEdgeIsShiftEdge);
1100 :
1101 0 : m_grid.detach_from_vertices(m_attInfoVecSudosTouchingVrtx);
1102 :
1103 0 : m_grid.detach_from_vertices(m_attVrtxIsMidPtOfShiftVrtx);
1104 :
1105 0 : m_grid.detach_from_vertices(m_attMidPtVrtxOfShiftVrtx);
1106 :
1107 0 : m_grid.detach_from_vertices(m_attCenterVrtxOfShiftVrtx);
1108 :
1109 0 : m_grid.detach_from_edges(m_attEdgeCanBeRemoved);
1110 :
1111 0 : m_grid.detach_from_vertices(m_attNewSudoOfVrtx);
1112 :
1113 0 : m_grid.detach_from_edges( m_attCentersCutPts );
1114 :
1115 0 : m_grid.detach_from_edges(m_attDiamCentrVrtx);
1116 :
1117 0 : return true;
1118 : }
1119 :
1120 : ////////////////////////////////////////////////////////////////////////////
1121 :
1122 0 : bool DiamondsEstablish3D::assignBasicAtts()
1123 : {
1124 :
1125 0 : for( ElemGroupVrtx2BQuenched4Diams egv2bq : m_vecElemGroupVrtx2BQuenched )
1126 : {
1127 : Vertex * centerVrtx;
1128 : egv2bq.spuckOrigCenterVertex(centerVrtx);
1129 :
1130 0 : m_sel.select(centerVrtx);
1131 :
1132 0 : m_attAccsElmGrpVrtx2BQnchd[centerVrtx] = egv2bq;
1133 0 : m_attAccsVrtxIsCenterVrtx[centerVrtx] = true;
1134 : IndxVec sudoList = egv2bq.spuckSudoList();
1135 0 : m_attAccsInfoVecSudosTouchingVrtx[centerVrtx] = sudoList;
1136 :
1137 0 : m_sel.select( m_grid.associated_edges_begin(centerVrtx), m_grid.associated_edges_end(centerVrtx) );
1138 0 : m_sel.select( m_grid.associated_faces_begin(centerVrtx), m_grid.associated_faces_end(centerVrtx) );
1139 0 : m_sel.select( m_grid.associated_volumes_begin(centerVrtx), m_grid.associated_volumes_end(centerVrtx) );
1140 :
1141 0 : bool sudoAlreadyKnown = generateNewDiamSudos(centerVrtx, sudoList);
1142 :
1143 0 : }
1144 :
1145 0 : return true;
1146 : }
1147 :
1148 : ////////////////////////////////////////////////////////////////////////////
1149 :
1150 0 : bool DiamondsEstablish3D::trafoCollectedInfo2Attachments()
1151 : {
1152 : //for(VertexIterator iterV = m_sel.vertices_begin(); iterV != m_sel.vertices_end(); iterV++)
1153 0 : for( ElemGroupVrtx2BQuenched4Diams & egv2bQ : m_vecElemGroupVrtx2BQuenched )
1154 : {
1155 : Vertex * centerV;
1156 :
1157 : egv2bQ.spuckOrigCenterVertex(centerV);
1158 :
1159 : VecElems2BQuenched ve2bq;
1160 : egv2bQ.spuckVecElems2BQuenched4Diams(ve2bq);
1161 :
1162 0 : for( Elems2BQuenched & e2bq : ve2bq )
1163 : {
1164 : VecVolumeElementFaceQuintuplet vvef5;
1165 :
1166 : e2bq.spuckVecFullLowDimManifQuintuplet(vvef5);
1167 :
1168 0 : for( VolumeElementFaceQuintuplet & vef5 : vvef5 )
1169 : {
1170 0 : if( ! trafoQuintupleInfo2Attachments(vef5))
1171 : {
1172 : UG_LOG("Quintuples do not get to attachments" << std::endl);
1173 : return false;
1174 : }
1175 : }
1176 : // debugE2bQ(e2bq);
1177 :
1178 : Vertex * midPtVrtx;
1179 : e2bq.spuckMidPointOfShiftVrtcs(midPtVrtx);
1180 : m_attAccsVrtxIsMidPtOfShiftVrtcs[midPtVrtx] = true;
1181 :
1182 0 : m_sel.select(midPtVrtx);
1183 :
1184 : EdgePair ep;
1185 : e2bq.spuckPairLowDimElem(ep);
1186 : m_sel.select(ep.first);
1187 : m_sel.select(ep.second);
1188 :
1189 : VrtxPair vp;
1190 : e2bq.spuckShiftVrtcs(vp);
1191 : Vertex * shiftVrtxOne = vp.first;
1192 : Vertex * shiftVrtxTwo = vp.second;
1193 : m_sel.select(shiftVrtxOne);
1194 : m_sel.select(shiftVrtxTwo);
1195 :
1196 0 : m_attAccsMidPtVrtxOfShiftVrtx[shiftVrtxOne] = midPtVrtx;
1197 0 : m_attAccsMidPtVrtxOfShiftVrtx[shiftVrtxTwo] = midPtVrtx;
1198 :
1199 0 : m_attAccsCenterVrtxOfShiftVrtx[shiftVrtxOne] = centerV;
1200 0 : m_attAccsCenterVrtxOfShiftVrtx[shiftVrtxTwo] = centerV;
1201 :
1202 : // TODO FIXME hier die sudo Liste auch noch an den Center Vertex hängen als attachment
1203 :
1204 0 : }
1205 0 : }
1206 :
1207 :
1208 :
1209 : return true;
1210 : }
1211 :
1212 :
1213 : //////////////////////////////////////////////////////////////////////////////////
1214 :
1215 0 : bool DiamondsEstablish3D::trafoQuintupleInfo2Attachments(VolumeElementFaceQuintuplet & vef5 )
1216 : {
1217 : VrtxPair vp;
1218 : vef5.spuckShiftVrtcs(vp);
1219 :
1220 : m_attAccsVrtxIsShiftVrtx[vp.first] = true;
1221 : m_attAccsVrtxIsShiftVrtx[vp.second] = true;
1222 :
1223 : PairVolumeEdgeTwin pvet;
1224 : vef5.spuckPairFullLowDimTwin(pvet);
1225 :
1226 : VolumeEdgeTwin vetOne = pvet.first;
1227 : VolumeEdgeTwin vetTwo = pvet.second;
1228 :
1229 : Volume * volOne;
1230 : Volume * volTwo;
1231 : Edge * edgOne;
1232 : Edge * edgTwo;
1233 :
1234 : vetOne.spuckFullDimElem(volOne);
1235 : vetOne.spuckLowDimElem(edgOne);
1236 :
1237 : vetTwo.spuckFullDimElem(volTwo);
1238 : vetTwo.spuckLowDimElem(edgTwo);
1239 :
1240 : m_attAccsVolGetsShrinked[volOne] = true;
1241 : m_attAccsVolGetsShrinked[volTwo] = true;
1242 :
1243 : m_attAccsEdgeIsShiftEdge[edgOne] = true;
1244 : m_attAccsEdgeIsShiftEdge[edgTwo] = true;
1245 :
1246 0 : m_sel.select(volOne);
1247 : m_sel.select(volTwo);
1248 :
1249 : m_sel.select(edgOne);
1250 : m_sel.select(edgTwo);
1251 :
1252 : Face * fac;
1253 :
1254 : vef5.spuckManifElem(fac);
1255 :
1256 : m_attAccsFacIsTwinFac[fac] = true;
1257 :
1258 : m_sel.select(fac);
1259 :
1260 0 : return true;
1261 : }
1262 :
1263 : //////////////////////////////////////////////////////////////////////////////////
1264 :
1265 0 : bool DiamondsEstablish3D::assignMidPointOfShiftVrtcs(Elems2BQuenched & e2bq)
1266 : {
1267 : VrtxPair shiftVrtcs;
1268 :
1269 : e2bq.spuckShiftVrtcs(shiftVrtcs);
1270 :
1271 : Vertex * midPtVrtx;
1272 :
1273 0 : if( !createMidVrtx(shiftVrtcs, midPtVrtx))
1274 : {
1275 : UG_LOG("mid vertex already created " << std::endl);
1276 : // return false;
1277 : }
1278 :
1279 : e2bq.assignMidPointOfShiftVrtcs(midPtVrtx);
1280 :
1281 0 : return true;
1282 : }
1283 :
1284 : //////////////////////////////////////////////////////////////////////////////////
1285 :
1286 0 : bool DiamondsEstablish3D::createMidVrtx( VrtxPair const & shiftVrtcs, Vertex * & midPtVrtx )
1287 : {
1288 :
1289 : bool pairAlreadyCombined = false;
1290 :
1291 0 : for( CombiShiftVrtxMidVrtx & csvmv : m_vecCombiShiftVrtxMidVrtx )
1292 : {
1293 : VrtxPair prNoSwap;
1294 :
1295 : csvmv.spuckShiftVrtxPair(prNoSwap);
1296 :
1297 : if( prNoSwap == shiftVrtcs )
1298 : {
1299 : pairAlreadyCombined = true;
1300 : csvmv.spuckSinglVrtx(midPtVrtx);
1301 :
1302 : break;
1303 : }
1304 : else
1305 : {
1306 : std::swap( prNoSwap.first, prNoSwap.second );
1307 :
1308 : if( prNoSwap == shiftVrtcs )
1309 : {
1310 : pairAlreadyCombined = true;
1311 : csvmv.spuckSinglVrtx(midPtVrtx);
1312 :
1313 : break;
1314 : }
1315 :
1316 : if( pairAlreadyCombined )
1317 : break;
1318 :
1319 : }
1320 : }
1321 :
1322 : if( ! pairAlreadyCombined )
1323 : {
1324 :
1325 : Vertex * vrtxOne;
1326 : Vertex * vrtxTwo;
1327 :
1328 0 : vrtxOne = shiftVrtcs.first;
1329 0 : vrtxTwo = shiftVrtcs.second;
1330 :
1331 :
1332 : vector3 posOne = m_aaPos[vrtxOne];
1333 : vector3 posTwo = m_aaPos[vrtxTwo];
1334 :
1335 : vector3 sum;
1336 :
1337 : VecAdd(sum,posOne,posTwo);
1338 :
1339 : vector3 scal;
1340 :
1341 : VecScale(scal,sum,0.5);
1342 :
1343 0 : midPtVrtx = *m_grid.create<RegularVertex>();
1344 :
1345 : m_aaPos[midPtVrtx] = scal;
1346 :
1347 : CombiShiftVrtxMidVrtx csvmvNew( shiftVrtcs, midPtVrtx );
1348 :
1349 0 : m_vecCombiShiftVrtxMidVrtx.push_back(csvmvNew);
1350 : }
1351 :
1352 0 : return (! pairAlreadyCombined);
1353 :
1354 : }
1355 :
1356 : //////////////////////////////////////////////////////////////////////////////////
1357 :
1358 0 : bool DiamondsEstablish3D::createConditionForNewVrtcs()
1359 : {
1360 :
1361 : // iterate over all surrounding volumes to enable volume changes, this loop taken from SR but shortened
1362 0 : for(VolumeIterator iterSurrVol = m_sel.volumes_begin(); iterSurrVol != m_sel.volumes_end(); iterSurrVol++ )
1363 : {
1364 : Volume * sv = *iterSurrVol;
1365 :
1366 : std::vector<Vertex*> & newVrts = m_attAccsVrtVecVol[sv];
1367 0 : newVrts.resize(sv->num_vertices());
1368 :
1369 0 : for(size_t iVrt = 0; iVrt < sv->num_vertices(); iVrt++ )
1370 : {
1371 0 : newVrts[iVrt] = nullptr;
1372 : }
1373 : // erstmal so tun, als ob keine neuen Vertizes erzeugt werden an den alten Vertizes
1374 : }
1375 :
1376 : // for( FaceIterator iterFac = m_sel.faces_begin(); iterFac != m_sel.faces_end(); iterFac++)
1377 : // {
1378 : // Face * fac = *iterFac;
1379 : //
1380 : // std::vector<Vertex*> & newVrts = m_attAccsVrtVecFace[fac];
1381 : //
1382 : // newVrts.resize(fac->num_vertices());
1383 : //
1384 : // for(size_t iVrt = 0; iVrt < sv->num_vertices(); iVrt++ )
1385 : // {
1386 : // newVrts[iVrt] = nullptr;
1387 : // }
1388 : //
1389 : // }
1390 :
1391 0 : return true;
1392 : }
1393 :
1394 :
1395 : //////////////////////////////////////////////////////////////////////////////////
1396 :
1397 0 : bool DiamondsEstablish3D::distributeInfosForShrinkingVols()
1398 : {
1399 0 : for( ElemGroupVrtx2BQuenched4Diams & egv2bq : m_vecElemGroupVrtx2BQuenched )
1400 : {
1401 : VecElems2BQuenched ve2bq;
1402 :
1403 : egv2bq.spuckVecElems2BQuenched4Diams(ve2bq);
1404 :
1405 0 : for( Elems2BQuenched & e2bq : ve2bq )
1406 : {
1407 : VecVolumeElementFaceQuintuplet vvef5;
1408 :
1409 : e2bq.spuckVecFullLowDimManifQuintuplet(vvef5);
1410 :
1411 : Vertex * centerVrtx;
1412 :
1413 : egv2bq.spuckOrigCenterVertex(centerVrtx);
1414 :
1415 : Vertex * newMidVrtx;
1416 :
1417 : e2bq.spuckMidPointOfShiftVrtcs(newMidVrtx);
1418 :
1419 0 : for( VolumeElementFaceQuintuplet & vef5 : vvef5 )
1420 : {
1421 : PairVolumeEdgeTwin pvet;
1422 :
1423 : vef5.spuckPairFullLowDimTwin(pvet);
1424 :
1425 : VolumeEdgeTwin vetOne, vetTwo;
1426 :
1427 : vetOne = pvet.first;
1428 : vetTwo = pvet.second;
1429 :
1430 : Volume * volOne;
1431 : Volume * volTwo;
1432 :
1433 : vetOne.spuckFullDimElem(volOne);
1434 : vetTwo.spuckFullDimElem(volTwo);
1435 :
1436 0 : if( !teachMidVrtx2Vol(volOne, centerVrtx, newMidVrtx) || ! teachMidVrtx2Vol(volTwo,centerVrtx,newMidVrtx) )
1437 : {
1438 : UG_LOG("not taughtable" << std::endl);
1439 0 : return false;
1440 : }
1441 : }
1442 0 : }
1443 0 : }
1444 :
1445 : // TODO FIXME FACES detektieren, die verschoben werden müssen!!!
1446 :
1447 : return true;
1448 : }
1449 :
1450 :
1451 : //////////////////////////////////////////////////////////////////////////////////
1452 :
1453 0 : bool DiamondsEstablish3D::teachMidVrtx2Vol( Volume * const & vol, Vertex * const & origVrtx, Vertex * const & midVrtx )
1454 : {
1455 : IndexType taught = 0;
1456 :
1457 0 : std::vector<Vertex*> & newVrts4Fac = m_attAccsVrtVecVol[ vol ];
1458 :
1459 0 : for( IndexType indVrt = 0; indVrt < (vol)->num_vertices(); indVrt++ )
1460 : {
1461 0 : Vertex* volVrt = (vol)->vertex(indVrt);
1462 :
1463 0 : if( origVrtx == volVrt )
1464 : {
1465 0 : newVrts4Fac[ indVrt ] = midVrtx;
1466 0 : taught++;
1467 : }
1468 : }
1469 :
1470 0 : if( taught != 1 )
1471 : {
1472 : UG_LOG("strange problem with the volume " << taught << std::endl);
1473 0 : m_sh.assign_subset(vol, m_sh.num_subsets());
1474 0 : return false;
1475 : }
1476 :
1477 : return true;
1478 : }
1479 :
1480 : //////////////////////////////////////////////////////////////////////////////////
1481 :
1482 0 : bool DiamondsEstablish3D::shrinkVolumes()
1483 : {
1484 : // holds local side vertex indices
1485 : std::vector<size_t> locVrtInds;
1486 :
1487 : // first we create new edges from selected ones which are connected to
1488 : // inner vertices. This allows to preserve old subsets.
1489 : // Since we have to make sure that we use the right vertices,
1490 : // we have to iterate over the selected volumes and perform all actions on the edges
1491 : // of those volumes.
1492 :
1493 : int volNum = 0;
1494 :
1495 0 : for(VolumeIterator iter_sv = m_sel.volumes_begin(); iter_sv != m_sel.volumes_end(); ++iter_sv)
1496 : {
1497 : volNum++;
1498 :
1499 : Volume* sv = *iter_sv;
1500 :
1501 0 : UG_LOG("Diamond entering volume to create new elems " << CalculateCenter(sv, m_aaPos) << std::endl);
1502 :
1503 : // check for each edge whether it has to be copied.
1504 0 : for(size_t i_edge = 0; i_edge < sv->num_edges(); ++i_edge)
1505 : {
1506 0 : Edge* e = m_grid.get_edge(sv, i_edge);
1507 :
1508 0 : if(m_sel.is_selected(e))
1509 : {
1510 : // check the associated vertices through the volumes aaVrtVecVol attachment.
1511 : // If at least one has an associated new vertex and if no edge between the
1512 : // new vertices already exists, we'll create the new edge.
1513 : size_t ind0, ind1;
1514 0 : sv->get_vertex_indices_of_edge(ind0, ind1, i_edge);
1515 0 : Vertex* nv0 = (m_attAccsVrtVecVol[sv])[ind0];
1516 0 : Vertex* nv1 = (m_attAccsVrtVecVol[sv])[ind1];
1517 :
1518 0 : if(nv0 || nv1)
1519 : {
1520 : // if one vertex has no associated new one, then we use the vertex itself
1521 0 : if(!nv0)
1522 0 : nv0 = sv->vertex(ind0);
1523 0 : if(!nv1)
1524 0 : nv1 = sv->vertex(ind1);
1525 :
1526 : // create the new edge if it not already exists.
1527 0 : if( ! m_grid.get_edge(nv0, nv1))
1528 0 : m_grid.create_by_cloning(e, EdgeDescriptor(nv0, nv1), e);
1529 : }
1530 : }
1531 : }
1532 : }
1533 :
1534 :
1535 :
1536 : UG_LOG("Vol enter clone finished " << std::endl);
1537 :
1538 : // now we create new faces from selected ones which are connected to
1539 : // inner vertices. This allows to preserve old subsets.
1540 : // Since we have to make sure that we use the right vertices,
1541 : // we have to iterate over the selected volumes and perform all actions on the side-faces
1542 : // of those volumes.
1543 :
1544 0 : FaceDescriptor fd;
1545 :
1546 : int volNum2 = 0;
1547 :
1548 0 : for(VolumeIterator iter_sv = m_sel.volumes_begin(); iter_sv != m_sel.volumes_end(); ++iter_sv)
1549 : {
1550 : volNum2++;
1551 :
1552 : Volume* sv = *iter_sv;
1553 : // check for each face whether it has to be copied.
1554 :
1555 0 : UG_LOG("Diamond Face descriptor for vol " << CalculateCenter(sv, m_aaPos) << std::endl);
1556 :
1557 0 : for(size_t i_face = 0; i_face < sv->num_faces(); ++i_face)
1558 : {
1559 0 : Face* sf = m_grid.get_face(sv, i_face);
1560 :
1561 0 : if( m_sel.is_selected(sf))
1562 : {
1563 : // check the associated vertices through the volumes aaVrtVecVol attachment.
1564 : // If no face between the new vertices already exists, we'll create the new face.
1565 0 : sv->get_vertex_indices_of_face(locVrtInds, i_face);
1566 0 : fd.set_num_vertices(sf->num_vertices());
1567 :
1568 0 : for(size_t i = 0; i < sf->num_vertices(); ++i)
1569 : {
1570 0 : Vertex* nVrt = (m_attAccsVrtVecVol[sv])[locVrtInds[i]];
1571 :
1572 0 : if(nVrt)
1573 0 : fd.set_vertex(i, nVrt);
1574 : else
1575 0 : fd.set_vertex(i, sv->vertex(locVrtInds[i]));
1576 : }
1577 :
1578 : // if the new face does not already exist, we'll create it
1579 0 : if(!m_grid.get_face(fd))
1580 0 : m_grid.create_by_cloning(sf, fd, sf);
1581 : }
1582 : }
1583 : }
1584 :
1585 : UG_LOG("Face descriptor left" << std::endl);
1586 :
1587 : // Expand all faces.
1588 : // Since volumes are replaced on the fly, we have to take care with the iterator.
1589 : // record all new volumes in a vector. This will help to adjust positions later on.
1590 :
1591 : std::vector<Volume*> newFractureVolumes;
1592 : std::vector<IndexType> subsOfNewVolumes;
1593 :
1594 : // create alternative volumes where there are ending crossing clefts
1595 :
1596 0 : VolumeDescriptor vd;
1597 :
1598 : // beim Volumen fängt das Abfangen an, es muss abgefragt werden, ob das Volumen aus
1599 : // einem Segment ist, wo bisher falsch expandiert wird
1600 : // erst beim Face ab zu fangen ist viel zu spät TODO FIXME
1601 : // nicht zu vergessen, die Edges ordentlich sortiert zu sammeln, die zwei endende Vertizes enthalten,
1602 : // und dann zu splitten, wenn alle Attachements entfernt sind, dann Neustart anfordern mit neuer Geometrie
1603 :
1604 : int sudoNewVols = m_sh.num_subsets();
1605 :
1606 : // std::vector<Volume*> newVolsTwoCut;
1607 : // std::vector<Volume*> newVolsThreeCut;
1608 :
1609 :
1610 0 : for(VolumeIterator iter_sv = m_sel.volumes_begin(); iter_sv != m_sel.volumes_end();)
1611 : {
1612 :
1613 0 : Volume* sv = *iter_sv;
1614 : ++iter_sv;
1615 :
1616 :
1617 0 : UG_LOG("Diamond Volume new creation try at " << CalculateCenter(sv, m_aaPos) << std::endl);
1618 :
1619 :
1620 : // now expand the fracture faces of sv to volumes.
1621 0 : for(IndexType iSideVol = 0; iSideVol < sv->num_sides(); iSideVol++)
1622 : {
1623 : // get the local vertex indices of the side of the volume
1624 0 : sv->get_vertex_indices_of_face(locVrtInds, iSideVol);
1625 :
1626 0 : Face* sideFace = m_grid.get_side(sv, iSideVol);
1627 :
1628 0 : Volume * shiftVol = nullptr;
1629 :
1630 0 : if( sideFace )
1631 : {
1632 :
1633 0 : if( m_attAccsFacIsShiftFac[sideFace] )
1634 : {
1635 :
1636 0 : if( m_attAccsFacIsShiftQuadriliteralFac[sideFace] )
1637 : {
1638 :
1639 : std::vector<Vertex*> centerVrtcs;
1640 : std::vector<Vertex*> shiftVrtcs;
1641 : std::vector<Vertex*> midPtVrtcs;
1642 :
1643 :
1644 0 : for( IndexType vrtxInd = 0; vrtxInd < sideFace->num_vertices(); vrtxInd++ )
1645 : {
1646 0 : Vertex * sideVrtx = sideFace->vertex(vrtxInd);
1647 :
1648 0 : if( m_attAccsVrtxIsCenterVrtx[sideVrtx] )
1649 : {
1650 0 : centerVrtcs.push_back(sideVrtx);
1651 : }
1652 0 : else if( m_attAccsVrtxIsShiftVrtx[sideVrtx] )
1653 : {
1654 0 : shiftVrtcs.push_back(sideVrtx);
1655 : }
1656 : else
1657 : {
1658 0 : UG_LOG("not center not shift but shift face? "
1659 : << CalculateCenter( sideFace, m_aaPos ) << std::endl );
1660 0 : return false;
1661 : }
1662 : }
1663 :
1664 0 : if( centerVrtcs.size() != 2 || shiftVrtcs.size() != 2 )
1665 : {
1666 0 : UG_LOG("strange face behaviour for shift face of fac and vol "
1667 : << CalculateCenter( sideFace, m_aaPos )
1668 : << " ---- "
1669 : << CalculateCenter( sv, m_aaPos )
1670 : << std::endl);
1671 0 : return false;
1672 : }
1673 :
1674 : std::swap( shiftVrtcs[0], shiftVrtcs[1] );
1675 :
1676 0 : if( ! findShiftFaceVertices( sv, centerVrtcs, midPtVrtcs ))
1677 : {
1678 0 : UG_LOG("vertices of shift face strange "
1679 : << CalculateCenter( sideFace, m_aaPos ) << std::endl);
1680 0 : return false;
1681 : }
1682 :
1683 0 : UG_LOG("Diamentenerzeugung fuer " << CalculateCenter( sideFace, m_aaPos ) << std::endl);
1684 :
1685 : // centers.push_back(centerVrtcs);
1686 : // shifts.push_back(shiftVrtcs);
1687 : // midPts.push_back(midPtVrtcs);
1688 :
1689 : // // sehr sonderbar, wieso folgende Reihenfolge notwendig, damit es klappt:
1690 : //
1691 : // shiftVol = *m_grid.create<Prism>(
1692 : // PrismDescriptor( centerVrtcs[0], midPtVrtcs[0], shiftVrtcs[1],
1693 : // centerVrtcs[1], midPtVrtcs[1], shiftVrtcs[0]
1694 : // )
1695 : // );
1696 : // // sehr sonderbar, wieso folgende Reihenfolge notwendig, damit es klappt:
1697 :
1698 0 : shiftVol = *m_grid.create<Prism>(
1699 0 : PrismDescriptor( centerVrtcs[0], midPtVrtcs[0], shiftVrtcs[0],
1700 : centerVrtcs[1], midPtVrtcs[1], shiftVrtcs[1]
1701 : )
1702 : );
1703 :
1704 : int isThreeCross = -1;
1705 : IndexType foundThreeCross = 0;
1706 :
1707 : VrtxIndxCombi vrtxSudosCombiCenters;
1708 :
1709 0 : for( int i = 0; i < 2; i++ )
1710 : {
1711 0 : Vertex * cntrVrtx = centerVrtcs[i];
1712 :
1713 0 : IndxVec sudoVec = m_attAccsInfoVecSudosTouchingVrtx[cntrVrtx];
1714 :
1715 0 : IndexType sudoNums = sudoVec.size();
1716 :
1717 : VrtxIndxPair centerVrtxSudos( cntrVrtx, sudoVec );
1718 :
1719 0 : vrtxSudosCombiCenters.push_back(centerVrtxSudos);
1720 :
1721 0 : if( sudoNums < 2 )
1722 : {
1723 : UG_LOG("no cross but new prism " << std::endl);
1724 : return false;
1725 : }
1726 0 : else if( sudoNums == 2 )
1727 : {
1728 : ;
1729 : }
1730 0 : else if( sudoNums == 3 )
1731 : {
1732 : isThreeCross = i;
1733 0 : foundThreeCross++;
1734 : }
1735 : else
1736 : {
1737 : UG_LOG("strange sudo number " << std::endl);
1738 : return false;
1739 : }
1740 0 : }
1741 :
1742 0 : if( foundThreeCross > 0 )
1743 : {
1744 0 : if( foundThreeCross != 1 )
1745 : {
1746 : UG_LOG("three cross vertices too much " << std::endl);
1747 0 : return false;
1748 : }
1749 :
1750 0 : CombiNewVolsProps cnvp( shiftVol, sideFace, vrtxSudosCombiCenters, shiftVrtcs, midPtVrtcs, isThreeCross );
1751 :
1752 0 : m_vecCombiNewVolsThreeCross.push_back(cnvp);
1753 0 : }
1754 : else
1755 : {
1756 0 : CombiNewVolsProps cnvp( shiftVol, sideFace, vrtxSudosCombiCenters, shiftVrtcs, midPtVrtcs );
1757 :
1758 0 : m_vecCombiNewVolsTwoCross.push_back(cnvp);
1759 0 : }
1760 :
1761 : //m_attAccsInfoVecSudosTouchingVrtx
1762 :
1763 : // newVols.push_back(shiftVol);
1764 : // numFacs++;
1765 :
1766 : // if( numFacs == 1 )
1767 : // return true;
1768 : // Volume * shiftVol2 = *m_grid.create<Prism>(
1769 : // PrismDescriptor(
1770 : // shiftVrtcs[1], midPtVrtcs[1], centerVrtcs[1],
1771 : // shiftVrtcs[0], midPtVrtcs[0], centerVrtcs[0]
1772 : // )
1773 : // );
1774 :
1775 : // m_sh.assign_subset(shiftVol,sudoNewVols);
1776 : // m_sh.assign_subset(shiftVol2,m_sh.num_subsets()+2);
1777 :
1778 0 : UG_LOG("Diamentenerzeugung geschafft " << CalculateCenter( sideFace, m_aaPos ) << std::endl);
1779 :
1780 : // return true;
1781 :
1782 0 : }
1783 :
1784 : // if( shiftVol && shiftVol2 )
1785 : // {
1786 : // m_sh.assign_subset(shiftVol, sudoNewVols);
1787 : // m_sh.assign_subset(shiftVol2, sudoNewVols+1);
1788 : //
1789 : // m_sh.assign_subset( centerVrtcs[0], m_sh.num_subsets() );
1790 : // m_sh.assign_subset( centerVrtcs[1], m_sh.num_subsets() );
1791 : // m_sh.assign_subset( midPtVrtcs[0], m_sh.num_subsets() );
1792 : // m_sh.assign_subset( midPtVrtcs[1], m_sh.num_subsets() );
1793 : // m_sh.assign_subset( shiftVrtcs[0], m_sh.num_subsets() );
1794 : // m_sh.assign_subset( shiftVrtcs[1], m_sh.num_subsets() );
1795 : //
1796 : //
1797 : // return true;
1798 : // }
1799 : }
1800 : }
1801 :
1802 : // if( shiftVol )
1803 : // {
1804 : // m_sh.assign_subset(shiftVol,m_sh.num_subsets());
1805 : // }
1806 :
1807 : }
1808 :
1809 : // now set up a new volume descriptor and replace the volume.
1810 0 : if(vd.num_vertices() != sv->num_vertices())
1811 0 : vd.set_num_vertices(sv->num_vertices());
1812 :
1813 0 : for(size_t i_vrt = 0; i_vrt < sv->num_vertices(); ++i_vrt)
1814 : {
1815 0 : if( (m_attAccsVrtVecVol[sv])[i_vrt] )
1816 0 : vd.set_vertex(i_vrt, (m_attAccsVrtVecVol[sv])[i_vrt]);
1817 : else
1818 0 : vd.set_vertex(i_vrt, sv->vertex(i_vrt));
1819 : }
1820 :
1821 0 : m_grid.create_by_cloning(sv, vd, sv);
1822 0 : m_grid.erase(sv);
1823 :
1824 : }
1825 :
1826 : // for( Volume * v : newVols )
1827 : // {
1828 : // m_sh.assign_subset(v,m_sh.num_subsets());
1829 : // }
1830 :
1831 : UG_LOG("Volumes erzeugt " << std::endl);
1832 :
1833 0 : for( EdgeIterator itEdg = m_sel.begin<Edge>(); itEdg != m_sel.end<Edge>(); )
1834 : {
1835 : Edge * edg = *itEdg;
1836 : itEdg++;
1837 :
1838 0 : if( m_attAccsEdgeCanBeRemoved[edg] )
1839 : {
1840 0 : m_grid.erase(edg);
1841 : }
1842 : }
1843 :
1844 : UG_LOG("unnoetige Ecken entfernt " << std::endl);
1845 :
1846 0 : for(FaceIterator iter = m_sel.begin<Face>(); iter != m_sel.end<Face>();)
1847 : {
1848 : Face* fac = *iter;
1849 : ++iter;
1850 :
1851 0 : if( ! m_attAccsFacIsShiftFac[fac] )
1852 0 : m_grid.erase(fac);
1853 : }
1854 :
1855 : UG_LOG("Gesichter entfernt " << std::endl);
1856 :
1857 :
1858 : // int suse = m_sh.num_subsets();
1859 : //
1860 : // for( int i = 0; i < centers.size(); i++ )
1861 : // {
1862 : // std::vector<Vertex*> centerVrtcs = centers[i];
1863 : // std::vector<Vertex*> shiftVrtcs = shifts[i];
1864 : // std::vector<Vertex*> midPtVrtcs = midPts[i];
1865 : //
1866 : // UG_LOG("versuche neues Volumen nr " << i << std::endl);
1867 : //
1868 : // Volume * shiftVol = *m_grid.create<Prism>(
1869 : // PrismDescriptor( centerVrtcs[0], shiftVrtcs[0], midPtVrtcs[0],
1870 : // centerVrtcs[1], shiftVrtcs[1], midPtVrtcs[1]
1871 : // )
1872 : // );
1873 : //
1874 : // m_sh.assign_subset(shiftVol, suse);
1875 : //
1876 : // UG_LOG("neue Volumen " << i << CalculateCenter(shiftVol,m_aaPos) << std::endl);
1877 : //
1878 : // }
1879 :
1880 : UG_LOG("neue Volumen erzeugt " << std::endl);
1881 :
1882 : return true;
1883 0 : }
1884 :
1885 : //////////////////////////////////////////////////////////////////////////////////
1886 :
1887 0 : bool DiamondsEstablish3D::determineShiftFaces()
1888 : {
1889 :
1890 0 : for(FaceIterator iterF = m_sel.begin<Face>(); iterF != m_sel.end<Face>();)
1891 : {
1892 : Face* fac = *iterF;
1893 :
1894 : IndexType numShiftVrcs = 0;
1895 : IndexType numCentrVrtcs = 0;
1896 :
1897 0 : if( typeid(*fac) == typeid(Triangle) )
1898 : {
1899 0 : UG_LOG("we have a triangle " << CalculateCenter( fac, m_aaPos ) << " -> " << typeid(*fac).name() << std::endl);
1900 : }
1901 0 : else if( typeid(*fac) == typeid(Quadrilateral) )
1902 : {
1903 0 : UG_LOG("we have a quadriliteral " << CalculateCenter( fac, m_aaPos ) << " -> " << typeid(*fac).name() << std::endl);
1904 : }
1905 : else
1906 : {
1907 0 : UG_LOG("we have whatever at " << CalculateCenter( fac, m_aaPos ) << " -> " << typeid(*fac).name() << std::endl);
1908 : }
1909 :
1910 0 : for( IndexType iVrt = 0; iVrt < fac->num_vertices(); iVrt++ )
1911 : {
1912 0 : Vertex * vrt = fac->vertex(iVrt);
1913 :
1914 0 : if( m_attAccsVrtxIsCenterVrtx[vrt] )
1915 : {
1916 0 : numCentrVrtcs++;
1917 : }
1918 0 : else if( m_attAccsVrtxIsShiftVrtx[vrt] )
1919 : {
1920 0 : numShiftVrcs++;
1921 : }
1922 : }
1923 :
1924 0 : if( numShiftVrcs == 2 && numCentrVrtcs == 2 )
1925 : {
1926 : m_attAccsFacIsShiftFac[fac] = true;
1927 : m_attAccsFacIsShiftQuadriliteralFac[fac] = true;
1928 :
1929 0 : if( typeid(*fac) != typeid(Quadrilateral) )
1930 : {
1931 0 : UG_LOG("4 vertices but not Quadriliteral? " << CalculateCenter( fac, m_aaPos ) << " -> " << typeid(*fac).name() << std::endl);
1932 : }
1933 :
1934 : }
1935 0 : else if( numShiftVrcs == 2 && numCentrVrtcs == 1 )
1936 : {
1937 : m_attAccsFacIsShiftFac[fac] = true;
1938 : m_attAccsFacIsShiftTriangleFac[fac] = true;
1939 0 : UG_LOG("was ist das für ein Typ " << CalculateCenter( fac, m_aaPos ) << " -> " << typeid(*fac).name() << std::endl );
1940 : }
1941 :
1942 0 : if( numCentrVrtcs > 2 || numShiftVrcs > 2 )
1943 : {
1944 : UG_LOG("too much shift or center vertices in one face " << std::endl);
1945 0 : UG_LOG("weia was ist das für ein Typ " << CalculateCenter( fac, m_aaPos ) << " -> " << typeid(*fac).name() << std::endl );
1946 : return false;
1947 : }
1948 :
1949 : iterF++;
1950 : }
1951 :
1952 : return true;
1953 : }
1954 :
1955 : //////////////////////////////////////////////////////////////////////////////////
1956 :
1957 : //bool DiamondsEstablish3D::findShiftFaceVertices(
1958 : // std::vector<Vertex*> & centerVrtcs,
1959 : // std::vector<Vertex*> & shiftVrtcs,
1960 : // std::vector<Vertex*> & midPtVrtcs
1961 : //)
1962 : //{
1963 : //
1964 : // // figure out to which center vertex which shift vertex belongs
1965 : //
1966 : // if( centerVrtcs.size() != 2 || shiftVrtcs.size() != 2 )
1967 : // {
1968 : // UG_LOG("strange face behaviour for shift face" << std::endl);
1969 : //// << CalculateCenter( fac, m_aaPos ) << std::endl);
1970 : // return false;
1971 : // }
1972 : //
1973 : // if( ! checkAttsOfShiftFaceVrtcs( centerVrtcs, shiftVrtcs) )
1974 : // {
1975 : // std::swap( shiftVrtcs[0], shiftVrtcs[1] );
1976 : //
1977 : // if( ! checkAttsOfShiftFaceVrtcs( centerVrtcs, shiftVrtcs) )
1978 : // {
1979 : // UG_LOG("shift and swap vertices do not agree independent of ordering" << std::endl);
1980 : // return false;
1981 : // }
1982 : // }
1983 : //
1984 : // for( Vertex * & spv : shiftVrtcs )
1985 : // {
1986 : // midPtVrtcs.push_back( m_attAccsMidPtVrtxOfShiftVrtx[spv] );
1987 : // }
1988 : //
1989 : // return true;
1990 : //}
1991 :
1992 : ////////////////////////////////////////////////////////////////////////////////////
1993 :
1994 0 : bool DiamondsEstablish3D::checkAttsOfShiftFaceVrtcs( std::vector<Vertex*> const & centerVrtcs, std::vector<Vertex*> const & shiftVrtcs )
1995 : {
1996 : Vertex * const & centerVrtxOne = centerVrtcs[0];
1997 : Vertex * const & centerVrtxTwo = centerVrtcs[1];
1998 :
1999 : Vertex * const & shiftVrtxOne = shiftVrtcs[0];
2000 : Vertex * const & shiftVrtxTwo = shiftVrtcs[1];
2001 :
2002 0 : if( m_attAccsCenterVrtxOfShiftVrtx[shiftVrtxOne] == centerVrtxOne )
2003 : {
2004 0 : if( m_attAccsCenterVrtxOfShiftVrtx[shiftVrtxTwo] == centerVrtxTwo )
2005 : {
2006 : return true;
2007 : }
2008 : else
2009 : {
2010 : // UG_LOG("shift and center vertex do not fit together " << std::endl);
2011 0 : return false;
2012 : }
2013 : }
2014 :
2015 : return false;
2016 : }
2017 :
2018 : ////////////////////////////////////////////////////////////////////////////////////
2019 :
2020 0 : bool DiamondsEstablish3D::findShiftFaceVertices( Volume * & vol,
2021 : std::vector<Vertex*> & centerVrtcs,
2022 : std::vector<Vertex*> & midPtVrtcs
2023 : )
2024 : {
2025 0 : midPtVrtcs = std::vector<Vertex*>(2);
2026 :
2027 0 : std::vector<Vertex*> & newVrts4Fac = m_attAccsVrtVecVol[ vol ];
2028 :
2029 0 : for( IndexType indVrt = 0; indVrt < (vol)->num_vertices(); indVrt++ )
2030 : {
2031 0 : Vertex* volVrt = (vol)->vertex(indVrt);
2032 :
2033 0 : for( int i = 0; i < 2; i++ )
2034 : {
2035 0 : Vertex * ceVe = centerVrtcs[i];
2036 :
2037 0 : if( ceVe == volVrt )
2038 : {
2039 0 : midPtVrtcs[i] = newVrts4Fac[indVrt];
2040 : }
2041 : }
2042 : }
2043 :
2044 :
2045 0 : return true;
2046 : }
2047 :
2048 : ////////////////////////////////////////////////////////////////////////////////////
2049 :
2050 0 : bool DiamondsEstablish3D::assignSudos2DiamsPreform()
2051 : {
2052 0 : for( CombiNewVolsProps & cnvp : m_vecCombiNewVolsTwoCross )
2053 : {
2054 : VrtxIndxCombi viCombiCenter;
2055 : cnvp.spuckCenterVrtcsSudos(viCombiCenter);
2056 : Volume * vol;
2057 : cnvp.spuckFulldimElem(vol);
2058 :
2059 : VrtxIndxPair vrtxIndxPr( viCombiCenter[0] );
2060 :
2061 : Vertex * center = vrtxIndxPr.first;
2062 0 : IndxVec sudoVec = vrtxIndxPr.second;
2063 :
2064 0 : for( CombiCntrVrtxSudo & ccvs: m_vecCombiCntrVrtxSudo )
2065 : {
2066 0 : if( ccvs.spuckSudoVec() == sudoVec )
2067 : {
2068 : // m_sh.assign_subset(vol, ccvs.spuckNewSudo());
2069 0 : if( ! assignSudoOfNewVols2VolAndSubElems(vol,ccvs.spuckNewSudo() ) )
2070 : {
2071 : UG_LOG("sudos not assignable" << std::endl);
2072 : return false;
2073 : }
2074 : break;
2075 : }
2076 : }
2077 :
2078 0 : }
2079 :
2080 0 : for( CombiNewVolsProps & cnvp : m_vecCombiNewVolsThreeCross )
2081 : {
2082 : VrtxIndxCombi viCombiCenter;
2083 : cnvp.spuckCenterVrtcsSudos(viCombiCenter);
2084 : Volume * vol;
2085 : cnvp.spuckFulldimElem(vol);
2086 :
2087 0 : IndexType i3c = cnvp.spuckThreeCrossIndex();
2088 :
2089 0 : VrtxIndxPair vrtxIndxPr( viCombiCenter[i3c] );
2090 :
2091 : Vertex * center = vrtxIndxPr.first;
2092 0 : IndxVec sudoVec = vrtxIndxPr.second;
2093 :
2094 0 : for( CombiCntrVrtxSudo & ccvs: m_vecCombiCntrVrtxSudo )
2095 : {
2096 0 : if( ccvs.spuckSudoVec() == sudoVec )
2097 : {
2098 : // m_sh.assign_subset(vol, ccvs.spuckNewSudo());
2099 0 : if( ! assignSudoOfNewVols2VolAndSubElems(vol,ccvs.spuckNewSudo() ) )
2100 : {
2101 : UG_LOG("sudos not assignable" << std::endl);
2102 : return false;
2103 : }
2104 : break;
2105 : }
2106 : }
2107 :
2108 0 : }
2109 :
2110 : return true;
2111 :
2112 : }
2113 :
2114 : /////////////////////////////////////////////////////////////////
2115 :
2116 0 : bool DiamondsEstablish3D::splitDiamsPreform2Diams()
2117 : {
2118 : UG_LOG("try to split diam vols " << std::endl);
2119 :
2120 0 : for( CombiNewVolsProps & cnvp : m_vecCombiNewVolsThreeCross )
2121 : {
2122 0 : if( ! determineSplitPts(cnvp) )
2123 : {
2124 : UG_LOG("split points determination not possible");
2125 : return false;
2126 : }
2127 : }
2128 :
2129 0 : for( Edge * const & edg : m_edges2BDeletedAtLastStep )
2130 : {
2131 0 : if( ! determineCutPtAtEdge(edg))
2132 : {
2133 : UG_LOG("edge point not determinable" << std::endl);
2134 : return false;
2135 : }
2136 : }
2137 :
2138 0 : for( CombiNewVolsProps & cnvp : m_vecCombiNewVolsThreeCross )
2139 : {
2140 0 : if( ! splitThreeCrossLargeDiams(cnvp) )
2141 : {
2142 : UG_LOG("splitting not possible");
2143 : return false;
2144 : }
2145 : }
2146 :
2147 :
2148 0 : for( Face * fac : m_faces2BDeletedAtLastStep )
2149 : {
2150 0 : if( ! fac )
2151 : {
2152 : UG_LOG("want to delete null fac " << std::endl);
2153 : return false;
2154 : }
2155 :
2156 0 : m_grid.erase(fac);
2157 : }
2158 :
2159 0 : for( Edge * edg : m_edges2BDeletedAtLastStep )
2160 : {
2161 0 : if( ! edg )
2162 : {
2163 : UG_LOG("want to delete edge unexisting at end " << std::endl);
2164 : return false;
2165 : }
2166 :
2167 0 : m_grid.erase(edg);
2168 : }
2169 :
2170 :
2171 : UG_LOG("all diam vols splitted " << std::endl);
2172 :
2173 : // else
2174 : // {
2175 : // UG_LOG("No split of preform diamonds, remaining at preform " << std::endl);
2176 : // }
2177 :
2178 0 : return true;
2179 : }
2180 :
2181 : ////////////////////////////////////////////////////////////////////////////////////
2182 :
2183 0 : bool DiamondsEstablish3D::determineCutPtAtEdge( Edge * const & edg )
2184 : {
2185 :
2186 0 : std::vector<vector3> cutPtVec = m_attAccsCentersCutPts[edg];
2187 :
2188 0 : IndexType cutPtVecSize = cutPtVec.size();
2189 :
2190 0 : if( cutPtVec.size() == 0 )
2191 : {
2192 : UG_LOG("no cut vectors " << std::endl);
2193 : return false;
2194 : }
2195 :
2196 : // average the cut points
2197 :
2198 : vector3 vecSum;
2199 :
2200 0 : for( vector3 const & vec2Add : cutPtVec )
2201 : {
2202 : vector3 tempSum = vecSum;
2203 :
2204 : VecAdd(vecSum, tempSum, vec2Add );
2205 : }
2206 :
2207 0 : number scale = 1. / static_cast<number>(cutPtVecSize);
2208 :
2209 : vector3 avrgVec;
2210 :
2211 : VecScale(avrgVec, vecSum, scale);
2212 :
2213 0 : Vertex * cutVrtx = *m_grid.create<RegularVertex>();
2214 :
2215 : m_aaPos[cutVrtx] = avrgVec;
2216 :
2217 0 : m_attAccsDiamCentrVrtx[edg] = cutVrtx;
2218 :
2219 0 : return true;
2220 0 : }
2221 :
2222 : ////////////////////////////////////////////////////////////////////////////////////
2223 :
2224 :
2225 0 : bool DiamondsEstablish3D::generateNewDiamSudos(Vertex * & centerV, IndxVec sudoList )
2226 : {
2227 0 : bool sudoCombiUnKnown = addElem(m_sudosTable, sudoList);
2228 :
2229 0 : if( ! sudoCombiUnKnown )
2230 : {
2231 0 : for( CombiCntrVrtxSudo & ccvs : m_vecCombiCntrVrtxSudo )
2232 : {
2233 0 : if( ccvs.spuckSudoVec() == sudoList )
2234 : {
2235 : ccvs.schluckVertex(centerV);
2236 0 : m_sh.assign_subset(centerV, ccvs.spuckNewSudo());
2237 0 : m_attAccsNewSudoOfVrtx[centerV] = ccvs.spuckNewSudo();
2238 0 : break;
2239 : }
2240 : }
2241 : }
2242 : else
2243 : {
2244 0 : IndexType newSudo = m_sh.num_subsets();
2245 0 : m_sh.assign_subset(centerV,newSudo);
2246 :
2247 0 : std::string sudoName = std::string("diamond_");
2248 0 : for( IndexType & sd : sudoList )
2249 : {
2250 0 : sudoName += std::string("_") + std::string( const_cast<char*>( m_sh.get_subset_name( sd ) ) );
2251 : }
2252 :
2253 0 : m_sh.set_subset_name(sudoName.c_str(), newSudo);
2254 :
2255 0 : m_attAccsNewSudoOfVrtx[centerV] = newSudo;
2256 : CombiCntrVrtxSudo ccvs( sudoList, newSudo );
2257 0 : m_vecCombiCntrVrtxSudo.push_back(ccvs);
2258 : UG_LOG("NNNNNNNNNNNNNNNNNNNN" << std::endl);
2259 : UG_LOG("creating new sudo of " << std::endl);
2260 0 : for( auto & i : sudoList )
2261 : {
2262 0 : UG_LOG("list part " << i << std::endl);
2263 : }
2264 : UG_LOG("DDDDDDDDDDDDDDDDDDDD" << std::endl);
2265 :
2266 : }
2267 :
2268 0 : return ( sudoCombiUnKnown );
2269 : }
2270 :
2271 : ////////////////////////////////////////////////////////////////////////////////////
2272 :
2273 0 : bool DiamondsEstablish3D::assignSudoOfNewVols2VolAndSubElems(Volume * & vol, IndexType sudo)
2274 : {
2275 0 : if( ! vol )
2276 : {
2277 : UG_LOG("vol to assign sudo null " << std::endl);
2278 0 : return false;
2279 : }
2280 :
2281 0 : m_sh.assign_subset(vol, sudo);
2282 :
2283 0 : for(IndexType iFace = 0; iFace < vol->num_faces(); ++iFace)
2284 : {
2285 0 : Face * fac = m_grid.get_face(vol, iFace);
2286 :
2287 0 : if( ! fac )
2288 : {
2289 : UG_LOG("face null sudo " << std::endl);
2290 0 : return false;
2291 : }
2292 :
2293 0 : m_sh.assign_subset( fac, sudo );
2294 :
2295 : }
2296 :
2297 0 : for(IndexType iEdge = 0; iEdge < vol->num_edges(); ++iEdge)
2298 : {
2299 0 : Edge* edg = m_grid.get_edge(vol, iEdge);
2300 :
2301 0 : if( ! edg )
2302 : {
2303 : UG_LOG("edge null sudo " << std::endl);
2304 0 : return false;
2305 : }
2306 :
2307 0 : m_sh.assign_subset( edg, sudo );
2308 :
2309 : }
2310 :
2311 0 : for( IndexType iVrt = 0; iVrt < vol->num_vertices(); iVrt++ )
2312 : {
2313 0 : Vertex * vrt = vol->vertex(iVrt);
2314 :
2315 0 : if( !vrt )
2316 : {
2317 : UG_LOG("vertex null sudo " << std::endl);
2318 0 : return false;
2319 : }
2320 :
2321 0 : m_sh.assign_subset( vrt, sudo );
2322 : }
2323 :
2324 : return true;
2325 : }
2326 :
2327 :
2328 : ///////////////////////////////////////////////////////////////////////////////////
2329 :
2330 0 : bool DiamondsEstablish3D::detectRemovableEdges()
2331 : {
2332 :
2333 0 : for( FaceIterator itFac = m_sel.begin<Face>(); itFac != m_sel.end<Face>(); itFac++)
2334 : {
2335 : Face * fac = *itFac;
2336 :
2337 : if( ! fac )
2338 : {
2339 : UG_LOG("face does not exist " << std::endl);
2340 : return false;
2341 : }
2342 :
2343 : UG_LOG("looking into twin facs try " << std::endl);
2344 :
2345 0 : if( m_attAccsFacIsTwinFac[fac] )
2346 : {
2347 : UG_LOG("inside twin facs try " << std::endl);
2348 :
2349 0 : for( IndexType iEdg = 0; iEdg < fac->num_edges(); iEdg++)
2350 : {
2351 : UG_LOG("inside edges of twins try " << std::endl);
2352 :
2353 0 : Edge * edg = m_grid.get_edge(fac,iEdg);
2354 :
2355 0 : if( ! edg )
2356 : {
2357 : UG_LOG("edge does not exist " << std::endl);
2358 0 : return false;
2359 : }
2360 :
2361 : IndexType numAssoCenterVrcs = 0;
2362 :
2363 0 : for( IndexType iV = 0; iV < 2; iV++ )
2364 : {
2365 0 : Vertex * vrt = edg->vertex(iV);
2366 :
2367 0 : if( ! vrt )
2368 : {
2369 : UG_LOG("vertex does not exist " << std::endl);
2370 0 : return false;
2371 : }
2372 :
2373 0 : if( m_attAccsVrtxIsCenterVrtx[vrt] )
2374 : {
2375 0 : numAssoCenterVrcs++;
2376 : }
2377 : }
2378 :
2379 : UG_LOG("number of associatec center vertices is " << numAssoCenterVrcs << std::endl);
2380 :
2381 0 : if( numAssoCenterVrcs == 1 )
2382 : {
2383 : m_attAccsEdgeCanBeRemoved[edg] = true;
2384 : }
2385 :
2386 : }
2387 : }
2388 : }
2389 :
2390 0 : return true;
2391 : }
2392 :
2393 : ///////////////////////////////////////////////////////////////////////////////////
2394 :
2395 0 : bool DiamondsEstablish3D::determineSplitPts( CombiNewVolsProps & combiNewVolsProps )
2396 : {
2397 : VrtxIndxCombi viCombiCenter;
2398 : combiNewVolsProps.spuckCenterVrtcsSudos(viCombiCenter);
2399 : Volume * vol;
2400 : combiNewVolsProps.spuckFulldimElem(vol);
2401 :
2402 0 : IndexType ind3CrossSide = combiNewVolsProps.spuckThreeCrossIndex();
2403 :
2404 : IndexType indOtherSide = ( ind3CrossSide + 1 ) % 2;
2405 :
2406 0 : VrtxIndxPair diamCenterVrtxPr( viCombiCenter[ind3CrossSide] );
2407 :
2408 0 : Vertex * centerDiam = diamCenterVrtxPr.first;
2409 :
2410 0 : VrtxIndxPair diamOtherVrtxPr( viCombiCenter[ indOtherSide ] );
2411 :
2412 0 : Vertex * centerOther = diamOtherVrtxPr.first;
2413 :
2414 : std::vector<Vertex*> midPtVrtcs;
2415 : combiNewVolsProps.spuckMidPtVrtcs(midPtVrtcs);
2416 :
2417 0 : Vertex * midPtDiam = midPtVrtcs[ind3CrossSide];
2418 : Vertex * midPtOther = midPtVrtcs[indOtherSide];
2419 :
2420 0 : vector3 centerDiamVec3 = m_aaPos[centerDiam];
2421 : vector3 midPtDiamVec3 = m_aaPos[midPtDiam];
2422 0 : vector3 centerOtherVec3 = m_aaPos[centerOther];
2423 :
2424 : vector3 cutPtVec3;
2425 :
2426 0 : DropAPerpendicular(cutPtVec3, midPtDiamVec3, centerDiamVec3, centerOtherVec3);
2427 :
2428 : // figure out the edge to which this point belongs
2429 :
2430 0 : Edge * edgOfVrtcs = nullptr;
2431 :
2432 0 : if( ! figureOutEdgeOfVrtcs( edgOfVrtcs, centerDiam, centerOther ) )
2433 : {
2434 : UG_LOG("vertices do not belong to an edge " << std::endl);
2435 : return false;
2436 : }
2437 :
2438 0 : if( edgOfVrtcs != nullptr )
2439 : {
2440 : // addElem(m_edgesCut4Diams, edgOfVrtcs);
2441 0 : addElem(m_edges2BDeletedAtLastStep,edgOfVrtcs);
2442 :
2443 0 : m_attAccsCentersCutPts[edgOfVrtcs].push_back( cutPtVec3 );
2444 :
2445 : if( ! combiNewVolsProps.schluckLowdimElem( edgOfVrtcs ) )
2446 : {
2447 : UG_LOG("edge already known" << std::endl);
2448 : return false;
2449 : }
2450 : }
2451 : else
2452 : {
2453 : UG_LOG("no edge found " << std::endl);
2454 : return false;
2455 : }
2456 :
2457 :
2458 : return true;
2459 :
2460 0 : }
2461 :
2462 : ///////////////////////////////////////////////////////////////////////////////////
2463 :
2464 0 : bool DiamondsEstablish3D::figureOutEdgeOfVrtcs( Edge * & edgOfVrtcs, Vertex * const & vertexOne, Vertex * const & vertexTwo )
2465 : {
2466 :
2467 0 : for( Grid::AssociatedEdgeIterator itEdgOne = m_grid.associated_edges_begin( vertexOne );
2468 0 : itEdgOne != m_grid.associated_edges_end( vertexOne );
2469 : itEdgOne++
2470 : )
2471 : {
2472 0 : Edge * edgOne = *itEdgOne;
2473 :
2474 0 : if( edgOne )
2475 : {
2476 :
2477 :
2478 0 : for( Grid::AssociatedEdgeIterator itEdgTwo = m_grid.associated_edges_begin( vertexTwo );
2479 0 : itEdgTwo != m_grid.associated_edges_end( vertexTwo );
2480 : itEdgTwo++
2481 : )
2482 : {
2483 0 : Edge * edgTwo = *itEdgTwo;
2484 :
2485 0 : if( edgTwo )
2486 : {
2487 0 : if( edgOne == edgTwo )
2488 : {
2489 0 : edgOfVrtcs = edgOne;
2490 : return true;
2491 : }
2492 : }
2493 : }
2494 : }
2495 : }
2496 :
2497 : return false;
2498 : }
2499 :
2500 : ///////////////////////////////////////////////////////////////////////////////////
2501 :
2502 :
2503 0 : bool DiamondsEstablish3D::splitThreeCrossLargeDiams( CombiNewVolsProps & combiNewVolsProps )
2504 : {
2505 : VrtxIndxCombi viCombiCenter;
2506 : combiNewVolsProps.spuckCenterVrtcsSudos(viCombiCenter);
2507 : Volume * vol;
2508 : combiNewVolsProps.spuckFulldimElem(vol);
2509 :
2510 : Edge * centersEdg;
2511 : combiNewVolsProps.spuckLowdimElem( centersEdg );
2512 :
2513 : Vertex * cutVrtx;
2514 0 : cutVrtx = m_attAccsDiamCentrVrtx[centersEdg];
2515 :
2516 : combiNewVolsProps.schluckNewSplitVrtx(cutVrtx);
2517 :
2518 : std::vector<Vertex*> shiftVrtcs;
2519 : combiNewVolsProps.spuckShiftVrtcs(shiftVrtcs);
2520 :
2521 0 : IndexType ind3CrossSide = combiNewVolsProps.spuckThreeCrossIndex();
2522 :
2523 : IndexType indOtherSide = ( ind3CrossSide + 1 ) % 2;
2524 :
2525 0 : VrtxIndxPair diamCenterVrtxPr( viCombiCenter[ind3CrossSide] );
2526 :
2527 0 : Vertex * centerDiam = diamCenterVrtxPr.first;
2528 :
2529 0 : VrtxIndxPair diamOtherVrtxPr( viCombiCenter[ indOtherSide ] );
2530 :
2531 0 : Vertex * centerOther = diamOtherVrtxPr.first;
2532 :
2533 : std::vector<Vertex*> midPtVrtcs;
2534 : combiNewVolsProps.spuckMidPtVrtcs(midPtVrtcs);
2535 :
2536 0 : Vertex * midPtDiam = midPtVrtcs[ind3CrossSide];
2537 0 : Vertex * midPtOther = midPtVrtcs[indOtherSide];
2538 :
2539 0 : Vertex * shiftVrtxDiam = shiftVrtcs[ind3CrossSide];
2540 0 : Vertex * shiftVrtxOther = shiftVrtcs[ indOtherSide];
2541 :
2542 0 : Volume * splitVolPrism = *m_grid.create<Prism>(
2543 0 : PrismDescriptor( cutVrtx, midPtDiam, shiftVrtxDiam,
2544 : centerOther, midPtOther, shiftVrtxOther
2545 : )
2546 0 : );
2547 :
2548 0 : assignSudoOfNewVols2VolAndSubElems(splitVolPrism, m_attAccsNewSudoOfVrtx[centerOther]);
2549 : // m_sh.assign_subset( splitVolPrism, m_attAccsNewSudoOfVrtx[centerOther] );
2550 :
2551 0 : Volume * splitVolTetra = *m_grid.create<Tetrahedron>(
2552 0 : TetrahedronDescriptor( shiftVrtxDiam, midPtDiam, cutVrtx, centerDiam ) );
2553 : //
2554 : //// m_sh.assign_subset( splitVolTetra, m_attAccsNewSudoOfVrtx[centerDiam] );
2555 0 : assignSudoOfNewVols2VolAndSubElems(splitVolTetra, m_attAccsNewSudoOfVrtx[centerDiam]);
2556 : //
2557 : // figure out the faces and the edges which later should be erased
2558 : //
2559 : // IndexType numEdgesFound = 0;
2560 : //
2561 : // Edge* centersEdg = nullptr;
2562 : //
2563 : // for( IndexType iEdg = 0; iEdg < vol->num_edges(); iEdg++)
2564 : // {
2565 : // Edge * edg = m_grid.get_edge(vol, iEdg);
2566 : //
2567 : // if( EdgeContains(edg, centerDiam ) && EdgeContains(edg, centerOther))
2568 : // {
2569 : // centersEdg = edg;
2570 : // numEdgesFound++;
2571 : // }
2572 : // }
2573 : //
2574 : // if( numEdgesFound != 1 )
2575 : // {
2576 : // UG_LOG("edges found at end " << numEdgesFound << std::endl);
2577 : // return false;
2578 : // }
2579 : //
2580 : // addElem(m_edges2BDeletedAtLastStep,centersEdg);
2581 : //
2582 : std::vector<Face*> faces2Del;
2583 :
2584 0 : for( IndexType iFac = 0; iFac < vol->num_faces(); iFac++)
2585 : {
2586 0 : Face * fac = m_grid.get_face(vol, iFac);
2587 :
2588 0 : if( FaceContains(fac, centersEdg) )
2589 : {
2590 0 : faces2Del.push_back(fac);
2591 : }
2592 : }
2593 :
2594 0 : if( faces2Del.size() != 2 )
2595 : {
2596 : UG_LOG("strange number of faces to del at fin " << faces2Del.size() << std::endl);
2597 : return false;
2598 : }
2599 :
2600 0 : for( Face * f : faces2Del )
2601 : {
2602 0 : addElem(m_faces2BDeletedAtLastStep, f);
2603 : }
2604 :
2605 0 : m_grid.erase(vol);
2606 :
2607 :
2608 : return true;
2609 0 : }
2610 :
2611 :
2612 : ///////////////////////////////////////////////////////////////////////////////////
2613 :
2614 :
2615 : } /* namespace diamonds */
2616 :
2617 : } /* namespace arte */
2618 :
2619 : } /* namespace ug */
|
