diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 89f2f13e2fdba3f8700dc8419abc1eb9788d12de..865f019b25c5f38c4e99cdaed964b05a75737b72 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -392,7 +392,7 @@ bool GFaceCompound::checkFolding(std::vector<MVertex*> &ordered) const
//the same normal orientation
bool GFaceCompound::checkOrientation(int iter) const
{
-
+
//Only check orientation for stl files (1 patch)
// if(_compound.size() > 1.0) return true;
@@ -433,7 +433,8 @@ bool GFaceCompound::checkOrientation(int iter) const
return checkOrientation(iter+1);
}
else if (oriented && iter < iterMax){
- Msg::Debug("Parametrization is 1 to 1 :-)");
+ Msg::Info("Parametrization is 1 to 1 :-)");
+ //printStuff();
}
return oriented;
@@ -507,7 +508,6 @@ bool GFaceCompound::parametrize() const
//Multiscale Laplace parametrization
//-----------------
else if (_mapping == MULTISCALE){
- // printf("multiscale exit \n");
std::vector<MElement*> _elements;
for( std::list<GFace*>::const_iterator itt = _compound.begin(); itt != _compound.end(); ++itt)
for(unsigned int i = 0; i < (*itt)->triangles.size(); ++i)
@@ -522,15 +522,13 @@ bool GFaceCompound::parametrize() const
bool withoutFolding = parametrize_conformal_spectral() ;
//bool withoutFolding = parametrize_conformal();
if ( withoutFolding == false ){
- printStuff(); //exit(1);
+ //buildOct(); exit(1);
Msg::Warning("$$$ Parametrization switched to harmonic map");
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
}
}
- printStuff();
-
buildOct();
if (!checkOrientation(0)){
@@ -1001,12 +999,28 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
return;
}
+ //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];
if(step == ITERU) myAssembler.fixVertex(v, 0, 1, cos(theta));
else if(step == ITERV) myAssembler.fixVertex(v, 0, 1, sin(theta));
}
+
+ //pin down two vertices
+ // MVertex *v1 = ordered[0];
+ // MVertex *v2 = ordered[(int)ceil((double)ordered.size()/2.)];
+ // if(step == ITERU){
+ // myAssembler.fixVertex(v1, 0, 1, 0.);
+ // myAssembler.fixVertex(v2, 0, 1, 1.);
+ // }
+ // else if(step == ITERV){
+ // myAssembler.fixVertex(v1, 0, 1, 0.);
+ // myAssembler.fixVertex(v2, 0, 1, 0.);
+ // }
+ // printf("Pinned vertex %g %g %g / %g %g %g \n", v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z());
+ //exit(1);
+
}
else if(_type == SQUARE){
if(step == ITERU){
@@ -1250,6 +1264,11 @@ bool GFaceCompound::parametrize_conformal() const
MVertex *v1 = ordered[0];
MVertex *v2 = ordered[(int)ceil((double)ordered.size()/2.)];
+ myAssembler.fixVertex(v1, 0, 1, 0.);
+ myAssembler.fixVertex(v1, 0, 2, 0.);
+ myAssembler.fixVertex(v2, 0, 1, 1.);
+ myAssembler.fixVertex(v2, 0, 2, 0.);
+//printf("Pinned vertex %g %g %g / %g %g %g \n", v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z());
// MVertex *v2 ;
// double maxSize = 0.0;
@@ -1261,13 +1280,6 @@ bool GFaceCompound::parametrize_conformal() const
// maxSize = dist;
// }
// }
-
- myAssembler.fixVertex(v1, 0, 1, 0);//0
- myAssembler.fixVertex(v1, 0, 2, 0);//0
- myAssembler.fixVertex(v2, 0, 1, 1);//1
- myAssembler.fixVertex(v2, 0, 2, 0);//0
-
- //printf("Pinned vertex %g %g %g / %g %g %g \n", v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z());
std::list<GFace*>::const_iterator it = _compound.begin();
for( ; it != _compound.end(); ++it){
@@ -1513,11 +1525,11 @@ GPoint GFaceCompound::point(double par1, double par2) const
Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
{
+
if(!oct) parametrize();
-
- if(trivial()){
+
+ if(trivial())
return (*(_compound.begin()))->firstDer(param);
- }
double U, V;
GFaceCompoundTriangle *lt;
@@ -1590,7 +1602,7 @@ void GFaceCompound::secondDer(const SPoint2 ¶m,
//printf("der second dvdv = %g %g %g \n", dvdv->x(), dvdv->y(), dvdv->z());
//printf("der second dudv = %g %g %g \n", dudv->x(), dudv->y(), dudv->z());
- Msg::Error("Computation of the second derivatives is not implemented for compound faces");
+ Msg::Warning("Computation of the second derivatives is not implemented for compound faces");
}
@@ -1790,6 +1802,9 @@ void GFaceCompound::buildOct() const
}
nbT = count;
Octree_Arrange(oct);
+
+ printStuff();
+
}
bool GFaceCompound::checkTopology() const
@@ -2003,6 +2018,7 @@ void GFaceCompound::printStuff() const
char name0[256], name1[256], name2[256], name3[256];
char name4[256], name5[256], name6[256];
+ char name7[256], name8[256], name9[256];
sprintf(name0, "UVAREA-%d.pos", (*it)->tag());
sprintf(name1, "UVX-%d.pos", (*it)->tag());
sprintf(name2, "UVY-%d.pos", (*it)->tag());
@@ -2011,6 +2027,10 @@ void GFaceCompound::printStuff() const
sprintf(name5, "XYZV-%d.pos", (*it)->tag());
sprintf(name6, "XYZC-%d.pos", (*it)->tag());
+ sprintf(name7, "UVME-%d.pos", (*it)->tag());
+ sprintf(name8, "UVMF-%d.pos", (*it)->tag());
+ sprintf(name9, "UVMG-%d.pos", (*it)->tag());
+
FILE * uva = fopen(name0,"w");
FILE * uvx = fopen(name1,"w");
FILE * uvy = fopen(name2,"w");
@@ -2018,6 +2038,10 @@ void GFaceCompound::printStuff() const
FILE * xyzu = fopen(name4,"w");
FILE * xyzv = fopen(name5,"w");
FILE * xyzc = fopen(name6,"w");
+ FILE * uvme = fopen(name7,"w");
+ FILE * uvmf = fopen(name8,"w");
+ FILE * uvmg = fopen(name9,"w");
+
fprintf(uva,"View \"\"{\n");
fprintf(uvx,"View \"\"{\n");
@@ -2026,6 +2050,10 @@ void GFaceCompound::printStuff() const
fprintf(xyzu,"View \"\"{\n");
fprintf(xyzv,"View \"\"{\n");
fprintf(xyzc,"View \"\"{\n");
+ fprintf(uvme,"View \"\"{\n");
+ fprintf(uvmf,"View \"\"{\n");
+ fprintf(uvmg,"View \"\"{\n");
+
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
@@ -2063,12 +2091,43 @@ void GFaceCompound::printStuff() const
double q1[3] = {it1->second.x(), it1->second.y(), 0.0};
double q2[3] = {it2->second.x(), it2->second.y(), 0.0};
double a_2D = fabs(triangle_area(q0, q1, q2));
- double area = a_2D/a_3D;
+ double area = (a_3D/a_2D); //*(a_3D/a_2D);
+
+ Pair<SVector3, SVector3> der = this->firstDer(SPoint2( it0->second.x(), it0->second.y()));
+ double metric0e = dot(der.first(), der.first());
+ double metric0f = dot(der.second()*(1./norm(der.second())), der.first()*(1./norm(der.first())));
+ double metric0g = dot(der.second(), der.second());
+ Pair<SVector3, SVector3> der1 = this->firstDer(SPoint2( it1->second.x(), it1->second.y()));
+ double metric1e = dot(der1.first(), der1.first());
+ double metric1f = dot(der1.second()*(1/norm(der1.second())), der1.first()*(1./norm(der1.first())));
+ double metric1g = dot(der1.second(), der1.second());
+ Pair<SVector3, SVector3> der2 = this->firstDer(SPoint2( it2->second.x(), it2->second.y()));
+ double metric2e = dot(der2.first(), der2.first());
+ double metric2f = dot(der2.second()*(1./norm(der2.second())), der2.first()*(1./norm(der2.first())));
+ double metric2g = dot(der2.second(), der2.second());
+
fprintf(uva,"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,
- area, area, area);
+ area, area, area);
+
+ fprintf(uvme,"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,
+ 0.3*(metric0e+metric1e+metric2e), 0.3*(metric0e+metric1e+metric2e), 0.3*(metric0e+metric1e+metric2e) );
+ fprintf(uvmf,"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,
+ 0.3*(metric0f+metric1f+metric2f), 0.3*(metric0f+metric1f+metric2f), 0.3*(metric0f+metric1f+metric2f) );
+ fprintf(uvmg,"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,
+ 0.3*(metric0g+metric1g+metric2g), 0.3*(metric0g+metric1g+metric2g), 0.3*(metric0g+metric1g+metric2g) );
+
fprintf(uvx,"ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E){%22.15E,%22.15E,%22.15E};\n",
it0->second.x(), it0->second.y(), 0.0,
it1->second.x(), it1->second.y(), 0.0,
@@ -2100,6 +2159,13 @@ void GFaceCompound::printStuff() const
fclose(xyzv);
fprintf(xyzc,"};\n");
fclose(xyzc);
+ fprintf(uvme,"};\n");
+ fclose(uvme);
+ fprintf(uvmf,"};\n");
+ fclose(uvmf);
+ fprintf(uvmg,"};\n");
+ fclose(uvmg);
+
}
#endif
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index e364b460b0fc9adc5aa5c1658914ac29f19f6955..2acb8f07add77083c5c3fc50c3e27d7a23952d97 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1218,7 +1218,9 @@ void GModel::createTopologyFromMesh()
for (int ifa = 0; ifa < nbFaces; ifa++){
std::vector<MElement*> myElements = conRegions[ifa];
- int numF = maxFaceNum()+1;
+ int numF;
+ if (nbFaces == 1) numF = (*itF)->tag();
+ else numF = maxFaceNum()+1;
discreteFace *f = new discreteFace(this, numF);
//printf("*** Created discreteFace %d \n", numF);
add(f);
@@ -1252,11 +1254,52 @@ void GModel::createTopologyFromMesh()
discFaces = newDiscFaces;
//------
- //EMI-FIX in case of createTopology for Volumes
- //all faces are set to each volume
- for (std::vector<discreteRegion*>::iterator it = discRegions.begin(); it != discRegions.end(); it++)
- (*it)->setBoundFaces();
-
+ //set boundary faces for each volume
+
+ // find boundary faces of each region and put them in a map_faces that
+ // associates the MElements with the tags of the adjacent regions
+ std::map<MFace, std::vector<int>, Less_Face > map_faces;
+ for (std::vector<discreteRegion*>::iterator it = discRegions.begin();
+ it != discRegions.end(); it++){
+ (*it)->findFaces(map_faces);
+ }
+
+ // create reverse map, for each region find set of MFaces that are
+ // candidate for new discrete Face
+ std::map<int, std::set<int> > region2Faces;
+
+ for (std::vector<discreteFace*>::iterator itF = discFaces.begin();
+ itF != discFaces.end(); itF++){
+ for (unsigned int i = 0; i < (*itF)->getNumMeshElements(); i++){
+ MFace mf = (*itF)->getMeshElement(i)->getFace(0);
+ std::map<MFace, std::vector<int>, Less_Face >::iterator itset = map_faces.find(mf);
+ if (itset != map_faces.end()){
+ std::vector<int> tagRegions = itset->second;
+ for (std::vector<int>::iterator itR =tagRegions.begin(); itR != tagRegions.end(); itR++){
+ std::map<int, std::set<int> >::iterator itmap = region2Faces.find(*itR);
+ if (itmap == region2Faces.end()){
+ std::set<int> oneFace; oneFace.insert((*itF)->tag());
+ region2Faces.insert(std::make_pair(*itR, oneFace));
+ }
+ else{
+ std::set<int> allFaces = itmap->second;
+ allFaces.insert(allFaces.begin(), (*itF)->tag());
+ itmap->second = allFaces;
+ }
+ }
+ }
+ }
+ }
+
+ for (std::vector<discreteRegion*>::iterator it = discRegions.begin(); it != discRegions.end(); it++){
+ std::map<int, std::set<int> >::iterator itr = region2Faces.find((*it)->tag());
+ if (itr != region2Faces.end()){
+ std::set<int> bcFaces = itr->second;
+ (*it)->setBoundFaces(bcFaces);
+ }
+ }
+
+ //create Topo for faces
createTopologyFromFaces(discFaces);
double t2 = Cpu();
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index e7243bbb8e98b47f9a855087d206d728ef7ee82d..783b666671ad05f0dffbef16c88b9146a05b6b87 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2594,6 +2594,7 @@ void ExtrudeShapes(int type, List_T *list_in,
&pv, e);
Surface *ps = FindSurface(top.Num);
top.Type = ps ? ps->Typ : 0;
+
List_Add(list_out, &top);
if(pv){
Shape body;
diff --git a/Geo/discreteRegion.cpp b/Geo/discreteRegion.cpp
index a61203cdda26d9d3260c79bbf6035ef6ae5ecb08..fda3cd105cd35e5527f3e2827f8d47455b957342 100644
--- a/Geo/discreteRegion.cpp
+++ b/Geo/discreteRegion.cpp
@@ -14,13 +14,52 @@ discreteRegion::discreteRegion(GModel *model, int num) : GRegion(model, num)
//Tree_Add(model->getGEOInternals()->Volumes, &v);
}
-void discreteRegion::setBoundFaces()
+void discreteRegion::setBoundFaces(std::set<int> tagFaces)
{
+
+ for (std::set<int>::iterator it = tagFaces.begin() ; it != tagFaces.end();++it){
+ GFace *face = model()->getFaceByTag(*it);
+ l_faces.push_back(face);
+ face->addRegion(this);
+ }
+
//in case discrete region already exist
//to modify to take into account appropriate faces
- for(GModel::fiter face = model()->firstFace(); face != model()->lastFace(); face++){
- l_faces.push_back(*face);
- (*face)->addRegion(this);
+ // for(GModel::fiter face = model()->firstFace(); face != model()->lastFace(); face++){
+ // l_faces.push_back(*face);
+ // (*face)->addRegion(this);
+ // }
+
+
+}
+void discreteRegion::findFaces(std::map<MFace, std::vector<int>, Less_Face> &map_faces)
+{
+
+ std::set<MFace, Less_Face> bound_faces;
+ for (unsigned int elem = 0; elem < getNumMeshElements() ; elem++) {
+ MElement *e = getMeshElement(elem);
+ for (int iFace = 0; iFace < e->getNumFaces(); iFace++) {
+ MFace tmp_face = e->getFace(iFace);
+ std::set<MFace, Less_Face >::iterator itset = bound_faces.find(tmp_face);
+ if (itset == bound_faces.end()) bound_faces.insert(tmp_face);
+ else bound_faces.erase(itset);
+ }
+ }
+
+ // for the boundary faces, associate the tag of the current discrete face
+ for (std::set<MFace, Less_Face>::iterator itv = bound_faces.begin();
+ itv != bound_faces.end(); ++itv){
+ std::map<MFace, std::vector<int>, Less_Face >::iterator itmap = map_faces.find(*itv);
+ if (itmap == map_faces.end()){
+ std::vector<int> tagRegions;
+ tagRegions.push_back(tag());
+ map_faces.insert(std::make_pair(*itv, tagRegions));
+ }
+ else{
+ std::vector<int> tagRegions = itmap->second;
+ tagRegions.push_back(tag());
+ itmap->second = tagRegions;
+ }
}
}
diff --git a/Geo/discreteRegion.h b/Geo/discreteRegion.h
index f78603099777bf51358367656682f9eed9726121..250f61154d2b6aa93e31dc1fefd18e986d370341 100644
--- a/Geo/discreteRegion.h
+++ b/Geo/discreteRegion.h
@@ -8,13 +8,15 @@
#include "GModel.h"
#include "GRegion.h"
+#include "MFace.h"
class discreteRegion : public GRegion {
public:
discreteRegion(GModel *model, int num);
virtual ~discreteRegion() {}
virtual GeomType geomType() const { return DiscreteVolume; }
- void setBoundFaces();
+ void setBoundFaces(std::set<int> tagFaces);
+ void findFaces(std::map<MFace, std::vector<int>, Less_Face> &map_faces);
};
#endif