dune-grid-glue  2.5.0
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 
72  ContactMerge(const T allowedOverlap=T(0),
73  std::function<WorldCoords(WorldCoords)> domainDirections=nullptr,
74  std::function<WorldCoords(WorldCoords)> targetDirections=nullptr,
76  : domainDirections_(domainDirections), targetDirections_(targetDirections),
77  overlap_(allowedOverlap), type_(type)
78  {}
79 
85  ContactMerge(const T allowedOverlap, ProjectionType type)
86  : overlap_(allowedOverlap),
87  type_(type)
88  {}
89 
98  inline
99  void setSurfaceDirections(std::function<WorldCoords(WorldCoords)> domainDirections,
100  std::function<WorldCoords(WorldCoords)> targetDirections)
101  {
102  domainDirections_ = domainDirections;
103  targetDirections_ = targetDirections;
104  this->valid = false;
105  }
106 
108  void setOverlap(T overlap)
109  {
110  overlap_ = overlap;
111  }
112 
114  T getOverlap() const
115  {
116  return overlap_;
117  }
118 
122  void minNormalAngle(T angle)
123  {
124  using std::cos;
125  maxNormalProduct_ = cos(angle);
126  }
127 
131  T minNormalAngle() const
132  {
133  using std::acos;
134  return acos(maxNormalProduct_);
135  }
136 
137 protected:
138  typedef typename StandardMerge<T,dimworld-1,dimworld-1,dimworld>::RemoteSimplicialIntersection RemoteSimplicialIntersection;
139 
140 private:
144  std::function<WorldCoords(WorldCoords)> domainDirections_;
145  std::vector<WorldCoords> nodalDomainDirections_;
146 
155  std::function<WorldCoords(WorldCoords)> targetDirections_;
156  std::vector<WorldCoords> nodalTargetDirections_;
157 
159  T overlap_;
160 
162  ProjectionType type_;
163 
167  T maxNormalProduct_ = T(-0.1);
168 
173  void computeIntersections(const Dune::GeometryType& grid1ElementType,
174  const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
175  std::bitset<(1<<dim)>& neighborIntersects1,
176  unsigned int grid1Index,
177  const Dune::GeometryType& grid2ElementType,
178  const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
179  std::bitset<(1<<dim)>& neighborIntersects2,
180  unsigned int grid2Index,
181  std::vector<RemoteSimplicialIntersection>& intersections);
182 
186 protected:
187  void build(const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
188  const std::vector<unsigned int>& grid1Elements,
189  const std::vector<Dune::GeometryType>& grid1ElementTypes,
190  const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
191  const std::vector<unsigned int>& grid2Elements,
192  const std::vector<Dune::GeometryType>& grid2ElementTypes)
193  {
194  std::cout<<"ContactMerge building grid!\n";
195  // setup the nodal direction vectors
196  setupNodalDirections(grid1Coords, grid1Elements, grid1ElementTypes,
197  grid2Coords, grid2Elements, grid2ElementTypes);
198 
199  Base::build(grid1Coords, grid1Elements, grid1ElementTypes,
200  grid2Coords, grid2Elements, grid2ElementTypes);
201 
202  }
203 
204 private:
205 
207  static LocalCoords localCornerCoords(int i, const Dune::GeometryType& gt)
208  {
209  const Dune::ReferenceElement<T,dim>& ref = Dune::ReferenceElements<T,dim>::general(gt);
210  return ref.position(i,dim);
211  }
212 
213 protected:
214 
216  void computeCyclicOrder(const std::vector<std::array<LocalCoords,2> >& polytopeCorners,
217  const LocalCoords& center, std::vector<int>& ordering) const;
218 
220  void setupNodalDirections(const std::vector<WorldCoords>& coords1,
221  const std::vector<unsigned int>& elements1,
222  const std::vector<Dune::GeometryType>& elementTypes1,
223  const std::vector<WorldCoords>& coords2,
224  const std::vector<unsigned int>& elements2,
225  const std::vector<Dune::GeometryType>& elementTypes2);
226 
228  void computeOuterNormalField(const std::vector<WorldCoords>& coords,
229  const std::vector<unsigned int>& elements,
230  const std::vector<Dune::GeometryType>& elementTypes,
231  std::vector<WorldCoords>& normals);
232 
234  void removeDoubles(std::vector<std::array<LocalCoords,2> >& polytopeCorners);
235 };
236 
237 } /* namespace GridGlue */
238 } /* namespace Dune */
239 
240 #include "contactmerge.cc"
241 
242 #endif // DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
Common base class for many merger implementations: produce pairs of entities that may intersect...
Definition: standardmerge.hh:53
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
ContactMerge(const T allowedOverlap, ProjectionType type)
Construct merger given overlap and type of the projection.
Definition: contactmerge.hh:85
StandardMerge< T, dimworld-1, dimworld-1, dimworld >::RemoteSimplicialIntersection RemoteSimplicialIntersection
Definition: contactmerge.hh:138
T ctype
the numeric type used in this interface
Definition: contactmerge.hh:55
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:271
Definition: gridglue.hh:31
ContactMerge(const T allowedOverlap=T(0), std::function< WorldCoords(WorldCoords)> domainDirections=nullptr, std::function< WorldCoords(WorldCoords)> targetDirections=nullptr, ProjectionType type=OUTER_NORMAL)
Construct merger given overlap and possible projection directions.
Definition: contactmerge.hh:72
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
T getOverlap() const
Get the allowed overlap of the surfaces.
Definition: contactmerge.hh:114
Merge two codimension-1 surfaces that may be a positive distance apart.
Definition: contactmerge.hh:41
ProjectionType
Type of the projection, closest point or outer normal projection.
Definition: contactmerge.hh:64
void minNormalAngle(T angle)
set minimum angle in radians between normals at x and Φ(x)
Definition: contactmerge.hh:122
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: contactmerge.hh:58
Dune::FieldVector< T, dim > LocalCoords
the coordinate type used in this interface
Definition: contactmerge.hh:61
void setOverlap(T overlap)
Set the allowed overlap of the surfaces.
Definition: contactmerge.hh:108
Definition: contactmerge.hh:64
Common base class for many merger implementations: produce pairs of entities that may intersect...
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:298
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:216
Definition: contactmerge.hh:64
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:187
void setSurfaceDirections(std::function< WorldCoords(WorldCoords)> domainDirections, std::function< WorldCoords(WorldCoords)> targetDirections)
Set surface direction functions.
Definition: contactmerge.hh:99
T minNormalAngle() const
get minimum angle in radians between normals at x and Φ(x)
Definition: contactmerge.hh:131
Central component of the module implementing the coupling of two grids.
void removeDoubles(std::vector< std::array< LocalCoords, 2 > > &polytopeCorners)
Remove all multiples.
Definition: contactmerge.cc:337