← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/function/DiscreteFunction.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/function/DiscreteFunction.cpp
1
 
// Copyright (C) 2007-2008 Anders Logg.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// Modified by Garth N. Wells, 2007.
5
 
//
6
 
// First added:  2007-04-02
7
 
// Last changed: 2008-01-03
8
 
 
9
 
#include <dolfin/dolfin_log.h>
10
 
#include <dolfin/Mesh.h>
11
 
#include <dolfin/Form.h>
12
 
#include <dolfin/DofMap.h>
13
 
#include <dolfin/Vertex.h>
14
 
#include <dolfin/Cell.h>
15
 
#include <dolfin/UFCMesh.h>
16
 
#include <dolfin/UFCCell.h>
17
 
#include <dolfin/ElementLibrary.h>
18
 
#include <dolfin/SubFunction.h>
19
 
#include <dolfin/SubSystem.h>
20
 
#include <dolfin/DiscreteFunction.h>
21
 
 
22
 
using namespace dolfin;
23
 
 
24
 
//-----------------------------------------------------------------------------
25
 
DiscreteFunction::DiscreteFunction(Mesh& mesh, Vector& x, Form& form, uint i)
26
 
  : GenericFunction(mesh),
27
 
    x(&x), finite_element(0), dof_map(0), dofs(0),
28
 
    local_vector(0), local_dof_map(0)
29
 
{
30
 
  // Update dof maps
31
 
  form.updateDofMaps(mesh);
32
 
  dof_map = &form.dofMaps()[i];
33
 
 
34
 
  // Initialise function
35
 
  init(mesh, x, form.form(), i);
36
 
}
37
 
//-----------------------------------------------------------------------------
38
 
DiscreteFunction::DiscreteFunction(Mesh& mesh, Vector& x, DofMap& dof_map,
39
 
                                   const ufc::form& form, uint i)
40
 
  : GenericFunction(mesh),
41
 
    x(&x), finite_element(0), dof_map(&dof_map), dofs(0),
42
 
    local_vector(0), local_dof_map(0)
43
 
{
44
 
  init(mesh, x, form, i);
45
 
}
46
 
//-----------------------------------------------------------------------------
47
 
DiscreteFunction::DiscreteFunction(Mesh& mesh, Vector& x,
48
 
                                   std::string finite_element_signature,
49
 
                                   std::string dof_map_signature)
50
 
  : GenericFunction(mesh),
51
 
    x(&x), finite_element(0), dof_map(0), dofs(0),
52
 
    local_vector(0), local_dof_map(0)
53
 
{
54
 
  // Create finite element
55
 
  finite_element = ElementLibrary::create_finite_element(finite_element_signature);
56
 
  if (!finite_element)
57
 
  {
58
 
    error("Unable to find finite element in library: \"%s\".",
59
 
                  finite_element_signature.c_str());
60
 
  }
61
 
 
62
 
  // Create dof map from signature
63
 
  dof_map = new DofMap(dof_map_signature, mesh);
64
 
 
65
 
  // Check size of vector
66
 
  if (x.size() != dof_map->global_dimension())
67
 
    error("Size of vector does not match global dimension of finite element space.");
68
 
 
69
 
  // Initialize local array for mapping of dofs
70
 
  dofs = new uint[dof_map->local_dimension()];
71
 
  for (uint i = 0; i < dof_map->local_dimension(); i++)
72
 
    dofs[i] = 0;
73
 
 
74
 
  // Assume responsibility for data
75
 
  local_dof_map = dof_map;
76
 
}
77
 
//-----------------------------------------------------------------------------
78
 
DiscreteFunction::DiscreteFunction(SubFunction& sub_function)
79
 
  : GenericFunction(sub_function.f->mesh),
80
 
    x(0), finite_element(0), dof_map(0), dofs(0),
81
 
    local_vector(0), local_dof_map(0)
82
 
