Galerkin FEM for elliptic PDEs
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

FEM4EllipticPDE.h 9.8KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273
  1. // ===========================================================================
  2. //
  3. // Filename: FEM4EllipticPDE.h
  4. //
  5. // Created: 08/16/12 18:19:41
  6. // Compiler: Tested with g++, icpc, and MSVC 2010
  7. //
  8. // Author: Trevor Irons (ti)
  9. //
  10. // Organisation: Colorado School of Mines (CSM)
  11. // United States Geological Survey (USGS)
  12. //
  13. // Email: tirons@mines.edu, tirons@usgs.gov
  14. //
  15. // This program is free software: you can redistribute it and/or modify
  16. // it under the terms of the GNU General Public License as published by
  17. // the Free Software Foundation, either version 3 of the License, or
  18. // (at your option) any later version.
  19. //
  20. // This program is distributed in the hope that it will be useful,
  21. // but WITHOUT ANY WARRANTY; without even the implied warranty of
  22. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  23. // GNU General Public License for more details.
  24. //
  25. // You should have received a copy of the GNU General Public License
  26. // along with this program. If not, see <http://www.gnu.org/licenses/>.
  27. //
  28. // ===========================================================================
  29. /**
  30. @file
  31. @author Trevor Irons
  32. @date 08/16/12
  33. @version $Rev$
  34. **/
  35. #ifndef FEM4ELLIPTICPDE_INC
  36. #define FEM4ELLIPTICPDE_INC
  37. #include "LemmaObject.h"
  38. #include "vtkImplicitFunction.h"
  39. #include "vtkDataSet.h"
  40. #include "vtkXMLDataSetWriter.h"
  41. #include "vtkSmartPointer.h"
  42. #include "vtkIdList.h"
  43. #include "vtkCell.h"
  44. #include "vtkDoubleArray.h"
  45. #include "vtkPointData.h"
  46. #include "vtkDataSetSurfaceFilter.h"
  47. #include "vtkCellArray.h"
  48. #include "vtkCellData.h"
  49. #include <Eigen/IterativeLinearSolvers>
  50. #include "DCSurvey.h"
  51. namespace Lemma {
  52. /** @defgroup FEM4EllipticPDE
  53. @brief Finite element solver for elliptic PDE's
  54. @details Uses Galerkin method.
  55. */
  56. /** \addtogroup FEM4EllipticPDE
  57. @{
  58. */
  59. // ===================================================================
  60. // Class: FEM4EllipticPDE
  61. /**
  62. @class FEM4EllipticPDE
  63. @brief Galerkin FEM solver for Elliptic PDE problems
  64. @details Solves problems of the type
  65. \f[ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
  66. \f]
  67. */
  68. class FEM4EllipticPDE : public LemmaObject {
  69. friend std::ostream &operator<<(std::ostream &stream,
  70. const FEM4EllipticPDE &ob);
  71. public:
  72. // ==================== LIFECYCLE =======================
  73. /** Returns a pointer to a new object of type FEM4EllipticPDE.
  74. * It allocates all necessary memory.
  75. */
  76. static FEM4EllipticPDE* New();
  77. /** Deletes this object. Delete also disconnects any
  78. * attachments to this object.
  79. */
  80. void Delete();
  81. // ==================== OPERATORS =======================
  82. // ==================== OPERATIONS =======================
  83. /** Sets the vtkImplictFunction that is used as sigma in
  84. \f$ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
  85. \f$. If this is not set, sigma evaluates as 1.
  86. */
  87. void SetSigmaFunction(vtkImplicitFunction* sigma);
  88. /** Sets the step to be used in applying boundary conditions.
  89. */
  90. void SetBoundaryStep(const Real& h);
  91. /** Sets the vtkImplictFunction that is used
  92. */
  93. void SetGFunction(vtkImplicitFunction* G);
  94. /** Uses a function pointer instead of a vtkImplicitFunction. For functions that can
  95. * be written in closed form, this can be much faster and more accurate then interpolating
  96. * off a grid.
  97. */
  98. void SetGFunction( Real (*pFcn3)(const Real&, const Real&, const Real&) );
  99. /** Sets the vtkDataSet that defines the calculation grid.
  100. */
  101. void SetGrid(vtkDataSet* Grid);
  102. /** Sets up a DC problem with a survey
  103. * @param[in] ij is the injection index
  104. */
  105. void SetupDC(DCSurvey* Survey, const int& ij);
  106. /** Sets up the potential problem from the VTK file
  107. */
  108. void SetupPotential();
  109. /** Performs solution */
  110. void Solve(const std::string& fname);
  111. /** Performs solution */
  112. void SolveOLD(const std::string& fname);
  113. // ==================== ACCESS =======================
  114. // ==================== INQUIRY =======================
  115. protected:
  116. // ==================== LIFECYCLE =======================
  117. /// Default protected constructor.
  118. FEM4EllipticPDE (const std::string& cname);
  119. /** Default protected constructor. */
  120. ~FEM4EllipticPDE ();
  121. /**
  122. * @copybrief LemmaObject::Release()
  123. * @copydetails LemmaObject::Release()
  124. */
  125. void Release();
  126. // ==================== OPERATIONS =======================
  127. /** Used internally to construct stiffness matrix, A
  128. */
  129. void ConstructAMatrix();
  130. /* Used internally to construct the load vector, g
  131. */
  132. //void ConstructLoadVector();
  133. // This function is copyright (C) 2012 Joseph Capriotti
  134. /**
  135. Returns a vtkIdList of points in member vtkGrid connected to point ido.
  136. @param[in] ido is the point of interest
  137. @param[in] grid is a pointer to a vtkDataSet
  138. */
  139. vtkSmartPointer<vtkIdList> GetConnectedPoints(const int& id0);
  140. /** Uses Simpson's rule to numerically integrate
  141. * (Shamelessly adapted from http://en.wikipedia.org/wiki/Simpson's_rule
  142. */
  143. Real CompositeSimpsons(vtkImplicitFunction* f, Real r0[3], Real r1[3], int numInt);
  144. /** Uses Simpson's rule to numerically integrate
  145. * (Shamelessly adapted from http://en.wikipedia.org/wiki/Simpson's_rule
  146. * this method is for a static f
  147. */
  148. Real CompositeSimpsons(const Real& f, Real r0[3], Real r1[3], int numInt);
  149. /**
  150. * Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
  151. * r0 and r1
  152. */
  153. Real CompositeSimpsons2(vtkImplicitFunction* f, Real r0[3], Real r1[3], int numInt);
  154. /**
  155. * Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
  156. * r0 and r1, uses the Hat function and the function pointer fPtr
  157. */
  158. Real CompositeSimpsons2( Real (*fPtr)(const Real&, const Real&, const Real&), Real r0[3], Real r1[3], int numInt);
  159. /**
  160. * Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
  161. * r0 and r1, uses the Hat function. Assumes a constand function value f
  162. */
  163. Real CompositeSimpsons2( const Real& f, Real r0[3], Real r1[3], int numInt);
  164. /**
  165. * Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
  166. * r0 and r1, uses the Hat function. Assumes a linear function from f0 to f1.
  167. */
  168. Real CompositeSimpsons2( const Real& f0, const Real& f1, Real r0[3], Real r1[3], int numInt);
  169. /**
  170. * Uses Simpson's rule to numerically integrate three functions in 1 dimension over the points
  171. * r0 and r1, uses two Hat functions. Assumes a linear function from f0 to f1.
  172. */
  173. Real CompositeSimpsons3( const Real& f0, const Real& f1, Real r0[3], Real r1[3], int numInt);
  174. /** Hat function for use in Galerkin FEM method, actually a 1/2 Hat.
  175. */
  176. //Real hat(Real r[3], Real ri[3], Real rm1[3] );
  177. Real Hat(const Vector3r &r, const Vector3r& ri, const Vector3r& rm1 );
  178. /** Computes distance between r0 and r1
  179. */
  180. Real dist(Real r0[3], Real r1[3]);
  181. /** Computes distance between r0 and r1
  182. */
  183. Real dist(const Vector3r& r0, const Vector3r& r1);
  184. // ==================== DATA MEMBERS =========================
  185. /** The effective step at the boundary condition. Defaults to 1
  186. */
  187. Real BndryH;
  188. /** The sigma value at the boundary condition. Defaults to 1
  189. */
  190. Real BndrySigma;
  191. /** Implicit function defining the \f$\sigma\f$ term in
  192. \f$ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
  193. \f$.
  194. If set to 1, the Poisson equation is solved. Any type of vtkImplicitFunction may
  195. be used, including vtkImplicitDataSet.
  196. */
  197. vtkSmartPointer<vtkImplicitFunction> vtkSigma;
  198. /** Implicit function defining the \f$g\f$ term in
  199. \f$ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
  200. \f$.
  201. Typically this is used to define the source term.
  202. */
  203. vtkSmartPointer<vtkImplicitFunction> vtkG;
  204. /** Any type of vtkDataSet may be used as a calculation Grid */
  205. vtkSmartPointer<vtkDataSet> vtkGrid;
  206. /** Function pointer to function describing g */
  207. Real (*gFcn3)(const Real&, const Real&, const Real&);
  208. Eigen::SparseMatrix<Real> A;
  209. VectorXr g;
  210. private:
  211. }; // ----- end of class FEM4EllipticPDE -----
  212. /*! @} */
  213. } // ----- end of Lemma name -----
  214. #endif // ----- #ifndef FEM4ELLIPTICPDE_INC -----