[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)
- 2 read 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):
- 1 CG Conjugate Gradient
- 2 CGNR CG (Normal Residual equation)
- 3 BCG Bi-Conjugate Gradient
- 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