12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514 |
- // ===========================================================================
- //
- // Filename: EMSchur3DBase.cpp
- //
- // Created: 09/20/2013 04:53:40 PM
- // Compiler: Tested with g++, icpc, and MSVC 2010
- //
- // Author: Trevor Irons (ti)
- //
- // Organisation: University of Utah
- //
- // Email: tirons@egi.utah.edu
- //
- // ===========================================================================
-
- /**
- @file
- @author Trevor Irons
- @date 09/20/2013
- @version $Id$
- THis is a port of the procedural code developed at School of Mines
- **/
-
- #include "EMSchur3DBase.h"
-
- typedef Eigen::Triplet<Lemma::Complex> Tc;
- typedef Eigen::Triplet<Lemma::Real> Tr;
-
- #define UPPER 0
- #define LOWER 1 // 1=true, 0=false
-
- namespace Lemma {
-
-
- // ==================== FRIEND METHODS =====================
-
- std::ostream &operator<<(std::ostream &stream, const EMSchur3DBase &ob) {
- stream << ob.Serialize() << "\n";
- return stream;
- }
-
-
- // ==================== LIFECYCLE =======================
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: EMSchur3DBase
- // Description: constructor (protected)
- //--------------------------------------------------------------------------------------
- // std::shared_ptr<EMSchur3DBase> EMSchur3DBase::NewSP() {
- // return std::make_shared<EMSchur3DBase>( ctor_key() );
- // }
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: EMSchur3DBase
- // Description: constructor (protected)
- //--------------------------------------------------------------------------------------
- EMSchur3DBase::EMSchur3DBase ( const ctor_key& key ) : LemmaObject( key ),
- Grid(nullptr),
- //Survey(nullptr),
- LayModel(nullptr), Cvec(nullptr),
- ResFile("source"), sigma(nullptr), sigmap(nullptr)
- {
- } // ----- end of method EMSchur3DBase::EMSchur3DBase (constructor) -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: EMSchur3DBase
- // Description: constructor (protected)
- //--------------------------------------------------------------------------------------
- EMSchur3DBase::EMSchur3DBase ( const YAML::Node& node, const ctor_key& key ) : LemmaObject( key ),
- Grid(nullptr),
- //Survey(nullptr),
- LayModel(nullptr), Cvec(nullptr),
- ResFile("source"), sigma(nullptr), sigmap(nullptr)
- {
- } // ----- end of method EMSchur3DBase::~EMSchur3DBase (destructor) -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: ~EMSchur3DBase
- // Description: destructor (protected)
- //--------------------------------------------------------------------------------------
- EMSchur3DBase::~EMSchur3DBase ( ) {
- if (sigma) Delete3DScalar(sigma);
- if (sigmap) Delete3DScalar(sigmap);
- if (Cvec) delete [] Cvec;
- //if (CvecRe) delete [] CvecRe;
- } // ----- end of method EMSchur3DBase::~EMSchur3DBase (destructor) -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3D
- // Method: Serialize
- //--------------------------------------------------------------------------------------
- YAML::Node EMSchur3DBase::Serialize ( ) const {
- YAML::Node node = LemmaObject::Serialize();
- node.SetTag( this->GetName() );
- node["EMSchur3D_VERSION"] = EMSCHUR3D_VERSION;
- return node;
- } // ----- end of method EMSchur3D::Serialize -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: BuildC
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::BuildC ( Real*** sigmax, Real*** sigmay, Real*** sigmaz, const int& iw) {
-
- Cvec[iw].resize( unx+uny+unz , unx+uny+unz );
- Cvec[iw].reserve(Eigen::VectorXi::Constant(unx+uny+unz, 4)); // Upper/Lower
- //Cvec[iw].reserve(Eigen::VectorXi::Constant(unx+uny+unz, 7)); // Whole
-
- //Cvec_s.resize( idx. )
- //CMMvec[iw].resize( unx+uny+unz , unx+uny+unz );
- //CMMvec[iw].reserve(Eigen::VectorXi::Constant(unx+uny+unz, 1));
-
- Real omega = Omegas[iw];
- std::cout << "Building C_" << iw << std::endl;
-
- //Complex hsig(0);
- Real hsig(0);
- Real EPS(1e-24);
- Real EPS2(1e-18);
-
- // book keeping
- int ir = 0;
- int ic = 0;
-
- // LAPL{Ax} + i omega mu sigmax
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx+1; ++ix) {
-
- // Calculate 1/2 step values on staggered grid
- Real alo_ihz2(0);
- Real ahi_ihz2(0);
- if (iz == 0) {
- ahi_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
- alo_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
- } else if (iz == nz-1) {
- alo_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
- ahi_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
- } else {
- alo_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
- ahi_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
- }
-
- Real alo_ihy2(0); // half step low jump in y
- Real ahi_ihy2(0); // half step high jump in y
- if (iy == 0) {
- ahi_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
- alo_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
- } else if (iy == ny-1) {
- alo_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
- ahi_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
- } else {
- alo_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
- ahi_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
- }
-
- /////////////////////////////////////////////////////////////////////////
- // Diagonal entry
- // Harmonically average sigmax
- if (ix == 0) {
- hsig = sigma[ix][iy][iz];
- } else if (ix == nx) {
- hsig = sigma[ix-1][iy][iz];
- } else {
- hsig = ((hx[ix]+hx[ix-1])/2.) *
- (1. / ( (hx[ix ] / (2.*sigma[ix ][iy][iz] + EPS)) +
- ( hx[ix-1] / (2.*sigma[ix-1][iy][iz] + EPS)) ) );
- }
- if (std::abs(hsig) < EPS2) hsig = 0;
- //hsig += Complex(0, omega*EPSILON0);
-
- Real Sum(0);
- if (ix == 0) {
- // Dirichlet bdry condition on A
- // Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + 2.*ihx2[ix ];
- // Neumann bdry
- Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + ihx2[ix ];
- } else if (ix == nx) {
- // Dirichlet bdry condition on A
- // Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + 2.*ihx2[ix-1];
- // Neumann bdry
- Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + ihx2[ix-1];
- } else {
- Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + ihx2[ix] + ihx2[ix-1];
- }
-
- /////////////////////////////////////////
- // minus side off diagonal terms
- #if LOWER
- // Third Off Diagonal
- if (iz!=0) Cvec[iw].insert(ir, ic-ny*(nx+1)) = Complex(-ahi_ihz2, 0);
-
- // Second Off Diagonal
- if (iy!=0) Cvec[iw].insert(ir, ic-(nx+1)) = Complex(-ahi_ihy2, 0);
-
- // First Off Diagonal
- if (ix!=0) Cvec[iw].insert(ir, ic-1) = Complex(-ihx2[ix-1], 0);
- #endif
- ////////////////////////////////////////////
- // Diagonal Term
- Cvec[iw].insert(ir,ic) = Complex(Sum, 0) + Complex(0,1)*omega*MU0*hsig;
- //Cvec[iw].insert(ir,ic) = Complex(Sum, omega*MU0*hsig);
- //CMMvec[iw].insert(ir,ic) = 1./ Complex(Sum, omega*MU0*hsig);
-
- ////////////////////////////////////////////////////////////////////////
- // plus side off diagonal terms
- #if UPPER
- // First Off Diagonal
- if (ix!=nx) Cvec[iw].insert(ir, ic+1) = Complex(-ihx2[ix], 0);
-
- // Second Off Diagonal
- if (iy!=ny-1) Cvec[iw].insert(ir, ic+(nx+1)) = Complex(-ahi_ihy2, 0);
-
- // Third Off Diagonal
- if (iz!=nz-1) Cvec[iw].insert(ir, ic+ny*(nx+1)) = Complex(-ahi_ihz2, 0);
- #endif
- ++ir;
- ++ic;
- }
- }
- }
- assert(ic == unx);
-
- // LAPL{Ay} + i omega mu sigmay
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny+1; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
-
- // Calculate 1/2 step values on staggered grid
- Real alo_ihz2(0);
- Real ahi_ihz2(0);
- if (iz == 0) {
- ahi_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
- alo_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
- } else if (iz == nz-1) {
- alo_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
- ahi_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
- } else {
- alo_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
- ahi_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
- }
-
- // Calculate 1/2 step values on staggered grid
- Real alo_ihx2(0);
- Real ahi_ihx2(0);
- if (ix == 0) {
- ahi_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
- alo_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
- } else if (ix == nx-1) {
- alo_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
- ahi_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
- } else {
- alo_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
- ahi_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
- }
-
- // Harmonically average sigmay
- if (iy == 0) {
- hsig = sigma[ix][iy][iz];
- } else if (iy == ny) {
- hsig = sigma[ix][iy-1][iz];
- } else {
- hsig = ((hy[iy]+hy[iy-1])/2.) *
- (1. / ( (hy[iy ] / (2.*sigma[ix][iy ][iz] + EPS)) +
- ( hy[iy-1] / (2.*sigma[ix][iy-1][iz] + EPS)) ) );
- }
- if (std::abs(hsig) < EPS2) hsig = 0;
- //hsig += Complex(0, omega*EPSILON0);
-
- Real Sum(0);
- if (iy == 0) {
- // Dirichlet bdry condition on A
- //Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + 2.*ihy2[iy ];
- // Neumann bdry
- Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + ihy2[iy ];
- } else if (iy == ny) {
- // Dirichlet bdry condition on A
- //Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + 2.*ihy2[iy-1];
- // Neumann bdry
- Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + ihy2[iy-1];
- } else {
- Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + ihy2[iy] + ihy2[iy-1];
- }
- #if LOWER
- // Third Off Diagonal
- if (iz!=0) Cvec[iw].insert(ir,ic-(ny+1)*(nx)) = Complex(-ahi_ihz2, 0);
-
- // Second Off Diagonal
- if (iy!=0) Cvec[iw].insert(ir,ic-(nx)) = Complex(-ihy2[iy-1], 0);
-
- // First Off Diagonal
- if (ix!=0) Cvec[iw].insert(ir,ic-1) = Complex(-ahi_ihx2, 0);
- #endif
- // Diagonal Term
- Cvec[iw].insert(ir,ic) = Complex(Sum, 0) + Complex(0,1)*omega*MU0*hsig;
- //Cvec[iw].insert(ir,ic) = Complex(Sum, omega*MU0*hsig);
- //CMMvec[iw].insert(ir,ic) = 1./ Complex(Sum, omega*MU0*hsig);
- #if UPPER
- // First Off Diagonal
- if (ix!=nx-1) Cvec[iw].insert(ir,ic+1) = Complex(-ahi_ihx2, 0);
-
- // Second Off Diagonal
- if (iy!=ny) Cvec[iw].insert(ir,ic+(nx)) = Complex(-ihy2[iy], 0);
-
- // Third Off Diagonal
- if (iz!=nz-1) Cvec[iw].insert(ir,ic+(ny+1)*(nx)) = Complex(-ahi_ihz2, 0);
- #endif
- ++ir;
- ++ic;
- }
- }
- }
- assert(ic == unx+uny);
-
- // LAPL{Az} + i omega mu sigmaz
- for (int iz=0; iz<nz+1; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
-
- // Calculate 1/2 step values on staggered grid
- Real alo_ihx2(0);
- Real ahi_ihx2(0);
- if (ix == 0) {
- ahi_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
- alo_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
- } else if (ix == nx-1) {
- alo_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
- ahi_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
- } else {
- alo_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
- ahi_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
- }
-
- // Calculate 1/2 step values on staggered grid
- Real alo_ihy2(0);
- Real ahi_ihy2(0);
- if (iy == 0) {
- ahi_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
- alo_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
- } else if (iy == ny-1) {
- alo_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
- ahi_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
- } else {
- alo_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
- ahi_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
- }
-
- // Harmonically average sigmaz
- if (iz==0) {
- hsig = sigma[ix][iy][iz];
- } else if (iz == nz) {
- hsig = sigma[ix][iy][iz-1];
- } else {
- hsig = ((hz[iz]+hz[iz-1])/2.) *
- (1. / ( (hz[iz ] / (2.*sigma[ix][iy][iz ] + EPS)) +
- ( hz[iz-1] / (2.*sigma[ix][iy][iz-1] + EPS)) ) );
- }
- if (std::abs(hsig) < EPS2) hsig = 0;
- //hsig += Complex(0, omega*EPSILON0);
-
- Real Sum(0);
- if (iz == 0) {
- // Dirichlet bdry condition on A
- //Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + 2.*ihz2[iz ];
- // Neumann bdry
- Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + ihz2[iz ];
- } else if (iz == nz) {
- // Dirichlet bdry condition on A
- Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + 2.*ihz2[iz-1];
- // Neumann bdry
- //Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + ihz2[iz-1];
- } else {
- Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + ihz2[iz] + ihz2[iz-1];
- }
- #if LOWER
- // Third Off Diagonal
- if (iz!=0) Cvec[iw].insert(ir,ic-ny*(nx)) = Complex(-ihz2[iz-1], 0);
-
- // Second Off Diagonal
- if (iy!=0) Cvec[iw].insert(ir,ic-(nx)) = Complex(-ahi_ihy2, 0);
-
- // First Off Diagonal
- if (ix!=0) Cvec[iw].insert(ir,ic-1) = Complex(-ahi_ihx2, 0);
- #endif
- // Diagonal Term
- Cvec[iw].insert(ir,ic) = Complex(Sum, 0) + Complex(0,1)*omega*MU0*hsig;
- //Cvec[iw].insert(ir,ic) = Complex(Sum, omega*MU0*hsig);
- //CMMvec[iw].insert(ir,ic) = 1. / Complex(Sum, omega*MU0*hsig);
- #if UPPER
- // First Off Diagonal
- if (ix!=nx-1) Cvec[iw].insert(ir,ic+1) = Complex(-ahi_ihx2, 0);
-
- // Second Off Diagonal
- if (iy!=ny-1) Cvec[iw].insert(ir,ic+(nx)) = Complex(-ahi_ihy2, 0);
-
- // Third Off Diagonal
- if (iz!=nz) Cvec[iw].insert(ir,ic+ny*(nx)) = Complex(-ihz2[iz], 0);
- #endif
- ++ir;
- ++ic;
- }
- }
- }
- assert(ic == unx+uny+unz);
- //Cvec[iw].makeCompressed();
- Cvec[iw].makeCompressed();
- //CMMvec[iw].makeCompressed();
- //CsymVec[iw] = Cvec[iw].selfadjointView<Eigen::Upper>();
- return ;
- } // ----- end of method EMSchur3DBase::BuildC -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: BuildCReal
- //--------------------------------------------------------------------------------------
- /*
- void EMSchur3DBase::BuildCReal ( Real*** sigmax, Real*** sigmay, Real*** sigmaz, const int& iw) {
-
- int CI = unx+uny+unz;
- CvecRe[iw].resize( 2*CI, 2*CI ); // twice the size, but Real
-
- std::vector<Tr> triplets;
- triplets.reserve( 9*CI );
-
- Real omega = Omegas[iw];
- std::cout << "Building real C_" << iw << std::endl;
-
- // LAPL{Ax} + i omega mu sigmax
- int ir = 0;
- int ic = 0;
- Real hsig(0);
- Real EPS(1e-24);
- Real EPS2(1e-18);
-
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx+1; ++ix) {
-
- // Calculate 1/2 step values on staggered grid
- Real alo_ihz2(0);
- Real ahi_ihz2(0);
- if (iz == 0) {
- ahi_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
- alo_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
- } else if (iz == nz-1) {
- alo_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
- ahi_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
- } else {
- alo_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
- ahi_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
- }
-
- Real alo_ihy2(0); // half step low jump in y
- Real ahi_ihy2(0); // half step high jump in y
- if (iy == 0) {
- ahi_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
- alo_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
- } else if (iy == ny-1) {
- alo_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
- ahi_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
- } else {
- alo_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
- ahi_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
- }
-
- /////////////////////////////////////////////////////////////////////////
- // Diagonal entry
- // Harmonically average sigmax
- if (ix == 0) {
- hsig = sigma[ix][iy][iz];
- } else if (ix == nx) {
- hsig = sigma[ix-1][iy][iz];
- } else {
- hsig = ((hx[ix]+hx[ix-1])/2.) *
- (1. / ( (hx[ix ] / (2.*sigma[ix ][iy][iz] + EPS)) +
- ( hx[ix-1] / (2.*sigma[ix-1][iy][iz] + EPS)) ) );
- }
- if (std::abs(hsig) < EPS2) hsig = 0;
- //hsig += Complex(0, omega*EPSILON0);
-
- Real Sum(0);
- if (ix == 0) {
- // Dirichlet bdry condition on A
- // Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + 2.*ihx2[ix ];
- // Neumann bdry
- Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + ihx2[ix ];
- } else if (ix == nx) {
- // Dirichlet bdry condition on A
- // Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + 2.*ihx2[ix-1];
- // Neumann bdry
- Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + ihx2[ix-1];
- } else {
- Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + ihx2[ix] + ihx2[ix-1];
- }
-
- ////////////////////////////////////////////
- // Diagonal Term
- triplets.push_back( Tr(ir, ic , Sum) );
- triplets.push_back( Tr(ir+CI, ic+CI, -Sum) );
-
- ////////////////////////////////////////////////////////////////////////
- // plus side off diagonal terms
- // First Off Diagonal
- if (ix!=nx) {
- triplets.push_back( Tr(ir, ic+1 , -ihx2[ix]) );
- triplets.push_back( Tr(ir+CI, ic+1+CI, ihx2[ix]) );
- }
- // Second Off Diagonal
- if (iy!=ny-1) {
- triplets.push_back( Tr(ir , ic+(nx+1) , -ahi_ihy2) );
- triplets.push_back( Tr(ir+CI, ic+(nx+1)+CI, ahi_ihy2) );
- }
- // Third Off Diagonal
- if (iz!=nz-1) {
- triplets.push_back( Tr(ir , ic+ny*(nx+1) , -ahi_ihz2) );
- triplets.push_back( Tr(ir+CI, ic+ny*(nx+1)+CI, ahi_ihz2) );
- }
-
- ////////////////////////////////////////////////////////////////////////
- // imaginary part
- triplets.push_back( Tr(ir, ic+CI, omega*MU0*hsig) );
-
- ++ir;
- ++ic;
- }
- }
- }
- assert(ic == unx);
-
- // LAPL{Ay} + i omega mu sigmay
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny+1; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
-
- // Calculate 1/2 step values on staggered grid
- Real alo_ihz2(0);
- Real ahi_ihz2(0);
- if (iz == 0) {
- ahi_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
- alo_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
- } else if (iz == nz-1) {
- alo_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
- ahi_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
- } else {
- alo_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
- ahi_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
- }
-
- // Calculate 1/2 step values on staggered grid
- Real alo_ihx2(0);
- Real ahi_ihx2(0);
- if (ix == 0) {
- ahi_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
- alo_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
- } else if (ix == nx-1) {
- alo_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
- ahi_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
- } else {
- alo_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
- ahi_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
- }
-
- // Harmonically average sigmay
- if (iy == 0) {
- hsig = sigma[ix][iy][iz];
- } else if (iy == ny) {
- hsig = sigma[ix][iy-1][iz];
- } else {
- hsig = ((hy[iy]+hy[iy-1])/2.) *
- (1. / ( (hy[iy ] / (2.*sigma[ix][iy ][iz] + EPS)) +
- ( hy[iy-1] / (2.*sigma[ix][iy-1][iz] + EPS)) ) );
- }
- if (std::abs(hsig) < EPS2) hsig = 0;
- //hsig += Complex(0, omega*EPSILON0);
-
- Real Sum(0);
- if (iy == 0) {
- // Dirichlet bdry condition on A
- //Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + 2.*ihy2[iy ];
- // Neumann bdry
- Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + ihy2[iy ];
- } else if (iy == ny) {
- // Dirichlet bdry condition on A
- //Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + 2.*ihy2[iy-1];
- // Neumann bdry
- Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + ihy2[iy-1];
- } else {
- Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + ihy2[iy] + ihy2[iy-1];
- }
-
- //////////////////////////////////////////////////////////////////////////////
- // Diagonal Term
- triplets.push_back( Tr(ir, ic , Sum) );
- triplets.push_back( Tr(ir+CI, ic+CI, -Sum) );
-
- // First Off Diagonal
- if (ix!=nx-1) {
- triplets.push_back( Tr(ir, ic+1 , -ahi_ihx2));
- triplets.push_back( Tr(ir+CI, ic+1+CI, ahi_ihx2));
- }
-
- // Second Off Diagonal
- if (iy!=ny) {
- triplets.push_back( Tr(ir , ic+(nx) , -ihy2[iy]) );
- triplets.push_back( Tr(ir+CI, ic+(nx)+CI, ihy2[iy]) );
- }
-
- // Third Off Diagonal
- if (iz!=nz-1) {
- triplets.push_back( Tr(ir , ic+(ny+1)*(nx) , -ahi_ihz2));
- triplets.push_back( Tr(ir+CI, ic+(ny+1)*(nx)+CI, ahi_ihz2));
- }
-
- ///////////////////////////////////////////////////////////////////////
- // imaginary part
- triplets.push_back( Tr(ir, ic+CI, omega*MU0*hsig) );
-
- ++ir;
- ++ic;
- }
- }
- }
- assert(ic == unx+uny);
-
- // LAPL{Az} + i omega mu sigmaz
- for (int iz=0; iz<nz+1; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
-
- // Calculate 1/2 step values on staggered grid
- Real alo_ihx2(0);
- Real ahi_ihx2(0);
- if (ix == 0) {
- ahi_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
- alo_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
- } else if (ix == nx-1) {
- alo_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
- ahi_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
- } else {
- alo_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
- ahi_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
- }
-
- // Calculate 1/2 step values on staggered grid
- Real alo_ihy2(0);
- Real ahi_ihy2(0);
- if (iy == 0) {
- ahi_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
- alo_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
- } else if (iy == ny-1) {
- alo_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
- ahi_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
- } else {
- alo_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
- ahi_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
- }
-
- // Harmonically average sigmaz
- if (iz==0) {
- hsig = sigma[ix][iy][iz];
- } else if (iz == nz) {
- hsig = sigma[ix][iy][iz-1];
- } else {
- hsig = ((hz[iz]+hz[iz-1])/2.) *
- (1. / ( (hz[iz ] / (2.*sigma[ix][iy][iz ] + EPS)) +
- ( hz[iz-1] / (2.*sigma[ix][iy][iz-1] + EPS)) ) );
- }
- if (std::abs(hsig) < EPS2) hsig = 0;
- //hsig += Complex(0, omega*EPSILON0);
-
- Real Sum(0);
- if (iz == 0) {
- // Dirichlet bdry condition on A
- //Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + 2.*ihz2[iz ];
- // Neumann bdry
- Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + ihz2[iz ];
- } else if (iz == nz) {
- // Dirichlet bdry condition on A
- //Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + 2.*ihz2[iz-1];
- // Neumann bdry
- Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + ihz2[iz-1];
- } else {
- Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + ihz2[iz] + ihz2[iz-1];
- }
-
- //////////////////////////////////////////////////////////////////
- // Diagonal Term
- triplets.push_back( Tr(ir , ic , Sum) );
- triplets.push_back( Tr(ir+CI, ic+CI, -Sum) );
-
- // First Off Diagonal
- if (ix!=nx-1) {
- triplets.push_back( Tr(ir , ic+1 , -ahi_ihx2) );
- triplets.push_back( Tr(ir+CI, ic+1+CI, ahi_ihx2) );
- }
- // Second Off Diagonal
- if (iy!=ny-1) {
- triplets.push_back( Tr(ir , ic+(nx) , -ahi_ihy2) );
- triplets.push_back( Tr(ir+CI, ic+(nx)+CI, ahi_ihy2) );
- }
- // Third Off Diagonal
- if (iz!=nz) {
- triplets.push_back( Tr(ir , ic+ny*(nx) , -ihz2[iz]) );
- triplets.push_back( Tr(ir+CI, ic+ny*(nx)+CI, ihz2[iz]) );
- }
-
- //////////////////////////////////////////////////////
- // imaginary part
- triplets.push_back( Tr(ir, ic+CI, omega*MU0*hsig) );
-
- ++ir;
- ++ic;
- }
- }
- }
- assert(ic == unx+uny+unz);
- CvecRe[iw].setFromTriplets(triplets.begin(), triplets.end()) ;
- CvecRe[iw].makeCompressed();
-
- return ;
- } // ----- end of method EMSchur3DBase::BuildCReal -----
- */
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: SetResFileName
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::SetResFileName ( const std::string& fname ) {
- ResFile = fname;
- return ;
- } // ----- end of method EMSchur3DBase::SetResFileName -----
-
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: BuildCPreconditioner
- //--------------------------------------------------------------------------------------
- // void EMSchur3DBase::BuildCPreconditioner ( const int& iw ) {
- // std::cout << "Building preconditioner for C_" << iw << std::endl;
- // CMvec[iw].setDroptol(1e-4); // 1e-4 previously (tested value)
- // CMvec[iw].setFillfactor(10);
- // Eigen::SparseMatrix<Complex> Csym;
- // Csym = Cvec[iw].selfadjointView<Eigen::Upper>();
- // CMvec[iw].compute( Csym );
- // return ;
- // } // ----- end of method EMSchur3DBase::BuildCPreconditioner -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: BuildD
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::BuildD ( ) {
- D.resize(uns, unx+uny+unz );
- //D.reserve(Eigen::VectorXi::Constant(unx+uny+unz, 6));
-
- std::vector<Tc> triplets;
- triplets.reserve( 2*unx + 2*uny + 2*unz );
-
- std::cout << "Building D...";
- int ic = 0;
- int ir = 0;
-
- ////////////////////////////////
- // Ax dx
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- ++ic; // TRICKY
- for (int ix=0; ix<nx; ++ix) {
- //D.insert(ir, ic-1) = Complex(-ihx[ix], 0);
- //D.insert(ir, ic ) = Complex( ihx[ix], 0);
- triplets.push_back( Tc(ir, ic-1, Complex(-ihx[ix], 0)) );
- triplets.push_back( Tc(ir, ic , Complex( ihx[ix], 0)) );
- ++ir;
- ++ic;
- }
- }
- }
- assert (ic==unx);
-
- ///////////////////////////////
- // Ay dy
- ir = 0;
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
- //D.insert(ir, ic ) = Complex(-ihy[iy], 0);
- //D.insert(ir, ic+nx) = Complex( ihy[iy], 0);
- triplets.push_back( Tc(ir, ic , Complex(-ihy[iy], 0)) );
- triplets.push_back( Tc(ir, ic+nx, Complex( ihy[iy], 0)) );
- ++ir;
- ++ic;
- }
- }
- ic += nx;
- }
- assert (ic==unx+uny);
-
- ///////////////////////////////
- // Az dz
- ir = 0;
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
- //D.insert(ir, ic ) = Complex(-ihz[iz], 0);
- //D.insert(ir, ic+nx*ny) = Complex( ihz[iz], 0);
- triplets.push_back( Tc(ir, ic , Complex(-ihz[iz], 0)) );
- triplets.push_back( Tc(ir, ic+nx*ny, Complex( ihz[iz], 0)) );
- ++ir;
- ++ic;
- }
- }
- }
- assert (ic+nx*ny==unx+uny+unz);
-
- D.setFromTriplets(triplets.begin(), triplets.end());
- D.makeCompressed();
- std::cout << "Done!" << std::endl;
- return ;
- } // ----- end of method EMSchur3DBase::BuildD -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: SetRectilinearGrid
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::SetRectilinearGrid ( std::shared_ptr<RectilinearGrid> GridPtr ) {
-
- Grid = GridPtr;
-
- // Do some convienience stuff, we call these so often it's nice to not have to use
- // the accessors every time. This is the only place these are set.
- nx = Grid->GetNx();
- ny = Grid->GetNy();
- nz = Grid->GetNz();
-
- hx = Grid->GetDx();
- hy = Grid->GetDy();
- hz = Grid->GetDz();
-
- unx = (nx+1)*ny*nz; // Number of Vectors x component
- uny = nx*(ny+1)*nz; // Number of Vectors y component
- unz = nx*ny*(nz+1); // Number of Vectors z component
- uns = nx*ny*nz; // On dual grid, number of scalars
-
- ihx = 1. / hx.array();
- ihy = 1. / hy.array();
- ihz = 1. / hz.array();
-
- ihx2 = 1. / (hx.array() * hx.array());
- ihy2 = 1. / (hy.array() * hy.array());
- ihz2 = 1. / (hz.array() * hz.array());
-
- return ;
- } // ----- end of method EMSchur3DBase::SetRectilinearGrid -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: SetRectilinearGrid
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::SetRectilinearGrid ( vtkRectilinearGrid* rGridPtr ) {
-
- // Convert to Lemma here
-
-
- } // ----- end of method EMSchur3DBase::SetRectilinearGrid -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: Setup
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::Setup ( ) {
-
- ////////////////////////////////////
- // First set up grid stuff
- if (Cvec) delete [] Cvec;
- //if (CSolver) delete [] CSolver;
- Cvec = new Eigen::SparseMatrix<Complex, Eigen::RowMajor> [Omegas.size()];
-
- //#ifdef LEMMAUSEOMP
- //#pragma omp parallel for schedule(static, 1)
- //#endif
- for (int iw=0; iw<Omegas.size(); ++iw) {
- std::cout << "Build C " << std::endl;
- BuildC(sigma, sigma, sigma, iw);
- // std::cout << "Build Re(C) " << std::endl;
- // BuildCReal(sigma, sigma, sigma, iw);
- }
-
- BuildCDirectSolver( ); // templated specialisation
- BuildD();
-
- return ;
- } // ----- end of method EMSchur3DBase::Setup -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: SetLayeredEarthEM
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::SetLayeredEarthEM ( std::shared_ptr<LayeredEarthEM> LayModelin ) {
- LayModel = LayModelin;
- return ;
- } // ----- end of method EMSchur3DBase::SetLayeredEarthEM -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: PrimaryField
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::PrimaryField ( std::shared_ptr<DipoleSource> source, std::shared_ptr<FieldPoints> dpoint ) {
-
- auto EmEarth = Lemma::EMEarth1D::NewSP();
- EmEarth->AttachDipoleSource(source);
- EmEarth->AttachLayeredEarthEM(LayModel);
- EmEarth->AttachFieldPoints(dpoint);
- EmEarth->SetFieldsToCalculate(Lemma::E);
- EmEarth->SetHankelTransformMethod(ANDERSON801);
- std::cout << "Primary Field Calculation" << std::endl;
- EmEarth->MakeCalc3();
- return ;
- } // ----- end of method EMSchur3DBase::PrimaryField -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: FillPoints
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::FillPoints ( std::shared_ptr<FieldPoints> dpoint) {
-
- dpoint->SetNumberOfPoints(unx + uny + unz);
- int ip(0);
- Real tx, ty, tz;
- Real ox = Grid->GetOx();
- Real oy = Grid->GetOy();
- Real oz = Grid->GetOz();
- tz = oz;
- for (int iz=0; iz<nz; ++iz) {
- ty = oy;
- for (int iy=0; iy<ny; ++iy) {
- tx = ox - hx(0)/2.;
- for (int ix=0; ix<nx+1; ++ix) {
- dpoint->SetLocation(ip, tx, ty, tz);
- if (ix < nx) tx += hx(ix);
- ++ip;
- }
- if (iy<ny-1) ty += .5*hy(iy) + .5*hy(iy+1);
- }
- if (iz<nz-1) tz += .5*hz(iz) + .5*hz(iz+1);
- }
-
- tz = oz;
- for (int iz=0; iz<nz; ++iz) {
- ty = oy - hy(0)/2.;
- for (int iy=0; iy<ny+1; ++iy) {
- tx = ox;
- for (int ix=0; ix<nx; ++ix) {
- dpoint->SetLocation(ip, tx, ty, tz);
- if (ix<nx-1) tx += .5*hx(ix) + .5*hx(ix+1);
- ++ip;
- }
- if (iy < ny) ty += hy(iy);
- }
- if (iz<nz-1) tz += .5*hz(iz) + .5*hz(iz+1);
- }
-
- tz = oz - hz(0)/2.;
- for (int iz=0; iz<nz+1; ++iz) {
- ty = oy;
- for (int iy=0; iy<ny; ++iy) {
- tx = ox;
- for (int ix=0; ix<nx; ++ix) {
- dpoint->SetLocation(ip, tx, ty, tz);
- if (ix<nx-1) tx += .5*hx(ix) + .5*hx(ix+1);
- ++ip;
- }
- if (iy<ny-1) ty += .5*hy(iy) + .5*hy(iy+1);
- }
- if (iz < nz) tz += hz(iz);
- }
-
- return ;
- } // ----- end of method EMSchur3DBase::FillPoints -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: FillSourceTerms
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::FillSourceTerms ( Eigen::Ref<VectorXcr> ms,
- Eigen::Ref<VectorXcr> Se, Eigen::Ref<VectorXcr> E0,
- std::shared_ptr<FieldPoints> dpoint, const Real& omega ) {
-
- //Complex hsig(0);
- //Complex hsigp(0);
- Real hsig(0);
- Real hsigp(0);
- Real EPS(0); //1e-24);
- Real EPS2(1e-18);
- int iv = 0;
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx+1; ++ix) {
-
- // Harmonically average sigma
- if (ix == 0) {
- hsig = sigma[ix][iy][iz];
- hsigp = sigmap[ix][iy][iz];
- } else if (ix == nx) {
- hsig = sigma[ix-1][iy][iz];
- hsigp = sigmap[ix-1][iy][iz];
- } else {
- hsig = ((hx[ix]+hx[ix-1])/2.) *
- (1. / ( (hx[ix ] / (2.*sigma[ix ][iy][iz] + EPS)) +
- ( hx[ix-1] / (2.*sigma[ix-1][iy][iz] + EPS)) ) );
- hsigp = ((hx[ix]+hx[ix-1])/2.) *
- (1. / ( (hx[ix ] / (2.*sigmap[ix ][iy][iz] + EPS)) +
- ( hx[ix-1] / (2.*sigmap[ix-1][iy][iz] + EPS)) ) );
- }
- if (std::abs(hsig) < EPS2) hsig = 0;
- if (std::abs(hsigp) < EPS2) hsigp = 0;
- //hsig += Complex(0, omega*EPSILON0);
- //ioms(iv) = Complex(0, omega*MU0*hsig);
- //Se(iv) = ioms(iv)*dpoint->GetEfield(0,iv)(0);
- ms(iv) = MU0*hsig;
- Se(iv) = MU0*hsigp*
- dpoint->GetEfield(0,iv)(0);
- E0(iv) = dpoint->GetEfield(0,iv)(0);
- ++iv;
- }
- }
- }
-
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny+1; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
-
- // Harmonically average sigma
- if (iy == 0) {
- hsig = sigma[ix][iy][iz];
- hsigp = sigmap[ix][iy][iz];
- } else if (iy == ny) {
- hsig = sigma[ix][iy-1][iz];
- hsigp = sigmap[ix][iy-1][iz];
- } else {
- hsig = ((hy[iy]+hy[iy-1])/2.) *
- (1. / ( (hy[iy ] / (2.*sigma[ix][iy ][iz] + EPS)) +
- ( hy[iy-1] / (2.*sigma[ix][iy-1][iz] + EPS)) ) );
- hsigp = ((hy[iy]+hy[iy-1])/2.) *
- (1. / ( (hy[iy ] / (2.*sigmap[ix][iy ][iz] + EPS)) +
- ( hy[iy-1] / (2.*sigmap[ix][iy-1][iz] + EPS)) ) );
- }
- if (std::abs(hsig) < EPS2) hsig = 0;
- if (std::abs(hsigp) < EPS2) hsigp = 0;
- //hsig += Complex(0, omega*EPSILON0);
- //ioms(iv) = Complex(0, omega*MU0*hsig);
- //Se(iv) = ioms(iv) * dpoint->GetEfield(0,iv)(1);
- ms(iv) = MU0*hsig;
- Se(iv) = MU0*hsigp* // TODO was hsigp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- dpoint->GetEfield(0,iv)(1);
- E0(iv) = dpoint->GetEfield(0,iv)(1);
- ++iv;
- }
- }
- }
-
- for (int iz=0; iz<nz+1; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
-
- // Harmonically average sigma
- if (iz==0) {
- hsig = sigma[ix][iy][iz];
- hsigp = sigma[ix][iy][iz];
- } else if (iz == nz) {
- hsig = sigma[ix][iy][nz-1];
- hsigp = sigmap[ix][iy][nz-1];
- } else {
- hsig = ((hz[iz]+hz[iz-1])/2.) *
- (1. / ( (hz[iz ] / (2.*sigma[ix][iy][iz ] + EPS)) +
- ( hz[iz-1] / (2.*sigma[ix][iy][iz-1] + EPS)) ) );
- hsigp = ((hz[iz]+hz[iz-1])/2.) *
- (1. / ( (hz[iz ] / (2.*sigmap[ix][iy][iz ] + EPS)) +
- ( hz[iz-1] / (2.*sigmap[ix][iy][iz-1] + EPS)) ) );
- }
- if (std::abs(hsig) < EPS2) hsig = 0;
- if (std::abs(hsigp) < EPS2) hsigp = 0;
- //hsig += Complex(0, omega*EPSILON0);
- //ioms(iv) = Complex(0, omega*MU0*hsig);
- //Se(iv) = ioms(iv) * dpoint->GetEfield(0, iv)(2);
- ms(iv) = MU0*hsig;
- Se(iv) = MU0*hsig
- *dpoint->GetEfield(0, iv)(2);
- E0(iv) = dpoint->GetEfield(0, iv)(2);
- ++iv;
- }
- }
- }
-
- return ;
- } // ----- end of method EMSchur3DBase::FillSourceTerms -----
-
- VectorXcr EMSchur3DBase::StaggeredGridCurl( Eigen::Ref<VectorXcr> A ) {
- VectorXcr B = VectorXcr::Zero(3*nx*ny*nz);
- VectorXcr dzdy = VectorXcr::Zero( nx*ny*nz);
- VectorXcr dydz = VectorXcr::Zero( nx*ny*nz);
- VectorXcr dxdz = VectorXcr::Zero( nx*ny*nz);
- VectorXcr dzdx = VectorXcr::Zero( nx*ny*nz);
- VectorXcr dydx = VectorXcr::Zero( nx*ny*nz);
- VectorXcr dxdy = VectorXcr::Zero( nx*ny*nz);
- /*
- Curl_x = dzdy - dydz
- Curl_y = dxdz - dzdx
- Curl_z = dydx - dxdy
- */
- ////////////////////////////////
- // X component
- int ii = 0;
- int iix = 0;
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
- dxdz(iix) = ( A(ii+(nx+1)*ny) - A(ii) ) / hz[iz];
- dxdy(iix) = ( A(ii+(nx+1) ) - A(ii) ) / hy[iy];
- ++ii;
- ++iix;
- }
- // Jump around boundary
- ++ii;
- }
- }
-
- ///////////////////////////////
- // Y component
- int iiy = 0;
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
- dydz(iiy) = ( A(ii+(nx)*(ny+1)) - A(ii) ) / hz[iz];
- dydx(iiy) = ( A(ii+1 ) - A(ii) ) / hx[ix];
- ++iiy;
- ++ii;
- }
- }
- // Jump around boundary
- ii += nx;
- }
-
- ////////////////////////////
- // Z component
- int iiz = 0;
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
- //Az(iiz) = 0.5 * ( A(ii) + A(ii+nx*ny));
- dzdy(iiz) = ( A(ii+nx ) - A(ii) ) / hy[iy];
- dzdx(iiz) = ( A(ii+1 ) - A(ii) ) / hx[ix];
- ++iiz;
- ++ii;
- }
- //++iiz;
- //ii += nx*ny;
- }
- }
-
- B.segment( 0, nx*ny*nz) = dzdy - dydz;
- B.segment( nx*ny*nz, nx*ny*nz) = dxdz - dzdx;
- B.segment(2*nx*ny*nz, nx*ny*nz) = dydx - dxdy;
-
- return B;
- }
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: WriteVTKResults
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::WriteVTKResults ( const std::string& fname, Eigen::Ref<VectorXcr> A,
- Eigen::Ref<VectorXcr> Se, Eigen::Ref<VectorXcr> E0, Eigen::Ref<VectorXcr> E,
- Eigen::Ref<VectorXcr> Phi, Eigen::Ref<VectorXcr> ADiv, Eigen::Ref<VectorXcr> ADiv2,
- Eigen::Ref<VectorXcr> B ) {
-
- auto VTKGridExporter = RectilinearGridVTKExporter::NewSP();
- VTKGridExporter->SetGrid(Grid);
-
- // VTK arrays to fill
- // scalars
- vtkDoubleArray *sigmaArray = vtkDoubleArray::New();
- vtkDoubleArray *sigmapArray = vtkDoubleArray::New();
- vtkDoubleArray *phiRealArray = vtkDoubleArray::New();
- vtkDoubleArray *phiImagArray = vtkDoubleArray::New();
- vtkDoubleArray *ADivRealArray = vtkDoubleArray::New();
- vtkDoubleArray *ADivImagArray = vtkDoubleArray::New();
- vtkDoubleArray *ADiv2RealArray = vtkDoubleArray::New();
- vtkDoubleArray *ADiv2ImagArray = vtkDoubleArray::New();
-
- // vectors
- vtkDoubleArray *AReal = vtkDoubleArray::New();
- vtkDoubleArray *AImag = vtkDoubleArray::New();
- vtkDoubleArray *BReal = vtkDoubleArray::New();
- vtkDoubleArray *BImag = vtkDoubleArray::New();
- vtkDoubleArray *SeReal = vtkDoubleArray::New();
- vtkDoubleArray *SeImag = vtkDoubleArray::New();
- vtkDoubleArray *E0Real = vtkDoubleArray::New();
- vtkDoubleArray *E0Imag = vtkDoubleArray::New();
- vtkDoubleArray *EReal = vtkDoubleArray::New();
- vtkDoubleArray *EImag = vtkDoubleArray::New();
-
- AReal->SetNumberOfComponents(3);
- AImag->SetNumberOfComponents(3);
- BReal->SetNumberOfComponents(3);
- BImag->SetNumberOfComponents(3);
- SeReal->SetNumberOfComponents(3);
- SeImag->SetNumberOfComponents(3);
- E0Real->SetNumberOfComponents(3);
- E0Imag->SetNumberOfComponents(3);
- EReal->SetNumberOfComponents(3);
- EImag->SetNumberOfComponents(3);
-
- // Interpolate
- Eigen::VectorXcd Ax(nx*ny*nz); // A
- Eigen::VectorXcd Ay(nx*ny*nz);
- Eigen::VectorXcd Az(nx*ny*nz);
- Eigen::VectorXcd Bx(nx*ny*nz); // B
- Eigen::VectorXcd By(nx*ny*nz);
- Eigen::VectorXcd Bz(nx*ny*nz);
- Eigen::VectorXcd Sex(nx*ny*nz); // Se
- Eigen::VectorXcd Sey(nx*ny*nz);
- Eigen::VectorXcd Sez(nx*ny*nz);
- Eigen::VectorXcd E0x(nx*ny*nz); // E0
- Eigen::VectorXcd E0y(nx*ny*nz);
- Eigen::VectorXcd E0z(nx*ny*nz);
- Eigen::VectorXcd Ex(nx*ny*nz); // E
- Eigen::VectorXcd Ey(nx*ny*nz);
- Eigen::VectorXcd Ez(nx*ny*nz);
-
- ////////////////////////////////
- // X component
- int ii = 0;
- int iix = 0;
- int ic = 0;
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
- Ax(iix) = 0.5 * ( A(ii) + A(ii+1));
- Sex(iix) = 0.5 * (Se(ii) + Se(ii+1));
- E0x(iix) = 0.5 * (E0(ii) + E0(ii+1));
- Ex(iix) = 0.5 * ( E(ii) + E(ii+1));
- Bx(iix) = B(ic);
- ++ii;
- ++iix;
- ++ic;
- }
- // Jump around boundary
- ++ii;
- }
- }
-
- ///////////////////////////////
- // Y component
- int iiy = 0;
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
- Ay(iiy) = 0.5 * ( A(ii) + A(ii+nx));
- Sey(iiy) = 0.5 * (Se(ii) + Se(ii+nx));
- E0y(iiy) = 0.5 * (E0(ii) + E0(ii+nx));
- Ey(iiy) = 0.5 * ( E(ii) + E(ii+nx));
- By(iiy) = B(ic);
- ++iiy;
- ++ii;
- ++ic;
- }
- }
- // Jump around boundary
- ii += nx;
- }
-
- ////////////////////////////
- // Z component
- int iiz = 0;
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
- Az(iiz) = 0.5 * ( A(ii) + A(ii+nx*ny));
- Sez(iiz) = 0.5 * (Se(ii) + Se(ii+nx*ny));
- E0z(iiz) = 0.5 * (E0(ii) + E0(ii+nx*ny));
- Ez(iiz) = 0.5 * ( E(ii) + E(ii+nx*ny));
- Bz(iiz) = B(ic);
- ++iiz;
- ++ii;
- ++ic;
- }
- //++iiz;
- //ii += nx*ny;
- }
- }
-
- //VTKGridExporter->GetVTKGrid();
- int i = 0;
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
- // scalars
- //sigmaArray->InsertTuple1(i, std::real(sigma[ix][iy][iz]) );
- //sigmapArray->InsertTuple1(i, std::real(sigmap[ix][iy][iz]) );
- sigmaArray->InsertTuple1(i, sigma[ix][iy][iz] );
- sigmapArray->InsertTuple1(i, sigmap[ix][iy][iz] );
-
- phiRealArray->InsertTuple1(i, std::real(Phi(i)));
- phiImagArray->InsertTuple1(i, std::imag(Phi(i)));
- ADivRealArray->InsertTuple1(i, std::real(ADiv(i)));
- ADivImagArray->InsertTuple1(i, std::imag(ADiv(i)));
- ADiv2RealArray->InsertTuple1(i, std::real(ADiv2(i)));
- ADiv2ImagArray->InsertTuple1(i, std::imag(ADiv2(i)));
-
- // vectors
- AReal->InsertTuple3(i, std::real(Ax(i)), std::real(Ay(i)), std::real(Az(i)));
- AImag->InsertTuple3(i, std::imag(Ax(i)), std::imag(Ay(i)), std::imag(Az(i)));
- BReal->InsertTuple3(i, std::real(Bx(i)), std::real(By(i)), std::real(Bz(i)));
- BImag->InsertTuple3(i, std::imag(Bx(i)), std::imag(By(i)), std::imag(Bz(i)));
- SeReal->InsertTuple3(i, std::real(Sex(i)), std::real(Sey(i)), std::real(Sez(i)));
- SeImag->InsertTuple3(i, std::imag(Sex(i)), std::imag(Sey(i)), std::imag(Sez(i)));
- E0Real->InsertTuple3(i, std::real(E0x(i)), std::real(E0y(i)), std::real(E0z(i)));
- E0Imag->InsertTuple3(i, std::imag(E0x(i)), std::imag(E0y(i)), std::imag(E0z(i)));
- EReal-> InsertTuple3(i, std::real(Ex(i)), std::real(Ey(i)), std::real(Ez(i)));
- EImag-> InsertTuple3(i, std::imag(Ex(i)), std::imag(Ey(i)), std::imag(Ez(i)));
-
- ++i;
- }
- }
- }
-
- AReal->SetName("A_real");
- AImag->SetName("A_imag");
- BReal->SetName("B_real");
- BImag->SetName("B_imag");
- SeReal->SetName("Se_real");
- SeImag->SetName("Se_imag");
- E0Real->SetName("E0_real");
- E0Imag->SetName("E0_imag");
- EReal->SetName("E_real");
- EImag->SetName("E_imag");
-
- phiRealArray->SetName("phi Real");
- phiImagArray->SetName("phi imag");
- ADivRealArray->SetName("ini div A real");
- ADivImagArray->SetName("ini div A imag");
- ADiv2RealArray->SetName("div A real");
- ADiv2ImagArray->SetName("div A imag");
-
- // Add scalar fields
- sigmaArray->SetName("sigma");
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(sigmaArray);
-
- sigmapArray->SetName("sigma prime");
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(sigmapArray);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(phiRealArray);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(phiImagArray);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(ADivRealArray);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(ADivImagArray);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(ADiv2RealArray);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(ADiv2ImagArray);
-
- // Add Vector Arrays
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(AReal);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(AImag);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(SeReal);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(SeImag);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(E0Real);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(E0Imag);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(EReal);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(EImag);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(BReal);
- VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(BImag);
-
- // Write
-
- vtkXMLRectilinearGridWriter *gridWrite = vtkXMLRectilinearGridWriter::New();
- gridWrite->SetInputData( VTKGridExporter->GetVTKGrid() );
- gridWrite->SetFileName( (fname + std::string(".vtr")).c_str());
- gridWrite->Update();
- gridWrite->Write();
-
- sigmaArray->Delete();
- gridWrite->Delete();
-
- phiRealArray->Delete();
- phiImagArray->Delete();
- ADivRealArray->Delete();
- ADivImagArray->Delete();
- ADiv2RealArray->Delete();
- ADiv2ImagArray->Delete();
-
- // TODO delete vectors too!
- AReal->Delete();
- AImag->Delete();
- BReal->Delete();
- BImag->Delete();
- SeReal->Delete();
- SeImag->Delete();
- E0Real->Delete();
- E0Imag->Delete();
- EReal->Delete();
- EImag->Delete();
-
- return ;
- } // ----- end of method EMSchur3DBase::WriteVTKResults -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: LoadMeshToolsConductivityModel
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::LoadMeshToolsConductivityModel ( const std::string& fname ) {
-
- if (sigma) Delete3DScalar(sigma);
- //Allocate3DScalar(sigma, Complex(0));
- Allocate3DScalar(sigma, 0.);
-
- // TODO, there are smarter ways to do this, since the bacground is 1D. But for now, we just have
- // them the same. Eases logic in harmonic averaging. And it's just in one thread, so not so
- // hateful
- if (!LayModel) {
- std::cerr << "In EMSchur3DBase::LoadMeshToolsConductivityModel, you need to set your bacground 1D model\n";
- exit(EXIT_FAILURE);
- }
-
- if (sigmap) Delete3DScalar(sigmap);
- //Allocate3DScalar(sigmap, Complex(0));
- Allocate3DScalar(sigmap, 0.);
-
-
- auto Parser = ASCIIParser::NewSP();
- Parser->Open(fname);
- std::vector<Real> sigmod = Parser->ReadReals( -1 );
- if ( static_cast<int>(sigmod.size()) != nx*ny*nz) {
- std::cerr << "In EMSchur3DBase::LoadMeshToolsConductivityModel model sizes are not in agreement\n";
- exit(EXIT_FAILURE);
- }
- Parser->Close();
-
-
- int ic(0);
- for (int ix=0; ix<nx; ++ix) {
- for (int iy=0; iy<ny; ++iy) {
- Real pz = Grid->GetOz();
- for (int iz=0; iz<nz; ++iz) {
- if ( sigmod[ic] != sigmod[ic] ) {
- std::cout << "sigmod problems " << ic << std::endl;
- }
- sigma[ix][iy][iz] = sigmod[ic];
- //sigmap[ix][iy][iz] = sigmod[ic] - LayModel->GetLayerConductivity(LayModel->GetLayerAtThisDepth(pz));
- sigmap[ix][iy][iz] = sigmod[ic] - LayModel->GetLayerConductivity(LayModel->GetLayerAtThisDepth(pz)).real() ;
-
- ++ic;
- pz += hz[iz];
- }
- }
- }
-
- // Marker and cell
- MAC = VectorXi::Zero(nx*ny*nz);
- idx.clear();
- ic = 0;
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
- if (sigma[ix][iy][iz] > 0.) {
- MAC(ic) = 1;
- idx.push_back(ic);
- }
- ++ic;
- }
- }
- }
-
- return ;
- }
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: SetAEMSurvey
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::SetAEMSurvey ( std::shared_ptr<AEMSurvey> SurveyIn ) {
- Survey = SurveyIn;
- // determine how many frequencies we are dealing with
- Omegas = 2.*PI*Survey->GetFrequencies();
- return ;
- } // ----- end of method EMSchur3DBase::SetAEMSurvey -----// ----- end of method EMSchur3DBase::LoadMeshToolsConductivityModel -----
-
- //--------------------------------------------------------------------------------------
- // Class: EMSchur3DBase
- // Method: Solve
- //--------------------------------------------------------------------------------------
- void EMSchur3DBase::Solve ( ) {
-
- // TODO, these should throw exceptions for a GUI to catch
- if (!LayModel or !Survey or !Grid or !sigma) {
- std::cerr << "You need to set something. EMSChur3D::Solve()";
- exit(EXIT_FAILURE);
- }
-
- //BuildD(); // strangely necessary, bug, track down.
- Setup();
-
- std::cout << "Entering parallel solve" << std::endl;
- //#ifdef LEMMAUSEOMP
- //#pragma omp parallel for schedule(static,1)
- //#endif
- for (int isource=0; isource<Survey->GetNumberOfSources(); ++isource) {
- SolveSource(Survey->GetSource(isource), isource);
- }
-
- return ;
- } // ----- end of method EMSchur3DBase::Solve -----
-
- } // ----- end of Lemma name -----
|