# [Getdp] given current => calculate the current density

Peter Kaufmann mail at peter.in-berlin.de
Fri Apr 18 17:42:35 CEST 2008

```Hello!

My problem sounds simple (and I think it have to be simple), but I did
not manage to solve it.

I want to calculate the current density in a 3D Volume. The current on
two surfaces are given (incoming and outgoing), the other surfaces are
isolated:
______
/_____/|
incoming |    | | outgoing
=>      |    | | =>
current  |____|/  current

My first attempt (testA.pro) was to set two Dirichlet Constraints
(incoming and outgoing current), use the equation of continuity and
calculate the current density direct via

Galerkin { [  Dof{j} ,  {d j} ] ;
In v_conductor ; Jacobian JVol ; Integration I1 ; }

... but I only get a current density on the two surfaces.

My second attempt was to calculate first the scalar potential and then
via PostProcessing the current density (testB.pro). The problem here is -
as far as I understand the problem - the Neumann Constraint on the

Generally

j[] = - sigma[] * {Grad{ phi }

holds. So on the surfaces I could write the Neumann Constraint using the
following Galerkin equation:

Galerkin { [ j[] + sigma[] * {Grad phi} , Dof {phi} ] ;
In s_current_in ; Jacobian JSur ; Integration I1 ; }

... but running getdp I only get the error message "Wrong Initial
Constraints: remaining Dof(s) with non-fixed initial conditions".

Maybe some one can give me a hint or an example?

Is there a way to calculate the size of a surface (for example integrate
1 over the surface) in the "Function" section and use it in the
Constraint conditions?

Kind regards,

Peter Kaufmann

-------------- next part --------------
/*
test case for the calculation of a current density

works on a gold cube (1mm x 1 mm x 1mm)
*/

Group {

v_conductor = Region[38];
s_current_in = Region[40];
s_current_out = Region[39];
c_all = Region[{v_conductor, s_current_in, s_current_out}];
}

Function {

I = 1;	// current in Ampere
f = 1; // surface in mm??
j0[ c_all ] = I/f; // incoming current density

resistivity_au = 2.44e-6; // electrical resistivity in Ohm*mm
conductance_au = 1/resistivity_au; // electrical conductance in 1/(Ohm*mm)

sigma [ c_all ] = conductance_au; // electrical conductance

}

Constraint {

// set the incoming current density
{ Name jIn ; Type Assign;
Case{
//{ Region s_current_in ; Value Vector[j0[], 0, 0] ; }
{ Region s_current_in ; Value Vector[j0[]/2, 0, 0] ; }
{ Region s_current_out ; Value Vector[j0[], 0, 0] ; }
}
}

}

Jacobian {
{ Name JVol ;
Case {
{ Region All ;        Jacobian Vol ; }
}
}
{ Name JSur ;
Case {
{ Region All ;        Jacobian Sur ; }
}
}
}

Integration {
{ Name I1 ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point       ; NumberOfPoints  1 ; }
{ GeoElement Line        ; NumberOfPoints  3 ; }
{ GeoElement Triangle    ; NumberOfPoints  4 ; }
{ GeoElement Quadrangle  ; NumberOfPoints  4 ; }
{ GeoElement Tetrahedron ; NumberOfPoints  4 ; }
{ GeoElement Hexahedron  ; NumberOfPoints  6 ; }
{ GeoElement Prism       ; NumberOfPoints  6 ; }
}
}
}
}
}

FunctionSpace {
{ Name Jspace ; Type Form0 ;
BasisFunction {
{ Name sn ; NameOfCoef tn ; Function BF_Node ;
Support Region[ v_conductor ] ; Entity NodesOf[ All ] ; }
}
Constraint {
{ NameOfCoef tn ; EntityType NodesOf ; NameOfConstraint jIn; } // Dirichlet constraint
}
}
}

Formulation {
{ Name StaticJ ; Type FemEquation ;
Quantity {
{ Name j ; Type Local ; NameOfSpace Jspace ; }
}
Equation {
Galerkin { [  Dof{j} ,  {d j} ] ;  In v_conductor ; Jacobian JVol ; Integration I1 ; }
}
}
}

Resolution {
{ Name StaticJSolution ;
System {
{ Name StaticJEquations ; NameOfFormulation StaticJ ; }
}
Operation {
Generate[StaticJEquations] ; Solve[StaticJEquations] ; SaveSolution[StaticJEquations] ;
}
}
}

PostProcessing {
{ Name  StaticJPP; NameOfFormulation StaticJ ;
Quantity {
{ Name j ; Value { Local { [ {j} ]; In v_conductor ; Jacobian JVol ; } } }
}
}
}

PostOperation {

{ Name coils ; NameOfPostProcessing StaticJPP;
Operation {
Print[ j, OnElementsOf v_conductor, File "j.pos"] ;
}
}

}
-------------- next part --------------
/*
test case for the calculation of a current density
*/

Group {

v_conductor = Region[38];
s_current_in = Region[40];
s_current_out = Region[39];
c_all = Region[{v_conductor, s_current_in, s_current_out}];
}

Function {

I = 1;	// current in Ampere
f_in = 2; // surface in mm??
f_out = 1; // surface in mm??
j[ s_current_in ] = Vector[I/f_in, 0, 0]; // incoming current density
j[ s_current_out ] = Vector[I/f_out, 0, 0]; // outgoing current density

resistivity_au = 2.44e-6; // electrical resistivity in Ohm*mm
conductance_au = 1/resistivity_au; // electrical conductance in 1/(Ohm*mm)

sigma [ c_all ] = conductance_au; // electrical conductance

}

Constraint {

// set the electrical potential on the drain side to zero (only the
// potential difference matters!)

{ Name phiScalarPotential ; Type Assign;
Case{
{ Region s_current_out ; Value 0.0 ; }
}
}
}

Jacobian {
{ Name JVol ;
Case {
{ Region All ;        Jacobian Vol ; }
}
}
{ Name JSur ;
Case {
{ Region All ;        Jacobian Sur ; }
}
}
}

Integration {
{ Name I1 ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point       ; NumberOfPoints  1 ; }
{ GeoElement Line        ; NumberOfPoints  3 ; }
{ GeoElement Triangle    ; NumberOfPoints  4 ; }
{ GeoElement Quadrangle  ; NumberOfPoints  4 ; }
{ GeoElement Tetrahedron ; NumberOfPoints  4 ; }
{ GeoElement Hexahedron  ; NumberOfPoints  6 ; }
{ GeoElement Prism       ; NumberOfPoints  6 ; }
}
}
}
}
}

FunctionSpace {
{ Name phi ; Type Form0 ;
BasisFunction {
{ Name sn ; NameOfCoef tn ; Function BF_Node ;
Support Region[ v_conductor ] ; Entity NodesOf[ All ] ; }
}
Constraint {
{ NameOfCoef tn ; EntityType NodesOf ; NameOfConstraint phiScalarPotential; } // Dirichlet constraint
}
}
}

Formulation {
{ Name StaticPotential ; Type FemEquation ;
Quantity {
{ Name phi ; Type Local ; NameOfSpace phi ; }
}
Equation {
Galerkin { [ sigma[] *  Dof {Grad phi} ,  {Grad phi} ] ;
In v_conductor ; Jacobian JVol ; Integration I1 ; }

Galerkin { [ j[] + sigma[] * {Grad phi} , Dof {phi} ] ;	// Neumann Constraint
In s_current_in ; Jacobian JSur ; Integration I1 ; }
Galerkin { [ j[] + sigma[] * {Grad phi} , Dof {phi} ] ;	// Neumann Constraint
In s_current_out ; Jacobian JSur ; Integration I1 ; }
}
}
}

Resolution {
{ Name StaticPotentialSolution ;
System {
{ Name StaticPotentialEquations ; NameOfFormulation StaticPotential ; }
}
Operation {
Generate[StaticPotentialEquations] ; Solve[StaticPotentialEquations] ; SaveSolution[StaticPotentialEquations] ;
}
}
}

PostProcessing {
{ Name  StaticPotentialPP; NameOfFormulation StaticPotential ;
Quantity {
{ Name phi ; Value { Local { [ {phi} ]; In v_conductor ; Jacobian JVol ; } } }
{ Name j; Value { Local { [ -sigma[] * {Grad phi} ]; In v_conductor ; Jacobian JVol ; } } }
}
}
}

PostOperation {
{ Name coils ; NameOfPostProcessing StaticPotentialPP;
Operation {
Print[ phi, OnElementsOf v_conductor, File "phi.pos"] ;
Print[ j, OnElementsOf v_conductor, File "j.pos"] ;
}
}
}
-------------- next part --------------
a=1;
Point(1) = {0, -0.5, 0, 0.5*a};
Point(2) = {0, 1.5, 0, 0.5*a};
Point(3) = {1, 1, 0, 0.5*a};
Point(4) = {1, 0, 0, 0.5*a};
Line (1) = {1, 2};
Line (2) = {2, 3};
Line (3) = {3, 4};
Line (4) = {4, 1};
Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};
Extrude {0, 0, 1} {
Surface{6};
}

Line Loop(29) = {3,22,-10,-18};
Plane Surface(30) = {29};
Line Loop(31) = {8,-14,-1,13};

Plane Surface(32) = {31};
Line Loop(33) = {9,10,11,8};
Plane Surface(34) = {33};
Plane Surface(35) = {5};
Surface Loop(36) = {28,23,6,19,27,15};
Volume(37) = {36};
Physical Volume(38) = {1}; // volume
Physical Surface(39) = {23}; // right surface
Physical Surface(40) = {15}; // left surface
```