Surface NMR forward modelling
Nelze vybrat více než 25 témat Téma musí začínat písmenem nebo číslem, může obsahovat pomlčky („-“) a může být dlouhé až 35 znaků.

KernelV0.cpp 28KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614
  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:25 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. #include "KernelV0.h"
  18. #include "FieldPoints.h"
  19. namespace Lemma {
  20. // ==================== FRIEND METHODS =====================
  21. std::ostream &operator << (std::ostream &stream, const KernelV0 &ob) {
  22. stream << ob.Serialize() << "\n---\n"; // End of doc ---
  23. return stream;
  24. }
  25. // ==================== LIFECYCLE =======================
  26. //--------------------------------------------------------------------------------------
  27. // Class: KernelV0
  28. // Method: KernelV0
  29. // Description: constructor (locked)
  30. //--------------------------------------------------------------------------------------
  31. KernelV0::KernelV0 (const ctor_key&) : LemmaObject( ) {
  32. } // ----- end of method KernelV0::KernelV0 (constructor) -----
  33. //--------------------------------------------------------------------------------------
  34. // Class: KernelV0
  35. // Method: KernelV0
  36. // Description: DeSerializing constructor (locked)
  37. //--------------------------------------------------------------------------------------
  38. KernelV0::KernelV0 (const YAML::Node& node, const ctor_key&) : LemmaObject(node) {
  39. } // ----- end of method KernelV0::KernelV0 (constructor) -----
  40. //--------------------------------------------------------------------------------------
  41. // Class: KernelV0
  42. // Method: NewSP()
  43. // Description: public constructor returing a shared_ptr
  44. //--------------------------------------------------------------------------------------
  45. std::shared_ptr< KernelV0 > KernelV0::NewSP() {
  46. return std::make_shared< KernelV0 >( ctor_key() );
  47. }
  48. //--------------------------------------------------------------------------------------
  49. // Class: KernelV0
  50. // Method: ~KernelV0
  51. // Description: destructor (protected)
  52. //--------------------------------------------------------------------------------------
  53. KernelV0::~KernelV0 () {
  54. } // ----- end of method KernelV0::~KernelV0 (destructor) -----
  55. //--------------------------------------------------------------------------------------
  56. // Class: KernelV0
  57. // Method: Serialize
  58. //--------------------------------------------------------------------------------------
  59. YAML::Node KernelV0::Serialize ( ) const {
  60. YAML::Node node = LemmaObject::Serialize();
  61. node.SetTag( GetName() );
  62. // Coils Transmitters & Receivers
  63. for ( auto txm : TxRx) {
  64. node[txm.first] = txm.second->Serialize();
  65. }
  66. // LayeredEarthEM
  67. node["SigmaModel"] = SigmaModel->Serialize();
  68. node["Larmor"] = Larmor;
  69. node["Temperature"] = Temperature;
  70. node["tol"] = tol;
  71. node["minLevel"] = minLevel;
  72. node["maxLevel"] = maxLevel;
  73. node["Taup"] = Taup;
  74. node["PulseI"] = PulseI;
  75. node["Interfaces"] = Interfaces;
  76. for ( int ilay=0; ilay<Interfaces.size()-1; ++ilay ) {
  77. node["Kern-" + to_string(ilay) ] = static_cast<VectorXcr>(Kern.row(ilay));
  78. }
  79. return node;
  80. } // ----- end of method KernelV0::Serialize -----
  81. //--------------------------------------------------------------------------------------
  82. // Class: KernelV0
  83. // Method: DeSerialize
  84. //--------------------------------------------------------------------------------------
  85. std::shared_ptr<KernelV0> KernelV0::DeSerialize ( const YAML::Node& node ) {
  86. if (node.Tag() != "KernelV0" ) {
  87. throw DeSerializeTypeMismatch( "KernelV0", node.Tag());
  88. }
  89. return std::make_shared< KernelV0 > ( node, ctor_key() );
  90. } // ----- end of method KernelV0::DeSerialize -----
  91. //--------------------------------------------------------------------------------------
  92. // Class: KernelV0
  93. // Method: AlignWithAkvoDataset
  94. //--------------------------------------------------------------------------------------
  95. void KernelV0::AlignWithAkvoDataset( const YAML::Node& node ) {
  96. if (node["processed"].as<std::string>().substr(0,4) == "Akvo") {
  97. std::cout << "Akvo file read\n";
  98. std::cout << node["processed"] << std::endl;
  99. }
  100. if (node["pulseType"].as<std::string>() == "FID") {
  101. PulseI = node["Pulse 1"]["current"].as<VectorXr>();
  102. this->SetPulseDuration( node["pulseLength"].as<double>() );
  103. } else {
  104. std::cerr << "Pulse Type " << node["PulseType"] << "is not supported\n";
  105. }
  106. }
  107. //--------------------------------------------------------------------------------------
  108. // Class: KernelV0
  109. // Method: DeSerialize
  110. //--------------------------------------------------------------------------------------
  111. void KernelV0::CalculateK0 (const std::vector< std::string>& Tx, const std::vector<std::string >& Rx,
  112. bool vtkOutput ) {
  113. // Set up
  114. Larmor = SigmaModel->GetMagneticFieldMagnitude()*GAMMA; // in rad 2246.*2.*PI;
  115. // All EM calculations will share same field points
  116. cpoints = FieldPoints::NewSP();
  117. cpoints->SetNumberOfPoints(8);
  118. for (auto tx : Tx) {
  119. // Set up EMEarth
  120. EMEarths[tx] = EMEarth1D::NewSP();
  121. EMEarths[tx]->AttachWireAntenna(TxRx[tx]);
  122. EMEarths[tx]->AttachLayeredEarthEM(SigmaModel);
  123. EMEarths[tx]->AttachFieldPoints( cpoints );
  124. EMEarths[tx]->SetFieldsToCalculate(H);
  125. // TODO query for method, altough with flat antennae, this is fastest
  126. //EMEarths[tx]->SetHankelTransformMethod(FHTKEY201);
  127. EMEarths[tx]->SetHankelTransformMethod(ANDERSON801);
  128. EMEarths[tx]->SetTxRxMode(TX);
  129. TxRx[tx]->SetCurrent(1.);
  130. }
  131. for (auto rx : Rx) {
  132. if (EMEarths.count(rx)) {
  133. EMEarths[rx]->SetTxRxMode(TXRX);
  134. } else {
  135. EMEarths[rx] = EMEarth1D::NewSP();
  136. EMEarths[rx]->AttachWireAntenna(TxRx[rx]);
  137. EMEarths[rx]->AttachLayeredEarthEM(SigmaModel);
  138. EMEarths[rx]->AttachFieldPoints( cpoints );
  139. EMEarths[rx]->SetFieldsToCalculate(H);
  140. // TODO query for method, altough with flat antennae, this is fastest
  141. //EMEarths[rx]->SetHankelTransformMethod(FHTKEY201);
  142. EMEarths[rx]->SetHankelTransformMethod(ANDERSON801);
  143. EMEarths[rx]->SetTxRxMode(RX);
  144. TxRx[rx]->SetCurrent(1.);
  145. }
  146. }
  147. std::cout << "Calculating K0 kernel\n";
  148. Kern = MatrixXcr::Zero( Interfaces.size()-1, PulseI.size() );
  149. for (ilay=0; ilay<Interfaces.size()-1; ++ilay) {
  150. std::cout << "Layer " << ilay << "\tfrom " << Interfaces(ilay) <<" to "<< Interfaces(ilay+1) << std::endl;
  151. Size(2) = Interfaces(ilay+1) - Interfaces(ilay);
  152. Origin(2) = Interfaces(ilay);
  153. IntegrateOnOctreeGrid( vtkOutput );
  154. }
  155. std::cout << "\nFinished KERNEL\n";
  156. }
  157. //--------------------------------------------------------------------------------------
  158. // Class: KernelV0
  159. // Method: IntegrateOnOctreeGrid
  160. //--------------------------------------------------------------------------------------
  161. void KernelV0::IntegrateOnOctreeGrid( bool vtkOutput) {
  162. Vector3r cpos = Origin + Size/2.;
  163. VOLSUM = 0;
  164. nleaves = 0;
  165. if (!vtkOutput) {
  166. EvaluateKids( Size, 0, cpos, VectorXcr::Ones(PulseI.size()) );
  167. } else {
  168. #ifdef LEMMAUSEVTK
  169. vtkHyperOctree* oct = vtkHyperOctree::New();
  170. oct->SetDimension(3);
  171. oct->SetOrigin( Origin(0), Origin(1), Origin(2) );
  172. oct->SetSize( Size(0), Size(1), Size(2) );
  173. vtkHyperOctreeCursor* curse = oct->NewCellCursor();
  174. curse->ToRoot();
  175. EvaluateKids2( Size, 0, cpos, VectorXcr::Ones(PulseI.size()), oct, curse );
  176. for (int iq=0; iq<PulseI.size(); ++iq) {
  177. // Fill in leaf data
  178. vtkDoubleArray* kr = vtkDoubleArray::New();
  179. kr->SetNumberOfComponents(1);
  180. kr->SetName("Re($K_0$)");
  181. kr->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  182. vtkDoubleArray* ki = vtkDoubleArray::New();
  183. ki->SetNumberOfComponents(1);
  184. ki->SetName("Im($K_0$)");
  185. ki->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  186. vtkDoubleArray* km = vtkDoubleArray::New();
  187. km->SetNumberOfComponents(1);
  188. km->SetName("mod($K_0$)");
  189. km->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  190. vtkIntArray* kid = vtkIntArray::New();
  191. kid->SetNumberOfComponents(1);
  192. kid->SetName("ID");
  193. kid->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  194. vtkIntArray* kerr = vtkIntArray::New();
  195. kerr->SetNumberOfComponents(1);
  196. kerr->SetName("nleaf");
  197. //Real LeafVol(0);
  198. for (auto leaf : LeafDict) {
  199. kr->InsertTuple1( leaf.first, std::real(leaf.second(iq)) );
  200. ki->InsertTuple1( leaf.first, std::imag(leaf.second(iq)) );
  201. km->InsertTuple1( leaf.first, std::abs(leaf.second(iq)) );
  202. kid->InsertTuple1( leaf.first, leaf.first );
  203. //LeafVol += std::real(leaf.second);
  204. }
  205. //std::cout << "\n\nLeafVol=" << LeafVol << std::endl;
  206. for (auto leaf : LeafDictIdx) {
  207. kerr->InsertTuple1( leaf.first, leaf.second );
  208. }
  209. auto kri = oct->GetLeafData()->AddArray(kr);
  210. auto kii = oct->GetLeafData()->AddArray(ki);
  211. auto kmi = oct->GetLeafData()->AddArray(km);
  212. auto kidi = oct->GetLeafData()->AddArray(kid);
  213. auto keri = oct->GetLeafData()->AddArray(kerr);
  214. auto write = vtkXMLHyperOctreeWriter::New();
  215. //write.SetDataModeToAscii()
  216. write->SetInputData(oct);
  217. std::string fname = std::string("octree-") + to_string(ilay)
  218. + std::string("-") + to_string(iq) + std::string(".vto");
  219. write->SetFileName(fname.c_str());
  220. write->Write();
  221. write->Delete();
  222. oct->GetLeafData()->RemoveArray( kri );
  223. oct->GetLeafData()->RemoveArray( kii );
  224. oct->GetLeafData()->RemoveArray( kmi );
  225. oct->GetLeafData()->RemoveArray( kidi );
  226. oct->GetLeafData()->RemoveArray( keri );
  227. kerr->Delete();
  228. kid->Delete();
  229. kr->Delete();
  230. ki->Delete();
  231. km->Delete();
  232. }
  233. curse->Delete();
  234. oct->Delete();
  235. #else
  236. throw std::runtime_error("IntegrateOnOctreeGrid with vtkOutput requires Lemma with VTK support");
  237. #endif
  238. }
  239. std::cout << "\nVOLSUM=" << VOLSUM << "\tActual=" << Size(0)*Size(1)*Size(2)
  240. << "\tDifference=" << VOLSUM - (Size(0)*Size(1)*Size(2)) << std::endl;
  241. }
  242. //--------------------------------------------------------------------------------------
  243. // Class: KernelV0
  244. // Method: f
  245. //--------------------------------------------------------------------------------------
  246. #if 1
  247. VectorXcr KernelV0::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr ) {
  248. // Compute the elliptic fields
  249. Vector3r B0hat = SigmaModel->GetMagneticFieldUnitVector();
  250. Vector3r B0 = SigmaModel->GetMagneticField();
  251. // Elliptic representation
  252. EllipticB EBT = EllipticFieldRep(MU0*Ht, B0hat);
  253. EllipticB EBR = EllipticFieldRep(MU0*Hr, B0hat);
  254. // Compute Mn0
  255. Vector3r Mn0 = ComputeMn0(1.0, B0);
  256. //std::cout << "Mn0\t" << Mn0.transpose() << std::endl;
  257. Real Mn0Abs = Mn0.norm();
  258. // Compute phase delay
  259. // TODO add transmiiter current phase and delay induced apparent time phase!
  260. Complex PhaseTerm = EBR.bhat.dot(EBT.bhat) + Complex(0, (B0hat.dot(EBR.bhat.cross(EBT.bhat))));
  261. Complex ejztr = std::exp(Complex(0, EBR.zeta + EBT.zeta));
  262. // Calcuate vector of all responses
  263. VectorXcr F = VectorXcr::Zero( PulseI.size() );
  264. for (int iq=0; iq<PulseI.size(); ++iq) {
  265. // Compute the tipping angle
  266. Real sintheta = std::sin(0.5*GAMMA*PulseI(iq)*Taup*(EBT.alpha-EBT.beta));
  267. F(iq) = -volume*Complex(0,Larmor)*Mn0Abs*(EBR.alpha+EBR.beta)*ejztr*sintheta*PhaseTerm;
  268. //TODO TEST FOR ASYMETRY
  269. //Real sintheta = std::sin(0.5*GAMMA*PulseI(iq)*Taup*(EBT.alpha-EBT.beta));
  270. //F(iq) = volume * Complex(EBT.Bperp.real().norm(), EBT.Bperp.imag().norm()); //Complex(sintheta, EBT.Bperp.norm() );
  271. //F(iq) = volume * Complex(EBT.alpha, EBT.beta);
  272. //F(iq) = volume * MU0*Hr.norm();
  273. //F(iq) = volume * EBT.err;
  274. //F(iq) = volume * sintheta;
  275. }
  276. return F;
  277. }
  278. #endif
  279. #if 0
  280. VectorXcr KernelV0::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr ) {
  281. VectorXcr F = VectorXcr::Ones( PulseI.size() );
  282. F.array() *= volume * Complex(Ht.norm(), Hr.norm()); //*Ht.dot(Hr);
  283. return F;
  284. }
  285. #endif
  286. // //--------------------------------------------------------------------------------------
  287. // // Class: KernelV0
  288. // // Method: ComputeV0Cell
  289. // //--------------------------------------------------------------------------------------
  290. // Complex KernelV0::ComputeV0Cell(const EllipticB& EBT, const EllipticB& EBR,
  291. // const Real& sintheta, const Real& phase, const Real& Mn0Abs,
  292. // const Real& vol) {
  293. // // earth response of receiver adjoint field
  294. // Vector3r B0hat = SigmaModel->GetMagneticFieldUnitVector();
  295. // Complex ejztr = std::exp(Complex(0, EBR.zeta + EBT.zeta));
  296. // Complex PhaseTerm = EBR.bhat.dot(EBT.bhat) + (B0hat.dot(EBR.bhat.cross(EBT.bhat) ));
  297. // return -vol*Complex(0,Larmor)*Mn0Abs*(EBR.alpha+EBR.beta)*ejztr*sintheta*PhaseTerm;
  298. // }
  299. //--------------------------------------------------------------------------------------
  300. // Class: KernelV0
  301. // Method: ComputeV0Cell
  302. //--------------------------------------------------------------------------------------
  303. Vector3r KernelV0::ComputeMn0(const Real& Porosity, const Vector3r& B0) {
  304. Real chi_n = NH2O*((GAMMA*GAMMA*HBAR*HBAR)/(4.*KB*Temperature));
  305. return chi_n*Porosity*B0;
  306. }
  307. //--------------------------------------------------------------------------------------
  308. // Class: KernelV0
  309. // Method: ComputeV0Cell
  310. //--------------------------------------------------------------------------------------
  311. EllipticB KernelV0::EllipticFieldRep (const Vector3cr& B, const Vector3r& B0hat) {
  312. // This all follows Weichman et al., 2000.
  313. // There are some numerical stability issues that arise when the two terms in the beta
  314. // formulation are nearly equivalent. The current formulation will result in a null-valued
  315. // beta, although this does not entirely recreate the true value of B perp.
  316. // Reformulating may be welcome
  317. EllipticB ElipB = EllipticB();
  318. Vector3cr Bperp = B - B0hat.dot(B)*B0hat; // Eigen is OK with this
  319. //Vector3r Bperpr = B.real() - B0hat.dot(B.real())*B0hat;
  320. //Vector3r Bperpi = B.imag() - B0hat.dot(B.imag())*B0hat;
  321. //Vector3cr Bperp = Bperpr + Complex(0,1)*Bperpi;
  322. //ElipB.BperpdotB = Bperp.dot(B0hat); // TODO remove
  323. Real BperpNorm = Bperp.norm();
  324. //Complex Bp2 = Bperp.transpose() * Bperp;
  325. Complex Bp2 = Bperp.conjugate().dot(Bperp);
  326. VectorXcr iB0 = Complex(0,1)*B0hat.cast<Complex>().array();
  327. ElipB.eizt = std::sqrt(Bp2 / std::abs(Bp2));
  328. ElipB.alpha = INVSQRT2*std::sqrt(BperpNorm*BperpNorm + std::abs(Bp2));
  329. ElipB.beta = std::copysign(1, std::real(iB0.dot( Bperp.cross(Bperp.conjugate())) )) *
  330. (INVSQRT2*std::sqrt(BperpNorm*BperpNorm - std::abs(Bp2)));
  331. ElipB.bhat = ((Real)1./ElipB.alpha)*(((Real)1./ElipB.eizt)*Bperp.array()).real().array();
  332. ElipB.bhatp = B0hat.cross(ElipB.bhat);
  333. ElipB.zeta = std::real(std::log(ElipB.eizt)/Complex(0,1));
  334. /* use as an error check decomposed field - computed actual */
  335. // Vector3cr Bperp2 = ElipB.eizt * (ElipB.alpha * ElipB.bhat
  336. // + (Complex(0,1) * ElipB.beta * ElipB.bhatp) );
  337. // ElipB.err = (Bperp-Bperp2).norm();
  338. // if (ElipB.err > .01*Bperp.norm() ) {
  339. // std::cout << "Elip error\n";
  340. // Real Beta2 = sgn( std::real(iB0.dot( Bperp.cross(Bperp.conjugate())) )) *
  341. // (INVSQRT2*std::sqrt(BperpNorm*BperpNorm - std::abs(Bp2)));
  342. // Vector3cr Bperp3 = ElipB.eizt * (ElipB.alpha * ElipB.bhat
  343. // + (Complex(0,1) * Beta2 * ElipB.bhatp) );
  344. // std::cout << "Beta term0\t" << (INVSQRT2*std::sqrt(BperpNorm*BperpNorm - std::abs(Bp2))) << std::endl;
  345. // std::cout << "Beta term1\t" << BperpNorm*BperpNorm << "\t" << std::abs(Bp2) << std::endl;
  346. // std::cout << "Beta \t" << ElipB.beta << std::endl;
  347. // std::cout << "Beta2 \t" << Beta2 << std::endl;
  348. // std::cout << "Bperp \t" << Bperp.transpose() << std::endl;
  349. // std::cout << "Bperp2\t" << Bperp2.transpose() << std::endl;
  350. // std::cout << "Bperp3\t" << Bperp3.transpose() << std::endl;
  351. // std::cout << "err \t" << ElipB.err << std::endl;
  352. // }
  353. return ElipB;
  354. }
  355. //--------------------------------------------------------------------------------------
  356. // Class: KernelV0
  357. // Method: EvaluateKids
  358. //--------------------------------------------------------------------------------------
  359. void KernelV0::EvaluateKids( const Vector3r& size, const int& level, const Vector3r& cpos,
  360. const VectorXcr& parentVal ) {
  361. std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
  362. //std::cout.flush();
  363. // Next level step, interested in one level below
  364. // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
  365. Vector3r step = size.array() / (Real)(1 << (level+1) );
  366. Real vol = (step(0)*step(1)*step(2)); // volume of each child
  367. Vector3r pos = cpos - step/2.;
  368. Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
  369. 0, 0, 0,
  370. step[0], 0, 0,
  371. 0, step[1], 0,
  372. step[0], step[1], 0,
  373. 0, 0, step[2],
  374. step[0], 0, step[2],
  375. 0, step[1], step[2],
  376. step[0], step[1], step[2] ).finished();
  377. MatrixXcr kvals(8, PulseI.size()); // individual kernel vals
  378. cpoints->ClearFields();
  379. for (int ichild=0; ichild<8; ++ichild) {
  380. Vector3r cp = pos; // Eigen complains about combining these
  381. cp += posadd.row(ichild);
  382. cpoints->SetLocation( ichild, cp );
  383. }
  384. Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
  385. Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
  386. for ( auto EMCalc : EMEarths ) {
  387. EMCalc.second->GetFieldPoints()->ClearFields();
  388. EMCalc.second->CalculateWireAntennaFields();
  389. switch (EMCalc.second->GetTxRxMode()) {
  390. case TX:
  391. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  392. break;
  393. case RX:
  394. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  395. break;
  396. case TXRX:
  397. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  398. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  399. break;
  400. default:
  401. break;
  402. }
  403. }
  404. for (int ichild=0; ichild<8; ++ichild) {
  405. Vector3r cp = pos; // Eigen complains about combining these
  406. cp += posadd.row(ichild);
  407. kvals.row(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
  408. }
  409. VectorXcr ksum = kvals.colwise().sum(); // Kernel sum
  410. // Evaluate whether or not furthur splitting is needed
  411. if ( (((ksum - parentVal).array().abs() > tol).any() && level<maxLevel) || level < minLevel ) {
  412. // Not a leaf dive further in
  413. for (int ichild=0; ichild<8; ++ichild) {
  414. Vector3r cp = pos; // Eigen complains about combining these
  415. cp += posadd.row(ichild);
  416. EvaluateKids( size, level+1, cp, kvals.row(ichild) );
  417. }
  418. return; // not leaf
  419. }
  420. // implicit else, is a leaf
  421. Kern.row(ilay) += ksum;
  422. VOLSUM += 8.*vol;
  423. nleaves += 8; // reflects the number of kernel evaluations
  424. return; // is leaf
  425. }
  426. #ifdef LEMMAUSEVTK
  427. //--------------------------------------------------------------------------------------
  428. // Class: KernelV0
  429. // Method: EvaluateKids2 -- same as Evaluate Kids, but include VTK octree generation
  430. //--------------------------------------------------------------------------------------
  431. void KernelV0::EvaluateKids2( const Vector3r& size, const int& level, const Vector3r& cpos,
  432. const VectorXcr& parentVal, vtkHyperOctree* oct, vtkHyperOctreeCursor* curse) {
  433. std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
  434. //std::cout.flush();
  435. // Next level step, interested in one level below
  436. // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
  437. Vector3r step = size.array() / (Real)(1 << (level+1) );
  438. Real vol = (step(0)*step(1)*step(2)); // volume of each child
  439. Vector3r pos = cpos - step/2.;
  440. Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
  441. 0, 0, 0,
  442. step[0], 0, 0,
  443. 0, step[1], 0,
  444. step[0], step[1], 0,
  445. 0, 0, step[2],
  446. step[0], 0, step[2],
  447. 0, step[1], step[2],
  448. step[0], step[1], step[2] ).finished();
  449. MatrixXcr kvals(8, PulseI.size()); // individual kernel vals
  450. cpoints->ClearFields();
  451. for (int ichild=0; ichild<8; ++ichild) {
  452. Vector3r cp = pos; // Eigen complains about combining these
  453. cp += posadd.row(ichild);
  454. cpoints->SetLocation( ichild, cp );
  455. }
  456. Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
  457. Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
  458. for ( auto EMCalc : EMEarths ) {
  459. //EMCalc->GetFieldPoints()->ClearFields();
  460. EMCalc.second->CalculateWireAntennaFields();
  461. switch (EMCalc.second->GetTxRxMode()) {
  462. case TX:
  463. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  464. break;
  465. case RX:
  466. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  467. break;
  468. case TXRX:
  469. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  470. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  471. break;
  472. default:
  473. break;
  474. }
  475. }
  476. for (int ichild=0; ichild<8; ++ichild) {
  477. Vector3r cp = pos; // Eigen complains about combining these
  478. cp += posadd.row(ichild);
  479. kvals.row(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
  480. }
  481. VectorXcr ksum = kvals.colwise().sum(); // Kernel sum
  482. // Evaluate whether or not furthur splitting is needed
  483. if ( (((ksum - parentVal).array().abs() > tol).any() && level<maxLevel) || level < minLevel ) {
  484. oct->SubdivideLeaf(curse);
  485. for (int ichild=0; ichild<8; ++ichild) {
  486. curse->ToChild(ichild);
  487. Vector3r cp = pos; // Eigen complains about combining these
  488. cp += posadd.row(ichild);
  489. /* Test for position via alternative means */
  490. /*
  491. Real p[3];
  492. GetPosition(curse, p);
  493. if ( (Vector3r(p) - cp).norm() > 1e-8 ) {
  494. std::cout << "ERROR @ nleaves" << nleaves << "\n" << cp[0] << "\t" << p[0] << "\t" << cp[1] << "\t" << p[1]
  495. << "\t" << cp[2] << "\t" << p[2] << "\t" << vol<< std::endl;
  496. throw std::runtime_error("doom");
  497. }
  498. */
  499. /* End of position test */
  500. EvaluateKids2( size, level+1, cp, kvals.row(ichild), oct, curse );
  501. curse->ToParent();
  502. }
  503. return; // not a leaf
  504. }
  505. /* just stuff with sum of the kids and don't subdivide */
  506. /*
  507. LeafDict[curse->GetLeafId()] = ksum/(8.*vol);
  508. LeafDictIdx[curse->GetLeafId()] = nleaves;
  509. */
  510. /* Alternatively, subdivide the VTK octree here and stuff the children to make better
  511. * visuals, but also 8x the storage...
  512. */
  513. oct->SubdivideLeaf(curse);
  514. for (int ichild=0; ichild<8; ++ichild) {
  515. curse->ToChild(ichild);
  516. LeafDict[curse->GetLeafId()] = ksum/(8.*vol);
  517. LeafHt[curse->GetLeafId()] = Ht.col(ichild);
  518. LeafHr[curse->GetLeafId()] = Hr.col(ichild);
  519. LeafDictIdx[curse->GetLeafId()] = nleaves;
  520. curse->ToParent();
  521. }
  522. Kern.row(ilay) += ksum;
  523. VOLSUM += 8*vol;
  524. nleaves += 8; // good reason to say 1 or 8 here...8 sounds better and reflects kernel evaluations
  525. return; // is a leaf
  526. }
  527. //--------------------------------------------------------------------------------------
  528. // Class: KernelV0
  529. // Method: GetPosition
  530. //--------------------------------------------------------------------------------------
  531. void KernelV0::GetPosition( vtkHyperOctreeCursor* Cursor, Real* p ) {
  532. Real ratio=1.0/(1<<(Cursor->GetCurrentLevel()));
  533. //step = ((Size).array() / std::pow(2.,Cursor->GetCurrentLevel()));
  534. p[0]=(Cursor->GetIndex(0)+.5)*ratio*this->Size[0]+this->Origin[0] ;//+ .5*step[0];
  535. p[1]=(Cursor->GetIndex(1)+.5)*ratio*this->Size[1]+this->Origin[1] ;//+ .5*step[1];
  536. p[2]=(Cursor->GetIndex(2)+.5)*ratio*this->Size[2]+this->Origin[2] ;//+ .5*step[2];
  537. }
  538. #endif
  539. } // ---- end of namespace Lemma ----
  540. /* vim: set tabstop=4 expandtab */
  541. /* vim: set filetype=cpp */