dune-grid-glue 2.5-git
contactmerge.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
6
7namespace Dune {
8namespace GridGlue {
9
10template<int dimworld, typename T>
11void ContactMerge<dimworld, T>::computeIntersections(const Dune::GeometryType& grid1ElementType,
12 const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
13 std::bitset<(1<<dim)>& neighborIntersects1,
14 unsigned int grid1Index,
15 const Dune::GeometryType& grid2ElementType,
16 const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
17 std::bitset<(1<<dim)>& neighborIntersects2,
18 unsigned int grid2Index,
19 std::vector<RemoteSimplicialIntersection>& intersections)
20{
21 using std::get;
22
23 std::vector<std::array<LocalCoords,2> > polytopeCorners;
24
25 // Initialize
26 neighborIntersects1.reset();
27 neighborIntersects2.reset();
28
29 const int nCorners1 = grid1ElementCorners.size();
30 const int nCorners2 = grid2ElementCorners.size();
31
32 if (nCorners1 != dimworld)
33 DUNE_THROW(Dune::Exception, "element1 must have " << dimworld << " corners, but has " << nCorners1);
34 if (nCorners2 != dimworld)
35 DUNE_THROW(Dune::Exception, "element2 must have " << dimworld << " corners, but has " << nCorners2);
36
37 // The grid1 projection directions
38 std::vector<WorldCoords> directions1(nCorners1);
39 for (size_t i=0; i<directions1.size(); i++)
40 directions1[i] = nodalDomainDirections_[this->grid1ElementCorners_[grid1Index][i]];
41
42 // The grid2 projection directions
43 // These are only needed in a check for validity of the projection and they should be chosen
44 // to be some outer pointing directions, e.g. outer normals
45 std::vector<WorldCoords> directions2(nCorners2);
46 for (size_t i=0; i<directions2.size(); i++)
47 directions2[i] = nodalTargetDirections_[this->grid2ElementCorners_[grid2Index][i]];
48
50 // Compute all corners of the intersection polytope
52
53 const auto corners = std::tie(grid1ElementCorners, grid2ElementCorners);
54 const auto normals = std::tie(directions1, directions2);
55 Projection<WorldCoords> p(overlap_, maxNormalProduct_);
56 p.project(corners, normals);
57
58 /* projection */
59 {
60 const auto& success = get<0>(p.success());
61 const auto& images = get<0>(p.images());
62 for (unsigned i = 0; i < dimworld; ++i) {
63 if (success[i]) {
64 std::array<LocalCoords, 2> corner;
65 corner[0] = localCornerCoords(i, grid1ElementType);
66 for (unsigned j = 0; j < dim; ++j)
67 corner[1][j] = images[i][j];
68 polytopeCorners.push_back(corner);
69 }
70 }
71 }
72
73 /* inverse projection */
74 {
75 const auto& success = get<1>(p.success());
76 const auto& preimages = get<1>(p.images());
77 for (unsigned i = 0; i < dimworld; ++i) {
78 if (success[i]) {
79 std::array<LocalCoords, 2> corner;
80 for (unsigned j = 0; j < dim; ++j)
81 corner[0][j] = preimages[i][j];
82 corner[1] = localCornerCoords(i, grid2ElementType);
83 polytopeCorners.push_back(corner);
84 }
85 }
86 }
87
88 /* edge intersections */
89 {
90 for (unsigned i = 0; i < p.numberOfEdgeIntersections(); ++i) {
91 std::array<LocalCoords, 2> corner;
92 const auto& local = p.edgeIntersections()[i].local;
93 for (unsigned j = 0; j < dim; ++j) {
94 corner[0][j] = local[0][j];
95 corner[1][j] = local[1][j];
96 }
97 polytopeCorners.push_back(corner);
98 }
99 }
100
101 // check which neighbors might also intersect
102 const Dune::ReferenceElement<T,dim>& ref2 = Dune::ReferenceElements<T,dim>::general(grid2ElementType);
103 for (int i=0; i<ref2.size(1); i++) {
104
105 // if all face corners hit the the other element then
106 // the neighbor might also intersect
107
108 bool intersects(true);
109 for (int k=0; k<ref2.size(i,1,dim); k++)
110 intersects &= get<1>(p.success())[ref2.subEntity(i,1,k,dim)];
111
112 if (intersects)
113 neighborIntersects2[i] = true;
114 }
115
116 const Dune::ReferenceElement<T,dim>& ref1 = Dune::ReferenceElements<T,dim>::general(grid1ElementType);
117 for (int i=0; i<ref1.size(1); i++) {
118
119 // if all face corners hit the the other element then
120 // the neighbor might also intersect
121
122 bool intersects(true);
123 for (int k=0; k<ref1.size(i,1,dim); k++)
124 intersects &= get<0>(p.success())[ref1.subEntity(i,1,k,dim)];
125
126 if (intersects)
127 neighborIntersects1[i] = true;
128 }
129
130 // Compute the edge intersections
131 for (unsigned i = 0; i < p.numberOfEdgeIntersections(); ++i) {
132 const auto& edge = p.edgeIntersections()[i].edge;
133 neighborIntersects1[edge[0]] = true;
134 neighborIntersects2[edge[1]] = true;
135 }
136
137 // remove possible doubles
138 removeDoubles(polytopeCorners);
139
140 // Compute an interior point of the polytope
141 int nPolyCorners = polytopeCorners.size();
142
143 // If the polytope is degenerated then there is no intersection
144 if (nPolyCorners<dimworld)
145 return;
146
147 // If the polytope is a simplex return it
148 if (nPolyCorners==dim+1) {
149
150 // std::cout<<"Add intersection: 1\n";
151 typename Base::RemoteSimplicialIntersection intersect;
152 intersect.grid1Entities_[0] = grid1Index;
153 intersect.grid2Entities_[0] = grid2Index;
154
155 for (int j=0;j<dim+1; j++) {
156 intersect.grid1Local_[0][j]=polytopeCorners[j][0];
157 intersect.grid2Local_[0][j]=polytopeCorners[j][1];
158 }
159 intersections.push_back(intersect);
160
161 return;
162 }
163
164 // At this point we must have dimworld>=3
165
167 // Compute a point in the middle of the polytope and order all corners cyclic
169
170 std::array<LocalCoords,2> center;
171 center[0] = 0; center[1] = 0;
172 for (int i=0; i<nPolyCorners; i++) {
173 center[0].axpy(1.0/nPolyCorners,polytopeCorners[i][0]);
174 center[1].axpy(1.0/nPolyCorners,polytopeCorners[i][1]);
175 }
176
177 // Order cyclic
178 std::vector<int> ordering;
179 computeCyclicOrder(polytopeCorners,center[0],ordering);
180
182 // Add intersections
184
185 for (size_t i=1; i<polytopeCorners.size()-1; i++) {
186
187 typename Base::RemoteSimplicialIntersection intersect;
188 intersect.grid1Entities_[0] = grid1Index;
189 intersect.grid2Entities_[0] = grid2Index;
190
191 for (int j=0;j<dim; j++) {
192 intersect.grid1Local_[0][j]=polytopeCorners[ordering[i+j]][0];
193 intersect.grid2Local_[0][j]=polytopeCorners[ordering[i+j]][1];
194 }
195
196 // last corner is the first for all intersections
197 intersect.grid1Local_[0][dim]=polytopeCorners[ordering[0]][0];
198 intersect.grid2Local_[0][dim]=polytopeCorners[ordering[0]][1];
199
200 intersections.push_back(intersect);
201 }
202}
203
204template<int dimworld, typename T>
205void ContactMerge<dimworld, T>::computeCyclicOrder(const std::vector<std::array<LocalCoords,2> >& polytopeCorners,
206 const LocalCoords& center, std::vector<int>& ordering) const
207{
208 ordering.resize(polytopeCorners.size());
209
210 for (size_t k=0; k<ordering.size(); k++)
211 ordering[k] = k;
212
213 //TODO Do I have to order triangles to get some correct orientation?
214 if (polytopeCorners.size()<=3)
215 return;
216
217 // compute angles inside the polygon plane w.r.t to this axis
218 LocalCoords edge0 = polytopeCorners[1][0] - polytopeCorners[0][0];
219
220 // Compute a vector that is perpendicular to the edge but lies in the polytope plane
221 // So we have a unique ordering
222 LocalCoords edge1 = polytopeCorners[2][0] - polytopeCorners[0][0];
223 LocalCoords normal0 = edge1;
224 normal0.axpy(-(edge0*edge1),edge0);
225
226 std::vector<T> angles(polytopeCorners.size());
227
228 for (size_t i=0; i<polytopeCorners.size(); i++) {
229
230 LocalCoords edge = polytopeCorners[i][0] - center;
231
232 T x(edge*edge0);
233 T y(edge*normal0);
234
235 angles[i] = std::atan2(y, x);
236 if (angles[i]<0)
237 angles[i] += 2*M_PI;
238 }
239
240 // bubblesort
241
242 for (int i=polytopeCorners.size(); i>1; i--){
243 bool swapped = false;
244
245 for (int j=0; j<i-1; j++){
246
247 if (angles[j] > angles[j+1]){
248 swapped = true;
249 std::swap(angles[j], angles[j+1]);
250 std::swap(ordering[j], ordering[j+1]);
251 }
252 }
253
254 if (!swapped)
255 break;
256 }
257}
258
259template<int dimworld, typename T>
260void ContactMerge<dimworld, T>::setupNodalDirections(const std::vector<WorldCoords>& coords1,
261 const std::vector<unsigned int>& elements1,
262 const std::vector<Dune::GeometryType>& elementTypes1,
263 const std::vector<WorldCoords>& coords2,
264 const std::vector<unsigned int>& elements2,
265 const std::vector<Dune::GeometryType>& elementTypes2)
266{
267 if (domainDirections_) {
268
269 // Sample the provided analytical contact direction field
270 nodalDomainDirections_.resize(coords1.size());
271 for (size_t i=0; i<coords1.size(); i++)
272 nodalDomainDirections_[i] = domainDirections_(coords1[i]);
273 } else
274 computeOuterNormalField(coords1,elements1,elementTypes1, nodalDomainDirections_);
275
276 if (targetDirections_) {
277
278 // Sample the provided analytical target direction field
279 nodalTargetDirections_.resize(coords2.size());
280 for (size_t i=0; i<coords2.size(); i++)
281 nodalTargetDirections_[i] = targetDirections_(coords2[i]);
282 } else
283 computeOuterNormalField(coords2,elements2,elementTypes2, nodalTargetDirections_);
284}
285
286template<int dimworld, typename T>
287void ContactMerge<dimworld, T>::computeOuterNormalField(const std::vector<WorldCoords>& coords,
288 const std::vector<unsigned int>& elements,
289 const std::vector<Dune::GeometryType>& elementTypes,
290 std::vector<WorldCoords>& normals)
291{
292 normals.assign(coords.size(),WorldCoords(0));
293
294
295 int offset = 0;
296
297 for (size_t i=0; i<elementTypes.size(); i++) {
298
299 int nCorners = Dune::ReferenceElements<T,dim>::general(elementTypes[i]).size(dim);
300
301 // For segments 1, for triangles or quadrilaterals take the first 2
302 std::array<WorldCoords, dim> edges;
303 for (int j=1; j<=dim; j++)
304 edges[j-1] = coords[elements[offset + j]] - coords[elements[offset]];
305
306 WorldCoords elementNormal;
307
308 if (dim==1) {
309 elementNormal[0] = edges[0][1]; elementNormal[1] = -edges[0][0];
310 } else
311 elementNormal = crossProduct(edges[0], edges[1]);
312
313 elementNormal /= elementNormal.two_norm();
314
315 for (int j=0; j<nCorners;j++)
316 normals[elements[offset + j]] += elementNormal;
317
318 offset += nCorners;
319 }
320
321 for (size_t i=0; i<coords.size(); i++)
322 normals[i] /= normals[i].two_norm();
323}
324
325template<int dimworld, typename T>
326void ContactMerge<dimworld, T>::removeDoubles(std::vector<std::array<LocalCoords,2> >& polytopeCorners)
327{
328
329 size_t counter(1);
330 for (size_t i=1; i<polytopeCorners.size(); i++) {
331 bool contained = false;
332 for (size_t j=0; j<counter; j++)
333 if ( (polytopeCorners[j][0]-polytopeCorners[i][0]).two_norm()<1e-10) {
334 assert((polytopeCorners[j][1]-polytopeCorners[i][1]).two_norm()<1e-10);
335 contained = true;
336 break;
337 }
338
339 if (!contained) {
340 if (counter < i)
341 polytopeCorners[counter] = polytopeCorners[i];
342 counter++;
343 }
344 }
345 polytopeCorners.resize(counter);
346}
347
348} /* namespace GridGlue */
349} /* namespace Dune */
Definition gridglue.hh:30
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
static Dune::FieldVector< T, dim > crossProduct(const Dune::FieldVector< T, dim > &a, const Dune::FieldVector< T, dim > &b)
compute cross product
Definition crossproduct.hh:13
Coordinate corner(unsigned c)
Definition projection_impl.hh:22
void computeCyclicOrder(const std::vector< std::array< LocalCoords, 2 > > &polytopeCorners, const LocalCoords &center, std::vector< int > &ordering) const
Order the corners of the intersection polytope in cyclic order.
Definition contactmerge.cc:205
void removeDoubles(std::vector< std::array< LocalCoords, 2 > > &polytopeCorners)
Remove all multiples.
Definition contactmerge.cc:326
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition contactmerge.hh:58
void setupNodalDirections(const std::vector< WorldCoords > &coords1, const std::vector< unsigned int > &elements1, const std::vector< Dune::GeometryType > &elementTypes1, const std::vector< WorldCoords > &coords2, const std::vector< unsigned int > &elements2, const std::vector< Dune::GeometryType > &elementTypes2)
Setup the direction vectors containing the directions for each vertex.
Definition contactmerge.cc:260
void computeOuterNormalField(const std::vector< WorldCoords > &coords, const std::vector< unsigned int > &elements, const std::vector< Dune::GeometryType > &elementTypes, std::vector< WorldCoords > &normals)
If no direction field was specified compute the outer normal field.
Definition contactmerge.cc:287
Dune::FieldVector< T, dim > LocalCoords
the coordinate type used in this interface
Definition contactmerge.hh:61