123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199 |
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- #include <time.h>
- #include "Lemma"
- #include "banner.h"
-
- using namespace Lemma;
-
- #ifdef LEMMAUSEVTK
- #include "matplot.h"
- using namespace matplot;
- #endif
-
- int main() {
-
- clock_t launch = clock();
-
-
- std::string name;
- std::string version;
- std::string usage;
- name = "FEM Forward Modeller - 1D ";
- version = "1.0alpha";
- usage = "utfemforward [inputfile]";
- banner(name,version,usage);
-
- time_t curr = time(0);
- std::cout << std::endl << " Start time: " << ctime(&curr) << std::endl;
-
-
- DipoleSource* Trans = DipoleSource::New();
- ReceiverPoints* Receivers = ReceiverPoints::New();
- LayeredEarthEM* Earth = LayeredEarthEM::New();
-
-
-
-
-
-
- DataFEM* modelledData = DataFEM::New();
-
-
- InstrumentFem* theinstrument = InstrumentFem::New();
-
-
- Earth->SetNumberOfLayers(5);
- Earth->SetLayerConductivity((VectorXcr(5) << 0.,1.e-4,1.e-2,
- 1.e-4,1.e-6).finished());
- Earth->SetLayerThickness((VectorXr(3) << 20.,5.,50.).finished());
-
-
- Real momtemp;
- momtemp = 1;
- Real freq;
- freq = 391;
- Trans->SetType(MAGNETICDIPOLE);
- Trans->SetLocation(0,0,-1.e-3);
- Trans->SetPolarisation(ZPOLARISATION);
-
- Trans->SetNumberOfFrequencies(1);
- Trans->SetFrequency(0,freq);
- Trans->SetMoment(momtemp);
-
-
- Vector3r loc1;
- Vector3r loc2;
- Receivers->SetNumberOfReceivers(1);
- loc1 << 9,0,-1.e-3;
- Receivers->SetLocation(0,loc1);
-
-
-
- theinstrument->SetReceiverPoints(Receivers);
- theinstrument->SetDipoleSource(Trans);
- theinstrument->EMEarthModel(Earth);
- theinstrument->SetOutputData(modelledData);
-
-
- theinstrument->MakeCalculation();
-
-
- std::cout << Receivers->GetHfield(0,0)(2) << std::endl;
-
-
-
-
- theinstrument->Delete();
- modelledData->Delete();
-
-
-
-
- InstrumentFem* theinstrument2 = InstrumentFem::New();
- DataFEM* modelledData2 = DataFEM::New();
-
-
-
- std::string datfile;
- datfile = "126.obs";
-
- DataReaderFemUBC* Reader = DataReaderFemUBC::New();
- DataFEM* inpdata = DataFEM::New();
- Reader->SetDataFEM(inpdata);
- try {
- Reader->ReadData(datfile,1);
- } catch(std::exception& e) {
- exit(EXIT_FAILURE);
- }
- std::cout << *inpdata << std::endl;
- theinstrument2->SetOutputData(modelledData2);
- theinstrument2->AlignWithData(inpdata);
-
-
- inpdata->Delete();
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- theinstrument2->Delete();
- Reader->Delete();
- Earth->Delete();
- Receivers->Delete();
- Trans->Delete();
-
-
- Real time;
- clock_t done = clock();
- time = (done-launch)/CLOCKS_PER_SEC;
- std::cout << " Execution time: " << time << " [CPU] seconds."
- << std::endl;
-
- return EXIT_SUCCESS;
- }
|