Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
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 "adaption_surface_grid_function.h"
34 : #include "lib_disc/dof_manager/orientation.h"
35 : #include "lib_disc/function_spaces/local_transfer_interface.h"
36 :
37 : namespace ug{
38 :
39 : template <typename TDomain>
40 0 : AdaptionSurfaceGridFunction<TDomain>::
41 : AdaptionSurfaceGridFunction(SmartPtr<TDomain> spDomain, bool bObserveStorage)
42 : : m_spDomain(spDomain),
43 : m_spGrid(spDomain->grid()),
44 0 : m_bObserveStorage(bObserveStorage),
45 0 : m_aValue(true)
46 0 : {}
47 :
48 : template <typename TDomain>
49 0 : void AdaptionSurfaceGridFunction<TDomain>::
50 : prolongate(const GridMessage_Adaption& msg)
51 : {
52 : // \todo: this temporary - handle more flexible
53 : m_vpProlong.clear();
54 0 : for(size_t fct = 0; fct < m_spDDInfo->num_fct(); ++fct)
55 0 : m_vpProlong.push_back(GetStandardElementProlongation<TDomain>(m_spDDInfo->lfeid(fct)));
56 :
57 : const GridObjectCollection& goc = msg.affected_elements();
58 0 : for(size_t lvl = 0; lvl < goc.num_levels(); ++lvl)
59 : {
60 0 : prolongate<Vertex>(msg, lvl);
61 0 : prolongate<Edge>(msg, lvl);
62 0 : prolongate<Face>(msg, lvl);
63 0 : prolongate<Volume>(msg, lvl);
64 : }
65 0 : }
66 :
67 : template <typename TDomain>
68 : template <typename TBaseElem>
69 0 : void AdaptionSurfaceGridFunction<TDomain>::
70 : prolongate(const GridMessage_Adaption& msg, const size_t lvl)
71 : {
72 : const GridBaseObjectId gbo = (GridBaseObjectId)TBaseElem::BASE_OBJECT_ID;
73 :
74 : // extract and init prolongations working on the base element type
75 : std::vector<SmartPtr<IElemProlongation<TDomain> > > vpProlong;
76 0 : for(size_t fct = 0; fct < m_vpProlong.size(); ++fct){
77 :
78 0 : if(!m_vpProlong[fct]->perform_prolongation_on(gbo)) continue;
79 :
80 0 : m_vpProlong[fct]->init(m_spDomain,
81 0 : make_sp(new ValueAccessor(*this, fct)),
82 0 : make_sp(new ValueAccessor(*this, fct)));
83 :
84 0 : vpProlong.push_back(m_vpProlong[fct]);
85 : }
86 :
87 : // check that something to do
88 0 : if(vpProlong.empty()) return;
89 :
90 : // iterators
91 : const GridObjectCollection& goc = msg.affected_elements();
92 : typedef typename GridObjectCollection::traits<TBaseElem>::const_iterator const_iterator;
93 :
94 :
95 : const_iterator iter = goc.begin<TBaseElem>(lvl);
96 : const_iterator iterEnd = goc.end<TBaseElem>(lvl);
97 :
98 : // loop base element type
99 0 : for( ; iter != iterEnd; ++iter)
100 : {
101 : // get parent element, that has been refined
102 : TBaseElem* parent = *iter;
103 :
104 : // call implementations
105 0 : for(size_t f = 0; f < vpProlong.size(); ++f){
106 0 : vpProlong[f]->prolongate(parent);
107 : }
108 : }
109 0 : }
110 :
111 : //template <typename TElem, typename TSubElem>
112 : //static void select_associated(MGSelector& sel, TElem* elem)
113 : //{
114 : // MultiGrid* mg = sel.multi_grid();
115 : //
116 : // typename Grid::traits<TSubElem>::secure_container vSubElem;
117 : // mg->associated_elements(vSubElem, elem);
118 : // for(size_t i = 0; i < vSubElem.size(); ++i)
119 : // sel.select(vSubElem[i]);
120 : //}
121 : //
122 : //template <typename TElem>
123 : //static void select_associated(MGSelector& sel, TElem* elem)
124 : //{
125 : // // add all subelements
126 : // switch(TElem::dim){
127 : // case 3: select_associated<TElem, Face>(sel, elem);
128 : // case 2: select_associated<TElem, Edge>(sel, elem);
129 : // case 1: select_associated<TElem, Vertex>(sel, elem);
130 : // default: ;
131 : // }
132 : //}
133 : //
134 : //static void select_associated(MGSelector& sel, GridObject* elem)
135 : //{
136 : // switch(elem->base_object_id()){
137 : // case VERTEX: return; // no sub elems
138 : // case EDGE: select_associated(sel, static_cast<Edge*>(elem)); return;
139 : // case FACE: select_associated(sel, static_cast<Face*>(elem)); return;
140 : // case VOLUME: select_associated(sel, static_cast<Volume*>(elem)); return;
141 : // default: UG_THROW("Dim not supported.");
142 : // }
143 : //}
144 :
145 :
146 : template <typename TDomain>
147 : template <typename TBaseElem>
148 0 : void AdaptionSurfaceGridFunction<TDomain>::
149 : select_parents(MGSelector& sel, const GridMessage_Adaption& msg)
150 : {
151 : // iterators
152 : const GridObjectCollection& goc = msg.affected_elements();
153 : typedef typename GridObjectCollection::traits<TBaseElem>::const_iterator const_iterator;
154 :
155 0 : for(size_t lvl = 0; lvl < goc.num_levels(); ++lvl)
156 : {
157 : const_iterator iter = goc.begin<TBaseElem>(lvl);
158 : const_iterator iterEnd = goc.end<TBaseElem>(lvl);
159 :
160 : // loop base element type
161 0 : for( ; iter != iterEnd; ++iter)
162 : {
163 : // get parent element, that is marked for coarsening
164 : TBaseElem* elem = *iter;
165 :
166 : // exclude if not in surface or no shadow
167 0 : if(m_aaValue[elem].size() != m_spDDInfo->num_fct())
168 0 : continue;
169 :
170 : // get parent
171 : GridObject* parent = sel.multi_grid()->get_parent(elem);
172 :
173 : // check whether a parent exist and if it is of the same type as the
174 : // elements considered. Parents of a different type are handled in
175 : // other calls to this method
176 0 : if(!parent || (parent->base_object_id() != TBaseElem::BASE_OBJECT_ID))
177 0 : continue;
178 :
179 : // select parent
180 0 : sel.select(parent);
181 : }
182 : }
183 0 : }
184 :
185 : template <typename TDomain>
186 0 : void AdaptionSurfaceGridFunction<TDomain>::
187 : do_restrict(const GridMessage_Adaption& msg)
188 : {
189 : // \todo: this temporary - handle more flexible
190 : m_vpRestrict.clear();
191 0 : for(size_t fct = 0; fct < m_spDDInfo->num_fct(); ++fct)
192 0 : m_vpRestrict.push_back(GetStandardElementRestriction<TDomain>(m_spDDInfo->lfeid(fct)));
193 :
194 0 : MGSelector sel(*m_spGrid);
195 :
196 0 : if(m_aaValue.is_valid_vertex_accessor())
197 0 : select_parents<Vertex>(sel, msg);
198 0 : if(m_aaValue.is_valid_edge_accessor())
199 0 : select_parents<Edge>(sel, msg);
200 0 : if(m_aaValue.is_valid_face_accessor())
201 0 : select_parents<Face>(sel, msg);
202 0 : if(m_aaValue.is_valid_volume_accessor())
203 0 : select_parents<Volume>(sel, msg);
204 :
205 0 : do_restrict<Vertex>(sel, msg);
206 0 : do_restrict<Edge>(sel, msg);
207 0 : do_restrict<Face>(sel, msg);
208 0 : do_restrict<Volume>(sel, msg);
209 0 : }
210 :
211 : template <typename TDomain>
212 : template <typename TBaseElem>
213 0 : void AdaptionSurfaceGridFunction<TDomain>::
214 : do_restrict(const MGSelector& sel, const GridMessage_Adaption& msg)
215 : {
216 : const GridBaseObjectId gbo = (GridBaseObjectId)TBaseElem::BASE_OBJECT_ID;
217 :
218 : // extract and init prolongations working on the base element type
219 : std::vector<SmartPtr<IElemRestriction<TDomain> > > vpRestrict;
220 0 : for(size_t fct = 0; fct < m_vpRestrict.size(); ++fct){
221 :
222 0 : if(!m_vpRestrict[fct]->perform_restriction_on(gbo)) continue;
223 :
224 0 : m_vpRestrict[fct]->init(m_spDomain,
225 0 : make_sp(new ValueAccessor(*this, fct)),
226 0 : make_sp(new ValueAccessor(*this, fct)));
227 :
228 0 : vpRestrict.push_back(m_vpRestrict[fct]);
229 : }
230 :
231 : // check that something to do
232 0 : if(vpRestrict.empty()) return;
233 :
234 : // check that correct attachments used
235 0 : if(m_spDDInfo->max_dofs(gbo) == 0) return;
236 :
237 : // iterators
238 : typedef typename Selector::traits<TBaseElem>::const_level_iterator const_iterator;
239 :
240 0 : for(int lvl = sel.num_levels() - 1; lvl >= 0; --lvl)
241 : {
242 : const_iterator iter = sel.begin<TBaseElem>(lvl);
243 : const_iterator iterEnd = sel.end<TBaseElem>(lvl);
244 :
245 : // loop base element type
246 0 : for( ; iter != iterEnd; ++iter)
247 : {
248 : // get parent element, that has been refined
249 : TBaseElem* parent = *iter;
250 :
251 : // add storage for parent
252 0 : obj_created(parent);
253 :
254 : // call implementations
255 0 : for(size_t f = 0; f < vpRestrict.size(); ++f){
256 0 : vpRestrict[f]->do_restrict(parent);
257 : }
258 : }
259 : }
260 0 : }
261 :
262 :
263 : template <typename TDomain>
264 0 : void AdaptionSurfaceGridFunction<TDomain>::
265 : attach_entries(ConstSmartPtr<DoFDistributionInfo> spDDInfo)
266 : {
267 : GFUNCADAPT_PROFILE_FUNC();
268 0 : if(m_spGrid.invalid()) UG_THROW("Grid missing")
269 :
270 0 : m_spDDInfo = spDDInfo;
271 :
272 : // get required elem types
273 0 : const bool vrt = (spDDInfo->max_dofs(VERTEX) > 0);
274 0 : const bool edge = (spDDInfo->max_dofs(EDGE) > 0);
275 0 : const bool face = (spDDInfo->max_dofs(FACE) > 0);
276 0 : const bool vol = (spDDInfo->max_dofs(VOLUME) > 0);
277 :
278 : // attach storage arrays
279 0 : if(vrt) m_spGrid->attach_to<Vertex>(m_aValue);
280 0 : if(edge) m_spGrid->attach_to<Edge>(m_aValue);
281 0 : if(face) m_spGrid->attach_to<Face>(m_aValue);
282 0 : if(vol) m_spGrid->attach_to<Volume>(m_aValue);
283 :
284 : // prepare access
285 0 : m_aaValue.access(*m_spGrid, m_aValue, vrt, edge, face, vol);
286 :
287 0 : if(m_bObserveStorage)
288 : {
289 : // get type of observer
290 : int type = OT_GRID_OBSERVER;
291 0 : if(vrt) type |= OT_VERTEX_OBSERVER;
292 0 : if(edge) type |= OT_EDGE_OBSERVER;
293 0 : if(face) type |= OT_FACE_OBSERVER;
294 0 : if(vol) type |= OT_VOLUME_OBSERVER;
295 :
296 : // register observer
297 0 : m_spGrid->register_observer(this, type);
298 : }
299 0 : }
300 :
301 : template <typename TDomain>
302 : template <typename TElem>
303 0 : void AdaptionSurfaceGridFunction<TDomain>::detach_entries()
304 : {
305 0 : if(m_spGrid.invalid()) UG_THROW("Grid missing")
306 :
307 0 : if(m_spGrid->has_attachment<TElem>(m_aValue))
308 0 : m_spGrid->detach_from<TElem>(m_aValue);
309 0 : }
310 :
311 : template <typename TDomain>
312 0 : void AdaptionSurfaceGridFunction<TDomain>::detach_entries()
313 : {
314 : GFUNCADAPT_PROFILE_FUNC();
315 : // remove attachments
316 0 : if(m_aaValue.is_valid_vertex_accessor()) detach_entries<Vertex>();
317 0 : if(m_aaValue.is_valid_edge_accessor()) detach_entries<Edge>();
318 0 : if(m_aaValue.is_valid_face_accessor()) detach_entries<Face>();
319 0 : if(m_aaValue.is_valid_volume_accessor()) detach_entries<Volume>();
320 0 : m_aaValue = MultiElementAttachmentAccessor<AValues>();
321 :
322 : // stop observing grid
323 0 : if(m_bObserveStorage)
324 0 : m_spGrid->unregister_observer(this);
325 0 : }
326 :
327 : ////////////////////////////////////////////////////////////////////////////////
328 : // Value accessor
329 : ////////////////////////////////////////////////////////////////////////////////
330 :
331 : template <typename TDomain>
332 0 : AdaptionSurfaceGridFunction<TDomain>::ValueAccessor::
333 : ValueAccessor(AdaptionSurfaceGridFunction<TDomain>& rASGF,
334 : size_t fct)
335 0 : : m_rASGF(rASGF), m_fct(fct)
336 : {
337 0 : for(int g = 0; g < NUM_GEOMETRIC_BASE_OBJECTS; ++g)
338 : {
339 : const GridBaseObjectId gbo = (GridBaseObjectId)g;
340 0 : m_HasDoFs[gbo] = (m_rASGF.m_spDDInfo->max_fct_dofs(fct, gbo) > 0);
341 : }
342 0 : }
343 :
344 : template <typename TDomain>
345 : void
346 0 : AdaptionSurfaceGridFunction<TDomain>::ValueAccessor::
347 : access_inner(GridObject* elem)
348 : {
349 : UG_ASSERT(m_fct < m_rASGF.m_aaValue[elem].size(), "Only storage for "
350 : <<m_rASGF.m_aaValue[elem].size()<<" fcts, but fct-cmp "
351 : <<m_fct<<" requested on "<< ElementDebugInfo(*m_rASGF.m_spDomain->grid(), elem))
352 0 : std::vector<number>& vVal = m_rASGF.m_aaValue[elem][m_fct];
353 0 : m_Val.resize(vVal.size());
354 0 : for(size_t i = 0; i < vVal.size(); ++i)
355 0 : m_Val[i] = &vVal[i];
356 0 : }
357 :
358 : template <typename TDomain>
359 : template <typename TBaseElem, typename TSubBaseElem>
360 : void
361 0 : AdaptionSurfaceGridFunction<TDomain>::ValueAccessor::
362 : access_closure(TBaseElem* elem)
363 : {
364 :
365 : typename Grid::traits<TSubBaseElem>::secure_container vSubElem;
366 0 : m_rASGF.m_spGrid->associated_elements_sorted(vSubElem, elem);
367 :
368 : std::vector<size_t> vOrientOffset;
369 :
370 0 : for(size_t i = 0; i < vSubElem.size(); ++i)
371 : {
372 : // get subelement
373 : TSubBaseElem* subElem = vSubElem[i];
374 : UG_ASSERT(m_fct < m_rASGF.m_aaValue[subElem].size(), "Only storage for "
375 : <<m_rASGF.m_aaValue[subElem].size()<<" fcts, but fct-cmp "
376 : <<m_fct<<" requested on "<<ElementDebugInfo(*m_rASGF.m_spDomain->grid(), subElem))
377 0 : std::vector<number>& vVal = m_rASGF.m_aaValue[subElem][m_fct];
378 :
379 : // get the orientation for this subelement
380 0 : ComputeOrientationOffset(vOrientOffset, elem, subElem, i,
381 : m_rASGF.m_spDDInfo->lfeid(m_fct));
382 :
383 : UG_ASSERT(vOrientOffset.size() == vVal.size() ||
384 : vOrientOffset.empty(), "Orientation wrong");
385 :
386 : // cache access
387 0 : if(vOrientOffset.empty()){
388 0 : for(size_t j = 0; j < vVal.size(); ++j)
389 0 : m_Val.push_back(&vVal[ j ]);
390 : }else{
391 0 : for(size_t j = 0; j < vVal.size(); ++j)
392 0 : m_Val.push_back(&vVal[ vOrientOffset[j] ]);
393 : }
394 : }
395 0 : }
396 :
397 : template <typename TDomain>
398 : template <typename TBaseElem>
399 : void
400 0 : AdaptionSurfaceGridFunction<TDomain>::ValueAccessor::
401 : access_closure(TBaseElem* elem)
402 : {
403 : m_Val.clear();
404 0 : if(m_HasDoFs[VERTEX]) access_closure<TBaseElem, Vertex>(elem);
405 0 : if(m_HasDoFs[EDGE]) access_closure<TBaseElem, Edge>(elem);
406 0 : if(m_HasDoFs[FACE]) access_closure<TBaseElem, Face>(elem);
407 0 : if(m_HasDoFs[VOLUME]) access_closure<TBaseElem, Volume>(elem);
408 0 : }
409 :
410 : template <typename TDomain>
411 : void
412 0 : AdaptionSurfaceGridFunction<TDomain>::ValueAccessor::
413 : access_closure(GridObject* elem)
414 : {
415 0 : switch(elem->base_object_id()){
416 0 : case VERTEX: access_closure(static_cast<Vertex*>(elem)); return;
417 0 : case EDGE: access_closure(static_cast<Edge*>(elem)); return;
418 0 : case FACE: access_closure(static_cast<Face*>(elem)); return;
419 0 : case VOLUME: access_closure(static_cast<Volume*>(elem)); return;
420 0 : default: UG_THROW("Base object id not found.")
421 : }
422 : }
423 :
424 :
425 : #ifdef UG_DIM_1
426 : template class AdaptionSurfaceGridFunction<Domain1d>;
427 : #endif
428 : #ifdef UG_DIM_2
429 : template class AdaptionSurfaceGridFunction<Domain2d>;
430 : #endif
431 : #ifdef UG_DIM_3
432 : template class AdaptionSurfaceGridFunction<Domain3d>;
433 : #endif
434 :
435 : } // end namespace ug
|
