462
460
options["ZMin"] = new FieldOptionDouble
463
461
(z_min, "Minimum Z coordinate of the box");
464
462
options["ZMax"] = new FieldOptionDouble
465
(z_max, "Maximum Z coordinate of the box");
463
(z_max, "Maximum Z coordinate of the box");
467
465
const char *getName()
471
469
double operator() (double x, double y, double z, GEntity *ge=0)
473
471
return (x <= x_max && x >= x_min && y <= y_max && y >= y_min && z <= z_max
474
472
&& z >= z_min) ? v_in : v_out;
476
class CylinderField : public Field
484
std::string getDescription()
486
return "The value of this field is VIn inside a frustrated cylinder, VOut outside.\n"
487
"The cylinder is given by\n\n"
488
" ||dX||^2 < R^2 &&\n"
489
" (X-X0).A < ||A||^2\n"
490
" dX = (X - X0) - ((X - X0).A)/(||A||^2) . A";
494
v_in = v_out = xc = yc = zc = xa = ya = R = 0;
497
options["VIn"] = new FieldOptionDouble
498
(v_in, "Value inside the cylinder");
499
options["VOut"] = new FieldOptionDouble
500
(v_out, "Value outside the cylinder");
502
options["XCenter"] = new FieldOptionDouble
503
(xc, "X coordinate of the cylinder center");
504
options["YCenter"] = new FieldOptionDouble
505
(yc, "Y coordinate of the cylinder center");
506
options["ZCenter"] = new FieldOptionDouble
507
(zc, "Z coordinate of the cylinder center");
510
options["XAxis"] = new FieldOptionDouble
511
(xa, "X component of the cylinder axis");
512
options["YAxis"] = new FieldOptionDouble
513
(ya, "Y component of the cylinder axis");
514
options["ZAxis"] = new FieldOptionDouble
515
(za, "Z component of the cylinder axis");
517
options["Radius"] = new FieldOptionDouble
521
const char *getName()
525
double operator() (double x, double y, double z, GEntity *ge=0)
531
double adx = (xa * dx + ya * dy + za * dz)/(xa*xa + ya*ya + za*za);
537
return ((dx*dx + dy*dy + dz*dz < R*R) && fabs(adx) < 1) ? v_in : v_out;
478
542
class ThresholdField : public Field
481
546
double dmin, dmax, lcmin, lcmax;
482
547
bool sigmoid, stopAtDistMax;
484
const char *getName()
549
virtual const char *getName()
486
551
return "Threshold";
488
std::string getDescription()
553
virtual std::string getDescription()
490
555
return "F = LCMin if Field[IField] <= DistMin,\n"
491
556
"F = LCMax if Field[IField] >= DistMax,\n"
606
class BoundaryLayerField : public ThresholdField {
610
virtual bool isotropic () const {return false;}
611
virtual const char *getName()
613
return "BoundaryLayer";
615
virtual std::string getDescription()
617
return "F = LCMin if Field[IField] <= DistMin,\n"
618
"F = LCMax if Field[IField] >= DistMax,\n"
619
"F = interpolation between LcMin and LcMax if DistMin < Field[IField] < DistMax";
621
virtual void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0)
623
Field *field = GModel::current()->getFields()->get(iField);
625
metr(0,0) = 1/(MAX_LC*MAX_LC);
626
metr(1,1) = 1/(MAX_LC*MAX_LC);
627
metr(2,2) = 1/(MAX_LC*MAX_LC);
628
metr(0,1) = metr(0,2) = metr(1,2) = 0;
631
double dist = (*field) (x, y, z);
633
double r = (dist - dmin) / (dmax - dmin);
634
r = std::max(std::min(r, 1.), 0.);
636
if(stopAtDistMax && r >= 1.){
640
double s = exp(12. * r - 6.) / (1. + exp(12. * r - 6.));
641
lc = lcmin * (1. - s) + lcmax * s;
644
lc = lcmin * (1 - r) + lcmax * r;
647
double delta = std::min(CTX::instance()->lc / 1e4, dist);
649
((*field) (x + delta / 2, y, z) -
650
(*field) (x - delta / 2, y, z)) / delta;
652
((*field) (x, y + delta / 2, z) -
653
(*field) (x, y - delta / 2, z)) / delta;
655
((*field) (x, y, z + delta / 2) -
656
(*field) (x, y, z - delta / 2)) / delta;
658
SVector3 g(gx,gy,gz);
662
if (fabs(gx) < fabs(gy) && fabs(gx) < fabs(gz))
663
t1 = SVector3(1,0,0);
664
else if (fabs(gy) < fabs(gx) && fabs(gy) < fabs(gz))
665
t1 = SVector3(0,1,0);
667
t1 = SVector3(0,0,1);
669
t2 = crossprod(g,t1);
671
t1 = crossprod(t2,g);
673
metr = SMetric3(1./(lc*lc),
541
680
class GradientField : public Field
543
682
int iField, kind;
730
870
Field *field = GModel::current()->getFields()->get(iField);
731
871
if(!field) return MAX_LC;
732
872
return ((*field) (x + delta , y, z)+ (*field) (x - delta , y, z)
733
+(*field) (x, y + delta , z)+ (*field) (x, y - delta , z)
734
+(*field) (x, y, z + delta )+ (*field) (x, y, z - delta )
735
-6* (*field) (x , y, z)) / (delta*delta);
873
+(*field) (x, y + delta , z)+ (*field) (x, y - delta , z)
874
+(*field) (x, y, z + delta )+ (*field) (x, y, z - delta )
875
-6* (*field) (x , y, z)) / (delta*delta);
766
906
Field *field = GModel::current()->getFields()->get(iField);
767
907
if(!field) return MAX_LC;
768
908
return ((*field) (x + delta , y, z) + (*field) (x - delta, y, z)
769
+ (*field) (x, y + delta, z) + (*field) (x, y - delta, z)
770
+ (*field) (x, y, z + delta) + (*field) (x, y, z - delta)
771
+ (*field) (x, y, z)) / 7;
909
+ (*field) (x, y + delta, z) + (*field) (x, y - delta, z)
910
+ (*field) (x, y, z + delta) + (*field) (x, y, z - delta)
911
+ (*field) (x, y, z)) / 7;
775
#if defined(HAVE_MATH_EVAL)
776
915
class MathEvalExpression
779
std::list<Field*> *list;
785
std::string function;
786
char *c_str_function;
919
std::set<int> _fields;
788
double evaluate(double x, double y, double z)
921
MathEvalExpression() : _f(0) {}
922
~MathEvalExpression(){ if(_f) delete _f; }
923
bool set_function(const std::string &f)
792
for(int i = 0; i < nvalues; i++){
793
switch (evaluators_id[i]) {
805
Field *f = GModel::current()->getFields()->get(evaluators_id[i]);
806
values[i] = f ? (*f) (x, y, z) : MAX_LC;
925
// get id numbers of fields appearing in the function
932
while(i + 1 + j < f.size() && f[i + 1 + j] >= '0' && f[i + 1 + j] <= '9'){
936
_fields.insert(atoi(id.c_str()));
810
return evaluator_evaluate(eval, nvalues, names, values);
819
bool set_function(const std::string & f)
822
error_status = false;
823
c_str_function = strdup(f.c_str());
824
eval = evaluator_create(c_str_function);
940
std::vector<std::string> expressions(1), variables(3 + _fields.size());
946
for(std::set<int>::iterator it = _fields.begin(); it != _fields.end(); it++){
947
std::ostringstream sstream;
948
sstream << "F" << *it;
949
variables[i++] = sstream.str();
952
_f = new mathEvaluator(expressions, variables);
953
if(expressions.empty()) {
829
evaluator_get_variables(eval, &names, &nvalues);
830
values = new double[nvalues];
831
evaluators_id = new int[nvalues];
832
for(int i = 0; i < nvalues; i++) {
834
if(!strcmp("x", names[i]))
835
evaluators_id[i] = -1;
836
else if(!strcmp("y", names[i]))
837
evaluators_id[i] = -2;
838
else if(!strcmp("z", names[i]))
839
evaluators_id[i] = -3;
840
else if(sscanf(names[i], "F%i", &id) == 1)
841
evaluators_id[i] = id;
843
Msg::Error("Unknown matheval argument \"%s\"\n", names[i]);
853
free(c_str_function);
855
evaluator_destroy(eval);
859
delete evaluators_id;
861
~MathEvalExpression()
960
double evaluate(double x, double y, double z)
962
if(!_f) return MAX_LC;
963
std::vector<double> values(3 + _fields.size()), res(1);
968
for(std::set<int>::iterator it = _fields.begin(); it != _fields.end(); it++){
969
Field *field = GModel::current()->getFields()->get(*it);
970
values[i++] = field ? (*field)(x, y, z) : MAX_LC;
972
if(_f->eval(values, res))
928
1040
for(int i = 0; i < 3; i++) {
929
1041
if(!expr[i].set_function(f[i]))
930
1042
Msg::Error("Field %i : Invalid matheval expression \"%s\"",
931
this->id, f[i].c_str());
1043
this->id, f[i].c_str());
933
1045
update_needed = false;
935
1047
Field *field = GModel::current()->getFields()->get(ifield);
936
1048
if(!field) return MAX_LC;
937
1049
return (*field)(expr[0].evaluate(x, y, z),
938
expr[1].evaluate(x, y, z),
939
expr[2].evaluate(x, y, z));
1050
expr[1].evaluate(x, y, z),
1051
expr[2].evaluate(x, y, z));
941
1053
const char *getName()
948
1059
#if !defined(HAVE_NO_POST)
949
1060
class PostViewField : public Field
1198
1309
Surface *s = FindSurface(*it);
1200
1311
for(int i = 0; i < n_nodes_by_edge; i++) {
1201
for(int j = 0; j < n_nodes_by_edge; j++) {
1202
double u = (double)i / (n_nodes_by_edge - 1);
1203
double v = (double)j / (n_nodes_by_edge - 1);
1204
Vertex V = InterpolateSurface(s, u, v, 0, 0);
1205
zeronodes[k][0] = V.Pos.X;
1206
zeronodes[k][1] = V.Pos.Y;
1207
zeronodes[k++][2] = V.Pos.Z;
1312
for(int j = 0; j < n_nodes_by_edge; j++) {
1313
double u = (double)i / (n_nodes_by_edge - 1);
1314
double v = (double)j / (n_nodes_by_edge - 1);
1315
Vertex V = InterpolateSurface(s, u, v, 0, 0);
1316
zeronodes[k][0] = V.Pos.X;
1317
zeronodes[k][1] = V.Pos.Y;
1318
zeronodes[k++][2] = V.Pos.Z;
1212
1323
GFace *f = GModel::current()->getFaceByTag(*it);
1214
1325
for(int i = 0; i < n_nodes_by_edge; i++) {
1215
for(int j = 0; j < n_nodes_by_edge; j++) {
1216
double u = (double)i / (n_nodes_by_edge - 1);
1217
double v = (double)j / (n_nodes_by_edge - 1);
1218
Range<double> b1 = ge->parBounds(0);
1219
Range<double> b2 = ge->parBounds(1);
1220
double t1 = b1.low() + u * (b1.high() - b1.low());
1221
double t2 = b2.low() + v * (b2.high() - b2.low());
1222
GPoint gp = f->point(t1, t2);
1223
zeronodes[k][0] = gp.x();
1224
zeronodes[k][1] = gp.y();
1225
zeronodes[k++][2] = gp.z();
1326
for(int j = 0; j < n_nodes_by_edge; j++) {
1327
double u = (double)i / (n_nodes_by_edge - 1);
1328
double v = (double)j / (n_nodes_by_edge - 1);
1329
Range<double> b1 = ge->parBounds(0);
1330
Range<double> b2 = ge->parBounds(1);
1331
double t1 = b1.low() + u * (b1.high() - b1.low());
1332
double t2 = b2.low() + v * (b2.high() - b2.low());
1333
GPoint gp = f->point(t1, t2);
1334
zeronodes[k][0] = gp.x();
1335
zeronodes[k][1] = gp.y();
1336
zeronodes[k++][2] = gp.z();
1309
1421
for(int ele = 0; ele < data->getNumElements(0, ent); ele++){
1310
1422
if(data->skipElement(0, ent, ele)) continue;
1311
1423
for(int nod = 0; nod < data->getNumNodes(0, ent, ele); nod++){
1313
data->getNode(0, ent, ele, nod, x, y, z);
1425
data->getNode(0, ent, ele, nod, x, y, z);
1314
1426
double val = (*this)(x, y, z);
1315
for(int comp = 0; comp < data->getNumComponents(0, ent, ele); comp++)
1316
data->setValue(0, ent, ele, nod, comp, val);
1427
for(int comp = 0; comp < data->getNumComponents(0, ent, ele); comp++)
1428
data->setValue(0, ent, ele, nod, comp, val);