Line data Source code
1 : /*
2 : * Copyright (c) 2011-2018: G-CSC, Goethe University Frankfurt
3 : * Author: Markus Breit
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 "marking_utils.h"
34 :
35 : #include "lib_disc/domain.h"
36 : #include "lib_disc/domain_traits.h"
37 : #include "lib_grid/algorithms/geom_obj_util/anisotropy_util.h"
38 : #include "lib_grid/tools/grid_level.h"
39 : #include "lib_grid/tools/subset_group.h"
40 : #include "lib_grid/tools/surface_view.h"
41 :
42 :
43 : namespace ug {
44 :
45 :
46 :
47 : template <typename TDomain>
48 0 : void MarkGlobal(SmartPtr<IRefiner> refiner, SmartPtr<TDomain> domain)
49 : {
50 : typedef typename domain_traits<TDomain::dim>::element_type elem_type;
51 : typedef typename SurfaceView::traits<elem_type>::const_iterator const_iterator;
52 :
53 : // get surface view
54 0 : SurfaceView sv(domain->subset_handler());
55 :
56 : // loop elements for marking
57 0 : const_iterator iter = sv.begin<elem_type>(GridLevel(), SurfaceView::ALL_BUT_SHADOW_COPY);
58 0 : const_iterator iterEnd = sv.end<elem_type>(GridLevel(), SurfaceView::ALL_BUT_SHADOW_COPY);
59 0 : for (; iter != iterEnd; ++iter)
60 0 : refiner->mark(*iter, RM_FULL);
61 0 : }
62 :
63 :
64 : template <typename TDomain>
65 0 : void MarkSubsets
66 : (
67 : SmartPtr<IRefiner> refiner,
68 : SmartPtr<TDomain> domain,
69 : const std::vector<std::string>& vSubset
70 : )
71 : {
72 : typedef typename domain_traits<TDomain::dim>::element_type elem_type;
73 : typedef typename SurfaceView::traits<elem_type>::const_iterator const_iterator;
74 :
75 : // get subset handler
76 : SmartPtr<MGSubsetHandler> sh = domain->subset_handler();
77 :
78 : // transform subset names to indices
79 0 : std::vector<bool> contained(sh->num_subsets(), false);
80 : {
81 0 : SubsetGroup ssg(sh);
82 0 : try {ssg.add(vSubset);}
83 0 : UG_CATCH_THROW("MarkSubsets failed to add subsets.");
84 :
85 : const size_t nSs = ssg.size();
86 0 : for (size_t i = 0; i < nSs; ++i)
87 : contained[ssg[i]] = true;
88 : }
89 :
90 : // get surface view
91 0 : SurfaceView sv(sh);
92 :
93 : // loop elements for marking
94 0 : const_iterator iter = sv.begin<elem_type>(GridLevel(), SurfaceView::ALL_BUT_SHADOW_COPY);
95 0 : const_iterator iterEnd = sv.end<elem_type>(GridLevel(), SurfaceView::ALL_BUT_SHADOW_COPY);
96 0 : for (; iter != iterEnd; ++iter)
97 0 : if (contained[sh->get_subset_index(*iter)])
98 0 : refiner->mark(*iter, RM_FULL);
99 0 : }
100 :
101 :
102 : template <typename TDomain>
103 0 : void MarkAlongSurface
104 : (
105 : SmartPtr<IRefiner> refiner,
106 : SmartPtr<TDomain> domain,
107 : const std::vector<std::string>& surfaceSubsets,
108 : const std::vector<std::string>& volumeSubsets
109 : )
110 : {
111 : typedef typename domain_traits<TDomain::dim>::element_type elem_type;
112 : typedef typename MultiGrid::traits<elem_type>::secure_container elem_list_type;
113 : typedef typename SurfaceView::traits<Vertex>::const_iterator iter_type;
114 :
115 : const size_t nSurf = surfaceSubsets.size();
116 0 : UG_COND_THROW(volumeSubsets.size() != nSurf, "Same number of surface and volume subsets required.");
117 :
118 0 : Grid* grid = refiner->grid();
119 : SmartPtr<MGSubsetHandler> spSH = domain->subset_handler();
120 0 : const SurfaceView sv(spSH);
121 :
122 0 : for (size_t i = 0; i < nSurf; ++i)
123 : {
124 : int surf_si;
125 : int vol_si;
126 : try
127 : {
128 0 : surf_si = spSH->get_subset_index(surfaceSubsets[i].c_str());
129 0 : vol_si = spSH->get_subset_index(volumeSubsets[i].c_str());
130 : }
131 0 : UG_CATCH_THROW("Cannot convert subset names " << surfaceSubsets[i] << " and "
132 : << volumeSubsets[i] << " to indices.");
133 :
134 : // loop elements for marking
135 0 : iter_type iter = sv.begin<Vertex>(surf_si, GridLevel(), SurfaceView::MG_ALL);
136 0 : iter_type iterEnd = sv.end<Vertex>(surf_si, GridLevel(), SurfaceView::MG_ALL);
137 0 : for (; iter != iterEnd; ++iter)
138 : {
139 : Vertex* v = *iter;
140 : elem_list_type el;
141 : grid->associated_elements(el, v);
142 : size_t el_sz = el.size();
143 0 : for (size_t j = 0; j < el_sz; ++j)
144 : {
145 : int elem_si = spSH->get_subset_index(el[j]);
146 0 : if (elem_si == vol_si && !sv.surface_state(el[j]).contains(SurfaceView::MG_SHADOW_PURE))
147 : {
148 : // mark the element for anisotropic refinement
149 0 : refiner->mark(el[j], RM_REFINE);
150 : }
151 : }
152 : }
153 : }
154 0 : }
155 :
156 :
157 : template <typename TDomain>
158 0 : void MarkAnisotropic
159 : (
160 : SmartPtr<IRefiner> refiner,
161 : SmartPtr<TDomain> domain,
162 : number thresholdRatio
163 : )
164 : {
165 : typedef typename domain_traits<TDomain::dim>::element_type elem_type;
166 : typedef typename domain_traits<TDomain::dim>::side_type side_type;
167 : typedef typename SurfaceView::traits<elem_type>::const_iterator const_iterator;
168 :
169 0 : Grid& grid = *refiner->grid();
170 : typename TDomain::position_accessor_type aaPos = domain->position_accessor();
171 :
172 : // get surface view and prepare loop over all surface elements
173 0 : SurfaceView sv(domain->subset_handler());
174 0 : const_iterator iter = sv.begin<elem_type>(GridLevel(), SurfaceView::ALL_BUT_SHADOW_COPY);
175 0 : const_iterator iterEnd = sv.end<elem_type>(GridLevel(), SurfaceView::ALL_BUT_SHADOW_COPY);
176 :
177 : // loop elements for marking
178 : std::vector<Edge*> longEdges;
179 0 : for (; iter != iterEnd; ++iter)
180 : {
181 : longEdges.clear();
182 0 : AnisotropyState state = long_edges_of_anisotropic_elem(*iter, grid, aaPos, thresholdRatio, longEdges);
183 0 : if (state != ISOTROPIC)
184 : {
185 : // mark elem
186 0 : refiner->mark(*iter, RM_CLOSURE);
187 :
188 : // mark long edges
189 : const size_t nEdges = longEdges.size();
190 0 : for (size_t e = 0; e < nEdges; ++e)
191 0 : refiner->mark(longEdges[e], RM_FULL);
192 :
193 : // mark all sides
194 : typename Grid::traits<side_type>::secure_container sl;
195 : grid.associated_elements(sl, *iter);
196 : const size_t slSz = sl.size();
197 0 : for (size_t s = 0; s < slSz; ++s)
198 0 : if (refiner->get_mark(sl[s]) != RM_FULL)
199 0 : refiner->mark(sl[s], RM_CLOSURE);
200 : }
201 : }
202 0 : }
203 :
204 :
205 : template <typename TDomain>
206 0 : void MarkAnisotropicOnlyX
207 : (
208 : SmartPtr<IRefiner> refiner,
209 : SmartPtr<TDomain> domain,
210 : number thresholdRatio
211 : )
212 : {
213 : typedef typename domain_traits<TDomain::dim>::element_type elem_type;
214 : typedef typename domain_traits<TDomain::dim>::side_type side_type;
215 : typedef typename SurfaceView::traits<elem_type>::const_iterator const_iterator;
216 :
217 0 : Grid& grid = *refiner->grid();
218 : typename TDomain::position_accessor_type aaPos = domain->position_accessor();
219 :
220 : // get surface view and prepare loop over all surface elements
221 0 : SurfaceView sv(domain->subset_handler());
222 0 : const_iterator iter = sv.begin<elem_type>(GridLevel(), SurfaceView::ALL_BUT_SHADOW_COPY);
223 0 : const_iterator iterEnd = sv.end<elem_type>(GridLevel(), SurfaceView::ALL_BUT_SHADOW_COPY);
224 :
225 : // loop elements for marking
226 : std::vector<Edge*> longEdges;
227 0 : for (; iter != iterEnd; ++iter)
228 : {
229 : longEdges.clear();
230 0 : AnisotropyState state = long_edges_of_anisotropic_elem(*iter, grid, aaPos, thresholdRatio, longEdges);
231 0 : if (state == ISOTROPIC)
232 0 : continue;
233 :
234 0 : UG_COND_THROW(!longEdges.size(), "Element is anisotropic, but no long edges present.");
235 :
236 : // check whether edges point in x direction
237 0 : Edge* longEdge = longEdges[0];
238 : MathVector<TDomain::dim> dir;
239 0 : VecSubtract(dir, aaPos[longEdge->vertex(1)], aaPos[longEdge->vertex(0)]);
240 0 : VecNormalize(dir, dir);
241 0 : if (fabs(dir[0]) > 0.9)
242 : {
243 : // mark elem
244 0 : refiner->mark(*iter, RM_CLOSURE);
245 :
246 : // mark long edges
247 : const size_t nEdges = longEdges.size();
248 0 : for (size_t e = 0; e < nEdges; ++e)
249 0 : refiner->mark(longEdges[e], RM_FULL);
250 :
251 : // mark all sides
252 : typename Grid::traits<side_type>::secure_container sl;
253 : grid.associated_elements(sl, *iter);
254 : const size_t slSz = sl.size();
255 0 : for (size_t s = 0; s < slSz; ++s)
256 0 : if (refiner->get_mark(sl[s]) != RM_FULL)
257 0 : refiner->mark(sl[s], RM_CLOSURE);
258 : }
259 : }
260 0 : }
261 :
262 :
263 :
264 : // explicit template specializations
265 : #ifdef UG_DIM_1
266 : template void MarkGlobal<Domain1d>(SmartPtr<IRefiner>, SmartPtr<Domain1d>);
267 : template void MarkSubsets<Domain1d>(SmartPtr<IRefiner>, SmartPtr<Domain1d>, const std::vector<std::string>&);
268 : template void MarkAlongSurface(SmartPtr<IRefiner>, SmartPtr<Domain1d>, const std::vector<std::string>&,
269 : const std::vector<std::string>&);
270 : template void MarkAnisotropic<Domain1d>(SmartPtr<IRefiner>, SmartPtr<Domain1d>, number);
271 : template void MarkAnisotropicOnlyX<Domain1d>(SmartPtr<IRefiner>, SmartPtr<Domain1d>, number);
272 : #endif
273 : #ifdef UG_DIM_2
274 : template void MarkGlobal<Domain2d>(SmartPtr<IRefiner>, SmartPtr<Domain2d>);
275 : template void MarkSubsets<Domain2d>(SmartPtr<IRefiner>, SmartPtr<Domain2d>, const std::vector<std::string>&);
276 : template void MarkAlongSurface(SmartPtr<IRefiner>, SmartPtr<Domain2d>, const std::vector<std::string>&,
277 : const std::vector<std::string>&);
278 : template void MarkAnisotropic<Domain2d>(SmartPtr<IRefiner>, SmartPtr<Domain2d>, number);
279 : template void MarkAnisotropicOnlyX<Domain2d>(SmartPtr<IRefiner>, SmartPtr<Domain2d>, number);
280 : #endif
281 : #ifdef UG_DIM_3
282 : template void MarkGlobal<Domain3d>(SmartPtr<IRefiner>, SmartPtr<Domain3d>);
283 : template void MarkSubsets<Domain3d>(SmartPtr<IRefiner>, SmartPtr<Domain3d>, const std::vector<std::string>&);
284 : template void MarkAlongSurface(SmartPtr<IRefiner>, SmartPtr<Domain3d>, const std::vector<std::string>&,
285 : const std::vector<std::string>&);
286 : template void MarkAnisotropic<Domain3d>(SmartPtr<IRefiner>, SmartPtr<Domain3d>, number);
287 : template void MarkAnisotropicOnlyX<Domain3d>(SmartPtr<IRefiner>, SmartPtr<Domain3d>, number);
288 : #endif
289 :
290 :
291 : } // end namespace ug
|
