Skip to content
Snippets Groups Projects
Commit ea6d4ad3 authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

*** empty log message ***

parent 07c95f33
No related branches found
No related tags found
No related merge requests found
......@@ -1321,12 +1321,17 @@ void GFaceCompound::buildOct() const
//Verify topological conditions for computing the harmonic map
bool GFaceCompound::checkTopology() const
{
// FIXME!!! I think those things are wrong with cross-patch reparametrization
if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)return true;
bool correctTopo = true;
int Nb = _interior_loops.size();
int G = genusGeom() ;
if( G != 0 || Nb < 1) correctTopo = false;
// printf("%d %d\n",Nb,G);
if (!correctTopo){
Msg::Info("--> Wrong topology: Genus=%d and N boundary=%d", G, Nb);
nbSplit = std::max(G+1, 2);
......
......@@ -718,18 +718,25 @@ void signedDistancesPointsTriangle (std::vector<double>&distances,
SVector3 t3 = p3-p2;
SVector3 n = crossprod(t1,t2);
n.normalize();
double mat[3][3] = {{t1.x(),t2.x(),-n.x()},
{t1.y(),t2.y(),-n.y()},
{t1.z(),t2.z(),-n.z()}};
double inv[3][3];
inv3x3(mat,inv);
double det = inv3x3(mat,inv);
distances.clear();
distances.resize(pts.size());
for (int i=0; i<pts.size();i++)
distances[i] = 1.e22;
if (det = 0.0)return;
const double n2t1 = dot(t1,t1);
const double n2t2 = dot(t2,t3);
const double n2t3 = dot(t3,t3);
distances.clear();
distances.resize(pts.size());
double u,v,d;
for (int i=0; i<pts.size();i++){
......@@ -739,7 +746,7 @@ void signedDistancesPointsTriangle (std::vector<double>&distances,
v = (inv[1][0] * pp1.x() + inv[1][1] * pp1.y() + inv[1][2] * pp1.z());
d = (inv[2][0] * pp1.x() + inv[2][1] * pp1.y() + inv[2][2] * pp1.z());
double sign = (d>0) ? 1.0 : -1.0;
if (d == 0) sign = 1.e22;
if (d == 0) sign = 1.e10;
if (u >= 0 && v >=0 && 1.-u-v >= 0.0){
distances[i] = d;
......@@ -749,7 +756,7 @@ void signedDistancesPointsTriangle (std::vector<double>&distances,
const double t13 = dot(pp1,t2) / n2t2;
SVector3 pp2 = p - p2;
const double t23 = dot(pp2,t3) / n2t3;
d = 1.e22;
d = 1.e10;
bool found = false;
if (t12 >=0 && t12 <= 1.0){
d = sign * std::min(fabs(d),p.distance(p1+(p2-p1)*t12));
......
......@@ -63,7 +63,7 @@ private:
}
typename std::map<INDEXTYPE,scalar>::const_iterator begin() const {return _nodal_values.begin();}
typename std::map<INDEXTYPE,scalar>::const_iterator ens() const {return _nodal_values.endn();}
typename std::map<INDEXTYPE,scalar>::const_iterator end() const {return _nodal_values.end();}
// add that in the ann search tool
void insert_point (double x, double y, double z){
......@@ -72,19 +72,27 @@ private:
// compute distance
double distance (double x, double y, double z) const;
// set the value
void setValue(INDEXTYPE i, scalar s){
_nodal_values[i] = s;
}
inline INDEXTYPE index_of_element (double x, double y, double z) const{
// P = P_0 + xi * _vdx + eta * _vdy + zeta *vdz
// DP = P-P_0 * _vdx = xi
SVector3 DP (x-_X,y-_Y,z-_Z);
const double xi = dot(DP,_xiAxis);
const double eta = dot(DP,_etaAxis);
const double zeta = dot(DP,_zetaAxis);
double xi = dot(DP,_xiAxis);
double eta = dot(DP,_etaAxis);
double zeta = dot(DP,_zetaAxis);
int i = xi/_dxi * _Nxi;
int j = eta/_deta * _Neta;
int k = zeta/_dzeta * _Nzeta;
const int i = xi/_dxi * _Nxi;
const int j = eta/_deta * _Neta;
const int k = zeta/_dzeta * _Nzeta;
if (i < 0) i = 0;if (i >= _Nxi) i = _Nxi-1;
if (j < 0) j = 0;if (j >= _Neta) j = _Neta-1;
if (k < 0) k = 0;if (k >= _Nzeta) k = _Nzeta-1;
return element_index(i,j,k);
}
......@@ -179,7 +187,8 @@ private:
typename std::map<INDEXTYPE,scalar>::const_iterator it = _nodal_values.begin();
for ( ; it!=_nodal_values.end();++it){
SPoint3 p = coordinates_of_node(it->first);
fprintf(f,"%d %g\n",it->first,distance(p.x(),p.y(),p.z()));
// fprintf(f,"%d %g\n",it->first,distance(p.x(),p.y(),p.z()));
fprintf(f,"%d %g\n",it->first,it->second);
}
fprintf(f,"$EndNodeData\n");
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment