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

*** empty log message ***

parent 3eb1bfd5
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/GFaceCompound.cpp
+ 91
− 66
View file @ 86a6ce05
...@@ -66,12 +66,22 @@ void GFaceCompound::parametrize() const ...@@ -66,12 +66,22 @@ void GFaceCompound::parametrize() const
void GFaceCompound::getBoundingEdges() void GFaceCompound::getBoundingEdges()
{ {
if (_U0.size() && !_U1.size()){ if (_U0.size() && !_V1.size()){
std::list<GEdge*> :: const_iterator it = _U0.begin(); std::list<GEdge*> :: const_iterator it = _U0.begin();
for ( ; it != _U0.end() ; ++it){ for ( ; it != _U0.end() ; ++it){
l_edges.push_back(*it); l_edges.push_back(*it);
(*it)->addFace(this); (*it)->addFace(this);
} }
it = _U1.begin();
for ( ; it != _U1.end() ; ++it){
l_edges.push_back(*it);
(*it)->addFace(this);
}
it = _V0.begin();
for ( ; it != _V0.end() ; ++it){
l_edges.push_back(*it);
(*it)->addFace(this);
}
return; return;
} }
...@@ -107,6 +117,11 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound, ...@@ -107,6 +117,11 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
: GFace(m, tag), _compound(compound), _U0(U0), _U1(U1), _V0(V0), _V1(V1), oct(0) : GFace(m, tag), _compound(compound), _U0(U0), _U1(U1), _V0(V0), _V1(V1), oct(0)
{ {
getBoundingEdges(); getBoundingEdges();
if (!_U0.size()) _type = UNITCIRCLE;
else if (!_V0.size()) _type = UNITCIRCLE;
else if (!_U1.size()) _type = CYLINDER;
else if (!_V1.size()) _type = BIFURCATION;
else _type = SQUARE;
} }
GFaceCompound::~GFaceCompound() GFaceCompound::~GFaceCompound()
...@@ -193,7 +208,7 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const ...@@ -193,7 +208,7 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
gmshGradientBasedDiffusivity diffusivity(coordinates); gmshGradientBasedDiffusivity diffusivity(coordinates);
if (ITER > 0) diffusivity.setComponent(_isU); if (ITER > 0) diffusivity.setComponent(_isU);
if (!_U1.size()){ if (_type == UNITCIRCLE){
// maps the boundary onto a circle // maps the boundary onto a circle
std::vector<MVertex*> ordered; std::vector<MVertex*> ordered;
std::vector<double> coords; std::vector<double> coords;
...@@ -209,6 +224,19 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const ...@@ -209,6 +224,19 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
else myAssembler.fixVertex(v, 0, 1, sin(theta)); else myAssembler.fixVertex(v, 0, 1, sin(theta));
} }
} }
else if (_type == CYLINDER){
// maps the boundary onto an annulus
std::vector<MVertex*> ordered;
std::vector<double> coords;
bool success = orderVertices (_U0, ordered, coords);
if (!success)throw;
for (unsigned int i = 0; i < ordered.size(); i++){
MVertex *v = ordered[i];
const double theta = 2 * M_PI * coords[i];
if (_isU) myAssembler.fixVertex(v, 0, 1, cos(theta));
else myAssembler.fixVertex(v, 0, 1, sin(theta));
}
}
else{ else{
if (_isU){ if (_isU){
std::list<GEdge*> :: const_iterator it = _U0.begin(); std::list<GEdge*> :: const_iterator it = _U0.begin();
...@@ -258,14 +286,19 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const ...@@ -258,14 +286,19 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
Msg::Debug("Assembly done"); Msg::Debug("Assembly done");
lsys.systemSolve(); lsys.systemSolve();
Msg::Debug("System solved"); Msg::Debug("System solved");
it = _compound.begin(); it = _compound.begin();
std::set<MVertex *>allNodes;
for ( ; it != _compound.end() ; ++it){ for ( ; it != _compound.end() ; ++it){
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){ for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i]; MTriangle *t = (*it)->triangles[i];
double uu[3], vv[3];
for (int j = 0; j < 3; j++){ for (int j = 0; j < 3; j++){
MVertex *v = t->getVertex(j); allNodes.insert(t->getVertex(j));
}
}
}
for (std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
double value = myAssembler.getDofValue(v,0,1); double value = myAssembler.getDofValue(v,0,1);
std::map<MVertex*,SPoint2>::iterator itf = coordinates.find(v); std::map<MVertex*,SPoint2>::iterator itf = coordinates.find(v);
if (itf == coordinates.end()){ if (itf == coordinates.end()){
...@@ -275,10 +308,6 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const ...@@ -275,10 +308,6 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
else{ else{
if(_isU) itf->second[0]= value; if(_isU) itf->second[0]= value;
else itf->second[1]= value; else itf->second[1]= value;
uu[j] = itf->second[0];
vv[j] = itf->second[1];
}
}
} }
} }
} }
...@@ -317,7 +346,19 @@ double GFaceCompound::curvature(const SPoint2 &param) const ...@@ -317,7 +346,19 @@ double GFaceCompound::curvature(const SPoint2 &param) const
if (!lt){ if (!lt){
return 0.0; return 0.0;
} }
return curvature(lt->t);
// if (lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
// SPoint2 pv1, pv2, pv3;
// bool ok = reparamMeshVertexOnFace(lt->t->getVertex(0), lt->gf, pv1);
// ok |= reparamMeshVertexOnFace(lt->t->getVertex(1), lt->gf, pv2);
// ok |= reparamMeshVertexOnFace(lt->t->getVertex(2), lt->gf, pv3);
// if (ok){
// SPoint2 pv = pv1*(1.-U-V) + pv2*U + pv3*V;
// return lt->gf->curvature(pv));
// }
// }
// return curvature(lt->t);
} }
double GFaceCompound::curvature(MTriangle *t) const double GFaceCompound::curvature(MTriangle *t) const
...@@ -328,7 +369,7 @@ double GFaceCompound::curvature(MTriangle *t) const ...@@ -328,7 +369,7 @@ double GFaceCompound::curvature(MTriangle *t) const
double val[9] = {n1.x(),n2.x(),n3.x(), double val[9] = {n1.x(),n2.x(),n3.x(),
n1.y(),n2.y(),n3.y(), n1.y(),n2.y(),n3.y(),
n1.z(),n2.z(),n3.z()}; n1.z(),n2.z(),n3.z()};
return fabs(t->interpolateDiv (val,0,0,0.)); // return fabs(t->interpolateDiv (val,0.,0.,0.));
} }
GPoint GFaceCompound::point(double par1, double par2) const GPoint GFaceCompound::point(double par1, double par2) const
...@@ -339,10 +380,15 @@ GPoint GFaceCompound::point(double par1, double par2) const ...@@ -339,10 +380,15 @@ GPoint GFaceCompound::point(double par1, double par2) const
getTriangle (par1, par2, &lt, U,V); getTriangle (par1, par2, &lt, U,V);
SPoint3 p(0, 0, 0); SPoint3 p(0, 0, 0);
if (!lt){ if (!lt){
Msg::Warning("Re-Parametrized face %d --> point (%g %g) lies outside the domain", Msg::Warning("Re-Parametrized face %d --> point (%g %g) lies outside the domain", tag(), par1,par2);
tag(), par1,par2);
return GPoint(p.x(),p.y(),p.z(),this); return GPoint(p.x(),p.y(),p.z(),this);
} }
if (0 && lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
return lt->gf->point(pv.x(),pv.y());
}
p = lt->v1*(1.-U-V) + lt->v2*U + lt->v3*V; p = lt->v1*(1.-U-V) + lt->v2*U + lt->v3*V;
double par[2] = {par1,par2}; double par[2] = {par1,par2};
return GPoint(p.x(),p.y(),p.z(),this,par); return GPoint(p.x(),p.y(),p.z(),this,par);
...@@ -353,12 +399,7 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 &param) const ...@@ -353,12 +399,7 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 &param) const
parametrize(); parametrize();
double U,V,UDU,UDV,VDU,VDV; double U,V,UDU,UDV,VDU,VDV;
GFaceCompoundTriangle *lt; GFaceCompoundTriangle *lt;
// GFaceCompoundTriangle *ltdu;
// GFaceCompoundTriangle *ltdv;
getTriangle (param.x(), param.y(), &lt, U,V); getTriangle (param.x(), param.y(), &lt, U,V);
// getTriangle (param.x()+1.e-4, param.y(), &ltdu, UDU,VDU);
// getTriangle (param.x(), param.y()+1.e-4, &ltdv, UDV,VDV);
if (!lt) 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));
...@@ -367,30 +408,12 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 &param) const ...@@ -367,30 +408,12 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 &param) const
double inv[2][2]; double inv[2][2];
inv2x2(mat,inv); inv2x2(mat,inv);
// 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);
// 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 dXdxi (lt->v2 - lt->v1);
SVector3 dXdeta (lt->v3 - 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]);
SVector3 dXdu (dXdxi*inv[0][0]+dXdeta*inv[1][0]); SVector3 dXdu (dXdxi*inv[0][0]+dXdeta*inv[1][0]);
SVector3 dXdv (dXdxi*inv[0][1]+dXdeta*inv[1][1]); SVector3 dXdv (dXdxi*inv[0][1]+dXdeta*inv[1][1]);
// SVector3 dXduFD(SVector3(pdu - p) * 1.e4);
// SVector3 dXdvFD(SVector3(pdv - p) * 1.e4);
return Pair<SVector3, SVector3>(dXdu,dXdv); return Pair<SVector3, SVector3>(dXdu,dXdv);
} }
...@@ -520,13 +543,15 @@ void GFaceCompound::buildOct() const ...@@ -520,13 +543,15 @@ void GFaceCompound::buildOct() const
_gfct[count].p1 = it0->second; _gfct[count].p1 = it0->second;
_gfct[count].p2 = it1->second; _gfct[count].p2 = it1->second;
_gfct[count].p3 = it2->second; _gfct[count].p3 = it2->second;
_gfct[count].t = t; if ((*it)->geomType() != GEntity::DiscreteSurface){
_gfct[count].v1 = SPoint3 reparamMeshVertexOnFace(t->getVertex(0), *it, _gfct[count].gfp1);
(t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z()); reparamMeshVertexOnFace(t->getVertex(1), *it, _gfct[count].gfp2);
_gfct[count].v2 = SPoint3 reparamMeshVertexOnFace(t->getVertex(2), *it, _gfct[count].gfp3);
(t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z()); }
_gfct[count].v3 = SPoint3 _gfct[count].v1 = SPoint3(t->getVertex(0)->x(),t->getVertex(0)->y(),t->getVertex(0)->z());
(t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->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].gf = *it;
fprintf(xyzu,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n", 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(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(), t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
......
Geo/GFaceCompound.h
+ 16
− 6
View file @ 86a6ce05
...@@ -26,11 +26,15 @@ The compound can therefore be re-meshed using any surface mesh ...@@ -26,11 +26,15 @@ The compound can therefore be re-meshed using any surface mesh
generator of gmsh! generator of gmsh!
*/ */
typedef struct { class GFaceCompoundTriangle {
public:
SPoint2 p1, p2, p3; SPoint2 p1, p2, p3;
SPoint2 gfp1, gfp2, gfp3;
SPoint3 v1, v2, v3; SPoint3 v1, v2, v3;
MTriangle *t; GFace *gf;
} GFaceCompoundTriangle; GFaceCompoundTriangle () : gf(0)
{}
} ;
class Octree; class Octree;
...@@ -51,9 +55,13 @@ class GFaceCompound : public GFace { ...@@ -51,9 +55,13 @@ class GFaceCompound : public GFace {
double &_u, double &_v) const; double &_u, double &_v) const;
virtual double curvature(MTriangle *t) const; virtual double curvature(MTriangle *t) const;
public: public:
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound, typedef enum {UNITCIRCLE, CYLINDER, BIFURCATION, SQUARE} typeOfIsomorphism;
std::list<GEdge*> &U0, std::list<GEdge*> &U1, GFaceCompound(GModel *m, int tag,
std::list<GEdge*> &V0, std::list<GEdge*> &V1); std::list<GFace*> &compound,
std::list<GEdge*> &U0,
std::list<GEdge*> &U1,
std::list<GEdge*> &V0,
std::list<GEdge*> &V1);
virtual ~GFaceCompound(); virtual ~GFaceCompound();
Range<double> parBounds(int i) const { return Range<double>(0, 1); } Range<double> parBounds(int i) const { return Range<double>(0, 1); }
virtual GPoint point(double par1, double par2) const; virtual GPoint point(double par1, double par2) const;
...@@ -64,6 +72,8 @@ public: ...@@ -64,6 +72,8 @@ public:
SPoint2 getCoordinates(MVertex *v) const { parametrize() ; return coordinates[v]; } SPoint2 getCoordinates(MVertex *v) const { parametrize() ; return coordinates[v]; }
virtual bool buildRepresentationCross(){ return false; } virtual bool buildRepresentationCross(){ return false; }
virtual double curvature(const SPoint2 &param) const; virtual double curvature(const SPoint2 &param) const;
private:
typeOfIsomorphism _type;
}; };
#endif #endif
Geo/MVertex.cpp
+ 9
− 0
View file @ 86a6ce05
...@@ -208,6 +208,14 @@ MVertex::linearSearch(std::set<MVertex*, MVertexLessThanLexicographic> &pos) ...@@ -208,6 +208,14 @@ MVertex::linearSearch(std::set<MVertex*, MVertexLessThanLexicographic> &pos)
static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> &params) static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> &params)
{ {
params.clear(); params.clear();
if (gf->geomType() == GEntity::CompoundSurface &&
v->onWhat()->dim() < 2){
GFaceCompound *gfc = (GFaceCompound*) gf;
params.push_back(gfc->getCoordinates(v));
return;
}
if(v->onWhat()->dim() == 0){ if(v->onWhat()->dim() == 0){
GVertex *gv = (GVertex*)v->onWhat(); GVertex *gv = (GVertex*)v->onWhat();
std::list<GEdge*> ed = gv->edges(); std::list<GEdge*> ed = gv->edges();
...@@ -250,6 +258,7 @@ static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> &params ...@@ -250,6 +258,7 @@ static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> &params
bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf, bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
SPoint2 &param1, SPoint2 &param2) SPoint2 &param1, SPoint2 &param2)
{ {
std::vector<SPoint2> p1, p2; std::vector<SPoint2> p1, p2;
getAllParameters(v1, gf, p1); getAllParameters(v1, gf, p1);
getAllParameters(v2, gf, p2); getAllParameters(v2, gf, p2);
......
Mesh/Refine.cpp
+ 14
− 6
View file @ 86a6ce05
...@@ -86,13 +86,21 @@ static void Subdivide(GFace *gf, bool splitTrianglesIntoQuads) ...@@ -86,13 +86,21 @@ static void Subdivide(GFace *gf, bool splitTrianglesIntoQuads)
MTriangle *t = gf->triangles[i]; MTriangle *t = gf->triangles[i];
if(t->getNumVertices() == 6){ if(t->getNumVertices() == 6){
SPoint2 pt, temp; SPoint2 pt, temp;
SPoint3 ptx; t->pnt(0.5,0.5,0,ptx);
bool reparamOK = true;
for(int k = 0; k<6; k++){ for(int k = 0; k<6; k++){
reparamMeshVertexOnFace(t->getVertex(k), gf, temp); reparamOK &= reparamMeshVertexOnFace(t->getVertex(k), gf, temp);
pt[0] += temp[0]/6.; pt[0] += temp[0]/6.;
pt[1] += temp[1]/6.; pt[1] += temp[1]/6.;
} }
MVertex *newv;
if (reparamOK){
GPoint gp = gf->point(pt); GPoint gp = gf->point(pt);
MFaceVertex *newv = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, pt[0], pt[1]); newv = new MFaceVertex (gp.x(),gp.y(),gp.z(),gf,pt[0],pt[1]);
}
else {
newv = new MVertex (ptx.x(),ptx.y(),ptx.z(),gf);
}
gf->mesh_vertices.push_back(newv); gf->mesh_vertices.push_back(newv);
quadrangles2.push_back quadrangles2.push_back
(new MQuadrangle(t->getVertex(0), t->getVertex(3), newv, t->getVertex(5))); (new MQuadrangle(t->getVertex(0), t->getVertex(3), newv, t->getVertex(5)));
......
benchmarks/3d/Torus.geo
+ 1
− 0
View file @ 86a6ce05
...@@ -13,3 +13,4 @@ Line Loop(5) = {4,1,2,3}; ...@@ -13,3 +13,4 @@ Line Loop(5) = {4,1,2,3};
Plane Surface(6) = {5}; Plane Surface(6) = {5};
Extrude Surface{6, {0.0,1,0}, {0,0.0,0.0}, 1*3.14159/2}; Extrude Surface{6, {0.0,1,0}, {0,0.0,0.0}, 1*3.14159/2};
Recombine Surface {6, 27, 23, 15, 19, 28};
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment