dune-grid-glue  2.4.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 
20 #include <dune/common/fvector.hh>
21 #include <dune/common/function.hh>
22 #include <dune/common/exceptions.hh>
23 #include <dune/common/bitsetvector.hh>
24 
25 #include <dune/grid/common/grid.hh>
26 
29 
30 namespace Dune {
31 namespace GridGlue {
32 
38 template<int dimworld, typename T = double>
40 : public StandardMerge<T,dimworld-1,dimworld-1,dimworld>
41 {
42  enum {dim = dimworld-1};
43 
44  static_assert( dim==1 || dim==2,
45  "ContactMerge yet only handles the cases dim==1 and dim==2!");
46 
48 public:
49 
50  /* 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 */
51 
53  typedef T ctype;
54 
56  typedef Dune::FieldVector<T, dimworld> WorldCoords;
57 
59  typedef Dune::FieldVector<T, dim> LocalCoords;
60 
61  ContactMerge(const T allowedOverlap=T(0),
62  const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections = NULL,
63  const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections = NULL)
64  : domainDirections_(domainDirections), targetDirections_(targetDirections),
65  overlap_(allowedOverlap)
66  {}
67 
76  inline
77  void setSurfaceDirections(const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections,
78  const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections)
79  {
80  domainDirections_ = domainDirections;
81  targetDirections_ = targetDirections;
82  this->valid = false;
83  }
84 
86  void setOverlap(T overlap)
87  {
88  overlap_ = overlap;
89  }
90 
92  T getOverlap() const
93  {
94  return overlap_;
95  }
96 
100  void minNormalAngle(T angle)
101  {
102  using std::cos;
103  maxNormalProduct_ = cos(angle);
104  }
105 
109  T minNormalAngle() const
110  {
111  using std::acos;
112  return acos(maxNormalProduct_);
113  }
114 
115 protected:
116  typedef typename StandardMerge<T,dimworld-1,dimworld-1,dimworld>::RemoteSimplicialIntersection RemoteSimplicialIntersection;
117 
118 private:
122  const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections_;
123  std::vector<WorldCoords> nodalDomainDirections_;
124 
133  const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections_;
134  std::vector<WorldCoords> nodalTargetDirections_;
135 
137  T overlap_;
138 
142  T maxNormalProduct_ = T(-0.1);
143 
148  void computeIntersections(const Dune::GeometryType& grid1ElementType,
149  const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
150  std::bitset<(1<<dim)>& neighborIntersects1,
151  unsigned int grid1Index,
152  const Dune::GeometryType& grid2ElementType,
153  const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
154  std::bitset<(1<<dim)>& neighborIntersects2,
155  unsigned int grid2Index,
156  std::vector<RemoteSimplicialIntersection>& intersections);
157 
161 protected:
162  void build(const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
163  const std::vector<unsigned int>& grid1Elements,
164  const std::vector<Dune::GeometryType>& grid1ElementTypes,
165  const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
166  const std::vector<unsigned int>& grid2Elements,
167  const std::vector<Dune::GeometryType>& grid2ElementTypes)
168  {
169  std::cout<<"ContactMerge building grid!\n";
170  // setup the nodal direction vectors
171  setupNodalDirections(grid1Coords, grid1Elements, grid1ElementTypes,
172  grid2Coords, grid2Elements, grid2ElementTypes);
173 
174  Base::build(grid1Coords, grid1Elements, grid1ElementTypes,
175  grid2Coords, grid2Elements, grid2ElementTypes);
176 
177  }
178 
179 private:
180 
182  static LocalCoords localCornerCoords(int i, const Dune::GeometryType& gt)
183  {
184  const Dune::ReferenceElement<T,dim>& ref = Dune::ReferenceElements<T,dim>::general(gt);
185  return ref.position(i,dim);
186  }
187 
188 protected:
189 
191  void computeCyclicOrder(const std::vector<std::array<LocalCoords,2> >& polytopeCorners,
192  const LocalCoords& center, std::vector<int>& ordering) const;
193 
195  void setupNodalDirections(const std::vector<WorldCoords>& coords1,
196  const std::vector<unsigned int>& elements1,
197  const std::vector<Dune::GeometryType>& elementTypes1,
198  const std::vector<WorldCoords>& coords2,
199  const std::vector<unsigned int>& elements2,
200  const std::vector<Dune::GeometryType>& elementTypes2);
201 
203  void computeOuterNormalField(const std::vector<WorldCoords>& coords,
204  const std::vector<unsigned int>& elements,
205  const std::vector<Dune::GeometryType>& elementTypes,
206  std::vector<WorldCoords>& normals);
207 
209  void removeDoubles(std::vector<std::array<LocalCoords,2> >& polytopeCorners);
210 };
211 
212 } /* namespace GridGlue */
213 } /* namespace Dune */
214 
215 #include "contactmerge.cc"
216 
217 #endif // DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
void minNormalAngle(T angle)
set minimum angle in radians between normals at x and Φ(x)
Definition: contactmerge.hh:100
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: contactmerge.hh:56
Definition: gridglue.hh:33
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:279
void setSurfaceDirections(const Dune::VirtualFunction< WorldCoords, WorldCoords > *domainDirections, const Dune::VirtualFunction< WorldCoords, WorldCoords > *targetDirections)
Set surface direction functions.
Definition: contactmerge.hh:77
StandardMerge< T, dimworld-1, dimworld-1, dimworld >::RemoteSimplicialIntersection RemoteSimplicialIntersection
Definition: contactmerge.hh:116
T getOverlap() const
Get the allowed overlap of the surfaces.
Definition: contactmerge.hh:92
Central component of the module implementing the coupling of two grids.
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
void removeDoubles(std::vector< std::array< LocalCoords, 2 > > &polytopeCorners)
Remove all multiples.
Definition: contactmerge.cc:345
T minNormalAngle() const
get minimum angle in radians between normals at x and Φ(x)
Definition: contactmerge.hh:109
Dune::FieldVector< T, dim > LocalCoords
the coordinate type used in this interface
Definition: contactmerge.hh:59
T ctype
the numeric type used in this interface
Definition: contactmerge.hh:53
Common base class for many merger implementations: produce pairs of entities that may intersect...
Definition: standardmerge.hh:54
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
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:224
Common base class for many merger implementations: produce pairs of entities that may intersect...
Merge two codimension-1 surfaces that may be a positive distance apart.
Definition: contactmerge.hh:39
ContactMerge(const T allowedOverlap=T(0), const Dune::VirtualFunction< WorldCoords, WorldCoords > *domainDirections=NULL, const Dune::VirtualFunction< WorldCoords, WorldCoords > *targetDirections=NULL)
Definition: contactmerge.hh:61
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:162
void setOverlap(T overlap)
Set the allowed overlap of the surfaces.
Definition: contactmerge.hh:86
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:306