12template<
int dimworld,
typename T>
13void ContactMerge<dimworld, T>::computeIntersections(
const Dune::GeometryType& grid1ElementType,
14 const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
15 std::bitset<(1<<dim)>& neighborIntersects1,
16 unsigned int grid1Index,
17 const Dune::GeometryType& grid2ElementType,
18 const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
19 std::bitset<(1<<dim)>& neighborIntersects2,
20 unsigned int grid2Index,
25 std::vector<std::array<LocalCoords,2> > polytopeCorners;
28 neighborIntersects1.reset();
29 neighborIntersects2.reset();
31 const int nCorners1 = grid1ElementCorners.size();
32 const int nCorners2 = grid2ElementCorners.size();
34 if (nCorners1 != dimworld)
35 DUNE_THROW(Dune::Exception,
"element1 must have " << dimworld <<
" corners, but has " << nCorners1);
36 if (nCorners2 != dimworld)
37 DUNE_THROW(Dune::Exception,
"element2 must have " << dimworld <<
" corners, but has " << nCorners2);
40 std::vector<WorldCoords> directions1(nCorners1);
41 for (
size_t i=0; i<directions1.size(); i++)
42 directions1[i] = nodalDomainDirections_[this->grid1ElementCorners_[grid1Index][i]];
45 std::vector<WorldCoords> directions2(nCorners2);
46 for (
size_t i=0; i<directions2.size(); i++)
47 directions2[i] = nodalTargetDirections_[this->grid2ElementCorners_[grid2Index][i]];
53 std::array<
decltype(std::cref(grid1ElementCorners)),2> cornersRef ={std::cref(grid1ElementCorners), std::cref(grid2ElementCorners)};
54 std::array<
decltype(std::cref(directions1)),2> directionsRef ={std::cref(directions1), std::cref(directions2)};
55 std::array<Dune::GeometryType,2> elementTypes = {grid1ElementType, grid2ElementType};
58 const size_t domGrid = (type_==ProjectionType::OUTER_NORMAL) ? 0 : 1;
59 const size_t tarGrid = (type_==ProjectionType::OUTER_NORMAL) ? 1 : 0;
65 const auto corners = std::tie(cornersRef[domGrid].get(),cornersRef[tarGrid].get());
66 const auto normals = std::tie(directionsRef[domGrid].get(), directionsRef[tarGrid].get());
67 Projection<WorldCoords> p(overlap_, maxNormalProduct_);
68 p.project(corners, normals);
72 const auto& success = get<0>(p.success());
73 const auto& images = get<0>(p.images());
74 for (
unsigned i = 0; i < dimworld; ++i) {
76 std::array<LocalCoords, 2>
corner;
77 corner[domGrid] = localCornerCoords(i, elementTypes[domGrid]);
78 for (
unsigned j = 0; j < dim; ++j)
79 corner[tarGrid][j] = images[i][j];
80 polytopeCorners.push_back(corner);
87 const auto& success = get<1>(p.success());
88 const auto& preimages = get<1>(p.images());
89 for (
unsigned i = 0; i < dimworld; ++i) {
91 std::array<LocalCoords, 2>
corner;
92 for (
unsigned j = 0; j < dim; ++j)
93 corner[domGrid][j] = preimages[i][j];
94 corner[tarGrid] = localCornerCoords(i, elementTypes[tarGrid]);
95 polytopeCorners.push_back(corner);
102 for (
unsigned i = 0; i < p.numberOfEdgeIntersections(); ++i) {
103 std::array<LocalCoords, 2>
corner;
104 const auto& local = p.edgeIntersections()[i].local;
105 for (
unsigned j = 0; j < dim; ++j) {
106 corner[domGrid][j] = local[0][j];
107 corner[tarGrid][j] = local[1][j];
109 polytopeCorners.push_back(corner);
114 std::array<
decltype(std::ref(neighborIntersects1)),2> neighborIntersectsRef = {std::ref(neighborIntersects1), std::ref(neighborIntersects2)};
115 const auto& refTar = Dune::ReferenceElements<T,dim>::general(elementTypes[tarGrid]);
116 for (
int i=0; i<refTar.size(1); i++) {
121 bool intersects(
true);
122 for (
int k=0; k<refTar.size(i,1,dim); k++)
123 intersects &= get<1>(p.success())[refTar.subEntity(i,1,k,dim)];
126 neighborIntersectsRef[tarGrid].get()[i] =
true;
129 const auto& refDom = Dune::ReferenceElements<T,dim>::general(elementTypes[domGrid]);
130 for (
int i=0; i<refDom.size(1); i++) {
135 bool intersects(
true);
136 for (
int k=0; k<refDom.size(i,1,dim); k++)
137 intersects &= get<0>(p.success())[refDom.subEntity(i,1,k,dim)];
140 neighborIntersectsRef[domGrid].get()[i] =
true;
144 for (
unsigned i = 0; i < p.numberOfEdgeIntersections(); ++i) {
145 const auto& edge = p.edgeIntersections()[i].edge;
146 neighborIntersects1[edge[domGrid]] =
true;
147 neighborIntersects2[edge[tarGrid]] =
true;
151 removeDoubles(polytopeCorners);
154 int nPolyCorners = polytopeCorners.size();
157 if (nPolyCorners<dimworld)
161 if (nPolyCorners==dim+1) {
164 typename Base::SimplicialIntersection intersect(grid1Index, grid2Index);
166 for (
int j=0;j<dim+1; j++) {
167 intersect.corners0[0][j]=polytopeCorners[j][0];
168 intersect.corners1[0][j]=polytopeCorners[j][1];
181 std::array<LocalCoords,2> center;
182 center[0] = 0; center[1] = 0;
183 for (
int i=0; i<nPolyCorners; i++) {
184 center[0].axpy(1.0/nPolyCorners,polytopeCorners[i][0]);
185 center[1].axpy(1.0/nPolyCorners,polytopeCorners[i][1]);
189 std::vector<int> ordering;
190 computeCyclicOrder(polytopeCorners,center[0],ordering);
196 for (
size_t i=1; i<polytopeCorners.size()-1; i++) {
198 typename Base::SimplicialIntersection intersect(grid1Index, grid2Index);
200 for (
int j=0;j<dim; j++) {
201 intersect.corners0[0][j]=polytopeCorners[ordering[i+j]][0];
202 intersect.corners1[0][j]=polytopeCorners[ordering[i+j]][1];
206 intersect.corners0[0][dim]=polytopeCorners[ordering[0]][0];
207 intersect.corners1[0][dim]=polytopeCorners[ordering[0]][1];
213template<
int dimworld,
typename T>
215 const LocalCoords& center, std::vector<int>& ordering)
const
217 ordering.resize(polytopeCorners.size());
219 for (
size_t k=0; k<ordering.size(); k++)
223 if (polytopeCorners.size()<=3)
227 LocalCoords edge0 = polytopeCorners[1][0] - polytopeCorners[0][0];
231 LocalCoords edge1 = polytopeCorners[2][0] - polytopeCorners[0][0];
233 normal0.axpy(-(edge0*edge1),edge0);
235 std::vector<T> angles(polytopeCorners.size());
237 for (
size_t i=0; i<polytopeCorners.size(); i++) {
244 angles[i] = std::atan2(y, x);
251 for (
int i=polytopeCorners.size(); i>1; i--){
252 bool swapped =
false;
254 for (
int j=0; j<i-1; j++){
256 if (angles[j] > angles[j+1]){
258 std::swap(angles[j], angles[j+1]);
259 std::swap(ordering[j], ordering[j+1]);
268template<
int dimworld,
typename T>
270 const std::vector<unsigned int>& elements1,
271 const std::vector<Dune::GeometryType>& elementTypes1,
272 const std::vector<WorldCoords>& coords2,
273 const std::vector<unsigned int>& elements2,
274 const std::vector<Dune::GeometryType>& elementTypes2)
276 if (domainDirections_) {
279 nodalDomainDirections_.resize(coords1.size());
280 for (
size_t i=0; i<coords1.size(); i++)
281 nodalDomainDirections_[i] = domainDirections_(coords1[i]);
283 computeOuterNormalField(coords1,elements1,elementTypes1, nodalDomainDirections_);
285 if (targetDirections_) {
288 nodalTargetDirections_.resize(coords2.size());
289 for (
size_t i=0; i<coords2.size(); i++)
290 nodalTargetDirections_[i] = targetDirections_(coords2[i]);
292 computeOuterNormalField(coords2,elements2,elementTypes2, nodalTargetDirections_);
295template<
int dimworld,
typename T>
297 const std::vector<unsigned int>& elements,
298 const std::vector<Dune::GeometryType>& elementTypes,
299 std::vector<WorldCoords>& normals)
306 for (
size_t i=0; i<elementTypes.size(); i++) {
308 int nCorners = Dune::ReferenceElements<T,dim>::general(elementTypes[i]).size(dim);
311 std::array<WorldCoords, dim> edges;
312 for (
int j=1; j<=dim; j++)
313 edges[j-1] = coords[elements[offset + j]] - coords[elements[offset]];
318 elementNormal[0] = edges[0][1]; elementNormal[1] = -edges[0][0];
322 elementNormal /= elementNormal.two_norm();
324 for (
int j=0; j<nCorners;j++)
325 normals[elements[offset + j]] += elementNormal;
330 for (
size_t i=0; i<coords.size(); i++)
331 normals[i] /= normals[i].two_norm();
334template<
int dimworld,
typename T>
339 for (
size_t i=1; i<polytopeCorners.size(); i++) {
340 bool contained =
false;
341 for (
size_t j=0; j<counter; j++)
342 if ( (polytopeCorners[j][0]-polytopeCorners[i][0]).two_norm()<1e-10) {
343 assert((polytopeCorners[j][1]-polytopeCorners[i][1]).two_norm()<1e-10);
350 polytopeCorners[counter] = polytopeCorners[i];
354 polytopeCorners.resize(counter);
Definition gridglue.hh:37
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:15
Coordinate corner(unsigned c)
Definition projection_impl.hh:24
void computeCyclicOrder(const std::vector< std::array< LocalCoords, 2 > > &polytopeCorners, const LocalCoords ¢er, std::vector< int > &ordering) const
Order the corners of the intersection polytope in cyclic order.
Definition contactmerge.cc:214
void removeDoubles(std::vector< std::array< LocalCoords, 2 > > &polytopeCorners)
Remove all multiples.
Definition contactmerge.cc:335
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition contactmerge.hh:59
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:269
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:296
Dune::FieldVector< T, dim > LocalCoords
the coordinate type used in this interface
Definition contactmerge.hh:62