NGSolve
4.9
|
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