NAME
FieldLattice - finite difference scalar field
INCLUDE
include "FieldLattice.h"
SYNTAX
class FieldLattice : public FieldWithPtValues
{
// Pointwise "finite difference" scalar field over a uniform regular grid
int nsd;
Handle(GridLattice) mesh; // finite difference grid (const. partition)
Handle(ArrayGenSel(NUMT)) vec; // used for point values
int fict_boundary;
bool reallocate
(
const GridLattice& grid,
const ArrayGen(NUMT)* coefficients =NULL,
const char* fieldname = NULL
);
// this one allows staggered grid and ghost boundary:
bool reallocate
(
GridLattice& grid,
const ArrayGenSel(NUMT)* coefficients = NULL,
GridLattice::Component comp = GridLattice:: none,
int fict_boundary = 0,
const char* fieldname = NULL,
const bool change_staggered_grid = true
);
public:
FieldLattice ();
// constructors are defined in terms of the redim functions below
FieldLattice (const GridLattice& grid,
const char* fieldname);
FieldLattice (const GridLattice& grid, ArrayGen(NUMT)& coefficients,
const char* fieldname);
FieldLattice (const GridLattice& grid, ArrayGenSel(NUMT)& coefficients,
const char* fieldname);
FieldLattice (const GridLattice& grid, const int fict_boundary,
const char* fieldname);
FieldLattice (const GridLattice& grid, GridLattice::Component comp,
const char* fieldname, const bool change_staggered_grid=true);
FieldLattice (const GridLattice& grid, GridLattice::Component comp,
const int fict_boundary, const char* fieldname,
const bool change_staggered_grid = true);
~FieldLattice ();
bool redim (const GridLattice& grid, const char* fieldname);
// memory is allocated
bool redim (const GridLattice& grid, const ArrayGen(NUMT)& coefficients,
const char* fieldname); // copies coefficients into vec
bool redim (const GridLattice& grid, const ArrayGenSel(NUMT)& coeff,
const char* fieldname); // vec points to coeff
bool redim (const GridLattice& grid, const int fict_boundary,
const char* fieldname); // allows a ghost boundary
bool redim (const GridLattice& grid, GridLattice::Component comp,
const char* fieldname,
const bool change_staggered_grid = true);
// allows staggered grid
bool redim (const GridLattice& grid, GridLattice::Component comp,
const int fict_boundary, const char* fieldname,
const bool change_staggered_grid = true);
// staggered grid and ghost boundary
bool redim (const FieldLattice& f, const char* fieldname);
GridLattice& grid () { return mesh.getRef(); }
const GridLattice& grid () const { return mesh.getRef(); }
ArrayGenSel(NUMT)& values () { return vec.getRef(); }
const ArrayGenSel(NUMT)& values () const { return vec.getRef(); }
ArrayGen(NUMT)& values0 (); // not inline!
const ArrayGen(NUMT)& values0 () const; // not inline!
bool ghostBoundary () const;
void exchangeValues (FieldLattice& f); // this <-> f
bool empty () const;
virtual bool ok () const;
virtual void minmax (NUMT& min, NUMT& max, GridWithPts* grid =NULL) const;
// The valueIndex(i,j) functions are inefficient and gives a warning,
// but are kept for backward compatibility. Use values()(i,j) instead!
NUMT& valueIndex (int i);
NUMT& valueIndex (int i, int j);
NUMT& valueIndex (int i, int j, int k);
NUMT valueIndex (int i) const;
NUMT valueIndex (int i, int j) const;
NUMT valueIndex (int i, int j, int k) const;
// convenient subscripting for staggered grids
// (these are never inlined; just for rapid prototyping!)
NUMT& valueIndex (double i);
NUMT& valueIndex (double i, double j);
NUMT& valueIndex (double i, double j, double k);
NUMT valueIndex (double i) const;
NUMT valueIndex (double i, double j) const;
NUMT valueIndex (double i, double j, double k) const;
int fictBoundary () const { return fict_boundary; }
// subscripting with an n-tuple index:
NUMT& valueIndex (const Ptv(int)& index);
NUMT valueIndex (const Ptv(int)& index) const;
void value (const ArrayGenSel(NUMT)& values);
virtual NUMT valuePt (const Ptv(real)& x, real t = DUMMY);
virtual Ptv(NUMT) derivativePt (const Ptv(real)& x, real t = DUMMY);
virtual NUMT valueFEM (const FiniteElement& fe, real t = DUMMY);
virtual void derivativeFEM (Ptv(NUMT)& d, const FiniteElement& fe,
real t = DUMMY);
void operator = (const FieldLattice& field);
void operator = (const FieldFunc& func);
void fill (const FieldFunc& func, real time = DUMMY);
void copyWithoutGhostBoundary (const FieldLattice& f);
virtual void print (Os os) const;
void scan (Is is);
virtual void fill (NUMT value);
virtual void add (NUMT value);
virtual void mult (NUMT value);
virtual void apply (Func(NUMT) f);
virtual void add (Field& field, int power, NUMT front_factor);
virtual Field& scale ();
virtual Field& unscale ();
virtual void unloadData (Os os) const;
virtual void loadData (Is is);
virtual void attach (Grid& grid);
virtual Grid* getGridBase();
// virtual functions from FieldWithPtValues:
virtual int getNoPoints () const // cannot be made inline
{ return mesh->getNoPoints(); }
virtual int getNoValues () const // cannot be made inline
{ return vec->size(); }
virtual NUMT& valuePoint (int point_no);
virtual NUMT valuePoint (int point_no) const;
virtual Ptv(real) getPt (int point_no) const
{ return mesh->getPt(mesh->single2multiple(point_no)); }
virtual GridWithPts& getGridWithPts ()
{ return mesh.getRef(); }
virtual const GridWithPts& getGridWithPts () const
{ return mesh.getRef(); }
virtual Vec(NUMT)& valuesVec ()
{ return vec.getRef(); }
virtual const Vec(NUMT)& valuesVec () const
{ return vec.getRef(); }
CLASS_INFO
VIRTUAL_CAST(FieldLattice)
// temporary:
COPY_CONSTRUCTOR(FieldLattice);
void pack (CharPack& package);
void unpack (CharPack& package, int package_size);
};
KEYWORDS
finite differences, scalar field
DESCRIPTION
The class implements a scalar field suited for finite difference
methods. The field is defined in terms of a GridLattice object
and an array containing the point values at the grid points.
The number of space dimensions can be arbitrary (not limited to
1D, 2D or 3D). However, in high dimensions the amount of data is
significant and subscripting is performed in terms of Ptv(int)
objects, which is not very efficient (it is then wise to design a
class with particular emphasis on efficient data handling).
This version of the FieldLattice class can operate on staggered
grids, with or without a ghost boundary.
If you use this FieldLattice heavily in finite difference pro
gramming, we strongly recommend to study the code in the member
functions, and the ArrayGenSel object and its base classes, to
get a complete knowledge of the tool. This is important for cre
ating efficient code.
CONSTRUCTORS AND INITIALIZATION
As for the other field objects, "initialization" here refers to
allocation of internal data structures for the field values and
setting a pointer to a valid grid object. Hence, initialization
of field values is not implied.
There are several types of redim functions and corresponding con
structors (constructors and redim functions call reallocate in
the same way).
(1) The constructor function without arguments performs no
actions. redim must be called later to allocate data structures
and bind the object to a grid.
(2) The redim function that takes a GridLattice object allocates
a corresponding field, according to the dimensions and lower and
upper values of the various indices in the grid.
(3) The redim function that takes a GridLattice object and an
ArrayGenSel(NUMT) object binds the mesh pointer to the grid
object and the vec pointer to the supplied array. Hence, there is
no allocation of data and the field object relies on external
data structures. This functionality is well suited for cases
where several fields can share the same fundamental array of
field values (transforming one field format to another is an
example). Note that when this constructor is used, the field val
ues are shared with the incoming array. That is, if that array is
changed outside the "FieldLattice object, the changes also affect
the FieldLattice object. Moreover, the constructor performs ini
tialization of field values in addition to the "initialization as
defined above.
(4) As point (3) above, but the array is of ArrayGen(NUMT) type.
In this case, the vec pointer cannot simply be bound to the sup
plied array. Instead, we allocate an ArrayGenSel(NUMT) array
(without any ghost boundary of course) and copy the values from
the supplied array. Notice that this seems to be a waste of mem
ory, but there is no way to "extend" a base class object (here
ArrayGen) to a subclass object (here ArrayGenSel).
(5) Another redim function takes a GridLattice object and an
integer representing the size of the fictitious (ghost) boundary.
The fict_boundary does not affect the GridLattice object, but the
array allocated can in addition to holding the values at the grid
points, also hold one or more extra points of points around grid.
This is often convenient when programming finite difference meth
ods. If the base value of the grid is i0, and the field was
allocated with fict_boundary=1, one can access and place values
in the fictitious point i0-1.
(6) A redim function takes a GridLattice object, an enumerator
comp_in_staggered_grid and a bool change_staggered_grid. The
enumerator is defined in class GridLattice and must have one of
the values GridLattice::p, GridLattice::u, GridLattice::v or
GridLattice::w, where p is the scalar value in the computations,
and u, v and w hold the values of the different velocity compo
nents. The underlying grid will now not be GridLattice, but one
of its subclasses GridLatticeB and GridLatticeC. See the documen
tation in these classes for more information on using the stag
gered grid utilities. The optional last argument change_stag
gered_grid is true by default, but might be set to false if the
staggered grid has been constructed in its right size before the
call to this redim or constructor. This functionality is normally
only used in parallel programming, where a manager is in charge
of dividing global grids into local ones and distributing these
to subdomain solvers.
(7) A redim function takes a GridLattice grid, an enumerator
comp_in_staggered_grid, an integer representing the size of the
fictitious boundary, and a bool variable change_staggered_grid.
The constructor is just a combination of the constructors (5) and
(6), see their documentation for more information.
MEMBER FUNCTIONS
redim - documented under "Construction and initialization".
grid - gives access to the underlying grid.
values - gives access to the array of field values (at the grid
points).
values0 - gives access to the array of field values (at the grid
points) as an ArrayGen(NUMT) array. This can have undesired side
effects if the underlying ArrayGenSel(NUMT) object, used in this
field, has a ghost boundary. You can use ghostBoundary() to check
if there are fictitious pints outside the grid. Notice that the
values0 function is not inline. One should therefore never use
it for subscripting line ...values0()(i,j). Use values0 to ini
tialize an ArrayGen(NUMT)& reference and perform subscripting
based on this reference!
valueIndex - gives access to a single field value. For example,
valueIndex(j,k) is identical to values()(j,k), but the former
construction is the "most protected" and hence the recommended
one. If one creates a loop over the field values, the loop
indices are given by the grid object, cf. the example below.
value - assigns numbers to all the field values based on the
array given as argument to the function. The internal data struc
tures in *this field must be allocated before the function can be
called. A copy of the array given as argument is taken.
The other functions are virtual and documented in the classes
FieldWithPtValues and Field.
EXAMPLES
First we show how a loop over the indices is constructed.
void exampleOnLoop (FieldLattice& f) {
FieldLattice f2; f2.redim(f); // f2 has same grid and size as
f
GridLattice& g = f.grid();
if (g.getNoSpaceDim() != 2) errorFP("exampleOnLoop","grid is
not 2D!");
int i0 = g.getBase(1); int j0 = g.getBase(2); // loop start
int in = g.getMaxI(1); int jn = g.getMaxI(2); // loop limits
int i,j;
for (i = i0; i <= in; i++)
for (j = j0; j <= jn; j++) {
f2(i,j) = 2*f(i,j) + 1; }
Now, this loop can be written more compactly using other member
functions. In the more compact notation, it is also easy to
write code that is independent of the number of space dimensions.
Here is an example:
void exampleOnLoop (FieldLattice& f) {
FieldLattice f2; f2.redim(f); // f2 has same grid and size as
f
f2 = f; // f2's values are a copy of f's
values
f2.mult(2.0); f2.add(1.0);
Next we discuss how f2 can be set equal to f. Here are some
alternatives:
FieldLattice f2; f2.redim(f); // f2 has same grid and size as
f
f2.value (f.values()); // f2 takes a copy of f.values
FieldLattice f3 (f.grid(), f.values(), "f3"); // f3 refer
ences f.values()
FieldLattice f4; f4.redim (f.grid(), f.values(), "f4"); //
reference
For f, f3 and f4, the arrays of field values in the various
fields are the same physical array. Therefore, if we set f3.val
ueIndex(1,3)=4.0 this affects the f and f4 fields as well. To
ensure that a field has its own local copy of the array values,
use the constructors or redim functions that takes the grid
object and not the array of field values as arguments, or use
values as we did above for f2.
SEE ALSO
class Field, class GridLattice, class ArrayGenSel(NUMT)
DEVELOPED BY
SINTEF Applied Mathematics, Oslo, Norway, and University of Oslo,
Dept. of Mathematics, Norway
AUTHOR
Hans Petter Langtangen, SINTEF/UiO extensions by Elizabeth Ack
lam, UiO