Line data Source code
1 : /*
2 : * Copyright (c) 2009-2015: G-CSC, Goethe University Frankfurt
3 : * Authors: Martin Stepniewski, 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 : #ifndef __H__UG__GRID_STATISTICS__
34 : #define __H__UG__GRID_STATISTICS__
35 :
36 : #include <algorithm>
37 : #include <limits>
38 : #include <vector>
39 : #include "lib_grid/lg_base.h"
40 :
41 : #ifdef UG_PARALLEL
42 : #include "lib_grid/parallelization/distributed_grid.h"
43 : #include "pcl/pcl_process_communicator.h"
44 : #endif
45 :
46 : namespace ug
47 : {
48 :
49 : //**********************************************************************
50 : // declarations
51 : //**********************************************************************
52 :
53 : ////////////////////////////////////////////////////////////////////////
54 : // AssignTetrahedronAttributesByAspectRatio - mstepnie
55 : /// assigns tetrahedral elements of a grid to subsets respecting their aspect ratio
56 : bool AssignTetrahedronAttributesByAspectRatio(Grid& grid,
57 : SubsetHandler& shVolume,
58 : AInt& aTetrahedronAspectRatioClass,
59 : std::vector<double>& offsets,
60 : Grid::VertexAttachmentAccessor<APosition>& aaPos);
61 :
62 : ////////////////////////////////////////////////////////////////////////
63 : /// assigns a subset based on the quality of the given element.
64 : /**
65 : * Currently only faces are supported.
66 : *
67 : * \param intervals contains the intervals which define into which subset
68 : * an element goes. Numbers have to be sorted, starting at
69 : * 0 and ending at 1 (0 and 1 should be contained in intervals).
70 : */
71 : template <class TIterator>
72 : bool AssignSubsetsByQuality(Grid& grid, SubsetHandler& sh,
73 : TIterator elemsBegin, TIterator elemsEnd,
74 : std::vector<number> intervals)
75 : {
76 : if(intervals.empty()){
77 : sh.assign_subset(elemsBegin, elemsEnd, 0);
78 : return true;
79 : }
80 :
81 : // access position
82 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPosition);
83 :
84 : // iterate over all elements
85 : for(TIterator iter = elemsBegin; iter != elemsEnd; ++iter)
86 : {
87 : typename TIterator::value_type elem = *iter;
88 : number quality = FaceQuality(elem, aaPos);
89 :
90 : size_t newInd = -1;
91 : for(size_t i = 0; i < intervals.size(); ++i){
92 : if(intervals[i] < quality)
93 : newInd = i;
94 : else
95 : break;
96 : }
97 :
98 : sh.assign_subset(elem, newInd);
99 : }
100 : return true;
101 : }
102 :
103 :
104 : template <class TIterator, class TAAPos>
105 0 : void PrintElementEdgeRatios(Grid& grid, TIterator elemsBegin, TIterator elemsEnd,
106 : TAAPos aaPos)
107 : {
108 : using namespace std;
109 : typedef typename PtrToValueType<typename TIterator::value_type>::base_type elem_t;
110 : UG_COND_THROW(elem_t::BASE_OBJECT_ID == VERTEX || elem_t::BASE_OBJECT_ID == EDGE,
111 : "Can't evaluate anisotropy statistics for vertices or edges.");
112 :
113 : Grid::edge_traits::secure_container edges;
114 :
115 0 : number minRatio = 1;
116 0 : number maxRatio = 0;
117 : number avRatio = 0;
118 : vector<number> ratios;
119 0 : for(TIterator i_elem = elemsBegin; i_elem != elemsEnd; ++i_elem){
120 : elem_t* elem = *i_elem;
121 :
122 : #ifdef UG_PARALLEL
123 : if(grid.distributed_grid_manager()->is_ghost(elem))
124 : continue;
125 : #endif
126 :
127 : grid.associated_elements(edges, elem);
128 0 : number shortest = numeric_limits<double>::max();
129 0 : number longest = 0;
130 0 : for_each_in_vec(Edge* e, edges){
131 0 : number l = EdgeLength(e, aaPos);
132 0 : shortest = min(shortest, l);
133 0 : longest = max(longest, l);
134 : }end_for;
135 :
136 0 : number ratio = 0;
137 0 : if(longest > 0)
138 0 : ratio = shortest / longest;
139 :
140 0 : minRatio = min(minRatio, ratio);
141 0 : maxRatio = max(maxRatio, ratio);
142 0 : avRatio += ratio;
143 0 : ratios.push_back(ratio);
144 : }
145 :
146 0 : int num = (int)ratios.size();
147 : #ifdef UG_PARALLEL
148 : pcl::ProcessCommunicator com;
149 : minRatio = com.allreduce(minRatio, PCL_RO_MIN);
150 : maxRatio = com.allreduce(maxRatio, PCL_RO_MAX);
151 : num = com.allreduce(num, PCL_RO_SUM);
152 : avRatio = com.allreduce(avRatio, PCL_RO_SUM);
153 : #endif
154 :
155 0 : if(num == 0){
156 : UG_LOG("---\n");
157 : }
158 : else{
159 0 : avRatio /= (number)num;
160 0 : UG_LOG("min: " << minRatio << ", max: " << maxRatio << ", av: " << avRatio);
161 :
162 0 : if(num > 1){
163 : number sdSum = 0;
164 0 : for(size_t i = 0; i < ratios.size(); ++i)
165 0 : sdSum += sq(avRatio - ratios[i]);
166 :
167 : #ifdef UG_PARALLEL
168 : sdSum = com.allreduce(sdSum, PCL_RO_SUM);
169 : #endif
170 :
171 0 : number sd = sqrt(sdSum / ((number)num - 1));
172 : UG_LOG(", sd: " << sd);
173 : }
174 : UG_LOG(endl);
175 : }
176 0 : }
177 :
178 : }// end of namespace
179 :
180 : #endif
|
