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

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493
  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. // debugging stuff
  50. #include <vtkDataSetMapper.h>
  51. #include <vtkSelectionNode.h>
  52. #include <vtkExtractSelection.h>
  53. #include <vtkSelection.h>
  54. #include <vtkVertexGlyphFilter.h>
  55. #define RENDERTEST
  56. //#define CONNECTTEST
  57. #ifdef RENDERTEST
  58. #include "vtkRectilinearGridGeometryFilter.h"
  59. #include "vtkRectilinearGridOutlineFilter.h"
  60. #include "vtkPolyDataMapper.h"
  61. #include "vtkActor.h"
  62. #include "vtkRenderWindowInteractor.h"
  63. #include "vtkRenderer.h"
  64. #include "vtkRenderWindow.h"
  65. #include "vtkCamera.h"
  66. #include "vtkProperty.h"
  67. #include <vtkDataSetMapper.h>
  68. #include <vtkSelectionNode.h>
  69. #include <vtkSelection.h>
  70. #include <vtkExtractSelection.h>
  71. #include <vtkExtractEdges.h>
  72. #include <vtkVertexGlyphFilter.h>
  73. #include <vtkTriangleFilter.h>
  74. #include <vtkImplicitHalo.h>
  75. #endif
  76. #ifdef CONNECTTEST
  77. #include <vtkVersion.h>
  78. #include <vtkSmartPointer.h>
  79. #include <vtkPolyData.h>
  80. #include <vtkCellData.h>
  81. #include <vtkDoubleArray.h>
  82. #include <vtkDataSet.h>
  83. #include <vtkSphereSource.h>
  84. #include <vtkTriangleFilter.h>
  85. #include <vtkExtractEdges.h>
  86. #include <vtkDataSetMapper.h>
  87. #include <vtkActor.h>
  88. #include <vtkRenderWindow.h>
  89. #include <vtkRenderer.h>
  90. #include <vtkRenderWindowInteractor.h>
  91. #include <vtkSelectionNode.h>
  92. #include <vtkSelection.h>
  93. #include <vtkExtractSelection.h>
  94. #include <vtkProperty.h>
  95. #include <vtkVertexGlyphFilter.h>
  96. #endif
  97. vtkSmartPointer<vtkIdList> GetConnectedVertices(vtkSmartPointer<vtkDataSet> mesh, int id);
  98. vtkSmartPointer<vtkIdList> GetConnectedVertices2(vtkSmartPointer<vtkPolyData> mesh, int id);
  99. // This function is copyright (C) 2012 Joeseph Cappriotti
  100. vtkSmartPointer<vtkIdList> GetConnectedPoints(int id0, vtkSmartPointer<vtkDataSet> grid);
  101. using namespace Lemma;
  102. int main() {
  103. // Define Sigma term
  104. // Defaults to 1, (Poisson equation)
  105. // Define source (G) term
  106. //vtkSphere *Sphere = vtkSphere::New();
  107. vtkImplicitHalo *Sphere = vtkImplicitHalo::New();
  108. Sphere->SetCenter (0, 0, 0);
  109. Sphere->SetRadius (60);
  110. Sphere->SetFadeOut(.95);
  111. ////////////////////////////////////////////
  112. // Define solution mesh
  113. // Create a rectilinear grid by defining three arrays specifying the
  114. // coordinates in the x-y-z directions.
  115. int nx = 160;
  116. double dx = 2.5;
  117. double ox = -200;
  118. vtkFloatArray *xCoords = vtkFloatArray::New();
  119. for (int i=0; i<nx+1; i++) xCoords->InsertNextValue(ox + i*dx);
  120. int ny = 160;
  121. double dy = 2.5;
  122. double oy = -200;
  123. vtkFloatArray *yCoords = vtkFloatArray::New();
  124. for (int i=0; i<ny+1; i++) yCoords->InsertNextValue(oy + i*dy);
  125. int nz = 160;
  126. double dz = 2.5;
  127. double oz = -200;
  128. vtkFloatArray *zCoords = vtkFloatArray::New();
  129. for (int i=0; i<nz+1; i++) zCoords->InsertNextValue(oz + i*dz);
  130. vtkRectilinearGrid *rGrid = vtkRectilinearGrid::New();
  131. rGrid->SetDimensions(nx+1, ny+1, nz+1);
  132. rGrid->SetExtent(0, nx, 0, ny, 0, nz);
  133. rGrid->SetXCoordinates(xCoords);
  134. rGrid->SetYCoordinates(yCoords);
  135. rGrid->SetZCoordinates(zCoords);
  136. vtkRectilinearGrid *sigGrid = vtkRectilinearGrid::New();
  137. sigGrid->SetDimensions(nx+1, ny+1, nz+1);
  138. sigGrid->SetExtent(0, nx, 0, ny, 0, nz);
  139. sigGrid->SetXCoordinates(xCoords);
  140. sigGrid->SetYCoordinates(yCoords);
  141. sigGrid->SetZCoordinates(zCoords);
  142. vtkDoubleArray *sigArray = vtkDoubleArray::New();
  143. sigArray->SetNumberOfComponents(1);
  144. for (int ii=0; ii<(nx+1)*(ny+1)*(nz+1); ++ii) {
  145. sigArray->InsertTuple1(ii, 1.);
  146. }
  147. sigArray->SetName("sigma");
  148. sigGrid->GetPointData()->AddArray(sigArray);
  149. sigGrid->GetPointData()->SetScalars(sigArray);
  150. vtkImplicitDataSet* implSigma = vtkImplicitDataSet::New();
  151. implSigma->SetDataSet(sigGrid);
  152. implSigma->SetOutValue(0.);
  153. implSigma->SetOutGradient(0., 0., 0.);
  154. FEM4EllipticPDE *Solver = FEM4EllipticPDE::New();
  155. Solver->SetGFunction(Sphere);
  156. //Solver->SetSigmaFunction(implSigma);
  157. Solver->SetBoundaryStep(.5);
  158. Solver->SetGrid(rGrid);
  159. Solver->Solve("FEM_results.vtr");
  160. vtkXMLDataSetWriter *Writer = vtkXMLDataSetWriter::New();
  161. Writer->SetInputData(sigGrid);
  162. std::string fname("sigma.vtr");
  163. //Writer->Update();
  164. //fname.append( Writer->GetDefaultFileExtension() );
  165. Writer->SetFileName(fname.c_str());
  166. Writer->Write();
  167. Writer->Delete();
  168. // Clean up
  169. Sphere->Delete();
  170. Solver->Delete();
  171. rGrid->Delete();
  172. xCoords->Delete();
  173. yCoords->Delete();
  174. zCoords->Delete();
  175. }
  176. // This function is copyright (C) 2012 Joeseph Cappriotti
  177. vtkSmartPointer<vtkIdList> GetConnectedPoints(int id0, vtkSmartPointer<vtkDataSet> grid){
  178. vtkSmartPointer<vtkIdList> pointIds = vtkSmartPointer<vtkIdList>::New();
  179. vtkSmartPointer<vtkIdList> cellList = vtkSmartPointer<vtkIdList>::New();
  180. grid->GetPointCells(id0, cellList);
  181. for(int i=0;i<cellList->GetNumberOfIds();++i){
  182. vtkCell* cell = grid->GetCell(cellList->GetId(i));
  183. if(cell->GetNumberOfEdges() > 0){
  184. for(int j=0; j<cell->GetNumberOfEdges(); ++j){
  185. vtkCell* edge = cell->GetEdge(j);
  186. vtkIdList* edgePoints=edge->GetPointIds();
  187. if(edgePoints->GetId(0)==id0){
  188. pointIds->InsertUniqueId(edgePoints->GetId(1));
  189. }else if(edgePoints->GetId(1)==id0){
  190. pointIds->InsertUniqueId(edgePoints->GetId(0));
  191. }
  192. }
  193. }
  194. }
  195. return pointIds;
  196. }
  197. vtkSmartPointer<vtkIdList> GetConnectedVertices(vtkSmartPointer<vtkDataSet> mesh, int id) {
  198. std::cout << "number of points " << mesh->GetNumberOfPoints() << std::endl;
  199. std::cout << "number of cells " << mesh->GetNumberOfCells() << std::endl;
  200. vtkSmartPointer<vtkIdList> connectedVertices = vtkSmartPointer<vtkIdList>::New();
  201. //get all cells that vertex 'id' is a part of
  202. vtkSmartPointer<vtkIdList> cellIdList = vtkSmartPointer<vtkIdList>::New();
  203. mesh->GetPointCells(id, cellIdList); // cells using point
  204. //mesh->GetCellPoints(id, cellIdList); // points defining cell
  205. cout << "Vertex " << id << " is used in cell(s) ";
  206. for(vtkIdType i = 0; i < cellIdList->GetNumberOfIds(); i++) {
  207. cout << cellIdList->GetId(i) << ", ";
  208. }
  209. cout << endl;
  210. // TODO this is where things don't work
  211. // for each cell, figure out what vertices are connected
  212. for(vtkIdType i = 0; i < cellIdList->GetNumberOfIds(); i++) {
  213. cout << "\tcell id " << i << " : " << cellIdList->GetId(i) << endl;
  214. std::cout << "\tcell has " << mesh->GetCell(i)->GetNumberOfPoints() << " points\n";
  215. vtkExtractCells *Cell = vtkExtractCells::New();
  216. Cell->SetInputData(mesh);
  217. Cell->AddCellRange(i, i);
  218. Cell->Update();
  219. //std::cout << *Cell->GetOutput() << std::endl;
  220. // extract surface
  221. vtkDataSetSurfaceFilter *surf = vtkDataSetSurfaceFilter::New();
  222. surf->UseStripsOn();
  223. surf->SetInputData(Cell->GetOutput());
  224. surf->Update();
  225. //std::cout << *surf->GetOutput() << std::endl;
  226. vtkSmartPointer<vtkTriangleFilter> triangleFilter =
  227. vtkSmartPointer<vtkTriangleFilter>::New();
  228. triangleFilter->SetInputData( surf->GetOutput() );
  229. triangleFilter->Update();
  230. vtkSmartPointer<vtkExtractEdges> extractEdges =
  231. vtkSmartPointer<vtkExtractEdges>::New();
  232. extractEdges->SetInputConnection(triangleFilter->GetOutputPort());
  233. extractEdges->Update();
  234. vtkSmartPointer<vtkPolyData> pmesh = extractEdges->GetOutput();
  235. //vtkSmartPointer<vtkPolyData> pmesh = triangleFilter->GetOutput();
  236. //vtkSmartPointer<vtkIdList> pconnectedVertices = GetConnectedVertices2(surf->GetOutput(), id);
  237. //
  238. // //vtkSmartPointer<vtkIdTypeArray> ids =
  239. // //vtkSmartPointer<vtkIdTypeArray>::New();
  240. // //ids->SetNumberOfComponents(1);
  241. //
  242. vtkSmartPointer<vtkIdList> pointIdList = vtkSmartPointer<vtkIdList>::New();
  243. //mesh->GetPointCells(cellIdList->GetId(i), pointIdList); // returns cells using point
  244. mesh->GetCellPoints(cellIdList->GetId(i), pointIdList); // points defining cell
  245. //vtkSmartPointer<vtkIdList> pconnectedVertices = GetConnectedVertices2(pmesh, id);
  246. vtkSmartPointer<vtkIdList> pconnectedVertices = GetConnectedVertices2(surf->GetOutput(), id);
  247. //vtkSmartPointer<vtkIdList> pconnectedVertices = GetConnectedVertices2(triangleFilter->GetOutput(), 1);
  248. std::cout << "\tPoly Connected vertices: ";
  249. for(vtkIdType ip = 0; ip < pconnectedVertices->GetNumberOfIds(); ip++) {
  250. std::cout << pconnectedVertices->GetId(ip) << " ";
  251. //ids->InsertNextValue(connectedVertices->GetId(i));
  252. connectedVertices->InsertNextId( pointIdList->GetId( pconnectedVertices->GetId(ip) ) );
  253. }
  254. std::cout << std::endl;
  255. cout << "\tcell " << i << " contains these points ";
  256. for(vtkIdType ii = 0; ii < pointIdList->GetNumberOfIds(); ii++) {
  257. cout << pointIdList->GetId(ii) << ", ";
  258. }
  259. cout << endl;
  260. /*
  261. cout << "\tEnd points are " << pointIdList->GetId(0) << " and " << pointIdList->GetId(1) << endl;
  262. if(pointIdList->GetId(0) != id) {
  263. //cout << "Connected to " << pointIdList->GetId(0) << endl;
  264. connectedVertices->InsertNextId(pointIdList->GetId(0));
  265. } else {
  266. //cout << "Connected to " << pointIdList->GetId(1) << endl;
  267. connectedVertices->InsertNextId(pointIdList->GetId(1));
  268. }
  269. */
  270. #ifdef RENDERTEST
  271. vtkSmartPointer<vtkDataSetMapper> sphereMapper =
  272. vtkSmartPointer<vtkDataSetMapper>::New();
  273. sphereMapper->SetInputConnection(extractEdges->GetOutputPort());
  274. vtkSmartPointer<vtkActor> sphereActor =
  275. vtkSmartPointer<vtkActor>::New();
  276. sphereActor->SetMapper(sphereMapper);
  277. vtkSmartPointer<vtkIdTypeArray> ids =
  278. vtkSmartPointer<vtkIdTypeArray>::New();
  279. ids->SetNumberOfComponents(1);
  280. std::cout << "RENDERTEST Connected vertices: ";
  281. for(vtkIdType i = 0; i < connectedVertices->GetNumberOfIds(); i++) {
  282. std::cout << connectedVertices->GetId(i) << " ";
  283. ids->InsertNextValue(connectedVertices->GetId(i));
  284. }
  285. std::cout << std::endl;
  286. vtkSmartPointer<vtkDataSetMapper> connectedVertexMapper =
  287. vtkSmartPointer<vtkDataSetMapper>::New();
  288. {
  289. vtkSmartPointer<vtkSelectionNode> selectionNode =
  290. vtkSmartPointer<vtkSelectionNode>::New();
  291. selectionNode->SetFieldType(vtkSelectionNode::POINT);
  292. selectionNode->SetContentType(vtkSelectionNode::INDICES);
  293. selectionNode->SetSelectionList(ids);
  294. vtkSmartPointer<vtkSelection> selection =
  295. vtkSmartPointer<vtkSelection>::New();
  296. selection->AddNode(selectionNode);
  297. vtkSmartPointer<vtkExtractSelection> extractSelection =
  298. vtkSmartPointer<vtkExtractSelection>::New();
  299. extractSelection->SetInputConnection(0, extractEdges->GetOutputPort());
  300. extractSelection->SetInputData(1, selection);
  301. extractSelection->Update();
  302. vtkSmartPointer<vtkVertexGlyphFilter> glyphFilter =
  303. vtkSmartPointer<vtkVertexGlyphFilter>::New();
  304. glyphFilter->SetInputConnection(extractSelection->GetOutputPort());
  305. glyphFilter->Update();
  306. connectedVertexMapper->SetInputConnection(glyphFilter->GetOutputPort());
  307. }
  308. vtkSmartPointer<vtkActor> connectedVertexActor =
  309. vtkSmartPointer<vtkActor>::New();
  310. connectedVertexActor->SetMapper(connectedVertexMapper);
  311. connectedVertexActor->GetProperty()->SetColor(1,0,0);
  312. connectedVertexActor->GetProperty()->SetPointSize(5);
  313. vtkSmartPointer<vtkDataSetMapper> queryVertexMapper =
  314. vtkSmartPointer<vtkDataSetMapper>::New();
  315. {
  316. vtkSmartPointer<vtkIdTypeArray> ids =
  317. vtkSmartPointer<vtkIdTypeArray>::New();
  318. ids->SetNumberOfComponents(1);
  319. ids->InsertNextValue(id);
  320. vtkSmartPointer<vtkSelectionNode> selectionNode =
  321. vtkSmartPointer<vtkSelectionNode>::New();
  322. selectionNode->SetFieldType(vtkSelectionNode::POINT);
  323. selectionNode->SetContentType(vtkSelectionNode::INDICES);
  324. selectionNode->SetSelectionList(ids);
  325. vtkSmartPointer<vtkSelection> selection =
  326. vtkSmartPointer<vtkSelection>::New();
  327. selection->AddNode(selectionNode);
  328. vtkSmartPointer<vtkExtractSelection> extractSelection =
  329. vtkSmartPointer<vtkExtractSelection>::New();
  330. extractSelection->SetInputConnection(0, extractEdges->GetOutputPort());
  331. extractSelection->SetInputData(1, selection);
  332. extractSelection->Update();
  333. vtkSmartPointer<vtkVertexGlyphFilter> glyphFilter =
  334. vtkSmartPointer<vtkVertexGlyphFilter>::New();
  335. glyphFilter->SetInputConnection(extractSelection->GetOutputPort());
  336. glyphFilter->Update();
  337. queryVertexMapper->SetInputConnection(glyphFilter->GetOutputPort());
  338. }
  339. vtkSmartPointer<vtkActor> queryVertexActor =
  340. vtkSmartPointer<vtkActor>::New();
  341. queryVertexActor->SetMapper(queryVertexMapper);
  342. queryVertexActor->GetProperty()->SetColor(0,1,0);
  343. queryVertexActor->GetProperty()->SetPointSize(5);
  344. // Create the usual rendering stuff.
  345. vtkRenderer *renderer = vtkRenderer::New();
  346. vtkRenderWindow *renWin = vtkRenderWindow::New();
  347. renWin->AddRenderer(renderer);
  348. vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
  349. iren->SetRenderWindow(renWin);
  350. renderer->AddActor(queryVertexActor);
  351. renderer->AddActor(connectedVertexActor);
  352. renderer->AddActor(sphereActor);
  353. //renderer->AddActor(wireActor);
  354. //renderer->AddActor(owireActor);
  355. renderer->SetBackground(.3,.2,.1);
  356. renderer->ResetCamera();
  357. renderer->GetActiveCamera()->Elevation(-150.0);
  358. renderer->GetActiveCamera()->Azimuth(0.0);
  359. renderer->GetActiveCamera()->Roll(0.0);
  360. //renderer->GetActiveCamera()->Pitch(-45.0);
  361. renderer->GetActiveCamera()->Zoom(1.0);
  362. renWin->SetSize(300,300);
  363. // interact with data
  364. renWin->Render();
  365. iren->Start();
  366. // clean up rendering stuff
  367. //plane->Delete();
  368. //rgridMapper->Delete();
  369. //wireActor->Delete();
  370. //owireActor->Delete();
  371. //outline->Delete();
  372. //outlineMapper->Delete();
  373. renderer->Delete();
  374. renWin->Delete();
  375. iren->Delete();
  376. #endif // RENDERTEST
  377. surf->Delete();
  378. Cell->Delete();
  379. }
  380. return connectedVertices;
  381. }
  382. // from http://www.vtk.org/Wiki/VTK/Examples/Cxx/PolyData/VertexConnectivity
  383. vtkSmartPointer<vtkIdList> GetConnectedVertices2(vtkSmartPointer<vtkPolyData> mesh, int id) {
  384. std::cout << "mesh points " << mesh->GetNumberOfPoints() << std::endl;
  385. vtkSmartPointer<vtkIdList> connectedVertices =
  386. vtkSmartPointer<vtkIdList>::New();
  387. //get all cells that vertex 'id' is a part of
  388. vtkSmartPointer<vtkIdList> cellIdList =
  389. vtkSmartPointer<vtkIdList>::New();
  390. mesh->GetPointCells(id, cellIdList);
  391. /*
  392. for(vtkIdType i = 0; i < cellIdList->GetNumberOfIds(); i++) {
  393. cout << cellIdList->GetId(i) << ", ";
  394. }
  395. cout << endl;
  396. */
  397. for(vtkIdType i = 0; i < cellIdList->GetNumberOfIds(); i++) {
  398. //cout << "id " << i << " : " << cellIdList->GetId(i) << endl;
  399. vtkSmartPointer<vtkIdList> pointIdList =
  400. vtkSmartPointer<vtkIdList>::New();
  401. mesh->GetCellPoints(cellIdList->GetId(i), pointIdList);
  402. //cout << "End points are " << pointIdList->GetId(0) << " and " << pointIdList->GetId(1) << endl;
  403. if(pointIdList->GetId(0) != id) {
  404. //cout << "Connected to " << pointIdList->GetId(0) << endl;
  405. connectedVertices->InsertNextId(pointIdList->GetId(0));
  406. } else {
  407. //cout << "Connected to " << pointIdList->GetId(1) << endl;
  408. connectedVertices->InsertNextId(pointIdList->GetId(1));
  409. }
  410. }
  411. return connectedVertices;
  412. }