Line |
Branch |
Exec |
Source |
1 |
|
|
// ______ ______ _ _ _____ ______ |
2 |
|
|
// | ____| ____| | (_)/ ____| | ____| |
3 |
|
|
// | |__ | |__ | | _| (___ ___| |__ |
4 |
|
|
// | __| | __| | | | |\___ \ / __| __| |
5 |
|
|
// | | | |____| |____| |____) | (__| |____ |
6 |
|
|
// |_| |______|______|_|_____/ \___|______| |
7 |
|
|
// Finite Elements for Life Sciences and Engineering |
8 |
|
|
// |
9 |
|
|
// License: LGL2.1 License |
10 |
|
|
// FELiScE default license: LICENSE in root folder |
11 |
|
|
// |
12 |
|
|
// Main authors: J. Foulon & J-F. Gerbeau & V. Martin |
13 |
|
|
// |
14 |
|
|
|
15 |
|
|
#ifndef SURFACEINTERPOLATOR_HPP |
16 |
|
|
#define SURFACEINTERPOLATOR_HPP |
17 |
|
|
|
18 |
|
|
// System includes |
19 |
|
|
|
20 |
|
|
// External includes |
21 |
|
|
|
22 |
|
|
// Project includes |
23 |
|
|
#include "Core/felisce.hpp" |
24 |
|
|
#include "Solver/linearProblem.hpp" |
25 |
|
|
#include "DegreeOfFreedom/dofBoundary.hpp" |
26 |
|
|
#include "PETScInterface/petscMatrix.hpp" |
27 |
|
|
#include "Geometry/geometricMeshRegion.hpp" |
28 |
|
|
|
29 |
|
|
namespace felisce |
30 |
|
|
{ |
31 |
|
|
class SurfaceInterpolator { |
32 |
|
|
public: |
33 |
|
|
typedef GeometricMeshRegion::ElementType ElementType; |
34 |
|
|
|
35 |
|
✗ |
SurfaceInterpolator()= default;; |
36 |
|
|
void initSurfaceInterpolator(LinearProblem* pb, DofBoundary* dofbd) { |
37 |
|
|
std::vector<DofBoundary*> tmp(1, dofbd); |
38 |
|
|
this->initSurfaceInterpolator(pb,tmp); |
39 |
|
|
} |
40 |
|
✗ |
void initSurfaceInterpolator(LinearProblem* pb, std::vector<DofBoundary*> dofbd) { |
41 |
|
✗ |
m_pb = pb; |
42 |
|
✗ |
m_dofBD = dofbd; |
43 |
|
✗ |
m_nrows.resize(dofbd.size()); |
44 |
|
✗ |
m_ncols.resize(dofbd.size()); |
45 |
|
✗ |
m_mesh = pb->mesh().get(); |
46 |
|
|
} |
47 |
|
|
void writeInterfaceFile(felInt iunknown, felInt numComp, std::string filename, std::string folder, std::size_t iBD=0 ) const; |
48 |
|
|
void readInterfaceFile(std::string filename, std::vector< std::pair< Point*, felInt > >& pointsAndCorrespondingSurfLabel, std::vector<int>& rowNumbering, std::string folder = FelisceParam::instance().resultDir); |
49 |
|
|
void buildInterpolator( PetscMatrix& interpolator, |
50 |
|
|
const std::vector< std::pair< Point*, felInt > >& pointsAndCorrespondingSurfLabel, |
51 |
|
|
const std::vector<int> rowNumbering, |
52 |
|
|
const Variable& var, felInt ivar, std::size_t iBD=0); |
53 |
|
|
void buildInterpolator( PetscMatrix& interpolator, |
54 |
|
|
const std::vector< std::pair< Point*, felInt > >& pointsAndCorrespondingSurfLabel, |
55 |
|
|
const std::vector<int> rowNumbering, |
56 |
|
|
const Variable& var, felInt ivar, std::size_t iBD, const std::vector<int> alternateLabels, std::string logFile); |
57 |
|
|
|
58 |
|
✗ |
int nrows(std::size_t iBD=0) {return m_nrows[iBD];} |
59 |
|
✗ |
int ncols(std::size_t iBD=0) {return m_ncols[iBD];} |
60 |
|
|
|
61 |
|
|
void initializeSurfaceInterpolator(); |
62 |
|
|
|
63 |
|
|
/** |
64 |
|
|
* @param thePoint [in] the point to be projected |
65 |
|
|
* @param surfaceLabel [in] the label of the surface where we are looking for the point |
66 |
|
|
* @param eltTypeReturn [out] if the point is found into one of the boundary element this is his type |
67 |
|
|
* @param elemPointReturn [out] if the point is found into one of the boundary element this are his points |
68 |
|
|
* @param projection [out] if the point is found into one of the boundary element this is its projection onto the element plane |
69 |
|
|
* @param ielSupportDof [out] ielSupportDof of the element |
70 |
|
|
* @param iel2ScoreByEdge [out] if the point is NOT found into one of the boundary element it contains a std::unordered_map from an ielSupportDofByType to the best edge distance, if the point is found on the element this values are incomplete and should not be used |
71 |
|
|
* @param distance [in/out] smaller distant from a single point, this value is updated only if it is better than the previous results |
72 |
|
|
* @param idElOfClosestPoint [in/out] the id of the el by typt closest point, this value is updated only if it is better than the previous results |
73 |
|
|
*/ |
74 |
|
|
int getClosestSurfaceElement( |
75 |
|
|
Point* thePoint, |
76 |
|
|
felInt surfaceLabel, ElementType& eltTypeReturn, |
77 |
|
|
std::vector<Point*>& elemPointReturn, |
78 |
|
|
Point& projection, felInt& ielSupportDof, |
79 |
|
|
std::map<int,double>& iel2ScoreByEdge, |
80 |
|
|
double& distance, |
81 |
|
|
int& idElOfClosestPoint |
82 |
|
|
); |
83 |
|
|
|
84 |
|
|
int lookIntoStar( |
85 |
|
|
int idClosestPoint, |
86 |
|
|
int surfaceLabel, |
87 |
|
|
Point* thePoint, |
88 |
|
|
ElementType& eltTypeReturn, |
89 |
|
|
std::vector<Point*>& elemPointReturn, |
90 |
|
|
Point& projection, |
91 |
|
|
felInt& ielSupportDof, |
92 |
|
|
felInt& depth, |
93 |
|
|
double& score, |
94 |
|
|
std::map<int,double>& iel2ScoreByEdge |
95 |
|
|
); |
96 |
|
|
|
97 |
|
|
void displayDataForSurfaceInterpolator(); |
98 |
|
|
|
99 |
|
|
ElementType getETandAssert2D(std::set<ElementType> setOfET); |
100 |
|
|
|
101 |
|
|
inline const bool& surfaceInterpolatorInitialized() const { |
102 |
|
|
return m_surfaceInterpolatorInitialized; |
103 |
|
|
} |
104 |
|
|
inline bool& surfaceInterpolatorInitialized() { |
105 |
|
|
return m_surfaceInterpolatorInitialized; |
106 |
|
|
} |
107 |
|
|
|
108 |
|
|
private: |
109 |
|
|
void checkInterpolatorStatistics(int foundInTria, int foundInTriaAltLab, int foundInEdge, int foundIntPoint, int pointNotFound, int nPoints ) const; |
110 |
|
|
|
111 |
|
|
|
112 |
|
|
// Is D inside ABC? |
113 |
|
|
bool isPointInTriangle(const double *A, const double *B,const double *C, const double *D, double& alpha, double& beta); |
114 |
|
|
bool isPointInEdgeAndIfYesHowFar( const double *p0,const double *p1,const double* thePoint,double& distance, double& t); |
115 |
|
|
|
116 |
|
|
void checkEdge(const double *p0,const double *p1,Point* thePoint,double& bestDistance); |
117 |
|
|
bool checkTriangle( const Point& projection, const std::vector<Point*>& elemPoint, double& alpha, double& beta, Point* thePoint, double& bestEdgeDistance ); |
118 |
|
|
bool checkQuadrangle( const Point& projection, const std::vector<Point*>& elemPoint, double& alpha, double& beta, Point* thePoint, double& bestEdgeDistance ); |
119 |
|
|
|
120 |
|
|
void projectOnElementPlane( const std::vector<Point*>& elemPoint, Point* thePoint, Point& projection ); |
121 |
|
|
void projectOnBestEdge( const std::vector<Point*>& elemPoint, Point* thePoint, Point& projection); |
122 |
|
|
void projectOnBestPoint( const std::vector<Point*>& elemPoint, Point* thePoint, Point& projection); |
123 |
|
|
void projectOnEdge(const double *p0,const double *p1,const double* thePoint, double& t); |
124 |
|
|
|
125 |
|
|
double computeScore(double alpha, double beta); |
126 |
|
|
void updateBestScore( double alpha, double beta, double& bestScore, double& alphabest, double& betabest); |
127 |
|
|
int closestPoint(const std::set< std::pair< Point*, felInt > >& setOfPointsAndIds, const Point* thePoint, Point*& closest, double& distance); |
128 |
|
|
|
129 |
|
|
|
130 |
|
|
LinearProblem* m_pb; |
131 |
|
|
std::vector<DofBoundary*> m_dofBD; |
132 |
|
|
|
133 |
|
|
std::vector<int> m_nrows; |
134 |
|
|
std::vector<int> m_ncols; |
135 |
|
|
|
136 |
|
|
GeometricMeshRegion* m_mesh; |
137 |
|
|
|
138 |
|
|
std::map< /*surfLabel*/ felInt, std::set< std::pair<Point*, /*idPoint*/ felInt> > > m_surfLab2SetOfPointsAndIds; |
139 |
|
|
std::map< /*surfLabel*/ felInt, std::map< /*IdPoint*/ felInt, /*star*/std::vector< /*idByType*/ felInt> > > m_surfLabelAndIdPoint2Star; |
140 |
|
|
|
141 |
|
|
bool m_surfaceInterpolatorInitialized = false; |
142 |
|
|
}; |
143 |
|
|
} |
144 |
|
|
#endif |
145 |
|
|
|