diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp
index fee79d26615752dc0e9970a8e27fba05dfec5947..fc5571630baa3d21e300f64d5481aeebd3cc74cf 100644
--- a/Geo/STensor3.cpp
+++ b/Geo/STensor3.cpp
@@ -48,6 +48,27 @@ SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
return iv;
}
+SMetric3 intersection_alauzet (const SMetric3 &m1, const SMetric3 &m2)
+{
+ SMetric3 im1 = m1.invert();
+ fullMatrix<double> V(3,3);
+ fullVector<double> S(3);
+ im1 *= m2;
+ im1.eig(V,S,true);
+ SVector3 v0(V(0,0),V(1,0),V(2,0));
+ SVector3 v1(V(0,1),V(1,1),V(2,1));
+ SVector3 v2(V(0,2),V(1,2),V(2,2));
+ // is this required??
+ v0.normalize();
+ v1.normalize();
+ v2.normalize();
+ double l0 = std::max(dot(v0,m1,v0),dot(v0,m2,v0));
+ double l1 = std::max(dot(v1,m1,v1),dot(v1,m2,v1));
+ double l2 = std::max(dot(v2,m1,v2),dot(v2,m2,v2));
+ SMetric3 iv(l0,l1,l2,v0,v1,v2);
+ return iv;
+}
+
// preserve orientation of m1 !!!
SMetric3 intersection_conserveM1 (const SMetric3 &m1, const SMetric3 &m2)
{
diff --git a/Geo/STensor3.h b/Geo/STensor3.h
index 2be7571fa02879748f52a74f66ec925263eb5a9d..220fde14dfa4a636cd1353e3f61a1a5f11fc1a1c 100644
--- a/Geo/STensor3.h
+++ b/Geo/STensor3.h
@@ -173,6 +173,8 @@ SMetric3 intersection_conserve_mostaniso_2d (const SMetric3 &m1, const SMetric3
// compute the largest inscribed ellipsoid...
SMetric3 intersection (const SMetric3 &m1,
const SMetric3 &m2);
+SMetric3 intersection_alauzet (const SMetric3 &m1,
+ const SMetric3 &m2);
SMetric3 interpolation (const SMetric3 &m1,
const SMetric3 &m2,
const double t);
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index dccd86205228e47a713b294d210cd0cc708977b8..988ed8048d05dc2ad72c6c9317e71e4ce3cda966 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -485,6 +485,55 @@ class CylinderField : public Field
}
};
+class SphereField : public Field
+{
+ double v_in, v_out;
+ double xc,yc,zc;
+ double R;
+
+ public:
+ std::string getDescription()
+ {
+ return "The value of this field is VIn inside a sphere, VOut outside. "
+ "The sphere is given by\n\n"
+ " ||dX||^2 < R^2 &&\n"
+ " dX = (X - XC)^2 + (Y-YC)^2 + (Z-ZC)^2";
+ }
+ SphereField()
+ {
+ v_in = v_out = xc = yc = zc = R = 0;
+
+ options["VIn"] = new FieldOptionDouble
+ (v_in, "Value inside the sphere");
+ options["VOut"] = new FieldOptionDouble
+ (v_out, "Value outside the sphere");
+
+ options["XCenter"] = new FieldOptionDouble
+ (xc, "X coordinate of the sphere center");
+ options["YCenter"] = new FieldOptionDouble
+ (yc, "Y coordinate of the sphere center");
+ options["ZCenter"] = new FieldOptionDouble
+ (zc, "Z coordinate of the sphere center");
+
+
+ options["Radius"] = new FieldOptionDouble
+ (R,"Radius");
+
+ }
+ const char *getName()
+ {
+ return "Sphere";
+ }
+ double operator() (double x, double y, double z, GEntity *ge=0)
+ {
+ double dx = x-xc;
+ double dy = y-yc;
+ double dz = z-zc;
+
+ return ( (dx*dx + dy*dy + dz*dz < R*R) ) ? v_in : v_out;
+ }
+};
+
class FrustumField : public Field
{
double x1,y1,z1;
@@ -1026,15 +1075,17 @@ class MathEvalField : public Field
std::string f;
public:
- void myAction(){
- printf("doing sthg \n");
+ void myAction()
+ {
+ printf("doing sthg\n");
}
MathEvalField()
{
options["F"] = new FieldOptionString
(f, "Mathematical function to evaluate.", &update_needed);
f = "F2 + Sin(z)";
- callbacks["test"] = new FieldCallbackGeneric<MathEvalField>(this, &MathEvalField::myAction, "description blabla");
+ callbacks["test"] = new FieldCallbackGeneric<MathEvalField>
+ (this, &MathEvalField::myAction, "description blabla");
}
double operator() (double x, double y, double z, GEntity *ge=0)
{
@@ -1337,6 +1388,76 @@ class MinAnisoField : public Field
}
};
+class IntersectAnisoField : public Field
+{
+ std::list<int> idlist;
+ public:
+ IntersectAnisoField()
+ {
+ options["FieldsList"] = new FieldOptionList
+ (idlist, "Field indices", &update_needed);
+ }
+ virtual bool isotropic () const {return false;}
+ std::string getDescription()
+ {
+ return "Take the intersection of 2 anisotropic fields according to Alauzet.";
+ }
+ virtual void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0)
+ {
+ // check if idlist contains 2 elements other error message
+ SMetric3 v;
+ for(std::list<int>::iterator it = idlist.begin(); it != idlist.end(); it++) {
+ Field *f = (GModel::current()->getFields()->get(*it));
+ SMetric3 ff;
+ if(f && *it != id) {
+ if (f->isotropic()){
+ double l = (*f) (x, y, z, ge);
+ ff = SMetric3(1./(l*l));
+ }
+ else{
+ (*f) (x, y, z, ff, ge);
+ }
+ if (it == idlist.begin())
+ v = ff;
+ else
+ v = intersection_alauzet(v,ff);
+ }
+ }
+ metr = v;
+ }
+ double operator() (double x, double y, double z, GEntity *ge=0)
+ {
+ // check if idlist contains 2 elements other error message
+ SMetric3 metr;
+ for(std::list<int>::iterator it = idlist.begin(); it != idlist.end(); it++) {
+ Field *f = (GModel::current()->getFields()->get(*it));
+ SMetric3 m;
+ if(f && *it != id){
+ if (!f->isotropic()){
+ (*f)(x, y, z, m, ge);
+ }
+ else {
+ double L = (*f)(x, y, z, ge);
+ for (int i = 0; i < 3; i++)
+ m(i,i) = 1. / (L*L);
+ }
+ }
+ if (it == idlist.begin())
+ metr = m;
+ else
+ metr = intersection_alauzet(metr,m);
+ }
+ fullMatrix<double> V(3,3);
+ fullVector<double> S(3);
+ metr.eig(V, S, 1);
+ return sqrt(1./S(2)); //S(2) is largest eigenvalue
+ }
+ const char *getName()
+ {
+ return "IntersectAniso";
+ }
+};
+
class MinField : public Field
{
std::list<int> idlist;
@@ -1355,7 +1476,18 @@ class MinField : public Field
double v = MAX_LC;
for(std::list<int>::iterator it = idlist.begin(); it != idlist.end(); it++) {
Field *f = (GModel::current()->getFields()->get(*it));
- if(f && *it != id) v = std::min(v, (*f) (x, y, z, ge));
+ if(f && *it != id) {
+ if (f->isotropic())
+ v = std::min(v, (*f) (x, y, z, ge));
+ else{
+ SMetric3 ff;
+ (*f) (x, y, z, ff, ge);
+ fullMatrix<double> V(3,3);
+ fullVector<double> S(3);
+ ff.eig(V, S, 1);
+ v = std::min(v, sqrt(1./S(2))); //S(2) is largest eigenvalue
+ }
+ }
}
return v;
}
@@ -1383,7 +1515,18 @@ class MaxField : public Field
double v = -MAX_LC;
for(std::list<int>::iterator it = idlist.begin(); it != idlist.end(); it++) {
Field *f = (GModel::current()->getFields()->get(*it));
- if(f && *it != id) v = std::max(v, (*f) (x, y, z, ge));
+ if(f && *it != id) {
+ if (f->isotropic())
+ v = std::max(v, (*f) (x, y, z, ge));
+ else{
+ SMetric3 ff;
+ (*f) (x, y, z, ff, ge);
+ fullMatrix<double> V(3,3);
+ fullVector<double> S(3);
+ ff.eig(V, S, 1);
+ v = std::max(v, sqrt(1./S(0))); //S(0) is smallest eigenvalue
+ }
+ }
}
return v;
}
@@ -2157,6 +2300,7 @@ FieldManager::FieldManager()
#endif
map_type_name["Box"] = new FieldFactoryT<BoxField>();
map_type_name["Cylinder"] = new FieldFactoryT<CylinderField>();
+ map_type_name["Sphere"] = new FieldFactoryT<SphereField>();
map_type_name["Frustum"] = new FieldFactoryT<FrustumField>();
map_type_name["LonLat"] = new FieldFactoryT<LonLatField>();
#if defined(HAVE_POST)
@@ -2166,6 +2310,7 @@ FieldManager::FieldManager()
map_type_name["Restrict"] = new FieldFactoryT<RestrictField>();
map_type_name["Min"] = new FieldFactoryT<MinField>();
map_type_name["MinAniso"] = new FieldFactoryT<MinAnisoField>();
+ map_type_name["IntersectAniso"] = new FieldFactoryT<IntersectAnisoField>();
map_type_name["Max"] = new FieldFactoryT<MaxField>();
map_type_name["Laplacian"] = new FieldFactoryT<LaplacianField>();
map_type_name["Mean"] = new FieldFactoryT<MeanField>();