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.

PolygonalWireAntenna.cpp 11KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326
  1. /* This file is part of Lemma, a geophysical modelling and inversion API */
  2. /* This Source Code Form is subject to the terms of the Mozilla Public
  3. * License, v. 2.0. If a copy of the MPL was not distributed with this
  4. * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
  5. /**
  6. @file
  7. @author Trevor Irons
  8. @date 05/18/2010
  9. **/
  10. #include "PolygonalWireAntenna.h"
  11. namespace Lemma {
  12. std::ostream &operator << (std::ostream &stream, const PolygonalWireAntenna &ob) {
  13. stream << ob.Serialize() << "\n";
  14. return stream;
  15. }
  16. // ==================== LIFECYCLE =======================
  17. PolygonalWireAntenna::PolygonalWireAntenna( const ctor_key& key ) :
  18. WireAntenna( key ), minDipoleRatio(.15),
  19. minDipoleMoment(1e-6), maxDipoleMoment(1e1), rRepeat(1e10,1e10,1e10) {
  20. Points.setZero();
  21. //rRepeat.setOnes();
  22. }
  23. PolygonalWireAntenna::PolygonalWireAntenna( const YAML::Node& node, const ctor_key& key) : WireAntenna(node, key ) {
  24. minDipoleRatio = node["minDipoleRatio"].as<Real>();
  25. maxDipoleMoment = node["maxDipoleMoment"].as<Real>();
  26. minDipoleMoment = node["minDipoleMoment"].as<Real>();
  27. }
  28. PolygonalWireAntenna::~PolygonalWireAntenna() {
  29. }
  30. //--------------------------------------------------------------------------------------
  31. // Class: PolygonalWireAntenna
  32. // Method: Serialize
  33. //--------------------------------------------------------------------------------------
  34. YAML::Node PolygonalWireAntenna::Serialize ( ) const {
  35. YAML::Node node = WireAntenna::Serialize();
  36. node.SetTag( this->GetName() );
  37. node["minDipoleRatio"] = minDipoleRatio;
  38. node["maxDipoleMoment"] = maxDipoleMoment;
  39. node["minDipoleMoment"] = minDipoleMoment;
  40. return node;
  41. } // ----- end of method PolygonalWireAntenna::Serialize -----
  42. //--------------------------------------------------------------------------------------
  43. // Class: WireAntenna
  44. // Method: DeSerialize
  45. //--------------------------------------------------------------------------------------
  46. std::shared_ptr<PolygonalWireAntenna> PolygonalWireAntenna::DeSerialize ( const YAML::Node& node ) {
  47. if (node.Tag() != "PolygonalWireAntenna") {
  48. throw DeSerializeTypeMismatch( "PolygonalWireAntenna", node.Tag());
  49. }
  50. return std::make_shared<PolygonalWireAntenna> ( node, ctor_key() );
  51. } // ----- end of method WireAntenna::DeSerialize -----
  52. std::shared_ptr<PolygonalWireAntenna> PolygonalWireAntenna::NewSP() {
  53. return std::make_shared<PolygonalWireAntenna>( ctor_key() );
  54. }
  55. std::shared_ptr<WireAntenna> PolygonalWireAntenna::Clone() const {
  56. auto copy = PolygonalWireAntenna::NewSP();
  57. copy->minDipoleRatio = this->minDipoleRatio;
  58. copy->minDipoleMoment = this->minDipoleMoment;
  59. copy->maxDipoleMoment = this->maxDipoleMoment;
  60. copy->NumberOfPoints = this->NumberOfPoints;
  61. copy->Freqs = this->Freqs;
  62. copy->Current = this->Current;
  63. copy->NumberOfTurns = this->NumberOfTurns;
  64. copy->Points = this->Points;
  65. //copy->Dipoles = this->Dipoles; // no, disaster
  66. return copy;
  67. }
  68. std::shared_ptr<PolygonalWireAntenna> PolygonalWireAntenna::ClonePA() const {
  69. auto copy = PolygonalWireAntenna::NewSP();
  70. copy->minDipoleRatio = this->minDipoleRatio;
  71. copy->minDipoleMoment = this->minDipoleMoment;
  72. copy->maxDipoleMoment = this->maxDipoleMoment;
  73. copy->NumberOfPoints = this->NumberOfPoints;
  74. copy->Freqs = this->Freqs;
  75. copy->Current = this->Current;
  76. copy->NumberOfTurns = this->NumberOfTurns;
  77. copy->Points = this->Points;
  78. //copy->Dipoles = this->Dipoles; // no, disaster
  79. return copy;
  80. }
  81. //--------------------------------------------------------------------------------------
  82. // Class: PolygonalWireAntenna
  83. // Method: GetName
  84. // Description: Class identifier
  85. //--------------------------------------------------------------------------------------
  86. inline std::string PolygonalWireAntenna::GetName ( ) const {
  87. return CName;
  88. } // ----- end of method PolygonalWireAntenna::GetName -----
  89. void PolygonalWireAntenna::SetMinDipoleRatio (const Real& ratio) {
  90. minDipoleRatio = ratio;
  91. }
  92. void PolygonalWireAntenna::SetMinDipoleMoment (const Real& m) {
  93. minDipoleMoment = m;
  94. }
  95. void PolygonalWireAntenna::SetMaxDipoleMoment (const Real& m) {
  96. maxDipoleMoment = m;
  97. }
  98. // ==================== OPERATIONS =======================
  99. void PolygonalWireAntenna::ApproximateWithElectricDipoles(const Vector3r &rp) {
  100. // Only resplit if necessary. Save a few cycles if repeated
  101. if ( (rRepeat-rp).norm() > 1e-16 ) {
  102. Dipoles.clear();
  103. // loop over all segments
  104. for (int iseg=0; iseg<NumberOfPoints-1; ++iseg) {
  105. InterpolateLineSegment(Points.col(iseg), Points.col(iseg+1), rp);
  106. }
  107. rRepeat = rp;
  108. } else {
  109. for (unsigned int id=0; id<Dipoles.size(); ++id) {
  110. Dipoles[id]->SetFrequencies(Freqs);
  111. }
  112. }
  113. }
  114. Vector3r PolygonalWireAntenna::ClosestPointOnLine(const Vector3r &p1,
  115. const Vector3r &p2, const Vector3r &tp) {
  116. Vector3r v1 = p2 - p1;
  117. Vector3r v2 = p1 - tp;
  118. Vector3r v3 = p1 - p2;
  119. Vector3r v4 = p2 - tp;
  120. Real dot1 = v2.dot(v1);
  121. Real dot2 = v1.dot(v1);
  122. Real dot3 = v4.dot(v3);
  123. Real dot4 = v3.dot(v3);
  124. Real t1 = -1.*dot1/dot2;
  125. Real t2 = -1.*dot3/dot4;
  126. Vector3r pos = p1+v1*t1 ;
  127. // check if on line
  128. // else give back the closest end point
  129. if ( t1>=0 && t2>=0. ) {
  130. return pos;
  131. } else if (t1<0) {
  132. return p1;
  133. } else {
  134. return p2;
  135. }
  136. }
  137. void PolygonalWireAntenna::PushXYZDipoles(const Vector3r &step,
  138. const Vector3r &cp, const Vector3r &dir,
  139. std::vector< std::shared_ptr<DipoleSource> > &xDipoles) {
  140. Real scale = (Real)(NumberOfTurns)*Current;
  141. auto tx = DipoleSource::NewSP();
  142. tx->SetLocation(cp);
  143. tx->SetType(UNGROUNDEDELECTRICDIPOLE);
  144. tx->SetPolarisation(dir);
  145. tx->SetFrequencies(Freqs);
  146. tx->SetMoment(scale*step.norm());
  147. xDipoles.push_back(tx);
  148. }
  149. void PolygonalWireAntenna::CorrectOverstepXYZDipoles(const Vector3r &step,
  150. const Vector3r &cp, const Vector3r &dir,
  151. std::vector< std::shared_ptr<DipoleSource> > &xDipoles ) {
  152. Real scale = (Real)(NumberOfTurns)*Current;
  153. // X oriented dipoles
  154. if (step.norm() > minDipoleMoment) {
  155. xDipoles[xDipoles.size()-1]->SetLocation(cp);
  156. xDipoles[xDipoles.size()-1]->SetMoment(scale*step.norm());
  157. }
  158. }
  159. void PolygonalWireAntenna::InterpolateLineSegment(const Vector3r &p1,
  160. const Vector3r &p2, const Vector3r & tp) {
  161. Vector3r phat = (p1-p2).array() / (p1-p2).norm();
  162. Vector3r c = this->ClosestPointOnLine(p1, p2, tp);
  163. Real dtp = (tp-c).norm(); // distance to point at c
  164. Real dc1 = (p1-c).norm(); // distance to c from p1
  165. Real dc2 = (p2-c).norm(); // distance to c from p1
  166. // unit vector
  167. Vector3r cdir = (p2-p1).array() / (p2-p1).norm();
  168. ///////////////////
  169. // dipoles for this segment
  170. std::vector< std::shared_ptr<DipoleSource> > xDipoles;
  171. // go towards p1
  172. if ( ((c-p1).array().abs() > minDipoleMoment).any() ) {
  173. // cp = current pos, lp = last pos
  174. Vector3r cp = c + phat*(dtp*minDipoleRatio)*.5;
  175. Vector3r lp = c;
  176. Real dist = (cp-p1).norm();
  177. Real dist_old = dist+1.;
  178. // check that we didn't run past the end, or that we aren't starting at
  179. // the end, or that initial step runs over!
  180. Vector3r dir = (p1-cp).array() / (p1-cp).norm(); // direction of movement
  181. Vector3r step = phat*(dtp*minDipoleRatio);
  182. Vector3r stepold = Vector3r::Zero();
  183. // (dir-phat) just shows if we are stepping towards or away from p1
  184. while (dist < dist_old && (dir-phat).norm() < 1e-8) {
  185. PushXYZDipoles(step, cp, cdir, xDipoles);
  186. // Make 1/2 of previous step, 1/2 of this step, store this step
  187. stepold = step;
  188. step = phat*( (cp-tp).norm()*minDipoleRatio );
  189. while ( (step.array().abs() > maxDipoleMoment).any() ) {
  190. step *= .5;
  191. }
  192. lp = cp;
  193. cp += .5*stepold + .5*step;
  194. dist = (cp-p1).norm();
  195. dir = (p1-cp).array() / (p1-cp).norm();
  196. }
  197. // cp now points to end last of dipole moments
  198. cp -= .5*step;
  199. // Fix last dipole position, so that entire wire is represented,
  200. // and no more
  201. Real distLastSeg = (c - cp).norm();
  202. if (distLastSeg + minDipoleMoment < dc1) {
  203. // case 1: understep, add dipole
  204. step = (p1-cp).array();
  205. cp += .5*step;
  206. PushXYZDipoles(step, cp, cdir, xDipoles);
  207. } else if (distLastSeg > dc1 + minDipoleMoment) {
  208. // case 2: overstep, reposition dipole and size
  209. step = (p1 - (lp-.5*stepold));
  210. cp = (lp-.5*stepold) + (.5*step);
  211. CorrectOverstepXYZDipoles(step, cp, cdir, xDipoles);
  212. }
  213. // else case 0: nearly 'perfect' fit do nothing
  214. }
  215. // go towards p2
  216. if ( ( (c-p2).array().abs() > minDipoleMoment).any() ) {
  217. // cp = current pos, lp = last pos
  218. Vector3r step = -phat*(dtp*minDipoleRatio);
  219. while ( (step.array().abs() > maxDipoleMoment).any() ) {
  220. step *= .5;
  221. }
  222. Vector3r cp = c + step*.5;
  223. Vector3r lp = c;
  224. Real dist = (p2-cp).norm();
  225. Real dist_old = dist+1e3;
  226. // check that we didn't run past the end, or that we aren't starting at
  227. // the end, or that initial step runs over!
  228. Vector3r dir = (p2-cp).array() / (p2-cp).norm(); // direction of movement
  229. Vector3r stepold = Vector3r::Zero();
  230. // (dir-phat) just shows if we are stepping towards or away from p1
  231. while (dist < dist_old && (dir+phat).norm() < 1e-8) {
  232. PushXYZDipoles(step, cp, cdir, xDipoles);
  233. // Make 1/2 of previous step, 1/2 of this step, store this step
  234. stepold = step;
  235. step = -phat*( (cp-tp).norm()*minDipoleRatio );
  236. while ( (step.array().abs() > maxDipoleMoment).any() ) {
  237. step *= .5;
  238. }
  239. lp = cp;
  240. cp += .5*stepold + .5*step;
  241. dist = (cp-p2).norm();
  242. dir = (p2-cp).array() / (p2-cp).norm();
  243. }
  244. // cp now points to end last of dipole moments
  245. cp -= .5*step;
  246. // Fix last dipole position, so that entire wire is represented,
  247. // and no more
  248. Real distLastSeg = (c - cp).norm();
  249. if (distLastSeg + minDipoleMoment < dc2) {
  250. // case 1: understep, add dipole
  251. step = (p2-cp).array();
  252. cp += .5*step;
  253. PushXYZDipoles(step, cp, cdir, xDipoles);
  254. } else if (distLastSeg > dc2 + minDipoleMoment) {
  255. // case 2: overstep, reposition dipole and size
  256. step = (p2 - (lp-.5*stepold));
  257. cp = (lp-.5*stepold) + (.5*step);
  258. CorrectOverstepXYZDipoles(step, cp, cdir, xDipoles);
  259. }
  260. // else case 0: nearly 'perfect' fit do nothing
  261. }
  262. Dipoles.insert(Dipoles.end(), xDipoles.begin(), xDipoles.end());
  263. }
  264. }