NGSolve  4.9
fem/pml.hpp
00001 #ifndef FILE_PML
00002 #define FILE_PML
00003 
00004 /*********************************************************************/
00005 /* File:   pml.hpp                                                   */
00006 /* Author: Joachim Schoeberl                                         */
00007 /* Date:   1. Jul. 2004                                              */
00008 /*********************************************************************/
00009 
00010 namespace ngfem
00011 {
00012   extern Complex alpha;
00013   extern double pml_r;
00014   extern double pml_x;
00015 
00016   extern double pml_xmin[3];
00017   extern double pml_xmax[3]; 
00018 
00019   /*
00020     rect_pml = 0 .... circular pml with radius pml_r
00021     rect_pml = 1 .... square pml on square (-pml_x, pml_x)^d
00022     rect_pml = 2 .... rectangular pml on (pml_xmin, pml_xmax) x (pml_ymin, pml_ymax) x (pml_zmin, pml_zmax)
00023   */
00024   extern int rect_pml;
00025   
00026   extern bool apply_deriv_alpha;
00027 
00028 
00029 
00030   
00031 
00032 
00033   
00034   template <>
00035   MappedIntegrationPoint<2,2,Complex> :: 
00036   MappedIntegrationPoint (const IntegrationPoint & aip,
00037                           const ElementTransformation & aeltrans);
00038   // LocalHeap & lh);
00039 
00040   
00041   template <>
00042   MappedIntegrationPoint<2,2,AutoDiff<1,Complex> > :: 
00043   MappedIntegrationPoint (const IntegrationPoint & aip,
00044                           const ElementTransformation & aeltrans);
00045   // LocalHeap & lh);
00046 
00047 
00048 
00049   extern void SetPMLParameters();
00050 
00051 
00052   template <class DIFFOP, class DMATOP, class FEL = FiniteElement>
00053   class PML_BDBIntegrator : public T_BDBIntegrator<DIFFOP,DMATOP,FEL>
00054   {
00055   public:
00056   
00057     enum { DIM_SPACE   = DIFFOP::DIM_SPACE };
00058     enum { DIM_ELEMENT = DIFFOP::DIM_ELEMENT };
00059     enum { DIM_DMAT    = DIFFOP::DIM_DMAT };
00060     enum { DIM         = DIFFOP::DIM };
00061   
00063     PML_BDBIntegrator (const DMATOP & admat)
00064       : T_BDBIntegrator<DIFFOP,DMATOP,FEL> (admat)
00065     { 
00066       SetPMLParameters();
00067 
00068     }
00069 
00071     virtual ~PML_BDBIntegrator ()
00072     { ; }
00073 
00074 
00076     virtual void
00077     CalcElementMatrix (const FiniteElement & bfel, 
00078                            const ElementTransformation & eltrans, 
00079                            FlatMatrix<double> & elmat,
00080                            LocalHeap & locheap) const
00081     {
00082       throw Exception ("PML cannot generate real matrices");
00083     }
00084 
00086     virtual void
00087     CalcElementMatrix (const FiniteElement & bfel, 
00088                            const ElementTransformation & eltrans, 
00089                            FlatMatrix<Complex> & elmat,
00090                            LocalHeap & locheap) const
00091     {
00092       try
00093         {
00094           const FEL & fel = static_cast<const FEL&> (bfel);
00095           int ndof = fel.GetNDof();
00096 
00097           elmat = 0;
00098         
00099           FlatMatrixFixHeight<DIM_DMAT, Complex> bmat (ndof * DIM, locheap);
00100           FlatMatrixFixHeight<DIM_DMAT, Complex> dbmat (ndof * DIM, locheap);
00101 
00102           Mat<DIM_DMAT,DIM_DMAT, Complex> dmat;
00103 
00104           const IntegrationRule & ir = this->GetIntegrationRule (fel);
00105 
00106           for (int i = 0; i < ir.GetNIP(); i++)
00107             {
00108               HeapReset hr (locheap);
00109 
00110               MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE,Complex> 
00111                 sip(ir[i], eltrans);
00112               MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE,double> 
00113                 sip_real(ir[i], eltrans);
00114 
00115               DIFFOP::GenerateMatrix (fel, sip, bmat, locheap);
00116               this->dmatop.GenerateMatrix (fel, sip_real, dmat, locheap);
00117 
00118               Complex fac = sip.GetJacobiDet() * sip.IP().Weight();
00119 
00120               dmat *= fac;
00121               dbmat = dmat * bmat;
00122               // elmat += Trans (bmat) * dbmat;
00123               if (DMATOP::SYMMETRIC)
00124                 FastMat<DIM_DMAT> (elmat.Height(), &dbmat(0,0), &bmat(0,0), &elmat(0,0));
00125               else
00126                 elmat += Trans (bmat) * dbmat;
00127             } 
00128         }
00129 
00130       catch (Exception & e)
00131         {
00132           e.Append ("in CalcElementMatrix, type = ");
00133           e.Append (typeid(*this).name());
00134           e.Append ("\n");
00135           throw;
00136         }
00137       catch (exception & e)
00138         {
00139           Exception e2(e.what());
00140           e2.Append ("\nin CalcElementMatrix, type = ");
00141           e2.Append (typeid(*this).name());
00142           e2.Append ("\n");
00143           throw e2;
00144         }
00145     }
00146 
00147 
00148 
00149     virtual void 
00150     ApplyElementMatrix (const FiniteElement & bfel, 
00151                         const ElementTransformation & eltrans, 
00152                         const FlatVector<double> & elx, 
00153                         FlatVector<double> & ely,
00154                         void * precomputed,
00155                         LocalHeap & locheap) const
00156     { ; }
00157 
00158 
00159     virtual void 
00160     ApplyElementMatrix (const FiniteElement & bfel, 
00161                         const ElementTransformation & eltrans, 
00162                         const FlatVector<Complex> & elx, 
00163                         FlatVector<Complex> & ely,
00164                         void * precomputed,
00165                         LocalHeap & locheap) const
00166     {
00167       const FEL & fel = static_cast<const FEL&> (bfel);
00168       int ndof = fel.GetNDof ();
00169     
00170       ely = 0;
00171 
00172 
00173       if (!apply_deriv_alpha)
00174         {
00175           Vec<DIM_DMAT,Complex> hv1;
00176           Vec<DIM_DMAT,Complex> hv2;
00177         
00178           FlatVector<Complex> hely (ndof*DIM, locheap);
00179           const IntegrationRule & ir = this->GetIntegrationRule (fel);
00180         
00181         
00182           for (int i = 0; i < ir.GetNIP(); i++)
00183             {
00184               MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE,Complex> 
00185                 sip(ir[i], eltrans);
00186 
00187               MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE,double> 
00188                 sip_real (ir[i], eltrans);
00189             
00190               DIFFOP::Apply (fel, sip, elx, hv1, locheap);
00191               this -> dmatop.Apply (fel, sip_real, hv1, hv2, locheap);
00192               DIFFOP::ApplyTrans (fel, sip, hv2, hely, locheap);
00193             
00194               Complex fac = sip.GetJacobiDet() * sip.IP().Weight();
00195               ely += fac * hely;
00196             }     
00197         }
00198       else
00199         {
00200           Vec<DIM_DMAT, AutoDiff<1,Complex> > hv1;
00201           Vec<DIM_DMAT, AutoDiff<1,Complex> > hv2;
00202         
00203           FlatVector<AutoDiff<1,Complex> > hely (ndof*DIM, locheap);
00204           FlatVector<AutoDiff<1,Complex> > helx (ndof*DIM, locheap);
00205           const IntegrationRule & ir = this->GetIntegrationRule (fel);
00206         
00207           for (int j = 0; j < helx.Size(); j++)
00208             helx(j) = elx(j);
00209         
00210           for (int i = 0; i < ir.GetNIP(); i++)
00211             {
00212               MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE,AutoDiff<1,Complex> >
00213                 sip(ir[i], eltrans);
00214 
00215               MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE,double> 
00216                 sip_real (ir[i], eltrans);
00217             
00218               DIFFOP::Apply (fel, sip, helx, hv1, locheap);
00219               this -> dmatop.Apply (fel, sip_real, hv1, hv2, locheap);
00220               DIFFOP::ApplyTrans (fel, sip, hv2, hely, locheap);
00221             
00222               AutoDiff<1,Complex> fac = sip.GetJacobiDet() * sip.IP().Weight();
00223 
00224               for (int j = 0; j < hely.Size(); j++)
00225                 ely(j) += (fac * hely(j)).DValue(0);
00226             }     
00227         }
00228     }
00229 
00230 
00232     virtual int GetDimension () const { return DIM; }
00234     virtual string Name () const { return "PML-BDB integrator"; }
00235   };
00236 
00237 }
00238 
00239 #endif