Commit 8d2147bb by Christophe Geuzaine

Merge branch 'new_reparametrization' into 'master'

New reparametrization See merge request !98
parents 4558e2c9 1de85893
Pipeline #1382 failed with stage
in 38 minutes 13 seconds
......@@ -55,7 +55,7 @@ windows64_msvc_ci:
script:
- md build
- cd build
- cmake -DGMSH_EXTRA_VERSION=%EXTRA_VERSION:~0,13% ..
- cmake -DGMSH_EXTRA_VERSION=%EXTRA_VERSION:~0,13% -DENABLE_HXT=0 ..
- msbuild package.vcxproj
tags:
- windows64
......
......@@ -52,6 +52,7 @@ opt(GETDP "Enable GetDP solver (as a directly linked library)" ${DEFAULT})
opt(GMM "Enable GMM linear solvers (simple alternative to PETSc)" ${DEFAULT})
opt(GMP "Enable GMP for Kbipack (advanced)" ON)
opt(GRAPHICS "Enable building graphics lib even without GUI (advanced)" OFF)
opt(HXT "Enable HXT library" ${DEFAULT})
opt(INTERNAL_DEVELOPER_API "Enable internal developer API" OFF)
opt(KBIPACK "Enable Kbipack (neeeded by homology solver)" ${DEFAULT})
opt(MATHEX "Enable math expression parser (used by plugins and options)" ${DEFAULT})
......@@ -1047,6 +1048,14 @@ if(HAVE_SOLVER)
endif(ENABLE_GETDP)
endif(HAVE_SOLVER)
if (ENABLE_HXT)
set_config_option(HAVE_HXT "Hxt")
add_subdirectory(contrib/hxt)
include_directories(BEFORE contrib/hxt)
endif(ENABLE_HXT)
if(ENABLE_OCC)
set(OCC_MINIMAL_VERSION "6.9.1")
if(WIN32 OR CYGWIN)
......
......@@ -21,6 +21,7 @@
#cmakedefine HAVE_FOURIER_MODEL
#cmakedefine HAVE_GMM
#cmakedefine HAVE_GMP
#cmakedefine HAVE_HXT
#cmakedefine HAVE_KBIPACK
#cmakedefine HAVE_LAPACK
#cmakedefine HAVE_LIBCGNS
......
......@@ -14,7 +14,7 @@ set(SRC
gmshSurface.cpp
OCCVertex.cpp OCCEdge.cpp OCCFace.cpp OCCRegion.cpp
GenericVertex.cpp GenericEdge.cpp GenericFace.cpp GenericRegion.cpp
discreteEdge.cpp discreteFace.cpp discreteDiskFace.cpp discreteRegion.cpp
discreteEdge.cpp discreteFace.cpp discreteRegion.cpp
fourierEdge.cpp fourierFace.cpp fourierProjectionFace.cpp
ACISVertex.cpp ACISEdge.cpp ACISFace.cpp
GModel.cpp
......
......@@ -15,13 +15,14 @@
#include "GaussLegendre1D.h"
#include "Context.h"
#include "closestPoint.h"
#include "discreteEdge.h"
#if defined(HAVE_MESH)
#include "meshGEdge.h"
#endif
GEdge::GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1)
: GEntity(model, tag), _length(0.), _tooSmall(false), _cp(0),
v0(_v0), v1(_v1), masterOrientation(0)
v0(_v0), v1(_v1), masterOrientation(0), compound_edge(NULL)
{
if(v0) v0->addEdge(this);
if(v1 && v1 != v0) v1->addEdge(this);
......@@ -663,10 +664,47 @@ void GEdge::discretize(double tol, std::vector<SPoint3> &dpts, std::vector<doubl
}
}
#if defined(HAVE_MESH)
static void meshCompound(GEdge* ge)
{
std::vector<MLine*> lines;
for (unsigned int i = 0; i < ge->_compound.size(); i++){
GEdge *c = (GEdge*)ge->_compound[i];
for (unsigned int j = 0; j<c->lines.size(); j++){
lines.push_back(new MLine(c->lines[j]->getVertex(0),
c->lines[j]->getVertex(1)));
}
}
discreteEdge *de = new discreteEdge(ge->model(), ge->tag() + 100000, NULL, NULL);
ge->model()->add(de);
de->lines = lines;
de->createGeometry();
de->mesh(false);
ge->compound_edge = de;
}
#endif
void GEdge::mesh(bool verbose)
{
#if defined(HAVE_MESH)
meshGEdge mesher;
mesher(this);
if(_compound.size()){ // Some faces are meshed together
if(_compound[0] == this){ // I'm the one that makes the compound job
bool ok = true;
for(unsigned int i = 0; i < _compound.size(); i++){
GEdge *ge = (GEdge*)_compound[i];
ok &= (ge->meshStatistics.status == GEdge::DONE);
}
if(!ok){
meshStatistics.status = GEdge::PENDING;
}
else{
meshCompound(this);
meshStatistics.status = GEdge::DONE;
return;
}
}
}
#endif
}
......@@ -216,6 +216,10 @@ class GEdge : public GEntity {
std::vector<MLine*> lines;
// when a compound of edges is created, both meshes should be kept alive this
// is due to Gmsh's flow and it only applies to model edges
GEdge *compound_edge;
void addLine(MLine *line){ lines.push_back(line); }
void addElement(int type, MElement *e);
void removeElement(int type, MElement *e);
......
......@@ -118,7 +118,7 @@ int GFace::delEdge(GEdge* edge)
pos++;
}
l_edges.erase(it);
std::list<int>::iterator itOri;
int posOri = 0;
int orientation = 0;
......@@ -130,7 +130,7 @@ int GFace::delEdge(GEdge* edge)
posOri++;
}
l_dirs.erase(itOri);
return orientation;
}
......@@ -1353,34 +1353,77 @@ bool GFace::fillPointCloud(double maxDist,
#if defined(HAVE_MESH)
static void meshCompound(GFace* gf, bool verbose)
{
discreteFace *df = new discreteFace(gf->model(), gf->tag() + 100000, true);
std::set<int> ec;
discreteFace *df = new discreteFace(gf->model(), gf->tag() + 100000);
std::vector< GFace* > triangles_tag;
std::vector< SPoint2 > triangles_uv;
for (unsigned int i = 0; i < gf->_compound.size(); i++){
GFace *c = (GFace*)gf->_compound[i];
std::list<GEdge*> edges = c->edges();
for (std::list<GEdge*>::iterator it = edges.begin() ; it != edges.end(); ++it){
std::set<int>::iterator found = ec.find((*it)->tag());
if (found == ec.end())ec.insert((*it)->tag());
else ec.erase(found);
}
df->triangles.insert(df->triangles.end(), c->triangles.begin(),
c->triangles.end());
c->triangles.end());
df->mesh_vertices.insert(df->mesh_vertices.end(), c->mesh_vertices.begin(),
c->mesh_vertices.end());
for (unsigned int j = 0; j < c->triangles.size(); j++)
df->_CAD.push_back(c);
c->mesh_vertices.end());
for (unsigned int j=0;j<c->triangles.size();j++){
triangles_tag.push_back(c);
for (int k = 0; k < 3; k++){
SPoint2 param;
reparamMeshVertexOnFace(c->triangles[j]->getVertex(k), c, param);
triangles_uv.push_back(param);
}
}
c->triangles.clear();
c->mesh_vertices.clear();
}
std::vector<int> cedges;
cedges.insert(cedges.begin(), ec.begin(), ec.end());
df->setBoundEdges(cedges);
df->createGeometry();
df->mesh(verbose);
for (int i = 0; i < df->mesh_vertices.size(); i++){
double u,v;
df->mesh_vertices[i]->getParameter(0, u);
df->mesh_vertices[i]->getParameter(1, v);
double U,V;
int position = df->trianglePosition(u, v, U, V);
if(position != -1) {
triangles_tag[position]->mesh_vertices.push_back(df->mesh_vertices[i]);
df->mesh_vertices[i]->setEntity(triangles_tag[position]);
if (0 && triangles_tag[position]->geomType() != GEntity::DiscreteSurface){
SPoint2 p0 = triangles_uv[3*position + 0];
SPoint2 p1 = triangles_uv[3*position + 1];
SPoint2 p2 = triangles_uv[3*position + 2];
SPoint2 p = p0 *(1-U-V) + p1 * U + p2 * V;
GPoint gp = triangles_tag[position]->point(p);
df->mesh_vertices[i]->setParameter(0,p.x());
df->mesh_vertices[i]->setParameter(1,p.y());
df->mesh_vertices[i]->x() = gp.x();
df->mesh_vertices[i]->y() = gp.y();
df->mesh_vertices[i]->z() = gp.z();
}
}
else {
df->mesh_vertices.push_back(df->mesh_vertices[i]);
df->mesh_vertices[i]->setEntity(gf);
}
}
for (int i = 0; i < df->triangles.size(); i++){
MTriangle *t = df->triangles[i];
if (t->getVertex(0)->onWhat()->dim() == 2)
((GFace*)t->getVertex(0)->onWhat())->triangles.push_back(t);
else if (t->getVertex(1)->onWhat()->dim() == 2)
((GFace*)t->getVertex(1)->onWhat())->triangles.push_back(t);
else if (t->getVertex(2)->onWhat()->dim() == 2)
((GFace*)t->getVertex(2)->onWhat())->triangles.push_back(t);
else gf->triangles.push_back(t); // FIXME could be better!
}
// gf->triangles = df->triangles;
df->triangles.clear();
df->mesh_vertices.clear();
delete df;
}
#endif
......@@ -1391,7 +1434,7 @@ void GFace::mesh(bool verbose)
meshGFace mesher;
mesher(this, verbose);
if(_compound.size()){ // Some faces are meshed together
if(_compound[0] == this){ // I'm the one that makes the compound job
if(_compound[0] == this){ // I'm the one that makes the compound job
bool ok = true;
for(unsigned int i = 0; i < _compound.size(); i++){
GFace *gf = (GFace*)_compound[i];
......@@ -1402,6 +1445,7 @@ void GFace::mesh(bool verbose)
}
else{
meshCompound(this, verbose);
meshStatistics.status = GFace::DONE;
return;
}
}
......@@ -1437,10 +1481,7 @@ void GFace::replaceEdges(std::list<GEdge*> &new_edges)
void GFace::moveToValidRange(SPoint2 &pt) const
{
// printf("coucou %8d %12.5E %12.5E %d %d\n",
// tag(),pt.x(),pt.y(),
// periodic(0),periodic(1));
for(int i=0; i < 2; i++){
for(int i = 0; i < 2; i++){
if(periodic(i)){
Range<double> range = parBounds(i);
double tol = 1e-6*(range.high()-range.low());
......@@ -1580,7 +1621,8 @@ void GFace::setMeshMaster(GFace* master, const std::vector<double>& tfo)
"for periodic connection of surface %d to %d "
"using the specified transformation"
"Minimum distance is %g with a tolerance of %g",
m_vertex->tag(),master->tag(),tag(),dist_min, CTX::instance()->geom.tolerance * CTX::instance()->lc);
m_vertex->tag(), master->tag(), tag(), dist_min,
CTX::instance()->geom.tolerance * CTX::instance()->lc);
return;
}
gVertexCounterparts[l_vertex] = m_vertex;
......@@ -1605,7 +1647,7 @@ void GFace::setMeshMaster(GFace* master, const std::vector<double>& tfo)
std::pair<GVertex*,GVertex*> mPair(gVertexCounterparts[lPair.first],
gVertexCounterparts[lPair.second]);
int sign = 1;
int sign = 1;
std::map<std::pair<GVertex*,GVertex*>,GEdge*>::iterator mv2eIter = m_vtxToEdge.find(mPair);
if (mv2eIter == m_vtxToEdge.end()) {
sign *= -1;
......@@ -1660,6 +1702,7 @@ struct myPlane {
return n.x() * x + n.y() * y + n.z() * z + a;
}
};
struct myLine {
SPoint3 p;
SVector3 t;
......
......@@ -1593,16 +1593,8 @@ void GModel::_createGeometryOfDiscreteEntities(bool force)
createTopologyFromMeshNew();
exportDiscreteGEOInternals();
}
if(force || CTX::instance()->meshDiscrete){
Msg::Info("Creating the geometry of discrete curves");
for(eiter it = firstEdge(); it != lastEdge(); ++it){
if((*it)->geomType() == GEntity::DiscreteCurve) {
discreteEdge *de = dynamic_cast<discreteEdge*> (*it);
if(de) de->createGeometry();
}
}
}
if(CTX::instance()->meshDiscrete){
// return;
if (CTX::instance()->meshDiscrete){
Msg::Info("Creating the geometry of discrete surfaces");
for(fiter it = firstFace(); it != lastFace(); ++it){
if((*it)->geomType() == GEntity::DiscreteSurface) {
......@@ -1611,6 +1603,15 @@ void GModel::_createGeometryOfDiscreteEntities(bool force)
}
}
}
if (force || CTX::instance()->meshDiscrete){
Msg::Info("Creating the geometry of discrete curves");
for(eiter it = firstEdge(); it != lastEdge(); ++it){
if((*it)->geomType() == GEntity::DiscreteCurve) {
discreteEdge *de = dynamic_cast<discreteEdge*> (*it);
if(de) de->createGeometry();
}
}
}
}
void GModel::_associateEntityWithMeshVertices()
......
......@@ -7,6 +7,8 @@
#include "MEdge.h"
#include "Numeric.h"
// FIXME
// remove that when MElementCut is removed
bool MEdge::isInside(MVertex *v) const
{
double tol = MVertexLessThanLexicographic::getTolerance();
......@@ -62,3 +64,57 @@ bool MEdge::isInside(MVertex *v) const
return true;
}
bool SortEdgeConsecutive(const std::vector<MEdge> &e,
std::vector<std::vector<MVertex*> >&vs)
{
if(e.empty()) return true;
std::map<MVertex*, std::pair<MVertex*, MVertex*> > c;
for (size_t i = 0; i<e.size();i++){
MVertex *v0 = e[i].getVertex(0);
MVertex *v1 = e[i].getVertex(1);
std::map<MVertex*, std::pair<MVertex*,MVertex*> >::iterator it0 = c.find(v0);
std::map<MVertex*, std::pair<MVertex*,MVertex*> >::iterator it1 = c.find(v1);
if (it0 == c.end()) c[v0] = std::make_pair(v1,(MVertex*)NULL);
else it0->second.second = v1;
if (it1 == c.end()) c[v1] = std::make_pair(v0,(MVertex*)NULL);
else it1->second.second = v0;
}
while (!c.empty()){
std::vector<MVertex*> v;
MVertex *start = NULL;
{
std::map<MVertex*, std::pair<MVertex*,MVertex*> >::iterator it = c.begin();
start = it->first;
for (; it != c.end(); ++it) {
if (it->second.second == NULL){
start = it->first;
break;
}
}
}
std::map<MVertex*, std::pair<MVertex*,MVertex*> >::iterator it = c.find(start);
MVertex *prev = (it->second.second == start) ? it->second.first : it->second.second;
MVertex *current = start;
do {
v.push_back(current);
std::map<MVertex*, std::pair<MVertex*,MVertex*> >::iterator it = c.find(current);
c.erase(it);
MVertex *v1 = it->second.first;
MVertex *v2 = it->second.second;
MVertex *temp = current;
if (v1 == prev)current = v2;
else if (v2 == prev)current = v1;
else {
break;
}
prev = temp;
if (current == start) {
v.push_back(current);
}
} while (current != start && current != NULL);
vs.push_back(v);
}
return true;
}
......@@ -39,7 +39,7 @@ class MEdge {
inline MVertex *getSortedVertex(const int i) const { return _v[int(_si[i])]; }
inline MVertex *getMinVertex() const { return _v[int(_si[0])]; }
inline MVertex *getMaxVertex() const { return _v[int(_si[1])]; }
inline int computeCorrespondence(MEdge& other) {
if (other.getVertex(0) == _v[0] && other.getVertex(1) == _v[1]) return 1;
else if (other.getVertex(0) == _v[1] && other.getVertex(1) == _v[0]) return -1;
......@@ -123,4 +123,8 @@ struct Less_Edge : public std::binary_function<MEdge, MEdge, bool> {
}
};
// assume a set of MEdge, give consecutive list of vertices
bool SortEdgeConsecutive(const std::vector<MEdge> &,
std::vector<std::vector<MVertex*> >&vs);
#endif
......@@ -119,6 +119,21 @@ double MTriangle::gammaShapeMeasure()
void MTriangle::xyz2uvw(double xyz[3], double uvw[3]) const
{
/*
double M[2][2], R[2];
const SPoint2 p0 (getVertex(0)->x(),getVertex(0)->y());
const SPoint2 p1 (getVertex(1)->x(),getVertex(1)->y());
const SPoint2 p2 (getVertex(2)->x(),getVertex(2)->y());
M[0][0] = p1.x() - p0.x();
M[0][1] = p2.x() - p0.x();
M[1][0] = p1.y() - p0.y();
M[1][1] = p2.y() - p0.y();
R[0] = (xyz[0] - p0.x());
R[1] = (xyz[1] - p0.y());
sys2x2(M, R, uvw);
return;*/
const double O[3] = {_v[0]->x(), _v[0]->y(), _v[0]->z()};
const double d[3] = {xyz[0] - O[0], xyz[1] - O[1], xyz[2] - O[2]};
const double d1[3] = {_v[1]->x() - O[0], _v[1]->y() - O[1], _v[1]->z() - O[2]};
......
......@@ -10,7 +10,6 @@
#include "GVertex.h"
#include "GEdge.h"
#include "GFace.h"
#include "discreteDiskFace.h"
#include "GmshMessage.h"
#include "StringUtils.h"
......@@ -459,9 +458,8 @@ static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> &params
params.clear();
#if defined(HAVE_ANN) && defined(HAVE_SOLVER)
if (gf->geomType() == GEntity::DiscreteDiskSurface ){
discreteDiskFace *gfc = (discreteDiskFace*) gf;
params.push_back(gfc->parFromVertex(v));
if (gf->geomType() == GEntity::DiscreteSurface ){
params.push_back(gf->parFromPoint(SPoint3(v->x(),v->y(),v->z())));
return ;
}
#endif
......@@ -563,13 +561,12 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
bool reparamMeshVertexOnFace(MVertex *v, const GFace *gf, SPoint2 &param,
bool onSurface)
{
#if defined(HAVE_ANN) && defined(HAVE_SOLVER)
if (gf->geomType() == GEntity::DiscreteDiskSurface ){
discreteDiskFace *gfc = (discreteDiskFace*) gf;
param = gfc->parFromVertex(v);
return true;
}
#endif
//#if defined(HAVE_ANN) && defined(HAVE_SOLVER)
// if (gf->geomType() == GEntity::DiscreteSurface ){
// param=gf->parFromPoint(SPoint3(v->x(),v->y(),v->z()));
// return true;
// }
//#endif
if(v->onWhat()->geomType() == GEntity::DiscreteCurve ||
v->onWhat()->geomType() == GEntity::BoundaryLayerCurve){
......@@ -578,6 +575,11 @@ bool reparamMeshVertexOnFace(MVertex *v, const GFace *gf, SPoint2 &param,
}
if(v->onWhat()->dim() == 0){
if (gf->geomType() == GEntity::DiscreteSurface ){
param=gf->parFromPoint(SPoint3(v->x(),v->y(),v->z()));
return true;
}
GVertex *gv = (GVertex*)v->onWhat();
// hack for bug in periodic curves
if (gv->getNativeType() == GEntity::GmshModel && gf->geomType() == GEntity::Plane)
......@@ -590,6 +592,12 @@ bool reparamMeshVertexOnFace(MVertex *v, const GFace *gf, SPoint2 &param,
if((*it)->isSeam(gf)) return false;
}
else if(v->onWhat()->dim() == 1){
if (gf->geomType() == GEntity::DiscreteSurface ){
param=gf->parFromPoint(SPoint3(v->x(),v->y(),v->z()));
return true;
}
GEdge *ge = (GEdge*)v->onWhat();
double t;
v->getParameter(0, t);
......
......@@ -12,16 +12,14 @@
class discreteEdge : public GEdge {
protected:
std::map<MVertex*,MVertex*> v2v;
std::vector<double> _pars;
std::vector<int> _orientation;// ?
std::map<MVertex*, MLine*> boundv;
bool createdTopo;
std::vector<MVertex*> discrete_vertices;
std::vector<MLine*> discrete_lines;
void orderMLines();
bool getLocalParameter(const double &t, int &iEdge, double &tLoc) const;
public:
discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1);
virtual ~discreteEdge() {}
virtual ~discreteEdge();
virtual GeomType geomType() const { return DiscreteCurve; }
virtual GPoint point(double p) const;
virtual SVector3 firstDer(double par) const;
......@@ -30,25 +28,10 @@ class discreteEdge : public GEdge {
double *curvMax, double *curvMin) const;
virtual bool haveParametrization(){ return !_pars.empty(); }
virtual Range<double> parBounds(int) const;
bool getLocalParameter(const double &t, int &iEdge, double &tLoc) const;
void interpolateInGeometry (MVertex *v, MVertex **v1, MVertex **v2, double &xi) const;
void parametrize(std::map<GFace*, std::map<MVertex*, MVertex*,
std::less<MVertex*> > > &face2Verts,
std::map<GRegion*, std::map<MVertex*, MVertex*,
std::less<MVertex*> > > &region2Vert);
void parametrize(std::map<MVertex*,MVertex*>& old2New);
void orderMLines();
void setBoundVertices();
void setTopo(std::vector<MLine*>);
void createTopo();
void createGeometry();
void computeNormals () const;
virtual void mesh(bool verbose) ;
void writeGEO(FILE *fp);
int minimumDrawSegments() const {return 2*_pars.size();}
MVertex * getGeometricalVertex (MVertex *v);
};
#endif
......@@ -6,6 +6,7 @@
#ifndef _DISCRETE_FACE_H_
#define _DISCRETE_FACE_H_
#include "GmshConfig.h"
#include "GModel.h"
#include "GFace.h"
#include "discreteEdge.h"
......@@ -14,29 +15,63 @@
#include "MEdge.h"
#include "MLine.h"
#if defined(HAVE_HXT)
#include <algorithm>
#include "rtree.h"
class MElementOctree;
extern "C" {
#include "hxt_mesh.h"
#include "hxt_parametrization.h"
#include "hxt_linear_system.h"
}
class hxt_reparam_surf {
public:
MElementOctree *oct;
mutable RTree< std::pair<MTriangle*,MTriangle*> *,double,3> rtree3d;
std::vector<MVertex> v2d;
std::vector<MVertex> v3d;
std::vector<MTriangle> t2d;
std::vector<MTriangle> t3d;
std::vector<discreteEdge*> bnd;
std::vector<discreteEdge*> emb;
hxt_reparam_surf() : oct(NULL) {}
~hxt_reparam_surf();
};
#else
class discreteDiskFace;
class triangulation;
#endif
class discreteFace : public GFace {
private:
bool _meshable;
void checkAndFixOrientation();
#ifdef HAVE_HXT
int _current_parametrization;
std::vector< hxt_reparam_surf > _parametrizations;
std::vector<discreteEdge*> e_internals;
std::vector<discreteVertex*> v_internals;
HXTStatus reparametrize_through_hxt();
bool compute_topology_of_partition(int nbColors,
int *colors,
int *nNodes,
int *nodes,
double *uv,
double angle_threshold,
std::vector<MVertex*> &c2v,
std::vector<std::vector<MEdge> > &boundaries,
std::vector<std::vector<MEdge> > &internals);
#endif
public:
discreteFace(GModel *model, int num, bool meshable=false);
discreteFace(GModel *model, int num);
virtual ~discreteFace() {}
void checkAndFixOrientation();
void setupDiscreteVertex(GVertex*,MVertex*,std::set<MVertex*>*);
void setupDiscreteEdge(discreteEdge*,std::vector<MLine*>,std::set<MVertex*>*);
void splitDiscreteEdge(GEdge*,GVertex*,discreteEdge*[2]);
void updateTopology(std::vector<triangulation*>&);
void split(triangulation*,std::vector<triangulation*>&,int);
void fillHoles(triangulation*);
virtual void addTriangle(triangulation*,MTriangle*);
using GFace::addTriangle;
void complex_crossField();
void crossField();
using GFace::point;
GPoint point(double par1, double par2) const;
SPoint2 parFromPoint(const SPoint3 &p, bool onSurface=true) const;
GPoint closestPoint(const SPoint3 &queryPoint, double maxDistance,
SVector3 *normal = NULL) const;
GPoint closestPoint(const SPoint3 &queryPoint, const double initialGuess[2]) const;
SVector3 normal(const SPoint2 &param) const;
double curvatureMax(const SPoint2 &param) const;
double curvatures(const SPoint2 &param, SVector3 *dirMax, SVector3 *dirMin,
......@@ -45,17 +80,16 @@ class discreteFace : public GFace {
virtual Pair<SVector3, SVector3> firstDer(const SPoint2 &param) const;
virtual void secondDer(const SPoint2 &param,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const;
void setBoundEdges(const std::vector<int> &tagEdges);
void setBoundEdges(const std::vector<int> &tagEdges,
const std::vector<int> &signEdges);
void findEdges(std::map<MEdge, std::vector<int>, Less_Edge > &map_edges);
void writeGEO(FILE *fp);
void createGeometry();
void gatherMeshes();
virtual void mesh(bool verbose);
std::vector<discreteDiskFace*> _atlas;
std::vector<GFace*> _CAD;
std::map<MEdge,std::vector<int>,Less_Edge> allEdg2Tri;
void setBoundEdges(const std::vector<int> &tagEdges);
void setBoundEdges(const std::vector<int> &tagEdges,
const std::vector<int> &signEdges);
int trianglePosition(double par1, double par2, double &u, double &v) const;
GPoint intersectionWithCircle(const SVector3 &n1, const SVector3 &n2,
const SVector3 &p, const double &R,
double uv[2]) ;
};
#endif
......@@ -12,6 +12,7 @@
#include "Numeric.h"
#include "BDS.h"
#include "GFace.h"
#include "discreteFace.h"
#include "meshGFaceDelaunayInsertion.h"
#include "qualityMeasures.h"
......@@ -634,12 +635,29 @@ bool BDS_Mesh::split_edge(BDS_Edge *e, BDS_Point *mid)
// p1,op2,mid +
*/
BDS_Point *op[2];
BDS_Point *p1 = e->p1;
BDS_Point *p2 = e->p2;
e->oppositeof(op);
// double _p1 [2] = {p1->u,p1->v};
// double _p2 [2] = {p2->u,p2->v};
// double _op1[2] = {op[0]->u,op[0]->v};
// double _op2[2] = {op[1]->u,op[1]->v};
// double _mid[2] = {mid->u,mid->v};
// double ori1 = robustPredicates::orient2d(_mid, _p1, _op2);
// double ori2 = robustPredicates::orient2d(_mid, _op2, _p2);
// double ori3 = robustPredicates::orient2d(_mid, _p2, _op1);