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