[Getdp] RectShell xform patch
Martin Devera
devik at cdi.cz
Sat Sep 27 09:20:36 CEST 2008
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) {