瀏覽代碼

Fixes for python wrapper with ABC inheritence. Also fixed EM1D to not clone antenna in tight loop

master
Trevor Irons 4 年之前
父節點
當前提交
a84f85654c

+ 5
- 5
Modules/FDEM1D/examples/inp/points.inp 查看文件

@@ -1,5 +1,5 @@
1
-20 20 20          // number of points (nx, ny, nz) in (northing, easting, depth) z pos down
2
--100 -100 1e-1    // offset in x, y, z (location of first calculation)
3
-5 4 3 2 1 .5 .5 .5 .5 .5 .5 .5 .5 .5 1 2 3 4 5 // x direction point spacing, in metres, (nx-1) 
4
-5 4 3 2 1 .5 .5 .5 .5 .5 .5 .5 .5 .5 1 2 3 4 5 // y direction point spacing, in metres, (ny-1) 
5
-.1 .1 .1 2 2 2 3 3 3 4 4 4 4 4 4 5 5 5 5       // z direction point spacing, in metres, (nz-1) 
1
+40 40 40          // number of points (nx, ny, nz) in (northing, easting, depth) z pos down
2
+32.5 32.5 1e-1    // offset in x, y, z (location of first calculation)
3
+5 5 5 5 5 5 5 5 5 5 5 4 3 2 1 .5 .5 .5 .5 .5 .5 .5 .5 .5 1 2 3 4 5 5 5 5 5 5 5 5 5 5 5 // x direction point spacing, in metres, (nx-1) 
4
+5 5 5 5 5 5 5 5 5 5 5 4 3 2 1 .5 .5 .5 .5 .5 .5 .5 .5 .5 1 2 3 4 5 5 5 5 5 5 5 5 5 5 5 // y direction point spacing, in metres, (nx-1) 
5
+.1 .1 .1 1 1 1 1 1 1 1 1 1 1 2 2 2 3 3 3 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5       // z direction point spacing, in metres, (nz-1) 

+ 2
- 0
Modules/FDEM1D/include/DipoleSource.h 查看文件

@@ -16,6 +16,7 @@
16 16
 #include <memory>
17 17
 #include "LemmaObject.h"
18 18
 #include "LayeredEarthEM.h"
19
+//#include "PolygonalWireAntenna.h"
19 20
 
20 21
 #ifdef LEMMAUSEVTK
21 22
 #include "vtkActor.h"
@@ -51,6 +52,7 @@ namespace Lemma {
51 52
         friend std::ostream &operator<<(std::ostream &stream, const DipoleSource &ob);
52 53
 
53 54
         friend class EMEarth1D;
55
+        friend class PolygonalWireAntenna;
54 56
 
55 57
         public:
56 58
 

+ 1
- 3
Modules/FDEM1D/python/pyFDEM1D.cpp 查看文件

@@ -139,8 +139,7 @@ PYBIND11_MODULE(FDEM1D, m) {
139 139
             "Sets all frequencies, argument is numpy array of frequencies")
140 140
         ;
141 141
 
142
-    py::class_<Lemma::LayeredEarthEM, std::shared_ptr<Lemma::LayeredEarthEM> >
143
-        LayeredEarthEM(m, "LayeredEarthEM");
142
+    py::class_<Lemma::LayeredEarthEM, std::shared_ptr<Lemma::LayeredEarthEM>, Lemma::EarthModel > LayeredEarthEM(m, "LayeredEarthEM");
144 143
 
145 144
         // lifecycle
146 145
         LayeredEarthEM.def(py::init(&Lemma::LayeredEarthEM::NewSP))
@@ -226,7 +225,6 @@ PYBIND11_MODULE(FDEM1D, m) {
226 225
         .def("SetLayerBreathPermitivity", &Lemma::LayeredEarthEM::SetLayerBreathPermitivity,
227 226
             "Sets the permitivity breath for Cole-COle model")
228 227
 
229
-
230 228
         // methods
231 229
         .def("EvaluateColeColeModel", &Lemma::LayeredEarthEM::EvaluateColeColeModel,
232 230
             "Calculates complex resistivity based on cole-cole parameters")

+ 6
- 6
Modules/FDEM1D/src/DipoleSource.cpp 查看文件

@@ -275,6 +275,7 @@ namespace Lemma {
275 275
         lays = Earth->GetLayerAtThisDepth(Location[2]);
276 276
         layr = Earth->GetLayerAtThisDepth(Receivers->GetLocation(irec)[2]);
277 277
 
278
+        // TODO, avoid smart pointer here maybe?
278 279
         KernelManager = KernelEM1DManager::NewSP();
279 280
 
280 281
             KernelManager->SetEarth(Earth);
@@ -955,8 +956,9 @@ namespace Lemma {
955 956
                             }
956 957
                             break;
957 958
                         case H:
958
-                            f(5) = Hankel->Zgauss(5, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[5]));
959
-                            f(6) = Hankel->Zgauss(6, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[6]));
959
+                            // TI - Comparisons with Amira show slight better agreement when neglecting these grounding terms
960
+                            f(5) = 0; //Hankel->Zgauss(5, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[5]));
961
+                            f(6) = 0; //Hankel->Zgauss(6, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[6]));
960 962
                             f(7) = Hankel->Zgauss(7, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[7]))*KernelManager->GetRAWKernel(ik[7])->GetZs()/KernelManager->GetRAWKernel(ik[7])->GetZm();
