[Getdp] OneTooth Example with non linear material

Matthias Lang lang_matthias at gmx.de
Tue Oct 22 23:55:37 CEST 2013


Hello ererybody,

I took the 3D magnetostatic example OneTooth.pro
and included nonlinear material. (found in Material.pro)
But unfortunatly I get the following Error Message:

.
.
Info    : GetDP - E n d   P r e - P r o c e s s i n g
Info    : GetDP - P r o c e s s i n g . . .
Info    : GetDP - IterativeLoop ...
Info    : GetDP - Non Linear Iteration Relaxation 1
Info    : GetDP - GenerateJac[Sys_Mag]
Error   : GetDP - Function _Exp_Fct_17 called with too few arguments
Info    : Done running 'GetDP'

Where is the problem? - inside or in front of the computer?

Regards

Matthias Lang
-------------- next part --------------
/* -------------------------------------------------------------------
    File "MagSta_a_3D_nl.pro"
    Magnetostatics - Magnetic vector potential a formulation (3D)
    -------------------------------------------------------------------
    I N P U T
    ---------
    GlobalGroup : (Extension Mag is for Magnetic problem)
    -----------
    Domain_Mag Whole magnetic domain
    DomainS_Mag Inductor regions (Source)
    Function :
    --------
    nu[] Magnetic reluctivity
    Constraint :
    ----------
*/
Group {
    DefineGroup[ Domain_Mag, DomainS_Mag, Domain_NL];
}


Function {
    DefineFunction[ nu, dhdb_NL ];
    DefineConstant[
	Nb_max_iter        = 30,
	relaxation_factor  = 1,
	stop_criterion     = 1e-5,
	NbT        	   = 1,
	Steel_1020_dhdb_NL
    ];
}



FunctionSpace {
    // Magnetic vector potential a (b = curl a)
    { Name Hcurl_a_Mag_3D; Type Form1;
      BasisFunction {
        // a = a s
        //  e   e
        { Name se; NameOfCoef ae; Function BF_Edge;
          Support Domain_Mag; Entity EdgesOf[ All ]; }
    }
    Constraint {
        { NameOfCoef ae; EntityType EdgesOf;
            NameOfConstraint MagneticVectorPotential_3D; }
        { NameOfCoef ae; EntityType EdgesOfTreeIn; EntitySubType StartingOn; 
            NameOfConstraint Gauge; }
        }
    }

}

Formulation {
    { Name Magnetostatics_a_3D; Type FemEquation;
      Quantity {
        { Name a ; Type Local; NameOfSpace Hcurl_a_Mag_3D; }
      }
      Equation {
            Galerkin { [ nu[] * Dof{Curl a} , {Curl a} ]; In Domain_Mag;
                        Jacobian Vol; Integration CurlCurl; }
        Galerkin { JacNL [ dhdb_NL[{d a}] * Dof{d a} , {d a} ] ;
        		In Domain_NL ; Jacobian Vol ; Integration CurlCurl ; }

    Galerkin { [ - js[] , {a} ]; In DomainS_Mag;
                        Jacobian Vol; Integration CurlCurl; }
	
      }
    }
}

Resolution {
    {   Name MagSta_a_3D;
        System {
            { Name Sys_Mag; NameOfFormulation Magnetostatics_a_3D; }
        }
        Operation {
//             Generate[Sys_Mag]; Solve[Sys_Mag]; //linear case
	IterativeLoop[Nb_max_iter, stop_criterion, relaxation_factor]{
            GenerateJac[Sys_Mag] ; SolveJac[Sys_Mag] ; }	
	SaveSolution[Sys_Mag];
        }
    }
}

