[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[ff3d-users] problem
From: |
Robert Li |
Subject: |
[ff3d-users] problem |
Date: |
Mon, 2 Feb 2004 09:32:31 -0800 (PST) |
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/


- [ff3d-users] problem,
Robert Li <=