Line data Source code
1 : /*
2 : * Copyright (c) 2014-2021: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Stepniewski
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 :
34 : #include "subdivision_volumes.h"
35 :
36 : namespace ug
37 : {
38 :
39 :
40 : /// global boundary refinement rule variable for switching between linear and Subdivision Loop refinement
41 : static GlobalBoundaryRefinementRule g_boundaryRefinementRule = LINEAR;
42 :
43 : ////////////////////////////////////////////////////////////////////////////////
44 0 : void SetBoundaryRefinementRule(GlobalBoundaryRefinementRule refRule)
45 : {
46 0 : g_boundaryRefinementRule = refRule;
47 0 : }
48 :
49 :
50 : ////////////////////////////////////////////////////////////////////////////////
51 0 : GlobalBoundaryRefinementRule GetBoundaryRefinementRule()
52 : {
53 0 : return g_boundaryRefinementRule;
54 : }
55 :
56 :
57 : ////////////////////////////////////////////////////////////////////////////////
58 0 : void CheckValences(MultiGrid& mg, MGSubsetHandler& markSH, const char* filename)
59 : {
60 : // Output file SubsetHandler
61 0 : SubsetHandler sh(mg);
62 :
63 : // Attachments for associated number of elements
64 : AInt aEdgeValence;
65 : AInt aVertexValence_toc;
66 : AInt aVertexValence_prism;
67 : AInt aVertexValence_hex;
68 :
69 : // attach previously declared attachments with initial value 0
70 0 : mg.attach_to_edges_dv(aEdgeValence, 0);
71 0 : mg.attach_to_vertices_dv(aVertexValence_toc, 0);
72 0 : mg.attach_to_vertices_dv(aVertexValence_prism, 0);
73 0 : mg.attach_to_vertices_dv(aVertexValence_hex, 0);
74 :
75 : // Define attachment accessors
76 : Grid::EdgeAttachmentAccessor<AInt> aaEdgeValence(mg, aEdgeValence);
77 : Grid::VertexAttachmentAccessor<AInt> aaVertexValence_toc(mg, aVertexValence_toc);
78 : Grid::VertexAttachmentAccessor<AInt> aaVertexValence_prism(mg, aVertexValence_prism);
79 : Grid::VertexAttachmentAccessor<AInt> aaVertexValence_hex(mg, aVertexValence_hex);
80 :
81 : // Calculate edge valence
82 0 : for(VolumeIterator vIter = mg.begin<Volume>(mg.top_level()); vIter != mg.end<Volume>(mg.top_level()); ++vIter)
83 : {
84 : Volume* v = *vIter;
85 :
86 : Grid::edge_traits::secure_container edges;
87 : mg.associated_elements(edges, v);
88 :
89 0 : for(size_t i = 0; i < edges.size(); ++i)
90 : {
91 0 : if(markSH.get_subset_index(edges[i]) != -1)
92 0 : continue;
93 :
94 0 : ++aaEdgeValence[edges[i]];
95 : }
96 : }
97 :
98 : // Calculate vertex valence
99 0 : CalculateNumElemsVertexAttachmentInTopLevel(mg, aVertexValence_toc, aVertexValence_prism, aVertexValence_hex);
100 :
101 : // Assign subset indices by valence
102 0 : for(EdgeIterator eIter = mg.begin<Edge>(mg.top_level()); eIter != mg.end<Edge>(mg.top_level()); ++eIter)
103 : {
104 : Edge* e = *eIter;
105 :
106 0 : if(markSH.get_subset_index(e) != -1)
107 0 : continue;
108 :
109 0 : if(aaEdgeValence[e] < 4)
110 : {
111 0 : sh.assign_subset(e, aaEdgeValence[e]);
112 : }
113 : }
114 :
115 : // Assign subset indices by valence
116 0 : for(VertexIterator vIter = mg.begin<Vertex>(mg.top_level()); vIter != mg.end<Vertex>(mg.top_level()); ++vIter)
117 : {
118 : Vertex* v = *vIter;
119 :
120 0 : if(markSH.get_subset_index(v) != -1)
121 0 : continue;
122 :
123 0 : if(aaVertexValence_toc[v] < 4)
124 : {
125 0 : sh.assign_subset(v, aaVertexValence_toc[v]);
126 : }
127 : }
128 :
129 : // Output grid with valence SubsetHandler
130 0 : SaveGridToFile(mg, sh, filename);
131 :
132 : // Clean up
133 : mg.detach_from_edges(aEdgeValence);
134 : mg.detach_from_vertices(aVertexValence_toc);
135 : mg.detach_from_vertices(aVertexValence_prism);
136 : mg.detach_from_vertices(aVertexValence_hex);
137 0 : }
138 :
139 :
140 : ////////////////////////////////////////////////////////////////////////////////
141 0 : void PrintSubdivisionVolumesRefinementMask()
142 : {
143 : // size_t n = 3;
144 : //
145 : // double lin[n][n][n];
146 : // double m[2*n-1][2*n-1][2*n-1];
147 : //
148 : // for(size_t i = 0; i < n; ++i)
149 : // {
150 : // for(size_t j = 0; j < n; ++j)
151 : // {
152 : // for(size_t k = 0; k < n; ++k)
153 : // {
154 : // lin[i][j][k] = 0.0;
155 : // }
156 : // }
157 : // }
158 : //
159 : // for(size_t i = 0; i < 2*n-1; ++i)
160 : // {
161 : // for(size_t j = 0; j < 2*n-1; ++j)
162 : // {
163 : // for(size_t k = 0; k < 2*n-1; ++k)
164 : // {
165 : // m[i][j][k] = 0.0;
166 : // }
167 : // }
168 : // }
169 : //
170 : // lin[0][2][0] = 1.0;
171 : // lin[1][1][0] = 3.0;
172 : // lin[1][2][0] = 3.0;
173 : // lin[2][0][0] = 1.0;
174 : // lin[2][1][0] = 3.0;
175 : // lin[2][2][0] = 1.0;
176 : //
177 : // lin[0][1][1] = 3.0;
178 : // lin[0][2][1] = 3.0;
179 : // lin[1][0][1] = 3.0;
180 : // lin[1][1][1] = 6.0;
181 : // lin[1][2][1] = 3.0;
182 : // lin[2][0][1] = 3.0;
183 : // lin[2][1][1] = 3.0;
184 : //
185 : // lin[0][0][2] = 1.0;
186 : // lin[0][1][2] = 3.0;
187 : // lin[0][2][2] = 1.0;
188 : // lin[1][0][2] = 3.0;
189 : // lin[1][1][2] = 3.0;
190 : // lin[2][0][2] = 1.0;
191 : //
192 : // for(size_t i = 0; i < n; ++i)
193 : // {
194 : // for(size_t j = 0; j < n; ++j)
195 : // {
196 : // for(size_t k = 0; k < n; ++k)
197 : // {
198 : // for(size_t a = 0; a < n; ++a)
199 : // {
200 : // for(size_t b = 0; b < n; ++b)
201 : // {
202 : // for(size_t c = 0; c < n; ++c)
203 : // {
204 : // m[i+a][j+b][k+c] += lin[i][j][k]*lin[a][b][c];
205 : // }
206 : // }
207 : // }
208 : // }
209 : // }
210 : // }
211 : //
212 : // for(size_t i = 0; i < 2*n-1; ++i)
213 : // {
214 : // for(size_t k = 0; k < 2*n-1; ++k)
215 : // {
216 : // for(size_t j = 0; j < 2*n-1; ++j)
217 : // {
218 : // UG_LOG(m[i][j][k] << "\t ");
219 : // }
220 : //
221 : // UG_LOG(" | ");
222 : // }
223 : //
224 : // UG_LOG(std::endl);
225 : // }
226 :
227 : const size_t n = 3;
228 :
229 : double lin[n][n];
230 : double m[2*n-1][2*n-1];
231 :
232 0 : for(size_t i = 0; i < n; ++i)
233 : {
234 0 : for(size_t j = 0; j < n; ++j)
235 : {
236 0 : lin[i][j] = 0.0;
237 : }
238 : }
239 :
240 0 : for(size_t i = 0; i < 2*n-1; ++i)
241 : {
242 0 : for(size_t j = 0; j < 2*n-1; ++j)
243 : {
244 0 : m[i][j] = 0.0;
245 : }
246 : }
247 :
248 0 : lin[0][1] = 1.0;
249 0 : lin[0][2] = 1.0;
250 0 : lin[1][0] = 1.0;
251 0 : lin[1][1] = 2.0;
252 0 : lin[1][2] = 1.0;
253 0 : lin[2][0] = 1.0;
254 0 : lin[2][1] = 1.0;
255 :
256 0 : for(size_t i = 0; i < n; ++i)
257 : {
258 0 : for(size_t j = 0; j < n; ++j)
259 : {
260 0 : for(size_t k = 0; k < n; ++k)
261 : {
262 0 : for(size_t l = 0; l < n; ++l)
263 : {
264 0 : m[i+k][j+l] += lin[i][j]*lin[k][l];
265 : }
266 : }
267 : }
268 : }
269 :
270 0 : for(size_t i = 0; i < 2*n-1; ++i)
271 : {
272 0 : for(size_t j = 0; j < 2*n-1; ++j)
273 : {
274 0 : UG_LOG(m[i][j] << "\t");
275 : }
276 :
277 : UG_LOG(std::endl);
278 : }
279 0 : }
280 :
281 :
282 : ////////////////////////////////////////////////////////////////////////////////
283 0 : void SplitOctahedronToTetrahedrons( Grid& grid, Octahedron* oct, Volume* parentVol,
284 : std::vector<Tetrahedron*>& vTetsOut, int bestDiag)
285 : {
286 : // Position attachment management
287 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPosition);
288 :
289 : // Determine the shortest diagonal to split upon the octahedron
290 0 : if(bestDiag != 0 && bestDiag != 1 && bestDiag != 2)
291 : {
292 : bestDiag = 2;
293 :
294 0 : number d05 = VecDistanceSq(aaPos[oct->vertex(1)], aaPos[oct->vertex(3)]);
295 0 : number d13 = VecDistanceSq(aaPos[oct->vertex(0)], aaPos[oct->vertex(5)]);
296 0 : number d = VecDistanceSq(aaPos[oct->vertex(2)], aaPos[oct->vertex(4)]);
297 :
298 0 : if(d13 < d){
299 : bestDiag = 1;
300 : d = d13;
301 : }
302 0 : if(d05 < d){
303 : bestDiag = 0;
304 : }
305 : }
306 :
307 : Tetrahedron* tet1;
308 : Tetrahedron* tet2;
309 : Tetrahedron* tet3;
310 : Tetrahedron* tet4;
311 :
312 0 : switch(bestDiag){
313 :
314 0 : case 0:// diag: 0-5
315 : // Remark: element creation without father element specification
316 0 : tet1 = *grid.create<Tetrahedron>(TetrahedronDescriptor( oct->vertex(1),
317 0 : oct->vertex(0),
318 0 : oct->vertex(4),
319 0 : oct->vertex(3)), parentVol);
320 :
321 0 : tet2 = *grid.create<Tetrahedron>(TetrahedronDescriptor( oct->vertex(0),
322 0 : oct->vertex(2),
323 0 : oct->vertex(3),
324 0 : oct->vertex(1)), parentVol);
325 :
326 0 : tet3 = *grid.create<Tetrahedron>(TetrahedronDescriptor( oct->vertex(4),
327 0 : oct->vertex(3),
328 0 : oct->vertex(5),
329 0 : oct->vertex(1)), parentVol);
330 :
331 0 : tet4 = *grid.create<Tetrahedron>(TetrahedronDescriptor( oct->vertex(1),
332 0 : oct->vertex(5),
333 0 : oct->vertex(2),
334 0 : oct->vertex(3)), parentVol);
335 :
336 0 : vTetsOut.push_back(tet1);
337 0 : vTetsOut.push_back(tet2);
338 0 : vTetsOut.push_back(tet3);
339 0 : vTetsOut.push_back(tet4);
340 :
341 : // compare case 0 tetrahedron pattern from tetrahedron_rules.cpp:
342 : /*
343 : inds[fi++] = 4;
344 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 1;
345 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 5;
346 :
347 : inds[fi++] = 4;
348 : inds[fi++] = NUM_VERTICES + 1; inds[fi++] = NUM_VERTICES + 4;
349 : inds[fi++] = NUM_VERTICES + 5; inds[fi++] = NUM_VERTICES + 0;
350 :
351 : inds[fi++] = 4;
352 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 5;
353 : inds[fi++] = NUM_VERTICES + 3; inds[fi++] = NUM_VERTICES + 0;
354 :
355 : inds[fi++] = 4;
356 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 3;
357 : inds[fi++] = NUM_VERTICES + 4; inds[fi++] = NUM_VERTICES + 5;
358 : */
359 :
360 : break;
361 :
362 0 : case 1:// diag: 1-3
363 0 : tet1 = *grid.create<Tetrahedron>(TetrahedronDescriptor( oct->vertex(1),
364 0 : oct->vertex(0),
365 0 : oct->vertex(4),
366 0 : oct->vertex(5)), parentVol);
367 :
368 0 : tet2 = *grid.create<Tetrahedron>(TetrahedronDescriptor( oct->vertex(0),
369 0 : oct->vertex(2),
370 0 : oct->vertex(3),
371 0 : oct->vertex(5)), parentVol);
372 :
373 0 : tet3 = *grid.create<Tetrahedron>(TetrahedronDescriptor( oct->vertex(4),
374 0 : oct->vertex(3),
375 0 : oct->vertex(5),
376 0 : oct->vertex(0)), parentVol);
377 :
378 0 : tet4 = *grid.create<Tetrahedron>(TetrahedronDescriptor( oct->vertex(1),
379 0 : oct->vertex(5),
380 0 : oct->vertex(2),
381 0 : oct->vertex(0)), parentVol);
382 :
383 0 : vTetsOut.push_back(tet1);
384 0 : vTetsOut.push_back(tet2);
385 0 : vTetsOut.push_back(tet3);
386 0 : vTetsOut.push_back(tet4);
387 :
388 : // compare case 1 tetrahedron pattern from tetrahedron_rules.cpp:
389 : /*
390 : inds[fi++] = 4;
391 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 1;
392 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 3;
393 :
394 : inds[fi++] = 4;
395 : inds[fi++] = NUM_VERTICES + 1; inds[fi++] = NUM_VERTICES + 4;
396 : inds[fi++] = NUM_VERTICES + 5; inds[fi++] = NUM_VERTICES + 3;
397 :
398 : inds[fi++] = 4;
399 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 5;
400 : inds[fi++] = NUM_VERTICES + 3; inds[fi++] = NUM_VERTICES + 1;
401 :
402 : inds[fi++] = 4;
403 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 3;
404 : inds[fi++] = NUM_VERTICES + 4; inds[fi++] = NUM_VERTICES + 1;
405 : */
406 :
407 : break;
408 :
409 0 : case 2:// diag 2-4
410 0 : tet1 = *grid.create<Tetrahedron>(TetrahedronDescriptor( oct->vertex(1),
411 0 : oct->vertex(4),
412 0 : oct->vertex(5),
413 0 : oct->vertex(2)), parentVol);
414 :
415 0 : tet2 = *grid.create<Tetrahedron>(TetrahedronDescriptor( oct->vertex(0),
416 0 : oct->vertex(4),
417 0 : oct->vertex(1),
418 0 : oct->vertex(2)), parentVol);
419 :
420 0 : tet3 = *grid.create<Tetrahedron>(TetrahedronDescriptor( oct->vertex(4),
421 0 : oct->vertex(5),
422 0 : oct->vertex(2),
423 0 : oct->vertex(3)), parentVol);
424 :
425 0 : tet4 = *grid.create<Tetrahedron>(TetrahedronDescriptor( oct->vertex(2),
426 0 : oct->vertex(0),
427 0 : oct->vertex(4),
428 0 : oct->vertex(3)), parentVol);
429 :
430 0 : vTetsOut.push_back(tet1);
431 0 : vTetsOut.push_back(tet2);
432 0 : vTetsOut.push_back(tet3);
433 0 : vTetsOut.push_back(tet4);
434 :
435 : // compare case 2 tetrahedron pattern from tetrahedron_rules.cpp:
436 : /*
437 : inds[fi++] = 4;
438 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 2;
439 : inds[fi++] = NUM_VERTICES + 3; inds[fi++] = NUM_VERTICES + 4;
440 :
441 : inds[fi++] = 4;
442 : inds[fi++] = NUM_VERTICES + 1; inds[fi++] = NUM_VERTICES + 2;
443 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 4;
444 :
445 : inds[fi++] = 4;
446 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 3;
447 : inds[fi++] = NUM_VERTICES + 4; inds[fi++] = NUM_VERTICES + 5;
448 :
449 : inds[fi++] = 4;
450 : inds[fi++] = NUM_VERTICES + 4; inds[fi++] = NUM_VERTICES + 1;
451 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 5;
452 : */
453 :
454 : break;
455 : }
456 0 : }
457 :
458 :
459 : ////////////////////////////////////////////////////////////////////////////////
460 0 : void TetrahedralizeHybridTetOctGrid(MultiGrid& mg, int bestDiag)
461 : {
462 : PROFILE_FUNC_GROUP("subdivision_volumes");
463 :
464 : #ifdef UG_PARALLEL
465 : DistributedGridManager* dgm = mg.distributed_grid_manager();
466 : #endif
467 :
468 0 : if(bestDiag != 0 && bestDiag != 1 && bestDiag != 2)
469 : {
470 : bestDiag = -1;
471 : }
472 :
473 : // Position attachment management
474 : Grid::VertexAttachmentAccessor<APosition> aaPos(mg, aPosition);
475 :
476 : std::vector<Tetrahedron*> vTetsOut;
477 :
478 : // Loop over all levels and split octahedrons
479 0 : for(size_t i = mg.num_levels()-1; i > 0; --i)
480 : {
481 : // Loop over all octahedrons on each level
482 0 : for(VolumeIterator octIter = mg.begin<Octahedron>(i); octIter != mg.end<Octahedron>(i); ++octIter)
483 : {
484 0 : Octahedron* oct = dynamic_cast<Octahedron*>(*octIter);
485 0 : Volume* parentVol = dynamic_cast<Volume*>(mg.get_parent(oct));
486 :
487 : #ifdef UG_PARALLEL
488 : /*
489 : * In case of parallel distribution e.g. in level 2
490 : * there can be subsequently be elements which
491 : * locally don’t have a parent (esp. V_SLAVES),
492 : * i.e. parentVol = NULL. And elements with
493 : * parent = NULL are associated to mg.level = 0.
494 : * Therefore, if(dgm->is_ghost(oct)) is not
495 : * sufficient.
496 : */
497 : if(dgm->contains_status(oct, ES_V_SLAVE))
498 : continue;
499 : #endif
500 :
501 0 : SplitOctahedronToTetrahedrons(mg, oct, parentVol, vTetsOut, bestDiag);
502 : }
503 : }
504 :
505 : // Erase all octahedrons in multigrid
506 0 : while(mg.begin<Octahedron>() != mg.end<Octahedron>())
507 : {
508 : Octahedron* oct = *mg.begin<Octahedron>();
509 0 : mg.erase(oct);
510 : }
511 0 : }
512 :
513 :
514 : ////////////////////////////////////////////////////////////////////////////////
515 : template <class TAPosition>
516 0 : void ProjectHierarchyToSubdivisionLimit(MultiGrid& mg, TAPosition& aPos)
517 : {
518 : PROFILE_FUNC_GROUP("subdivision_volumes");
519 :
520 : if(TAPosition::ValueType::Size == 1){
521 0 : UG_THROW("Error in ProjectHierarchyToSubdivisionLimit:\n"
522 : "Currently only dimensions 2 and 3 are supported.\n");
523 : }
524 :
525 : // Catch use of procedure for MultiGrids with just one level
526 : size_t numLevels = mg.num_levels();
527 :
528 : #ifdef UG_PARALLEL
529 : if(pcl::NumProcs() > 1){
530 : pcl::ProcessCommunicator pc;
531 : numLevels = pc.allreduce(numLevels, PCL_RO_MAX);
532 : }
533 : #endif
534 :
535 0 : if(numLevels == 1)
536 : {
537 0 : UG_THROW("Error in ProjectHierarchyToSubdivisionLimit:\n"
538 : "Procedure only to be used for MultiGrids with more than one level.");
539 : }
540 :
541 : #ifdef UG_PARALLEL
542 : // Attachment communication policies COPY
543 : ComPol_CopyAttachment<VertexLayout, TAPosition> comPolCopyAPosition(mg, aPos);
544 :
545 : // Interface communicators and distributed domain manager
546 : pcl::InterfaceCommunicator<VertexLayout> com;
547 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
548 : GridLayoutMap& glm = dgm.grid_layout_map();
549 : #endif
550 :
551 : Grid::VertexAttachmentAccessor<TAPosition> aaPos(mg, aPos);
552 :
553 : // AttachmentCopy VSLAVE->VMASTER on mg.top_level() in case not all VMasters in toplevel don't have correct position already
554 : #ifdef UG_PARALLEL
555 : // copy v_slaves to ghosts = VMASTER
556 : com.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, comPolCopyAPosition);
557 : com.communicate();
558 : #endif
559 :
560 : // Loop all levels from toplevel down to base level
561 0 : for(int lvl = (int)mg.top_level(); lvl > 0; --lvl)
562 : {
563 : // Loop all vertices of current level and submit positions to parent vertices
564 0 : for(VertexIterator vrtIter = mg.begin<Vertex>(lvl); vrtIter != mg.end<Vertex>(lvl); ++vrtIter)
565 : {
566 : Vertex* v = *vrtIter;
567 0 : Vertex* parent = dynamic_cast<Vertex*>(mg.get_parent(v));
568 :
569 : // Only, if parent vertex exists
570 0 : if(parent)
571 : {
572 : aaPos[parent] = aaPos[v];
573 : }
574 : }
575 :
576 : #ifdef UG_PARALLEL
577 : // copy v_slaves to ghosts = VMASTER
578 : com.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, comPolCopyAPosition);
579 : com.communicate();
580 : #endif
581 : }
582 0 : }
583 :
584 :
585 : ////////////////////////////////////////////////////////////////////////////////
586 : template <class TAPosition>
587 0 : void CalculateSmoothCreaseManifoldPosInParentLevelLoopScheme(MultiGrid& mg, TAPosition& aPos, MGSubsetHandler& markSH,
588 : MGSubsetHandler& linearManifoldSH,
589 : TAPosition& aSmoothBndPosEvenVrt,
590 : TAPosition& aSmoothBndPosOddVrt,
591 : AInt& aNumManifoldEdges)
592 : {
593 : /*
594 : * Scheme reference:
595 : *
596 : * C. Loop, Smooth subdivision surfaces based on triangles,
597 : * master’s thesis, University of Utah, 1987.
598 : */
599 :
600 : if(TAPosition::ValueType::Size == 1){
601 : UG_THROW("Error in CalculateSmoothCreaseManifoldPosInParentLevelLoopScheme:\n"
602 : "Currently only dimensions 2 and 3 are supported.\n");
603 : }
604 :
605 : // Catch use of procedure for MultiGrids with just one level
606 : size_t numLevels = mg.num_levels();
607 :
608 : #ifdef UG_PARALLEL
609 : if(pcl::NumProcs() > 1){
610 : pcl::ProcessCommunicator pc;
611 : numLevels = pc.allreduce(numLevels, PCL_RO_MAX);
612 : }
613 : #endif
614 :
615 0 : if(numLevels == 1)
616 : {
617 0 : UG_THROW("Error in CalculateSmoothCreaseManifoldPosInParentLevelLoopScheme: "
618 : "Procedure only to be used for MultiGrids with more than one level.");
619 : }
620 :
621 : // Define attachment accessors
622 : Grid::VertexAttachmentAccessor<TAPosition> aaPos(mg, aPos);
623 : Grid::VertexAttachmentAccessor<TAPosition> aaSmoothBndPosEvenVrt(mg, aSmoothBndPosEvenVrt);
624 : Grid::EdgeAttachmentAccessor<TAPosition> aaSmoothBndPosOddVrt(mg, aSmoothBndPosOddVrt);
625 : Grid::VertexAttachmentAccessor<AInt> aaNumManifoldEdges(mg, aNumManifoldEdges);
626 :
627 : #ifdef UG_PARALLEL
628 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
629 : #endif
630 :
631 : // Declare centroid coordinate vector
632 : typedef typename TAPosition::ValueType position_type;
633 : position_type p;
634 : VecSet(p, 0);
635 :
636 : // Load subdivision surfaces rules
637 0 : SubdivRules_PLoop& subdiv = SubdivRules_PLoop::inst();
638 :
639 : /*
640 : * Smoothing of even vertices x
641 : *
642 : * Weights of center vertex : (1 - k*b[k]), k = vertex valence
643 : * Weights of edge adjacent vertices: b[k] = (0.625 - a*a)/k,
644 : * a = 0.375 + 0.25 * cos((2.0*PI)/k)
645 : *
646 : * Closed surface weights Polyline weights
647 : * for non-closed surfaces
648 : *
649 : (b) (b) (b)
650 : \ | /
651 : ...
652 : (b) -- (1-kb) -- (b) 0.125 -- 0.75 -- 0.125
653 :
654 : / | \
655 : (b) (b) (b)
656 : *
657 : */
658 :
659 : // Calculate smooth position for EVEN vertices
660 0 : for(VertexIterator vrtIter = mg.begin<Vertex>(mg.top_level()-1); vrtIter != mg.end<Vertex>(mg.top_level()-1); ++vrtIter)
661 : {
662 : VecSet(p, 0);
663 :
664 : Vertex* vrt = *vrtIter;
665 :
666 : // Skip ghost vertices
667 : #ifdef UG_PARALLEL
668 : if(dgm.is_ghost(vrt))
669 : continue;
670 : #endif
671 :
672 : // In case of vertices, which do not belong to the user-specified linear boundary manifold subsets,
673 : // and activated subdivision Loop refinement calculate subdivision surfaces smooth position
674 0 : if(linearManifoldSH.get_subset_index(vrt) == -1)
675 : {
676 : // perform loop subdivision on even manifold vertices
677 : // first get neighbored manifold vertices
678 0 : for(Grid::AssociatedEdgeIterator eIter = mg.associated_edges_begin(vrt); eIter != mg.associated_edges_end(vrt); ++eIter)
679 : {
680 0 : Edge* e = *eIter;
681 :
682 : // In case of a polyline surface boundary vertex (marked in non-closed surfaces) exclude non-boundary neighbor vertices
683 0 : if(markSH.get_subset_index(vrt) != -1){
684 0 : if(markSH.get_subset_index(e) == -1)
685 0 : continue;
686 : }
687 :
688 : // Exclude ghost and horizontal slave neighbor vertices from contributing to centroid
689 : #ifdef UG_PARALLEL
690 : if(dgm.is_ghost(e))
691 : continue;
692 :
693 : if(dgm.contains_status(e, ES_H_SLAVE))
694 : continue;
695 : #endif
696 :
697 0 : VecAdd(p, p, aaPos[GetConnectedVertex(e, vrt)]);
698 : }
699 :
700 : number centerWght;
701 : number nbrWght;
702 :
703 : // Polyline surface boundary vertex (marked in non-closed surfaces)
704 0 : if(markSH.get_subset_index(vrt) != -1){
705 : // Polyline weights in case of boundary vertices of non-closed surfaces
706 :
707 : // Prototype Butterfly interpolation
708 : // centerWght = 1.0;
709 : // nbrWght = 0.0;
710 :
711 : centerWght = 0.75;
712 : nbrWght = 0.125;
713 : }
714 : // Inner vertex
715 : else{
716 : // Regular subdivision weights by C. Loop for closed triangular surfaces
717 0 : centerWght = subdiv.ref_even_inner_center_weight(aaNumManifoldEdges[vrt]);
718 0 : nbrWght = subdiv.ref_even_inner_nbr_weight(aaNumManifoldEdges[vrt]);
719 : }
720 :
721 : // Exclude horizontal slaves of the currently smoothed vertex to avoid multiple contributions to centroid
722 : #ifdef UG_PARALLEL
723 : if(dgm.is_ghost(vrt))
724 : continue;
725 :
726 : if(dgm.contains_status(vrt, ES_H_SLAVE))
727 : {
728 : VecScaleAppend(aaSmoothBndPosEvenVrt[vrt], nbrWght, p);
729 : continue;
730 : }
731 : #endif
732 :
733 : VecScaleAppend(aaSmoothBndPosEvenVrt[vrt], centerWght, aaPos[vrt], nbrWght, p);
734 : }
735 : }
736 :
737 : /*
738 : * Smoothing of odd vertices x
739 : *
740 : * Weights of face adjacent vertices to x: 1/8
741 : * Weights of edge adjacent vertices to x: 3/8
742 : *
743 : * Closed surface weights Polyline weights
744 : * for non-closed surfaces
745 : 1/8
746 : / \
747 : 3/8--x--3/8 0.5--x--0.5
748 : \ /
749 : 1/8
750 : *
751 : */
752 :
753 : // Calculate smooth position for ODD vertices
754 0 : for(EdgeIterator eIter = mg.begin<Edge>(mg.top_level()-1); eIter != mg.end<Edge>(mg.top_level()-1); ++eIter)
755 : {
756 : VecSet(p, 0);
757 :
758 : Edge* e = *eIter;
759 :
760 : // Skip ghost edges
761 : #ifdef UG_PARALLEL
762 : if(dgm.is_ghost(e))
763 : continue;
764 : #endif
765 :
766 : // In case of manifold edges, which do not belong to the user-specified linear boundary manifold subsets,
767 : // and activated subdivision Loop refinement calculate subdivision surfaces smooth position
768 0 : if(linearManifoldSH.get_subset_index(e) == -1)
769 : {
770 : // In case of non-closed surfaces check if the current edge belongs to the polyline surface boundary (marked in non-closed surfaces)
771 0 : if(markSH.get_subset_index(e) != -1)
772 : {
773 : // Exclude ghost and horizontal slaves of the parent edge vertices of the currently smoothed vertex
774 : // to avoid multiple contributions to the centroid of the edge adjacent vertices
775 : #ifdef UG_PARALLEL
776 : if(dgm.is_ghost(e) || dgm.contains_status(e, ES_H_SLAVE))
777 : continue;
778 : #endif
779 :
780 : // Prototype Butterfly interpolation
781 : // std::vector<Edge*> vEdges;
782 : // CollectEdges(vEdges, mg, e->vertex(0), true);
783 : // for(size_t i = 0; i < vEdges.size(); ++i){
784 : // Edge* nbrEdge = vEdges[i];
785 : // if(nbrEdge != e && aaNumManifoldFaces[nbrEdge] == 1){
786 : // VecAdd(p, p, aaPos[GetConnectedVertex(nbrEdge, e->vertex(0))]);
787 : // }
788 : // }
789 : //
790 : // CollectEdges(vEdges, mg, e->vertex(1), true);
791 : // for(size_t i = 0; i < vEdges.size(); ++i){
792 : // Edge* nbrEdge = vEdges[i];
793 : // if(nbrEdge != e && aaNumManifoldFaces[nbrEdge] == 1){
794 : // VecAdd(p, p, aaPos[GetConnectedVertex(nbrEdge, e->vertex(1))]);
795 : // }
796 : // }
797 : //
798 : // VecScaleAppend(aaSmoothBndPosOddVrt[e], 9.0/16, aaPos[e->vertex(0)], 9.0/16, aaPos[e->vertex(1)], -1.0/16, p);
799 :
800 0 : VecScaleAppend(aaSmoothBndPosOddVrt[e], 0.5, aaPos[e->vertex(0)], 0.5, aaPos[e->vertex(1)]);
801 : }
802 : else
803 : {
804 : // perform loop subdivision on odd manifold vertices
805 : // get the neighbored manifold triangles
806 : std::vector<Face*> associatedFaces;
807 : std::vector<Face*> associatedManifoldFaces;
808 :
809 0 : CollectAssociated(associatedFaces, mg, e);
810 :
811 0 : for(size_t i = 0; i < associatedFaces.size(); ++i)
812 : {
813 : // Exclude ghost and horizontal slave manifold faces
814 : #ifdef UG_PARALLEL
815 : if(dgm.is_ghost(associatedFaces[i]))
816 : continue;
817 :
818 : if(dgm.contains_status(associatedFaces[i], ES_H_SLAVE))
819 : continue;
820 : #endif
821 :
822 0 : associatedManifoldFaces.push_back(associatedFaces[i]);
823 : }
824 :
825 0 : if(associatedManifoldFaces.size() <= 2)
826 : {
827 : // Check, if all faces are triangles
828 0 : for(size_t i = 0; i < associatedManifoldFaces.size(); ++i)
829 : {
830 0 : if(associatedManifoldFaces[i]->num_vertices() != 3)
831 : {
832 0 : UG_THROW("ERROR in CalculateSmoothCreaseManifoldPosInParentLevelLoopScheme: "
833 : "Non triangular faces included in grid: " << ElementDebugInfo(mg, associatedManifoldFaces[i]));
834 : }
835 : }
836 :
837 : // Summate centroid of face adjacent vertices
838 0 : for(size_t i = 0; i < associatedManifoldFaces.size(); ++i)
839 : {
840 0 : VecAdd(p, p, aaPos[GetConnectedVertex(e, associatedManifoldFaces[i])]);
841 : }
842 :
843 : // Exclude ghost and horizontal slaves of the parent edge vertices of the currently smoothed vertex
844 : // to avoid multiple contributions to the centroid of the edge adjacent vertices
845 : #ifdef UG_PARALLEL
846 : if(dgm.is_ghost(e))
847 : continue;
848 :
849 : if(dgm.contains_status(e, ES_H_SLAVE))
850 : {
851 : VecScaleAppend(aaSmoothBndPosOddVrt[e], 0.125, p);
852 : continue;
853 : }
854 : #endif
855 :
856 0 : VecScaleAppend(aaSmoothBndPosOddVrt[e], 0.375, aaPos[e->vertex(0)], 0.375, aaPos[e->vertex(1)], 0.125, p);
857 : }
858 : else
859 0 : UG_THROW("ERROR in CalculateSmoothCreaseManifoldPosInParentLevelLoopScheme: numAssociatedManifoldFaces > 2.");
860 0 : }
861 : }
862 : }
863 :
864 : // Manage vertex and edge attachment communication in parallel case -> COMMUNICATE aSmoothBndPosEvenVrt, aSmoothBndPosOddVrt
865 : #ifdef UG_PARALLEL
866 : // Reduce add operations:
867 : // sum up h_slaves into h_masters
868 :
869 : // Copy operations:
870 : // copy h_masters to h_slaves for consistency
871 : AttachmentAllReduce<Vertex>(mg, aSmoothBndPosEvenVrt, PCL_RO_SUM);
872 : AttachmentAllReduce<Edge>(mg, aSmoothBndPosOddVrt, PCL_RO_SUM);
873 : #endif
874 0 : }
875 :
876 :
877 : ////////////////////////////////////////////////////////////////////////////////
878 : template <class TAPosition>
879 0 : void CalculateSmoothManifoldPosInParentLevelLoopScheme(MultiGrid& mg, TAPosition& aPos, MGSubsetHandler& markSH,
880 : MGSubsetHandler& linearManifoldSH,
881 : TAPosition& aSmoothBndPosEvenVrt,
882 : TAPosition& aSmoothBndPosOddVrt,
883 : AInt& aNumManifoldEdges)
884 : {
885 : /*
886 : * Scheme reference:
887 : *
888 : * C. Loop, Smooth subdivision surfaces based on triangles,
889 : * master’s thesis, University of Utah, 1987.
890 : */
891 :
892 : if(TAPosition::ValueType::Size == 1){
893 : UG_THROW("Error in CalculateSmoothManifoldPosInParentLevelLoopScheme:\n"
894 : "Currently only dimensions 2 and 3 are supported.\n");
895 : }
896 :
897 : // Catch use of procedure for MultiGrids with just one level
898 : size_t numLevels = mg.num_levels();
899 :
900 : #ifdef UG_PARALLEL
901 : if(pcl::NumProcs() > 1){
902 : pcl::ProcessCommunicator pc;
903 : numLevels = pc.allreduce(numLevels, PCL_RO_MAX);
904 : }
905 : #endif
906 :
907 0 : if(numLevels == 1)
908 : {
909 0 : UG_THROW("Error in CalculateSmoothManifoldPosInParentLevelLoopScheme: "
910 : "Procedure only to be used for MultiGrids with more than one level.");
911 : }
912 :
913 : // Define attachment accessors
914 : Grid::VertexAttachmentAccessor<TAPosition> aaPos(mg, aPos);
915 : Grid::VertexAttachmentAccessor<TAPosition> aaSmoothBndPosEvenVrt(mg, aSmoothBndPosEvenVrt);
916 : Grid::EdgeAttachmentAccessor<TAPosition> aaSmoothBndPosOddVrt(mg, aSmoothBndPosOddVrt);
917 : Grid::VertexAttachmentAccessor<AInt> aaNumManifoldEdges(mg, aNumManifoldEdges);
918 :
919 : #ifdef UG_PARALLEL
920 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
921 : #endif
922 :
923 : // Declare centroid coordinate vector
924 : typedef typename TAPosition::ValueType position_type;
925 : position_type p;
926 : VecSet(p, 0);
927 :
928 : // Load subdivision surfaces rules
929 0 : SubdivRules_PLoop& subdiv = SubdivRules_PLoop::inst();
930 :
931 : // Calculate smooth position for EVEN vertices
932 0 : for(VertexIterator vrtIter = mg.begin<Vertex>(mg.top_level()-1); vrtIter != mg.end<Vertex>(mg.top_level()-1); ++vrtIter)
933 : {
934 : VecSet(p, 0);
935 :
936 : Vertex* vrt = *vrtIter;
937 :
938 : // Skip ghost vertices
939 : #ifdef UG_PARALLEL
940 : if(dgm.is_ghost(vrt))
941 : continue;
942 : #endif
943 :
944 : // In case of marked manifold vertices, which do not belong to the user-specified linear boundary manifold subsets,
945 : // and activated subdivision Loop refinement calculate subdivision surfaces smooth position
946 0 : if(markSH.get_subset_index(vrt) != -1 && linearManifoldSH.get_subset_index(vrt) == -1)
947 : {
948 : // perform loop subdivision on even manifold vertices
949 : // first get neighbored manifold vertices
950 0 : for(Grid::AssociatedEdgeIterator eIter = mg.associated_edges_begin(vrt); eIter != mg.associated_edges_end(vrt); ++eIter)
951 : {
952 0 : Edge* e = *eIter;
953 :
954 : // Only consider associated edges, which are marked as manifold edges
955 0 : if(markSH.get_subset_index(e) != -1)
956 : {
957 : // Exclude ghost and horizontal slave neighbor vertices from contributing to centroid
958 : #ifdef UG_PARALLEL
959 : if(dgm.is_ghost(e))
960 : continue;
961 :
962 : if(dgm.contains_status(e, ES_H_SLAVE))
963 : continue;
964 : #endif
965 :
966 0 : VecAdd(p, p, aaPos[GetConnectedVertex(e, vrt)]);
967 : }
968 : }
969 :
970 0 : number centerWgt = subdiv.ref_even_inner_center_weight(aaNumManifoldEdges[vrt]);
971 0 : number nbrWgt = subdiv.ref_even_inner_nbr_weight(aaNumManifoldEdges[vrt]);
972 :
973 : // Exclude horizontal slaves of the currently smoothed vertex to avoid multiple contributions to centroid
974 : #ifdef UG_PARALLEL
975 : if(dgm.is_ghost(vrt))
976 : continue;
977 :
978 : if(dgm.contains_status(vrt, ES_H_SLAVE))
979 : {
980 : VecScaleAppend(aaSmoothBndPosEvenVrt[vrt], nbrWgt, p);
981 : continue;
982 : }
983 : #endif
984 :
985 : VecScaleAppend(aaSmoothBndPosEvenVrt[vrt], centerWgt, aaPos[vrt], nbrWgt, p);
986 : }
987 : }
988 :
989 : /*
990 : * Smoothing of odd vertices x
991 : *
992 : * Weights of face adjacent vertices to x: 1/8
993 : * Weights of edge adjacent vertices to x: 3/8
994 : *
995 : 1/8
996 : / \
997 : 3/8--x--3/8
998 : \ /
999 : 1/8
1000 : *
1001 : */
1002 :
1003 : // Calculate smooth position for ODD vertices
1004 0 : for(EdgeIterator eIter = mg.begin<Edge>(mg.top_level()-1); eIter != mg.end<Edge>(mg.top_level()-1); ++eIter)
1005 : {
1006 : VecSet(p, 0);
1007 :
1008 : Edge* e = *eIter;
1009 :
1010 : // Skip ghost edges
1011 : #ifdef UG_PARALLEL
1012 : if(dgm.is_ghost(e))
1013 : continue;
1014 : #endif
1015 :
1016 : // In case of marked manifold edges, which do not belong to the user-specified linear boundary manifold subsets,
1017 : // and activated subdivision Loop refinement calculate subdivision surfaces smooth position
1018 0 : if(markSH.get_subset_index(e) != -1 && linearManifoldSH.get_subset_index(e) == -1)
1019 : {
1020 : // perform loop subdivision on odd manifold vertices
1021 : // get the neighbored manifold triangles
1022 : std::vector<Face*> associatedFaces;
1023 : std::vector<Face*> associatedManifoldFaces;
1024 :
1025 0 : CollectAssociated(associatedFaces, mg, e);
1026 :
1027 0 : for(size_t i = 0; i < associatedFaces.size(); ++i)
1028 : {
1029 :
1030 : // Only consider associated faces, which are marked as manifold faces
1031 0 : if(markSH.get_subset_index(associatedFaces[i]) != -1)
1032 : {
1033 : // Exclude ghost and horizontal slave manifold faces
1034 : #ifdef UG_PARALLEL
1035 : if(dgm.is_ghost(associatedFaces[i]))
1036 : continue;
1037 :
1038 : if(dgm.contains_status(associatedFaces[i], ES_H_SLAVE))
1039 : continue;
1040 : #endif
1041 :
1042 0 : associatedManifoldFaces.push_back(associatedFaces[i]);
1043 : }
1044 : }
1045 :
1046 0 : if(associatedManifoldFaces.size() <= 2)
1047 : {
1048 : // Check, if all faces are triangles
1049 0 : for(size_t i = 0; i < associatedManifoldFaces.size(); ++i)
1050 : {
1051 0 : if(associatedManifoldFaces[i]->num_vertices() != 3)
1052 : {
1053 0 : UG_THROW("ERROR in CalculateSmoothManifoldPosInParentLevelLoopScheme: "
1054 : "Non triangular faces included in grid: " << ElementDebugInfo(mg, associatedManifoldFaces[i]));
1055 : }
1056 : }
1057 :
1058 : // Summate centroid of face adjacent vertices
1059 0 : for(size_t i = 0; i < associatedManifoldFaces.size(); ++i)
1060 : {
1061 0 : VecAdd(p, p, aaPos[GetConnectedVertex(e, associatedManifoldFaces[i])]);
1062 : }
1063 :
1064 : // Exclude ghost and horizontal slaves of the parent edge vertices of the currently smoothed vertex
1065 : // to avoid multiple contributions to the centroid of the edge adjacent vertices
1066 : #ifdef UG_PARALLEL
1067 : if(dgm.is_ghost(e))
1068 : {
1069 : continue;
1070 : }
1071 :
1072 : if(dgm.contains_status(e, ES_H_SLAVE))
1073 : {
1074 : VecScaleAppend(aaSmoothBndPosOddVrt[e], 0.125, p);
1075 : continue;
1076 : }
1077 : #endif
1078 :
1079 0 : VecScaleAppend(aaSmoothBndPosOddVrt[e], 0.375, aaPos[e->vertex(0)], 0.375, aaPos[e->vertex(1)], 0.125, p);
1080 : }
1081 : else
1082 0 : UG_THROW("ERROR in CalculateSmoothManifoldPosInParentLevelLoopScheme: numAssociatedManifoldFaces > 2.");
1083 0 : }
1084 : }
1085 :
1086 : // Manage vertex and edge attachment communication in parallel case -> COMMUNICATE aSmoothBndPosEvenVrt, aSmoothBndPosOddVrt
1087 : #ifdef UG_PARALLEL
1088 : // Reduce add operations:
1089 : // sum up h_slaves into h_masters
1090 :
1091 : // Copy operations:
1092 : // copy h_masters to h_slaves for consistency
1093 : AttachmentAllReduce<Vertex>(mg, aSmoothBndPosEvenVrt, PCL_RO_SUM);
1094 : AttachmentAllReduce<Edge>(mg, aSmoothBndPosOddVrt, PCL_RO_SUM);
1095 : #endif
1096 0 : }
1097 :
1098 :
1099 : ////////////////////////////////////////////////////////////////////////////////
1100 : template <class TAPosition>
1101 0 : void CalculateSmoothManifoldPosInTopLevelAveragingScheme(MultiGrid& mg, TAPosition& aPos, MGSubsetHandler& markSH,
1102 : MGSubsetHandler& linearManifoldSH,
1103 : TAPosition& aSmoothBndPos_tri,
1104 : TAPosition& aSmoothBndPos_quad)
1105 : {
1106 : /*
1107 : * Scheme reference:
1108 : *
1109 : * J. Warren and H. Weimer, Subdivision Methods for Geometric Design: A Constructive Approach,
1110 : * Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1st ed., 2001.
1111 : */
1112 :
1113 : if(TAPosition::ValueType::Size == 1){
1114 : UG_THROW("Error in CalculateSmoothManifoldPosInTopLevelAveragingScheme:\n"
1115 : "Currently only dimensions 2 and 3 are supported.\n");
1116 : }
1117 :
1118 : // Define attachment accessors
1119 : Grid::VertexAttachmentAccessor<TAPosition> aaPos(mg, aPos);
1120 : Grid::VertexAttachmentAccessor<TAPosition> aaSmoothBndPos_tri(mg, aSmoothBndPos_tri);
1121 : Grid::VertexAttachmentAccessor<TAPosition> aaSmoothBndPos_quad(mg, aSmoothBndPos_quad);
1122 :
1123 : #ifdef UG_PARALLEL
1124 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
1125 : #endif
1126 :
1127 : // Declare centroid coordinate vector
1128 : typedef typename TAPosition::ValueType position_type;
1129 : position_type p;
1130 : VecSet(p, 0);
1131 :
1132 : // Loop all manifold faces of top_level
1133 0 : for(FaceIterator fIter = mg.begin<Face>(mg.top_level()); fIter != mg.end<Face>(mg.top_level()); ++fIter)
1134 : {
1135 : Face* f = *fIter;
1136 :
1137 : // Skip ghost volumes
1138 : #ifdef UG_PARALLEL
1139 : if(dgm.is_ghost(f))
1140 : continue;
1141 : #endif
1142 :
1143 : // In case of marked manifold faces, which do not belong to the user-specified linear boundary manifold subsets,
1144 : // and activated Averaging scheme calculate subdivision surfaces smooth position
1145 0 : if(markSH.get_subset_index(f) != -1 && linearManifoldSH.get_subset_index(f) == -1)
1146 : {
1147 : // TRIANGLE CASE
1148 0 : if(f->reference_object_id() == ROID_TRIANGLE)
1149 : {
1150 : // Iterate over all face vertices, calculate and apply local centroid masks
1151 0 : for(size_t i = 0; i < f->num_vertices(); ++i)
1152 : {
1153 : // Init
1154 0 : Vertex* vrt = f->vertex(i);
1155 : VecSet(p, 0);
1156 :
1157 : // Summate coordinates of neighbor vertices to vrt inside face
1158 0 : for(size_t j = 0; j < f->num_vertices(); ++j)
1159 : {
1160 0 : if(j != i)
1161 : {
1162 0 : VecAdd(p, p, aaPos[f->vertex(j)]);
1163 : }
1164 : }
1165 :
1166 : // Smooth vertex position
1167 : VecScaleAppend(aaSmoothBndPos_tri[vrt], 2.0/8, aaPos[vrt], 3.0/8, p);
1168 : }
1169 : }
1170 :
1171 : // QUADRILATERAL CASE
1172 0 : else if(f->reference_object_id() == ROID_QUADRILATERAL)
1173 : {
1174 : // Iterate over all face vertices, calculate and apply local centroid masks
1175 0 : for(size_t i = 0; i < f->num_vertices(); ++i)
1176 : {
1177 : // Init
1178 0 : Vertex* vrt = f->vertex(i);
1179 : VecSet(p, 0);
1180 :
1181 : // Summate coordinates of neighbor vertices to vrt inside face
1182 0 : for(size_t j = 0; j < f->num_vertices(); ++j)
1183 : {
1184 0 : if(j != i)
1185 : {
1186 0 : VecAdd(p, p, aaPos[f->vertex(j)]);
1187 : }
1188 : }
1189 :
1190 : // Smooth vertex position
1191 : VecScaleAppend(aaSmoothBndPos_quad[vrt], 1.0/4, aaPos[vrt], 1.0/4, p);
1192 : }
1193 : }
1194 :
1195 : // UNSUPPORTED MANIFOLD ELEMENT CASE
1196 : else
1197 0 : UG_THROW("ERROR in CalculateSmoothManifoldPosInTopLevelAveragingScheme: Non triangular/quadrilateral faces included in grid.");
1198 : }
1199 : }
1200 :
1201 : // Manage vertex attachment communication in parallel case -> COMMUNICATE aaSmoothBndPos
1202 : #ifdef UG_PARALLEL
1203 : // Reduce add operations:
1204 : // sum up h_slaves into h_masters
1205 :
1206 : // Copy operations:
1207 : // copy h_masters to h_slaves for consistency
1208 : AttachmentAllReduce<Vertex>(mg, aSmoothBndPos_tri, PCL_RO_SUM);
1209 : AttachmentAllReduce<Vertex>(mg, aSmoothBndPos_quad, PCL_RO_SUM);
1210 : #endif
1211 0 : }
1212 :
1213 :
1214 : ////////////////////////////////////////////////////////////////////////////////
1215 : template <class TAPosition>
1216 0 : void CalculateSmoothManifoldPosInParentLevelButterflyScheme(MultiGrid& mg, TAPosition& aPos,
1217 : MGSubsetHandler& markSH,
1218 : MGSubsetHandler& linearManifoldSH,
1219 : TAPosition& aSmoothBndPosOddVrt,
1220 : AInt& aNumManifoldEdges)
1221 : {
1222 : /*
1223 : * Scheme reference:
1224 : *
1225 : * D. N. Zorin, Interpolating Subdivision for Meshes with Arbitrary Topology,
1226 : * SIGGRAPH '96 Proceedings of the 23rd annual conference on Computer graphics
1227 : * and interactive techniques, 1996.
1228 : */
1229 :
1230 : if(TAPosition::ValueType::Size == 1){
1231 : UG_THROW("Error in CalculateSmoothManifoldPosInParentLevelButterflyScheme:\n"
1232 : "Currently only dimensions 2 and 3 are supported.\n");
1233 : }
1234 :
1235 : // WARNING: Parallel implementation has to be fixed
1236 : #ifdef UG_PARALLEL
1237 : UG_LOG("WARNING: CalculateSmoothManifoldPosInParentLevelButterflyScheme: "
1238 : "Parallel implementation has to be fixed." << std::endl);
1239 : #endif
1240 :
1241 : // Catch use of procedure for MultiGrids with just one level
1242 : size_t numLevels = mg.num_levels();
1243 :
1244 : #ifdef UG_PARALLEL
1245 : if(pcl::NumProcs() > 1){
1246 : pcl::ProcessCommunicator pc;
1247 : numLevels = pc.allreduce(numLevels, PCL_RO_MAX);
1248 : }
1249 : #endif
1250 :
1251 0 : if(numLevels == 1)
1252 : {
1253 0 : UG_THROW("Error in CalculateSmoothManifoldPosInParentLevelButterflyScheme: "
1254 : "Procedure only to be used for MultiGrids with more than one level.");
1255 : }
1256 :
1257 : // Define attachment accessors
1258 : Grid::VertexAttachmentAccessor<TAPosition> aaPos(mg, aPos);
1259 : Grid::EdgeAttachmentAccessor<TAPosition> aaSmoothBndPosOddVrt(mg, aSmoothBndPosOddVrt);
1260 : Grid::VertexAttachmentAccessor<AInt> aaNumManifoldEdges(mg, aNumManifoldEdges);
1261 :
1262 : #ifdef UG_PARALLEL
1263 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
1264 : #endif
1265 :
1266 : // Declare centroid coordinate vector
1267 : typedef typename TAPosition::ValueType position_type;
1268 : position_type p;
1269 : position_type q;
1270 : VecSet(p, 0);
1271 : VecSet(q, 0);
1272 :
1273 : /*
1274 : * Smoothing of odd vertices x
1275 : *
1276 : -1/16 1/8 -1/16
1277 : \ / \ /
1278 : 1/2--x--1/2
1279 : / \ / \
1280 : -1/16 1/8 -1/16
1281 : *
1282 : */
1283 :
1284 : // Calculate smooth position for ODD vertices (EVEN vertices will be interpolated by this scheme)
1285 0 : for(EdgeIterator eIter = mg.begin<Edge>(mg.top_level()-1); eIter != mg.end<Edge>(mg.top_level()-1); ++eIter)
1286 : {
1287 : // Reset centroids
1288 : VecSet(p, 0);
1289 : VecSet(q, 0);
1290 :
1291 : Edge* e = *eIter;
1292 :
1293 : // Skip ghost edges
1294 : #ifdef UG_PARALLEL
1295 : if(dgm.is_ghost(e))
1296 : continue;
1297 : #endif
1298 :
1299 : // In case of marked manifold edges, which do not belong to the user-specified linear boundary manifold subsets,
1300 : // and activated subdivision Butterfly refinement calculate subdivision surfaces smooth position
1301 0 : if(markSH.get_subset_index(e) != -1 && linearManifoldSH.get_subset_index(e) == -1)
1302 : {
1303 : // REGULAR CASE: both edge vertices are of valence 6
1304 0 : if(aaNumManifoldEdges[e->vertex(0)] == 6 && aaNumManifoldEdges[e->vertex(1)] == 6)
1305 : {
1306 : // perform Butterfly subdivision on odd manifold vertices
1307 : // get the neighbored manifold triangles
1308 : std::vector<Face*> associatedFaces;
1309 : std::vector<Face*> associatedButterflyFaces;
1310 : std::vector<Face*> associatedManifoldFaces;
1311 : std::vector<Face*> associatedButterflyManifoldFaces;
1312 :
1313 0 : CollectAssociated(associatedFaces, mg, e);
1314 :
1315 0 : for(size_t i = 0; i < associatedFaces.size(); ++i)
1316 : {
1317 : // Only consider associated faces, which are marked as manifold faces
1318 0 : if(markSH.get_subset_index(associatedFaces[i]) != -1)
1319 : {
1320 : // Exclude ghost and horizontal slave manifold faces
1321 : #ifdef UG_PARALLEL
1322 : if(dgm.is_ghost(associatedFaces[i]))
1323 : continue;
1324 :
1325 : if(dgm.contains_status(associatedFaces[i], ES_H_SLAVE))
1326 : continue;
1327 : #endif
1328 :
1329 0 : if(associatedManifoldFaces.size() < 2)
1330 : {
1331 0 : associatedManifoldFaces.push_back(associatedFaces[i]);
1332 : }
1333 : }
1334 : }
1335 :
1336 : // THROW, if more then 2 associated manifold faces have been found
1337 0 : if(associatedManifoldFaces.size() <= 2)
1338 : {
1339 : // Check, if all faces are triangles
1340 0 : for(size_t i = 0; i < associatedManifoldFaces.size(); ++i)
1341 : {
1342 0 : if(associatedManifoldFaces[i]->num_vertices() != 3)
1343 : {
1344 0 : UG_THROW("ERROR in CalculateSmoothManifoldPosInParentLevelButterflyScheme3d: "
1345 : "Non triangular faces included in grid: " << ElementDebugInfo(mg, associatedManifoldFaces[i]));
1346 : }
1347 : }
1348 :
1349 : // Summate centroid of face adjacent vertices (with corresponding weights 1/8)
1350 0 : for(size_t i = 0; i < associatedManifoldFaces.size(); ++i)
1351 : {
1352 0 : VecAdd(p, p, aaPos[GetConnectedVertex(e, associatedManifoldFaces[i])]);
1353 : }
1354 :
1355 : // Extend original "Loop's neighborhood diamond" to include 'BUTTERFLY VERTICES' (with corresponding weights -1/16)
1356 0 : for(size_t i = 0; i < associatedManifoldFaces.size(); ++i)
1357 : {
1358 : // Iterate over edges of original "Loop's neighborhood diamond" to extend to Butterfly neighborhood
1359 0 : for(size_t j = 0; j < associatedManifoldFaces[i]->num_edges(); ++j)
1360 : {
1361 : // Clear face container
1362 : associatedButterflyFaces.clear();
1363 : associatedButterflyManifoldFaces.clear();
1364 :
1365 : // Exclude edge e currently being edited
1366 0 : if(j != (size_t)GetEdgeIndex(associatedManifoldFaces[i], e))
1367 : {
1368 : // Collect associated Butterfly face adjacent to edge j
1369 0 : GetNeighbours(associatedButterflyFaces, mg, associatedManifoldFaces[i], j);
1370 :
1371 0 : for(size_t k = 0; k < associatedButterflyFaces.size(); ++k)
1372 : {
1373 : // Only consider associated butterfly faces, which are marked as manifold faces
1374 0 : if(markSH.get_subset_index(associatedButterflyFaces[k]) != -1)
1375 : {
1376 : // Exclude ghost and horizontal slave manifold faces
1377 : #ifdef UG_PARALLEL
1378 : if(dgm.is_ghost(associatedButterflyFaces[k]))
1379 : continue;
1380 :
1381 : if(dgm.contains_status(associatedButterflyFaces[k], ES_H_SLAVE))
1382 : continue;
1383 : #endif
1384 :
1385 0 : if(associatedButterflyFaces[k]->num_vertices() != 3)
1386 : {
1387 0 : UG_THROW("ERROR in CalculateSmoothManifoldPosInParentLevelButterflyScheme: "
1388 : "Non triangular faces included in grid: " << ElementDebugInfo(mg, associatedButterflyFaces[k]));
1389 : }
1390 :
1391 0 : associatedButterflyManifoldFaces.push_back(associatedButterflyFaces[k]);
1392 : }
1393 : }
1394 :
1395 0 : if(associatedButterflyManifoldFaces.size() != 1)
1396 0 : UG_THROW("ERROR in CalculateSmoothManifoldPosInParentLevelButterflyScheme: number of edge associated Butterfly Manifold faces != 1.");
1397 :
1398 : // Summate centroid of butterly face adjacent vertex
1399 0 : VecAdd(q, q, aaPos[GetConnectedVertex(mg.get_edge(associatedManifoldFaces[i]->edge_desc(j)), associatedButterflyManifoldFaces[0])]);
1400 : }
1401 : }
1402 : }
1403 :
1404 : // Exclude ghost and horizontal slaves of the parent edge vertices of the currently smoothed vertex
1405 : // to avoid multiple contributions to the centroid of the edge adjacent vertices
1406 : #ifdef UG_PARALLEL
1407 : if(dgm.is_ghost(e))
1408 : {
1409 : continue;
1410 : }
1411 :
1412 : if(dgm.contains_status(e, ES_H_SLAVE))
1413 : {
1414 : VecScaleAppend(aaSmoothBndPosOddVrt[e], 0.125, p, -1.0/16, q);
1415 : continue;
1416 : }
1417 : #endif
1418 :
1419 0 : VecScaleAppend(aaSmoothBndPosOddVrt[e], 0.5, aaPos[e->vertex(0)], 0.5, aaPos[e->vertex(1)], 0.125, p, -1.0/16, q);
1420 : }
1421 : else
1422 0 : UG_THROW("ERROR in CalculateSmoothManifoldPosInParentLevelButterflyScheme: numAssociatedManifoldFaces > 2.");
1423 0 : }
1424 :
1425 : // IRREGULAR CASE: at least one edge vertex irregular
1426 0 : if(aaNumManifoldEdges[e->vertex(0)] != 6 || aaNumManifoldEdges[e->vertex(1)] != 6)
1427 : {
1428 : // Number of centroids to calculate (1 or 2, depending on the valence of the vertices of e)
1429 : size_t numLoops;
1430 :
1431 0 : if((aaNumManifoldEdges[e->vertex(0)] != 6 && aaNumManifoldEdges[e->vertex(1)] == 6) ||
1432 0 : (aaNumManifoldEdges[e->vertex(0)] == 6 && aaNumManifoldEdges[e->vertex(1)] != 6))
1433 : {
1434 : numLoops = 1;
1435 : }
1436 : // case aaNumManifoldEdges[e->vertex(0)] != 6 && aaNumManifoldEdges[e->vertex(1)] != 6
1437 : else
1438 : numLoops = 2;
1439 :
1440 : // Loop centroids to calculate
1441 0 : for(size_t n = 0; n < numLoops; ++n)
1442 : {
1443 : // Reset centroids
1444 : VecSet(p, 0);
1445 : VecSet(q, 0);
1446 :
1447 : Vertex* vrt;
1448 : Vertex* butterflyVertex;
1449 :
1450 : std::vector<Face*> associatedFaces;
1451 : std::vector<Face*> associatedManifoldFaces;
1452 :
1453 : // Determine smoothing case
1454 0 : if(aaNumManifoldEdges[e->vertex(0)] != 6 && aaNumManifoldEdges[e->vertex(1)] == 6)
1455 : {
1456 0 : vrt = e->vertex(0);
1457 :
1458 : // Push back e->vertex(1) as vertex s_0
1459 0 : butterflyVertex = e->vertex(1);
1460 : }
1461 0 : else if (aaNumManifoldEdges[e->vertex(0)] == 6 && aaNumManifoldEdges[e->vertex(1)] != 6)
1462 : {
1463 0 : vrt = e->vertex(1);
1464 :
1465 : // Push back e->vertex(0) as vertex s_0
1466 0 : butterflyVertex = e->vertex(0);
1467 : }
1468 : // case aaNumManifoldEdges[e->vertex(0)] != 6 && aaNumManifoldEdges[e->vertex(1)] != 6
1469 : else
1470 : {
1471 0 : vrt = e->vertex(n);
1472 :
1473 : // Push back the other vertex as vertex s_0
1474 0 : butterflyVertex = e->vertex(1 - (n % 2));
1475 : }
1476 :
1477 : // perform Butterfly subdivision on odd manifold vertices
1478 : // get the neighbored manifold triangles
1479 0 : CollectAssociated(associatedFaces, mg, e);
1480 :
1481 0 : for(size_t i = 0; i < associatedFaces.size(); ++i)
1482 : {
1483 : // Only consider associated faces, which are marked as manifold faces
1484 0 : if(markSH.get_subset_index(associatedFaces[i]) != -1)
1485 : {
1486 : // Exclude ghost and horizontal slave manifold faces
1487 : #ifdef UG_PARALLEL
1488 : if(dgm.is_ghost(associatedFaces[i]))
1489 : continue;
1490 :
1491 : if(dgm.contains_status(associatedFaces[i], ES_H_SLAVE))
1492 : continue;
1493 : #endif
1494 :
1495 0 : if(associatedManifoldFaces.size() < 2)
1496 : {
1497 0 : associatedManifoldFaces.push_back(associatedFaces[i]);
1498 : }
1499 : }
1500 : }
1501 :
1502 0 : if(associatedManifoldFaces.size() <= 2)
1503 : {
1504 : // Check, if all faces are triangles
1505 0 : for(size_t i = 0; i < associatedManifoldFaces.size(); ++i)
1506 : {
1507 0 : if(associatedManifoldFaces[i]->num_vertices() != 3)
1508 : {
1509 0 : UG_THROW("ERROR in CalculateSmoothManifoldPosInParentLevelButterflyScheme3d: "
1510 : "Non triangular faces included in grid: " << ElementDebugInfo(mg, associatedManifoldFaces[i]));
1511 : }
1512 : }
1513 :
1514 : // Start with one triangle of "Loop's neighborhood diamond"
1515 0 : Face* f = associatedManifoldFaces[0];
1516 0 : size_t k = (size_t)aaNumManifoldEdges[vrt];
1517 :
1518 : number butterflyWeight = 0.0;
1519 : std::vector<number> butterflyWeights;
1520 :
1521 0 : if(k == 3)
1522 : {
1523 0 : butterflyWeights.push_back(5.0/12);
1524 0 : butterflyWeights.push_back(-1.0/12);
1525 0 : butterflyWeights.push_back(-1.0/12);
1526 : }
1527 :
1528 0 : if(k == 4)
1529 : {
1530 0 : butterflyWeights.push_back(3.0/8);
1531 0 : butterflyWeights.push_back(0.0);
1532 0 : butterflyWeights.push_back(-1.0/8);
1533 0 : butterflyWeights.push_back(0.0);
1534 : }
1535 :
1536 : // Special parallel treatment for s_0 in case e is ghost or horizontal slave
1537 0 : if(k != 3 && k != 4)
1538 0 : butterflyWeight = 1.0/k * 7.0/4;
1539 : else
1540 0 : butterflyWeight = butterflyWeights[0];
1541 :
1542 0 : VecScaleAppend(p, butterflyWeight, aaPos[butterflyVertex]);
1543 :
1544 : // Ordered traversing of associated edges of currently considered irregular vertex vrt
1545 0 : for(size_t i = 1; i < k; ++i)
1546 : {
1547 : // Clear face containers
1548 : associatedFaces.clear();
1549 : associatedManifoldFaces.clear();
1550 :
1551 : // Get connecting edge of next butterfly vertex s_i in line to vrt
1552 0 : if(f->get_opposing_object(butterflyVertex).first != EDGE)
1553 : {
1554 0 : UG_THROW("ERROR in CalculateSmoothManifoldPosInParentLevelButterflyScheme3d: "
1555 : "Opposing object of butterfly vertex in manifold face is not an edge: "
1556 : << ElementDebugInfo(mg, f))
1557 : }
1558 :
1559 0 : Edge* egdeOfNextFace = mg.get_edge(f, f->get_opposing_object(butterflyVertex).second);
1560 :
1561 : // Store butterfly vertex s_i
1562 0 : egdeOfNextFace->get_opposing_side(vrt, &butterflyVertex);
1563 :
1564 : // butterfly weight s_i
1565 0 : if(k != 3 && k != 4)
1566 0 : butterflyWeight = 1.0/k * (1.0/4 + cos(2*i*PI/k) + 1.0/2*cos(4*i*PI/k));
1567 : else
1568 0 : butterflyWeight = butterflyWeights[i];
1569 :
1570 : // Summate centroid of butterly vertices s_i
1571 0 : VecScaleAppend(q, butterflyWeight, aaPos[butterflyVertex]);
1572 :
1573 : // Get next face to traverse
1574 0 : GetNeighbours(associatedFaces, mg, f, GetEdgeIndex(f, egdeOfNextFace));
1575 :
1576 0 : for(size_t j = 0; j < associatedFaces.size(); ++j)
1577 : {
1578 : // Only consider associated faces, which are marked as manifold faces
1579 0 : if(markSH.get_subset_index(associatedFaces[j]) != -1)
1580 : {
1581 : // Exclude ghost and horizontal slave manifold faces
1582 : #ifdef UG_PARALLEL
1583 : if(dgm.is_ghost(associatedFaces[j]))
1584 : continue;
1585 :
1586 : if(dgm.contains_status(associatedFaces[j], ES_H_SLAVE))
1587 : continue;
1588 : #endif
1589 :
1590 0 : associatedManifoldFaces.push_back(associatedFaces[j]);
1591 : }
1592 : }
1593 :
1594 0 : if(associatedManifoldFaces.size() != 1)
1595 0 : UG_THROW("ERROR in CalculateSmoothManifoldPosInParentLevelButterflyScheme: number of edge associated Butterfly Manifold faces != 1.");
1596 :
1597 : // Store next face in line to traverse in next iteration
1598 0 : f = associatedManifoldFaces[0];
1599 : }
1600 :
1601 : // Exclude ghost and horizontal slaves of the parent edge vertices of the currently smoothed vertex
1602 : // to avoid multiple contributions to the centroid of the edge adjacent vertices
1603 : #ifdef UG_PARALLEL
1604 : if(dgm.is_ghost(e))
1605 : {
1606 : continue;
1607 : }
1608 :
1609 : if(dgm.contains_status(e, ES_H_SLAVE))
1610 : {
1611 : VecScaleAppend(aaSmoothBndPosOddVrt[e], 1.0/numLoops, q);
1612 : continue;
1613 : }
1614 : #endif
1615 :
1616 0 : VecScaleAppend(aaSmoothBndPosOddVrt[e], 1.0/numLoops*3.0/4, aaPos[vrt], 1.0/numLoops, p, 1.0/numLoops, q);
1617 0 : }
1618 : else
1619 0 : UG_THROW("ERROR in CalculateSmoothManifoldPosInParentLevelButterflyScheme: numAssociatedManifoldFaces > 2.");
1620 : }
1621 : }
1622 : }
1623 : }
1624 :
1625 : // Manage vertex and edge attachment communication in parallel case -> COMMUNICATE aSmoothBndPosEvenVrt, aSmoothBndPosOddVrt
1626 : #ifdef UG_PARALLEL
1627 : // Reduce add operations:
1628 : // sum up h_slaves into h_masters
1629 :
1630 : // Copy operations:
1631 : // copy h_masters to h_slaves for consistency
1632 : AttachmentAllReduce<Edge>(mg, aSmoothBndPosOddVrt, PCL_RO_SUM);
1633 : #endif
1634 0 : }
1635 :
1636 :
1637 : ////////////////////////////////////////////////////////////////////////////////
1638 0 : void CalculateSmoothVolumePosInTopLevel(MultiGrid& mg, MGSubsetHandler& markSH,
1639 : APosition& aSmoothVolPos_toc,
1640 : APosition& aSmoothVolPos_prism,
1641 : APosition& aSmoothVolPos_hex)
1642 : {
1643 : /*
1644 : * Scheme references:
1645 : *
1646 : * S. Schaefer, J. Hakenberg, and J. Warren, Smooth subdivision of tetrahedral meshes,
1647 : * Proceedings of the 2004 Eurographics/ACM Symposium on Geometry Processing.
1648 : *
1649 : * J. Hakenberg, Smooth Subdivision for Mixed Volumetric Meshes, thesis, 2004
1650 : */
1651 :
1652 : #ifdef UG_PARALLEL
1653 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
1654 : #endif
1655 :
1656 : // Define attachment accessors
1657 : Grid::VertexAttachmentAccessor<APosition> aaPos(mg, aPosition);
1658 : Grid::VertexAttachmentAccessor<APosition> aaSmoothVolPos_toc(mg, aSmoothVolPos_toc);
1659 : Grid::VertexAttachmentAccessor<APosition> aaSmoothVolPos_prism(mg, aSmoothVolPos_prism);
1660 : Grid::VertexAttachmentAccessor<APosition> aaSmoothVolPos_hex(mg, aSmoothVolPos_hex);
1661 :
1662 : // Declare volume centroid coordinate vector
1663 : typedef APosition::ValueType pos_type;
1664 : pos_type p;
1665 :
1666 : // Loop all volumes of top_level
1667 0 : for(VolumeIterator vIter = mg.begin<Volume>(mg.top_level()); vIter != mg.end<Volume>(mg.top_level()); ++vIter)
1668 : {
1669 : Volume* vol = *vIter;
1670 :
1671 : // Skip ghost volumes
1672 : #ifdef UG_PARALLEL
1673 : if(dgm.is_ghost(vol))
1674 : continue;
1675 : #endif
1676 :
1677 : // Iterate over all volume vertices, calculate and apply local centroid masks
1678 0 : for(size_t i = 0; i < vol->num_vertices(); ++i)
1679 : {
1680 : // Init
1681 0 : Vertex* vrt = vol->vertex(i);
1682 : VecSet(p, 0);
1683 :
1684 : // In case of linear or subdivision Loop boundary manifold refinement:
1685 : // handle vertices of separating manifolds separately
1686 0 : if(markSH.get_subset_index(vrt) != -1 && g_boundaryRefinementRule != SUBDIV_VOL)
1687 : {
1688 0 : continue;
1689 : }
1690 :
1691 : // TETRAHEDRON CASE
1692 0 : if(vol->reference_object_id() == ROID_TETRAHEDRON)
1693 : {
1694 : // Summate coordinates of neighbor vertices to vrt inside tetrahedron
1695 0 : for(size_t j = 0; j < vol->num_vertices(); ++j)
1696 : {
1697 0 : if(j != i)
1698 : {
1699 0 : VecAdd(p, p, aaPos[vol->vertex(j)]);
1700 : }
1701 : }
1702 :
1703 : // Smooth vertex position
1704 : VecScaleAppend(aaSmoothVolPos_toc[vrt], -1.0/16, aaPos[vrt], 17.0/48, p);
1705 : }
1706 :
1707 : // OCTAHEDRON CASE
1708 0 : else if(vol->reference_object_id() == ROID_OCTAHEDRON)
1709 : {
1710 : // Get cell-adjacent vertex
1711 0 : Vertex* oppVrt = vol->vertex(vol->get_opposing_object(vrt).second);
1712 :
1713 : // Summate coordinates of DIRECT neighbor vertices to vrt inside octahedron
1714 0 : for(size_t j = 0; j < vol->num_vertices(); ++j)
1715 : {
1716 0 : if(GetVertexIndex(vol, oppVrt) == -1)
1717 : {
1718 0 : UG_THROW("ERROR in CalculateSmoothVolumePosInTopLevel: identified opposing vertex actually not included in current volume.");
1719 : }
1720 :
1721 0 : if(j != i && j != (size_t)GetVertexIndex(vol, oppVrt))
1722 : {
1723 0 : VecAdd(p, p, aaPos[vol->vertex(j)]);
1724 : }
1725 : }
1726 :
1727 : // Smooth vertex position
1728 : VecScaleAppend(aaSmoothVolPos_toc[vrt], 3.0/8, aaPos[vrt], 1.0/12, p, 7.0/24, aaPos[oppVrt]);
1729 : }
1730 :
1731 : // PRISM CASE
1732 0 : else if(vol->reference_object_id() == ROID_PRISM)
1733 : {
1734 : // Get cell-adjacent vertex
1735 : Vertex* oppVrt;
1736 :
1737 : if(i == 0)
1738 0 : oppVrt = vol->vertex(3);
1739 : else if(i == 1)
1740 0 : oppVrt = vol->vertex(4);
1741 : else if(i == 2)
1742 0 : oppVrt = vol->vertex(5);
1743 : else if(i == 3)
1744 0 : oppVrt = vol->vertex(0);
1745 : else if(i == 4)
1746 0 : oppVrt = vol->vertex(1);
1747 : else
1748 0 : oppVrt = vol->vertex(2);
1749 :
1750 : // Summate coordinates of DIRECT neighbor vertices to vrt inside octahedron
1751 0 : for(size_t j = 0; j < vol->num_vertices(); ++j)
1752 : {
1753 0 : if(j != i && j != (size_t)GetVertexIndex(vol, oppVrt))
1754 : {
1755 0 : VecAdd(p, p, aaPos[vol->vertex(j)]);
1756 : }
1757 : }
1758 :
1759 : // Smooth vertex position
1760 : VecScaleAppend(aaSmoothVolPos_prism[vrt], 1.0/8, aaPos[vrt], 3.0/16, p, 1.0/8, aaPos[oppVrt]);
1761 : }
1762 :
1763 : // HEXAHEDRON CASE
1764 0 : else if(vol->reference_object_id() == ROID_HEXAHEDRON)
1765 : {
1766 : // Summate coordinates of neighbor vertices to vrt inside hexahedron
1767 0 : for(size_t j = 0; j < vol->num_vertices(); ++j)
1768 : {
1769 0 : if(j != i)
1770 : {
1771 0 : VecAdd(p, p, aaPos[vol->vertex(j)]);
1772 : }
1773 : }
1774 :
1775 : // Smooth vertex position
1776 : VecScaleAppend(aaSmoothVolPos_hex[vrt], 1.0/8, aaPos[vrt], 1.0/8, p);
1777 : }
1778 :
1779 : // UNSUPPORTED VOLUME ELEMENT CASE
1780 : else
1781 : {
1782 0 : UG_THROW("ERROR in CalculateSmoothVolumePosInTopLevel: Volume type not supported for subdivision volumes refinement.");
1783 : }
1784 : }
1785 : }
1786 :
1787 :
1788 : // Manage vertex attachment communication in parallel case -> COMMUNICATE aSmoothVolPos
1789 : #ifdef UG_PARALLEL
1790 : // Reduce add operations:
1791 : // sum up h_slaves into h_masters
1792 :
1793 : // Copy operations:
1794 : // copy h_masters to h_slaves for consistency
1795 : AttachmentAllReduce<Vertex>(mg, aSmoothVolPos_toc, PCL_RO_SUM);
1796 : AttachmentAllReduce<Vertex>(mg, aSmoothVolPos_prism, PCL_RO_SUM);
1797 : AttachmentAllReduce<Vertex>(mg, aSmoothVolPos_hex, PCL_RO_SUM);
1798 : #endif
1799 0 : }
1800 :
1801 :
1802 : ////////////////////////////////////////////////////////////////////////////////
1803 0 : void CalculateConstrainedSmoothVolumePosInTopLevel(MultiGrid& mg, MGSubsetHandler& markSH,
1804 : APosition& aSmoothVolPos_toc)
1805 : {
1806 : #ifdef UG_PARALLEL
1807 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
1808 : #endif
1809 :
1810 : // Define attachment accessors
1811 : Grid::VertexAttachmentAccessor<APosition> aaPos(mg, aPosition);
1812 : Grid::VertexAttachmentAccessor<APosition> aaSmoothVolPos_toc(mg, aSmoothVolPos_toc);
1813 :
1814 : // Declare volume centroid coordinate vector
1815 : typedef APosition::ValueType pos_type;
1816 : pos_type p;
1817 :
1818 : // boundary neighbor counter
1819 : size_t bndNbrCnt;
1820 :
1821 : // Loop all volumes of top_level
1822 0 : for(VolumeIterator vIter = mg.begin<Volume>(mg.top_level()); vIter != mg.end<Volume>(mg.top_level()); ++vIter)
1823 : {
1824 : Volume* vol = *vIter;
1825 :
1826 : // Skip ghost volumes
1827 : #ifdef UG_PARALLEL
1828 : if(dgm.is_ghost(vol))
1829 : continue;
1830 : #endif
1831 :
1832 : // Iterate over all volume vertices, calculate and apply local centroid masks
1833 0 : for(size_t i = 0; i < vol->num_vertices(); ++i)
1834 : {
1835 : // Init
1836 0 : Vertex* vrt = vol->vertex(i);
1837 : VecSet(p, 0);
1838 : bndNbrCnt = 0;
1839 :
1840 : // In case of linear or subdivision Loop boundary manifold refinement:
1841 : // handle vertices of separating manifolds separately
1842 0 : if(markSH.get_subset_index(vrt) != -1 && g_boundaryRefinementRule != SUBDIV_VOL)
1843 : {
1844 0 : continue;
1845 : }
1846 :
1847 : // TETRAHEDRON CASE
1848 0 : if(vol->reference_object_id() == ROID_TETRAHEDRON)
1849 : {
1850 : // Summate coordinates of neighbor vertices to vrt inside tetrahedron
1851 0 : for(size_t j = 0; j < vol->num_vertices(); ++j)
1852 : {
1853 0 : if(j != i)
1854 : {
1855 : // Only consider non-manifold neighbors
1856 0 : if(markSH.get_subset_index(vol->vertex(j)) == -1)
1857 : {
1858 0 : VecAdd(p, p, aaPos[vol->vertex(j)]);
1859 : }
1860 : else
1861 : {
1862 0 : bndNbrCnt++;
1863 : }
1864 : }
1865 : }
1866 :
1867 : // Smooth vertex position
1868 : //VecScaleAppend(aaSmoothVolPos_toc[vrt], -1.0/16 + bndNbrCnt*17.0/48, aaPos[vrt], 17.0/48, p);
1869 :
1870 0 : if(bndNbrCnt == 3)
1871 : VecScaleAppend(aaSmoothVolPos_toc[vrt], 1.0, aaPos[vrt]);
1872 0 : else if(bndNbrCnt == 0)
1873 : VecScaleAppend(aaSmoothVolPos_toc[vrt], -1.0/16, aaPos[vrt], 17.0/48, p);
1874 0 : else if(bndNbrCnt == 1)
1875 : VecScaleAppend(aaSmoothVolPos_toc[vrt], -1.0/16, aaPos[vrt], 17.0/48 + 17.0/(48*2), p);
1876 : else
1877 : VecScaleAppend(aaSmoothVolPos_toc[vrt], -1.0/16, aaPos[vrt], 51.0/48, p);
1878 : }
1879 :
1880 : // OCTAHEDRON CASE
1881 0 : else if(vol->reference_object_id() == ROID_OCTAHEDRON)
1882 : {
1883 : // Get cell-adjacent vertex
1884 0 : Vertex* oppVrt = vol->vertex(vol->get_opposing_object(vrt).second);
1885 :
1886 : // Summate coordinates of DIRECT neighbor vertices to vrt inside octahedron
1887 0 : for(size_t j = 0; j < vol->num_vertices(); ++j)
1888 : {
1889 0 : if(GetVertexIndex(vol, oppVrt) == -1)
1890 : {
1891 0 : UG_THROW("ERROR in CalculateConstrainedSmoothVolumePosInTopLevel: identified opposing vertex actually not included in current volume.");
1892 : }
1893 :
1894 0 : if(j != i && j != (size_t)GetVertexIndex(vol, oppVrt))
1895 : {
1896 : // Only consider non-manifold direct neighbors
1897 0 : if(markSH.get_subset_index(vol->vertex(j)) == -1)
1898 : {
1899 0 : VecAdd(p, p, aaPos[vol->vertex(j)]);
1900 : }
1901 : else
1902 : {
1903 0 : bndNbrCnt++;
1904 : }
1905 : }
1906 : }
1907 :
1908 : // Smooth vertex position
1909 :
1910 : // if opposing vertex is a non-manifold vertex
1911 0 : if(markSH.get_subset_index(oppVrt) == -1)
1912 : {
1913 : //VecScaleAppend(aaSmoothVolPos_toc[vrt], 3.0/8 + bndNbrCnt*1.0/12, aaPos[vrt], 1.0/12, p, 7.0/24, aaPos[oppVrt]);
1914 :
1915 0 : if(bndNbrCnt == 4)
1916 : VecScaleAppend(aaSmoothVolPos_toc[vrt], 3.0/8, aaPos[vrt], 15.0/24, aaPos[oppVrt]);
1917 0 : else if(bndNbrCnt == 0)
1918 : VecScaleAppend(aaSmoothVolPos_toc[vrt], 3.0/8, aaPos[vrt], 1.0/12, p, 7.0/24, aaPos[oppVrt]);
1919 0 : else if(bndNbrCnt == 1)
1920 : VecScaleAppend(aaSmoothVolPos_toc[vrt], 3.0/8, aaPos[vrt], 1.0/12 + 1.0/(12*3), p, 7.0/24, aaPos[oppVrt]);
1921 0 : else if(bndNbrCnt == 2)
1922 : VecScaleAppend(aaSmoothVolPos_toc[vrt], 3.0/8, aaPos[vrt], 2.0/12, p, 7.0/24, aaPos[oppVrt]);
1923 : else
1924 : VecScaleAppend(aaSmoothVolPos_toc[vrt], 3.0/8, aaPos[vrt], 4.0/12, p, 7.0/24, aaPos[oppVrt]);
1925 : }
1926 : else
1927 : {
1928 : //VecScaleAppend(aaSmoothVolPos_toc[vrt], 3.0/8 + bndNbrCnt*1.0/12 + 7.0/24, aaPos[vrt], 1.0/12, p);
1929 :
1930 0 : if(bndNbrCnt == 4)
1931 : VecScaleAppend(aaSmoothVolPos_toc[vrt], 1.0, aaPos[vrt]);
1932 0 : else if(bndNbrCnt == 0)
1933 : VecScaleAppend(aaSmoothVolPos_toc[vrt], 3.0/8, aaPos[vrt], 1.0/12 + 7.0/(24*4), p);
1934 0 : else if(bndNbrCnt == 1)
1935 : VecScaleAppend(aaSmoothVolPos_toc[vrt], 3.0/8, aaPos[vrt], 1.0/12 + 1.0/(12*3) + 7.0/(24*3), p);
1936 0 : else if(bndNbrCnt == 2)
1937 : VecScaleAppend(aaSmoothVolPos_toc[vrt], 3.0/8, aaPos[vrt], 2.0/12 + 7.0/(24*2), p);
1938 : else
1939 : VecScaleAppend(aaSmoothVolPos_toc[vrt], 3.0/8, aaPos[vrt], 4.0/12 + 7.0/24, p);
1940 : }
1941 : }
1942 :
1943 : // UNSUPPORTED VOLUME ELEMENT CASE
1944 : else
1945 : {
1946 0 : UG_THROW("ERROR in CalculateConstrainedSmoothVolumePosInTopLevel: Volume type not supported for subdivision volumes refinement, only tetrahedra and octahedra.");
1947 : }
1948 : }
1949 : }
1950 :
1951 : // Manage vertex attachment communication in parallel case -> COMMUNICATE aSmoothVolPos
1952 : #ifdef UG_PARALLEL
1953 : // Reduce add operations:
1954 : // sum up h_slaves into h_masters
1955 :
1956 : // Copy operations:
1957 : // copy h_masters to h_slaves for consistency
1958 : AttachmentAllReduce<Vertex>(mg, aSmoothVolPos_toc, PCL_RO_SUM);
1959 : #endif
1960 0 : }
1961 :
1962 :
1963 : ////////////////////////////////////////////////////////////////////////////////
1964 0 : void CalculateNumElemsVertexAttachmentInTopLevel(MultiGrid& mg, AInt& aNumElems_toc, AInt& aNumElems_prism, AInt& aNumElems_hex)
1965 : {
1966 : // Define attachment accessor
1967 : Grid::VertexAttachmentAccessor<AInt> aaNumElems_toc(mg, aNumElems_toc);
1968 : Grid::VertexAttachmentAccessor<AInt> aaNumElems_prism(mg, aNumElems_prism);
1969 : Grid::VertexAttachmentAccessor<AInt> aaNumElems_hex(mg, aNumElems_hex);
1970 :
1971 : // Manage vertex attachment communication in parallel case:
1972 : // - Setup communication policy for the above attachment
1973 : // - Setup interface communicator
1974 : // - Setup distributed grid manager
1975 : // - Setup grid layout map
1976 : #ifdef UG_PARALLEL
1977 : // Attachment communication policies COPY
1978 : ComPol_CopyAttachment<VertexLayout, AInt> comPolCopyNumElems_toc(mg, aNumElems_toc);
1979 : ComPol_CopyAttachment<VertexLayout, AInt> comPolCopyNumElems_prism(mg, aNumElems_prism);
1980 : ComPol_CopyAttachment<VertexLayout, AInt> comPolCopyNumElems_hex(mg, aNumElems_hex);
1981 :
1982 : // Interface communicators and distributed domain manager
1983 : pcl::InterfaceCommunicator<VertexLayout> com;
1984 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
1985 : GridLayoutMap& glm = dgm.grid_layout_map();
1986 : #endif
1987 :
1988 : // Loop all volumes of top level and calculate number of volumes each vertex is contained by
1989 0 : for(VolumeIterator vIter = mg.begin<Volume>(mg.top_level()); vIter != mg.end<Volume>(mg.top_level()); ++vIter)
1990 : {
1991 : Volume* vol = *vIter;
1992 :
1993 : // Skip ghosts
1994 : #ifdef UG_PARALLEL
1995 : if(dgm.is_ghost(vol))
1996 : continue;
1997 : #endif
1998 :
1999 0 : for(size_t i = 0; i < vol->num_vertices(); ++i)
2000 : {
2001 0 : if(vol->reference_object_id() == ROID_TETRAHEDRON || vol->reference_object_id() == ROID_OCTAHEDRON)
2002 0 : ++aaNumElems_toc[vol->vertex(i)];
2003 0 : else if(vol->reference_object_id() == ROID_PRISM)
2004 0 : ++aaNumElems_prism[vol->vertex(i)];
2005 0 : else if(vol->reference_object_id() == ROID_HEXAHEDRON)
2006 0 : ++aaNumElems_hex[vol->vertex(i)];
2007 : else
2008 0 : UG_THROW("ERROR in CalculateNumElemsVertexAttachmentInTopLevel: not supported element type included in grid.");
2009 : }
2010 : }
2011 :
2012 : // Manage vertex attachment communication in parallel case -> COMMUNICATE aNumElems
2013 : #ifdef UG_PARALLEL
2014 : // Reduce add operations:
2015 : // sum up h_slaves into h_masters
2016 :
2017 : // Copy operations:
2018 : // copy h_masters to h_slaves for consistency
2019 : AttachmentAllReduce<Vertex>(mg, aNumElems_toc, PCL_RO_SUM);
2020 : AttachmentAllReduce<Vertex>(mg, aNumElems_prism, PCL_RO_SUM);
2021 : AttachmentAllReduce<Vertex>(mg, aNumElems_hex, PCL_RO_SUM);
2022 :
2023 : // copy v_slaves to ghosts = VMASTER
2024 : com.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, comPolCopyNumElems_toc);
2025 : com.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, comPolCopyNumElems_prism);
2026 : com.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, comPolCopyNumElems_hex);
2027 : com.communicate();
2028 : #endif
2029 0 : }
2030 :
2031 :
2032 : ////////////////////////////////////////////////////////////////////////////////
2033 0 : void CalculateNumManifoldEdgesVertexAttachmentInParentLevel(MultiGrid& mg, MGSubsetHandler& markSH,
2034 : AInt& aNumManifoldEdges, bool bCreaseSurf)
2035 : {
2036 : // Define attachment accessor
2037 : Grid::VertexAttachmentAccessor<AInt> aaNumManifoldEdges(mg, aNumManifoldEdges);
2038 :
2039 : // Manage vertex attachment communication in parallel case:
2040 : // - Setup communication policy for the above attachment
2041 : // - Setup interface communicator
2042 : // - Setup distributed grid manager
2043 : // - Setup grid layout map
2044 : #ifdef UG_PARALLEL
2045 : // Attachment communication policies COPY
2046 : ComPol_CopyAttachment<VertexLayout, AInt> comPolCopyNumManifoldEdges(mg, aNumManifoldEdges);
2047 :
2048 : // Interface communicators and distributed domain manager
2049 : pcl::InterfaceCommunicator<VertexLayout> com;
2050 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
2051 : GridLayoutMap& glm = dgm.grid_layout_map();
2052 : #endif
2053 :
2054 : // Catch use of procedure for MultiGrids with just one level
2055 : size_t numLevels = mg.num_levels();
2056 :
2057 : #ifdef UG_PARALLEL
2058 : if(pcl::NumProcs() > 1){
2059 : pcl::ProcessCommunicator pc;
2060 : numLevels = pc.allreduce(numLevels, PCL_RO_MAX);
2061 : }
2062 : #endif
2063 :
2064 0 : if(numLevels == 1)
2065 0 : UG_THROW("CalculateNumManifoldEdgesVertexAttachmentInParentLevel: method may not be used in base level 0.");
2066 :
2067 : // Loop all edges of parent level and calculate number of associated manifold edges of each vertex
2068 0 : for(EdgeIterator eIter = mg.begin<Edge>(mg.top_level()-1); eIter != mg.end<Edge>(mg.top_level()-1); ++eIter)
2069 : {
2070 : Edge* e = *eIter;
2071 :
2072 : // Smooth surfaces:
2073 : // Count only marked manifold edges
2074 : // -> bool bSwitch = (markSH.get_subset_index(e) != -1);
2075 :
2076 : // Crease surfaces;
2077 : // Count all edges
2078 : // -> bSwitch = true;
2079 :
2080 : bool bSwitch;
2081 :
2082 0 : if(bCreaseSurf)
2083 : bSwitch = true;
2084 : else
2085 0 : bSwitch = (markSH.get_subset_index(e) != -1);
2086 :
2087 : // Check, if edge is contained in subset with marked manifold elements
2088 : // if (markSH.get_subset_index(e) != -1)
2089 0 : if(bSwitch)
2090 : {
2091 : // Skip ghost and horizontal slave edges
2092 : #ifdef UG_PARALLEL
2093 : if(dgm.is_ghost(e))
2094 : continue;
2095 :
2096 : if(dgm.contains_status(e, ES_H_SLAVE))
2097 : continue;
2098 : #endif
2099 :
2100 0 : ++aaNumManifoldEdges[e->vertex(0)];
2101 0 : ++aaNumManifoldEdges[e->vertex(1)];
2102 : }
2103 : }
2104 :
2105 : // Manage vertex attachment communication in parallel case -> COMMUNICATE aNumManifoldEdges
2106 : #ifdef UG_PARALLEL
2107 : // Reduce add operations:
2108 : // sum up h_slaves into h_masters
2109 :
2110 : // Copy operations:
2111 : // copy h_masters to h_slaves for consistency
2112 : AttachmentAllReduce<Vertex>(mg, aNumManifoldEdges, PCL_RO_SUM);
2113 :
2114 : // copy v_slaves to ghosts = VMASTER
2115 : com.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, comPolCopyNumManifoldEdges);
2116 : com.communicate();
2117 : #endif
2118 0 : }
2119 :
2120 :
2121 : ////////////////////////////////////////////////////////////////////////////////
2122 0 : void CalculateNumManifoldFacesVertexAttachmentInTopLevel(MultiGrid& mg, MGSubsetHandler& markSH, AInt& aNumManifoldFaces_tri, AInt& aNumManifoldFaces_quad)
2123 : {
2124 : // Define attachment accessor
2125 : Grid::VertexAttachmentAccessor<AInt> aaNumManifoldFaces_tri(mg, aNumManifoldFaces_tri);
2126 : Grid::VertexAttachmentAccessor<AInt> aaNumManifoldFaces_quad(mg, aNumManifoldFaces_quad);
2127 :
2128 : // Manage vertex attachment communication in parallel case:
2129 : // - Setup communication policy for the above attachment
2130 : // - Setup interface communicator
2131 : // - Setup distributed grid manager
2132 : // - Setup grid layout map
2133 : #ifdef UG_PARALLEL
2134 : // Attachment communication policies COPY
2135 : ComPol_CopyAttachment<VertexLayout, AInt> comPolCopyNumManifoldFaces_tri(mg, aNumManifoldFaces_tri);
2136 : ComPol_CopyAttachment<VertexLayout, AInt> comPolCopyNumManifoldFaces_quad(mg, aNumManifoldFaces_quad);
2137 :
2138 : // Interface communicators and distributed domain manager
2139 : pcl::InterfaceCommunicator<VertexLayout> com;
2140 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
2141 : GridLayoutMap& glm = dgm.grid_layout_map();
2142 : #endif
2143 :
2144 : // Loop all manifold faces of top level and calculate number of faces each vertex is contained by
2145 0 : for(FaceIterator fIter = mg.begin<Face>(mg.top_level()); fIter != mg.end<Face>(mg.top_level()); ++fIter)
2146 : {
2147 : Face* f = *fIter;
2148 :
2149 : // Only consider boundary manifold faces
2150 0 : if(markSH.get_subset_index(f) != -1)
2151 : {
2152 : // Skip ghosts
2153 : #ifdef UG_PARALLEL
2154 : if(dgm.is_ghost(f))
2155 : continue;
2156 : #endif
2157 :
2158 0 : for(size_t i = 0; i < f->num_vertices(); ++i)
2159 : {
2160 0 : if(f->reference_object_id() == ROID_TRIANGLE)
2161 0 : ++aaNumManifoldFaces_tri[f->vertex(i)];
2162 0 : else if(f->reference_object_id() == ROID_QUADRILATERAL)
2163 0 : ++aaNumManifoldFaces_quad[f->vertex(i)];
2164 : else
2165 0 : UG_THROW("ERROR in CalculateNumManifoldFacesVertexAttachmentInTopLevel: Non triangular/quadrilateral faces included in grid.");
2166 : }
2167 : }
2168 : }
2169 :
2170 : // Manage vertex attachment communication in parallel case -> COMMUNICATE aNumElems
2171 : #ifdef UG_PARALLEL
2172 : // Reduce add operations:
2173 : // sum up h_slaves into h_masters
2174 :
2175 : // Copy operations:
2176 : // copy h_masters to h_slaves for consistency
2177 : AttachmentAllReduce<Vertex>(mg, aNumManifoldFaces_tri, PCL_RO_SUM);
2178 : AttachmentAllReduce<Vertex>(mg, aNumManifoldFaces_quad, PCL_RO_SUM);
2179 :
2180 : // copy v_slaves to ghosts = VMASTER
2181 : com.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, comPolCopyNumManifoldFaces_tri);
2182 : com.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, comPolCopyNumManifoldFaces_quad);
2183 : com.communicate();
2184 : #endif
2185 0 : }
2186 :
2187 :
2188 : ////////////////////////////////////////////////////////////////////////////////
2189 0 : void InitLinearManifoldSubsetHandler(MultiGrid& mg, MGSubsetHandler& sh,
2190 : MGSubsetHandler& linearManifoldSH,
2191 : const char* linearManifoldSubsets)
2192 : {
2193 0 : if(linearManifoldSubsets[0] == '\0')
2194 0 : return;
2195 :
2196 : // tokenize user input
2197 0 : std::vector<std::string> linearManifoldSubsetsString = TokenizeString(linearManifoldSubsets);
2198 :
2199 : // remove white space
2200 0 : for(size_t i = 0; i < linearManifoldSubsetsString.size(); ++i)
2201 : {
2202 0 : RemoveWhitespaceFromString(linearManifoldSubsetsString[i]);
2203 : }
2204 :
2205 : // if no subset passed, clear subsets
2206 0 : if(linearManifoldSubsetsString.size() == 1 && linearManifoldSubsetsString[0].empty())
2207 : {
2208 : linearManifoldSubsetsString.clear();
2209 : }
2210 :
2211 : // if subsets passed with separator, but not all tokens filled, throw error
2212 0 : for(size_t i = 0; i < linearManifoldSubsetsString.size(); ++i)
2213 : {
2214 0 : if(linearManifoldSubsetsString.empty())
2215 : {
2216 0 : UG_THROW( "ERROR in InitLinearManifoldSubsetHandler: "
2217 : "linear boundary manifold subsets string passed lacks a "
2218 : "subset specification at position "<<i<<"(of "
2219 : <<linearManifoldSubsetsString.size()-1<<")");
2220 : }
2221 : }
2222 :
2223 : // assign all user specified vertices to linear boundary manifold SubsetHandler
2224 0 : for(size_t i = 0; i < linearManifoldSubsetsString.size(); ++i)
2225 : {
2226 0 : int j = sh.get_subset_index(linearManifoldSubsetsString[i].c_str());
2227 0 : UG_COND_THROW(j < 0, "Linear manifold subset named '"
2228 : << linearManifoldSubsetsString[i] << "' could not be identified.");
2229 :
2230 0 : for(VertexIterator vrtIter = sh.begin<Vertex>(j, mg.top_level()); vrtIter != sh.end<Vertex>(j, mg.top_level()); ++vrtIter)
2231 : {
2232 : Vertex* vrt = *vrtIter;
2233 0 : linearManifoldSH.assign_subset(vrt, 0);
2234 : }
2235 :
2236 : /* In case InitLinearManifoldSubsetHandler is not used as part of
2237 : * the GlobalSubdivisionDomainRefiner factory solely in the base level */
2238 0 : if(mg.num_levels() > 1){
2239 0 : for(VertexIterator vrtIter = sh.begin<Vertex>(j, mg.top_level()-1); vrtIter != sh.end<Vertex>(j, mg.top_level()-1); ++vrtIter)
2240 : {
2241 : Vertex* vrt = *vrtIter;
2242 0 : linearManifoldSH.assign_subset(vrt, 0);
2243 : }
2244 : }
2245 : }
2246 :
2247 : // assign all user specified edges to linear boundary manifold SubsetHandler
2248 0 : for(size_t i = 0; i < linearManifoldSubsetsString.size(); ++i)
2249 : {
2250 0 : int j = sh.get_subset_index(linearManifoldSubsetsString[i].c_str());
2251 :
2252 0 : for(EdgeIterator eIter = sh.begin<Edge>(j, mg.top_level()); eIter != sh.end<Edge>(j, mg.top_level()); ++eIter)
2253 : {
2254 : Edge* e = *eIter;
2255 0 : linearManifoldSH.assign_subset(e, 0);
2256 : }
2257 :
2258 : /* In case InitLinearManifoldSubsetHandler is not used as part of
2259 : * the GlobalSubdivisionDomainRefiner factory solely in the base level */
2260 0 : if(mg.num_levels() > 1){
2261 0 : for(EdgeIterator eIter = sh.begin<Edge>(j, mg.top_level()-1); eIter != sh.end<Edge>(j, mg.top_level()-1); ++eIter)
2262 : {
2263 : Edge* e = *eIter;
2264 0 : linearManifoldSH.assign_subset(e, 0);
2265 : }
2266 : }
2267 : }
2268 :
2269 : // assign all user specified faces to linear boundary manifold SubsetHandler
2270 0 : for(size_t i = 0; i < linearManifoldSubsetsString.size(); ++i)
2271 : {
2272 0 : int j = sh.get_subset_index(linearManifoldSubsetsString[i].c_str());
2273 :
2274 0 : for(FaceIterator fIter = sh.begin<Face>(j, mg.top_level()); fIter != sh.end<Face>(j, mg.top_level()); ++fIter)
2275 : {
2276 : Face* f = *fIter;
2277 0 : linearManifoldSH.assign_subset(f, 0);
2278 : }
2279 :
2280 : /* In case InitLinearManifoldSubsetHandler is not used as part of
2281 : * the GlobalSubdivisionDomainRefiner factory solely in the base level */
2282 0 : if(mg.num_levels() > 1){
2283 0 : for(FaceIterator fIter = sh.begin<Face>(j, mg.top_level()-1); fIter != sh.end<Face>(j, mg.top_level()-1); ++fIter)
2284 : {
2285 : Face* f = *fIter;
2286 0 : linearManifoldSH.assign_subset(f, 0);
2287 : }
2288 : }
2289 : }
2290 :
2291 : // Debug log
2292 : // UG_LOG("InitLinearManifoldSubsetHandler:" << std::endl);
2293 : // UG_LOG(">> Applying linear subdivision on the following boundary manifold subsets:" << std::endl);
2294 : //
2295 : // for(size_t i = 0; i < linearManifoldSubsetsString.size(); ++i)
2296 : // {
2297 : // UG_LOG("Subset # " << sh.get_subset_index(linearManifoldSubsetsString[i].c_str()) << ": " << linearManifoldSubsetsString[i] << std::endl);
2298 : // }
2299 0 : }
2300 :
2301 :
2302 : ////////////////////////////////////////////////////////////////////////////////
2303 : template <class TAPosition>
2304 0 : void ApplySmoothManifoldPosToTopLevelLoopScheme(MultiGrid& mg, TAPosition& aPos, MGSubsetHandler& markSH,
2305 : MGSubsetHandler& linearManifoldSH, bool bCreaseSurf)
2306 : {
2307 : if(TAPosition::ValueType::Size == 1){
2308 : UG_THROW("Error in ApplySmoothManifoldPosToTopLevelLoopScheme:\n"
2309 : "Currently only dimensions 2 and 3 are supported.\n");
2310 : }
2311 :
2312 : // Catch use of procedure for MultiGrids with just one level
2313 : size_t numLevels = mg.num_levels();
2314 :
2315 : #ifdef UG_PARALLEL
2316 : if(pcl::NumProcs() > 1){
2317 : pcl::ProcessCommunicator pc;
2318 : numLevels = pc.allreduce(numLevels, PCL_RO_MAX);
2319 : }
2320 : #endif
2321 :
2322 0 : if(numLevels == 1)
2323 : {
2324 0 : UG_THROW("Error in ApplySmoothManifoldPosToTopLevelLoopScheme: "
2325 : "Procedure only to be used for MultiGrids with more than one level.");
2326 : }
2327 :
2328 :
2329 : /*****************************************
2330 : *
2331 : * (1) SETUP
2332 : *
2333 : *****************************************/
2334 :
2335 : // Position attachment value type
2336 : typedef typename TAPosition::ValueType position_type;
2337 :
2338 : // Vertex attachments for associated number of manifold edges and smooth position
2339 : // (distinguish between volume and boundary smooth vertex positions
2340 : // and in case of boundary between EVEN and ODD smooth vertex positions)
2341 : AInt aNumManifoldEdges;
2342 : TAPosition aSmoothBndPosEvenVrt;
2343 : TAPosition aSmoothBndPosOddVrt;
2344 :
2345 : // attach previously declared vertex attachments with initial value 0
2346 0 : mg.attach_to_vertices_dv(aNumManifoldEdges, 0);
2347 0 : mg.attach_to_vertices_dv(aSmoothBndPosEvenVrt, position_type(0));
2348 0 : mg.attach_to_edges_dv(aSmoothBndPosOddVrt, position_type(0));
2349 :
2350 : // Define attachment accessors
2351 : Grid::VertexAttachmentAccessor<TAPosition> aaPos(mg, aPos);
2352 : Grid::VertexAttachmentAccessor<TAPosition> aaSmoothBndPosEvenVrt(mg, aSmoothBndPosEvenVrt);
2353 : Grid::EdgeAttachmentAccessor<TAPosition> aaSmoothBndPosOddVrt(mg, aSmoothBndPosOddVrt);
2354 :
2355 : // Manage vertex attachment communication in parallel case:
2356 : // - Setup communication policy for the above attachment aPosition
2357 : // - Setup interface communicator
2358 : // - Setup distributed grid manager
2359 : // - Setup grid layout map
2360 : #ifdef UG_PARALLEL
2361 : // Attachment communication policies COPY
2362 : ComPol_CopyAttachment<VertexLayout, TAPosition> comPolCopyAPosition(mg, aPos);
2363 :
2364 : // Interface communicators and distributed domain manager
2365 : pcl::InterfaceCommunicator<VertexLayout> com;
2366 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
2367 : GridLayoutMap& glm = dgm.grid_layout_map();
2368 : #endif
2369 :
2370 :
2371 : /*****************************************
2372 : *
2373 : * (2) DETERMINE aNumManifoldEdges
2374 : *
2375 : *****************************************/
2376 :
2377 0 : CalculateNumManifoldEdgesVertexAttachmentInParentLevel(mg, markSH, aNumManifoldEdges, bCreaseSurf);
2378 :
2379 :
2380 : /*****************************************
2381 : *
2382 : * (3) CALCULATE aSmoothBndPosEvenVrt,
2383 : * aSmoothBndPosOddVrt
2384 : *
2385 : *****************************************/
2386 :
2387 : // Calculate aSmoothBndPosEvenVrt, aSmoothBndPosOddVrt
2388 0 : if(bCreaseSurf)
2389 0 : CalculateSmoothCreaseManifoldPosInParentLevelLoopScheme(mg, aPos, markSH, linearManifoldSH,
2390 : aSmoothBndPosEvenVrt, aSmoothBndPosOddVrt, aNumManifoldEdges);
2391 : else
2392 0 : CalculateSmoothManifoldPosInParentLevelLoopScheme(mg, aPos, markSH, linearManifoldSH,
2393 : aSmoothBndPosEvenVrt, aSmoothBndPosOddVrt, aNumManifoldEdges);
2394 :
2395 :
2396 : /*****************************************
2397 : *
2398 : * (4) APPLY
2399 : *
2400 : *****************************************/
2401 :
2402 : // Loop all vertices of top_level
2403 0 : for(VertexIterator vrtIter = mg.begin<Vertex>(mg.top_level()); vrtIter != mg.end<Vertex>(mg.top_level()); ++vrtIter)
2404 : {
2405 : Vertex* vrt = *vrtIter;
2406 :
2407 : // Catch vertices without parent
2408 0 : if(mg.get_parent(vrt) == NULL)
2409 0 : continue;
2410 :
2411 : // Smooth surfaces:
2412 : // In case of marked manifold vertices, which do not belong to the user-specified linear boundary manifold subsets,
2413 : // and activated Loop scheme refinement apply subdivision surfaces smoothing, else linear refinement
2414 : // -> bool bSwitch = (markSH.get_subset_index(vrt) != -1 && linearManifoldSH.get_subset_index(vrt) == -1);
2415 :
2416 : // Crease surfaces;
2417 : // In case of vertices, which do not belong to the user-specified linear boundary manifold subsets,
2418 : // and activated Loop scheme refinement apply subdivision surfaces smoothing, else linear refinement
2419 : // -> bSwitch = (linearManifoldSH.get_subset_index(vrt) == -1);
2420 :
2421 : bool bSwitch;
2422 :
2423 0 : if(bCreaseSurf)
2424 0 : bSwitch = (linearManifoldSH.get_subset_index(vrt) == -1);
2425 : else
2426 0 : bSwitch = (markSH.get_subset_index(vrt) != -1 && linearManifoldSH.get_subset_index(vrt) == -1);
2427 :
2428 0 : if(bSwitch)
2429 : {
2430 : // EVEN VERTEX
2431 0 : if(mg.get_parent(vrt)->reference_object_id() == ROID_VERTEX)
2432 : {
2433 : // Get parent vertex
2434 : Vertex* parentVrt = static_cast<Vertex*>(mg.get_parent(vrt));
2435 :
2436 : aaPos[vrt] = aaSmoothBndPosEvenVrt[parentVrt];
2437 : }
2438 :
2439 : // ODD VERTEX
2440 0 : else if(mg.get_parent(vrt)->reference_object_id() == ROID_EDGE)
2441 : {
2442 : // Get parent edge
2443 : Edge* parentEdge = static_cast<Edge*>(mg.get_parent(vrt));
2444 :
2445 : aaPos[vrt] = aaSmoothBndPosOddVrt[parentEdge];
2446 : }
2447 : }
2448 : }
2449 :
2450 :
2451 : /*****************************************
2452 : *
2453 : * (5) COMMUNICATE VERTICALLY
2454 : * AFTER SUBDIVISION SURFACES
2455 : *
2456 : *****************************************/
2457 :
2458 : // Communicate aPosition in parallel case
2459 : #ifdef UG_PARALLEL
2460 : // copy ghosts = VMASTER to v_slaves
2461 : com.exchange_data(glm, INT_V_MASTER, INT_V_SLAVE, comPolCopyAPosition);
2462 : com.communicate();
2463 : #endif
2464 :
2465 :
2466 : /*****************************************
2467 : *
2468 : * (6) CLEAN UP
2469 : *
2470 : *****************************************/
2471 :
2472 : // detach vertex attachments
2473 : mg.detach_from_vertices(aNumManifoldEdges);
2474 : mg.detach_from_vertices(aSmoothBndPosEvenVrt);
2475 : mg.detach_from_edges(aSmoothBndPosOddVrt);
2476 0 : }
2477 :
2478 :
2479 : ////////////////////////////////////////////////////////////////////////////////
2480 : template <class TAPosition>
2481 0 : void ApplySmoothManifoldPosToTopLevelButterflyScheme(MultiGrid& mg, TAPosition& aPos, MGSubsetHandler& markSH,
2482 : MGSubsetHandler& linearManifoldSH)
2483 : {
2484 : if(TAPosition::ValueType::Size == 1){
2485 : UG_THROW("Error in ApplySmoothManifoldPosToTopLevelButterflyScheme:\n"
2486 : "Currently only dimensions 2 and 3 are supported.\n");
2487 : }
2488 :
2489 : // Catch use of procedure for MultiGrids with just one level
2490 : size_t numLevels = mg.num_levels();
2491 :
2492 : #ifdef UG_PARALLEL
2493 : if(pcl::NumProcs() > 1){
2494 : pcl::ProcessCommunicator pc;
2495 : numLevels = pc.allreduce(numLevels, PCL_RO_MAX);
2496 : }
2497 : #endif
2498 :
2499 0 : if(numLevels == 1)
2500 : {
2501 0 : UG_THROW("Error in ApplySmoothManifoldPosToTopLevelButterflyScheme: "
2502 : "Procedure only to be used for MultiGrids with more than one level.");
2503 : }
2504 :
2505 :
2506 : /*****************************************
2507 : *
2508 : * (1) SETUP
2509 : *
2510 : *****************************************/
2511 :
2512 : // Position attachment value type
2513 : typedef typename TAPosition::ValueType position_type;
2514 :
2515 : // Vertex attachments for associated number of manifold edges and smooth position
2516 : // (distinguish between volume and boundary smooth vertex positions
2517 : // and in case of boundary between EVEN and ODD smooth vertex positions)
2518 : AInt aNumManifoldEdges;
2519 : TAPosition aSmoothBndPosOddVrt;
2520 :
2521 : // attach previously declared vertex attachments with initial value 0
2522 0 : mg.attach_to_vertices_dv(aNumManifoldEdges, 0);
2523 0 : mg.attach_to_edges_dv(aSmoothBndPosOddVrt, position_type(0));
2524 :
2525 : // Define attachment accessors
2526 : Grid::VertexAttachmentAccessor<TAPosition> aaPos(mg, aPos);
2527 : Grid::EdgeAttachmentAccessor<TAPosition> aaSmoothBndPosOddVrt(mg, aSmoothBndPosOddVrt);
2528 :
2529 : // Manage vertex attachment communication in parallel case:
2530 : // - Setup communication policy for the above attachment aPosition
2531 : // - Setup interface communicator
2532 : // - Setup distributed grid manager
2533 : // - Setup grid layout map
2534 : #ifdef UG_PARALLEL
2535 : // Attachment communication policies COPY
2536 : ComPol_CopyAttachment<VertexLayout, TAPosition> comPolCopyAPosition(mg, aPos);
2537 :
2538 : // Interface communicators and distributed domain manager
2539 : pcl::InterfaceCommunicator<VertexLayout> com;
2540 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
2541 : GridLayoutMap& glm = dgm.grid_layout_map();
2542 : #endif
2543 :
2544 :
2545 : /*****************************************
2546 : *
2547 : * (2) DETERMINE aNumManifoldEdges
2548 : *
2549 : *****************************************/
2550 :
2551 0 : CalculateNumManifoldEdgesVertexAttachmentInParentLevel(mg, markSH, aNumManifoldEdges, false);
2552 :
2553 :
2554 : /*****************************************
2555 : *
2556 : * (3) CALCULATE aSmoothBndPosEvenVrt,
2557 : * aSmoothBndPosOddVrt
2558 : *
2559 : *****************************************/
2560 :
2561 : // Calculate aSmoothBndPosOddVrt
2562 0 : CalculateSmoothManifoldPosInParentLevelButterflyScheme(mg, aPos, markSH, linearManifoldSH, aSmoothBndPosOddVrt, aNumManifoldEdges);
2563 :
2564 :
2565 : /*****************************************
2566 : *
2567 : * (4) APPLY
2568 : *
2569 : *****************************************/
2570 :
2571 : // Loop all vertices of top_level
2572 0 : for(VertexIterator vrtIter = mg.begin<Vertex>(mg.top_level()); vrtIter != mg.end<Vertex>(mg.top_level()); ++vrtIter)
2573 : {
2574 : Vertex* vrt = *vrtIter;
2575 :
2576 : // Catch vertices without parent
2577 0 : if(mg.get_parent(vrt) == NULL)
2578 0 : continue;
2579 :
2580 : // In case of marked manifold vertices, which do not belong to the user-specified linear boundary manifold subsets,
2581 : // and activated Loop scheme refinement apply subdivision surfaces smoothing, else linear refinement
2582 0 : if(markSH.get_subset_index(vrt) != -1 && linearManifoldSH.get_subset_index(vrt) == -1)
2583 : {
2584 : // ODD VERTEX
2585 0 : if(mg.get_parent(vrt)->reference_object_id() == ROID_EDGE)
2586 : {
2587 : // Get parent edge
2588 : Edge* parentEdge = static_cast<Edge*>(mg.get_parent(vrt));
2589 :
2590 : aaPos[vrt] = aaSmoothBndPosOddVrt[parentEdge];
2591 : }
2592 : }
2593 : }
2594 :
2595 :
2596 : /*****************************************
2597 : *
2598 : * (5) COMMUNICATE VERTICALLY
2599 : * AFTER SUBDIVISION SURFACES
2600 : *
2601 : *****************************************/
2602 :
2603 : // Communicate aPosition in parallel case
2604 : #ifdef UG_PARALLEL
2605 : // copy ghosts = VMASTER to v_slaves
2606 : com.exchange_data(glm, INT_V_MASTER, INT_V_SLAVE, comPolCopyAPosition);
2607 : com.communicate();
2608 : #endif
2609 :
2610 :
2611 : /*****************************************
2612 : *
2613 : * (6) CLEAN UP
2614 : *
2615 : *****************************************/
2616 :
2617 : // detach vertex attachments
2618 : mg.detach_from_vertices(aNumManifoldEdges);
2619 : mg.detach_from_edges(aSmoothBndPosOddVrt);
2620 0 : }
2621 :
2622 :
2623 : ////////////////////////////////////////////////////////////////////////////////
2624 : template <class TAPosition>
2625 0 : void ApplySmoothManifoldPosToTopLevelAveragingScheme(MultiGrid& mg, TAPosition& aPos, MGSubsetHandler& markSH,
2626 : MGSubsetHandler& linearManifoldSH)
2627 : {
2628 : /*****************************************
2629 : *
2630 : * (1) SETUP
2631 : *
2632 : *****************************************/
2633 :
2634 : if(TAPosition::ValueType::Size == 1){
2635 : UG_THROW("Error in ApplySmoothManifoldPosToTopLevelAveragingScheme:\n"
2636 : "Currently only dimensions 2 and 3 are supported.\n");
2637 : }
2638 :
2639 : // Position attachment value type
2640 : typedef typename TAPosition::ValueType position_type;
2641 :
2642 : // Vertex attachments for associated number of manifold faces and smooth position
2643 : AInt aNumManifoldFaces_tri;
2644 : AInt aNumManifoldFaces_quad;
2645 : TAPosition aSmoothBndPos_tri;
2646 : TAPosition aSmoothBndPos_quad;
2647 :
2648 : // attach previously declared vertex attachments with initial value 0
2649 0 : mg.attach_to_vertices_dv(aNumManifoldFaces_tri, 0);
2650 0 : mg.attach_to_vertices_dv(aNumManifoldFaces_quad, 0);
2651 0 : mg.attach_to_vertices_dv(aSmoothBndPos_tri, position_type(0));
2652 0 : mg.attach_to_vertices_dv(aSmoothBndPos_quad, position_type(0));
2653 :
2654 : // Define attachment accessors
2655 : Grid::VertexAttachmentAccessor<TAPosition> aaPos(mg, aPos);
2656 : Grid::VertexAttachmentAccessor<AInt> aaNumManifoldFaces_tri(mg, aNumManifoldFaces_tri);
2657 : Grid::VertexAttachmentAccessor<AInt> aaNumManifoldFaces_quad(mg, aNumManifoldFaces_quad);
2658 : Grid::VertexAttachmentAccessor<TAPosition> aaSmoothBndPos_tri(mg, aSmoothBndPos_tri);
2659 : Grid::VertexAttachmentAccessor<TAPosition> aaSmoothBndPos_quad(mg, aSmoothBndPos_quad);
2660 :
2661 : // Manage vertex attachment communication in parallel case:
2662 : // - Setup communication policy for the above attachment aPosition
2663 : // - Setup interface communicator
2664 : // - Setup distributed grid manager
2665 : // - Setup grid layout map
2666 : #ifdef UG_PARALLEL
2667 : // Attachment communication policies COPY
2668 : ComPol_CopyAttachment<VertexLayout, TAPosition> comPolCopyAPosition(mg, aPos);
2669 :
2670 : // Interface communicators and distributed domain manager
2671 : pcl::InterfaceCommunicator<VertexLayout> com;
2672 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
2673 : GridLayoutMap& glm = dgm.grid_layout_map();
2674 : #endif
2675 :
2676 :
2677 : /*****************************************
2678 : *
2679 : * (2) DETERMINE aNumManifoldEdges
2680 : *
2681 : *****************************************/
2682 :
2683 0 : CalculateNumManifoldFacesVertexAttachmentInTopLevel(mg, markSH, aNumManifoldFaces_tri, aNumManifoldFaces_quad);
2684 :
2685 :
2686 : /*****************************************
2687 : *
2688 : * (3) CALCULATE
2689 : *
2690 : *****************************************/
2691 :
2692 : // Calculate aSmoothBndPosEvenVrt, aSmoothBndPosOddVrt
2693 0 : CalculateSmoothManifoldPosInTopLevelAveragingScheme(mg, aPos, markSH, linearManifoldSH, aSmoothBndPos_tri, aSmoothBndPos_quad);
2694 :
2695 :
2696 : /*****************************************
2697 : *
2698 : * (4) APPLY
2699 : *
2700 : *****************************************/
2701 :
2702 : // Move manifold vertices to their smoothed position
2703 0 : for(VertexIterator vrtIter = mg.begin<Vertex>(mg.top_level()); vrtIter != mg.end<Vertex>(mg.top_level()); ++vrtIter)
2704 : {
2705 : Vertex* vrt = *vrtIter;
2706 :
2707 : // In case of marked manifold vertices, which do not belong to the user-specified linear boundary manifold subsets,
2708 : // and activated averaging scheme apply subdivision surfaces smoothing, else linear refinement
2709 0 : if(markSH.get_subset_index(vrt) != -1 && linearManifoldSH.get_subset_index(vrt) == -1)
2710 : {
2711 0 : if(aaNumManifoldFaces_tri[vrt] == 0 && aaNumManifoldFaces_quad[vrt] == 0)
2712 0 : UG_THROW("ERROR in ApplySmoothManifoldPosToTopLevelAveragingScheme: grid contains manifold vertex not contained in any manifold face.");
2713 :
2714 : // Scale smooth vertex position by the number of associated volume elements (SubdivisionVolumes smoothing)
2715 0 : VecScale(aaSmoothBndPos_tri[vrt], aaSmoothBndPos_tri[vrt], 1.0/6.0/(aaNumManifoldFaces_tri[vrt]/6.0 + aaNumManifoldFaces_quad[vrt]/4.0));
2716 0 : VecScale(aaSmoothBndPos_quad[vrt], aaSmoothBndPos_quad[vrt], 1.0/4.0/(aaNumManifoldFaces_tri[vrt]/6.0 + aaNumManifoldFaces_quad[vrt]/4.0));
2717 : VecScaleAdd(aaPos[vrt], 1.0, aaSmoothBndPos_tri[vrt], 1.0, aaSmoothBndPos_quad[vrt]);
2718 : }
2719 : }
2720 :
2721 :
2722 : /*****************************************
2723 : *
2724 : * (5) COMMUNICATE VERTICALLY
2725 : * AFTER SUBDIVISION SURFACES
2726 : *
2727 : *****************************************/
2728 :
2729 : // Communicate aPosition in parallel case
2730 : #ifdef UG_PARALLEL
2731 : // copy v_slaves to ghosts = VMASTER
2732 : com.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, comPolCopyAPosition);
2733 : com.communicate();
2734 : #endif
2735 :
2736 :
2737 : /*****************************************
2738 : *
2739 : * (6) CLEAN UP
2740 : *
2741 : *****************************************/
2742 :
2743 : // detach vertex attachments
2744 : mg.detach_from_vertices(aNumManifoldFaces_tri);
2745 : mg.detach_from_vertices(aNumManifoldFaces_quad);
2746 : mg.detach_from_vertices(aSmoothBndPos_tri);
2747 : mg.detach_from_vertices(aSmoothBndPos_quad);
2748 0 : }
2749 :
2750 :
2751 : ////////////////////////////////////////////////////////////////////////////////
2752 0 : void ApplySmoothVolumePosToTopLevel(MultiGrid& mg, MGSubsetHandler& markSH,
2753 : bool bConstrained)
2754 : {
2755 : /*****************************************
2756 : *
2757 : * (1) SETUP
2758 : *
2759 : *****************************************/
2760 :
2761 : // Vertex attachments for associated number of elements and smooth position
2762 : AInt aNumElems_toc;
2763 : AInt aNumElems_prism;
2764 : AInt aNumElems_hex;
2765 : APosition aSmoothVolPos_toc;
2766 : APosition aSmoothVolPos_prism;
2767 : APosition aSmoothVolPos_hex;
2768 :
2769 : // attach previously declared vertex attachments with initial value 0
2770 0 : mg.attach_to_vertices_dv(aNumElems_toc, 0);
2771 0 : mg.attach_to_vertices_dv(aNumElems_prism, 0);
2772 0 : mg.attach_to_vertices_dv(aNumElems_hex, 0);
2773 0 : mg.attach_to_vertices_dv(aSmoothVolPos_toc, vector3(0, 0, 0));
2774 0 : mg.attach_to_vertices_dv(aSmoothVolPos_prism, vector3(0, 0, 0));
2775 0 : mg.attach_to_vertices_dv(aSmoothVolPos_hex, vector3(0, 0, 0));
2776 :
2777 : // Define attachment accessors
2778 : Grid::VertexAttachmentAccessor<APosition> aaPos(mg, aPosition);
2779 : Grid::VertexAttachmentAccessor<AInt> aaNumElems_toc(mg, aNumElems_toc);
2780 : Grid::VertexAttachmentAccessor<AInt> aaNumElems_prism(mg, aNumElems_prism);
2781 : Grid::VertexAttachmentAccessor<AInt> aaNumElems_hex(mg, aNumElems_hex);
2782 : Grid::VertexAttachmentAccessor<APosition> aaSmoothVolPos_toc(mg, aSmoothVolPos_toc);
2783 : Grid::VertexAttachmentAccessor<APosition> aaSmoothVolPos_prism(mg, aSmoothVolPos_prism);
2784 : Grid::VertexAttachmentAccessor<APosition> aaSmoothVolPos_hex(mg, aSmoothVolPos_hex);
2785 :
2786 : // Manage vertex attachment communication in parallel case:
2787 : // - Setup communication policy for the above attachment aPosition
2788 : // - Setup interface communicator
2789 : // - Setup distributed grid manager
2790 : // - Setup grid layout map
2791 : #ifdef UG_PARALLEL
2792 : // Attachment communication policies COPY
2793 : ComPol_CopyAttachment<VertexLayout, AVector3> comPolCopyAPosition(mg, aPosition);
2794 :
2795 : // Interface communicators and distributed domain manager
2796 : pcl::InterfaceCommunicator<VertexLayout> com;
2797 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
2798 : GridLayoutMap& glm = dgm.grid_layout_map();
2799 : #endif
2800 :
2801 :
2802 : /*****************************************
2803 : *
2804 : * (2) DETERMINE aNumElems
2805 : *
2806 : *****************************************/
2807 :
2808 0 : CalculateNumElemsVertexAttachmentInTopLevel(mg, aNumElems_toc, aNumElems_prism, aNumElems_hex);
2809 :
2810 :
2811 : /*****************************************
2812 : *
2813 : * (3) CALCULATE
2814 : *
2815 : *****************************************/
2816 :
2817 : // Calculate aSmoothVolPos
2818 0 : if(bConstrained)
2819 0 : CalculateConstrainedSmoothVolumePosInTopLevel(mg, markSH, aSmoothVolPos_toc);
2820 : else
2821 0 : CalculateSmoothVolumePosInTopLevel(mg, markSH, aSmoothVolPos_toc, aSmoothVolPos_prism, aSmoothVolPos_hex);
2822 :
2823 :
2824 : /*****************************************
2825 : *
2826 : * (4) APPLY
2827 : *
2828 : *****************************************/
2829 :
2830 : // Move vertices to their smoothed position
2831 0 : for(VertexIterator vrtIter = mg.begin<Vertex>(mg.top_level()); vrtIter != mg.end<Vertex>(mg.top_level()); ++vrtIter)
2832 : {
2833 : Vertex* vrt = *vrtIter;
2834 :
2835 0 : if(aaNumElems_toc[vrt] == 0 && aaNumElems_prism[vrt] == 0 && aaNumElems_hex[vrt] == 0)
2836 0 : UG_THROW("ERROR in ApplySmoothVolumePosToTopLevel: grid contains vertex not contained in any volume.");
2837 :
2838 0 : if(g_boundaryRefinementRule == SUBDIV_VOL)
2839 : {
2840 : // Scale smooth vertex position by the number of associated volume elements (SubdivisionVolumes smoothing)
2841 0 : VecScale(aaSmoothVolPos_toc[vrt], aaSmoothVolPos_toc[vrt], 1.0/14.0/(aaNumElems_toc[vrt]/14.0 + aaNumElems_prism[vrt]/12.0 + aaNumElems_hex[vrt]/8.0));
2842 0 : VecScale(aaSmoothVolPos_prism[vrt], aaSmoothVolPos_prism[vrt], 1.0/12.0/(aaNumElems_toc[vrt]/14.0 + aaNumElems_prism[vrt]/12.0 + aaNumElems_hex[vrt]/8.0));
2843 0 : VecScale(aaSmoothVolPos_hex[vrt], aaSmoothVolPos_hex[vrt], 1.0/8.0/(aaNumElems_toc[vrt]/14.0 + aaNumElems_prism[vrt]/12.0 + aaNumElems_hex[vrt]/8.0));
2844 : VecScaleAdd(aaPos[vrt], 1.0, aaSmoothVolPos_toc[vrt], 1.0, aaSmoothVolPos_prism[vrt], 1.0, aaSmoothVolPos_hex[vrt]);
2845 : }
2846 : else
2847 : {
2848 : // Only in case of inner vertices
2849 0 : if(markSH.get_subset_index(vrt) == -1)
2850 : {
2851 : // Scale smooth vertex position by the number of associated volume elements
2852 0 : VecScale(aaSmoothVolPos_toc[vrt], aaSmoothVolPos_toc[vrt], 1.0/14.0/(aaNumElems_toc[vrt]/14.0 + aaNumElems_prism[vrt]/12.0 + aaNumElems_hex[vrt]/8.0));
2853 0 : VecScale(aaSmoothVolPos_prism[vrt], aaSmoothVolPos_prism[vrt], 1.0/12.0/(aaNumElems_toc[vrt]/14.0 + aaNumElems_prism[vrt]/12.0 + aaNumElems_hex[vrt]/8.0));
2854 0 : VecScale(aaSmoothVolPos_hex[vrt], aaSmoothVolPos_hex[vrt], 1.0/8.0/(aaNumElems_toc[vrt]/14.0 + aaNumElems_prism[vrt]/12.0 + aaNumElems_hex[vrt]/8.0));
2855 : VecScaleAdd(aaPos[vrt], 1.0, aaSmoothVolPos_toc[vrt], 1.0, aaSmoothVolPos_prism[vrt], 1.0, aaSmoothVolPos_hex[vrt]);
2856 : }
2857 : }
2858 : }
2859 :
2860 :
2861 : /*****************************************
2862 : *
2863 : * (5) COMMUNICATE VERTICALLY
2864 : * AFTER SUBDIVISION VOLUMES
2865 : *
2866 : *****************************************/
2867 :
2868 : // Communicate aPosition in parallel case
2869 : #ifdef UG_PARALLEL
2870 : // copy v_slaves to ghosts = VMASTER
2871 : com.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, comPolCopyAPosition);
2872 : com.communicate();
2873 : #endif
2874 :
2875 :
2876 : /*****************************************
2877 : *
2878 : * (6) CLEAN UP
2879 : *
2880 : *****************************************/
2881 :
2882 : // detach vertex attachments
2883 : mg.detach_from_vertices(aNumElems_toc);
2884 : mg.detach_from_vertices(aNumElems_prism);
2885 : mg.detach_from_vertices(aNumElems_hex);
2886 : mg.detach_from_vertices(aSmoothVolPos_toc);
2887 : mg.detach_from_vertices(aSmoothVolPos_prism);
2888 : mg.detach_from_vertices(aSmoothVolPos_hex);
2889 0 : }
2890 :
2891 :
2892 : ////////////////////////////////////////////////////////////////////////////////
2893 : template <class TAPosition>
2894 0 : void ApplySmoothSubdivisionSurfacesToTopLevel(MultiGrid& mg, TAPosition& aPos, MGSubsetHandler& sh,
2895 : MGSubsetHandler& markSH, MGSubsetHandler& linearManifoldSH,
2896 : bool bCreaseSurf)
2897 : {
2898 : /*****************************************
2899 : *
2900 : * (1) SETUP
2901 : *
2902 : *****************************************/
2903 :
2904 : PROFILE_FUNC_GROUP("subdivision_volumes");
2905 :
2906 : if(TAPosition::ValueType::Size == 1){
2907 0 : UG_THROW("Error in ApplySmoothSubdivisionSurfacesToTopLevel:\n"
2908 : "Currently only dimensions 2 and 3 are supported.\n");
2909 : }
2910 :
2911 : // Catch use of procedure for MultiGrids with just one level
2912 : size_t numLevels = mg.num_levels();
2913 :
2914 : #ifdef UG_PARALLEL
2915 : if(pcl::NumProcs() > 1){
2916 : pcl::ProcessCommunicator pc;
2917 : numLevels = pc.allreduce(numLevels, PCL_RO_MAX);
2918 : }
2919 : #endif
2920 :
2921 0 : if(numLevels == 1)
2922 : {
2923 0 : UG_THROW("Error in ApplySmoothSubdivisionSurfacesToTopLevel: "
2924 : "Procedure only to be used for MultiGrids with more than one level.");
2925 : }
2926 :
2927 :
2928 : /*****************************************
2929 : *
2930 : * (2) SUBDIVISION SURFACES
2931 : *
2932 : *****************************************/
2933 :
2934 0 : if(bCreaseSurf){
2935 0 : if(!(g_boundaryRefinementRule == SUBDIV_SURF_LOOP_SCHEME || g_boundaryRefinementRule == LINEAR)){
2936 0 : UG_THROW("Error in ApplySmoothSubdivisionSurfacesToTopLevel: "
2937 : "Crease surfaces are currently only supported with 'subdiv_surf_loop_scheme' or 'linear' refinement.");
2938 : }
2939 : }
2940 :
2941 0 : if(g_boundaryRefinementRule == SUBDIV_SURF_LOOP_SCHEME)
2942 0 : ApplySmoothManifoldPosToTopLevelLoopScheme(mg, aPos, markSH, linearManifoldSH, bCreaseSurf);
2943 : else if(g_boundaryRefinementRule == SUBDIV_SURF_AVERAGING_SCHEME)
2944 0 : ApplySmoothManifoldPosToTopLevelAveragingScheme(mg, aPos, markSH, linearManifoldSH);
2945 : else if(g_boundaryRefinementRule == SUBDIV_SURF_BUTTERFLY_SCHEME)
2946 0 : ApplySmoothManifoldPosToTopLevelButterflyScheme(mg, aPos, markSH, linearManifoldSH);
2947 : else if(g_boundaryRefinementRule == SUBDIV_VOL){}
2948 : else if(g_boundaryRefinementRule == LINEAR){}
2949 : else
2950 0 : UG_THROW("ERROR in ApplySmoothSubdivisionSurfacesToTopLevel: Unknown boundary refinement rule. Known rules are 'subdiv_surf_loop_scheme', 'subdiv_surf_averaging_scheme', 'subdiv_surf_butterfly_scheme' or 'linear'.");
2951 0 : }
2952 :
2953 :
2954 : ////////////////////////////////////////////////////////////////////////////////
2955 0 : void ApplySmoothSubdivisionVolumesToTopLevel(MultiGrid& mg, MGSubsetHandler& sh, MGSubsetHandler& markSH,
2956 : MGSubsetHandler& linearManifoldSH, bool bConstrained)
2957 : {
2958 : /*****************************************
2959 : *
2960 : * (1) SETUP
2961 : *
2962 : *****************************************/
2963 :
2964 : PROFILE_FUNC_GROUP("subdivision_volumes");
2965 :
2966 : // Ensure, that hybrid tet-/oct refinement is used as refinement rule for tetrahedrons
2967 0 : if(tet_rules::GetRefinementRule() != tet_rules::HYBRID_TET_OCT)
2968 0 : UG_THROW("ERROR in ApplySubdivisionVolumesToTopLevel: Set necessary refinement rule by SetTetRefinementRule('hybrid_tet_oct').");
2969 :
2970 : // Catch use of procedure for MultiGrids with just one level
2971 : size_t numLevels = mg.num_levels();
2972 :
2973 : #ifdef UG_PARALLEL
2974 : if(pcl::NumProcs() > 1){
2975 : pcl::ProcessCommunicator pc;
2976 : numLevels = pc.allreduce(numLevels, PCL_RO_MAX);
2977 : }
2978 : #endif
2979 :
2980 0 : if(numLevels == 1)
2981 : {
2982 0 : UG_THROW("Error in ApplySmoothSubdivisionToTopLevel: "
2983 : "Procedure only to be used for MultiGrids with more than one level.");
2984 : }
2985 :
2986 :
2987 : /*****************************************
2988 : *
2989 : * (2) SUBDIVISION SURFACES
2990 : *
2991 : *****************************************/
2992 :
2993 0 : if(g_boundaryRefinementRule == SUBDIV_SURF_LOOP_SCHEME)
2994 0 : ApplySmoothManifoldPosToTopLevelLoopScheme(mg, aPosition, markSH, linearManifoldSH, false);
2995 : else if(g_boundaryRefinementRule == SUBDIV_SURF_AVERAGING_SCHEME)
2996 0 : ApplySmoothManifoldPosToTopLevelAveragingScheme(mg, aPosition, markSH, linearManifoldSH);
2997 : else if(g_boundaryRefinementRule == SUBDIV_SURF_BUTTERFLY_SCHEME)
2998 0 : ApplySmoothManifoldPosToTopLevelButterflyScheme(mg, aPosition, markSH, linearManifoldSH);
2999 : else if(g_boundaryRefinementRule == SUBDIV_VOL){}
3000 : else if(g_boundaryRefinementRule == LINEAR){}
3001 : else
3002 0 : UG_THROW("ERROR in ApplySubdivisionVolumesToTopLevel: Unknown boundary refinement rule. Known rules are 'subdiv_surf_loop_scheme', 'subdiv_surf_averaging_scheme', 'subdiv_surf_butterfly_scheme', 'linear' or 'subdiv_vol'.");
3003 :
3004 :
3005 : /*****************************************
3006 : *
3007 : * (3) SUBDIVISION VOLUMES
3008 : *
3009 : *****************************************/
3010 :
3011 0 : ApplySmoothVolumePosToTopLevel(mg, markSH, bConstrained);
3012 0 : }
3013 :
3014 :
3015 : //////////////////////////////////////////////////////////////////////////////
3016 : // Wrapper procedures
3017 : template <class TAPosition>
3018 0 : void ApplySmoothSubdivisionSurfacesToTopLevel(MultiGrid& mg, TAPosition& aPos, MGSubsetHandler& sh,
3019 : MGSubsetHandler& markSH, const char* linearManifoldSubsets,
3020 : bool bCreaseSurf)
3021 : {
3022 0 : MGSubsetHandler linearManifoldSH(mg);
3023 0 : InitLinearManifoldSubsetHandler(mg, sh, linearManifoldSH, linearManifoldSubsets);
3024 :
3025 0 : ApplySmoothSubdivisionSurfacesToTopLevel(mg, aPos, sh, markSH, linearManifoldSH, bCreaseSurf);
3026 0 : }
3027 :
3028 0 : void ApplySmoothSubdivisionVolumesToTopLevel(MultiGrid& mg, MGSubsetHandler& sh, MGSubsetHandler& markSH,
3029 : const char* linearManifoldSubsets)
3030 : {
3031 0 : MGSubsetHandler linearManifoldSH(mg);
3032 0 : InitLinearManifoldSubsetHandler(mg, sh, linearManifoldSH, linearManifoldSubsets);
3033 :
3034 0 : ApplySmoothSubdivisionVolumesToTopLevel(mg, sh, markSH, linearManifoldSH, false);
3035 0 : }
3036 :
3037 0 : void ApplyConstrainedSmoothSubdivisionVolumesToTopLevel(MultiGrid& mg, MGSubsetHandler& sh, MGSubsetHandler& markSH,
3038 : const char* linearManifoldSubsets)
3039 : {
3040 0 : MGSubsetHandler linearManifoldSH(mg);
3041 0 : InitLinearManifoldSubsetHandler(mg, sh, linearManifoldSH, linearManifoldSubsets);
3042 :
3043 0 : ApplySmoothSubdivisionVolumesToTopLevel(mg, sh, markSH, linearManifoldSH, true);
3044 0 : }
3045 :
3046 :
3047 : //////////////////////////////////////////////////////////////////////////////
3048 : // Explicit instantiations
3049 : template void ApplySmoothSubdivisionSurfacesToTopLevel<APosition1>(MultiGrid& mg, APosition1& aPos, MGSubsetHandler& sh,
3050 : MGSubsetHandler& markSH, const char* linearManifoldSubsets,
3051 : bool bCreaseSurf);
3052 : template void ApplySmoothSubdivisionSurfacesToTopLevel<APosition2>(MultiGrid& mg, APosition2& aPos, MGSubsetHandler& sh,
3053 : MGSubsetHandler& markSH, const char* linearManifoldSubsets,
3054 : bool bCreaseSurf);
3055 : template void ApplySmoothSubdivisionSurfacesToTopLevel<APosition>(MultiGrid& mg, APosition& aPos, MGSubsetHandler& sh,
3056 : MGSubsetHandler& markSH, const char* linearManifoldSubsets,
3057 : bool bCreaseSurf);
3058 :
3059 : template void ProjectHierarchyToSubdivisionLimit(MultiGrid& mg, APosition1& aPos);
3060 : template void ProjectHierarchyToSubdivisionLimit(MultiGrid& mg, APosition2& aPos);
3061 : template void ProjectHierarchyToSubdivisionLimit(MultiGrid& mg, APosition& aPos);
3062 :
3063 :
3064 :
3065 : }// end of namespace
|
