Skip to content
Snippets Groups Projects
Commit a2f6433e authored by Thomas Toulorge's avatar Thomas Toulorge
Browse files

Changed determination of blobs for high-order mesh optimization (GUI still to be adapted)

parent d18285c5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+ 267
− 392
View file @ a2f6433e
...@@ -22,6 +22,8 @@ ...@@ -22,6 +22,8 @@
#include "FlGui.h" #include "FlGui.h"
#endif #endif
void OptHomMessage (const char *s, ...) { void OptHomMessage (const char *s, ...) {
char str[1024]; char str[1024];
va_list args; va_list args;
...@@ -39,16 +41,17 @@ void OptHomMessage (const char *s, ...) { ...@@ -39,16 +41,17 @@ void OptHomMessage (const char *s, ...) {
} }
} }
else else
fprintf(stdout,"%s\n",str); fprintf(stdout,"%s\n", str);
#else #else
fprintf(stdout,"%s\n",str); fprintf(stdout,"%s\n", str);
#endif #endif
} }
// get all elements that are neighbors
double distMaxStraight (MElement *el){ double distMaxStraight (MElement *el)
{
const polynomialBasis *lagrange = el->getFunctionSpace(); const polynomialBasis *lagrange = el->getFunctionSpace();
const polynomialBasis *lagrange1 = el->getFunctionSpace(1); const polynomialBasis *lagrange1 = el->getFunctionSpace(1);
int nV = lagrange->points.size1(); int nV = lagrange->points.size1();
...@@ -68,23 +71,76 @@ double distMaxStraight (MElement *el){ ...@@ -68,23 +71,76 @@ double distMaxStraight (MElement *el){
for (int iV = nV1; iV < nV; iV++) { for (int iV = nV1; iV < nV; iV++) {
SVector3 d = el->getVertex(iV)->point()-sxyz[iV]; SVector3 d = el->getVertex(iV)->point()-sxyz[iV];
double dx = d.norm(); double dx = d.norm();
if (dx > maxdx)maxdx=dx; if (dx > maxdx) maxdx = dx;
} }
return maxdx; return maxdx;
}
SVector3 vecMaxStraight (MElement *el)
{
const polynomialBasis *lagrange = el->getFunctionSpace();
const polynomialBasis *lagrange1 = el->getFunctionSpace(1);
int nV = lagrange->points.size1();
int nV1 = lagrange1->points.size1();
SPoint3 sxyz[256];
for (int i = 0; i < nV1; ++i) {
sxyz[i] = el->getVertex(i)->point();
}
for (int i = nV1; i < nV; ++i) {
double f[256];
lagrange1->f(lagrange->points(i, 0), lagrange->points(i, 1), lagrange->points(i, 2), f);
for (int j = 0; j < nV1; ++j)
sxyz[i] += sxyz[j] * f[j];
}
double maxdx = 0.0;
SVector3 maxVec;
for (int iV = nV1; iV < nV; iV++) {
SVector3 d = el->getVertex(iV)->point()-sxyz[iV];
double dx = d.norm();
if (dx > maxdx) {
maxdx = dx;
maxVec = d;
}
}
return maxVec;
}
static void PrintBlob (std::set<MElement*> &eles, int iTag){
char name[256];
sprintf(name,"BLOB%d.pos", iTag);
FILE *f = fopen (name,"w");
fprintf(f,"View\"%s\"{\n", name);
for (std::set<MElement*>::iterator it = eles.begin(); it != eles.end();++it){
(*it)->writePOS(f, true, false, false, false, false, false,1.0, iTag);
}
fprintf(f,"};\n");
fclose(f);
} }
void exportMeshToDassault (GModel *gm, const std::string &fn, int dim){ void exportMeshToDassault (GModel *gm, const std::string &fn, int dim){
FILE *f = fopen(fn.c_str(),"w"); FILE *f = fopen(fn.c_str(),"w");
int numVertices = gm->indexMeshVertices(true); int numVertices = gm->indexMeshVertices(true);
std::vector<GEntity*> entities; std::vector<GEntity*> entities;
gm->getEntities(entities); gm->getEntities(entities);
fprintf(f,"%d %d\n",numVertices,dim); fprintf(f,"%d %d\n", numVertices, dim);
for(unsigned int i = 0; i < entities.size(); i++) for(unsigned int i = 0; i < entities.size(); i++)
for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){ for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
MVertex *v = entities[i]->mesh_vertices[j]; MVertex *v = entities[i]->mesh_vertices[j];
if (dim == 2)fprintf(f,"%d %22.15E %22.15E\n",v->getIndex(),v->x(),v->y()); if (dim == 2)fprintf(f,"%d %22.15E %22.15E\n", v->getIndex(), v->x(), v->y());
else if (dim == 3)fprintf(f,"%d %22.15E %22.15E %22.5E\n",v->getIndex(),v->x(),v->y(),v->z()); else if (dim == 3)fprintf(f,"%d %22.15E %22.15E %22.5E\n", v->getIndex(), v->x(), v->y(), v->z());
} }
if (dim == 2){ if (dim == 2){
...@@ -95,15 +151,15 @@ void exportMeshToDassault (GModel *gm, const std::string &fn, int dim){ ...@@ -95,15 +151,15 @@ void exportMeshToDassault (GModel *gm, const std::string &fn, int dim){
nt += tris.size(); nt += tris.size();
if (tris.size())order = tris[0]->getPolynomialOrder(); if (tris.size())order = tris[0]->getPolynomialOrder();
} }
fprintf(f,"%d %d\n",nt,(order+1)*(order+2)/2); fprintf(f,"%d %d\n", nt,(order+1)*(order+2)/2);
int count = 1; int count = 1;
for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf){ for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf){
std::vector<MTriangle*> &tris = (*itf)->triangles; std::vector<MTriangle*> &tris = (*itf)->triangles;
for (size_t i=0;i<tris.size();i++){ for (size_t i=0;i<tris.size();i++){
MTriangle *t = tris[i]; MTriangle *t = tris[i];
fprintf(f,"%d ",count++); fprintf(f,"%d ", count++);
for (int j=0;j<t->getNumVertices();j++){ for (int j=0;j<t->getNumVertices();j++){
fprintf(f,"%d ",t->getVertex(j)->getIndex()); fprintf(f,"%d ", t->getVertex(j)->getIndex());
} }
fprintf(f,"\n"); fprintf(f,"\n");
} }
...@@ -113,15 +169,15 @@ void exportMeshToDassault (GModel *gm, const std::string &fn, int dim){ ...@@ -113,15 +169,15 @@ void exportMeshToDassault (GModel *gm, const std::string &fn, int dim){
std::vector<MLine*> &l = (*ite)->lines; std::vector<MLine*> &l = (*ite)->lines;
ne += l.size(); ne += l.size();
} }
fprintf(f,"%d %d\n",ne,(order+1)); fprintf(f,"%d %d\n", ne,(order+1));
count = 1; count = 1;
for (GModel::eiter ite = gm->firstEdge(); ite != gm->lastEdge(); ++ite){ for (GModel::eiter ite = gm->firstEdge(); ite != gm->lastEdge(); ++ite){
std::vector<MLine*> &l = (*ite)->lines; std::vector<MLine*> &l = (*ite)->lines;
for (size_t i=0;i<l.size();i++){ for (size_t i=0;i<l.size();i++){
MLine *t = l[i]; MLine *t = l[i];
fprintf(f,"%d ",count++); fprintf(f,"%d ", count++);
for (int j=0;j<t->getNumVertices();j++){ for (int j=0;j<t->getNumVertices();j++){
fprintf(f,"%d ",t->getVertex(j)->getIndex()); fprintf(f,"%d ", t->getVertex(j)->getIndex());
} }
fprintf(f,"%d \n",(*ite)->tag()); fprintf(f,"%d \n",(*ite)->tag());
} }
...@@ -132,242 +188,136 @@ void exportMeshToDassault (GModel *gm, const std::string &fn, int dim){ ...@@ -132,242 +188,136 @@ void exportMeshToDassault (GModel *gm, const std::string &fn, int dim){
static void getTopologicalNeighbors(int nbLayers, static std::set<MVertex *> getBndVertices(std::set<MElement*> &elements, std::map<MVertex *, std::vector<MElement*> > &vertex2elements)
const std::vector<MElement*> &all, {
const std::vector<MElement*> &elements, std::set<MVertex*> bnd;
std::set<MElement*> &result){ for (std::set<MElement*>::iterator itE = elements.begin(); itE != elements.end(); ++itE) {
for (int i = 0; i < (*itE)->getNumPrimaryVertices(); ++i) {
const std::vector<MElement*> &neighbours = vertex2elements[(*itE)->getVertex(i)];
std::set<MVertex*> touched; for (size_t k = 0; k < neighbours.size(); ++k) {
if (elements.find(neighbours[k]) == elements.end()) {
for (int i = 0; i < elements.size(); i++) for (int j = 0; j < neighbours[k]->getNumVertices(); ++j) {
for (int j=0;j<elements[i]->getNumVertices();j++) bnd.insert(neighbours[k]->getVertex(j));
touched.insert(elements[i]->getVertex(j)); }
for (int layer = 0; layer < nbLayers; layer++) {
for (int i = 0; i < all.size(); i++) {
MElement *t = all[i];
bool add_ = false;
for (int j = 0; j < t->getNumVertices(); j++) {
if (touched.find(t->getVertex(j)) != touched.end()) {
add_ = true;
} }
} }
if (add_) result.insert(t);
}
for (std::set<MElement*>::iterator it = result.begin(); it != result.end(); ++it) {
for (int j = 0; j < (*it)->getNumVertices(); j++) {
touched.insert((*it)->getVertex(j));
}
} }
} }
// printf("%d %d %d %d\n",all.size(),elements.size(), result.size(),nbLayers); return bnd;
// exit(1);
} }
static void getGeometricalNeighbors (MElement *e, const std::vector<MElement*> &all, double F, std::set<MElement*> &result) {
double R = distMaxStraight (e);
SPoint3 p = e->barycenter_infty(); static std::set<MElement*> getSurroundingBlob(MElement *el, int depth, std::map<MVertex *, std::vector<MElement*> > &vertex2elements, const double distFactor)
for (int i = 0; i < all.size(); i++) { {
MElement *e = all[i];
double dist = p.distance(e->barycenter_infty());
if (dist < R*F)result.insert(e);
}
}
static void intersection (const std::set<MElement*> &a, const std::set<MElement*> &b, std::set<MElement*> &result){ const SPoint3 p = el->barycenter_infty();
for (std::set<MElement*>::const_iterator it = a.begin() ; it != a.end() ; ++it ){ const double limDist = distMaxStraight(el)*distFactor;
if (b.find(*it) != b.end()){
result.insert(*it); std::set<MElement*> blob;
} std::list<MElement*> currentLayer, lastLayer;
}
}
static void computeVertices (const std::set<MElement*> &a, std::set<MVertex*> &result){ blob.insert(el);
for (std::set<MElement*>::const_iterator it = a.begin() ; it != a.end() ; ++it ){ lastLayer.push_back(el);
for (int i=0;i<(*it)->getNumVertices();++i){ for (int d = 0; d < depth; ++d) {
result.insert((*it)->getVertex(i)); currentLayer.clear();
for (std::list<MElement*>::iterator it = lastLayer.begin(); it != lastLayer.end(); ++it) {
for (int i = 0; i < (*it)->getNumPrimaryVertices(); ++i) {
const std::vector<MElement*> &neighbours = vertex2elements[(*it)->getVertex(i)];
for (std::vector<MElement*>::const_iterator itN = neighbours.begin(); itN != neighbours.end(); ++itN)
if (p.distance((*itN)->barycenter_infty()) < limDist)
if (blob.insert(*itN).second) currentLayer.push_back(*itN); // Assume that if an el is too far, its neighbours are too far as well
}
} }
lastLayer = currentLayer;
} }
}
static MElement * compare_worst (MElement *a, MElement *b){ return blob;
if (!a)return b;
if (!b)return a;
if (a->distoShapeMeasure() < b->distoShapeMeasure()) return a;
return b;
}
template <class ITERATOR>
MElement* getTheWorstElementDown (const ITERATOR &beg, const ITERATOR &end, double &q) {
MElement *worst = 0;
q = 1.e22;
ITERATOR it = beg;
for (;it != end;++it){
MElement *t = *it;
double jmin,jmax;
t->scaledJacRange(jmin,jmax);
if (jmin < q) {
worst = t;q = jmin;
}
}
return worst;
} }
template <class ITERATOR>
MElement* getTheWorstElementUp (const ITERATOR &beg, const ITERATOR &end, double &q) {
MElement *worst = 0;
q = -1.e22;
ITERATOR it = beg;
for (;it != end;++it){
MElement *t = *it;
double jmin,jmax;
t->scaledJacRange(jmin,jmax);
if (jmax > q) {
worst = t;q = jmax;
}
}
return worst;
}
static std::set<MVertex*> filterSimple(GEntity &ge, std::set<MElement*> &badElements, int nbLayers, std::set<MElement*> &result)
// Returns max. dist. between primary vertices of an element
static double elementSize(MElement* el)
{ {
std::vector<MElement*> allElements;
for (int i = 0; i < ge.getNumMeshElements(); ++i)
allElements.push_back(ge.getMeshElement(i));
OptHomMessage("%d bad elements", badElements.size());
getTopologicalNeighbors(nbLayers, allElements, std::vector<MElement*> (badElements.begin(), badElements.end()),result);
std::set<MVertex*> vs;
for (int i = 0; i < allElements.size(); i++) {
if (result.find(allElements[i]) == result.end()) {
for (int j = 0; j < allElements[i]->getNumVertices(); j++) {
vs.insert(allElements[i]->getVertex(j));
}
}
}
return vs;
}
std::set<MVertex*> filter2D_boundaryLayer(GFace *gf, std::vector<SPoint3> vert;
int nbLayers, vert.reserve(el->getNumPrimaryVertices());
double _qual_min, for (int i = 0; i < el->getNumPrimaryVertices(); i++) vert.push_back(el->getVertex(i)->point());
double _qual_max,
double F ,
std::set<MElement*> & badasses,
std::set<MElement*> & result
) {
double jmin, jmax;
MElement *worstDown = getTheWorstElementDown(badasses.begin(),badasses.end(),jmin);
MElement *worstUp = getTheWorstElementUp(badasses.begin(),badasses.end(),jmax);
MElement *worst;
if (jmin < _qual_min)worst = worstDown;
else worst = worstUp;
// MElement *worst = compare_worst (getTheWorstElement(gf->triangles),
// getTheWorstElement(gf->quadrangles));
std::vector<MElement*> vworst; vworst.push_back(worst);
std::vector<MElement*> all;
all.insert(all.begin(),gf->triangles.begin(),gf->triangles.end());
all.insert(all.begin(),gf->quadrangles.begin(),gf->quadrangles.end());
std::set<MElement*> result1;
getTopologicalNeighbors(nbLayers, all, vworst,result1);
std::set<MElement*> result2;
getGeometricalNeighbors (worst, all, F, result2);
intersection (result1,result2,result);
// printf("intsersection(%d,%d) = %d\n",result1.size(),result2.size(),result.size());
std::set<MVertex*> vs;
for (int i = 0; i < all.size(); i++) {
if (result.find(all[i]) == result.end()) {
for (int j = 0; j < all[i]->getNumVertices(); j++) {
vs.insert(all[i]->getVertex(j));
}
}
}
return vs;
}
double maxDistSq = 0;
for (int i = 0; i < el->getNumPrimaryVertices()-1; i++)
for (int j = i+1; j < el->getNumPrimaryVertices(); j++)
maxDistSq = std::max(maxDistSq, SVector3(vert[j]-vert[i]).normSq());
std::set<MVertex*> filter3D_boundaryLayer(GRegion *gr, return sqrt(maxDistSq);
int nbLayers,
double _qual_min,
double _qual_max,
double F,
std::set<MElement *> &result) {
double jmin,jmax;
MElement *worst = compare_worst (getTheWorstElementDown(gr->tetrahedra.begin(),gr->tetrahedra.end(),jmin),
getTheWorstElementDown(gr->prisms.begin(),gr->prisms.end(),jmin));
worst = compare_worst (worst,getTheWorstElementDown(gr->hexahedra.begin(),gr->hexahedra.end(),jmin));
std::vector<MElement*> vworst; vworst.push_back(worst);
std::vector<MElement*> all;
all.insert(all.begin(),gr->tetrahedra.begin(),gr->tetrahedra.end());
all.insert(all.begin(),gr->prisms.begin(),gr->prisms.end());
all.insert(all.begin(),gr->hexahedra.begin(),gr->hexahedra.end());
std::set<MElement*> result1;
getTopologicalNeighbors(nbLayers, all, vworst,result1);
std::set<MElement*> result2;
getGeometricalNeighbors (worst, all, F, result2);
intersection (result1,result2,result);
// printf("intsersection(%d,%d) = %d\n",result1.size(),result2.size(),result.size());
std::set<MVertex*> vs;
for (int i = 0; i < all.size(); i++) {
if (result.find(all[i]) == result.end()) {
for (int j = 0; j < all[i]->getNumVertices(); j++) {
vs.insert(all[i]->getVertex(j));
}
}
}
return vs;
}
static std::set<MVertex *> getBndVertices(std::set<MElement*> &elements, std::map<MVertex *, std::vector<MElement*> > &vertex2elements)
{
std::set<MVertex*> bnd;
for (std::set<MElement*>::iterator itE = elements.begin(); itE != elements.end(); ++itE) {
for (int i = 0; i < (*itE)->getNumPrimaryVertices(); ++i) {
const std::vector<MElement*> &neighbours = vertex2elements[(*itE)->getVertex(i)];
for (size_t k = 0; k < neighbours.size(); ++k) {
if (elements.find(neighbours[k]) == elements.end()) {
for (int j = 0; j < neighbours[k]->getNumVertices(); ++j) {
bnd.insert(neighbours[k]->getVertex(j));
}
}
}
}
}
return bnd;
} }
static std::set<MElement*> getSurroundingBlob(MElement *el, int depth, std::map<MVertex *, std::vector<MElement*> > &vertex2elements)
static std::set<MElement*> getSurroundingBlob_BL(MElement *el, int depth, std::map<MVertex *, std::vector<MElement*> > &vertex2elements, const double distFactor)
{ {
std::set<MElement *> blob;
std::set<MElement *> currentLayer, lastLayer; const SPoint3 p = el->barycenter_infty();
const double halfElSize = 0.5*elementSize(el);
//std::cout << "El. " << el->getNum() << " -> halfElSize = " << halfElSize << "\n";
SVector3 dir = vecMaxStraight(el);
dir.normalize();
//std::cout << "El. " << el->getNum() << " -> direction (" << dir(0) << ", " << dir(1) << ", " << dir(2) << ")\n";
std::set<MElement*> blob;
std::list<MElement*> currentLayer, lastLayer;
blob.insert(el); blob.insert(el);
lastLayer.insert(el); lastLayer.push_back(el);
for (int d = 0; d < depth; ++d) { for (int d = 0; d < depth; ++d) {
for (std::set<MElement *>::iterator it = lastLayer.begin(); it != lastLayer.end(); ++it) { currentLayer.clear();
for (int i = 0; i < (*it)->getNumPrimaryVertices(); ++i){ for (std::list<MElement*>::iterator it = lastLayer.begin(); it != lastLayer.end(); ++it) {
for (int i = 0; i < (*it)->getNumPrimaryVertices(); ++i) {
const std::vector<MElement*> &neighbours = vertex2elements[(*it)->getVertex(i)]; const std::vector<MElement*> &neighbours = vertex2elements[(*it)->getVertex(i)];
for (size_t k = 0; k < neighbours.size(); ++k) { for (std::vector<MElement*>::const_iterator itN = neighbours.begin(); itN != neighbours.end(); ++itN) {
if (blob.find(neighbours[k]) != blob.end()) const SVector3 barBar = (*itN)->barycenter_infty()-p;
continue; const SVector3 barProj = dir*dot(dir, barBar);
blob.insert(neighbours[k]); double projDist = (barBar-barProj).norm();
currentLayer.insert(neighbours[k]); //std::cout << "El. " << el->getNum() << ", d = " << d << ", i = " << i << ", projDist = " << projDist << "\n";
double elSize2 = elementSize(*itN);
if (projDist <= halfElSize+elSize2)
if (blob.insert(*itN).second) currentLayer.push_back(*itN); // Assume that if an el is too far, its neighbours are too far as well
} }
} }
} }
lastLayer = currentLayer; lastLayer = currentLayer;
} }
//std::cout << " -> Blob of " << blob.size() << " elements\n";
return blob; return blob;
}
static void addBlobChaintoGroup(std::set<int> &group, const std::vector<std::set<int> > &groupConnect, std::vector<bool> &todoPB, int iB)
{
if (todoPB[iB]) {
todoPB[iB] = false;
group.insert(iB);
const std::set<int> &connect = groupConnect[iB];
for (std::set<int>::iterator itBC = connect.begin(); itBC != connect.end(); ++itBC)
addBlobChaintoGroup(group, groupConnect, todoPB, *itBC);
}
} }
static std::vector<std::pair<std::set<MElement*> , std::set<MVertex*> > > getConnectedBlobs(GEntity &entity, const std::set<MElement*> &badElements, int depth)
static std::vector<std::pair<std::set<MElement*> , std::set<MVertex*> > > getConnectedBlobs(GEntity &entity, const std::set<MElement*> &badElements, int depth, const double distFactor)
{ {
// Compute vertex -> element connectivity
// std::cout << "Computing vertex -> element connectivity...\n";
std::map<MVertex*, std::vector<MElement *> > vertex2elements; std::map<MVertex*, std::vector<MElement *> > vertex2elements;
for (size_t i = 0; i < entity.getNumMeshElements(); ++i) { for (size_t i = 0; i < entity.getNumMeshElements(); ++i) {
MElement &element = *entity.getMeshElement(i); MElement &element = *entity.getMeshElement(i);
...@@ -375,45 +325,69 @@ static std::vector<std::pair<std::set<MElement*> , std::set<MVertex*> > > getCon ...@@ -375,45 +325,69 @@ static std::vector<std::pair<std::set<MElement*> , std::set<MVertex*> > > getCon
vertex2elements[element.getVertex(j)].push_back(&element); vertex2elements[element.getVertex(j)].push_back(&element);
} }
} }
std::vector<std::pair<std::set<MElement *>, std::set<MVertex*> > > result;
std::set<MElement *> badElementsDone; // Contruct primary blobs
// std::cout << "Constructing " << badElements.size() << " primary blobs...\n";
std::vector<std::set<MElement*> > primBlobs;
primBlobs.reserve(badElements.size());
for (std::set<MElement*>::const_iterator it = badElements.begin(); it != badElements.end(); ++it) { for (std::set<MElement*>::const_iterator it = badElements.begin(); it != badElements.end(); ++it) {
if(badElementsDone.find(*it) != badElementsDone.end()) primBlobs.push_back(getSurroundingBlob(*it, depth, vertex2elements, distFactor));
continue; // blobs.push_back(getSurroundingBlob_BL(*it, depth, vertex2elements, distFactor));
std::stack<MElement*> todo; }
todo.push(*it);
result.resize(result.size() + 1);
while (!todo.empty()) { // Compute blob connectivity
MElement *el = todo.top(); // std::cout << "Computing blob connectivity...\n";
todo.pop(); std::map<MElement*, std::set<int> > tags;
if(badElementsDone.find(el) == badElementsDone.end()) { std::vector<std::set<int> > blobConnect(primBlobs.size());
std::set<MElement*> blob = getSurroundingBlob(el, depth, vertex2elements); for (int iB = 0; iB < primBlobs.size(); ++iB) {
for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE) { std::set<MElement*> &blob = primBlobs[iB];
if(badElements.find(*itE) != badElements.end()) for(std::set<MElement*>::iterator itEl = blob.begin(); itEl != blob.end(); ++itEl) {
todo.push(*itE); std::set<int> &blobInd = tags[*itEl];
} if (!blobInd.empty()) {
badElementsDone.insert(el); for (std::set<int>::iterator itBS = blobInd.begin(); itBS != blobInd.end(); ++itBS) blobConnect[*itBS].insert(iB);
result.back().first.insert(blob.begin(), blob.end()); blobConnect[iB].insert(blobInd.begin(), blobInd.end());
} }
blobInd.insert(iB);
} }
result.back().second = getBndVertices(result.back().first, vertex2elements);
} }
return result;
}
// Identify groups of connected blobs
// std::cout << "Identifying groups of primary blobs...\n";
std::list<std::set<int> > groups;
std::vector<bool> todoPB(primBlobs.size(), true);
for (int iB = 0; iB < primBlobs.size(); ++iB)
if (todoPB[iB]) {
std::set<int> group;
addBlobChaintoGroup(group, blobConnect, todoPB, iB);
groups.push_back(group);
}
static void PrintBlob (std::set<MElement*> &eles, int iTag){ // Merge primary blobs according to groups
char name[256]; // std::cout << "Merging primary blobs into " << groups.size() << " blobs...\n";
sprintf(name,"BLOB%d.pos",iTag); std::list<std::set<MElement*> > blobs;
FILE *f = fopen (name,"w"); for (std::list<std::set<int> >::iterator itG = groups.begin(); itG != groups.end(); ++itG) {
fprintf(f,"View\"%s\"{\n",name); blobs.push_back(std::set<MElement*>());
for (std::set<MElement*>::iterator it = eles.begin(); it != eles.end();++it){ for (std::set<int>::iterator itB = itG->begin(); itB != itG->end(); ++itB) {
(*it)->writePOS(f,true,false,false,false,false,false,1.0,iTag); std::set<MElement*> primBlob = primBlobs[*itB];
blobs.back().insert(primBlob.begin(), primBlob.end());
}
} }
fprintf(f,"};\n");
fclose(f); // Store and compute blob boundaries
// std::cout << "Storing and computing boundaries for " << blobs.size() << " blobs...\n";
std::vector<std::pair<std::set<MElement *>, std::set<MVertex*> > > result;
for (std::list<std::set<MElement*> >::iterator itB = blobs.begin(); itB != blobs.end(); ++itB)
result.push_back(std::pair<std::set<MElement*>, std::set<MVertex*> >(*itB, getBndVertices(*itB, vertex2elements)));
// std::cout << "Done with blobs\n";
return result;
} }
void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p) void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
{ {
...@@ -421,9 +395,6 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p) ...@@ -421,9 +395,6 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
int samples = 30; int samples = 30;
// int method = Mesh::METHOD_RELAXBND | Mesh::METHOD_PHYSCOORD | Mesh::METHOD_PROJJAC;
// int method = Mesh::METHOD_FIXBND | Mesh::METHOD_PHYSCOORD | Mesh::METHOD_PROJJAC;
// int method = Mesh::METHOD_FIXBND | Mesh::METHOD_SURFCOORD | Mesh::METHOD_PROJJAC;
int method = 0; int method = 0;
if (p.method == 0) if (p.method == 0)
method = Mesh::METHOD_FIXBND | Mesh::METHOD_PROJJAC; method = Mesh::METHOD_FIXBND | Mesh::METHOD_PROJJAC;
...@@ -434,161 +405,65 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p) ...@@ -434,161 +405,65 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
else if(p.method < 0) else if(p.method < 0)
method = -p.method; method = -p.method;
// printf("p.method = %d\n",p.method);
// int method = Mesh::METHOD_SURFCOORD | Mesh::METHOD_PROJJAC;
// int method = Mesh::METHOD_RELAXBND | Mesh::METHOD_SURFCOORD | Mesh::METHOD_PROJJAC;
// int method = Mesh::METHOD_PROJJAC;
SVector3 BND(gm->bounds().max(), gm->bounds().min()); SVector3 BND(gm->bounds().max(), gm->bounds().min());
OptHomMessage("High order mesh optimizer starts"); OptHomMessage("High order mesh optimizer starts");
std::vector<GEntity*> entities;
gm->getEntities(entities);
if (p.filter == 0) { for (int i = 0; i < entities.size(); ++i) {
std::vector<GEntity*> entities;
gm->getEntities(entities);
for (int i = 0; i < entities.size(); ++i) {
double tf1 = Cpu();;
GEntity &entity = *entities[i];
if (entity.dim() != p.dim || (p.onlyVisible && !entity.getVisibility()))
continue;
OptHomMessage("Optimizing Model entity %4d...", entity.tag());
std::set<MElement*> badasses;
for (int i = 0; i < entity.getNumMeshElements();i++){
double jmin,jmax;
//entity.getMeshElement(i)->scaledJacRange(jmin,jmax);
jmin = jmax = entity.getMeshElement(i)->distoShapeMeasure();
if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX)
badasses.insert(entity.getMeshElement(i));
}
if (badasses.empty())
continue;
std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > toOptimize = getConnectedBlobs(entity, badasses, p.nbLayers);
//#pragma omp parallel for schedule(dynamic, 1)
p.SUCCESS = 1;
for (int i = 0; i < toOptimize.size(); ++i) {
PrintBlob (toOptimize[i].first, i+1);
OptHomMessage("Optimizing a blob %i/%i composed of %4d elements", i+1, toOptimize.size(), toOptimize[i].first.size());
fflush(stdout);
OptHOM temp(&entity, toOptimize[i].first, toOptimize[i].second, method);
int success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
OptHomMessage("jacobian optimization succeed, starting svd optimization");
success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC, p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax);
}
temp.mesh.updateGEntityPositions();
if (success <= 0) {
std::ostringstream ossI2;
ossI2 << "final_" << entity.tag() << "ITER_" << i << ".msh";
temp.mesh.writeMSH(ossI2.str().c_str());
}
//#pragma omp critical
p.SUCCESS = std::min(p.SUCCESS, success);
}
double DTF = Cpu()-tf1;
if (p.SUCCESS == 1)
OptHomMessage("Optimization succeeded for entity %d (CPU %g sec)", entity.tag(), DTF);
else if (p.SUCCESS == 0)
OptHomMessage("Warning : Model entity %4d has all jacobians positive but not all in the range CPU %g",entity.tag(),DTF);
else if (p.SUCCESS == -1)
OptHomMessage("Error : Model entity %4d has still negative jacobians",entity.tag());
}
}
else {
double distMaxBND, distAvgBND, minJac, maxJac;
if (p.dim == 2) {
double tf1 = Cpu();;
for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf) {
if (p.onlyVisible && !(*itf)->getVisibility())continue;
int ITER = 0;
OptHomMessage("Optimizing Model face %4d...",(*itf)->tag());
p.SUCCESS = 1;
std::set<MElement*> badasses;
for (int i=0;i<(*itf)->getNumMeshElements();i++){
double jmin,jmax;
(*itf)->getMeshElement(i)->scaledJacRange(jmin,jmax);
if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX)badasses.insert((*itf)->getMeshElement(i));
}
// printf("START WITH %d bad asses\n",badasses.size());
if (badasses.size() == 0)continue;
while (1){
std::set<MVertex*> toFix;
std::set<MElement*> toOptimize;
toFix = filter2D_boundaryLayer(*itf, p.nbLayers, p.BARRIER_MIN, p.BARRIER_MAX, p.DistanceFactor, badasses, toOptimize);
OptHOM temp(*itf, toOptimize, toFix, method);
temp.recalcJacDist();
temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
// if (minJac < 1.e2)OptHomMessage("Optimizing a blob of %4d elements minJ %12.5E -- maxJ %12.5E",(*itf)->getNumMeshElements(), minJac, maxJac);
std::ostringstream ossI;
ossI << "initial_" << (*itf)->tag() << "ITER_" << ITER << ".msh";
temp.mesh.writeMSH(ossI.str().c_str());
if (minJac > p.BARRIER_MIN && maxJac < p.BARRIER_MAX) break;
p.SUCCESS = std::min(p.SUCCESS,temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax));
// temp.recalcJacDist();
// temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
temp.mesh.updateGEntityPositions();
if (p.SUCCESS == -1) break;
ITER ++;
if (p.filter == 1 && ITER > badasses.size() * 2)break;
// std::ostringstream ossF;
// ossF << "final_" << (*itf)->tag() << ".msh";
// temp.mesh.writeMSH(ossF.str().c_str());
PrintBlob (toOptimize, ITER);
}
double DTF = Cpu()-tf1;
if (p.SUCCESS == 1){
OptHomMessage("Optimization succeeded (CPU %g sec)",DTF);
}
else if (p.SUCCESS == 0)
OptHomMessage("Warning : Model face %4d has all jacobians positive but not all in the range CPU %g",(*itf)->tag(),DTF);
else if (p.SUCCESS == -1)
OptHomMessage("Error : Model face %4d has still negative jacobians",(*itf)->tag());
Msg::Info("---------------------------------------------------------------------------------------------------------------"); double tf1 = Cpu();
} GEntity &entity = *entities[i];
exportMeshToDassault (gm,gm->getName() + "_opti.das", 2);
if (entity.dim() != p.dim || (p.onlyVisible && !entity.getVisibility())) continue;
OptHomMessage("Optimizing Model entity %4d...", entity.tag());
std::set<MElement*> badasses;
for (int i = 0; i < entity.getNumMeshElements();i++) {
double jmin, jmax;
entity.getMeshElement(i)->scaledJacRange(jmin, jmax);
if (p.BARRIER_MIN_METRIC > 0) jmax = jmin;
if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX)
badasses.insert(entity.getMeshElement(i));
} }
else if (p.dim == 3) { if (badasses.empty()) continue;
for (GModel::riter itr = gm->firstRegion(); itr != gm->lastRegion(); ++itr) {
if (p.onlyVisible && !(*itr)->getVisibility())continue; std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > toOptimize = getConnectedBlobs(entity, badasses, p.nbLayers, p.DistanceFactor);
int ITER = 0;
Msg::Info("Model region %4d is considered",(*itr)->tag()); //#pragma omp parallel for schedule(dynamic, 1)
p.SUCCESS = true; p.SUCCESS = 1;
while (1){ for (int i = 0; i < toOptimize.size(); ++i) {
std::set<MVertex*> toFix; //PrintBlob (toOptimize[i].first, i+1);
std::set<MElement*> toOptimize; OptHomMessage("Optimizing a blob %i/%i composed of %4d elements", i+1, toOptimize.size(), toOptimize[i].first.size());
fflush(stdout);
toFix = filter3D_boundaryLayer(*itr, p.nbLayers, p.BARRIER_MIN, p.BARRIER_MAX, p.DistanceFactor, toOptimize); OptHOM temp(&entity, toOptimize[i].first, toOptimize[i].second, method);
int success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
// if ((*itr)->tetrahedra.size() > 0) { if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
OptHOM temp(*itr, toOptimize, toFix, method); OptHomMessage("Jacobian optimization succeed, starting svd optimization");
double distMaxBND, distAvgBND, minJac, maxJac; success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC, p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax);
temp.recalcJacDist();
temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
if (minJac < 1.e2)Msg::Info("Optimizing a blob of %4d elements minJ %12.5E -- maxJ %12.5E",(*itr)->getNumMeshElements(), minJac, maxJac);
std::ostringstream ossI;
ossI << "initial_" << (*itr)->tag() << "ITER_" << ITER << ".msh";
temp.mesh.writeMSH(ossI.str().c_str());
if (minJac > p.BARRIER_MIN && maxJac < p.BARRIER_MAX) break;
p.SUCCESS = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
temp.recalcJacDist();
temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
temp.mesh.updateGEntityPositions();
if (!p.SUCCESS) break;
// }
ITER ++;
}
Msg::Info("Model region %4d is done (%d subdomains were computed) SUCCESS %1d",(*itr)->tag(),ITER,p.SUCCESS);
Msg::Info("----------------------------------------------------------------");
// temp.mesh.writeMSH("final.msh");
} }
temp.mesh.updateGEntityPositions();
if (success <= 0) {
std::ostringstream ossI2;
ossI2 << "final_" << entity.tag() << "ITER_" << i << ".msh";
temp.mesh.writeMSH(ossI2.str().c_str());
}
//#pragma omp critical
p.SUCCESS = std::min(p.SUCCESS, success);
} }
double DTF = Cpu()-tf1;
if (p.SUCCESS == 1)
OptHomMessage("Optimization succeeded for entity %d (CPU %g sec)", entity.tag(), DTF);
else if (p.SUCCESS == 0)
OptHomMessage("Warning : Model entity %4d has all jacobians positive but not all in the range (CPU %g sec)", entity.tag(), DTF);
else if (p.SUCCESS == -1)
OptHomMessage("Error : Model entity %4d has still negative jacobians (CPU %g sec)", entity.tag(), DTF);
} }
double t2 = Cpu(); double t2 = Cpu();
p.CPU = t2-t1; p.CPU = t2-t1;
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment