33 const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
34 std::bitset<(1<<dim1)>& neighborIntersects1,
35 unsigned int grid1Index,
36 const Dune::GeometryType& grid2ElementType,
37 const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
38 std::bitset<(1<<dim2)>& neighborIntersects2,
39 unsigned int grid2Index,
50 const Dune::ReferenceElement<T,dim1>& refElement1 = Dune::ReferenceElements<T,dim1>::general(grid1ElementType);
51 const Dune::ReferenceElement<T,dim2>& refElement2 = Dune::ReferenceElements<T,dim2>::general(grid2ElementType);
54 assert((
unsigned int)(refElement1.size(dim1)) == grid1ElementCorners.size());
55 assert((
unsigned int)(refElement2.size(dim2)) == grid2ElementCorners.size());
60 typedef MultiLinearGeometry<T,dim1,dimworld> Geometry1;
61 typedef MultiLinearGeometry<T,dim2,dimworld> Geometry2;
63 Geometry1 grid1Geometry(grid1ElementType, grid1ElementCorners);
64 Geometry2 grid2Geometry(grid2ElementType, grid2ElementCorners);
67 std::vector<Dune::FieldVector<T,dimworld> > scaledGrid1ElementCorners(grid1ElementCorners.size());
68 std::vector<Dune::FieldVector<T,dimworld> > scaledGrid2ElementCorners(grid2ElementCorners.size());
72 T scaling = min((grid1ElementCorners[0] - grid1ElementCorners[1]).two_norm(),
73 (grid2ElementCorners[0] - grid2ElementCorners[1]).two_norm());
76 for (
unsigned int i = 0; i < grid1ElementCorners.size(); ++i) {
77 scaledGrid1ElementCorners[i] = grid1ElementCorners[i];
78 scaledGrid1ElementCorners[i] *= (1.0/scaling);
80 for (
unsigned int i = 0; i < grid2ElementCorners.size(); ++i) {
81 scaledGrid2ElementCorners[i] = grid2ElementCorners[i];
82 scaledGrid2ElementCorners[i] *= (1.0/scaling);
86 typedef typename std::vector<Empty>::size_type size_type;
88 const int dimis = dim1 < dim2 ? dim1 : dim2;
89 const size_type n_intersectionnodes = dimis+1;
92 std::vector<FieldVector<T,dimworld> > scaledP(0), P(0);
93 std::vector<std::vector<int> > H,SX(1<<dim1),SY(1<<dim2);
94 FieldVector<T,dimworld> centroid;
96 std::vector<FieldVector<T,dim1> > g1local(n_intersectionnodes);
98 std::vector<FieldVector<T,dim2> > g2local(n_intersectionnodes);
102 scaledGrid2ElementCorners,
106 P.resize(scaledP.size());
107 for (
unsigned int i = 0; i < scaledP.size(); ++i) {
112 for (size_type i = 0; i < neighborIntersects1.size(); ++i) {
114 neighborIntersects1[i] = (SX[i].size() > 0);
116 neighborIntersects1[i] =
false;
118 for (size_type i = 0; i < neighborIntersects2.size(); ++i) {
120 neighborIntersects2[i] = (SY[i].size() > 0);
122 neighborIntersects2[i] =
false;
126 if (P.size() == n_intersectionnodes) {
128 for (i = 0; i < n_intersectionnodes; ++i) {
129 g1local[i] = grid1Geometry.local(P[i]);
130 g2local[i] = grid2Geometry.local(P[i]);
134 for (i = 0; i < n_intersectionnodes; ++i) {
139 }
else if (P.size() > n_intersectionnodes) {
142 std::vector<FieldVector<T,dimworld> > global(n_intersectionnodes);
146 for (size_type i=0; i < P.size(); i++)
148 centroid /=
static_cast<T
>(P.size()) ;
155 for (size_type i=0; i < H.size(); i++) {
156 int hs = H[i].size();
162 for ( size_type j=0 ; j < dimis; ++j) {
163 global[j]= P[H[i][j]];
167 global[dimis]=centroid;
170 for (size_type j = 0; j < n_intersectionnodes; ++j) {
171 g1local[j] = grid1Geometry.local(global[j]);
172 g2local[j] = grid2Geometry.local(global[j]);
176 for (size_type j = 0; j < n_intersectionnodes; ++j) {
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< RemoteSimplicialIntersection > &intersections)
Compute the intersection between two overlapping elements.
Definition overlappingmerge.cc:32