Line data Source code
1 : /*
2 : * Copyright (c) 2011-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include <vector>
34 : #include <string>
35 : #include <sstream>
36 :
37 : // include bridge
38 : #include "bindings/lua/lua_user_data.h"
39 : #include "bridge/bridge.h"
40 : #include "bridge/util.h"
41 : #include "bridge/util_domain_dependent.h"
42 : #include "lib_disc/domain.h"
43 : #include "lib_disc/domain_traits.h"
44 : #include "lib_grid/lib_grid.h"
45 : //todo: include this in algorithms.h
46 : #include "lib_grid/algorithms/refinement_mark_util.h"
47 : #include "lib_grid/refinement/adaptive_regular_mg_refiner.h"
48 : #include "lib_grid/refinement/global_fractured_media_refiner.h"
49 : #include "lib_grid/refinement/global_multi_grid_refiner.h"
50 : #include "lib_grid/refinement/global_subdivision_multi_grid_refiner.h"
51 : #include "lib_grid/refinement/hanging_node_refiner_multi_grid.h"
52 : #include "lib_grid/refinement/ref_mark_adjusters/horizontal_anisotropy_adjuster.h"
53 : #include "lib_grid/refinement/ref_mark_adjusters/shadow_copy_adjuster.h"
54 : #include "lib_grid/algorithms/subdivision/subdivision_volumes.h"
55 : #include "lib_grid/grid_objects/tetrahedron_rules.h"
56 : #ifdef UG_PARALLEL
57 : #include "lib_grid/parallelization/util/attachment_operations.hpp"
58 : #include "lib_grid/parallelization/parallel_refinement/parallel_refinement.h"
59 : #endif
60 :
61 : #include "common/math/misc/math_util.h"
62 :
63 : #include "lib_grid/refinement/refiner_factory.hpp"
64 :
65 : using namespace std;
66 :
67 : namespace ug{
68 :
69 : /**
70 : * \defgroup refinement_bridge Refinement Bridge
71 : * \ingroup domain_bridge
72 : * \{
73 : */
74 :
75 : ////////////////////////////////////////////////////////////////////////////////
76 : /// Returns a refinement mark, given a string describing it.
77 : /** Valid parameters are: "refine", "coarsen"*/
78 0 : static RefinementMark StringToRefinementMark(std::string markType)
79 : {
80 0 : TrimString(markType);
81 : std::transform(markType.begin(), markType.end(), markType.begin(), ::tolower);
82 0 : if(markType == "closure") return RM_CLOSURE;
83 0 : if(markType == "local") return RM_LOCAL;
84 0 : if(markType == "refine") return RM_REFINE;
85 0 : if(markType == "coarsen") return RM_COARSEN;
86 :
87 :
88 0 : UG_THROW("StringToRefinementMark: non-supported type: "<<markType);
89 : return RM_NONE;
90 : }
91 :
92 :
93 : ////////////////////////////////////////////////////////////////////////////////
94 : /// Factory method, that creates a global subdivision domain refiner.
95 : /** Automatically chooses whether a parallel refiner is required.*/
96 : template <class TDomain>
97 : static SmartPtr<GlobalSubdivisionMultiGridRefiner<typename TDomain::position_attachment_type> >
98 0 : GlobalSubdivisionDomainRefiner(TDomain* dom, const char* linearManifoldSubsets, bool constrained)
99 : {
100 : typedef typename TDomain::position_attachment_type position_attachment_type;
101 :
102 : // Check if mandatory manifold marks exist in the given domain grid
103 : try{
104 0 : dom->additional_subset_handler("markSH");
105 : }
106 0 : UG_CATCH_THROW("Error in CreateGlobalSubdivisionDomainRefiner:\n"
107 : "Given multigrid does not contain surface manifold mark SubsetHandler 'markSH'.\n"
108 : "Please mark all surface manifold elements in the current multigrid.");
109 :
110 : // Linear manifold subset handler management: check for optional specification by user
111 : MGSubsetHandler* linearManifoldSH = NULL;
112 0 : if(linearManifoldSubsets[0] != '\0'){
113 : try{
114 0 : dom->additional_subset_handler("linearManifoldSH");
115 : }
116 : // Linear manifold subsets specified but no corresponding SubsetHandler user-provided
117 0 : UG_CATCH_THROW("Error in CreateGlobalSubdivisionDomainRefiner:\n"
118 : "Linear manifold subsets specified, but no corresponding SubsetHandler 'linearManifoldSH' provided.");
119 : // Assignment of user-provided linear manifold SubsetHandler
120 0 : linearManifoldSH = &(*dom->additional_subset_handler("linearManifoldSH"));
121 : }
122 : else{
123 : // Assignment of empty dummy for linear manifold SubsetHandler
124 0 : dom->create_additional_subset_handler("linearManifoldSH");
125 0 : linearManifoldSH = &(*dom->additional_subset_handler("linearManifoldSH"));
126 : }
127 :
128 : // Instantiation of GlobalSubdivisionMultiGridRefiner object
129 : GlobalSubdivisionMultiGridRefiner<position_attachment_type>* ref = NULL;
130 : #ifdef UG_PARALLEL
131 : if(pcl::NumProcs() > 1){
132 : ref = new ParallelGlobalSubdivisionRefiner<position_attachment_type>(
133 : *dom->distributed_grid_manager(),
134 : dom->refinement_projector());
135 :
136 : ref->assign_subset_handler(*dom->subset_handler());
137 : ref->assign_mark_subset_handler(*dom->additional_subset_handler("markSH"));
138 : ref->assign_position_attachment(dom->position_attachment());
139 : }
140 : #endif
141 :
142 : if(!ref){
143 0 : ref = new GlobalSubdivisionMultiGridRefiner<position_attachment_type>(
144 : *dom->grid(), dom->position_attachment(),
145 0 : *dom->subset_handler(),
146 0 : *dom->additional_subset_handler("markSH"),
147 : dom->refinement_projector());
148 : }
149 :
150 0 : ref->set_linear_manifold_subsets(*linearManifoldSH, linearManifoldSubsets);
151 : ref->set_constrained_subdivision(constrained);
152 :
153 0 : return SmartPtr<GlobalSubdivisionMultiGridRefiner<position_attachment_type> >(ref);
154 : }
155 :
156 : ////////////////////////////////////////////////////////////////////////////////
157 : /// Creates an adaptive hanging node domain refiner.
158 : /** Automatically chooses whether a parallel refiner is required.*/
159 : template <typename TDomain>
160 0 : static SmartPtr<IRefiner> HangingNodeDomainRefiner(TDomain* dom)
161 : {
162 0 : if(!dom->is_adaptive()){
163 0 : UG_THROW("Can't create an adaptive refiner for the given domain. "
164 : "Construct the domain with isAdaptive enabled.");
165 : }
166 :
167 : //todo: support normal grids, too!
168 : #ifdef UG_PARALLEL
169 : if(pcl::NumProcs() > 1){
170 : return SmartPtr<IRefiner>(
171 : new ParallelHangingNodeRefiner_MultiGrid(
172 : *dom->distributed_grid_manager(),
173 : dom->refinement_projector()));
174 : }
175 : #endif
176 :
177 : return SmartPtr<IRefiner>(
178 0 : new HangingNodeRefiner_MultiGrid(
179 : *dom->grid(),
180 0 : dom->refinement_projector()));
181 : }
182 :
183 : ////////////////////////////////////////////////////////////////////////////////
184 : /// Creates an adaptive regular domain refiner.
185 : /** Automatically chooses whether a parallel refiner is required.*/
186 : template <typename TDomain>
187 0 : static SmartPtr<IRefiner> CreateAdaptiveRegularDomainRefiner(TDomain* dom)
188 : {
189 0 : if(!dom->is_adaptive()){
190 0 : UG_THROW("Can't create an adaptive refiner for the given domain. "
191 : "Construct the domain with isAdaptive enabled.");
192 : }
193 :
194 : //todo: support normal grids, too!
195 : //todo: support parallelism, too!
196 : /*#ifdef UG_PARALLEL
197 : if(pcl::NumProcs() > 1){
198 : return SmartPtr<IRefiner>(new ParallelHangingNodeRefiner_MultiGrid(*dom->distributed_grid_manager()));
199 : }
200 : #endif*/
201 :
202 : return SmartPtr<IRefiner>(
203 0 : new AdaptiveRegularRefiner_MultiGrid(
204 : *dom->grid(),
205 0 : dom->refinement_projector()));
206 : }
207 :
208 : ////////////////////////////////////////////////////////////////////////////////
209 : /// Creates a global fractured domain refiner.
210 : template <class TDomain>
211 : static SmartPtr<GlobalFracturedMediaRefiner>
212 0 : CreateGlobalFracturedDomainRefiner(TDomain* dom)
213 : {
214 0 : if(!dom->is_adaptive()){
215 0 : UG_THROW("Can't create an fractured domain refiner for the given domain. "
216 : "Construct the domain with isAdaptive enabled.");
217 : }
218 :
219 : GlobalFracturedMediaRefiner* ref = NULL;
220 : #ifdef UG_PARALLEL
221 : if(pcl::NumProcs() > 1){
222 : ref = new ParallelGlobalFracturedMediaRefiner(
223 : *dom->distributed_grid_manager(),
224 : dom->refinement_projector());
225 : }
226 : #endif
227 :
228 : if(!ref)
229 0 : ref = new GlobalFracturedMediaRefiner(
230 : *dom->grid(),
231 : dom->refinement_projector());
232 :
233 0 : ref->set_subset_handler(*dom->subset_handler());
234 :
235 0 : return SmartPtr<GlobalFracturedMediaRefiner>(ref);
236 : }
237 :
238 :
239 : /// Adds a horizontal-anisotropy-adjuster to the given refiner.
240 : /** Adds a horizontal-anisotropy-adjuster to the given refiner. The provided
241 : * domain is used to obtain the used position attachment.
242 : *
243 : * \note ref has to be derived from HangingNodeRefinerBase*/
244 : template <class TDomain>
245 : static SPIRefMarkAdjuster
246 0 : AddHorizontalAnisotropyAdjuster(IRefiner* ref, TDomain* dom)
247 : {
248 0 : HangingNodeRefiner_MultiGrid* href = dynamic_cast<HangingNodeRefiner_MultiGrid*>(ref);
249 0 : UG_COND_THROW(!href, "A horizontal anisotropy adjuster can only be added to an instance of HangingNodeRefiner_MultiGrid.");
250 :
251 : SmartPtr<HorizontalAnisotropyAdjuster<
252 : typename TDomain::position_attachment_type> >
253 0 : haa = make_sp(new HorizontalAnisotropyAdjuster<
254 : typename TDomain::position_attachment_type>
255 : (dom->position_attachment()));
256 :
257 0 : href->add_ref_mark_adjuster(haa);
258 0 : return haa;
259 : }
260 :
261 :
262 : /**
263 : * @brief Adds a ShadowCopyAdjuster to the given refiner
264 : * @param ref the refiner; has to be derived from HangingNodeRefiner_MultiGrid
265 : */
266 : static SPIRefMarkAdjuster
267 0 : AddShadowCopyAdjuster(IRefiner* ref)
268 : {
269 0 : HangingNodeRefiner_MultiGrid* href = dynamic_cast<HangingNodeRefiner_MultiGrid*>(ref);
270 0 : UG_COND_THROW(!href, "A ShadowCopyAdjuster can only be added to an instance of HangingNodeRefiner_MultiGrid.");
271 :
272 0 : SmartPtr<ShadowCopyAdjuster> sca = make_sp(new ShadowCopyAdjuster());
273 :
274 0 : href->add_ref_mark_adjuster(sca);
275 0 : return sca;
276 : }
277 :
278 :
279 :
280 : // ////////////////////////////////////////////////////////////////////////////////
281 : // /// marks face for anisotropic refinement, if it contains edges below given sizeRatio
282 : // /**
283 : // * marks face and for anisotropic refinement, if it contains edges
284 : // * below given sizeRatio. These edges are also marked.
285 : // * @return true, if face has been marked for anisotropic refinement.
286 : // * This is the case, when one of its edges has been marked for refinement.*/
287 : // template <class TAAPos> bool MarkIfAnisotropic(Face* f, IRefiner* refiner, number sizeRatio, TAAPos& aaPos)
288 : // {
289 : // bool marked = false;
290 : // uint num_edges = f->num_edges();
291 : // vector<Edge*> edges(num_edges);
292 : // // collect associated edges
293 : // CollectAssociated(edges, *refiner->grid(), f);
294 :
295 : // // find the shortest edge
296 : // Edge* minEdge = FindShortestEdge(edges.begin(), edges.end(), aaPos);
297 : // UG_ASSERT(minEdge,
298 : // "Associated edges of each face have to exist at this point.");
299 : // number minLen = EdgeLength(minEdge, aaPos);
300 :
301 : // // compare all associated edges of f against minEdge (even minEdge itself,
302 : // // if somebody sets edgeRatio to 1 or higher)
303 : // for (uint i_edge = 0; i_edge < num_edges; ++i_edge) {
304 : // Edge* e = edges[i_edge];
305 : // number len = EdgeLength(e, aaPos);
306 : // // to avoid division by 0, we only consider edges with a length > 0
307 : // if(len > 0) {
308 : // if(minLen / len <= sizeRatio) {
309 : // // the edge will be refined
310 : // refiner->mark(e);
311 : // marked = true;
312 : // }
313 : // }
314 : // }
315 :
316 : // if(marked) {
317 : // // if a vertex has been marked, also mark the face, or else just a hanging
318 : // // node would be inserted.
319 : // // We do not mark other associated objects here since this would
320 : // // cause creation of a closure and would also behave differently
321 : // // in a parallel environment, compared to a serial environment.
322 : // // By using RM_ANISOTROPIC, we avoid that associated edges of
323 : // // the marked face will automatically be marked, too.
324 : // refiner->mark(f, RM_ANISOTROPIC);
325 : // }
326 :
327 : // return marked;
328 : // }
329 :
330 : ////////////////////////////////////////////////////////////////////////////////
331 : /// Marks all elements from refinement.
332 : /** If used in a parallel environment only elements on the calling procs
333 : * are marked.
334 : */
335 0 : static void MarkForRefinement_All(SmartPtr<IRefiner> ref)
336 : {
337 : PROFILE_FUNC_GROUP("grid");
338 0 : Grid* g = ref->get_associated_grid();
339 0 : if(!g){
340 : UG_LOG("Refiner is not registered at a grid. Aborting.\n");
341 0 : return;
342 : }
343 : ref->mark(g->vertices_begin(), g->vertices_end());
344 : ref->mark(g->edges_begin(), g->edges_end());
345 : ref->mark(g->faces_begin(), g->faces_end());
346 : ref->mark(g->volumes_begin(), g->volumes_end());
347 : }
348 :
349 : template <class TElem>
350 0 : static void MarkForRefinement_AllAnisotropic(SmartPtr<IRefiner> ref)
351 : {
352 : PROFILE_FUNC_GROUP("grid");
353 0 : Grid* g = ref->get_associated_grid();
354 0 : if(!g){
355 : UG_LOG("Refiner is not registered at a grid. Aborting.\n");
356 0 : return;
357 : }
358 : ref->mark(g->begin<TElem>(), g->end<TElem>(), RM_ANISOTROPIC);
359 : }
360 :
361 : ////////////////////////////////////////////////////////////////////////////////
362 : /// Marks all vertices in the given d-dimensional sphere.
363 : template <class TDomain>
364 0 : void MarkForAdaption_VerticesInSphereMaxLvl(TDomain& dom, SmartPtr<IRefiner> refiner,
365 : const typename TDomain::position_type& center,
366 : number radius, std::string markType, int maxLvl)
367 : {
368 : PROFILE_FUNC_GROUP("grid");
369 : typedef typename TDomain::position_accessor_type position_accessor_type;
370 :
371 0 : RefinementMark rmMark = StringToRefinementMark(markType);
372 :
373 : // make sure that the refiner was created for the given domain
374 0 : if(refiner->get_associated_grid() != dom.grid().get()){
375 0 : throw(UGError("ERROR in MarkForRefinement_VerticesInSphere: "
376 : "Refiner was not created for the specified domain. Aborting."));
377 : }
378 :
379 : typedef typename TDomain::grid_type TGrid;
380 :
381 0 : Grid& grid = *refiner->get_associated_grid();
382 0 : TGrid& g = *dom.grid();
383 : position_accessor_type& aaPos = dom.position_accessor();
384 :
385 : // we'll compare against the square radius.
386 0 : number radiusSq = radius * radius;
387 :
388 : // we'll store associated edges, faces and volumes in those containers
389 : vector<Edge*> vEdges;
390 : vector<Face*> vFaces;
391 : vector<Volume*> vVols;
392 :
393 : // iterate over all vertices of the grid. If a vertex is inside the given sphere,
394 : // then we'll mark all associated elements.
395 : for(VertexIterator iter = grid.begin<Vertex>();
396 0 : iter != grid.end<Vertex>(); ++iter)
397 : {
398 : int lvl = g.get_level(*iter);
399 0 : if(maxLvl >= 0 && lvl >= maxLvl)
400 0 : continue;
401 :
402 0 : if(VecDistanceSq(center, aaPos[*iter]) <= radiusSq){
403 : CollectAssociated(vEdges, grid, *iter);
404 : CollectAssociated(vFaces, grid, *iter);
405 : CollectAssociated(vVols, grid, *iter);
406 :
407 : refiner->mark(vEdges.begin(), vEdges.end(), rmMark);
408 : refiner->mark(vFaces.begin(), vFaces.end(), rmMark);
409 : refiner->mark(vVols.begin(), vVols.end(), rmMark);
410 : }
411 : }
412 0 : }
413 :
414 : template <class TDomain>
415 0 : void MarkForRefinement_VerticesInSphere(TDomain& dom, SmartPtr<IRefiner> refiner,
416 : const typename TDomain::position_type& center,
417 : number radius)
418 : {
419 0 : MarkForAdaption_VerticesInSphereMaxLvl(dom, refiner, center, radius, "refine", -1);
420 0 : }
421 :
422 : template <class TDomain>
423 0 : void MarkForAdaption_VerticesInSphere(TDomain& dom, SmartPtr<IRefiner> refiner,
424 : const typename TDomain::position_type& center,
425 : number radius, std::string markType)
426 : {
427 0 : MarkForAdaption_VerticesInSphereMaxLvl(dom, refiner, center, radius, markType, -1);
428 0 : }
429 :
430 : ////////////////////////////////////////////////////////////////////////////////
431 : /// Marks all elements which lie completely in the given d-dimensional sphere.
432 : /** Valid parameters for TElem are Edge, Face, Volume.*/
433 : template <class TDomain, class TElem>
434 0 : void MarkForRefinement_ElementsInSphere(TDomain& dom, SmartPtr<IRefiner> refiner,
435 : const typename TDomain::position_type& center,
436 : number radius)
437 : {
438 : PROFILE_FUNC_GROUP("grid");
439 : typedef typename TDomain::position_accessor_type position_accessor_type;
440 : typedef typename geometry_traits<TElem>::iterator ElemIter;
441 :
442 : // make sure that the refiner was created for the given domain
443 0 : if(refiner->get_associated_grid() != dom.grid().get()){
444 0 : throw(UGError("ERROR in MarkForRefinement_VerticesInCube: "
445 : "Refiner was not created for the specified domain. Aborting."));
446 : }
447 :
448 0 : Grid& grid = *refiner->get_associated_grid();
449 : position_accessor_type& aaPos = dom.position_accessor();
450 :
451 : // we'll compare against the square radius.
452 0 : number radiusSq = radius * radius;
453 :
454 : // we'll store associated edges, faces and volumes in those containers
455 : vector<Edge*> vEdges;
456 : vector<Face*> vFaces;
457 : vector<Volume*> vVols;
458 :
459 : // iterate over all elements of the grid. If all vertices of an element are inside
460 : // the given sphere, then we'll mark those elements.
461 : for(ElemIter iter = grid.begin<TElem>();
462 0 : iter != grid.end<TElem>(); ++iter)
463 : {
464 : // get element
465 : TElem* elem = *iter;
466 :
467 : // bool flag to check whether all vertices are in the sphere
468 : bool bInSphere = true;
469 :
470 : // loop all vertices
471 0 : for(size_t i = 0; i < elem->num_vertices(); ++i)
472 : {
473 : // check if vertex is in sphere
474 0 : if(VecDistanceSq(center, aaPos[elem->vertex(i)]) > radiusSq)
475 : bInSphere = false;
476 : }
477 :
478 : // mark the element
479 0 : if(bInSphere)
480 0 : refiner->mark(elem);
481 : }
482 0 : }
483 :
484 : ////////////////////////////////////////////////////////////////////////////////
485 : /// Marks all elements which have vertices in the given d-dimensional cube.
486 : /** Make sure that TAPos is an attachment of vector_t position types.*/
487 : template <class TDomain>
488 0 : void MarkForAdaption_VerticesInCube(TDomain& dom, SmartPtr<IRefiner> refiner,
489 : const typename TDomain::position_type& min,
490 : const typename TDomain::position_type& max,
491 : std::string markType)
492 : {
493 : PROFILE_FUNC_GROUP("grid");
494 : typedef typename TDomain::position_type position_type;
495 : typedef typename TDomain::position_accessor_type position_accessor_type;
496 :
497 0 : RefinementMark rmMark = StringToRefinementMark(markType);
498 :
499 : // make sure that the refiner was created for the given domain
500 0 : if(refiner->get_associated_grid() != dom.grid().get()){
501 0 : throw(UGError("ERROR in MarkForRefinement_VerticesInCube: "
502 : "Refiner was not created for the specified domain. Aborting."));
503 : }
504 :
505 0 : Grid& grid = *refiner->get_associated_grid();
506 : position_accessor_type& aaPos = dom.position_accessor();
507 :
508 : // we'll store associated edges, faces and volumes in those containers
509 : vector<Edge*> vEdges;
510 : vector<Face*> vFaces;
511 : vector<Volume*> vVols;
512 :
513 : // iterate over all vertices of the grid. If a vertex is inside the given cube,
514 : // then we'll mark all associated elements.
515 : for(VertexIterator iter = grid.begin<Vertex>();
516 0 : iter != grid.end<Vertex>(); ++iter)
517 : {
518 : // Position
519 : position_type& pos = aaPos[*iter];
520 :
521 : // check flag
522 : bool bRefine = true;
523 :
524 : // check node
525 0 : for(size_t d = 0; d < pos.size(); ++d)
526 0 : if(pos[d] < min[d] || max[d] < pos[d])
527 : bRefine = false;
528 :
529 0 : if(bRefine)
530 : {
531 : CollectAssociated(vEdges, grid, *iter);
532 : CollectAssociated(vFaces, grid, *iter);
533 : CollectAssociated(vVols, grid, *iter);
534 :
535 : refiner->mark(vEdges.begin(), vEdges.end(), rmMark);
536 : refiner->mark(vFaces.begin(), vFaces.end(), rmMark);
537 : refiner->mark(vVols.begin(), vVols.end(), rmMark);
538 : }
539 : }
540 0 : }
541 :
542 : template <class TDomain>
543 0 : void MarkForRefinement_VerticesInCube(TDomain& dom, SmartPtr<IRefiner> refiner,
544 : const typename TDomain::position_type& min,
545 : const typename TDomain::position_type& max)
546 : {
547 0 : MarkForAdaption_VerticesInCube(dom, refiner, min, max, "refine");
548 0 : }
549 :
550 :
551 : /// Marks all elements for anisotropic refienment and also marks all edges > minLen.
552 : template <class TDomain>
553 0 : void MarkAnisotropic_LongEdges(TDomain& dom, IRefiner& refiner, number minLen)
554 : {
555 : UG_ASSERT(dom.grid().get() == refiner.get_associated_grid(),
556 : "Grids in domain and in refiner have to match!");
557 :
558 : typedef typename domain_traits<TDomain::dim>::element_type elem_t;
559 : typedef typename MultiGrid::traits<elem_t>::iterator iter_t;
560 :
561 : typename TDomain::position_accessor_type aaPos = dom.position_accessor();
562 0 : MultiGrid& mg = *dom.grid();
563 : MultiGrid::edge_traits::secure_container edges;
564 : MultiGrid::face_traits::secure_container faces;
565 :
566 : number minLenSq = sq(minLen);
567 :
568 0 : for(iter_t e_iter = mg.begin<elem_t>(); e_iter != mg.end<elem_t>(); ++e_iter){
569 : elem_t* elem = *e_iter;
570 0 : if(mg.has_children(elem))
571 0 : continue;
572 :
573 : // we'll mark all volumes as anisotropic, since we currently have to use
574 : // copy elements during anisotropic refinement.
575 0 : refiner.mark(elem, RM_ANISOTROPIC);
576 :
577 0 : mg.associated_elements(edges, elem);
578 :
579 0 : for(size_t i = 0; i < edges.size(); ++i){
580 0 : if(EdgeLengthSq(edges[i], aaPos) >= minLenSq){
581 0 : refiner.mark(edges[i]);
582 : }
583 : }
584 : }
585 0 : }
586 :
587 :
588 : // ////////////////////////////////////////////////////////////////////////////////
589 : // /// Marks the long edges in anisotropic faces and faces with a big area in anisotropic volumes.
590 : // /**
591 : // * THE NOT HIGHLY SUCCESSFUL IMPLEMENTATION OF THIS METHOD LEAD TO NEW INSIGHT.
592 : // * A NEW ANISOTROPIC REFINEMENT STRATEGY WILL BE IMPLEMENTED, WHICH WILL ALSO
593 : // * LEAD TO A REIMPLEMENTATION OF THIS METHOD. BEST NOT TO USE IT UNTIL THEN.
594 : // *
595 : // * The sizeRatio determines Whether a face or a volume is considered anisotropic.
596 : // * Make sure that the sizeRatio is in the interval [0, 1]. If the
597 : // * ratio of the shortest edge to an other edge falls below the given threshold,
598 : // * then the associated face is considered anisotropic and the longer edge will be
599 : // * marked. The face itself will then be marked for anisotropic refinement.
600 : // * The same technique is applied to volumes, this time however the ratio between
601 : // * face-areas is considered when judging whether a volume is anisotropic.
602 : // *
603 : // * VOLUME MARKS ARE CURRENTLY DISABLED
604 : // *
605 : // * Note that this algorithm only really works for a serial environment.
606 : // * \todo improve it so that it works in parallel, too. A specialization of
607 : // * HangingNodeRefiner may be required for this.
608 : // *
609 : // * \todo activate and improve volume marks
610 : // **/
611 : // template <class TDomain>
612 : // void MarkForRefinement_AnisotropicElements(TDomain& dom, SmartPtr<IRefiner> refiner,
613 : // number sizeRatio)
614 : // {
615 : // PROFILE_FUNC_GROUP("grid");
616 : // typedef typename TDomain::position_accessor_type position_accessor_type;
617 :
618 : // // make sure that the refiner was created for the given domain
619 : // if(refiner->get_associated_grid() != dom.grid().get()){
620 : // throw(UGError("ERROR in MarkForRefinement_VerticesInCube: "
621 : // "Refiner was not created for the specified domain. Aborting."));
622 : // }
623 :
624 : // // access the grid and the position attachment
625 : // Grid& grid = *refiner->get_associated_grid();
626 : // position_accessor_type& aaPos = dom.position_accessor();
627 :
628 : // // If the grid is a multigrid, we want to avoid marking of elements, which do
629 : // // have children
630 : // MultiGrid* pmg = dynamic_cast<MultiGrid*>(&grid);
631 :
632 : // // make sure that the grid automatically generates sides for each element
633 : // if(!grid.option_is_enabled(GRIDOPT_AUTOGENERATE_SIDES)){
634 : // UG_LOG("WARNING in MarkForRefinement_AnisotropicElements: "
635 : // "Enabling GRIDOPT_AUTOGENERATE_SIDES.\n");
636 : // grid.enable_options(GRIDOPT_AUTOGENERATE_SIDES);
637 : // }
638 :
639 : // // we'll store associated edges and faces in those containers
640 : // vector<Edge*> edges;
641 : // vector<Face*> faces;
642 :
643 : // // iterate over all faces of the grid.
644 : // for(FaceIterator iter = grid.begin<Face>();
645 : // iter != grid.end<Face>(); ++iter)
646 : // {
647 : // Face* f = *iter;
648 :
649 : // // ignore faces with children
650 : // if(pmg && pmg->has_children(f))
651 : // continue;
652 :
653 : // // collect associated edges
654 : // CollectAssociated(edges, grid, f);
655 :
656 : // // find the shortest edge
657 : // Edge* minEdge = FindShortestEdge(edges.begin(), edges.end(), aaPos);
658 : // UG_ASSERT(minEdge, "Associated edges of each face have to exist at this point.");
659 : // number minLen = EdgeLength(minEdge, aaPos);
660 :
661 : // // compare all associated edges of f against minEdge (even minEdge itself,
662 : // // if somebody sets edgeRatio to 1 or higher)
663 : // for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
664 : // Edge* e = edges[i_edge];
665 : // number len = EdgeLength(e, aaPos);
666 : // // to avoid division by 0, we only consider edges with a length > 0
667 : // if(len > 0){
668 : // if(minLen / len <= sizeRatio){
669 : // // the edge will be refined
670 : // refiner->mark(e);
671 :
672 : // // we'll also mark the current face, or else just a hanging
673 : // // node would be inserted.
674 : // // We do not mark other associated objects here since this would
675 : // // cause creation of a closure and would also behave differently
676 : // // in a parallel environment, compared to a serial environment.
677 : // // By using RM_ANISOTROPIC, we avoid that associated edges of
678 : // // the marked face will automatically be marked, too.
679 : // refiner->mark(f, RM_ANISOTROPIC);
680 : // }
681 : // }
682 : // }
683 : // }
684 :
685 : // // iterate over all faces again. We have to make sure that faces which have
686 : // // a marked short edge are refined regular.
687 : // // first push all marked faces into a queue
688 : // // we're using grid::mark to make sure that each face lies only once on the queue.
689 : // // Grid::mark has nothing to do with refinement. It is just a helper for us.
690 : // grid.begin_marking();
691 :
692 : // queue<Face*> queFaces;
693 : // for(FaceIterator iter = grid.begin<Face>(); iter != grid.end<Face>(); ++iter){
694 : // queFaces.push(*iter);
695 : // grid.mark(*iter);
696 : // }
697 :
698 : // while(!queFaces.empty()){
699 : // Face* f = queFaces.front();
700 : // queFaces.pop();
701 :
702 : // if(pmg && pmg->has_children(f)){
703 : // grid.unmark(f);
704 : // continue;
705 : // }
706 :
707 : // // collect associated edges
708 : // CollectAssociated(edges, grid, f);
709 : // /*
710 : // if(refiner->get_mark(f) == RM_NONE){
711 : // bool gotOne = false;
712 : // for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
713 : // Edge* e = edges[i_edge];
714 : // if(refiner->get_mark(e) != RM_NONE){
715 : // gotOne = true;
716 : // break;
717 : // }
718 : // }
719 :
720 : // if(gotOne){
721 : // for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
722 : // Edge* e = edges[i_edge];
723 : // if(refiner->get_mark(e) == RM_NONE){
724 : // refiner->mark(e);
725 : // CollectFaces(faces, grid, e);
726 : // for(size_t i_face = 0; i_face < faces.size(); ++i_face){
727 : // Face* nbr = faces[i_face];
728 : // if(!grid.is_marked(nbr)
729 : // && (refiner->get_mark(nbr) == RM_ANISOTROPIC))
730 : // {
731 : // grid.mark(nbr);
732 : // queFaces.push(nbr);
733 : // }
734 : // }
735 : // }
736 : // }
737 : // refiner->mark(f);
738 : // }
739 : // }
740 : // else */if(refiner->get_mark(f) == RM_ANISOTROPIC){
741 : // // find the shortest edge
742 : // Edge* minEdge = FindShortestEdge(edges.begin(), edges.end(), aaPos);
743 : // UG_ASSERT(minEdge, "Associated edges of each face have to exist at this point.");
744 : // number minLen = EdgeLength(minEdge, aaPos);
745 :
746 : // // check if a short edge and a long edge is selected
747 : // bool longEdgeSelected = false;
748 : // bool shortEdgeSelected = false;
749 :
750 : // for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
751 : // Edge* e = edges[i_edge];
752 : // if(refiner->get_mark(e) == RM_NONE)
753 : // continue;
754 :
755 : // number len = EdgeLength(e, aaPos);
756 : // // to avoid division by 0, we only consider edges with a length > 0
757 : // if(len > 0){
758 : // if(minLen / len <= sizeRatio){
759 : // longEdgeSelected = true;
760 : // }
761 : // else{
762 : // shortEdgeSelected = true;
763 : // }
764 : // }
765 : // }
766 :
767 : // // if a short edge and a long edge was selected, we'll have to mark all
768 : // // edges and push associated faces with anisotropic mark to the queue
769 : // if(longEdgeSelected && shortEdgeSelected){
770 : // for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
771 : // Edge* e = edges[i_edge];
772 : // if(refiner->get_mark(e) == RM_NONE){
773 : // // mark it and push associated anisotropic faces to the queue
774 : // refiner->mark(e);
775 : // //!!!
776 :
777 : // if(ConstrainingEdge::type_match(e)){
778 : // UG_LOG("CONSTRAINING EDGE MARKED (2)\n");
779 : // }
780 :
781 : // CollectFaces(faces, grid, e);
782 : // for(size_t i_face = 0; i_face < faces.size(); ++i_face){
783 : // Face* nbr = faces[i_face];
784 : // if(!grid.is_marked(nbr)
785 : // && (refiner->get_mark(nbr) == RM_ANISOTROPIC))
786 : // {
787 : // grid.mark(nbr);
788 : // queFaces.push(nbr);
789 : // }
790 : // }
791 : // }
792 : // }
793 : // }
794 : // }
795 : // // now unmark the face
796 : // grid.unmark(f);
797 : // }
798 :
799 : // grid.end_marking();
800 :
801 : // // mark unmarked neighbors of marked edges for regular refinement
802 : // /*
803 : // for(FaceIterator iter = grid.begin<Face>();
804 : // iter != grid.end<Face>(); ++iter)
805 : // {
806 : // Face* f = *iter;
807 :
808 : // // if it is already marked, leave it as it is
809 : // if(refiner->get_mark(f) != RM_NONE)
810 : // continue;
811 :
812 : // // ignore faces with children
813 : // if(pmg && pmg->has_children(f))
814 : // continue;
815 :
816 : // for(size_t i = 0; i < f->num_edges(); ++i){
817 : // if(refiner->get_mark(grid.get_side(f, i)) != RM_NONE)
818 : // refiner->mark(f);
819 : // }
820 : // }*/
821 :
822 : // /*
823 : // // now that all faces are marked, we can process volumes. We consider a
824 : // // volume which has an anisotropic side as an anisotropic volume
825 : // for(VolumeIterator iter = grid.begin<Volume>();
826 : // iter != grid.end<Volume>(); ++iter)
827 : // {
828 : // Volume* v = *iter;
829 :
830 : // // collect associated faces
831 : // CollectAssociated(faces, grid, v);
832 :
833 : // // find the smallest face
834 : // Face* minFace = FindSmallestFace(faces.begin(), faces.end(), aaPos);
835 : // UG_ASSERT(minFace, "Associated faces of each volume have to exist at this point.");
836 : // number minArea = FaceArea(minFace, aaPos);
837 :
838 : // // compare all associated faces of v against minArea
839 : // for(size_t i_face = 0; i_face < faces.size(); ++i_face){
840 : // Face* f = faces[i_face];
841 : // number area = FaceArea(f, aaPos);
842 : // // avoid division by 0
843 : // if(area > 0){
844 : // if(minArea / area <= sizeRatio){
845 : // // the face will be refined. If it is already marked, we'll
846 : // // leave it at that, to keep the anisotropy.
847 : // // If it wasn't marked, we'll mark it for full refinement
848 : // // (all anisotropic faces have already been marked).
849 : // if(refiner->get_mark(f) == RM_NONE)
850 : // refiner->mark(f);
851 :
852 : // // the volume itself now has to be marked, too.
853 : // refiner->mark(v, RM_ANISOTROPIC);
854 : // }
855 : // }
856 : // }
857 : // }*/
858 : // }
859 :
860 : // ///////
861 : // /**
862 : // *
863 : // */
864 : // template <class TDomain>
865 : // void MarkForRefinement_AnisotropicElements2(TDomain& dom, SmartPtr<IRefiner> refiner,
866 : // number sizeRatio)
867 : // {
868 : // PROFILE_FUNC_GROUP("grid");
869 : // typedef typename TDomain::position_accessor_type position_accessor_type;
870 :
871 : // // make sure that the refiner was created for the given domain
872 : // if(refiner->get_associated_grid() != dom.grid().get()){
873 : // throw(UGError("ERROR in MarkForRefinement_VerticesInCube: "
874 : // "Refiner was not created for the specified domain. Aborting."));
875 : // }
876 :
877 : // // access the grid and the position attachment
878 : // Grid& grid = *refiner->get_associated_grid();
879 : // position_accessor_type& aaPos = dom.position_accessor();
880 : // IRefiner* ref = refiner.get_nonconst();
881 :
882 : // // If the grid is a multigrid, we want to avoid marking of elements, which do
883 : // // have children
884 : // MultiGrid* pmg = dynamic_cast<MultiGrid*>(&grid);
885 :
886 : // // make sure that the grid automatically generates sides for each element
887 : // if(!grid.option_is_enabled(GRIDOPT_AUTOGENERATE_SIDES)){
888 : // UG_LOG("WARNING in MarkForRefinement_AnisotropicElements: "
889 : // "Enabling GRIDOPT_AUTOGENERATE_SIDES.\n");
890 : // grid.enable_options(GRIDOPT_AUTOGENERATE_SIDES);
891 : // }
892 :
893 : // // we'll store associated edges, faces and volumes in those containers
894 : // vector<Edge*> edges;
895 : // vector<Face*> faces;
896 : // vector<Volume*> volumes;
897 : // // iterate over all faces of the grid.
898 : // for(FaceIterator iter = grid.begin<Face>();
899 : // iter != grid.end<Face>(); ++iter)
900 : // {
901 : // Face* f = *iter;
902 : // // ignore faces with children
903 : // if(pmg && pmg->has_children(f))
904 : // continue;
905 :
906 : // // if faces has been marked, store it for later marking of its neighbours
907 : // if(MarkIfAnisotropic(f, ref, sizeRatio, aaPos))
908 : // faces.push_back(f);
909 : // else // fixme: mark for regular refinement should not be needed!
910 : // refiner->mark(f);
911 : // }
912 :
913 : // // if a face is marked for anisotropic refinement (AR),
914 : // // also mark associated volumes for AR
915 : // for(vector<Face*>::iterator iter = faces.begin(); iter != faces.end(); ++iter)
916 : // CollectVolumes(volumes, grid, *iter, false);
917 :
918 : // refiner->mark(volumes.begin(), volumes.end(), RM_ANISOTROPIC);
919 : // }
920 :
921 : ////////////////////////////////////////////////////////////////////////////////
922 : /// Marks all elements which have vertices in the given d-dimensional cube.
923 : /** Make sure that TAPos is an attachment of vector_t position types.*/
924 : template <class TDomain>
925 0 : void AssignSubset_VerticesInCube(TDomain& dom,
926 : const typename TDomain::position_type& min,
927 : const typename TDomain::position_type& max, int si)
928 : {
929 : typedef typename TDomain::position_type position_type;
930 : typedef typename TDomain::position_accessor_type position_accessor_type;
931 :
932 0 : Grid& grid = *dom.grid();
933 : position_accessor_type& aaPos = dom.position_accessor();
934 0 : MGSubsetHandler& sh = *dom.subset_handler();
935 :
936 : // iterate over all vertices of the grid. If a vertex is inside the given cube,
937 : // then we'll mark all associated elements.
938 : for(VertexIterator iter = grid.begin<Vertex>();
939 0 : iter != grid.end<Vertex>(); ++iter)
940 : {
941 : // Position
942 : position_type& pos = aaPos[*iter];
943 :
944 : // check flag
945 : bool bInside = true;
946 :
947 : // check node
948 0 : for(size_t d = 0; d < pos.size(); ++d)
949 0 : if(pos[d] < min[d] || max[d] < pos[d])
950 : bInside = false;
951 :
952 0 : if(bInside)
953 : {
954 0 : sh.assign_subset(*iter, si);
955 : }
956 : }
957 0 : }
958 :
959 : ////////////////////////////////////////////////////////////////////////////////
960 : /// Marks all elements which have vertices in the given d-dimensional cube.
961 : /** Make sure that TAPos is an attachment of vector_t position types.*/
962 : template <class TDomain>
963 0 : void AssignSubset_VerticesInSphere(TDomain& dom,
964 : const typename TDomain::position_type& center,
965 : const number radius, int si)
966 : {
967 : typedef typename TDomain::position_type position_type;
968 : typedef typename TDomain::position_accessor_type position_accessor_type;
969 :
970 0 : Grid& grid = *dom.grid();
971 : position_accessor_type& aaPos = dom.position_accessor();
972 0 : MGSubsetHandler& sh = *dom.subset_handler();
973 :
974 : // iterate over all vertices of the grid. If a vertex is inside the given cube,
975 : // then we'll mark all associated elements.
976 : for(VertexIterator iter = grid.begin<Vertex>();
977 0 : iter != grid.end<Vertex>(); ++iter)
978 : {
979 : // Position
980 : position_type& pos = aaPos[*iter];
981 :
982 : // check flag
983 : bool bInside = true;
984 :
985 : // check node
986 : number normSq = 0.0;
987 0 : for(int d = 0; d < TDomain::dim; ++d){
988 0 : normSq += (pos[d] - center[d])*(pos[d] - center[d]);
989 : }
990 0 : number dist = sqrt( normSq );
991 :
992 0 : if(dist >= radius)
993 : bInside = false;
994 :
995 : if(bInside)
996 : {
997 0 : sh.assign_subset(*iter, si);
998 : }
999 : }
1000 0 : }
1001 :
1002 :
1003 : #ifdef UG_DIM_1
1004 0 : static number DistanceToSurfaceDomain(MathVector<1>& tpos, const Vertex* vrt,
1005 : Grid::VertexAttachmentAccessor<Attachment<MathVector<1> > >& aaPosSurf)
1006 : {
1007 0 : UG_THROW("Not implemented");
1008 : }
1009 : #endif
1010 :
1011 : #ifdef UG_DIM_2
1012 0 : static number DistanceToSurfaceDomain(MathVector<2>& tpos, const Edge* edge,
1013 : Grid::VertexAttachmentAccessor<Attachment<MathVector<2> > >& aaPosSurf)
1014 : {
1015 0 : return DistancePointToLine(tpos , aaPosSurf[edge->vertex(0)], aaPosSurf[edge->vertex(1)]);
1016 : }
1017 : #endif
1018 :
1019 : #ifdef UG_DIM_3
1020 0 : static number DistanceToSurfaceDomain(MathVector<3>& tpos, const Face* face,
1021 : Grid::VertexAttachmentAccessor<Attachment<MathVector<3> > >& aaPosSurf)
1022 : {
1023 :
1024 : MathVector<3> a, b, n;
1025 0 : VecSubtract(a, aaPosSurf[face->vertex(1)], aaPosSurf[face->vertex(0)]);
1026 0 : VecSubtract(b, aaPosSurf[face->vertex(2)], aaPosSurf[face->vertex(0)]);
1027 0 : VecCross(n, a,b);
1028 :
1029 : number bc1Out, bc2Out;
1030 0 : return DistancePointToTriangle( a, bc1Out, bc2Out,
1031 : tpos,
1032 0 : aaPosSurf[face->vertex(0)], aaPosSurf[face->vertex(1)], aaPosSurf[face->vertex(2)], n);
1033 : }
1034 : #endif
1035 :
1036 : template <class TDomain>
1037 0 : void MarkForRefinement_CloseToSurface(TDomain& dom, SmartPtr<IRefiner> refiner,
1038 : double time, size_t maxLvl,
1039 : const TDomain& surfaceDomain)
1040 : {
1041 : PROFILE_FUNC();
1042 : typedef typename TDomain::grid_type TGrid;
1043 : //typedef typename TDomain::subset_handler_type TSubsetHandler;
1044 : typedef typename TDomain::position_type TPos;
1045 : typedef typename TDomain::position_accessor_type TAAPos;
1046 : typedef typename domain_traits<TDomain::dim>::element_type TElem;
1047 : typedef typename domain_traits<TDomain::dim>::side_type TSide;
1048 : typedef typename TGrid::template traits<TElem>::iterator TIter;
1049 :
1050 0 : TGrid& g = *dom.grid();
1051 : //TSubsetHandler& sh = *dom.subset_handler();
1052 : TAAPos aaPos = dom.position_accessor();
1053 : TAAPos aaPosSurf = surfaceDomain.position_accessor();
1054 0 : const TGrid& surface = *surfaceDomain.grid();
1055 : Grid::edge_traits::secure_container edges;
1056 :
1057 0 : for(TIter iter = g.template begin<TElem>(); iter != g.template end<TElem>(); ++iter)
1058 : {
1059 : TElem* e = *iter;
1060 0 : size_t lvl = g.get_level(e);
1061 0 : if(lvl >= maxLvl)
1062 0 : continue;
1063 :
1064 0 : if(!g.has_children(e)){
1065 : // calculate the element center and the length of the longest edge
1066 0 : TPos tpos = CalculateCenter(e, aaPos);
1067 : vector3 pos;
1068 : VecCopy(pos, tpos, 0);
1069 :
1070 :
1071 : // the minimal edge defines h for this element
1072 : number h = numeric_limits<number>::max();
1073 0 : g.associated_elements(edges, e);
1074 0 : for(size_t i = 0; i < edges.size(); ++i){
1075 0 : number l = EdgeLengthSq(edges[i], aaPos);
1076 0 : if(l < h){
1077 : h = l;
1078 : }
1079 : }
1080 0 : h = sqrt(h);
1081 :
1082 :
1083 : int refine = 0;
1084 :
1085 : for( typename TGrid::template traits<TSide>::const_iterator it = surface.template begin<TSide>();
1086 0 : it != surface.template end<TSide>(); ++it){
1087 : const TSide* side = *it;
1088 :
1089 0 : number distance = DistanceToSurfaceDomain(tpos, side, aaPosSurf);
1090 :
1091 0 : if(distance < h)
1092 : refine = 1;
1093 : }
1094 :
1095 0 : if(refine)
1096 0 : refiner->mark(e);
1097 : }
1098 : }
1099 0 : }
1100 :
1101 :
1102 : template <class TDomain>
1103 0 : void MarkForRefinement_ContainsSurfaceNode(TDomain& dom, SmartPtr<IRefiner> refiner,
1104 : double time, size_t maxLvl,
1105 : const TDomain& surfaceDomain)
1106 : {
1107 : PROFILE_FUNC();
1108 : typedef typename TDomain::grid_type TGrid;
1109 : //typedef typename TDomain::subset_handler_type TSubsetHandler;
1110 : //typedef typename TDomain::position_type TPos;
1111 : typedef typename TDomain::position_accessor_type TAAPos;
1112 : typedef typename domain_traits<TDomain::dim>::element_type TElem;
1113 : typedef typename domain_traits<TDomain::dim>::side_type TSide;
1114 : typedef typename TGrid::template traits<TElem>::iterator TIter;
1115 :
1116 0 : TGrid& g = *dom.grid();
1117 : //TSubsetHandler& sh = *dom.subset_handler();
1118 : TAAPos aaPos = dom.position_accessor();
1119 : TAAPos aaPosSurf = surfaceDomain.position_accessor();
1120 0 : const TGrid& surface = *surfaceDomain.grid();
1121 : Grid::edge_traits::secure_container edges;
1122 :
1123 0 : for(TIter iter = g.template begin<TElem>(); iter != g.template end<TElem>(); ++iter)
1124 : {
1125 : TElem* e = *iter;
1126 0 : size_t lvl = g.get_level(e);
1127 0 : if(lvl >= maxLvl)
1128 0 : continue;
1129 :
1130 0 : if(!g.has_children(e)){
1131 : // calculate the element center and the length of the longest edge
1132 :
1133 : int refine = 0;
1134 :
1135 : for( typename TGrid::template traits<TSide>::const_iterator it = surface.template begin<TSide>();
1136 0 : it != surface.template end<TSide>(); ++it){
1137 : TSide* side = const_cast<TSide*>(*it);
1138 :
1139 0 : for(size_t co = 0; co < NumVertices(side); ++co){
1140 :
1141 0 : if(ContainsPoint(e, aaPosSurf[GetVertex(side,co)], aaPos))
1142 : refine = 1;
1143 : }
1144 : }
1145 :
1146 0 : if(refine)
1147 0 : refiner->mark(e);
1148 : }
1149 : }
1150 0 : }
1151 :
1152 :
1153 :
1154 : template <class TDomain>
1155 0 : void MarkForRefinement_ElementsByLuaCallback(TDomain& dom, SmartPtr<IRefiner> refiner,
1156 : double time, size_t maxLvl,
1157 : const char* luaCallbackName)
1158 : {
1159 : PROFILE_FUNC();
1160 : typedef typename TDomain::grid_type TGrid;
1161 : typedef typename TDomain::subset_handler_type TSubsetHandler;
1162 : typedef typename TDomain::position_type TPos;
1163 : typedef typename TDomain::position_accessor_type TAAPos;
1164 : typedef typename domain_traits<TDomain::dim>::element_type TElem;
1165 : typedef typename TGrid::template traits<TElem>::iterator TIter;
1166 :
1167 0 : TGrid& g = *dom.grid();
1168 0 : TSubsetHandler& sh = *dom.subset_handler();
1169 : TAAPos aaPos = dom.position_accessor();
1170 : // Grid::edge_traits::secure_container edges;
1171 :
1172 0 : LuaFunction<int, number> callback;
1173 : // we'll pass the following arguments: x, y, z, /*h,*/ lvl, si, time
1174 0 : callback.set_lua_callback(luaCallbackName, /*7*/6);
1175 :
1176 0 : for(TIter iter = g.template begin<TElem>(); iter != g.template end<TElem>(); ++iter)
1177 : {
1178 : TElem* e = *iter;
1179 0 : size_t lvl = g.get_level(e);
1180 0 : if(lvl >= maxLvl)
1181 0 : continue;
1182 :
1183 0 : if(!g.has_children(e)){
1184 : // calculate the element center and the length of the longest edge
1185 0 : TPos tpos = CalculateCenter(e, aaPos);
1186 : vector3 pos;
1187 : VecCopy(pos, tpos, 0);
1188 :
1189 : // note: this part had been inserted, but the lua callbacks rely on 6 arguments
1190 : // if you need the mesh size in the lua callback you may add it AT THE END of the parameter list
1191 : // HOWEVER, please contact the author of this file BEFORE inserting it
1192 : /*
1193 : // the longest edge defines h for this element
1194 : number h = numeric_limits<number>::max();
1195 : g.associated_elements(edges, e);
1196 : for(size_t i = 0; i < edges.size(); ++i){
1197 : number l = EdgeLengthSq(edges[i], aaPos);
1198 : if(l < h){
1199 : h = l;
1200 : }
1201 : }
1202 : h = sqrt(h);
1203 : */
1204 :
1205 0 : int refine = 0;
1206 0 : callback(refine, 6, pos.x(), pos.y(), pos.z(), /*h,*/ (number)lvl,
1207 0 : (number)sh.get_subset_index(e), (number)time);
1208 0 : if(refine)
1209 0 : refiner->mark(e);
1210 : }
1211 : }
1212 0 : }
1213 :
1214 : template <class TDomain>
1215 0 : void MarkForCoarsen_ElementsByLuaCallback(TDomain& dom, SmartPtr<IRefiner> refiner,
1216 : double time, const char* luaCallbackName)
1217 : {
1218 : PROFILE_FUNC();
1219 0 : if(!refiner->coarsening_supported()){
1220 : UG_LOG("WARNING in MarkForCoarsen_ElementsByLuaCallback: "
1221 : "Refiner doesn't support coarsening!\n");
1222 0 : return;
1223 : }
1224 :
1225 : typedef typename TDomain::grid_type TGrid;
1226 : typedef typename TDomain::subset_handler_type TSubsetHandler;
1227 : typedef typename TDomain::position_type TPos;
1228 : typedef typename TDomain::position_accessor_type TAAPos;
1229 : typedef typename domain_traits<TDomain::dim>::element_type TElem;
1230 : typedef typename TGrid::template traits<TElem>::iterator TIter;
1231 :
1232 0 : TGrid& g = *dom.grid();
1233 0 : TSubsetHandler& sh = *dom.subset_handler();
1234 : TAAPos aaPos = dom.position_accessor();
1235 :
1236 0 : LuaFunction<int, number> callback;
1237 : // we'll pass the following arguments: x, y, z, lvl, si, time
1238 0 : callback.set_lua_callback(luaCallbackName, 6);
1239 :
1240 0 : for(TIter iter = g.template begin<TElem>(); iter != g.template end<TElem>(); ++iter)
1241 : {
1242 : TElem* e = *iter;
1243 0 : if(!g.has_children(e)){
1244 0 : int coarsen = 0;
1245 0 : TPos tpos = CalculateCenter(e, aaPos);
1246 : vector3 pos;
1247 : VecCopy(pos, tpos, 0);
1248 0 : callback(coarsen, 6, pos.x(), pos.y(), pos.z(), (number)g.get_level(e),
1249 0 : (number)sh.get_subset_index(e), (number)time);
1250 0 : if(coarsen){
1251 0 : refiner->mark(e, RM_COARSEN);
1252 : }
1253 : }
1254 : }
1255 : }
1256 :
1257 :
1258 :
1259 : template <class TDomain, class TSubsetHandler, class TElem>
1260 0 : void MarkForAdaption_ElementsInSubset_(TDomain& dom, IRefiner& refiner,
1261 : TSubsetHandler& sh, int subsetIndex, RefinementMark rmMark)
1262 : {
1263 : PROFILE_FUNC();
1264 : typedef typename GridObjectCollection::traits<TElem>::iterator iterator_t;
1265 :
1266 0 : GridObjectCollection goc = sh.get_grid_objects_in_subset(subsetIndex);
1267 :
1268 0 : for(size_t lvl = 0; lvl < goc.num_levels(); ++lvl){
1269 : for(iterator_t iter = goc.begin<TElem>(lvl);
1270 0 : iter != goc.end<TElem>(lvl); ++iter)
1271 : {
1272 0 : refiner.mark(*iter, rmMark);
1273 : }
1274 : }
1275 0 : }
1276 :
1277 : template <class TDomain, class TSubsetHandler, class TElem>
1278 0 : void MarkForAdaption_ElementsInSubset(TDomain& dom, IRefiner& refiner,
1279 : TSubsetHandler& sh, int subsetIndex, std::string markType)
1280 : {
1281 0 : RefinementMark rmMark = StringToRefinementMark(markType);
1282 0 : MarkForAdaption_ElementsInSubset_<TDomain,TSubsetHandler,TElem>(dom, refiner, sh, subsetIndex, rmMark);
1283 0 : }
1284 :
1285 :
1286 : template <class TDomain, class TSubsetHandler, class TElem>
1287 0 : void MarkForRefinement_ElementsInSubset(TDomain& dom, IRefiner& refiner,
1288 : TSubsetHandler& sh, int subsetIndex)
1289 : {
1290 0 : MarkForAdaption_ElementsInSubset_<TDomain,TSubsetHandler,TElem>(dom, refiner, sh, subsetIndex, RM_REFINE);
1291 0 : }
1292 :
1293 :
1294 : template <class TDomain, class TElem>
1295 : void MarkForAdaption_ElementsContainingPoint(TDomain& dom, IRefiner& refiner,
1296 : number x, number y, number z,
1297 : std::string markType)
1298 : {
1299 : PROFILE_FUNC();
1300 : typedef typename TDomain::grid_type TGrid;
1301 : typedef typename TDomain::position_type TPos;
1302 : typedef typename TDomain::position_accessor_type TAAPos;
1303 : typedef typename TGrid::template traits<TElem>::iterator TIter;
1304 :
1305 : RefinementMark rmMark = StringToRefinementMark(markType);
1306 :
1307 : vector3 p3(x, y, z);
1308 : TPos p;
1309 : VecCopy(p, p3, 0);
1310 :
1311 : TGrid& g = *dom.grid();
1312 : TAAPos aaPos = dom.position_accessor();
1313 :
1314 : for(TIter iter = g.template begin<TElem>();
1315 : iter != g.template end<TElem>(); ++iter)
1316 : {
1317 : if(ContainsPoint(*iter, p, aaPos))
1318 : refiner.mark(*iter, rmMark);
1319 : }
1320 : }
1321 :
1322 :
1323 : template <class TDomain>
1324 0 : void MarkForAdaption_ElementsTouchingSubset(TDomain& dom, IRefiner& refiner,
1325 : ISubsetHandler& sh, int subsetIndex,
1326 : std::string markType)
1327 : {
1328 : PROFILE_FUNC();
1329 : typedef typename TDomain::grid_type TGrid;
1330 : typedef typename domain_traits<TDomain::dim>::element_type TElem;
1331 : typedef typename TGrid::template traits<TElem>::iterator TIter;
1332 :
1333 0 : RefinementMark rmMark = StringToRefinementMark(markType);
1334 :
1335 0 : TGrid& g = *dom.grid();
1336 :
1337 : Grid::vertex_traits::secure_container vrts;
1338 : Grid::edge_traits::secure_container edges;
1339 : Grid::face_traits::secure_container faces;
1340 :
1341 : for(TIter iter = g.template begin<TElem>();
1342 0 : iter != g.template end<TElem>(); ++iter)
1343 : {
1344 : TElem* e = *iter;
1345 : bool gotOne = false;
1346 : if(VERTEX < TElem::BASE_OBJECT_ID){
1347 0 : g.associated_elements(vrts, e);
1348 0 : for(size_t i = 0; i < vrts.size(); ++i){
1349 0 : if(sh.get_subset_index(vrts[i]) == subsetIndex){
1350 0 : refiner.mark(e, rmMark);
1351 : gotOne = true;
1352 : break;
1353 : }
1354 : }
1355 0 : if(gotOne)
1356 0 : continue;
1357 : }
1358 : if(EDGE < TElem::BASE_OBJECT_ID){
1359 : g.associated_elements(edges, e);
1360 0 : for(size_t i = 0; i < edges.size(); ++i){
1361 0 : if(sh.get_subset_index(edges[i]) == subsetIndex){
1362 0 : refiner.mark(e, rmMark);
1363 : gotOne = true;
1364 : break;
1365 : }
1366 : }
1367 0 : if(gotOne)
1368 0 : continue;
1369 : }
1370 : if(FACE < TElem::BASE_OBJECT_ID){
1371 : g.associated_elements(faces, e);
1372 0 : for(size_t i = 0; i < faces.size(); ++i){
1373 0 : if(sh.get_subset_index(faces[i]) == subsetIndex){
1374 0 : refiner.mark(e, rmMark);
1375 : gotOne = true;
1376 : break;
1377 : }
1378 : }
1379 0 : if(gotOne)
1380 0 : continue;
1381 : }
1382 : }
1383 0 : }
1384 :
1385 : template <class TDomain>
1386 0 : void MarkForAdaption_ElementsTouchingSubsets(TDomain& dom, IRefiner& refiner,
1387 : const char* subsets,
1388 : std::string markType)
1389 : {
1390 0 : ISubsetHandler& sh = *dom.subset_handler();
1391 0 : vector<string> vSubsets = TokenizeTrimString(std::string(subsets));
1392 0 : for(size_t i = 0; i < vSubsets.size(); ++i){
1393 0 : const int si = sh.get_subset_index(vSubsets[i].c_str());
1394 0 : if(si < 0)
1395 0 : UG_THROW("MarkForAdaption_ElementsTouchingSubsets: Subset '"<<vSubsets[i]<<"' not found");
1396 :
1397 0 : MarkForAdaption_ElementsTouchingSubset(dom, refiner, sh, si, markType);
1398 : }
1399 0 : }
1400 :
1401 :
1402 : template <class elem_t>
1403 0 : void MarkForRefinement_SubsetInterfaceElements(IRefiner& refiner, ISubsetHandler& sh)
1404 : {
1405 : // ALGORITHM
1406 : // - store max subset index of associated higher-dim elements in elem_t.
1407 : // - AllReduce (max)
1408 : // - Check if an associated higher-dim elem exists with a different subset index
1409 : // - if so: set attachment value to -2
1410 : // - AllReduce (min)
1411 : // - Mark all elems with attachment value -2
1412 : //
1413 : // NOTE: a serial implementation could do without the attachment by iterating
1414 : // over associated higher-dim elements of elem_t and directly checking
1415 : // for different subsets!
1416 :
1417 : AInt aSI;
1418 0 : UG_COND_THROW(!refiner.grid(), "A grid has to be associated with the given refiner.");
1419 0 : UG_COND_THROW(refiner.grid() != sh.grid(), "Grid of supplied refiner and subset handler have to match");
1420 0 : Grid& grid = *refiner.grid();
1421 :
1422 0 : grid.attach_to<elem_t> (aSI);
1423 : Grid::AttachmentAccessor<elem_t, AInt> aaSI(grid, aSI);
1424 :
1425 : SetAttachmentValues (aaSI, grid.begin<elem_t>(), grid.end<elem_t>(), -1);
1426 :
1427 : typedef typename elem_t::sideof hdelem_t;
1428 : typedef typename Grid::traits<elem_t>::iterator iter_t;
1429 : typedef typename Grid::traits<hdelem_t>::iterator hditer_t;
1430 :
1431 : typename Grid::traits<elem_t>::secure_container sides;
1432 :
1433 : for(hditer_t i_hd = grid.begin<hdelem_t>();
1434 0 : i_hd != grid.end<hdelem_t>(); ++i_hd)
1435 : {
1436 : hdelem_t* elem = *i_hd;
1437 0 : int si = sh.get_subset_index(elem);
1438 : grid.associated_elements(sides, elem);
1439 :
1440 0 : for(size_t i_side = 0; i_side < sides.size(); ++i_side){
1441 : elem_t* side = sides[i_side];
1442 0 : aaSI[side] = max(aaSI[side], si);
1443 : }
1444 : }
1445 :
1446 : #ifdef UG_PARALLEL
1447 : AttachmentAllReduce<elem_t>(grid, aSI, PCL_RO_MAX);
1448 : #endif
1449 :
1450 : for(hditer_t i_hd = grid.begin<hdelem_t>();
1451 0 : i_hd != grid.end<hdelem_t>(); ++i_hd)
1452 : {
1453 : hdelem_t* elem = *i_hd;
1454 : int si = sh.get_subset_index(elem);
1455 : grid.associated_elements(sides, elem);
1456 :
1457 0 : for(size_t i_side = 0; i_side < sides.size(); ++i_side){
1458 : elem_t* side = sides[i_side];
1459 0 : if(aaSI[side] != si){
1460 0 : aaSI[side] = -2;
1461 : }
1462 : }
1463 : }
1464 :
1465 : #ifdef UG_PARALLEL
1466 : AttachmentAllReduce<elem_t>(grid, aSI, PCL_RO_MIN);
1467 : #endif
1468 :
1469 : for(iter_t i_elem = grid.begin<elem_t>();
1470 0 : i_elem != grid.end<elem_t>(); ++i_elem)
1471 : {
1472 0 : if(aaSI[*i_elem] == -2)
1473 0 : refiner.mark(*i_elem);
1474 : }
1475 :
1476 : grid.detach_from<elem_t> (aSI);
1477 0 : }
1478 :
1479 : template <class TDomain, class TElem>
1480 0 : void MarkForRefinement_SubsetInterfaceElements(TDomain& dom, IRefiner& refiner)
1481 : {
1482 0 : MarkForRefinement_SubsetInterfaceElements<TElem>(refiner, *dom.subset_handler());
1483 0 : }
1484 :
1485 :
1486 : template <class TDomain>
1487 0 : void MarkForRefinement_AnisotropicElements(TDomain& dom, IRefiner& refiner,
1488 : number minEdgeRatio)
1489 : {
1490 : typedef typename domain_traits<TDomain::dim>::element_type TElem;
1491 :
1492 0 : Grid& grid = *dom.grid();
1493 0 : MarkForAnisotropicRefinement(grid, refiner, minEdgeRatio,
1494 : grid.begin<TElem>(), grid.end<TElem>(),
1495 : dom.position_accessor());
1496 0 : }
1497 :
1498 :
1499 0 : void MarkNeighborsForFullRefinement(IRefiner& refiner, bool sideNbrsOnly)
1500 : {
1501 0 : refiner.mark_neighborhood(1, RM_REFINE, sideNbrsOnly);
1502 0 : }
1503 :
1504 0 : void MarkNeighborsForAnisotropicRefinement(IRefiner& refiner, bool sideNbrsOnly)
1505 : {
1506 0 : refiner.mark_neighborhood(1, RM_ANISOTROPIC, sideNbrsOnly);
1507 0 : }
1508 :
1509 0 : void MarkNeighborsForLocalRefinement(IRefiner& refiner, bool sideNbrsOnly)
1510 : {
1511 0 : refiner.mark_neighborhood(1, RM_LOCAL, sideNbrsOnly);
1512 0 : }
1513 :
1514 :
1515 : template <class TDomain>
1516 0 : void MarkForRefinement_AnisotropicDirection_ (
1517 : TDomain& dom,
1518 : IRefiner& refiner,
1519 : const MathVector<TDomain::dim>& dir,
1520 : number minEdgeRatio,
1521 : RefinementMark elem_default_mark)
1522 :
1523 : {
1524 : using std::min;
1525 : using std::max;
1526 :
1527 : typedef MathVector<TDomain::dim> vector_t;
1528 : typedef typename domain_traits<TDomain::dim>::element_type TElem;
1529 :
1530 : // (normalized) user direction d
1531 : vector_t ndir;
1532 0 : VecNormalize(ndir, dir);
1533 :
1534 : typename TDomain::position_accessor_type aaPos = dom.position_accessor();
1535 0 : MultiGrid& mg = *dom.grid();
1536 :
1537 : MultiGrid::edge_traits::secure_container assEdges;
1538 :
1539 : vector<Edge*> anisoEdges;
1540 : vector<Edge*> normalEdges;
1541 :
1542 0 : lg_for_each_template(TElem, elem, mg){
1543 0 : if(mg.has_children(elem))
1544 0 : continue;
1545 :
1546 0 : number shortestAnisoEdgeSq = numeric_limits<number>::max();
1547 0 : number longestNormalEdgeSq = 0;
1548 :
1549 : // we'll mark all elements as closure since we use copy elements anyway
1550 0 : refiner.mark(elem, elem_default_mark);
1551 :
1552 : anisoEdges.clear();
1553 : normalEdges.clear();
1554 :
1555 0 : mg.associated_elements(assEdges, elem);
1556 0 : for_each_in_vec(Edge* e, assEdges){
1557 :
1558 : // normalized edge direction e
1559 : vector_t edgeDir;
1560 0 : VecSubtract(edgeDir, aaPos[e->vertex(0)], aaPos[e->vertex(1)]);
1561 0 : VecNormalize(edgeDir, edgeDir);
1562 :
1563 : // e*d = cos(alpha)
1564 : number dot = VecDot(edgeDir, ndir);
1565 0 : if((dot + SMALL >= 1.0) || (dot - SMALL <= -1.0)){
1566 : // measure shortest edge in direction
1567 0 : anisoEdges.push_back(e);
1568 0 : shortestAnisoEdgeSq = min(shortestAnisoEdgeSq, EdgeLengthSq(e, aaPos));
1569 : }
1570 : else{
1571 : // measure longest normal edge
1572 0 : normalEdges.push_back(e);
1573 0 : longestNormalEdgeSq = max(longestNormalEdgeSq, EdgeLengthSq(e, aaPos));
1574 : }
1575 : }end_for;
1576 :
1577 0 : if(longestNormalEdgeSq > 0){
1578 0 : if(shortestAnisoEdgeSq / longestNormalEdgeSq <= sq(minEdgeRatio)){
1579 : // mark all normal edges for full refinement
1580 0 : for_each_in_vec(Edge* e, normalEdges){
1581 0 : refiner.mark(e, RM_REFINE);
1582 : }end_for;
1583 : }
1584 : }
1585 : }lg_end_for;
1586 0 : }
1587 :
1588 :
1589 : template <class TDomain>
1590 0 : void MarkForRefinement_AnisotropicDirection (
1591 : TDomain& dom,
1592 : IRefiner& refiner,
1593 : const MathVector<TDomain::dim>& dir,
1594 : number minEdgeRatio)
1595 : {
1596 : MarkForRefinement_AnisotropicDirection_<TDomain>
1597 0 : (dom, refiner, dir, minEdgeRatio, RM_CLOSURE);
1598 0 : }
1599 :
1600 :
1601 : template <class TDomain>
1602 0 : void MarkForRefinement_AnisotropicDirection2 (
1603 : TDomain& dom,
1604 : IRefiner& refiner,
1605 : const MathVector<TDomain::dim>& dir,
1606 : number minEdgeRatio, std::string markType)
1607 : {
1608 : MarkForRefinement_AnisotropicDirection_<TDomain>
1609 0 : (dom, refiner, dir, minEdgeRatio, StringToRefinementMark(markType));
1610 0 : }
1611 :
1612 :
1613 :
1614 :
1615 : template <class TDomain>
1616 0 : void MarkForRefinement_EdgeDirection (
1617 : TDomain& dom,
1618 : IRefiner& refiner,
1619 : MathVector<TDomain::dim>& dir,
1620 : number minDeviationAngle,
1621 : number maxDeviationAngle,
1622 : bool selectFlipped)
1623 : {
1624 0 : MultiGrid& mg = *dom.grid();
1625 0 : MarkForRefinementByDirection( refiner,
1626 : dom.position_accessor(),
1627 : mg.begin<Edge>(),
1628 : mg.end<Edge>(),
1629 : dir,
1630 : minDeviationAngle,
1631 : maxDeviationAngle,
1632 : selectFlipped);
1633 0 : }
1634 :
1635 : // end group refinement_bridge
1636 : /// \}
1637 :
1638 : ////////////////////////////////////////////////////////////////////////////////
1639 : // REFINEMENT PROJECTORS
1640 : // template <class TDomain>
1641 : // SmartPtr<ProjectionHandler>
1642 : // DomainProjectionHandler (TDomain& dom)
1643 : // {
1644 : // return make_sp(new ProjectionHandler(dom.geometry3d(), dom.subset_handler()));
1645 : // }
1646 :
1647 : // template <class TDomain>
1648 : // SmartPtr<IRefinementCallback>
1649 : // LinearProjectorFactory(TDomain* dom)
1650 : // {
1651 : // return make_sp(new LinearProjector(
1652 : // make_sp(new Geometry<3, TDomain::dim>(
1653 : // *dom->grid(), dom->position_attachment()))));
1654 : // }
1655 :
1656 : // template <class vector_t>
1657 : // static
1658 : // vector_t StdVecToMathVec(const std::vector<number>& v)
1659 : // {
1660 : // vector_t mv;
1661 : // VecSet(mv, 0);
1662 : // for(size_t i = 0; (i < vector_t::Size) && (i < v.size()); ++i)
1663 : // mv[i] = v[i];
1664 : // return mv;
1665 : // }
1666 :
1667 : // /// Creates a refinement projector which projects new vertices onto a sphere
1668 : // /** Specify a domain, the center of the sphere (cx, cy, cz), and the sphere's radius.
1669 : // */
1670 : // template <class TDomain>
1671 : // SmartPtr<IRefinementCallback>
1672 : // SphereProjectorFactory(TDomain* dom, std::vector<number> center)
1673 : // {
1674 : // typedef SphereProjector<typename TDomain::position_attachment_type> TRefProj;
1675 : // return SmartPtr<TRefProj>(
1676 : // new TRefProj(*dom->grid(), dom->position_attachment(),
1677 : // StdVecToMathVec<typename TDomain::position_type>(center)));
1678 : // }
1679 :
1680 : // /// Creates a refinement projector which projects new vertices onto a sphere
1681 : // /** An outer radius can also be specified. Vertices outside this outer radius will
1682 : // * be projected linear.
1683 : // * Specify a domain, the center of the sphere (cx, cy, cz), an inner and an outer radius.
1684 : // */
1685 : // template <class TDomain>
1686 : // SmartPtr<IRefinementCallback>
1687 : // SphericalFalloffProjectorFactory(TDomain* dom, std::vector<number> center,
1688 : // number innerRadius, number outerRadius)
1689 : // {
1690 : // typedef SphericalFalloffProjector<typename TDomain::position_attachment_type> TRefProj;
1691 : // return SmartPtr<TRefProj>(
1692 : // new TRefProj(*dom->grid(), dom->position_attachment(),
1693 : // StdVecToMathVec<typename TDomain::position_type>(center),
1694 : // innerRadius, outerRadius));
1695 : // }
1696 :
1697 : // /// Creates a refinement projector which projects new vertices onto a cylinder
1698 : // /** Specify a domain, a point on the cylinder's axis c, the direction
1699 : // * of the axis
1700 : // */
1701 : // template <class TDomain>
1702 : // SmartPtr<IRefinementCallback>
1703 : // CylinderProjectorFactory(TDomain* dom, std::vector<number> c, std::vector<number> axis)
1704 : // {
1705 : // typedef CylinderProjector<typename TDomain::position_attachment_type> TRefProj;
1706 : // return SmartPtr<TRefProj>(
1707 : // new TRefProj(*dom->grid(), dom->position_attachment(),
1708 : // StdVecToMathVec<typename TDomain::position_type>(c),
1709 : // StdVecToMathVec<typename TDomain::position_type>(axis)));
1710 : // }
1711 :
1712 : // template <class TDomain>
1713 : // SmartPtr<IRefinementCallback>
1714 : // CylindricalFalloffProjectorFactory(TDomain* dom, std::vector<number> c,
1715 : // std::vector<number> a,
1716 : // number innerRadius, number outerRadius)
1717 : // {
1718 : // typedef CylindricalFalloffProjector<typename TDomain::position_attachment_type> TRefProj;
1719 : // return SmartPtr<TRefProj>(
1720 : // new TRefProj(*dom->grid(), dom->position_attachment(),
1721 : // StdVecToMathVec<typename TDomain::position_type>(c),
1722 : // StdVecToMathVec<typename TDomain::position_type>(a),
1723 : // innerRadius, outerRadius));
1724 : // }
1725 :
1726 : // template <class TDomain>
1727 : // SmartPtr<IRefinementCallback>
1728 : // SubdivisionLoopProjectorFactory(TDomain* dom)
1729 : // {
1730 : // typedef SubdivisionLoopProjector<typename TDomain::position_attachment_type> TRefProj;
1731 : // return SmartPtr<TRefProj>(
1732 : // new TRefProj(*dom->grid(), dom->position_attachment(),
1733 : // dom->position_attachment()));
1734 : // }
1735 :
1736 :
1737 : namespace domain_wrappers {
1738 :
1739 : }// end of namespace
1740 :
1741 :
1742 : /// Setting procedure for global refinement rule variable
1743 : /** The global refinement rule information switches between regular
1744 : * and subdivision volume refinement using hybrid tetra-/octahedral
1745 : * splitting.
1746 : */
1747 0 : void SetTetRefinementRule(std::string ruleName)
1748 : {
1749 0 : ruleName = ToLower(ruleName);
1750 0 : if(ruleName.compare("standard") == 0)
1751 0 : tet_rules::SetRefinementRule(tet_rules::STANDARD);
1752 0 : else if(ruleName.compare("hybrid_tet_oct") == 0)
1753 0 : tet_rules::SetRefinementRule(tet_rules::HYBRID_TET_OCT);
1754 : else{
1755 0 : UG_THROW("ERROR in SetTetRefinementRule:\n"
1756 : "Unknown refinement rule! Known rules are: standard, hybrid_tet_oct");
1757 : }
1758 0 : }
1759 :
1760 : /// Setting procedure for global boundary refinement rule variable
1761 : /** The global boundary refinement rule information switches between
1762 : * regular and a collection of subdivision surface refinement schemes.
1763 : */
1764 0 : void SetSmoothSubdivisionVolumesBoundaryRefinementRule(std::string bndRefRule)
1765 : {
1766 0 : bndRefRule = ToLower(bndRefRule);
1767 0 : if(bndRefRule.compare("linear") == 0)
1768 0 : SetBoundaryRefinementRule(LINEAR);
1769 0 : else if(bndRefRule.compare("subdiv_surf_loop_scheme") == 0)
1770 0 : SetBoundaryRefinementRule(SUBDIV_SURF_LOOP_SCHEME);
1771 0 : else if(bndRefRule.compare("subdiv_surf_averaging_scheme") == 0)
1772 0 : SetBoundaryRefinementRule(SUBDIV_SURF_AVERAGING_SCHEME);
1773 0 : else if(bndRefRule.compare("subdiv_surf_butterfly_scheme") == 0)
1774 0 : SetBoundaryRefinementRule(SUBDIV_SURF_BUTTERFLY_SCHEME);
1775 0 : else if(bndRefRule.compare("subdiv_vol") == 0)
1776 0 : SetBoundaryRefinementRule(SUBDIV_VOL);
1777 : else
1778 0 : UG_THROW("ERROR in SetBoundaryRefinementRule: Unknown boundary refinement rule!\n"
1779 : "Known rules are: 'linear', 'subdiv_surf_loop_scheme', 'subdiv_surf_averaging_scheme',"
1780 : " 'subdiv_surf_butterfly_scheme' or 'subdiv_vol'.");
1781 0 : }
1782 :
1783 : /// Helper function for subdivision volumes multigrid smoothing
1784 : /** Returns the default position attachment
1785 : * from the given domain needed for hierarchy smoothing
1786 : * using subdivision volumes schemes.
1787 : */
1788 : template <class TDomain>
1789 : typename TDomain::position_attachment_type&
1790 0 : GetDomainPositionAttachment(TDomain& dom)
1791 : {
1792 0 : return dom.position_attachment();
1793 : }
1794 :
1795 :
1796 : ////////////////////////////////////////////////////////////////////////////////
1797 : ////////////////////////////////////////////////////////////////////////////////
1798 : namespace bridge{
1799 : namespace Refinement{
1800 :
1801 : /// \addtogroup refinement_bridge
1802 : /// \{
1803 :
1804 : /**
1805 : * Class exporting the functionality. All functionality that is to
1806 : * be used in scripts or visualization must be registered here.
1807 : */
1808 : struct Functionality
1809 : {
1810 :
1811 : /**
1812 : * Function called for the registration of Domain and Algebra independent parts.
1813 : * All Functions and Classes not depending on Domain and Algebra
1814 : * are to be placed here when registering.
1815 : *
1816 : * @param reg registry
1817 : * @param parentGroup group for sorting of functionality
1818 : */
1819 1 : static void Common(Registry& reg, string grp)
1820 : {
1821 : // register domain independent mark methods
1822 3 : reg.add_function("MarkForRefinement_All",
1823 : &MarkForRefinement_All, grp, "", "ref");
1824 3 : reg.add_function("MarkForRefinement_AllVolumesAnisotropic",
1825 : &MarkForRefinement_AllAnisotropic<Volume>, grp, "", "ref");
1826 : // register refinement rule switch function
1827 3 : reg.add_function("SetTetRefinementRule",
1828 : &SetTetRefinementRule, grp, "", "refRuleName",
1829 : "Sets the refinement rule which is used to refine tetrahedrons. Possible parameters: 'standard', 'hybrid_tet_oct");
1830 : // register boundary refinement rule switch function for Subdivision Volumes smoothing
1831 3 : reg.add_function("SetSmoothSubdivisionVolumesBoundaryRefinementRule",
1832 : &SetSmoothSubdivisionVolumesBoundaryRefinementRule, grp, "", "bndRefRule",
1833 : "Sets the boundary refinement rule used during Subdivision Volumes smoothing. Possible parameters: 'linear', 'subdiv_surf_loop_scheme', 'subdiv_surf_averaging_scheme' or 'subdiv_vol'.");
1834 : // smooth volume/surface subdivision
1835 3 : reg.add_function("ApplySmoothSubdivisionVolumesToTopLevel",
1836 : (void (*)(ug::MultiGrid&, ug::MGSubsetHandler&, ug::MGSubsetHandler&, const char*)) (&ug::ApplySmoothSubdivisionVolumesToTopLevel), grp);
1837 3 : reg.add_function("ApplyConstrainedSmoothSubdivisionVolumesToTopLevel",
1838 : &ApplyConstrainedSmoothSubdivisionVolumesToTopLevel, grp);
1839 3 : reg.add_function("TetrahedralizeHybridTetOctGrid",
1840 : &TetrahedralizeHybridTetOctGrid, grp);
1841 3 : reg.add_function("CheckValences",
1842 : &CheckValences, grp);
1843 :
1844 3 : reg.add_function("MarkNeighborsForFullRefinement",
1845 : &MarkNeighborsForFullRefinement,
1846 : grp, "", "refiner#sideNeighborsOnly")
1847 3 : .add_function("MarkNeighborsForAnisotropicRefinement",
1848 : &MarkNeighborsForAnisotropicRefinement,
1849 : grp, "", "refiner#sideNeighborsOnly")
1850 3 : .add_function("MarkNeighborsForLocalRefinement",
1851 : &MarkNeighborsForLocalRefinement,
1852 : grp, "", "refiner#sideNeighborsOnly");
1853 :
1854 3 : reg.add_function("AddShadowCopyAdjuster", &AddShadowCopyAdjuster, grp, "", "refiner");
1855 :
1856 1 : }
1857 :
1858 : /**
1859 : * Function called for the registration of Domain dependent parts.
1860 : * All Functions and Classes depending on the Domain
1861 : * are to be placed here when registering. The method is called for all
1862 : * available Domain types, based on the current build options.
1863 : *
1864 : * @param reg registry
1865 : * @param parentGroup group for sorting of functionality
1866 : */
1867 : template <typename TDomain>
1868 3 : static void Domain(Registry& reg, string grp)
1869 : {
1870 : typedef TDomain domain_type;
1871 : typedef typename TDomain::position_attachment_type apos_type;
1872 :
1873 : string suffix = GetDomainSuffix<TDomain>();
1874 : string tag = GetDomainTag<TDomain>();
1875 :
1876 : // refiner factory-method registration
1877 : // Note that the refiners themselves have already been registered in lib_grid_bridge.
1878 9 : reg.add_function("GlobalDomainRefiner",
1879 : &GlobalDomainRefiner<domain_type>, grp, "GlobalDomainRefiner", "dom");
1880 9 : reg.add_function("HangingNodeDomainRefiner",
1881 : &HangingNodeDomainRefiner<domain_type>, grp, "HangingNodeDomainRefiner", "dom");
1882 9 : reg.add_function("GlobalFracturedDomainRefiner",
1883 : &CreateGlobalFracturedDomainRefiner<domain_type>, grp, "GlobalFracturedDomainRefiner", "dom");
1884 9 : reg.add_function("AdaptiveRegularDomainRefiner",
1885 : &CreateAdaptiveRegularDomainRefiner<domain_type>, grp, "AdaptiveRegularDomainRefiner", "dom");
1886 9 : reg.add_function("AddHorizontalAnisotropyAdjuster",
1887 : &AddHorizontalAnisotropyAdjuster<domain_type>, grp, "", "refiner # dom");
1888 9 : reg.add_function("GlobalSubdivisionDomainRefiner",
1889 : &GlobalSubdivisionDomainRefiner<domain_type>, grp, "GlobalSubdivisionDomainRefiner", "dom");
1890 9 : reg.add_function("ApplySmoothSubdivisionSurfacesToTopLevel",
1891 : (void (*)(ug::MultiGrid&, apos_type&, ug::MGSubsetHandler&, ug::MGSubsetHandler&, const char*, bool))
1892 : (&ApplySmoothSubdivisionSurfacesToTopLevel<apos_type>), grp);
1893 9 : reg.add_function("ProjectHierarchyToSubdivisionLimit",
1894 : &ProjectHierarchyToSubdivisionLimit<apos_type>, grp);
1895 9 : reg.add_function("GetDomainPositionAttachment",
1896 : &GetDomainPositionAttachment<domain_type>, grp);
1897 :
1898 : // register domain dependent mark methods
1899 9 : reg.add_function("MarkForRefinement_VerticesInSphere",
1900 : &MarkForRefinement_VerticesInSphere<domain_type>, grp,
1901 : "", "dom#refiner#center#radius")
1902 9 : .add_function("MarkForAdaption_VerticesInSphere",
1903 : &MarkForAdaption_VerticesInSphere<domain_type>, grp,
1904 : "", "dom#refiner#center#radius#adaption_type")
1905 9 : .add_function("MarkForAdaption_VerticesInSphereMaxLvl",
1906 : &MarkForAdaption_VerticesInSphereMaxLvl<domain_type>, grp,
1907 : "", "dom#refiner#center#radius#adaption_type")
1908 9 : .add_function("MarkForRefinement_EdgesInSphere",
1909 : &MarkForRefinement_ElementsInSphere<domain_type, Edge>, grp,
1910 : "", "dom#refiner#center#radius")
1911 9 : .add_function("MarkForRefinement_FacesInSphere",
1912 : &MarkForRefinement_ElementsInSphere<domain_type, Face>, grp,
1913 : "", "dom#refiner#center#radius")
1914 9 : .add_function("MarkForRefinement_VolumesInSphere",
1915 : &MarkForRefinement_ElementsInSphere<domain_type, Volume>, grp,
1916 : "", "dom#refiner#center#radius")
1917 9 : .add_function("MarkForRefinement_VerticesInCube",
1918 : &MarkForRefinement_VerticesInCube<domain_type>, grp,
1919 : "", "dom#refiner#min#max")
1920 9 : .add_function("MarkForAdaption_VerticesInCube",
1921 : &MarkForAdaption_VerticesInCube<domain_type>, grp,
1922 : "", "dom#refiner#min#max#adaption_type")
1923 9 : .add_function("MarkAnisotropic_LongEdges",
1924 : &MarkAnisotropic_LongEdges<domain_type>, grp,
1925 : "", "dom#refiner#maxEdgeLen")
1926 : // .add_function("MarkForRefinement_AnisotropicElements",
1927 : // &MarkForRefinement_AnisotropicElements<domain_type>, grp,
1928 : // "", "dom#refiner#sizeRatio")
1929 : // .add_function("MarkForRefinement_AnisotropicElements2",
1930 : // &MarkForRefinement_AnisotropicElements2<domain_type>, grp,
1931 : // "", "dom#refiner#sizeRatio")
1932 9 : .add_function("MarkForRefinement_ElementsByLuaCallback",
1933 : &MarkForRefinement_ElementsByLuaCallback<domain_type>, grp,
1934 : "", "dom#refiner#time#callbackName")
1935 9 : .add_function("MarkForCoarsen_ElementsByLuaCallback",
1936 : &MarkForCoarsen_ElementsByLuaCallback<domain_type>, grp,
1937 : "", "dom#refiner#time#callbackName")
1938 9 : .add_function("MarkForRefinement_ElementsInSubset",
1939 : &MarkForRefinement_ElementsInSubset<domain_type, MGSubsetHandler,
1940 : typename domain_traits<TDomain::dim>::element_type>,
1941 : grp, "", "dom#refiner#subsetHandler#subsetIndex")
1942 9 : .add_function("MarkForAdaption_ElementsInSubset",
1943 : &MarkForAdaption_ElementsInSubset<domain_type, MGSubsetHandler,
1944 : typename domain_traits<TDomain::dim>::element_type>,
1945 : grp, "", "dom#refiner#subsetHandler#subsetIndex#mark")
1946 9 : .add_function("MarkForRefinement_VerticesInSubset",
1947 : &MarkForRefinement_ElementsInSubset<domain_type, MGSubsetHandler, Vertex>,
1948 : grp, "", "dom#refiner#subsetHandler#subsetIndex")
1949 9 : .add_function("MarkForRefinement_EdgesInSubset",
1950 : &MarkForRefinement_ElementsInSubset<domain_type, MGSubsetHandler, Edge>,
1951 : grp, "", "dom#refiner#subsetHandler#subsetIndex")
1952 9 : .add_function("MarkForRefinement_FacesInSubset",
1953 : &MarkForRefinement_ElementsInSubset<domain_type, MGSubsetHandler, Face>,
1954 : grp, "", "dom#refiner#subsetHandler#subsetIndex")
1955 9 : .add_function("MarkForRefinement_VolumesInSubset",
1956 : &MarkForRefinement_ElementsInSubset<domain_type, MGSubsetHandler, Volume>,
1957 : grp, "", "dom#refiner#subsetHandler#subsetIndex")
1958 9 : .add_function("MarkForAdaption_ElementsTouchingSubset",
1959 : &MarkForAdaption_ElementsTouchingSubset<domain_type>,
1960 : grp, "", "dom#refiner#subsetHandler#subsetIndex#strMark")
1961 9 : .add_function("MarkForAdaption_ElementsTouchingSubsets",
1962 : &MarkForAdaption_ElementsTouchingSubsets<domain_type>,
1963 : grp, "", "dom#refiner#subsetHandler#subsetNames#strMark")
1964 9 : .add_function("MarkForRefinement_SubsetInterfaceVertices",
1965 : &MarkForRefinement_SubsetInterfaceElements<TDomain, Vertex>,
1966 : grp, "dom#refiner")
1967 9 : .add_function("MarkForRefinement_SubsetInterfaceEdges",
1968 : &MarkForRefinement_SubsetInterfaceElements<TDomain, Edge>,
1969 : grp, "dom#refiner")
1970 9 : .add_function("MarkForRefinement_SubsetInterfaceFaces",
1971 : &MarkForRefinement_SubsetInterfaceElements<TDomain, Face>,
1972 : grp, "dom#refiner")
1973 9 : .add_function("MarkForRefinement_AnisotropicElements",
1974 : &MarkForRefinement_AnisotropicElements<domain_type>,
1975 : grp, "", "dom#refiner#minEdgeRatio")
1976 9 : .add_function("MarkForRefinement_AnisotropicDirection",
1977 : &MarkForRefinement_AnisotropicDirection<domain_type>,
1978 : grp, "", "dom#refiner#dir#minEdgeRatio")
1979 9 : .add_function("MarkForRefinement_AnisotropicDirection2",
1980 : &MarkForRefinement_AnisotropicDirection2<domain_type>,
1981 : grp, "", "dom#refiner#dir#minEdgeRatio")
1982 9 : .add_function("MarkForRefinement_EdgeDirection",
1983 : &MarkForRefinement_EdgeDirection<domain_type>,
1984 : grp, "", "dom#refiner#dir#minDeviationAngle#maxDeviationAngle#selectFlipped")
1985 9 : .add_function("MarkForRefinement_CloseToSurface",
1986 : &MarkForRefinement_CloseToSurface<domain_type>,
1987 : grp, "", "")
1988 9 : .add_function("MarkForRefinement_ContainsSurfaceNode",
1989 : &MarkForRefinement_ContainsSurfaceNode<domain_type>,
1990 : grp, "", "")
1991 9 : .add_function("AssignSubset_VerticesInCube",
1992 : &AssignSubset_VerticesInCube<domain_type>,
1993 : grp, "", "")
1994 9 : .add_function("AssignSubset_VerticesInSphere",
1995 : &AssignSubset_VerticesInSphere<domain_type>,
1996 : grp, "", "");
1997 :
1998 :
1999 : // .add_function("MarkForAdaption_EdgesContainingPoint",
2000 : // &MarkForAdaption_ElementsContainingPoint<domain_type, Edge>,
2001 : // grp, "", "dom#refiner#x#y#z#markType")
2002 : // .add_function("MarkForAdaption_FacesContainingPoint",
2003 : // &MarkForAdaption_ElementsContainingPoint<domain_type, Face>,
2004 : // grp, "", "dom#refiner#x#y#z#markType");
2005 : // .add_function("MarkForRefinement_VolumesContainingPoint",
2006 : // &MarkForAdaption_ElementsContainingPoint<domain_type, Volume>,
2007 : // grp, "", "dom#refiner#x#y#z#markType");
2008 :
2009 :
2010 : // register refinement projection handler and factories
2011 : // {
2012 : // typedef RefinementProjectionHandler<apos_type> T;
2013 : // string name = string("RefinementProjectionHandler").append(suffix);
2014 : // reg.add_class_<T, IRefinementCallback>(name, grp)
2015 : // .add_method("set_default_callback", &T::set_default_callback, grp)
2016 : // .add_method("set_callback",
2017 : // static_cast<void (T::*)(int, SmartPtr<IRefinementCallback>) >
2018 : // (&T::set_callback), grp)
2019 : // .add_method("set_callback",
2020 : // static_cast<void (T::*)(std::string, SmartPtr<IRefinementCallback>) >
2021 : // (&T::set_callback), grp);
2022 : // reg.add_class_to_group(name, "RefinementProjectionHandler", tag);
2023 : // }
2024 :
2025 : // reg.add_function("DomainRefinementProjectionHandler",
2026 : // &DomainRefinementProjectionHandler<TDomain>, grp,
2027 : // "RefinementProjectionHandler", "domain")
2028 : // .add_function("LinearProjector", &LinearProjectorFactory<TDomain>, grp,
2029 : // "IRefinementCallback", "domain")
2030 : // .add_function("SphereProjector", &SphereProjectorFactory<TDomain>, grp,
2031 : // "IRefinementCallback", "domain#centerX#centerY#centerZ#radius")
2032 : // .add_function("SphericalFalloffProjector", &SphericalFalloffProjectorFactory<TDomain>, grp,
2033 : // "IRefinementCallback", "domain#centerX#centerY#centerZ#innerRadius#outerRadius")
2034 : // .add_function("CylinderProjector", &CylinderProjectorFactory<TDomain>, grp,
2035 : // "IRefinementCallback", "domain#centerX#centerY#centerZ#axisX#axisY#axisZ")
2036 : // .add_function("CylindricalFalloffProjector", &CylindricalFalloffProjectorFactory<TDomain>, grp,
2037 : // "IRefinementCallback", "domain#centerX#centerY#centerZ#axisX#axisY#axisZ#innerRadius#outerRadius")
2038 : // .add_function("SubdivisionLoopProjector", &SubdivisionLoopProjectorFactory<TDomain>, grp,
2039 : // "IRefinementCallback", "domain");
2040 3 : }
2041 :
2042 : }; // end Functionality
2043 :
2044 : // end group refinement_bridge
2045 : /// \}
2046 :
2047 : }// end Refinement
2048 :
2049 : /// \addtogroup refinement_bridge
2050 1 : void RegisterBridge_Refinement(Registry& reg, string grp)
2051 : {
2052 1 : grp.append("/Refinement");
2053 : typedef Refinement::Functionality Functionality;
2054 :
2055 : try{
2056 2 : RegisterCommon<Functionality>(reg,grp);
2057 1 : RegisterDomainDependent<Functionality>(reg,grp);
2058 : }
2059 0 : UG_REGISTRY_CATCH_THROW(grp);
2060 1 : }
2061 :
2062 : }// end of namespace bridge
2063 : }// end of namespace ug
|
