Main Lemma Repository
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.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345
  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 01/29/2014 09:53:49 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 <iostream>
  18. #include <fstream>
  19. #include "dipolesource.h"
  20. #include "layeredearthem.h"
  21. #include "receiverpoints.h"
  22. #include "emearth1d.h"
  23. #ifdef LEMMAUSEOMP
  24. #include "omp.h"
  25. #endif
  26. // For testing purposes disable VTK and run scale.sh
  27. //#undef LEMMAUSEVTK
  28. //#undef KIHALEE_EM1D
  29. #ifdef LEMMAUSEVTK
  30. #include "vtkRenderer.h"
  31. #include "vtkRenderWindow.h"
  32. #include "vtkRenderWindowInteractor.h"
  33. #include "vtkDoubleArray.h"
  34. #include "vtkXYPlotActor.h"
  35. #include "vtkXYPlotWidget.h"
  36. #include "vtkProperty2D.h"
  37. #endif
  38. #include "timer.h"
  39. using namespace Lemma;
  40. int main() {
  41. // Do calculation
  42. jsw_timer timer;
  43. DipoleSource *dipole = DipoleSource::New();
  44. dipole->SetType(MAGNETICDIPOLE);
  45. dipole->SetPolarisation(0., 0.0, 1.0);
  46. //dipole->SetPolarisation(XPOLARISATION);
  47. //dipole->SetPolarisation(YPOLARISATION);
  48. //dipole->SetPolarisation(ZPOLARISATION);
  49. //dipole->SetMoment(1);
  50. //dipole->SetLocation(1,1,-.0100328);
  51. dipole->SetLocation(0., 0., -0.001);
  52. //dipole->SetLocation(-2.5,1.25,0);
  53. int iif=10000;
  54. dipole->SetNumberOfFrequencies(iif);
  55. Real w0 = 10;
  56. for (int iw=0; iw<iif; ++iw) {
  57. dipole->SetFrequency(iw, w0);
  58. w0 += 1523.;
  59. }
  60. VectorXcr sigma(8);
  61. sigma << 0., 3e-4, 5e-4, 1e-2, .5, 5e-6, .03, .04;
  62. VectorXr thick(6);
  63. thick << 10 , 10 , 10 , 10 , 10 , 10;
  64. LayeredEarthEM *earth = LayeredEarthEM::New();
  65. earth->SetNumberOfLayers(8);
  66. earth->SetLayerConductivity(sigma);
  67. earth->SetLayerThickness(thick);
  68. ReceiverPoints *receivers = ReceiverPoints::New();
  69. Vector3r loc;
  70. receivers->SetNumberOfReceivers(1);
  71. loc << 150, 175, -.01;
  72. receivers->SetLocation(0, loc);
  73. // EmEarth
  74. EMEarth1D *EmEarth = EMEarth1D::New();
  75. EmEarth->AttachDipoleSource(dipole);
  76. EmEarth->AttachLayeredEarthEM(earth);
  77. EmEarth->AttachReceiverPoints(receivers);
  78. //EmEarth->SetFieldsToCalculate(BOTH);
  79. EmEarth->SetFieldsToCalculate(H);
  80. //EmEarth->SetHankelTransformMethod(GAUSSIANQUADRATURE);
  81. //EmEarth->SetHankelTransformMethod(DIGITALFILTERING);
  82. //EmEarth->SetHankelTransformMethod(FHTKEY201);
  83. EmEarth->SetHankelTransformMethod(QWEKEY);
  84. timer.begin();
  85. EmEarth->MakeCalc3();
  86. //std::cout << "H\n" << receivers->GetHfield(0)
  87. // << "End of H" << std::endl;
  88. Real lemmaTime = timer.end();
  89. //std::cout << receivers->GetHfield() << std::endl;
  90. std::cout << lemmaTime << "\t";
  91. #ifdef LEMMAUSEVTK
  92. vtkDataObject *_dataObjectxr = receivers->GetVtkDataObjectFreq(HFIELDREAL, 0,
  93. 0, iif, XCOMPONENT, dipole->GetFrequencies());
  94. vtkDataObject *_dataObjectxi = receivers->GetVtkDataObjectFreq(HFIELDIMAG, 0,
  95. 0, iif, XCOMPONENT, dipole->GetFrequencies());
  96. vtkDataObject *_dataObjectyr= receivers->GetVtkDataObjectFreq(HFIELDREAL, 0,
  97. 0, iif, YCOMPONENT, dipole->GetFrequencies());
  98. vtkDataObject *_dataObjectyi = receivers->GetVtkDataObjectFreq(HFIELDIMAG, 0,
  99. 0, iif, YCOMPONENT, dipole->GetFrequencies());
  100. vtkDataObject *_dataObjectzr = receivers->GetVtkDataObjectFreq(HFIELDREAL, 0,
  101. 0, iif, ZCOMPONENT, dipole->GetFrequencies());
  102. vtkDataObject *_dataObjectzi = receivers->GetVtkDataObjectFreq(HFIELDIMAG, 0,
  103. 0, iif, ZCOMPONENT, dipole->GetFrequencies());
  104. #endif
  105. // //dipole->SetPolarisation(XPOLARISATION);
  106. // //std::cout << "Fortran Calc" << std::endl;
  107. // #ifdef KIHALEE_EM1D
  108. // timer.begin();
  109. // EmEarth->MakeCalc();
  110. // Real fortranTime = timer.end();
  111. // std::cout << fortranTime << "\t";
  112. // std::cout << fortranTime/lemmaTime << std::endl;
  113. // #endif
  114. // #ifdef LEMMAUSEVTK
  115. // vtkDataObject *_fdataObjectxr = receivers->GetVtkDataObjectFreq(HFIELDREAL, 0,
  116. // 0, receivers->GetNumberOfReceivers(), XCOMPONENT, dipole->GetFrequencies());
  117. //
  118. // vtkDataObject *_fdataObjectxi = receivers->GetVtkDataObjectFreq(HFIELDIMAG, 0,
  119. // 0, receivers->GetNumberOfReceivers(), XCOMPONENT, dipole->GetFrequencies());
  120. //
  121. // vtkDataObject *_fdataObjectyr = receivers->GetVtkDataObjectFreq(HFIELDREAL, 0,
  122. // 0, receivers->GetNumberOfReceivers(), YCOMPONENT, dipole->GetFrequencies());
  123. //
  124. // vtkDataObject *_fdataObjectyi = receivers->GetVtkDataObjectFreq(HFIELDIMAG, 0,
  125. // 0, receivers->GetNumberOfReceivers(), YCOMPONENT, dipole->GetFrequencies());
  126. //
  127. // vtkDataObject *_fdataObjectzr = receivers->GetVtkDataObjectFreq(HFIELDREAL, 0,
  128. // 0, receivers->GetNumberOfReceivers(), ZCOMPONENT, dipole->GetFrequencies());
  129. //
  130. // vtkDataObject *_fdataObjectzi = receivers->GetVtkDataObjectFreq(HFIELDIMAG, 0,
  131. // 0, receivers->GetNumberOfReceivers(), ZCOMPONENT, dipole->GetFrequencies());
  132. // #endif
  133. #ifdef LEMMAUSEVTK
  134. vtkRenderer *_ren = vtkRenderer::New();
  135. vtkRenderWindow *_renwin = vtkRenderWindow::New();
  136. _renwin->AddRenderer(_ren);
  137. ///////////////////////////////////////////////
  138. // X Component Plot
  139. vtkXYPlotActor *_xplot = vtkXYPlotActor::New();
  140. //vtkXYPlotWidget *_wplot = vtkXYPlotWidget::New();
  141. // Add the datasets
  142. _xplot->AddDataObjectInput( _dataObjectxr);
  143. _xplot->AddDataObjectInput( _dataObjectxi);
  144. //_xplot->AddDataObjectInput(_fdataObjectxr);
  145. //_xplot->AddDataObjectInput(_fdataObjectxi);
  146. _xplot->SetTitle("X field");
  147. _xplot->SetXTitle("Depth");
  148. _xplot->SetYTitle("E Field");
  149. _xplot->SetXValuesToValue();
  150. //_plot->SetXValuesToIndex();
  151. // set which parts of the data object are to be used for which axis
  152. _xplot->SetDataObjectXComponent(0,1);
  153. _xplot->SetDataObjectYComponent(0,0);
  154. _xplot->SetDataObjectXComponent(1,1);
  155. _xplot->SetDataObjectYComponent(1,0);
  156. _xplot->SetDataObjectXComponent(2,1);
  157. _xplot->SetDataObjectYComponent(2,0);
  158. _xplot->SetDataObjectXComponent(3,1);
  159. _xplot->SetDataObjectYComponent(3,0);
  160. _xplot->SetNumberOfXLabels(3);
  161. _xplot->SetPlotColor(0,1,1,1);
  162. _xplot->SetPlotColor(1,1,1,1);
  163. _xplot->SetPlotColor(2,1,0,0);
  164. _xplot->SetPlotColor(3,1,0,0);
  165. _xplot->GetProperty()->SetLineWidth(2);
  166. _xplot->GetProperty()->SetPointSize(4);
  167. // _plot->SetPlotPoints(0,0);
  168. // _plot->SetPlotPoints(1,0);
  169. // _plot->SetPlotPoints(2,0);
  170. // _plot->SetPlotPoints(3,0);
  171. // _plot->SetPlotLines(0,0);
  172. // _plot->SetPlotLines(1,0);
  173. // _plot->SetPlotLines(2,0);
  174. // _plot->SetPlotLines(3,0);
  175. _xplot->PlotCurvePointsOff();
  176. _xplot->PlotCurveLinesOff();
  177. _xplot->GetPositionCoordinate()->SetValue(0.0, 0.67, 0);
  178. _xplot->GetPosition2Coordinate()->SetValue(1.0, 0.33, 0);
  179. //_plot->SetReverseYAxis(1); // Just flips axis, not data!
  180. ///////////////////////////////////////////////
  181. // Y Component Plot
  182. vtkXYPlotActor *_yplot = vtkXYPlotActor::New();
  183. //vtkXYPlotWidget *_wplot = vtkXYPlotWidget::New();
  184. // Add the datasets
  185. _yplot->AddDataObjectInput( _dataObjectyr);
  186. _yplot->AddDataObjectInput( _dataObjectyi);
  187. //_yplot->AddDataObjectInput(_fdataObjectyr);
  188. //_yplot->AddDataObjectInput(_fdataObjectyi);
  189. _yplot->SetTitle("Y field");
  190. _yplot->SetXTitle("Depth");
  191. _yplot->SetYTitle("E Field");
  192. _yplot->SetXValuesToValue();
  193. //_yplot->SetXValuesToIndex();
  194. // set which parts of the data object are to be used for which axis
  195. _yplot->SetDataObjectXComponent(0,1);
  196. _yplot->SetDataObjectYComponent(0,0);
  197. _yplot->SetDataObjectXComponent(1,1);
  198. _yplot->SetDataObjectYComponent(1,0);
  199. _yplot->SetDataObjectXComponent(2,1);
  200. _yplot->SetDataObjectYComponent(2,0);
  201. _yplot->SetDataObjectXComponent(3,1);
  202. _yplot->SetDataObjectYComponent(3,0);
  203. _yplot->SetNumberOfXLabels(3);
  204. _yplot->SetPlotColor(0,1,1,1);
  205. _yplot->SetPlotColor(1,1,1,1);
  206. _yplot->SetPlotColor(2,1,0,0);
  207. _yplot->SetPlotColor(3,1,0,0);
  208. _yplot->GetProperty()->SetLineWidth(2);
  209. _yplot->GetProperty()->SetPointSize(4);
  210. _yplot->GetPositionCoordinate()->SetValue(0.00, 0.33, 0);
  211. _yplot->GetPosition2Coordinate()->SetValue(1.0, 0.33, 0);
  212. // _plot->SetPlotPoints(0,0);
  213. // _plot->SetPlotPoints(1,0);
  214. // _plot->SetPlotPoints(2,0);
  215. // _plot->SetPlotPoints(3,0);
  216. // _plot->SetPlotLines(0,0);
  217. // _plot->SetPlotLines(1,0);
  218. // _plot->SetPlotLines(2,0);
  219. // _plot->SetPlotLines(3,0);
  220. _yplot->PlotCurvePointsOff();
  221. _yplot->PlotCurveLinesOff();
  222. //_plot->SetReverseYAxis(1)
  223. ///////////////////////////////////////////////
  224. // Z Component Plot
  225. vtkXYPlotActor *_zplot = vtkXYPlotActor::New();
  226. //vtkXYPlotWidget *_wplot = vtkXYPlotWidget::New();
  227. // Add the datasets
  228. _zplot->AddDataObjectInput( _dataObjectzr);
  229. _zplot->AddDataObjectInput( _dataObjectzi);
  230. //_zplot->AddDataObjectInput(_fdataObjectzr);
  231. //_zplot->AddDataObjectInput(_fdataObjectzi);
  232. _zplot->SetTitle("Z field");
  233. _zplot->SetXTitle("Depth");
  234. _zplot->SetYTitle("E Field");
  235. _zplot->SetXValuesToValue();
  236. //_plot->SetXValuesToIndex();
  237. // set which parts of the data object are to be used for which axis
  238. _zplot->SetDataObjectXComponent(0,1);
  239. _zplot->SetDataObjectYComponent(0,0);
  240. _zplot->SetDataObjectXComponent(1,1);
  241. _zplot->SetDataObjectYComponent(1,0);
  242. _zplot->SetDataObjectXComponent(2,1);
  243. _zplot->SetDataObjectYComponent(2,0);
  244. _zplot->SetDataObjectXComponent(3,1);
  245. _zplot->SetDataObjectYComponent(3,0);
  246. _zplot->SetNumberOfXLabels(3);
  247. _zplot->SetPlotColor(0,1,1,1);
  248. _zplot->SetPlotColor(1,1,1,1);
  249. _zplot->SetPlotColor(2,1,0,0);
  250. _zplot->SetPlotColor(3,1,0,0);
  251. _zplot->GetProperty()->SetLineWidth(2);
  252. _zplot->GetProperty()->SetPointSize(4);
  253. _zplot->GetPositionCoordinate()->SetValue(0.0, 0.0, 0);
  254. _zplot->GetPosition2Coordinate()->SetValue(1.0, 0.33, 0);
  255. // _plot->SetPlotPoints(0,0);
  256. // _plot->SetPlotPoints(1,0);
  257. // _plot->SetPlotPoints(2,0);
  258. // _plot->SetPlotPoints(3,0);
  259. // _plot->SetPlotLines(0,0);
  260. // _plot->SetPlotLines(1,0);
  261. // _plot->SetPlotLines(2,0);
  262. // _plot->SetPlotLines(3,0);
  263. _zplot->PlotCurvePointsOff();
  264. _zplot->PlotCurveLinesOff();
  265. //_plot->SetReverseYAxis(1)
  266. _ren->AddActor2D(_xplot);
  267. _ren->AddActor2D(_yplot);
  268. _ren->AddActor2D(_zplot);
  269. vtkRenderWindowInteractor *_iren = vtkRenderWindowInteractor::New();
  270. _iren->SetRenderWindow(_renwin);
  271. _iren->Initialize();
  272. _renwin->Render();
  273. //_wplot->SetXYPlotActor(_plot);
  274. //_wplot->SetInteractor(_iren);
  275. //_wplot->EnabledOn();
  276. _iren->Start();
  277. // TODO clean up you dirty whore
  278. // _iren->Delete();
  279. // _plot->Delete();
  280. // _dataObject->Delete();
  281. // _fieldData->Delete();
  282. // _renwin->Delete();
  283. // _ren->Delete();
  284. // _xdatareal->Delete();
  285. // _ydata->Delete();
  286. // _depth->Delete();
  287. #endif
  288. EmEarth->Delete();
  289. receivers->Delete();
  290. earth->Delete();
  291. dipole->Delete();
  292. return 0;
  293. }