← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/ode/Partition.cpp

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/kernel/ode/Partition.cpp
1
 
// Copyright (C) 2004-2008 Johan Jansson and Anders Logg.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// First added:  2004
5
 
// Last changed: 2008-12-11
6
 
 
7
 
#include <algorithm>
8
 
#include <cmath>
9
 
 
10
 
#include <dolfin/dolfin_log.h>
11
 
#include <dolfin/dolfin_math.h>
12
 
#include <dolfin/parameters.h>
13
 
#include <dolfin/MultiAdaptivity.h>
14
 
#include <dolfin/Partition.h>
15
 
 
16
 
using namespace dolfin;
17
 
 
18
 
//-----------------------------------------------------------------------------
19
 
Partition::Partition(uint N) : indices(N)
20
 
{
21
 
  // Get parameter for threshold
22
 
  threshold = dolfin_get("ODE partitioning threshold");
23
 
 
24
 
  // Reset all indices
25
 
  for (uint i = 0; i < N; i++)
26
 
    indices[i] = i;
27
 
}
28
 
//-----------------------------------------------------------------------------
29
 
Partition::~Partition()
30
 
{
31
 
  // Do nothing
32
 
}
33
 
//-----------------------------------------------------------------------------
34
 
dolfin::uint Partition::size() const
35
 
{
36
 
  return indices.size();
37
 
}
38
 
//-----------------------------------------------------------------------------
39
 
dolfin::uint Partition::index(uint i) const
40
 
{
41
 
  dolfin_assert(i < indices.size());  
42
 
  return indices[i];
43
 
}
44
 
//-----------------------------------------------------------------------------
45
 
real Partition::update(uint offset, uint& end, MultiAdaptivity& adaptivity,
46
 
                       real K)
47
 
{
48
 
  // Compute time steps for partition. We partition the components into two
49
 
  // groups, one group with k < Kpivot and one with k >= Kpivot.
50
 
 
51
 
  // Compute time step for partitioning
52
 
  real Kpivot = threshold * maximum(offset, adaptivity);
53
 
 
54
 
  // Comparison operator
55
 
  Less less(Kpivot, adaptivity);
56
 
 
57
 
  // Partition using std::partition
58
 
  Array<uint>::iterator start = indices.begin();
59
 
  std::advance(start, offset);
60
 
  Array<uint>::iterator middle = std::partition(start, indices.end(), less);
61
 
  
62
 
  // Compute pivot index
63
 
  end = std::distance(indices.begin(), middle);
64
 
 
65
 
  // Modify time step to the smallest k such that k >= Kpivot.
66
 
  Kpivot = minimum(offset, end, adaptivity);
67
 
 
68
 
  // Modify time step so K is a multiple of Kpivot
69
 
  //Kpivot = K / ceil(K / Kpivot);
70
 
 
71
 
  /*
72
 
  for (uint i = offset; i < end ; i++)
73
 
  {
74
 
    uint index = indices[i];
75
 
    cout << "i = " << index << ": " << adaptivity.timestep(index)
76
 
         << " --> " << Kpivot << endl;
77
 
  }
78
 
  */
79
 
 
80
 
  return Kpivot;
81
 
}
82
 
//-----------------------------------------------------------------------------
83
 
void Partition::debug(uint offset, uint end) const
84
 
{
85
 
  cout << "Partition:";
86
 
 
87
 
  for (uint i = 0; i < indices.size(); i++)
88
 
  {
89
 
    if ( i == offset || i == end )
90
 
      cout << " |";
91
 
 
92
 
    cout << " " << indices[i];
93
 
  }
94
 
 
95
 
  cout << endl;
96
 
}
97
 
//-----------------------------------------------------------------------------
98
 
real Partition::maximum(uint offset, MultiAdaptivity& adaptivity) const
99
 
{
100
 
  real K = 0.0;
101
 
 
102
 
  for (uint i = offset; i < indices.size(); i++)
103
 
    K = std::max(adaptivity.timestep(indices[i]), K);
104
 
 
105
 
  return K;
106
 
}
107
 
//-----------------------------------------------------------------------------
108
 
real Partition::minimum(uint offset, uint end,
109
 
                           MultiAdaptivity& adaptivity) const
110
 
{
111
 
  real k = adaptivity.timestep(indices[offset]);
112
 
 
113
 
  for (uint i = offset + 1; i < end; i++)
114
 
    k = std::min(adaptivity.timestep(indices[i]), k);
115
 
 
116
 
  return k;
117
 
}
118
 
//-----------------------------------------------------------------------------
119
 
Partition::Less::Less(real& K, MultiAdaptivity& adaptivity)
120
 
  : K(K), adaptivity(adaptivity)
121
 
{
122
 
  // Do nothing
123
 
}
124
 
//-----------------------------------------------------------------------------
125
 
bool Partition::Less::operator()(uint index) const
126
 
{
127
 
  return adaptivity.timestep(index) >= K;
128
 
}
129
 
//-----------------------------------------------------------------------------