← Back to branch summary

~ubuntu-branches/ubuntu/quantal/python-demgengeo/quantal

« back to all changes in this revision

Viewing changes to src/CircMNTable2D.cc

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2011-11-18 21:47:18 UTC
  • Revision ID: package-import@ubuntu.com-20111118214718-4ysqm3dhpqwdd7gd
Tags: upstream-0.99~bzr106
Import upstream version 0.99~bzr106

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/CircMNTable2D.cc
 
1
/////////////////////////////////////////////////////////////
 
2
//                                                         //
 
3
// Copyright (c) 2007-2011 by The University of Queensland //
 
4
// Earth Systems Science Computational Centre (ESSCC)      //
 
5
// http://www.uq.edu.au/esscc                              //
 
6
//                                                         //
 
7
// Primary Business: Brisbane, Queensland, Australia       //
 
8
// Licensed under the Open Software License version 3.0    //
 
9
// http://www.opensource.org/licenses/osl-3.0.php          //
 
10
//                                                         //
 
11
/////////////////////////////////////////////////////////////
 
12
 
 
13
#include "CircMNTable2D.h"
 
14
 
 
15
// --- System includes ---
 
16
#include <cmath>
 
17
 
 
18
// --- IO includes ---
 
19
#include <iostream>
 
20
 
 
21
using std::endl;
 
22
 
 
23
using std::floor;
 
24
 
 
25
CircMNTable2D::CircMNTable2D()
 
26
{}
 
27
 
 
28
/*!
 
29
  Construct CircMNTable2D. 
 
30
 
 
31
  \param MinPt minimum point (z component ignored)
 
32
  \param MaxPt maximum point (z component ignored)
 
33
  \param cd cell dimension
 
34
  \param ngroups initial number of particle groups
 
35
*/
 
36
CircMNTable2D::CircMNTable2D(const Vector3& MinPt,const Vector3& MaxPt,double cd,unsigned int ngroups):
 
37
  MNTable2D(MinPt, MaxPt, cd, ngroups)
 
38
{
 
39
  set_x_circ();
 
40
  // check if grid spacing fits size in circular direction:
 
41
  double nx=(MaxPt-MinPt).X()/m_celldim;
 
42
  // error message if not
 
43
  if(nx!=floor(nx)){
 
44
    std::cerr << "WARNING! grid spacing " << m_celldim << " doesn't fit periodic x-dimension " << (MaxPt-MinPt).X() << std::endl;
 
45
  }
 
46
  double shift_x=(m_nx-2)*m_celldim;
 
47
  m_shift_vec_x=Vector3(shift_x,0.0,0.0);
 
48
}
 
49
 
 
50
/*!
 
51
  Destruct CircMNTable2D. Just calls MNTable2D destructor.
 
52
*/
 
53
CircMNTable2D::~CircMNTable2D()
 
54
{}
 
55
 
 
56
/*!
 
57
  set circularity of x-dimension to 1
 
58
*/
 
59
void CircMNTable2D::set_x_circ()
 
60
{
 
61
  m_x_periodic=1;
 
62
} 
 
63
 
 
64
int CircMNTable2D::getXIndex(const Vector3& Pos) const
 
65
{
 
66
  return int(floor((Pos.x()-m_x0)/m_celldim));
 
67
}
 
68
 
 
69
int CircMNTable2D::getYIndex(const Vector3& Pos) const
 
70
{
 
71
  return int(floor((Pos.y()-m_y0)/m_celldim));
 
72
}
 
73
 
 
74
 
 
75
int CircMNTable2D::getFullIndex(const Vector3& Pos) const
 
76
{
 
77
  int ix=int(floor((Pos.x()-m_x0)/m_celldim));
 
78
  int iy=int(floor((Pos.y()-m_y0)/m_celldim));
 
79
 
 
80
  return idx(ix,iy);
 
81
}
 
82
 
 
83
/*!
 
84
  get the cell index for a given position
 
85
 
 
86
  \param Pos the position
 
87
  \return the cell index if Pos is inside the table, -1 otherwise
 
88
*/
 
89
int CircMNTable2D::getIndex(const Vector3& Pos) const
 
