[ff3d-users] Problem with elliptic hollow cylinder section

[Tobias Nähring]

Note:
A short summary of the problem with the corresponding pictures and
files is available at: http://www.tn-home.de/Soft/Fem3d/index.html.



I want to solve the Laplacian equation in the domain

 G := {(x,y,z)|   1 < (x/rx)^2 + (y/ry)^2 < a^2;   0 < z < L }

which resembles a section of a hollow cylinder with elliptical
cross-section. The constants are: rx=1, ry=0.5, a=2, L=5.

The boundary condition is

 u(x,y,z) = (|x,y|-1)*(a-|x,y|)*(L-z)/L

with

 |x,y|:=sqrt( (x/rx)^2 + (y/ry)^2 ).

That means effectively that u(x,y,z)=0 for z=L, and for (x,y,z) on the
cylinder mantles. Only at z=0 one has some nice non-zero boundary
condition.

Because of the symmetry of the boundary condition the
computational domain can be confined to the quarter

Q:={(x,y,z) in G| x>0, y>0}

of G with natural boundary conditions for u on the boundaries at x=0
and y=0.

I tried to solve this problem with FreeFEM3D without much luck.

The results of my attempts are shown at
http://www.tn-home.de/Tobias/Soft/Fem3d/index.html

ff3d accepts the input file and the corresponding pov file (included below)
and the solver options are set to the save values

cg(epsilon=1e-8) and cg(maxiter=2000)

but the output is rather erroneous. I have played
already with other solver settings, such as bicg(...) and
krylov(type=bicg) without any luck. The grid size is as large as
possible on my computer with 512MB RAM. (If I enlarge the grid size
then ff3d freezes during "Marching cube ...".

With the current grid size I get the error message:

 .
 warning : meshing object plane {
 (0, 1, 0),2
 }
  a problem occured ...
 Check that:
 - the object is not to small compare to the background mesh
 - the object intersects the background mesh

I don't know how to handle that.

The ff3d computed solution of the problem are graphically represented
in the figures:

http://www.tn-home.de/Tobias/Soft/Fem3d/

Is there any chance to get a better solution?

One more remark: There would be a way to specify the mantle surface of
a cylinder section: the `open' directive of the povray object
cylinder (as described in the povray manual).

The input file for ff3d:
//////////////////////////////////////////////////////////////////////
// Model of a cable with ellipsoidal cross-section.


vector Grid=(70,60,30);
double rx = 1.0; // radius of the inner conductor in x-direction
double ry = 0.5; // radius of the inner conductor in y-direction
double a=2.0; // The cross-section boundary of the inner conductor
         // scaled by a gives the cross-section boundary of the
         // shielding.
double CableW = a*rx;
double CableH = a*ry;

double CableL = 5.0; // length of the cable
double Border = 0.1; // Border added around the modelled piece of the
            // cable to build a fictitious domain

scene myScene = pov("cable-ellipse-annulus-num.pov");

mesh myMesh = structured( Grid, (-Border,-Border,-2.0*Border), 
(a*rx+Border,a*ry+Border,CableL+2.0*Border));
// Comment: in z-direcion we provided a bit more space for the 
fictitious domain (just guessing).

domain myDomain = domain(myScene, inside(<1,0,0>));
mesh myCylMesh = surface(myDomain,myMesh);

function norm2d=sqrt((x/rx)^2+(y/ry)^2);

// Realizing the boundary conditions with method-type penalty doesnot
// work (no convergence).
solve(u) in myDomain by myMesh
 cg(epsilon=1e-8),
 cg(maxiter=2000)
{
 pde(u) div(grad(u))=0;
 u=(norm2d-1)*(a-norm2d)*(CableL-z)/CableL on myCylMesh;
}

//////////////////////////////////////////////////////////////////////
// Results output for postprocessing via gnuplot:

double myx; double myy; double myz;
double mydx = CableW/x(Grid);
double mydy = CableH/y(Grid);
double mydz = CableL/z(Grid);
double myxmax=CableW+0.5*mydx;
double myymax=CableH+0.5*mydy;
double myzmax=CableL+0.5*mydz;

for( myz=0.0; myz <= myzmax; myz = myz + mydz ) {
 for( myx=0.0; myx <= myxmax; myx = myx + mydx ) {
   for( myy=0.0; myy <= myymax; myy = myy + mydy ) {
     cout << "@" << myx << ", " << myy << ", " << myz << ", " << 
u(myx,myy,myz) << "\n";
   }
   cout << "@\n";
 }
 cout << "@\n";
}
//////////////////////////////////////////////////////////////////////

The pov file (named cable-ellipse-annulus-num.pov):
//////////////////////////////////////////////////////////////////////
// Model of a cable with ellipsoidal cross-section.
// rx=1.0 radius of the inner conductor in x-direction
// ry=0.5 radius of the inner conductor in y-direction
// a=2.0 the inner conductor scaled by a gives the shielding
// CableL=5 length of the cable.
// Since we impose symmetric boundary conditions we need only
// to model a quarter of the cable. (0 < x < a*rx, 0 < y < a*ry)
intersection {
 difference {
   cylinder {
     <0,0,0>,
     <0,0,1>,
     2.0 // a==2.0
   }
   cylinder{
     <0,0,-0.01>,
     <0,0,1.01>,
     1.0
   }
   plane{
     <0,1,0>, 0
   }
   plane{
     <1,0,0>, 0
   }
 }
 box{
   <0,0,0>,
   <2,2,1>
 }
 scale <1.0,0.5,5.0> // <rx,ry,CableL>
 pigment { color rgb <1,0,0> }
}