961 963
                             f(8) = Hankel->Zgauss(8, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[8]))*KernelManager->GetRAWKernel(ik[8])->GetZs()/KernelManager->GetRAWKernel(ik[8])->GetZm();
962 964
                             f(9) = Hankel->Zgauss(9, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[9]))*KernelManager->GetRAWKernel(ik[9])->GetZs()/KernelManager->GetRAWKernel(ik[9])->GetZm();
@@ -965,14 +967,12 @@ namespace Lemma {
965 967
                                         Pol[1]*QPI*(sps*f(5)+c2p*f(6)/rho-cps*f(7)+c2p*f(8)/rho)*Moment,
966 968
                                         Pol[1]*QPI*scp*(-f(5)+(Real)(2.)*f(6)/rho-f(7)+(Real)(2.)*f(8)/rho)*Moment,
967 969
                                         -Pol[1]*QPI*cp*f(9)*Moment );
968
-                                // Analytic whole space solution could go here
969 970
                             }
970 971
                             if (std::abs(Pol[0]) > 0) {
971 972
                                 this->Receivers->AppendHfield(ifreq, irec,
972 973
                                     Pol[0]*Moment*QPI*scp*(f(5)-(Real)(2.)*f(6)/rho+f(7)-(Real)(2.)*f(8)/rho),
973 974
                                     Pol[0]*Moment*QPI*(-cps*f(5)+c2p*f(6)/rho+sps*f(7)+c2p*f(8)/rho),
974 975
                                     Pol[0]*Moment*QPI*sp*f(9) );
975
-                                // Analytic whole space solution
976 976
                             }
977 977
                             break;
978 978
 
@@ -982,8 +982,8 @@ namespace Lemma {
982 982
                             f(4) = 0;
983 983
                             f(2) = Hankel->Zgauss(2, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[2])) * KernelManager->GetRAWKernel(0)->GetZs();
984 984
                             f(3) = Hankel->Zgauss(3, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[3])) * KernelManager->GetRAWKernel(1)->GetZs();
985
-                            f(5) = Hankel->Zgauss(5, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[5]));
986
-                            f(6) = Hankel->Zgauss(6, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[6]));
985
+                            f(5) = 0;//Hankel->Zgauss(5, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[5]));
986
+                            f(6) = 0;//Hankel->Zgauss(6, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[6]));
987 987
                             f(7) = Hankel->Zgauss(7, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[7]))*KernelManager->GetRAWKernel(ik[7])->GetZs()/KernelManager->GetRAWKernel(ik[7])->GetZm();
988 988
                             f(8) = Hankel->Zgauss(8, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[8]))*KernelManager->GetRAWKernel(ik[8])->GetZs()/KernelManager->GetRAWKernel(ik[8])->GetZm();
989 989
                             f(9) = Hankel->Zgauss(9, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[9]))*KernelManager->GetRAWKernel(ik[9])->GetZs()/KernelManager->GetRAWKernel(ik[9])->GetZm();

+ 3
- 8
Modules/FDEM1D/src/EMEarth1D.cpp 查看文件

