dune-grid-glue  2.5-git
contactmerge.hh
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:
8 #ifndef DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
9 #define DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
10 
11 
12 #include <iostream>
13 #include <fstream>
14 #include <iomanip>
15 #include <vector>
16 #include <algorithm>
17 #include <limits>
18 #include <memory>
19 #include <functional>
20 
21 #include <dune/common/fvector.hh>
22 #include <dune/common/function.hh>
23 #include <dune/common/exceptions.hh>
24 #include <dune/common/bitsetvector.hh>
25 #include <dune/common/deprecated.hh>
26 
27 #include <dune/grid/common/grid.hh>
28 
31 
32 namespace Dune {
33 namespace GridGlue {
34 
40 template<int dimworld, typename T = double>
42 : public StandardMerge<T,dimworld-1,dimworld-1,dimworld>
43 {
44  enum {dim = dimworld-1};
45 
46  static_assert( dim==1 || dim==2,
47  "ContactMerge yet only handles the cases dim==1 and dim==2!");
48 
50 public:
51 
52  /* E X P O R T E D T Y P E S A N D C O N S T A N T S */
53 
55  typedef T ctype;
56 
58  typedef Dune::FieldVector<T, dimworld> WorldCoords;
59 
61  typedef Dune::FieldVector<T, dim> LocalCoords;
62 
70  DUNE_DEPRECATED_MSG("Please use a std::function<FieldVector(FieldVector)> to prescribe non-default projections")
71  ContactMerge(const T allowedOverlap,
72  const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections,
73  const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections)
74  : overlap_(allowedOverlap)
75  {
76  if (domainDirections) {
77  domainDirections_ = [domainDirections](const WorldCoords& in) {
78  WorldCoords out;
79  domainDirections->evaluate(in,out);
80  return out;
81  };
82  }
83  if (targetDirections) {
84  targetDirections_ = [targetDirections](const WorldCoords& in) {
85  WorldCoords out;
86  targetDirections->evaluate(in,out);
87  return out;
88  };
89  }
90  }
91 
99  ContactMerge(const T allowedOverlap=T(0),
100  std::function<WorldCoords(WorldCoords)> domainDirections=nullptr,
101  std::function<WorldCoords(WorldCoords)> targetDirections=nullptr)
102  : domainDirections_(domainDirections), targetDirections_(targetDirections),
103  overlap_(allowedOverlap)
104  {}
105 
114  inline
115  void setSurfaceDirections(std::function<WorldCoords(WorldCoords)> domainDirections,
116  std::function<WorldCoords(WorldCoords)> targetDirections)
117  {
118  domainDirections_ = domainDirections;
119  targetDirections_ = targetDirections;
120  this->valid = false;
121  }
122 
131  DUNE_DEPRECATED_MSG("Please use a std::function<FieldVector(FieldVector)> to prescribe non-default projections")
132  inline
133  void setSurfaceDirections(const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections,
134  const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections)
135  {
136  domainDirections_ = [domainDirections](const WorldCoords& in) {
137  WorldCoords out;
138  domainDirections->evaluate(in,out);
139  return out;
140  };
141 
142  targetDirections_ = [targetDirections](const WorldCoords& in) {
143  WorldCoords out;
144  targetDirections->evaluate(in,out);
145  return out;
146  };
147  this->valid = false;
148  }
149 
151  void setOverlap(T overlap)
152  {
153  overlap_ = overlap;
154  }
155 
157  T getOverlap() const
158  {
159  return overlap_;
160  }
161 
165  void minNormalAngle(T angle)
166  {
167  using std::cos;
168  maxNormalProduct_ = cos(angle);
169  }
170 
174  T minNormalAngle() const
175  {
176  using std::acos;
177  return acos(maxNormalProduct_);
178  }
179 
180 protected:
181  typedef typename StandardMerge<T,dimworld-1,dimworld-1,dimworld>::RemoteSimplicialIntersection RemoteSimplicialIntersection;
182 
183 private:
187  std::function<WorldCoords(WorldCoords)> domainDirections_;
188  std::vector<WorldCoords> nodalDomainDirections_;
189 
198  std::function<WorldCoords(WorldCoords)> targetDirections_;
199  std::vector<WorldCoords> nodalTargetDirections_;
200 
202  T overlap_;
203 
207  T maxNormalProduct_ = T(-0.1);
208 
213  void computeIntersections(const Dune::GeometryType& grid1ElementType,
214  const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
215  std::bitset<(1<<dim)>& neighborIntersects1,
216  unsigned int grid1Index,
217  const Dune::GeometryType& grid2ElementType,
218  const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
219  std::bitset<(1<<dim)>& neighborIntersects2,
220  unsigned int grid2Index,
221  std::vector<RemoteSimplicialIntersection>& intersections);
222 
226 protected:
227  void build(const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
228  const std::vector<unsigned int>& grid1Elements,
229  const std::vector<Dune::GeometryType>& grid1ElementTypes,
230  const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
231  const std::vector<unsigned int>& grid2Elements,
232  const std::vector<Dune::GeometryType>& grid2ElementTypes)
233  {
234  std::cout<<"ContactMerge building grid!\n";
235  // setup the nodal direction vectors
236  setupNodalDirections(grid1Coords, grid1Elements, grid1ElementTypes,
237  grid2Coords, grid2Elements, grid2ElementTypes);
238 
239  Base::build(grid1Coords, grid1Elements, grid1ElementTypes,
240  grid2Coords, grid2Elements, grid2ElementTypes);
241 
242  }
243 
244 private:
245 
247  static LocalCoords localCornerCoords(int i, const Dune::GeometryType& gt)
248  {
249  const Dune::ReferenceElement<T,dim>& ref = Dune::ReferenceElements<T,dim>::general(gt);
250  return ref.position(i,dim);
251  }
252 
253 protected:
254 
256  void computeCyclicOrder(const std::vector<std::array<LocalCoords,2> >& polytopeCorners,
257  const LocalCoords& center, std::vector<int>& ordering) const;
258 
260  void 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 
268  void computeOuterNormalField(const std::vector<WorldCoords>& coords,
269  const std::vector<unsigned int>& elements,
270  const std::vector<Dune::GeometryType>& elementTypes,
271  std::vector<WorldCoords>& normals);
272 
274  void removeDoubles(std::vector<std::array<LocalCoords,2> >& polytopeCorners);
275 };
276 
277 } /* namespace GridGlue */
278 } /* namespace Dune */
279 
280 #include "contactmerge.cc"
281 
282 #endif // DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
Central component of the module implementing the coupling of two grids.
Common base class for many merger implementations: produce pairs of entities that may intersect.
Definition: gridglue.hh:30
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
Merge two codimension-1 surfaces that may be a positive distance apart.
Definition: contactmerge.hh:43
void build(const std::vector< Dune::FieldVector< T, dimworld > > &grid1Coords, const std::vector< unsigned int > &grid1Elements, const std::vector< Dune::GeometryType > &grid1ElementTypes, const std::vector< Dune::FieldVector< T, dimworld > > &grid2Coords, const std::vector< unsigned int > &grid2Elements, const std::vector< Dune::GeometryType > &grid2ElementTypes)
builds the merged grid
Definition: contactmerge.hh:227
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
void setOverlap(T overlap)
Set the allowed overlap of the surfaces.
Definition: contactmerge.hh:151
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 minNormalAngle(T angle)
set minimum angle in radians between normals at x and Φ(x)
Definition: contactmerge.hh:165
T ctype
the numeric type used in this interface
Definition: contactmerge.hh:55
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
T getOverlap() const
Get the allowed overlap of the surfaces.
Definition: contactmerge.hh:157
ContactMerge(const T allowedOverlap=T(0), std::function< WorldCoords(WorldCoords)> domainDirections=nullptr, std::function< WorldCoords(WorldCoords)> targetDirections=nullptr)
Construct merger given overlap and possible projection directions.
Definition: contactmerge.hh:99
void setSurfaceDirections(std::function< WorldCoords(WorldCoords)> domainDirections, std::function< WorldCoords(WorldCoords)> targetDirections)
Set surface direction functions.
Definition: contactmerge.hh:115
StandardMerge< T, dimworld-1, dimworld-1, dimworld >::RemoteSimplicialIntersection RemoteSimplicialIntersection
Definition: contactmerge.hh:181
T minNormalAngle() const
get minimum angle in radians between normals at x and Φ(x)
Definition: contactmerge.hh:174
Dune::FieldVector< T, dim > LocalCoords
the coordinate type used in this interface
Definition: contactmerge.hh:61
Common base class for many merger implementations: produce pairs of entities that may intersect.
Definition: standardmerge.hh:55
virtual void build(const std::vector< Dune::FieldVector< T, dimworld > > &grid1_Coords, const std::vector< unsigned int > &grid1_elements, const std::vector< Dune::GeometryType > &grid1_element_types, const std::vector< Dune::FieldVector< T, dimworld > > &grid2_coords, const std::vector< unsigned int > &grid2_elements, const std::vector< Dune::GeometryType > &grid2_element_types)
builds the merged grid
Definition: standardmerge.hh:465