← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/ode/Method.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/Method.cpp
1
 
// Copyright (C) 2003-2005 Anders Logg.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// Modified by Garth N. Wells, 2006.
5
 
//
6
 
// First added:  2003
7
 
// Last changed: 2006-03-27
8
 
 
9
 
#include <dolfin/dolfin_log.h>
10
 
#include <dolfin/dolfin_math.h>
11
 
#include <dolfin/Lagrange.h>
12
 
#include <dolfin/Method.h>
13
 
 
14
 
using namespace dolfin;
15
 
 
16
 
//-----------------------------------------------------------------------------
17
 
Method::Method(unsigned int q, unsigned int nq, unsigned int nn)
18
 
{
19
 
  dolfin_assert(nq > 0);
20
 
  dolfin_assert(nn > 0);
21
 
 
22
 
  this->q = q;
23
 
  this->nq = nq;
24
 
  this->nn = nn;
25
 
 
26
 
  // Allocate quadrature points
27
 
  qpoints  = new real[nq];
28
 
  for (unsigned int i = 0; i < nq; i++)
29
 
    qpoints[i] = 0.0;
30
 
 
31
 
  // Allocate nodal points
32
 
  npoints  = new real[nn];
33
 
  for (unsigned int i = 0; i < nn; i++)
34
 
    npoints[i] = 0.0;
35
 
 
36
 
  // Allocate quadrature weights
37
 
  qweights = new real[nq];
38
 
  for (unsigned int i = 0; i < nq; i++)
39
 
    qweights[i] = 0.0;
40
 
 
41
 
  // Allocate weights
42
 
  nweights = new real*[nn];
43
 
  for (unsigned int i = 0; i < nn; i++)
44
 
  {
45
 
    nweights[i] = new real[nq];
46
 
    for (unsigned int j = 0; j < nq; j++)
47
 
      nweights[i][j] = 0.0;
48
 
  }
49
 
 
50
 
  // Allocate derivatives
51
 
  derivatives = new real[nq];
52
 
  for (unsigned int i = 0; i < nq; i++)
53
 
    derivatives[i] = 0.0;
54
 
 
55
 
  trial = 0;
56
 
  test = 0;
57
 
 
58
 
  _type = none;
59
 
}
60
 
//-----------------------------------------------------------------------------
61
 
Method::~Method()
62
 
{
63
 
  if ( qpoints ) delete [] qpoints;
64
 
  if ( npoints ) delete [] npoints;
65
 
  if ( qweights ) delete [] qweights;
66
 
 
67
 
  if ( nweights )
68
 
  {
69
 
    for (unsigned int i = 0; i < nn; i++)
70
 
      delete [] nweights[i];
71
 
    delete [] nweights;
72
 
  }
73
 
 
74
 
  if ( derivatives ) delete [] derivatives;
75
 
 
76
 
  if ( trial ) delete trial;
77
 
  if ( test ) delete test;
78
 
}
79
 
//-----------------------------------------------------------------------------
80
 
void Method::init()
81
 
{
82
 
  computeQuadrature();
83
 
  computeBasis();
84
 
  computeWeights();
85
 
  computeDerivatives();
86
 
}
87
 
//-----------------------------------------------------------------------------
88
 
void Method::update(real x0, real f[], real k, real values[]) const
89
 
{
90
 
  // Update values
91
 
  for (uint i = 0; i < nn; i++)
92
 
  {
93
 
    real sum = 0.0;
94
 
    for (uint j = 0; j < nq; j++)
95
 
      sum += nweights[i][j] * f[j];
96
 
    values[i] = x0 + k*sum;
97
 
  }
98
 
}
99
 
//-----------------------------------------------------------------------------
100
 
void Method::update(real x0, real f[], real k, real values[], real alpha) const
101
 
{
102
 
  // Update values
103
 
  for (uint i = 0; i < nn; i++)
104
 
  {
105
 
    real sum = 0.0;
106
 
    for (uint j = 0; j < nq; j++)
107
 
      sum += nweights[i][j] * f[j];
108
 
    values[i] += alpha*(x0 + k*sum - values[i]);
109
 
  }
110
 
}
111
 
//-----------------------------------------------------------------------------
112
 
void Method::computeDerivatives()
113
 
{
114
 
  for (unsigned int i = 0; i < nq; i++)
115
 
    derivatives[i] = trial->ddx(i, 1.0);
116
 
}
117
 
//-----------------------------------------------------------------------------