90
{
 
91
  int ret;
 
92
 
 
93
  int ix=int(floor((Pos.x()-m_x0)/m_celldim));
 
94
  int iy=int(floor((Pos.y()-m_y0)/m_celldim));
 
95
 
 
96
  // check if pos is in table excl. padding
 
97
  if((ix>=0) && (ix<=m_nx-1) && (iy>0) && (iy<m_ny-1)){
 
98
    ret=idx(ix,iy);
 
99
  } else {
 
100
    ret=-1;
 
101
  }
 
102
  
 
103
  return ret;
 
104
}
 
105
 
 
106
/*!
 
107
  Insert sphere. Insert clone into other side of Table.
 
108
 
 
109
  \param S the Sphere
 
110
  \param gid the group id
 
111
*/
 
112
bool CircMNTable2D::insert(const Sphere& S,unsigned int gid)
 
113
{
 
114
  bool res;
 
115
 
 
116
  int id=getIndex(S.Center());
 
117
  int idx=getXIndex(S.Center());
 
118
 
 
119
  if((id!=-1) && (idx!=0) && (idx!=m_nx-1) && (gid<m_ngroups)){ // valid index
 
120
    // insert sphere
 
121
    m_data[id].insert(S,gid);
 
122
    res=true;
 
123
    int xidx=getXIndex(S.Center());
 
124
    // insert clone
 
125
    if (xidx==1){
 
126
      Sphere SClone=S;
 
127
      SClone.shift(m_shift_vec_x);
 
128
      int clone_id=getFullIndex(SClone.Center());
 
129
      m_data[clone_id].insert(SClone,gid);
 
130
    } else if  (xidx==m_nx-2){
 
131
      Sphere SClone=S;
 
132
      SClone.shift(-1.0*m_shift_vec_x);
 
133
      int clone_id=getFullIndex(SClone.Center());
 
134
      m_data[clone_id].insert(SClone,gid);
 
135
    } 
 
136
  } else {
 
137
    res=false;
 
138
  }
 
139
 
 
140
  return res;
 
141
}
 
142
 
 
143
/*!
 
144
  Insert sphere if it doesn't collide with other spheres. Insert clone into other side of Table.
 
145
 
 
146
  \param S the Sphere
 
147
  \param gid the group id
 
148
*/
 
149
bool CircMNTable2D::insertChecked(const Sphere& S,unsigned int gid,double tol)
 
150
{
 
151
  bool res;
 
152
  int id=getIndex(S.Center());
 
153
  int idx=getXIndex(S.Center());
 
154
  
 
155
//   std::cerr << "CircMNTable2D::insertChecked(" << S << ") id=" << id << " idx= " << idx << std::endl;
 
156
 
 
157
  tol+=s_small_value;
 
158
  if((id!=-1) && (idx!=0) && (idx!=m_nx-1) && (gid<m_ngroups)){
 
159
    // insert original
 
160
    multimap<double,const Sphere*> close_spheres=getSpheresFromGroupNear(S.Center(),S.Radius()-tol,gid);
 
161
    if(close_spheres.size()==0){
 
162
      m_data[id].insert(S,gid);
 
163
      res=true;
 
164
    } else res=false;
 
165
    int xidx=getXIndex(S.Center());
 
166
    // insert clone
 
167
    if (xidx==1){
 
168
      Sphere SClone=S;
 
169
      SClone.shift(m_shift_vec_x);
 
170
      multimap<double,const Sphere*> close_spheres=getSpheresFromGroupNear(SClone.Center(),SClone.Radius()-tol,gid);
 
171
      if(close_spheres.size()==0){
 
172
        int clone_id=getFullIndex(SClone.Center());
 
173
        m_data[clone_id].insert(SClone,gid);
 
174
      }
 
175
    } else if  (xidx==m_nx-2){
 
176
      Sphere SClone=S;
 
177
      SClone.shift(-1.0*m_shift_vec_x);
 
178
      multimap<double,const Sphere*> close_spheres=getSpheresFromGroupNear(SClone.Center(),SClone.Radius()-tol,gid);
 
179
      if(close_spheres.size()==0){
 
180
        int clone_id=getFullIndex(SClone.Center());
 
181
        m_data[clone_id].insert(SClone,gid);
 
182
      }
 
183
    } 
 
184
  } else {
 
185
    res=false;
 
186
  }
 
187
 
 
188
  return res;
 
189
}
 
190
 
 
191
/*!
 
192
  insert sphere if it doesn't collide with other spheres
 
193
 
 
194
  \param S the Sphere
 
195
  \param gid the group id
 
196
*/
 
