/*
$Id: materialstr.cc,v 1.4 1996/11/19 10:02:53 bzferdma Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include "materialstr.h"
#include "cmdpars.h"
extern CmdPars Cmd;
//-------------------------------------------------------------------------
UserTransMaterial:: UserTransMaterial(const char* fileName, int spaceDim0)
: DefaultMaterial(fileName, spaceDim0)
{
setMaterialType();
}
//-------------------------------------------------------------------------
StepSource:: StepSource(const char* fileName, int spaceDim0)
: DefaultMaterial(fileName, spaceDim0)
{
setMaterialType();
}
Real StepSource:: S(int type, Vector<Real>* /*x*/, Real time)
{
if (type == 1 && time <= 10.0) return 1.0;
else return 0.0;
}
//-------------------------------------------------------------------------
TransPeakSource:: TransPeakSource(const char* fileName, int spaceDim0)
: PeakSource(fileName, spaceDim0)
{
setMaterialType();
}
Real TransPeakSource:: S(int /*type*/, Vector<Real>* x, Real time)
{
Vector<Real>& xx = *x;
Real val;
switch(spaceDim)
{
case 1: val = S1d(xx[1],time); return val;
case 2: val = S2d(xx[1], xx[2],time); return val;
case 3: val = S3d(xx[1], xx[2], xx[3],time); return val;
default: cout.flush(); abort(); return 0;
}
}
Real TransPeakSource:: S1d(Real x, Real time)
{
Real xm, t = time, x0;
Real f, fx, fxx, ft;
Real e, ex, exx;
Real uxx, p=100.0;
Real pi = 3.141592653589;
x0 = 0.5 + 0.375*sin(2.0*pi*t/20.0);
xm = x-x0;
e = exp(-p*(xm*xm));
ex = -2.0*p*xm*e;
exx = -2.0*p*e + (-2.0*p*xm)*ex;
f = x*(x-1.);
fx = (2.0*x-1.0);
fxx = 2.0;
ft = f*ex*(pi/10.0*cos(pi*t/10.0));
uxx = ft - ( fxx*e + 2.0*fx*ex + f*exx );
return uxx;
}
Real TransPeakSource:: S2d(Real x, Real y, Real time)
{
Real xm, ym;
Real f, fx, fy, fxx, fyy, ft;
Real e, ex, exx, ey, eyy;
Real uxx, uyy, p=100.0;
xm=x-0.5;
ym=y-0.5;
Real r = 0.375, t = time;
Real pi = 3.1414;
xm = 0.5 + r * cos(2.0*pi*t/20.0);
ym = 0.5 + r * sin(2.0*pi*t/20.0);
xm = x-xm;
ym = y-ym;
e = exp(-p*(xm*xm+ym*ym));
ex = -2.0*p*xm*e;
exx = -2.0*p*e + (-2.0*p*xm)*ex;
ey = -2.0*p*ym*e;
eyy = -2.0*p*e + (-2.0*p*ym)*ey;
f = x*(x-1.) * y*(y-1.);
fx = (2.0*x-1.0)*y*(y-1.);
fxx = 2.0*y*(y-1.);
fy = (2.0*y-1.0)*x*(x-1.);
fyy = 2.0*x*(x-1.);
uxx = fxx*e + 2.0*fx*ex + f*exx;
uyy = fyy*e + 2.0*fy*ey + f*eyy;
ft = -2.0/10.0*p*e*f*r*pi*( - xm*sin(pi*t/10.0) + ym*cos(pi*t/10.0) );
return (ft - uxx - uyy);
}
/* constant peak
Real TransPeakSource:: S2d(Real x, Real y, Real time)
{
Real xm, ym;
Real f, fx, fy, fxx, fyy;
Real e, ex, exx, ey, eyy;
Real uxx, uyy, p=100.0;
xm=x-0.5;
ym=y-0.5;
e = exp(-p*(xm*xm+ym*ym));
ex = -2.0*p*xm*e;
exx = -2.0*p*e + (-2.0*p*xm)*ex;
ey = -2.0*p*ym*e;
eyy = -2.0*p*e + (-2.0*p*ym)*ey;
f = x*(x-1.) * y*(y-1.);
fx = (2.0*x-1.0)*y*(y-1.);
fxx = 2.0*y*(y-1.);
fy = (2.0*y-1.0)*x*(x-1.);
fyy = 2.0*x*(x-1.);
uxx = fxx*e + 2.0*fx*ex + f*exx;
uyy = fyy*e + 2.0*fy*ey + f*eyy;
return (uxx + uyy);
}
*/
Real TransPeakSource:: S3d(Real x, Real y, Real z, Real /*time*/)
{
Real a, b, c, p;
Real xm1,ym1,zm1,xm1x,ym1y,zm1z,xma,ymb,zmc;
Real sum,divisor;
a = 0.5;
b = 0.5;
c = 0.5;
p = 150.0;
xma = x-a;
ymb = y-b;
zmc = z-c;
xm1 = x-1.0;
ym1 = y-1.0;
zm1 = z-1.0;
xm1x = x*xm1;
ym1y = y*ym1;
zm1z = z*zm1;
divisor = exp(p*(xma*xma + ymb*ymb + zmc*zmc));
sum = 2.0*(xm1x*ym1y+xm1x*zm1z+ym1y*zm1z);
sum = sum - 6.0*p*(xm1x*ym1y*zm1z);
sum = sum - 4.0*p*(xm1*xma*ym1y*zm1z + x*xma*ym1y*zm1z);
sum = sum + 4.0*p*p*(xm1x*xma*xma*ym1y*zm1z);
sum = sum - 4.0*p*(xm1x*ym1*ymb*zm1z + xm1x*y*ymb*zm1z);
sum = sum + 4.0*p*p*(xm1x*ym1y*ymb*ymb*zm1z);
sum = sum - 4.0*p*(xm1x*ym1y*zm1*zmc + xm1x*ym1y*z*zmc);
sum = sum + 4.0*p*p*(xm1x*ym1y*zm1z*zmc*zmc);
return -100.0*sum/divisor;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
CastingMaterial:: CastingMaterial(const char* fileName, int spaceDim0)
: DefaultMaterial(fileName, spaceDim0)
{
qMold = 1.0; Cmd.get("qMold",&qMold);
qSpray = 0.5; Cmd.get("qSpray",&qSpray);
setMaterialType();
}
Real CastingMaterial:: Neumann(int /*type*/, Vector<Real>* /*x*/, Real time)
{
switch(spaceDim)
{
case 1: cout.flush(); abort(); return 0.0;
case 2: if ( time < 12.5 ) return -qMold*0.142;
else return -0.1403*qSpray; // - 0.0328;
case 3: cout.flush(); abort(); return 0;
default: cout.flush(); abort(); return 0;
}
}
Real CastingMaterial:: Cauchy(int /*type*/, Vector<Real>* /*x*/, Real time)
{
switch(spaceDim)
{
case 1: cout.flush(); abort(); return 0;
case 2: if ( time < 12.5 ) return qMold * 2.0/10.0; // U <= 0 recommended!
else return qSpray*2.0/10.0;
// if ( (time > 12.5) && (time < 121) ) return 1.38*2.0/7.0;
case 3: cout.flush(); abort(); return 0;
default: cout.flush(); abort(); return 0;
}
}
//-------------------------------------------------------------------------
StefanSource:: StefanSource(const char* fileName, int spaceDim0)
: VarSourceMaterial(fileName, spaceDim0)
{
kappa0 = 1.0; Cmd.get("kappa0", &kappa0);
kappa1 = 1.0; Cmd.get("kappa1", &kappa1);
c0 = 1.0; Cmd.get("c0", &c0);
c1 = 1.0; Cmd.get("c1", &c1);
s0 = 0.0; Cmd.get("s0", &s0);
s1 = 0.0; Cmd.get("s1", &s1);
theta1 = 0.0; Cmd.get("StefanTheta1", &theta1);
setMaterialType();
}
//-------------------------------------------------------------------------
Real StefanSource:: S(int /*type*/, Vector<Real>* x, Real time)
{
Vector<Real>& xx = *x;
switch(spaceDim)
{
case 1: cout.flush(); abort(); return 0;;
case 2: return S2d(xx[1], xx[2], time);
case 3: cout.flush(); abort(); return 0;
default: cout.flush(); abort(); return 0;
}
// return 0;
}
//-------------------------------------------------------------------------
Real StefanSource:: S2d(Real x, Real y, Real t)
{
Real f=0.0, theta;
theta = (x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)-exp(-4.*t)/4.;
if (theta < theta1) f = c0 *exp(-4.*t) - 4*kappa0;
if (theta == theta1) f =(c0+c1)*exp(-4.*t)/2. - 2*(kappa0+kappa1);
if (theta > theta1) f = c1 *exp(-4.*t) - 4*kappa1;
return f;
}
syntax highlighted by Code2HTML, v. 0.9.1