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_bhmag.cpp 5.3KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176
  1. // ===========================================================================
  2. //
  3. // Filename: utFEM4EllipticPDE.cpp
  4. //
  5. // Created: 08/16/12 19:49:10
  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 0.0
  34. **/
  35. #include "FEM4EllipticPDE.h"
  36. #include "vtkSphere.h"
  37. #include "vtkRectilinearGrid.h"
  38. #include "vtkFloatArray.h"
  39. #include <vtkIdList.h>
  40. #include <vtkIdTypeArray.h>
  41. #include <vtkCell.h>
  42. #include <vtkTriangleFilter.h>
  43. #include <vtkDataSetSurfaceFilter.h>
  44. #include <vtkExtractCells.h>
  45. #include <vtkGeometryFilter.h>
  46. #include <vtkUnstructuredGrid.h>
  47. #include <vtkExtractEdges.h>
  48. #include <vtkImplicitDataSet.h>
  49. #include <vtkXMLUnstructuredGridReader.h>
  50. // debugging stuff
  51. #include <vtkDataSetMapper.h>
  52. #include <vtkSelectionNode.h>
  53. #include <vtkExtractSelection.h>
  54. #include <vtkSelection.h>
  55. #include <vtkVertexGlyphFilter.h>
  56. #define RENDERTEST
  57. //#define CONNECTTEST
  58. #ifdef RENDERTEST
  59. #include "vtkRectilinearGridGeometryFilter.h"
  60. #include "vtkRectilinearGridOutlineFilter.h"
  61. #include "vtkPolyDataMapper.h"
  62. #include "vtkActor.h"
  63. #include "vtkRenderWindowInteractor.h"
  64. #include "vtkRenderer.h"
  65. #include "vtkRenderWindow.h"
  66. #include "vtkCamera.h"
  67. #include "vtkProperty.h"
  68. #include <vtkDataSetMapper.h>
  69. #include <vtkSelectionNode.h>
  70. #include <vtkSelection.h>
  71. #include <vtkExtractSelection.h>
  72. #include <vtkExtractEdges.h>
  73. #include <vtkVertexGlyphFilter.h>
  74. #include <vtkTriangleFilter.h>
  75. #include <vtkImplicitHalo.h>
  76. #endif
  77. using namespace Lemma;
  78. Real Magnet(const Real& x, const Real& y, const Real& z);
  79. int main(int argc, char**argv) {
  80. Eigen::initParallel();
  81. if (argc < 4) {
  82. std::cout << "usage: ./utFEMEllipticPDE_bhmag <g.vtu> <mesh.vtu> <results.vtu>" << std::endl;
  83. //std::cout << "usage: ./utFEMEllipticPDE_bhmag <mesh.vtu> <results.vtu>" << std::endl;
  84. exit(EXIT_SUCCESS);
  85. }
  86. int nx = 80;
  87. //double dx = .0021875 ; // 160
  88. double dx = .004375 ; // 160
  89. double ox = -.175;
  90. vtkFloatArray *xCoords = vtkFloatArray::New();
  91. for (int i=0; i<nx+1; i++) xCoords->InsertNextValue(ox + i*dx);
  92. int ny = 80;
  93. double dy = .004375;
  94. double oy = -.175;
  95. vtkFloatArray *yCoords = vtkFloatArray::New();
  96. for (int i=0; i<ny+1; i++) yCoords->InsertNextValue(oy + i*dy);
  97. int nz = 343; // 685; // 160
  98. double dz = .004375;
  99. double oz = 9.75;
  100. vtkFloatArray *zCoords = vtkFloatArray::New();
  101. for (int i=0; i<nz+1; i++) zCoords->InsertNextValue(oz + i*dz);
  102. vtkRectilinearGrid *rGrid = vtkRectilinearGrid::New();
  103. rGrid->SetDimensions(nx+1, ny+1, nz+1);
  104. rGrid->SetExtent(0, nx, 0, ny, 0, nz);
  105. rGrid->SetXCoordinates(xCoords);
  106. rGrid->SetYCoordinates(yCoords);
  107. rGrid->SetZCoordinates(zCoords);
  108. ////////////////////////////////////////////
  109. // Define Sigma/mu term
  110. // TODO
  111. ////////////////////////////////////////////
  112. // Define source (G) term
  113. vtkXMLUnstructuredGridReader* GReader = vtkXMLUnstructuredGridReader::New();
  114. std::cout << "Reading G file " << argv[1] << std::endl;
  115. GReader->SetFileName(argv[1]);
  116. GReader->Update();
  117. vtkImplicitDataSet* implG = vtkImplicitDataSet::New();
  118. implG->SetDataSet(GReader->GetOutput());
  119. implG->SetOutValue(0.);
  120. implG->SetOutGradient(0., 0., 0.);
  121. ////////////////////////////////////////////
  122. // Define solution mesh
  123. vtkXMLUnstructuredGridReader* MeshReader = vtkXMLUnstructuredGridReader::New();
  124. std::cout << "Reading Mesh file " << argv[2] << std::endl;
  125. MeshReader->SetFileName(argv[2]);
  126. MeshReader->Update();
  127. ////////////////////////////////////////////
  128. // Solve
  129. FEM4EllipticPDE *Solver = FEM4EllipticPDE::New();
  130. Solver->SetGFunction(implG);
  131. //Solver->SetGFunction(Magnet);
  132. //Solver->SetSigmaFunction(implSigma);
  133. Solver->SetBoundaryStep(.05);
  134. Solver->SetGrid(MeshReader->GetOutput());
  135. //Solver->SetupPotential();
  136. //Solver->SetGrid(rGrid);
  137. Solver->Solve(argv[3]);
  138. // Clean up
  139. Solver->Delete();
  140. GReader->Delete();
  141. implG->Delete();
  142. exit(EXIT_SUCCESS);
  143. }
  144. Real Magnet(const Real& x, const Real& y, const Real& z) {
  145. if (z < 10 || z > 11 ) {
  146. return 0.;
  147. } else if ( std::sqrt(x*x + y*y) <= 0.05 ) {
  148. if (x>0.) return 1.;
  149. if (x<0.) return -1.;
  150. }
  151. return 0.;
  152. }