197
bool CircMNTable2D::checkInsertable(const Sphere& S,unsigned int gid)
 
198
{
 
199
  bool res;
 
200
 
 
201
  int id=getIndex(S.Center());
 
202
  int idx=getXIndex(S.Center());
 
203
  
 
204
 if((id!=-1) && (idx!=0) && (idx!=m_nx-1) && (gid<m_ngroups)){ 
 
205
    multimap<double,const Sphere*> close_spheres=getSpheresFromGroupNear(S.Center(),S.Radius()-s_small_value,gid);
 
206
    if(close_spheres.size()==0){ 
 
207
      res=true;
 
208
    } else res=false;
 
209
  } else {
 
210
    res=false;
 
211
  }
 
212
 
 
213
  return res;
 
214
}
 
215
 
 
216
 
 
217
/*!
 
218
  Generate bonds between particles of a group. Takes cloned particles into account.
 
219
 
 
220
  \param gid the group ID
 
221
  \param tol max. difference between bond length and equilibrium dist.
 
222
  \param btag bond tag
 
223
*/
 
224
void CircMNTable2D::generateBonds(int gid,double tol,int btag)
 
225
{
 
226
  std::cerr << "CircMNTable2D::generateBonds( " << gid << " , " << tol << " , " << btag << " )" << std::endl;
 
227
  // loop over all cells 
 
228
  for(int i=0;i<m_nx-1;i++){
 
229
    for(int j=1;j<m_ny-1;j++){
 
230
      int id=idx(i,j);
 
231
      // loop over "upper" neighbors of each cell
 
232
      for(int ii=-1;ii<=1;ii++){
 
233
        for(int jj=-1;jj<=1;jj++){
 
234
          int id2=idx(i+ii,j+jj);
 
235
          vector<pair<int,int> > bonds;
 
236
          if((id2==id) && (i!=0)){ // intra-cell, not for boundary 
 
237
            //std::cerr << id << " - " << id << std::endl;
 
238
            bonds=m_data[id].getBonds(gid,tol);
 
239
          } else if(id2>id){ // inter-cell
 
240
            //std::cerr << id << " - " << id2 << std::endl;
 
241
            bonds=m_data[id].getBonds(gid,tol,m_data[id2]);
 
242
          }
 
243
          for(vector<pair<int,int> >::iterator iter=bonds.begin();
 
244
              iter!=bonds.end();
 
245
              iter++){
 
246
            //std::cerr << iter->first << " | " << iter->second << "   "; 
 
247
            m_bonds[btag].insert(*iter);
 
248
          }
 
249
          //std::cerr << std::endl;
 
250
        }
 
251
      }
 
252
      //std::cerr << std::endl;
 
253
    }
 
254
  }
 
255
}
 
256
 
 
257
/*!
 
258
  Generate bonds between particles of a group with identical particle tag. Takes cloned particles into account.
 
259
 
 
260
  \param gid the group ID
 
261
  \param tol max. difference between bond length and equilibrium dist.
 
262
  \param btag bond tag
 
263
  \param ptag the particle tag
 
264
  \param mask the mask for the particle tag
 
265
*/
 
266
void CircMNTable2D::generateBondsWithMask(int gid,double tol,int btag, int ptag, int mask)
 
267
{
 
268
  std::cerr << "CircMNTable2D::generateBondsWithMask( " << gid << " , " << tol << " , " << btag << " " 
 
269
            << ptag <<" " << mask<<" )" << std::endl;
 
270
  // loop over all cells 
 
271
  for(int i=0;i<m_nx-1;i++){
 
272
    for(int j=1;j<m_ny-1;j++){
 
273
      int id=idx(i,j);
 
274
      // loop over "upper" neighbors of each cell
 
275
      for(int ii=-1;ii<=1;ii++){
 
276
        for(int jj=-1;jj<=1;jj++){
 
277
          int id2=idx(i+ii,j+jj);
 
278
          vector<pair<int,int> > bonds;
 
279
          if((id2==id) && (i!=0)){ // intra-cell, not for boundary 
 
280
            //std::cerr << id << " - " << id << std::endl;
 
281
            bonds=m_data[id].getBonds(gid,tol,ptag,mask);
 
282
          } else if(id2>id){ // inter-cell
 
283
            //std::cerr << id << " - " << id2 << std::endl;
 
284
            bonds=m_data[id].getBonds(gid,tol,m_data[id2],ptag,mask);
 
285
          }
 
286
          for(vector<pair<int,int> >::iterator iter=bonds.begin();
 
287
              iter!=bonds.end();
 
288
              iter++){
 
289
            //std::cerr << iter->first << " | " << iter->second << "   "; 
 
290
            m_bonds[btag].insert(*iter);
 
291
          }
 
292
          //std::cerr << std::endl;
 
293
        }
 
294
      }
 
295
      //std::cerr << std::endl;
 
296
    }
 
297
  }
 
298
}
 
