/*
$Id: materialsA.cc,v 1.3 1996/11/19 09:57:55 bzferdma Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include "materialsA.h"
#include "physics.h"
#include "utils.h"
#include "cmdpars.h"
extern CmdPars Cmd;
//-------------------------------------------------------------------------
DefaultMaterial:: DefaultMaterial(int spaceDim0)
: Material(spaceDim0),
eVal(1), mVal(1), pVal(1), cVal(1), sVal(1), neumannVal(1), cauchyVal(1),
eTerm(1), mTerm(1), pTerm(1), cTerm(1), sTerm(1),
neumannTerm(1), cauchyTerm(1)
{
int i;
FORALL(eVal,i) eVal[i] = 0;
FORALL(mVal,i) mVal[i] = 0;
FORALL(pVal,i) pVal[i] = 0;
FORALL(cVal,i) cVal[i] = 0;
FORALL(neumannVal,i) neumannVal[i] = 0;
FORALL(cauchyVal,i) cauchyVal[i] = 0;
setMaterialType();
}
//-------------------------------------------------------------------------
DefaultMaterial:: DefaultMaterial(const char* fileName, int spaceDim0)
: Material(spaceDim0),
eVal(1), mVal(1), pVal(1), cVal(1), sVal(1),
neumannVal(1), cauchyVal(1),
eTerm(1), mTerm(1), pTerm(1), cTerm(1), sTerm(1),
neumannTerm(1), cauchyTerm(1)
{
int i, j, k;
char s[MaxLineLength];
if (Cmd.get("matFile",s))
{
if (!strchr(s,'.')) strcat(s,".mat");
}
else
{
strcpy(s,fileName);
if (strchr(s,'.')) strcpy(strchr(s,'.'),".mat");
else strcat(s,".mat");
}
BufferedParser parser(s, "end", CommentFlag);
readMaterials(parser);
// setMaterialType(); // not used: definition of matType in readMaterials
readBCValues (parser);
// -- read transform factors and apply them to the input values:
Real eFact=1.0, mFact=1.0, pFact=1.0, cFact=1.0, sFact=1.0,
neumannFact=1.0, cauchyFact=1.0;
readFactors (parser, &eFact, &mFact, &pFact, &cFact, &sFact,
&neumannFact, &cauchyFact);
if (spaceDim == 1)
{
FORALL(eVal,i)
{
FORALL_ROWS(*eVal[i],j)
FORALL_COLUMNS(*eVal[i],k)
{
(*eVal[i])(j,k) *= eFact;
(*mVal[i])(j,k) *= mFact;
(*pVal[i])(j,k) *= pFact;
}
}
FORALL(cVal,i) FORALL(*cVal[i],j) (*cVal[i])[j] *= cFact;
FORALL(sVal,i) sVal[i] *= sFact;
FORALL(neumannVal,i) neumannVal[i] *= neumannFact;
FORALL(cauchyVal,i) cauchyVal[i] *= cauchyFact;
}
else if (spaceDim == 2)
{
FORALL(eVal,i)
{
FORALL_ROWS(*eVal[i],j)
FORALL_COLUMNS(*eVal[i],k)
{
(*eVal[i])(j,k) *= eFact;
(*mVal[i])(j,k) *= mFact;
(*pVal[i])(j,k) *= pFact;
}
}
FORALL(cVal,i) FORALL(*cVal[i],j) (*cVal[i])[j] *= cFact;
FORALL(sVal,i) sVal[i] *= sFact;
FORALL(neumannVal,i) neumannVal[i] *= neumannFact;
FORALL(cauchyVal,i) cauchyVal[i] *= cauchyFact;
}
else if (spaceDim == 3)
{
FORALL(eVal,i)
{
FORALL_ROWS(*eVal[i],j)
FORALL_COLUMNS(*eVal[i],k)
{
(*eVal[i])(j,k) *= eFact;
(*mVal[i])(j,k) *= mFact;
(*pVal[i])(j,k) *= pFact;
}
}
FORALL(cVal,i) FORALL(*cVal[i],j) (*cVal[i])[j] *= cFact;
FORALL(sVal,i) sVal[i] *= sFact;
FORALL(neumannVal,i) neumannVal[i] *= neumannFact;
FORALL(cauchyVal,i) cauchyVal[i] *= cauchyFact;
}
}
//-------------------------------------------------------------------------
DefaultMaterial:: ~DefaultMaterial()
{
int i;
FORALL(eVal,i) delete eVal[i];
FORALL(mVal,i) delete mVal[i];
FORALL(pVal,i) delete pVal[i];
FORALL(cVal,i) delete cVal[i];
}
//-------------------------------------------------------------------------
void DefaultMaterial:: setMaterialType()
{
// Do nothing, matType is defined in readMaterials
// cout << "In DefaultMaterial::setMaterialType: ";
// switch (matType)
// {
// case anisotropic: cout << "anisotropic" << endl; break;
// case isotropic: cout << "isotropic" << endl; break;
// default : cout << "undefinedmaterial" << endl;
// }
}
//-------------------------------------------------------------------------
void DefaultMaterial:: readMaterials(Parser& parser)
{
int i, j, k, dim=0, classA, maxClass=0, minClass=1;
Real val;
Matrix<Real>* mat=0;
char s[MaxLineLength];
parser.rewind();
if (parser.findWord("Mat") != parser.ok)
{
cout <<"\n*** No Materials encountered in " << s << "\n";
cout.flush(); abort();
}
long position = parser.ftell();
parser.getString(s);
strToLower(s);
if (strstr(s,"ani")) { matType = anisotropic; dim = spaceDim; }
else if (strstr(s,"iso")) { matType = isotropic; dim = 1; }
else {
cout << "\n*** Materials: unknown material type read: " << s << "\n";
cout.flush(); abort();
}
// look at the entries:
while (parser.getValue(&classA) == parser.ok)
{
if (classA <= 0) {
cout << "\n*** Class Material: identifiers for Materials"
<< " must be > 0 (read: " << classA << ")\n";
cout.flush(); abort();
}
if (classA > maxClass) maxClass = classA;
while (parser.getString(s) != parser.end) ; // skip the rest
}
// ... allocate:
eTerm.resize(minClass, maxClass);
mTerm.resize(minClass, maxClass);
pTerm.resize(minClass, maxClass);
cTerm.resize(minClass, maxClass);
sTerm.resize(minClass, maxClass);
FORALL(eTerm,i) eTerm[i]=mTerm[i]=pTerm[i]=cTerm[i]=sTerm[i]=False;
eVal.resize(minClass, maxClass);
mVal.resize(minClass, maxClass);
pVal.resize(minClass, maxClass);
cVal.resize(minClass, maxClass);
sVal.resize(minClass, maxClass);
FORALL(eVal,i) eVal[i] = new Matrix<Real>(dim,dim);
FORALL(mVal,i) mVal[i] = new Matrix<Real>(dim,dim);
FORALL(pVal,i) pVal[i] = new Matrix<Real>(dim,dim);
FORALL(cVal,i) cVal[i] = new Vector<Real>(dim);
FORALL(eVal,i)
{
FORALL_ROWS(*eVal[i],j)
FORALL_COLUMNS(*eVal[i],k)
{
(*eVal[i])(j,k) = 0.0;
(*mVal[i])(j,k) = 0.0;
(*pVal[i])(j,k) = 0.0;
}
}
FORALL(cVal,i) { FORALL(*cVal[i],j) (*cVal[i])[j] = 0.0; }
FORALL(sVal,i) sVal[i] = 0.0;
// ... and read:
parser.fseek(position);
parser.getString(s); // skip iso / anisotropic
while (parser.getValue(&classA) == parser.ok)
{
while (parser.getString(s) != parser.end) // read entries E, M, C, ...
{
strToLower(s);
if (s[0]=='m' || s[0]=='e' || s[0]=='p')
{
if (s[0]=='m') { mat = mVal[classA]; mTerm[classA] = True; }
else if (s[0]=='e') { mat = eVal[classA]; eTerm[classA] = True; }
else if (s[0]=='p') { mat = pVal[classA]; pTerm[classA] = True; }
for (i=1; i<=dim; ++i)
for (k=1; k<=dim; ++k)
{
parser.getValue(&val);
(*mat)(i,k) = val;
}
}
else if (s[0]=='c') // convection term
{
for (k=1; k<=dim; ++k)
{
parser.getValue(&val);
(*cVal[classA])[k] = val;
}
cTerm[classA] = True;
}
else if (s[0]=='s') // source term
{
parser.getValue(&val);
sVal[classA] = val;
sTerm[classA] = True;
}
else {
cout <<"\n*** Materials: unknown key word for material: "
<< s << "\n";
cout.flush(); abort();
}
}
}
}
//----------------------------------------------------------------------------
void DefaultMaterial:: readBCValues(Parser& parser)
{
int i, classA, maxClass = 0, minClass = 1;
Real val;
char s[MaxLineLength], BCType;
//Real huge = machMax(Real(0.0));
FORALL(neumannVal,i) neumannVal[i] = 0.0;
FORALL(cauchyVal, i) cauchyVal[i] = 0.0;
FORALL(neumannVal,i) neumannTerm[i] = cauchyTerm[i] = False;
parser.rewind();
if (parser.findWord("Bound") != parser.ok)
{
cout <<"\n*** No boundary conditions encountered in " << s << "\n";
cout.flush(); abort();
}
long position = parser.ftell();
// scan the entries:
while (parser.getValue(&classA) == parser.ok)
{
parser.getValue(&BCType);
BCType = toupper(BCType);
parser.getValue(&val);
if (BCType == 'D') continue; // not treated in materials
if (classA <= 0)
{
cout << "\n*** Boundary condition class must be > 0 (read: "
<< classA << ")\n";
cout.flush(); abort();
}
if (classA > maxClass) maxClass = classA;
if (BCType == 'C') parser.getValue(&val); // one additional term
// for Cauchy BCs
}
if (maxClass == 0) return;
// ... allocate and read:
neumannVal.resize (minClass, maxClass);
cauchyVal.resize (minClass, maxClass);
neumannTerm.resize(minClass, maxClass);
cauchyTerm.resize (minClass, maxClass);
FORALL(neumannVal,i) neumannVal[i] = 0.0;
FORALL(cauchyVal, i) cauchyVal[i] = 0.0;
FORALL(neumannVal,i) neumannTerm[i] = cauchyTerm[i] = False;
parser.fseek(position);
while (parser.getValue(&classA) == parser.ok)
{
parser.getValue(&BCType);
BCType = toupper(BCType);
parser.getValue(&val);
if (BCType == 'D') { } // will be read by DirichletBCs
else if (BCType == 'N')
{
if (neumannTerm[classA]) doubleID(classA);
neumannVal [classA] = val;
neumannTerm[classA] = True;
}
else if (BCType == 'C')
{
if (cauchyTerm [classA]) doubleID(classA);
if (neumannTerm[classA]) doubleID(classA);
cauchyVal [classA] = val;
cauchyTerm[classA] = True;
parser.getValue(&val);
neumannVal [classA] = val;
neumannTerm[classA] = True;
}
}
}
//-------------------------------------------------------------------------
void DefaultMaterial:: readFactors(Parser& parser, Real* eFact, Real* mFact, Real* pFact,
Real* cFact, Real* sFact,
Real* neumannFact, Real* cauchyFact)
{
char s[MaxLineLength];
Real fact;
parser.rewind();
if (parser.findWord("Fac") != parser.ok) return;
while (parser.getString(s) == parser.ok)
{
parser.getValue(&fact);
strToLower(s);
if (s[0]=='e') *eFact = fact;
else if (s[0]=='m') *mFact = fact;
else if (s[0]=='p') *pFact = fact;
else if (s[0]=='c')
{
if (s[1]=='o') *cFact = fact; // convection term
else *cauchyFact = fact;
}
else if (s[0]=='s') *sFact = fact;
else if (s[0]=='n') *neumannFact = fact;
else if (s[0]=='d') ;
else {
cout << "\n*** Class Material: unknown identifier for factor: "
<< s << "\n";
cout.flush(); abort();
}
}
}
//-------------------------------------------------------------------------
Real DefaultMaterial:: Neumann(int type, Vector<Real>* /*x*/, Real /*time*/)
{
if (!neumannTerm[type]) unknownID(type);
return neumannVal[type];
}
Real DefaultMaterial:: Inner(int /*type*/, Vector<Real>* x, Real /*time*/)
{
Vector<Real>& xx = *x;
switch(spaceDim)
{
case 2:
if (!eTerm[1]) unknownID(1);
if (!eTerm[2]) unknownID(2);
return (eVal[1]-eVal[2])*(xx[1] + xx[2]);
case 3:
if (!eTerm[1]) unknownID(1);
if (!eTerm[2]) unknownID(2);
return (eVal[1]-eVal[2])*(xx[1] + xx[2] + xx[3]);
default: cout.flush(); abort(); return 0;
}
}
Real DefaultMaterial:: Cauchy(int type, Vector<Real>* /*x*/, Real /*time*/)
{
if (!cauchyTerm[type]) unknownID(type);
return cauchyVal[type];
}
//-------------------------------------------------------------------------
VarSourceMaterial:: VarSourceMaterial(const char* fileName, int spaceDim0)
: DefaultMaterial(fileName, spaceDim0)
{
setMaterialType();
}
//-------------------------------------------------------------------------
VarEllipticMaterial:: VarEllipticMaterial(const char* fileName, int spaceDim0)
: DefaultMaterial(fileName, spaceDim0)
{
setMaterialType();
}
//-------------------------------------------------------------------------
VarPoissonMaterial:: VarPoissonMaterial(const char* fileName, int spaceDim0)
: DefaultMaterial(fileName, spaceDim0)
{
setMaterialType();
}
//-------------------------------------------------------------------------
LinMaterial:: LinMaterial(const char* fileName, int spaceDim0)
: DefaultMaterial(fileName, spaceDim0)
{
setMaterialType();
}
Real LinMaterial:: S(int /*type*/, Vector<Real>* x, Real /*time*/)
{
Vector<Real>& xx = *x;
switch(spaceDim)
{
case 1: return xx[1];
case 2: return xx[1] + xx[2];
case 3: return 0.0;
default: cout.flush(); abort(); return 0;
}
}
Real LinMaterial:: M(int /*type*/, Vector<Real>* x, int /*i*/, int /*j*/)
{
Vector<Real>& xx = *x;
switch(spaceDim)
{
case 1: return xx[1];
case 2: return 1.0;//xx[1] + xx[2];
case 3: return 0.0;
default: cout.flush(); abort(); return 0;
}
}
Real LinMaterial:: Neumann(int /*type*/, Vector<Real>* x, Real /*time*/)
{
Vector<Real>& xx = *x;
switch(spaceDim)
{
case 1: return 1.0 + xx[1];
case 2: return 1.0 + xx[1]+xx[2];
case 3: return -1.0 + xx[1]+xx[2]+xx[3];
default: cout.flush(); abort(); return 0;
}
}
Real LinMaterial:: Cauchy(int /*type*/, Vector<Real>* /*x*/, Real /*time*/)
{
// Vector<Real>& xx = *x;
switch(spaceDim)
{
case 1: return 1.0;
case 2: return 1.0;
case 3: return 1.0;
default: cout.flush(); abort(); return 0;
}
}
Num LinMaterial:: trueSolInPoint(const Vector<Real>& x, const Real /*time*/)
{
switch(spaceDim)
{
case 1: return x[1];
case 2: return x[1]+x[2];
case 3: return x[1]+x[2]+x[3];
default: cout.flush(); abort(); return 0;
}
}
//-------------------------------------------------------------------------
PeakPoisson:: PeakPoisson(const char* fileName, int spaceDim0)
: VarPoissonMaterial(fileName, spaceDim0)
{
setMaterialType();
}
Real PeakPoisson:: S(int /*type*/, Vector<Real>* x, Real /*time*/)
{
Vector<Real>& xx = *x;
switch(spaceDim)
{
case 1: return S1d(xx[1]);
case 2: return S2d(xx[1], xx[2]);
case 3: return S3d(xx[1], xx[2], xx[3]);
default: cout.flush(); abort(); return 0;
}
}
Real PeakPoisson:: S1d(Real x)
{
Real xm=x-0.5;
Real f, fx, fxx;
Real e, ex, exx;
Real uxx, p=100.0;
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;
uxx = fxx*e + 2.0*fx*ex + f*exx;
return uxx;
}
Real PeakPoisson:: S2d(Real x, Real y)
{
Real xm=x-0.5, ym=y-0.5;
Real f, fx, fy, fxx, fyy;
Real e, ex, exx, ey, eyy;
Real uxx, uyy, p=100.0;
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 PeakPoisson:: S3d(Real x, Real y, Real z)
{
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;
}
//-------------------------------------------------------------------------
PeakSource:: PeakSource(const char* fileName, int spaceDim0)
: VarSourceMaterial(fileName, spaceDim0)
{
setMaterialType();
}
Real PeakSource:: S(int /*type*/, Vector<Real>* x, Real /*time*/)
{
Vector<Real>& xx = *x;
switch(spaceDim)
{
case 1: return S1d(xx[1]);
case 2: return S2d(xx[1], xx[2]);
case 3: return S3d(xx[1], xx[2], xx[3]);
default: cout.flush(); abort(); return 0;
}
}
Real PeakSource:: S1d(Real x)
{
Real xm=x-0.5;
Real f, fx, fxx;
Real e, ex, exx;
Real uxx, p=100.0;
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;
uxx = fxx*e + 2.0*fx*ex + f*exx;
return uxx;
// return exp(-p*(xm*xm));
}
Real PeakSource:: S2d(Real x, Real y)
{
Real xm=x-0.5, ym=y-0.5;
Real f, fx, fy, fxx, fyy;
Real e, ex, exx, ey, eyy;
Real uxx, uyy, p=100.0;
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);
// return exp(-p*(xm*xm+ym*ym));
}
Real PeakSource:: S3d(Real x, Real y, Real z)
{
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;
//p = -100.; return exp(p*(xma*xma + ymb*ymb + zmc*zmc));
}
Num PeakSource:: trueSolInPoint(const Vector<Real>& x, const Real /*time*/)
{
Real val = 0.0;
switch(spaceDim)
{
case 1: // peak-1d.cmd
{
val = x[1]*(x[1]-1.0)*exp(-100.0*(x[1]-0.5)*(x[1]-0.5));
return -val;
}
case 2: // peak-2d.cmd
{
Real xm=x[1]-0.5, ym=x[2]-0.5;
Real preF= x[1]*(x[1]-1.0)*x[2]*(x[2]-1.0);
val = preF*exp(-100.0*(xm*xm+ym*ym));
return val;
}
case 3: // peak-3d.cmd
{
Real xm=x[1]-0.5, ym=x[2]-0.5, zm=x[3]-0.5;
Real preF= x[1]*(x[1]-1.0)*x[2]*(x[2]-1.0)*x[3]*(x[3]-1.0);
val = preF*exp(-150.0*(xm*xm+ym*ym+zm*zm));
return 100.0*val;
}
default: cout.flush(); abort(); return 0;
}
}
//-------------------------------------------------------------------------
MultiPeakSource:: MultiPeakSource(const char* fileName, int spaceDim0)
: VarSourceMaterial(fileName, spaceDim0)
{
setMaterialType();
}
Real MultiPeakSource:: S(int /*type*/, Vector<Real>* x, Real /*time*/)
{
Vector<Real>& xx = *x;
switch(spaceDim)
{
case 1: return S1d(xx[1]);
case 2: return S2d(xx[1], xx[2]);
case 3: return S3d(xx[1], xx[2], xx[3]);
default: cout.flush(); abort(); return 0;
}
}
Real MultiPeakSource:: S1d(Real x)
{
Real xm;
Real f, fx, fxx;
Real e, ex, exx;
Real uxx, p, factor;
// peak 1
xm=x-0.5;
factor = 1.0;
p = 100.0;
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;
uxx = factor*(fxx*e + 2.0*fx*ex + f*exx);
// peak 2
xm=x-0.2;
factor = 1.0;
p = 150.0;
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;
uxx += factor*(fxx*e + 2.0*fx*ex + f*exx);
// peak 3
xm=x-0.9;
factor = 1.0;
p = 120.0;
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;
uxx += factor*(fxx*e + 2.0*fx*ex + f*exx);
return uxx;
}
Real MultiPeakSource:: S2d(Real x, Real y)
{
Real xm, ym;
Real f, fx, fy, fxx, fyy;
Real e, ex, exx, ey, eyy;
Real uxx, uyy, p, factor, sum=0.0;
// peak 1
xm=x-0.5;
ym=y-0.5;
p=100.0;
factor = 100.0;
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;
sum = factor*(uxx + uyy);
// peak 2
xm=x-0.9;
ym=y-0.5;
p=80.0;
factor = 50.0;
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;
sum += factor*(uxx + uyy);
// peak 3
xm = x-0.2;
ym = y-0.5;
p = 150.0;
factor = 50.0;
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;
sum += factor*(uxx + uyy);
return sum;
}
Real MultiPeakSource:: S3d(Real x, Real y, Real z)
{
Real a, b, c, p, factor;
Real xm1,ym1,zm1,xm1x,ym1y,zm1z,xma,ymb,zmc;
Real sum,divisor, fVal;
a = 0.5;
b = 0.5;
c = 0.5;
p = 150.0;
factor = 100.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);
fVal = -factor*sum/divisor; /* g peak no.1 */
a = 0.7;
b = 0.6;
c = 0.5;
p = 50.0;
factor = 180.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);
fVal += -factor*sum/divisor; /* g peak no.2 */
a = 0.3;
b = 0.6;
c = 0.5;
p = 150.0;
factor = 120.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);
fVal += -factor*sum/divisor; /* g peak no.3 */
return fVal;
}
//-------------------------------------------------------------------------
UserStaticMaterial:: UserStaticMaterial(const char* /*fileName*/, int spaceDim0)
: DefaultMaterial(spaceDim0)
{
setMaterialType();
}
//-------------------------------------------------------------------------
CylindricCoord:: CylindricCoord(const char* fileName, int spaceDim0)
: UserStaticMaterial(fileName,spaceDim0)
{ }
//-------------------------------------------------------------------------
UserVarSource:: UserVarSource(const char* fileName, int spaceDim0)
: VarSourceMaterial(fileName, spaceDim0)
{
setMaterialType();
}
Real UserVarSource:: S(int /*type*/, Vector<Real>* x, Real /*time*/)
{
Vector<Real>& xx = *x;
return xx[1]*(xx[1]-1.0)*xx[2]*(xx[2]-1.0)*1000.0;
}
//-------------------------------------------------------------------------
syntax highlighted by Code2HTML, v. 0.9.1