Surface NMR forward modelling
Ви не можете вибрати більше 25 тем Теми мають розпочинатися з літери або цифри, можуть містити дефіси (-) і не повинні перевищувати 35 символів.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645
  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 11/11/2016 01:47:34 PM
  11. * @author Trevor Irons (ti)
  12. * @email tirons@egi.utah.edu
  13. * @copyright Copyright (c) 2016, University of Utah
  14. * @copyright Copyright (c) 2016, Lemma Software, LLC
  15. * @copyright Copyright (c) 2008, Colorado School of Mines
  16. */
  17. #ifndef LOOPINTERACTIONS
  18. #define LOOPINTERACTIONS
  19. #pragma once
  20. #include "LemmaObject.h"
  21. #include "LayeredEarthEM.h"
  22. #include "PolygonalWireAntenna.h"
  23. #include "EMEarth1D.h"
  24. #include "FieldPoints.h"
  25. #ifdef LEMMAUSEVTK
  26. #include "vtkHyperOctree.h"
  27. #include "vtkHyperOctreeCursor.h"
  28. #include "vtkXMLHyperOctreeWriter.h"
  29. #include "vtkDoubleArray.h"
  30. #endif
  31. namespace Lemma {
  32. enum INTERACTION {COUPLING, INTERFERENCE, PHASE};
  33. /// convert enums to string saves repeated code useful for YAML serializing
  34. std::string enum2String(const INTERACTION& type);
  35. /**
  36. * \ingroup Merlin
  37. * \brief
  38. * \details
  39. */
  40. template< INTERACTION Type >
  41. class LoopInteractions : public LemmaObject {
  42. friend std::ostream &operator << (std::ostream &stream, const LoopInteractions &ob) {
  43. stream << ob.Serialize() << "\n---\n"; // End of doc ---
  44. return stream;
  45. }
  46. protected:
  47. /*
  48. * This key is used to lock the constructor. It is protected so that inhereted
  49. * classes also have the key to contruct their base class.
  50. */
  51. struct ctor_key {};
  52. public:
  53. // ==================== LIFECYCLE =======================
  54. /**
  55. * Default constructor.
  56. * @note This method is locked, and cannot be called directly.
  57. * The reason that the method is public is to enable the use
  58. * of make_shared whilst enforcing the use of shared_ptr,
  59. * in c++-17, this curiosity may be resolved.
  60. * @see LoopInteractions::NewSP
  61. */
  62. explicit LoopInteractions ( const ctor_key& ) : LemmaObject () { }
  63. /**
  64. * DeSerializing constructor.
  65. * @note This method is locked, and cannot be called directly.
  66. * The reason that the method is public is to enable the use
  67. * of make_shared whilst enforcing the use of shared_ptr,
  68. * in c++-17, this curiosity may be resolved.
  69. * @see LoopInteractions::DeSerialize
  70. */
  71. LoopInteractions ( const YAML::Node& node, const ctor_key& ) : LemmaObject(node) { }
  72. /**
  73. * Default destructor.
  74. * @note This method should never be called due to the mandated
  75. * use of smart pointers. It is necessary to keep the method
  76. * public in order to allow for the use of the more efficient
  77. * make_shared constructor.
  78. */
  79. virtual ~LoopInteractions () { }
  80. /**
  81. * Uses YAML to serialize this object.
  82. * @return a YAML::Node
  83. * @see LoopInteractions::DeSerialize
  84. */
  85. virtual YAML::Node Serialize() const {
  86. YAML::Node node = LemmaObject::Serialize();
  87. node.SetTag( GetName() );
  88. // Coils Transmitters & Receivers
  89. for ( auto txm : TxRx) {
  90. node[txm.first] = txm.second->Serialize();
  91. }
  92. // LayeredEarthEM
  93. node["SigmaModel"] = SigmaModel->Serialize();
  94. node["tol"] = tol;
  95. node["minLevel"] = minLevel;
  96. node["maxLevel"] = maxLevel;
  97. return node;
  98. }
  99. /*
  100. * Factory method for generating concrete class.
  101. * @return a std::shared_ptr of type LoopInteractions
  102. */
  103. static std::shared_ptr< LoopInteractions > NewSP() {
  104. return std::make_shared< LoopInteractions >( ctor_key() );
  105. }
  106. /**
  107. * Constructs an LoopInteractions object from a YAML::Node.
  108. * @see LoopInteractions::Serialize
  109. */
  110. static std::shared_ptr<LoopInteractions> DeSerialize(const YAML::Node& node) {
  111. if (node.Tag() != "LoopInteractions" ) {
  112. throw DeSerializeTypeMismatch( "LoopInteractions", node.Tag());
  113. }
  114. return std::make_shared< LoopInteractions > ( node, ctor_key() );
  115. } // ----- end of method LoopInteractions::DeSerialize -----
  116. // ==================== OPERATORS =======================
  117. // ==================== OPERATIONS =======================
  118. /**
  119. * @return std::shared_ptr<LayeredEarthEM>
  120. */
  121. inline std::shared_ptr<LayeredEarthEM> GetSigmaModel ( ) {
  122. return SigmaModel;
  123. } // ----- end of method LoopInteractions::get_SigmaModel -----
  124. /**
  125. * @param[in] value the 1D-EM model used for calculations
  126. */
  127. inline void SetLayeredEarthEM ( std::shared_ptr< LayeredEarthEM > value ) {
  128. SigmaModel = value;
  129. return ;
  130. } // ----- end of method LoopInteractions::set_SigmaModel -----
  131. /**
  132. *
  133. */
  134. inline void SetIntegrationSize ( const Vector3r& size ) {
  135. Size = size;
  136. return ;
  137. } // ----- end of method LoopInteractions::SetIntegrationSize -----
  138. /**
  139. *
  140. */
  141. inline void SetIntegrationOrigin ( const Vector3r& origin ) {
  142. Origin = origin;
  143. return ;
  144. } // ----- end of method LoopInteractions::SetIntegrationOrigin -----
  145. /**
  146. * Assign transmiter coils
  147. */
  148. inline void PushCoil( const std::string& label, std::shared_ptr<PolygonalWireAntenna> ant ) {
  149. TxRx[label] = ant;
  150. }
  151. /**
  152. * Calculates a single imaging kernel, however, phased arrays are supported
  153. * so that more than one transmitter and/or receiver can be specified.
  154. * @param[in] tx is the list of transmitters to use for a kernel, use the same labels as
  155. * used in PushCoil.
  156. * @param[in] rx is the list of receivers to use for a kernel, use the same labels as
  157. * used in PushCoil. @see PushCoil
  158. * @param[in] vtkOutput generates a VTK hyperoctree file as well, useful for visualization.
  159. * requires compilation of Lemma with VTK.
  160. */
  161. Complex Calculate (const std::vector< std::string >& tx, const std::vector< std::string >& rx,
  162. bool vtkOutput=false );
  163. /**
  164. * Sets the tolerance to use for making the adaptive mesh
  165. *
  166. */
  167. inline void SetTolerance(const Real& ttol) {
  168. tol = ttol;
  169. }
  170. inline void SetPulseDuration(const Real& taup) {
  171. Taup = taup;
  172. }
  173. // ==================== INQUIRY =======================
  174. /**
  175. * Returns the name of the underlying class, similiar to Python's type
  176. * @return string of class name
  177. */
  178. virtual inline std::string GetName() const {
  179. return CName;
  180. }
  181. protected:
  182. // ==================== LIFECYCLE =======================
  183. /** Copy is disabled */
  184. LoopInteractions( const LoopInteractions& ) = delete;
  185. private:
  186. /**
  187. * Returns the kernel value for an input prism
  188. */
  189. virtual Complex f( const Vector3r& r, const Real& volume , const Vector3cr& Ht, const Vector3cr& Hr);
  190. void IntegrateOnOctreeGrid( bool vtkOutput=false );
  191. /**
  192. * Recursive call to integrate a function on an adaptive Octree Grid.
  193. * For efficiency's sake the octree grid is not stored, as only the
  194. * integral (sum) is of interest. The logic for grid refinement is based
  195. * on an Octree representation of the domain. If an Octree representation
  196. * of the kernel is desired, call alternative version @see EvaluateKids2
  197. * @param[in] size gives the domain size, in metres
  198. * @param[in] level gives the current level of the octree grid, call with 0 initially
  199. * @param[in] cpos is the centre position of the parent cuboid
  200. */
  201. void EvaluateKids( const Vector3r& size, const int& level, const Vector3r& cpos,
  202. const Complex& parentVal );
  203. #ifdef LEMMAUSEVTK
  204. /**
  205. * Same functionality as @see EvaluateKids, but includes generation of a VTK
  206. * HyperOctree, which is useful for visualization.
  207. */
  208. void EvaluateKids2( const Vector3r& size, const int& level, const Vector3r& cpos,
  209. const Complex& parentVal, vtkHyperOctree* octree, vtkHyperOctreeCursor* curse );
  210. void GetPosition( vtkHyperOctreeCursor* Cursor, Real* p );
  211. #endif
  212. // ==================== DATA MEMBERS =========================
  213. int ilay;
  214. int nleaves;
  215. int minLevel=4;
  216. int maxLevel=8;
  217. Real VOLSUM;
  218. Real tol=1e-11;
  219. Real Taup = .020; // Sec
  220. Complex SUM;
  221. Vector3r Size;
  222. Vector3r Origin;
  223. std::shared_ptr< LayeredEarthEM > SigmaModel = nullptr;
  224. std::shared_ptr< FieldPoints > cpoints;
  225. std::map< std::string , std::shared_ptr< PolygonalWireAntenna > > TxRx;
  226. std::map< std::string , std::shared_ptr< EMEarth1D > > EMEarths;
  227. #ifdef LEMMAUSEVTK
  228. std::map< int, Complex > LeafDict;
  229. std::map< int, int > LeafDictIdx;
  230. std::map< int, Real > LeafDictErr;
  231. #endif
  232. /** ASCII string representation of the class name */
  233. static constexpr auto CName = "LoopInteractions";
  234. }; // ----- end of class LoopInteractions -----
  235. ///////////////////////////////////////////////////////////////
  236. // Implimentation of non specialized args -- templated class //
  237. ///////////////////////////////////////////////////////////////
  238. // forward declare specs
  239. template <>
  240. Complex LoopInteractions<COUPLING>::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr );
  241. template <>
  242. Complex LoopInteractions<INTERFERENCE>::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr );
  243. template <>
  244. Complex LoopInteractions<PHASE>::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr );
  245. //--------------------------------------------------------------------------------------
  246. // Class: LoopInteractions
  247. // Method: DeSerialize
  248. //--------------------------------------------------------------------------------------
  249. template< INTERACTION Type >
  250. Complex LoopInteractions<Type>::Calculate (const std::vector< std::string>& Tx, const std::vector<std::string >& Rx,
  251. bool vtkOutput ) {
  252. static bool first = false; // a little hackish
  253. if (!first) {
  254. // All EM calculations will share same field points
  255. cpoints = FieldPoints::NewSP();
  256. cpoints->SetNumberOfPoints(8);
  257. }
  258. first = true;
  259. for (auto tx : Tx) {
  260. // Set up EMEarth
  261. EMEarths[tx] = EMEarth1D::NewSP();
  262. EMEarths[tx]->AttachWireAntenna(TxRx[tx]);
  263. EMEarths[tx]->AttachLayeredEarthEM(SigmaModel);
  264. EMEarths[tx]->AttachFieldPoints( cpoints );
  265. EMEarths[tx]->SetFieldsToCalculate(H);
  266. // TODO query for method, altough with flat antennae, this is fastest
  267. EMEarths[tx]->SetHankelTransformMethod(ANDERSON801);
  268. EMEarths[tx]->SetTxRxMode(TX);
  269. TxRx[tx]->SetCurrent(1.);
  270. }
  271. for (auto rx : Rx) {
  272. if (EMEarths.count(rx)) {
  273. EMEarths[rx]->SetTxRxMode(TXRX);
  274. } else {
  275. EMEarths[rx] = EMEarth1D::NewSP();
  276. EMEarths[rx]->AttachWireAntenna(TxRx[rx]);
  277. EMEarths[rx]->AttachLayeredEarthEM(SigmaModel);
  278. EMEarths[rx]->AttachFieldPoints( cpoints );
  279. EMEarths[rx]->SetFieldsToCalculate(H);
  280. // TODO query for method, altough with flat antennae, this is fastest
  281. EMEarths[rx]->SetHankelTransformMethod(ANDERSON801);
  282. EMEarths[rx]->SetTxRxMode(RX);
  283. TxRx[rx]->SetCurrent(1.);
  284. }
  285. }
  286. SUM = 0;
  287. IntegrateOnOctreeGrid( vtkOutput );
  288. std::cout << "\nFinished KERNEL\n";
  289. EMEarths.clear();
  290. return SUM;
  291. }
  292. //--------------------------------------------------------------------------------------
  293. // Class: LoopInteractions
  294. // Method: IntegrateOnOctreeGrid
  295. //--------------------------------------------------------------------------------------
  296. template< INTERACTION Type >
  297. void LoopInteractions<Type>::IntegrateOnOctreeGrid( bool vtkOutput ) {
  298. static int count = 0;
  299. Vector3r cpos = Origin + Size/2.;
  300. VOLSUM = 0;
  301. nleaves = 0;
  302. if (!vtkOutput) {
  303. EvaluateKids( Size, 0, cpos, Complex(100.));
  304. } else {
  305. #ifdef LEMMAUSEVTK
  306. vtkHyperOctree* oct = vtkHyperOctree::New();
  307. oct->SetDimension(3);
  308. oct->SetOrigin( Origin(0), Origin(1), Origin(2) );
  309. oct->SetSize( Size(0), Size(1), Size(2) );
  310. vtkHyperOctreeCursor* curse = oct->NewCellCursor();
  311. curse->ToRoot();
  312. EvaluateKids2( Size, 0, cpos, Complex(100.0), oct, curse );
  313. // Fill in leaf data
  314. vtkDoubleArray* kr = vtkDoubleArray::New();
  315. kr->SetNumberOfComponents(1);
  316. kr->SetName("Re($\\sum$)");
  317. kr->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  318. vtkDoubleArray* ki = vtkDoubleArray::New();
  319. ki->SetNumberOfComponents(1);
  320. ki->SetName("Im($\\sum$)");
  321. ki->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  322. vtkDoubleArray* km = vtkDoubleArray::New();
  323. km->SetNumberOfComponents(1);
  324. km->SetName("mod($\\sum$)");
  325. km->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  326. vtkIntArray* kid = vtkIntArray::New();
  327. kid->SetNumberOfComponents(1);
  328. kid->SetName("ID");
  329. kid->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  330. vtkIntArray* kerr = vtkIntArray::New();
  331. kerr->SetNumberOfComponents(1);
  332. kerr->SetName("nleaf");
  333. //Real LeafVol(0);
  334. for (auto leaf : LeafDict) {
  335. kr->InsertTuple1( leaf.first, std::real(leaf.second) );
  336. ki->InsertTuple1( leaf.first, std::imag(leaf.second) );
  337. km->InsertTuple1( leaf.first, std::abs(leaf.second) );
  338. kid->InsertTuple1( leaf.first, leaf.first );
  339. //LeafVol += std::real(leaf.second);
  340. }
  341. //std::cout << "\n\nLeafVol=" << LeafVol << std::endl;
  342. for (auto leaf : LeafDictIdx) {
  343. kerr->InsertTuple1( leaf.first, leaf.second );
  344. }
  345. auto kri = oct->GetLeafData()->AddArray(kr);
  346. auto kii = oct->GetLeafData()->AddArray(ki);
  347. auto kmi = oct->GetLeafData()->AddArray(km);
  348. auto kidi = oct->GetLeafData()->AddArray(kid);
  349. auto keri = oct->GetLeafData()->AddArray(kerr);
  350. auto write = vtkXMLHyperOctreeWriter::New();
  351. //write.SetDataModeToAscii()
  352. write->SetInputData(oct);
  353. std::string fname = std::string("octree-") + enum2String(Type) + std::string("-")
  354. + to_string(count) + std::string(".vto");
  355. write->SetFileName(fname.c_str());
  356. write->Write();
  357. write->Delete();
  358. oct->GetLeafData()->RemoveArray( kri );
  359. oct->GetLeafData()->RemoveArray( kii );
  360. oct->GetLeafData()->RemoveArray( kmi );
  361. oct->GetLeafData()->RemoveArray( kidi );
  362. oct->GetLeafData()->RemoveArray( keri );
  363. kerr->Delete();
  364. kid->Delete();
  365. kr->Delete();
  366. ki->Delete();
  367. km->Delete();
  368. curse->Delete();
  369. oct->Delete();
  370. #else
  371. throw std::runtime_error("IntegrateOnOctreeGrid with vtkOutput requires Lemma with VTK support");
  372. #endif
  373. }
  374. std::cout << "\nVOLSUM=" << VOLSUM << "\tActual=" << Size(0)*Size(1)*Size(2)
  375. << "\tDifference=" << VOLSUM - (Size(0)*Size(1)*Size(2)) << std::endl;
  376. count += 1;
  377. }
  378. //--------------------------------------------------------------------------------------
  379. // Class: LoopInteractions
  380. // Method: EvaluateKids
  381. //--------------------------------------------------------------------------------------
  382. template< INTERACTION Type >
  383. void LoopInteractions<Type>::EvaluateKids( const Vector3r& size, const int& level, const Vector3r& cpos,
  384. const Complex& parentVal ) {
  385. std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
  386. std::cout.flush();
  387. // Next level step, interested in one level below
  388. // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
  389. Vector3r step = size.array() / (Real)(1 << (level+1) );
  390. Real vol = (step(0)*step(1)*step(2)); // volume of each child
  391. Vector3r pos = cpos - step/2.;
  392. Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
  393. 0, 0, 0,
  394. step[0], 0, 0,
  395. 0, step[1], 0,
  396. step[0], step[1], 0,
  397. 0, 0, step[2],
  398. step[0], 0, step[2],
  399. 0, step[1], step[2],
  400. step[0], step[1], step[2] ).finished();
  401. VectorXcr kvals(8); // individual kernel vals
  402. cpoints->ClearFields();
  403. for (int ichild=0; ichild<8; ++ichild) {
  404. Vector3r cp = pos; // Eigen complains about combining these
  405. cp += posadd.row(ichild);
  406. cpoints->SetLocation( ichild, cp );
  407. }
  408. Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
  409. Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
  410. for ( auto EMCalc : EMEarths ) {
  411. EMCalc.second->GetFieldPoints()->ClearFields();
  412. EMCalc.second->CalculateWireAntennaFields();
  413. switch (EMCalc.second->GetTxRxMode()) {
  414. case TX:
  415. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  416. break;
  417. case RX:
  418. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  419. break;
  420. case TXRX:
  421. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  422. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  423. break;
  424. default:
  425. break;
  426. }
  427. }
  428. for (int ichild=0; ichild<8; ++ichild) {
  429. Vector3r cp = pos; // Eigen complains about combining these
  430. cp += posadd.row(ichild);
  431. kvals(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
  432. }
  433. Complex ksum = kvals.sum(); // Kernel sum
  434. // Evaluate whether or not furthur splitting is needed
  435. if ( (std::abs(ksum-parentVal) > tol && level < maxLevel) || level < minLevel ) {
  436. // Not a leaf dive further in
  437. for (int ichild=0; ichild<8; ++ichild) {
  438. Vector3r cp = pos; // Eigen complains about combining these
  439. cp += posadd.row(ichild);
  440. EvaluateKids( size, level+1, cp, kvals(ichild) );
  441. }
  442. return; // not leaf
  443. }
  444. // implicit else, is a leaf
  445. SUM += ksum;
  446. VOLSUM += 8.*vol;
  447. nleaves += 1; // could say += 8 just as fairly
  448. return; // is leaf
  449. }
  450. #ifdef LEMMAUSEVTK
  451. //--------------------------------------------------------------------------------------
  452. // Class: LoopInteractions
  453. // Method: EvaluateKids2 -- same as Evaluate Kids, but include VTK octree generation
  454. //--------------------------------------------------------------------------------------
  455. template< INTERACTION Type >
  456. void LoopInteractions<Type>::EvaluateKids2( const Vector3r& size, const int& level, const Vector3r& cpos,
  457. const Complex& parentVal, vtkHyperOctree* oct, vtkHyperOctreeCursor* curse) {
  458. std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
  459. std::cout.flush();
  460. // Next level step, interested in one level below
  461. // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
  462. Vector3r step = size.array() / (Real)(1 << (level+1) );
  463. Real vol = (step(0)*step(1)*step(2)); // volume of each child
  464. Vector3r pos = cpos - step/2.;
  465. Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
  466. 0, 0, 0,
  467. step[0], 0, 0,
  468. 0, step[1], 0,
  469. step[0], step[1], 0,
  470. 0, 0, step[2],
  471. step[0], 0, step[2],
  472. 0, step[1], step[2],
  473. step[0], step[1], step[2] ).finished();
  474. VectorXcr kvals(8); // individual kernel vals
  475. cpoints->ClearFields();
  476. for (int ichild=0; ichild<8; ++ichild) {
  477. Vector3r cp = pos; // Eigen complains about combining these
  478. cp += posadd.row(ichild);
  479. cpoints->SetLocation( ichild, cp );
  480. }
  481. Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
  482. Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
  483. for ( auto EMCalc : EMEarths ) {
  484. //EMCalc->GetFieldPoints()->ClearFields();
  485. EMCalc.second->CalculateWireAntennaFields();
  486. switch (EMCalc.second->GetTxRxMode()) {
  487. case TX:
  488. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  489. break;
  490. case RX:
  491. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  492. break;
  493. case TXRX:
  494. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  495. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  496. break;
  497. default:
  498. break;
  499. }
  500. }
  501. for (int ichild=0; ichild<8; ++ichild) {
  502. Vector3r cp = pos; // Eigen complains about combining these
  503. cp += posadd.row(ichild);
  504. kvals(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
  505. }
  506. Complex ksum = kvals.sum(); // Kernel sum
  507. // Evaluate whether or not furthur splitting is needed
  508. if ( (std::abs(ksum-parentVal) > tol && level < maxLevel) || level < minLevel ) {
  509. oct->SubdivideLeaf(curse);
  510. for (int ichild=0; ichild<8; ++ichild) {
  511. curse->ToChild(ichild);
  512. Vector3r cp = pos; // Eigen complains about combining these
  513. cp += posadd.row(ichild);
  514. /* Test for position via alternative means */
  515. /*
  516. Real p[3];
  517. GetPosition(curse, p);
  518. if ( (Vector3r(p) - cp).norm() > 1e-8 ) {
  519. std::cout << "ERROR @ nleaves" << nleaves << "\n" << cp[0] << "\t" << p[0] << "\t" << cp[1] << "\t" << p[1]
  520. << "\t" << cp[2] << "\t" << p[2] << "\t" << vol<< std::endl;
  521. throw std::runtime_error("doom");
  522. }
  523. */
  524. /* End of position test */
  525. EvaluateKids2( size, level+1, cp, kvals(ichild), oct, curse );
  526. curse->ToParent();
  527. }
  528. return; // not a leaf
  529. }
  530. LeafDict[curse->GetLeafId()] = ksum/(8.*vol);
  531. LeafDictIdx[curse->GetLeafId()] = nleaves;
  532. SUM += ksum;
  533. VOLSUM += 8*vol;
  534. nleaves += 1;
  535. return; // is a leaf
  536. }
  537. //--------------------------------------------------------------------------------------
  538. // Class: LoopInteractions
  539. // Method: GetPosition
  540. //--------------------------------------------------------------------------------------
  541. template< INTERACTION Type >
  542. void LoopInteractions<Type>::GetPosition( vtkHyperOctreeCursor* Cursor, Real* p ) {
  543. Real ratio=1.0/(1<<(Cursor->GetCurrentLevel()));
  544. //step = ((Size).array() / std::pow(2.,Cursor->GetCurrentLevel()));
  545. p[0]=(Cursor->GetIndex(0)+.5)*ratio*this->Size[0]+this->Origin[0] ;//+ .5*step[0];
  546. p[1]=(Cursor->GetIndex(1)+.5)*ratio*this->Size[1]+this->Origin[1] ;//+ .5*step[1];
  547. p[2]=(Cursor->GetIndex(2)+.5)*ratio*this->Size[2]+this->Origin[2] ;//+ .5*step[2];
  548. }
  549. #endif
  550. } // ----- end of namespace Lemma ----
  551. /* vim: set tabstop=4 expandtab */
  552. /* vim: set filetype=cpp */
  553. #endif // ----- #ifndef LOOPINTERACTIONS -----