dune-grid-glue  2.4-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 
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 
30 
36 template<int dimworld, typename T = double>
38 : public StandardMerge<T,dimworld-1,dimworld-1,dimworld>
39 {
40  enum {dim = dimworld-1};
41 
42  static_assert( dim==1 || dim==2,
43  "ContactMerge yet only handles the cases dim==1 and dim==2!");
44 
46 public:
47 
48  /* 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 */
49 
51  typedef T ctype;
52 
54  typedef Dune::FieldVector<T, dimworld> WorldCoords;
55 
57  typedef Dune::FieldVector<T, dim> LocalCoords;
58 
59  ContactMerge(const T allowedOverlap=T(0),
60  const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections = NULL,
61  const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections = NULL)
62  : domainDirections_(domainDirections), targetDirections_(targetDirections),
63  overlap_(allowedOverlap)
64  {}
65 
74  inline
75  void setSurfaceDirections(const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections,
76  const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections)
77  {
78  domainDirections_ = domainDirections;
79  targetDirections_ = targetDirections;
80  this->valid = false;
81  }
82 
84  void setOverlap(T overlap)
85  {
86  overlap_ = overlap;
87  }
88 
89 protected:
90  typedef typename StandardMerge<T,dimworld-1,dimworld-1,dimworld>::RemoteSimplicialIntersection RemoteSimplicialIntersection;
91 
92 private:
96  const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections_;
97  std::vector<WorldCoords> nodalDomainDirections_;
98 
107  const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections_;
108  std::vector<WorldCoords> nodalTargetDirections_;
109 
111  T overlap_;
112 
117  void computeIntersections(const Dune::GeometryType& grid1ElementType,
118  const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
119  std::bitset<(1<<dim)>& neighborIntersects1,
120  unsigned int grid1Index,
121  const Dune::GeometryType& grid2ElementType,
122  const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
123  std::bitset<(1<<dim)>& neighborIntersects2,
124  unsigned int grid2Index,
125  std::vector<RemoteSimplicialIntersection>& intersections);
126 
130 protected:
131  void build(const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
132  const std::vector<unsigned int>& grid1Elements,
133  const std::vector<Dune::GeometryType>& grid1ElementTypes,
134  const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
135  const std::vector<unsigned int>& grid2Elements,
136  const std::vector<Dune::GeometryType>& grid2ElementTypes)
137  {
138  std::cout<<"ContactMerge building grid!\n";
139  // setup the nodal direction vectors
140  setupNodalDirections(grid1Coords, grid1Elements, grid1ElementTypes,
141  grid2Coords, grid2Elements, grid2ElementTypes);
142 
143  Base::build(grid1Coords, grid1Elements, grid1ElementTypes,
144  grid2Coords, grid2Elements, grid2ElementTypes);
145 
146  }
147 
148 private:
149 
151  static LocalCoords localCornerCoords(int i, const Dune::GeometryType& gt)
152  {
153  const Dune::ReferenceElement<T,dim>& ref = Dune::ReferenceElements<T,dim>::general(gt);
154  return ref.position(i,dim);
155  }
156 
158  static int isCorner(const Dune::GeometryType& gt, const LocalCoords& local)
159  {
160  const Dune::ReferenceElement<T,dim>& ref = Dune::ReferenceElements<T,dim>::general(gt);
161 
162  for (int i=0; i<ref.size(dim); i++)
163  if ((ref.position(i,dim)-local).two_norm()<1e-8)
164  return i;
165  return -1;
166  }
167 
169  static std::vector<T> localToBarycentric(const Dune::GeometryType& gt, const LocalCoords& local)
170  {
171  std::vector<T> bar(Dune::ReferenceElements<T,dim>::general(gt).size(dim));
172  if (bar.size()<4) {
173  bar[0] = 1.0;
174  for (int i=0; i<dim; i++) {
175  bar[i+1] = local[i];
176  bar[0] -= local[i];
177  }
178  } else {
179  // quadrilateral case
180  for (int i=0; i<4; i++) {
181  bar[i] = 1;
182 
183  for (int j=0; j<dim; j++)
184  // if j-th bit of i is set multiply with in[j], else with 1-in[j]
185  bar[i] *= (i & (1<<j)) ? local[j] : 1-local[j];
186  }
187  }
188  return bar;
189  }
190 
191 protected:
193  static WorldCoords interpolate(const std::vector<WorldCoords>& p,
194  const Dune::GeometryType& gt, const LocalCoords& local)
195  {
196  // Compute barycentric coordinates
197  std::vector<T> bar = localToBarycentric(gt, local);
198 
199  WorldCoords global(0);
200  for (size_t i=0; i<p.size(); i++)
201  global.axpy(bar[i],p[i]);
202 
203  return global;
204  }
205 
217  bool isFeasibleProjection(const std::vector<WorldCoords>& elementCorners,
218  const std::vector<WorldCoords>& contactDirections,
219  const Dune::GeometryType& gt,
220  const WorldCoords& preImage,
221  const WorldCoords& preContactDirection,
222  const LocalCoords& localCoords)
223  {
224  T eps = 1e-6;
225 
226  // Compute the global coordinates of the projected point and the corr. contact direction
227  WorldCoords projPoint = interpolate(elementCorners, gt, localCoords);
228  WorldCoords projDirection = interpolate(contactDirections, gt, localCoords);
229 
230  WorldCoords segment = preImage - projPoint;
231  T distance = segment.two_norm();
232  segment /= distance;
233 
234  // The 'segment' vector always points from the projection to the preimage, thus
235  // if the contact direction of the preimage and the segment point in the same direction,
236  // then one of the faces must be on the backside of the body - invalid.
237  T angleSegPreDir = segment*preContactDirection;
238 
239  // The same holds if the segment and the contact direction at the projection.
240  T angleSegProjDir = segment*projDirection;
241 
242  // If both invalid conditions hold, then the intersection can be valid, if the two underlying
243  // bodies overlap a bit, which is often the case in the large deformation framework
244  // This intersection is then valid if the distance is smaller than an allowed overlap.
245 
246  // If both invalid-criterions hold
247  if (angleSegPreDir > -eps && angleSegProjDir < eps)
248  {
249  // then the projection is valid if the overlap is small
250  if (distance > overlap_)
251  return false;
252 
253  // if only one holds then the projection is invalid
254  } else if (angleSegPreDir > -eps || angleSegProjDir < eps)
255  return false;
256 
257  return true;
258  }
259 
261  void computeCyclicOrder(const std::vector<std::array<LocalCoords,2> >& polytopeCorners,
262  const LocalCoords& center, std::vector<int>& ordering) const;
263 
265  void setupNodalDirections(const std::vector<WorldCoords>& coords1,
266  const std::vector<unsigned int>& elements1,
267  const std::vector<Dune::GeometryType>& elementTypes1,
268  const std::vector<WorldCoords>& coords2,
269  const std::vector<unsigned int>& elements2,
270  const std::vector<Dune::GeometryType>& elementTypes2);
271 
273  void computeOuterNormalField(const std::vector<WorldCoords>& coords,
274  const std::vector<unsigned int>& elements,
275  const std::vector<Dune::GeometryType>& elementTypes,
276  std::vector<WorldCoords>& normals);
277 
279  void removeDoubles(std::vector<std::array<LocalCoords,2> >& polytopeCorners);
280 };
281 
282 #include "contactmerge.cc"
283 
284 #endif // DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
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:261
Central component of the module implementing the coupling of two grids.
void setSurfaceDirections(const Dune::VirtualFunction< WorldCoords, WorldCoords > *domainDirections, const Dune::VirtualFunction< WorldCoords, WorldCoords > *targetDirections)
Set surface direction functions.
Definition: contactmerge.hh:75
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:288
ContactMerge(const T allowedOverlap=T(0), const Dune::VirtualFunction< WorldCoords, WorldCoords > *domainDirections=NULL, const Dune::VirtualFunction< WorldCoords, WorldCoords > *targetDirections=NULL)
Definition: contactmerge.hh:59
void removeDoubles(std::vector< std::array< LocalCoords, 2 > > &polytopeCorners)
Remove all multiples.
Definition: contactmerge.cc:327
void setOverlap(T overlap)
Set the allowed overlap of the surfaces.
Definition: contactmerge.hh:84
Merge two codimension-1 surfaces that may be a positive distance apart.
Definition: contactmerge.hh:37
IteratorRange<...> intersections(const GridGlue<...> &glue, const Reverse<...> &reverse=!reversed)
Iterate over all intersections of a GridGlue.
Common base class for many merger implementations: produce pairs of entities that may intersect...
Definition: standardmerge.hh:52
static WorldCoords interpolate(const std::vector< WorldCoords > &p, const Dune::GeometryType &gt, const LocalCoords &local)
Compute global coordinates of a local point using linear interpolation of the corners.
Definition: contactmerge.hh:193
T ctype
the numeric type used in this interface
Definition: contactmerge.hh:51
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:206
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: contactmerge.hh:54
Base class for predicates selecting the part of a grid to be extracted.
StandardMerge< T, dimworld-1, dimworld-1, dimworld >::RemoteSimplicialIntersection RemoteSimplicialIntersection
Definition: contactmerge.hh:90
bool isFeasibleProjection(const std::vector< WorldCoords > &elementCorners, const std::vector< WorldCoords > &contactDirections, const Dune::GeometryType &gt, const WorldCoords &preImage, const WorldCoords &preContactDirection, const LocalCoords &localCoords)
Check if the projection of a point and an opposing face is valid.
Definition: contactmerge.hh:217
Dune::FieldVector< T, dim > LocalCoords
the coordinate type used in this interface
Definition: contactmerge.hh:57
Common base class for many merger implementations: produce pairs of entities that may intersect...
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 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:131
bool valid
Definition: standardmerge.hh:74