Skip to content
Snippets Groups Projects
Commit 99f025f0 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

add StopAtDistMax for Treshold field

tmp fix build
parent f27147bd
No related branches found
No related tags found
No related merge requests found
......@@ -451,7 +451,7 @@ class ThresholdField : public Field
{
int iField;
double dmin, dmax, lcmin, lcmax;
bool sigmoid;
bool sigmoid, stopAtDistMax;
public:
const char *get_name()
{
......@@ -465,6 +465,7 @@ class ThresholdField : public Field
lcmin = 0.1;
lcmax = 1;
sigmoid = false;
stopAtDistMax = false;
options["IField"] = new FieldOptionInt(iField, "Field index");
options["DistMin"] = new FieldOptionDouble(dmin, "Distance from entity up to which "
"element size will be LcMin");
......@@ -475,6 +476,8 @@ class ThresholdField : public Field
options["Sigmoid"] = new FieldOptionBool(sigmoid, "True to interpolate between LcMin "
"and LcMax using a sigmoid, false to "
"interpolate linearly");
options["StopAtDistMax"] = new FieldOptionBool(stopAtDistMax, "True to not impose element "
"size outside DistMax");
}
double operator() (double x, double y, double z)
{
......@@ -482,7 +485,10 @@ class ThresholdField : public Field
double r = ((*field) (x, y, z) - dmin) / (dmax - dmin);
r = std::max(std::min(r, 1.), 0.);
double lc;
if(sigmoid){
if(stopAtDistMax && r >= 1.){
lc = MAX_LC;
}
else if(sigmoid){
double s = exp(12. * r - 6.) / (1. + exp(12. * r - 6.));
lc = lcmin * (1. - s) + lcmax * s;
}
......
......@@ -861,21 +861,23 @@ void getRegionVertices(GRegion *gr,
const double t2 = points(k, 1);
const double t3 = points(k, 2);
SPoint3 pos;
incomplete->pnt(t1,t2,t3,pos);
// FIXME: KOEN - I had to comment this out (MElement does not have
// pnt() member) -- CG
v = new MVertex(pos.x(),pos.y(),pos.z(),gr);
// SPoint3 pos;
// incomplete->pnt(t1,t2,t3,pos);
// v = new MVertex(pos.x(),pos.y(),pos.z(),gr);
// double X(0),Y(0),Z(0);
// for (int j=0; j<incomplete->getNumVertices(); j++){
// double sf ; incomplete->getShapeFunction(j,t1,t2,t3,sf);
// MVertex *vt = incomplete->getVertex(j);
// X += sf * vt->x();
// Y += sf * vt->y();
// Z += sf * vt->z();
// }
//
// v = new MVertex(X,Y,Z, gr);
double X(0),Y(0),Z(0);
for (int j=0; j<incomplete->getNumVertices(); j++){
double sf ; incomplete->getShapeFunction(j,t1,t2,t3,sf);
MVertex *vt = incomplete->getVertex(j);
X += sf * vt->x();
Y += sf * vt->y();
Z += sf * vt->z();
}
v = new MVertex(X,Y,Z, gr);
gr->mesh_vertices.push_back(v);
vr.push_back(v);
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment