Main Lemma Repository
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

femforward.cpp 5.0KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199
  1. // ===========================================================================
  2. //
  3. // Filename: utfemforward.cpp
  4. //
  5. // Description:
  6. //
  7. // Version: 0.0
  8. // Created: 01/15/2013 12:11:34 PM
  9. // Revision: none
  10. // Compiler: Tested with g++
  11. //
  12. // Author: M. Andy Kass (MAK)
  13. //
  14. // Organisation: Broken Spoke Development, LLC
  15. //
  16. //
  17. // Email: mkass@numericalgeo.com
  18. //
  19. // This program is free software: you can redistribute it and/or modify
  20. // it under the terms of the GNU General Public License as published by
  21. // the Free Software Foundation, either version 3 of the License, or
  22. // (at your option) any later version.
  23. //
  24. // This program is distributed in the hope that it will be useful,
  25. // but WITHOUT ANY WARRANTY; without even the implied warranty of
  26. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  27. // GNU General Public License for more details.
  28. //
  29. // You should have received a copy of the GNU General Public License
  30. // along with this program. If not, see <http://www.gnu.org/licenses/>.
  31. //
  32. // ===========================================================================
  33. #include <time.h>
  34. #include "Lemma"
  35. #include "banner.h"
  36. using namespace Lemma;
  37. #ifdef LEMMAUSEVTK
  38. #include "matplot.h"
  39. using namespace matplot;
  40. #endif
  41. int main() {
  42. clock_t launch = clock();
  43. // Banner display
  44. std::string name;
  45. std::string version;
  46. std::string usage;
  47. name = "FEM Forward Modeller - 1D ";
  48. version = "1.0alpha";
  49. usage = "utfemforward [inputfile]";
  50. banner(name,version,usage);
  51. time_t curr = time(0);
  52. std::cout << std::endl << " Start time: " << ctime(&curr) << std::endl;
  53. // Define some objects
  54. DipoleSource* Trans = DipoleSource::New();
  55. ReceiverPoints* Receivers = ReceiverPoints::New();
  56. LayeredEarthEM* Earth = LayeredEarthEM::New();
  57. // SINGLE SOUNDING
  58. // Create new data object to receive forward model
  59. DataFEM* modelledData = DataFEM::New();
  60. // Create instrument object
  61. InstrumentFem* theinstrument = InstrumentFem::New();
  62. // Set up EMEarthModel
  63. Earth->SetNumberOfLayers(5);
  64. Earth->SetLayerConductivity((VectorXcr(5) << 0.,1.e-4,1.e-2,
  65. 1.e-4,1.e-6).finished());
  66. Earth->SetLayerThickness((VectorXr(3) << 20.,5.,50.).finished());
  67. // Set up transmitter
  68. Real momtemp;
  69. momtemp = 1;
  70. Real freq;
  71. freq = 391;
  72. Trans->SetType(MAGNETICDIPOLE);
  73. Trans->SetLocation(0,0,-1.e-3);
  74. Trans->SetPolarisation(ZPOLARISATION);
  75. // as a test
  76. Trans->SetNumberOfFrequencies(1);
  77. Trans->SetFrequency(0,freq);
  78. Trans->SetMoment(momtemp);
  79. // Set up receivers
  80. Vector3r loc1;
  81. Vector3r loc2;
  82. Receivers->SetNumberOfReceivers(1);
  83. loc1 << 9,0,-1.e-3;
  84. Receivers->SetLocation(0,loc1);
  85. // Attach to instrument
  86. theinstrument->SetReceiverPoints(Receivers);
  87. theinstrument->SetDipoleSource(Trans);
  88. theinstrument->EMEarthModel(Earth);
  89. theinstrument->SetOutputData(modelledData);
  90. // Calculate
  91. theinstrument->MakeCalculation();
  92. // See what comes back...
  93. std::cout << Receivers->GetHfield(0,0)(2) << std::endl;
  94. // Need to convert H field to ppm
  95. theinstrument->Delete();
  96. modelledData->Delete();
  97. // MULTIPLE SOUNDING MODEL USING READ-IN DATA
  98. InstrumentFem* theinstrument2 = InstrumentFem::New();
  99. DataFEM* modelledData2 = DataFEM::New();
  100. // Using the UBC data reader to read in the survey params
  101. // this will be used for the multiple sounding testing
  102. std::string datfile;
  103. datfile = "126.obs";
  104. DataReaderFemUBC* Reader = DataReaderFemUBC::New();
  105. DataFEM* inpdata = DataFEM::New();
  106. Reader->SetDataFEM(inpdata);
  107. try {
  108. Reader->ReadData(datfile,1);
  109. } catch(std::exception& e) {
  110. exit(EXIT_FAILURE);
  111. }
  112. std::cout << *inpdata << std::endl;
  113. theinstrument2->SetOutputData(modelledData2);
  114. theinstrument2->AlignWithData(inpdata);
  115. // Now, how do I want to pass it an earth model for each sounding?
  116. //
  117. inpdata->Delete();
  118. // Need to loop over transmitters
  119. // The relative positions of the transmitter and receiver is fixed
  120. // For testing, do some kind of random realisation of Earth model?
  121. // int nObs;
  122. // int nFreqs;
  123. // nObs = inpdata->GetnObs();
  124. // nFreqs = inpdata->GetnFreq();
  125. // Trans->SetType(MAGNETICDIPOLE);
  126. // Trans->SetLocation(0,0,-1e-3);
  127. // Vector3r loc;
  128. // Receivers->SetNumberOfReceivers(1);
  129. // loc << 0,0,-1e-3;
  130. // Receivers->SetLocation(0,loc);
  131. // Earth->SetNumberOfLayers(5); //Including the two halfspaces
  132. // Earth->SetLayerConductivity((VectorXcr(5) << 0.,1.e-4,1.e-2,
  133. // 1.e-4,1.e-6).finished());
  134. // Earth->SetLayerThickness((VectorXr(3) << 20.,5.,50.).finished());
  135. // for (int ii=0;ii<nObs;ii++) {
  136. // for (int jj=0;jj<nFreqs;jj++) {
  137. // Trans->SetMoment(inpdata->GetTxMom().coeff(jj));
  138. //
  139. // }
  140. // }
  141. //This is a bunch of testing stuff
  142. //theinstrument->AlignWithData(inpdata);
  143. //theinstrument->MakeCalculation();
  144. //inpdata->Delete();
  145. //theinstrument->MakeCalculation();
  146. theinstrument2->Delete();
  147. Reader->Delete();
  148. Earth->Delete();
  149. Receivers->Delete();
  150. Trans->Delete();
  151. // Timey-wimey stuff
  152. Real time;
  153. clock_t done = clock();
  154. time = (done-launch)/CLOCKS_PER_SEC;
  155. std::cout << " Execution time: " << time << " [CPU] seconds."
  156. << std::endl;
  157. return EXIT_SUCCESS;
  158. }