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
32namespace Dune {
33namespace GridGlue {
34
40template<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
50public:
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
175 {
176 using std::acos;
177 return acos(maxNormalProduct_);
178 }
179
180protected:
181 typedef typename StandardMerge<T,dimworld-1,dimworld-1,dimworld>::RemoteSimplicialIntersection RemoteSimplicialIntersection;
182
183private:
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
226protected:
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
244private:
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
253protected:
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)
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)
Definition standardmerge.hh:465
bool valid
Definition standardmerge.hh:75