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__VERTEX_UTIL_IMPL__
34 : #define __H__LIB_GRID__VERTEX_UTIL_IMPL__
35 :
36 : #include "vertex_util.h"
37 : #include "edge_util.h"
38 : #include "face_util.h"
39 : #include "../trees/kd_tree_static.h"
40 : #include "misc_util.h"
41 :
42 : namespace ug
43 : {
44 :
45 : ////////////////////////////////////////////////////////////////////////
46 : template <class TAAPos>
47 : number VertexDistanceSq(Vertex* v0, Vertex* v1, TAAPos& aaPos)
48 : {
49 : return VecDistanceSq(aaPos[v0], aaPos[v1]);
50 : }
51 :
52 : ////////////////////////////////////////////////////////////////////////
53 : template <class TAAPos>
54 0 : number VertexDistance(Vertex* v0, Vertex* v1, TAAPos& aaPos)
55 : {
56 0 : return VecDistance(aaPos[v0], aaPos[v1]);
57 : }
58 :
59 : ////////////////////////////////////////////////////////////////////////
60 : template <class TAAPosVRT>
61 0 : void CalculateVertexNormal(vector3& nOut, Grid& grid, Vertex* vrt, TAAPosVRT& aaPos)
62 : {
63 : // set all normal to zero
64 : nOut = vector3(0, 0, 0);
65 :
66 : // loop through all associated faces, calculate their normal and add them to thee normal
67 0 : Grid::AssociatedFaceIterator iterEnd = grid.associated_faces_end(vrt);
68 0 : for(Grid::AssociatedFaceIterator iter = grid.associated_faces_begin(vrt);
69 0 : iter != iterEnd; iter++)
70 : {
71 : vector3 vN;
72 0 : CalculateNormal(vN, *iter, aaPos);
73 : VecAdd(nOut, nOut, vN);
74 : }
75 :
76 0 : VecNormalize(nOut, nOut);
77 0 : }
78 :
79 : ////////////////////////////////////////////////////////////////////////
80 : template <class TAAPosVRT>
81 0 : void CalculateBoundaryVertexNormal2D(typename TAAPosVRT::ValueType& nOut,
82 : Grid& grid, Vertex* vrt,
83 : TAAPosVRT& aaPos)
84 : {
85 : // The algorithm is a little cumbersome. However, through this setup, we
86 : // make sure that the orientation of the normal indeed points outwards,
87 : // based only on the topology.
88 :
89 : // set nOut to 0
90 : VecSet(nOut, 0);
91 :
92 : // iterate over associated faces
93 : std::vector<Face*> faces;
94 : CollectAssociated(faces, grid, vrt);
95 :
96 0 : EdgeDescriptor ed;
97 0 : for(size_t i_face = 0; i_face < faces.size(); ++i_face){
98 0 : Face* f = faces[i_face];
99 : // check for each side of f whether it is a boundary edge
100 0 : for(size_t i_side = 0; i_side < f->num_sides(); ++i_side){
101 0 : if(IsBoundaryEdge2D(grid, grid.get_edge(f, i_side))){
102 0 : f->edge_desc(i_side, ed);
103 : // make sure that e contains the specified vertex
104 0 : if(!EdgeContains(&ed, vrt))
105 0 : continue;
106 :
107 : // the normal pointing outwards is clearly defined from the
108 : // orientation of the edge descriptor
109 0 : nOut.x() += aaPos[ed.vertex(1)].y() - aaPos[ed.vertex(0)].y();
110 0 : nOut.y() += aaPos[ed.vertex(0)].x() - aaPos[ed.vertex(1)].x();
111 : }
112 : }
113 : }
114 :
115 0 : VecNormalize(nOut, nOut);
116 0 : }
117 :
118 : ////////////////////////////////////////////////////////////////////////
119 : template <class TAAPosVRT>
120 0 : void CalculateBoundaryVertexNormal3D(vector3& nOut, Grid& grid, Vertex* vrt,
121 : TAAPosVRT& aaPos)
122 : {
123 : // The algorithm is a little cumbersome. However, through this setup, we
124 : // make sure that the orientation of the normal indeed points outwards,
125 : // based only on the topology.
126 :
127 : // set nOut to 0
128 : VecSet(nOut, 0);
129 :
130 : // iterate over associated volumes
131 : std::vector<Volume*> vols;
132 : CollectAssociated(vols, grid, vrt);
133 :
134 0 : FaceDescriptor fd;
135 0 : for(size_t i_vol = 0; i_vol < vols.size(); ++i_vol){
136 0 : Volume* v = vols[i_vol];
137 : // check for each side of f whether it is a boundary edge
138 0 : for(size_t i_side = 0; i_side < v->num_sides(); ++i_side){
139 0 : if(IsBoundaryFace3D(grid, grid.get_face(v, i_side))){
140 0 : v->face_desc(i_side, fd);
141 :
142 : // make sure that fd contains the given vertex
143 0 : if(!FaceContains(&fd, vrt))
144 0 : continue;
145 :
146 : // the normal pointing outwards is clearly defined from the
147 : // orientation of the face descriptor
148 : vector3 n;
149 0 : CalculateNormal(n, &fd, aaPos);
150 : VecAdd(nOut, nOut, n);
151 : }
152 : }
153 : }
154 :
155 0 : VecNormalize(nOut, nOut);
156 0 : }
157 :
158 : ////////////////////////////////////////////////////////////////////////
159 : template <class TVrtIter, class TAPosition>
160 : void
161 : CalculateBoundingBox(typename TAPosition::ValueType& vMinOut,
162 : typename TAPosition::ValueType& vMaxOut,
163 : TVrtIter vrtsBegin, TVrtIter vrtsEnd,
164 : Grid::AttachmentAccessor<Vertex, TAPosition>& aaPos)
165 : {
166 : size_t dim = TAPosition::ValueType::Size;
167 :
168 : if(vrtsBegin != vrtsEnd)
169 : {
170 : vMinOut = aaPos[*vrtsBegin];
171 : vMaxOut = vMinOut;
172 :
173 : for(TVrtIter iter = vrtsBegin; iter != vrtsEnd; ++iter)
174 : {
175 : for(size_t i = 0; i < dim; ++i){
176 : vMinOut[i] = std::min(vMinOut[i], aaPos[*iter][i]);
177 : vMaxOut[i] = std::max(vMaxOut[i], aaPos[*iter][i]);
178 : }
179 : }
180 : }
181 : }
182 :
183 : ////////////////////////////////////////////////////////////////////////
184 : template<class TVertexPositionAttachmentAccessor>
185 : inline
186 : typename TVertexPositionAttachmentAccessor::ValueType
187 : CalculateCenter(const Vertex* v, TVertexPositionAttachmentAccessor& aaPosVRT)
188 : {
189 : return aaPosVRT[v];
190 : }
191 :
192 : template<class TAAPosVRT, class TAAWeightVRT>
193 : UG_API
194 : typename TAAPosVRT::ValueType
195 : CalculateCenter(const Vertex* v, TAAPosVRT& aaPosVRT, TAAWeightVRT&)
196 : {
197 : return aaPosVRT[v];
198 : }
199 :
200 : ////////////////////////////////////////////////////////////////////////
201 : template <class TVrtIter, class TAPosition>
202 : typename TAPosition::ValueType
203 : CalculateCenter(TVrtIter vrtsBegin, TVrtIter vrtsEnd,
204 : Grid::AttachmentAccessor<Vertex, TAPosition>& aaPos)
205 : {
206 : typename TAPosition::ValueType vMin, vMax;
207 : CalculateBoundingBox(vMin, vMax, vrtsBegin, vrtsEnd, aaPos);
208 : typename TAPosition::ValueType vRet;
209 : VecScaleAdd(vRet, 0.5, vMin, 0.5, vMax);
210 : return vRet;
211 : }
212 :
213 : ////////////////////////////////////////////////////////////////////////
214 : template <class TVrtIter, class TAPosition>
215 : typename TAPosition::ValueType
216 : CalculateBarycenter(TVrtIter vrtsBegin, TVrtIter vrtsEnd,
217 : Grid::VertexAttachmentAccessor<TAPosition>& aaPos)
218 : {
219 : typename TAPosition::ValueType v;
220 : VecSet(v, 0);
221 : int counter = 0;
222 : for(TVrtIter iter = vrtsBegin; iter != vrtsEnd; ++iter)
223 : {
224 : VecAdd(v,v,aaPos[*iter]);
225 : counter++;
226 : }
227 :
228 : if(counter>0)
229 : VecScale(v,v,1.f/counter);
230 : return v;
231 : }
232 :
233 : ////////////////////////////////////////////////////////////////////////
234 : template <class TVrtIterator>
235 : Vertex* MergeMultipleVertices(Grid& grid, TVrtIterator vrtsBegin,
236 : TVrtIterator vrtsEnd)
237 : {
238 : if(vrtsBegin == vrtsEnd)
239 : return NULL;
240 :
241 : Vertex* v = *vrtsBegin;
242 : ++vrtsBegin;
243 : while(vrtsBegin != vrtsEnd){
244 : Vertex* v2 = *vrtsBegin;
245 : ++vrtsBegin;
246 : MergeVertices(grid, v, v2);
247 : }
248 : return v;
249 : }
250 :
251 : ////////////////////////////////////////////////////////////////////////
252 : //TODO: replace KDTreeStatic by a dynamic kd-tree.
253 : //TODO: Better support for various iterators.
254 : template <int dim, class TVrtIterator>
255 0 : void RemoveDoubles(Grid& grid, const TVrtIterator& iterBegin,
256 : const TVrtIterator& iterEnd, Attachment<MathVector<dim> >& aPos,
257 : number threshold)
258 : {
259 0 : if(!grid.has_vertex_attachment(aPos))
260 : return;
261 :
262 : typedef Attachment<MathVector<dim> > attachment_type;
263 : Grid::VertexAttachmentAccessor<attachment_type> aaPos(grid, aPos);
264 :
265 0 : RemoveDoubles<dim>(grid, iterBegin, iterEnd, aaPos, threshold);
266 : }
267 :
268 : template <int dim, class TVrtIterator, class TAAPos>
269 0 : void RemoveDoubles(Grid& grid, const TVrtIterator& iterBegin,
270 : const TVrtIterator& iterEnd,
271 : TAAPos aaPos,
272 : number threshold)
273 : {
274 : typedef Attachment<MathVector<dim> > attachment_type;
275 : KDTreeStatic<attachment_type, dim, MathVector<dim> > kdTree;
276 0 : kdTree.create_from_grid(grid, iterBegin, iterEnd, aaPos, 20, 10, KDSD_LARGEST);
277 :
278 : // we need temporary attachments:
279 : // a vector<Vertex*> attachment, that stores for each vertex all other vertices
280 : // closer than threshold, which have higher attachment data index.
281 : typedef Attachment<std::list<Vertex*> > AVertexList;
282 : AVertexList aVertexList;
283 : grid.attach_to_vertices(aVertexList);
284 : Grid::VertexAttachmentAccessor<AVertexList> aaVL(grid, aVertexList);
285 :
286 : // we'll store in this attachment whether a vertex will be merged or not.
287 : AInt aInt;
288 : grid.attach_to_vertices(aInt);
289 : Grid::VertexAttachmentAccessor<AInt> aaInt(grid, aInt);
290 : {
291 0 : for(TVrtIterator iter = iterBegin; iter != iterEnd; ++iter)
292 0 : aaInt[*iter] = 0;
293 : }
294 :
295 : // compare squares.
296 0 : threshold *= threshold;
297 : std::vector<Vertex*> neighbours;
298 : // iterate over all vertices and collect all that have aInt == 0 and are within range.
299 0 : for(TVrtIterator iter = iterBegin; iter != iterEnd; ++iter)
300 : {
301 : Vertex* v = *iter;
302 0 : if(aaInt[v] == 0)
303 : {// the vertex will not be removed during merge
304 : // find all vertices closer than threshold
305 : uint numClosest = 3;
306 0 : while(numClosest < grid.num_vertices())
307 : {
308 : neighbours.clear();
309 0 : kdTree.get_neighbourhood(neighbours, aaPos[v], numClosest);
310 :
311 0 : if(VecDistanceSq(aaPos[neighbours.back()], aaPos[v]) < threshold)
312 0 : numClosest *= 2;
313 : else
314 : break;
315 : }
316 :
317 : // store them in the vertexVec attachment
318 0 : if(!neighbours.empty())
319 : {
320 : for(std::vector<Vertex*>::iterator nIter = neighbours.begin();
321 0 : nIter != neighbours.end(); ++nIter)
322 : {
323 0 : Vertex* nv = *nIter;
324 0 : if(aaInt[nv] == 0)
325 : {
326 0 : if(nv != v)
327 : {
328 0 : if(VecDistanceSq(aaPos[v], aaPos[nv]) < threshold)
329 : {
330 : aaVL[v].push_back(nv);
331 0 : aaInt[nv] = 1;
332 : }
333 : else
334 : break;
335 : }
336 : }
337 : }
338 : }
339 : }
340 : }
341 :
342 : // iterate over all vertices again and merge collected ones
343 : // This iteration only works, if the iterators stem from a class
344 : // like Selector or SubsetHandler or Grid etc, which automatically
345 : // handle removed elements.
346 : //TODO: This should be improved somehow!
347 : {
348 : TVrtIterator iter = iterBegin;
349 0 : while(iter != iterEnd)
350 : {
351 : Vertex* v = *iter;
352 0 : if(!aaVL[v].empty())
353 : {
354 : std::list<Vertex*>::iterator nIter = aaVL[v].begin();
355 0 : while(nIter != aaVL[v].end())
356 : {
357 0 : Vertex* delVrt = *nIter;
358 : nIter++;
359 0 : MergeVertices(grid, v, delVrt);
360 : }
361 : }
362 :
363 : ++iter;
364 : }
365 : }
366 :
367 : grid.detach_from_vertices(aVertexList);
368 : grid.detach_from_vertices(aInt);
369 0 : }
370 :
371 :
372 : ////////////////////////////////////////////////////////////////////////
373 : template<class TAAPos> inline
374 : void TransformVertex(Vertex* vrt, matrix33& m, TAAPos& aaPos)
375 : {
376 : // todo: avoid unnecessary copy
377 : vector3 oldPos = aaPos[vrt];
378 : MatVecMult(aaPos[vrt], m, oldPos);
379 : }
380 :
381 : ////////////////////////////////////////////////////////////////////////
382 : template<class TIterator, class TAAPos>
383 : void TransformVertices(TIterator vrtsBegin, TIterator vrtsEnd,
384 : matrix33& m, TAAPos& aaPos)
385 : {
386 : for(TIterator iter = vrtsBegin; iter != vrtsEnd; ++iter)
387 : TransformVertex(*iter, m, aaPos);
388 : }
389 :
390 : ////////////////////////////////////////////////////////////////////////
391 : template<class TIterator, class TAAPos> inline
392 : void MoveVertices(TIterator vrtsBegin, TIterator vrtsEnd, TAAPos aaPos,
393 : const typename TAAPos::ValueType& offset)
394 : {
395 : for(TIterator iter = vrtsBegin; iter != vrtsEnd; ++iter)
396 : aaPos[*iter] += offset;
397 : }
398 :
399 : ////////////////////////////////////////////////////////////////////////
400 : template <class vector_t, class TAAPos>
401 : UG_API bool
402 : ContainsPoint(const Vertex* v, const vector_t& p, TAAPos aaPos)
403 : {
404 : const vector_t& pv = aaPos[v];
405 0 : for(size_t i = 0; i < vector_t::Size; ++i){
406 0 : if(pv[i] != p[i])
407 : return false;
408 : }
409 : return true;
410 : }
411 :
412 : template <size_t dim>
413 : int FindVertexByCoordinate(const MathVector<dim>& coord, size_t ncoords, const MathVector<dim> vCoords[])
414 : {
415 :
416 : if (ncoords <= 0) return -1;
417 :
418 : size_t bestVrt = 0;
419 : number bestDistSq = VecDistanceSq(coord, vCoords[0]);
420 :
421 : for (size_t i=1; i<ncoords; ++i)
422 : {
423 : number distSq = VecDistance(coord, vCoords[i]);
424 : if(distSq < bestDistSq)
425 : {
426 : bestDistSq = distSq;
427 : bestVrt = i;
428 : }
429 : }
430 : return bestVrt;
431 : }
432 :
433 : }// end of namespace
434 :
435 : #endif
|