{
83
 
  // Create sub system
84
 
  SubSystem sub_system(sub_function.i);
85
 
 
86
 
  // Extract sub element
87
 
  finite_element = sub_system.extractFiniteElement(*sub_function.f->finite_element);
88
 
 
89
 
  // Extract sub dof map and offset
90
 
  uint offset = 0;
91
 
  dof_map = sub_function.f->dof_map->extractDofMap(sub_system.array(), offset);
92
 
 
93
 
  // Create vector of dofs and copy values
94
 
  const uint n = dof_map->global_dimension();
95
 
  x = new Vector(n);
96
 
  real* values = new real[n];
97
 
  uint* get_rows = new uint[n];
98
 
  uint* set_rows = new uint[n];
99
 
  for (uint i = 0; i < n; i++)
100
 
  {
101
 
    get_rows[i] = offset + i;
102
 
    set_rows[i] = i;
103
 
  }
104
 
  sub_function.f->x->get(values, n, get_rows);
105
 
  x->set(values, n, set_rows);
106
 
  x->apply();
107
 
  delete [] values;
108
 
  delete [] get_rows;
109
 
  delete [] set_rows;
110
 
 
111
 
  // Initialize local array for mapping of dofs
112
 
  dofs = new uint[dof_map->local_dimension()];
113
 
  for (uint i = 0; i < dof_map->local_dimension(); i++)
114
 
    dofs[i] = 0;
115
 
 
116
 
  // Assume responsibility for vector and dof map
117
 
  local_vector = x;
118
 
  local_dof_map = dof_map;
119
 
}
120
 
//-----------------------------------------------------------------------------
121
 
DiscreteFunction::DiscreteFunction(const DiscreteFunction& f)
122
 
  : GenericFunction(f.mesh),
123
 
    x(0), finite_element(0), dof_map(0), dofs(0),
124
 
    local_vector(0), local_dof_map(0)
125
 
{
126
 
  cout << "Copy constructor for discrete function" << endl;
127
 
 
128
 
  // FIXME: Why don't we just copy the finite_element?
129
 
  // Create finite element
130
 
  finite_element = ElementLibrary::create_finite_element(f.finite_element->signature());
131
 
  if (!finite_element)
132
 
    error("Unable to find finite element in library: \"%s\".",
133
 
                  f.finite_element->signature());
134
 
 
135
 
  // Create dof map
136
 
  dof_map = new DofMap(f.dof_map->signature(), mesh); 
137
 
 
138
 
  // Create vector and copy values
139
 
  x  = new Vector(dof_map->global_dimension());
140
 
  *x = *f.x;
141
 
 
142
 
  // Initialize local array for mapping of dofs
143
 
  dofs = new uint[dof_map->local_dimension()];
144
 
  for (uint i = 0; i < dof_map->local_dimension(); i++)
145
 
    dofs[i] = 0;
146
 
 
147
 
  // Assume responsibility for vector and dof map
148
 
  local_vector = x;
149
 
  local_dof_map = dof_map;
150
 
}
151
 
//-----------------------------------------------------------------------------
152
 
DiscreteFunction::~DiscreteFunction()
153
 
{
154
 
  if (finite_element)
155
 
    delete finite_element;
156
 
      
157
 
  if (dofs)
158
 
    delete [] dofs;
159
 
 
160
 
  if (local_vector)
161
 
    delete local_vector;
162
 
 
163
 
  if (local_dof_map)
164
 
    delete local_dof_map;
165
 
}
166
 
//-----------------------------------------------------------------------------
167
 
dolfin::uint DiscreteFunction::rank() const
168
 
{
169
 
  dolfin_assert(finite_element);
170
 
  return finite_element->value_rank();
171
 
}
172
 
//-----------------------------------------------------------------------------
173
 
dolfin::uint DiscreteFunction::dim(uint i) const
174
 
{
175
 
  dolfin_assert(finite_element);
176
 
  return finite_element->value_dimension(i);
177
 
}
178
 
//-----------------------------------------------------------------------------
179
 
dolfin::uint DiscreteFunction::numSubFunctions() const
180
 
{
181
 
  dolfin_assert(finite_element);
182
 
  return finite_element->num_sub_elements();
183
 
}
184
 
//-----------------------------------------------------------------------------
185
 
const DiscreteFunction& DiscreteFunction::operator= (const DiscreteFunction& f)
186
 
{
187
 
  cout << "Assignment in discrete function" << endl;
188
 
 
189
 
  // Check that data matches
190
 
  if (strcmp(finite_element->signature(), f.finite_element->signature()) != 0 ||
191
 
       strcmp(dof_map->signature(), f.dof_map->signature()) != 0              ||
192
 
      x->size() != f.x->size())
193
 
  {
194
 
    error("Assignment of discrete function failed. Finite element spaces or dimensions don't match.");
195
 
  }
196
 
 
197
 
  // Copy vector
198
 
  *x = *f.x;
199
 
 
200
 
  return *this;
201
 
}
202
 
