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