Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Stepniewski
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 "lib_grid/algorithms/geom_obj_util/geom_obj_util.h"
34 : #include "subdivision_loop.h"
35 :
36 : using namespace std;
37 :
38 : namespace ug
39 : {
40 : ////////////////////////////////////////////////////////////////////////
41 : // ProjectToLimitLoop
42 : /// projects surface vertices to their limit subdivision surface position
43 0 : bool ProjectToLimitLoop(Grid& grid, APosition& aProjPos)
44 : {
45 : // grid management
46 : Grid::VertexAttachmentAccessor<AVector3> aaPos(grid,aPosition);
47 : Grid::VertexAttachmentAccessor<APosition> aaProjPos(grid, aProjPos);
48 :
49 : // needed variables
50 : double const pi = 3.14159265;
51 : double x = 0;
52 : double y = 0;
53 : double z = 0;
54 : int valence = 0;
55 : const int numPrecalculated = 10;
56 : double beta[numPrecalculated];
57 : double b = 0;
58 : double chi = 0;
59 :
60 : // calculate weights for subdivision mask
61 0 : for(int i = 1; i < numPrecalculated; ++i)
62 : {
63 0 : double tmp = 0.375 + 0.25 * cos( (2.0 * pi) / (float)i );
64 0 : beta[i] = ( 0.625 - tmp * tmp ) / (float)i ;
65 : }
66 :
67 : beta[0] = 0;
68 0 : beta[6] = 0.0625;
69 :
70 : // iterate through all vertices, evaluate their limit positions and save them in their projection attachment
71 0 : for(VertexIterator vIter = grid.vertices_begin(); vIter != grid.vertices_end(); ++vIter)
72 : {
73 : Vertex* v = *vIter;
74 : valence = 0;
75 : x = 0;
76 : y = 0;
77 : z = 0;
78 :
79 0 : for(Grid::AssociatedEdgeIterator eIter = grid.associated_edges_begin(v); eIter != grid.associated_edges_end(v); ++eIter)
80 : {
81 0 : Edge* e = *eIter;
82 0 : valence++;
83 :
84 0 : if(valence >= numPrecalculated)
85 : {
86 0 : double tmp = 0.375 + 0.25 * cos( (2.0*pi) / (float)valence );
87 0 : b = (0.625 - tmp*tmp) / (float)valence;
88 : }
89 :
90 : else
91 0 : b = beta[valence];
92 :
93 0 : chi = 1.0 / (0.375 / b + valence);
94 :
95 0 : if(aaPos[v].x() == aaPos[e->vertex(0)].x() && aaPos[v].y() == aaPos[e->vertex(0)].y() && aaPos[v].z() == aaPos[e->vertex(0)].z())
96 : {
97 0 : x += aaPos[e->vertex(1)].x();
98 0 : y += aaPos[e->vertex(1)].y();
99 0 : z += aaPos[e->vertex(1)].z();
100 : }
101 :
102 : else
103 : {
104 0 : x += aaPos[e->vertex(0)].x();
105 0 : y += aaPos[e->vertex(0)].y();
106 0 : z += aaPos[e->vertex(0)].z();
107 : }
108 : }
109 :
110 0 : x*=chi;
111 0 : y*=chi;
112 0 : z*=chi;
113 :
114 0 : aaProjPos[v].x() = aaPos[v].x() * (1.0 - (float)valence * chi);
115 0 : aaProjPos[v].y() = aaPos[v].y() * (1.0 - (float)valence * chi);
116 0 : aaProjPos[v].z() = aaPos[v].z() * (1.0 - (float)valence * chi);
117 :
118 0 : aaProjPos[v].x() += x;
119 0 : aaProjPos[v].y() += y;
120 0 : aaProjPos[v].z() += z;
121 : }
122 :
123 0 : return true;
124 : }
125 :
126 : }// end of namespace
127 :
|
