1
// Copyright (C) 2007-2008 Anders Logg.
2
// Licensed under the GNU LGPL Version 2.1.
4
// Modified by Garth N. Wells, 2007.
6
// First added: 2007-04-02
7
// Last changed: 2008-01-03
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>
22
using namespace dolfin;
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)
31
form.updateDofMaps(mesh);
32
dof_map = &form.dofMaps()[i];
34
// Initialise function
35
init(mesh, x, form.form(), i);
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)
44
init(mesh, x, form, i);
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)
54
// Create finite element
55
finite_element = ElementLibrary::create_finite_element(finite_element_signature);
58
error("Unable to find finite element in library: \"%s\".",
59
finite_element_signature.c_str());
62
// Create dof map from signature
63
dof_map = new DofMap(dof_map_signature, mesh);
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.");
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++)
74
// Assume responsibility for data
75
local_dof_map = dof_map;
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)
84
SubSystem sub_system(sub_function.i);
86
// Extract sub element
87
finite_element = sub_system.extractFiniteElement(*sub_function.f->finite_element);
89
// Extract sub dof map and offset
91
dof_map = sub_function.f->dof_map->extractDofMap(sub_system.array(), offset);
93
// Create vector of dofs and copy values
94
const uint n = dof_map->global_dimension();
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++)
101
get_rows[i] = offset + i;
104
sub_function.f->x->get(values, n, get_rows);
105
x->set(values, n, set_rows);
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++)
116
// Assume responsibility for vector and dof map
118
local_dof_map = dof_map;
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)
126
cout << "Copy constructor for discrete function" << endl;
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());
132
error("Unable to find finite element in library: \"%s\".",
133
f.finite_element->signature());
136
dof_map = new DofMap(f.dof_map->signature(), mesh);
138
// Create vector and copy values
139
x = new Vector(dof_map->global_dimension());
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++)
147
// Assume responsibility for vector and dof map
149
local_dof_map = dof_map;
151
//-----------------------------------------------------------------------------
152
DiscreteFunction::~DiscreteFunction()
155
delete finite_element;
164
delete local_dof_map;
166
//-----------------------------------------------------------------------------
167
dolfin::uint DiscreteFunction::rank() const
169
dolfin_assert(finite_element);
170
return finite_element->value_rank();
172
//-----------------------------------------------------------------------------
173
dolfin::uint DiscreteFunction::dim(uint i) const
175
dolfin_assert(finite_element);
176
return finite_element->value_dimension(i);
178
//-----------------------------------------------------------------------------
179
dolfin::uint DiscreteFunction::numSubFunctions() const
181
dolfin_assert(finite_element);
182
return finite_element->num_sub_elements();
184
//-----------------------------------------------------------------------------
185
const DiscreteFunction& DiscreteFunction::operator= (const DiscreteFunction& f)
187
cout << "Assignment in discrete function" << endl;
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())
194
error("Assignment of discrete function failed. Finite element spaces or dimensions don't match.");
202
//-----------------------------------------------------------------------------
203
void DiscreteFunction::interpolate(real* values)
205
dolfin_assert(values);
206
dolfin_assert(finite_element);
207
dolfin_assert(dof_map);
210
// Compute size of value (number of entries in tensor value)
212
for (uint i = 0; i < finite_element->value_rank(); i++)
213
size *= finite_element->value_dimension(i);
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()];
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)
226
// Update to current cell
227
ufc_cell.update(*cell);
230
dof_map->tabulate_dofs(dofs, ufc_cell);
232
// Pick values from global vector
233
x->get(dof_values, dof_map->local_dimension(), dofs);
235
// Interpolate values at the vertices
236
finite_element->interpolate_vertex_values(vertex_values, dof_values, ufc_cell);
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()];
245
delete [] vertex_values;
246
delete [] dof_values;
248
//-----------------------------------------------------------------------------
249
void DiscreteFunction::interpolate(real* coefficients,
250
const ufc::cell& cell,
251
const ufc::finite_element& finite_element)
253
dolfin_assert(coefficients);
254
dolfin_assert(this->finite_element);
255
dolfin_assert(this->dof_map);
256
dolfin_assert(this->dofs);
258
// FIXME: Better test here, compare against the local element
261
if (finite_element.space_dimension() != dof_map->local_dimension())
262
error("Finite element does not match for interpolation of discrete function.");
265
dof_map->tabulate_dofs(dofs, cell);
267
// Pick values from global vector
268
x->get(coefficients, dof_map->local_dimension(), dofs);
270
//-----------------------------------------------------------------------------
271
Vector& DiscreteFunction::vector() const
274
error("Vector associated with DiscreteFunction has not been initialised.");
278
//-----------------------------------------------------------------------------
279
void DiscreteFunction::init(Mesh& mesh, Vector& x, const ufc::form& form, uint i)
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.",
287
// Create finite element
288
finite_element = form.create_finite_element(i);
291
if (x.size() != dof_map->global_dimension())
292
x.init(dof_map->global_dimension());
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++)
299
//-----------------------------------------------------------------------------