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 24KB

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