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 11KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309
  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 <Eigen/SparseLU>
  51. #include "DCSurvey.h"
  52. namespace Lemma {
  53. /** @defgroup FEM4EllipticPDE
  54. @brief Finite element solver for elliptic PDE's
  55. @details Uses Galerkin method.
  56. */
  57. /** \addtogroup FEM4EllipticPDE
  58. @{
  59. */
  60. // ===================================================================
  61. // Class: FEM4EllipticPDE
  62. /**
  63. @class FEM4EllipticPDE
  64. @brief Galerkin FEM solver for Elliptic PDE problems
  65. @details Solves problems of the type
  66. \f[ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
  67. \f]
  68. */
  69. class FEM4EllipticPDE : public LemmaObject {
  70. friend std::ostream &operator<<(std::ostream &stream,
  71. const FEM4EllipticPDE &ob);
  72. public:
  73. // ==================== LIFECYCLE =======================
  74. /** Returns a pointer to a new object of type FEM4EllipticPDE.
  75. * It allocates all necessary memory.
  76. */
  77. static FEM4EllipticPDE* New();
  78. /** Deletes this object. Delete also disconnects any
  79. * attachments to this object.
  80. */
  81. void Delete();
  82. #ifdef HAVE_YAMLCPP
  83. /**
  84. * Uses YAML to serialize this object.
  85. * @return a YAML::Node
  86. */
  87. YAML::Node Serialize() const;
  88. /**
  89. * Constructs an object from a YAML::Node.
  90. */
  91. static FEM4EllipticPDE* DeSerialize(const YAML::Node& node);
  92. #endif
  93. // ==================== OPERATORS =======================
  94. // ==================== OPERATIONS =======================
  95. /** Sets the vtkImplictFunction that is used as sigma in
  96. \f$ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
  97. \f$. If this is not set, sigma evaluates as 1.
  98. */
  99. void SetSigmaFunction(vtkImplicitFunction* sigma);
  100. /** Sets the step to be used in applying boundary conditions.
  101. */
  102. void SetBoundaryStep(const Real& h);
  103. /** Sets the vtkImplictFunction that is used
  104. */
  105. void SetGFunction(vtkImplicitFunction* G);
  106. /** Uses a function pointer instead of a vtkImplicitFunction. For functions that can
  107. * be written in closed form, this can be much faster and more accurate then interpolating
  108. * off a grid.
  109. */
  110. void SetGFunction( Real (*pFcn3)(const Real&, const Real&, const Real&) );
  111. /** Sets the vtkDataSet that defines the calculation grid.
  112. */
  113. void SetGrid(vtkDataSet* Grid);
  114. /** Sets up a DC problem with a survey
  115. * @param[in] ij is the injection index
  116. */
  117. void SetupDC(DCSurvey* Survey, const int& ij);
  118. /** Sets up the potential problem from the VTK file
  119. */
  120. void SetupPotential();
  121. /** Performs solution */
  122. void Solve(const std::string& fname);
  123. /** Performs solution */
  124. void SolveOLD2(const std::string& fname);
  125. // ==================== ACCESS =======================
  126. // ==================== INQUIRY =======================
  127. protected:
  128. // ==================== LIFECYCLE =======================
  129. /// Default protected constructor.
  130. FEM4EllipticPDE (const std::string& cname);
  131. #ifdef HAVE_YAMLCPP
  132. /** Protected DeDerializing constructor, use factory DeSerialize method*/
  133. FEM4EllipticPDE (const YAML::Node& node);
  134. #endif
  135. /** Default protected constructor. */
  136. ~FEM4EllipticPDE ();
  137. /**
  138. * @copybrief LemmaObject::Release()
  139. * @copydetails LemmaObject::Release()
  140. */
  141. void Release();
  142. // ==================== OPERATIONS =======================
  143. /** Used internally to construct stiffness matrix, A
  144. */
  145. void ConstructAMatrix();
  146. /* Used internally to construct the load vector, g
  147. */
  148. //void ConstructLoadVector();
  149. // This function is copyright (C) 2012 Joseph Capriotti
  150. /**
  151. Returns a vtkIdList of points in member vtkGrid connected to point ido.
  152. @param[in] ido is the point of interest
  153. @param[in] grid is a pointer to a vtkDataSet
  154. */
  155. vtkSmartPointer<vtkIdList> GetConnectedPoints(const int& id0);
  156. /** Uses Simpson's rule to numerically integrate
  157. * (Shamelessly adapted from http://en.wikipedia.org/wiki/Simpson's_rule
  158. */
  159. Real CompositeSimpsons(vtkImplicitFunction* f, Real r0[3], Real r1[3], int numInt);
  160. /** Uses Simpson's rule to numerically integrate
  161. * (Shamelessly adapted from http://en.wikipedia.org/wiki/Simpson's_rule
  162. * this method is for a static f
  163. */
  164. Real CompositeSimpsons(const Real& f, Real r0[3], Real r1[3], int numInt);
  165. /**
  166. * Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
  167. * r0 and r1
  168. */
  169. Real CompositeSimpsons2(vtkImplicitFunction* f, Real r0[3], Real r1[3], int numInt);
  170. /**
  171. * Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
  172. * r0 and r1, uses the Hat function and the function pointer fPtr
  173. */
  174. Real CompositeSimpsons2( Real (*fPtr)(const Real&, const Real&, const Real&), Real r0[3], Real r1[3], int numInt);
  175. /**
  176. * Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
  177. * r0 and r1, uses the Hat function. Assumes a constand function value f
  178. */
  179. Real CompositeSimpsons2( const Real& f, Real r0[3], Real r1[3], int numInt);
  180. /**
  181. * Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
  182. * r0 and r1, uses the Hat function. Assumes a linear function from f0 to f1.
  183. */
  184. Real CompositeSimpsons2( const Real& f0, const Real& f1, Real r0[3], Real r1[3], int numInt);
  185. /**
  186. * Uses Simpson's rule to numerically integrate three functions in 1 dimension over the points
  187. * r0 and r1, uses two Hat functions. Assumes a linear function from f0 to f1.
  188. */
  189. Real CompositeSimpsons3( const Real& f0, const Real& f1, Real r0[3], Real r1[3], int numInt);
  190. /** Hat function for use in Galerkin FEM method, actually a 1/2 Hat.
  191. */
  192. //Real hat(Real r[3], Real ri[3], Real rm1[3] );
  193. Real Hat(const Vector3r &r, const Vector3r& ri, const Vector3r& rm1 );
  194. /** Computes distance between r0 and r1
  195. */
  196. Real dist(Real r0[3], Real r1[3]);
  197. /** Computes distance between r0 and r1
  198. */
  199. Real dist(const Vector3r& r0, const Vector3r& r1);
  200. /** Sets up the potential problem from the VTK file
  201. */
  202. void SetupLineSourcePotential();
  203. /** Sets up the potential problem from the VTK file
  204. */
  205. void SetupSurfaceSourcePotential();
  206. /** Sets up the potential problem from the VTK file
  207. */
  208. void SetupPointSourcePotential();
  209. /** Sets up the potential problem from the VTK file
  210. */
  211. void SetupVolumeSourcePotential();
  212. // ==================== DATA MEMBERS =========================
  213. /** The effective step at the boundary condition. Defaults to 1
  214. */
  215. Real BndryH;
  216. /** The sigma value at the boundary condition. Defaults to 1
  217. */
  218. Real BndrySigma;
  219. /** Implicit function defining the \f$\sigma\f$ term in
  220. \f$ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
  221. \f$.
  222. If set to 1, the Poisson equation is solved. Any type of vtkImplicitFunction may
  223. be used, including vtkImplicitDataSet.
  224. */
  225. vtkSmartPointer<vtkImplicitFunction> vtkSigma;
  226. /** Implicit function defining the \f$g\f$ term in
  227. \f$ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
  228. \f$.
  229. Typically this is used to define the source term.
  230. */
  231. vtkSmartPointer<vtkImplicitFunction> vtkG;
  232. /** Any type of vtkDataSet may be used as a calculation Grid */
  233. vtkSmartPointer<vtkDataSet> vtkGrid;
  234. /** Function pointer to function describing g */
  235. Real (*gFcn3)(const Real&, const Real&, const Real&);
  236. Eigen::SparseMatrix<Real> A;
  237. VectorXr g;
  238. private:
  239. }; // ----- end of class FEM4EllipticPDE -----
  240. /*! @} */
  241. } // ----- end of Lemma name -----
  242. #endif // ----- #ifndef FEM4ELLIPTICPDE_INC -----