diff --git a/Common/Context.h b/Common/Context.h
index 8b43003b95de128629e48eda5a5f908615779484..e9121396c3ead20d5a6057fc4162cfd22a45ba16 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -236,8 +236,10 @@ class CTX {
int mouseSelection, mouseHoverMeshes, pickElements;
// disable some warnings for expert users?
int expertMode;
+#if defined(HAVE_VISUDEV)
// Enable heavy visualization capabilities (for development purpose)
int heavyVisu;
+#endif
// dynamic: equal to 1 while gmsh is printing
int printing;
// hide all unselected entities?
diff --git a/Common/Options.h b/Common/Options.h
index 53a35ad8bb72b56ac60a329035d6ee7f62334b34..12b9a3e5e2c366be1bf9cad79154e25d3234dbed 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -285,7 +285,9 @@ double opt_general_zoom_factor(OPT_ARGS_NUM);
double opt_general_expert_mode(OPT_ARGS_NUM);
double opt_general_stereo_mode(OPT_ARGS_NUM);
double opt_general_camera_mode(OPT_ARGS_NUM);
+#if defined(HAVE_VISUDEV)
double opt_general_heavy_visualization(OPT_ARGS_NUM);
+#endif
double opt_general_eye_sep_ratio(OPT_ARGS_NUM);
double opt_general_focallength_ratio(OPT_ARGS_NUM);
double opt_general_camera_aperture(OPT_ARGS_NUM);
diff --git a/Common/gmsh.cpp b/Common/gmsh.cpp
index b28ca2aa7ed376ccd194ea51b31a9c1aed213b27..7a342eb79b3be2ad021cfeb695bb274f060e66e0 100644
--- a/Common/gmsh.cpp
+++ b/Common/gmsh.cpp
@@ -430,16 +430,14 @@ GMSH_API gmshModelGetMeshElements(const int dim, const int tag,
GMSH_API gmshModelSetMeshSize(const vector_pair &dimTags, const double size)
{
if(!isInitialized()) return GMSH_ERROR(-1);
- bool ok = true;
for(unsigned int i = 0; i < dimTags.size(); i++){
int dim = dimTags[i].first, tag = dimTags[i].second;
- if(dim != 1) ok = false;
- GVertex *gv = GModel::current()->getVertexByTag(tag);
- if(gv) gv->setPrescribedMeshSizeAtVertex(size);
+ if(dim == 0){
+ GVertex *gv = GModel::current()->getVertexByTag(tag);
+ if(gv) gv->setPrescribedMeshSizeAtVertex(size);
+ }
}
- if(ok)
- return GMSH_OK;
- return GMSH_ERROR(1);
+ return GMSH_OK;
}
GMSH_API gmshModelSetTransfiniteLine(const int tag, const int nPoints,
diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index eff590d77ee18f76474cff015e5057ce9fab9308..758b8ccb18831c104859d3a14b31fcf873cf0a39 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -377,6 +377,13 @@ bool discreteDiskFace::parametrize() const
double t2 = Cpu();
Msg::Debug("Assembly done in %8.3f seconds", t2 - t1);
+
+
+ if(!myAssemblerU.sizeOfR() || !myAssemblerV.sizeOfR()){
+ Msg::Warning("No unknowns in systems to compute parametrization - skipping");
+ return false;
+ }
+
lsys_u->systemSolve();
lsys_v->systemSolve();
Msg::Debug("Systems solved");
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index d333e2a286cecd925fe0625959b3f18d7da3931f..72fbcfa2f73189279ce29b69ce2bf449b22dc1b1 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -405,7 +405,13 @@ static void filterPoints(GEdge*ge, int nMinimumPoints)
v->getParameter(0,t);
if (i != 0){
double t0;
- v0->getParameter(0,t0);
+ if (v0->onWhat()->dim() == 0){
+ // Vertex is begin point
+ t0 = ge->parFromPoint(SPoint3(v0->x(), v0->y(), v0->z()));
+ }
+ else
+ v0->getParameter(0, t0);
+
t=0.5*(t+t0);
}
double lc = F_LcB(ge, t);
diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp
index cbc5d6b0d9380b71bebbc56e2285180c0a86df8d..b73bdc3a22b7692e7883df92e01c58ecc0141182 100644
--- a/Mesh/qualityMeasuresJacobian.cpp
+++ b/Mesh/qualityMeasuresJacobian.cpp
@@ -403,6 +403,32 @@ double minICNMeasure(MElement *el,
}
void sampleIGEMeasure(MElement *el, int deg, double &min, double &max)
+{
+ fullVector<double> ige;
+ sampleIGEMeasure(el, deg, ige);
+
+ min = std::numeric_limits<double>::infinity();
+ max = -min;
+ for (int i = 0; i < ige.size(); ++i) {
+ min = std::min(min, ige(i));
+ max = std::max(max, ige(i));
+ }
+}
+
+void sampleJacobian(MElement *el, int deg, fullVector<double> &jac,
+ const fullMatrix<double> *normals)
+{
+ FuncSpaceData sampleSpace = FuncSpaceData(el, deg);
+ const JacobianBasis *jacBasis = BasisFactory::getJacobianBasis(sampleSpace);
+
+ fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
+ el->getNodesCoord(nodesXYZ);
+
+ jac.resize(jacBasis->getNumJacNodes());
+ jacBasis->getSignedJacobian(nodesXYZ, jac, normals);
+}
+
+void sampleIGEMeasure(MElement *el, int deg, fullVector<double> &ige)
{
fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
el->getNodesCoord(nodesXYZ);
@@ -425,7 +451,7 @@ void sampleIGEMeasure(MElement *el, int deg, double &min, double &max)
jacDetSpace = FuncSpaceData(el, true, deg-1, 1, &serendipFalse);
break;
default:
- Msg::Error("ICN not implemented for type of element %d", el->getType());
+ Msg::Error("IGE not implemented for type of element %d", el->getType());
return;
}
@@ -437,21 +463,12 @@ void sampleIGEMeasure(MElement *el, int deg, double &min, double &max)
fullMatrix<double> coeffMatLag(gradBasis->getNumSamplingPoints(), 9);
gradBasis->getAllGradientsFromNodes(nodesXYZ, coeffMatLag);
+ if (el->getDim() == 2) coeffMatLag.resize(coeffMatLag.size1(), 6, false);
fullMatrix<double> v;
computeCoeffLengthVectors_(coeffMatLag, v, type);
- fullVector<double> ige;
computeIGE_(coeffDeterminant, v, ige, type);
-
- min = std::numeric_limits<double>::infinity();
- max = -min;
- for (int i = 0; i < ige.size(); ++i) {
- min = std::min(min, ige(i));
- max = std::max(max, ige(i));
- }
-
- return;
}
double minSampledICNMeasure(MElement *el, int deg)//fordebug
diff --git a/Mesh/qualityMeasuresJacobian.h b/Mesh/qualityMeasuresJacobian.h
index 6c3b62992a1df10b80eed5c812d02207a1c11103..1bb4fc2a35d7dc8fcf8df15649c6fecf4e46dced 100644
--- a/Mesh/qualityMeasuresJacobian.h
+++ b/Mesh/qualityMeasuresJacobian.h
@@ -24,6 +24,10 @@ double minICNMeasure(MElement *el,
bool knownValid = false,
bool reversedOk = false);
void sampleIGEMeasure(MElement *el, int order, double &min, double &max);
+void sampleJacobian(MElement *el, int order, fullVector<double> &jac,
+ const fullMatrix<double> *normals = NULL);
+void sampleIGEMeasure(MElement *el, int order, fullVector<double> &ige);
+void sampleICNMeasure(MElement *el, int order, fullVector<double> &icn);
double minSampledICNMeasure(MElement *el, int order);//fordebug
double minSampledIGEMeasure(MElement *el, int order);//fordebug
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index c798ab731f8d45ea8a26b9f7c94ef87afa4643af..cc33b71d9966b15b7edbcae1eedab29aba0c26c3 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -7,6 +7,7 @@
#if defined(HAVE_MESH)
+
#include "AnalyseCurvedMesh.h"
#include "OS.h"
#include "Context.h"
@@ -19,6 +20,9 @@
#include <sstream>
#include <fstream>
#include "qualityMeasuresJacobian.h"
+#if defined(HAVE_VISUDEV)
+#include "BasisFactory.h"
+#endif
class bezierBasis;
@@ -30,6 +34,9 @@ StringXNumber CurvedMeshOptions_Number[] = {
{GMSH_FULLRC, "Draw PView", NULL, 0},
{GMSH_FULLRC, "Recompute", NULL, 0},
{GMSH_FULLRC, "Dimension of elements", NULL, -1}
+#if defined(HAVE_VISUDEV)
+ ,{GMSH_FULLRC, "Element to draw quality", NULL, 0}
+#endif
};
extern "C"
@@ -97,13 +104,25 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
{
_m = GModel::current();
- bool computeJac = static_cast<int>(CurvedMeshOptions_Number[0].def);
- bool computeIGE = static_cast<int>(CurvedMeshOptions_Number[1].def);
- bool computeICN = static_cast<int>(CurvedMeshOptions_Number[2].def);
- _threshold = static_cast<double>(CurvedMeshOptions_Number[3].def);
- bool drawPView = static_cast<int>(CurvedMeshOptions_Number[4].def);
- bool recompute = static_cast<bool>(CurvedMeshOptions_Number[5].def);
- int askedDim = static_cast<int>(CurvedMeshOptions_Number[6].def);
+ int computeJac = static_cast<int>(CurvedMeshOptions_Number[0].def);
+ int computeIGE = static_cast<int>(CurvedMeshOptions_Number[1].def);
+ int computeICN = static_cast<int>(CurvedMeshOptions_Number[2].def);
+ _threshold = CurvedMeshOptions_Number[3].def;
+ bool drawPView = static_cast<bool>(CurvedMeshOptions_Number[4].def);
+ bool recompute = static_cast<bool>(CurvedMeshOptions_Number[5].def);
+ int askedDim = static_cast<int>(CurvedMeshOptions_Number[6].def);
+
+#if defined(HAVE_VISUDEV)
+ _pwJac = computeJac/2;
+ _pwIGE = computeIGE/2;
+ _pwICN = computeICN/2;
+
+ _numElementToScan = static_cast<int>(CurvedMeshOptions_Number[7].def);
+ _viewOrder = 0;
+ _dataPViewJac.clear();
+ _dataPViewIGE.clear();
+ _dataPViewICN.clear();
+#endif
if (askedDim < 0 || askedDim > 4) askedDim = _m->getDim();
@@ -163,6 +182,10 @@ PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
if (printStatS) _printStatIGE();
if (printStatI) _printStatICN();
+#if defined(HAVE_VISUDEV)
+ _createPViewPointwise();
+#endif
+
// Create PView
if (drawPView)
for (int dim = 1; dim <= 3; ++dim) {
@@ -342,6 +365,10 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
_data.push_back(data_elementMinMax(el, min, max));
if (min < 0 && max < 0) ++cntInverted;
progress.next();
+
+#if defined(HAVE_VISUDEV)
+ _computePointwiseQuantities(el, normals);
+#endif
}
delete normals;
}
@@ -515,4 +542,84 @@ void GMSH_AnalyseCurvedMeshPlugin::_printStatICN()
infminI, avgminI, supminI);
}
+#if defined(HAVE_VISUDEV)
+void GMSH_AnalyseCurvedMeshPlugin::_computePointwiseQuantities(MElement *el,
+ const fullMatrix<double> *normals) {
+ if (_numElementToScan != -1 && el->getNum() != _numElementToScan) return;
+
+ if (!_type2tag[el->getType()])
+ _type2tag[el->getType()] = el->getTypeForMSH();
+ else
+ // Skip if element tag is different from previous elements of same type
+ if (_type2tag[el->getType()] != el->getTypeForMSH()) return;
+
+ static fullVector<double> tmpVector;
+ int num = el->getNum();
+
+ if (!_viewOrder) {
+// _viewOrder = std::min(10, 2 * el->getPolynomialOrder());
+ _viewOrder = 9;
+ }
+
+ if (_pwJac) {
+ jacobianBasedQuality::sampleJacobian(el, _viewOrder, tmpVector, normals);
+ std::vector<double> &vec = _dataPViewJac[num];
+ for (int j = 0; j < tmpVector.size(); ++j)
+ vec.push_back(tmpVector(j));
+ }
+
+ if (_pwIGE) {
+ jacobianBasedQuality::sampleIGEMeasure(el, _viewOrder, tmpVector);
+ std::vector<double> &vec = _dataPViewIGE[num];
+ for (int j = 0; j < tmpVector.size(); ++j)
+ vec.push_back(tmpVector(j));
+ }
+
+ if (_pwICN) {
+// jacobianBasedQuality::sampleICNMeasure(el, _viewOrder, tmpVector);
+// std::vector<double> &vec = _dataPViewICN[num];
+// for (int j = 0; j < tmpVector.size(); ++j)
+// vec.push_back(tmpVector(j));
+ }
+}
+
+void GMSH_AnalyseCurvedMeshPlugin::_setInterpolationMatrices(PView *view)
+{
+ PViewData *viewData = view->getData();
+ for (int type = 0; type < 20; ++type) {
+ if (!_type2tag[type]) continue;
+ viewData->deleteInterpolationMatrices(type);
+ const nodalBasis *fsE = BasisFactory::getNodalBasis(_type2tag[type]);
+ const polynomialBasis *pfsE = dynamic_cast<const polynomialBasis *>(fsE);
+ const int hoTag = ElementType::getTag(type, _viewOrder);
+ const nodalBasis *fsV = BasisFactory::getNodalBasis(hoTag);
+ const polynomialBasis *pfsV = dynamic_cast<const polynomialBasis *>(fsV);
+ viewData->setInterpolationMatrices(type,
+ pfsV->coefficients, pfsV->monomials,
+ pfsE->coefficients, pfsE->monomials);
+ }
+}
+
+void GMSH_AnalyseCurvedMeshPlugin::_createPViewPointwise()
+{
+ if (_pwJac) {
+ _setInterpolationMatrices(
+ new PView("Pointwise Jacobian", "ElementNodeData", _m, _dataPViewJac, 0, 1)
+ );
+ }
+
+ if (_pwIGE) {
+ _setInterpolationMatrices(
+ new PView("Pointwise IGE", "ElementNodeData", _m, _dataPViewIGE, 0, 1)
+ );
+ }
+
+ if (_pwICN) {
+ _setInterpolationMatrices(
+ new PView("Pointwise ICN", "ElementNodeData", _m, _dataPViewICN, 0, 1)
+ );
+ }
+}
+#endif
+
#endif
diff --git a/Plugin/AnalyseCurvedMesh.h b/Plugin/AnalyseCurvedMesh.h
index d1aac032a5306c5db3438f49739f82779a8e130d..cddb2796fcbfb1086125fd1c20b0588b447be9f9 100644
--- a/Plugin/AnalyseCurvedMesh.h
+++ b/Plugin/AnalyseCurvedMesh.h
@@ -22,11 +22,9 @@ private:
double _minJac, _maxJac, _minIGE, _minICN;
public:
data_elementMinMax(MElement *e,
- double minJ = 2,
- double maxJ = 0,
- double minIGE = -1,
- double minICN = -1)
- : _el(e), _minJac(minJ), _maxJac(maxJ), _minIGE(minIGE), _minICN(minICN) {}
+ double minJ = 2, double maxJ = 0,
+ double minIGE = -1, double minICN = -1)
+ : _el(e), _minJac(minJ), _maxJac(maxJ), _minIGE(minIGE), _minICN(minICN) {}
void setMinS(double r) { _minIGE = r; }
void setMinI(double r) { _minICN = r; }
MElement* element() { return _el; }
@@ -42,6 +40,17 @@ private :
GModel *_m;
double _threshold;
+#if defined(HAVE_VISUDEV)
+ // Pointwise data
+ int _numElementToScan;
+ bool _pwJac, _pwIGE, _pwICN;
+ std::map<int, std::vector<double>> _dataPViewJac;
+ std::map<int, std::vector<double>> _dataPViewIGE;
+ std::map<int, std::vector<double>> _dataPViewICN;
+ int _type2tag[20] = {0};
+ int _viewOrder = 0;
+#endif
+
// for 1d, 2d, 3d
bool _computedJac[3], _computedIGE[3], _computedICN[3];
bool _PViewJac[3], _PViewIGE[3], _PViewICN[3];
@@ -81,6 +90,12 @@ private :
void _printStatJacobian();
void _printStatIGE();
void _printStatICN();
+
+#if defined(HAVE_VISUDEV)
+ void _computePointwiseQuantities(MElement *, const fullMatrix<double> *normals);
+ void _createPViewPointwise();
+ void _setInterpolationMatrices(PView *);
+#endif
};
#endif
diff --git a/Post/PViewData.cpp b/Post/PViewData.cpp
index 59d56e5155b7a33ff0feac502e4407cc291d5402..0473d24ec2da431bc3d1266b4d80f29a5edf348e 100644
--- a/Post/PViewData.cpp
+++ b/Post/PViewData.cpp
@@ -184,6 +184,11 @@ bool PViewData::haveInterpolationMatrices(int type)
return _interpolation.count(type) ? true : false;
}
+void PViewData::deleteInterpolationMatrices(int type)
+{
+ _interpolation.erase(type);
+}
+
void PViewData::removeInterpolationScheme(const std::string &name)
{
std::map<std::string, interpolationMatrices>::iterator it = _interpolationSchemes.find(name);
diff --git a/Post/PViewData.h b/Post/PViewData.h
index 9b5449ec1227eb610707e2891c246d09dccdd6de..922083f18097416c4b61b40962f604d8b01bfd35 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -225,6 +225,7 @@ class PViewData {
const fullMatrix<double> &expGeo);
int getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p);
bool haveInterpolationMatrices(int type=0);
+ void deleteInterpolationMatrices(int type=0);
// access to global interpolation schemes
static void removeInterpolationScheme(const std::string &name);
diff --git a/demos/api/t4.cpp b/demos/api/t4.cpp
index a8e985a53da03f494e6cc72e73245e6ddf989902..7a47f21c99a91a1f81f9dcaa09351dafad7ede12 100644
--- a/demos/api/t4.cpp
+++ b/demos/api/t4.cpp
@@ -22,7 +22,6 @@ int main(int argc, char **argv)
double ccos = (-h5*R1 + e2 * hypot(h5, hypot(e2, R1))) / (h5*h5 + e2*e2);
double ssin = sqrt(1 - ccos*ccos);
- int o;
gmshModelGeoAddPoint(1, -e1-e2, 0 , 0, Lc1);
gmshModelGeoAddPoint(2, -e1-e2, h1 , 0, Lc1);
gmshModelGeoAddPoint(3, -e3-r , h1 , 0, Lc2);
diff --git a/demos/api/t5.cpp b/demos/api/t5.cpp
index 8ce20129a0bab54b0ac5198def54fa8f6784764c..9b7cbacf818641efe25e6fefee9e157673934683 100644
--- a/demos/api/t5.cpp
+++ b/demos/api/t5.cpp
@@ -62,7 +62,6 @@ int main(int argc, char **argv)
double lcar2 = .0005;
double lcar3 = .055;
- int o;
gmshModelGeoAddPoint(1, 0.5,0.5,0.5, lcar2);
gmshModelGeoAddPoint(2, 0.5,0.5,0, lcar1);
gmshModelGeoAddPoint(3, 0,0.5,0.5, lcar1);
@@ -134,7 +133,7 @@ int main(int argc, char **argv)
t, x, y, z, r, volumes.back());
}
- int v = gmshModelGeoAddVolume(186, shells)[1];
+ gmshModelGeoAddVolume(186, shells);
gmshModelAddPhysicalGroup(3, 10, {186});
gmshModelGeoSynchronize();