[ff3d-users] problem

[Robert Li]

Dear Stephane Del Pino,
  When I tested my sample, I found that the result of
FEM and FDM is similiar, but is different. It seems
that there is scale between them.
  When I traced the source,I found that there are
differnet Vertex coordinate value to same Verter
No(eg. I=78 or i=78) in 
template <typename SurfaceMeshType> 
   void eval(const SurfaceMeshType& surfMesh) const
and 
void InstructionSaveFunction::execute().
Is this right?


Robert


>>>>>>>>>>>>>>>>>>>>>>
BoundaryConditionDiscretizationElimination.hpp

  template <typename SurfaceMeshType> 
   void eval(const SurfaceMeshType& surfMesh) const 
   {  

typedef
FiniteElementTraits<SurfaceMeshType::ElementGeometry> 
           FEType; 
 typedef typename FEType::Transformation
ConformTransformation;
 for (typename SurfaceMeshType::const_iterator
icell(surfMesh);
        !(icell.end()); ++icell) {
const typename SurfaceMeshType::ElementGeometry& cell
= *icell;
const ConformTransformation CT(cell); 
    
           for (size_t k=0;
k<CT.numberOfQuadratureVertices(); k++) { 
             TinyVector<3> q = CT.quadratureVertex(k);

    
             // Index of the vertex of the mesh the
closer to q. 
//   #warning Use of Index is not necessary 
             const Index& iv =
__bc.mesh().vertexIndex(q); 
             const Vertex& V = __bc.mesh().vertex(iv);

    
           const real_t GValue = __dirichlet.g(V);  
2004/01/28 li comment to revise boundary conditions

             const size_t I =
__bc.__degreeOfFreedomSet(__equationNumber, 
                                                      
 __bc.mesh().vertexNumber(V)); 




>>>>>>>>>>>>>>>>>>>>
Instruction.cpp

void InstructionSaveFunction::execute()
{
  (*__mesh).execute();
  (*__function).execute();
  (*__fileName).execute();
  Information::instance().setMesh(__mesh);

  const std::string CR = (*__fileDescriptor).cr();

  const Mesh& mesh = static_cast<const
Mesh&>(*(*__mesh).mesh());

  switch((*__fileDescriptor).format()) {
  case (FileDescriptor::dataExplorer): {
    std::ofstream fout((*__fileName).value());
    if (fout.bad()) {
      fferr(0) << "error: cannot open file '" <<
(*__fileName).value() << "'\n";
      exit(1);
    }

    for (size_t i=0; i<mesh.numberOfVertices(); ++i) {
      const TinyVector<3>& X = mesh.vertex(i);
      fout << (*__function).value(X[0], X[1], X[2]) <<
CR;
    }
    break;
  }
  case (FileDescriptor::medit): {
    //! plots the solution
    std::string sol = (*__fileName).value();
    sol += ".bb";
    std::ofstream fsol(sol.c_str());
        if (!(fsol)) { //2004/01/24 li add
         fferr(0) << "error: cannot open file '" <<
sol << "'\n";
         exit(1); 
       } 

        

    fsol << "3 1 " << mesh.numberOfVertices() << " 2"
<< CR;
    ReferenceCounting<FunctionExpression> __f;

    switch ((*(*__function).value()).type()) {
    case FunctionExpression::fem: {
      FunctionExpressionFEM& f
        =
static_cast<FunctionExpressionFEM&>(*(*__function).value());
      if ((*f.mesh()).mesh() == ((*__mesh).mesh())) {
        for(size_t i=0; i<mesh.numberOfVertices(); ++i) {
          fsol << f.value(i) << CR;
          ///2004/01/27 li add debug start
                const TinyVector<3>& X = mesh.vertex(i);
                
        ///end
    }


>>>>>>>>>>>>>>.
heat5.txt
//This example is  dpendent for time t(transient
state)
//dt(H) !=0:
//k=1920000.0 density=3000000.0  Cap=0.85 
//g!=0:
//materials=1:

vector n = (20,20,20);
vector a = (-3,-3,-2);
vector b = (3,3,2);

scene S = pov("Heat5.pov"); // the pov-ray file for
the geometry

mesh M = structured(n,a,b);

domain O = domain(S, (inside(<1,0,0>) and
inside(<0,1,0>) and inside(<0,0,1>) and
                                inside(<1,1,0>) and inside(<1,0,1>) and
inside(<0,1,1>) ) /*or
                                (inside(<0.9,0,0>) and inside(<0,0.9,0>) and
inside(<0,0,0.9>) and
                                inside(<0.9,0.9,0>) and inside(<0.9,0,0.9>) and
inside(<0,0.9,0.9>) )*/
                                );

mesh Omesh = surface(O,M);
mesh s1  = extract(Omesh, 1);   //xmin of bottom box
mesh s2  = extract(Omesh, 2);   //xmax of bottom box
mesh s3  = extract(Omesh, 3);   //ymin of bottom box
mesh s4  = extract(Omesh, 4);   //ymax of bottom box
mesh s5  = extract(Omesh, 5);   //zmin of bottom box
mesh s6  = extract(Omesh, 6);   //zmax of bottom box
//mesh s7  = extract(Omesh, 7); //xmin of middle box
//mesh s8  = extract(Omesh, 8); //xmax of middle box
//mesh s9  = extract(Omesh, 9); //ymin of middle box
//mesh s10 = extract(Omesh, 10);        //ymax of middle box
//mesh s11 = extract(Omesh, 11);        //zmin of middle box
//mesh s12 = extract(Omesh, 12);        //zmax of middle box

femfunction Tn(M) = 0;  //T0; when t=0; T
double  i = 0;                  //loop times
double  dt = 0.1;               //time step
double  Kx=1920000.0;   //conductivity in x direction
double  Ky=1920000.0;   //conductivity in y direction
double  Kz=920000.0;    //conductivity in z direction
double  Cu=4184000.0;   //unfrozen volumetric heat
capacity(1000*4184)
double  Cf=1932000.0;   //frozen volumetric heat
capacity (920*2100)
double  Tpc=0.0;                //phase change temperature
function Cp=Cu*one(Tn>Tpc) +
Cf*one(Tn<=Tpc);//volumetric heat capacity
double  Lf=334000000.0; //Latent heat of fusion of
water
double  Theat=0.8;              //Normalized volumetric water
content(0-1)
function Theatu1=0.04*one(Tn<=-2) +
(0.5*Tn+1.05)*one(Tn<=-0.5 and Tn>-2) +
                                 (1.6*Tn)*one(Tn<=0 and Tn>-0.5) + 0 * 
one(Tn>0);
function factor=dt/(Cp+(Lf)*(Theat)*(Theatu1));
do {
         solve(T) in O by M {
                 pde(T)
                        T-(
dx(((Kx)*(factor))*dx(T))+dy(((Ky)*(factor))*dy(T))+dz(((Kz)*(factor))*dz(T))
) = Tn;
//                      T = -10 on s12;
//                  T = 0 on s11;
                    T = -10 on s6;
                    T = 0 on s5;
//                      T = -10 on M zmax;  
//                  T = 0 on M zmin;
         }

        //Save result 
        save (medit , "heat5.00".i , T, M, dos ) ;
        save (medit , "heat5.00".i,  M, dos ) ;

        // compute deltaT here
        Tn = T;
        i=i+1;

} while (i<=10);


>>>>>>>>>heat5.pov
plane  { <-1,0,0>, 2 pigment { color rgb <1,0,0> }}
plane  { <1,0,0>, 2 pigment { color rgb <0,1,0> }}
plane  { <0,-1,0>, 2 pigment { color rgb <0,0,1> }}
plane  { <0,1,0>, 2 pigment { color rgb <1,1,0> }}
plane  { <0,0,-1>, 2 pigment { color rgb <1,0,1> }}
plane  { <0,0,1>, 2 pigment { color rgb <0,1,1> }}

>>>>>>>>>heat6.txt
//This example is  dpendent for time t(transient
state)
//dt(H) !=0:
//k=1920000.0 density=3000000.0  Cap=0.85 
//g!=0:
//materials=1:

vector n = (20,20,20);
vector a = (-2,-2,-2);
vector b = (2,2,2);

scene S = pov("void.pov"); // the pov-ray file for the
geometry

mesh M = structured(n,a,b);

//domain O = domain(S, (inside(<1,0,0>) and
inside(<0,1,0>) and inside(<0,0,1>) and
//                              inside(<1,1,0>) and inside(<1,0,1>) and
inside(<0,1,1>) ) /*or
//                              (inside(<0.9,0,0>) and inside(<0,0.9,0>) and
inside(<0,0,0.9>) and
//                              inside(<0.9,0.9,0>) and inside(<0.9,0,0.9>) and
inside(<0,0.9,0.9>) )*/
//                              );
/*
mesh Omesh = surface(O,M);
mesh s1  = extract(Omesh, 1);   //xmin of bottom box
mesh s2  = extract(Omesh, 2);   //xmax of bottom box
mesh s3  = extract(Omesh, 3);   //ymin of bottom box
mesh s4  = extract(Omesh, 4);   //ymax of bottom box
mesh s5  = extract(Omesh, 5);   //zmin of bottom box
mesh s6  = extract(Omesh, 6);   //zmax of bottom box
//mesh s7  = extract(Omesh, 7); //xmin of middle box
//mesh s8  = extract(Omesh, 8); //xmax of middle box
//mesh s9  = extract(Omesh, 9); //ymin of middle box
//mesh s10 = extract(Omesh, 10);        //ymax of middle box
//mesh s11 = extract(Omesh, 11);        //zmin of middle box
//mesh s12 = extract(Omesh, 12);        //zmax of middle box
*/

//femfunction Tn(M) = 0;        //T0; when t=0; T
function Tn = 0;                        //T0; when t=0; T
double  i = 0;                  //loop times
double  dt = 0.1;               //time step
double  Kx=1920000.0;   //conductivity in x direction
double  Ky=1920000.0;   //conductivity in y direction
double  Kz=1920000.0;   //conductivity in z direction
double  Cu=4184000.0;   //unfrozen volumetric heat
capacity(1000*4184)
double  Cf=1932000.0;   //frozen volumetric heat
capacity (920*2100)
double  Tpc=0.0;                //phase change temperature
function Cp=Cu*one(Tn>Tpc) +
Cf*one(Tn<=Tpc);//volumetric heat capacity
double  Lf=334000000.0; //Latent heat of fusion of
water
double  Theat=0.8;              //Normalized volumetric water
content(0-1)
function Theatu1=0.04*one(Tn<=-2) +
(0.5*Tn+1.05)*one(Tn<=-0.5 and Tn>-2) +
                                 (1.6*Tn)*one(Tn<=0 and Tn>-0.5) + 0 * 
one(Tn>0);
function factor=dt/(Cp+(Lf)*(Theat)*(Theatu1));
do {
         solve(T) in M {
                 pde(T)
                        T-(
dx(((Kx)*(factor))*dx(T))+dy(((Ky)*(factor))*dy(T))+dz(((Kz)*(factor))*dz(T))
) = Tn;
//                      T = -10 on s12;
//                  T = 0 on s11;
//                  T = 0 on s6;
//                  T = 3 on s5;
                        T = -10 on M zmax;  
                    T = 0 on M zmin;
         }

        //Save result 
        save (medit , "heat6.00".i , T, M, dos ) ;
        save (medit , "heat6.00".i,  M, dos ) ;

        // compute deltaT here
        Tn = T;
        i=i+1;

} while (i<=10);

__________________________________
Do you Yahoo!?
Yahoo! SiteBuilder - Free web site building tool. Try it!
http://webhosting.yahoo.com/ps/sb/