Galerkin FEM for elliptic PDEs
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495
  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 <vtkXMLUnstructuredGridReader.h>
  21. #include <vtkXMLUnstructuredGridWriter.h>
  22. #include <vtkDoubleArray.h>
  23. #include <vtkPointData.h>
  24. #include <cmath>
  25. int main(int argc, char**argv) {
  26. std::cout << "Mesh processing routine\n";
  27. if (argc<3) {
  28. std::cout << "usage:\n" << "./utVTKEdge inMesh.vtu outG.vtu" << std::endl;
  29. }
  30. vtkXMLUnstructuredGridReader* Reader = vtkXMLUnstructuredGridReader::New();
  31. std::cout << "Reading" << argv[1] << std::endl;
  32. Reader->SetFileName(argv[1]);
  33. Reader->Update();
  34. int nc = Reader->GetOutput()->GetNumberOfCells() ;
  35. int nn = Reader->GetOutput()->GetNumberOfPoints() ;
  36. std::cout << "Processing grid with nodes=" << nn << "\t cells=" << nc << "\n";
  37. // Input magnet size TODO get from file or .geo file even?
  38. double R = .25; // radius of magnet
  39. double D0 = 10; // Top of magnet
  40. double D1 = 11; // Bottom of magnet
  41. double eps = 1e-2; // epsilon for on the point
  42. // Pol = double[3];
  43. // Declare new array for the data
  44. vtkDoubleArray* G = vtkDoubleArray::New();
  45. G->SetNumberOfComponents(1);
  46. G->SetNumberOfTuples(nn);
  47. G->SetName("G");
  48. // Loop over nodes, look at position, set with G if on boundary
  49. double point[3];
  50. for (int in=0; in<nn; ++in) {
  51. Reader->GetOutput()->GetPoint(in, &point[0]);
  52. //std::cout << point[0] << "\t" <<point[1] << "\t" << point[2] << std::endl;
  53. double raddist = sqrt( point[0]*point[0] + point[1]*point[1] );
  54. if ( raddist > R - eps && raddist < R + eps && point[2] < D1 && point[2] > D0 ) {
  55. if ( point[1] < 0 ) {
  56. G->InsertTuple1( in, 1 );
  57. } else if (point[1] > 0) {
  58. G->InsertTuple1( in, -1 );
  59. } else {
  60. G->InsertTuple1( in, 0 );
  61. }
  62. } else {
  63. G->InsertTuple1( in, 0 );
  64. }
  65. }
  66. // Add new data
  67. Reader->GetOutput()->GetPointData()->SetScalars( G );
  68. // Write out new file
  69. vtkXMLUnstructuredGridWriter* Writer = vtkXMLUnstructuredGridWriter::New();
  70. Writer->SetInputConnection(Reader->GetOutputPort());
  71. Writer->SetFileName( argv[2] );
  72. Writer->Write();
  73. Reader->Delete();
  74. }