PostProcessing {
    {   Name MagSta_a_3D; NameOfFormulation Magnetostatics_a_3D;
        Quantity {
            { Name a;
                Value {
                    Local { [ {a} ]; In Domain_Mag; Jacobian Vol; }
                }
            }
			{ Name a_norm ; 
				Value { 
					Local { [ Norm[{a}] ] ; In Domain_Mag ; Jacobian Vol ; }
					} 
			} 
            { Name b;
                Value {
                    Local { [ {d a} ]; In Domain_Mag; Jacobian Vol; }
                }
            }
            { Name h;
                Value {
                    Local { [ nu[] * {d a} ]; In Domain_Mag; Jacobian Vol; }
                }
            }
            { Name js;
                Value {
                    Local { [js[]] ; In DomainS_Mag; Jacobian Vol; }
                }
            }
                
        }
    }
}
-------------- next part --------------
Function{
/*--------------------------------------------------------------------------------------*/
  // nu = 100. + 10. * exp ( 1.8 * b * b )
  // analytical
  nu_1a[] = 100. + 10. * Exp[1.8*SquNorm[$1]] ;
  dnudb2_1a[] = 18. * Exp[1.8*SquNorm[$1]] ;
  h_1a[] = nu_1a[$1]*$1 ;
  dhdb_1a[] = TensorDiag[1,1,1] * nu_1a[$1#1] + 2*dnudb2_1a[#1] * SquDyadicProduct[#1]  ;
  dhdb_1a_NL[] = 2*dnudb2_1a[$1#1] * SquDyadicProduct[#1]  ;
/*--------------------------------------------------------------------------------------*/

 
/*--------------------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------------------*/

  // BH data EIcore

  MatEIcore_h = {0.,
    2.1827e+02,3.8272e+02,5.2036e+02,7.1167e+02,
    8.4921e+02,1.0405e+03,1.4741e+03,2.0409e+03,
    3.0109e+03,5.0572e+03,7.5335e+03,1.0037e+04,
    1.2486e+04,1.5015e+04,1.7464e+04,2.0021e+04,
    2.2040e+04,2.5000e+04} ;

  MatEIcore_b = {0.,
    2.0329e-01,4.0287e-01,6.0986e-01,8.0575e-01,
    1.0053e+00,1.1975e+00,1.4008e+00,1.5154e+00,
    1.5967e+00,1.6706e+00,1.7076e+00,1.7335e+00,
    1.7446e+00,1.7593e+00,1.7667e+00,1.7815e+00,
    1.7889e+00,1.7963e+00} ;

  MatEIcore_b2 = List[MatEIcore_b]^2 ;
  MatEIcore_nu = List[MatEIcore_h]/List[MatEIcore_b] ;
  MatEIcore_nu(0) = MatEIcore_nu(1);

  MatEIcore_nu_b2  = ListAlt[MatEIcore_b2, MatEIcore_nu] ;
  //nu_EIcore[] = InterpolationLinear[SquNorm[$1]]{List[MatEIcore_nu_b2]} ; // used later in *.pro file
  MatEIcore_interp_nu[] = InterpolationLinear[SquNorm[$1]]{List[MatEIcore_nu_b2]} ; // used later in *.pro file
  //dnudb2_EIcore[] = dInterpolationLinear[SquNorm[$1]]{List[MatEIcore_nu_b2]} ;
  MatEIcore_interp_dnudb2[] = dInterpolationLinear[SquNorm[$1]]{List[MatEIcore_nu_b2]} ;
//h_EIcore[] = nu_EIcore[$1] * $1 ;
  MatEIcore_interp_h[] = MatEIcore_interp_nu[$1] * $1 ;
//dhdb_EIcore[] = TensorDiag[1,1,1]*nu_EIcore[$1#1] + 2*dnudb2_EIcore[#1] * SquDyadicProduct[#1] ;
  MatEIcore_dhdb[] = TensorDiag[1,1,1]*MatEIcore_interp_nu[$1#1] + 2*MatEIcore_interp_dnudb2[#1] * SquDyadicProduct[#1] ;
//dhdb_EIcore_NL[] = 2*dnudb2_EIcore[$1] * SquDyadicProduct[$1] ; // used later in *.pro file
  MatEIcore_dhdb_NL[] = 2*MatEIcore_interp_dnudb2[$1] * SquDyadicProduct[$1] ; // used later in *.pro file
/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/

  // BH data for: 1020 Steel

  Steel_1020_h = {0.000000,
    79.577472,    100.182101,    126.121793,    158.777930,
   199.889571,    251.646061,    316.803620,    398.832128,
   502.099901,    632.106325,    795.774715,   1001.821011,
  1261.217929,   1587.779301,   1998.895710,   2516.460605,
  3168.036204,   3988.321282,   5020.999013,   6321.063250,
  7957.747155,  10018.210114,  12612.179293,  15877.793010,
 19988.957103,  25164.606052,  31680.362037,  39883.212823,
 50209.990127,  63210.632497,  79577.471546, 100182.101136,
126121.792926, 158777.930096, 199889.571030, 251646.060522,
316803.620370};

  Steel_1020_b = {0.000000,
	0.166252, 0.208725, 0.261635, 0.327146,
	0.407495, 0.504604, 0.619366, 0.750502, 
	0.893195, 1.038258, 1.173276, 1.286731, 
	1.373533, 1.437658, 1.489189, 1.538344, 
	1.591321, 1.649723, 1.711927, 1.774712, 
	1.834976, 1.891442, 1.945262, 1.998705, 
	2.052852, 2.106184, 2.155149, 2.196375, 
	2.229008, 2.255450, 2.280053, 2.307136, 
	2.339975, 2.381044, 2.432708, 2.497748, 
	2.579627}; 

	
  Steel_1020_b2 = List[Steel_1020_b]^2 ;
  Steel_1020_nu = List[Steel_1020_h]/List[Steel_1020_b] ;
  Steel_1020_nu(0) = Steel_1020_nu(1);

  Steel_1020_nu_b2  = ListAlt[Steel_1020_b2, Steel_1020_nu] ;
  //nu_EIcore[] = InterpolationLinear[SquNorm[$1]]{List[Steel_1020_nu_b2]} ; // used later in *.pro file
  Steel_1020_interp_nu[] = InterpolationLinear[SquNorm[$1]]{List[Steel_1020_nu_b2]} ; // used later in *.pro file
  //dnudb2_EIcore[] = dInterpolationLinear[SquNorm[$1]]{List[Steel_1020_nu_b2]} ;
  Steel_1020_interp_dnudb2[] = dInterpolationLinear[SquNorm[$1]]{List[Steel_1020_nu_b2]} ;
//h_EIcore[] = nu_EIcore[$1] * $1 ;
  Steel_1020_interp_h[] = Steel_1020_interp_nu[$1] * $1 ;
//dhdb_EIcore[] = TensorDiag[1,1,1]*nu_EIcore[$1#1] + 2*dnudb2_EIcore[#1] * SquDyadicProduct[#1] ;
  Steel_1020_dhdb[] = TensorDiag[1,1,1]*Steel_1020_interp_nu[$1#1] + 2*Steel_1020_interp_dnudb2[#1] * SquDyadicProduct[#1] ;
//dhdb_EIcore_NL[] = 2*dnudb2_EIcore[$1] * SquDyadicProduct[$1] ; // used later in *.pro file
  Steel_1020_dhdb_NL[] = 2*Steel_1020_interp_dnudb2[$1] * SquDyadicProduct[$1] ; // used later in *.pro file
/*--------------------------------------------------------------------------------------*/
}






-------------- next part --------------
A non-text attachment was scrubbed...
Name: OneTooth.geo
Type: application/vnd.dynageo
Size: 7369 bytes
Desc: not available
URL: <http://www.geuz.org/pipermail/getdp/attachments/20131022/186a4a61/attachment.geo>
-------------- next part --------------
// Fichier "OneTooth.pro"
// Magnétostatique

Include "OneTooth-dat.pro";
Include "Materials.pro"; // nonlinear BH caracteristic of magnetic material

Group {
    Air          = #AIR;
    Core         = #CORE;
    Ind          = #COIL;
    Plate        = #PLATE;
    
	SurfaceGe0   = #CUT_SYMMETRY ;
    SurfaceGInf  = #AIR_INF;
    SkinCore     = #SKIN_CORE;
    SkinPlate    = #SKIN_PLATE;

    DomainS_Mag = Region[{Ind}];
    Domain_Mag  = Region[{Core, Air, Ind, Plate}];
    Domain_NL  = Region[{Core, Plate}];
}

Function {
    mu0 = 4.e-7 * Pi;  
    murSteel = 100.;
    nu [ Region[{Air, Ind}] ] = 1. / mu0;
//linear Case
//    nu [ Region[{Core}]] = 1. / (murSteel*mu0);
//    nu [ Region[{Plate}]] = 1. / (murSteel*mu0);
//Nonlinear Case	
//	nu [ Domain_NL ] = MatEIcore_interp_nu[$1] ;
//	dhdb_NL [ Domain_NL ] = MatEIcore_dhdb_NL[$1];
	nu [ Domain_NL ] = Steel_1020_interp_nu[$1] ;
	dhdb_NL [ Domain_NL ] = Steel_1020_dhdb_NL[$1];
    
    Sc          = th_coil*ly_coil;  //cross section of coil
    Ns          = 100;		//Number of turns of the inductor Ind
    Val_Current = 1.0;		//current in the inductor

    //calculation of current vectors
    xc = 0.;  // center of the coil
    zc = 0.;   // center of the coil
    
    js[ Ind ] = Ns * Val_Current * Vector[-Sin[Atan2[Z[]-zc,X[]-xc]],
											  0.,
											  Cos[Atan2[Z[]-zc,X[]-xc]]] ;
}

Constraint {
    { Name MagneticVectorPotential_3D; Type Assign;
        Case {
            { Region SurfaceGInf ; Value 0.; }
			{ Region SurfaceGe0 ; Value 0.; }
        }
    }

    { Name Gauge; Type Assign;
         Case {
               { Region Domain_Mag; SubRegion #{SurfaceGInf, SurfaceGe0, 
                  SkinPlate, SkinCore}; Value 0.; }
   }

}
}
Jacobian {
  { Name Vol ;
    Case {
      { Region All ; Jacobian Vol ; }
    }
  }
}

Integration {
    { Name CurlCurl ;
        Case { {Type Gauss ;
            Case {  
                    { GeoElement Triangle ; NumberOfPoints 4 ; }
                    { GeoElement Tetrahedron ; NumberOfPoints 4 ; }
                 }
            }
        }
    }
}

Include "Magsta_a_3D_nl.pro"

PostOperation {
    { Name Map_a; NameOfPostProcessing MagSta_a_3D;
        Operation {
			Print[ a, OnElementsOf Domain_Mag, File "Magsta_a_3D_a.pos" ];
			Print[ a_norm, OnElementsOf Domain_Mag, File "Magsta_a_3D_a_norm.pos" ];
            Print[ b,  OnElementsOf Domain_Mag, File "Magsta_a_3D_b.pos" ];
            Print[ js, OnElementsOf DomainS_Mag, File "Magsta_a_3D_js.pos" ];
			Print[ js, OnElementsOf DomainS_Mag, Format Table,File  "js.txt" ];
        }
    }
}

// EOF
-------------- next part --------------
// Fichier "OneTooth-dat.pro"

mm = 1.e-3; // mm -> m

/*
  x -> longueur		-- gauche, droite
  y -> hauteur		-- haut, bas
  z -> ?paisseur	-- avant, arri?re
*/

// core
lx_core = 30. * mm;
ly_core = 100. * mm;
lz_core = 30. * mm;

// coil
th_coil = 10. * mm;
ri_coil = (lx_core^2. + lz_core^2.)^0.5 / 2. + 1. * mm;
re_coil = ri_coil + th_coil;
ly_coil = 50. * mm;

// plate
airgap = 2. * mm;
th_plate = 20. * mm;
ri_plate = ((lx_core/2.)^2. + (ly_core/2.)^2.)^0.5 + airgap;
re_plate = ri_plate+th_plate;
lz_plate = lz_core + 10. * mm;

// la bo?te ? air
bordure = 20. * mm;
axl = -re_plate - bordure; //left of airbox 
axr =  re_plate + bordure; //right of airbox 
ayb = -re_plate - bordure; //bottom of airbox
ayt =  re_plate + bordure; //top of airbox
azb =  0.; //back of airbox
azf =  re_coil +  bordure; //front of airbox

// les entit?s physiques
AIR    = 9001;
CORE   = 9002;
COIL   = 9003;
PLATE  = 9004;

AIR_INF      = 9011;
SKIN_PLATE   = 9012;
SKIN_CORE    = 9013;
CUT_SYMMETRY = 9014;

// EOF