diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp index a4acb2df9200be146c6fabdd83fcdead059d3e4a..726ab5289dcc461e99201740166e4d8095286e1e 100644 --- a/Geo/GFaceCompound.cpp +++ b/Geo/GFaceCompound.cpp @@ -13,45 +13,41 @@ class gmshGradientBasedDiffusivity : public gmshFunction { + private: MElement *_current; int _iComp; - mutable std::map<MVertex*,SPoint2> _coordinates; -public: + mutable std::map<MVertex*, SPoint2> _coordinates; + public: gmshGradientBasedDiffusivity (std::map<MVertex*,SPoint2> &coordinates) : _current (0), _iComp(-1),_coordinates(coordinates){} - void setCurrent (MElement *current){_current=current;} - void setComponent (int iComp){_iComp=iComp;} - virtual double operator () (double x, double y, double z) const { - if (_iComp == -1)return 1.0; - double xyz[3]={x,y,z},uvw[3]; - _current->xyz2uvw(xyz,uvw); + void setCurrent (MElement *current){ _current = current; } + void setComponent (int iComp){ _iComp = iComp; } + virtual double operator () (double x, double y, double z) const + { + if(_iComp == -1) return 1.0; + double xyz[3] = {x, y, z}, uvw[3]; + _current->xyz2uvw(xyz, uvw); double value1[256]; double value2[256]; - for (int i=0;i<_current->getNumVertices();i++){ + for (int i = 0; i < _current->getNumVertices(); i++){ const SPoint2 p = _coordinates[_current->getVertex(i)]; value1[i] = p[0]; value2[i] = p[1]; } - double g1[3],g2[3]; - _current->interpolateGrad(value1,uvw[0],uvw[1],uvw[2],g1); - _current->interpolateGrad(value2,uvw[0],uvw[1],uvw[2],g2); - // printf("%g %g %g -- %g %g %g \n",g[0],g[1],g[2],value[0],value[1],value[2]); - return sqrt(g1[0]*g1[0]+g1[1]*g1[1]+g1[2]*g1[2] + - g2[0]*g2[0]+g2[1]*g2[1]+g2[2]*g2[2]); + double g1[3], g2[3]; + _current->interpolateGrad(value1, uvw[0], uvw[1], uvw[2], g1); + _current->interpolateGrad(value2, uvw[0], uvw[1], uvw[2], g2); + return sqrt(g1[0] * g1[0] + g1[1] * g1[1] + g1[2] * g1[2] + + g2[0] * g2[0] + g2[1] * g2[1] + g2[2] * g2[2]); } }; - -static void fixEdgeToValue (GEdge *ed, - double value, - gmshAssembler &myAssembler) +static void fixEdgeToValue (GEdge *ed, double value, gmshAssembler &myAssembler) { - myAssembler.fixVertex(ed->getBeginVertex()->mesh_vertices[0], - 0, 1, value); - myAssembler.fixVertex(ed->getEndVertex()->mesh_vertices[0], - 0, 1, value); - for (unsigned int i = 0 ; i < ed->mesh_vertices.size(); i++){ - myAssembler.fixVertex(ed->mesh_vertices[i],0, 1, value); + myAssembler.fixVertex(ed->getBeginVertex()->mesh_vertices[0], 0, 1, value); + myAssembler.fixVertex(ed->getEndVertex()->mesh_vertices[0], 0, 1, value); + for (unsigned int i = 0; i < ed->mesh_vertices.size(); i++){ + myAssembler.fixVertex(ed->mesh_vertices[i], 0, 1, value); } } @@ -59,17 +55,16 @@ void GFaceCompound::parametrize() const { if (!oct){ coordinates.clear(); - parametrize(0,0); - parametrize(1,0); - // parametrize(0,1); - // parametrize(1,1); - computeNormals () ; + parametrize(0, 0); + parametrize(1, 0); + // parametrize(0,1); + // parametrize(1,1); + computeNormals(); buildOct(); } } - -void GFaceCompound::getBoundingEdges () +void GFaceCompound::getBoundingEdges() { if (_U0.size() && !_U1.size()){ std::list<GEdge*> :: const_iterator it = _U0.begin(); @@ -82,86 +77,76 @@ void GFaceCompound::getBoundingEdges () std::set<GEdge*> _unique; std::multiset<GEdge*> _touched; - std::list<GFace*> :: iterator it = _compound.begin(); - for ( ; it != _compound.end() ; ++it){ + std::list<GFace*>::iterator it = _compound.begin(); + for ( ; it != _compound.end(); ++it){ std::list<GEdge*> ed = (*it)->edges(); std::list<GEdge*> :: iterator ite = ed.begin(); - for ( ; ite != ed.end() ; ++ite){ + for ( ; ite != ed.end(); ++ite){ _touched.insert(*ite); } } it = _compound.begin(); - for ( ; it != _compound.end() ; ++it){ + for ( ; it != _compound.end(); ++it){ std::list<GEdge*> ed = (*it)->edges(); - std::list<GEdge*> :: iterator ite = ed.begin(); - for ( ; ite != ed.end() ; ++ite){ - if (!(*ite)->degenerate(0) && _touched.count(*ite) == 1)_unique.insert(*ite); - } - } - - + std::list<GEdge*>::iterator ite = ed.begin(); + for ( ; ite != ed.end(); ++ite){ + if (!(*ite)->degenerate(0) && _touched.count(*ite) == 1) _unique.insert(*ite); + } + } std::set<GEdge*>::iterator itf = _unique.begin(); - for ( ; itf != _unique.end() ; ++itf){ + for ( ; itf != _unique.end(); ++itf){ l_edges.push_back(*itf); (*itf)->addFace(this); } - } GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound, - std::list<GEdge*> &U0, - std::list<GEdge*> &V0, - std::list<GEdge*> &U1, - std::list<GEdge*> &V1) : - GFace(m,tag),_compound(compound),_U0(U0),_U1(U1),_V0(V0),_V1(V1), oct(0) + std::list<GEdge*> &U0, std::list<GEdge*> &V0, + std::list<GEdge*> &U1, std::list<GEdge*> &V1) + : GFace(m, tag), _compound(compound), _U0(U0), _U1(U1), _V0(V0), _V1(V1), oct(0) { getBoundingEdges(); } GFaceCompound::~GFaceCompound() { - if (oct){ + if(oct){ Octree_Delete(oct); delete [] _gfct; } } -static bool orderVertices (const std::list<GEdge*> &e , - std::vector<MVertex*> &l, - std::vector<double> &coord){ +static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l, + std::vector<double> &coord) +{ l.clear(); coord.clear(); - std::list<GEdge*> :: const_iterator it = e.begin(); + std::list<GEdge*>::const_iterator it = e.begin(); std::list<MLine*> temp; double tot_length = 0; - for ( ; it != e.end() ; ++it ){ - // printf("GLine %d\n",(*it)->tag()); - for (unsigned int i = 0 ; i < (*it)->lines.size(); i++ ){ - temp.push_back((*it)->lines[i]); + for ( ; it != e.end(); ++it ){ + for (unsigned int i = 0; i < (*it)->lines.size(); i++ ){ + temp.push_back((*it)->lines[i]); MVertex *v0 = (*it)->lines[i]->getVertex(0); MVertex *v1 = (*it)->lines[i]->getVertex(1); - // printf("UNSORTED %d -> %d \n",v0->getNum(),v1->getNum()); - const double length = sqrt ((v0->x() - v1->x()) * (v0->x() - v1->x()) + - (v0->y() - v1->y()) * (v0->y() - v1->y()) + - (v0->z() - v1->z()) * (v0->z() - v1->z()) ); + const double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) + + (v0->y() - v1->y()) * (v0->y() - v1->y()) + + (v0->z() - v1->z()) * (v0->z() - v1->z())); tot_length += length; } } - MVertex *first_v = (*temp.begin())->getVertex(0); + MVertex *first_v = (*temp.begin())->getVertex(0); MVertex *current_v = (*temp.begin())->getVertex(1); l.push_back(first_v); coord.push_back(0.0); temp.erase(temp.begin()); - - // printf("SORTED %d -> %d \n",first_v->getNum(),current_v->getNum()); - while (temp.size()){ bool found = 0; - for (std::list<MLine*>::iterator itl = temp.begin() ; itl !=temp.end() ; ++itl){ + for (std::list<MLine*>::iterator itl = temp.begin(); itl != temp.end(); ++itl){ MLine *ll = *itl; MVertex *v0 = ll->getVertex(0); MVertex *v1 = ll->getVertex(1); @@ -170,11 +155,10 @@ static bool orderVertices (const std::list<GEdge*> &e , l.push_back(current_v); current_v = v1; temp.erase(itl); - const double length = sqrt ((v0->x() - v1->x()) * (v0->x() - v1->x()) + - (v0->y() - v1->y()) * (v0->y() - v1->y()) + - (v0->z() - v1->z()) * (v0->z() - v1->z()) ); - coord.push_back(coord[coord.size()-1] + length/tot_length); - // printf("SORTED %d -> %d (%d,%g)\n",v0->getNum(),v1->getNum(),temp.size(),coord[coord.size()-1]); + const double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) + + (v0->y() - v1->y()) * (v0->y() - v1->y()) + + (v0->z() - v1->z()) * (v0->z() - v1->z())); + coord.push_back(coord[coord.size()-1] + length / tot_length); break; } else if (v1 == current_v){ @@ -182,23 +166,19 @@ static bool orderVertices (const std::list<GEdge*> &e , l.push_back(current_v); current_v = v0; temp.erase(itl); - const double length = sqrt ((v0->x() - v1->x()) * (v0->x() - v1->x()) + - (v0->y() - v1->y()) * (v0->y() - v1->y()) + - (v0->z() - v1->z()) * (v0->z() - v1->z()) ); - coord.push_back(coord[coord.size()-1] + length/tot_length); - // printf("SORTED %d -> %d (%d,%g)\n",v1->getNum(),v0->getNum(),temp.size(),coord[coord.size()-1]); + const double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) + + (v0->y() - v1->y()) * (v0->y() - v1->y()) + + (v0->z() - v1->z()) * (v0->z() - v1->z())); + coord.push_back(coord[coord.size()-1] + length / tot_length); break; } } - if (!found)return false; + if(!found) return false; } - // printf(" DONE\n"); - return true; } - -void GFaceCompound::parametrize (bool _isU, int ITER) const +void GFaceCompound::parametrize(bool _isU, int ITER) const { Msg::Info("Parametrizing Surface %d coordinate %d (ITER %d)", tag(), _isU, ITER); @@ -210,23 +190,21 @@ void GFaceCompound::parametrize (bool _isU, int ITER) const gmshLinearSystemFull lsys; #endif gmshAssembler myAssembler(&lsys); - gmshGradientBasedDiffusivity diffusivity (coordinates); + gmshGradientBasedDiffusivity diffusivity(coordinates); if (ITER > 0) diffusivity.setComponent(_isU); if (!_U1.size()){ // maps the boundary onto a circle std::vector<MVertex*> ordered; std::vector<double> coords; - // Msg::Info("%d edges on the contour", l_edges.size()); - // for (std::list<GEdge*>::const_iterator it = l_edges.begin(); - // it !=l_edges.end();++it)printf("%d ",(*it)->tag()); - // printf("\n"); - bool success = orderVertices (l_edges, ordered, coords); - if (!success)throw; + bool success = orderVertices(l_edges, ordered, coords); + if (!success){ + Msg::Error("Could not order vertices on boundary"); + return; + } for (unsigned int i = 0; i < ordered.size(); i++){ MVertex *v = ordered[i]; const double theta = 2 * M_PI * coords[i]; - // printf("fixing %d to %g\n",v->getIndex(),theta); if (_isU) myAssembler.fixVertex(v, 0, 1, cos(theta)); else myAssembler.fixVertex(v, 0, 1, sin(theta)); } @@ -234,52 +212,55 @@ void GFaceCompound::parametrize (bool _isU, int ITER) const else{ if (_isU){ std::list<GEdge*> :: const_iterator it = _U0.begin(); - for ( ; it != _U0.end() ; ++it){ - fixEdgeToValue (*it, 0.0, myAssembler); + for ( ; it != _U0.end(); ++it){ + fixEdgeToValue(*it, 0.0, myAssembler); } it = _U1.begin(); - for ( ; it != _U1.end() ; ++it){ - fixEdgeToValue (*it, 1.0, myAssembler); + for ( ; it != _U1.end(); ++it){ + fixEdgeToValue(*it, 1.0, myAssembler); } } else{ std::list<GEdge*> :: const_iterator it = _V0.begin(); - for ( ; it != _V0.end() ; ++it){ + for ( ; it != _V0.end(); ++it){ fixEdgeToValue (*it, 0.0, myAssembler); } it = _V1.begin(); - for ( ; it != _V1.end() ; ++it){ + for ( ; it != _V1.end(); ++it){ fixEdgeToValue (*it, 1.0, myAssembler); } } } - std::list<GFace*> :: const_iterator it = _compound.begin(); - for ( ; it != _compound.end() ; ++it){ + std::list<GFace*>::const_iterator it = _compound.begin(); + for ( ; it != _compound.end(); ++it){ for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){ MTriangle *t = (*it)->triangles[i]; - myAssembler.numberVertex(t->getVertex(0),0, 1); - myAssembler.numberVertex(t->getVertex(1),0, 1); - myAssembler.numberVertex(t->getVertex(2),0, 1); + myAssembler.numberVertex(t->getVertex(0), 0, 1); + myAssembler.numberVertex(t->getVertex(1), 0, 1); + myAssembler.numberVertex(t->getVertex(2), 0, 1); } } - // printf("%d %d %d %d\n",_U0.size(),_U1.size(),_V0.size(),_V1.size()); - // printf("creating term %d dofs numbered %d fixed\n", myAssembler.sizeOfR(),myAssembler.sizeOfF()); + + Msg::Debug("Creating term %d dofs numbered %d fixed", + myAssembler.sizeOfR(), myAssembler.sizeOfF()); + gmshLaplaceTerm laplace (model(), &diffusivity, 1); it = _compound.begin(); for ( ; it != _compound.end() ; ++it){ for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){ MTriangle *t = (*it)->triangles[i]; diffusivity.setCurrent(t); - laplace.addToMatrix (myAssembler, t); + laplace.addToMatrix(myAssembler, t); } } - // printf("assembly done\n"); + + Msg::Debug("Assembly done"); lsys.systemSolve(); - // printf("system souleved\n"); + Msg::Debug("System solved"); it = _compound.begin(); - for ( ; it != _compound.end() ; ++it){ + for ( ; it != _compound.end(); ++it){ for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){ MTriangle *t = (*it)->triangles[i]; double uu[3], vv[3]; @@ -288,7 +269,7 @@ void GFaceCompound::parametrize (bool _isU, int ITER) const double value = myAssembler.getDofValue(v,0,1); std::map<MVertex*,SPoint2>::iterator itf = coordinates.find(v); if (itf == coordinates.end()){ - SPoint2 p (_isU ? value:0.0,_isU ? 0.0:value); + SPoint2 p (_isU ? value : 0.0, _isU ? 0.0 : value); coordinates[v] = p; } else{ @@ -297,12 +278,11 @@ void GFaceCompound::parametrize (bool _isU, int ITER) const uu[j] = itf->second[0]; vv[j] = itf->second[1]; } - } + } } } } - void GFaceCompound::computeNormals () const { _normals.clear(); @@ -312,7 +292,7 @@ void GFaceCompound::computeNormals () const for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){ MTriangle *t = (*it)->triangles[i]; t->getJacobian(0, 0, 0, J); - // SVector3 n (J[2][0],J[2][1],J[2][2]); + // SVector3 n (J[2][0],J[2][1],J[2][2]); SVector3 d1(J[0][0], J[0][1], J[0][2]); SVector3 d2(J[1][0], J[1][1], J[1][2]); SVector3 n = crossprod(d1, d2); @@ -351,79 +331,66 @@ double GFaceCompound::curvature(MTriangle *t) const return fabs(t->interpolateDiv (val,0,0,0.)); } - GPoint GFaceCompound::point(double par1, double par2) const { parametrize(); double U,V; GFaceCompoundTriangle *lt; getTriangle (par1, par2, <, U,V); - SPoint3 p(0,0,0); + SPoint3 p(0, 0, 0); if (!lt){ - Msg::Warning("Re-Parametrized face %d --> point (%g %g) lies outside the domain", tag(), par1,par2); - - return GPoint(p.x(),p.y(),p.z(),this); + Msg::Warning("Re-Parametrized face %d --> point (%g %g) lies outside the domain", + tag(), par1,par2); + return GPoint(p.x(), p.y(), p.z(), this); } - p = lt->v1*(1.-U-V) + lt->v2*U + lt->v3*V; - // lt->t->pnt(U,V,0,p); - double par[2] = {par1,par2}; - // printf("coucou %g %g -> %g %g %g\n",par1,par2,p.x(),p.y(),p.z()); - return GPoint(p.x(),p.y(),p.z(),this,par); + p = lt->v1 * (1. - U - V) + lt->v2 * U + lt->v3 * V; + double par[2] = {par1, par2}; + return GPoint(p.x(), p.y(), p.z(), this, par); } -/* - computing dX/du and dX/dv -*/ - Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const { parametrize(); double U,V,UDU,UDV,VDU,VDV; GFaceCompoundTriangle *lt; - // GFaceCompoundTriangle *ltdu; - // GFaceCompoundTriangle *ltdv; - getTriangle (param.x(), param.y(), <, U,V); - // getTriangle (param.x()+1.e-4, param.y(), <du, UDU,VDU); - // getTriangle (param.x(), param.y()+1.e-4, <dv, UDV,VDV); + // GFaceCompoundTriangle *ltdu; + // GFaceCompoundTriangle *ltdv; + getTriangle(param.x(), param.y(), <, U,V); + // getTriangle (param.x()+1.e-4, param.y(), <du, UDU,VDU); + // getTriangle (param.x(), param.y()+1.e-4, <dv, UDV,VDV); if (!lt) - return Pair<SVector3, SVector3>(SVector3(1,0,0),SVector3(0,1,0)); + return Pair<SVector3, SVector3>(SVector3(1, 0, 0), SVector3(0, 1, 0)); - double mat[2][2] = {{lt->p2.x()-lt->p1.x(),lt->p3.x()-lt->p1.x()}, - {lt->p2.y()-lt->p1.y(),lt->p3.y()-lt->p1.y()}}; + double mat[2][2] = {{lt->p2.x() - lt->p1.x(), lt->p3.x() - lt->p1.x()}, + {lt->p2.y() - lt->p1.y(), lt->p3.y() - lt->p1.y()}}; double inv[2][2]; inv2x2(mat,inv); - // MVertex *v0 = lt->t->getVertex(0); - // MVertex *v1 = lt->t->getVertex(1); - // MVertex *v2 = lt->t->getVertex(2); + // MVertex *v0 = lt->t->getVertex(0); + // MVertex *v1 = lt->t->getVertex(1); + // MVertex *v2 = lt->t->getVertex(2); - // SPoint3 p(0,0,0); - // SPoint3 pdu(0,0,0); - // SPoint3 pdv(0,0,0); + // SPoint3 p(0,0,0); + // SPoint3 pdu(0,0,0); + // SPoint3 pdv(0,0,0); - // lt->t->pnt(U,V,0,p); - // ltdu->t->pnt(UDU,VDU,0,pdu); - // ltdv->t->pnt(UDV,VDV,0,pdv); + // lt->t->pnt(U,V,0,p); + // ltdu->t->pnt(UDU,VDU,0,pdu); + // ltdv->t->pnt(UDV,VDV,0,pdv); - SVector3 dXdxi (lt->v2 - lt->v1); SVector3 dXdeta (lt->v3 - lt->v1); - // return Pair<SVector3, SVector3>(dXdxi*inv[0][0]+dXdeta*inv[1][0], - // dXdxi*inv[0][1]+dXdeta*inv[1][1]); + // return Pair<SVector3, SVector3>(dXdxi * inv[0][0] + dXdeta * inv[1][0], + // dXdxi * inv[0][1] + dXdeta * inv[1][1]); - SVector3 dXdu (dXdxi*inv[0][0]+dXdeta*inv[1][0]); - SVector3 dXdv (dXdxi*(inv[0][1])+dXdeta*(inv[1][1])); + SVector3 dXdu (dXdxi * inv[0][0] + dXdeta * inv[1][0]); + SVector3 dXdv (dXdxi * inv[0][1] + dXdeta * inv[1][1]); - // SVector3 dXduFD (SVector3(pdu-p)*1.e4); - // SVector3 dXdvFD (SVector3(pdv-p)*1.e4); - -// printf("FD %g %g %g / %g %g %g AN %g %g %g / %g %g %g \n", -// dXdu.x(),dXdu.y(),dXdu.z(), -// dXdv.x(),dXdv.y(),dXdv.z(), -// dXduFD.x(),dXduFD.y(),dXduFD.z(), -// dXdvFD.x(),dXdvFD.y(),dXdvFD.z()); + // SVector3 dXduFD(SVector3(pdu - p) * 1.e4); + // SVector3 dXdvFD(SVector3(pdv - p) * 1.e4); + return Pair<SVector3, SVector3>(dXdu,dXdv); } @@ -448,8 +415,8 @@ static void GFaceCompoundBB(void *a, double*mmin, double*mmax) static void GFaceCompoundCentroid(void *a, double*c) { GFaceCompoundTriangle *t = (GFaceCompoundTriangle *)a; - c[0] = (t->p1.x() + t->p2.x() + t->p3.x())/3.0; - c[1] = (t->p1.y() + t->p2.y() + t->p3.y())/3.0; + c[0] = (t->p1.x() + t->p2.x() + t->p3.x()) / 3.0; + c[1] = (t->p1.y() + t->p2.y() + t->p3.y()) / 3.0; c[2] = 0.0; } @@ -474,15 +441,15 @@ static int GFaceCompoundInEle(void *a, double*c) return 0; } -void GFaceCompound::getTriangle (double u, double v, - GFaceCompoundTriangle **lt, - double &_u, double &_v) const +void GFaceCompound::getTriangle(double u, double v, + GFaceCompoundTriangle **lt, + double &_u, double &_v) const { parametrize(); - double uv[3] = {u,v,0}; + double uv[3] = {u, v, 0}; *lt = (GFaceCompoundTriangle*)Octree_Search(uv, oct); - if (!(*lt)){return;} + if (!(*lt)) return; double M[2][2],X[2],R[2]; const SPoint2 p0 = (*lt)->p1; @@ -497,7 +464,6 @@ void GFaceCompound::getTriangle (double u, double v, sys2x2(M,R,X); _u = X[0]; _v = X[1]; - return; } void GFaceCompound::buildOct() const @@ -505,7 +471,7 @@ void GFaceCompound::buildOct() const SBoundingBox3d bb; int count = 0; std::list<GFace*> :: const_iterator it = _compound.begin(); - // printf("building octree %d coords \n",coordinates.size()); + for ( ; it != _compound.end() ; ++it){ for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){ MTriangle *t = (*it)->triangles[i]; @@ -518,19 +484,19 @@ void GFaceCompound::buildOct() const } } _gfct = new GFaceCompoundTriangle[count]; - double origin[3] = {bb.min().x(),bb.min().y(),bb.min().z()}; - double ssize[3] = {bb.max().x()-bb.min().x(), - bb.max().y()-bb.min().y(), - bb.max().z()-bb.min().z()}; + double origin[3] = {bb.min().x(), bb.min().y(), bb.min().z()}; + double ssize[3] = {bb.max().x() - bb.min().x(), + bb.max().y() - bb.min().y(), + bb.max().z() - bb.min().z()}; oct = Octree_Create(10, origin, ssize, GFaceCompoundBB, GFaceCompoundCentroid, GFaceCompoundInEle); it = _compound.begin(); count = 0; - FILE * uvx = fopen ("UVX.pos","w"); - FILE * uvy = fopen ("UVY.pos","w"); - FILE * uvz = fopen ("UVZ.pos","w"); + FILE * uvx = fopen("UVX.pos","w"); + FILE * uvy = fopen("UVY.pos","w"); + FILE * uvz = fopen("UVZ.pos","w"); FILE * xyzu = fopen("XYZU.pos","w"); FILE * xyzv = fopen("XYZV.pos","w"); FILE * xyzc = fopen("XYZC.pos","w"); @@ -555,59 +521,45 @@ void GFaceCompound::buildOct() const _gfct[count].p2 = it1->second; _gfct[count].p3 = it2->second; _gfct[count].t = t; - _gfct[count].v1 = SPoint3(t->getVertex(0)->x(),t->getVertex(0)->y(),t->getVertex(0)->z()); - _gfct[count].v2 = SPoint3(t->getVertex(1)->x(),t->getVertex(1)->y(),t->getVertex(1)->z()); - _gfct[count].v3 = SPoint3(t->getVertex(2)->x(),t->getVertex(2)->y(),t->getVertex(2)->z()); + _gfct[count].v1 = SPoint3 + (t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z()); + _gfct[count].v2 = SPoint3 + (t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z()); + _gfct[count].v3 = SPoint3 + (t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z()); fprintf(xyzu,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n", - t->getVertex(0)->x(), - t->getVertex(0)->y(), - t->getVertex(0)->z(), - t->getVertex(1)->x(), - t->getVertex(1)->y(), - t->getVertex(1)->z(), - t->getVertex(2)->x(), - t->getVertex(2)->y(), - t->getVertex(2)->z(), + t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(), + t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(), + t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(), it0->second.x(),it1->second.x(),it2->second.x()); fprintf(xyzv,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n", - t->getVertex(0)->x(), - t->getVertex(0)->y(), - t->getVertex(0)->z(), - t->getVertex(1)->x(), - t->getVertex(1)->y(), - t->getVertex(1)->z(), - t->getVertex(2)->x(), - t->getVertex(2)->y(), - t->getVertex(2)->z(), + t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(), + t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(), + t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(), it0->second.y(),it1->second.y(),it2->second.y()); const double K = fabs(curvature (t)); fprintf(xyzc,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n", - t->getVertex(0)->x(), - t->getVertex(0)->y(), - t->getVertex(0)->z(), - t->getVertex(1)->x(), - t->getVertex(1)->y(), - t->getVertex(1)->z(), - t->getVertex(2)->x(), - t->getVertex(2)->y(), - t->getVertex(2)->z(),K,K,K); + t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(), + t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(), + t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(), + K, K, K); fprintf(uvx,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n", - it0->second.x(),it0->second.y(),0.0, - it1->second.x(),it1->second.y(),0.0, - it2->second.x(),it2->second.y(),0.0, - t->getVertex(0)->x(),t->getVertex(1)->x(),t->getVertex(2)->x()); + it0->second.x(), it0->second.y(), 0.0, + it1->second.x(), it1->second.y(), 0.0, + it2->second.x(), it2->second.y(), 0.0, + t->getVertex(0)->x(), t->getVertex(1)->x(), t->getVertex(2)->x()); fprintf(uvy,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n", - it0->second.x(),it0->second.y(),0.0, - it1->second.x(),it1->second.y(),0.0, - it2->second.x(),it2->second.y(),0.0, - t->getVertex(0)->y(),t->getVertex(1)->y(),t->getVertex(2)->y()); + it0->second.x(), it0->second.y(), 0.0, + it1->second.x(), it1->second.y(), 0.0, + it2->second.x(), it2->second.y(), 0.0, + t->getVertex(0)->y(), t->getVertex(1)->y(), t->getVertex(2)->y()); fprintf(uvz,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n", - it0->second.x(),it0->second.y(),0.0, - it1->second.x(),it1->second.y(),0.0, - it2->second.x(),it2->second.y(),0.0, - t->getVertex(0)->z(),t->getVertex(1)->z(),t->getVertex(2)->z()); - + it0->second.x(), it0->second.y(), 0.0, + it1->second.x(), it1->second.y(), 0.0, + it2->second.x(), it2->second.y(), 0.0, + t->getVertex(0)->z(), t->getVertex(1)->z(), t->getVertex(2)->z()); + Octree_Insert(&_gfct[count], oct); count ++; } @@ -626,4 +578,3 @@ void GFaceCompound::buildOct() const fclose(xyzc); Octree_Arrange(oct); } - diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h index d7db76c451ccb6a7f24e086ac39275ceadf02a8d..2e72276e176a1c7c39ec7ed06ae8754e83e5a26c 100644 --- a/Geo/GFaceCompound.h +++ b/Geo/GFaceCompound.h @@ -47,7 +47,8 @@ class GFaceCompound : public GFace { void parametrize() const ; void computeNormals () const; void getBoundingEdges(); - void getTriangle(double u, double v, GFaceCompoundTriangle **lt, double &_u, double &_v) const; + void getTriangle(double u, double v, GFaceCompoundTriangle **lt, + double &_u, double &_v) const; virtual double curvature(MTriangle *t) const; public: GFaceCompound(GModel *m, int tag, @@ -66,7 +67,6 @@ public: SPoint2 getCoordinates(MVertex *v) const { parametrize() ; return coordinates[v]; } virtual bool buildRepresentationCross(){ return false; } virtual double curvature(const SPoint2 ¶m) const; - }; #endif