[Getdp] RectShell xform patch

[Martin Devera]

Hello,

I started to learn about Maxwell FEM using GetDP and while trying
to implement VolAxiRectShell I found inconvenient to split the geometry
of shell to X/Y axes.
I patched code to autodetect rectangular shell's axis using Axis=0 code:

Jacobian VolAxiRectShell{Val_Rint, Val_Rext,0}

the patch follows.
Also when already here, I tested in on a simple magnetostatic problem
(to determine strenght of B field of aircored helical cone coil) and
found some numerical inconsistencies.
When you see coil.jpg (B strength as Norm[{d a}]) you can clearly see
random colored triangles near symmetry axis and in bottom shell space.
Any clue ?

best regards, Martin

PS: the patch contains "fix", change of 1e-2 to 1e-12 - seemed as typo

--- Main/Get_Geometry.c.old     2008-09-27 09:05:11.000000000 +0200
+++ Main/Get_Geometry.c 2008-09-27 09:05:17.000000000 +0200
@@ -444,6 +444,13 @@ double  Transformation (int Dim, int Typ
      dRdz = (Z-Cz)/R ;
    }
    else{
+    if (!Axis) {
+      /* Autodetect axis to allow rectangular slab */
+      double aB = fabs(B) + 1.e-12*fabs(B), aA = fabs(A) - 1.e-12*fabs(A);
+      if (fabs(X-Cx) < aB && fabs(X-Cx) > aA) Axis = 1;
+      else if (fabs(Y-Cy) < aB && fabs(Y-Cy) > aA) Axis = 2;
+      else if (fabs(Z-Cz) < aB && fabs(Z-Cz) > aA) Axis = 3;
+    }
      switch(Axis) {
      case 1: R = (X-Cx) ; dRdx =1.0 ; dRdy =0.0 ; dRdz =0.0 ; break;
      case 2: R = (Y-Cy) ; dRdx =0.0 ; dRdy =1.0 ; dRdz =0.0 ; break;
@@ -453,7 +460,7 @@ double  Transformation (int Dim, int Typ
    }

    if ( (fabs(R) > fabs(B) + 1.e-12*fabs(B)) ||
-       (fabs(R) < fabs(A) - 1.e-2*fabs(A)) )
+       (fabs(R) < fabs(A) - 1.e-12*fabs(A)) )
      Msg(GERROR, "Bad parameters for transformation Jacobian: %g not in [%g,%g]", R, A, B) ;

    if (B == R) {