[ff3d-users] Re: Negative values, memory=none

[Joerg R. Weimar]

Dear Stephane,

I did not know that you have another job and just maintain ff3d on your 
spare time.

The CVS-server savannah is down for the rest of the week, so I don't 
expect any changes.

The version that I got
is from Dec 2. 2003
(example:
//  $Id: Mesh.hpp,v 1.7 2003/12/02 21:37:57 delpinux Exp $
)

In this version, the problem with negative values is not yet fixed.
I can work around that by using an offset.

I removed the box which is equal to the mesh box in my pov-files, and 
ignore the inaccuracy coming from the inaccurate volume calculations. 
(Once you have calculated the marching cube intersection, should not an 
accurate calculation of the intersection volume be fairly 
straightforward, by summing up the volumes of a few tetrahedrons (1/6*(a 
x b).c )?)


I now have the task of  using the solution of one PDE as boundary 
condition for a second PDE.
I use the following input, which does give me the desired value of u, 
but some problems persist with u2: according to the implied von-Neumann 
boundary conditions, I would expect an almost constant solution for u2 
below the object <0,0,1> and above the object <1,0,0>, which I do not 
get. It seems related to the point at the end of section 4.2.3 of the 
old freeFEM3d documentation, where it says "Add other remarks and 
Robin/Neumann boundar conditions description".

Best regards, Jörg Weimar.

--
PD Dr. Joerg R. Weimar, Inst. f. Wissensch. Rechnen, TU-Braunschweig
address@hidden, http://www.tu-bs.de/institute/WiR/weimar
Tel. +49-531-391-3006  Mail: D-38092 Braunschweig

vector n = (41, 41, 41);
vertex a = (-150, -150, -100);
vertex b = ( 150, 150, 400);
mesh M = structured(n,a,b);
 
scene S = pov("jrwtest3.pov"); // the pov-ray file for the geometry
domain filament = domain(S,outside(<0,0,1>) && outside(<1,0,0>) 
                    && inside(<0,1,0>));
domain air = domain(S,outside(<0,0,1>) && outside(<1,0,0>) 
                    && outside(<0,1,0>));


function u0 = (10.0);   
function u1 = u0 + 10.0;
function um = (u1+u0)/2;   

function filamentind = one(filament);   

function conductivity = 1.0 ;   
function permittivity = 1.0 ;   

solve(u) in filament by M 
//    memory(matrix=none),
//    method(type=penalty)
{
   pde(u)
     - div(conductivity*grad(u)) = 0;
     u = u0  on <0,0,1>;
     u = u1  on <1,0,0>;
     u = um on M xmin;
     u = um on M xmax;
};

double I=int(M)(u);
double J=int(M)(1);
cout << I << " " << J << " " <<I/J << "\n";
save(opendx,"u.dat",u,M);

solve(u2) in air by M 
//    memory(matrix=none),
//    method(type=penalty)
{
   pde(u2)
     - div(permittivity*grad(u2)) = 0;
     u2 = u   on <0,1,0>;
     u2 = u0  on <0,0,1>;
     u2 = u1  on <1,0,0>;
//     u2 = um on M xmin;
//     u2 = um on M xmax;
};

double I2=int(M)(u2);
double J2=int(M)(1);

cout << I2 << " " << J2 << " Average: " <<I2/J2 << "\n";

save(opendx,"u2.dat",u2,M);
save(opendx,"uu.dat",u*filamentind+u2*(1-filamentind),M);

/*
box {
   // Domain
   <-200, -200, -99>, <200, 200, 400>
   pigment {
      color rgb <1, 1, 0>
   }
}
*/

sphere {
      // UpperElectrode
      <0, 0, 300>, 30
   pigment {
      color rgb <1, 0, 0>
   }
}
   
cylinder {
   // LowerElectrode
   <0, 0, -100>, <0, 0, 0>, 50
   
   pigment {
      color rgb <0, 0, 1>
   }
}

union {
    // filament
    cylinder {
        <0, 0, -50>, <0, 0, 120>, 15
    }
    sphere {
      <0, 0, 120>, 20
    }
    cylinder {
        <0, 0, 120>, <100, 0, 50>, 15
    }
    sphere {
      <100, 0, 50>, 20
    }
    cylinder {
        <100, 0, 50>, <100, 0, 250>, 15
    }
    sphere {
      <100, 0, 250>, 20
    }
    cylinder {
        <100, 0, 250>, <0, 0, 200>, 15
    }
    sphere {
      <0, 0, 200>, 20
    }
    cylinder {
        <0, 0, 200>, <0, 0, 300>, 15
    }
   
   pigment {
      color rgb <0, 1, 0>
   }
}