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.

VTKEdgeGsphere.cpp 5.9KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173
  1. /* This file is part of Lemma, a geophysical modelling and inversion API.
  2. * More information is available at http://lemmasoftware.org
  3. */
  4. /* This Source Code Form is subject to the terms of the Mozilla Public
  5. * License, v. 2.0. If a copy of the MPL was not distributed with this
  6. * file, You can obtain one at http://mozilla.org/MPL/2.0/.
  7. */
  8. /**
  9. * @file
  10. * @date 08/07/2014 04:16:25 PM
  11. * @version $Id$
  12. * @author Trevor Irons (ti)
  13. * @email Trevor.Irons@xri-geo.com
  14. * @copyright Copyright (c) 2014, XRI Geophysics, LLC
  15. * @copyright Copyright (c) 2014, Trevor Irons
  16. */
  17. #include <vtkUnstructuredGrid.h>
  18. #include <vtkUnstructuredGridReader.h>
  19. #include <vtkUnstructuredGridWriter.h>
  20. #include <vtkGenericDataObjectReader.h>
  21. #include <vtkXMLUnstructuredGridReader.h>
  22. #include <vtkXMLUnstructuredGridWriter.h>
  23. #include <vtkDoubleArray.h>
  24. #include <vtkPointData.h>
  25. #include <cmath>
  26. #include <Eigen/Core>
  27. const double PI = 4.0*atan(1.0);
  28. int main(int argc, char**argv) {
  29. std::cout << "Mesh processing routine\n";
  30. if (argc<4) {
  31. std::cout << "usage:\n" << "./utVTKEdge inMesh.vtu <radius> <problemType> outG.vtu" << std::endl;
  32. exit(EXIT_SUCCESS);
  33. }
  34. vtkUnstructuredGrid* uGrid = vtkUnstructuredGrid::New();
  35. std::string fn = argv[1];
  36. if(fn.substr(fn.find_last_of(".") + 1) == "vtk") {
  37. vtkGenericDataObjectReader* Reader = vtkGenericDataObjectReader::New();
  38. Reader->SetFileName(argv[1]);
  39. Reader->Update();
  40. if(Reader->IsFileUnstructuredGrid()) {
  41. std::cout << "Found Unstructured grid legacy file" << std::endl;
  42. uGrid = Reader->GetUnstructuredGridOutput();
  43. } else {
  44. std::cerr << "Unknown legacy format";
  45. exit(EXIT_FAILURE);
  46. }
  47. } else {
  48. vtkXMLUnstructuredGridReader* Reader = vtkXMLUnstructuredGridReader::New();
  49. std::cout << "Reading" << argv[1] << std::endl;
  50. Reader->SetFileName(argv[1]);
  51. Reader->Update();
  52. uGrid = Reader->GetOutput();
  53. }
  54. int nc = uGrid->GetNumberOfCells() ;
  55. int nn = uGrid->GetNumberOfPoints() ;
  56. std::cout << "Processing grid with nodes=" << nn << "\t cells=" << nc << "\n";
  57. // Input magnet size TODO get from file or .geo file even?
  58. double R = atof(argv[2]); //.25; // radius of magnet
  59. double eps = 1e-4; // epsilon for on the point
  60. // Pol = double[3];
  61. // Declare new array for the data
  62. vtkDoubleArray* G = vtkDoubleArray::New();
  63. G->SetNumberOfComponents(1);
  64. G->SetNumberOfTuples(nn);
  65. G->SetName("G");
  66. vtkDoubleArray* phi = vtkDoubleArray::New();
  67. phi->SetNumberOfComponents(1);
  68. phi->SetNumberOfTuples(nn);
  69. phi->SetName("analytic_phi");
  70. // Loop over nodes, look at position, set with G if on boundary
  71. double point[3];
  72. Eigen::Vector3d M; M << 0,0,1;
  73. for (int in=0; in<nn; ++in) {
  74. uGrid->GetPoint(in, &point[0]);
  75. //std::cout << point[0] << "\t" <<point[1] << "\t" << point[2] << std::endl;
  76. double raddist = sqrt( point[0]*point[0] + point[1]*point[1] + point[2]*point[2] );
  77. //double rho = sqrt( point[1]*point[1] + point[2]*point[2] );
  78. // Calculate RHS
  79. if ( std::string(argv[3]) == "gravity") {
  80. if ( raddist < R + eps) {
  81. G->InsertTuple1( in, 1 );
  82. } else {
  83. G->InsertTuple1( in, 0 );
  84. }
  85. }
  86. else if ( std::string(argv[3]) == "magnet") {
  87. if ( raddist > R - eps && raddist < R + eps ) {
  88. // \rho = \nabla \cdot \mathbf{M}
  89. //
  90. if ( point[2] < -eps ) {
  91. G->InsertTuple1( in, 1 );
  92. } else if (point[2] > eps) {
  93. G->InsertTuple1( in, -1 );
  94. } else {
  95. G->InsertTuple1( in, 0 );
  96. }
  97. //
  98. // Above is inaccurate and unstable
  99. // Use divergence theorm --> \mahtbf{M} \cdot \hat{n}
  100. Eigen::Vector3d n;
  101. n << point[0],point[1],point[2];
  102. n.array() /= raddist;
  103. //double sigma = n.dot(M);
  104. G->InsertTuple1(in, n.dot(M) );
  105. } else {
  106. G->InsertTuple1( in, 0 );
  107. }
  108. }
  109. // Insert analytic phi part
  110. /* magnetics problem p. 198 Jackson */
  111. if (std::string(argv[3]) == "magnet") {
  112. if (raddist < R) {
  113. phi->InsertTuple1( in, (1./3.)*point[2] );
  114. } else {
  115. phi->InsertTuple1( in, 0);
  116. double costheta = point[2]/raddist ;
  117. //phi->InsertTuple1( in, -(1./3.)*(R*R*R) * ( costheta / (rho*rho) ) );
  118. phi->InsertTuple1( in, (1./3.) * (R*R*R) * (costheta / (raddist*raddist)) );
  119. }
  120. } else if (std::string(argv[3]) == "gravity") {
  121. double mass = 4./3. * PI * (R*R*R); // volume * density (1)
  122. if (raddist < R) {
  123. phi->InsertTuple1( in, mass * (( 3*(R*R) - raddist*raddist )/(2*(R*R*R))) ); // (1./3.)*point[2] );
  124. //phi->InsertTuple1( in, (2*PI/3)*(3*R*R-raddist*raddist) ); // (1./3.)*point[2] );
  125. } else {
  126. //phi->InsertTuple1( in, 0);
  127. //double costheta = point[2]/raddist ;
  128. //phi->InsertTuple1( in, -(1./3.)*(R*R*R) * ( costheta / (rho*rho) ) );
  129. phi->InsertTuple1( in, mass/raddist );
  130. }
  131. }
  132. }
  133. // Add new data
  134. uGrid->GetPointData()->AddArray( phi );
  135. uGrid->GetPointData()->SetScalars( G );
  136. // Write out new file
  137. vtkXMLUnstructuredGridWriter* Writer = vtkXMLUnstructuredGridWriter::New();
  138. //Writer->SetInputConnection(Reader->GetOutputPort());
  139. Writer->SetInputData( uGrid );
  140. Writer->SetFileName( argv[4] );
  141. Writer->Write();
  142. //Reader->Delete();
  143. }