diff --git a/Fltk/highOrderToolsWindow.cpp b/Fltk/highOrderToolsWindow.cpp index 9205b992a8c53141460efb05f5ff5bcf4f135c52..e39cfdebc5fb6989a4016d18ae8e91615191bd0c 100644 --- a/Fltk/highOrderToolsWindow.cpp +++ b/Fltk/highOrderToolsWindow.cpp @@ -96,7 +96,6 @@ static void chooseopti_cb(Fl_Widget *w, void *data) if (elastic == 1){ o->choice[0]->deactivate(); - o->choice[1]->deactivate(); for (int i=3;i<=6;i++) o->value[i]->deactivate(); // o->push[1]->deactivate(); @@ -104,7 +103,6 @@ static void chooseopti_cb(Fl_Widget *w, void *data) else { o->push[0]->activate(); o->choice[0]->activate(); - o->choice[1]->activate(); for (int i=3;i<=6;i++) o->value[i]->activate(); // o->push[1]->activate(); @@ -137,7 +135,6 @@ static void highordertools_runelas_cb(Fl_Widget *w, void *data) p.weightFree = o->value[6]->value(); p.DistanceFactor = o->value[7]->value(); p.method = o->CAD ? (int)o->choice[0]->value() : 2; - p.filter = (int)o->choice[1]->value(); HighOrderMeshOptimizer (GModel::current(),p); printf("CPU TIME = %4f seconds\n",p.CPU); } @@ -295,12 +292,6 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize) choice[0]->menu(menu_objf); choice[0]->align(FL_ALIGN_RIGHT); - y += BH; - choice[1] = new Fl_Choice - (x,y, IW, BH, "Blob strategy"); - choice[1]->menu(menu_strategy); - choice[1]->align(FL_ALIGN_RIGHT); - y += BH; value[5] = new Fl_Value_Input (x,y, IW, BH, "W fixed"); @@ -358,7 +349,10 @@ void highOrderToolsWindow::show(bool redrawOnly) else { value[0]->value(meshOrder); butt[0]->value(!complete); - if (CAD) output[0]->value("Available"); + if (CAD) { + output[0]->value("Available"); + choice[0]->value(1); + } else { output[0]->value("Not Available"); choice[0]->deactivate(); diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp index 35b6fd192c49748d69b07a0d1fba1afd176b73f5..52ce0aaf19f3a220ffb4c25eeb6e926a374b8632 100644 --- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp +++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp @@ -80,55 +80,6 @@ double distMaxStraight (MElement *el) -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){ FILE *f = fopen(fn.c_str(),"w"); @@ -238,66 +189,6 @@ static std::set<MElement*> getSurroundingBlob(MElement *el, int depth, std::map< -// Returns max. dist. between primary vertices of an element -static double elementSize(MElement* el) -{ - - std::vector<SPoint3> vert; - vert.reserve(el->getNumPrimaryVertices()); - for (int i = 0; i < el->getNumPrimaryVertices(); i++) vert.push_back(el->getVertex(i)->point()); - - 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()); - - return sqrt(maxDistSq); - -} - - - -static std::set<MElement*> getSurroundingBlob_BL(MElement *el, int depth, std::map<MVertex *, std::vector<MElement*> > &vertex2elements, const double distFactor) -{ - - 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); - lastLayer.push_back(el); - for (int d = 0; d < depth; ++d) { - 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) { - const SVector3 barBar = (*itN)->barycenter_infty()-p; - const SVector3 barProj = dir*dot(dir, barBar); - double projDist = (barBar-barProj).norm(); -//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; - } - -//std::cout << " -> Blob of " << blob.size() << " elements\n"; - return blob; - -} - - - static void addBlobChaintoGroup(std::set<int> &group, const std::vector<std::set<int> > &groupConnect, std::vector<bool> &todoPB, int iB) { @@ -316,8 +207,10 @@ static void addBlobChaintoGroup(std::set<int> &group, const std::vector<std::set static std::vector<std::pair<std::set<MElement*> , std::set<MVertex*> > > getConnectedBlobs(GEntity &entity, const std::set<MElement*> &badElements, int depth, const double distFactor) { + OptHomMessage("Starting blob generation from %i bad elements...",badElements.size()); + // Compute vertex -> element connectivity -// std::cout << "Computing vertex -> element connectivity...\n"; + std::cout << "Computing vertex -> element connectivity...\n"; std::map<MVertex*, std::vector<MElement *> > vertex2elements; for (size_t i = 0; i < entity.getNumMeshElements(); ++i) { MElement &element = *entity.getMeshElement(i); @@ -327,17 +220,15 @@ static std::vector<std::pair<std::set<MElement*> , std::set<MVertex*> > > getCon } // Contruct primary blobs -// std::cout << "Constructing " << badElements.size() << " primary blobs...\n"; + 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) primBlobs.push_back(getSurroundingBlob(*it, depth, vertex2elements, distFactor)); -// blobs.push_back(getSurroundingBlob_BL(*it, depth, vertex2elements, distFactor)); - } // Compute blob connectivity -// std::cout << "Computing blob connectivity...\n"; + std::cout << "Computing blob connectivity...\n"; std::map<MElement*, std::set<int> > tags; std::vector<std::set<int> > blobConnect(primBlobs.size()); for (int iB = 0; iB < primBlobs.size(); ++iB) { @@ -353,7 +244,7 @@ static std::vector<std::pair<std::set<MElement*> , std::set<MVertex*> > > getCon } // Identify groups of connected blobs -// std::cout << "Identifying groups of primary blobs...\n"; + 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) @@ -364,7 +255,7 @@ static std::vector<std::pair<std::set<MElement*> , std::set<MVertex*> > > getCon } // Merge primary blobs according to groups -// std::cout << "Merging primary blobs into " << groups.size() << " blobs...\n"; + std::cout << "Merging primary blobs into " << groups.size() << " blobs...\n"; std::list<std::set<MElement*> > blobs; for (std::list<std::set<int> >::iterator itG = groups.begin(); itG != groups.end(); ++itG) { blobs.push_back(std::set<MElement*>()); @@ -375,12 +266,12 @@ static std::vector<std::pair<std::set<MElement*> , std::set<MVertex*> > > getCon } // Store and compute blob boundaries -// std::cout << "Storing and computing boundaries for " << blobs.size() << " blobs...\n"; + std::cout << "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"; + OptHomMessage("Generated %i blobs",blobs.size()); return result; @@ -433,8 +324,9 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p) //#pragma omp parallel for schedule(dynamic, 1) p.SUCCESS = 1; + p.minJac = 1.e100; + p.maxJac = -1.e100; 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); @@ -443,6 +335,11 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p) 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); } + double minJac, maxJac, distMaxBND, distAvgBND; + temp.recalcJacDist(); + temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND); + p.minJac = std::min(p.minJac,minJac); + p.maxJac = std::max(p.maxJac,maxJac); temp.mesh.updateGEntityPositions(); if (success <= 0) { std::ostringstream ossI2; diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h index 29bf7bd3f9720d0588b353a2cd5344075b638e05..1e1c3a082c12031bd23014888a5a1c081dcb39e5 100644 --- a/contrib/HighOrderMeshOptimizer/OptHomRun.h +++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h @@ -19,7 +19,6 @@ struct OptHomParameters { double DistanceFactor; // filter elements such that no elements further away than // DistanceFactor times the max distance to straight sided version of an element are optimized int method ; // how jacobians are computed and if points can move on boundaries - int filter ; // 0--> standard 1--> boundary layers // OUTPUT ------> int SUCCESS ; // 0 --> success , 1 --> Not converged double minJac, maxJac; // after optimization, range of jacobians @@ -31,7 +30,7 @@ struct OptHomParameters { OptHomParameters () // default values : BARRIER_MIN_METRIC (-1.), BARRIER_MIN (0.1), BARRIER_MAX (2.0) , weightFixed (1.e6), weightFree (1.e2), - nbLayers (6) , dim(3) , itMax(10000), onlyVisible(true), DistanceFactor(12), method(1), filter(1) + nbLayers (6) , dim(3) , itMax(10000), onlyVisible(true), DistanceFactor(12), method(1) { } };