Skip to content
Snippets Groups Projects
Commit 3e69fc3b authored by Amaury Johnen's avatar Amaury Johnen
Browse files

plugin now available for everybody

parent 15687030
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Plugin/AnalyseCurvedMesh.cpp
+ 698
− 479
View file @ 3e69fc3b
......@@ -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,6 +234,7 @@ 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());
......@@ -251,9 +258,7 @@ static void setJacobian(MElement *el, const JacobianBasis *jfs, fullVector<doubl
jacobian.add(interm1, -1);
break;
case 3 :
nodesY.resize(numVertices,numEl);
nodesZ.resize(numVertices,numEl);
interm1.resize(jacobian.size1(),jacobian.size2());
......@@ -308,9 +313,11 @@ static void setJacobian(MElement *el, const JacobianBasis *jfs, fullVector<doubl
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);
......@@ -680,202 +698,127 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax()
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 (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->getNum());
return 0;
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();
fullVector<double> jacobian(numSamplingPt);
fullVector<double> jacBez(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 k = 0; k < numEl; ++k) {
setJacobian(el[k], jfs, jacobian);
for (int i = 0; i < jacobian.size(); i++) {
if (jacobian(i) <= 0.) return -1;
int i;
for (i = 0; i < jacobian.size() && jacobian(i) > _jacBreak; ++i);
if (i < jacobian.size()) {
invalids.push_back(el[k]);
++_numInvalid;
continue;
}
if (_maxDepth < 1) {
invalids.push_back(el[k]);
++_numUncertain;
continue;
}
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;
for (i = 0; i < jacBez.size() && jacBez(i) > _bezBreak; ++i);
if (i >= jacBez.size()) {
++_numValid;
continue;
}
if (allPtPositive) return 1;
if (depth <= 1) {
return 0;
if (_maxDepth < 2) {
invalids.push_back(el[k]);
++_numUncertain;
continue;
}
else {
int tag = division(bb, jacBez, depth-1);
if (tag < 0)
return tag - 1;
if (tag > 0)
return tag + 1;
return tag;
int result = subDivision(bb, jacBez, _maxDepth-1);
if (result < 0) {
invalids.push_back(el[k]);
++_numInvalid;
continue;
}
if (result > 0) {
++_numValid;
continue;
}
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;
invalids.push_back(el[k]);
++_numUncertain;
continue;
}
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;
}
int GMSH_AnalyseCurvedMeshPlugin::
subDivision(const bezierBasis *bb, const fullVector<double> &jacobian, int depth)
{
fullVector<double> newJacobian(bb->subDivisor.size1());
bb->subDivisor.mult(jacobian, newJacobian);
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;
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{
int tag = division(bb, jacBez, depth-1);
fullVector<double> subJacobian;
std::vector<int> negTag, posTag;
bool zeroTag = false;
for (int i = 0; i < bb->numDivisions; i++) {
subJacobian.setAsProxy(newJacobian, i * jacobian.size(), jacobian.size());
int tag = subDivision(bb, subJacobian, depth-1);
if (tag < 0)
return tag - 1;
if (tag > 0)
return tag + 1;
return tag;
negTag.push_back(tag);
else if (tag > 0)
posTag.push_back(tag);
else
zeroTag = true;
}
if (negTag.size() > 0)
return *std::min_element(negTag.begin(), negTag.end()) - 1;
if (zeroTag) return 0;
return *std::max_element(posTag.begin(), posTag.end()) + 1;
}
}
void GMSH_AnalyseCurvedMeshPlugin::method_2_2
(MElement *const *el, std::vector<int> &tags, int depth)
void GMSH_AnalyseCurvedMeshPlugin::
computeMinMax(MElement *const*el, int numEl)
{
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());
......@@ -885,147 +828,423 @@ void GMSH_AnalyseCurvedMeshPlugin::method_2_2
int numSamplingPt = bb->points.size1();
int dim = bb->points.size2();
int numEl = tags.size();
fullMatrix<double> jacobian(numSamplingPt, numEl);
fullMatrix<double> jacBez(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; // dfinir suivant le temps de chaque opration
if ((double) numBad / (double) numEl < critere) {// il vaut mieux calculer les points de controle de chaque lment
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 lments 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);
}
fullVector<double> prox(numSamplingPt);
int numGood = 0;
for (int k = 0; k < numEl; ++k) {
prox.setAsProxy(jacBez, k);
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++;
}
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();
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);
//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 < jfs->numDivisions; i++) {
for (int j = 0; j < jfs->numLagPts; j++)
if (newJacobian(i * jfs->points.size1() + j) <= 0.) return -1;
//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;
//}
}
bool allPtPositive = true;
for (int i = 0; i < newJacobian.size(); i++) {
if (newJacobian(i) <= 0.) allPtPositive = false;
}
if (allPtPositive) return 1;
//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; // dfinir suivant le temps de chaque opration
// if ((double) numBad / (double) numEl < critere) {// il vaut mieux calculer les points de controle de chaque lment
// 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 lments 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;
}
}
// 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();
......
Plugin/AnalyseCurvedMesh.h
+ 1
− 1
View file @ 3e69fc3b
......@@ -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);
......
Plugin/PluginManager.cpp
+ 1
− 1
View file @ 3e69fc3b
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment