diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp
index d0098daf3f31f08ac23fea9d8c6322331a3bfda6..b84169170aa8d4085806eacaa59de64fb423ef18 100644
--- a/Common/OpenFile.cpp
+++ b/Common/OpenFile.cpp
@@ -340,6 +340,9 @@ int MergeFile(std::string fileName, bool warnIfMissing)
else if(ext == ".ply2"){
status = GModel::current()->readPLY2(fileName);
}
+ else if(ext == ".ply"){
+ status = GModel::current()->readPLY(fileName);
+ }
else {
CTX::instance()->geom.draw = 1;
if(!strncmp(header, "$PTS", 4) || !strncmp(header, "$NO", 3) ||
diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index 1f00a940b4e83e83c1c54ea3f3aa76dda209d9d0..546439453fe39fb445b54e0462a8b4a431a86ce2 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -239,7 +239,8 @@ void statisticsWindow::compute(bool elementQuality)
dMax = std::max(nbE, dMax);
if (nbE == 4) d4 += 1;
}
- printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n", dMin, dMax, (double)d4/nbElems);
+ if (nbElems > 0)
+ printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n", dMin, dMax, (double)d4/nbElems);
//end emi hack
int num = 0;
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index ee14cf5c7f9feee62a62d2bd7571690a4c2d34e7..002f5232ff7b8eeb5beb10688b124be48c845420 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -664,19 +664,18 @@ bool GFaceCompound::parametrize() const
fillNeumannBCS();
bool noOverlap = parametrize_conformal_spectral() ;
if (!noOverlap){
- //buildOct(); exit(1);
Msg::Warning("!!! Overlap: parametrization switched to 'FE conformal' map");
noOverlap = parametrize_conformal();
}
// if (!noOverlap) {
- // Msg::Warning("!!! Overlap: parametrization switched to 'nonLIN conformal' map");
- // noOverlap = parametrize_conformal_nonLinear() ;
- // }
- if ( !noOverlap || !checkOrientation(0) ){
- Msg::Warning("$$$ Overlap: parametrization switched to 'harmonic' map");
- parametrize(ITERU,HARMONIC);
- parametrize(ITERV,HARMONIC);
- }
+ // Msg::Warning("!!! Overlap: parametrization switched to 'nonLIN conformal' map");
+ // noOverlap = parametrize_conformal_nonLinear() ;
+ //}
+ if ( !noOverlap || !checkOrientation(0) ){
+ Msg::Warning("$$$ Overlap: parametrization switched to 'harmonic' map");
+ parametrize(ITERU,HARMONIC);
+ parametrize(ITERV,HARMONIC);
+ }
}
buildOct();
@@ -1079,9 +1078,22 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
{
dofManager<double> myAssembler(_lsys);
+
+ // EMI-test for paper Dong: fix only 2 vertices
+ // MVertex *v1;
+ // MVertex *v2;
+ // double zmin = 1.e6;
+ // double zmax = 0.0;
+ // for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+ // MVertex *v = *itv;
+ // double z = v->z();
+ // if (z > zmax){ zmax = z; v2 = v;}
+ // if (z < zmin){ zmin = z; v1 = v;}
+ // }
+ // myAssembler.fixVertex(v1, 0, 1, 0.0);
+ // myAssembler.fixVertex(v2, 0, 1, 1.0);
if(_type == UNITCIRCLE){
- //map to a unit circle
for(unsigned int i = 0; i < _ordered.size(); i++){
MVertex *v = _ordered[i];
const double theta = 2 * M_PI * _coords[i];
@@ -1954,58 +1966,6 @@ void GFaceCompound::getTriangle(double u, double v,
_v = X[1];
}
-void GFaceCompound::partitionFaceCM()
-{
-
- if(!oct) parametrize();
-
- double CMu = 0.0;
- double sumArea = 0.0;
-
- 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];
- std::map<MVertex*,SPoint3>::const_iterator it0 = coordinates.find(t->getVertex(0));
- std::map<MVertex*,SPoint3>::const_iterator it1 = coordinates.find(t->getVertex(1));
- std::map<MVertex*,SPoint3>::const_iterator it2 = coordinates.find(t->getVertex(2));
- double q0[3] = {it0->second.x(), it0->second.y(), 0.0};
- double q1[3] = {it1->second.x(), it1->second.y(), 0.0};
- double q2[3] = {it2->second.x(), it2->second.y(), 0.0};
- double area = 1/fabs(triangle_area(q0, q1, q2));
- double cg_u = (q0[0]+q1[0]+q2[0])/3.;
- CMu += cg_u*area;
- sumArea += area;
- }
- CMu /= sumArea;
- }
-
- //printf("min size partition =%d \n", (int)allNodes.size()/2);
- model()->setMinPartitionSize((int)allNodes.size()/2);
- model()->setMaxPartitionSize((int)allNodes.size()/2+1);
-
- it = _compound.begin();
- for( ; it != _compound.end() ; ++it){
- for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
- MTriangle *t = (*it)->triangles[i];
- std::map<MVertex*,SPoint3>::const_iterator it0 = coordinates.find(t->getVertex(0));
- std::map<MVertex*,SPoint3>::const_iterator it1 = coordinates.find(t->getVertex(1));
- std::map<MVertex*,SPoint3>::const_iterator it2 = coordinates.find(t->getVertex(2));
- double cg_u = (it0->second.x()+it1->second.x()+it2->second.x())/3.;
- if (cg_u <= CMu)
- t->setPartition(1);
- else
- t->setPartition(2);
- }
- }
-
- model()->recomputeMeshPartitions();
-
- CreateOutputFile("toto.msh", CTX::instance()->mesh.fileFormat);
- Msg::Exit(1);
-
- return;
-}
void GFaceCompound::buildOct() const
{
@@ -2079,7 +2039,6 @@ bool GFaceCompound::checkTopology() const
// FIXME!!! I think those things are wrong with cross-patch reparametrization
//if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)return true;
-
//TODO: smthg to exit here for lloyd remeshing
bool correctTopo = true;
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 47304e85e72b5c53240a2f7b12563798d7dff90b..2613de7a4d31e39768a715208d40585a3218a979 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -131,7 +131,6 @@ class GFaceCompound : public GFace {
bool parametrize() const ;
void coherenceNormals();
void coherencePatches() const;
- void partitionFaceCM();
virtual std::list<GFace*> getCompounds() const { return _compound; }
mutable int nbSplit;
int getNbSplit() const { return nbSplit; }
diff --git a/Geo/GModel.h b/Geo/GModel.h
index f4445d0def335bfb0ced1e4d2ad14fb0b72530c9..8a1297ea358f7e0bda994fcde11720086cfb54d0 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -491,7 +491,8 @@ class GModel
int writeSTL(const std::string &name, bool binary=false,
bool saveAll=false, double scalingFactor=1.0);
- // PLY2 format (asciii text format)
+ // PLY(2) format (ascii text format)
+ int readPLY(const std::string &name);
int readPLY2(const std::string &name);
int writePLY2(const std::string &name);
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 77b6043e99d8033dc5358eb7a4dc990919cb9d82..718d3f1374895adb3cb97ce90668fdcbecd19d45 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -29,6 +29,9 @@
#include "discreteFace.h"
#include "discreteRegion.h"
#include "MVertexPositionSet.h"
+#include "PView.h"
+#include "PViewData.h"
+#include "PViewDataList.h"
void GModel::_storePhysicalTagsInEntities(int dim,
std::map<int, std::map<int, std::string> > &map)
@@ -53,6 +56,14 @@ void GModel::_storePhysicalTagsInEntities(int dim,
}
}
+ static void replaceCommaByDot(const std::string name){
+ char myCommand[1000], myCommand2[1000];
+ sprintf(myCommand, "sed 's/,/./g' %s > temp.txt", name.c_str());
+ system(myCommand);
+ sprintf(myCommand2, "mv temp.txt %s ", name.c_str());
+ system(myCommand2);
+ }
+
static bool getVertices(int num, int *indices, std::map<int, MVertex*> &map,
std::vector<MVertex*> &vertices)
{
@@ -80,6 +91,19 @@ static bool getVertices(int num, int *indices, std::vector<MVertex*> &vec,
}
return true;
}
+static bool getProperties(int num, int *indices, std::vector<double> &vec,
+ std::vector<double> &properties)
+{
+ for(int i = 0; i < num; i++){
+ if(indices[i] < 0 || indices[i] > (int)(vec.size() - 1)){
+ Msg::Error("Wrong vertex index %d", indices[i]);
+ return false;
+ }
+ else
+ properties.push_back(vec[indices[i]]);
+ }
+ return true;
+}
static MElement *createElementMSH(GModel *m, int num, int typeMSH, int physical,
int reg, int part, std::vector<MVertex*> &v,
@@ -1226,6 +1250,131 @@ static int readElementsVRML(FILE *fp, std::vector<MVertex*> &vertexVector, int r
return 1;
}
+int GModel::readPLY(const std::string &name)
+{
+
+ replaceCommaByDot(name);
+
+ FILE *fp = fopen(name.c_str(), "r");
+ if(!fp){
+ Msg::Error("Unable to open file '%s'", name.c_str());
+ return 0;
+ }
+
+
+ std::vector<MVertex*> vertexVector;
+ std::map<int, std::vector<MElement*> > elements[5];
+ std::map<int, std::vector<double> > properties;
+
+ char buffer[256], str[256], str2[256], str3[256];
+ std::string s1;
+ int elementary = getMaxElementaryNumber(-1) + 1;
+ int nbv, nbf;
+ int nbprop = 0;
+ int nbView = 0;
+ std::vector<std::string> propName;
+ while(!feof(fp)) {
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ if(buffer[0] != '#'){ // skip comments
+ sscanf(buffer, "%s %s", str, str2);
+ if(!strcmp(str, "element") && !strcmp(str2, "vertex")){
+ sscanf(buffer, "%s %s %d", str, str2, &nbv);
+ }
+ if(!strcmp(str, "format") && strcmp(str2, "ascii")){
+ Msg::Error("Only reading of ascii PLY files implemented");
+ return 0;
+ }
+ if(!strcmp(str, "property") && strcmp(str2, "list")){
+ nbprop++;
+ sscanf(buffer, "%s %s %s", str, str2, str3);
+ if (nbprop > 3) propName.push_back(s1+str3);
+ }
+ if(!strcmp(str, "element") && !strcmp(str2, "face")){
+ sscanf(buffer, "%s %s %d", str, str2, &nbf);
+ }
+ if(!strcmp(str, "end_header")){
+ nbView = nbprop -3;
+ Msg::Info("%d elements", nbv);
+ Msg::Info("%d triangles", nbf);
+ Msg::Info("%d properties", nbView);
+
+ //printf("*********READING VERTEX \n");
+ vertexVector.resize(nbv);
+ for(int i = 0; i < nbv; i++) {
+ double x,y,z;
+ char line[10000], *pEnd, *pEnd2, *pEnd3;
+ fgets(line, sizeof(line), fp);
+ x = strtod(line, &pEnd);
+ y = strtod(pEnd, &pEnd2);
+ z = strtod(pEnd2, &pEnd3);
+ //printf("xyz=%g %g %g \n", x,y,z);
+ vertexVector[i] = new MVertex(x, y, z);
+
+ pEnd = pEnd3;
+ double prop[nbView];
+ for (int k=0; k< nbView; k++){
+ prop[k]=strtod(pEnd, &pEnd2);
+ pEnd = pEnd2;
+ properties[k].push_back(prop[k]);
+ }
+ }
+
+ //printf("*********READING ELEMS \n");
+ for(int i = 0; i < nbf; i++) {
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int n[3], nbe;
+ int nb = sscanf(buffer, "%d %d %d %d", &nbe, &n[0], &n[1], &n[2]);
+ std::vector<MVertex*> vertices;
+ if(!getVertices(3, n, vertexVector, vertices)) return 0;
+ elements[0][elementary].push_back(new MTriangle(vertices));
+ }
+
+ }
+
+ }
+ }
+
+ for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
+ _storeElementsInEntities(elements[i]);
+ _associateEntityWithMeshVertices();
+ _storeVerticesInEntities(vertexVector);
+
+
+ //Create PViews here
+ std::vector<GEntity*> _entities;
+ getEntities(_entities);
+ for (int iV=0; iV< nbView; iV++){
+ PView *view = new PView();
+ PViewDataList *data = dynamic_cast<PViewDataList*>(view->getData());
+ for(unsigned int ii = 0; ii < _entities.size(); ii++){
+ for(unsigned int i = 0; i < _entities[ii]->getNumMeshElements(); i++){
+ MElement *e = _entities[ii]->getMeshElement(i);
+ int numNodes = e->getNumVertices();
+ std::vector<double> x(numNodes), y(numNodes), z(numNodes);
+ std::vector<double> *out = data->incrementList(1, e->getType());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->x());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->y());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->z());
+ std::vector<double> props;
+ int n[3];
+ n[0] = e->getVertex(0)->getNum()-1;
+ n[1] = e->getVertex(1)->getNum()-1;
+ n[2] = e->getVertex(2)->getNum()-1;
+ if(!getProperties(3, n, properties[iV], props)) return 0;
+ for(int nod = 0; nod < numNodes; nod++) out->push_back(props[nod]);
+ }
+ }
+ data->setName(propName[iV]);
+ data->Time.push_back(0);
+ data->setFileName("property.pos");
+ data->finalize();
+ }
+
+ fclose(fp);
+
+ return 1;
+}
+
int GModel::readPLY2(const std::string &name)
{
FILE *fp = fopen(name.c_str(), "r");
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index b336104087a0a7973eea00a9359889de653e8c98..8e88e8cf462103c4f5c7c7b014b77afb892a66d5 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -460,7 +460,7 @@ int PartitionGraph(Graph &graph, meshPartitionOptions &options)
}
break;
case 3: // Nodal weight
- printf("METIS with weights\n");
+ //printf("METIS with weights\n");
metisOptions[0] = 1;
metisOptions[1] = options.edge_matching;
metisOptions[2] = 1;
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index 96751985b04a1b3d7c44a026d25d568571a16c51..fed14108c18b73cca93c9848dc793951930266cc 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -5,6 +5,7 @@
#include "MEdge.h"
#include "MElement.h"
#include "multiscaleLaplace.h"
+#include "GFaceCompound.h"
#include "Numeric.h"
#include "Context.h"
@@ -299,6 +300,11 @@ multiscalePartition::multiscalePartition(std::vector<MElement *> &elements,
options.global_method = 1;// 1 Multilevel-KL, 2 Spectral
options.mesh_dims[0] = nbParts;
}
+ else if (options.partitioner == 2){
+ options.algorithm = 2;//1 recursive, 2=kway, 3=nodal weights
+ options.refine_algorithm=2;
+ options.edge_matching = 3;
+ }
partitionLevel *level = new partitionLevel;
level->elements.insert(level->elements.begin(),elements.begin(),elements.end());
@@ -364,7 +370,7 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts,
nextLevel->recur,nextLevel->region, genus, AR, nbParts);
partition(*nextLevel, nbParts, MULTILEVEL);
}
- else if (genus == 0 && AR > 5 ){//|| genus == 0 && NB > 1){
+ else if (genus == 0 && AR > AR_MAX || genus == 0 && NB > 1){
int nbParts = 2;
if(!onlyMultilevel){
Msg::Info("Mesh partition: level (%d-%d) is ZERO-GENUS (AR=%d NB=%d) ---> LAPLACIAN partition %d parts",
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 3779c6716aa73fcba9622598b5ec0401a4acf0d1..a3fb4edf1be5bb7f3a00290f1c496525bcf10da6 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -12,6 +12,8 @@
#include "fullMatrix.h"
#include <iostream>
+#define SQU(a) ((a)*(a))
+
inline double pow_int(const double &a, const int &n)
{
switch (n) {