Main Lemma Repository
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

PolygonalWireAntenna.cpp 11KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338
  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. ///////////////////
  104. // dipole array, this has impoved performance over directly pushing
  105. std::vector< std::shared_ptr<DipoleSource> > xDipoles;
  106. // loop over all segments
  107. int iseg;
  108. for (iseg=0; iseg<NumberOfPoints-1; ++iseg) {
  109. InterpolateLineSegment(Points.col(iseg), Points.col(iseg+1), rp, xDipoles);
  110. }
  111. // Check to see if the loop is closed, if not, assume its grounded on ends,
  112. if ( (Points.col(0)-Points.col(iseg)).norm() > 1e-3 ) {
  113. xDipoles[0]->SetType(GROUNDEDELECTRICDIPOLE);
  114. xDipoles.back()->SetType(GROUNDEDELECTRICDIPOLE);
  115. }
  116. Dipoles = std::move(xDipoles);
  117. rRepeat = rp;
  118. } else {
  119. for (unsigned int id=0; id<Dipoles.size(); ++id) {
  120. Dipoles[id]->SetFrequencies(Freqs);
  121. }
  122. }
  123. }
  124. Vector3r PolygonalWireAntenna::ClosestPointOnLine(const Vector3r &p1,
  125. const Vector3r &p2, const Vector3r &tp) {
  126. Vector3r v1 = p2 - p1;
  127. Vector3r v2 = p1 - tp;
  128. Vector3r v3 = p1 - p2;
  129. Vector3r v4 = p2 - tp;
  130. Real dot1 = v2.dot(v1);
  131. Real dot2 = v1.dot(v1);
  132. Real dot3 = v4.dot(v3);
  133. Real dot4 = v3.dot(v3);
  134. Real t1 = -1.*dot1/dot2;
  135. Real t2 = -1.*dot3/dot4;
  136. Vector3r pos = p1+v1*t1 ;
  137. // check if on line
  138. // else give back the closest end point
  139. if ( t1>=0 && t2>=0. ) {
  140. return pos;
  141. } else if (t1<0) {
  142. return p1;
  143. } else {
  144. return p2;
  145. }
  146. }
  147. void PolygonalWireAntenna::PushXYZDipoles(const Vector3r &step,
  148. const Vector3r &cp, const Vector3r &dir,
  149. std::vector< std::shared_ptr<DipoleSource> > &xDipoles) {
  150. Real scale = (Real)(NumberOfTurns)*Current;
  151. auto tx = DipoleSource::NewSP();
  152. tx->SetLocation(cp);
  153. tx->SetType(UNGROUNDEDELECTRICDIPOLE);
  154. tx->SetPolarisation(dir);
  155. tx->SetFrequencies(Freqs);
  156. tx->SetMoment(scale*step.norm());
  157. xDipoles.push_back(tx);
  158. }
  159. void PolygonalWireAntenna::CorrectOverstepXYZDipoles(const Vector3r &step,
  160. const Vector3r &cp, const Vector3r &dir,
  161. std::vector< std::shared_ptr<DipoleSource> > &xDipoles ) {
  162. Real scale = (Real)(NumberOfTurns)*Current;
  163. // X oriented dipoles
  164. if (step.norm() > minDipoleMoment) {
  165. xDipoles[xDipoles.size()-1]->SetLocation(cp);
  166. xDipoles[xDipoles.size()-1]->SetMoment(scale*step.norm());
  167. }
  168. }
  169. void PolygonalWireAntenna::InterpolateLineSegment(const Vector3r &p1,
  170. const Vector3r &p2, const Vector3r & tp,
  171. std::vector< std::shared_ptr<DipoleSource> > &xDipoles ) {
  172. Vector3r phat = (p1-p2).array() / (p1-p2).norm();
  173. Vector3r c = this->ClosestPointOnLine(p1, p2, tp);
  174. Real dtp = (tp-c).norm(); // distance to point at c
  175. Real dc1 = (p1-c).norm(); // distance to c from p1
  176. Real dc2 = (p2-c).norm(); // distance to c from p1
  177. // unit vector
  178. Vector3r cdir = (p2-p1).array() / (p2-p1).norm();
  179. // go towards p1
  180. if ( ((c-p1).array().abs() > minDipoleMoment).any() ) {
  181. // cp = current pos, lp = last pos
  182. Vector3r cp = c + phat*(dtp*minDipoleRatio)*.5;
  183. Vector3r lp = c;
  184. Real dist = (cp-p1).norm();
  185. Real dist_old = dist+1.;
  186. // check that we didn't run past the end, or that we aren't starting at
  187. // the end, or that initial step runs over!
  188. Vector3r dir = (p1-cp).array() / (p1-cp).norm(); // direction of movement
  189. Vector3r step = phat*(dtp*minDipoleRatio);
  190. Vector3r stepold = Vector3r::Zero();
  191. // (dir-phat) just shows if we are stepping towards or away from p1
  192. while (dist < dist_old && (dir-phat).norm() < 1e-8) {
  193. PushXYZDipoles(step, cp, cdir, xDipoles);
  194. // Make 1/2 of previous step, 1/2 of this step, store this step
  195. stepold = step;
  196. step = phat*( (cp-tp).norm()*minDipoleRatio );
  197. while ( (step.array().abs() > maxDipoleMoment).any() ) {
  198. step *= .5;
  199. }
  200. lp = cp;
  201. cp += .5*stepold + .5*step;
  202. dist = (cp-p1).norm();
  203. dir = (p1-cp).array() / (p1-cp).norm();
  204. }
  205. // cp now points to end last of dipole moments
  206. cp -= .5*step;
  207. // Fix last dipole position, so that entire wire is represented,
  208. // and no more
  209. Real distLastSeg = (c - cp).norm();
  210. if (distLastSeg + minDipoleMoment < dc1) {
  211. // case 1: understep, add dipole
  212. step = (p1-cp).array();
  213. cp += .5*step;
  214. PushXYZDipoles(step, cp, cdir, xDipoles);
  215. } else if (distLastSeg > dc1 + minDipoleMoment) {
  216. // case 2: overstep, reposition dipole and size
  217. step = (p1 - (lp-.5*stepold));
  218. cp = (lp-.5*stepold) + (.5*step);
  219. CorrectOverstepXYZDipoles(step, cp, cdir, xDipoles);
  220. }
  221. // else case 0: nearly 'perfect' fit do nothing
  222. }
  223. // go towards p2
  224. if ( ( (c-p2).array().abs() > minDipoleMoment).any() ) {
  225. // cp = current pos, lp = last pos
  226. Vector3r step = -phat*(dtp*minDipoleRatio);
  227. while ( (step.array().abs() > maxDipoleMoment).any() ) {
  228. step *= .5;
  229. }
  230. Vector3r cp = c + step*.5;
  231. Vector3r lp = c;
  232. Real dist = (p2-cp).norm();
  233. Real dist_old = dist+1e3;
  234. // check that we didn't run past the end, or that we aren't starting at
  235. // the end, or that initial step runs over!
  236. Vector3r dir = (p2-cp).array() / (p2-cp).norm(); // direction of movement
  237. Vector3r stepold = Vector3r::Zero();
  238. // (dir-phat) just shows if we are stepping towards or away from p1
  239. while (dist < dist_old && (dir+phat).norm() < 1e-8) {
  240. PushXYZDipoles(step, cp, cdir, xDipoles);
  241. // Make 1/2 of previous step, 1/2 of this step, store this step
  242. stepold = step;
  243. step = -phat*( (cp-tp).norm()*minDipoleRatio );
  244. while ( (step.array().abs() > maxDipoleMoment).any() ) {
  245. step *= .5;
  246. }
  247. lp = cp;
  248. cp += .5*stepold + .5*step;
  249. dist = (cp-p2).norm();
  250. dir = (p2-cp).array() / (p2-cp).norm();
  251. }
  252. // cp now points to end last of dipole moments
  253. cp -= .5*step;
  254. // Fix last dipole position, so that entire wire is represented,
  255. // and no more
  256. Real distLastSeg = (c - cp).norm();
  257. if (distLastSeg + minDipoleMoment < dc2) {
  258. // case 1: understep, add dipole
  259. step = (p2-cp).array();
  260. cp += .5*step;
  261. PushXYZDipoles(step, cp, cdir, xDipoles);
  262. } else if (distLastSeg > dc2 + minDipoleMoment) {
  263. // case 2: overstep, reposition dipole and size
  264. step = (p2 - (lp-.5*stepold));
  265. cp = (lp-.5*stepold) + (.5*step);
  266. CorrectOverstepXYZDipoles(step, cp, cdir, xDipoles);
  267. }
  268. // else case 0: nearly 'perfect' fit do nothing
  269. }
  270. //Dipoles.insert(Dipoles.end(), xDipoles.begin(), xDipoles.end());
  271. //Dipoles = xDipoles;
  272. //xDipoles.clear();
  273. }
  274. }