| 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 |
|
|
// System includes |
| 16 |
|
|
|
| 17 |
|
|
// External includes |
| 18 |
|
|
|
| 19 |
|
|
// Project includes |
| 20 |
|
|
#include "Geometry/surfaceInterpolator.hpp" |
| 21 |
|
|
#include "Geometry/point.hpp" |
| 22 |
|
|
#include "Geometry/geo_utilities.hpp" |
| 23 |
|
|
#include "Tools/math_utilities.hpp" |
| 24 |
|
|
|
| 25 |
|
|
namespace felisce { |
| 26 |
|
|
//default parameter for alternate labels: an empty std::vector. |
| 27 |
|
✗ |
void SurfaceInterpolator::buildInterpolator( PetscMatrix& interpolator, |
| 28 |
|
|
const std::vector< std::pair< Point*, felInt > >& pointsAndCorrespondingSurfLabel, |
| 29 |
|
|
const std::vector<int> rowNumbering, |
| 30 |
|
|
const Variable& var, felInt ivar, std::size_t iBD ) { |
| 31 |
|
✗ |
std::vector<int> alternativeLabels; |
| 32 |
|
✗ |
std::string logFile="logFile"; |
| 33 |
|
✗ |
this->buildInterpolator( interpolator, |
| 34 |
|
|
pointsAndCorrespondingSurfLabel, |
| 35 |
|
|
rowNumbering, |
| 36 |
|
|
var,ivar,iBD, |
| 37 |
|
|
alternativeLabels, |
| 38 |
|
|
logFile); |
| 39 |
|
|
} |
| 40 |
|
|
|
| 41 |
|
|
namespace { |
| 42 |
|
✗ |
void printLog(std::ostream& log, int id, Point* thePoint, Point& pointProjection, double distance) { |
| 43 |
|
✗ |
log<<id<<'\t'; |
| 44 |
|
✗ |
thePoint->printWithTab(log); |
| 45 |
|
✗ |
log<<'\t'; |
| 46 |
|
✗ |
pointProjection.printWithTab(log); |
| 47 |
|
✗ |
log<<'\t'<<distance<<std::endl; |
| 48 |
|
|
} |
| 49 |
|
|
} |
| 50 |
|
|
|
| 51 |
|
|
/* |
| 52 |
|
|
* Let f_1 be a finite-element function defined on the surface on the current domain (1). |
| 53 |
|
|
* f_1(x) = sum_j phi_{1,j}(x) f_{1,j} |
| 54 |
|
|
* where phi_{1,j} are the basis function of the surface. |
| 55 |
|
|
* |
| 56 |
|
|
* We want to pass this information to another domain (2) that share a surface, but has not conforming mesh. |
| 57 |
|
|
* |
| 58 |
|
|
* f_2(x) = sum_i phi_{2,i}(x) f_{2,i} |
| 59 |
|
|
* where f_{2,i} = f_1(x_i), x_i are the points coming from the second domain. |
| 60 |
|
|
* f_{2,i}=f_1(x_i) = sum_j phi_{1,j}(x_i) f_{1,j}. |
| 61 |
|
|
* |
| 62 |
|
|
* From this we get the matrix representation of the relationship between the two finite-element representation. |
| 63 |
|
|
* f_{2,i}= A_{i,j} f_{1,j} |
| 64 |
|
|
* where A_{i,j} = phi_{1,j}(x_i) is the interpolator that we want to build in this function! |
| 65 |
|
|
* |
| 66 |
|
|
* A will be parallel with respect to the column! |
| 67 |
|
|
*/ |
| 68 |
|
✗ |
void SurfaceInterpolator::buildInterpolator( PetscMatrix& interpolator, |
| 69 |
|
|
const std::vector< std::pair< Point*, felInt > >& pointsAndCorrespondingSurfLabel, |
| 70 |
|
|
const std::vector<int> rowNumbering, |
| 71 |
|
|
const Variable& var, felInt ivar, std::size_t iBD, |
| 72 |
|
|
const std::vector<int> alternativeLabels, std::string logFile) |
| 73 |
|
|
{ |
| 74 |
|
✗ |
if ( MpiInfo::rankProc() == 0 ) { |
| 75 |
|
✗ |
std::ofstream log; |
| 76 |
|
✗ |
log.open(logFile.c_str()); |
| 77 |
|
|
// First we extract the information. |
| 78 |
|
|
// the points in pointsAndCorrespondingSurfLabel are the x_i. |
| 79 |
|
|
|
| 80 |
|
|
// Warning: if the x_i belongs to more than one label, it does not matter. |
| 81 |
|
|
// in fact, we assume that the edges of the surface are geometrically conforming. Which means that they can have, as 1D manifolds, different finite-element |
| 82 |
|
|
// representation, but they overlap. |
| 83 |
|
|
|
| 84 |
|
|
// Number of points coming from the other mesh. |
| 85 |
|
✗ |
felInt nPoints = pointsAndCorrespondingSurfLabel.size(); // physical points that can support more than one dofs. Think about std::vector variables. |
| 86 |
|
✗ |
if ( nPoints == 0 ) |
| 87 |
|
✗ |
FEL_ERROR("Zero points in surface interpolator."); |
| 88 |
|
|
// Number of dofs of the other mesh, which is also the number of rows |
| 89 |
|
✗ |
m_nrows[iBD] = rowNumbering.size(); |
| 90 |
|
|
// We assume that all the points support the same number of dofs |
| 91 |
|
✗ |
FEL_ASSERT( m_nrows[iBD] % nPoints == 0 ); |
| 92 |
|
|
// We compute this number |
| 93 |
|
✗ |
felInt nComp = m_nrows[iBD]/nPoints; // 1 if we are passing scalars, 2 or 3 if we are passing vectors |
| 94 |
|
|
|
| 95 |
|
|
// We first save the interpolator in this data structure |
| 96 |
|
|
// We assume that rowNumbering contains all and only the numbers from 0 to m_nrows-1; |
| 97 |
|
✗ |
std::vector< std::vector < std::pair < /*column numbering*/ felInt, /*value*/ double > > > temporaryInterpolator; |
| 98 |
|
✗ |
temporaryInterpolator.resize(m_nrows[iBD]); |
| 99 |
|
|
|
| 100 |
|
|
// eltType will be used in the loop to distinguish between triangles and quadrangles |
| 101 |
|
|
// previousEltType will be used to store the previous value of eltType and check if it has changed |
| 102 |
|
|
// It is initialized to Tetra4 because Tetra4 can not be an element of the boundary, the idea is to std::set this element type to "NULL" |
| 103 |
|
✗ |
ElementType eltType, previousEltType=GeometricMeshRegion::Tetra4; |
| 104 |
|
|
// Vector to store vertex of the current element |
| 105 |
|
✗ |
std::vector< Point*> elemPoint; |
| 106 |
|
|
// Curvilinear finite element on the boundary, to access basis functions |
| 107 |
|
✗ |
CurvilinearFiniteElement* currentFe=nullptr; |
| 108 |
|
|
// Number of degrees of freedom of the current element |
| 109 |
|
✗ |
felInt nDofEle = 0; |
| 110 |
|
|
// Global index of the element |
| 111 |
|
|
felInt ielSupportDof; |
| 112 |
|
|
// Id of the column |
| 113 |
|
|
felInt idColumn; |
| 114 |
|
|
// Auxiliary vectors to temporary store the values of the matrix |
| 115 |
|
✗ |
std::vector<double> values; |
| 116 |
|
|
// projected point |
| 117 |
|
✗ |
Point pointProjection(0.0); |
| 118 |
|
|
|
| 119 |
|
✗ |
int idForLog = 0; //0 inside elem, 1 inside elem alternative lab, 2 edge, 3 point 4 not found |
| 120 |
|
|
|
| 121 |
|
✗ |
const int iMesh = var.idMesh(); |
| 122 |
|
|
|
| 123 |
|
|
// For each x_i |
| 124 |
|
✗ |
int pointNotFound(0); |
| 125 |
|
✗ |
int foundInsideAnElement(0), foundInsideAnElementAltLab(0), foundInEdge(0),foundInPoint(0); |
| 126 |
|
✗ |
for ( std::size_t cPoint(0); cPoint< (std::size_t) nPoints; ++cPoint) { |
| 127 |
|
|
|
| 128 |
|
|
// 1) IDENTIFY THE ELEMENT, ON THIS MESH, CLOSEST TO THE POINT. |
| 129 |
|
✗ |
std::map<int,double> iel2ScoreByEdge; |
| 130 |
|
✗ |
double pointDistance(1.e20); |
| 131 |
|
✗ |
int idOfElOfClosestPoint,ielSupportDofByType(0); |
| 132 |
|
|
int notFound; |
| 133 |
|
✗ |
Point* thePoint = pointsAndCorrespondingSurfLabel[cPoint].first; |
| 134 |
|
|
// 1a) check if the projection can be done on a triangle on the given label. |
| 135 |
|
|
// at the same time, best point and best edges are computed for the current label |
| 136 |
|
✗ |
notFound = getClosestSurfaceElement(thePoint, /* [in] the point*/ |
| 137 |
|
✗ |
pointsAndCorrespondingSurfLabel[cPoint].second,/* [in] the surface label where to look for*/ |
| 138 |
|
|
eltType, /* [out] type of the element */ |
| 139 |
|
|
elemPoint, /* [out] verteces of the element */ |
| 140 |
|
|
pointProjection, /* [out] the projection of the given point on the element*/ |
| 141 |
|
|
ielSupportDof, /* [out] the id of the element*/ |
| 142 |
|
|
iel2ScoreByEdge, /* [out] updated std::unordered_map from ielByType to the best edge score*/ |
| 143 |
|
|
pointDistance, /* [out] updated best distance from a point of the mesh*/ |
| 144 |
|
|
idOfElOfClosestPoint); /* [out] updated id of the closest point of the mesh*/ |
| 145 |
|
✗ |
foundInsideAnElement += (1-notFound); |
| 146 |
|
✗ |
if ( ! notFound ) { |
| 147 |
|
✗ |
idForLog=0; |
| 148 |
|
|
} |
| 149 |
|
✗ |
if ( notFound ) { |
| 150 |
|
|
// 1b) check if the projection can be done on a triangle on the other possible labels. |
| 151 |
|
|
// at the same time, best point and best edges are possibily updated for the current label |
| 152 |
|
✗ |
for ( std::size_t iL(0); iL<alternativeLabels.size() && notFound; ++iL ) { |
| 153 |
|
✗ |
if ( alternativeLabels[iL] != pointsAndCorrespondingSurfLabel[cPoint].second ) { |
| 154 |
|
✗ |
notFound = getClosestSurfaceElement(thePoint, /* [in] the point*/ |
| 155 |
|
✗ |
alternativeLabels[iL],/* [in] the surface label where to look for*/ |
| 156 |
|
|
eltType, /* [out] type of the element */ |
| 157 |
|
|
elemPoint, /* [out] verteces of the element */ |
| 158 |
|
|
pointProjection,/* [out] the projection of the given point on the element*/ |
| 159 |
|
|
ielSupportDof,iel2ScoreByEdge,pointDistance,idOfElOfClosestPoint); /* [out] the id of the element*/ |
| 160 |
|
✗ |
foundInsideAnElementAltLab+= (1-notFound); |
| 161 |
|
✗ |
if ( ! notFound ) { |
| 162 |
|
✗ |
idForLog=1; |
| 163 |
|
|
} |
| 164 |
|
|
} |
| 165 |
|
|
} |
| 166 |
|
|
// 1c) It was not possible to project it on a triangle, we have to project it either on an edge or on a point |
| 167 |
|
✗ |
if ( notFound ) { |
| 168 |
|
|
// 1d) we look for the best edge |
| 169 |
|
✗ |
double best=100; |
| 170 |
|
✗ |
for ( auto it = iel2ScoreByEdge.begin(); it != iel2ScoreByEdge.end(); ++it ) { |
| 171 |
|
✗ |
if ( it->second < best ) { |
| 172 |
|
✗ |
ielSupportDofByType = it->first; |
| 173 |
|
✗ |
best = it->second; |
| 174 |
|
|
} |
| 175 |
|
|
} |
| 176 |
|
|
// 1e) we check if it is better to project on a point instead |
| 177 |
|
✗ |
if ( best > pointDistance ) { |
| 178 |
|
✗ |
foundInPoint++; |
| 179 |
|
✗ |
ielSupportDofByType = idOfElOfClosestPoint; |
| 180 |
|
✗ |
idForLog=3; |
| 181 |
|
|
} else { |
| 182 |
|
✗ |
foundInEdge++; |
| 183 |
|
✗ |
idForLog=2; |
| 184 |
|
|
} |
| 185 |
|
|
// We hope the eltType of the element is the same on all the alternative labels..todo!! |
| 186 |
|
✗ |
ElementType eltType2 = getETandAssert2D(m_pb->mesh(iMesh)->intRefToEnum().at(pointsAndCorrespondingSurfLabel[cPoint].second)); |
| 187 |
|
✗ |
ielSupportDof=0;//very important |
| 188 |
|
✗ |
m_pb->mesh(iMesh)->getIdElemFromTypeElemAndIdByType(eltType,ielSupportDofByType,ielSupportDof); |
| 189 |
|
✗ |
std::vector<int> elemIdPoint( GeometricMeshRegion::m_numPointsPerElt[eltType2] ); |
| 190 |
|
✗ |
elemPoint.resize( GeometricMeshRegion::m_numPointsPerElt[eltType2] ); |
| 191 |
|
|
|
| 192 |
|
✗ |
m_pb->mesh(iMesh)->getOneElement(eltType2, ielSupportDofByType, elemIdPoint, 0); |
| 193 |
|
|
|
| 194 |
|
✗ |
for (int iPoint = 0; iPoint < GeometricMeshRegion::m_numPointsPerElt[eltType2]; iPoint++) { |
| 195 |
|
✗ |
elemPoint[iPoint] = &m_pb->mesh(iMesh)->listPoints()[elemIdPoint[iPoint]]; |
| 196 |
|
|
} |
| 197 |
|
|
|
| 198 |
|
✗ |
if ( best > pointDistance ) {//case of the point |
| 199 |
|
✗ |
projectOnBestPoint ( elemPoint,thePoint,pointProjection); |
| 200 |
|
|
} else {//case of the edge |
| 201 |
|
✗ |
projectOnBestEdge ( elemPoint,thePoint,pointProjection); |
| 202 |
|
|
} |
| 203 |
|
|
} |
| 204 |
|
|
} |
| 205 |
|
|
// Final assert! |
| 206 |
|
✗ |
double projectionDistance = thePoint->dist(pointProjection); |
| 207 |
|
✗ |
double epsAssert(0.1);//TODO |
| 208 |
|
✗ |
if ( projectionDistance > epsAssert ) { |
| 209 |
|
✗ |
std::cout<<"=========================The point is quite far from the projection."<<std::endl; |
| 210 |
|
✗ |
idForLog=4; |
| 211 |
|
✗ |
pointNotFound++; |
| 212 |
|
|
} |
| 213 |
|
✗ |
printLog(log, idForLog, thePoint, pointProjection, projectionDistance); |
| 214 |
|
|
// 2) IF NEEDED, INITIALIZATION OF THE CURVILINEAR FINITE ELEMENT |
| 215 |
|
✗ |
if ( eltType != previousEltType ) { |
| 216 |
|
✗ |
if ( currentFe ) |
| 217 |
|
✗ |
delete currentFe; |
| 218 |
|
|
|
| 219 |
|
|
// Get standard information for the element |
| 220 |
|
✗ |
int typeOfFiniteElement = var.finiteElementType(); |
| 221 |
|
✗ |
const GeoElement* geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second; |
| 222 |
|
✗ |
const RefElement* refEle = geoEle->defineFiniteEle(eltType,typeOfFiniteElement, *m_pb->mesh(iMesh)); |
| 223 |
|
|
|
| 224 |
|
|
// Initialization |
| 225 |
|
✗ |
currentFe = new CurvilinearFiniteElement(*refEle,*geoEle,var.degreeOfExactness()); |
| 226 |
|
|
|
| 227 |
|
✗ |
nDofEle= currentFe->numDof(); |
| 228 |
|
✗ |
values.resize(nDofEle); |
| 229 |
|
|
|
| 230 |
|
|
// Update of previous type of element |
| 231 |
|
✗ |
previousEltType=eltType; |
| 232 |
|
|
} |
| 233 |
|
|
|
| 234 |
|
|
// 3) UPDATE OF THE BASIS FUNCTIONS IN CURVILINEAR FINITE ELEMENT |
| 235 |
|
✗ |
currentFe->updateMeas(0,elemPoint); |
| 236 |
|
|
|
| 237 |
|
|
// 4) EVALUATION OF THE BASIS FUNCTION IN THE POINT |
| 238 |
|
✗ |
currentFe->evaluateBasisFunctionsInPoint(&pointProjection,values); |
| 239 |
|
|
|
| 240 |
|
|
// 5) FILLING THE DATA STRUCTURE |
| 241 |
|
✗ |
for ( std::size_t iComp(0); iComp < static_cast<std::size_t>(nComp); ++iComp ) { |
| 242 |
|
✗ |
for ( std::size_t iDof(0); iDof < static_cast<std::size_t>(nDofEle); ++iDof ) { |
| 243 |
|
|
// get the global application numbering volume |
| 244 |
|
✗ |
m_pb->dof().loc2glob(ielSupportDof,iDof,ivar,iComp, idColumn); |
| 245 |
|
|
// Tranforming it into globalPetscOrdering volume |
| 246 |
|
✗ |
AOApplicationToPetsc(m_pb->ao(), 1, & idColumn); |
| 247 |
|
|
// Moving to global application ordering bd |
| 248 |
|
✗ |
idColumn = m_dofBD[iBD]->petscVol2ApplicationBD( idColumn ); |
| 249 |
|
|
// Saving the values in the temporary interpolator matrix |
| 250 |
|
✗ |
temporaryInterpolator[ rowNumbering[ cPoint*nComp + iComp ] ].push_back( std::make_pair( idColumn, values[iDof] ) ); |
| 251 |
|
|
} |
| 252 |
|
|
} |
| 253 |
|
|
} |
| 254 |
|
✗ |
this->checkInterpolatorStatistics(foundInsideAnElement, foundInsideAnElementAltLab, foundInEdge, foundInPoint, pointNotFound, nPoints ); |
| 255 |
|
|
|
| 256 |
|
|
// 6) Build the pattern and the data in csr format |
| 257 |
|
✗ |
std::vector<felInt> icsr(m_nrows[iBD]+1); |
| 258 |
|
✗ |
std::vector<felInt> jcsr; |
| 259 |
|
✗ |
std::vector<double> data; |
| 260 |
|
✗ |
icsr[0]=0; |
| 261 |
|
✗ |
for ( felInt iRow(0); iRow< m_nrows[iBD]; ++iRow ) { |
| 262 |
|
✗ |
icsr[iRow+1]=icsr[iRow] + temporaryInterpolator[iRow].size(); |
| 263 |
|
✗ |
if ( temporaryInterpolator[iRow].size () > 5 ) |
| 264 |
|
✗ |
std::cout<<"iRow: "<<iRow<<" size "<<temporaryInterpolator[iRow].size()<<std::endl; |
| 265 |
|
✗ |
for ( std::size_t j(0); j<temporaryInterpolator[iRow].size(); ++j) { |
| 266 |
|
✗ |
jcsr.push_back(temporaryInterpolator[iRow][j].first); |
| 267 |
|
✗ |
data.push_back(temporaryInterpolator[iRow][j].second); |
| 268 |
|
|
} |
| 269 |
|
|
} |
| 270 |
|
|
// 7) Allocating the matrix |
| 271 |
|
✗ |
std::vector<felInt> nnz( m_nrows[iBD] ); |
| 272 |
|
✗ |
for ( std::size_t nRow = 0; nRow < (std::size_t) m_nrows[iBD]; nRow++ ) { |
| 273 |
|
✗ |
nnz[nRow] = icsr[nRow+1]-icsr[nRow]; |
| 274 |
|
|
} |
| 275 |
|
✗ |
m_ncols[iBD] = m_dofBD[iBD]->numGlobalDofInterface(); |
| 276 |
|
✗ |
interpolator.createSeqAIJ(MPI_COMM_SELF, m_nrows[iBD], m_ncols[iBD], 0, nnz.data()); |
| 277 |
|
✗ |
interpolator.setOption(MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); |
| 278 |
|
✗ |
interpolator.setFromOptions(); |
| 279 |
|
✗ |
for ( felInt iRow(0); iRow< m_nrows[iBD]; ++iRow ) { |
| 280 |
|
✗ |
interpolator.setValues(1,&iRow,nnz[iRow],&jcsr[icsr[iRow]],&data[icsr[iRow]],INSERT_VALUES); |
| 281 |
|
|
} |
| 282 |
|
✗ |
interpolator.assembly(MAT_FINAL_ASSEMBLY); |
| 283 |
|
|
}//endif |
| 284 |
|
|
} |
| 285 |
|
|
|
| 286 |
|
✗ |
void SurfaceInterpolator::writeInterfaceFile(felInt iunknown, felInt numComp, std::string filename, std::string folder, std::size_t iBD) const |
| 287 |
|
|
{ |
| 288 |
|
✗ |
if ( MpiInfo::rankProc() == 0 ) { |
| 289 |
|
|
//All the infos are in dofBD |
| 290 |
|
✗ |
m_dofBD[iBD]->writeSurfaceInterpolatorFile(iunknown,numComp,filename,folder); |
| 291 |
|
|
} |
| 292 |
|
|
} |
| 293 |
|
|
|
| 294 |
|
✗ |
void SurfaceInterpolator::readInterfaceFile( std::string filename, std::vector< std::pair< Point*, felInt > >& pointsAndCorrespondingSurfLabel, std::vector<int>& rowNumbering, std::string folder) |
| 295 |
|
|
{ |
| 296 |
|
✗ |
if ( MpiInfo::rankProc() == 0 ) { |
| 297 |
|
|
|
| 298 |
|
|
// 0) Make sure we got empty output variable |
| 299 |
|
✗ |
pointsAndCorrespondingSurfLabel.clear(); |
| 300 |
|
✗ |
rowNumbering.clear(); |
| 301 |
|
|
|
| 302 |
|
|
// 1) Opening the file |
| 303 |
|
✗ |
std::stringstream fname; |
| 304 |
|
✗ |
fname<<folder<<"/"<<filename; |
| 305 |
|
✗ |
std::ifstream iFile; |
| 306 |
|
✗ |
iFile.open(fname.str().c_str()); |
| 307 |
|
✗ |
if ( ! iFile.good() ) { |
| 308 |
|
✗ |
iFile.close(); |
| 309 |
|
✗ |
std::cout<<fname.str()<<std::endl; |
| 310 |
|
✗ |
FEL_ERROR("Problem in opening the file"); |
| 311 |
|
|
} |
| 312 |
|
|
// 2) Get the number of points in the file |
| 313 |
|
|
felInt nPoints; |
| 314 |
|
✗ |
iFile>>nPoints; |
| 315 |
|
✗ |
if ( FelisceParam::verbose() > 2 ) { |
| 316 |
|
✗ |
std::cout<<"Found "<<nPoints<<" points in "<<fname.str().c_str()<<std::endl; |
| 317 |
|
|
} |
| 318 |
|
|
|
| 319 |
|
|
// 3) Temporary variable declaration |
| 320 |
|
|
// Number of components |
| 321 |
|
|
felInt nComp; |
| 322 |
|
|
// The pointer to the Point |
| 323 |
|
|
Point* newPoint; |
| 324 |
|
|
// The surface label |
| 325 |
|
|
felInt label; |
| 326 |
|
|
// The if of the dof |
| 327 |
|
|
felInt id; |
| 328 |
|
|
|
| 329 |
|
|
// 4) Reading |
| 330 |
|
✗ |
for ( int curPoint(0); curPoint<nPoints; ++curPoint) { |
| 331 |
|
|
// get the number of components |
| 332 |
|
✗ |
iFile>>nComp; |
| 333 |
|
|
|
| 334 |
|
|
// get the point |
| 335 |
|
✗ |
newPoint = new Point; |
| 336 |
|
✗ |
iFile>>newPoint->x(); |
| 337 |
|
✗ |
iFile>>newPoint->y(); |
| 338 |
|
✗ |
iFile>>newPoint->z(); |
| 339 |
|
|
|
| 340 |
|
|
// get the label |
| 341 |
|
✗ |
iFile>>label; |
| 342 |
|
|
|
| 343 |
|
|
// verify label correctness, TODO improve this check |
| 344 |
|
|
// especially in case 1 label is mapped into more labels or viceversa. |
| 345 |
|
✗ |
if ( FelisceParam::instance().withCVG ) { |
| 346 |
|
✗ |
int otherLabel=label; |
| 347 |
|
✗ |
FEL_ASSERT( FelisceParam::instance().interfaceLabels.size() == FelisceParam::instance().correspondingInterfaceLabels.size() ); |
| 348 |
|
✗ |
for ( std::size_t k(0); k<FelisceParam::instance().interfaceLabels.size(); ++k) { |
| 349 |
|
✗ |
if ( FelisceParam::instance().correspondingInterfaceLabels[k] == otherLabel ) { |
| 350 |
|
✗ |
label = FelisceParam::instance().interfaceLabels[k]; |
| 351 |
|
|
} |
| 352 |
|
|
} |
| 353 |
|
|
} |
| 354 |
|
|
// save point and label in data structure |
| 355 |
|
✗ |
pointsAndCorrespondingSurfLabel.emplace_back(newPoint,label); |
| 356 |
|
|
|
| 357 |
|
|
// get the id and save them |
| 358 |
|
✗ |
for ( int iComp(0); iComp<nComp; ++iComp) { |
| 359 |
|
✗ |
iFile>>id; |
| 360 |
|
✗ |
rowNumbering.push_back(id); |
| 361 |
|
|
} |
| 362 |
|
|
} |
| 363 |
|
|
} |
| 364 |
|
|
} |
| 365 |
|
|
|
| 366 |
|
✗ |
void SurfaceInterpolator::checkInterpolatorStatistics(int foundInsideAnElement, int foundInsideAnElementAltLab, int foundInEdge, int foundInPoint, int pointNotFound, int nPoints ) const |
| 367 |
|
|
{ |
| 368 |
|
✗ |
if ( FelisceParam::verbose() > 1 ) { |
| 369 |
|
✗ |
std::stringstream msg; |
| 370 |
|
✗ |
double percentageInsideElement = (double) foundInsideAnElement / (double) nPoints; |
| 371 |
|
✗ |
double percentageInsideElementAltLab = (double) foundInsideAnElementAltLab / (double) nPoints; |
| 372 |
|
✗ |
double percentageEdge = (double) foundInEdge / (double) nPoints; |
| 373 |
|
✗ |
double percentagePoint= (double) foundInPoint / (double) nPoints; |
| 374 |
|
|
|
| 375 |
|
✗ |
int prec = 2; |
| 376 |
|
✗ |
int oldPrec=std::cout.precision(); |
| 377 |
|
✗ |
if ( percentageInsideElement*100 < 0.01 ) |
| 378 |
|
✗ |
prec = 4; |
| 379 |
|
✗ |
else if (percentageInsideElement*100 > 1 ) |
| 380 |
|
✗ |
prec=oldPrec; |
| 381 |
|
✗ |
std::cout<<"Points projected on an element: "<<foundInsideAnElement<<"/" |
| 382 |
|
|
<<nPoints<<", percentage=" |
| 383 |
|
✗ |
<<std::setprecision(prec) <<percentageInsideElement*100 |
| 384 |
|
✗ |
<<std::setprecision(oldPrec) <<"%"<<std::endl; |
| 385 |
|
✗ |
if ( percentageInsideElement*100 < 0.01 ) |
| 386 |
|
✗ |
prec = 4; |
| 387 |
|
✗ |
else if (percentageInsideElement*100 > 1 ) |
| 388 |
|
✗ |
prec=oldPrec; |
| 389 |
|
✗ |
std::cout<<"Points projected on an element on a different label than the one specified: "<<foundInsideAnElementAltLab<<"/" |
| 390 |
|
|
<<nPoints<<", percentage=" |
| 391 |
|
✗ |
<<std::setprecision(prec) <<percentageInsideElementAltLab*100 |
| 392 |
|
✗ |
<<std::setprecision(oldPrec) <<"%"<<std::endl; |
| 393 |
|
|
|
| 394 |
|
✗ |
if ( percentageEdge*100 < 0.01 ) |
| 395 |
|
✗ |
prec = 4; |
| 396 |
|
✗ |
std::cout<<"Points projected on an edge: "<<foundInEdge<<"/" |
| 397 |
|
|
<<nPoints<<", percentage=" |
| 398 |
|
✗ |
<<std::setprecision(prec) <<percentageEdge*100 |
| 399 |
|
✗ |
<<std::setprecision(oldPrec) <<"%"<<std::endl; |
| 400 |
|
✗ |
if ( percentagePoint*100 < 0.01 ) |
| 401 |
|
✗ |
prec = 4; |
| 402 |
|
✗ |
std::cout<<"Points projected on the closest point: "<<foundInPoint<<"/" |
| 403 |
|
|
<<nPoints<<", percentage=" |
| 404 |
|
✗ |
<<std::setprecision(prec) <<percentagePoint*100 |
| 405 |
|
✗ |
<<std::setprecision(oldPrec) <<"%"<<std::endl; |
| 406 |
|
|
} |
| 407 |
|
✗ |
if ( pointNotFound > 0 ) { |
| 408 |
|
✗ |
std::stringstream msg; |
| 409 |
|
✗ |
double percentage = (double) pointNotFound/ (double) nPoints; |
| 410 |
|
✗ |
int prec = 2; |
| 411 |
|
✗ |
int oldPrec=std::cout.precision(); |
| 412 |
|
✗ |
if ( percentage*100 < 0.01 ) |
| 413 |
|
✗ |
prec = 4; |
| 414 |
|
✗ |
msg<<"For "<<pointNotFound<<" point(s) the algorithm (total nb of points=" |
| 415 |
|
|
<<nPoints<<", percentage=" |
| 416 |
|
✗ |
<<std::setprecision(prec) <<percentage*100 |
| 417 |
|
✗ |
<<std::setprecision(oldPrec) <<"%), with threshold " |
| 418 |
|
✗ |
<<FelisceParam::instance().interpolatorThreshold<<", was not able to find an element containing the point projection."<<std::endl; |
| 419 |
|
✗ |
FEL_ERROR(msg.str().c_str()); |
| 420 |
|
|
} |
| 421 |
|
|
} |
| 422 |
|
|
|
| 423 |
|
✗ |
void SurfaceInterpolator::initializeSurfaceInterpolator() |
| 424 |
|
|
{ |
| 425 |
|
✗ |
std::vector<Point*> elemPoint; |
| 426 |
|
✗ |
std::vector<felInt> elemIdPoint; |
| 427 |
|
|
// Counter, necessary to compute the id of element by type; |
| 428 |
|
|
felInt numElementPerEltType; |
| 429 |
|
|
|
| 430 |
|
|
// For each type of element present on the boundary: i.e. triangles, quadrangles |
| 431 |
|
✗ |
for (std::size_t cElemType = 0; cElemType < m_mesh->bagElementTypeDomainBoundary().size(); ++cElemType) { |
| 432 |
|
✗ |
ElementType eltType = m_mesh->bagElementTypeDomainBoundary()[cElemType]; |
| 433 |
|
✗ |
elemPoint.resize(GeometricMeshRegion::m_numPointsPerElt[eltType]); |
| 434 |
|
✗ |
elemIdPoint.resize(GeometricMeshRegion::m_numPointsPerElt[eltType]); |
| 435 |
|
✗ |
numElementPerEltType=0; |
| 436 |
|
|
// We loop over the associated labels |
| 437 |
|
✗ |
for (auto itLabelMesh = m_mesh->intRefToBegEndMaps[eltType].begin(); itLabelMesh != m_mesh->intRefToBegEndMaps[eltType].end(); itLabelMesh++) { |
| 438 |
|
|
// We get the label |
| 439 |
|
✗ |
felInt cLabel = itLabelMesh->first; |
| 440 |
|
|
// We get the number of elements present on this label |
| 441 |
|
✗ |
felInt numElemsPerLabel = itLabelMesh->second.second; |
| 442 |
|
|
// For each of this element... |
| 443 |
|
✗ |
for ( felInt iel = 0; iel < numElemsPerLabel; iel++) { |
| 444 |
|
|
// Even in the case of P2, I just need the verteces of the element. |
| 445 |
|
|
// The reason is that I need them to find the element in the future! |
| 446 |
|
|
// Than of course, the interpolation will depend on the type of finite-element |
| 447 |
|
|
// used. |
| 448 |
|
✗ |
m_mesh->getOneElement(eltType, numElementPerEltType+iel, elemIdPoint, 0); |
| 449 |
|
✗ |
for (int iPoint = 0; iPoint < GeometricMeshRegion::m_numPointsPerElt[eltType]; iPoint++) { |
| 450 |
|
✗ |
elemPoint[iPoint]= &m_mesh->listPoints()[elemIdPoint[iPoint]]; |
| 451 |
|
|
// For each point we save the index of the elements of the part of the star |
| 452 |
|
|
// on this label. |
| 453 |
|
✗ |
m_surfLabelAndIdPoint2Star[cLabel][elemIdPoint[iPoint]].push_back(numElementPerEltType+iel); // this is an index local to this proc, pay attention when going parallel |
| 454 |
|
✗ |
m_surfLab2SetOfPointsAndIds[cLabel].insert(std::make_pair(elemPoint[iPoint],elemIdPoint[iPoint])); |
| 455 |
|
|
} |
| 456 |
|
|
} |
| 457 |
|
✗ |
numElementPerEltType+=numElemsPerLabel; |
| 458 |
|
|
} |
| 459 |
|
|
} |
| 460 |
|
✗ |
m_surfaceInterpolatorInitialized = true; |
| 461 |
|
|
} |
| 462 |
|
|
|
| 463 |
|
✗ |
int SurfaceInterpolator::getClosestSurfaceElement( |
| 464 |
|
|
Point* thePoint, |
| 465 |
|
|
felInt surfaceLabel, |
| 466 |
|
|
ElementType& eltTypeReturn, |
| 467 |
|
|
std::vector<Point*>& elemPointReturn, |
| 468 |
|
|
Point& projection, |
| 469 |
|
|
felInt& ielSupportDof, |
| 470 |
|
|
std::map<int,double>& iel2ScoreByEdge, |
| 471 |
|
|
double& distance, |
| 472 |
|
|
int& idElOfClosestPoint |
| 473 |
|
|
) |
| 474 |
|
|
{ |
| 475 |
|
✗ |
if ( ! m_surfaceInterpolatorInitialized ) { |
| 476 |
|
✗ |
FEL_ERROR("You have to call initializeSurfaceInterpolator"); |
| 477 |
|
|
} |
| 478 |
|
✗ |
if ( FelisceParam::verbose() > 2 ) { |
| 479 |
|
✗ |
std::cout<<"Point to be projected: "; thePoint->print(10); |
| 480 |
|
|
} |
| 481 |
|
|
|
| 482 |
|
|
// We extract the std::set of points and the corresponding ids for this label |
| 483 |
|
✗ |
auto setOfPointsAndIds = m_surfLab2SetOfPointsAndIds.at(surfaceLabel); |
| 484 |
|
|
|
| 485 |
|
|
// We start looking for the closest point |
| 486 |
|
|
Point* closest; |
| 487 |
|
|
double curDist; |
| 488 |
|
✗ |
felInt curIdClosestPoint = closestPoint( setOfPointsAndIds, thePoint, closest, curDist ); |
| 489 |
|
✗ |
if ( curDist < distance ) { |
| 490 |
|
✗ |
distance = curDist; |
| 491 |
|
✗ |
std::vector<felInt> star = m_surfLabelAndIdPoint2Star.at(surfaceLabel).at(curIdClosestPoint); |
| 492 |
|
✗ |
idElOfClosestPoint = star[0]; |
| 493 |
|
|
} |
| 494 |
|
✗ |
if ( FelisceParam::verbose() > 2 ) { |
| 495 |
|
✗ |
std::cout<<"The closest point on label "<<surfaceLabel<<" is: "; closest->print(10); |
| 496 |
|
|
} |
| 497 |
|
✗ |
int depth(0); |
| 498 |
|
|
double score; |
| 499 |
|
✗ |
return lookIntoStar(curIdClosestPoint, surfaceLabel, thePoint, eltTypeReturn, elemPointReturn, projection, ielSupportDof, depth, score,iel2ScoreByEdge); |
| 500 |
|
|
} |
| 501 |
|
|
|
| 502 |
|
✗ |
int SurfaceInterpolator::lookIntoStar( |
| 503 |
|
|
int idClosestPoint, |
| 504 |
|
|
int surfaceLabel, |
| 505 |
|
|
Point* thePoint, |
| 506 |
|
|
ElementType & eltTypeReturn, |
| 507 |
|
|
std::vector<Point*>& elemPointReturn, |
| 508 |
|
|
Point& projection, felInt& ielSupportDof, |
| 509 |
|
|
felInt& depth, |
| 510 |
|
|
double& minScore, |
| 511 |
|
|
std::map<int,double>& iel2distFromEdge |
| 512 |
|
|
) |
| 513 |
|
|
{ |
| 514 |
|
|
// We start looking for the closest element |
| 515 |
|
✗ |
felInt idClosestElement=-1; |
| 516 |
|
|
// Points of the element |
| 517 |
|
✗ |
std::vector<Point*> elemPoint; |
| 518 |
|
|
// ..and corresponding ids. |
| 519 |
|
✗ |
std::vector<felInt> elemIdPoint; |
| 520 |
|
|
|
| 521 |
|
|
// We assume that the closest element is in the star of this point |
| 522 |
|
|
// We extract the star |
| 523 |
|
✗ |
std::vector<felInt> star = m_surfLabelAndIdPoint2Star.at(surfaceLabel).at(idClosestPoint); |
| 524 |
|
|
// We loop over its element, until we find the closest one |
| 525 |
|
✗ |
const std::size_t star_size = star.size(); |
| 526 |
|
✗ |
if ( FelisceParam::verbose() > 2 ) { |
| 527 |
|
✗ |
std::cout<<"Size of the star: "<<star_size<<std::endl; |
| 528 |
|
|
} |
| 529 |
|
✗ |
std::vector<double> alpha(star_size); |
| 530 |
|
✗ |
std::vector<double> beta(star_size); |
| 531 |
|
✗ |
for( std::size_t k(0); k<star_size && idClosestElement<0; ++k) { |
| 532 |
|
|
// ==================EXTRACTING ELEMENT INFORMATION==================== // |
| 533 |
|
|
// We get the type of element |
| 534 |
|
✗ |
ElementType eltType=getETandAssert2D(m_mesh->intRefToEnum().at(surfaceLabel)); |
| 535 |
|
|
// ..and we resize this two vectors accoringly |
| 536 |
|
✗ |
elemPoint.resize(GeometricMeshRegion::m_numPointsPerElt[eltType]); |
| 537 |
|
✗ |
elemIdPoint.resize(GeometricMeshRegion::m_numPointsPerElt[eltType]); |
| 538 |
|
|
// We extract the element and its points |
| 539 |
|
✗ |
m_mesh->getOneElement(eltType, star[k], elemIdPoint, 0); |
| 540 |
|
|
// We assume that we have only triangles or quadrangles |
| 541 |
|
✗ |
FEL_ASSERT(elemIdPoint.size()>=3 && elemIdPoint.size()<=4); |
| 542 |
|
|
|
| 543 |
|
|
// We get the pointers of the points trough their ids |
| 544 |
|
✗ |
for (int iPoint = 0; iPoint < GeometricMeshRegion::m_numPointsPerElt[eltType]; iPoint++) { |
| 545 |
|
✗ |
elemPoint[iPoint] = &m_mesh->listPoints()[elemIdPoint[iPoint]]; |
| 546 |
|
|
} |
| 547 |
|
|
|
| 548 |
|
|
// ==================PLANE PROJETION=================================== // |
| 549 |
|
✗ |
projectOnElementPlane(elemPoint,thePoint,projection); |
| 550 |
|
✗ |
if ( FelisceParam::verbose() > 3 ) { |
| 551 |
|
✗ |
std::cout<<" its idByType: "<<star[k]<<std::endl; |
| 552 |
|
✗ |
std::cout<<" its position in the star: "<<k<<std::endl; |
| 553 |
|
✗ |
std::cout<<" its points: "<<std::endl; |
| 554 |
|
✗ |
for ( int iPoint = 0; iPoint < GeometricMeshRegion::m_numPointsPerElt[eltType]; iPoint++) { |
| 555 |
|
✗ |
elemPoint[iPoint]->print(10); |
| 556 |
|
|
} |
| 557 |
|
✗ |
std::cout<<"TheProjection "; |
| 558 |
|
✗ |
projection.print(10); |
| 559 |
|
|
} |
| 560 |
|
|
|
| 561 |
|
|
// ==================CHECK IF THE THE PROJECTION BELONGS TO THE ELEMENT AND COMPUTE BEST EDGE DISTANCE==================== // |
| 562 |
|
✗ |
bool pointInElement=false; |
| 563 |
|
✗ |
if ( elemPoint.size() == 3 ) { |
| 564 |
|
✗ |
double edgeDist=1e20; |
| 565 |
|
✗ |
pointInElement = checkTriangle(projection,elemPoint,alpha[k],beta[k],thePoint,edgeDist); |
| 566 |
|
|
// saving bestEdge information |
| 567 |
|
✗ |
iel2distFromEdge[star[k]]=edgeDist; |
| 568 |
|
✗ |
} else if ( elemPoint.size() == 4 ) { |
| 569 |
|
✗ |
double edgeDist=1.e20; |
| 570 |
|
✗ |
pointInElement = checkQuadrangle(projection,elemPoint,alpha[k],beta[k],thePoint,edgeDist); |
| 571 |
|
|
// saving bestEdge information |
| 572 |
|
✗ |
iel2distFromEdge[star[k]]=edgeDist; |
| 573 |
|
|
} |
| 574 |
|
|
// We check if the projection is inside or outside the element. |
| 575 |
|
✗ |
if ( pointInElement ) { |
| 576 |
|
|
// The correct element has been found |
| 577 |
|
✗ |
idClosestElement = star[k]; |
| 578 |
|
✗ |
ielSupportDof=0; |
| 579 |
|
✗ |
m_mesh->getIdElemFromTypeElemAndIdByType(eltType,idClosestElement,ielSupportDof); |
| 580 |
|
✗ |
eltTypeReturn = eltType; |
| 581 |
|
✗ |
elemPointReturn = elemPoint; |
| 582 |
|
✗ |
if ( FelisceParam::verbose() > 2 ) { |
| 583 |
|
✗ |
std::cout<<" its idByType: "<<idClosestElement<<std::endl; |
| 584 |
|
✗ |
std::cout<<" its position in the star: "<<k<<std::endl; |
| 585 |
|
✗ |
std::cout<<" its points: "<<std::endl; |
| 586 |
|
✗ |
for ( int iPoint = 0; iPoint < GeometricMeshRegion::m_numPointsPerElt[eltType]; iPoint++) { |
| 587 |
|
✗ |
elemPoint[iPoint]->print(10); |
| 588 |
|
|
} |
| 589 |
|
✗ |
std::cout<<"TheProjection "; projection.print(10); |
| 590 |
|
|
} |
| 591 |
|
|
} |
| 592 |
|
|
} |
| 593 |
|
|
// If found it so simply we just skip all this part and return 0 (success); |
| 594 |
|
|
// If we have not found it |
| 595 |
|
✗ |
if (idClosestElement<0) { |
| 596 |
|
✗ |
if ( FelisceParam::verbose() > 3 ) { |
| 597 |
|
✗ |
std::cout<<"Point: "; |
| 598 |
|
✗ |
thePoint->print(10); |
| 599 |
|
|
} |
| 600 |
|
|
// We compute the score of each element in the star! |
| 601 |
|
✗ |
minScore=1e20; |
| 602 |
|
✗ |
for ( std::size_t k(0); k<star.size(); ++k ) { |
| 603 |
|
✗ |
double score=computeScore(alpha[k],beta[k]); |
| 604 |
|
✗ |
if ( score < minScore) { |
| 605 |
|
✗ |
minScore = score; |
| 606 |
|
✗ |
idClosestElement = star[k]; |
| 607 |
|
|
} |
| 608 |
|
|
} |
| 609 |
|
✗ |
if ( FelisceParam::verbose() > 3 ) |
| 610 |
|
✗ |
std::cout<<"SPECIAL CASE OF GET CLOSEST SURFACE ELEMENT: Score of current star: "<<minScore<<std::endl; |
| 611 |
|
|
|
| 612 |
|
|
// For now this element is our best! |
| 613 |
|
|
// We extract the return parameters |
| 614 |
|
✗ |
ElementType eltType=getETandAssert2D(m_mesh->intRefToEnum().at(surfaceLabel)); |
| 615 |
|
✗ |
eltTypeReturn = eltType; |
| 616 |
|
✗ |
ielSupportDof=0; |
| 617 |
|
✗ |
m_mesh->getIdElemFromTypeElemAndIdByType(eltType,idClosestElement,ielSupportDof); |
| 618 |
|
✗ |
m_mesh->getOneElement(eltType, idClosestElement, elemIdPoint, 0); |
| 619 |
|
|
|
| 620 |
|
✗ |
elemPointReturn.resize( GeometricMeshRegion::m_numPointsPerElt[eltType] ); |
| 621 |
|
✗ |
for (int iPoint = 0; iPoint < GeometricMeshRegion::m_numPointsPerElt[eltType]; iPoint++) { |
| 622 |
|
✗ |
elemPointReturn[iPoint] = &m_mesh->listPoints()[elemIdPoint[iPoint]]; |
| 623 |
|
|
} |
| 624 |
|
|
|
| 625 |
|
|
// We also look into the adjacent stars. |
| 626 |
|
✗ |
if ( depth >= FelisceParam::instance().maxDepth ) { |
| 627 |
|
✗ |
if ( FelisceParam::verbose() > 3 ) |
| 628 |
|
✗ |
std::cout<<"Max depth reached, iterations concluded"<<std::endl; |
| 629 |
|
|
} else { |
| 630 |
|
✗ |
depth++; |
| 631 |
|
|
int result; |
| 632 |
|
✗ |
double score=100; |
| 633 |
|
|
ElementType eltTypeReturnTmp; |
| 634 |
|
✗ |
std::vector<Point*> elemPointReturnTmp; |
| 635 |
|
✗ |
Point projectionTmp(0.0); |
| 636 |
|
|
felInt ielSupportDofTmp; |
| 637 |
|
|
// I look in the other stars close to this element. |
| 638 |
|
✗ |
for ( std::size_t k(0); k<elemIdPoint.size(); ++k) { |
| 639 |
|
|
// I skip the current star. |
| 640 |
|
✗ |
if ( elemIdPoint[k] != idClosestPoint ) { |
| 641 |
|
✗ |
if ( FelisceParam::verbose() > 3 ) |
| 642 |
|
✗ |
std::cout<<"Descending into the star of point number "<<k<<std::endl; |
| 643 |
|
✗ |
result = this->lookIntoStar(/*where to look*/elemIdPoint[k], surfaceLabel, thePoint, eltTypeReturnTmp, elemPointReturnTmp, projectionTmp, ielSupportDofTmp, depth , score, iel2distFromEdge ); |
| 644 |
|
|
// We found it in this star! |
| 645 |
|
✗ |
if ( result == 0 ) { |
| 646 |
|
✗ |
if ( FelisceParam::verbose() > 3 ) |
| 647 |
|
✗ |
std::cout<<"found it in the current star!"<<std::endl; |
| 648 |
|
|
// We do not go any further! we update the parameters and return success. |
| 649 |
|
✗ |
eltTypeReturn = eltTypeReturnTmp; |
| 650 |
|
✗ |
elemPointReturn = elemPointReturnTmp; |
| 651 |
|
✗ |
projection = projectionTmp; |
| 652 |
|
✗ |
ielSupportDof = ielSupportDofTmp; |
| 653 |
|
✗ |
return 0;//success |
| 654 |
|
|
} else { |
| 655 |
|
|
// We have not found it, but maybe we found something slightly better |
| 656 |
|
|
|
| 657 |
|
|
// is this better than what we have? |
| 658 |
|
|
// if not we ignore it |
| 659 |
|
✗ |
if ( score < minScore ) { |
| 660 |
|
✗ |
if ( FelisceParam::verbose() > 3 ) |
| 661 |
|
✗ |
std::cout<<"Not found, but score improved!"<<std::endl; |
| 662 |
|
✗ |
minScore = score; |
| 663 |
|
|
//update current best |
| 664 |
|
✗ |
eltTypeReturn = eltTypeReturnTmp; |
| 665 |
|
✗ |
elemPointReturn = elemPointReturnTmp; |
| 666 |
|
✗ |
projection = projectionTmp; |
| 667 |
|
✗ |
ielSupportDof = ielSupportDofTmp; |
| 668 |
|
|
} |
| 669 |
|
|
} |
| 670 |
|
|
} |
| 671 |
|
|
} |
| 672 |
|
|
} |
| 673 |
|
✗ |
return 1; // not found |
| 674 |
|
|
} else { |
| 675 |
|
|
//standard result |
| 676 |
|
✗ |
return 0; // found |
| 677 |
|
|
} |
| 678 |
|
|
} |
| 679 |
|
|
|
| 680 |
|
✗ |
void SurfaceInterpolator::displayDataForSurfaceInterpolator() |
| 681 |
|
|
{ |
| 682 |
|
|
int rank; |
| 683 |
|
|
int numproc; |
| 684 |
|
✗ |
MPI_Comm_rank(PETSC_COMM_WORLD, &rank); |
| 685 |
|
✗ |
MPI_Comm_size(PETSC_COMM_WORLD, &numproc); |
| 686 |
|
|
|
| 687 |
|
✗ |
for ( int cProc(0); cProc<numproc; ++cProc ) { |
| 688 |
|
✗ |
if( rank == cProc ) { |
| 689 |
|
✗ |
std::cout << "---------- Proc numero: " << rank << std::endl; |
| 690 |
|
✗ |
for (auto i = m_surfLab2SetOfPointsAndIds.begin(); i != m_surfLab2SetOfPointsAndIds.end(); ++i) { |
| 691 |
|
✗ |
std::cout << "----- " << "Label: " << i->first << std::endl; |
| 692 |
|
✗ |
for (auto j = i->second.begin(); j != i->second.end(); ++j) { |
| 693 |
|
✗ |
std::cout << "P: "; j->first->print(10); |
| 694 |
|
✗ |
std::cout << "--> Number of Elements in the star: " << m_surfLabelAndIdPoint2Star[i->first][j->second].size() << std::endl; |
| 695 |
|
✗ |
std::cout << "----> Elements ids: "; |
| 696 |
|
✗ |
for(std::size_t k = 0; k < m_surfLabelAndIdPoint2Star[i->first][j->second].size(); ++k) |
| 697 |
|
✗ |
std::cout << m_surfLabelAndIdPoint2Star[i->first][j->second][k] << " "; |
| 698 |
|
✗ |
std::cout << std::endl; |
| 699 |
|
|
} |
| 700 |
|
|
} |
| 701 |
|
|
} |
| 702 |
|
✗ |
MPI_Barrier(MpiInfo::petscComm()); |
| 703 |
|
|
} |
| 704 |
|
|
} |
| 705 |
|
|
|
| 706 |
|
✗ |
GeometricMeshRegion::ElementType SurfaceInterpolator::getETandAssert2D(std::set<ElementType> setOfET){ |
| 707 |
|
✗ |
int nTypeElPerLab(0); |
| 708 |
|
✗ |
ElementType eltType = static_cast<ElementType>(0); |
| 709 |
|
✗ |
for ( std::set<ElementType>::iterator et=setOfET.begin(); |
| 710 |
|
✗ |
et!=setOfET.end(); |
| 711 |
|
✗ |
++et ) { |
| 712 |
|
✗ |
int first3DElement = 12; //TODO! |
| 713 |
|
✗ |
int first2DElement = 4; //TODO! |
| 714 |
|
✗ |
ElementType eltType2 = *et; |
| 715 |
|
✗ |
if( (int) eltType2 < first3DElement && (int)eltType2 >= first2DElement ) { |
| 716 |
|
✗ |
nTypeElPerLab++; |
| 717 |
|
✗ |
eltType=*et; |
| 718 |
|
|
} |
| 719 |
|
|
} |
| 720 |
|
|
// We assume that each label can have only one type of element. |
| 721 |
|
✗ |
FEL_ASSERT(nTypeElPerLab==1); |
| 722 |
|
✗ |
FEL_ASSERT(eltType!=0); |
| 723 |
|
✗ |
return eltType; |
| 724 |
|
|
} |
| 725 |
|
|
|
| 726 |
|
|
|
| 727 |
|
|
|
| 728 |
|
|
|
| 729 |
|
|
|
| 730 |
|
|
|
| 731 |
|
✗ |
bool SurfaceInterpolator::isPointInTriangle(const double* A,const double* B,const double* C,const double* D, double& alpha, double& beta) |
| 732 |
|
|
{ |
| 733 |
|
|
|
| 734 |
|
|
double V1[3],V2[3],H[3]; |
| 735 |
|
|
|
| 736 |
|
✗ |
V1[0]=B[0]-A[0]; |
| 737 |
|
✗ |
V1[1]=B[1]-A[1]; |
| 738 |
|
✗ |
V1[2]=B[2]-A[2]; |
| 739 |
|
|
|
| 740 |
|
✗ |
V2[0]=C[0]-A[0]; |
| 741 |
|
✗ |
V2[1]=C[1]-A[1]; |
| 742 |
|
✗ |
V2[2]=C[2]-A[2]; |
| 743 |
|
|
|
| 744 |
|
✗ |
H[0]=D[0]-A[0]; |
| 745 |
|
✗ |
H[1]=D[1]-A[1]; |
| 746 |
|
✗ |
H[2]=D[2]-A[2]; |
| 747 |
|
|
|
| 748 |
|
|
double det,det1,det2; |
| 749 |
|
✗ |
det = V1[0]*V2[1] - V2[0]*V1[1]; |
| 750 |
|
✗ |
det1 = V1[0]*V2[2] - V2[0]*V1[2]; |
| 751 |
|
✗ |
det2 = V1[1]*V2[2] - V2[1]*V1[2]; |
| 752 |
|
|
double fdet,fdet1,fdet2; |
| 753 |
|
✗ |
fdet=std::fabs(det); |
| 754 |
|
✗ |
fdet1=std::fabs(det1); |
| 755 |
|
✗ |
fdet2=std::fabs(det2); |
| 756 |
|
|
|
| 757 |
|
✗ |
alpha = -1; |
| 758 |
|
✗ |
beta = -1; |
| 759 |
|
✗ |
if ( fdet> fdet1 && fdet > fdet2 ) { |
| 760 |
|
✗ |
alpha=(V2[1]*H[0] - V2[0]*H[1])/det; |
| 761 |
|
✗ |
beta =(V1[0]*H[1] - V1[1]*H[0])/det; |
| 762 |
|
✗ |
} else if ( fdet1> fdet && fdet1 > fdet2 ){ |
| 763 |
|
✗ |
alpha=(V2[2]*H[0] - V2[0]*H[2])/det1; |
| 764 |
|
✗ |
beta =(V1[0]*H[2] - V1[2]*H[0])/det1; |
| 765 |
|
✗ |
} else if ( fdet2> fdet && fdet2 > fdet1 ){ |
| 766 |
|
✗ |
alpha=(V2[2]*H[1] - V2[1]*H[2])/det2; |
| 767 |
|
✗ |
beta =(V1[1]*H[2] - V1[2]*H[1])/det2; |
| 768 |
|
|
} |
| 769 |
|
✗ |
double eps = FelisceParam::instance().interpolatorThreshold; |
| 770 |
|
✗ |
return alpha>=-eps && beta >= -eps && alpha+beta<=1+eps; |
| 771 |
|
|
} |
| 772 |
|
|
|
| 773 |
|
✗ |
void SurfaceInterpolator::checkEdge(const double *p0,const double *p1,Point* thePoint,double& bestDistance) |
| 774 |
|
|
{ |
| 775 |
|
|
double distance,t; |
| 776 |
|
✗ |
bool inside = isPointInEdgeAndIfYesHowFar( p0, p1, thePoint->coor(), distance, t); |
| 777 |
|
✗ |
if ( inside && distance < bestDistance ) { |
| 778 |
|
✗ |
bestDistance = distance; |
| 779 |
|
|
} |
| 780 |
|
|
} |
| 781 |
|
|
|
| 782 |
|
✗ |
bool SurfaceInterpolator::checkTriangle( const Point& projection, const std::vector<Point*>& elemPoint, double& alpha, double& beta, Point* thePoint, double& bestEdgeDistance ) |
| 783 |
|
|
{ |
| 784 |
|
|
// Check if it is in the triangle |
| 785 |
|
✗ |
bool pointInElement(false); |
| 786 |
|
✗ |
pointInElement = isPointInTriangle(elemPoint[0]->coor(),elemPoint[1]->coor(),elemPoint[2]->coor(),projection.coor(),alpha,beta); |
| 787 |
|
✗ |
if ( pointInElement ){ |
| 788 |
|
✗ |
return pointInElement; //no need to check the edges; |
| 789 |
|
|
} |
| 790 |
|
✗ |
checkEdge(elemPoint[0]->coor(),elemPoint[1]->coor(),thePoint,bestEdgeDistance); |
| 791 |
|
✗ |
checkEdge(elemPoint[1]->coor(),elemPoint[2]->coor(),thePoint,bestEdgeDistance); |
| 792 |
|
✗ |
checkEdge(elemPoint[2]->coor(),elemPoint[0]->coor(),thePoint,bestEdgeDistance); |
| 793 |
|
✗ |
return pointInElement; |
| 794 |
|
|
} |
| 795 |
|
|
|
| 796 |
|
✗ |
bool SurfaceInterpolator::checkQuadrangle( const Point& projection, const std::vector<Point*>& elemPoint, double& alpha, double& beta, Point* thePoint, double& bestEdgeDistance ) |
| 797 |
|
|
{ |
| 798 |
|
✗ |
bool pointInElement(false); |
| 799 |
|
|
double alphaTmp,betaTmp, bestScore; |
| 800 |
|
✗ |
bestScore=1e20; |
| 801 |
|
✗ |
pointInElement = isPointInTriangle(elemPoint[0]->coor(),elemPoint[1]->coor(),elemPoint[2]->coor(),projection.coor(),alphaTmp,betaTmp); |
| 802 |
|
✗ |
updateBestScore(alphaTmp,betaTmp,bestScore,alpha,beta); |
| 803 |
|
✗ |
pointInElement = isPointInTriangle(elemPoint[1]->coor(),elemPoint[2]->coor(),elemPoint[3]->coor(),projection.coor(),alphaTmp,betaTmp) || pointInElement; |
| 804 |
|
✗ |
updateBestScore(alphaTmp,betaTmp,bestScore,alpha,beta); |
| 805 |
|
✗ |
pointInElement = isPointInTriangle(elemPoint[2]->coor(),elemPoint[3]->coor(),elemPoint[0]->coor(),projection.coor(),alphaTmp,betaTmp) || pointInElement; |
| 806 |
|
✗ |
updateBestScore(alphaTmp,betaTmp,bestScore,alpha,beta); |
| 807 |
|
✗ |
pointInElement = isPointInTriangle(elemPoint[3]->coor(),elemPoint[0]->coor(),elemPoint[1]->coor(),projection.coor(),alphaTmp,betaTmp) || pointInElement; |
| 808 |
|
✗ |
updateBestScore(alphaTmp,betaTmp,bestScore,alpha,beta); |
| 809 |
|
✗ |
if ( pointInElement ){ |
| 810 |
|
✗ |
return pointInElement; //no need to check the edges; |
| 811 |
|
|
} |
| 812 |
|
✗ |
checkEdge(elemPoint[0]->coor(),elemPoint[1]->coor(),thePoint,bestEdgeDistance); |
| 813 |
|
✗ |
checkEdge(elemPoint[1]->coor(),elemPoint[2]->coor(),thePoint,bestEdgeDistance); |
| 814 |
|
✗ |
checkEdge(elemPoint[2]->coor(),elemPoint[3]->coor(),thePoint,bestEdgeDistance); |
| 815 |
|
✗ |
checkEdge(elemPoint[3]->coor(),elemPoint[0]->coor(),thePoint,bestEdgeDistance); |
| 816 |
|
✗ |
return pointInElement; |
| 817 |
|
|
} |
| 818 |
|
|
|
| 819 |
|
✗ |
bool SurfaceInterpolator::isPointInEdgeAndIfYesHowFar( const double *p0,const double *p1,const double* thePoint,double& distance, double&t) |
| 820 |
|
|
{ |
| 821 |
|
✗ |
projectOnEdge(p0,p1,thePoint,t); |
| 822 |
|
✗ |
double eps = FelisceParam::instance().interpolatorThreshold; |
| 823 |
|
✗ |
if ( t < 1 + eps && t > -eps ) { |
| 824 |
|
✗ |
double dist(0); |
| 825 |
|
✗ |
for ( int i(0); i<3; ++i ) { |
| 826 |
|
✗ |
dist += ((1-t)*p0[i] +t*p1[i]-thePoint[i])*((1-t)*p0[i] +t*p1[i]-thePoint[i]); |
| 827 |
|
|
} |
| 828 |
|
✗ |
distance = std::sqrt(dist); |
| 829 |
|
✗ |
return true; |
| 830 |
|
|
} else { |
| 831 |
|
✗ |
distance = 1.e10; |
| 832 |
|
✗ |
return false; |
| 833 |
|
|
} |
| 834 |
|
|
} |
| 835 |
|
|
|
| 836 |
|
|
// elemPoint [in] points defining the plane, only the first three are used |
| 837 |
|
|
// thePoint [in] the point to be projected |
| 838 |
|
|
// projection [out] the projection of this point on the plane |
| 839 |
|
✗ |
void SurfaceInterpolator::projectOnElementPlane( const std::vector<Point*>& elemPoint, Point* thePoint, Point& projection ) |
| 840 |
|
|
{ |
| 841 |
|
|
Point normal; // Normal to the element plane |
| 842 |
|
✗ |
MathUtilities::CrossProduct(normal.getCoor(), (*elemPoint[0] - *elemPoint[1]).getCoor(), (*elemPoint[0] - *elemPoint[2]).getCoor()); |
| 843 |
|
|
|
| 844 |
|
✗ |
double nnorm = normal.norm(); |
| 845 |
|
✗ |
for( int coor(0); coor<3; coor++) { |
| 846 |
|
✗ |
normal.coor(coor)/=nnorm; |
| 847 |
|
|
} |
| 848 |
|
|
// We compute the intercept of the plane containing the element |
| 849 |
|
✗ |
double intercept = (*elemPoint[0]) * normal; |
| 850 |
|
|
|
| 851 |
|
|
// P_proj = P + s*normal |
| 852 |
|
|
// We compute s, by imposing that P_proj is on the plane. |
| 853 |
|
✗ |
double s = intercept - (*thePoint) * normal; |
| 854 |
|
✗ |
projection = s * normal + (*thePoint); |
| 855 |
|
|
} |
| 856 |
|
|
|
| 857 |
|
✗ |
void SurfaceInterpolator::projectOnEdge(const double *p0,const double *p1,const double* thePoint, double& t) |
| 858 |
|
|
{ |
| 859 |
|
✗ |
double denom(0); |
| 860 |
|
✗ |
double num(0); |
| 861 |
|
✗ |
for ( int i(0); i<3; ++i ) { |
| 862 |
|
✗ |
num += - p0[i]*p1[i]-thePoint[i]*p0[i]+thePoint[i]*p1[i] + p0[i]*p0[i]; |
| 863 |
|
✗ |
denom += -2*p0[i]*p1[i] +p1[i]*p1[i] + p0[i]*p0[i]; |
| 864 |
|
|
} |
| 865 |
|
✗ |
t = num/denom; |
| 866 |
|
|
} |
| 867 |
|
|
|
| 868 |
|
✗ |
void SurfaceInterpolator::projectOnBestEdge(const std::vector<Point*>& elemPoint, Point* thePoint, Point& projection) |
| 869 |
|
|
{ |
| 870 |
|
✗ |
double bestDist(1e20),bestT(0.0),bestEdgeA(0.0),bestEdgeB(0.0); |
| 871 |
|
|
double d,t; |
| 872 |
|
✗ |
for (std::size_t iPoint(0); iPoint < elemPoint.size(); iPoint++) { |
| 873 |
|
|
// We assume that the point describing the element are in an order such that |
| 874 |
|
|
// the edges are described by consecutive points |
| 875 |
|
✗ |
isPointInEdgeAndIfYesHowFar(elemPoint[iPoint]->coor(),elemPoint[( iPoint + 1 ) % elemPoint.size() ]->coor(),thePoint->coor(), d,t); |
| 876 |
|
✗ |
if ( d < bestDist ) { |
| 877 |
|
✗ |
bestDist = d; |
| 878 |
|
✗ |
bestT = t; |
| 879 |
|
✗ |
bestEdgeA = iPoint; |
| 880 |
|
✗ |
bestEdgeB = (iPoint+1)%elemPoint.size(); |
| 881 |
|
|
} |
| 882 |
|
|
} |
| 883 |
|
✗ |
for( int coor(0); coor<3; coor++) { |
| 884 |
|
✗ |
projection.coor()[coor] = (1-bestT)*elemPoint[bestEdgeA]->coor()[coor]+bestT*elemPoint[bestEdgeB]->coor()[coor]; |
| 885 |
|
|
} |
| 886 |
|
|
} |
| 887 |
|
|
|
| 888 |
|
✗ |
void SurfaceInterpolator::projectOnBestPoint(const std::vector<Point*>& elemPoint, Point* thePoint, Point& projection) |
| 889 |
|
|
{ |
| 890 |
|
✗ |
double best=1e20; |
| 891 |
|
✗ |
for (std::size_t iPoint = 0; iPoint < elemPoint.size(); iPoint++) { |
| 892 |
|
✗ |
double dist = elemPoint[iPoint]->dist( *thePoint); |
| 893 |
|
✗ |
if ( dist<best ) { |
| 894 |
|
✗ |
best=dist; |
| 895 |
|
✗ |
projection=*elemPoint[iPoint]; |
| 896 |
|
|
} |
| 897 |
|
|
} |
| 898 |
|
|
} |
| 899 |
|
|
|
| 900 |
|
✗ |
double SurfaceInterpolator::computeScore(double alpha, double beta) |
| 901 |
|
|
{ |
| 902 |
|
✗ |
double score(0); |
| 903 |
|
✗ |
if ( alpha < 0 ) score += -alpha; |
| 904 |
|
✗ |
if ( alpha > 1 ) score += alpha-1; |
| 905 |
|
✗ |
if ( beta < 0 ) score += -beta; |
| 906 |
|
✗ |
if ( beta > 1 ) score += beta-1; |
| 907 |
|
✗ |
if ( alpha+beta < 0 ) score += -(alpha+beta); |
| 908 |
|
✗ |
if ( alpha+beta > 1 ) score += alpha+beta-1; |
| 909 |
|
✗ |
return score; |
| 910 |
|
|
} |
| 911 |
|
|
|
| 912 |
|
✗ |
void SurfaceInterpolator::updateBestScore( double alpha, double beta, double& bestScore, double& alphabest, double& betabest) |
| 913 |
|
|
{ |
| 914 |
|
✗ |
double score = computeScore(alpha,beta); |
| 915 |
|
✗ |
if( score < bestScore ) { |
| 916 |
|
✗ |
bestScore = score; |
| 917 |
|
✗ |
alphabest = alpha; |
| 918 |
|
✗ |
betabest = beta; |
| 919 |
|
|
} |
| 920 |
|
|
} |
| 921 |
|
|
|
| 922 |
|
✗ |
int SurfaceInterpolator::closestPoint(const std::set< std::pair< Point*, felInt > >& setOfPointsAndIds, const Point* thePoint, Point*& closest, double& distance) |
| 923 |
|
|
{ |
| 924 |
|
|
// We start looking for the closest point |
| 925 |
|
✗ |
felInt idClosestPoint=-1; |
| 926 |
|
✗ |
double dist(0), distMin=1.e30; |
| 927 |
|
|
// We compute the distance between each point and our point |
| 928 |
|
✗ |
for(auto cPoint=setOfPointsAndIds.begin(); |
| 929 |
|
✗ |
cPoint!=setOfPointsAndIds.end(); ++cPoint ) { |
| 930 |
|
|
// Distance |
| 931 |
|
✗ |
dist = 0; |
| 932 |
|
✗ |
for( int coor(0); coor<3; coor++) { |
| 933 |
|
✗ |
dist += std::pow((thePoint->coor(coor)-cPoint->first->coor(coor)),2); |
| 934 |
|
|
} |
| 935 |
|
|
// Check if it is closer than the previous best |
| 936 |
|
✗ |
if ( dist < distMin ) { |
| 937 |
|
✗ |
distMin=dist; |
| 938 |
|
✗ |
idClosestPoint=cPoint->second; |
| 939 |
|
✗ |
closest = cPoint->first; |
| 940 |
|
|
} |
| 941 |
|
|
} |
| 942 |
|
✗ |
distance = std::sqrt(distMin); |
| 943 |
|
✗ |
return idClosestPoint; |
| 944 |
|
|
} |
| 945 |
|
|
|
| 946 |
|
|
} |
| 947 |
|
|
|
| 948 |
|
|
|