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 "debug_util.h"
34 : #include "attachment_util.h"
35 :
36 : #ifdef UG_PARALLEL
37 : #include "lib_grid/parallelization/distributed_grid.h"
38 : #include "lib_grid/parallelization/util/compol_copy_attachment.h"
39 : #include "lib_grid/parallelization/util/compol_attachment_reduce.h"
40 : #include "pcl/pcl.h"
41 : #endif
42 :
43 : using namespace std;
44 :
45 : namespace ug{
46 :
47 0 : void PrintElementNumbers(const GridObjectCollection& goc)
48 : {
49 : UG_LOG("grid element numbers:\n");
50 0 : for(size_t i = 0; i < goc.num_levels(); ++i)
51 : {
52 0 : if(goc.num_levels() > 1){
53 : UG_LOG("level " << i << endl);
54 : }
55 : UG_LOG(" vertices total:\t" << goc.num<Vertex>(i) << endl);
56 0 : if(goc.num<Vertex>() > 0){
57 : UG_LOG(" normal vrts:\t" << goc.num<RegularVertex>(i) << endl);
58 : UG_LOG(" hanging vrts:\t" << goc.num<ConstrainedVertex>(i) << endl);
59 : }
60 :
61 : UG_LOG(" edges total:\t\t" << goc.num<Edge>(i) << endl);
62 0 : if(goc.num<Edge>() > 0){
63 : UG_LOG(" normal edges:\t" << goc.num<RegularEdge>(i) << endl);
64 : UG_LOG(" constraining edges:\t" << goc.num<ConstrainingEdge>(i) << endl);
65 : UG_LOG(" constrained edges:\t" << goc.num<ConstrainedEdge>(i) << endl);
66 : }
67 :
68 : UG_LOG(" faces total:\t\t" << goc.num<Face>(i) << endl);
69 0 : if(goc.num<Face>() > 0){
70 : UG_LOG(" normal triangles:\t" << goc.num<Triangle>(i) << endl);
71 : UG_LOG(" constraining tris:\t" << goc.num<ConstrainingTriangle>(i) << endl);
72 : UG_LOG(" constrained tris:\t" << goc.num<ConstrainedTriangle>(i) << endl);
73 :
74 : UG_LOG(" normal quads:\t" << goc.num<Quadrilateral>(i) << endl);
75 : UG_LOG(" constraining quads:\t" << goc.num<ConstrainingQuadrilateral>(i) << endl);
76 : UG_LOG(" constrained quads:\t" << goc.num<ConstrainedQuadrilateral>(i) << endl);
77 : }
78 :
79 : UG_LOG(" volumes total:\t" << goc.num<Volume>(i) << endl);
80 0 : if(goc.num<Volume>() > 0){
81 : UG_LOG(" tetrahedrons:\t" << goc.num<Tetrahedron>(i) << endl);
82 : UG_LOG(" pyramids:\t" << goc.num<Pyramid>(i) << endl);
83 : UG_LOG(" prisms:\t" << goc.num<Prism>(i) << endl);
84 : UG_LOG(" hexahedrons:\t" << goc.num<Hexahedron>(i) << endl);
85 : }
86 :
87 : UG_LOG(endl);
88 : }
89 0 : }
90 :
91 0 : void PrintGridElementNumbers(Grid& grid)
92 : {
93 0 : PrintElementNumbers(grid.get_grid_objects());
94 0 : }
95 :
96 0 : void PrintGridElementNumbers(MultiGrid& mg)
97 : {
98 0 : PrintElementNumbers(mg.get_grid_objects());
99 0 : }
100 :
101 0 : void PrintGridElementNumbers(GridSubsetHandler& sh)
102 : {
103 0 : PrintElementNumbers(sh.get_grid_objects());
104 0 : }
105 :
106 : template <class TGeomObj>
107 0 : void PrintAttachmentInfo(Grid& grid)
108 : {
109 : typedef typename Grid::traits<TGeomObj>::AttachmentPipe AttachmentPipe;
110 : typedef typename AttachmentPipe::ConstAttachmentEntryIterator AttIter;
111 :
112 : // iterate over all attachments of the grid
113 : AttachmentPipe& pipe = grid.get_attachment_pipe<TGeomObj>();
114 :
115 : int counter = 1;
116 : size_t totalSize = 0;
117 0 : for(AttIter iter = pipe.attachments_begin();
118 0 : iter != pipe.attachments_end(); ++iter, ++counter)
119 : {
120 : // name
121 0 : IAttachment* att = iter->m_pAttachment;
122 0 : UG_LOG("Attachment " << counter << " (" << att->get_name() << "): ");
123 :
124 : // size
125 0 : IAttachmentDataContainer* con = iter->m_pContainer;
126 0 : UG_LOG(con->occupied_memory() << " bytes\n");
127 0 : totalSize += con->occupied_memory();
128 : }
129 :
130 0 : UG_LOG(counter - 1 << " attachments with a total size of "
131 : << totalSize << " bytes.\n");
132 0 : }
133 :
134 0 : void PrintAttachmentInfo(Grid& grid)
135 : {
136 : UG_LOG("Vertex Attachments:\n");
137 0 : PrintAttachmentInfo<Vertex>(grid);
138 :
139 : UG_LOG("\nEdge Attachments:\n");
140 0 : PrintAttachmentInfo<Edge>(grid);
141 :
142 : UG_LOG("\nFace Attachments:\n");
143 0 : PrintAttachmentInfo<Face>(grid);
144 :
145 : UG_LOG("\nVolume Attachments:\n");
146 0 : PrintAttachmentInfo<Volume>(grid);
147 0 : }
148 :
149 : template <class TElem>
150 0 : static void CheckMultiGridConsistencyImpl(MultiGrid& mg)
151 : {
152 : #ifdef UG_PARALLEL
153 : DistributedGridManager* dgm = mg.distributed_grid_manager();
154 : #endif
155 :
156 0 : for(size_t lvl = 0; lvl < mg.num_levels(); ++lvl){
157 0 : for(typename MultiGrid::traits<TElem>::iterator iter = mg.begin<TElem>(lvl);
158 0 : iter != mg.end<TElem>(lvl); ++iter)
159 : {
160 : TElem* e = *iter;
161 : // make sure that all children have the local element as parent
162 0 : for(size_t i_child = 0; i_child < mg.num_children<Vertex>(e); ++i_child)
163 : {
164 0 : if(e != mg.get_parent(mg.get_child<Vertex>(e, i_child))){
165 0 : UG_THROW("parent is not referenced by children!");
166 : }
167 : }
168 :
169 0 : for(size_t i_child = 0; i_child < mg.num_children<Edge>(e); ++i_child)
170 : {
171 0 : if(e != mg.get_parent(mg.get_child<Edge>(e, i_child))){
172 0 : UG_THROW("parent is not referenced by children!");
173 : }
174 : }
175 :
176 0 : for(size_t i_child = 0; i_child < mg.num_children<Face>(e); ++i_child)
177 : {
178 0 : if(e != mg.get_parent(mg.get_child<Face>(e, i_child))){
179 0 : UG_THROW("parent is not referenced by children!");
180 : }
181 : }
182 :
183 0 : for(size_t i_child = 0; i_child < mg.num_children<Volume>(e); ++i_child)
184 : {
185 0 : if(e != mg.get_parent(mg.get_child<Volume>(e, i_child))){
186 0 : UG_THROW("parent is not referenced by children!");
187 : }
188 : }
189 :
190 : // also make sure that each element with a parent is contained in the
191 : // children list of its parent
192 : GridObject* parent = mg.get_parent(e);
193 0 : if(parent){
194 : bool gotIt = false;
195 0 : for(size_t i = 0; i < mg.num_children<TElem>(parent); ++i){
196 0 : if(mg.get_child<TElem>(parent, i) == e){
197 : gotIt = true;
198 : break;
199 : }
200 : }
201 :
202 0 : if(!gotIt){
203 0 : UG_THROW("Element not contained in child list of its parent");
204 : }
205 : }
206 : else{
207 0 : if(lvl > 0){
208 : #ifdef UG_PARALLEL
209 : if(!dgm->contains_status(e, INT_V_SLAVE))
210 : {
211 : UG_THROW("Each element in a higher level has to have a parent"
212 : " unless it is a v-slave!");
213 : }
214 : #else
215 0 : UG_THROW("Each element in a higher level has to have a parent"
216 : " unless it is a v-slave!");
217 : #endif
218 : }
219 : }
220 : }
221 : }
222 0 : }
223 :
224 0 : void CheckMultiGridConsistency(MultiGrid& mg)
225 : {
226 0 : CheckMultiGridConsistencyImpl<Vertex>(mg);
227 0 : CheckMultiGridConsistencyImpl<Edge>(mg);
228 0 : CheckMultiGridConsistencyImpl<Face>(mg);
229 0 : CheckMultiGridConsistencyImpl<Volume>(mg);
230 0 : }
231 :
232 : // the following define is only used by the CheckHangingNodeConsistency methods
233 : // and is highly specialized for them.
234 : #define FAILED_CHECK(elem, msg)\
235 : {UG_ERR_LOG("CheckHangingNodeConsistency: " << msg << endl);\
236 : UG_ERR_LOG(" at: " << GetGridObjectCenter(g, elem));\
237 : if(MultiGrid* mg = dynamic_cast<MultiGrid*>(&g))\
238 : {UG_ERR_LOG(" on level " << mg->get_level(elem));}\
239 : UG_ERR_LOG(endl);\
240 : isConsistent = false;\
241 : }//continue;}
242 :
243 :
244 : template<typename T>
245 0 : void CheckHangingNodeConstrainingFace(bool isConsistent, Grid &g, T iter, T end)
246 : {
247 : Grid::edge_traits::secure_container edges;
248 :
249 0 : for(; iter != end; ++iter)
250 : {
251 : ConstrainingFace* cf = *iter;
252 :
253 : // make sure that all associated edges are constraining, too
254 : g.associated_elements(edges, cf);
255 0 : for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
256 0 : if(!edges[i_edge]->is_constraining()){
257 0 : FAILED_CHECK(edges[i_edge], "Non constraining side of constraining face detected!");
258 : }
259 : }
260 : }
261 0 : }
262 :
263 :
264 0 : bool CheckHangingNodeConsistency(Grid& g)
265 : {
266 : bool isConsistent = true;
267 :
268 : // iterate over all hanging nodes and check whether the associated parent
269 : // contains the node in its list of constraiend objects
270 : for(Grid::traits<ConstrainedVertex>::iterator iter = g.begin<ConstrainedVertex>();
271 0 : iter != g.end<ConstrainedVertex>(); ++iter)
272 : {
273 : ConstrainedVertex* hnode = *iter;
274 : GridObject* constrObj = hnode->get_constraining_object();
275 0 : if(!constrObj){
276 0 : FAILED_CHECK(hnode, "Hanging Vertex has no constraining object!");
277 : }
278 :
279 0 : if(ConstrainingEdge* ce = dynamic_cast<ConstrainingEdge*>(constrObj)){
280 : // check whether hnode is a constrained object of ce
281 0 : if(!ce->is_constrained_object(hnode)){
282 0 : FAILED_CHECK(hnode, "Hanging Vertex is not constrained by parent edge!");
283 : }
284 : }
285 0 : else if(ConstrainingFace* cf = dynamic_cast<ConstrainingFace*>(constrObj)){
286 : // check whether hnode is a constraiend object of cf
287 0 : if(!cf->is_constrained_object(hnode)){
288 0 : FAILED_CHECK(hnode, "Hanging Vertex is not constrained by parent face!");
289 : }
290 : }
291 : else{
292 0 : FAILED_CHECK(hnode, "Parent of Hanging Vertex is not a constraining object!");
293 : }
294 : }
295 :
296 : // iterate over all constrained edges and check whether the associated constraining
297 : // object contains the edge in its list of constraiend objects
298 : for(Grid::traits<ConstrainedEdge>::iterator iter = g.begin<ConstrainedEdge>();
299 0 : iter != g.end<ConstrainedEdge>(); ++iter)
300 : {
301 : ConstrainedEdge* e = *iter;
302 : GridObject* parent = e->get_constraining_object();
303 :
304 0 : if(!parent){
305 0 : FAILED_CHECK(e, "Constrained Edge has no parent!");
306 : }
307 :
308 0 : if(!parent->is_constraining()){
309 0 : FAILED_CHECK(e, "Parent of Constrained Edge is not a constraining object!");
310 : }
311 :
312 0 : if(ConstrainingEdge* ce = dynamic_cast<ConstrainingEdge*>(parent)){
313 : // check whether e is a constraiend object of ce
314 0 : if(!ce->is_constrained_object(e)){
315 0 : FAILED_CHECK(e, "Constrained Edge is not constrained by parent edge!");
316 : }
317 :
318 : // since the constraining object is an edge, we'll make sure, that the constrained edge
319 : // is connected to exactly one constrained vertex.
320 0 : if(e->vertex(0)->is_constrained() && e->vertex(1)->is_constrained()){
321 0 : FAILED_CHECK(e, "Constrained Edge (constrained by edge) "
322 : " may not be connected to two constrained vertices!");
323 : }
324 : }
325 0 : else if(ConstrainingFace* cf = dynamic_cast<ConstrainingFace*>(parent)){
326 : // check whether hnode is a constraiend object of cf
327 0 : if(!cf->is_constrained_object(e)){
328 0 : FAILED_CHECK(e, "Constrained Edge is not constrained by parent face!");
329 : }
330 : }
331 : else{
332 0 : FAILED_CHECK(e, "Unknown parent type of Constrained Edge!");
333 : }
334 :
335 : // make sure that the edge connected to at least one constrained vertex
336 : // we check above, that an edge which is constrained by an edge is not
337 : // connected to two constrained vertices.
338 0 : if(!(e->vertex(0)->is_constrained() || e->vertex(1)->is_constrained())){
339 0 : FAILED_CHECK(e, "Constrained Edge has to be connected to a constrained vertex!");
340 : }
341 : }
342 :
343 : // CHECK FACES
344 : // fix since assign does not seem to work on some machines.
345 0 : CheckHangingNodeConstrainingFace(isConsistent, g, g.begin<ConstrainingTriangle>(), g.end<ConstrainingTriangle>());
346 0 : CheckHangingNodeConstrainingFace(isConsistent, g, g.begin<ConstrainingQuadrilateral>(), g.end<ConstrainingQuadrilateral>());
347 :
348 0 : return isConsistent;
349 : }
350 :
351 0 : bool CheckHangingNodeConsistency(MultiGrid& mg)
352 : {
353 0 : Grid& g = mg;
354 0 : bool isConsistent = CheckHangingNodeConsistency(g);
355 :
356 : // iterate over all hanging nodes and check whether the associated parent
357 : // matches the constraining object. Other checks have already been performed!
358 : for(Grid::traits<ConstrainedVertex>::iterator iter = g.begin<ConstrainedVertex>();
359 0 : iter != g.end<ConstrainedVertex>(); ++iter)
360 : {
361 : ConstrainedVertex* hnode = *iter;
362 : GridObject* co = hnode->get_constraining_object();
363 : GridObject* parent = mg.get_parent(hnode);
364 :
365 0 : if(co != parent){
366 0 : FAILED_CHECK(hnode, "Hanging Vertex: Constraining object and parent do not match!");
367 : }
368 :
369 0 : if(co && ((mg.get_level(co) + 1) != mg.get_level(hnode))){
370 0 : FAILED_CHECK(hnode, "Hanging Vertex has to be one level higher than its constraining object!");
371 : }
372 : }
373 :
374 : // iterate over all edges if an edge has children, then collect its ass. faces.
375 : // If not all such faces have children, then the edge has to be a constraining
376 : // edge and its children have to be constrained.
377 : vector<Face*> faces;
378 : for(EdgeIterator iter = mg.begin<Edge>();
379 0 : iter != mg.end<Edge>(); ++iter)
380 : {
381 : Edge* e = *iter;
382 :
383 0 : if(mg.get_parent(e) && (mg.get_level(mg.get_parent(e)) + 1 != mg.get_level(e))){
384 0 : FAILED_CHECK(e, "Edge and parent are not in consecutive levels!");
385 : }
386 :
387 0 : for(size_t i = 0; i < 2; ++i){
388 0 : if(mg.get_level(e) != mg.get_level(e->vertex(i))){
389 0 : FAILED_CHECK(e, "Edge-Vertex is not on the same level as the edge itself!");
390 0 : UG_ERR_LOG(" Vertex at " << GetGridObjectCenter(mg, e->vertex(i))
391 : << " on level " << mg.get_level(e->vertex(i)) << endl);
392 : }
393 : }
394 :
395 : if(mg.has_children<Edge>(e)){
396 : // e may not be a constrained edge
397 0 : if(e->is_constrained()){
398 0 : FAILED_CHECK(e, "Constrained Edge may not have children!");
399 : }
400 : /*
401 : // check whether all ass. faces have children.
402 : bool allNbrsAreParents= true;
403 : bool hasConstrainingNbrFaces = false;
404 : CollectAssociated(faces, mg, e);
405 : for(size_t i = 0; i < faces.size(); ++i){
406 : if(!mg.has_children<Face>(faces[i]))
407 : allNbrsAreParents = false;
408 : if(faces[i]->is_constraining())
409 : hasConstrainingNbrFaces = true;
410 : }
411 : */
412 :
413 : // depending whether e is constrained or normal, either all nbrs
414 : // or only some have to have children.
415 0 : if(e->is_constraining()){
416 : /* This causes false positives for distributed grids!
417 : if(allNbrsAreParents && (!hasConstrainingNbrFaces)){
418 : // e should be a normal edge
419 : FAILED_CHECK(e, "At least one neighbor face of a"
420 : " constraining edge should not have children!");
421 : }
422 : else{*/
423 : // make sure that e has two constrained edge-children and a
424 : // constrained vertex child.
425 0 : ConstrainingEdge* ce = dynamic_cast<ConstrainingEdge*>(e);
426 : UG_ASSERT(ce, "Only ConstrainingEdges should return true in Edge::is_constrained()!");
427 0 : if(ce){
428 0 : if(ce->num_constrained_edges() != 2){
429 0 : FAILED_CHECK(e, "ConstrainingEdge has to constrain 2 edges not "
430 : << ce->num_constrained_edges() << "!");
431 : }
432 :
433 0 : if(!(ce->constrained_edge(0)->is_constrained())
434 0 : && ce->constrained_edge(1)->is_constrained())
435 : {
436 0 : FAILED_CHECK(e, "Both edges constrained by "
437 : " a constraining edge have to be ConstrainedEdges!");
438 : }
439 :
440 0 : if(ce->num_constrained_vertices() != 1){
441 0 : FAILED_CHECK(e, "ConstrainingEdge has to constrain 1 vertex, not "
442 : << ce->num_constrained_vertices() << "!");
443 : }
444 :
445 0 : if(!ce->constrained_vertex(0)->is_constrained()){
446 0 : FAILED_CHECK(e, "a vertex constrained by a constraining "
447 : "edge has to be a hangingVertex!");
448 : }
449 : }
450 : //}
451 : }
452 : else{
453 : // e is not constraining
454 : // all neighbors should be parents, too!
455 : /*this check is only valid if we don't create a closure
456 : if(!allNbrsAreParents){
457 : // e should be a normal edge
458 : FAILED_CHECK(e, "All neighbor faces of a normal edge should have children!");
459 : }*/
460 : }
461 : }
462 : }
463 : // UG_LOG("2\n");
464 : // check faces
465 0 : for(FaceIterator iter = g.begin<Face>(); iter != g.end<Face>(); ++iter){
466 : Face* f = *iter;
467 :
468 0 : if(f->is_constraining()){
469 : // UG_LOG("2.1/n");
470 0 : ConstrainingFace* cf = dynamic_cast<ConstrainingFace*>(f);
471 : UG_ASSERT(cf, "All constraining faces should derive from ConstrainingFace");
472 :
473 : // make sure that the face has the right number of children
474 0 : if(mg.num_children<Face>(f) != 4){
475 0 : FAILED_CHECK(f, "Face has bad number of child faces. "
476 : "4 required, " << mg.num_children<Face>(f) << " found.");
477 : }
478 :
479 0 : if(mg.num_children<Edge>(f) != f->num_vertices()){
480 0 : FAILED_CHECK(f, "Face has bad number of child edges. "
481 : << f->num_vertices() << " required, "
482 : << mg.num_children<Edge>(f) << " found.");
483 : }
484 :
485 0 : if((f->num_vertices() > 3) && mg.num_children<Vertex>(f) != 1){
486 0 : FAILED_CHECK(f, "Face has bad number of child vertices. "
487 : "1 required, " << mg.num_children<Vertex>(f) << " found.");
488 : }
489 :
490 : // make sure that number of children and number of constrained elements match
491 0 : if(mg.num_children<Face>(f) != cf->num_constrained_faces())
492 0 : FAILED_CHECK(f, "Number of child faces of constraining face does not match number of constrained faces.");
493 :
494 0 : if(mg.num_children<Edge>(f) != cf->num_constrained_edges())
495 0 : FAILED_CHECK(f, "Number of child edges of constraining face does not match number of constrained edges.");
496 :
497 0 : if(mg.num_children<Vertex>(f) != cf->num_constrained_vertices())
498 0 : FAILED_CHECK(f, "Number of child vertices of constraining face does not match number of constrained vertices.");
499 :
500 : // make sure that all children are constrained and that they are contained
501 : // in the list of constrained elements
502 0 : for(size_t i = 0; i < mg.num_children<Face>(f); ++i){
503 : Face* child = mg.get_child<Face>(f, i);
504 0 : if(!child->is_constrained()){
505 0 : FAILED_CHECK(f, "All child faces of a constraining face have to be constrained faces.");
506 : }
507 : else{
508 : // make sure that the child is contained in the list of constrained objects of cf
509 0 : if(!cf->is_constrained_object(child))
510 0 : FAILED_CHECK(f, "Child face of constraining face is not in list of constrained faces.");
511 : }
512 : }
513 0 : for(size_t i = 0; i < mg.num_children<Edge>(f); ++i){
514 : Edge* child = mg.get_child<Edge>(f, i);
515 0 : if(!child->is_constrained()){
516 0 : FAILED_CHECK(f, "All child edges of a constraining face have to be constrained edges.");
517 : }
518 : else{
519 : // make sure that the child is contained in the list of constrained objects of cf
520 0 : if(!cf->is_constrained_object(child))
521 0 : FAILED_CHECK(f, "Child edge of constraining face is not in list of constrained edges.");
522 : }
523 : }
524 0 : for(size_t i = 0; i < mg.num_children<Vertex>(f); ++i){
525 : Vertex* child = mg.get_child<Vertex>(f, i);
526 0 : if(!child->is_constrained()){
527 0 : FAILED_CHECK(f, "All child vertices of a constraining face have to be constrained vertices.");
528 : }
529 : else{
530 : // make sure that the child is contained in the list of constrained objects of cf
531 0 : if(!cf->is_constrained_object(child))
532 0 : FAILED_CHECK(f, "Child vertex of constraining face is not in list of constrained vertices.");
533 : }
534 : }
535 : }
536 :
537 0 : else if(f->is_constrained()){
538 : // UG_LOG("2.2.0/n");
539 0 : ConstrainedFace* cdf = dynamic_cast<ConstrainedFace*>(f);
540 : UG_ASSERT(cdf, "All constrained faces should derive from ConstrainedFace");
541 : // UG_LOG("2.2.1/n");
542 : // we don't have to check all interconnections, since we already checked a lot of
543 : // stuff for constraining faces. So just do the rest now.
544 0 : ConstrainingFace* cf = dynamic_cast<ConstrainingFace*>(cdf->get_constraining_object());
545 : // UG_LOG("2.2.2/n");
546 0 : if(!cf){
547 0 : FAILED_CHECK(cdf, "No constraining face found for given constrained face.");
548 : }
549 :
550 0 : if(cf != mg.get_parent(cdf)){
551 0 : FAILED_CHECK(cdf, "The parent of a constrained face should always be its constraining face!");
552 : }
553 : }
554 : }
555 : // UG_LOG("3\n");
556 0 : return isConsistent;
557 0 : }
558 :
559 :
560 : ////////////////////////////////////////////////////////////////////////////////
561 : ////////////////////////////////////////////////////////////////////////////////
562 : enum ConstraintTypes{
563 : CT_NONE = 0,
564 : CT_CONSTRAINING = 1,
565 : CT_CONSTRAINED = 1 << 1
566 : };
567 :
568 : template <class TElem>
569 0 : static bool CheckDistributedObjectConstraintTypes(MultiGrid& mg)
570 : {
571 : typedef typename Grid::traits<TElem>::iterator ElemIter;
572 : bool retVal = true;
573 :
574 : AInt aState;
575 0 : mg.attach_to<TElem>(aState);
576 : Grid::AttachmentAccessor<TElem, AInt> aaState(mg, aState);
577 :
578 : // set up local states
579 0 : for(ElemIter iter = mg.begin<TElem>(); iter != mg.end<TElem>(); ++iter){
580 : TElem* e = *iter;
581 0 : aaState[e] = CT_NONE;
582 0 : if(e->is_constraining())
583 0 : aaState[e] |= CT_CONSTRAINING;
584 0 : if(e->is_constrained())
585 0 : aaState[e] |= CT_CONSTRAINED;
586 : }
587 :
588 : // communicate states
589 : #ifdef UG_PARALLEL
590 : typedef typename GridLayoutMap::Types<TElem>::Layout Layout;
591 : pcl::InterfaceCommunicator<Layout> com;
592 : DistributedGridManager& distGridMgr = *mg.distributed_grid_manager();
593 : GridLayoutMap& glm = distGridMgr.grid_layout_map();
594 :
595 : ComPol_AttachmentReduce<Layout, AInt> compolOr(mg, aState, PCL_RO_BOR);
596 : com.exchange_data(glm, INT_H_SLAVE, INT_H_MASTER, compolOr);
597 : com.communicate();
598 :
599 : com.exchange_data(glm, INT_H_MASTER, INT_H_SLAVE, compolOr);
600 : com.communicate();
601 :
602 : com.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, compolOr);
603 : com.communicate();
604 :
605 : com.exchange_data(glm, INT_V_MASTER, INT_V_SLAVE, compolOr);
606 : com.communicate();
607 :
608 : com.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, compolOr);
609 : com.communicate();
610 : #endif
611 :
612 : // check whether communicated states match the actual element states on this proc.
613 0 : for(ElemIter iter = mg.begin<TElem>(); iter != mg.end<TElem>(); ++iter){
614 : TElem* e = *iter;
615 : int state = CT_NONE;
616 0 : if(e->is_constraining())
617 : state |= CT_CONSTRAINING;
618 0 : if(e->is_constrained())
619 0 : state |= CT_CONSTRAINED;
620 :
621 0 : if(state != aaState[e]){
622 0 : UG_LOG("ERROR: Distributed object has different constraint states on different procs. "
623 : << "At: " << GetGridObjectCenter(mg, e)
624 : << " on level " << mg.get_level(e) << endl);
625 : retVal = false;
626 : }
627 : }
628 :
629 : mg.detach_from<TElem>(aState);
630 0 : return retVal;
631 : }
632 :
633 0 : bool CheckDistributedObjectConstraintTypes(MultiGrid& mg)
634 : {
635 : bool retVal = true;
636 :
637 : // assign constraint states
638 : UG_LOG("Checking constraint types of VERTICES\n");
639 0 : retVal &= CheckDistributedObjectConstraintTypes<Vertex>(mg);
640 : UG_LOG("Checking constraint types of EDGES\n");
641 0 : retVal &= CheckDistributedObjectConstraintTypes<Edge>(mg);
642 : UG_LOG("Checking constraint types of FACES\n");
643 0 : retVal &= CheckDistributedObjectConstraintTypes<Face>(mg);
644 : UG_LOG("Checking constraint types DONE\n");
645 :
646 : // make sure that we return the same value globally
647 : #ifdef UG_PARALLEL
648 : retVal = pcl::AllProcsTrue(retVal);
649 : #endif
650 :
651 0 : return retVal;
652 : }
653 :
654 : #ifdef UG_PARALLEL
655 : template <class TLayout>
656 : class ComPol_CheckDistributedParentStates : public pcl::ICommunicationPolicy<TLayout>
657 : {
658 : public:
659 : typedef TLayout Layout;
660 : typedef typename Layout::Type GeomObj;
661 : typedef typename Layout::Element Element;
662 : typedef typename Layout::Interface Interface;
663 : typedef typename Interface::const_iterator InterfaceIter;
664 :
665 : ComPol_CheckDistributedParentStates(MultiGrid& mg) :
666 : m_mg(mg),
667 : m_dgm(mg.distributed_grid_manager()),
668 : m_comparisionFailed(false),
669 : m_performMasterCheck(false)
670 : {
671 : }
672 :
673 : virtual ~ComPol_CheckDistributedParentStates() {}
674 :
675 : virtual int
676 : get_required_buffer_size(const Interface& interface)
677 : {return interface.size() * sizeof(int);}
678 :
679 : virtual bool
680 : collect(ug::BinaryBuffer& buff, const Interface& interface)
681 : {
682 : for(InterfaceIter iter = interface.begin();
683 : iter != interface.end(); ++iter)
684 : {
685 : Element elem = interface.get_element(iter);
686 : int parentType = m_mg.parent_type(elem);
687 : Serialize(buff, parentType);
688 : }
689 : return true;
690 : }
691 :
692 : virtual bool
693 : extract(ug::BinaryBuffer& buff, const Interface& interface)
694 : {
695 : for(InterfaceIter iter = interface.begin();
696 : iter != interface.end(); ++iter)
697 : {
698 : Element elem = interface.get_element(iter);
699 : int parentType;
700 : Deserialize(buff, parentType);
701 :
702 : if(m_performMasterCheck){
703 : GridObject* parent = m_mg.get_parent(elem);
704 : if(parent && (parentType != parent->base_object_id())){
705 : UG_LOG(" PARENT-TYPE MISMATCH AT CHILD ELEMENT WITH OBJECT ID " << elem->base_object_id()
706 : << " at " << GetGridObjectCenter(m_mg, elem) << " on level " << m_mg.get_level(elem) << "\n");
707 : UG_LOG(" Parent object id: " << parent->base_object_id() << ", received id: " << parentType << "\n");
708 : m_comparisionFailed = true;
709 : }
710 : else if((!parent) && (parentType != -1)){
711 : UG_LOG(" PARENT-TYPE MISMATCH AT CHILD ELEMENT WITH OBJECT ID " << elem->base_object_id()
712 : << " at " << GetGridObjectCenter(m_mg, elem) << " on level " << m_mg.get_level(elem) << "\n");
713 : UG_LOG(" The element hasn't got a parent but received parent id is != -1. Received id: " << parentType << "\n");
714 : m_comparisionFailed = true;
715 : }
716 : }
717 : else if(parentType != m_mg.parent_type(elem)){
718 : UG_LOG(" PARENT-TYPE MISMATCH AT ELEMENT WITH OBJECT ID " << elem->base_object_id()
719 : << " at " << GetGridObjectCenter(m_mg, elem) << " on level " << m_mg.get_level(elem) << "\n");
720 : UG_LOG(" Parent object id: " << (int)m_mg.parent_type(elem) << ", received id: " << parentType << "\n");
721 : m_comparisionFailed = true;
722 : }
723 : }
724 : return true;
725 : }
726 :
727 : bool exchange_data()
728 : {
729 : pcl::InterfaceCommunicator<TLayout> com;
730 :
731 : m_comparisionFailed = false;
732 : GridLayoutMap& glm = m_dgm->grid_layout_map();
733 : m_performMasterCheck = false;
734 : if(glm.has_layout<GeomObj>(INT_H_SLAVE))
735 : com.send_data(glm.get_layout<GeomObj>(INT_H_SLAVE), *this);
736 : if(glm.has_layout<GeomObj>(INT_H_MASTER))
737 : com.receive_data(glm.get_layout<GeomObj>(INT_H_MASTER), *this);
738 : com.communicate();
739 :
740 : m_performMasterCheck = true;
741 : if(glm.has_layout<GeomObj>(INT_V_SLAVE))
742 : com.send_data(glm.get_layout<GeomObj>(INT_V_SLAVE), *this);
743 : if(glm.has_layout<GeomObj>(INT_V_MASTER))
744 : com.receive_data(glm.get_layout<GeomObj>(INT_V_MASTER), *this);
745 : com.communicate();
746 :
747 : return !m_comparisionFailed;
748 : }
749 :
750 : private:
751 : MultiGrid& m_mg;
752 : DistributedGridManager* m_dgm;
753 : bool m_comparisionFailed;
754 : bool m_performMasterCheck;
755 : };
756 : #endif
757 :
758 : template <class TElem>
759 : bool CheckLocalParentTypes(MultiGrid& mg)
760 : {
761 : bool success = true;
762 : typedef typename MultiGrid::traits<TElem>::iterator iter_t;
763 : for(iter_t iter = mg.begin<TElem>(); iter != mg.end<TElem>(); ++iter){
764 : TElem* e = *iter;
765 : GridObject* parent = mg.get_parent(e);
766 : if(parent){
767 : if(mg.parent_type(e) != parent->base_object_id()){
768 : UG_LOG(" LOCAL PARENT-TYPE MISMATCH AT ELEMENT WITH OBJECT ID " << e->base_object_id()
769 : << " at " << GetGridObjectCenter(mg, e) << " on level " << mg.get_level(e) << "\n");
770 : UG_LOG(" Stored parent id: " << (int)mg.parent_type(e) << ", actual parent id: " << parent->base_object_id() << "\n");
771 : success = false;
772 : }
773 : }
774 : }
775 : return success;
776 : }
777 :
778 0 : bool CheckDistributedParentTypes(MultiGrid& mg)
779 : {
780 : UG_LOG("DEBUG: Checking distributed parent types...\n");
781 : #ifdef UG_PARALLEL
782 : ComPol_CheckDistributedParentStates<VertexLayout> vrtChecker(mg);
783 : ComPol_CheckDistributedParentStates<EdgeLayout> edgeChecker(mg);
784 : ComPol_CheckDistributedParentStates<FaceLayout> faceChecker(mg);
785 : ComPol_CheckDistributedParentStates<VolumeLayout> volChecker(mg);
786 :
787 : bool success = true;
788 : UG_LOG(" checking vertices...\n");
789 : success &= CheckLocalParentTypes<Vertex>(mg);
790 : success &= vrtChecker.exchange_data();
791 : UG_LOG(" checking edges...\n");
792 : success &= CheckLocalParentTypes<Edge>(mg);
793 : success &= edgeChecker.exchange_data();
794 : UG_LOG(" checking faces...\n");
795 : success &= CheckLocalParentTypes<Face>(mg);
796 : success &= faceChecker.exchange_data();
797 : UG_LOG(" checking volumes...\n");
798 : success &= CheckLocalParentTypes<Volume>(mg);
799 : success &= volChecker.exchange_data();
800 : UG_LOG(" checking done with status ");
801 :
802 : success = pcl::AllProcsTrue(success);
803 :
804 : if(success){
805 : UG_LOG("SUCCESS\n");
806 : }
807 : else{
808 : UG_LOG("FAILURE\n");
809 : }
810 : return success;
811 : #else
812 0 : return true;
813 : #endif
814 : }
815 :
816 :
817 0 : bool CheckElementConsistency(MultiGrid& mg, Vertex* v)
818 : {
819 : bool success = true;
820 0 : UG_LOG("DEBUG: Checking vertex at " << GetGridObjectCenter(mg, v) << endl);
821 : UG_LOG(" vertex type:");
822 0 : if(v->is_constrained()){
823 : UG_LOG(" constrained");
824 : }
825 : else{
826 : UG_LOG(" normal");
827 : }
828 : UG_LOG("\n");
829 :
830 0 : if(v->is_constrained()){
831 : ConstrainedVertex* cdv = dynamic_cast<ConstrainedVertex*>(v);
832 : (void) cdv; // removes unused warning
833 : UG_ASSERT(cdv, "Bad type!");
834 : UG_ASSERT(cdv->get_constraining_object() == mg.get_parent(cdv),
835 : "ConstrainingObject / Parent mismatch");
836 : }
837 : else{
838 : }
839 :
840 0 : return success;
841 : }
842 :
843 0 : bool CheckElementConsistency(MultiGrid& mg, Edge* e)
844 : {
845 : bool success = true;
846 0 : UG_LOG("DEBUG: Checking edge at " << GetGridObjectCenter(mg, e) << endl);
847 : UG_LOG(" edge type:");
848 0 : if(e->is_constrained()){
849 : UG_LOG(" constrained");
850 : }
851 0 : else if(e->is_constraining()){
852 : UG_LOG(" constraining");
853 : }
854 : else{
855 : UG_LOG(" normal");
856 : }
857 : UG_LOG("\n");
858 :
859 : // constrained/constraining checks
860 0 : if(e->is_constrained()){
861 : ConstrainedEdge* cde = dynamic_cast<ConstrainedEdge*>(e);
862 : (void) cde; // removes unused warning
863 : UG_ASSERT(cde, "Bad type!");
864 : UG_ASSERT(cde->get_constraining_object() == mg.get_parent(cde),
865 : "ConstrainingObject / Parent mismatch");
866 :
867 : UG_ASSERT(cde->vertex(0)->is_constrained() || cde->vertex(1)->is_constrained(),
868 : "At least one corner of a constrained edge has to be a constrained vertex");
869 : }
870 0 : else if(e->is_constraining()){
871 0 : ConstrainingEdge* cge = dynamic_cast<ConstrainingEdge*>(e);
872 : UG_ASSERT(cge, "Bad type!");
873 : UG_ASSERT(cge->constrained_vertex(0), "A constrained vertex has to exist");
874 : UG_ASSERT(mg.get_child_vertex(cge) == cge->constrained_vertex(0),
875 : "Mismatch between vertex child and constrained vertex.");
876 : UG_ASSERT(cge->num_constrained_edges() == 2, "2 constrained edges have to exist!");
877 :
878 0 : for(size_t i1 = 0; i1 < cge->num_constrained_edges(); ++i1){
879 : bool constrainedEdgeMatch = false;
880 0 : for(size_t i2 = 0; i2 < mg.num_children<Edge>(cge); ++i2){
881 0 : if(mg.get_child<Edge>(cge, i2) == cge->constrained_edge(i1)){
882 : constrainedEdgeMatch = true;
883 : break;
884 : }
885 : }
886 0 : if(!constrainedEdgeMatch){
887 0 : UG_THROW("no matching child found to constrained edge");
888 : }
889 : }
890 :
891 0 : CheckElementConsistency(mg, cge->constrained_vertex(0));
892 : }
893 :
894 : // check vertices
895 0 : for(size_t i = 0; i < e->num_vertices(); ++i){
896 0 : success &= CheckElementConsistency(mg, e->vertex(i));
897 : }
898 :
899 0 : return success;
900 : }
901 :
902 0 : bool CheckElementConsistency(MultiGrid& mg, Face* f)
903 : {
904 : bool success = true;
905 0 : UG_LOG("DEBUG: Checking face at " << GetGridObjectCenter(mg, f) << endl);
906 :
907 : UG_LOG(" face type:");
908 0 : if(f->is_constrained()){
909 : UG_LOG(" constrained");
910 : }
911 0 : else if(f->is_constraining()){
912 : UG_LOG(" constraining");
913 : }
914 : else{
915 : UG_LOG(" normal");
916 : }
917 : UG_LOG("\n");
918 :
919 : // check sides
920 : Grid::edge_traits::secure_container edges;
921 : mg.associated_elements(edges, f);
922 :
923 0 : for(size_t i = 0; i < edges.size(); ++i){
924 0 : success &= CheckElementConsistency(mg, edges[i]);
925 : }
926 :
927 0 : return success;
928 : }
929 :
930 :
931 : template <class TElem>
932 : static
933 0 : std::string ElementDebugInfo_IMPL(const Grid& grid, TElem* e)
934 : {
935 0 : std::stringstream ss;
936 0 : if(!e)
937 0 : return std::string("invalid element");
938 :
939 0 : if(e->is_constrained())
940 0 : ss << "constrained ";
941 0 : else if(e->is_constraining())
942 0 : ss << "constraining ";
943 : else
944 0 : ss << "normal ";
945 :
946 0 : switch(e->base_object_id()){
947 0 : case VERTEX: ss << "vertex "; break;
948 0 : case EDGE: ss << "edge "; break;
949 0 : case FACE: ss << "face "; break;
950 0 : case VOLUME: ss << "volume "; break;
951 0 : default: ss << "UNKNOWN ELEMENT TYPE"; return ss.str();
952 : }
953 :
954 0 : ss << "at " << GetGridObjectCenter(*const_cast<Grid*>(&grid), e) << " ";
955 :
956 0 : if(const MultiGrid* mg = dynamic_cast<const MultiGrid*>(&grid)){
957 0 : ss << "on level " << mg->get_level(e);
958 0 : ss << " with assigned parent type " << (int)mg->parent_type(e) << " ";
959 : }
960 :
961 0 : if(grid.periodic_boundary_manager()){
962 : typedef typename TElem::grid_base_object base_t;
963 : base_t* b = e;
964 0 : const PeriodicBoundaryManager* pdm = grid.periodic_boundary_manager();
965 0 : if(pdm->is_slave(b)){
966 0 : ss << "-- periodic slave ";
967 0 : if(pdm->master(b)){
968 : ss << "with master at "
969 0 : << GetGridObjectCenter(*const_cast<Grid*>(&grid), pdm->master(b))
970 0 : << " ";
971 : }
972 : }
973 0 : else if(pdm->is_master(b)){
974 : typedef typename PeriodicBoundaryManager::Group<base_t>::SlaveContainer
975 : slave_container_t;
976 0 : ss << "-- periodic master ";
977 0 : slave_container_t* slaves = pdm->slaves(b);
978 0 : if(!slaves->empty()){
979 0 : ss << "with slaves at ";
980 0 : for(typename slave_container_t::iterator i = slaves->begin();
981 0 : i != slaves->end(); ++i)
982 : {
983 0 : ss << GetGridObjectCenter(*const_cast<Grid*>(&grid), *i) << " ";
984 : }
985 : }
986 : }
987 : }
988 : return ss.str();
989 0 : }
990 :
991 0 : std::string ElementDebugInfo(const Grid& grid, GridObject* e)
992 : {
993 0 : switch(e->base_object_id()){
994 0 : case VERTEX: return ElementDebugInfo_IMPL(grid, static_cast<Vertex*>(e));
995 0 : case EDGE: return ElementDebugInfo_IMPL(grid, static_cast<Edge*>(e));
996 0 : case FACE: return ElementDebugInfo_IMPL(grid, static_cast<Face*>(e));
997 0 : case VOLUME: return ElementDebugInfo_IMPL(grid, static_cast<Volume*>(e));
998 0 : default: return string("UNKNOWN ELEMENT TYPE");
999 : }
1000 : }
1001 :
1002 : }// end of namespace
|
