From 3e69fc3be19a52305b7440783961e2eabc066542 Mon Sep 17 00:00:00 2001 From: Amaury Johnan <amjohnen@gmail.com> Date: Mon, 24 Oct 2011 22:34:51 +0000 Subject: [PATCH] plugin now available for everybody --- Plugin/AnalyseCurvedMesh.cpp | 1177 ++++++++++++++++++++-------------- Plugin/AnalyseCurvedMesh.h | 2 +- Plugin/PluginManager.cpp | 2 +- 3 files changed, 700 insertions(+), 481 deletions(-) diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp index 2e6a1b0f9d..8b2e1adced 100644 --- a/Plugin/AnalyseCurvedMesh.cpp +++ b/Plugin/AnalyseCurvedMesh.cpp @@ -9,17 +9,21 @@ #include "JacobianBasis.h" #include "polynomialBasis.h" #include "AnalyseCurvedMesh.h" +#include "FlGui.h" +#include "Context.h" +#include "drawContext.h" +#include "OS.h" #define UNDEF_JAC_TAG -999 StringXNumber JacobianOptions_Number[] = { {GMSH_FULLRC, "Dim", NULL, -1}, - {GMSH_FULLRC, "Analysis", NULL, 1}, + /*{GMSH_FULLRC, "Analysis", NULL, 1},*/ {GMSH_FULLRC, "Effect", NULL, 6}, {GMSH_FULLRC, "MaxDepth", NULL, 10}, {GMSH_FULLRC, "JacBreak", NULL, .0}, {GMSH_FULLRC, "BezBreak", NULL, .0}, - {GMSH_FULLRC, "Tolerance", NULL, 1e-5} + /*{GMSH_FULLRC, "Tolerance", NULL, 1e-5}*/ }; extern "C" @@ -47,10 +51,11 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const //"Others are wrong or suppose to be wrong.\n\n" //"Plugin(AnalyseCurvedMesh) write in the console which elements are wrong. " //"(if labels of analysed type of elements are set visible, only wrong elements will be visible)\n\n" - "Analysis : 0 do nothing\n" + /*"Analysis : 0 do nothing\n" " +1 find invalid elements (*)\n" - " +2 compute min and max of all elements and print some statistics\n\n" - "Effect (for *) : 0 do nothing\n" + " +2 compute min and max of all elements and print some statistics\n\n"*/ + /*"Effect (for *) : 0 do nothing\n"*/ + "Effect : 0 do nothing\n" " +1 print a list of invalid elements\n" " +2 print some statistics\n" " +4 hide valid elements (for GUI)\n\n" @@ -62,7 +67,7 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const "the element is said to be invalid\n\n" "BezBreak = [0,JacBreak[ : if all Bezier coefficients are > 'BezBreak', " "the element is said to be valid\n\n" - "Tolerance = R+ , << 1 : tolerance (relatively to J_avg.) used during the computation of min and max"; + /*"Tolerance = R+ , << 1 : tolerance (relatively to J_avg.) used during the computation of min and max"*/; } // Miscellaneous method @@ -207,7 +212,8 @@ static void setJacobian(MElement *el, const JacobianBasis *jfs, fullVector<doubl } -/*void setJacobian(MElement *const *el, const JacobianBasis *jfs, fullMatrix<double> &jacobian) + +void setJacobian(MElement *const *el, const JacobianBasis *jfs, fullMatrix<double> &jacobian) { int numEl = jacobian.size2(); int numVertices = el[0]->getNumVertices(); @@ -228,37 +234,36 @@ static void setJacobian(MElement *el, const JacobianBasis *jfs, fullVector<doubl jfs->gradShapeMatX.mult(nodesX, jacobian); break; + case 2 : nodesY.resize(numVertices,numEl); interm1.resize(jacobian.size1(),jacobian.size2()); interm2.resize(jacobian.size1(),jacobian.size2()); - + for (int j = 0; j < numEl; j++) { for (int i = 0; i < numVertices; i++) { nodesX(i,j) = el[j]->getVertex(i)->x(); nodesY(i,j) = el[j]->getVertex(i)->y(); } } - + jfs->gradShapeMatX.mult(nodesX, jacobian); jfs->gradShapeMatY.mult(nodesY, interm2); jacobian.multTByT(interm2); - + jfs->gradShapeMatY.mult(nodesX, interm1); jfs->gradShapeMatX.mult(nodesY, interm2); interm1.multTByT(interm2); - + jacobian.add(interm1, -1); break; - - + case 3 : - nodesY.resize(numVertices,numEl); nodesZ.resize(numVertices,numEl); interm1.resize(jacobian.size1(),jacobian.size2()); interm2.resize(jacobian.size1(),jacobian.size2()); - + for (int j = 0; j < numEl; j++) { for (int i = 0; i < numVertices; i++) { nodesX(i,j) = el[j]->getVertex(i)->x(); @@ -266,51 +271,53 @@ static void setJacobian(MElement *el, const JacobianBasis *jfs, fullVector<doubl nodesZ(i,j) = el[j]->getVertex(i)->z(); } } - + jfs->gradShapeMatX.mult(nodesX, jacobian); jfs->gradShapeMatY.mult(nodesY, interm2); jacobian.multTByT(interm2); jfs->gradShapeMatZ.mult(nodesZ, interm2); jacobian.multTByT(interm2); - + jfs->gradShapeMatX.mult(nodesY, interm1); jfs->gradShapeMatY.mult(nodesZ, interm2); interm1.multTByT(interm2); jfs->gradShapeMatZ.mult(nodesX, interm2); interm1.multTByT(interm2); jacobian.add(interm1, 1); - + jfs->gradShapeMatX.mult(nodesZ, interm1); jfs->gradShapeMatY.mult(nodesX, interm2); interm1.multTByT(interm2); jfs->gradShapeMatZ.mult(nodesY, interm2); interm1.multTByT(interm2); jacobian.add(interm1, 1); - - + + jfs->gradShapeMatX.mult(nodesY, interm1); jfs->gradShapeMatY.mult(nodesX, interm2); interm1.multTByT(interm2); jfs->gradShapeMatZ.mult(nodesZ, interm2); interm1.multTByT(interm2); jacobian.add(interm1, -1); - + jfs->gradShapeMatX.mult(nodesZ, interm1); jfs->gradShapeMatY.mult(nodesY, interm2); interm1.multTByT(interm2); jfs->gradShapeMatZ.mult(nodesX, interm2); interm1.multTByT(interm2); jacobian.add(interm1, -1); - + jfs->gradShapeMatX.mult(nodesX, interm1); jfs->gradShapeMatY.mult(nodesZ, interm2); interm1.multTByT(interm2); jfs->gradShapeMatZ.mult(nodesY, interm2); interm1.multTByT(interm2); jacobian.add(interm1, -1); + + default : + break; } -}*/ - +} // Execution PView *GMSH_AnalyseCurvedMeshPlugin::execute(PView *v) @@ -319,159 +326,163 @@ PView *GMSH_AnalyseCurvedMeshPlugin::execute(PView *v) _dim = (int) JacobianOptions_Number[0].def; if (_dim < 0 || _dim > 3) _dim = _m->getDim(); - int analysis = (int) JacobianOptions_Number[1].def % 4; - int toDo = (int) JacobianOptions_Number[2].def % 8; - _maxDepth = (int) JacobianOptions_Number[3].def; - _jacBreak = (double) JacobianOptions_Number[4].def; - _bezBreak = (double) JacobianOptions_Number[5].def; - _tol = (double) JacobianOptions_Number[6].def; + /*int analysis = (int) JacobianOptions_Number[1].def % 2;*/ + int analysis = 1; + int toDo = (int) JacobianOptions_Number[1].def % 8; + _maxDepth = (int) JacobianOptions_Number[2].def; + _jacBreak = (double) JacobianOptions_Number[3].def; + _bezBreak = (double) JacobianOptions_Number[4].def; + _tol = 1;/*(double) JacobianOptions_Number[5].def;*/ if (analysis % 2) { + double t = Cpu(); + Msg::Info("Strating validity check..."); checkValidity(toDo); + Msg::Info("Done validity check (%fs)", Cpu()-t); } if (analysis / 2) { computeMinMax(); } } -{ - Msg::Info("AnalyseCurvedMeshPlugin : Starting analysis."); - int numBadEl = 0; - int numUncertain = 0; - int numAnalysedEl = 0; - std::vector<int> tag; - int method = (int)JacobianOptions_Number[0].def % 4; - int maxDepth = (int)JacobianOptions_Number[1].def; - - if (method < 1 || method > 2) { -#if defined(HAVE_BLAS) - method = 2; -#else - method = 1; -#endif - } - - GModel *m = GModel::current(); - switch (m->getDim()) { - - case 3 : - Msg::Info("Only 3D elements will be analyse."); - for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++) { - GRegion *r = *it; - - unsigned int numType[5] = {0, 0, 0, 0, 0}; - r->getNumMeshElements(numType); - - for(int type = 0; type < 5; type++) { - MElement *const *el = r->getStartElementType(type); - - int *a; - a = checkJacobian(el, numType[type], maxDepth, method); - numUncertain += a[0]; - numBadEl += a[1]; - numAnalysedEl += numType[type]; - delete[] a; - - /*for(unsigned int i = 0; i < numType[type]; i++) { - numAnalysedEl++; - if (checkJacobian(el[i], maxDepth) <= 0) numBadEl++; - }*/ - } - } - break; - - case 2 : - Msg::Info("Only 2D elements will be analyse."); - Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct."); - for (GModel::fiter it = m->firstFace(); it != m->lastFace(); it++) { - GFace *f = *it; - - unsigned int numType[3] = {0, 0, 0}; - f->getNumMeshElements(numType); - - for (int type = 0; type < 3; type++) { - MElement *const *el = f->getStartElementType(type); - - int *a; - a = checkJacobian(el, numType[type], maxDepth, method); - numUncertain += a[0]; - numBadEl += a[1]; - numAnalysedEl += numType[type]; - delete[] a; - - /*for (unsigned int i = 0; i < numType[type]; i++) { - numAnalysedEl++; - if (checkJacobian(el[i], maxDepth) <= 0) numBadEl++; - }*/ - } - } - break; - - case 1 : - Msg::Info("Only 1D elements will be analyse."); - Msg::Warning("1D elements must be on a y=cst & z=cst line ! If they aren't, results won't be correct."); - for (GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++) { - GEdge *e = *it; - - unsigned int numElement = e->getNumMeshElements(); - MElement *const *el = e->getStartElementType(0); - - int *a; - a = checkJacobian(el, numElement, maxDepth, method); - numUncertain += a[0]; - numBadEl += a[1]; - numAnalysedEl += numElement; - delete[] a; - - /*for (unsigned int i = 0; i < numElement; i++) { - numAnalysedEl++; - if (checkJacobian(el[i], maxDepth) <= 0) numBadEl++; - }*/ - } - break; - - default : - Msg::Error("I can't analyse any element."); - } - - - //Set all visibility of smaller dimension elements to 0 - switch (m->getDim()) { - case 2 : - for (GModel::fiter it = m->firstFace(); it != m->lastFace(); it++) { - GFace *f = *it; - - unsigned int numType[3] = {0, 0, 0}; - f->getNumMeshElements(numType); - - for (int type = 0; type < 3; type++) { - MElement *const *el = f->getStartElementType(type); - for (unsigned int i = 0; i < numType[type]; i++) { - el[i]->setVisibility(0); - } - } - } - case 1 : - for (GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++) { - GEdge *e = *it; - - unsigned int numElement = e->getNumMeshElements(); - MElement *const *el = e->getStartElementType(0); - - for (unsigned int i = 0; i < numElement; i++) { - el[i]->setVisibility(0); - } - } - break; - } - - Msg::Info("%d elements have been analysed.", numAnalysedEl); - Msg::Info("%d elements were bad.", numBadEl); - Msg::Info("%d elements were undetermined.", numUncertain); - Msg::Info("AnalyseCurvedMeshPlugin : Job finished."); - - return 0; -} +//{ +// Msg::Info("AnalyseCurvedMeshPlugin : Starting analysis."); +// int numBadEl = 0; +// int numUncertain = 0; +// int numAnalysedEl = 0; +// std::vector<int> tag; +// int method = (int)JacobianOptions_Number[0].def % 4; +// int maxDepth = (int)JacobianOptions_Number[1].def; +// +// if (method < 1 || method > 2) { +//#if defined(HAVE_BLAS) +// method = 2; +//#else +// method = 1; +//#endif +// } +// +// GModel *m = GModel::current(); +// switch (m->getDim()) { +// +// case 3 : +// Msg::Info("Only 3D elements will be analyse."); +// for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++) { +// GRegion *r = *it; +// +// unsigned int numType[5] = {0, 0, 0, 0, 0}; +// r->getNumMeshElements(numType); +// +// for(int type = 0; type < 5; type++) { +// MElement *const *el = r->getStartElementType(type); +// +// int *a; +// a = checkJacobian(el, numType[type], maxDepth, method); +// numUncertain += a[0]; +// numBadEl += a[1]; +// numAnalysedEl += numType[type]; +// delete[] a; +// +// /*for(unsigned int i = 0; i < numType[type]; i++) { +// numAnalysedEl++; +// if (checkJacobian(el[i], maxDepth) <= 0) numBadEl++; +// }*/ +// } +// } +// break; +// +// case 2 : +// Msg::Info("Only 2D elements will be analyse."); +// Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct."); +// for (GModel::fiter it = m->firstFace(); it != m->lastFace(); it++) { +// GFace *f = *it; +// +// unsigned int numType[3] = {0, 0, 0}; +// f->getNumMeshElements(numType); +// +// for (int type = 0; type < 3; type++) { +// MElement *const *el = f->getStartElementType(type); +// +// int *a; +// a = checkJacobian(el, numType[type], maxDepth, method); +// numUncertain += a[0]; +// numBadEl += a[1]; +// numAnalysedEl += numType[type]; +// delete[] a; +// +// /*for (unsigned int i = 0; i < numType[type]; i++) { +// numAnalysedEl++; +// if (checkJacobian(el[i], maxDepth) <= 0) numBadEl++; +// }*/ +// } +// } +// break; +// +// case 1 : +// Msg::Info("Only 1D elements will be analyse."); +// Msg::Warning("1D elements must be on a y=cst & z=cst line ! If they aren't, results won't be correct."); +// for (GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++) { +// GEdge *e = *it; +// +// unsigned int numElement = e->getNumMeshElements(); +// MElement *const *el = e->getStartElementType(0); +// +// int *a; +// a = checkJacobian(el, numElement, maxDepth, method); +// numUncertain += a[0]; +// numBadEl += a[1]; +// numAnalysedEl += numElement; +// delete[] a; +// +// /*for (unsigned int i = 0; i < numElement; i++) { +// numAnalysedEl++; +// if (checkJacobian(el[i], maxDepth) <= 0) numBadEl++; +// }*/ +// } +// break; +// +// default : +// Msg::Error("I can't analyse any element."); +// } +// +// +// //Set all visibility of smaller dimension elements to 0 +// switch (m->getDim()) { +// case 2 : +// for (GModel::fiter it = m->firstFace(); it != m->lastFace(); it++) { +// GFace *f = *it; +// +// unsigned int numType[3] = {0, 0, 0}; +// f->getNumMeshElements(numType); +// +// for (int type = 0; type < 3; type++) { +// MElement *const *el = f->getStartElementType(type); +// for (unsigned int i = 0; i < numType[type]; i++) { +// el[i]->setVisibility(0); +// } +// } +// } +// case 1 : +// for (GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++) { +// GEdge *e = *it; +// +// unsigned int numElement = e->getNumMeshElements(); +// MElement *const *el = e->getStartElementType(0); +// +// for (unsigned int i = 0; i < numElement; i++) { +// el[i]->setVisibility(0); +// } +// } +// break; +// } +// +// Msg::Info("%d elements have been analysed.", numAnalysedEl); +// Msg::Info("%d elements were bad.", numBadEl); +// Msg::Info("%d elements were undetermined.", numUncertain); +// Msg::Info("AnalyseCurvedMeshPlugin : Job finished."); +// +// return 0; +//} void GMSH_AnalyseCurvedMeshPlugin::checkValidity(int toDo) @@ -499,7 +510,7 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(int toDo) case 2 : Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct."); - for (GModel::fiter it = m->firstFace(); it != m->lastFace(); it++) { + for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++) { GFace *f = *it; unsigned int numType[3] = {0, 0, 0}; f->getNumMeshElements(numType); @@ -514,7 +525,7 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(int toDo) case 1 : Msg::Warning("1D elements must be on a y=cst & z=cst line ! If they aren't, results won't be correct."); - for (GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++) { + for (GModel::eiter it = _m->firstEdge(); it != _m->lastEdge(); it++) { GEdge *e = *it; unsigned int numElement = e->getNumMeshElements(); MElement *const *el = e->getStartElementType(0); @@ -535,7 +546,7 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(int toDo) } } if (toDo / 2 % 2) { - Msg::Info("Found %d invalid elements against %d valid", _numInvalid, _numValid); + Msg::Info("Found %d invalid elements and %d valid", _numInvalid, _numValid); if (_numUncertain) { Msg::Info("%d uncertain elements.", _numUncertain); if (_jacBreak < _bezBreak) { @@ -545,9 +556,16 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(int toDo) Msg::Info("Increase MaxDepth in order to decrease the number of uncertain elements."); } } + Msg::Info("%d elements analysed", _numAnalysedEl); } - if (toDo / 4 % 2) + if (toDo / 4 % 2) { + Msg::Info("Note: Valid elements are hidden. Change 'Effect' to disable this functionality."); + Msg::Info("(To revert visibility : Tools > Visibility > Interactive > Show All)"); hideValid_ShowInvalid(invalids); + CTX::instance()->mesh.changed = (ENT_ALL); + FlGui::instance()->check(); + drawContext::global()->draw(); + } } void GMSH_AnalyseCurvedMeshPlugin::hideValid_ShowInvalid(std::vector<MElement*> &invalids) @@ -577,7 +595,7 @@ void GMSH_AnalyseCurvedMeshPlugin::hideValid_ShowInvalid(std::vector<MElement*> break; case 2 : - for (GModel::fiter it = m->firstFace(); it != m->lastFace(); it++) { + for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++) { GFace *f = *it; unsigned int numType[3] = {0, 0, 0}; f->getNumMeshElements(numType); @@ -597,7 +615,7 @@ void GMSH_AnalyseCurvedMeshPlugin::hideValid_ShowInvalid(std::vector<MElement*> break; case 1 : - for (GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++) { + for (GModel::eiter it = _m->firstEdge(); it != _m->lastEdge(); it++) { GEdge *e = *it; unsigned int numElement = e->getNumMeshElements(); MElement *const *el = e->getStartElementType(0); @@ -645,7 +663,7 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax() case 2 : Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct."); - for (GModel::fiter it = m->firstFace(); it != m->lastFace(); it++) { + for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++) { GFace *f = *it; unsigned int numType[3] = {0, 0, 0}; f->getNumMeshElements(numType); @@ -660,7 +678,7 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax() case 1 : Msg::Warning("1D elements must be on a y=cst & z=cst line ! If they aren't, results won't be correct."); - for (GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++) { + for (GModel::eiter it = _m->firstEdge(); it != _m->lastEdge(); it++) { GEdge *e = *it; unsigned int numElement = e->getNumMeshElements(); MElement *const *el = e->getStartElementType(0); @@ -674,338 +692,108 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax() return; } Msg::Info("Extrema of J_avg : %g, %g", _min_Javg, _max_Javg); - Msg::Info("Minimum of min(distortion) : %g", _min_pJmin); - Msg::Info("Average of min(distortion) : %g", _avg_pJmin / _numAnalysedEl); - Msg::Info("Minimum of min(J) / max(J) : %g", _min_ratioJ); - Msg::Info("Average of min(J) / max(J) : %g", _avg_ratioJ / _numAnalysedEl); + Msg::Info("Minimum of min(distortion) : %g", _min_pJmin); + Msg::Info("Average of min(distortion) : %g", _avg_pJmin / _numAnalysedEl); + Msg::Info("Minimum of min(J) / max(J) : %g", _min_ratioJ); + Msg::Info("Average of min(J) / max(J) : %g", _avg_ratioJ / _numAnalysedEl); } -int *GMSH_AnalyseCurvedMeshPlugin::checkJacobian -(MElement *const *el, int numEl, int maxDepth, int method) -{ - int *a = new int[2]; - a[0] = a[1] = 0; - if (numEl <= 0) return a; - - switch (method) { - - case 1 : - for (int i = 0; i < numEl; i++) { - int tag = method_1_2(el[i], maxDepth); - if (tag < 0) { - a[1]++; - if (tag < -1) - Msg::Info("Bad element : %d (with tag %d)", el[i]->getNum(), tag); - else - Msg::Info("Bad element : %d", el[i]->getNum()); - } - else if (tag > 0) { - el[i]->setVisibility(0); - if (tag > 1) - Msg::Info("Good element : %d (with tag %d)", el[i]->getNum(), tag); - } - else { - a[0]++; - Msg::Info("Element %d may be bad", el[i]->getNum()); - } - } - return a; - - case 2 : - std::vector<int> tag(numEl, UNDEF_JAC_TAG); - method_2_2(el, tag, maxDepth); - - Msg::Info(" "); - Msg::Info("Bad elements :"); - for (unsigned int i = 0; i < tag.size(); i++) { - if (tag[i] < 0) { - if (tag[i] < -1) - Msg::Info("%d (with tag %d)", el[i]->getNum(), tag[i]); - else - Msg::Info("%d", el[i]->getNum()); - a[1]++; - } - } - Msg::Info(" "); - Msg::Info("Uncertain elements :"); - for (unsigned int i = 0; i < tag.size(); i++) { - if (tag[i] == 0) { - Msg::Info("%d", el[i]->getNum()); - a[0]++; - } - } - - Msg::Info(" "); - return a; - } -} -int GMSH_AnalyseCurvedMeshPlugin::method_1_1(MElement *el, int depth) +void GMSH_AnalyseCurvedMeshPlugin:: + checkValidity(MElement *const*el, int numEl, std::vector<MElement*> &invalids) { - const polynomialBasis *fs = el->getFunctionSpace(-1); - if (!fs) { - Msg::Error("Function space not implemented for type of element %d", el->getNum()); - return 0; - } - const JacobianBasis *jfs = el->getJacobianFuncSpace(-1); - if (!jfs) { - Msg::Error("Jacobian function space not implemented for type of element %d", el->getNum()); - return 0; - } - const bezierBasis *bb = jfs->bezier; - - int numSamplingPt = bb->points.size1(); - int dim = bb->points.size2(); - fullVector<double> jacobian(numSamplingPt); - - for (int i = 0; i < numSamplingPt; i++) { - double gsf[256][3]; - switch (dim) { - case 1 : - fs->df(bb->points(i,0),0,0, gsf); - break; - case 2 : - fs->df(bb->points(i,0),bb->points(i,1),0, gsf); - break; - case 3 : - fs->df(bb->points(i,0),bb->points(i,1),bb->points(i,2), gsf); - break; - default : - Msg::Error("Can't get the gradient for %dD elements.", dim); - return false; - } - double jac[3][3]; - jacobian(i) = getJacobian(gsf, jac, el); - } - - for (int i = 0; i < jacobian.size(); i++) { - if (jacobian(i) <= 0.) return -1; - } - - - fullVector<double> jacBez(jacobian.size()); - bb->matrixLag2Bez.mult(jacobian, jacBez); - - bool allPtPositive = true; - for (int i = 0; i < jacBez.size(); i++) { - if (jacBez(i) <= 0.) allPtPositive = false; - } - if (allPtPositive) return 1; - - - if (depth <= 1) { - return 0; - } - else{ - int tag = division(bb, jacBez, depth-1); - if (tag < 0) - return tag - 1; - if (tag > 0) - return tag + 1; - return tag; - } -} - - -int GMSH_AnalyseCurvedMeshPlugin::method_1_2(MElement *el, int depth) -{ - const polynomialBasis *fs = el->getFunctionSpace(-1); - if (!fs) { - Msg::Error("Function space not implemented for type of element %d", el->getNum()); - return 0; - } - const JacobianBasis *jfs = el->getJacobianFuncSpace(-1); - if (!jfs) { - Msg::Error("Jacobian function space not implemented for type of element %d", el->getNum()); - return 0; - } - const bezierBasis *bb = jfs->bezier; - - int numSamplingPt = bb->points.size1(); - int dim = bb->points.size2(); - fullVector<double> jacobian(numSamplingPt); - - - setJacobian(el, jfs, jacobian); - - //{ - // Msg::Info("Printing vector jac[%d]", el->getNum()); - // Msg::Info(" "); - // for(int I = 0; I < jacobian.size(); I++){ - // Msg::Info("%lf ", jacobian(I)); - // } - // Msg::Info(" "); - //} - - for (int i = 0; i < jacobian.size(); i++) { - if (jacobian(i) <= 0.) return -1; - } - - - fullVector<double> jacBez(jacobian.size()); - bb->matrixLag2Bez.mult(jacobian, jacBez); - - bool allPtPositive = true; - for (int i = 0; i < jacBez.size(); i++) { - if (jacBez(i) <= 0.) allPtPositive = false; - } - if (allPtPositive) return 1; - - - if (depth <= 1) { - return 0; - } - else{ - int tag = division(bb, jacBez, depth-1); - if (tag < 0) - return tag - 1; - if (tag > 0) - return tag + 1; - return tag; - } -} - - -void GMSH_AnalyseCurvedMeshPlugin::method_2_2 - (MElement *const *el, std::vector<int> &tags, int depth) -{ - const polynomialBasis *fs = el[0]->getFunctionSpace(-1); - if (!fs) { - Msg::Error("Function space not implemented for type of element %d", el[0]->getNum()); + if (numEl < 1) return; - } + const JacobianBasis *jfs = el[0]->getJacobianFuncSpace(-1); if (!jfs) { Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum()); return; } const bezierBasis *bb = jfs->bezier; - + int numSamplingPt = bb->points.size1(); int dim = bb->points.size2(); - int numEl = tags.size(); - - fullMatrix<double> jacobian(numSamplingPt, numEl); - setJacobian(el, jfs, jacobian); - + fullVector<double> jacobian(numSamplingPt); + fullVector<double> jacBez(numSamplingPt); - int numBad = 0; - - for (int i = 0; i < numEl; i++) { - bool isBad = false; - for (int j = 0; j < jacobian.size1(); j++) { - if (jacobian(j,i) <= 0.) isBad = true; + for (int k = 0; k < numEl; ++k) { + setJacobian(el[k], jfs, jacobian); + + int i; + for (i = 0; i < jacobian.size() && jacobian(i) > _jacBreak; ++i); + if (i < jacobian.size()) { + invalids.push_back(el[k]); + ++_numInvalid; + continue; } - if (isBad) { - tags[i] = -1; - numBad++; + + if (_maxDepth < 1) { + invalids.push_back(el[k]); + ++_numUncertain; + continue; } - } - - - fullMatrix<double> jacBez; - std::vector<int> correspondingEl; - - double critere = .2; // � d�finir suivant le temps de chaque op�ration - if ((double) numBad / (double) numEl < critere) {// il vaut mieux calculer les points de controle de chaque �l�ment - jacBez.resize(numSamplingPt, numEl); + bb->matrixLag2Bez.mult(jacobian, jacBez); - correspondingEl.resize(numEl); - for (int i = 0; i < numEl; i++) correspondingEl[i] = i; - } - else {// il vaut mieux ne calculer que les points de controle des �l�ments encore incertains - fullMatrix<double> newJac(numSamplingPt, numEl-numBad); - fullMatrix<double> prox1(numSamplingPt,1); - fullMatrix<double> prox2(numSamplingPt,1); - int index = 0; - correspondingEl.resize(numEl - numBad); - for (int i = 0; i < numEl; i++) { - if (tags[i] == UNDEF_JAC_TAG) { - correspondingEl[index] = i; - prox1.setAsProxy(newJac,index,1); - prox2.setAsProxy(jacobian,i,1); - prox1.copy(prox2); - } + + for (i = 0; i < jacBez.size() && jacBez(i) > _bezBreak; ++i); + if (i >= jacBez.size()) { + ++_numValid; + continue; } - jacBez.resize(numSamplingPt, numEl-numBad); - bb->matrixLag2Bez.mult(jacobian, jacBez); - } - - - int numGood = 0; - - for (int i = 0; i < jacBez.size2(); i++) { - bool allPtPositive = true; - for (int j = 0; j < jacBez.size1(); j++) { - if (jacBez(j,i) <= 0.) allPtPositive = false; + + if (_maxDepth < 2) { + invalids.push_back(el[k]); + ++_numUncertain; + continue; } - if (allPtPositive) { - tags[correspondingEl[i]] = 1; - numGood++; + else { + int result = subDivision(bb, jacBez, _maxDepth-1); + if (result < 0) { + invalids.push_back(el[k]); + ++_numInvalid; + continue; + } + if (result > 0) { + ++_numValid; + continue; + } + invalids.push_back(el[k]); + ++_numUncertain; + continue; } } - - - Msg::Warning("Not yet implemented for maxDepth > 1"); - depth = 1; - if (depth <= 1) { - for (int i = 0; i < jacBez.size2(); i++) - if (tags[correspondingEl[i]] == UNDEF_JAC_TAG) - tags[correspondingEl[i]] = 0; - return; - } - /*else{ - int tag = division(jfs, jacBez, depth-1); - if (tag < 0) - return tag - 1; - if (tag > 0) - return tag + 1; - return tag; - }*/ - } -int GMSH_AnalyseCurvedMeshPlugin::division - (const bezierBasis *jfs, const fullVector<double> &jacobian, int depth) +int GMSH_AnalyseCurvedMeshPlugin:: + subDivision(const bezierBasis *bb, const fullVector<double> &jacobian, int depth) { - if (jfs->subDivisor.size2() != jacobian.size()) { - Msg::Error("Wrong sizes in division : [%d,%d] * [%d]", - jfs->subDivisor.size1(), jfs->subDivisor.size2(), jacobian.size()); - Msg::Info(" "); - return 0; - } - - fullVector<double> newJacobian(jfs->subDivisor.size1()); - jfs->subDivisor.mult(jacobian, newJacobian); + fullVector<double> newJacobian(bb->subDivisor.size1()); + bb->subDivisor.mult(jacobian, newJacobian); - for (int i = 0; i < jfs->numDivisions; i++) { - for (int j = 0; j < jfs->numLagPts; j++) - if (newJacobian(i * jfs->points.size1() + j) <= 0.) return -1; - } - - bool allPtPositive = true; - for (int i = 0; i < newJacobian.size(); i++) { - if (newJacobian(i) <= 0.) allPtPositive = false; - } - if (allPtPositive) return 1; - - - // We don't know if the jacobian is positif everywhere or not - // We subdivide if we still can do it (if depth > 0). - if (depth <= 0) { + for (int i = 0; i < bb->numDivisions; i++) + for (int j = 0; j < bb->numLagPts; j++) + if (newJacobian(i * bb->points.size1() + j) <= _jacBreak) + return -1; + + int i = 0; + while (i < newJacobian.size() && newJacobian(i) > _bezBreak) + ++i; + if (i >= newJacobian.size()) return 1; + + if (depth <= 1) { return 0; } else{ fullVector<double> subJacobian; std::vector<int> negTag, posTag; bool zeroTag = false; - - for (int i = 0; i < jfs->numDivisions; i++) - { + + for (int i = 0; i < bb->numDivisions; i++) { subJacobian.setAsProxy(newJacobian, i * jacobian.size(), jacobian.size()); - int tag = division(jfs, subJacobian, depth-1); + int tag = subDivision(bb, subJacobian, depth-1); if (tag < 0) negTag.push_back(tag); @@ -1014,18 +802,449 @@ int GMSH_AnalyseCurvedMeshPlugin::division else zeroTag = true; } - - if (negTag.size() > 0) return max(negTag) - 1; + + if (negTag.size() > 0) + return *std::min_element(negTag.begin(), negTag.end()) - 1; if (zeroTag) return 0; - - return max(posTag) - 1; + + return *std::max_element(posTag.begin(), posTag.end()) + 1; } } - +void GMSH_AnalyseCurvedMeshPlugin:: + computeMinMax(MElement *const*el, int numEl) +{ + if (numEl < 1) + return; + + const JacobianBasis *jfs = el[0]->getJacobianFuncSpace(-1); + if (!jfs) { + Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum()); + return; + } + const bezierBasis *bb = jfs->bezier; + int numSamplingPt = bb->points.size1(); + int dim = bb->points.size2(); + + fullMatrix<double> jacobian(numSamplingPt, numEl); + fullMatrix<double> jacBez(numSamplingPt, numEl); + setJacobian(el, jfs, jacobian); + bb->matrixLag2Bez.mult(jacobian, jacBez); + + fullVector<double> prox(numSamplingPt); + + for (int k = 0; k < numEl; ++k) { + prox.setAsProxy(jacBez, k); + + double minJ, maxJ = minJ = jacobian(0,k);//, avgJac = 0; + for (int i = 1; i < jacobian.size1(); ++i) { + if (jacobian(i,k) < minJ) minJ = jacobian(i,k); + if (jacobian(i,k) > maxJ) maxJ = jacobian(i,k); + //avgJac += prox(i); + } + //avgJac /= jacobian.size1(); + + + + //int minJ, maxJ = minJ = prox(0); + // + //for (int i = 1; i < prox.size(); ++i) { + // if (prox(i) > maxJ) maxJ = prox(i); + // if (prox(i) < minJ) minJ = prox(i); + //} + + //for (int i = 0; i < jacobian.size(); i++) + //if (jacobian(i) <= _jacBreak) { + // invalids.push_back(el[k]); + // ++_numInvalid; + // continue; + //} + // + //if (_maxDepth < 1) { + // ++_numUncertain; + // continue; + //} + // + //bb->matrixLag2Bez.mult(jacobian, jacBez); + // + //int i = 0; + //while (i < jacBez.size() && jacBez(i) > _bezBreak) + // ++i; + //if (i >= jacBez.size()) { + // ++_numValid; + // continue; + //} + // + //if (_maxDepth < 2) { + // ++_numUncertain; + // continue; + //} + //else { + // int result = subDivision(bb, jacBez, _maxDepth-1); + // if (result < 0) { + // invalids.push_back(el[k]); + // ++_numInvalid; + // continue; + // } + // if (result > 0) { + // ++_numValid; + // continue; + // } + // ++_numUncertain; + // continue; + //} + } +} + +//int *GMSH_AnalyseCurvedMeshPlugin::checkJacobian +//(MElement *const *el, int numEl, int maxDepth, int method) +//{ +// int *a = new int[2]; +// a[0] = a[1] = 0; +// if (numEl <= 0) return a; +// +// switch (method) { +// +// case 1 : +// for (int i = 0; i < numEl; i++) { +// int tag = method_1_2(el[i], maxDepth); +// if (tag < 0) { +// a[1]++; +// if (tag < -1) +// Msg::Info("Bad element : %d (with tag %d)", el[i]->getNum(), tag); +// else +// Msg::Info("Bad element : %d", el[i]->getNum()); +// } +// else if (tag > 0) { +// el[i]->setVisibility(0); +// if (tag > 1) +// Msg::Info("Good element : %d (with tag %d)", el[i]->getNum(), tag); +// } +// else { +// a[0]++; +// Msg::Info("Element %d may be bad", el[i]->getNum()); +// } +// } +// return a; +// +// case 2 : +// std::vector<int> tag(numEl, UNDEF_JAC_TAG); +// method_2_2(el, tag, maxDepth); +// +// Msg::Info(" "); +// Msg::Info("Bad elements :"); +// for (unsigned int i = 0; i < tag.size(); i++) { +// if (tag[i] < 0) { +// if (tag[i] < -1) +// Msg::Info("%d (with tag %d)", el[i]->getNum(), tag[i]); +// else +// Msg::Info("%d", el[i]->getNum()); +// a[1]++; +// } +// } +// +// Msg::Info(" "); +// Msg::Info("Uncertain elements :"); +// for (unsigned int i = 0; i < tag.size(); i++) { +// if (tag[i] == 0) { +// Msg::Info("%d", el[i]->getNum()); +// a[0]++; +// } +// } +// +// Msg::Info(" "); +// return a; +// } +//} +// +// +//int GMSH_AnalyseCurvedMeshPlugin::method_1_1(MElement *el, int depth) +//{ +// const polynomialBasis *fs = el->getFunctionSpace(-1); +// if (!fs) { +// Msg::Error("Function space not implemented for type of element %d", el->getNum()); +// return 0; +// } +// const JacobianBasis *jfs = el->getJacobianFuncSpace(-1); +// if (!jfs) { +// Msg::Error("Jacobian function space not implemented for type of element %d", el->getNum()); +// return 0; +// } +// const bezierBasis *bb = jfs->bezier; +// +// int numSamplingPt = bb->points.size1(); +// int dim = bb->points.size2(); +// fullVector<double> jacobian(numSamplingPt); +// +// for (int i = 0; i < numSamplingPt; i++) { +// double gsf[256][3]; +// switch (dim) { +// case 1 : +// fs->df(bb->points(i,0),0,0, gsf); +// break; +// case 2 : +// fs->df(bb->points(i,0),bb->points(i,1),0, gsf); +// break; +// case 3 : +// fs->df(bb->points(i,0),bb->points(i,1),bb->points(i,2), gsf); +// break; +// default : +// Msg::Error("Can't get the gradient for %dD elements.", dim); +// return false; +// } +// double jac[3][3]; +// jacobian(i) = getJacobian(gsf, jac, el); +// } +// +// for (int i = 0; i < jacobian.size(); i++) { +// if (jacobian(i) <= 0.) return -1; +// } +// +// +// fullVector<double> jacBez(jacobian.size()); +// bb->matrixLag2Bez.mult(jacobian, jacBez); +// +// bool allPtPositive = true; +// for (int i = 0; i < jacBez.size(); i++) { +// if (jacBez(i) <= 0.) allPtPositive = false; +// } +// if (allPtPositive) return 1; +// +// +// if (depth <= 1) { +// return 0; +// } +// else{ +// int tag = division(bb, jacBez, depth-1); +// if (tag < 0) +// return tag - 1; +// if (tag > 0) +// return tag + 1; +// return tag; +// } +//} +// +// +//int GMSH_AnalyseCurvedMeshPlugin::method_1_2(MElement *el, int depth) +//{ +// const polynomialBasis *fs = el->getFunctionSpace(-1); +// if (!fs) { +// Msg::Error("Function space not implemented for type of element %d", el->getNum()); +// return 0; +// } +// const JacobianBasis *jfs = el->getJacobianFuncSpace(-1); +// if (!jfs) { +// Msg::Error("Jacobian function space not implemented for type of element %d", el->getNum()); +// return 0; +// } +// const bezierBasis *bb = jfs->bezier; +// +// int numSamplingPt = bb->points.size1(); +// int dim = bb->points.size2(); +// fullVector<double> jacobian(numSamplingPt); +// +// +// setJacobian(el, jfs, jacobian); +// +// { +// Msg::Info("Printing vector jac[%d]", el->getNum()); +// Msg::Info(" "); +// for(int I = 0; I < jacobian.size(); I++){ +// Msg::Info("%lf ", jacobian(I)); +// } +// Msg::Info(" "); +// } +// +// for (int i = 0; i < jacobian.size(); i++) { +// if (jacobian(i) <= 0.) return -1; +// } +// +// +// fullVector<double> jacBez(jacobian.size()); +// bb->matrixLag2Bez.mult(jacobian, jacBez); +// +// bool allPtPositive = true; +// for (int i = 0; i < jacBez.size(); i++) { +// if (jacBez(i) <= 0.) allPtPositive = false; +// } +// if (allPtPositive) return 1; +// +// +// if (depth <= 1) { +// return 0; +// } +// else{ +// int tag = division(bb, jacBez, depth-1); +// if (tag < 0) +// return tag - 1; +// if (tag > 0) +// return tag + 1; +// return tag; +// } +//} +// +// +//void GMSH_AnalyseCurvedMeshPlugin::method_2_2 +// (MElement *const *el, std::vector<int> &tags, int depth) +//{ +// const polynomialBasis *fs = el[0]->getFunctionSpace(-1); +// if (!fs) { +// Msg::Error("Function space not implemented for type of element %d", el[0]->getNum()); +// return; +// } +// const JacobianBasis *jfs = el[0]->getJacobianFuncSpace(-1); +// if (!jfs) { +// Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum()); +// return; +// } +// const bezierBasis *bb = jfs->bezier; +// +// int numSamplingPt = bb->points.size1(); +// int dim = bb->points.size2(); +// int numEl = tags.size(); +// +// fullMatrix<double> jacobian(numSamplingPt, numEl); +// setJacobian(el, jfs, jacobian); +// +// +// int numBad = 0; +// +// for (int i = 0; i < numEl; i++) { +// bool isBad = false; +// for (int j = 0; j < jacobian.size1(); j++) { +// if (jacobian(j,i) <= 0.) isBad = true; +// } +// if (isBad) { +// tags[i] = -1; +// numBad++; +// } +// } +// +// +// fullMatrix<double> jacBez; +// std::vector<int> correspondingEl; +// +// double critere = .2; // � d�finir suivant le temps de chaque op�ration +// if ((double) numBad / (double) numEl < critere) {// il vaut mieux calculer les points de controle de chaque �l�ment +// jacBez.resize(numSamplingPt, numEl); +// bb->matrixLag2Bez.mult(jacobian, jacBez); +// correspondingEl.resize(numEl); +// for (int i = 0; i < numEl; i++) correspondingEl[i] = i; +// } +// else {// il vaut mieux ne calculer que les points de controle des �l�ments encore incertains +// fullMatrix<double> newJac(numSamplingPt, numEl-numBad); +// fullMatrix<double> prox1(numSamplingPt,1); +// fullMatrix<double> prox2(numSamplingPt,1); +// int index = 0; +// correspondingEl.resize(numEl - numBad); +// for (int i = 0; i < numEl; i++) { +// if (tags[i] == UNDEF_JAC_TAG) { +// correspondingEl[index] = i; +// prox1.setAsProxy(newJac,index,1); +// prox2.setAsProxy(jacobian,i,1); +// prox1.copy(prox2); +// } +// } +// jacBez.resize(numSamplingPt, numEl-numBad); +// bb->matrixLag2Bez.mult(jacobian, jacBez); +// } +// +// +// int numGood = 0; +// +// for (int i = 0; i < jacBez.size2(); i++) { +// bool allPtPositive = true; +// for (int j = 0; j < jacBez.size1(); j++) { +// if (jacBez(j,i) <= 0.) allPtPositive = false; +// } +// if (allPtPositive) { +// tags[correspondingEl[i]] = 1; +// numGood++; +// } +// } +// +// +// Msg::Warning("Not yet implemented for maxDepth > 1"); +// depth = 1; +// if (depth <= 1) { +// for (int i = 0; i < jacBez.size2(); i++) +// if (tags[correspondingEl[i]] == UNDEF_JAC_TAG) +// tags[correspondingEl[i]] = 0; +// return; +// } +// /*else{ +// int tag = division(jfs, jacBez, depth-1); +// if (tag < 0) +// return tag - 1; +// if (tag > 0) +// return tag + 1; +// return tag; +// }*/ +// +//} +// +//int GMSH_AnalyseCurvedMeshPlugin::division +// (const bezierBasis *jfs, const fullVector<double> &jacobian, int depth) +//{ +// if (jfs->subDivisor.size2() != jacobian.size()) { +// Msg::Error("Wrong sizes in division : [%d,%d] * [%d]", +// jfs->subDivisor.size1(), jfs->subDivisor.size2(), jacobian.size()); +// Msg::Info(" "); +// return 0; +// } +// +// fullVector<double> newJacobian(jfs->subDivisor.size1()); +// jfs->subDivisor.mult(jacobian, newJacobian); +// +// for (int i = 0; i < jfs->numDivisions; i++) { +// for (int j = 0; j < jfs->numLagPts; j++) +// if (newJacobian(i * jfs->points.size1() + j) <= 0.) return -1; +// } +// +// bool allPtPositive = true; +// for (int i = 0; i < newJacobian.size(); i++) { +// if (newJacobian(i) <= 0.) allPtPositive = false; +// } +// if (allPtPositive) return 1; +// +// +// We don't know if the jacobian is positif everywhere or not +// We subdivide if we still can do it (if depth > 0). +// if (depth <= 0) { +// return 0; +// } +// else{ +// fullVector<double> subJacobian; +// std::vector<int> negTag, posTag; +// bool zeroTag = false; +// +// for (int i = 0; i < jfs->numDivisions; i++) +// { +// subJacobian.setAsProxy(newJacobian, i * jacobian.size(), jacobian.size()); +// int tag = division(jfs, subJacobian, depth-1); +// +// if (tag < 0) +// negTag.push_back(tag); +// else if (tag > 0) +// posTag.push_back(tag); +// else +// zeroTag = true; +// } +// +// if (negTag.size() > 0) return max(negTag) - 1; +// +// if (zeroTag) return 0; +// +// return max(posTag) - 1; +// } +//} +// +// +// +// /*{ Msg::Info("Printing matrix %s :", "nodesX"); int ni = nodesX.size1(); diff --git a/Plugin/AnalyseCurvedMesh.h b/Plugin/AnalyseCurvedMesh.h index 05da9831a3..ba3d966fc8 100644 --- a/Plugin/AnalyseCurvedMesh.h +++ b/Plugin/AnalyseCurvedMesh.h @@ -46,9 +46,9 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin private : void checkValidity(int toDo); void computeMinMax(); - void checkValidity(MElement *const *, int numEl, std::vector<MElement*> &invalids); void computeMinMax(MElement *const *, int numEl); void hideValid_ShowInvalid(std::vector<MElement*> &invalids); + int subDivision(const bezierBasis*, const fullVector<double>&, int depth); /*bool isJacPositive(MElement *); int method_1_1(MElement *, int depth); int method_1_2(MElement *, int depth); diff --git a/Plugin/PluginManager.cpp b/Plugin/PluginManager.cpp index 1c848ddab8..32260b3a14 100644 --- a/Plugin/PluginManager.cpp +++ b/Plugin/PluginManager.cpp @@ -172,7 +172,7 @@ void PluginManager::registerDefaultPlugins() ("Skin", GMSH_RegisterSkinPlugin())); allPlugins.insert(std::pair<std::string, GMSH_Plugin*> ("MathEval", GMSH_RegisterMathEvalPlugin())); -# if 0 // experimental (Amaury) +# if 1 // experimental (Amaury) allPlugins.insert(std::pair<std::string, GMSH_Plugin*> ("AnalyseCurvedMesh", GMSH_RegisterAnalyseCurvedMeshPlugin())); #endif -- GitLab