//-----------------------------------------------------------------------------
203
 
void DiscreteFunction::interpolate(real* values)
204
 
{
205
 
  dolfin_assert(values);
206
 
  dolfin_assert(finite_element);
207
 
  dolfin_assert(dof_map);
208
 
  dolfin_assert(dofs);
209
 
  
210
 
  // Compute size of value (number of entries in tensor value)
211
 
  uint size = 1;
212
 
  for (uint i = 0; i < finite_element->value_rank(); i++)
213
 
    size *= finite_element->value_dimension(i);
214
 
 
215
 
  // Local data for interpolation on each cell
216
 
  CellIterator cell(mesh);
217
 
  UFCCell ufc_cell(*cell);
218
 
  const uint num_cell_vertices = mesh.type().numVertices(mesh.topology().dim());
219
 
  real* vertex_values = new real[size*num_cell_vertices];
220
 
  real* dof_values = new real[finite_element->space_dimension()];
221
 
 
222
 
  // Interpolate vertex values on each cell and pick the last value
223
 
  // if two or more cells disagree on the vertex values
224
 
  for (; !cell.end(); ++cell)
225
 
  {
226
 
    // Update to current cell
227
 
    ufc_cell.update(*cell);
228
 
 
229
 
    // Tabulate dofs
230
 
    dof_map->tabulate_dofs(dofs, ufc_cell);
231
 
    
232
 
    // Pick values from global vector
233
 
    x->get(dof_values, dof_map->local_dimension(), dofs);
234
 
 
235
 
    // Interpolate values at the vertices
236
 
    finite_element->interpolate_vertex_values(vertex_values, dof_values, ufc_cell);
237
 
 
238
 
    // Copy values to array of vertex values
239
 
    for (VertexIterator vertex(*cell); !vertex.end(); ++vertex)
240
 
      for (uint i = 0; i < size; ++i)
241
 
        values[i*mesh.numVertices() + vertex->index()] = vertex_values[i*num_cell_vertices + vertex.pos()];
242
 
  }
243
 
 
244
 
  // Delete local data
245
 
  delete [] vertex_values;
246
 
  delete [] dof_values;
247
 
}
248
 
//-----------------------------------------------------------------------------
249
 
void DiscreteFunction::interpolate(real* coefficients,
250
 
                                   const ufc::cell& cell,
251
 
                                   const ufc::finite_element& finite_element)
252
 
{
253
 
  dolfin_assert(coefficients);
254
 
  dolfin_assert(this->finite_element);
255
 
  dolfin_assert(this->dof_map);
256
 
  dolfin_assert(this->dofs);
257
 
 
258
 
  // FIXME: Better test here, compare against the local element
259
 
 
260
 
  // Check dimension
261
 
  if (finite_element.space_dimension() != dof_map->local_dimension())
262
 
    error("Finite element does not match for interpolation of discrete function.");
263
 
 
264
 
  // Tabulate dofs
265
 
  dof_map->tabulate_dofs(dofs, cell);
266
 
  
267
 
  // Pick values from global vector
268
 
  x->get(coefficients, dof_map->local_dimension(), dofs);
269
 
}
270
 
//-----------------------------------------------------------------------------
271
 
Vector& DiscreteFunction::vector() const
272
 
{
273
 
  if( !x )
274
 
    error("Vector associated with DiscreteFunction has not been initialised.");
275
 
 
276
 
  return *x;
277
 
}
278
 
//-----------------------------------------------------------------------------
279
 
void DiscreteFunction::init(Mesh& mesh, Vector& x, const ufc::form& form, uint i)
280
 
{
281
 
  // Check argument
282
 
  const uint num_arguments = form.rank() + form.num_coefficients();
283
 
  if (i >= num_arguments)
284
 
    error("Illegal function index %d. Form only has %d arguments.",
285
 
                  i, num_arguments);
286
 
 
287
 
  // Create finite element
288
 
  finite_element = form.create_finite_element(i);
289
 
 
290
 
  // Initialize vector
291
 
  if (x.size() != dof_map->global_dimension())
292
 
    x.init(dof_map->global_dimension());
293
 
 
294
 
  // Initialize local array for mapping of dofs
295
 
  dofs = new uint[dof_map->local_dimension()];
296
 
  for (uint i = 0; i < dof_map->local_dimension(); i++)
297
 
    dofs[i] = 0;
298
 
}
299
 
//-----------------------------------------------------------------------------