Line data Source code
1 : /*
2 : * Copyright (c) 2010-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 : #ifndef __H__LIB_GRID__VOLUME_UTIL_IMPL__
34 : #define __H__LIB_GRID__VOLUME_UTIL_IMPL__
35 :
36 : #include <algorithm>
37 : #include <vector>
38 : #include "lib_grid/lg_base.h"
39 : #include "common/static_assert.h"
40 : #include "common/util/vec_for_each.h"
41 : #include "lib_grid/grid_objects/hexahedron_rules.h"
42 : #include "lib_grid/grid_objects/prism_rules.h"
43 : #include "lib_grid/grid_objects/pyramid_rules.h"
44 :
45 : namespace ug
46 : {
47 : ////////////////////////////////////////////////////////////////////////
48 : template <class TAAPos>
49 : bool
50 0 : ContainsPoint(Volume* vol, const vector3& p, TAAPos aaPos)
51 : {
52 : // iterate over face descriptors of the sides and check whether the point
53 : // lies inside or outside
54 0 : FaceDescriptor fd;
55 : vector3 n, dir;
56 :
57 : // to minimize rounding errors we'll compare against a small-constant which is
58 : // relative to the first edge of the examined volume.
59 0 : EdgeDescriptor ed;
60 0 : vol->edge_desc(0, ed);
61 0 : number len = EdgeLength(&ed, aaPos);
62 :
63 : // the constant should be relative to the same geometric measure as what it is
64 : // compared against later on, i.e. length*area, since otherwise problems arise
65 : // with geometries scaled to very small extensions;
66 : // which is why I changed sqrt(lenSq) to lenSq^1.5 (mbreit, 2015-05-11)
67 0 : const number locSmall = len * len * len * SMALL;
68 :
69 0 : for(size_t i = 0; i < vol->num_faces(); ++i){
70 0 : vol->face_desc(i, fd);
71 0 : CalculateNormalNoNormalize(n, &fd, aaPos);
72 : VecSubtract(dir, aaPos[fd.vertex(0)], p);
73 :
74 0 : if(VecDot(dir, n) < -locSmall)
75 : return false;
76 : }
77 : return true;
78 : }
79 :
80 : ////////////////////////////////////////////////////////////////////////
81 : // PointIsInsideTetrahedron
82 : inline bool
83 0 : PointIsInsideTetrahedron(const vector3& v, Tetrahedron* tet,
84 : Grid::VertexAttachmentAccessor<APosition>& aaPos)
85 : {
86 0 : return PointIsInsideTetrahedron(v, aaPos[tet->vertex(0)], aaPos[tet->vertex(1)],
87 0 : aaPos[tet->vertex(2)], aaPos[tet->vertex(3)]);
88 : }
89 :
90 : ////////////////////////////////////////////////////////////////////////
91 : // CalculateCenter
92 : template<class TVertexPositionAttachmentAccessor>
93 : typename TVertexPositionAttachmentAccessor::ValueType
94 0 : CalculateCenter(const VolumeVertices* vol, TVertexPositionAttachmentAccessor& aaPosVRT)
95 : {
96 : typename TVertexPositionAttachmentAccessor::ValueType v;
97 : // init v with 0.
98 : VecSet(v, 0);
99 :
100 0 : size_t numVrts = vol->num_vertices();
101 0 : VolumeVertices::ConstVertexArray vrts = vol->vertices();
102 :
103 : // sum up
104 0 : for(size_t i = 0; i < numVrts; ++i)
105 : {
106 0 : VecAdd(v, v, aaPosVRT[vrts[i]]);
107 : }
108 :
109 : // average
110 0 : if(numVrts > 0)
111 0 : VecScale(v, v, 1./(number)numVrts);
112 :
113 0 : return v;
114 : }
115 :
116 : ////////////////////////////////////////////////////////////////////////
117 : template<class TAAPosVRT, class TAAWeightVRT>
118 : UG_API
119 : typename TAAPosVRT::ValueType
120 : CalculateCenter(const VolumeVertices* vol, TAAPosVRT& aaPos, TAAWeightVRT& aaWeight)
121 : {
122 : typename TAAPosVRT::ValueType v;
123 : typedef typename TAAWeightVRT::ValueType weight_t;
124 :
125 : // init v with 0.
126 : VecSet(v, 0);
127 :
128 : size_t numVrts = vol->num_vertices();
129 : VolumeVertices::ConstVertexArray vrts = vol->vertices();
130 :
131 : // sum up
132 : weight_t totalWeight = 0;
133 : for(size_t i = 0; i < numVrts; ++i)
134 : {
135 : weight_t w = aaWeight[vrts[i]];
136 : VecScaleAppend(v, w, aaPos[vrts[i]]);
137 : totalWeight += w;
138 : }
139 :
140 : // average
141 : if(totalWeight != 0)
142 : VecScale(v, v, 1./(number)totalWeight);
143 :
144 : return v;
145 : }
146 :
147 :
148 : /// Can be used to compare vertices of their grids through their hash-value.
149 : template <class TElem>
150 : class CmpVrtsByHash{
151 : public:
152 : CmpVrtsByHash(TElem* e) : m_e(e) {};
153 : bool operator () (int i0, int i1) {
154 : return m_e->vertex(i0)->get_hash_value() < m_e->vertex(i1)->get_hash_value();
155 : }
156 : private:
157 : TElem* m_e;
158 : };
159 :
160 :
161 : template <class TVolIter>
162 : void ConvertToTetrahedra (
163 : Grid& grid,
164 : TVolIter volsBegin,
165 : TVolIter volsEnd)
166 : {
167 : using namespace std;
168 :
169 : // first we'll collect all quadrilaterals that are connected to selected
170 : // volumes. Avoid 'grid.mark' here, since it may be used by 'grid.associated_elements'
171 : Grid::face_traits::secure_container faces;
172 : vector<Face*> quads;
173 :
174 : for(TVolIter iv = volsBegin; iv != volsEnd; ++iv){
175 : grid.associated_elements(faces, *iv);
176 : for_each_in_vec(Face* f, faces){
177 : if(f->num_vertices() == 4)
178 : quads.push_back(f);
179 : }end_for;
180 : }
181 :
182 : // remove double entries
183 : grid.begin_marking();
184 : size_t offset = 0;
185 :
186 : for(size_t i = 0; i + offset < quads.size();){
187 : if(offset > 0)
188 : quads[i] = quads[i + offset];
189 :
190 : if(!grid.is_marked(quads[i])){
191 : grid.mark(quads[i]);
192 : ++i;
193 : }
194 : else{
195 : ++offset;
196 : }
197 : }
198 :
199 : if(offset > 0)
200 : quads.resize(quads.size() - offset);
201 :
202 : grid.end_marking();
203 :
204 : for_each_in_vec(Face* f, quads){
205 : //todo in a parallel environment, global id's should be compared here
206 : CmpVrtsByHash<Face> cmp(f);
207 : // get the smallest vertex of the face
208 : int smallest = 0;
209 : for(int i = 1; i < 4; ++i){
210 : if(cmp(i, smallest))
211 : smallest = i;
212 : }
213 :
214 : int i0 = smallest;
215 : int i1 = (smallest + 1) % 4;
216 : int i2 = (smallest + 2) % 4;
217 : int i3 = (smallest + 3) % 4;
218 : grid.create<Triangle>(TriangleDescriptor(f->vertex(i0), f->vertex(i1), f->vertex(i2)), f);
219 : grid.create<Triangle>(TriangleDescriptor(f->vertex(i2), f->vertex(i3), f->vertex(i0)), f);
220 : }end_for;
221 :
222 :
223 : // now convert the given volume-elements
224 : UG_STATIC_ASSERT((prism_rules::MAX_NUM_CONVERT_TO_TETS_INDS_OUT <
225 : hex_rules::MAX_NUM_CONVERT_TO_TETS_INDS_OUT) &&
226 : (pyra_rules::MAX_NUM_CONVERT_TO_TETS_INDS_OUT <
227 : hex_rules::MAX_NUM_CONVERT_TO_TETS_INDS_OUT),
228 : HEX_RULES_MAX_NUM_CONVERT_TO_TETS_INDS_OUT__considered_to_be_highest_among_prism_pyra_and_hex);
229 : static const int arrayLen = hex_rules::MAX_NUM_CONVERT_TO_TETS_INDS_OUT;
230 : int inds[arrayLen];
231 :
232 : vector<Volume*> volsToErase;
233 : for(TVolIter iv = volsBegin; iv != volsEnd; ++iv){
234 : Volume* vol = *iv;
235 : const ReferenceObjectID roid = vol->reference_object_id();
236 : CmpVrtsByHash<Volume> cmp(vol);
237 : size_t numEntries = 0;
238 :
239 : switch(roid){
240 : case ROID_PYRAMID:
241 : numEntries = pyra_rules::ConvertToTetrahedra(inds, cmp);
242 : break;
243 :
244 : case ROID_PRISM:
245 : numEntries = prism_rules::ConvertToTetrahedra(inds, cmp);
246 : break;
247 :
248 : case ROID_HEXAHEDRON:
249 : UG_THROW("ConvertToTetrahedra for hexahedra not yet implemented!");
250 : break;
251 : default:
252 : break;
253 : }
254 :
255 : if(numEntries > 0){
256 : volsToErase.push_back(vol);
257 : size_t i = 0;
258 : Volume::ConstVertexArray vrts = vol->vertices();
259 : while(i < numEntries){
260 : int goid = inds[i++];
261 : UG_COND_THROW(goid != GOID_TETRAHEDRON,
262 : "Only tetrahedra may result from ConvertToTetrahedra");
263 : int i0 = inds[i++];
264 : int i1 = inds[i++];
265 : int i2 = inds[i++];
266 : int i3 = inds[i++];
267 : grid.create<Tetrahedron>(
268 : TetrahedronDescriptor(vrts[i0], vrts[i1], vrts[i2], vrts[i3]),
269 : vol);
270 : }
271 : }
272 : }
273 :
274 : // finally erase all unused volumes and quadrilaterals
275 : for_each_in_vec(Volume* v, volsToErase){
276 : grid.erase(v);
277 : }end_for;
278 :
279 : Grid::volume_traits::secure_container assVols;
280 : for_each_in_vec(Face* f, quads){
281 : grid.associated_elements(assVols, f);
282 : if(assVols.empty())
283 : grid.erase(f);
284 : }end_for;
285 : }
286 :
287 : }// end of namespace
288 :
289 : #endif
|
