Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Authors: Sebastian Reiter, 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 : #ifndef __H__UG__SUBDIVISION_LOOP__
34 : #define __H__UG__SUBDIVISION_LOOP__
35 :
36 : #include "lib_grid/lg_base.h"
37 : #include "lib_grid/algorithms/geom_obj_util/geom_obj_util.h"
38 : #include "subdivision_rules_piecewise_loop.h"
39 :
40 : namespace ug
41 : {
42 :
43 : /// \addtogroup lib_grid_algorithms_refinement_subdivision
44 : /// @{
45 :
46 : ////////////////////////////////////////////////////////////////////////
47 : /// projects all vertices in the given grid to their limit-positions using the piecewise loop scheme.
48 : template <class TAVrtPos> void
49 0 : ProjectToLimitPLoop(Grid& grid, TAVrtPos aPos, TAVrtPos aProjPos)
50 : {
51 : // position type
52 : typedef typename TAVrtPos::ValueType pos_type;
53 :
54 : // access the subdivision weights
55 0 : SubdivRules_PLoop& subdiv = SubdivRules_PLoop::inst();
56 :
57 : // if aPos and aProjPos are equal, we'll have to create a temporary
58 : // attachment
59 : bool usingTmpAttachment = false;
60 : if(aPos == aProjPos){
61 : // create temporary attachment
62 : usingTmpAttachment = true;
63 : aProjPos = TAVrtPos();
64 : }
65 :
66 : // attach aProjPos if required
67 0 : if(!grid.has_vertex_attachment(aProjPos))
68 0 : grid.attach_to_vertices(aProjPos);
69 :
70 : // attachment accessors
71 : Grid::VertexAttachmentAccessor<TAVrtPos> aaPos(grid, aPos);
72 : Grid::VertexAttachmentAccessor<TAVrtPos> aaProjPos(grid, aProjPos);
73 :
74 : // vectors to hold temporary results
75 : typedef SubdivRules_PLoop::NeighborInfo NbrInfo;
76 : std::vector<NbrInfo> nbrInfos;
77 : std::vector<Vertex*> vrts;
78 : std::vector<number> nbrWgts;
79 :
80 : // if volumes are contained in the grid, we currently only perform projection
81 : // of boundary elements (no creases...)
82 0 : if(grid.num<Volume>() > 0){
83 : for(VertexIterator iter = grid.vertices_begin();
84 0 : iter != grid.vertices_end(); ++iter)
85 : {
86 : Vertex* vrt = *iter;
87 : // collect all surface neighbors of vrt in vrts
88 : vrts.clear();
89 0 : for(Grid::AssociatedEdgeIterator iter = grid.associated_edges_begin(vrt);
90 0 : iter != grid.associated_edges_end(vrt); ++iter)
91 : {
92 0 : if(IsBoundaryEdge3D(grid, *iter))
93 0 : vrts.push_back(GetConnectedVertex(*iter, vrt));
94 : }
95 :
96 : // default loop projection
97 : VecScale(aaProjPos[vrt], aaPos[vrt],
98 : subdiv.proj_inner_center_weight(vrts.size()));
99 :
100 0 : for(size_t i = 0; i < vrts.size(); ++i){
101 : VecScaleAdd(aaProjPos[vrt], 1.0, aaProjPos[vrt],
102 : subdiv.proj_inner_nbr_weight(vrts.size()),
103 0 : aaPos[vrts[i]]);
104 : }
105 : }
106 : }
107 : else{
108 : // iterate through all vertices
109 : for(VertexIterator iter = grid.vertices_begin();
110 0 : iter != grid.vertices_end(); ++iter)
111 : {
112 : Vertex* v = *iter;
113 :
114 : // check whether the vertex is a crease vertex or not
115 : //todo: this has to be more flexible
116 0 : if(IsBoundaryVertex2D(grid, v)){
117 : // project the crease vertex
118 : Edge* nbrs[2];
119 : size_t numNbrs = 0;
120 0 : for(Grid::AssociatedEdgeIterator iter = grid.associated_edges_begin(v);
121 0 : iter != grid.associated_edges_end(v); ++iter)
122 : {
123 0 : if(IsBoundaryEdge2D(grid, *iter)){
124 0 : nbrs[numNbrs] = *iter;
125 0 : ++numNbrs;
126 0 : if(numNbrs == 2)
127 : break;
128 : }
129 : }
130 :
131 0 : if(numNbrs == 2){
132 0 : pos_type& p0 = aaPos[GetConnectedVertex(nbrs[0], v)];
133 0 : pos_type& p1 = aaPos[GetConnectedVertex(nbrs[1], v)];
134 : vector3 w = subdiv.proj_crease_weights();
135 : VecScaleAdd(aaProjPos[v], w.x(), aaPos[v], w.y(), p0, w.z(), p1);
136 : }
137 : else
138 : aaProjPos[v] = aaPos[v];
139 : }
140 : else{
141 : // project an inner vertex
142 : // we have to calculate the valence and
143 : // we have to check whether the vertex is a neighbor to a crease.
144 : //todo: This could be given my some kind of mark.
145 : bool creaseNbr = false;
146 :
147 : // collect all neighbors in an ordered set.
148 : //todo: the order is only important if the node is indeed a crease neighbor.
149 0 : if(!CollectSurfaceNeighborsSorted(vrts, grid, v)){
150 : UG_LOG("WARNING in ProjectToLimitPLoop: surface is not regular.");
151 : UG_LOG(" Ignoring vertex...\n");
152 : aaProjPos[v] = aaPos[v];
153 0 : continue;
154 : }
155 :
156 0 : nbrInfos.resize(vrts.size());
157 0 : for(size_t i = 0; i < vrts.size(); ++i)
158 : {
159 0 : Vertex* nbrVrt = vrts[i];
160 : // we have to check whether nbrVrt is a crease edge. If it is,
161 : // we have to calculate its valence.
162 0 : if(IsBoundaryVertex2D(grid, nbrVrt)){
163 : creaseNbr = true;
164 : size_t creaseValence = 0;
165 0 : for(Grid::AssociatedEdgeIterator iter = grid.associated_edges_begin(nbrVrt);
166 0 : iter != grid.associated_edges_end(nbrVrt); ++iter)
167 : {
168 0 : ++creaseValence;
169 : }
170 :
171 0 : nbrInfos[i] = NbrInfo(nbrVrt, creaseValence);
172 : }
173 : else
174 0 : nbrInfos[i] = NbrInfo(nbrVrt, 0);
175 : }
176 :
177 : // if the vertex is a crease-nbr, we'll apply special weights.
178 : // if not, normal loop-projection-masks are used.
179 :
180 0 : if(creaseNbr){
181 : number cntrWgt;
182 0 : nbrWgts.resize(vrts.size());
183 0 : subdiv.proj_inner_crease_nbr_weights(cntrWgt, &nbrWgts.front(),
184 : &nbrInfos.front(), nbrInfos.size());
185 :
186 : // special crease neighbor projection
187 0 : VecScale(aaProjPos[v], aaPos[v], cntrWgt);
188 :
189 0 : for(size_t i = 0; i < nbrWgts.size(); ++i){
190 0 : VecScaleAdd(aaProjPos[v], 1.0, aaProjPos[v],
191 0 : nbrWgts[i], aaPos[nbrInfos[i].nbr]);
192 : }
193 : }
194 : else{
195 : // default loop projection
196 : VecScale(aaProjPos[v], aaPos[v],
197 : subdiv.proj_inner_center_weight(nbrInfos.size()));
198 :
199 0 : for(size_t i = 0; i < nbrInfos.size(); ++i){
200 : VecScaleAdd(aaProjPos[v], 1.0, aaProjPos[v],
201 : subdiv.proj_inner_nbr_weight(nbrInfos.size()),
202 0 : aaPos[nbrInfos[i].nbr]);
203 : }
204 : }
205 : }
206 : }
207 : }
208 :
209 : // clean up
210 0 : if(usingTmpAttachment){
211 : // swap attachment buffers
212 : aaPos.swap(aaProjPos);
213 :
214 : // detach it
215 0 : grid.detach_from_vertices(aProjPos);
216 : }
217 0 : }
218 :
219 : ////////////////////////////////////////////////////////////////////////
220 : /// projects all boundary vertices in the given grid to their limit-positions using the piecewise loop scheme.
221 : template <class TAVrtPos> void
222 : ProjectToLimitSubdivBoundary(Grid& grid, TAVrtPos aPos, TAVrtPos aProjPos)
223 : {
224 : // position type
225 : typedef typename TAVrtPos::ValueType pos_type;
226 :
227 : // access the subdivision weights
228 : SubdivRules_PLoop& subdiv = SubdivRules_PLoop::inst();
229 :
230 : // if aPos and aProjPos are equal, we'll have to create a temporary
231 : // attachment
232 : bool usingTmpAttachment = false;
233 : if(aPos == aProjPos){
234 : // create temporary attachment
235 : usingTmpAttachment = true;
236 : aProjPos = TAVrtPos();
237 : }
238 :
239 : // attach aProjPos if required
240 : if(!grid.has_vertex_attachment(aProjPos))
241 : grid.attach_to_vertices(aProjPos);
242 :
243 : // attachment accessors
244 : Grid::VertexAttachmentAccessor<TAVrtPos> aaPos(grid, aPos);
245 : Grid::VertexAttachmentAccessor<TAVrtPos> aaProjPos(grid, aProjPos);
246 :
247 : // iterate through all vertices
248 : for(VertexIterator iter = grid.vertices_begin();
249 : iter != grid.vertices_end(); ++iter)
250 : {
251 : Vertex* v = *iter;
252 :
253 : // check whether the vertex is a crease vertex or not
254 : //todo: this has to be more flexible
255 : if(IsBoundaryVertex2D(grid, v)){
256 : // project the crease vertex
257 : Edge* nbrs[2];
258 : size_t numNbrs = 0;
259 : for(Grid::AssociatedEdgeIterator iter = grid.associated_edges_begin(v);
260 : iter != grid.associated_edges_end(v); ++iter)
261 : {
262 : if(IsBoundaryEdge2D(grid, *iter)){
263 : nbrs[numNbrs] = *iter;
264 : ++numNbrs;
265 : if(numNbrs == 2)
266 : break;
267 : }
268 : }
269 :
270 : if(numNbrs == 2){
271 : pos_type& p0 = aaPos[GetConnectedVertex(nbrs[0], v)];
272 : pos_type& p1 = aaPos[GetConnectedVertex(nbrs[1], v)];
273 : vector3 w = subdiv.proj_crease_weights();
274 : VecScaleAdd(aaProjPos[v], w.x(), aaPos[v], w.y(), p0, w.z(), p1);
275 : }
276 : else
277 : aaProjPos[v] = aaPos[v];
278 : }
279 : else{
280 : // inner vertices are not moved
281 : aaProjPos[v] = aaPos[v];
282 : }
283 : }
284 :
285 : // clean up
286 : if(usingTmpAttachment){
287 : // swap attachment buffers
288 : aaPos.swap(aaProjPos);
289 :
290 : // detach it
291 : grid.detach_from_vertices(aProjPos);
292 : }
293 : }
294 :
295 : ////////////////////////////////////////////////////////////////////////
296 : // ProjectToLimitLoop
297 : /// projects surface vertices to their limit subdivision surface position
298 : bool ProjectToLimitLoop(Grid& grid, APosition& aProjPos);
299 :
300 : /// @} // end of add_to_group
301 :
302 : }// end of namespace
303 :
304 : #endif
|