@@ -231,7 +231,7 @@ namespace Lemma {
231 231
         if (Antenna->GetName() == std::string("PolygonalWireAntenna") || Antenna->GetName() == std::string("TEMTransmitter") ) {
232 232
             icalc += 1;
233 233
             // Check to see if they are all on a plane? If so we can do this fast
234
-            if (Antenna->IsHorizontallyPlanar() && ( HankelType == ANDERSON801 || HankelType== FHTKEY201 || HankelType==FHTKEY101 ||
234
+            if ( Antenna->IsHorizontallyPlanar() && ( HankelType == ANDERSON801 || HankelType== FHTKEY201 || HankelType==FHTKEY101 ||
235 235
                                                      HankelType == FHTKEY51 || HankelType == FHTKONG61 || HankelType == FHTKONG121 ||
236 236
                                                      HankelType == FHTKONG241 || HankelType == IRONS )) {
237 237
                 #ifdef HAVE_BOOST_PROGRESS
@@ -246,11 +246,11 @@ namespace Lemma {
246 246
                     {
247 247
                     #endif
248 248
                     auto Hankel = HankelTransformFactory::NewSP( HankelType );
249
+                    auto AntCopy = static_cast<PolygonalWireAntenna*>(Antenna.get())->ClonePA();
249 250
                     #ifdef LEMMAUSEOMP
250 251
                     #pragma omp for schedule(static, 1)
251 252
                     #endif
252 253
                     for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
253
-                        auto AntCopy = static_cast<PolygonalWireAntenna*>(Antenna.get())->ClonePA();
254 254
                         SolveLaggedTxRxPair(irec, Hankel.get(), wavef, ifreq, AntCopy.get());
255 255
                         #ifdef HAVE_BOOST_PROGRESS
256 256
                         if (progressbar) ++(*disp);
@@ -270,7 +270,6 @@ namespace Lemma {
270 270
                     disp = new boost::progress_display( Receivers->GetNumberOfPoints()*Antenna->GetNumberOfFrequencies() );
271 271
                 }
272 272
                 #endif
273
-
274 273
                 // parallelise across receivers
275 274
                 #ifdef LEMMAUSEOMP
276 275
                 #pragma omp parallel
@@ -452,8 +451,7 @@ namespace Lemma {
452 451
                                     //}
453 452
                                     auto tDipole = Antenna->GetDipoleSource(idip);
454 453
                                     // Propogation constant in free space
455
-                                    Real wavef   = tDipole->GetAngularFrequency(ifreq) *
456
-                                          std::sqrt(MU0*EPSILON0);
454
+                                    Real wavef   = tDipole->GetAngularFrequency(ifreq) * std::sqrt(MU0*EPSILON0);
457 455
                                     SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
458 456
                                 } // dipole loop
459 457
                             } // frequency loop
@@ -778,10 +776,7 @@ namespace Lemma {
778 776
     void EMEarth1D::SolveLaggedTxRxPair(const int &irec, HankelTransform* Hankel,
779 777
                     const Real &wavef, const int &ifreq, PolygonalWireAntenna* antenna) {
780 778
 
781
-        //std::cout << "SolveLaggedTxRxPair" << std::endl;
782
-
783 779
         antenna->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
784
-
785 780
         // Determine the min and max arguments
786 781
         Real rhomin = 1e9;
787 782
         Real rhomax = 1e-9;

+ 4
- 2
Modules/FDEM1D/testing/CMakeLists.txt 查看文件

@@ -6,8 +6,10 @@ target_link_libraries(unittest_FEM1D_GetNameCheck "lemmacore" "fdem1d" "yaml-cpp
6 6
 CXXTEST_ADD_TEST(unittest_FEM1D_SerializeCheck SerializeCheck.cc ${CMAKE_CURRENT_SOURCE_DIR}/SerializeCheck.h)
7 7
 target_link_libraries(unittest_FEM1D_SerializeCheck "lemmacore" "fdem1d" "yaml-cpp")
8 8
 
9
-#CXXTEST_ADD_TEST(benchKiHa BenchKiHa.cc ${CMAKE_CURRENT_SOURCE_DIR}/BenchKiHa.h)
10
-#target_link_libraries(benchKiHa "lemmacore" "fdem1d" "yaml-cpp")
9
+if(KIHA_EM1D)
10
+	CXXTEST_ADD_TEST(benchKiHa BenchKiHa.cc ${CMAKE_CURRENT_SOURCE_DIR}/BenchKiHa.h)
11
+	target_link_libraries(benchKiHa "lemmacore" "fdem1d" "yaml-cpp")
12
+endif()
11 13
 
12 14
 #set_target_properties(unittest_FEM1D_GetNameCheck
13 15
 #	              unittest_FEM1D_SerializeCheck 

+ 4
- 2
Modules/LemmaCore/include/EarthModel.h 查看文件

@@ -33,6 +33,10 @@ namespace Lemma {
33 33
         public:
34 34
 
35 35
             // ====================  LIFECYCLE     =======================
36
+
37
+            /// Default destructor
38
+            ~EarthModel ();
39
+
36 40
             /** YAML Serializing method
37 41
              */
38 42
             YAML::Node Serialize() const;
@@ -112,8 +116,6 @@ namespace Lemma {
112 116
             /** Deserialize constructor */
113 117
             EarthModel (const YAML::Node& node, const ctor_key&);
114 118
 
115
-            /// Default protected constructor.
116
-            ~EarthModel ();
117 119
 
118 120
             // ====================  DATA MEMBERS  =========================
119 121
 

+ 13
- 1
Modules/LemmaCore/python/pyLemmaCore.cpp 查看文件

@@ -226,7 +226,19 @@ PYBIND11_MODULE(LemmaCore, m) {
226 226
     ;
227 227
 
228 228
     // ABC
229
-    //py::class_<Lemma::EarthModel, std::shared_ptr<Lemma::EarthModel> >(m, "EarthModel")
229
+    py::class_<Lemma::EarthModel, std::shared_ptr<Lemma::EarthModel> > abcEarthModel (m, "EarthModel");
230
+
231
+    // Modifiers
232
+    abcEarthModel.def("SetMagneticFieldIncDecMag", &Lemma::EarthModel::SetMagneticFieldIncDecMag, "Sets the inducing field" );
233
+    abcEarthModel.def("SetMagneticFieldComponents", &Lemma::EarthModel::SetMagneticFieldComponents, "Sets the inducing field" );
234
+
235
+    // Accessors
236
+    abcEarthModel.def("GetMagneticField", &Lemma::EarthModel::GetMagneticField, "returns the inducing field in Tesla" );
237
+    abcEarthModel.def("GetMagneticFieldInGauss", &Lemma::EarthModel::GetMagneticFieldInGauss, "returns the inducing field in Gauss" );
238
+    abcEarthModel.def("GetMagneticFieldUnitVector", &Lemma::EarthModel::GetMagneticFieldUnitVector, "returns the inducing field unit vector" );
239
+    abcEarthModel.def("GetMagneticFieldMagnitude", &Lemma::EarthModel::GetMagneticFieldMagnitude, "returns the magnitude of the inducing field" );
240
+    abcEarthModel.def("GetMagneticFieldMagnitudeInGauss", &Lemma::EarthModel::GetMagneticFieldMagnitudeInGauss, "returns the magnitude of the inducing field" );
241
+
230 242
 
231 243
     /*
232 244
     py::class_<Lemma::LayeredEarth, std::shared_ptr<Lemma::LayeredEarth> >(m, "LayeredEarth")

+ 1
- 1
python/setup.py 查看文件

@@ -44,7 +44,7 @@ setup(
44 44
     #'pyLemma.pyFDEM1D': ['pyFDEM1D.*.so']
45 45
     #'pyLemma.pyFDEM1D': ['pyFDEM1D.*.so']
46 46
     #'': ['LemmaCore.*.so','FDEM1D.*.so']
47
-    '': ['LemmaCore.*','FDEM1D.*']
47
+    '': ['LemmaCore.*','FDEM1D.*','Merlin.*']
48 48
   },
49 49
   zip_safe=False,
50 50
   cmdclass={'install':InstallPlatlib},

+ 1
- 1
vim/c.vim 查看文件

@@ -17,7 +17,7 @@ syn keyword nspace Lemma YAML Eigen
17 17
 
18 18
 " Lemma types
19 19
 highlight leType ctermfg=Yellow guifg=Yellow
20
-syn keyword leType HankelTransform FHTAnderson801 FHTKey201 FHTKey101 FHTKey51 GQChave QWEKey KernelEM1D KernelEM1DSpec KernelEM1DBase KernelEM1DManager DipoleSource EarthModel LayeredEarth LayeredEarthEM TEMSurvey TEMSurveyLine TEMSurveyLineRecord TEMInductiveReceiver WireAntenna PolygonalWireAntenna TEMTransmitter TEMReceiver Instrument InstrumentTem LemmaObject FieldPoints DCIPElectrode TEMSurveyData TEMSurveyLineData TEMSurveyLineRecordData  ASCIIParser CubicSplineInterpolator RectilinearGrid GridReader RectilinearGridReader RectilinearGridVTKExporter Filter WindowFilter DEMParticle DEMGrain Data DataReader EMEarth1D 
20
+syn keyword leType HankelTransform FHTAnderson801 FHTKey201 FHTKey101 FHTKey51 GQChave QWEKey KernelEM1D KernelEM1DSpec KernelEM1DBase KernelEM1DManager DipoleSource EarthModel LayeredEarth LayeredEarthEM TEMSurvey TEMSurveyLine TEMSurveyLineRecord TEMInductiveReceiver WireAntenna PolygonalWireAntenna TEMTransmitter TEMReceiver Instrument InstrumentTem LemmaObject FieldPoints DCIPElectrode TEMSurveyData TEMSurveyLineData TEMSurveyLineRecordData  ASCIIParser CubicSplineInterpolator RectilinearGrid GridReader RectilinearGridReader RectilinearGridVTKExporter Filter WindowFilter DEMParticle DEMGrain Data DataReader EMEarth1D KernelV0 
21 21
 
22 22
 " Deprecated Lemma Types
23 23
 highlight dleType ctermfg=Blue guifg=Blue 

Loading…
取消
儲存