Line data Source code
1 : /*
2 : * Copyright (c) 2016: 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 : #include "subdivision_projector.h"
34 : #include "../../algorithms/subdivision/subdivision_rules_piecewise_loop.h"
35 :
36 : using namespace std;
37 :
38 : namespace ug{
39 :
40 0 : void SubdivisionProjector::
41 : refinement_begins(const ISubGrid* psg)
42 : {
43 : // PLEASE NOTE: The implementation below is able to perform subdivision on
44 : // manifolds which consist of sides of volume-elements, even
45 : // if volume-elements are on both sides of those manifolds.
46 0 : RefinementProjector::refinement_begins(psg);
47 : const ISubGrid& sg = *psg;
48 :
49 : m_newPositions.clear();
50 :
51 : // calculate new positions of all selected vertices and store them in m_newPositions
52 0 : SubdivRules_PLoop& subdiv = SubdivRules_PLoop::inst();
53 : Grid& g = geom().grid();
54 : Grid::edge_traits::secure_container edges;
55 :
56 0 : lg_for_each_const (Vertex, vrt, sg.goc())
57 : {
58 : g.associated_elements (edges, vrt);
59 :
60 : Edge* creaseNbrs[2];
61 0 : const size_t numCreaseNbrs = nbr_crease_edges(vrt, &edges, creaseNbrs);
62 0 : if(numCreaseNbrs == 2){
63 0 : vector3 p0 = pos(GetConnectedVertex(creaseNbrs[0], vrt));
64 0 : vector3 p1 = pos(GetConnectedVertex(creaseNbrs[1], vrt));
65 : vector3 w = subdiv.ref_even_crease_weights();
66 :
67 : vector3 p;
68 0 : VecScaleAdd(p, w.x(), pos(vrt),
69 : w.y(), p0, w.z(), p1);
70 0 : m_newPositions.push_back(make_pair(vrt, p));
71 : }
72 0 : else if(numCreaseNbrs > 0){
73 0 : m_newPositions.push_back(make_pair(vrt, pos(vrt)));
74 : }
75 : else{
76 : // perform loop subdivision on even vertices
77 : size_t valence = 0;
78 : vector3 p;
79 : VecSet(p, 0);
80 :
81 : Grid::face_traits::secure_container faces;
82 0 : for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
83 : Edge* e = edges[i_edge];
84 : g.associated_elements(faces, e);
85 0 : if(concerned_nbr_faces(e, &faces)){
86 0 : VecAdd(p, p, pos(GetConnectedVertex(e, vrt)));
87 0 : ++valence;
88 : }
89 : }
90 :
91 : number centerWgt = subdiv.ref_even_inner_center_weight(valence);
92 : number nbrWgt = subdiv.ref_even_inner_nbr_weight(valence);
93 :
94 : vector3 np;
95 0 : VecScaleAdd(np, centerWgt, pos(vrt), nbrWgt, p);
96 0 : m_newPositions.push_back(make_pair(vrt, np));
97 : }
98 : } lg_end_for;
99 0 : }
100 :
101 0 : void SubdivisionProjector::
102 : refinement_ends()
103 : {
104 : // we have to adjust positions of old vertices (no multigrid) or vertex children
105 : // (multigrid) here.
106 : Grid& g = geom().grid();
107 0 : MultiGrid* pmg = dynamic_cast<MultiGrid*>(&g);
108 :
109 0 : if(pmg){
110 : MultiGrid& mg = *pmg;
111 0 : for(new_pos_vec_t::iterator i = m_newPositions.begin();
112 0 : i != m_newPositions.end(); ++i)
113 : {
114 0 : Vertex* child = mg.get_child_vertex(i->first);
115 0 : if(child)
116 0 : set_pos(child, i->second);
117 : else
118 0 : set_pos(i->first, i->second);
119 : }
120 : }
121 : else{
122 0 : for(new_pos_vec_t::iterator i = m_newPositions.begin();
123 0 : i != m_newPositions.end(); ++i)
124 : {
125 0 : set_pos(i->first, i->second);
126 : }
127 : }
128 0 : }
129 :
130 0 : number SubdivisionProjector::
131 : new_vertex(Vertex* vrt, Edge* parent)
132 : {
133 : //todo: Adjust to new helper methods.
134 :
135 0 : const SubdivRules_PLoop& subdiv = SubdivRules_PLoop::inst();
136 :
137 : // check whether the parent edge lies inside a volume geometry. If it does,
138 : // perform linear refinement.
139 : Grid& g = geom().grid();
140 :
141 0 : const bool parentIsCrease = m_cbIsCrease(parent);
142 : Face* concernedNbrs[2];
143 : size_t numConcernedFaces = 0;
144 :
145 0 : if(!parentIsCrease)
146 0 : numConcernedFaces = concerned_nbr_faces (parent, NULL, concernedNbrs);
147 :
148 0 : if(parentIsCrease || numConcernedFaces != 2){
149 : vector2 wghts = subdiv.ref_odd_crease_weights();
150 : vector3 p;
151 0 : VecScaleAdd(p, wghts.x(), pos(parent->vertex(0)),
152 0 : wghts.y(), pos(parent->vertex(1)));
153 : set_pos(vrt, p);
154 : return 1;
155 : }
156 :
157 0 : if(concernedNbrs[0]->num_vertices() == 3 && concernedNbrs[1]->num_vertices() == 3){
158 : // the 4 vertices that are important for the calculation
159 : Vertex* v[4];
160 0 : v[0] = parent->vertex(0); v[1] = parent->vertex(1);
161 0 : v[2] = GetConnectedVertex(parent, concernedNbrs[0]);
162 0 : v[3] = GetConnectedVertex(parent, concernedNbrs[1]);
163 :
164 : vector4 wghts;
165 :
166 0 : bool isCrease0 = nbr_crease_edges(v[0]) > 0;
167 0 : bool isCrease1 = nbr_crease_edges(v[1]) > 0;
168 : // if exactly one of the two is a crease vertex, special
169 : // weighting has to be performed
170 : //todo: this does not yet work for inner creases.
171 0 : if((isCrease0 && !isCrease1) || (!isCrease0 && isCrease1))
172 : {
173 : // the crease vertex has to be in v[0]
174 0 : if(isCrease1)
175 : swap(v[0], v[1]);
176 :
177 : // get the number of edges that are connected to v[0]
178 : // todo: only count edges that are on the correct side of the crease.
179 : Grid::edge_traits::secure_container edges;
180 : g.associated_elements(edges, v[0]);
181 :
182 0 : wghts = subdiv.ref_odd_inner_weights(edges.size());
183 : }
184 : else{
185 : wghts = subdiv.ref_odd_inner_weights();
186 : }
187 :
188 : vector3 p;
189 0 : VecScaleAdd(p, wghts.x(), pos(v[0]), wghts.y(), pos(v[1]),
190 0 : wghts.z(), pos(v[2]), wghts.w(), pos(v[3]));
191 : set_pos(vrt, p);
192 : }
193 : else
194 : RefinementProjector::new_vertex(vrt, parent);
195 :
196 : return 1;
197 : }
198 :
199 0 : size_t SubdivisionProjector::
200 : nbr_crease_edges (Vertex* vrt,
201 : Grid::edge_traits::secure_container* assEdges,
202 : Edge* creaseEdgesOut[2])
203 : {
204 : Grid& grid = geom().grid();
205 :
206 : Grid::edge_traits::secure_container localAssEdges;
207 0 : if(!assEdges){
208 : grid.associated_elements(localAssEdges, vrt);
209 : assEdges = &localAssEdges;
210 : }
211 :
212 : size_t creaseEdges = 0;
213 : Grid::face_traits::secure_container faces;
214 0 : if(creaseEdgesOut){
215 0 : for(size_t i = 0; i < assEdges->size(); ++i){
216 : Edge* edge = (*assEdges)[i];
217 0 : if(m_cbIsCrease(edge)){
218 0 : if(creaseEdges < 2)
219 0 : creaseEdgesOut[creaseEdges] = edge;
220 0 : ++creaseEdges;
221 : }
222 : else{
223 : grid.associated_elements(faces, edge);
224 : size_t num = 0;
225 0 : if(faces.size() != 0)
226 0 : num = concerned_nbr_faces(edge, &faces);
227 0 : if(faces.size() == 0 || ((num > 0) && (num!= 2))){
228 0 : if(creaseEdges < 2)
229 0 : creaseEdgesOut[creaseEdges] = edge;
230 0 : ++creaseEdges;
231 : }
232 : }
233 : }
234 : }
235 : else{
236 0 : for(size_t i = 0; i < assEdges->size(); ++i){
237 : Edge* edge = (*assEdges)[i];
238 0 : if(m_cbIsCrease(edge))
239 0 : ++creaseEdges;
240 : else{
241 : grid.associated_elements(faces, edge);
242 0 : if(faces.size() == 0)
243 0 : ++creaseEdges;
244 : else{
245 0 : size_t num = concerned_nbr_faces(edge, &faces);
246 0 : creaseEdges += int((num > 0) && (num != 2));
247 : }
248 : }
249 : }
250 : }
251 :
252 0 : return creaseEdges;
253 : }
254 :
255 0 : size_t SubdivisionProjector::
256 : concerned_nbr_faces (Edge* edge,
257 : Grid::face_traits::secure_container* assFaces,
258 : Face* facesOut[2])
259 : {
260 : Grid::face_traits::secure_container localAssFaces;
261 0 : if(!assFaces){
262 : Grid& grid = geom().grid();
263 : grid.associated_elements(localAssFaces, edge);
264 : assFaces = &localAssFaces;
265 : }
266 :
267 : size_t numConcernedFaces = 0;
268 0 : if(facesOut){
269 0 : for(size_t i = 0; i < assFaces->size(); ++i){
270 0 : if(is_concerned((*assFaces)[i])){
271 0 : if(numConcernedFaces < 2)
272 0 : facesOut[numConcernedFaces] = (*assFaces)[i];
273 0 : ++numConcernedFaces;
274 : }
275 : }
276 : }
277 : else{
278 0 : for(size_t i = 0; i < assFaces->size(); ++i)
279 0 : numConcernedFaces += (size_t)is_concerned((*assFaces)[i]);
280 : }
281 :
282 0 : return numConcernedFaces;
283 : }
284 :
285 : }// end of namespace
|
