Surface NMR forward modelling
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.

LoopInteractions.h 26KB

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