Surface NMR forward modelling
Você não pode selecionar mais de 25 tópicos Os tópicos devem começar com uma letra ou um número, podem incluir traços ('-') e podem ter até 35 caracteres.

KernelV0.cpp 25KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562
  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(ANDERSON801);
  127. EMEarths[tx]->SetTxRxMode(TX);
  128. TxRx[tx]->SetCurrent(1.);
  129. }
  130. for (auto rx : Rx) {
  131. if (EMEarths.count(rx)) {
  132. EMEarths[rx]->SetTxRxMode(TXRX);
  133. } else {
  134. EMEarths[rx] = EMEarth1D::NewSP();
  135. EMEarths[rx]->AttachWireAntenna(TxRx[rx]);
  136. EMEarths[rx]->AttachLayeredEarthEM(SigmaModel);
  137. EMEarths[rx]->AttachFieldPoints( cpoints );
  138. EMEarths[rx]->SetFieldsToCalculate(H);
  139. // TODO query for method, altough with flat antennae, this is fastest
  140. EMEarths[rx]->SetHankelTransformMethod(ANDERSON801);
  141. EMEarths[rx]->SetTxRxMode(RX);
  142. TxRx[rx]->SetCurrent(1.);
  143. }
  144. }
  145. std::cout << "Calculating K0 kernel\n";
  146. Kern = MatrixXcr::Zero( Interfaces.size()-1, PulseI.size() );
  147. for (ilay=0; ilay<Interfaces.size()-1; ++ilay) {
  148. std::cout << "Layer " << ilay << "\tfrom " << Interfaces(ilay) <<" to "<< Interfaces(ilay+1) << std::endl;
  149. Size(2) = Interfaces(ilay+1) - Interfaces(ilay);
  150. Origin(2) = Interfaces(ilay);
  151. IntegrateOnOctreeGrid( vtkOutput );
  152. }
  153. std::cout << "\nFinished KERNEL\n";
  154. }
  155. //--------------------------------------------------------------------------------------
  156. // Class: KernelV0
  157. // Method: IntegrateOnOctreeGrid
  158. //--------------------------------------------------------------------------------------
  159. void KernelV0::IntegrateOnOctreeGrid( bool vtkOutput) {
  160. Vector3r cpos = Origin + Size/2.;
  161. VOLSUM = 0;
  162. nleaves = 0;
  163. if (!vtkOutput) {
  164. EvaluateKids( Size, 0, cpos, VectorXcr::Ones(PulseI.size()) );
  165. } else {
  166. #ifdef LEMMAUSEVTK
  167. vtkHyperOctree* oct = vtkHyperOctree::New();
  168. oct->SetDimension(3);
  169. oct->SetOrigin( Origin(0), Origin(1), Origin(2) );
  170. oct->SetSize( Size(0), Size(1), Size(2) );
  171. vtkHyperOctreeCursor* curse = oct->NewCellCursor();
  172. curse->ToRoot();
  173. EvaluateKids2( Size, 0, cpos, VectorXcr::Ones(PulseI.size()), oct, curse );
  174. for (int iq=0; iq<PulseI.size(); ++iq) {
  175. // Fill in leaf data
  176. vtkDoubleArray* kr = vtkDoubleArray::New();
  177. kr->SetNumberOfComponents(1);
  178. kr->SetName("Re($K_0$)");
  179. kr->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  180. vtkDoubleArray* ki = vtkDoubleArray::New();
  181. ki->SetNumberOfComponents(1);
  182. ki->SetName("Im($K_0$)");
  183. ki->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  184. vtkDoubleArray* km = vtkDoubleArray::New();
  185. km->SetNumberOfComponents(1);
  186. km->SetName("mod($K_0$)");
  187. km->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  188. vtkIntArray* kid = vtkIntArray::New();
  189. kid->SetNumberOfComponents(1);
  190. kid->SetName("ID");
  191. kid->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  192. vtkIntArray* kerr = vtkIntArray::New();
  193. kerr->SetNumberOfComponents(1);
  194. kerr->SetName("nleaf");
  195. //Real LeafVol(0);
  196. for (auto leaf : LeafDict) {
  197. kr->InsertTuple1( leaf.first, std::real(leaf.second(iq)) );
  198. ki->InsertTuple1( leaf.first, std::imag(leaf.second(iq)) );
  199. km->InsertTuple1( leaf.first, std::abs(leaf.second(iq)) );
  200. kid->InsertTuple1( leaf.first, leaf.first );
  201. //LeafVol += std::real(leaf.second);
  202. }
  203. //std::cout << "\n\nLeafVol=" << LeafVol << std::endl;
  204. for (auto leaf : LeafDictIdx) {
  205. kerr->InsertTuple1( leaf.first, leaf.second );
  206. }
  207. auto kri = oct->GetLeafData()->AddArray(kr);
  208. auto kii = oct->GetLeafData()->AddArray(ki);
  209. auto kmi = oct->GetLeafData()->AddArray(km);
  210. auto kidi = oct->GetLeafData()->AddArray(kid);
  211. auto keri = oct->GetLeafData()->AddArray(kerr);
  212. auto write = vtkXMLHyperOctreeWriter::New();
  213. //write.SetDataModeToAscii()
  214. write->SetInputData(oct);
  215. std::string fname = std::string("octree-") + to_string(ilay)
  216. + std::string("-") + to_string(iq) + std::string(".vto");
  217. write->SetFileName(fname.c_str());
  218. write->Write();
  219. write->Delete();
  220. oct->GetLeafData()->RemoveArray( kri );
  221. oct->GetLeafData()->RemoveArray( kii );
  222. oct->GetLeafData()->RemoveArray( kmi );
  223. oct->GetLeafData()->RemoveArray( kidi );
  224. oct->GetLeafData()->RemoveArray( keri );
  225. kerr->Delete();
  226. kid->Delete();
  227. kr->Delete();
  228. ki->Delete();
  229. km->Delete();
  230. }
  231. curse->Delete();
  232. oct->Delete();
  233. #else
  234. throw std::runtime_error("IntegrateOnOctreeGrid with vtkOutput requires Lemma with VTK support");
  235. #endif
  236. }
  237. std::cout << "\nVOLSUM=" << VOLSUM << "\tActual=" << Size(0)*Size(1)*Size(2)
  238. << "\tDifference=" << VOLSUM - (Size(0)*Size(1)*Size(2)) << std::endl;
  239. }
  240. //--------------------------------------------------------------------------------------
  241. // Class: KernelV0
  242. // Method: f
  243. //--------------------------------------------------------------------------------------
  244. VectorXcr KernelV0::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr ) {
  245. // Compute the elliptic fields
  246. Vector3r B0hat = SigmaModel->GetMagneticFieldUnitVector();
  247. Vector3r B0 = SigmaModel->GetMagneticField();
  248. // Elliptic representation
  249. EllipticB EBT = EllipticFieldRep(MU0*Ht, B0hat);
  250. EllipticB EBR = EllipticFieldRep(MU0*Hr, B0hat);
  251. // Compute Mn0
  252. Vector3r Mn0 = ComputeMn0(1.0, B0);
  253. Real Mn0Abs = Mn0.norm();
  254. // Compute phase delay
  255. // TODO add transmiiter current phase and delay induced apparent time phase!
  256. Complex PhaseTerm = EBR.bhat.dot(EBT.bhat) + Complex(0, (B0hat.dot(EBR.bhat.cross(EBT.bhat))));
  257. Complex ejztr = std::exp(Complex(0, EBR.zeta + EBT.zeta));
  258. // Calcuate vector of all responses
  259. VectorXcr F = VectorXcr::Zero( PulseI.size() );
  260. for (int iq=0; iq<PulseI.size(); ++iq) {
  261. // Compute the tipping angle
  262. Real sintheta = std::sin(0.5*GAMMA*PulseI(iq)*Taup*(EBT.alpha-EBT.beta)); // why std::abs
  263. F(iq) = -volume*Complex(0,Larmor)*Mn0Abs*(EBR.alpha+EBR.beta)*ejztr*sintheta*PhaseTerm;
  264. }
  265. return F;
  266. }
  267. // //--------------------------------------------------------------------------------------
  268. // // Class: KernelV0
  269. // // Method: ComputeV0Cell
  270. // //--------------------------------------------------------------------------------------
  271. // Complex KernelV0::ComputeV0Cell(const EllipticB& EBT, const EllipticB& EBR,
  272. // const Real& sintheta, const Real& phase, const Real& Mn0Abs,
  273. // const Real& vol) {
  274. // // earth response of receiver adjoint field
  275. // Vector3r B0hat = SigmaModel->GetMagneticFieldUnitVector();
  276. // Complex ejztr = std::exp(Complex(0, EBR.zeta + EBT.zeta));
  277. // Complex PhaseTerm = EBR.bhat.dot(EBT.bhat) + (B0hat.dot(EBR.bhat.cross(EBT.bhat) ));
  278. // return -vol*Complex(0,Larmor)*Mn0Abs*(EBR.alpha+EBR.beta)*ejztr*sintheta*PhaseTerm;
  279. // }
  280. //--------------------------------------------------------------------------------------
  281. // Class: KernelV0
  282. // Method: ComputeV0Cell
  283. //--------------------------------------------------------------------------------------
  284. Vector3r KernelV0::ComputeMn0(const Real& Porosity, const Vector3r& B0) {
  285. Real chi_n = NH2O*((GAMMA*GAMMA*HBAR*HBAR)/(4.*KB*Temperature));
  286. return chi_n*Porosity*B0;
  287. }
  288. //--------------------------------------------------------------------------------------
  289. // Class: KernelV0
  290. // Method: ComputeV0Cell
  291. //--------------------------------------------------------------------------------------
  292. EllipticB KernelV0::EllipticFieldRep (const Vector3cr& B, const Vector3r& B0hat) {
  293. EllipticB ElipB = EllipticB();
  294. Vector3cr Bperp = B.array() - B0hat.dot(B)*B0hat.array();
  295. Real BperpNorm = Bperp.norm();
  296. Complex Bp2 = Bperp.transpose() * Bperp;
  297. VectorXcr iB0 = Complex(0,1)*B0hat.cast<Complex>().array();
  298. ElipB.eizt = std::sqrt(Bp2 / std::abs(Bp2));
  299. ElipB.alpha = INVSQRT2*std::sqrt(BperpNorm*BperpNorm + std::abs(Bp2));
  300. ElipB.beta = sgn(std::real(iB0.dot(Bperp.cross(Bperp.conjugate())))) *
  301. (INVSQRT2)*std::sqrt(std::abs(BperpNorm*BperpNorm-std::abs(Bp2)));
  302. ElipB.bhat = ((Real)1./ElipB.alpha)*(((Real)1./ElipB.eizt)*Bperp.array()).real().array();
  303. ElipB.bhatp = B0hat.cross(ElipB.bhat);
  304. ElipB.zeta = std::real(std::log(ElipB.eizt)/Complex(0,1));
  305. return ElipB;
  306. }
  307. //--------------------------------------------------------------------------------------
  308. // Class: KernelV0
  309. // Method: EvaluateKids
  310. //--------------------------------------------------------------------------------------
  311. void KernelV0::EvaluateKids( const Vector3r& size, const int& level, const Vector3r& cpos,
  312. const VectorXcr& parentVal ) {
  313. std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
  314. //std::cout.flush();
  315. // Next level step, interested in one level below
  316. // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
  317. Vector3r step = size.array() / (Real)(1 << (level+1) );
  318. Real vol = (step(0)*step(1)*step(2)); // volume of each child
  319. Vector3r pos = cpos - step/2.;
  320. Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
  321. 0, 0, 0,
  322. step[0], 0, 0,
  323. 0, step[1], 0,
  324. step[0], step[1], 0,
  325. 0, 0, step[2],
  326. step[0], 0, step[2],
  327. 0, step[1], step[2],
  328. step[0], step[1], step[2] ).finished();
  329. MatrixXcr kvals(8, PulseI.size()); // individual kernel vals
  330. cpoints->ClearFields();
  331. for (int ichild=0; ichild<8; ++ichild) {
  332. Vector3r cp = pos; // Eigen complains about combining these
  333. cp += posadd.row(ichild);
  334. cpoints->SetLocation( ichild, cp );
  335. }
  336. Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
  337. Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
  338. for ( auto EMCalc : EMEarths ) {
  339. EMCalc.second->GetFieldPoints()->ClearFields();
  340. EMCalc.second->CalculateWireAntennaFields();
  341. switch (EMCalc.second->GetTxRxMode()) {
  342. case TX:
  343. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  344. break;
  345. case RX:
  346. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  347. break;
  348. case TXRX:
  349. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  350. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  351. break;
  352. default:
  353. break;
  354. }
  355. }
  356. for (int ichild=0; ichild<8; ++ichild) {
  357. Vector3r cp = pos; // Eigen complains about combining these
  358. cp += posadd.row(ichild);
  359. kvals.row(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
  360. }
  361. VectorXcr ksum = kvals.colwise().sum(); // Kernel sum
  362. // Evaluate whether or not furthur splitting is needed
  363. if ( (((ksum - parentVal).array().abs() > tol).any() && level<maxLevel) || level < minLevel ) {
  364. // Not a leaf dive further in
  365. for (int ichild=0; ichild<8; ++ichild) {
  366. Vector3r cp = pos; // Eigen complains about combining these
  367. cp += posadd.row(ichild);
  368. EvaluateKids( size, level+1, cp, kvals.row(ichild) );
  369. }
  370. return; // not leaf
  371. }
  372. // implicit else, is a leaf
  373. Kern.row(ilay) += ksum;
  374. VOLSUM += 8.*vol;
  375. nleaves += 8; // reflects the number of kernel evaluations
  376. return; // is leaf
  377. }
  378. #ifdef LEMMAUSEVTK
  379. //--------------------------------------------------------------------------------------
  380. // Class: KernelV0
  381. // Method: EvaluateKids2 -- same as Evaluate Kids, but include VTK octree generation
  382. //--------------------------------------------------------------------------------------
  383. void KernelV0::EvaluateKids2( const Vector3r& size, const int& level, const Vector3r& cpos,
  384. const VectorXcr& parentVal, vtkHyperOctree* oct, vtkHyperOctreeCursor* curse) {
  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. MatrixXcr kvals(8, PulseI.size()); // 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->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.row(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
  432. }
  433. VectorXcr ksum = kvals.colwise().sum(); // Kernel sum
  434. // Evaluate whether or not furthur splitting is needed
  435. if ( (((ksum - parentVal).array().abs() > tol).any() && level<maxLevel) || level < minLevel ) {
  436. oct->SubdivideLeaf(curse);
  437. for (int ichild=0; ichild<8; ++ichild) {
  438. curse->ToChild(ichild);
  439. Vector3r cp = pos; // Eigen complains about combining these
  440. cp += posadd.row(ichild);
  441. /* Test for position via alternative means */
  442. /*
  443. Real p[3];
  444. GetPosition(curse, p);
  445. if ( (Vector3r(p) - cp).norm() > 1e-8 ) {
  446. std::cout << "ERROR @ nleaves" << nleaves << "\n" << cp[0] << "\t" << p[0] << "\t" << cp[1] << "\t" << p[1]
  447. << "\t" << cp[2] << "\t" << p[2] << "\t" << vol<< std::endl;
  448. throw std::runtime_error("doom");
  449. }
  450. */
  451. /* End of position test */
  452. EvaluateKids2( size, level+1, cp, kvals.row(ichild), oct, curse );
  453. curse->ToParent();
  454. }
  455. return; // not a leaf
  456. }
  457. /* just stuff with sum of the kids and don't subdivide */
  458. /*
  459. LeafDict[curse->GetLeafId()] = ksum/(8.*vol);
  460. LeafDictIdx[curse->GetLeafId()] = nleaves;
  461. */
  462. /* Alternatively, subdivide the VTK octree here and stuff the children to make better
  463. * visuals, but also 8x the storage...
  464. */
  465. oct->SubdivideLeaf(curse);
  466. for (int ichild=0; ichild<8; ++ichild) {
  467. curse->ToChild(ichild);
  468. LeafDict[curse->GetLeafId()] = ksum/(8.*vol);
  469. LeafDictIdx[curse->GetLeafId()] = nleaves;
  470. curse->ToParent();
  471. }
  472. Kern.row(ilay) += ksum;
  473. VOLSUM += 8*vol;
  474. nleaves += 8; // good reason to say 1 or 8 here...8 sounds better and reflects kernel evaluations
  475. return; // is a leaf
  476. }
  477. //--------------------------------------------------------------------------------------
  478. // Class: KernelV0
  479. // Method: GetPosition
  480. //--------------------------------------------------------------------------------------
  481. void KernelV0::GetPosition( vtkHyperOctreeCursor* Cursor, Real* p ) {
  482. Real ratio=1.0/(1<<(Cursor->GetCurrentLevel()));
  483. //step = ((Size).array() / std::pow(2.,Cursor->GetCurrentLevel()));
  484. p[0]=(Cursor->GetIndex(0)+.5)*ratio*this->Size[0]+this->Origin[0] ;//+ .5*step[0];
  485. p[1]=(Cursor->GetIndex(1)+.5)*ratio*this->Size[1]+this->Origin[1] ;//+ .5*step[1];
  486. p[2]=(Cursor->GetIndex(2)+.5)*ratio*this->Size[2]+this->Origin[2] ;//+ .5*step[2];
  487. }
  488. #endif
  489. } // ---- end of namespace Lemma ----
  490. /* vim: set tabstop=4 expandtab */
  491. /* vim: set filetype=cpp */