/* $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* 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(dim,dim); FORALL(mVal,i) mVal[i] = new Matrix(dim,dim); FORALL(pVal,i) pVal[i] = new Matrix(dim,dim); FORALL(cVal,i) cVal[i] = new Vector(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* /*x*/, Real /*time*/) { if (!neumannTerm[type]) unknownID(type); return neumannVal[type]; } Real DefaultMaterial:: Inner(int /*type*/, Vector* x, Real /*time*/) { Vector& 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* /*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* x, Real /*time*/) { Vector& 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* x, int /*i*/, int /*j*/) { Vector& 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* x, Real /*time*/) { Vector& 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* /*x*/, Real /*time*/) { // Vector& 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& 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* x, Real /*time*/) { Vector& 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* x, Real /*time*/) { Vector& 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& 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* x, Real /*time*/) { Vector& 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* x, Real /*time*/) { Vector& xx = *x; return xx[1]*(xx[1]-1.0)*xx[2]*(xx[2]-1.0)*1000.0; } //-------------------------------------------------------------------------