299
 
 
300
/*!
 
301
  Output the content of a MNTable2D to an ostream. The output format depends on the value of MNTable2D::s_output_style (see  MNTable2D::SetOutputStyle). If it is set to 0, output suitable for debugging is generated, if it is set to 1 output in the esys .geo format is generated. If MNTable2D::s_output_style is set to 2, the output format is VTK-XML.
 
302
 
 
303
  \param ost the output stream
 
304
  \param T the table
 
305
*/
 
306
ostream& operator << (ostream& ost,const CircMNTable2D& T)
 
307
{
 
308
  if(CircMNTable2D::s_output_style==0){ // debug style
 
309
    // don't print padding
 
310
    MNTCell::SetOutputStyle(0);
 
311
    for(int i=0;i<T.m_nx;i++){
 
312
      for(int j=1;j<T.m_ny-1;j++){
 
313
        ost << "=== Cell " << i << " , " << j << " === " << endl;
 
314
        ost << T.m_data[T.idx(i,j)];
 
315
      }
 
316
    }
 
317
  } else if (CircMNTable2D::s_output_style==1){ // geo style
 
318
    int nparts=0;
 
319
    for(int i=1;i<T.m_nx-1;i++){
 
320
      for(int j=1;j<T.m_ny-1;j++){
 
321
        nparts+=T.m_data[T.idx(i,j)].NParts();
 
322
      }
 
323
    }
 
324
    // header
 
325
    ost << "LSMGeometry 1.2" << endl;
 
326
    ost << "BoundingBox " << T.m_x0+T.m_celldim << " " <<  T.m_y0+T.m_celldim << " 0.0 " << T.m_x0+double(T.m_nx-1)*T.m_celldim << " " <<  T.m_y0+double(T.m_ny-1)*T.m_celldim << " 0.0 " << endl;
 
327
    ost << "PeriodicBoundaries "  << T.m_x_periodic << " 0 0" << endl;
 
328
    ost << "Dimension 2D" << endl;
 
329
    // particles
 
330
    ost << "BeginParticles" << endl;
 
331
    ost << "Simple" << endl;
 
332
    ost << nparts << endl;
 
333
    MNTCell::SetOutputStyle(1);
 
334
     for(int i=1;i<T.m_nx-1;i++){
 
335
      for(int j=1;j<T.m_ny-1;j++){
 
336
        ost <<  T.m_data[T.idx(i,j)];
 
337
      }
 
338
    }
 
339
    ost << "EndParticles" << endl;
 
340
    // bonds
 
341
    ost << "BeginConnect" << endl;
 
342
    // sum bonds
 
343
    int nbonds=0;
 
344
    for(map<int,set<pair<int,int> > >::const_iterator iter=T.m_bonds.begin();
 
345
        iter!=T.m_bonds.end();
 
346
        iter++){
 
347
      nbonds+=iter->second.size();
 
348
    }
 
349
    ost << nbonds << endl;
 
350
    for(map<int,set<pair<int,int> > >::const_iterator iter=T.m_bonds.begin();
 
351
        iter!=T.m_bonds.end();
 
352
        iter++){
 
353
      for(set<pair<int,int> >::const_iterator v_iter=iter->second.begin();
 
354
          v_iter!=iter->second.end();
 
355
          v_iter++){
 
356
        ost << v_iter->first << " " << v_iter->second << " " << iter->first << endl;
 
357
      }
 
358
    }
 
359
    ost << "EndConnect" << endl;
 
360
  } else if (CircMNTable2D::s_output_style==2){ // vtk-xml
 
361
    T.WriteAsVtkXml(ost);
 
362
  }
 
363
 
 
364
  return ost;
 
365
}