dune-grid-glue 2.10
Loading...
Searching...
No Matches
overlappingmerge.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
6#ifndef DUNE_GRIDGLUE_OVERLAPPINGMERGE_CC
7#define DUNE_GRIDGLUE_OVERLAPPINGMERGE_CC
8//#include <algorithm>
9
10namespace Dune {
11namespace GridGlue {
12
13template<int dim1, int dim2, int dimworld, typename T>
14bool OverlappingMerge<dim1,dim2,dimworld, T>::inPlane(std::vector<FieldVector<T,dimworld> >& points) {
15
16 T eps = 1e-8;
17
18 assert(dim1 == 3 && dim2 == 3 && dimworld == 3);
19 assert(points.size() == 4);
20
21 FieldVector<T,dimworld> v1 = points[1]-points[0];
22 FieldVector<T,dimworld> v2 = points[2]-points[0];
23 FieldVector<T,dimworld> v3 = points[3]-points[0];
24
25 FieldVector<T,dimworld> v1xv2;
26 v1xv2[0] = v1[1]*v2[2] - v1[2]*v2[1];
27 v1xv2[1] = v1[2]*v2[0] - v1[0]*v2[2];
28 v1xv2[2] = v1[0]*v2[1] - v1[1]*v2[0];
29
30 return (std::abs(v3.dot(v1xv2)) < eps);
31}
32
33template<int dim1, int dim2, int dimworld, typename T>
34void OverlappingMerge<dim1,dim2,dimworld, T>::computeIntersections(const Dune::GeometryType& grid1ElementType,
35 const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
36 std::bitset<(1<<dim1)>& neighborIntersects1,
37 unsigned int grid1Index,
38 const Dune::GeometryType& grid2ElementType,
39 const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
40 std::bitset<(1<<dim2)>& neighborIntersects2,
41 unsigned int grid2Index,
42 std::vector<SimplicialIntersection>& intersections)
43{
44 using std::min;
45
46 this->counter++;
47 intersections.clear();
48
50
51#ifndef NDEBUG
52 const auto& refElement1 = Dune::ReferenceElements<T,dim1>::general(grid1ElementType);
53 const auto& refElement2 = Dune::ReferenceElements<T,dim2>::general(grid2ElementType);
54
55 // A few consistency checks
56 assert((unsigned int)(refElement1.size(dim1)) == grid1ElementCorners.size());
57 assert((unsigned int)(refElement2.size(dim2)) == grid2ElementCorners.size());
58#endif
59
60 // Make generic geometries representing the grid1- and grid2 element.
61 // this eases computation of local coordinates.
62 typedef MultiLinearGeometry<T,dim1,dimworld> Geometry1;
63 typedef MultiLinearGeometry<T,dim2,dimworld> Geometry2;
64
65 Geometry1 grid1Geometry(grid1ElementType, grid1ElementCorners);
66 Geometry2 grid2Geometry(grid2ElementType, grid2ElementCorners);
67
68 // Dirty workaround for the intersection computation scaling problem (part 1/2)
69 std::vector<Dune::FieldVector<T,dimworld> > scaledGrid1ElementCorners(grid1ElementCorners.size());
70 std::vector<Dune::FieldVector<T,dimworld> > scaledGrid2ElementCorners(grid2ElementCorners.size());
71
72 // the scaling parameter is the length minimum of the lengths of the first edge in the grid1 geometry
73 // and the first edge in the grid2 geometry
74 T scaling = min((grid1ElementCorners[0] - grid1ElementCorners[1]).two_norm(),
75 (grid2ElementCorners[0] - grid2ElementCorners[1]).two_norm());
76
77 // scale the coordinates according to scaling parameter
78 for (unsigned int i = 0; i < grid1ElementCorners.size(); ++i) {
79 scaledGrid1ElementCorners[i] = grid1ElementCorners[i];
80 scaledGrid1ElementCorners[i] *= (1.0/scaling);
81 }
82 for (unsigned int i = 0; i < grid2ElementCorners.size(); ++i) {
83 scaledGrid2ElementCorners[i] = grid2ElementCorners[i];
84 scaledGrid2ElementCorners[i] *= (1.0/scaling);
85 }
86
87 // get size_type for all the vectors we are using
88 typedef typename std::vector<Empty>::size_type size_type;
89
90 const int dimis = dim1 < dim2 ? dim1 : dim2;
91 const size_type n_intersectionnodes = dimis+1;
92 size_type i;
93
94 std::vector<FieldVector<T,dimworld> > scaledP(0), P(0);
95 std::vector<std::vector<int> > H,SX(1<<dim1),SY(1<<dim2);
96 FieldVector<T,dimworld> centroid;
97 // local grid1 coordinates of the intersections
98 std::vector<FieldVector<T,dim1> > g1local(n_intersectionnodes);
99 // local grid2 coordinates of the intersections
100 std::vector<FieldVector<T,dim2> > g2local(n_intersectionnodes);
101
102 // compute the intersection nodes
103 IntersectionComputation<CM>::computeIntersection(scaledGrid1ElementCorners,
104 scaledGrid2ElementCorners,
105 SX,SY,scaledP);
106
107 // Dirty workaround - rescale the result (part 2/2)
108 P.resize(scaledP.size());
109 for (unsigned int i = 0; i < scaledP.size(); ++i) {
110 P[i] = scaledP[i];
111 P[i] *= scaling;
112 }
113
114 for (size_type i = 0; i < neighborIntersects1.size(); ++i) {
115 if (i < SX.size())
116 neighborIntersects1[i] = (SX[i].size() > 0);
117 else
118 neighborIntersects1[i] = false;
119 }
120 for (size_type i = 0; i < neighborIntersects2.size(); ++i) {
121 if (i < SY.size())
122 neighborIntersects2[i] = (SY[i].size() > 0);
123 else
124 neighborIntersects2[i] = false;
125 }
126
127 // P is an simplex of dimension dimis
128 if (P.size() == n_intersectionnodes) {
129
130 for (i = 0; i < n_intersectionnodes; ++i) {
131 g1local[i] = grid1Geometry.local(P[i]);
132 g2local[i] = grid2Geometry.local(P[i]);
133 }
134
135 intersections.emplace_back(grid1Index, grid2Index);
136 for (i = 0; i < n_intersectionnodes; ++i) {
137 intersections.back().corners0[0][i] = g1local[i];
138 intersections.back().corners1[0][i] = g2local[i];
139 }
140
141 } else if (P.size() > n_intersectionnodes) { // P is a union of simplices of dimension dimis
142
143 assert(dimis != 1);
144 std::vector<FieldVector<T,dimworld> > global(n_intersectionnodes);
145
146 // Compute the centroid
147 centroid=0;
148 for (size_type i=0; i < P.size(); i++)
149 centroid += P[i] ;
150 centroid /= static_cast<T>(P.size()) ;
151
152 // order the points and get intersection face indices
153 H.clear() ;
154 IntersectionComputation<CM>::template orderPoints<dimis,dimworld>(centroid,SX,SY,P,H);
155
156 // Loop over all intersection elements
157 for (size_type i=0; i < H.size(); i++) {
158 int hs = H[i].size(); // number of nodes of the intersection
159
160 // if the intersection element is not degenerated
161 if (hs==dimis) {
162
163 // create the intersection geometry
164 for ( size_type j=0 ; j < dimis; ++j) {
165 global[j]= P[H[i][j]]; // get the intersection face
166 }
167
168 // intersection face + centroid = new element
169 global[dimis]=centroid;
170
171 // create local representation of the intersection
172 for (size_type j = 0; j < n_intersectionnodes; ++j) {
173 g1local[j] = grid1Geometry.local(global[j]);
174 g2local[j] = grid2Geometry.local(global[j]);
175 }
176
177 intersections.emplace_back(grid1Index,grid2Index);
178 for (size_type j = 0; j < n_intersectionnodes; ++j) {
179 intersections.back().corners0[0][j] = g1local[j];
180 intersections.back().corners1[0][j] = g2local[j];
181 }
182 }
183 }
184 }
185}
186
187} /* namespace Dune::GridGlue */
188} /* namespace Dune */
189
190#endif // DUNE_GRIDGLUE_OVERLAPPINGMERGE_CC
Definition gridglue.hh:37
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
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
void computeIntersections(const Dune::GeometryType &grid1ElementType, const std::vector< Dune::FieldVector< T, dimworld > > &grid1ElementCorners, std::bitset<(1<< dim1)> &neighborIntersects1, unsigned int grid1Index, const Dune::GeometryType &grid2ElementType, const std::vector< Dune::FieldVector< T, dimworld > > &grid2ElementCorners, std::bitset<(1<< dim2)> &neighborIntersects2, unsigned int grid2Index, std::vector< SimplicialIntersection > &intersections)
Compute the intersection between two overlapping elements.
Definition overlappingmerge.cc:34
Definition simplexintersection.cc:30