MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Panel.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 #include <cstddef>
26 #include <cmath>
27 #include <vector>
28 #include "Dipole.h"
29 #include "Panel.h"
30 #include "Math/NumericalVector.h"
31 #include "Math/Vector.h"
32 #include "Multipole.h"
33 #include "Source.h"
34 #include "Sphere.h"
35 #include "Box.h"
36 
38  int panelLevel,
39  Vec3D leftBound,
40  Vec3D rightBound,
41  std::vector<Source*> sources,
42  std::vector<Dipole*> dipoles,
43  NumericalVector<>* squaredFactorials,
44  Box* box) :
45  panelLevel_(panelLevel),
46  leftBound_(leftBound),
47  rightBound_(rightBound),
48  root_(root),
49  sources_(sources),
50  dipoles_(dipoles),
51  box_(box)
52 {
53  //Initialise the panel characteristics
54  // I am using Vec3D, so the size is always constant
55 /* if (leftBound.size() != rightBound.size())
56  {
57  std::cout << "Bounds are not of correct dimensions" << std::endl;
58  std::exit(-1);
59  }*/
60 
61  // This code only works for 3D problems, since the maths is significantly different from 2D and 1D
62  // The datastructure however allows for 1D and 2D structures,
63  // In future 2D support might be added (by a master student?)
64  dim_ = 3;
65 
66  // The panel is a cube with same size in all directions
67  // Note this is half the length of a cube side
68  size_ = 0.5 * (rightBound.getComponent(0) - leftBound.getComponent(0));
69 
70  // Determine the centre
71  for (int iD = 0; iD < dim_; iD++)
72  {
73  centre_.setComponent(iD, (leftBound.getComponent(iD) + size_));
74  }
75 
76  //Add panel to the complete box
77  box->addPanel(panelLevel, this);
78 
79  //Set multipole and local expansion coefficients
80  multipoleAroundCentre_ = new Multipole(box_->getNumberOfTerms(), squaredFactorials, centre_);
81  multipoleAroundCentre_->computeMultipoleExpansion(); // Set vector with zeros to correct length
83  partialLocalExpansionAroundCentre_->initialiseLocalExpansion(); // Set vector with zeros to correct length
85  localExpansionAroundCentre_->initialiseLocalExpansion(); // Set vector with zeros to correct length
86 
87  //Create childeren based on the given dimensions;
88  if (panelLevel < box->getMaxLevel()) // final level is obtained if this value is 0
89  {
90  //Note: For an adaptive grid, add a check to see if this panel contains enough sources to have childeren.
91 
92  Vec3D leftBoundChild;
93  Vec3D rightBoundChild;
94 
95  createPanels(1, sources, dipoles, leftBoundChild, rightBoundChild, squaredFactorials);
96  }
97 
98 }
99 
100 void Panel::createPanels(int dim,
101  std::vector<Source*>& sources,
102  std::vector<Dipole*>& dipoles,
103  Vec3D& leftBoundChild,
104  Vec3D& rightBoundChild,
105  NumericalVector<>* squaredFactorials)
106 {
107  std::vector<Source*> sourcesChild;
108  std::vector<Dipole*> dipolesChild;
109 
110 
111  for (int iSide = 0; iSide < 2; ++iSide)
112  {
113  //Determine the x boundaries
114  leftBoundChild.setComponent(dim - 1, (leftBound_.getComponent(dim - 1) + size_ * iSide));
115  rightBoundChild.setComponent(dim - 1, (leftBoundChild.getComponent(dim - 1) + size_));
116 
117  //Determine the sources in this domain
118  for (Source* iS : sources)
119  {
120  if ((iS->getLocation().getComponent(dim - 1) > leftBoundChild.getComponent(dim - 1)) &&
121  (iS->getLocation().getComponent(dim - 1) <= rightBoundChild.getComponent(dim - 1)))
122  {
123  sourcesChild.push_back(iS);
124  }
125  }
126 
127  //Determine the dipoles in this domain
128  for (Dipole* iD : dipoles)
129  {
130  if ((iD->getLocation().getComponent(dim - 1) > leftBoundChild.getComponent(dim - 1)) &&
131  (iD->getLocation().getComponent(dim - 1) <= rightBoundChild.getComponent(dim - 1)))
132  {
133  dipolesChild.push_back(iD);
134  }
135  }
136 
137  if (dim != dim_)
138  {
139  createPanels(dim + 1, sourcesChild, dipolesChild, leftBoundChild, rightBoundChild, squaredFactorials);
140  }
141  else
142  {
143  //Create childeren
144 /* std::cout << "Creating a child with bounds: " << std::endl;
145  std::cout << "In x-direction: " << std::endl;
146  std::cout << "Left bound: " << leftBoundChild[0] << " and right bound: " << rightBoundChild[0] << std::endl;
147  std::cout << "In y-direction: " << std::endl;
148  std::cout << "Left bound: " << leftBoundChild[1] << " and right bound: " << rightBoundChild[1] << std::endl;
149  std::cout << "In z-direction: " << std::endl;
150  std::cout << "Left bound: " << leftBoundChild[2] << " and right bound: " << rightBoundChild[2] << std::endl;
151  std::cout << "Child contains " << sourcesChild.size() << " sources" << std::endl;
152  std::cout << "==============" << std::endl;*/
153  Panel* child = new Panel(this, panelLevel_ + 1, leftBoundChild, rightBoundChild, sourcesChild, dipolesChild,
154  squaredFactorials, box_);
155  childeren_.push_back(child);
156 
157  //reset leftBoundChild and rightBoundChild
158  }
159  sourcesChild.clear();
160  }
161 }
162 
164 {
165  //Create neighbour list and interaction list
167 
168 }
169 
171 {
172  //Start upward pass creating multipole methods and shifting them on parent panels
173  box_->upwardPass();
174 
175  //Start the downward pass
176  box_->downwardPass();
177 }
178 
180 {
181  if (getRoot() == nullptr)
182  {
183  //nothing happens
184  }
185  else
186  {
187  //Crate the links between panels
188  this->setPanelInteractions();
189  }
190 
191  //Initialise child panels
192  for (Panel* iP : childeren_)
193  {
194  iP->findPanelInteractions();
195  }
196 }
197 
199 {
200  //Childeren of a parent are always a neighbour of eachother
201  //Note: by definition, a panel is always a neighbour of itself
202  for (Panel* iC : root_->getChilderen())
203  {
204  neighbours_.push_back(iC);
205  }
206 
207  //For all the neighbours of the parent of this panel
208  for (Panel* iN : root_->getNeighbours())
209  {
210  //For all their childeren, check if the childeren are a neighbour or a second neighbour of this panel
211  // Note that all second neighbour panels are also childeren of these paretn neighbours
212  for (Panel* iC : iN->getChilderen())
213  {
214  //Compute distance between the centre of the current panel and possible neighbour panel
215  Mdouble distance = Vec3D::getDistance(iC->getCentre(), this->getCentre());
216 
217  //Check if the possible neighbour panel is close enough
218  if (distance < 3.0 * size_) //warning: slight hack. I assume all sides are equally long
219  {
220  neighbours_.push_back(iC);
221  //std::cout << "A neighbour has been created with a centre around: (" << iC->getCentre()[0] << "," << iC->getCentre()[1] << ")" << std::endl;
222  //std::cout << "Current panel has centre: (" << this->getCentre()[0] << "," << this->getCentre()[1] << ")" << std::endl;
223  //std::cout << "===============" << std::endl;
224  }
225  else if (distance < 6.0 * size_)
226  {
227  secondNeighbours_.push_back(iC);
228  }
229  else
230  {
231  // If it is not a neighbour or nearest neighbour, it belongs to the interactionlist group
232  interactionList_.push_back(iC);
233  //std::cout << "A panel is added to the interactionlist with a centre around: (" << iC->getCentre()[0] << "," << iC->getCentre()[1] << ")" << std::endl;
234  //std::cout << "Current panel has centre: (" << this->getCentre()[0] << "," << this->getCentre()[1] << ")" << std::endl;
235  //std::cout << "===============" << std::endl;
236  }
237  }
238  }
239 
240  // All childeren of the second nearest neighbour of the current panel's parent, are in the interactionlist
241  for (Panel* iN : root_->getSecondNeighbours())
242  {
243  for (Panel* iC : iN->getChilderen())
244  {
245  secondNeighbours_.push_back(iC);
246  }
247  }
248 }
249 
251 {
252 
254  // todo: write this part of the code
255  // A source can be expanded around the centre of the panel without translations
256  //add(iS->computeMultipoleExpansion();)
257 
259  for (Dipole* iD : dipoles_)
260  {
261  //Convert to a multipole
262  //todo: This must be done when the actual dipole is made.
263  iD->computeMultipoleExpansion();
264 
265  //Translate multipole to centre of the panel
266  multipoleAroundCentre_->addMultipoleCoefficients(iD->TranslateMultipoleExpansionTo(centre_));
267 
268  }
269 
271  for (Multipole* iM : multipoles_)
272  {
273  multipoleAroundCentre_->addMultipoleCoefficients(iM->TranslateMultipoleExpansionTo(centre_));
274  }
275 
276 }
277 
279 {
280  // Extract the multipole of all the childeren
281  for (Panel* iC : childeren_)
282  {
283  //Translate the multipole in the center of the child to the centre of the current panel
284  //And then add it to the multipole in the centre of the current panel
286  iC->multipoleAroundCentre_->TranslateMultipoleExpansionTo(centre_)
287  );
288  }
289 }
290 
292 {
293  //Note: nothing has to be done here, the current implementation initialises all local expansions with 0.
294 }
295 
297 {
298  //For all panels in the interactionlist, add a local expansion of the multipole to the partialLocalExpansion
299  for (Panel* panel : interactionList_)
300  {
302  panel->multipoleAroundCentre_->convertMultipoleToLocal(centre_)
303  );
304  }
305 }
306 
308 {
311  );
312 }
313 
315 {
316  for (Panel* iC : childeren_)
317  {
318  Vec3D centreChild = iC->getCentre();
319  iC->localExpansionAroundCentre_->addLocalExpansionCoefficients(
321  );
322  }
323 }
324 
325 
326 
327 
328 
329 
330 
331 
332 
333 
std::vector< Multipole * > multipoles_
Definition: Panel.h:172
Definition: Panel.h:42
void translateMultipoleExpansion()
Definition: Panel.cc:278
Panel * getRoot()
Definition: Panel.h:88
virtual void computeMultipoleExpansion()
Definition: Multipole.cc:43
void findPanelInteractions()
Definition: Panel.cc:179
Panel * root_
Definition: Panel.h:163
void createPanels(int dim, std::vector< Source * > &sources, std::vector< Dipole * > &dipoles, Vec3D &leftBoundChild, Vec3D &rightBoundChild, NumericalVector<> *squaredFactorials)
Definition: Panel.cc:100
double Mdouble
Definition: GeneralDefine.h:34
Multipole * multipoleAroundCentre_
Definition: Panel.h:182
std::vector< Panel * > interactionList_
Definition: Panel.h:167
Mdouble dim_
Definition: Panel.h:156
void addMultipoleCoefficients(NumericalVector< std::complex< Mdouble >> multipoleExpansionCoefficients)
Adds multipole coefficients to an existing multipole todo: remove this function; it should not be req...
Definition: Multipole.cc:153
std::vector< Panel * > getChilderen()
Definition: Panel.h:93
Definition: Dipole.h:34
std::vector< Panel * > secondNeighbours_
Definition: Panel.h:166
LocalExpansion * localExpansionAroundCentre_
Definition: Panel.h:184
void computeMultipoleExpansion()
Definition: Panel.cc:250
void addLocalExpansionCoefficients(NumericalVector< std::complex< Mdouble >> localExpansionCoefficients)
void initialise()
Definition: Panel.cc:163
void computePartialLocalExpansion()
Definition: Panel.cc:296
void translateLocalExpansion()
Definition: Panel.cc:314
void upwardPass()
Definition: Box.cc:44
void setComponent(int index, double val)
Sets the requested component of this Vec3D to the requested value.
Definition: Vector.cc:217
const int panelLevel_
Definition: Panel.h:155
Mdouble getComponent(int index) const
Returns the requested component of this Vec3D.
Definition: Vector.cc:194
static Mdouble getDistance(const Vec3D &a, const Vec3D &b)
Calculates the distance between two Vec3D: .
Definition: Vector.cc:175
Panel(Panel *root, int maximumPanelLevel, Vec3D leftBound, Vec3D rightBound, std::vector< Source * > sources, std::vector< Dipole * > dipoles, NumericalVector<> *squaredFactorials, Box *box)
Definition: Panel.cc:37
LocalExpansion * partialLocalExpansionAroundCentre_
Definition: Panel.h:183
void setPanelInteractions()
Definition: Panel.cc:198
void computeCoefficients()
Definition: Panel.cc:170
Vec3D getCentre()
Definition: Panel.h:83
Definition: Box.h:31
std::vector< Panel * > childeren_
Definition: Panel.h:164
std::vector< Panel * > neighbours_
Definition: Panel.h:165
int getNumberOfTerms()
Definition: Box.h:55
void addPanel(int level, Panel *panel)
Definition: Box.cc:39
std::vector< Panel * > getSecondNeighbours()
Definition: Panel.h:103
Vec3D leftBound_
Definition: Panel.h:157
Vec3D centre_
Definition: Panel.h:160
void computeLocalExpansion()
Definition: Panel.cc:307
NumericalVector< std::complex< Mdouble > > getExpansionCoefficients()
void initialiseLocalExpansion()
Box * box_
Definition: Panel.h:178
void setLocalExpansionZero()
Definition: Panel.cc:291
NumericalVector< std::complex< Mdouble > > translateLocalExpansion(Vec3D location)
Definition: Vector.h:49
void downwardPass()
Definition: Box.cc:68
double size_
Definition: Panel.h:159
Definition: Source.h:33
std::vector< Dipole * > dipoles_
Definition: Panel.h:171
std::vector< Panel * > getNeighbours()
Definition: Panel.h:98