dune-grid-glue 2.5-git
computeintersection.cc
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4namespace Dune {
5namespace GridGlue {
6
7//****************************************************************************************
8// PUBLIC
9//****************************************************************************************
10
11template<class CM>
13 const std::vector<V> Y,
14 std::vector<std::vector<int> >& SX,
15 std::vector<std::vector<int> >& SY,
16 std::vector<V>& P) {
17
18 std::vector<std::vector<unsigned int> > subElementsX, subElementsY;
19 std::vector<std::vector<int> > faceIdsX, faceIdsY;
20 std::vector<V> subElementX(CM::grid1Dimension+1), subElementY(CM::grid2Dimension+1), sP;
21 std::vector<std::vector<int> > sSX, sSY;
22
23 CM::grid1_subdivisions(X,subElementsX,faceIdsX);
24 CM::grid2_subdivisions(Y,subElementsY,faceIdsY);
25
26 bool intersectionFound = false;
27
28 for (unsigned int i = 0; i < subElementsX.size(); ++i) { // iterate over all X subelements
29 for (unsigned int ki = 0; ki < subElementsX[i].size(); ++ki) // define the X subelement
30 subElementX[ki] = X[subElementsX[i][ki]];
31 for (unsigned int j = 0; j < subElementsY.size(); ++j) { // iterate over all Y subelemetns
32 for (unsigned int kj = 0; kj < subElementsY[j].size(); ++kj) // define the Y subleement
33 subElementY[kj] = Y[subElementsY[j][kj]];
34
35 sP.clear();
36
37 // compute the intersection
38 bool b = CM::computeIntersectionPoints(subElementX,subElementY,sSX,sSY,sP);
39 intersectionFound = intersectionFound || b;
40
41 // only insert points on outer faces
42 for (unsigned int ki = 0; ki < sSX.size(); ++ki) { // iterate over all faces
43 if (faceIdsX[i][ki] >= 0) {
44 for (unsigned int kii = 0; kii < sSX[ki].size(); ++kii) {
45 int k = insertPoint(sP[sSX[ki][kii]],P); // determine index in P
46 SX[faceIdsX[i][ki]].push_back(k);
47 }
48 }
49 }
50 for (unsigned int kj = 0; kj < sSY.size(); ++kj) { // iterate over all faces
51 if (faceIdsY[j][kj] >= 0) {
52 for (unsigned int kjj = 0; kjj < sSY[kj].size(); ++kjj) {
53 int k = insertPoint(sP[sSY[kj][kjj]],P); // determine index in P
54 SY[faceIdsY[j][kj]].push_back(k);
55 }
56 }
57 }
58 }
59 }
60
61 return intersectionFound;
62}
63
64//****************************************************************************************
65// PRIVATE
66//****************************************************************************************
67
68template<class CM>
69void IntersectionComputation<CM>::orderPoints_(std::integral_constant<int,3>,
70 std::integral_constant<int,3>,
71 const V centroid,
72 const std::vector<std::vector<int> > SX,
73 const std::vector<std::vector<int> > SY,
74 const std::vector<V> P,
75 std::vector<std::vector<int> >& H)
76{
77 int n_facesX = SX.size();
78 int n_facesY = SY.size();
79 int m;
80
81 std::vector<int> no,id,temp ;
82 std::vector<V> p ;
83 std::vector<std::vector<int> > tempH;
84
85 std::vector<int> faceOrderingX(n_facesX);
86 std::vector<int> faceOrderingY(n_facesY);
87
88 if (n_facesX==3) {
89 faceOrderingX[0] = 0; faceOrderingX[1] = 2; faceOrderingX[2] = 1;
90 } else {
91 faceOrderingX[0] = 0; faceOrderingX[1] = 3; faceOrderingX[2] = 2; faceOrderingX[3] = 1;
92 }
93 if (n_facesY==3) {
94 faceOrderingY[0] = 0; faceOrderingY[1] = 2; faceOrderingY[2] = 1;
95 } else {
96 faceOrderingY[0] = 0; faceOrderingY[1] = 3; faceOrderingY[2] = 2; faceOrderingY[3] = 1;
97 }
98
99 if (P.size() > 3) {
100 for (int i = 0; i < n_facesX; ++i) { // loop on faces of X
101 if (SX[i].size() > 0) {
102 no = SX[faceOrderingX[i]];
103 removeDuplicates(no);
104 m = no.size() ;
105 if ((m>=3) && newFace3D(no,tempH)) // don't compute degenerate polygons and check if face is new
106 {
107 for ( int l=0; l<m; ++l)
108 p.push_back(P[no[l]]);
109 orderPointsCC(std::integral_constant<int,3>(), centroid,id,p); // order points counter-clock-wise
110 for ( int l=0; l<m; ++l)
111 temp.push_back(no[id[l]]) ;
112 tempH.push_back(temp) ;
113 temp.clear();
114 p.clear();
115 id.clear(); // clean
116 }
117 no.clear() ; // clean
118 }
119 }
120 for (int i = 0; i < n_facesY; ++i) { // loop on faces of Y
121 if (SY[i].size() > 0) {
122 no = SY[faceOrderingY[i]];
123 removeDuplicates(no);
124 m = no.size() ;
125 if ((m>=3) && newFace3D(no,tempH)) // don't compute degenerate polygons and check if face is new
126 {
127 for ( int l=0; l<m; ++l)
128 p.push_back(P[no[l]]) ;
129 orderPointsCC(std::integral_constant<int,3>(),centroid,id,p); // order points counter-clock-wise
130 for ( int l=0; l<m; ++l)
131 temp.push_back(no[id[l]]) ;
132 tempH.push_back(temp) ;
133 temp.clear();
134 p.clear();
135 id.clear(); // clean
136 }
137 no.clear() ; // clean
138 }
139 }
140 }
141
142 for (int i = 0; i < tempH.size(); ++i) {
143 int hs = tempH[i].size();
144 if (hs >= 3) {
145 for (int j = 1; j <= hs-2;++j) {
146 temp.clear();
147 temp.push_back(tempH[i][0]);
148 for (int k = 0; k < 2; ++k)
149 temp.push_back(tempH[i][j+k]);
150 H.push_back(temp);
151 }
152 }
153 }
154}
155
156template<class CM>
157void IntersectionComputation<CM>::orderPoints_(std::integral_constant<int,2>,
158 std::integral_constant<int,2>,
159 const V centroid,
160 const std::vector<std::vector<int> > SX,
161 const std::vector<std::vector<int> > SY,
162 const std::vector<V> P,
163 std::vector<std::vector<int> >& H)
164{
165 H.clear();
166 std::vector<int> id, temp(2);
167
168 orderPointsCC(std::integral_constant<int,2>(),centroid,id,P);
169
170 for (std::size_t i = 0; i < id.size();++i) {
171 temp[0] = id[i];
172 temp[1] = id[(i+1)%(id.size())];
173 H.push_back(temp);
174 }
175}
176
177template<class CM>
178void IntersectionComputation<CM>::orderPoints_(std::integral_constant<int,2>,
179 std::integral_constant<int,3>,
180 const V centroid,
181 const std::vector<std::vector<int> > SX,
182 const std::vector<std::vector<int> > SY,
183 const std::vector<V> P,
184 std::vector<std::vector<int> >& H)
185{
186 H.clear();
187 std::vector<int> id, temp(2);
188
189 orderPointsCC(std::integral_constant<int,3>(),centroid,id,P);
190
191 for (int i = 0; i < id.size();++i) {
192 temp[0] = id[i];
193 temp[1] = id[(i+1)%(id.size())];
194 H.push_back(temp);
195 }
196}
197
198template<class CM>
199void IntersectionComputation<CM>::removeDuplicates(std::vector<int> & p)
200{
201 sort(p.begin(),p.end());
202 std::vector<int>::iterator it = std::unique(p.begin(),p.end());
203 p.erase(it,p.end());
204}
205
206template<class CM>
207bool IntersectionComputation<CM>::newFace3D(const std::vector<int> id,
208 const std::vector<std::vector<int> > H)
209{
210 // get size_type for all the vectors we are using
211 typedef typename std::vector<Empty>::size_type size_type;
212
213 int n = H.size() ;
214 int m = id.size() ;
215 std::vector<int> A ;
216 std::vector<int> B = id ;
217 sort(B.begin(),B.end()) ;
218 int i = 0 ;
219 bool b = true ;
220 double tp ;
221
222 while ( b && (i<n) )
223 {
224 if ((H[i].size())>=m)
225 {
226 A=H[i] ;
227 sort(A.begin(),A.end());
228 tp = 0 ;
229 for ( size_type j=0 ; j < m; j++)
230 tp += std::fabs(A[j]-B[j]) ;
231 b = (tp>0) ;
232 }
233
234 i += 1 ;
235 }
236
237 return b ;
238}
239
240
241template<class CM>
242void IntersectionComputation<CM>::orderPointsCC(std::integral_constant<int,3>,
243 const V centroid,
244 std::vector<int>& id,
245 const std::vector<V> P)
246{
247 typedef typename std::vector<Empty>::size_type size_type;
248
249 id.clear();
250
251 // get size_type for all the vectors we are using
252 V c,d1,d2,dr,dn,cross,d ;
253 std::vector<typename V::value_type> ai ;
254
255 d1 = P[1] - P[0] ; // two reference vectors
256 d2 = P[2] - P[0] ;
257
258 cross[0] = d1[1]*d2[2] - d1[2]*d2[1] ; // cross product
259 cross[1] = d1[2]*d2[0] - d1[0]*d2[2] ;
260 cross[2] = d1[0]*d2[1] - d1[1]*d2[0] ;
261
262 if (((centroid - P[0])*cross)<0) // good orientation ?
263 {
264 dr = d1 ;
265 dr /= dr.two_norm() ; // 'x-axis' unit vector
266 dn = dr ;
267 dn *= -(d2*dr) ;
268 dn += d2 ;
269 dn /= dn.two_norm() ; // 'y-axis' unit vector
270 }
271 else
272 {
273 dr = d2 ;
274 dr /= dr.two_norm() ; // 'y-axis' unit vector
275 dn = dr ;
276 dn *= -(d1*dr) ;
277 dn += d1 ;
278 dn /= dn.two_norm() ; // 'x-axis' unit vector
279 }
280
281 // definition of angles, using projection on the local reference, ie by scalarly multipliying by dr and dn resp.
282 for ( size_type j=1 ; j < P.size() ; j++)
283 {
284 ai.push_back(atan2((P[j]-P[0])*dn,(P[j]-P[0])*dr)) ;
285 id.push_back(j) ;
286 }
287
288 // sort according to increasing angles
289 for ( size_type j=1; j < ai.size(); j++) {
290 for ( size_type i=0; i < j; i++) {
291 if (ai[j]<ai[i]) {
292 std::swap<typename V::value_type>(ai[i],ai[j]) ;
293 std::swap<int>(id[i],id[j]) ;
294 }
295 }
296 }
297
298 id.insert(id.begin(),0);
299}
300
301template<class CM>
302void IntersectionComputation<CM>::orderPointsCC(std::integral_constant<int,2>,
303 const V centroid,
304 std::vector<int>& id,
305 const std::vector<V> P)
306{
307 typedef typename std::vector<Empty>::size_type size_type;
308
309 // get size_type for all the vectors we are using
310 typedef typename std::vector<Empty>::size_type size_type;
311
312 std::vector<typename V::value_type> ai(P.size());
313 id.resize(P.size());
314
315 // definition of angles
316 for ( size_type i=0; i < P.size(); i++) {
317 ai[i] = atan2(P[i][1]-centroid[1],P[i][0]-centroid[0]);
318 id[i] = i;
319 }
320
321 // sort according to increasing angles
322 for ( size_type j=1; j < ai.size(); j++) {
323 for ( size_type i=0; i < j; i++) if (ai[j]<ai[i]) {
324 std::swap<typename V::value_type>(ai[i],ai[j]);
325 std::swap<int>(id[i],id[j]);
326 }
327 }
328}
329
330} /* namespace Dune::GridGlue */
331} /* namespace Dune */
Definition gridglue.hh:30
int insertPoint(const V p, std::vector< V > &P)
Definition computeintersection.hh:162
Intersection computation method for two elements of arbitrary dimension.
Definition computeintersection.hh:37
static bool computeIntersection(const std::vector< V > X, const std::vector< V > Y, std::vector< std::vector< int > > &SX, std::vector< std::vector< int > > &SY, std::vector< V > &P)
Compute the intersection of two elements X and Y Compute the intersection of two elements X and Y,...
Definition computeintersection.cc:12