dune-grid-glue  2.9
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:
3 // SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-GPL-2.0-only-with-dune-grid-glue-exception
10 #ifndef DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
11 #define DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
12 
13 
14 #include <iostream>
15 #include <fstream>
16 #include <iomanip>
17 #include <vector>
18 #include <algorithm>
19 #include <limits>
20 #include <memory>
21 #include <functional>
22 
23 #include <dune/common/fvector.hh>
24 #include <dune/common/exceptions.hh>
25 #include <dune/common/bitsetvector.hh>
26 #include <dune/common/deprecated.hh>
27 
28 #include <dune/grid/common/grid.hh>
29 
32 
33 namespace Dune {
34 namespace GridGlue {
35 
41 template<int dimworld, typename T = double>
43 : public StandardMerge<T,dimworld-1,dimworld-1,dimworld>
44 {
45  static constexpr int dim = dimworld-1;
46 
47  static_assert( dim==1 || dim==2,
48  "ContactMerge yet only handles the cases dim==1 and dim==2!");
49 
51 public:
52 
53  /* 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 */
54 
56  typedef T ctype;
57 
59  typedef Dune::FieldVector<T, dimworld> WorldCoords;
60 
62  typedef Dune::FieldVector<T, dim> LocalCoords;
63 
73  ContactMerge(const T allowedOverlap=T(0),
74  std::function<WorldCoords(WorldCoords)> domainDirections=nullptr,
75  std::function<WorldCoords(WorldCoords)> targetDirections=nullptr,
77  : domainDirections_(domainDirections), targetDirections_(targetDirections),
78  overlap_(allowedOverlap), type_(type)
79  {}
80 
86  ContactMerge(const T allowedOverlap, ProjectionType type)
87  : overlap_(allowedOverlap),
88  type_(type)
89  {}
90 
99  inline
100  void setSurfaceDirections(std::function<WorldCoords(WorldCoords)> domainDirections,
101  std::function<WorldCoords(WorldCoords)> targetDirections)
102  {
103  domainDirections_ = domainDirections;
104  targetDirections_ = targetDirections;
105  this->valid = false;
106  }
107 
109  void setOverlap(T overlap)
110  {
111  overlap_ = overlap;
112  }
113 
115  T getOverlap() const
116  {
117  return overlap_;
118  }
119 
123  void minNormalAngle(T angle)
124  {
125  using std::cos;
126  maxNormalProduct_ = cos(angle);
127  }
128 
132  T minNormalAngle() const
133  {
134  using std::acos;
135  return acos(maxNormalProduct_);
136  }
137 
138 protected:
139  typedef typename StandardMerge<T,dimworld-1,dimworld-1,dimworld>::SimplicialIntersection SimplicialIntersection;
140 
141 private:
145  std::function<WorldCoords(WorldCoords)> domainDirections_;
146  std::vector<WorldCoords> nodalDomainDirections_;
147 
156  std::function<WorldCoords(WorldCoords)> targetDirections_;
157  std::vector<WorldCoords> nodalTargetDirections_;
158 
160  T overlap_;
161 
163  ProjectionType type_;
164 
168  T maxNormalProduct_ = T(-0.1);
169 
174  void computeIntersections(const Dune::GeometryType& grid1ElementType,
175  const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
176  std::bitset<(1<<dim)>& neighborIntersects1,
177  unsigned int grid1Index,
178  const Dune::GeometryType& grid2ElementType,
179  const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
180  std::bitset<(1<<dim)>& neighborIntersects2,
181  unsigned int grid2Index,
182  std::vector<SimplicialIntersection>& intersections) override;
183 
187 protected:
188  void build(const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
189  const std::vector<unsigned int>& grid1Elements,
190  const std::vector<Dune::GeometryType>& grid1ElementTypes,
191  const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
192  const std::vector<unsigned int>& grid2Elements,
193  const std::vector<Dune::GeometryType>& grid2ElementTypes) override
194  {
195  std::cout<<"ContactMerge building grid!\n";
196  // setup the nodal direction vectors
197  setupNodalDirections(grid1Coords, grid1Elements, grid1ElementTypes,
198  grid2Coords, grid2Elements, grid2ElementTypes);
199 
200  Base::build(grid1Coords, grid1Elements, grid1ElementTypes,
201  grid2Coords, grid2Elements, grid2ElementTypes);
202 
203  }
204 
205 private:
206 
208  static LocalCoords localCornerCoords(int i, const Dune::GeometryType& gt)
209  {
210  const auto& ref = Dune::ReferenceElements<T,dim>::general(gt);
211  return ref.position(i,dim);
212  }
213 
214 protected:
215 
217  void computeCyclicOrder(const std::vector<std::array<LocalCoords,2> >& polytopeCorners,
218  const LocalCoords& center, std::vector<int>& ordering) const;
219 
221  void setupNodalDirections(const std::vector<WorldCoords>& coords1,
222  const std::vector<unsigned int>& elements1,
223  const std::vector<Dune::GeometryType>& elementTypes1,
224  const std::vector<WorldCoords>& coords2,
225  const std::vector<unsigned int>& elements2,
226  const std::vector<Dune::GeometryType>& elementTypes2);
227 
229  void computeOuterNormalField(const std::vector<WorldCoords>& coords,
230  const std::vector<unsigned int>& elements,
231  const std::vector<Dune::GeometryType>& elementTypes,
232  std::vector<WorldCoords>& normals);
233 
235  void removeDoubles(std::vector<std::array<LocalCoords,2> >& polytopeCorners);
236 };
237 
238 } /* namespace GridGlue */
239 } /* namespace Dune */
240 
241 #include "contactmerge.cc"
242 
243 #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:37
Merge two codimension-1 surfaces that may be a positive distance apart.
Definition: contactmerge.hh:44
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:214
StandardMerge< T, dimworld-1, dimworld-1, dimworld >::SimplicialIntersection SimplicialIntersection
Definition: contactmerge.hh:139
void removeDoubles(std::vector< std::array< LocalCoords, 2 > > &polytopeCorners)
Remove all multiples.
Definition: contactmerge.cc:335
void setOverlap(T overlap)
Set the allowed overlap of the surfaces.
Definition: contactmerge.hh:109
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 minNormalAngle(T angle)
set minimum angle in radians between normals at x and Φ(x)
Definition: contactmerge.hh:123
T ctype
the numeric type used in this interface
Definition: contactmerge.hh:56
ProjectionType
Type of the projection, closest point or outer normal projection.
Definition: contactmerge.hh:65
@ CLOSEST_POINT
Definition: contactmerge.hh:65
@ OUTER_NORMAL
Definition: contactmerge.hh:65
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
T getOverlap() const
Get the allowed overlap of the surfaces.
Definition: contactmerge.hh:115
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:73
void setSurfaceDirections(std::function< WorldCoords(WorldCoords)> domainDirections, std::function< WorldCoords(WorldCoords)> targetDirections)
Set surface direction functions.
Definition: contactmerge.hh:100
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) override
builds the merged grid
Definition: contactmerge.hh:188
ContactMerge(const T allowedOverlap, ProjectionType type)
Construct merger given overlap and type of the projection.
Definition: contactmerge.hh:86
T minNormalAngle() const
get minimum angle in radians between normals at x and Φ(x)
Definition: contactmerge.hh:132
Dune::FieldVector< T, dim > LocalCoords
the coordinate type used in this interface
Definition: contactmerge.hh:62
Common base class for many merger implementations: produce pairs of entities that may intersect.
Definition: standardmerge.hh:58
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) override
builds the merged grid
Definition: standardmerge.hh:392