8 #ifndef DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
9 #define DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
19 #include <dune/common/fvector.hh>
20 #include <dune/common/function.hh>
21 #include <dune/common/exceptions.hh>
22 #include <dune/common/bitsetvector.hh>
24 #include <dune/grid/common/grid.hh>
35 template<
int dimworld,
typename T =
double>
39 enum {dim = dimworld-1};
41 static_assert( dim==1 || dim==2,
42 "ContactMerge yet only handles the cases dim==1 and dim==2!");
56 typedef Dune::FieldVector<T, dim> LocalCoords;
59 const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections = NULL,
60 const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections = NULL)
61 : domainDirections_(domainDirections), targetDirections_(targetDirections),
62 overlap_(allowedOverlap)
74 void setSurfaceDirections(
const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections,
75 const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections)
77 domainDirections_ = domainDirections;
78 targetDirections_ = targetDirections;
83 void setOverlap(T overlap)
89 typedef typename StandardMerge<T,dimworld-1,dimworld-1,dimworld>::RemoteSimplicialIntersection RemoteSimplicialIntersection;
95 const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections_;
96 std::vector<WorldCoords> nodalDomainDirections_;
106 const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections_;
107 std::vector<WorldCoords> nodalTargetDirections_;
117 const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
118 std::bitset<(1<<dim)>& neighborIntersects1,
119 unsigned int grid1Index,
120 const Dune::GeometryType& grid2ElementType,
121 const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
122 std::bitset<(1<<dim)>& neighborIntersects2,
123 unsigned int grid2Index,
124 std::vector<RemoteSimplicialIntersection>& intersections);
129 void build(
const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
130 const std::vector<unsigned int>& grid1Elements,
131 const std::vector<Dune::GeometryType>& grid1ElementTypes,
132 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
133 const std::vector<unsigned int>& grid2Elements,
134 const std::vector<Dune::GeometryType>& grid2ElementTypes)
136 std::cout<<
"ContactMerge building grid!\n";
138 setupNodalDirections(grid1Coords, grid1Elements, grid1ElementTypes,
139 grid2Coords, grid2Elements, grid2ElementTypes);
141 Base::build(grid1Coords, grid1Elements, grid1ElementTypes,
142 grid2Coords, grid2Elements, grid2ElementTypes);
149 static LocalCoords localCornerCoords(
int i,
const Dune::GeometryType& gt)
151 const Dune::ReferenceElement<T,dim>& ref = Dune::ReferenceElements<T,dim>::general(gt);
152 return ref.position(i,dim);
156 static int isCorner(
const Dune::GeometryType& gt,
const LocalCoords& local)
158 const Dune::ReferenceElement<T,dim>& ref = Dune::ReferenceElements<T,dim>::general(gt);
160 for (
int i=0; i<ref.size(dim); i++)
161 if ((ref.position(i,dim)-local).two_norm()<1e-8)
167 static std::vector<T> localToBarycentric(
const Dune::GeometryType& gt,
const LocalCoords& local)
169 std::vector<T> bar(Dune::ReferenceElements<T,dim>::general(gt).size(dim));
172 for (
int i=0; i<dim; i++) {
178 for (
int i=0; i<4; i++) {
181 for (
int j=0; j<dim; j++)
183 bar[i] *= (i & (1<<j)) ? local[j] : 1-local[j];
190 static WorldCoords interpolate(
const std::vector<WorldCoords>& p,
191 const Dune::GeometryType& gt,
const LocalCoords& local)
194 std::vector<T> bar = localToBarycentric(gt, local);
196 WorldCoords global(0);
197 for (
size_t i=0; i<p.size(); i++)
198 global.axpy(bar[i],p[i]);
204 void computeCyclicOrder(
const std::vector<Dune::array<LocalCoords,2> >& polytopeCorners,
205 const LocalCoords& center, std::vector<int>& ordering)
const;
208 void setupNodalDirections(
const std::vector<WorldCoords>& coords1,
209 const std::vector<unsigned int>& elements1,
210 const std::vector<Dune::GeometryType>& elementTypes1,
211 const std::vector<WorldCoords>& coords2,
212 const std::vector<unsigned int>& elements2,
213 const std::vector<Dune::GeometryType>& elementTypes2);
216 void computeOuterNormalField(
const std::vector<WorldCoords>& coords,
217 const std::vector<unsigned int>& elements,
218 const std::vector<Dune::GeometryType>& elementTypes,
219 std::vector<WorldCoords>& normals);
222 void removeDoubles(std::vector<Dune::array<LocalCoords,2> >& polytopeCorners);
227 #endif // DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
Central component of the module implementing the coupling of two grids.
Merge two codimension-1 surfaces that may be a positive distance apart.
Definition: contactmerge.hh:36
virtual void computeIntersections(const Dune::GeometryType &grid1ElementType, const std::vector< Dune::FieldVector< T, dimworld > > &grid1ElementCorners, std::bitset<(1<< grid1Dim)> &neighborIntersects1, unsigned int grid1Index, const Dune::GeometryType &grid2ElementType, const std::vector< Dune::FieldVector< T, dimworld > > &grid2ElementCorners, std::bitset<(1<< grid2Dim)> &neighborIntersects2, unsigned int grid2Index, std::vector< RemoteSimplicialIntersection > &intersections)=0
Compute the intersection between two overlapping elements.
T ctype
the numeric type used in this interface
Definition: standardmerge.hh:59
Common base class for many merger implementations: produce pairs of entities that may intersect...
Definition: standardmerge.hh:50
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
Common base class for many merger implementations: produce pairs of entities that may intersect...
Base class for predicates selecting the part of a grid to be extracted.
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: standardmerge.hh:68
bool valid
Definition: standardmerge.hh:72