# [Getdp] Extracting fluxes over an area

Stephen Brown beevesb at sover.net
Wed Jul 2 21:04:33 CEST 2008

```Hi I am doing some simulations of electrical conduction in a
heterogeneous volume with boundary conditions of fixed potentials over
small surface regions ("electrodes"). Many simulations will be done, and
the +- current electrodes will change from one to the next. I want to
calculate the total electric current injected by integrating the flux
(dn) over the electrode area. The electric potentials (n) will be
observed on the surface of the model (specifically an average over the
other unused electrodes). Is there a clever way to do this in
postprocessing, either in getdp or in gmsh?

Thanks for any  help or ideas.
Stephen Brown

Here is a model geometry and problem description:

**************************************************************
/* geometry.pro */
L = 15.0;
R = 30.5;
LEl = 15.24;
REl = 0.5;
Ip = 14461;
Im = 141461;

**************************************************************
/* cylinder.geo */
/* Cylinder with electrodes */

Include "geometry1.pro";

/* lc = (R/20);
lc2 = (R/160);
lc2 = (R/80); */
lc = (R/10);
lc2 = (R/20);

/* End faces */
Point(1) = {0.0, 0.0, 0.0, lc};
Point(2) = {R,   0.0, 0.0, lc};
Point(3) = {0.0, R,   0.0, lc};
Point(4) = {-R,  0.0, 0.0, lc};
Point(5) = {0.0, -R,  0.0, lc};

Point(11) = {0.0, 0.0, L, lc};
Point(12) = {R,   0.0, L, lc};
Point(13) = {0.0, R,   L, lc};
Point(14) = {-R,  0.0, L, lc};
Point(15) = {0.0, -R,  L, lc};

Circle(201) = {2,1,3};
Circle(202) = {3,1,4};
Circle(203) = {4,1,5};
Circle(204) = {5,1,2};

Circle(211) = {13,11,12};
Circle(212) = {14,11,13};
Circle(213) = {15,11,14};
Circle(214) = {12,11,15};

/* Sides */
Line(101) = {2, 12};
Line(102) = {3, 13};
Line(103) = {4, 14};
Line(104) = {5, 15};

Line(111) = {12, 2};
Line(112) = {13, 3};
Line(113) = {14, 4};
Line(114) = {15, 5};

/* Electrodes */
/* 46 56 66
45 55 65
44 54 64

Z=0 - 0
Z=L - 1    */

/* ******************************************* */
RCx = 0.0;
RCy = 0.0;
RCm = RCx - REl;
RCp = RCx + REl;
Point(551) = {RCx, RCy,  0, lc2};
Point(552) = {RCm, RCy,  0, lc2};
Point(553) = {RCp, RCy,  0, lc2};
Point(554) = {RCx, RCy-REl, 0, lc2};
Point(555) = {RCx, RCy+REl, 0, lc2};

Circle(551) = {552,551,554};
Circle(552) = {554,551,553};
Circle(553) = {553,551,555};
Circle(554) = {555,551,552};

RCx = -LEl;
RCy = 0.0;
RCm = RCx - REl;
RCp = RCx + REl;
Point(451) = {RCx, RCy,  0, lc2};
Point(452) = {RCm, RCy,  0, lc2};
Point(453) = {RCp, RCy,  0, lc2};
Point(454) = {RCx, RCy-REl, 0, lc2};
Point(455) = {RCx, RCy+REl, 0, lc2};

Circle(451) = {452,451,454};
Circle(452) = {454,451,453};
Circle(453) = {453,451,455};
Circle(454) = {455,451,452};

RCx = +LEl;
RCy = 0.0;
RCm = RCx - REl;
RCp = RCx + REl;
Point(651) = {RCx, RCy,  0, lc2};
Point(652) = {RCm, RCy,  0, lc2};
Point(653) = {RCp, RCy,  0, lc2};
Point(654) = {RCx, RCy-REl, 0, lc2};
Point(655) = {RCx, RCy+REl, 0, lc2};

Circle(651) = {652,651,654};
Circle(652) = {654,651,653};
Circle(653) = {653,651,655};
Circle(654) = {655,651,652};

RCx = 0.0;
RCy = 0.0;
RCm = RCx - REl;
RCp = RCx + REl;
Point(1551) = {RCx, RCy,  L, lc2};
Point(1552) = {RCm, RCy,  L, lc2};
Point(1553) = {RCp, RCy,  L, lc2};
Point(1554) = {RCx, RCy-REl, L, lc2};
Point(1555) = {RCx, RCy+REl, L, lc2};

Circle(1551) = {1552,1551,1554};
Circle(1552) = {1554,1551,1553};
Circle(1553) = {1553,1551,1555};
Circle(1554) = {1555,1551,1552};

RCx = -LEl;
RCy = 0.0;
RCm = RCx - REl;
RCp = RCx + REl;
Point(1451) = {RCx, RCy,  L, lc2};
Point(1452) = {RCm, RCy,  L, lc2};
Point(1453) = {RCp, RCy,  L, lc2};
Point(1454) = {RCx, RCy-REl, L, lc2};
Point(1455) = {RCx, RCy+REl, L, lc2};

Circle(1451) = {1452,1451,1454};
Circle(1452) = {1454,1451,1453};
Circle(1453) = {1453,1451,1455};
Circle(1454) = {1455,1451,1452};

RCx = +LEl;
RCy = 0.0;
RCm = RCx - REl;
RCp = RCx + REl;
Point(1651) = {RCx, RCy,  L, lc2};
Point(1652) = {RCm, RCy,  L, lc2};
Point(1653) = {RCp, RCy,  L, lc2};
Point(1654) = {RCx, RCy-REl, L, lc2};
Point(1655) = {RCx, RCy+REl, L, lc2};

Circle(1651) = {1652,1651,1654};
Circle(1652) = {1654,1651,1653};
Circle(1653) = {1653,1651,1655};
Circle(1654) = {1655,1651,1652};

/* ******************************************* */
RCx = 0.0;
RCy = +LEl;
RCm = RCx - REl;
RCp = RCx + REl;
Point(561) = {RCx, RCy,  0, lc2};
Point(562) = {RCm, RCy,  0, lc2};
Point(563) = {RCp, RCy,  0, lc2};
Point(564) = {RCx, RCy-REl, 0, lc2};
Point(565) = {RCx, RCy+REl, 0, lc2};

Circle(561) = {562,561,564};
Circle(562) = {564,561,563};
Circle(563) = {563,561,565};
Circle(564) = {565,561,562};

RCx = -LEl;
RCy = +LEl;
RCm = RCx - REl;
RCp = RCx + REl;
Point(461) = {RCx, RCy,  0, lc2};
Point(462) = {RCm, RCy,  0, lc2};
Point(463) = {RCp, RCy,  0, lc2};
Point(464) = {RCx, RCy-REl, 0, lc2};
Point(465) = {RCx, RCy+REl, 0, lc2};

Circle(461) = {462,461,464};
Circle(462) = {464,461,463};
Circle(463) = {463,461,465};
Circle(464) = {465,461,462};

RCx = +LEl;
RCy = +LEl;
RCm = RCx - REl;
RCp = RCx + REl;
Point(661) = {RCx, RCy,  0, lc2};
Point(662) = {RCm, RCy,  0, lc2};
Point(663) = {RCp, RCy,  0, lc2};
Point(664) = {RCx, RCy-REl, 0, lc2};
Point(665) = {RCx, RCy+REl, 0, lc2};

Circle(661) = {662,661,664};
Circle(662) = {664,661,663};
Circle(663) = {663,661,665};
Circle(664) = {665,661,662};

RCx = 0.0;
RCy = +LEl;
RCm = RCx - REl;
RCp = RCx + REl;
Point(1561) = {RCx, RCy,  L, lc2};
Point(1562) = {RCm, RCy,  L, lc2};
Point(1563) = {RCp, RCy,  L, lc2};
Point(1564) = {RCx, RCy-REl, L, lc2};
Point(1565) = {RCx, RCy+REl, L, lc2};

Circle(1561) = {1562,1561,1564};
Circle(1562) = {1564,1561,1563};
Circle(1563) = {1563,1561,1565};
Circle(1564) = {1565,1561,1562};

RCx = -LEl;
RCy = +LEl;
RCm = RCx - REl;
RCp = RCx + REl;
Point(1461) = {RCx, RCy,  L, lc2};
Point(1462) = {RCm, RCy,  L, lc2};
Point(1463) = {RCp, RCy,  L, lc2};
Point(1464) = {RCx, RCy-REl, L, lc2};
Point(1465) = {RCx, RCy+REl, L, lc2};

Circle(1461) = {1462,1461,1464};
Circle(1462) = {1464,1461,1463};
Circle(1463) = {1463,1461,1465};
Circle(1464) = {1465,1461,1462};

RCx = +LEl;
RCy = +LEl;
RCm = RCx - REl;
RCp = RCx + REl;
Point(1661) = {RCx, RCy,  L, lc2};
Point(1662) = {RCm, RCy,  L, lc2};
Point(1663) = {RCp, RCy,  L, lc2};
Point(1664) = {RCx, RCy-REl, L, lc2};
Point(1665) = {RCx, RCy+REl, L, lc2};

Circle(1661) = {1662,1661,1664};
Circle(1662) = {1664,1661,1663};
Circle(1663) = {1663,1661,1665};
Circle(1664) = {1665,1661,1662};

/* ******************************************* */
RCx = 0.0;
RCy = -LEl;
RCm = RCx - REl;
RCp = RCx + REl;
Point(541) = {RCx, RCy,  0, lc2};
Point(542) = {RCm, RCy,  0, lc2};
Point(543) = {RCp, RCy,  0, lc2};
Point(544) = {RCx, RCy-REl, 0, lc2};
Point(545) = {RCx, RCy+REl, 0, lc2};

Circle(541) = {542,541,544};
Circle(542) = {544,541,543};
Circle(543) = {543,541,545};
Circle(544) = {545,541,542};

RCx = -LEl;
RCy = -LEl;
RCm = RCx - REl;
RCp = RCx + REl;
Point(441) = {RCx, RCy,  0, lc2};
Point(442) = {RCm, RCy,  0, lc2};
Point(443) = {RCp, RCy,  0, lc2};
Point(444) = {RCx, RCy-REl, 0, lc2};
Point(445) = {RCx, RCy+REl, 0, lc2};

Circle(441) = {442,441,444};
Circle(442) = {444,441,443};
Circle(443) = {443,441,445};
Circle(444) = {445,441,442};

RCx = +LEl;
RCy = -LEl;
RCm = RCx - REl;
RCp = RCx + REl;
Point(641) = {RCx, RCy,  0, lc2};
Point(642) = {RCm, RCy,  0, lc2};
Point(643) = {RCp, RCy,  0, lc2};
Point(644) = {RCx, RCy-REl, 0, lc2};
Point(645) = {RCx, RCy+REl, 0, lc2};

Circle(641) = {642,641,644};
Circle(642) = {644,641,643};
Circle(643) = {643,641,645};
Circle(644) = {645,641,642};

RCx = 0.0;
RCy = -LEl;
RCm = RCx - REl;
RCp = RCx + REl;
Point(1541) = {RCx, RCy,  L, lc2};
Point(1542) = {RCm, RCy,  L, lc2};
Point(1543) = {RCp, RCy,  L, lc2};
Point(1544) = {RCx, RCy-REl, L, lc2};
Point(1545) = {RCx, RCy+REl, L, lc2};

Circle(1541) = {1542,1541,1544};
Circle(1542) = {1544,1541,1543};
Circle(1543) = {1543,1541,1545};
Circle(1544) = {1545,1541,1542};

RCx = -LEl;
RCy = -LEl;
RCm = RCx - REl;
RCp = RCx + REl;
Point(1441) = {RCx, RCy,  L, lc2};
Point(1442) = {RCm, RCy,  L, lc2};
Point(1443) = {RCp, RCy,  L, lc2};
Point(1444) = {RCx, RCy-REl, L, lc2};
Point(1445) = {RCx, RCy+REl, L, lc2};

Circle(1441) = {1442,1441,1444};
Circle(1442) = {1444,1441,1443};
Circle(1443) = {1443,1441,1445};
Circle(1444) = {1445,1441,1442};

RCx = +LEl;
RCy = -LEl;
RCm = RCx - REl;
RCp = RCx + REl;
Point(1641) = {RCx, RCy,  L, lc2};
Point(1642) = {RCm, RCy,  L, lc2};
Point(1643) = {RCp, RCy,  L, lc2};
Point(1644) = {RCx, RCy-REl, L, lc2};
Point(1645) = {RCx, RCy+REl, L, lc2};

Circle(1641) = {1642,1641,1644};
Circle(1642) = {1644,1641,1643};
Circle(1643) = {1643,1641,1645};
Circle(1644) = {1645,1641,1642};

/* hook it all together */
/* electrodes */
Line Loop(3551) = {553,554,551,552};
Plane Surface(4551) = {-3551};

Line Loop(3451) = {453,454,451,452};
Plane Surface(4451) = {-3451};

Line Loop(3651) = {653,654,651,652};
Plane Surface(4651) = {-3651};

Line Loop(31551) = {1553,1554,1551,1552};
Plane Surface(41551) = {31551};

Line Loop(31451) = {1453,1454,1451,1452};
Plane Surface(41451) = {31451};

Line Loop(31651) = {1653,1654,1651,1652};
Plane Surface(41651) = {31651};

Line Loop(3561) = {563,564,561,562};
Plane Surface(4561) = {-3561};

Line Loop(3461) = {463,464,461,462};
Plane Surface(4461) = {-3461};

Line Loop(3661) = {663,664,661,662};
Plane Surface(4661) = {-3661};

Line Loop(31561) = {1563,1564,1561,1562};
Plane Surface(41561) = {31561};

Line Loop(31461) = {1463,1464,1461,1462};
Plane Surface(41461) = {31461};

Line Loop(31661) = {1663,1664,1661,1662};
Plane Surface(41661) = {31661};

Line Loop(3541) = {543,544,541,542};
Plane Surface(4541) = {-3541};

Line Loop(3441) = {443,444,441,442};
Plane Surface(4441) = {-3441};

Line Loop(3641) = {643,644,641,642};
Plane Surface(4641) = {-3641};

Line Loop(31541) = {1543,1544,1541,1542};
Plane Surface(41541) = {31541};

Line Loop(31441) = {1443,1444,1441,1442};
Plane Surface(41441) = {31441};

Line Loop(31641) = {1643,1644,1641,1642};
Plane Surface(41641) = {31641};

/* sides */
Line Loop(301) = {111,201,-112,211};
Ruled Surface(401) = {301};
Line Loop(302) = {-113,212,112,202};
Ruled Surface(402) = {302};
Line Loop(303) = {113,203,-114,213};
Ruled Surface(403) = {303};
Line Loop(304) = {114,204,-111,214};
Ruled Surface(404) = {304};

/* end faces */
Line Loop(305) = {212,211,214,213};
Plane Surface(405) =
{-305,31461,31561,31661,31451,31551,31651,31441,31541,31641};

Line Loop(306) = {203,204,201,202};
Plane Surface(406) =
{-306,3461,3561,3661,3451,3551,3651,3441,3541,3641};

/* Entire outside surface */
Surface Loop(501) = {401,402,403,404,405,406,
41461,41561,41661,41451,41551,41651,41441,41541,41641,
4461,4561,4661,4451,4551,4651,4441,4541,4641};
Volume(601) = {501};

/* Surface of outer cylindrical shell */
Surface Loop(502) = {401,402,403,404};

Physical Surface(1001) = {501};
Physical Surface(1002) = {401,402,403,404};
Physical Surface(1003) = {401};
Physical Surface(1004) = {403};

Physical Surface(14461) = {4461};
Physical Surface(14561) = {4561};
Physical Surface(14661) = {4661};
Physical Surface(14451) = {4451};
Physical Surface(14551) = {4551};
Physical Surface(14651) = {4651};
Physical Surface(14441) = {4441};
Physical Surface(14541) = {4541};
Physical Surface(14641) = {4641};

Physical Surface(141461) = {41461};
Physical Surface(141561) = {41561};
Physical Surface(141661) = {41661};
Physical Surface(141451) = {41451};
Physical Surface(141551) = {41551};
Physical Surface(141651) = {41651};
Physical Surface(141441) = {41441};
Physical Surface(141541) = {41541};
Physical Surface(141641) = {41641};

Physical Volume(2001) = {601};

**************************************************************88
/* cylinder.pro *

Include "geometry1.pro"

Group {
Sur           = Region[1001];
OuterShell    = Region[1002];

ElectrodeIp    = Region[Ip];
ElectrodeIm    = Region[Im];

Vol           = Region[2001];
Tot = Region[{Vol,Sur}];
}

Function {
// not a matrix but just an isotropic,position-dependent constant

k[Vol]=5.0 + ((3000.0-5.0) * (1.0+Tanh[10.0*X[]])/2.0);
/*  k[Vol]=1.0; */

// An anisotropic Tensor
/*  k[Vol]=Tensor
[1.0,  0,  0,
0,  1.0,  0,
0,  0,  kzz];
*/

// A position dependent, anisotropic tensor
/*
k[Vol]=Tensor
[10.0+2.0*Tanh(X[]),  0,  0,
0,  10.0+2.0*Tanh(X[]),  0,
0,  0,  1.0];
*/

}

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

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 ; }
}
}
}
}
}

Constraint {
{ Name Sta_T ; Type Assign;
Case {
{ Region ElectrodeIp; Value 1.0; }
{ Region ElectrodeIm; Value 0.0; }
/*     { Region OuterShell; Value 0.0; } */
}
}
}

FunctionSpace {
{ Name spacen; Type Form0;
BasisFunction {
{ Name sn; NameOfCoef Tn; Function BF_Node; Support Tot;
Entity NodesOf[All]; }
}
Constraint {
{ NameOfCoef Tn; EntityType NodesOf ; NameOfConstraint Sta_T; }
}
}

}

Formulation {
{ Name f_n ; Type FemEquation;
Quantity {
{ Name n;  Type Local; NameOfSpace spacen; }
}
Equation {
Galerkin { [ k[]*Dof{d n} , {d n} ]; In Vol; Integration I1;
Jacobian JVol;  }
}
}
}

Resolution {
{ Name all;
System {
{ Name T; NameOfFormulation f_n; }
}
Operation {
Generate T ; Solve T ; SaveSolution T;
}
}
}

PostProcessing {
{ Name The; NameOfFormulation f_n;
Quantity {
{ Name n; Value { Local{ [ {n} ] ; In Vol; Jacobian JVol; } } }
{ Name dn; Value { Local{ [ {d n} ] ; In Vol; Jacobian JVol; } } }
/*
{ Name dnx ; Value { Term { [CompX[k[]*{d n}]]  ; In Vol; Jacobian
JVol; } } }
{ Name dny ; Value { Term { [CompY[k[]*{d n}]]  ; In Vol; Jacobian
JVol; } } }
{ Name dnz ; Value { Term { [CompZ[k[]*{d n}]]  ; In Vol; Jacobian
JVol; } } }
*/
}
}
}

PostOperation {
{ Name all ; NameOfPostProcessing The ;
Operation {
Print[ n, OnElementsOf Vol, File "n.pos"];
Print[ dn, OnElementsOf Vol, File "dn.pos"];
/*
Print [ dnx, OnElementsOf Vol, File "dnx.pos"] ;
Print [ dny, OnElementsOf Vol, File "dny.pos"] ;
Print [ dnz, OnElementsOf Vol, File "dnz.pos"] ;
*/
}
}
}

***************************************************************************
/* solver.par */

/*

Matrix_Format (Integer):
- 1  Sparse
- 2  Full
- default : 1

Matrix_Printing (Integer): Disk write ('fort.*')
- 1  matrix (csr)
- 2  preconditioner (msr)
- 3  both
- default : 0

Matrix_Storage (Integer): Disk Write or Read in internal format
- 0  none
- 1  write matrix (sparse)
- default : 0

Scaling (Integer): Scale system
- 0  no
- 1  on basis of diagonal elements  (no loss of possible symmetry)
- 2  on basis of inf. norm  of first rows and then columns
(asymmetric)
- 3  on basis of norm 1     of first rows and then columns
(asymmetric)
- 4  on basis of norm 2     of first rows and then columns
(asymmetric)
- default : 0

Renumbering_Technique (Integer):
- 0  No renumbering
- 1  Reverse Cuthill-Mc Kee
- default : 1

Preconditioner (Integer):
- 0  NONE     No Factorization
- 1  ILUT     Incomplete LU factorization with dual truncation
strategy
- 2  ILUTP    ILUT with column  pivoting
- 3  ILUD     ILU with single dropping + diagonal compensation
(~MILUT)
- 4  ILUDP    ILUD with column pivoting
- 5  ILUK     level-k ILU
- 6  ILU0     simple ILU(0) preconditioning
- 7  MILU0    MILU(0) preconditioning
- 8  DIAGONAL
- default : 2

Preconditioner_Position (Integer):
- 0  No Preconditioner
- 1  Left Preconditioner
- 2  Right Preconditioner
- 3  Both Left and Right Preconditioner
- default : 2

Nb_Fill (Integer):
- ILUT/ILUTP : maximum number of elements per line
of L and U (except diagonal element)
- ILUK : each element whose fill-in level is greater than NB_FILL
is dropped.
- default : 20

Permutation_Tolerance (Real): Tolerance for column permutation in
ILUTP/ILUDP.
At stage i, columns i and j are permuted if
abs(a(i,j))*PERMUTATION_TOLERANCE > abs(a(i,i)).
- 0  no permutations
- 0.001 -> 0.1  classical
- default : 0.05

Dropping_Tolerance (Real):
- ILUT/ILUTP/ILUK: a(i,j) is dropped if
abs(a(i,j)) < DROPPING_TOLERANCE * abs(diagonal element in U).
- ILUD/ILUDP : a(i,j) is dropped if
abs(a(i,j)) < DROPPING_TOLERANCE * [weighted norm of line i].
Weighted norm = 1-norm / number of nonzero elements on the line.
- default : 0

Diagonal_Compensation (Real): ILUD/ILUDP: the term
'DIAGONAL_COMPENSATION * (sum
of all dropped elements of the line)' is added to the diagonal
element in U
- 0  ~ ILU with threshold
1  ~ MILU with threshold.
- default : 0

Re_Use_ILU (Integer): Reuse ILU decomposition (and renumbering if any)
- 0  no
- 1  yes
- default : 0

Algorithm (Integer):
- 2  CGNR     CG (Normal Residual equation)
- 4  DBCG     BCG with partial pivoting
- 5  BCGSTAB  BCG stabilized
- 6  TFQMR    Transpose-Free Quasi-Minimum Residual
- 7  FOM      Full Orthogonalization Method
- 8  GMRES    Generalized Minimum RESidual
- 9  FGMRES   Flexible version of GMRES
- 10 DQGMRES  Direct versions of GMRES
- 11 LU       LU Factorization
- 12 PGMRES   Alternative version of GMRES
- default : 8

Krylov_Size (Integer): Krylov subspace size
- default : 40

IC_Acceleration (Real): IC accelerator
- default : 1

Re_Use_LU (Integer): Reuse LU decomposition
- 0  no
- 1  yes
- default : 0

Iterative_Improvement (Integer): Iterative improvement of the solution
obtained by a LU
- default : 0

Nb_Iter_Max (Integer): Maximum number of iterations
- default : 1000

Stopping_Test (Real): Target relative residual
- default : 1e-10

*/

Matrix_Format            1
Matrix_Printing            0
Matrix_Storage            0
Scaling            0
Renumbering_Technique            1
Preconditioner            2
Preconditioner_Position            2
Nb_Fill           20
Permutation_Tolerance         0.05
Dropping_Tolerance            0
Diagonal_Compensation            0
Re_Use_ILU            0
Algorithm            8
Krylov_Size           40
IC_Acceleration            1
Re_Use_LU            0
Iterative_Improvement            0
Nb_Iter_Max         1000
Stopping_Test        1e-10

--
Stephen Brown
beevesb at sover.net

```