Skip to content
Snippets Groups Projects
Commit 1a908da3 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

fix

parent 822a90f4
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Geo/discreteDiskFace.cpp
+ 79
− 124
View file @ 1a908da3
...@@ -180,7 +180,7 @@ static void discreteDiskFaceCentroid(void *a, double*c) ...@@ -180,7 +180,7 @@ static void discreteDiskFaceCentroid(void *a, double*c)
c[2] = 0.0; c[2] = 0.0;
} }
static int discreteDiskFaceInEle(void *a, double*c)// # mark static int discreteDiskFaceInEle(void *a, double*c)
{ {
discreteDiskFaceTriangle *t = (discreteDiskFaceTriangle *)a; discreteDiskFaceTriangle *t = (discreteDiskFaceTriangle *)a;
double Xi[2]; double Xi[2];
...@@ -196,55 +196,48 @@ static int discreteDiskFaceInEle(void *a, double*c)// # mark ...@@ -196,55 +196,48 @@ static int discreteDiskFaceInEle(void *a, double*c)// # mark
static bool orderVertices(const double &tot_length, const std::vector<MVertex*> &l, static bool orderVertices(const double &tot_length, const std::vector<MVertex*> &l,
std::vector<double> &coord) std::vector<double> &coord)
{ // called once by constructor ; organize the vertices for the linear system {
// called once by constructor ; organize the vertices for the linear system
// expressing the mapping // expressing the mapping
coord.clear(); coord.clear();
coord.push_back(0.); coord.push_back(0.);
MVertex* first = l[0]; MVertex* first = l[0];
for(unsigned int i=1; i < l.size(); i++){ for(unsigned int i=1; i < l.size(); i++){
MVertex* next = l[i]; MVertex* next = l[i];
const double length = sqrt( (next->x() - first->x()) * (next->x() - first->x()) + const double length = sqrt( (next->x() - first->x()) * (next->x() - first->x()) +
(next->y() - first->y()) * (next->y() - first->y()) + (next->y() - first->y()) * (next->y() - first->y()) +
(next->z() - first->z()) * (next->z() - first->z()) ); (next->z() - first->z()) * (next->z() - first->z()) );
coord.push_back(coord[coord.size()-1] + length / tot_length); coord.push_back(coord[coord.size()-1] + length / tot_length);
first = next; first = next;
} }
return true; return true;
} }
/*BUILDER*/
discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation, discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
int p, std::vector<GFace*> *CAD) : int p, std::vector<GFace*> *CAD)
GFace(gf->model(),diskTriangulation->idNum), _parent (gf), _ddft(NULL), oct(NULL) : GFace(gf->model(), diskTriangulation->idNum), _parent (gf), _ddft(NULL), oct(NULL)
{ {
initialTriangulation = diskTriangulation; initialTriangulation = diskTriangulation;
std::vector<MElement*> mesh = diskTriangulation->tri; std::vector<MElement*> mesh = diskTriangulation->tri;
_order = p; _order = p;
_n = (p+1)*(p+2)/2; _n = (p+1)*(p+2)/2;
discrete_triangles.resize(mesh.size()); discrete_triangles.resize(mesh.size());
std::map<MVertex*,MVertex*> v2v;// mesh vertex |-> face vertex std::map<MVertex*,MVertex*> v2v; // mesh vertex |-> face vertex
std::map<MEdge,MVertex*,Less_Edge> ed2nodes; // edge to interior node(s) std::map<MEdge,MVertex*,Less_Edge> ed2nodes; // edge to interior node(s)
for (unsigned int i=0;i<mesh.size();i++){ // triangle by triangle for (unsigned int i=0;i<mesh.size();i++){ // triangle by triangle
std::vector<MVertex*> vs; // MTriangle vertices std::vector<MVertex*> vs; // MTriangle vertices
for (unsigned int j=0; j<3; j++){ // loop over vertices AND edges of the current triangle // loop over vertices AND edges of the current triangle
MVertex *v = mesh[i]->getVertex(j);// firstly, edge vertices for (unsigned int j = 0; j < 3; j++){
MVertex *v = mesh[i]->getVertex(j); // firstly, edge vertices
if (v->onWhat()->dim() == 2) { if (v->onWhat()->dim() == 2) {
std::map<MVertex*,MVertex*> :: iterator it = v2v.find(v); std::map<MVertex*,MVertex*>::iterator it = v2v.find(v);
if (it == v2v.end()){ if (it == v2v.end()){
MFaceVertex *vv; MFaceVertex *vv;
if (!CAD) vv = new MFaceVertex ( v->x(), v->y(), v->z(), v->onWhat(), 0, 0); if(!CAD || (*CAD)[i] != v->onWhat()){
vv = new MFaceVertex(v->x(), v->y(), v->z(), v->onWhat(), 0, 0);
}
else{ else{
GFace *cad = (*CAD)[i]; double pu, pv; v->getParameter(0, pu); v->getParameter(1, pv);
if(cad != v->onWhat())
Msg::Fatal("Line %d FILE %s : erroneous cad list",__LINE__,__FILE__);
double pu,pv; v->getParameter(0,pu);v->getParameter(1,pv);
vv = new MFaceVertex ( v->x(), v->y(), v->z(), v->onWhat(), pu, pv); vv = new MFaceVertex ( v->x(), v->y(), v->z(), v->onWhat(), pu, pv);
} }
v2v[v] = vv; v2v[v] = vv;
...@@ -293,7 +286,6 @@ discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation, ...@@ -293,7 +286,6 @@ discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
//putOnView(gf->tag(),diskTriangulation->idNum,true,true); //putOnView(gf->tag(),diskTriangulation->idNum,true,true);
//printParamMesh(); //printParamMesh();
} }
/*end BUILDER*/
discreteDiskFace::~discreteDiskFace() discreteDiskFace::~discreteDiskFace()
{ {
...@@ -338,7 +330,6 @@ bool discreteDiskFace::parametrize() const ...@@ -338,7 +330,6 @@ bool discreteDiskFace::parametrize() const
linearSystem<double> * lsys_u, *lsys_v; linearSystem<double> * lsys_u, *lsys_v;
#ifdef HAVE_MUMPS #ifdef HAVE_MUMPS
lsys_u = new linearSystemMUMPS<double>; lsys_u = new linearSystemMUMPS<double>;
lsys_v = new linearSystemMUMPS<double>; lsys_v = new linearSystemMUMPS<double>;
...@@ -389,7 +380,8 @@ bool discreteDiskFace::parametrize() const ...@@ -389,7 +380,8 @@ bool discreteDiskFace::parametrize() const
lsys_u->systemSolve(); lsys_u->systemSolve();
lsys_v->systemSolve(); lsys_v->systemSolve();
Msg::Debug("Systems solved"); Msg::Debug("Systems solved");
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){ for(std::set<MVertex *>::iterator itv = allNodes.begin();
itv !=allNodes.end() ; ++itv){
MVertex *v = *itv; MVertex *v = *itv;
double value_U, value_V; double value_U, value_V;
myAssemblerU.getDofValue(v, 0, 1, value_U); myAssemblerU.getDofValue(v, 0, 1, value_U);
...@@ -415,15 +407,15 @@ bool discreteDiskFace::parametrize() const ...@@ -415,15 +407,15 @@ bool discreteDiskFace::parametrize() const
return true; return true;
} }
void discreteDiskFace::getTriangleUV(const double u,const double v, void discreteDiskFace::getTriangleUV(const double u,const double v,
discreteDiskFaceTriangle **mt, discreteDiskFaceTriangle **mt,
double &_xi, double &_eta)const{ double &_xi, double &_eta)const
{
double uv[3] = {u,v,0.}; double uv[3] = {u,v,0.};
*mt = (discreteDiskFaceTriangle*) Octree_Search(uv,oct); *mt = (discreteDiskFaceTriangle*) Octree_Search(uv,oct);
if (!(*mt)) { if (!(*mt)) {
for (unsigned int i=0; i<discrete_triangles.size()-geoTriangulation->fillingHoles.size(); i++){ for (unsigned int i = 0; i < discrete_triangles.size() -
geoTriangulation->fillingHoles.size(); i++){
discreteDiskFaceTriangle *ct = &_ddft[i]; discreteDiskFaceTriangle *ct = &_ddft[i];
double Xi[2]; double Xi[2];
int xxx = discreteDiskFaceInEle(ct, Xi); int xxx = discreteDiskFaceInEle(ct, Xi);
...@@ -434,7 +426,9 @@ void discreteDiskFace::getTriangleUV(const double u,const double v, ...@@ -434,7 +426,9 @@ void discreteDiskFace::getTriangleUV(const double u,const double v,
return; return;
} }
} }
Msg::Debug("discreteDiskFace::getTriangleUV(), didn't find the reference coordinate (xi;eta) for (u;v)=(%f;%f) among %d triangles",u,v,discrete_triangles.size()-geoTriangulation->fillingHoles.size()); Msg::Debug("discreteDiskFace::getTriangleUV(), didn't find the reference "
"coordinate (xi;eta) for (u;v)=(%f;%f) among %d triangles",
u,v,discrete_triangles.size()-geoTriangulation->fillingHoles.size());
return; return;
} }
...@@ -442,7 +436,8 @@ void discreteDiskFace::getTriangleUV(const double u,const double v, ...@@ -442,7 +436,8 @@ void discreteDiskFace::getTriangleUV(const double u,const double v,
double U[2] = {u,v}; double U[2] = {u,v};
bool pass = uv2xi(*mt,U,Xi); bool pass = uv2xi(*mt,U,Xi);
if (!pass){ if (!pass){
Msg::Error("discreteDiskFace::getTriangleUV(), didn't find the reference coordinate (xi;eta) for (u;v)=(%f;%f)",u,v); Msg::Error("discreteDiskFace::getTriangleUV(), didn't find the reference "
"coordinate (xi;eta) for (u;v)=(%f;%f)", u, v);
return; return;
} }
_xi = Xi[0]; _xi = Xi[0];
...@@ -453,7 +448,7 @@ bool discreteDiskFace::checkOrientationUV() ...@@ -453,7 +448,7 @@ bool discreteDiskFace::checkOrientationUV()
{ {
discreteDiskFaceTriangle *ct; discreteDiskFaceTriangle *ct;
if(_order==1){ if(_order == 1){
double current; // initial and current orientation double current; // initial and current orientation
ct = &_ddft[0]; ct = &_ddft[0];
double p1[2] = {ct->p[0].x(), ct->p[0].y()}; double p1[2] = {ct->p[0].x(), ct->p[0].y()};
...@@ -471,7 +466,8 @@ bool discreteDiskFace::checkOrientationUV() ...@@ -471,7 +466,8 @@ bool discreteDiskFace::checkOrientationUV()
else nbP++; else nbP++;
} }
if (nbP*nbM){ if (nbP*nbM){
Msg::Info("Map %d of the atlas : Triangles have different orientations (%d + / %d -)",tag(),nbP,nbM); Msg::Info("Map %d of the atlas: triangles have different orientations (%d + / %d -)",
tag(), nbP, nbM);
return false; return false;
} }
return true; return true;
...@@ -481,7 +477,8 @@ bool discreteDiskFace::checkOrientationUV() ...@@ -481,7 +477,8 @@ bool discreteDiskFace::checkOrientationUV()
double min, max; double min, max;
std::vector<MVertex*> localVertices; std::vector<MVertex*> localVertices;
localVertices.resize(_n); localVertices.resize(_n);
for(unsigned int i=0; i<discrete_triangles.size()-geoTriangulation->fillingHoles.size(); i++){ for(unsigned int i=0; i<discrete_triangles.size() -
geoTriangulation->fillingHoles.size(); i++){
ct = &_ddft[i]; ct = &_ddft[i];
for(int j=0; j<_n; j++) for(int j=0; j<_n; j++)
localVertices[j] = new MVertex(ct->p[j].x(),ct->p[j].y(),0.); localVertices[j] = new MVertex(ct->p[j].x(),ct->p[j].y(),0.);
...@@ -524,11 +521,14 @@ void discreteDiskFace::optimize() ...@@ -524,11 +521,14 @@ void discreteDiskFace::optimize()
// -- generation of parametric nodes // -- generation of parametric nodes
std::map<SPoint3,MVertex*> sp2mv; std::map<SPoint3,MVertex*> sp2mv;
std::vector<MElement*> paramTriangles; std::vector<MElement*> paramTriangles;
for(std::map<MVertex*,SPoint3>::iterator it=coordinates.begin(); it!= coordinates.end(); ++it) for(std::map<MVertex*,SPoint3>::iterator it=coordinates.begin();
it != coordinates.end(); ++it)
sp2mv[it->second] = new MVertex(it->second.x(),it->second.y(),0.); sp2mv[it->second] = new MVertex(it->second.x(),it->second.y(),0.);
// -- generation of parametric triangles // -- generation of parametric triangles
paramTriangles.resize(discrete_triangles.size() - geoTriangulation->fillingHoles.size()); paramTriangles.resize(discrete_triangles.size() -
for(unsigned int i=0; i<discrete_triangles.size() -geoTriangulation->fillingHoles.size(); i++){ geoTriangulation->fillingHoles.size());
for(unsigned int i = 0; i < discrete_triangles.size() -
geoTriangulation->fillingHoles.size(); i++){
discreteDiskFaceTriangle* ct = &_ddft[i]; discreteDiskFaceTriangle* ct = &_ddft[i];
std::vector<MVertex*> mv; std::vector<MVertex*> mv;
mv.resize(ct->tri->getNumVertices()); mv.resize(ct->tri->getNumVertices());
...@@ -551,36 +551,35 @@ void discreteDiskFace::optimize() ...@@ -551,36 +551,35 @@ void discreteDiskFace::optimize()
// ---- discrete vertex // ---- discrete vertex
std::set<GFace*, GEntityLessThan>::iterator it = paramDisk->firstFace(); std::set<GFace*, GEntityLessThan>::iterator it = paramDisk->firstFace();
GFace *dgf = *it; GFace *dgf = *it;
discreteVertex *dv = new discreteVertex(paramDisk,paramDisk->getMaxElementaryNumber(0)+1); discreteVertex *dv = new discreteVertex(paramDisk,
paramDisk->getMaxElementaryNumber(0) + 1);
sp2mv[coordinates[_U0[0]]]->setEntity(dv); sp2mv[coordinates[_U0[0]]]->setEntity(dv);
dv->mesh_vertices.push_back(sp2mv[coordinates[_U0[0]]]); dv->mesh_vertices.push_back(sp2mv[coordinates[_U0[0]]]);
todelete.insert(sp2mv[coordinates[_U0[0]]]); todelete.insert(sp2mv[coordinates[_U0[0]]]);
paramDisk->add(dv); paramDisk->add(dv);
// ---- discrete edge // ---- discrete edge
discreteEdge *de = new discreteEdge(paramDisk,paramDisk->getMaxElementaryNumber(1)+1,dv,dv); discreteEdge *de = new discreteEdge(paramDisk,
paramDisk->getMaxElementaryNumber(1)+1, dv, dv);
paramDisk->add(de); paramDisk->add(de);
std::vector<MLine*> lines; std::vector<MLine*> lines;
for(unsigned int i=1; i<_U0.size(); i++){ for(unsigned int i=1; i<_U0.size(); i++){
sp2mv[coordinates[_U0[i]]]->setEntity(de); sp2mv[coordinates[_U0[i]]]->setEntity(de);
de->mesh_vertices.push_back(sp2mv[coordinates[_U0[i]]]); de->mesh_vertices.push_back(sp2mv[coordinates[_U0[i]]]);
todelete.insert(sp2mv[coordinates[_U0[i]]]); todelete.insert(sp2mv[coordinates[_U0[i]]]);
lines.push_back(new MLine(sp2mv[coordinates[_U0[i-1]]],sp2mv[coordinates[_U0[i]]])); lines.push_back(new MLine(sp2mv[coordinates[_U0[i-1]]],
sp2mv[coordinates[_U0[i]]]));
} }
lines.push_back(new MLine(sp2mv[coordinates[_U0[_U0.size()-1]]],sp2mv[coordinates[_U0[0]]])); lines.push_back(new MLine(sp2mv[coordinates[_U0[_U0.size()-1]]],
sp2mv[coordinates[_U0[0]]]));
de->setTopo(lines); de->setTopo(lines);
de->createGeometry();// !!!! setTopo ... MLine's de->createGeometry();// !!!! setTopo ... MLine's
// optimization // optimization
if(_order >1) if(_order >1)
HighOrderMeshOptimizer(paramDisk, optParams); HighOrderMeshOptimizer(paramDisk, optParams);
else else
MeshQualityOptimizer(paramDisk,opt); MeshQualityOptimizer(paramDisk,opt);
// update the parametrization // update the parametrization
paramTriangles = e2e[0]; paramTriangles = e2e[0];
for(unsigned int i=0; i< paramTriangles.size(); i++){ for(unsigned int i=0; i< paramTriangles.size(); i++){
...@@ -606,7 +605,6 @@ void discreteDiskFace::optimize() ...@@ -606,7 +605,6 @@ void discreteDiskFace::optimize()
// cleaning // cleaning
delete paramDisk; delete paramDisk;
#endif #endif
} }
...@@ -647,27 +645,24 @@ SPoint2 discreteDiskFace::parFromVertex(MVertex *v) const ...@@ -647,27 +645,24 @@ SPoint2 discreteDiskFace::parFromVertex(MVertex *v) const
std::map<MVertex*,SPoint3>::iterator it = coordinates.find(v); std::map<MVertex*,SPoint3>::iterator it = coordinates.find(v);
if(it != coordinates.end()) return SPoint2(it->second.x(),it->second.y()); if(it != coordinates.end()) return SPoint2(it->second.x(),it->second.y());
// The 1D mesh has been re-done // The 1D mesh has been re-done
if (v->onWhat()->dim()==1){ if (v->onWhat()->dim() == 1){
if (v->onWhat()->geomType() == DiscreteCurve){ if (v->onWhat()->geomType() == DiscreteCurve){
discreteEdge *de = dynamic_cast<discreteEdge*> (v->onWhat()); discreteEdge *de = dynamic_cast<discreteEdge*> (v->onWhat());
if (de){ if(de){
MVertex *v1,*v2; MVertex *v1,*v2;
double xi; double xi;
de->interpolateInGeometry (v,&v1,&v2,xi); // modify de->interpolateInGeometry (v,&v1,&v2,xi); // modify
std::map<MVertex*,SPoint3>::iterator it1 = coordinates.find(v1); std::map<MVertex*,SPoint3>::iterator it1 = coordinates.find(v1);
std::map<MVertex*,SPoint3>::iterator it2 = coordinates.find(v2); std::map<MVertex*,SPoint3>::iterator it2 = coordinates.find(v2);
if(it1 == coordinates.end()) Msg::Fatal("FIXME TO DO %d %s",__LINE__,__FILE__); if(it1 != coordinates.end() && it2 != coordinates.end()){
if(it2 == coordinates.end()) Msg::Fatal("FIXME TO DO %d %s",__LINE__,__FILE__); return SPoint2(it1->second.x(),it1->second.y()) * (1.-xi) +
return SPoint2(it1->second.x(),it1->second.y()) * (1.-xi) + SPoint2(it2->second.x(),it2->second.y()) * xi; // modify
SPoint2(it2->second.x(),it2->second.y()) * xi; // modify }
} }
} }
Msg::Fatal("FIXME TO DO %d %s",__LINE__,__FILE__);
} }
else if (v->onWhat()->dim()==0)
Msg::Fatal("discreteDiskFace::parFromVertex vertex classified on a model "
"vertex that is not part of the face");
Msg::Error("discreteDiskFace::parFromVertex failed");
return SPoint2(0,0); return SPoint2(0,0);
} }
...@@ -678,15 +673,13 @@ SVector3 discreteDiskFace::normal(const SPoint2 &param) const ...@@ -678,15 +673,13 @@ SVector3 discreteDiskFace::normal(const SPoint2 &param) const
double discreteDiskFace::curvatureMax(const SPoint2 &param) const double discreteDiskFace::curvatureMax(const SPoint2 &param) const
{ {
throw; return 0;
return false;
} }
double discreteDiskFace::curvatures(const SPoint2 &param, SVector3 *dirMax, SVector3 *dirMin, double discreteDiskFace::curvatures(const SPoint2 &param, SVector3 *dirMax, SVector3 *dirMin,
double *curvMax, double *curvMin) const double *curvMax, double *curvMin) const
{ {
throw; return 0;
return false;
} }
Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 &param) const Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 &param) const
...@@ -709,8 +702,7 @@ Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 &param) const ...@@ -709,8 +702,7 @@ Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 &param) const
double dudxi[2][2] = {{0.,0.},{0.,0.}}; double dudxi[2][2] = {{0.,0.},{0.,0.}};
for (int io=0; io<_n; io++){ for (int io = 0; io < _n; io++){
double X = tri->getVertex(io)->x(); double X = tri->getVertex(io)->x();
double Y = tri->getVertex(io)->y(); double Y = tri->getVertex(io)->y();
double Z = tri->getVertex(io)->z(); double Z = tri->getVertex(io)->z();
...@@ -732,7 +724,6 @@ Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 &param) const ...@@ -732,7 +724,6 @@ Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 &param) const
dudxi[1][0] += V*df[io][0]; dudxi[1][0] += V*df[io][0];
dudxi[1][1] += V*df[io][1]; dudxi[1][1] += V*df[io][1];
} }
double dxidu[2][2]; double dxidu[2][2];
...@@ -755,19 +746,18 @@ Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 &param) const ...@@ -755,19 +746,18 @@ Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 &param) const
void discreteDiskFace::secondDer(const SPoint2 &param, void discreteDiskFace::secondDer(const SPoint2 &param,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
{ // cf Sauvage's thesis {
return; // cf Sauvage's thesis
} }
void discreteDiskFace::putOnView(int iFace, int iMap, bool Xu, bool Ux) void discreteDiskFace::putOnView(int iFace, int iMap, bool Xu, bool Ux)
{// #improveme using built-in methods {
// FIXME using built-in methods
char mybuffer [64]; char mybuffer [64];
FILE *view_u=NULL, *view_v=NULL, *UVx=NULL, *UVy=NULL, *UVz=NULL; FILE *view_u=NULL, *view_v=NULL, *UVx=NULL, *UVy=NULL, *UVz=NULL;
if(Xu){ if(Xu){
sprintf(mybuffer, "param_u_gface%d_part%d_order%d.pos", sprintf(mybuffer, "param_u_gface%d_part%d_order%d.pos",
iFace, iMap,_order); iFace, iMap,_order);
...@@ -803,7 +793,8 @@ void discreteDiskFace::putOnView(int iFace, int iMap, bool Xu, bool Ux) ...@@ -803,7 +793,8 @@ void discreteDiskFace::putOnView(int iFace, int iMap, bool Xu, bool Ux)
fprintf(UVy,"View \"y(U)\"{\n"); fprintf(UVy,"View \"y(U)\"{\n");
fprintf(UVz,"View \"z(U)\"{\n"); fprintf(UVz,"View \"z(U)\"{\n");
} }
for (unsigned int i=0; i<discrete_triangles.size()-geoTriangulation->fillingHoles.size(); i++){ for (unsigned int i = 0; i < discrete_triangles.size() -
geoTriangulation->fillingHoles.size(); i++){
discreteDiskFaceTriangle* my_ddft = &_ddft[i]; discreteDiskFaceTriangle* my_ddft = &_ddft[i];
if (_order == 1){ if (_order == 1){
if(Xu){ if(Xu){
...@@ -893,26 +884,13 @@ void discreteDiskFace::putOnView(int iFace, int iMap, bool Xu, bool Ux) ...@@ -893,26 +884,13 @@ void discreteDiskFace::putOnView(int iFace, int iMap, bool Xu, bool Ux)
} }
} }
// useful for mesh generators
// Intersection of a circle and a plane
//FILE *allProblems = NULL;
//void openProblemsON(void){
// allProblems = fopen ("op.pos","w");
//}
//void openProblemsOFF(void){
// fclose(allProblems);
// allProblems = NULL;
//}
GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVector3 &n2, GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVector3 &n2,
const SVector3 &p, const double &R, const SVector3 &p, const double &R,
double uv[2]) const double uv[2]) const
{ {
SVector3 n = crossprod(n1,n2); SVector3 n = crossprod(n1,n2);
n.normalize(); n.normalize();
// printf("n %g %g %g\n",n.x(), n.y(), n.z());
const int N = (int)(discrete_triangles.size()-geoTriangulation->fillingHoles.size()); const int N = (int)(discrete_triangles.size()-geoTriangulation->fillingHoles.size());
for (int i=-1;i<N;i++){ for (int i=-1;i<N;i++){
discreteDiskFaceTriangle *ct = NULL; discreteDiskFaceTriangle *ct = NULL;
...@@ -1001,32 +979,18 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto ...@@ -1001,32 +979,18 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
GPoint pp(0); GPoint pp(0);
pp.setNoSuccess(); pp.setNoSuccess();
Msg::Debug("ARGG no success intersection circle"); Msg::Debug("ARGG no success intersection circle");
// Msg::Info("ARGG no success intersection circle");
// printf("Point(1) = {%g,%g,%g};\n",p.x(),p.y(),p.z());
// printf("Point(2) = {%g,%g,%g};\n",p.x()+d*n1.x(),p.y()+d*n1.y(),p.z()+d*n1.z());
// printf("Point(3) = {%g,%g,%g};\n",p.x()+d*n2.x(),p.y()+d*n2.y(),p.z()+d*n2.z());
// // printf("Circle(4) = {2,1,3};\n");
// printf("{%g,%g,%g};\n",n1.x(),n1.y(),n1.z());
// printf("{%g,%g,%g};\n",n2.x(),n2.y(),n2.z());
// printf("coucou --> \n");
// if (allProblems){
// fprintf(allProblems,"VP(%g,%g,%g){%g,%g,%g};\n",p.x(),p.y(),p.z(),R*n2.x(),R*n2.y(),R*n2.z());
// }
// getchar();
return pp; return pp;
} }
GPoint discreteDiskFace::intersectionWithCircle2(const SVector3 &n1, const SVector3 &n2, GPoint discreteDiskFace::intersectionWithCircle2(const SVector3 &n1, const SVector3 &n2,
const SVector3 &p, const double &d, const SVector3 &p, const double &d,
double uv[2]) const double uv[2]) const
{ {
// n2 is exterior // n2 is exterior
SVector3 n = crossprod(n1,n2); SVector3 n = crossprod(n1,n2);
n.normalize(); n.normalize();
for (int i=-1;i<(int)(discrete_triangles.size()-geoTriangulation->fillingHoles.size());i++){ for (int i = -1; i < (int)(discrete_triangles.size() -
geoTriangulation->fillingHoles.size()); i++){
discreteDiskFaceTriangle *ct = NULL; discreteDiskFaceTriangle *ct = NULL;
double U,V; double U,V;
if (i == -1) getTriangleUV(uv[0],uv[1], &ct, U,V); if (i == -1) getTriangleUV(uv[0],uv[1], &ct, U,V);
...@@ -1118,7 +1082,8 @@ GPoint discreteDiskFace::intersectionWithCircle2(const SVector3 &n1, const SVect ...@@ -1118,7 +1082,8 @@ GPoint discreteDiskFace::intersectionWithCircle2(const SVector3 &n1, const SVect
// printf("{%g,%g,%g};\n",n2.x(),n2.y(),n2.z()); // printf("{%g,%g,%g};\n",n2.x(),n2.y(),n2.z());
// printf("coucou --> \n"); // printf("coucou --> \n");
// if (allProblems){ // if (allProblems){
// fprintf(allProblems,"VP(%g,%g,%g){%g,%g,%g};\n",p.x(),p.y(),p.z(),d*n2.x(),d*n2.y(),d*n2.z()); // fprintf(allProblems,"VP(%g,%g,%g){%g,%g,%g};\n",
// p.x(),p.y(),p.z(),d*n2.x(),d*n2.y(),d*n2.z());
// } // }
// getchar(); // getchar();
return pp; return pp;
...@@ -1126,7 +1091,8 @@ GPoint discreteDiskFace::intersectionWithCircle2(const SVector3 &n1, const SVect ...@@ -1126,7 +1091,8 @@ GPoint discreteDiskFace::intersectionWithCircle2(const SVector3 &n1, const SVect
// computes some kind of maximal distance in a mesh // computes some kind of maximal distance in a mesh
static void addTo (std::map<MVertex*, std::vector<MElement*> > &v2t, MVertex *v, MElement *t) static void addTo(std::map<MVertex*, std::vector<MElement*> > &v2t,
MVertex *v, MElement *t)
{ {
std::map<MVertex*, std::vector<MElement*> > :: iterator it = v2t.find(v); std::map<MVertex*, std::vector<MElement*> > :: iterator it = v2t.find(v);
if (it == v2t.end()){ if (it == v2t.end()){
...@@ -1145,26 +1111,16 @@ static void update(std::map<MVertex*,double> &Close, MVertex *v2, double d) ...@@ -1145,26 +1111,16 @@ static void update(std::map<MVertex*,double> &Close, MVertex *v2, double d)
} }
static MEdge getEdge (MElement *t, MVertex *v) static MEdge getEdge(MElement *t, MVertex *v)
{ {
for (int i=0;i<3;i++) for (int i=0;i<3;i++)
if (t->getVertex(i) == v) return t->getEdge((i+1)%3); if (t->getVertex(i) == v) return t->getEdge((i+1)%3);
return MEdge(); return MEdge();
} }
/* #warning inline double computeDistance(MVertex *v1, double d1, MVertex *v2, double d2,
static double computeDistanceLinePoint (MVertex *v1, MVertex *v2, MVertex *v){ MVertex *v)
{
SVector3 U = v2->point() - v1->point();
SVector3 BA = v2->point() - v->point();
SVector3 xx = crossprod(U,BA);
return xx.norm() / U.norm();
}
*/
inline double computeDistance (MVertex *v1, double d1, MVertex *v2, double d2, MVertex *v){
// o------------a // o------------a
// //
// //
...@@ -1179,9 +1135,6 @@ inline double computeDistance (MVertex *v1, double d1, MVertex *v2, double d2, M ...@@ -1179,9 +1135,6 @@ inline double computeDistance (MVertex *v1, double d1, MVertex *v2, double d2, M
return std::min(d2+v2->distance(v),d1+v1->distance(v)); return std::min(d2+v2->distance(v),d1+v1->distance(v));
double a = v2->distance(v1); double a = v2->distance(v1);
// center (seed) to compute the distance (put it down) // center (seed) to compute the distance (put it down)
...@@ -1276,12 +1229,14 @@ double triangulation::geodesicDistance () ...@@ -1276,12 +1229,14 @@ double triangulation::geodesicDistance ()
std::map<MVertex*,double>::iterator it1 = Fixed.find(ed.getVertex(1)); std::map<MVertex*,double>::iterator it1 = Fixed.find(ed.getVertex(1));
// printf("coucou %p %p\n",ed.getVertex(0),ed.getVertex(1)); // printf("coucou %p %p\n",ed.getVertex(0),ed.getVertex(1));
if (it0 != Fixed.end() && it1 == Fixed.end()){ if (it0 != Fixed.end() && it1 == Fixed.end()){
double d = computeDistance (it->first,it->second,it0->first,it0->second,ed.getVertex(1)); double d = computeDistance (it->first,it->second,it0->first,it0->second,
ed.getVertex(1));
// printf("neigh %d fixed 0 --> d = %g\n",i,d); // printf("neigh %d fixed 0 --> d = %g\n",i,d);
update(Close, ed.getVertex(1), d); update(Close, ed.getVertex(1), d);
} }
else if (it1 != Fixed.end() && it0 == Fixed.end()){ else if (it1 != Fixed.end() && it0 == Fixed.end()){
double d = computeDistance (it->first,it->second,it1->first,it1->second,ed.getVertex(0)); double d = computeDistance(it->first,it->second,it1->first,it1->second,
ed.getVertex(0));
// printf("neigh %d fixed 1 --> d = %g\n",i,d); // printf("neigh %d fixed 1 --> d = %g\n",i,d);
update(Close, ed.getVertex(0), d); update(Close, ed.getVertex(0), d);
} }
...@@ -1302,13 +1257,12 @@ double triangulation::geodesicDistance () ...@@ -1302,13 +1257,12 @@ double triangulation::geodesicDistance ()
fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n", fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
tri[i]->getVertex(0)->x(),tri[i]->getVertex(0)->y(),tri[i]->getVertex(0)->z(), tri[i]->getVertex(0)->x(),tri[i]->getVertex(0)->y(),tri[i]->getVertex(0)->z(),
tri[i]->getVertex(1)->x(),tri[i]->getVertex(1)->y(),tri[i]->getVertex(1)->z(), tri[i]->getVertex(1)->x(),tri[i]->getVertex(1)->y(),tri[i]->getVertex(1)->z(),
tri[i]->getVertex(2)->x(),tri[i]->getVertex(2)->y(),tri[i]->getVertex(2)->z(),d0,d1,d2); tri[i]->getVertex(2)->x(),tri[i]->getVertex(2)->y(),tri[i]->getVertex(2)->z(),
d0,d1,d2);
} }
fprintf(f,"};\n"); fprintf(f,"};\n");
fclose(f); fclose(f);
return CLOSEST; return CLOSEST;
} }
...@@ -1367,7 +1321,8 @@ void discreteDiskFace::printParamMesh() ...@@ -1367,7 +1321,8 @@ void discreteDiskFace::printParamMesh()
sprintf(buffer,"param_mesh%d.msh",tag()); sprintf(buffer,"param_mesh%d.msh",tag());
FILE* pmesh = fopen(buffer,"w"); FILE* pmesh = fopen(buffer,"w");
int count = 1; int count = 1;
fprintf(pmesh,"$MeshFormat\n2.2 0 8\n$EndMeshFormat\n$Nodes\n%u\n",(unsigned int)allNodes.size()); fprintf(pmesh,"$MeshFormat\n2.2 0 8\n$EndMeshFormat\n$Nodes\n%u\n",
(unsigned int)allNodes.size());
for(std::set<MVertex*>::iterator it = allNodes.begin(); it!=allNodes.end(); ++it){ for(std::set<MVertex*>::iterator it = allNodes.begin(); it!=allNodes.end(); ++it){
fprintf(pmesh,"%d %f %f 0\n",count,(coordinates[(*it)]).x(),(coordinates[(*it)]).y()); fprintf(pmesh,"%d %f %f 0\n",count,(coordinates[(*it)]).x(),(coordinates[(*it)]).y());
mv2int[*it] = count; mv2int[*it] = count;
......
benchmarks/2d/compound.geo
+ 1
− 1
View file @ 1a908da3
...@@ -23,4 +23,4 @@ Line(15) = {8, 7}; ...@@ -23,4 +23,4 @@ Line(15) = {8, 7};
Line Loop(100) = {13, -14, 15}; Line Loop(100) = {13, -14, 15};
Plane Surface(11) = {10,100}; Plane Surface(11) = {10,100};
Compound Surface(12)={9,11}; Compound Surface{9,11};
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment