Galerkin FEM for elliptic PDEs
Nevar pievienot vairāk kā 25 tēmas Tēmai ir jāsākas ar burtu vai ciparu, tā var saturēt domu zīmes ('-') un var būt līdz 35 simboliem gara.

VTKEdgeGsphere.cpp 7.2KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200
  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 <vtkCellData.h>
  26. #include <vtkCell.h>
  27. #include <cmath>
  28. #include <Eigen/Core>
  29. const double PI = 4.0*atan(1.0);
  30. int main(int argc, char**argv) {
  31. std::cout << "Mesh processing routine\n";
  32. if (argc<4) {
  33. std::cout << "usage:\n" << "./utVTKEdge inMesh.vtu <radius> <problemType> outG.vtu" << std::endl;
  34. exit(EXIT_SUCCESS);
  35. }
  36. vtkUnstructuredGrid* uGrid;// = vtkUnstructuredGrid::New();
  37. std::string fn = argv[1];
  38. if(fn.substr(fn.find_last_of(".") + 1) == "vtk") {
  39. vtkGenericDataObjectReader* Reader = vtkGenericDataObjectReader::New();
  40. Reader->SetFileName(argv[1]);
  41. Reader->Update();
  42. if(Reader->IsFileUnstructuredGrid()) {
  43. std::cout << "Found Unstructured grid legacy file" << std::endl;
  44. uGrid = Reader->GetUnstructuredGridOutput();
  45. } else {
  46. std::cerr << "Unknown legacy format";
  47. exit(EXIT_FAILURE);
  48. }
  49. } else {
  50. vtkXMLUnstructuredGridReader* Reader = vtkXMLUnstructuredGridReader::New();
  51. std::cout << "Reading" << argv[1] << std::endl;
  52. Reader->SetFileName(argv[1]);
  53. Reader->Update();
  54. uGrid = Reader->GetOutput();
  55. }
  56. int nc = uGrid->GetNumberOfCells() ;
  57. int nn = uGrid->GetNumberOfPoints() ;
  58. std::cout << "Processing grid with nodes=" << nn << "\t cells=" << nc << "\n";
  59. // Input magnet size TODO get from file or .geo file even?
  60. double R = atof(argv[2]); //.25; // radius of magnet
  61. double eps = 1e-4; // epsilon for on the point
  62. // Pol = double[3];
  63. // Declare new array for the data
  64. vtkDoubleArray* G = vtkDoubleArray::New();
  65. G->SetNumberOfComponents(1);
  66. G->SetNumberOfTuples(nn);
  67. //G->SetName("G");
  68. if ( std::string(argv[3]) == "mag") {
  69. G->SetName("kappa");
  70. G->SetNumberOfTuples(nc);
  71. }
  72. else {
  73. G->SetName("G");
  74. G->SetNumberOfTuples(nn);
  75. }
  76. vtkDoubleArray* phi = vtkDoubleArray::New();
  77. phi->SetNumberOfComponents(1);
  78. phi->SetNumberOfTuples(nn);
  79. phi->SetName("analytic_phi");
  80. double point[3];
  81. if ( std::string(argv[3]) == "mag") {
  82. // Loop over cells and decide if all of them lie within sphere
  83. for (int ic=0; ic<nc; ++ic) {
  84. assert ( uGrid->GetCell(ic)->GetNumberOfPoints() == 4 );
  85. // TODO, in production code we might not want to do this check here
  86. if ( uGrid->GetCell(ic)->GetNumberOfPoints() != 4 ) {
  87. throw std::runtime_error("Non-tetrahedral mesh encountered!");
  88. }
  89. // construct coordinate matrix C
  90. //Eigen::Matrix<Real, 4, 4> C = Eigen::Matrix<Real, 4, 4>::Zero() ;
  91. bool cellin(true);
  92. for (int ii=0; ii<4; ++ii) {
  93. double* point = uGrid->GetCell(ic)->GetPoints()->GetPoint(ii); //[ipc] ;
  94. if ( sqrt( point[0]*point[0] + point[1]*point[1] + point[2]*point[2] ) > R+eps ) {
  95. cellin = false;
  96. }
  97. }
  98. if (cellin) {
  99. G->InsertTuple1( ic, 1 );
  100. } else {
  101. G->InsertTuple1( ic, 0 );
  102. }
  103. }
  104. uGrid->GetCellData()->AddArray( G );
  105. } else {
  106. // Loop over nodes, look at position, set with G if on boundary
  107. Eigen::Vector3d M; M << 0,0,1;
  108. for (int in=0; in<nn; ++in) {
  109. uGrid->GetPoint(in, &point[0]);
  110. //std::cout << point[0] << "\t" <<point[1] << "\t" << point[2] << std::endl;
  111. double raddist = sqrt( point[0]*point[0] + point[1]*point[1] + point[2]*point[2] );
  112. //double rho = sqrt( point[1]*point[1] + point[2]*point[2] );
  113. // Calculate RHS
  114. if ( std::string(argv[3]) == "gravity") {
  115. if ( raddist < R + eps) {
  116. G->InsertTuple1( in, 1 );
  117. } else {
  118. G->InsertTuple1( in, 0 );
  119. }
  120. }
  121. if ( std::string(argv[3]) == "magnet") {
  122. if ( raddist > R - eps && raddist < R + eps ) {
  123. // \rho = \nabla \cdot \mathbf{M}
  124. // Use divergence theorm --> \mahtbf{M} \cdot \hat{n}
  125. Eigen::Vector3d n;
  126. n << point[0],point[1],point[2];
  127. n.array() /= raddist;
  128. G->InsertTuple1(in, n.dot(M) );
  129. } else {
  130. G->InsertTuple1( in, 0 );
  131. }
  132. }
  133. // Insert analytic phi part
  134. /* magnetics problem p. 198 Jackson */
  135. if (std::string(argv[3]) == "magnet") {
  136. if (raddist < R) {
  137. phi->InsertTuple1( in, (1./3.)*point[2] );
  138. } else {
  139. phi->InsertTuple1( in, 0);
  140. double costheta = point[2]/raddist ;
  141. //phi->InsertTuple1( in, -(1./3.)*(R*R*R) * ( costheta / (rho*rho) ) );
  142. phi->InsertTuple1( in, (1./3.) * (R*R*R) * (costheta / (raddist*raddist)) );
  143. }
  144. } else if (std::string(argv[3]) == "gravity") {
  145. double mass = 4./3. * PI * (R*R*R); // volume * density (1)
  146. if (raddist < R) {
  147. phi->InsertTuple1( in, mass * (( 3*(R*R) - raddist*raddist )/(2*(R*R*R))) ); // (1./3.)*point[2] );
  148. //phi->InsertTuple1( in, (2*PI/3)*(3*R*R-raddist*raddist) ); // (1./3.)*point[2] );
  149. } else {
  150. //phi->InsertTuple1( in, 0);
  151. //double costheta = point[2]/raddist ;
  152. //phi->InsertTuple1( in, -(1./3.)*(R*R*R) * ( costheta / (rho*rho) ) );
  153. phi->InsertTuple1( in, mass/raddist );
  154. }
  155. }
  156. }
  157. // Add new data
  158. uGrid->GetPointData()->AddArray( phi );
  159. uGrid->GetPointData()->SetScalars( G );
  160. }
  161. std::cout << "Writing file with ncells=" << uGrid->GetNumberOfCells() << "\tnpoints=" << uGrid->GetNumberOfPoints() << std::endl;
  162. // Write out new file
  163. vtkXMLUnstructuredGridWriter* Writer = vtkXMLUnstructuredGridWriter::New();
  164. //Writer->SetInputConnection(Reader->GetOutputPort());
  165. Writer->SetInputData( uGrid );
  166. Writer->SetFileName( argv[4] );
  167. Writer->Update();
  168. Writer->Write();
  169. //Reader->Delete();
  170. }