righthandsidebyequation.h

Go to the documentation of this file.
00001 #ifndef  __RightHandSideByEquation_h
00002 #define  __RightHandSideByEquation_h
00003 
00004 #include  "domainrighthandside.h"
00005 #include  "exactsolution.h"
00006 #include  "equation.h"
00007 
00008 /*-----------------------------------------*/
00009 
00010 
00011 namespace Gascoigne
00012 {
00013 class RightHandSideByEquation : public DomainRightHandSide
00014 {
00015 protected:
00016 
00017   const Equation*      _EQ;
00018   const ExactSolution* _ES;
00019 
00020 public:
00021 
00022   RightHandSideByEquation(const Equation* eq, const ExactSolution* es)
00023     : DomainRightHandSide(), _EQ(eq), _ES(es) { assert(es); assert(eq);}
00024   
00025   std::string GetName() const { return "RightHandSideByEquation";} 
00026   int GetNcomp() const { return _EQ->GetNcomp();}
00027 
00028   double operator()(int c, const Vertex2d& v)const 
00029     {
00030       int n = _EQ->GetNcomp();
00031       DoubleVector b(n,0.);
00032       FemFunction U(n);
00033       for (int i=0; i<n; i++)
00034         {
00035           U[i].n() = 0.;
00036           U[i].m() = (*_ES)(i,v);
00037           U[i].x() = _ES->x(i,v);
00038           U[i].y() = _ES->y(i,v);
00039           U[i].D() = _ES->xx(i,v)+_ES->yy(i,v);
00040         }
00041       _EQ->OperatorStrong(b,U);
00042       if (GetTimeStep()>0.)
00043         {
00044           double eps = 1.e-6;
00045           DoubleVector ut(n,0.);
00046           double time = GetTime();
00047          _ES->SetTime(time+0.5*eps,GetTimeStep());
00048           for (int i=0; i<n; i++)
00049             {     
00050               ut[i] = (*_ES)(i,v);
00051             }
00052           _ES->SetTime(time-0.5*eps,GetTimeStep());
00053           for (int i=0; i<n; i++)
00054             {     
00055               ut[i] -= (*_ES)(i,v);
00056             }
00057           _ES->SetTime(time,GetTimeStep());
00058           b.add(1./eps,ut);
00059         }
00060       return b[c];
00061     }
00062   double operator()(int c, const Vertex3d& v)const 
00063     {
00064       int n = _EQ->GetNcomp();
00065       DoubleVector b(n,0.);
00066       FemFunction U(n);
00067       for (int i=0; i<n; i++)
00068         {
00069           U[i].n() = 0.;
00070           U[i].m() = (*_ES)(i,v);
00071           U[i].x() = _ES->x(i,v);
00072           U[i].y() = _ES->y(i,v);
00073           U[i].z() = _ES->z(i,v);
00074           U[i].D() = _ES->xx(i,v)+_ES->yy(i,v)+_ES->zz(i,v);
00075         }
00076       _EQ->OperatorStrong(b,U);
00077       if (GetTimeStep()>0.)
00078         {
00079           double eps = 1.e-6;
00080           DoubleVector ut(n,0.);
00081           double time = GetTime();
00082          _ES->SetTime(time+0.5*eps,GetTimeStep());
00083           for (int i=0; i<n; i++)
00084             {     
00085               ut[i] = (*_ES)(i,v);
00086             }
00087           _ES->SetTime(time-0.5*eps,GetTimeStep());
00088           for (int i=0; i<n; i++)
00089             {     
00090               ut[i] -= (*_ES)(i,v);
00091             }
00092           _ES->SetTime(time,GetTimeStep());
00093           b.add(1./eps,ut);
00094         }
00095       return b[c];
00096     }
00097 
00098 };
00099 }
00100 
00101 #endif

Generated on Thu Sep 14 10:34:38 2006 for Gascoigne by  doxygen 1.4.7