Galerkin FEM for elliptic PDEs
Вы не можете выбрать более 25 тем Темы должны начинаться с буквы или цифры, могут содержать дефисы(-) и должны содержать не более 35 символов.

FEM4EllipticPDE.h 9.7KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269
  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. /** Performs solution */
  107. void Solve(const std::string& fname);
  108. /** Performs solution */
  109. void SolveOLD(const std::string& fname);
  110. // ==================== ACCESS =======================
  111. // ==================== INQUIRY =======================
  112. protected:
  113. // ==================== LIFECYCLE =======================
  114. /// Default protected constructor.
  115. FEM4EllipticPDE (const std::string& cname);
  116. /** Default protected constructor. */
  117. ~FEM4EllipticPDE ();
  118. /**
  119. * @copybrief LemmaObject::Release()
  120. * @copydetails LemmaObject::Release()
  121. */
  122. void Release();
  123. // ==================== OPERATIONS =======================
  124. /** Used internally to construct stiffness matrix, A
  125. */
  126. void ConstructAMatrix();
  127. /* Used internally to construct the load vector, g
  128. */
  129. //void ConstructLoadVector();
  130. // This function is copyright (C) 2012 Joseph Capriotti
  131. /**
  132. Returns a vtkIdList of points in member vtkGrid connected to point ido.
  133. @param[in] ido is the point of interest
  134. @param[in] grid is a pointer to a vtkDataSet
  135. */
  136. vtkSmartPointer<vtkIdList> GetConnectedPoints(const int& id0);
  137. /** Uses Simpson's rule to numerically integrate
  138. * (Shamelessly adapted from http://en.wikipedia.org/wiki/Simpson's_rule
  139. */
  140. Real CompositeSimpsons(vtkImplicitFunction* f, Real r0[3], Real r1[3], int numInt);
  141. /** Uses Simpson's rule to numerically integrate
  142. * (Shamelessly adapted from http://en.wikipedia.org/wiki/Simpson's_rule
  143. * this method is for a static f
  144. */
  145. Real CompositeSimpsons(const Real& f, Real r0[3], Real r1[3], int numInt);
  146. /**
  147. * Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
  148. * r0 and r1
  149. */
  150. Real CompositeSimpsons2(vtkImplicitFunction* f, Real r0[3], Real r1[3], int numInt);
  151. /**
  152. * Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
  153. * r0 and r1, uses the Hat function and the function pointer fPtr
  154. */
  155. Real CompositeSimpsons2( Real (*fPtr)(const Real&, const Real&, const Real&), Real r0[3], Real r1[3], int numInt);
  156. /**
  157. * Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
  158. * r0 and r1, uses the Hat function. Assumes a constand function value f
  159. */
  160. Real CompositeSimpsons2( const Real& f, Real r0[3], Real r1[3], int numInt);
  161. /**
  162. * Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
  163. * r0 and r1, uses the Hat function. Assumes a linear function from f0 to f1.
  164. */
  165. Real CompositeSimpsons2( const Real& f0, const Real& f1, Real r0[3], Real r1[3], int numInt);
  166. /**
  167. * Uses Simpson's rule to numerically integrate three functions in 1 dimension over the points
  168. * r0 and r1, uses two Hat functions. Assumes a linear function from f0 to f1.
  169. */
  170. Real CompositeSimpsons3( const Real& f0, const Real& f1, Real r0[3], Real r1[3], int numInt);
  171. /** Hat function for use in Galerkin FEM method, actually a 1/2 Hat.
  172. */
  173. //Real hat(Real r[3], Real ri[3], Real rm1[3] );
  174. Real Hat(const Vector3r &r, const Vector3r& ri, const Vector3r& rm1 );
  175. /** Computes distance between r0 and r1
  176. */
  177. Real dist(Real r0[3], Real r1[3]);
  178. /** Computes distance between r0 and r1
  179. */
  180. Real dist(const Vector3r& r0, const Vector3r& r1);
  181. // ==================== DATA MEMBERS =========================
  182. /** The effective step at the boundary condition. Defaults to 1
  183. */
  184. Real BndryH;
  185. /** The sigma value at the boundary condition. Defaults to 1
  186. */
  187. Real BndrySigma;
  188. /** Implicit function defining the \f$\sigma\f$ term in
  189. \f$ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
  190. \f$.
  191. If set to 1, the Poisson equation is solved. Any type of vtkImplicitFunction may
  192. be used, including vtkImplicitDataSet.
  193. */
  194. vtkSmartPointer<vtkImplicitFunction> vtkSigma;
  195. /** Implicit function defining the \f$g\f$ term in
  196. \f$ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
  197. \f$.
  198. Typically this is used to define the source term.
  199. */
  200. vtkSmartPointer<vtkImplicitFunction> vtkG;
  201. /** Any type of vtkDataSet may be used as a calculation Grid */
  202. vtkSmartPointer<vtkDataSet> vtkGrid;
  203. /** Function pointer to function describing g */
  204. Real (*gFcn3)(const Real&, const Real&, const Real&);
  205. Eigen::SparseMatrix<Real> A;
  206. VectorXr g;
  207. private:
  208. }; // ----- end of class FEM4EllipticPDE -----
  209. /*! @} */
  210. } // ----- end of Lemma name -----
  211. #endif // ----- #ifndef FEM4ELLIPTICPDE_INC -----