8 #ifndef DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
9 #define DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
20 #include <dune/common/fvector.hh>
21 #include <dune/common/function.hh>
22 #include <dune/common/exceptions.hh>
23 #include <dune/common/bitsetvector.hh>
25 #include <dune/grid/common/grid.hh>
36 template<
int dimworld,
typename T =
double>
40 enum {dim = dimworld-1};
42 static_assert( dim==1 || dim==2,
43 "ContactMerge yet only handles the cases dim==1 and dim==2!");
60 const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections = NULL,
61 const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections = NULL)
62 : domainDirections_(domainDirections), targetDirections_(targetDirections),
63 overlap_(allowedOverlap)
76 const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections)
78 domainDirections_ = domainDirections;
79 targetDirections_ = targetDirections;
96 const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections_;
97 std::vector<WorldCoords> nodalDomainDirections_;
107 const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections_;
108 std::vector<WorldCoords> nodalTargetDirections_;
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,
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)
138 std::cout<<
"ContactMerge building grid!\n";
141 grid2Coords, grid2Elements, grid2ElementTypes);
143 Base::build(grid1Coords, grid1Elements, grid1ElementTypes,
144 grid2Coords, grid2Elements, grid2ElementTypes);
151 static LocalCoords localCornerCoords(
int i,
const Dune::GeometryType& gt)
153 const Dune::ReferenceElement<T,dim>& ref = Dune::ReferenceElements<T,dim>::general(gt);
154 return ref.position(i,dim);
158 static int isCorner(
const Dune::GeometryType& gt,
const LocalCoords& local)
160 const Dune::ReferenceElement<T,dim>& ref = Dune::ReferenceElements<T,dim>::general(gt);
162 for (
int i=0; i<ref.size(dim); i++)
163 if ((ref.position(i,dim)-local).two_norm()<1e-8)
169 static std::vector<T> localToBarycentric(
const Dune::GeometryType& gt,
const LocalCoords& local)
171 std::vector<T> bar(Dune::ReferenceElements<T,dim>::general(gt).size(dim));
174 for (
int i=0; i<dim; i++) {
180 for (
int i=0; i<4; i++) {
183 for (
int j=0; j<dim; j++)
185 bar[i] *= (i & (1<<j)) ? local[j] : 1-local[j];
194 const Dune::GeometryType& gt,
const LocalCoords& local)
197 std::vector<T> bar = localToBarycentric(gt, local);
199 WorldCoords global(0);
200 for (
size_t i=0; i<p.size(); i++)
201 global.axpy(bar[i],p[i]);
218 const std::vector<WorldCoords>& contactDirections,
219 const Dune::GeometryType& gt,
220 const WorldCoords& preImage,
221 const WorldCoords& preContactDirection,
222 const LocalCoords& localCoords)
227 WorldCoords projPoint =
interpolate(elementCorners, gt, localCoords);
228 WorldCoords projDirection =
interpolate(contactDirections, gt, localCoords);
230 WorldCoords segment = preImage - projPoint;
231 T distance = segment.two_norm();
237 T angleSegPreDir = segment*preContactDirection;
240 T angleSegProjDir = segment*projDirection;
247 if (angleSegPreDir > -eps && angleSegProjDir < eps)
250 if (distance > overlap_)
254 }
else if (angleSegPreDir > -eps || angleSegProjDir < eps)
261 void computeCyclicOrder(
const std::vector<std::array<LocalCoords,2> >& polytopeCorners,
262 const LocalCoords& center, std::vector<int>& ordering)
const;
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);
274 const std::vector<unsigned int>& elements,
275 const std::vector<Dune::GeometryType>& elementTypes,
276 std::vector<WorldCoords>& normals);
279 void removeDoubles(std::vector<std::array<LocalCoords,2> >& polytopeCorners);
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 >, 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 ¢er, 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 >, 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