1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle 2 // 3 // See the LICENSE.txt file in the Gmsh root directory for license information. 4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues. 5 6 #ifndef MESH_METRIC_H 7 #define MESH_METRIC_H 8 9 #include <map> 10 #include <algorithm> 11 #include "STensor3.h" 12 #include "Field.h" 13 #include "meshGFaceOptimize.h" 14 15 template <class scalar> class simpleFunction; 16 class MVertex; 17 class gLevelset; 18 class MElementOctree; 19 class STensor3; 20 21 /**Anisotropic mesh size field based on a metric */ 22 class meshMetric : public Field { 23 public: 24 typedef enum { 25 LEVELSET = 1, 26 HESSIAN = 2, 27 FREY = 3, 28 EIGENDIRECTIONS = 4, 29 EIGENDIRECTIONS_LINEARINTERP_H = 5, 30 ISOTROPIC_LINEARINTERP_H = 6 31 } MetricComputationTechnique; 32 33 private: 34 // intersect all metrics added in "setOfMetrics", preserve eigendirections of 35 // the "most anisotropic" metric 36 void updateMetrics(); 37 int _dim; 38 double _epsilon, _e, _e_moins, _np; 39 bool needMetricUpdate; 40 bool hasAnalyticalMetric; 41 meshMetric::MetricComputationTechnique _technique; 42 double hmin, hmax; 43 simpleFunction<double> *_fct; 44 45 std::vector<MElement *> _elements; 46 v2t_cont _adj; 47 MElementOctree *_octree; 48 std::map<int, MVertex *> _vertexMap; 49 50 std::map<MVertex *, double> vals; 51 std::map<MVertex *, SVector3> grads; 52 std::map<MVertex *, SMetric3> hessians; 53 54 public: 55 typedef std::map<MVertex *, SMetric3> nodalMetricTensor; 56 typedef std::map<MVertex *, double> nodalField; 57 58 private: 59 nodalMetricTensor _nodalMetrics; 60 nodalField _nodalSizes, _detMetric; 61 62 std::map<int, nodalMetricTensor> setOfMetrics; 63 std::map<int, nodalField> setOfSizes; 64 std::map<int, bool> setOfRecomputeBoolean; 65 std::map<int, simpleFunction<double> *> setOfFcts; 66 std::map<int, std::vector<double> > setOfParameters; 67 std::map<int, int> setOfTechniques; 68 // std::map<int,nodalField> setOfDetMetric; 69 70 public: 71 meshMetric(std::vector<MElement *> elements); 72 meshMetric(GModel *gm); 73 74 ~meshMetric(); 75 76 // compute a new metric and add it to the set of metrics 77 // parameters[1] = lcmin (default : in global gmsh options 78 // CTX::instance()->mesh.lcMin) parameters[2] = lcmax (default : in global 79 // gmsh options CTX::instance()->mesh.lcMax) Available algorithms 80 // ("techniques"): 1: fct is a LS, metric based on Coupez technique 81 // parameters[0] = thickness of the interface (mandatory) 82 // 2: metric based on the hessian of fct 83 // parameters[0] = the final number of elements 84 // 3: fct is a LS, variant of (1) based on Frey technique (combines Coupez and 85 // curvature) 86 // parameters[0] = thickness of the interface (mandatory) 87 // parameters[3] = the required minimum number of elements to represent a 88 // circle - used for curvature-based metric (default: = 15) 89 // 4: fct is a LS, variant of (3), metric computed in LS eigendirections 90 // parameters[0] = thickness of the interface in the positive ls direction 91 // (mandatory) parameters[4] = thickness of the interface in the negative 92 // ls direction (default: =parameters[0] if not specified) parameters[3] = 93 // the required minimum number of elements to represent a circle - used for 94 // curvature-based metric (default: = 15) 95 // 5: same as 4, except that the transition in band E uses linear 96 // interpolation of h, instead of linear interpolation of metric 6: fct is a 97 // LS, metric is isotropic with linear interpolation of h in band E 7: metric 98 // based on the Hessian of fct, scaled so that the smallest element has size 99 // lcmin 100 void addMetric(int technique, simpleFunction<double> *fct, 101 const std::vector<double> ¶meters); 102 metricAtVertex(MVertex * v)103 inline SMetric3 metricAtVertex(MVertex *v) 104 { 105 if(needMetricUpdate) updateMetrics(); 106 return _nodalMetrics[v]; 107 } 108 // this function scales the mesh metric in order 109 // to reach a target number of elements 110 void scaleMetric(int nbElementsTarget, nodalMetricTensor &nmt); 111 112 void computeMetric(int metricNumber); 113 void computeMetricLevelSet(MVertex *ver, SMetric3 &hessian, SMetric3 &metric, 114 double &size, double x = 0.0, double y = 0.0, 115 double z = 0.0); 116 void computeMetricHessian(MVertex *ver, SMetric3 &hessian, SMetric3 &metric, 117 double &size, double x = 0.0, double y = 0.0, 118 double z = 0.0); 119 void computeMetricFrey(MVertex *ver, SMetric3 &hessian, SMetric3 &metric, 120 double &size, double x = 0.0, double y = 0.0, 121 double z = 0.0); 122 void computeMetricEigenDir(MVertex *ver, SMetric3 &hessian, SMetric3 &metric, 123 double &size, double x = 0.0, double y = 0.0, 124 double z = 0.0); 125 void computeMetricIsoLinInterp(MVertex *ver, SMetric3 &hessian, 126 SMetric3 &metric, double &size, double x = 0.0, 127 double y = 0.0, double z = 0.0); 128 129 void computeValues(); 130 void computeHessian(); 131 132 double getLaplacian(MVertex *v); 133 SVector3 getGradient(MVertex *v); isotropic()134 virtual bool isotropic() const { return false; } getName()135 virtual const char *getName() { return "metricField"; } getDescription()136 virtual std::string getDescription() 137 { 138 return "Anisotropic size field based on hessian of a given function"; 139 } 140 141 // get metric at point(x,y,z) (previously computes intersection of metrics if 142 // not done yet) 143 virtual double operator()(double x, double y, double z, 144 GEntity *ge = nullptr); 145 virtual void operator()(double x, double y, double z, SMetric3 &metr, 146 GEntity *ge = nullptr); 147 148 // export pos files of fct, fct gradients (fct is the lattest fct passed to 149 // meshMetric !!) and resulting metric (intersection of all computed metrics) 150 void exportInfo(const char *fileendname); 151 }; 152 153 #endif 154