Line data Source code
1 : /*
2 : * Copyright (c) 2012-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 : #ifndef __H__UG__LIB_DISC__DOMAIN_IMPL__
34 : #define __H__UG__LIB_DISC__DOMAIN_IMPL__
35 :
36 : #include "domain.h"
37 : #include "common/serialization.h"
38 : #include "common/profiler/profiler.h"
39 :
40 : #ifdef UG_PARALLEL
41 : #include "lib_grid/refinement/projectors/projectors.h"
42 : #include "common/boost_serialization_routines.h"
43 : #include "common/util/archivar.h"
44 : #include "common/util/factory.h"
45 : #include <boost/archive/text_oarchive.hpp>
46 : #include <boost/archive/text_iarchive.hpp>
47 : #endif
48 :
49 : #ifdef UG_PARALLEL
50 : #include "pcl/pcl_process_communicator.h"
51 : #endif
52 :
53 :
54 : namespace ug{
55 :
56 : ////////////////////////////////////////////////////////////////////////////////
57 : // IDomain
58 : ////////////////////////////////////////////////////////////////////////////////
59 :
60 : template <typename TGrid, typename TSubsetHandler>
61 0 : IDomain<TGrid,TSubsetHandler>::IDomain(bool isAdaptive)
62 : :
63 0 : m_spGrid(new TGrid(GRIDOPT_NONE)), // Note: actual options are set by the derived class (dimension dependent).
64 0 : m_spSH(new TSubsetHandler(*m_spGrid)),
65 0 : m_isAdaptive(isAdaptive),
66 0 : m_adaptionIsActive(false)
67 : {
68 : #ifdef UG_PARALLEL
69 : // the grid has to be prepared for parallelism
70 : m_spGrid->set_parallel(true);
71 : #endif
72 :
73 : // register function for grid adaption
74 0 : m_spGridAdaptionCallbackID =
75 : message_hub()->register_class_callback(this,
76 : &ug::IDomain<ug::MultiGrid, ug::MultiGridSubsetHandler>::grid_adaption_callback);
77 :
78 0 : m_spGridCreationCallbackID =
79 : message_hub()->register_class_callback(this,
80 : &ug::IDomain<ug::MultiGrid, ug::MultiGridSubsetHandler>::grid_creation_callback);
81 :
82 0 : m_spGridDistributionCallbackID =
83 : message_hub()->register_class_callback(this,
84 : &ug::IDomain<ug::MultiGrid, ug::MultiGridSubsetHandler>::grid_distribution_callback);
85 0 : }
86 :
87 : /// Destructor
88 : template <typename TGrid, typename TSubsetHandler>
89 0 : IDomain<TGrid,TSubsetHandler>::~IDomain()
90 : {
91 0 : }
92 :
93 :
94 : template <typename TGrid, typename TSubsetHandler>
95 0 : void IDomain<TGrid,TSubsetHandler>::
96 : update_subset_infos(int rootProc)
97 : {
98 : PROFILE_FUNC();
99 :
100 : TSubsetHandler& sh = *m_spSH;
101 0 : for(int i = 0; i < sh.num_subsets(); ++i){
102 : int dim = -1;
103 0 : if(sh.contains_volumes(i))
104 : dim = 3;
105 0 : else if(sh.contains_faces(i))
106 : dim = 2;
107 0 : else if(sh.contains_edges(i))
108 : dim = 1;
109 0 : else if(sh.contains_vertices(i))
110 : dim = 0;
111 :
112 0 : sh.subset_info(i).set_property("dim", dim);
113 : }
114 :
115 : // do not communicate if geom is created on all procs
116 : if (rootProc == -2) return;
117 :
118 : #ifdef UG_PARALLEL
119 : pcl::ProcessCommunicator procCom;
120 :
121 : // prepare the subset-info package, send it to all procs and extract the info again.
122 : BinaryBuffer buf;
123 : if(pcl::ProcRank() == rootProc){
124 : Serialize(buf, sh.num_subsets());
125 : for(int i = 0; i < sh.num_subsets(); ++i){
126 : Serialize(buf, sh.subset_info(i).name);
127 : Serialize(buf, sh.subset_info(i).get_property("dim").to_int());
128 : }
129 : }
130 :
131 : procCom.broadcast(buf, rootProc);
132 :
133 : if(pcl::ProcRank() != rootProc){
134 : int numSubsets;
135 : Deserialize(buf, numSubsets);
136 : if(numSubsets > 0)
137 : sh.subset_required(numSubsets - 1);
138 :
139 : for(int i = 0; i < numSubsets; ++i){
140 : Deserialize(buf, sh.subset_info(i).name);
141 : int dim;
142 : Deserialize(buf, dim);
143 : sh.subset_info(i).set_property("dim", dim);
144 : }
145 : }
146 :
147 : //todo: distribute projectors from rootProc to all other processors.
148 : //note: first check whether source-proc has a ProjectionHandler. If so,
149 : // create a local projection handler first and perform serialization
150 : // afterwards.
151 : set_refinement_projector(
152 : broadcast_refinement_projector(
153 : rootProc, procCom, geometry3d(), m_refinementProjector));
154 :
155 : #endif
156 :
157 : }
158 :
159 :
160 : template <typename TGrid, typename TSubsetHandler>
161 : inline
162 0 : void IDomain<TGrid,TSubsetHandler>::
163 : grid_adaption_callback(const GridMessage_Adaption& msg)
164 : {
165 0 : if(msg.adaption_begins())
166 0 : m_adaptionIsActive = true;
167 :
168 0 : else if(m_adaptionIsActive){
169 : //if(msg.adaptive()){
170 0 : if(msg.adaption_ends())
171 : {
172 0 : update_domain_info();
173 0 : m_adaptionIsActive = false;
174 : }
175 : //}
176 : }
177 :
178 : else{
179 0 : UG_THROW("Before any grid-adaption may be performed, the domain "
180 : "has to be informed that grid-adaption shall begin. "
181 : "You may use IRefiner::grid_adaption_begins() or schedule "
182 : "an appropriate message to the associated grids message-hub.");
183 : }
184 0 : }
185 :
186 : template <typename TGrid, typename TSubsetHandler>
187 : inline
188 0 : void IDomain<TGrid,TSubsetHandler>::
189 : grid_creation_callback(const GridMessage_Creation& msg)
190 : {
191 0 : if(msg.msg() == GMCT_CREATION_STOPS){
192 0 : if(msg.proc_id() != -1)
193 0 : update_subset_infos(msg.proc_id());
194 0 : update_domain_info();
195 : }
196 0 : }
197 :
198 : template <typename TGrid, typename TSubsetHandler>
199 : inline
200 0 : void IDomain<TGrid,TSubsetHandler>::
201 : grid_distribution_callback(const GridMessage_Distribution& msg)
202 : {
203 : // this is already handled in grid_creation_callback
204 : /*if(msg.msg() == GMDT_DISTRIBUTION_STOPS){
205 : update_domain_info();
206 : }*/
207 0 : }
208 :
209 :
210 : template <typename TGrid, typename TSubsetHandler>
211 0 : void IDomain<TGrid,TSubsetHandler>::
212 : update_domain_info()
213 : {
214 : PROFILE_FUNC();
215 :
216 : TGrid& mg = *m_spGrid;
217 : TSubsetHandler& sh = *m_spSH;
218 :
219 : GridBaseObjectId locElemType;
220 0 : if(mg.template num<Volume>() > 0)
221 : locElemType = VOLUME;
222 0 : else if(mg.template num<Face>() > 0)
223 : locElemType = FACE;
224 0 : else if(mg.template num<Edge>() > 0)
225 : locElemType = EDGE;
226 : else
227 : locElemType = VERTEX;
228 :
229 : GridBaseObjectId elemType;
230 : #ifdef UG_PARALLEL
231 : pcl::ProcessCommunicator commWorld;
232 :
233 : // all processes should have the same number of grid levels
234 : int numLevLocal = m_spGrid->num_levels();
235 : int numLevGlobal;
236 : commWorld.allreduce(&numLevLocal, &numLevGlobal, 1, PCL_DT_INT, PCL_RO_MAX);
237 :
238 : if(numLevGlobal > 0){
239 : m_spGrid->level_required(numLevGlobal-1);
240 : m_spSH->level_required(numLevGlobal-1);
241 : }
242 :
243 : // communicate the element type of highest dimension present in the global grid.
244 : elemType = (GridBaseObjectId) commWorld.allreduce((int)locElemType, PCL_RO_MAX);
245 : #else
246 : elemType = locElemType;
247 : #endif
248 :
249 : // the number of levels of the multi-grid is now equal on all processes
250 : typedef DomainInfo::int_t int_t;
251 :
252 : std::vector<int_t> subsetDims;
253 0 : subsetDims.reserve(sh.num_subsets());
254 0 : for(int i = 0; i < sh.num_subsets(); ++i)
255 0 : subsetDims.push_back(sh.subset_info(i).get_property("dim").to_int());
256 :
257 : std::vector<int_t> numLocalGhosts;
258 : #ifdef UG_PARALLEL
259 : switch(elemType){
260 : case VOLUME:
261 : count_ghosts<Volume>(numLocalGhosts);
262 : break;
263 : case FACE:
264 : count_ghosts<Face>(numLocalGhosts);
265 : break;
266 : case EDGE:
267 : count_ghosts<Edge>(numLocalGhosts);
268 : break;
269 : case VERTEX:
270 : count_ghosts<Vertex>(numLocalGhosts);
271 : break;
272 : default:
273 : UG_THROW("Unknown base object type");
274 : break;
275 : }
276 : #else
277 0 : numLocalGhosts.resize(mg.num_levels(), 0);
278 : #endif
279 :
280 : int_t numLocalSurfElems = 0;
281 : std::vector<int_t> numLocalElems;
282 0 : numLocalElems.reserve(mg.num_levels());
283 :
284 0 : switch(elemType){
285 : case VOLUME:
286 0 : for(size_t i = 0; i < mg.num_levels(); ++i)
287 0 : numLocalElems.push_back(mg.template num<Volume>(i) - numLocalGhosts[i]);
288 0 : numLocalSurfElems = count_local_unique_surface_elements<Volume>();
289 0 : break;
290 : case FACE:
291 0 : for(size_t i = 0; i < mg.num_levels(); ++i)
292 0 : numLocalElems.push_back(mg.template num<Face>(i) - numLocalGhosts[i]);
293 0 : numLocalSurfElems = count_local_unique_surface_elements<Face>();
294 0 : break;
295 : case EDGE:
296 0 : for(size_t i = 0; i < mg.num_levels(); ++i)
297 0 : numLocalElems.push_back(mg.template num<Edge>(i) - numLocalGhosts[i]);
298 0 : numLocalSurfElems = count_local_unique_surface_elements<Edge>();
299 0 : break;
300 : case VERTEX:
301 0 : for(size_t i = 0; i < mg.num_levels(); ++i)
302 0 : numLocalElems.push_back(mg.template num<Vertex>(i) - numLocalGhosts[i]);
303 : numLocalSurfElems = count_local_unique_surface_elements<Vertex>();
304 : break;
305 : default:
306 : UG_THROW("Unknown base object type");
307 : break;
308 : }
309 :
310 :
311 :
312 : std::vector<int_t> numGlobalElems, minNumLocalElems, maxNumLocalElems;
313 : int_t numGlobalSurfElems = 0;
314 : #ifdef UG_PARALLEL
315 : // we have to sum local element counts excluding ghosts.
316 : numGlobalElems.resize(numLocalElems.size());
317 : minNumLocalElems.resize(numLocalElems.size());
318 : maxNumLocalElems.resize(numLocalElems.size());
319 : commWorld.allreduce(numLocalElems, numGlobalElems, PCL_RO_SUM);
320 : commWorld.allreduce(numLocalElems, minNumLocalElems, PCL_RO_MIN);
321 : commWorld.allreduce(numLocalElems, maxNumLocalElems, PCL_RO_MAX);
322 : numGlobalSurfElems = commWorld.allreduce(numLocalSurfElems, PCL_RO_SUM);
323 : #else
324 0 : numGlobalElems = numLocalElems;
325 0 : minNumLocalElems = numLocalElems;
326 0 : maxNumLocalElems = numLocalElems;
327 : numGlobalSurfElems = numLocalSurfElems;
328 : #endif
329 :
330 0 : m_domainInfo.set_info(elemType, numGlobalElems, numLocalElems, minNumLocalElems,
331 : maxNumLocalElems, numLocalGhosts, subsetDims, numGlobalSurfElems);
332 0 : }
333 :
334 : template <typename TGrid, typename TSubsetHandler>
335 0 : bool IDomain<TGrid,TSubsetHandler>::
336 : create_additional_subset_handler(std::string name)
337 : {
338 0 : if(m_additionalSH[name].valid())
339 : return false;
340 :
341 0 : m_additionalSH[name] = SmartPtr<TSubsetHandler>(new TSubsetHandler(*m_spGrid));
342 0 : return true;
343 : }
344 :
345 : template <typename TGrid, typename TSubsetHandler>
346 0 : std::vector<std::string> IDomain<TGrid,TSubsetHandler>::
347 : additional_subset_handler_names() const
348 : {
349 : typedef typename std::map<std::string, SmartPtr<TSubsetHandler> >::const_iterator iterator_t;
350 : std::vector<std::string> names;
351 0 : for(iterator_t iter = m_additionalSH.begin(); iter != m_additionalSH.end(); ++iter){
352 0 : if(iter->second.valid())
353 0 : names.push_back(iter->first);
354 : }
355 0 : return names;
356 0 : }
357 :
358 : template <typename TGrid, typename TSubsetHandler>
359 0 : SmartPtr<TSubsetHandler> IDomain<TGrid,TSubsetHandler>::
360 : additional_subset_handler(std::string name)
361 : {
362 0 : SmartPtr<TSubsetHandler> sp = m_additionalSH[name];
363 0 : if(!sp.valid()){
364 0 : UG_THROW("Requested additional subset handler with name '" << name
365 : << "' doesn't exist in the given domain!");
366 : }
367 0 : return sp;
368 : }
369 :
370 : template <typename TGrid, typename TSubsetHandler>
371 : const ConstSmartPtr<TSubsetHandler> IDomain<TGrid,TSubsetHandler>::
372 : additional_subset_handler(std::string name) const
373 : {
374 : SmartPtr<TSubsetHandler> sp = m_additionalSH[name];
375 : if(!sp.valid()){
376 : UG_THROW("Requested additional subset handler with name '" << name
377 : << "' doesn't exist in the given domain!");
378 : }
379 : return sp;
380 : }
381 :
382 : template <typename TGrid, typename TSubsetHandler>
383 0 : void IDomain<TGrid, TSubsetHandler>::
384 : set_refinement_projector(SPRefinementProjector proj)
385 : {
386 0 : m_refinementProjector = proj;
387 0 : if(proj.valid())
388 0 : proj->set_geometry(geometry3d());
389 0 : }
390 :
391 : template <typename TGrid, typename TSubsetHandler>
392 0 : SPRefinementProjector IDomain<TGrid, TSubsetHandler>::
393 : refinement_projector() const
394 : {
395 0 : return m_refinementProjector;
396 : }
397 :
398 :
399 : #ifdef UG_PARALLEL
400 :
401 : template <typename TGrid, typename TSubsetHandler>
402 : void IDomain<TGrid, TSubsetHandler>::
403 : serialize_projector (
404 : BinaryBuffer& buf,
405 : SPRefinementProjector proj)
406 : {
407 : static Archivar<boost::archive::text_oarchive,
408 : RefinementProjector,
409 : ProjectorTypes>
410 : archivar;
411 :
412 : static Factory<RefinementProjector, ProjectorTypes> projFac;
413 :
414 : if(proj.invalid()){
415 : Serialize(buf, std::string(""));
416 : return;
417 : }
418 :
419 : Serialize(buf, projFac.class_name(*proj));
420 :
421 : std::stringstream ss;
422 : boost::archive::text_oarchive ar(ss, boost::archive::no_header);
423 : archivar.archive(ar, *proj);
424 : Serialize(buf, ss.str());
425 : }
426 :
427 : template <typename TGrid, typename TSubsetHandler>
428 : SPRefinementProjector IDomain<TGrid, TSubsetHandler>::
429 : deserialize_projector (BinaryBuffer& buf)
430 : {
431 : static Archivar<boost::archive::text_iarchive,
432 : RefinementProjector,
433 : ProjectorTypes>
434 : archivar;
435 :
436 : static Factory<RefinementProjector, ProjectorTypes> projFac;
437 :
438 : std::string name;
439 : Deserialize(buf, name);
440 :
441 : if(name.empty())
442 : return SPRefinementProjector();
443 :
444 : SPRefinementProjector proj = projFac.create(name);
445 :
446 : std::string data;
447 : Deserialize(buf, data);
448 : std::stringstream ss(data, std::ios_base::in);
449 : boost::archive::text_iarchive ar(ss, boost::archive::no_header);
450 : archivar.archive(ar, *proj);
451 : return proj;
452 : }
453 :
454 :
455 : template <typename TGrid, typename TSubsetHandler>
456 : SPRefinementProjector IDomain<TGrid, TSubsetHandler>::
457 : broadcast_refinement_projector (
458 : int rootProc,
459 : pcl::ProcessCommunicator& procCom,
460 : SPIGeometry3d geometry,
461 : SPRefinementProjector projector)
462 : {
463 : BinaryBuffer buf;
464 : const int magicNumber = 3243578;
465 : const bool isRoot = (pcl::ProcRank() == rootProc);
466 :
467 : if(isRoot){
468 : // if the specified projector is a projection handler, we'll perform a
469 : // special operation.
470 : ProjectionHandler* ph = NULL;
471 : int projectorType = -1;// -1: none, 0: normal projector, 1: projection handler
472 : if(projector.valid()){
473 : ph = dynamic_cast<ProjectionHandler*>(projector.get());
474 : if(ph)
475 : projectorType = 1;
476 : else
477 : projectorType = 0;
478 : }
479 : Serialize(buf, projectorType);
480 : if(ph){
481 : const ISubsetHandler* psh = ph->subset_handler();
482 : if (psh == subset_handler().get())
483 : Serialize(buf, std::string(""));
484 : else
485 : {
486 : typedef typename std::map<std::string, SmartPtr<TSubsetHandler> >::const_iterator map_it_t;
487 : map_it_t it = m_additionalSH.begin();
488 : map_it_t it_end = m_additionalSH.end();
489 : for (; it != it_end; ++it)
490 : {
491 : if (it->second.get() == psh)
492 : {
493 : Serialize(buf, it->first);
494 : break;
495 : }
496 : }
497 : UG_COND_THROW(it == it_end, "Subset handler for projection handler not found "
498 : "in list of available subset handlers.");
499 : }
500 :
501 :
502 : serialize_projector(buf, ph->default_projector());
503 :
504 : size_t numProjectors = ph->num_projectors();
505 :
506 : Serialize(buf, numProjectors);
507 : for(size_t iproj = 0; iproj < numProjectors; ++iproj){
508 : SPRefinementProjector proj = ph->projector(iproj);
509 : serialize_projector(buf, proj);
510 : }
511 : }
512 : else if(projector.valid()){
513 : serialize_projector(buf, projector);
514 : }
515 :
516 : Serialize(buf, magicNumber);
517 : }
518 :
519 : procCom.broadcast(buf, rootProc);
520 :
521 : if(!isRoot){
522 : int projectorType;
523 : Deserialize(buf, projectorType);
524 : if(projectorType == 1){
525 : std::string sh_name;
526 : Deserialize(buf, sh_name);
527 : ProjectionHandler* ph;
528 : if (sh_name == std::string(""))
529 : ph = new ProjectionHandler(geometry, subset_handler());
530 : else
531 : {
532 : typedef typename std::map<std::string, SmartPtr<TSubsetHandler> >::const_iterator map_it_t;
533 : map_it_t it = m_additionalSH.begin();
534 : map_it_t it_end = m_additionalSH.end();
535 : for (; it != it_end; ++it)
536 : {
537 : if (it->first == sh_name)
538 : {
539 : ph = new ProjectionHandler(geometry, it->second);
540 : break;
541 : }
542 : }
543 : UG_COND_THROW(it == it_end, "Subset handler name for projection handler not found "
544 : "in list of available names.");
545 : }
546 : SPProjectionHandler projHandler = make_sp(ph);
547 :
548 :
549 : ph->set_default_projector(deserialize_projector(buf));
550 :
551 : size_t numProjectors;
552 : Deserialize(buf, numProjectors);
553 : for(size_t iproj = 0; iproj < numProjectors; ++iproj){
554 : ph->set_projector(iproj, deserialize_projector(buf));
555 : }
556 :
557 : projector = projHandler;
558 : }
559 : else if(projectorType == 0){
560 : projector = deserialize_projector(buf);
561 : }
562 : else if(projectorType == -1){
563 : projector = SPNULL;
564 : }
565 : else{
566 : UG_THROW("Invalid projector type in 'BroadcastRefinementProjector': "
567 : << projectorType);
568 : }
569 :
570 : int tmp;
571 : Deserialize(buf, tmp);
572 : UG_COND_THROW(tmp != magicNumber, "Magic number mismatch in "
573 : "'BroadcastRefinementProjector'. Received "
574 : << tmp << ", but expected " << magicNumber);
575 : }
576 :
577 : return projector;
578 : }
579 :
580 : #endif
581 :
582 :
583 : #ifdef UG_PARALLEL
584 : template <typename TGrid, typename TSubsetHandler>
585 : template <class TElem>
586 : void IDomain<TGrid,TSubsetHandler>::
587 : count_ghosts(std::vector<DomainInfo::int_t>& numGhostsOnLvlOut)
588 : {
589 : TGrid& mg = *m_spGrid;
590 : DistributedGridManager& dgm = *mg.distributed_grid_manager();
591 : GridLayoutMap& glm = dgm.grid_layout_map();
592 :
593 : numGhostsOnLvlOut.clear();
594 : numGhostsOnLvlOut.resize(mg.num_levels(), 0);
595 :
596 : if(glm.has_layout<TElem>(INT_V_MASTER)){
597 : typedef typename GridLayoutMap::Types<TElem>::Layout Layout;
598 : typedef typename GridLayoutMap::Types<TElem>::Interface Interface;
599 : Layout& layout = glm.get_layout<TElem>(INT_V_MASTER);
600 : for(size_t lvl = 0; lvl < layout.num_levels(); ++lvl){
601 : if(lvl >= mg.num_levels())
602 : break;
603 :
604 : typename Layout::LevelLayout& lvlLayout = layout.layout_on_level(lvl);
605 : for(typename Layout::LevelLayout::iterator iiter = lvlLayout.begin();
606 : iiter != lvlLayout.end(); ++iiter)
607 : {
608 : Interface& intfc = lvlLayout.interface(iiter);
609 : for(typename Interface::iterator eiter = intfc.begin();
610 : eiter != intfc.end(); ++eiter)
611 : {
612 : if(dgm.is_ghost(intfc.get_element(eiter)))
613 : ++numGhostsOnLvlOut[lvl];
614 : }
615 : }
616 : }
617 : }
618 : }
619 : #endif
620 :
621 :
622 : template <typename TGrid, typename TSubsetHandler>
623 : template <class TElem>
624 0 : size_t IDomain<TGrid,TSubsetHandler>::
625 : count_local_unique_surface_elements()
626 : {
627 : TGrid& mg = *m_spGrid;
628 : # ifdef UG_PARALLEL
629 : DistributedGridManager* dgm = mg.distributed_grid_manager();
630 : # endif
631 :
632 : size_t numLoc = 0;
633 :
634 : typedef typename TGrid::template traits<TElem>::iterator iter_t;
635 0 : for(iter_t ielem = mg.template begin<TElem>(); ielem != mg.template end<TElem>(); ++ielem)
636 : {
637 : TElem* e = *ielem;
638 :
639 0 : if(!mg.has_children(e)){
640 : #ifdef UG_PARALLEL
641 : if(!(dgm->is_ghost(e) || dgm->contains_status(e, ES_H_SLAVE)))
642 : ++numLoc;
643 : #else
644 0 : ++numLoc;
645 : #endif
646 : }
647 : }
648 :
649 0 : return numLoc;
650 : }
651 :
652 :
653 : ////////////////////////////////////////////////////////////////////////////////
654 : // Domain
655 : ////////////////////////////////////////////////////////////////////////////////
656 :
657 : template <int d, typename TGrid, typename TSubsetHandler>
658 0 : Domain<d,TGrid,TSubsetHandler>::
659 0 : Domain(bool isAdaptive) : IDomain<TGrid, TSubsetHandler>(isAdaptive)
660 : {
661 : // Depending on the dimesion, we'll activeate different options.
662 : // Otherwise we probably would waste memory...
663 : // In any case, sides of elements should always be present
664 : uint gridOpts = GRIDOPT_AUTOGENERATE_SIDES;
665 :
666 : // Furthermore vertices should store associated elements.
667 : // This option depends on the dimension of the domain
668 : if(dim > 0)
669 : gridOpts |= VRTOPT_STORE_ASSOCIATED_EDGES;
670 : if(dim > 1)
671 : gridOpts |= VRTOPT_STORE_ASSOCIATED_FACES;
672 : if(dim > 2)
673 : gridOpts |= VRTOPT_STORE_ASSOCIATED_VOLUMES;
674 :
675 : // thats it for now. One could think about enabling
676 : // FACEOPT_STORE_ASSOCIATED_EDGES, VOLOPT_STORE_ASSOCIATED_EDGES
677 : // and VOLOPT_STORE_ASSOCIATED_FACES. However this costs considerably
678 : // more memory compared to the performance benefits.
679 : // Now set the options
680 0 : this->grid()->set_options(gridOpts);
681 :
682 : // get position attachment
683 : m_aPos = GetDefaultPositionAttachment<position_attachment_type>();
684 :
685 : // let position accessor access Vertex Coordinates
686 0 : if(!this->grid()->template has_attachment<Vertex>(m_aPos))
687 0 : this->grid()->template attach_to<Vertex>(m_aPos);
688 0 : m_aaPos.access(*(this->grid()), m_aPos);
689 :
690 0 : m_geometry3d = MakeGeometry3d(*(this->grid()), m_aPos);
691 0 : this->m_refinementProjector = make_sp(new RefinementProjector(m_geometry3d));
692 0 : }
693 :
694 :
695 : } // end namespace ug
696 :
697 : #endif /* __H__UG__LIB_DISC__DOMAIN_IMPL__ */
|
