[Gmsh] length scale dependence
ghather at uclink.berkeley.edu
ghather at uclink.berkeley.edu
Thu May 30 03:10:35 CEST 2002
So I'm using gmsh and getdp to find the polarization inside a dialetric object. The object is inside a uniform eletric field. As a test, I put a sphere inside a box and the polarization was correct to within 4%.
But i've noticed that the results change with the length scale of the problem. A 1m sphere in a 10m box gives a slightly different result than a 10m sphere in a 100m box (w/ characteristic lengths also scaled).
Also I tried to get a more accurate result by giving the sphere more points (so the sphere becomes more sphereical). Somehow, the polarization error became greater and the field inside the sphere became less uniform.
Off of the top of your head, can you think of any reasons for this?
--Greg Hather
(1st year at UC Berkeley)
HERE IS THE GEO
r = 1;
boxlength = 14*r;
cl1 = .1*boxlength;
cl2 = .2*r;
nth = 15;
nphi = nth;
//they must be equal!
/*
r = 1;
boxlength = 20*r;
cl1 = .05*boxlength;
cl2 = .1*r;
nth = 15;
nphi = nth;
//they must be equal!
*/
temp1 = nth*nphi;
//points, 1000s
For t In {1:temp1}
phinumber = Fmod(t-1,nth)+1;
thnumber = ((t-(Fmod(t-1,nth)+1))/nth) + 1;
phi = (2*Pi/(nphi))*(phinumber-1);
th = ((Pi)/(nth+1))*(thnumber);
Point(1000+t) = {r*Sin(th)*Cos(phi),r*Sin(th)*Sin(phi),r*Cos(th),cl2};
EndFor
//lines in phi dir, 2000s
For t In {1:temp1}
phinumber = Fmod(t-1,nth)+1;
thnumber = ((t-(Fmod(t-1,nth)+1))/nth) + 1;
tthp1 = (phinumber) + ((thnumber+1)-1)*nth;
tphip1 = (phinumber+1) + ((thnumber)-1)*nphi;
tphiback = (1) + ((thnumber)-1)*nphi;
If (phinumber < nphi)
Line(2000 + t) = {1000 + t,1000 + tphip1};
EndIf
If (phinumber == nphi)
Line(2000 + t) = {1000+t, 1000 + tphiback};
EndIf
EndFor
//lines in th dir, 3000s
For t In {1:temp1}
phinumber = Fmod(t-1,nth)+1;
thnumber = ((t-(Fmod(t-1,nth)+1))/nth) + 1;
tthp1 = (phinumber) + ((thnumber+1)-1)*nth;
tphip1 = (phinumber+1) + ((thnumber)-1)*nphi;
tphiback = (1) + ((thnumber)-1)*nphi;
If (thnumber < nth)
Line(3000 + t) = {1000 + t,1000 + tthp1};
EndIf
EndFor
//computer gen all surface but the ends
//line loops 4000, planes 5000
For t In {1:temp1}
phinumber = Fmod(t-1,nth)+1;
thnumber = ((t-(Fmod(t-1,nth)+1))/nth) + 1;
tthp1 = (phinumber) + ((thnumber+1)-1)*nth;
tphip1 = (phinumber+1) + ((thnumber)-1)*nphi;
tphiback = (1) + ((thnumber)-1)*nphi;
tphibackp1 = (2) + ((thnumber)-1)*nphi;
If ((thnumber < nth) && (phinumber < nphi))
Line Loop(4000 + t) = {3000 + t,2000 + tthp1, -(3000 + tphip1), -(2000 + t)};
Plane Surface(5000 + t) = {4000 + t};
EndIf
If ((thnumber < nth) && (phinumber == nphi))
Line Loop(4000 + t) = {3000 + t,2000 + tthp1, -(3000 + tphiback), -(2000 + t)};
Plane Surface(5000 + t) = {4000 + t};
EndIf
EndFor
//Extra stuff, like ends is 6000
tv1 = 2000 + nth;
tv2 = 2000 + (nth)*(nth-1) + 1;
tv3 = 2000 + nth^2;
tv4 = 5000 + (nth)*(nth-1);
Line Loop(6001) = {2001:tv1};
Line Loop(6002) = {tv2:tv3};
Plane Surface(6003) = {6001};
Plane Surface(6004) = {6002};
Surface Loop(6005) = {5001:tv4,6003,6004};
//oh, and a box to keep it in
Point(1) = {boxlength/2,-boxlength/2,-boxlength/2,cl1};
Point(2) = {boxlength/2,boxlength/2,-boxlength/2,cl1};
Point(3) = {-boxlength/2,boxlength/2,-boxlength/2,cl1};
Point(4) = {-boxlength/2,-boxlength/2,-boxlength/2,cl1};
Point(5) = {boxlength/2,-boxlength/2,boxlength/2,cl1};
Point(6) = {boxlength/2,boxlength/2,boxlength/2,cl1};
Point(7) = {-boxlength/2,boxlength/2,boxlength/2,cl1};
Point(8) = {-boxlength/2,-boxlength/2,boxlength/2,cl1};
Line(9) = {1,2};
Line(10) = {2,3};
Line(11) = {3,4};
Line(12) = {4,1};
Line(13) = {5,6};
Line(14) = {6,7};
Line(15) = {7,8};
Line(16) = {8,5};
Line(17) = {2,6};
Line(18) = {3,7};
Line(19) = {4,8};
Line(20) = {1,5};
Line Loop(21) = {9,17,-13,-20};
Line Loop(22) = {10,18,-14,-17};
Line Loop(23) = {11,19,-15,-18};
Line Loop(24) = {12,20,-16,-19};
Line Loop(25) = {9,10,11,12};
Line Loop(26) = {13,14,15,16};
Plane Surface(27) = {21};
Plane Surface(28) = {22};
Plane Surface(29) = {23};
Plane Surface(30) = {24};
Plane Surface(31) = {25};
Plane Surface(32) = {26};
Physical Surface(36) = {27};
Physical Surface(37) = {28};
Physical Surface(38) = {29};
Physical Surface(39) = {30};
Physical Surface(40) = {31};
Physical Surface(41) = {32};
Surface Loop(33) = {27,28,29,30,31,32};
Volume(34) = {33,6005};
Volume(6006) = {6005};
Physical Volume(35) = {34};
Physical Volume(6007) = {6006};
HERE IS THE PRO
Group {
face1 = Region[36];
face2 = Region[37];
face3 = Region[38];
face4 = Region[39];
face5 = Region[40];
face6 = Region[41];
air = Region[35];
die = Region[6007];
inside = Region[{air,die}];
}
Function{
//a0 = 1;
//b0 = 0;
c0 = 3;
epsr[die] = c0; //+ b0*Tanh[a0*Norm[$1]];
epsr[air] = 1; //+ 0*Norm[$1];
}
Constraint{
{Name pol1; Type Assign;
Case{
{Region face1; Value $X;}
{Region face2; Value $X;}
{Region face3; Value $X;}
{Region face4; Value $X;}
{Region face5; Value $X;}
{Region face6; Value $X;}
}}
{Name pol2; Type Assign;
Case{
{Region face1; Value $Y;}
{Region face2; Value $Y;}
{Region face3; Value $Y;}
{Region face4; Value $Y;}
{Region face5; Value $Y;}
{Region face6; Value $Y;}
}}
}
FunctionSpace{
{Name vspace; Type Form0;
BasisFunction{
{Name sn; NameOfCoef vn; Function BF_Node;
Support inside; Entity NodesOf[All];}
}
Constraint{
{NameOfCoef vn; EntityType NodesOf;
NameOfConstraint pol1;}
}}
}
Jacobian{
{Name jac1;
Case{
{Region inside; Jacobian Vol;}
}}
}
Integration{
{Name int1;
Case{
{Type Gauss;
Case{
{GeoElement Triangle; NumberOfPoints 4;}
{GeoElement Quadrangle; NumberOfPoints 4;}
{GeoElement Tetrahedron; NumberOfPoints 4;}
{GeoElement Hexahedron; NumberOfPoints 6;}
{GeoElement Prism; NumberOfPoints 9;}
{GeoElement Pyramid; NumberOfPoints 5;}
}}
}}
}
Formulation{
{Name electroform; Type FemEquation;
Quantity{
{Name v; Type Local; NameOfSpace vspace;}
}
Equation{
Galerkin{[epsr[]*Dof{Grad v}, {Grad v}]; In inside; Jacobian jac1; Integration int1;}
}}
}
Resolution{
{Name eres;
System{
{Name esyst; NameOfFormulation electroform;}
}
Operation{
Generate[esyst]; Solve[esyst];
//InitSolution[esyst];SaveSolution[esyst];
// max iterations, rel error, relaxation factor
//aperently, generatejac etc must be used with nonlinear iteration
//IterativeLoop[10, .001, 1]{
//Generate[esyst]; Solve[esyst];
// GenerateJac[esyst]; SolveJac[esyst];
//}
SaveSolution[esyst];
}}
}
PostProcessing{
{Name epostp; NameOfFormulation electroform;
Quantity{
//{Name vin; Value{ Local{[{v}]; In die;}} }
//{Name ein; Value{ Local{[{Grad v}]; In die;}}}
//{Name vout; Value{ Local{[{v}]; In air;}} }
//{Name eout; Value{ Local{[{Grad v}]; In air;}}}
//{Name vall; Value{ Local{[{v}]; In inside;}}}
{Name eall; Value{ Local{[Norm[{Grad v}]]; In inside;}}}
{Name nete; Value{ Integral{[{Grad v}]; In die; Integration int1;}}}
{Name netv; Value{ Integral{[1]; In die; Integration int1;}}}
{Name esnet; Value{ Integral{[SquNorm[{Grad v}]]; In die; Integration int1;}}}
//{Name epsin; Value{ Local{[epsr[{Grad v}]]; In die;}}}
//[c0 + b0*Tanh[a0*Norm[{Grad v}]]]
}}
}
PostOperation{
{Name eposto; NameOfPostProcessing epostp;
Operation{
//Print[vin, OnElementsOf die, File "vinlsphere.pos"];
//Print[ein, OnElementsOf die, File "einlsphere.pos"];
//Print[vout, OnElementsOf air, File "voutlsphere.pos"];
//Print[eout, OnElementsOf air, File "eoutlsphere.pos"];
//Print[vall, OnLine{{0,0,-.5} {0,0,.5}} {100}, Format Table, File "valllsphere.txt"];
Print[eall, OnLine{{0,0,-7} {0,0,7}} {100}, Format Table, File "ealllsphere.txt"];
Print[nete[die], OnRegion die, Format Table, File "dipollsphere.txt"];
//Print[enets[die], OnRegion die, Format Table, File "dipollsphere.txt"];
Print[netv[die], OnRegion die, Format Table, File >> "dipollsphere.txt"];
Print[esnet[die], OnRegion die, Format Table, File >> "dipollsphere.txt"];
//Print[epsin, OnElementsOf die, File "epsinsphere.pos"];
//Print[epsin, OnPoint{0, 0, 0}, File "epsinpsquare.pos"];
}}
}
--------------------------------------------------------------------
mail2web - Check your email from the web at
http://mail2web.com/ .