[Getdp] Help, Why does this work
Bernhard Kubicek
bernhard.kubicek at arsenal.ac.at
Tue Oct 30 12:32:21 CET 2007
Dear Christophe, Dear others,
I have a problem. I can calculate something in GetDP, and more or less
mostly fail at understanding why it works.
The goal is the following.
Calculate some electrical scalar potential in a domain. The domain is
separated into two halves. At the boundary connecting the halves there
should be a nonlinear potential jump, depending on the current flowing
through the halves.
-------- Potential boundary
| | ftop
|------| lborder, potential jump
| | fbottom
-------- Potential boundary
I therefore adopted my already strange "faceconduct2d" wiki example:
usual conductivity in the two subdomains
>Galerkin { [ sigma[]*Dof{d vtop} , {d vtop} ]; In ftop;}
>Galerkin { [ sigma[]*Dof{d vbottom} , {d vbottom} ]; In fbottom;}
on the border there is a strange interface condition:
>Galerkin { [-fsigma[]*Dof{vtop} , {vbottom} ]; In lborder;}
>Galerkin { [ fsigma[]*Dof{vbottom} , {vbottom} ]; In lborder; }
>Galerkin { [-fsigma[]* Dof{vbottom} , {vtop} ]; In lborder; }
>Galerkin { [ fsigma[]*Dof{vtop} , {vtop} ]; In lborder;}
>fsigma[]= border_sigma/border_virtual_thickness
All this leads to a potential jump like it would happen if an the
lborder there were a thick region having the given thickness and
conductivity.
To create a constant potential jump, I add the following two terms:
>Galerkin { [ UJUMP[]*fsigma[], {vtop} ]; In lborder;}
>Galerkin { [ -UJUMP[]*fsigma[] , {vbottom} ]; In lborder;}
with UJUMP[] being e.g. 5V. and set fsigma to a large value (=good
conducting, very thin)
I find it quite a miracle that this works.
Finally, to make the jump current dependent, it is even more tricky.
In lborder, there is no good definition of the current density because
of the potential discontinuity.
Therefore I use this very unelegant method:
I add another function space
>{ Name jline; Type Local; NameOfSpace fooo; }
defined on nodes of ftop and lborder
it is set equal to the current in ftop:
>Galerkin {[ Dof{jline}, { jline} ]; In ftop; }
>Galerkin {[ -CompY[sigma[]*{d vtop}] , { jline} ]; In ftop; }
the jump then of course has to be:
>Galerkin { [ UJUMP[{jline}]*fsigma[], {vtop} ]; In lborder;}
>Galerkin { [ -UJUMP[{jline}]*fsigma[] , {vbottom} ]; In lborder;}
with UJUMP[]=5+$1/100.;
This kind of works while solving using a iterative loop.
My problem is that I want to use it for bent "lborders", where I cannot
use CompY anymore. Also, the use of the jline space takes a lot of
unnecessary memory.
So my questions are:
*Is there a less insane way of doing this?
*Does somebody understand why this works?
My feeling is that this is only understandable at matrix level.
Attached are all necessary files.
very nice greetings,
and thank you very much
Bernhard
-------------- next part --------------
l=1.;
ml=l/80;
Point(1)={0,0,0,ml};
Point(2)={0,l,0,ml};
Point(3)={l,l,0,ml};
Point(4)={l,0,0,ml};
Point(5)={0,l/2.,0,ml};
Point(6)={l,l/2.,0,ml};
Line(1) = {1,4};
Line(2) = {4,6};
Line(3) = {6,3};
Line(4) = {3,2};
Line(5) = {2,5};
Line(6) = {5,1};
Line(7) = {5,6};
Line Loop(8) = {1,2,-7,6};
Plane Surface(9) = {8};
Line Loop(10) = {7,3,4,5};
Plane Surface(11) = {10};
Physical Line(100) = {1};
Physical Line(101) = {4};
Physical Line(200) = {7};
Physical Surface(500) = {9};
Physical Surface(501) = {11};
-------------- next part --------------
Group {
//lines
ltop=# 101 ;
lbottom=# 100 ;
lborder=# 200;
//faces
ftop= Region[500];
fbottom= Region[501];
//domains
vtopl = Region[{ftop,fbottom}] ;
Lines=Region[{ltop,lbottom,lborder}] ;
}
Function {
//conductivity sigma [S/m]
Curve[]=Exp[-(X[]-0.5)*(X[]-0.5)/0.15/0.15];
sigma[ftop]= Curve[];
sigmatop[lborder]= Curve[];
sigma[fbottom]= 10*Curve[];
sigmabot[lborder]= 10*Curve[];
border_sigma= 1;
//border_virtual_thickness= 0.5; //[meter]
border_virtual_thickness= 0.000001; //[meter]
//face conductivity [S/m^2]
fsigma[]= border_sigma/border_virtual_thickness;
//UJUMP[]=10-Fabs[CompY[$1]]/180.;
//UJUMP[]=20+X[]*15;
UJUMP[]=0.1+10* $1/100.;
}
Jacobian {
//elements extend indefinitely into z directions
{ Name JVol ;
Case {
{ Region Lines ; Jacobian Sur ; }
{ Region All ; Jacobian Vol ; }
}
}
}
Integration {
{ Name Int ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 3 ; }
{ GeoElement Triangle ; NumberOfPoints 4 ; }
}
}
}
}
}
Constraint {
{ Name Sta_V ; Type Assign;
Case {
{ Region ltop; Value 100.; }
{ Region lbottom; Value 0; }
}
}
}
FunctionSpace {
//spaces a and b both have independent nodal values at the border-line
{ Name spacea; Type Form0;
BasisFunction {
{ Name sn; NameOfCoef Tn; Function BF_Node; Support Region[{ftop,lborder}];
Entity NodesOf[All]; }
}
Constraint {
{ NameOfCoef Tn; EntityType NodesOf ; NameOfConstraint Sta_V; }
}
}
{ Name spaceb; Type Form0;
BasisFunction {
{ Name sn; NameOfCoef Tn; Function BF_Node; Support Region[{fbottom,lborder}];
Entity NodesOf[All]; }
}
Constraint {
{ NameOfCoef Tn; EntityType NodesOf ; NameOfConstraint Sta_V; }
}
}
{ Name fooo; Type Form0;
BasisFunction {
{ Name sn; NameOfCoef Tn; Function BF_Node; Support Region[{ftop,lborder}];
Entity NodesOf[All]; }
}
Constraint {
{ NameOfCoef Tn; EntityType NodesOf ; NameOfConstraint Foooo; }
}
}
}
Formulation {
{ Name f_v ; Type FemEquation;
Quantity {
{ Name vtop; Type Local; NameOfSpace spacea; }
{ Name vbottom; Type Local; NameOfSpace spaceb; }
{ Name jline; Type Local; NameOfSpace fooo; }
}
Equation {
//
Galerkin {[ sigma[]*Dof{d vtop} , {d vtop} ]; In ftop; Integration Int; Jacobian JVol;}
Galerkin { [ sigma[]*Dof{d vbottom} , {d vbottom} ]; In fbottom; Integration Int; Jacobian JVol; }
Galerkin { [-fsigma[]*Dof{vtop} , {vbottom} ]; In lborder;Jacobian JVol; Integration Int; }
Galerkin { [ fsigma[]*Dof{vbottom} , {vbottom} ]; In lborder;Jacobian JVol; Integration Int; }
Galerkin { [ -UJUMP[{jline}]*fsigma[] , {vbottom} ]; In lborder;Jacobian JVol; Integration Int; }
Galerkin { [-fsigma[]* Dof{vbottom} , {vtop} ]; In lborder;Jacobian JVol; Integration Int; }
Galerkin { [ fsigma[]*Dof{vtop} , {vtop} ]; In lborder;Jacobian JVol; Integration Int; }
Galerkin { [ UJUMP[{jline}]*fsigma[], {vtop} ]; In lborder;Jacobian JVol; Integration Int; }
Galerkin {[ Dof{jline}, { jline} ]; In ftop; Integration Int; Jacobian JVol;}
Galerkin {[ -CompY[sigma[]*{d vtop}] , { jline} ]; In ftop; Integration Int; Jacobian JVol;}
}
}
}
Resolution {
{ Name all;
System {
{ Name T; NameOfFormulation f_v; }
}
Operation {
InitSolution[T];
//IterativeLoop[10, 0, 0.75] {
IterativeLoop[30, 0, 0.3] {
GenerateJac[T]; SolveJac[T];
}
SaveSolution T ;
}
}
}
PostProcessing {
{ Name The; NameOfFormulation f_v;
Quantity {
{ Name v; Value{
Local{ [ {vtop} ] ; In ftop;Jacobian JVol; }
Local{ [ {vbottom} ] ; In fbottom;Jacobian JVol; }
} }
{ Name vt; Value{
Local{ [ {vtop} ] ; In ftop;Jacobian JVol; }
Local{ [ {vtop} ] ; In lborder;Jacobian JVol; }
} }
{ Name vb; Value{
Local{ [ {vbottom} ] ; In fbottom;Jacobian JVol; }
Local{ [ {vbottom} ] ; In lborder;Jacobian JVol; }
} }
{ Name dv; Value{
Local{ [ {vtop}-{vbottom} ] ; In lborder;Jacobian JVol; }
} }
{ Name j; Value{
Local{ [-sigma[]* {d vtop} ] ; In ftop;Jacobian JVol; }
Local{ [-sigma[]* {d vbottom} ] ; In fbottom;Jacobian JVol; }
} }
{ Name jl; Value{
Local{ [{jline} ] ; In ftop;Jacobian JVol; }
} }
}
}
}
PostOperation {
{ Name all ; NameOfPostProcessing The ;
Operation {
Print[ v, OnElementsOf vtopl , File "v.pos"];
Print[ vt, OnElementsOf Region[{ftop,lborder}], File "vt.pos"];
Print[ vb, OnElementsOf Region[{fbottom,lborder}], File "vb.pos"];
Print[ dv, OnElementsOf lborder, File "dv.pos"];
Print[ dv, OnElementsOf lborder, Format Table];
Print[ j, OnElementsOf vtopl , File "j.pos"];
Print[ jl, OnElementsOf ftop, File "jl.pos"];
}
}
}