diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 078e119e9d994816f34521816862c7856ce2a6e9..bc965d60e403cd1ee54d4126d128bd0c7af47895 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -72,18 +72,24 @@ static void printBound(std::vector<MVertex*> &l, int tag)
}
*/
-static std::vector<MVertex*> getBlob(unsigned int minNbPt, v2t_cont::iterator it, v2t_cont &adj)
+static std::vector<MVertex*> getBlob(unsigned int minNbPt, v2t_cont::iterator it,
+ v2t_cont &adj)
{
-
- std::vector<MVertex*> vv(1,it->first), bvv = vv; // Vector of vertices in blob and in boundary of blob
+ std::vector<MVertex*> vv(1,it->first), bvv = vv;
+ // Vector of vertices in blob and in boundary of blob
do {
- std::set<MVertex*> nbvv; // Set of vertices in new boundary
- for (std::vector<MVertex*>::iterator itBV = bvv.begin(); itBV != bvv.end(); itBV++) { // For each boundary vertex...
+ std::set<MVertex*> nbvv;
+ // Set of vertices in new boundary
+ for (std::vector<MVertex*>::iterator itBV = bvv.begin(); itBV != bvv.end(); itBV++) {
+ // For each boundary vertex...
std::vector<MElement*> &adjBV = adj[*itBV];
- for (std::vector<MElement*>::iterator itBVEl = adjBV.begin(); itBVEl != adjBV.end(); itBVEl++) {
- for (int iV=0; iV<(*itBVEl)->getNumVertices(); iV++){ // ... look for adjacent vertices...
+ for (std::vector<MElement*>::iterator itBVEl = adjBV.begin(); itBVEl != adjBV.end();
+ itBVEl++) {
+ for (int iV=0; iV<(*itBVEl)->getNumVertices(); iV++){
+ // ... look for adjacent vertices...
MVertex *v = (*itBVEl)->getVertex(iV);
- if (find(vv.begin(),vv.end(),v) == vv.end()) nbvv.insert(v); // ... and add them in the new boundary if they are not already in the blob
+ if (find(vv.begin(),vv.end(),v) == vv.end()) nbvv.insert(v);
+ // ... and add them in the new boundary if they are not already in the blob
}
}
}
@@ -93,10 +99,10 @@ static std::vector<MVertex*> getBlob(unsigned int minNbPt, v2t_cont::iterator it
vv.insert(vv.end(),nbvv.begin(),nbvv.end());
}
// } while (!bvv.empty());
- } while (vv.size() < minNbPt); // Repeat until min. number of points is reached
+ } while (vv.size() < minNbPt);
+ // Repeat until min. number of points is reached
return vv;
-
}
static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
@@ -984,8 +990,8 @@ bool GFaceCompound::parametrize() const
_type = UNITCIRCLE;
coordinates.clear();
Octree_Delete(oct);
- parametrize(ITERU,CONVEX);
- parametrize(ITERV,CONVEX);
+ parametrize(ITERU, CONVEX);
+ parametrize(ITERV, CONVEX);
checkOrientation(0);
buildOct();
}
@@ -1177,10 +1183,9 @@ void GFaceCompound::computeALoop(std::set<GEdge*> &_unique, std::list<GEdge*> &l
GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, typeOfCompound toc,
- int allowPartition,
- linearSystem<double>* lsys)
+ int allowPartition)
: GFace(m, tag), _compound(compound), _U0(U0), oct(0), octNew(0),
- _lsys(lsys), _toc(toc), _allowPartition(allowPartition)
+ _toc(toc), _allowPartition(allowPartition)
{
ONE = new simpleFunction<double>(1.0);
MONE = new simpleFunction<double>(-1.0);
@@ -1223,10 +1228,9 @@ 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,
typeOfCompound toc,
- int allowPartition,
- linearSystem<double>* lsys)
+ int allowPartition)
: GFace(m, tag), _compound(compound), _U0(U0), _V0(V0), _U1(U1), _V1(V1),
- oct(0), octNew(0), _lsys(lsys), _toc(toc), _allowPartition(allowPartition)
+ oct(0), octNew(0), _toc(toc), _allowPartition(allowPartition)
{
ONE = new simpleFunction<double>(1.0);
MONE = new simpleFunction<double>(-1.0);
@@ -1270,29 +1274,70 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
GFaceCompound::~GFaceCompound()
{
- for (unsigned int i = 0; i < myParamVert.size(); i++) delete myParamVert[i];
- for (unsigned int i = 0; i < myParamElems.size(); i++) delete myParamElems[i];
+ delete ONE;
+ delete MONE;
+ _deleteInternals();
+}
+
+void GFaceCompound::_deleteInternals()
+{
+ for (unsigned int i = 0; i < myParamVert.size(); i++)
+ delete myParamVert[i];
+ myParamVert.clear();
+ for (unsigned int i = 0; i < myParamElems.size(); i++)
+ delete myParamElems[i];
+ myParamElems.clear();
+ _3Dto2D.clear();
+ _2Dto3D.clear();
+ XYZoct.clear();
+ allNodes.clear();
+ adjv.clear();
+ coordinates.clear();
+ firstDerivatives.clear();
+ xuu.clear();
+ xvv.clear();
+ xuv.clear();
+ xu.clear();
+ xv.clear();
+ _coordPoints.clear();
+ _normals.clear();
+ fillTris.clear();
+ fillNodes.clear();
+ fillFaces.clear();
+ _ordered.clear();
+ _coords.clear();
+ _mapV.clear();
if(oct){
Octree_Delete(oct);
delete [] _gfct;
+ oct = 0;
+ }
+ if(octNew){
+ delete octNew;
+ octNew = 0;
}
- if(octNew) delete(octNew);
- if (_lsys)delete _lsys;
- delete ONE;
- delete MONE;
if(uv_kdtree) {
ANNpointArray uv_nodes = uv_kdtree->thePoints();
if(uv_nodes) annDeallocPts(uv_nodes);
delete uv_kdtree;
+ uv_kdtree = 0;
}
if(kdtree){
ANNpointArray nodes = kdtree->thePoints();
if(nodes) annDeallocPts(nodes);
delete kdtree;
+ kdtree = 0;
}
}
+
+void GFaceCompound::deleteMesh()
+{
+ _deleteInternals();
+ GFace::deleteMesh();
+}
+
SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
{
if(trivial()){
@@ -1389,33 +1434,31 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
{
- linearSystem<double> *lsys = 0;
- if (_lsys) lsys = _lsys;
- else{
- if (tom==HARMONIC){
+ linearSystem<double> *lsys;
+
+ if (tom==HARMONIC){
#if defined(HAVE_TAUCS)
- lsys = new linearSystemCSRTaucs<double>;
-#elif defined(HAVE_PETSC) && !defined(HAVE_TAUCS)
- lsys = new linearSystemPETSc<double>;
-#elif defined(HAVE_GMM) && !defined(HAVE_TAUCS)
- linearSystemGmm<double> *lsysb = new linearSystemGmm<double>;
- lsysb->setGmres(1);
- lsys = lsysb;
+ lsys = new linearSystemCSRTaucs<double>;
+#elif defined(HAVE_PETSC)
+ lsys = new linearSystemPETSc<double>;
+#elif defined(HAVE_GMM)
+ linearSystemGmm<double> *lsysb = new linearSystemGmm<double>;
+ lsysb->setGmres(1);
+ lsys = lsysb;
#else
- lsys = new linearSystemFull<double>;
+ lsys = new linearSystemFull<double>;
#endif
- }
- else if (tom==CONVEX){
+ }
+ else{
#if defined(HAVE_PETSC)
- lsys = new linearSystemPETSc<double>;
+ lsys = new linearSystemPETSc<double>;
#elif defined(HAVE_GMM)
- linearSystemGmm<double> *lsysb = new linearSystemGmm<double>;
- lsysb->setGmres(1);
- lsys = lsysb;
+ linearSystemGmm<double> *lsysb = new linearSystemGmm<double>;
+ lsysb->setGmres(1);
+ lsys = lsysb;
#else
- lsys = new linearSystemFull<double>;
+ lsys = new linearSystemFull<double>;
#endif
- }
}
dofManager<double> myAssembler(lsys);
@@ -1548,10 +1591,7 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
}
}
- if (_lsys)
- lsys->clear();
- else
- delete lsys;
+ delete lsys;
}
bool GFaceCompound::parametrize_conformal_spectral() const
@@ -1561,8 +1601,8 @@ bool GFaceCompound::parametrize_conformal_spectral() const
return parametrize_conformal(0, NULL, NULL);
#else
- linearSystem <double> *lsysA = new linearSystemPETSc<double>;
- linearSystem <double> *lsysB = new linearSystemPETSc<double>;
+ linearSystem <double> *lsysA = new linearSystemPETSc<double>;
+ linearSystem <double> *lsysB = new linearSystemPETSc<double>;
dofManager<double> myAssembler(lsysA, lsysB);
myAssembler.setCurrentMatrix("A");
@@ -1675,15 +1715,15 @@ bool GFaceCompound::parametrize_conformal_spectral() const
bool GFaceCompound::parametrize_conformal(int iter, MVertex *v1, MVertex *v2) const
{
- linearSystem<double> *lsys = 0;
-#if defined(HAVE_PETSC) && !defined(HAVE_TAUCS)
+ linearSystem<double> *lsys;
+#if defined(HAVE_TAUCS)
+ lsys = new linearSystemCSRTaucs<double>;
+#elif defined(HAVE_PETSC)
lsys = new linearSystemPETSc<double>;
-#elif defined(HAVE_GMM) && !defined(HAVE_TAUCS)
+#elif defined(HAVE_GMM)
linearSystemGmm<double> *lsysb = new linearSystemGmm<double>;
lsysb->setGmres(1);
lsys = lsysb;
-#elif defined(HAVE_TAUCS)
- lsys = new linearSystemCSRTaucs<double>;
#else
lsys = new linearSystemFull<double>;
#endif
@@ -2185,48 +2225,8 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
SVector3 dXdu = dXdu1*(1.-U-V) + dXdu2*U + dXdu3*V;
SVector3 dXdv = dXdv1*(1.-U-V) + dXdv2*U + dXdv3*V;
return Pair<SVector3, SVector3>(dXdu,dXdv);
-
}
-// Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
-// {
-// if(!oct) parametrize();
-
-// if(trivial())
-// return (*(_compound.begin()))->firstDer(param);
-
-// double U, V;
-// GFaceCompoundTriangle *lt;
-// getTriangle(param.x(), param.y(), <, U,V);
-// if(!lt && _mapping != RBF)
-// return Pair<SVector3, SVector3>(SVector3(1, 0, 0), SVector3(0, 1, 0));
-// else if (!lt && _mapping == RBF){
-// double x, y, z;
-// SVector3 dXdu, dXdv ;
-// _rbf->UVStoXYZ(param.x(), param.y(), x, y, z, dXdu, dXdv);
-// return Pair<SVector3, SVector3>(dXdu, dXdv);
-// }
-
-// 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];
-// double det = inv2x2(mat,inv);
-// if (!det && _mapping == RBF){
-// double x, y, z;
-// SVector3 dXdu, dXdv ;
-// _rbf->UVStoXYZ(param.x(), param.y(), x, y, z, dXdu, dXdv);
-// return Pair<SVector3, SVector3>(dXdu, dXdv);
-// }
-
-// SVector3 dXdxi(lt->v2 - lt->v1);
-// SVector3 dXdeta(lt->v3 - lt->v1);
-
-// SVector3 dXdu(dXdxi * inv[0][0] + dXdeta * inv[1][0]);
-// SVector3 dXdv(dXdxi * inv[0][1] + dXdeta * inv[1][1]);
-
-// return Pair<SVector3, SVector3>(dXdu, dXdv);
-// }
-
void GFaceCompound::secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
{
@@ -2270,7 +2270,8 @@ void GFaceCompound::computeHessianMapping() const
for (v2t_cont::iterator it = adjv.begin(); it != adjv.end(); it++) {
MVertex *ver = it->first;
std::vector<MVertex*> vv = getBlob(minNbPtBlob, it, adjv);
- fullMatrix<double> A(vv.size(),sysDim), ATAx(sysDim,sysDim), ATAy(sysDim,sysDim), ATAz(sysDim,sysDim);
+ fullMatrix<double> A(vv.size(),sysDim), ATAx(sysDim,sysDim);
+ fullMatrix<double> ATAy(sysDim,sysDim), ATAz(sysDim,sysDim);
fullVector<double> bx(vv.size()), ATbx(sysDim), coeffsx(sysDim);
fullVector<double> by(vv.size()), ATby(sysDim), coeffsy(sysDim);
fullVector<double> bz(vv.size()), ATbz(sysDim), coeffsz(sysDim);
@@ -2304,6 +2305,7 @@ void GFaceCompound::computeHessianMapping() const
#endif
}
+
static void GFaceCompoundBB(void *a, double*mmin, double*mmax)
{
GFaceCompoundTriangle *t = (GFaceCompoundTriangle *)a;
@@ -2381,7 +2383,6 @@ void GFaceCompound::getTriangle(double u, double v,
void GFaceCompound::buildOct() const
{
-
#if defined(HAVE_MESH)
SBoundingBox3d bb;
int count = 0;
@@ -2556,7 +2557,8 @@ bool GFaceCompound::checkTopology() const
}
else if (G == 0 && AR > AR_MAX){
correctTopo = false;
- Msg::Warning("Wrong topology: Aspect ratio is too high AR=%d (AR1=%d AR2=%d)", AR, AR1, AR2);
+ Msg::Warning("Wrong topology: Aspect ratio is too high AR=%d (AR1=%d AR2=%d)",
+ AR, AR1, AR2);
if (_allowPartition == 1){
nbSplit = -2;
Msg::Info("-----------------------------------------------------------");
@@ -2586,7 +2588,6 @@ bool GFaceCompound::checkTopology() const
double GFaceCompound::checkAspectRatio() const
{
-
if(allNodes.empty()) buildAllNodes();
double limit = 1.e-20;
@@ -2663,9 +2664,9 @@ void GFaceCompound::coherencePatches() const
mySet.insert(t);
it->second = mySet;
}
+ }
}
}
-}
std::set<MElement* , std::less<MElement*> > touched;
int iE, si, iE2, si2;
@@ -2863,7 +2864,6 @@ void GFaceCompound::printStuff(int iNewton) const
fclose(xyzu);
fprintf(xyzv,"};\n");
fclose(xyzv);
-
}
// useful for mesh generators ----------------------------------------
@@ -2895,7 +2895,8 @@ GPoint GFaceCompound::intersectionWithCircle(const SVector3 &n1, const SVector3
// Compute w = t1 x t2 = (a,b,c) and write a point of the plabe as
// X(x,y,z) = ax + by + cz - (v1_x a + v1_y b + v1_z * c) = 0
// Then
- // a (p_x + t n1_x) + a (p_y + t n1_y) + c (p_z + t n1_z) - (v1_x a + v1_y b + v1_z * c) = 0
+ // a (p_x + t n1_x) + a (p_y + t n1_y) + c (p_z + t n1_z) -
+ // (v1_x a + v1_y b + v1_z * c) = 0
// gives t
SVector3 w = crossprod(t1,t2);
@@ -2909,7 +2910,8 @@ GPoint GFaceCompound::intersectionWithCircle(const SVector3 &n1, const SVector3
// (x-p_x)^2 + (y-p_y)^2 + (z-p_z)^2 = d^2
// Inserting the parametric form of the line into the sphere gives
// (t m_x + q_x - p_x)^2 + (t m_y + q_y - p_y)^2 + (t m_z + q_z - p_z)^2 = d^2
- // t^2 (m_x^2 + m_y^2 + m_z^2) + 2 t (m_x (q_x - p_x) + m_y (q_y - p_z) + m_z (q_z - p_z)) + ((q_x - p_x)^2 + (q_y - p_y)^2 + (q_z - p_z)^2 - d^2) = 0
+ // t^2 (m_x^2 + m_y^2 + m_z^2) + 2 t (m_x (q_x - p_x) + m_y (q_y - p_z) +
+ // m_z (q_z - p_z)) + ((q_x - p_x)^2 + (q_y - p_y)^2 + (q_z - p_z)^2 - d^2) = 0
// t^2 m^2 + 2 t (m . (q-p)) + ((q-p)^2 - d^2) = 0
// Let us call ta and tb the two roots
// they correspond to two points s_1 = q + ta m and s_2 = q + tb m
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 0022b889823c9ddd5d59729f343b72e7de51aa5d..6e1b0e22f3c85ef94a5b72e43d571d2a10d10a66 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -80,7 +80,6 @@ class GFaceCompound : public GFace {
mutable std::map<int,SPoint3> XYZoct;
mutable std::set<MVertex*> allNodes;
mutable v2t_cont adjv;
- mutable bool mapv2Tri;
mutable std::map<MVertex*, SPoint3> coordinates;
mutable std::map<MVertex*, Pair<SVector3,SVector3> > firstDerivatives;
mutable std::map<MVertex*, SVector3> xuu;
@@ -96,7 +95,6 @@ class GFaceCompound : public GFace {
mutable std::vector<MVertex*> _ordered;
mutable std::vector<double> _coords;
mutable std::map<MVertex*, int> _mapV;
- linearSystem <double> *_lsys;
mutable ANNkd_tree *uv_kdtree;
mutable ANNkd_tree *kdtree;
void buildOct() const ;
@@ -136,15 +134,14 @@ class GFaceCompound : public GFace {
public:
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, typeOfCompound typ = HARMONIC_CIRCLE,
- int allowPartition=1,
- linearSystem<double>* lsys =0);
+ int allowPartition=1);
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,
typeOfCompound typ = HARMONIC_CIRCLE,
- int allowPartition=1,
- linearSystem<double>* lsys =0);
+ int allowPartition=1);
~GFaceCompound();
+ virtual void deleteMesh();
Range<double> parBounds(int i) const
{ return trivial() ? (*(_compound.begin()))->parBounds(i) : Range<double>(-1, 1); }
@@ -182,6 +179,7 @@ class GFaceCompound : public GFace {
const double &d, double uv[2]) const;
private:
+ void _deleteInternals();
mutable typeOfCompound _toc;
mutable typeOfMapping _mapping;
mutable typeOfIsomorphism _type;
@@ -200,8 +198,7 @@ class GFaceCompound : public GFace {
typedef enum {UNITCIRCLE, MEANPLANE, SQUARE, ALREADYFIXED,SPECTRAL, FE} typeOfIsomorphism;
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, typeOfCompound typ = HARMONIC_CIRCLE,
- int allowPartition=1,
- linearSystem<double>* lsys =0)
+ int allowPartition=1)
: GFace(m, tag)
{
Msg::Error("Gmsh has to be compiled with Solver and ANN support to use GFaceCompounds");
@@ -210,13 +207,13 @@ class GFaceCompound : public GFace {
std::list<GEdge*> &U0, std::list<GEdge*> &V0,
std::list<GEdge*> &U1, std::list<GEdge*> &V1,
typeOfCompound typ = HARMONIC_CIRCLE,
- int allowPartition=1,
- linearSystem<double>* lsys =0)
+ int allowPartition=1)
: GFace(m, tag)
{
Msg::Error("Gmsh has to be compiled with Solver and ANN support to use GFaceCompounds");
}
virtual ~GFaceCompound() {}
+ virtual void deleteMesh() {}
GPoint point(double par1, double par2) const { return GPoint(); }
Pair<SVector3, SVector3> firstDer(const SPoint2 ¶m) const
{