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