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.

FieldPoints.cpp 16KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541
  1. /* This file is part of Lemma, a geophysical modelling and inversion API */
  2. /* This Source Code Form is subject to the terms of the Mozilla Public
  3. * License, v. 2.0. If a copy of the MPL was not distributed with this
  4. * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
  5. /**
  6. @file
  7. @author Trevor Irons
  8. @date 12/02/2009
  9. @version $Id: receiverpoints.cpp 203 2015-01-09 21:19:04Z tirons $
  10. **/
  11. #include "FieldPoints.h"
  12. namespace Lemma {
  13. // ==================== FRIENDS ======================
  14. std::ostream &operator << (std::ostream &stream, const FieldPoints &ob) {
  15. stream << ob.Serialize() << "\n";
  16. return stream;
  17. }
  18. // ==================== LIFECYCLE ===================================
  19. FieldPoints::FieldPoints( const ctor_key& key ) : LemmaObject( key ),
  20. NumberOfPoints(0), NumberOfBinsE(0), NumberOfBinsH(0) {
  21. }
  22. //--------------------------------------------------------------------------------------
  23. // Class: FieldPoints
  24. // Method: FieldPoints
  25. // Description: constructor (protected)
  26. //--------------------------------------------------------------------------------------
  27. FieldPoints::FieldPoints (const YAML::Node& node, const ctor_key& key) : LemmaObject(node, key) {
  28. //DeSerialize
  29. NumberOfPoints = node["NumberOfPoints"].as<int>();
  30. NumberOfBinsE = node["NumberOfBinsE"].as<int>();
  31. NumberOfBinsH = node["NumberOfBinsH"].as<int>();
  32. Mask = node["Mask"].as<VectorXi>();
  33. Locations = node["Locations"].as<Vector3Xr>();
  34. } // ----- end of method FieldPoints::FieldPoints (constructor) -----
  35. FieldPoints::~FieldPoints() {
  36. }
  37. std::shared_ptr<FieldPoints> FieldPoints::NewSP() {
  38. return std::make_shared<FieldPoints> ( ctor_key() );
  39. }
  40. //--------------------------------------------------------------------------------------
  41. // Class: FieldPoints
  42. // Method: Serialize
  43. //--------------------------------------------------------------------------------------
  44. YAML::Node FieldPoints::Serialize ( ) const {
  45. YAML::Node node = LemmaObject::Serialize();
  46. node.SetTag( GetName() );
  47. // update here
  48. node["NumberOfPoints"] = NumberOfPoints;
  49. node["NumberOfBinsE"] = NumberOfBinsE;
  50. node["NumberOfBinsH"] = NumberOfBinsH;
  51. node["Mask"] = Mask;
  52. node["Locations"] = Locations;// Can be huge
  53. //std::cout << "Locations.data" << Locations.data()[0] << std::endl;
  54. return node;
  55. } // ----- end of method FieldPoints::Serialize -----
  56. //--------------------------------------------------------------------------------------
  57. // Class: FieldPoints
  58. // Method: DeSerialize
  59. //--------------------------------------------------------------------------------------
  60. std::shared_ptr<FieldPoints> FieldPoints::DeSerialize ( const YAML::Node& node ) {
  61. if (node.Tag() != "FieldPoints") {
  62. throw DeSerializeTypeMismatch( "FieldPoints", node.Tag());
  63. }
  64. return std::make_shared<FieldPoints> ( node, ctor_key() );
  65. } // ----- end of method FieldPoints::DeSerialize -----
  66. //--------------------------------------------------------------------------------------
  67. // Class: FieldPoints
  68. // Method: GetName
  69. // Description: Class identifier
  70. //--------------------------------------------------------------------------------------
  71. inline std::string FieldPoints::GetName ( ) const {
  72. return CName;
  73. } // ----- end of method FieldPoints::GetName -----
  74. // ==================== ACCESS ===================================
  75. void FieldPoints::SetNumberOfPoints(const int &nrec) {
  76. if (nrec > 0)
  77. this->NumberOfPoints = nrec;
  78. else
  79. throw std::runtime_error("NUMBER RECEIVERS LESS THAN 1");
  80. this->Locations.resize(Eigen::NoChange, nrec);
  81. Locations.setZero();
  82. this->Mask.resize(nrec);
  83. Mask.setZero();
  84. ResizeEField();
  85. ResizeHField();
  86. }
  87. void FieldPoints::ResizeEField() {
  88. Efield.clear();
  89. for (int i=0; i<NumberOfBinsE; ++i) {
  90. Eigen::Matrix<Complex, 3, Eigen::Dynamic> tempe;
  91. this->Efield.push_back(tempe);
  92. this->Efield[i].resize(Eigen::NoChange, NumberOfPoints);
  93. this->Efield[i].setZero();
  94. }
  95. }
  96. void FieldPoints::ResizeHField() {
  97. Hfield.clear();
  98. for (int i=0; i<NumberOfBinsH; ++i) {
  99. Eigen::Matrix<Complex, 3, Eigen::Dynamic> temph;
  100. this->Hfield.push_back(temph);
  101. this->Hfield[i].resize(Eigen::NoChange, NumberOfPoints);
  102. this->Hfield[i].setZero();
  103. }
  104. }
  105. void FieldPoints::SetNumberOfBinsE(const int& nbins) {
  106. NumberOfBinsE = nbins;
  107. ResizeEField();
  108. }
  109. void FieldPoints::SetNumberOfBinsH(const int& nbins) {
  110. NumberOfBinsH = nbins;
  111. ResizeHField();
  112. }
  113. void FieldPoints::SetLocation(const int&nrec,const Vector3r& loc) {
  114. this->Locations.col(nrec) = loc;
  115. }
  116. void FieldPoints::SetLocation(const int&nrec,const Real& xp,
  117. const Real& yp, const Real& zp) {
  118. this->Locations.col(nrec) << xp, yp, zp;
  119. }
  120. void FieldPoints::SetEfield(const int& nbin,
  121. const int& loc, const Complex &ex,
  122. const Complex &ey, const Complex &ez) {
  123. this->Efield[nbin].col(loc) << ex, ey, ez;
  124. }
  125. void FieldPoints::AppendEfield(const int&nbin, const int& loc,
  126. const Complex &ex,
  127. const Complex &ey, const Complex &ez) {
  128. #ifdef LEMMAUSEOMP
  129. #pragma omp critical
  130. #endif
  131. this->Efield[nbin].col(loc) += Vector3cr(ex, ey, ez); //temp;
  132. }
  133. void FieldPoints::SetHfield(const int &nbin, const int& loc,
  134. const Complex &hx, const Complex &hy,
  135. const Complex &hz) {
  136. this->Hfield[nbin].col(loc) << hx, hy, hz;
  137. }
  138. void FieldPoints::AppendHfield(const int &nbin, const int& loc,
  139. const Complex &hx, const Complex &hy,
  140. const Complex &hz) {
  141. // #ifdef LEMMAUSEOMP
  142. // #pragma omp atomic
  143. // std::real(this->Hfield[nbin].col(loc)[0]) += std::real(hx);
  144. // #pragma omp atomic
  145. // std::imag(this->Hfield[nbin].col(loc)[0]) += std::imag(hx);
  146. // #pragma omp atomic
  147. // std::real(this->Hfield[nbin].col(loc)[1]) += std::real(hy);
  148. // #pragma omp atomic
  149. // std::imag(this->Hfield[nbin].col(loc)[1]) += std::imag(hy);
  150. // #pragma omp atomic
  151. // std::real(this->Hfield[nbin].col(loc)[2]) += std::real(hz);
  152. // #pragma omp atomic
  153. // std::imag(this->Hfield[nbin].col(loc)[2]) += std::imag(hz);
  154. // #else
  155. //(critical sections are slow)
  156. #ifdef LEMMAUSEOMP
  157. #pragma omp critical
  158. #endif
  159. this->Hfield[nbin].col(loc) += Vector3cr(hx,hy,hz);
  160. //#endif
  161. }
  162. // ==================== INQUIRY ===================================
  163. Vector3Xr FieldPoints::GetLocations() {
  164. return this->Locations;
  165. }
  166. MatrixXr FieldPoints::GetLocationsMat() {
  167. return MatrixXr(this->Locations);
  168. }
  169. Vector3r FieldPoints::GetLocation(const int &nrec) {
  170. return this->Locations.col(nrec);
  171. }
  172. Real FieldPoints::GetLocationX(const int &nrec) {
  173. return this->Locations.col(nrec)[0];
  174. }
  175. Real FieldPoints::GetLocationY(const int &nrec) {
  176. return this->Locations.col(nrec)[1];
  177. }
  178. Real FieldPoints::GetLocationZ(const int &nrec) {
  179. return this->Locations.col(nrec)[2];
  180. }
  181. Vector3cr FieldPoints::GetEfield(const int &nfreq, const int&nrec) {
  182. return this->Efield[nfreq].col(nrec);
  183. }
  184. Vector3cr FieldPoints::GetHfield(const int &nfreq, const int&nrec) {
  185. return this->Hfield[nfreq].col(nrec);
  186. }
  187. std::vector<Vector3Xcr> FieldPoints::GetHfield( ) {
  188. return this->Hfield;
  189. }
  190. std::vector<Vector3Xcr> FieldPoints::GetEfield( ) {
  191. return this->Efield;
  192. }
  193. Vector3Xcr FieldPoints::GetEfield (const int &nfreq) {
  194. return this->Efield[nfreq];
  195. }
  196. MatrixXcr FieldPoints::GetEfieldMat (const int &nfreq) {
  197. return MatrixXcr(this->Efield[nfreq]);
  198. }
  199. MatrixXcr FieldPoints::GetHfieldMat (const int &nfreq) {
  200. return MatrixXcr(this->Hfield[nfreq]);
  201. }
  202. Vector3Xcr FieldPoints::GetHfield (const int &nfreq) {
  203. return this->Hfield[nfreq];
  204. }
  205. void FieldPoints::MaskPoint(const int& imask) {
  206. Mask(imask) = true;
  207. }
  208. void FieldPoints::UnMaskPoint(const int& imask) {
  209. Mask(imask) = false;
  210. }
  211. void FieldPoints::UnMaskAllPoints() {
  212. Mask.setZero();
  213. }
  214. int FieldPoints::GetMask(const int& i) {
  215. return Mask(i);
  216. }
  217. int FieldPoints::GetNumberOfPoints() {
  218. return this->NumberOfPoints;
  219. }
  220. void FieldPoints::ClearFields() {
  221. for (int i=0; i<NumberOfBinsE; ++i) {
  222. this->Efield[i].setZero();
  223. }
  224. for (int i=0; i<NumberOfBinsH; ++i) {
  225. this->Hfield[i].setZero();
  226. }
  227. }
  228. #ifdef LEMMAUSEVTK
  229. vtkActor* FieldPoints::GetVtkGlyphActor(const FIELDTYPE &ftype,
  230. const Real &clip, const Real &scale,
  231. const int &nfreq) {
  232. vtkArrowSource *vArrow;
  233. vtkGlyph3D *vGlyph;
  234. vtkActor *vActor;
  235. vtkPolyDataMapper *vPolyMapper;
  236. vtkPoints *vPoints;
  237. vtkDoubleArray *vVects;
  238. vtkPolyData *vPointSet;
  239. vActor = vtkActor::New();
  240. vGlyph = vtkGlyph3D::New();
  241. vArrow = vtkArrowSource::New();
  242. vPolyMapper = vtkPolyDataMapper::New();
  243. vPoints = vtkPoints::New();
  244. vPointSet = vtkPolyData::New();
  245. vVects = vtkDoubleArray::New();
  246. vVects->SetNumberOfComponents(3);
  247. // Make PointData
  248. for (int ic=0; ic<NumberOfPoints; ++ic) {
  249. Vector3r loc = this->GetLocation(ic);
  250. vPoints->InsertPoint(ic, loc[0], loc[1],
  251. loc[2]);
  252. Vector3r temp;
  253. temp.setZero();
  254. switch (ftype) {
  255. case HFIELDREAL:
  256. temp = scale*(GetHfield(nfreq,ic).real());
  257. break;
  258. case HFIELDIMAG:
  259. temp = scale*(GetHfield(nfreq,ic).imag());
  260. break;
  261. case EFIELDREAL:
  262. temp = scale*(GetEfield(nfreq, ic).real());
  263. break;
  264. case EFIELDIMAG:
  265. temp = scale*(GetEfield(nfreq, ic).imag());
  266. break;
  267. }
  268. if (temp.norm() > clip) {
  269. temp /= temp.norm(); // norm
  270. temp *= clip; // clip dimension
  271. }
  272. vVects->InsertTuple3(ic, temp(0), temp(1),
  273. temp(2));
  274. }
  275. vPointSet->SetPoints(vPoints);
  276. vPointSet->GetPointData()->SetVectors(vVects);
  277. vGlyph->ScalingOn();
  278. vGlyph->SetScaleModeToScaleByVector();
  279. vGlyph->SetSourceConnection(vArrow->GetOutputPort());
  280. vGlyph->SetVectorMode(true);
  281. vGlyph->SetVectorModeToUseVector();
  282. vGlyph->OrientOn();
  283. vGlyph->SetInputData(vPointSet);
  284. vPolyMapper->SetInputConnection(vGlyph->GetOutputPort());
  285. vActor->SetMapper(vPolyMapper);
  286. if (vArrow != NULL) {
  287. vArrow->Delete();
  288. vArrow = NULL;
  289. }
  290. //if (vActor != NULL) {
  291. // vActor->Delete();
  292. // vActor = NULL;
  293. //}
  294. if (vPolyMapper != NULL) {
  295. vPolyMapper->Delete();
  296. vPolyMapper = NULL;
  297. }
  298. if (vVects != NULL) {
  299. vVects->Delete();
  300. vVects = NULL;
  301. }
  302. if (vPointSet != NULL) {
  303. vPointSet->Delete();
  304. vPointSet = NULL;
  305. }
  306. if (vPoints != NULL) {
  307. vPoints->Delete();
  308. vPoints = NULL;
  309. }
  310. if (vGlyph != NULL) {
  311. vGlyph->Delete();
  312. vGlyph = NULL;
  313. }
  314. return vActor;
  315. }
  316. vtkDataObject* FieldPoints::GetVtkDataObject(const FIELDTYPE &ftype,
  317. const int& nbin,
  318. const int& start, const int& end,
  319. const FIELDCOMPONENT& fcomp,
  320. const SPATIALCOORDINANT &scord) {
  321. if (start < 0) throw 77;
  322. if (end > NumberOfPoints) throw 78;
  323. if (start > end) throw 79;
  324. int ifc(-1);
  325. switch (fcomp) {
  326. case XCOMPONENT:
  327. ifc = 0;
  328. break;
  329. case YCOMPONENT:
  330. ifc = 1;
  331. break;
  332. case ZCOMPONENT:
  333. ifc = 2;
  334. break;
  335. }
  336. int isc(-1);
  337. switch (scord) {
  338. case XCOORD:
  339. isc = 0;
  340. break;
  341. case YCOORD:
  342. isc = 1;
  343. break;
  344. case ZCOORD:
  345. isc = 2;
  346. break;
  347. }
  348. vtkDoubleArray *_data = vtkDoubleArray::New();
  349. _data->SetNumberOfComponents(1);
  350. vtkDoubleArray *_pos = vtkDoubleArray::New();
  351. _pos->SetNumberOfComponents(1);
  352. int id=0;
  353. for (int irec=start; irec<end; ++irec) {
  354. switch (ftype) {
  355. case HFIELDREAL:
  356. _data->InsertTuple1(id, std::real(Hfield[nbin](ifc, irec)));
  357. break;
  358. case HFIELDIMAG:
  359. _data->InsertTuple1(id, std::imag(Hfield[nbin](ifc, irec)));
  360. break;
  361. case EFIELDREAL:
  362. _data->InsertTuple1(id, std::real(Efield[nbin](ifc, irec)));
  363. break;
  364. case EFIELDIMAG:
  365. _data->InsertTuple1(id, std::imag(Efield[nbin](ifc, irec)));
  366. break;
  367. }
  368. _pos->InsertTuple1(id, Locations(isc, irec));
  369. ++id;
  370. }
  371. vtkFieldData *_fieldData = vtkFieldData::New();
  372. _fieldData->AllocateArrays(2);
  373. _fieldData->AddArray(_data);
  374. _fieldData->AddArray(_pos);
  375. vtkDataObject *_dataObject = vtkDataObject::New();
  376. _dataObject->SetFieldData(_fieldData);
  377. _data->Delete();
  378. _pos->Delete();
  379. _fieldData->Delete();
  380. return _dataObject;
  381. }
  382. vtkDataObject* FieldPoints::GetVtkDataObjectFreq(const FIELDTYPE &ftype,
  383. const int& nrec,
  384. const int& fstart, const int& fend,
  385. const FIELDCOMPONENT& fcomp,
  386. const VectorXr& Freqs) {
  387. if (fstart < 0) throw 77;
  388. //if (fend > NumberOfFrequencies) throw 78;
  389. if (fstart > fend) throw 79;
  390. int ifc(-1);
  391. switch (fcomp) {
  392. case XCOMPONENT:
  393. ifc = 0;
  394. break;
  395. case YCOMPONENT:
  396. ifc = 1;
  397. break;
  398. case ZCOMPONENT:
  399. ifc = 2;
  400. break;
  401. }
  402. vtkDoubleArray *_data = vtkDoubleArray::New();
  403. _data->SetNumberOfComponents(1);
  404. vtkDoubleArray *_pos = vtkDoubleArray::New();
  405. _pos->SetNumberOfComponents(1);
  406. int id=0;
  407. //std::cout.precision(12);
  408. for (int ifreq=fstart; ifreq<fend; ++ifreq) {
  409. switch (ftype) {
  410. case HFIELDREAL:
  411. _data->InsertTuple1(id, std::real(Hfield[ifreq](ifc, nrec)));
  412. //std::cout << Hfield[ifreq](ifc, nrec) << std::endl;
  413. break;
  414. case HFIELDIMAG:
  415. _data->InsertTuple1(id, std::imag(Hfield[ifreq](ifc, nrec)));
  416. break;
  417. case EFIELDREAL:
  418. _data->InsertTuple1(id, std::real(Efield[ifreq](ifc, nrec)));
  419. break;
  420. case EFIELDIMAG:
  421. _data->InsertTuple1(id, std::imag(Efield[ifreq](ifc, nrec)));
  422. break;
  423. }
  424. _pos->InsertTuple1(id, Freqs[ifreq]);
  425. ++id;
  426. }
  427. vtkFieldData *_fieldData = vtkFieldData::New();
  428. _fieldData->AllocateArrays(2);
  429. _fieldData->AddArray(_data);
  430. _fieldData->AddArray(_pos);
  431. vtkDataObject *_dataObject = vtkDataObject::New();
  432. _dataObject->SetFieldData(_fieldData);
  433. _data->Delete();
  434. _pos->Delete();
  435. _fieldData->Delete();
  436. return _dataObject;
  437. }
  438. #endif
  439. }