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 : #include <queue>
34 : #include <stack>
35 : #include "lib_grid/lg_base.h"
36 : #include "../geom_obj_util/geom_obj_util.h"
37 : #include "lib_grid/refinement/regular_refinement.h"
38 : #include "lib_grid/refinement/projectors/cylinder_cut_projector.h"
39 : #include "grid_adaption.h"
40 :
41 : using namespace std;
42 :
43 : namespace ug
44 : {
45 :
46 0 : bool AdaptSurfaceGridToCylinder(Selector& selOut, Grid& grid,
47 : Vertex* vrtCenter, const vector3& normal,
48 : number radius, number rimSnapThreshold, AInt& aInt,
49 : APosition& aPos)
50 : {
51 0 : if(!grid.has_vertex_attachment(aPos)){
52 0 : UG_THROW("Position attachment required!");
53 : }
54 :
55 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPos);
56 :
57 0 : if(rimSnapThreshold < 0)
58 : rimSnapThreshold = 0;
59 :
60 0 : if(rimSnapThreshold > (radius - SMALL))
61 : rimSnapThreshold = radius - SMALL;
62 :
63 0 : const number smallRadius = radius - rimSnapThreshold;
64 0 : const number smallRadiusSq = smallRadius * smallRadius;
65 0 : const number largeRadius = radius + rimSnapThreshold;
66 0 : const number largeRadiusSq = largeRadius * largeRadius;
67 :
68 : // the cylinder geometry
69 : vector3 axis;
70 0 : VecNormalize(axis, normal);
71 0 : vector3 center = aaPos[vrtCenter];
72 :
73 : // recursively select all vertices in the cylinder which can be reached from a
74 : // selected vertex by following an edge. Start with the given one.
75 : // We'll also select edges which connect inner with outer vertices. Note that
76 : // some vertices are considered rim-vertices (those with a distance between
77 : // smallRadius and largeRadius). Those are neither considered inner nor outer.
78 : Selector& sel = selOut;
79 0 : sel.clear();
80 0 : sel.select(vrtCenter);
81 :
82 : stack<Vertex*> vrtStack;
83 : vrtStack.push(vrtCenter);
84 :
85 : Grid::edge_traits::secure_container edges;
86 : Grid::face_traits::secure_container faces;
87 : vector<Quadrilateral*> quads;
88 :
89 0 : while(!vrtStack.empty()){
90 0 : Vertex* curVrt = vrtStack.top();
91 : vrtStack.pop();
92 :
93 : // we have to convert associated quadrilaterals to triangles.
94 : // Be careful not to alter the array of associated elements while we iterate
95 : // over it...
96 : quads.clear();
97 : grid.associated_elements(faces, curVrt);
98 0 : for(size_t i = 0; i < faces.size(); ++i){
99 0 : if(faces[i]->num_vertices() == 4){
100 0 : Quadrilateral* q = dynamic_cast<Quadrilateral*>(faces[i]);
101 0 : if(q)
102 0 : quads.push_back(q);
103 : }
104 : }
105 :
106 0 : for(size_t i = 0; i < quads.size(); ++i){
107 0 : Triangulate(grid, quads[i], &aaPos);
108 : }
109 :
110 : // now check whether edges leave the cylinder and mark them accordingly.
111 : // Perform projection of vertices to the cylinder rim for vertices which
112 : // lie in the threshold area.
113 : grid.associated_elements(edges, curVrt);
114 :
115 0 : for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
116 : Edge* e = edges[i_edge];
117 0 : Vertex* vrt = GetConnectedVertex(e, curVrt);
118 :
119 0 : if(sel.is_selected(vrt))
120 0 : continue;
121 :
122 : vector3 p = aaPos[vrt];
123 : vector3 proj;
124 0 : ProjectPointToRay(proj, p, center, axis);
125 0 : number distSq = VecDistanceSq(p, proj);
126 :
127 0 : if(distSq < smallRadiusSq){
128 : sel.select(vrt);
129 : vrtStack.push(vrt);
130 : }
131 0 : else if(distSq < largeRadiusSq){
132 : sel.select(vrt);
133 : // cut the ray from center through p with the cylinder hull to calculate
134 : // the new position of vrt.
135 : vector3 dir;
136 : VecSubtract(dir, p, center);
137 : number t0, t1;
138 0 : if(RayCylinderIntersection(t0, t1, center, dir, center, axis, radius))
139 0 : VecScaleAdd(aaPos[vrt], 1., center, t1, dir);
140 : }
141 : else{
142 : // the edge will be refined later on
143 : sel.select(e);
144 : }
145 : }
146 : }
147 :
148 : // refine selected edges and use a special refinement callback, which places
149 : // new vertices on edges which intersect a cylinder on the cylinders hull.
150 0 : CylinderCutProjector refCallback(MakeGeometry3d(grid, aPos),
151 0 : center, axis, radius);
152 0 : Refine(grid, sel, aInt, &refCallback);
153 :
154 : // finally select all triangles which lie in the cylinder
155 0 : sel.clear();
156 : vrtStack.push(vrtCenter);
157 : sel.select(vrtCenter);
158 :
159 0 : while(!vrtStack.empty()){
160 0 : Vertex* curVrt = vrtStack.top();
161 : vrtStack.pop();
162 : grid.associated_elements(faces, curVrt);
163 :
164 0 : for(size_t i_face = 0; i_face < faces.size(); ++i_face){
165 : Face* f = faces[i_face];
166 0 : if(sel.is_selected(f))
167 0 : continue;
168 :
169 : sel.select(f);
170 :
171 0 : for(size_t i = 0; i < f->num_vertices(); ++i){
172 0 : Vertex* vrt = f->vertex(i);
173 0 : if(!sel.is_selected(vrt)){
174 0 : number dist = DistancePointToRay(aaPos[vrt], center, axis);
175 0 : if(dist < (radius - SMALL)){
176 : sel.select(vrt);
177 : vrtStack.push(vrt);
178 : }
179 : }
180 : }
181 : }
182 : }
183 :
184 0 : sel.clear<Vertex>();
185 :
186 0 : return true;
187 0 : }
188 :
189 : ////////////////////////////////////////////////////////////////////////
190 0 : bool AdaptSurfaceGridToCylinder(Selector& selOut, Grid& grid,
191 : Vertex* vrtCenter, const vector3& normal,
192 : number radius, number rimSnapThreshold, APosition& aPos)
193 : {
194 : AInt aInt;
195 : grid.attach_to_edges(aInt);
196 0 : bool retVal = AdaptSurfaceGridToCylinder(selOut, grid, vrtCenter, normal,
197 : radius, rimSnapThreshold, aInt, aPos);
198 : grid.detach_from_edges(aInt);
199 0 : return retVal;
200 : }
201 :
202 : }// end of namespace
|
