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